Force constants comparison

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 the structures (primitive cell and training structures)

examples/compare_fcs_phono3py/1_prepare_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)

Secondly we compute the reference force constants using phono3py

examples/compare_fcs_phono3py/3_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'
                % (dim, dim, dim), shell=True)
os.chdir('..')

Training of force constant potentials

In order to obtain force constant matrices from hiPhive we use the generated training structures and compiled into a structure container.

examples/compare_fcs_phono3py/2_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.

examples/compare_fcs_phono3py/4_compute_hiphive_fcs.py
from ase.io import read
from hiphive import StructureContainer, ForceConstantPotential
from hiphive.fitting import Optimizer


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.cluster_space
cs3 = sc3.cluster_space
cs4 = sc4.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)
fcs2.write('model2.fcs')

fcs3 = fcp3.get_force_constants(supercell)
fcs3.write('model3.fcs')

fcs4 = fcp4.get_force_constants(supercell)
fcs4.write('model4.fcs')

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.

examples/compare_fcs_phono3py/5_compare_fcs.py
import numpy as np
from ase.io import read
from hiphive import ForceConstants


def compute_array_relative_error(array1, array2):
    """ Computes relative error between two arrays """
    return 100.0 * np.linalg.norm(array1 - array2) / np.linalg.norm(array2)


def compute_relative_errors(fcs1, fcs2):
    """ Return relative fc difference in percentage for second and third order """
    fc2_array1 = fcs1.get_fc_array(order=2)
    fc2_array2 = fcs2.get_fc_array(order=2)
    fc2_error = compute_array_relative_error(fc2_array1, fc2_array2)
    if 3 in fcs1.orders and 3 in fcs2.orders:
        fc3_array1 = fcs1.get_fc_array(order=3)
        fc3_array2 = fcs2.get_fc_array(order=3)
        fc3_error = compute_array_relative_error(fc3_array1, fc3_array2)
    else:
        fc3_error = None
    return fc2_error, fc3_error


# Read phono3py fcs
atoms = read('phono3py_calculation/SPOSCAR')
fcs2_phono3py = ForceConstants.read_phonopy(atoms, 'phono3py_calculation/fc2.hdf5')
fcs3_phono3py = ForceConstants.read_phono3py(atoms, 'phono3py_calculation/fc3.hdf5')
fcs_phono3py = fcs2_phono3py + fcs3_phono3py
omega_phonopy = fcs_phono3py.compute_gamma_frequencies()

# Read hiphive fcs
fcs_model2 = ForceConstants.read('model2.fcs')
omega_model2 = fcs_model2.compute_gamma_frequencies()

fcs_model3 = ForceConstants.read('model3.fcs')
omega_model3 = fcs_model3.compute_gamma_frequencies()

fcs_model4 = ForceConstants.read('model4.fcs')
omega_model4 = fcs_model4.compute_gamma_frequencies()


# Compare fcs
fc2_model2_error, _ = compute_relative_errors(fcs_model2, fcs_phono3py)
omega_model2_error = compute_array_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, fc3_model3_error = compute_relative_errors(fcs_model3, fcs_phono3py)
omega_model3_error = compute_array_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, fc3_model4_error = compute_relative_errors(fcs_model4, fcs_phono3py)
omega_model4_error = compute_array_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