Rotational sum rules

In this tutorial we will study the importance of rotational sum rules and how they can be enforced in hiphive. These sum rules can currently only be applied to the second order force constants.

The rotational sum rules are dependent on the translational sum rule being fulfilled, i.e.

\[\sum_i\Phi_{ij}^{\alpha\beta} = 0.\]

The first rotational sum rule we consider is the Born-Huang condition

\[\sum_j\Phi_{ij}^{\alpha\beta}r_j^\gamma = \sum_j\Phi_{ij}^{\alpha\gamma}r_j^\beta,\]

where \(r_j^\beta\) is the position vector of atom \(j\). Furthermore we consider the Huang invariances

\[\sum_{ij}\Phi_{ij}^{\alpha\beta}r_{ij}^\gamma r_{ij}^\delta = \sum_{ij}\Phi_{ij}^{\gamma\delta}r_{ij}^\alpha r_{ij}^\beta.\]

MoS2 monolayer

We will use a MoS2 monolayer to showcase the importance of correctly imposing rotational sum rules, which are necessary in order to correctly obtain a quadratic dispersion for the lowermost transverse acoustic mode in the vicinity of the \(\Gamma\)-point.

Forces from density functional theory (DFT) calculations are supplied in the MoS2_rattled_structures.extxyz file.

Phonon dispersion for MoS2 close to the gamma point with translational (T) sum rules as well as Born-Huang (BH) and Huang (H) invariances enforced.

It is evident from the figure that it is important to enforce both of the aforementioned rotational sum rules in order to obtain a correct dispersion.

Source code

Set up structure container
examples/advanced_topics/rotational_sum_rules/MoS2_monolayer/1_setup_structure_container.py

from ase.io import read
from hiphive import ClusterSpace, StructureContainer

# parameters
cutoffs = [8.0]
structures = read('MoS2_rattled_structures.extxyz@:')

# build clusterspace and structurecontainer
cs = ClusterSpace(structures[0], cutoffs)
sc = StructureContainer(cs)
for structure in structures:
    sc.add_structure(structure)
sc.write('structure_container.sc')

Construct force constant potentials while imposing different sum rules
examples/advanced_topics/rotational_sum_rules/MoS2_monolayer/2_construct_fcps.py

from hiphive import ForceConstantPotential, StructureContainer
from hiphive import enforce_rotational_sum_rules
from trainstation import Optimizer


# fit model
sc = StructureContainer.read('structure_container.sc')
cs = sc.cluster_space
opt = Optimizer(sc.get_fit_data(), train_size=1.0)
opt.train()
print(opt)

# apply sum rules
parameters = opt.parameters
parameters_huang = enforce_rotational_sum_rules(
    cs, parameters, ['Huang'])
parameters_bornhuang = enforce_rotational_sum_rules(
    cs, parameters, ['Born-Huang'])
parameters_rot = enforce_rotational_sum_rules(
    cs, parameters, ['Huang', 'Born-Huang'])

fcp_normal = ForceConstantPotential(cs, parameters)
fcp_huang = ForceConstantPotential(cs, parameters_huang)
fcp_bornhuang = ForceConstantPotential(cs, parameters_bornhuang)
fcp_rot = ForceConstantPotential(cs, parameters_rot)

fcp_normal.write('fcp_normal.fcp')
fcp_huang.write('fcp_huang.fcp')
fcp_bornhuang.write('fcp_bornhuang.fcp')
fcp_rot.write('fcp_rot.fcp')

Plot the phonon dispersion for each force constant potential
examples/advanced_topics/rotational_sum_rules/MoS2_monolayer/3_plot_dispersions.py

import numpy as np
import matplotlib.pyplot as plt
from ase import Atoms
from hiphive import ForceConstantPotential
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms


def get_band(q_start, q_stop, N=500):
    return np.array([q_start + (q_stop-q_start)*i/(N-1) for i in range(N)])


