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.

../_images/MoS2_phonon_dispersion.svg

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/advanced_topics/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/advanced_topics/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/advanced_topics/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')

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 extract parameters from already calculated force constants. Together with the functionality to enforce the rotational sum rules this opens up for the possibility to correct force constants calculated with e.g. phonopoy.

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

../_images/graphene_phonon_dispersion.svg

Phonon dispersion for graphene close to the gamma point.

The dispersion calcualted using the extrated parameters closely match the phonopy dispersion except at the gamma point where both dispersions have a non-quadratic imaginary pocket. After the parameters have been symmetrized with the rotational sum rules the quadratic mode emerges.

As a bonus, the same procedure is performed but instead the parameters are fitted directly to the single displaced supercell given by phonopy.

Source code

graphene example examples/advanced_topics/rotational_sum_rules/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, StructureContainer
from hiphive.fitting import Optimizer
from hiphive import enforce_rotational_sum_rules
from hiphive.utilities import extract_parameters

# parameters
cutoff = 8
dim = 7
mesh = [24] * 3
a0 = 2.466340583      # Lattice parameter DFT - gamma 5x5x1 kpoints
vacuum = 40


def THz_to_meV(THz):
    return THz * 4.13567


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

# 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_supercell.extxyz')

# Add forces to phonopy and calcualte force constants
phonopy.set_forces([supercell.arrays['forces']])
phonopy.produce_force_constants()
fcs_phonopy = ForceConstants.from_arrays(supercell, phonopy.get_force_constants())

# Define band path
N_band_pts = 100
G2K = np.linspace([0.0, 0.0, 0.0], [1/3, 1/3, 0.0], N_band_pts)
K2M = np.linspace([1/3, 1/3, 0.0], [1/2, 0.0, 0.0], N_band_pts)
M2G = np.linspace([1/2, 0.0, 0.0], [0.0, 0.0, 0.0], N_band_pts)
bands = [G2K, K2M, M2G]


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

    # get dispersion
    phonopy.set_band_structure(bands)
    _, qnorms, freqs, _ = phonopy.get_band_structure()
    energies = THz_to_meV(np.array(freqs))
    qnorms = np.array(qnorms)

    lines = plt.plot(qnorms.ravel(), energies.reshape(-1, 6), color=color)
    lines[0].set_label(label)


# plot phonopy dispersion
plot_dispersion(fcs_phonopy, 'tab:blue', 'phonopy')

# get hiphive dispersion
cs = ClusterSpace(prim, [cutoff])
parameters = extract_parameters(fcs_phonopy, cs)
fcp = ForceConstantPotential(cs, parameters)
fcs_hiphive = fcp.get_force_constants(supercell)
plot_dispersion(fcs_hiphive, 'tab:orange', 'hiphive')

# plot rotational dispersion
enforced_parameters = enforce_rotational_sum_rules(cs, parameters, ['Huang', 'Born-Huang'])
fcp_enforced = ForceConstantPotential(cs, enforced_parameters)
fcs_hiphive_enforced = fcp_enforced.get_force_constants(supercell)
plot_dispersion(fcs_hiphive_enforced, 'tab:green', 'hiphive rotational')

# BONUS!
sc = StructureContainer(cs)
sc.add_structure(supercell)
opt = Optimizer(sc.get_fit_data())
opt.train()
parameters_fitted = opt.parameters
fcp_fitted = ForceConstantPotential(cs, parameters_fitted)
fcs_fitted = fcp_fitted.get_force_constants(supercell)
plot_dispersion(fcs_fitted, 'tab:red', 'hiphive fitted')

# fit enforce
enforced_parameters = enforce_rotational_sum_rules(cs, parameters_fitted, ['Huang', 'Born-Huang'])
fcp_enforced = ForceConstantPotential(cs, enforced_parameters)
fcs_hiphive_enforced = fcp_enforced.get_force_constants(supercell)
plot_dispersion(fcs_hiphive_enforced, 'tab:purple', 'hiphive fitted rotational')

# show
plt.legend()
plt.xlim(0, 0.15)
plt.ylim(-5, 25)

plt.xlabel('Wave vector')
plt.ylabel('Energy (meV)')

plt.savefig('graphene_phonon_dispersion.svg')
plt.show()