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

[lishouhang123456]

Dear Djordje,

I can obtain reasonable thermal conductivity using phono3py with the force constants converted from SSCHA by following your instructions. Recently, I attempted to convert the phono3py HDF5 force constants to the ShengBTE format using hiphive to validate my results. However, I found that the size of the third-order force constants from SSCHA is [num_atom_unit_cell, num_atom_super_cell, num_atom_super_cell, 3, 3, 3], which differs from the official phono3py format with the size of [num_atom_super_cell, num_atom_super_cell, num_atom_super_cell, 3, 3, 3]. As a result, I encountered an error when trying to convert the HDF5 force constants to the ShengBTE format. I am confused about the following two questions:

1. Why does SSCHA provide a different format of HDF5 force constants, which still works well in phono3py?

2. How can I obtain the phono3py format of HDF5 force constants, specifically with the size of [num_atom_super_cell, num_atom_super_cell, num_atom_super_cell, 3, 3, 3]?

I hope my questions are clear. Thank you for your help!

Best regards,
Shouhang

[DjordjeDangic]

Hello Shouhang,

1. Because we have a crystal (ie. discrete translational symmetry) [num_atom_unit_cell, num_atom_super_cell, num_atom_super_cell, 3, 3, 3] format of third order force constants is enough to fully represent third order anharmonicity for supercell with num_atom_super_cell number of atoms. SSCHA has that format probably because whoever wrote it thought it is more simple to represent them like that.
2. You could build the entire [num_atom_super_cell, num_atom_super_cell, num_atom_super_cell, 3, 3, 3] tensor from the smaller one respecting translational symmetry. Also, there might be a phono3py routine that does that for you. I kinda remember there is a function distribute_forceconstants or something similar that does exactly this. You should ask that in Phono3py forum though.

[lishouhang123456]

Dear Djordje,

Got it. I found that Phono3py can automatically recognize whether the force constants are compact or full. That's why the compact SSCHA d3 tensor also works in Phono3py. Thanks for your explanation!

Best regards,
Shouhang