def compute_dispersion(prim, fcp):
    dim_matrix = np.diag([6, 6, 1])

    # setup phonopy
    atoms_phonopy = PhonopyAtoms(symbols=prim.get_chemical_symbols(),
                                 scaled_positions=prim.get_scaled_positions(),
                                 cell=prim.cell)
    phonopy = Phonopy(atoms_phonopy, supercell_matrix=dim_matrix,
                      primitive_matrix=None)

    # set force constants
    supercell = phonopy.supercell
    supercell = Atoms(cell=supercell.cell, numbers=supercell.numbers, pbc=True,
                      scaled_positions=supercell.scaled_positions)
    fcs = fcp.get_force_constants(supercell)
    phonopy.force_constants = fcs.get_fc_array(order=2)

    # get dispersion
    band1 = get_band(np.array([0, 0, 0]), np.array([0.5, 0, 0]))
    bands = [band1]

    # get phonon dispersion
    phonopy.run_band_structure(bands)
    result = phonopy.get_band_structure_dict()
    qnorms = result['distances'][0]
    freqs = result['frequencies'][0]
    return qnorms, freqs


# compute dispersions
fcp_normal = ForceConstantPotential.read('fcp_normal.fcp')
fcp_rot = ForceConstantPotential.read('fcp_rot.fcp')
fcp_huang = ForceConstantPotential.read('fcp_huang.fcp')
fcp_bornhuang = ForceConstantPotential.read('fcp_bornhuang.fcp')

prim = fcp_normal.primitive_structure

qvec_normal, freqs_normal = compute_dispersion(prim, fcp_normal)
qvec_rot, freqs_rot = compute_dispersion(prim, fcp_rot)
qvec_huang, freqs_huang = compute_dispersion(prim, fcp_huang)
qvec_bornhuang, freqs_bornhuang = compute_dispersion(prim, fcp_bornhuang)


# plot
lw = 3.0
fs = 14
fig = plt.figure(figsize=(7.0, 4.8))
ax1 = fig.add_subplot(1, 1, 1)

ax1.plot(qvec_normal, freqs_normal, '-', lw=lw, color='tab:blue')
ax1.plot(np.nan, np.nan, '-', lw=lw, color='tab:blue', label='T')


ax1.plot(qvec_huang, freqs_huang, '--', lw=lw, color='tab:red')
ax1.plot(np.nan, np.nan, '--', lw=lw, color='tab:red', label='T + H')

ax1.plot(qvec_bornhuang, freqs_bornhuang, '--', lw=lw, color='tab:green')
ax1.plot(np.nan, np.nan, '--', lw=lw, color='tab:green', label='T + BH')

ax1.plot(qvec_rot, freqs_rot, '--', lw=lw, color='tab:orange')
ax1.plot(np.nan, np.nan, '--', lw=lw, color='tab:orange', label='T + H + BH')

ax1.set_xlabel('Wave vector', fontsize=fs)
ax1.set_ylabel('Frequency (THz)', fontsize=fs)
ax1.tick_params(labelsize=fs)
ax1.legend(fontsize=fs)

# zoom on gamma point
ax1.set_xlim([0.0, 0.07])
ax1.set_ylim([-0.5, 3.0])

fig.tight_layout()
fig.savefig('MoS2_phonon_dispersion.svg')

The alpha parameter

This post processing method has a parameter \(\alpha\) which determines how much the parameters can change in order to enforce the sum rules. Below the phonon dispersion for various values of \(\alpha\) is shown.

The analysis shows that for small values of \(\alpha\) the rotational sum rules are correctly enforced as the parameters are allowed to change a lot. It is also apparent that the lowest optical branch is affected by the post-processing of the force constants. When \(\alpha\) increases the dispersion moves towards the dispersions without rotational sum rules.

Fitting with constraints

In the examples above the rotational constraints are imposed after the fitting/extraction of the parameters (force-constants) is completed. Another approach to enforcing rotational constraints is to include them during fitting. This task can be cast in an extended linear form

\[\left\Vert\boldsymbol{A}\boldsymbol{x} - \boldsymbol{f}\right\Vert^2_2 + \lambda \left\Vert\ \boldsymbol{A}_\text{rotational} \boldsymbol{x}\right\Vert^2_2\]

where the first term is the normal equation for minizing the error of the forces, \(\boldsymbol{A}_\text{rotational}\) is the rotational constraint matrix, and \(\lambda\) is the parameter that adjusts how strongly to enforce the rotational constraints. This integrated approach is advantageous to the post-processing approach as the enforcement of the rotational sum rules has a smaller effect on the overall dispersion as shown in the figure below. The quadratic dispersion is recovered and the lowest optical mode is not affected by this approach.

The source code for this analysis can be found in the hiphive-examples repository. The same scheme also allow one to enforce the translational sum rules as demonstrated by this example.

Graphene

