Wrong 2nd order force constants for silion

hiPhive
1

I am currently working on extracting force constants for Silicon using hiphive and have encountered an issue where the obtained force constants appear to be incorrect.

Here are the details of my workflow:

  • I generated rattled structures using TDEP’s canonical ensemble approach.
  • DFT calculations were performed on these displaced structures to obtain forces.
  • I constructed a StructureContainer in hiphive with the ideal Silicon structure and the corresponding displacements and forces from the DFT runs.
  • I then trained a ForceConstantPotential and extracted the second-order force constants.
  • However, the resulting force constants (and derived phonon dispersion via Phonopy) do not match the expected values for Silicon.
    This is my script:
    import os
    from ase.io import read, write
    from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential
    from hiphive.utilities import prepare_structures
    from trainstation import Optimizer

parent_dir = “rattled_structures” # parent directory containing disp-* folders
unitcell_file = “POSCAR” # primitive/unit cell POSCAR
supercell_file = “SPOSCAR” # ideal supercell POSCAR
prim = read(unitcell_file, format=“vasp”)
write(“prim.extxyz”, prim, format=“extxyz”)

atoms_ideal = read(supercell_file, format=“vasp”)
write(“supercell_ideal.extxyz”, atoms_ideal, format=“extxyz”)

def get_total_energy_from_outcar(outcar_file):
“”“Read the last TOTEN energy from VASP OUTCAR.”“”
energy = None
with open(outcar_file, ‘r’) as f:
for line in f:
if “TOTEN” in line:
parts = line.split()
energy = float(parts[4])
if energy is None:
raise ValueError(f"No TOTEN energy found in {outcar_file}")
return energy

folders = sorted(
[f for f in os.listdir(parent_dir)
if os.path.isdir(os.path.join(parent_dir, f)) and f.startswith(“disp-”)]
)

supercells =
for folder in folders:
folder_path = os.path.join(parent_dir, folder)

# structure
struct_file = os.path.join(folder_path, "CONTCAR")
if not os.path.exists(struct_file):
    struct_file = os.path.join(folder_path, "POSCAR")
atoms = read(struct_file, format="vasp")

# forces (from last ionic step of OUTCAR)
outcar_file = os.path.join(folder_path, "OUTCAR")
if not os.path.exists(outcar_file):
    raise FileNotFoundError(f"{outcar_file} not found")
ref_atoms = read(outcar_file, index=-1)  # last step
atoms.arrays["forces"] = ref_atoms.get_forces()

# total energy
atoms.info["energy"] = get_total_energy_from_outcar(outcar_file)

supercells.append(atoms)

write training set

write(“supercells_rattled.extxyz”, supercells, format=“extxyz”)
print(“:white_check_mark: Wrote prim.extxyz, supercell_ideal.extxyz, supercells_rattled.extxyz”)

cutoffs = [5.0, 4.0, 4.0]
cs = ClusterSpace(prim, cutoffs)
print(cs)
cs.print_orbits()

structure container

structures = prepare_structures(supercells, atoms_ideal)
sc = StructureContainer(cs)
for structure in structures:
sc.add_structure(structure)
print(sc)

optimizer

opt = Optimizer(sc.get_fit_data())
opt.train()
print(opt)

trained potential

fcp = ForceConstantPotential(cs, opt.parameters)
fcp.write(“Si_trained.fcp”)
print(“:white_check_mark: Trained FCP written to Si_trained.fcp”)

supercell = read(“supercell_ideal.extxyz”)
fcp = ForceConstantPotential.read(“Si_trained.fcp”)

fc_matrix = fcp.get_force_constants(supercell)
fc_array = fc_matrix.get_fc_array(order=2)
print(“Force constants array shape:”, fc_array.shape)

fc_matrix.write_to_phonopy(“FORCE_CONSTANTS_Si”, format=“text”)
print(“:white_check_mark: FORCE_CONSTANTS_Si written for Phonopy”)

2

Typically the problematic step is that phonopy uses another permutation of the supercell atoms. Try to read the supercell from phonopy and write the fcs using that one.