Compare force constants

This tutorial demonstrates that hiPhive reproduces second and third order force constants obtained using phonopy and/or phono3py. Please consult the respective websites for installation instructions for these codes. The present example furthermore employs the lammps calculator from ASE, which requires a working lammps installation.

Preparation of reference data

First we generate reference force constant matrices using phono3py

tutorial/advanced/compare_fcs_phono3py/compute_phono3py_fcs.py
import subprocess
import glob
import os
import shutil
from ase.io import write, read
from tools import get_single_calculator, write_vaspruns

# Parameters
dim = 3
potential_file = 'Si.tersoff'
calc = get_single_calculator(potential_file, 'Si', 'tersoff',
                             pair_coeff_tag='Si(B)')

# Phono3py calculation
work_dir = 'phono3py_calculation/'
if not os.path.exists(work_dir):
    os.makedirs(work_dir)
os.chdir(work_dir)

# read and write POSCAR
atoms_prim = read('../structures/POSCAR')
write('POSCAR', atoms_prim, format='vasp')

# Generate displaced supercells
subprocess.call('phono3py -d --dim=\"%d %d %d\"' % (dim, dim, dim), shell=True)
if not os.path.exists('poscars'):
    os.makedirs('poscars')
else:
    shutil.rmtree('poscars')
    os.makedirs('poscars')
subprocess.call('mv POSCAR-* poscars', shell=True)

# Compute forces and write vasprun-xxx.xml
poscars = glob.glob('poscars/POSCAR-*')
if not os.path.exists('vaspruns'):
    os.makedirs('vaspruns')
else:
    shutil.rmtree('vaspruns')
    os.makedirs('vaspruns')

write_vaspruns('vaspruns/', poscars, calc, translate=True)
subprocess.call('phono3py --cf3 vaspruns/vasprun.xml*', shell=True)

# Write fc2 and fc3
subprocess.call('phono3py --dim=\"%d %d %d\" --sym_fc3r --sym_fc2 --tsym'
                % (dim, dim, dim), shell=True)
os.chdir('..')

Training of force constant potentials

In order to obtain force constant matrices from hiPhive we then generate input structures with random displacements and compute the computing forces.

tutorial/advanced/compare_fcs_phono3py/prepare_training_structures.py
import os
from ase.io import write
from ase.build import bulk
from hiphive.structure_generation import generate_rattled_structures
from tools import get_single_calculator


# Parameters
a0 = 5.4323
dim = 3
rattle_amplitude = 0.02
number_of_structures = 5
potential_file = 'Si.tersoff'
calc = get_single_calculator(potential_file, 'Si', 'tersoff',
                             pair_coeff_tag='Si(B)')
primitive_fname = 'structures/POSCAR'
structures_fname = 'structures/rattled_structures.extxyz'

# Generate rattled structures
atoms_prim = bulk('Si', 'diamond', a=a0)
atoms_ideal = atoms_prim.repeat(dim)

rattled_structures = generate_rattled_structures(
    atoms_ideal, number_of_structures, rattle_amplitude)

structures = []
for structure in rattled_structures:
    structure.set_calculator(calc)
    forces = structure.get_forces()
    displacements = structure.positions - atoms_ideal.get_positions()
    structure.new_array('forces', forces)
    structure.new_array('displacements', displacements)
    structure.calc = None
    structure.positions = atoms_ideal.get_positions()
    structures.append(structure)

# save structures
if not os.path.isdir(os.path.dirname(structures_fname)):
    os.mkdir(os.path.dirname(structures_fname))
write(primitive_fname, atoms_prim)
write(structures_fname, structures)

Next these input structures are compiled into a structure container.

tutorial/advanced/compare_fcs_phono3py/setup_containers.py
from ase.io import read
from hiphive import StructureContainer
from hiphive import ClusterSpace


# read structures
rattled_structures = read('structures/rattled_structures.extxyz@:')

# Second order model
cs = ClusterSpace(rattled_structures[0], [4.0])
sc = StructureContainer(cs)
for structure in rattled_structures:
    sc.add_structure(structure)
sc.write('structure_container2')


# Third order model
cs = ClusterSpace(rattled_structures[0], [4.0, 4.0])
sc = StructureContainer(cs)
for structure in rattled_structures:
    sc.add_structure(structure)
sc.write('structure_container3')

# Fourth order model
cs = ClusterSpace(rattled_structures[0], [4.0, 4.0, 4.0])
sc = StructureContainer(cs)
for structure in rattled_structures:
    sc.add_structure(structure)
sc.write('structure_container4')

Afterwards force constant potentials are trained and used to set up force constant matrices.

tutorial/advanced/compare_fcs_phono3py/compute_hiphive_fcs.py
from ase.io import read
from hiphive import StructureContainer, ForceConstantPotential
from hiphive.fitting import Optimizer
from hiphive.io.phonopy import write_phonopy_fc2, write_phonopy_fc3