Graphene is another 2D-material where the rotational sum rules can make a huge difference in the dispersion near the gamma point.

Functionality exists in hiphive to project existing force constants onto a ClusterSpace, as described here. Together with the functionality to enforce the rotational sum rules this opens up for the possibility to correct force constants calculated with external tools such as e.g., phonopoy.

Extracting the harmonic force constants with phonopy for graphene only requires one supercell calculation. This structure together with forces from density functional theory (DFT) calculations are supplied here in the graphene_phonopy_supercell.extxyz file.

Phonon dispersion for graphene close to the gamma point.

The dispersions from phonopy and the one parametrized using the extracted parameters closely match either other. Both dispersions exhibit, however, a non-quadratic imaginary pocket in the vicinity of the Gamma point. This pocket vanishes and the correct quadratic dispersion is obtained by enforcing the rotational sum rules.

Source code

graphene example examples/advanced_topics/rotational_sum_rules/graphene/graphene.py

import numpy as np
from matplotlib import pyplot as plt
from ase.io import read
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
from ase.lattice.hexagonal import Graphene

from hiphive import ForceConstants, ClusterSpace, ForceConstantPotential
from hiphive import enforce_rotational_sum_rules
from hiphive.utilities import extract_parameters


THz_to_meV = 4.13567


def get_band(q_start, q_stop, N):
    """ Return path between q_start and q_stop """
    return np.array([q_start + (q_stop-q_start)*i/(N-1) for i in range(N)])


def plot_dispersion(fcs, color, label):
    # set force constants
    phonopy.force_constants = fcs.get_fc_array(order=2)

    # get dispersion
    phonopy.run_band_structure(bands)
    result = phonopy.get_band_structure_dict()
    qnorms = np.hstack(result['distances'])
    freqs = THz_to_meV * np.vstack(result['frequencies'])
    lines = ax1.plot(qnorms, freqs, color=color)
    lines[0].set_label(label)


# parameters
cutoff = 8
dim = 7
mesh = [24] * 3
a0 = 2.466340583
vacuum = 40

# set up prim structure
prim = Graphene(symbol='C', latticeconstant={'a': a0, 'c': vacuum})

# define band path
N_band_pts = 500
G2K = get_band(np.array([0.0, 0.0, 0.0]), np.array([1/3, 1/3, 0.0]), N_band_pts)
bands = [G2K]

# init phonopy
phonopy_prim = PhonopyAtoms(numbers=prim.numbers, positions=prim.positions, cell=prim.cell)
phonopy = Phonopy(phonopy_prim, supercell_matrix=np.diag([dim, dim, 1]), primitive_matrix=None)
phonopy.generate_displacements()

# load the precalculated single phonopy supercell with displacements and forces
supercell = read('graphene_phonopy_supercell.extxyz')

# Add forces to phonopy and calculate force constants
phonopy.forces = [supercell.get_forces()]
phonopy.produce_force_constants()
fcs_phonopy = ForceConstants.from_arrays(supercell, phonopy.force_constants)

# map fcs_phonopy onto ClusterSpace and produce fcs_hiphive
cs = ClusterSpace(prim, [cutoff])
parameters = extract_parameters(fcs_phonopy, cs)
fcp = ForceConstantPotential(cs, parameters)
fcs_hiphive = fcp.get_force_constants(supercell)

# enforce rotational sum rules and produce fcs_hiphive_rot
enforced_parameters = enforce_rotational_sum_rules(cs, parameters, ['Huang', 'Born-Huang'])
fcp_rot = ForceConstantPotential(cs, enforced_parameters)
fcs_hiphive_rot = fcp_rot.get_force_constants(supercell)


# plotting
fig = plt.figure(figsize=(5.5, 3.8))
ax1 = fig.add_subplot(111)

ax1.axhline(y=0.0, ls='-', c='k', lw=1.0)
plot_dispersion(fcs_phonopy, 'tab:blue', 'phonopy raw')
plot_dispersion(fcs_hiphive, 'tab:orange', 'phonopy-based FCP')
plot_dispersion(fcs_hiphive_rot, 'tab:green', 'phonopy-based FCP rotational invariant')

# show
ax1.legend(loc=3)
ax1.set_xlim(0, 0.1)
ax1.set_ylim(-8, 15)

ax1.set_xlabel('Wave vector')
ax1.set_ylabel('Energy (meV)')

fig.tight_layout()
fig.savefig('graphene_phonon_dispersion.svg')