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.\]

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 enfore both of the aforementioned rotational sum rules in order to obtain a correct dispersion.

Source code

setup structure container
examples/rotational_sum_rules/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 fcps with different sum rules imposed
examples/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
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 fcp
examples/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
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')