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

../_images/graphene_phonon_dispersion.svg

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.set_force_constants(fcs.get_fc_array(order=2))

    # get dispersion
    phonopy.set_band_structure(bands)
    _, qnorms, freqs, _ = phonopy.get_band_structure()
    qnorms = np.hstack(qnorms)
    freqs = THz_to_meV * np.vstack(freqs)
    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.set_forces([supercell.arrays['forces']])
phonopy.produce_force_constants()
fcs_phonopy = ForceConstants.from_arrays(supercell, phonopy.get_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')