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

$\Phi_{ij}^{\alpha\beta}r_j^\gamma = \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.$

MoS:sub:2 monolayer¶

We will use a MoS:sub:2 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.

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

Source code¶

setup structure container
tutorial/advanced/rotational_sum_rules/1_setup_structure_container.py
from ase.io import read
from hiphive import ClusterSpace, StructureContainer

# parameters
cutoffs = [8.0]

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

construct fcps with different sum rules imposed
tutorial/advanced/rotational_sum_rules/2_construct_fcps.py
from hiphive import ForceConstantPotential, StructureContainer
from hiphive import enforce_rotational_sum_rules
from hiphive.fitting import Optimizer

# fit model
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 fcp
tutorial/advanced/rotational_sum_rules/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.get_supercell()
supercell = Atoms(cell=supercell.cell, numbers=supercell.numbers, pbc=True,
scaled_positions=supercell.get_scaled_positions())
fcs = fcp.get_force_constants(supercell)
phonopy.set_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.set_band_structure(bands)
qvecs, qnorms, freqs, _ = phonopy.get_band_structure()
return qnorms[0], freqs[0]

# compute dispersions

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.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')