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)

import os
from ase.io import write
from ase.build import bulk
from hiphive.structure_generation import generate_rattled_structures
from hiphive.utilities import prepare_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)

structures = generate_rattled_structures(atoms_ideal, number_of_structures, rattle_amplitude)
structures = prepare_structures(structures, atoms_ideal, calc)

# 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

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.

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.

from ase.io import read
from hiphive import StructureContainer, ForceConstantPotential
from trainstation 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.

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