supercell = read('phono3py_calculation/SPOSCAR')
# Read structure containers and cluster spaces
sc2 = StructureContainer.read('structure_container2')
sc3 = StructureContainer.read('structure_container3')
sc4 = StructureContainer.read('structure_container4')
cs2 = sc2.get_cluster_space()
cs3 = sc3.get_cluster_space()
cs4 = sc4.get_cluster_space()

# Fit models
opt = Optimizer(sc2.get_fit_data())
opt.train()
print(opt)
fcp2 = ForceConstantPotential(cs2, opt.parameters)

opt = Optimizer(sc3.get_fit_data())
opt.train()
print(opt)
fcp3 = ForceConstantPotential(cs3, opt.parameters)

opt = Optimizer(sc4.get_fit_data())
opt.train()
print(opt)
fcp4 = ForceConstantPotential(cs4, opt.parameters)

# Write force constants
# this would look smoother if we had our own hdf5 format for the full fcs!
fcs2 = fcp2.get_force_constants(supercell)
write_phonopy_fc2('fc2_model2.hdf5', fcs2)

fcs3 = fcp3.get_force_constants(supercell)
write_phonopy_fc2('fc2_model3.hdf5', fcs3)
write_phonopy_fc3('fc3_model3.hdf5', fcs3)

fcs4 = fcp4.get_force_constants(supercell)
write_phonopy_fc2('fc2_model4.hdf5', fcs4)
write_phonopy_fc3('fc3_model4.hdf5', fcs4)

Comparison of force constant matrices

Finally, the force constant matrices from hiPhive and phono3py are compared using the Froebenius norm.

\[\Delta = ||\matrix{\Phi}_\text{hiPhive} - \matrix{\Phi}_\text{phono3py}||_2 / ||\matrix{\Phi}_\text{phono3py}||_2\]

We also compute the relative error for the (easily computed), \(\Gamma\) frequencies obtained using the different second order force constants.

tutorial/advanced/compare_fcs_phono3py/compare_fcs.py
import numpy as np
from ase.io import read
from hiphive.io.phonopy import read_phonopy_fc2, read_phonopy_fc3
from hiphive.core.utilities import compute_gamma_frequencies


def compute_relative_error(fcs1, fcs2):
    """ Return relative fc difference in percentage """
    return 100.0 * np.linalg.norm(fcs1 - fcs2) / np.linalg.norm(fcs2)


# Read phonopy fcs
atoms = read('phono3py_calculation/SPOSCAR')
fc2_phono3py = read_phonopy_fc2('phono3py_calculation/fc2.hdf5')
fc3_phono3py = read_phonopy_fc3('phono3py_calculation/fc3.hdf5')
omega_phonopy = compute_gamma_frequencies(fc2_phono3py, atoms.get_masses())

# Read hiphive fcs
fc2_model2 = read_phonopy_fc2('fc2_model2.hdf5')
omega_model2 = compute_gamma_frequencies(fc2_model2, atoms.get_masses())

fc2_model3 = read_phonopy_fc2('fc2_model3.hdf5')
fc3_model3 = read_phonopy_fc3('fc3_model3.hdf5')
omega_model3 = compute_gamma_frequencies(fc2_model3, atoms.get_masses())


fc2_model4 = read_phonopy_fc2('fc2_model4.hdf5')
fc3_model4 = read_phonopy_fc3('fc3_model4.hdf5')
omega_model4 = compute_gamma_frequencies(fc2_model4, atoms.get_masses())

# Compare fcs
fc2_model2_error = compute_relative_error(fc2_model2, fc2_phono3py)
omega_model2_error = compute_relative_error(omega_phonopy, omega_model2)
print('Model 2: FC2 error        {:.4f} %'.format(fc2_model2_error))
print('Model 2: Frequency error  {:.4f} %'.format(omega_model2_error))

print('')

fc2_model3_error = compute_relative_error(fc2_model3, fc2_phono3py)
fc3_model3_error = compute_relative_error(fc3_model3, fc3_phono3py)
omega_model3_error = compute_relative_error(omega_phonopy, omega_model3)
print('Model 3: FC2 error        {:.4f} %'.format(fc2_model3_error))
print('Model 3: FC3 error        {:.4f} %'.format(fc3_model3_error))
print('Model 3: Frequency error  {:.4f} %'.format(omega_model3_error))

print('')

fc2_model4_error = compute_relative_error(fc2_model4, fc2_phono3py)
fc3_model4_error = compute_relative_error(fc3_model4, fc3_phono3py)
omega_model4_error = compute_relative_error(omega_phonopy, omega_model4)
print('Model 4: FC2 error        {:.4f} %'.format(fc2_model4_error))
print('Model 4: FC3 error        {:.4f} %'.format(fc3_model4_error))
print('Model 4: Frequency error  {:.4f} %'.format(omega_model4_error))

The results are compiled in the following table, which illustrates that it is crucial to include higher-order terms in the expansion in order to recover the (lower-order) terms of interest.

Terms FC2 error (%) FC3 error (%) Frequency error (%)
2nd 1.6209   1.5564
2nd+3rd 0.8561 2.5278 0.4223
2nd+3rd+4th 0.1129 1.1812 0.0601