Third order force constants with external thermal calculators

[lishouhang123456]

I am trying to convert the third-order force constants with other external thermal calculators, like ShengBTE and Phonopy, as the direct solution for the phonon Boltzmann transport equation is not included in SSCHA. I followed the instructions on GitHub (https://github.com/orgs/SSCHAcode/discussions/104) and successfully obtained the force constants in ShengBTE form. However, there are large differences in the lattice thermal conductivity (RTA level) from SSCHA, ShengBTE, and phono3py.

I tested on the silicon example. The RTA thermal conductivity for silicon at 300 K from SSCHA, ShengBTE, and phono3py are 251.32 W/mK, 1.375, and 10.289 W/mK, respectively. I think there is something wrong with the unit of the third-order force constants. I paste the Python code to calculate the second- and third-order force constants in the following.

import cellconstructor.ForceTensor
import numpy as np
from quippy.potential import Potential
from ase import Atoms
import ase.io
from ase.eos import calculate_eos
from ase.units import kJ
from ase.phonons import Phonons as AsePhonons
from ase.constraints import ExpCellFilter
from ase.optimize import BFGS, QuasiNewton
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.Structure
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax

def get_starting_dynamical_matrices(structure_filename, potential, supercell):
atoms = ase.io.read(structure_filename)
atoms.set_calculator(potential)
ecf = ExpCellFilter(atoms)
qn = QuasiNewton(ecf)
qn.run(fmax=0.0005)
structure = cellconstructor.Structure.Structure()
structure.generate_from_ase_atoms(atoms, get_masses=True)
dyn = cellconstructor.Phonons.compute_phonons_finite_displacements(structure,potential,supercell=supercell)
dyn.Symmetrize()
dyn.ForcePositiveDefinite()
eos = calculate_eos(atoms)
v0, e0, B = eos.fit()
bulk = B / kJ * 1.0e24
return dyn, bulk

temperature = 300.0
nconf = 2000
max_pop = 1000

pot = Potential('IP Tersoff', param_filename='./ip.parms.Tersoff.xml')
supercell = tuple((4*np.ones(3,dtype=int)).tolist())
dyn, bulk = get_starting_dynamical_matrices('./POSCAR', pot, supercell)

ensemble = sscha.Ensemble.Ensemble(dyn, T0 = temperature, supercell=dyn.GetSupercell())
ensemble.generate(N=nconf)
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minimizer.min_step_dyn = 0.1
minimizer.kong_liu_ratio = 0.5
minimizer.meaningful_factor = 0.001

minimizer.max_ka = 4000

relax = sscha.Relax.SSCHA(minimizer,ase_calculator=pot, N_configs=nconf, max_pop = max_pop, save_ensemble=True)
relax.vc_relax(static_bulk_modulus=bulk, ensemble_loc="directory_of_the_ensemble")

new_ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0=temperature, supercell=relax.minim.dyn.GetSupercell())
new_ensemble.generate(N = nconf*5)
new_ensemble.compute_ensemble(pot, compute_stress=True, stress_numerical=False, cluster = None, verbose = True)

new_minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(new_ensemble)
new_minimizer.minim_struct = False
new_minimizer.set_minimization_step(0.1)
new_minimizer.meaningful_factor=0.001
new_minimizer.max_ka = 2000
new_minimizer.init()
new_minimizer.run()
new_minimizer.dyn.save_qe('final_dyn')

new_ensemble.update_weights(new_minimizer.dyn, temperature)

dyn_hessian, d3_tensor = new_ensemble.get_free_energy_hessian(include_v4=False, get_full_hessian=True, verbose = True, return_d3 = True)
np.save("d3.npy",d3_tensor)
dyn_hessian.save_qe('hessian_dyn')

I use this block to convert the third-order force constants

tensor3 = cellconstructor.ForceTensor.Tensor3(new_minimizer.dyn.structure,new_minimizer.dyn.structure.generate_supercell(supercell),supercell)
tensor3.SetupFromTensor(d3_tensor)
tensor3.tensor *= CC.Units.RY_TO_EV / CC.Units.BOHR_TO_ANGSTROM**3
tensor3.WriteOnFile(fname='FORCE_CONSTANTS_3RD',file_format='phonopy')

I believe the unit for the third-order force constants from SCCHA is RY/BOHR³, while ShengBTE uses eV/ų. I am unsure about the cause of the final discrepancy in thermal conductivities. If anyone has encountered a similar issue, could you please share your experience with converting third-order force constants for use in other thermal conductivity calculators?

[DjordjeDangic]

Hi, there are functions in the Methods class that translate SSCHA force constants to PHONOPY/PHONO3PY ones: tensor2_to_phonopy_fc2(SSCHA_tensor, phonon) and tensor3_to_phonopy_fc3(SSCHA_tensor, phonon). You can try them with phonopy (SSCHA.Phonons is dyn and d3 is 3rd order tensor which depending on how you save it might have a factor of two):

unitcell = PhonopyAtoms(symbols=dyn.structure.atoms, cell=dyn.structure.unit_cell, positions=dyn.structure.coords)
phonon = Phonopy(unitcell, dyn.GetSupercell(), primitive_matrix = np.eye(3).tolist())
tensor2 = CC.ForceTensor.Tensor2(dyn.structure, dyn.structure.generate_supercell(dyn.GetSupercell()), dyn.GetSupercell())
tensor2.SetupFromPhonons(dyn)
phonon.force_constants = CC.Methods.tensor2_to_phonopy_fc2(tensor2, phonon)
fc31 = CC.ForceTensor.Tensor3(dyn.structure, dyn.structure.generate_supercell(dyn.GetSupercell()), dyn.GetSupercell())
fc31.SetupFromTensor(d3)
force_constants_3rd_order = CC.Methods.tensor3_to_phonopy_fc3(fc31, phonon)
thermal_conductivity = phono3py.Phono3py(unitcell, dyn.GetSupercell(), primitive_matrix = np.eye(3).tolist())

[lishouhang123456]

Dear Djordje,

It works now! Thanks for your great answer.

Best regards,
Shouhang