In this article, we explain how to perform phonon calculations using "Phonopy," a package specifically designed for phonon analysis. With Phonopy, you can calculate phonon band structures and density of states (DOS), and use those results to derive thermodynamic properties.
When plotting phonon bands, it is convenient to use the sumo package to automatically determine the paths in reciprocal space based on the given cell geometry. In this guide, we will use sumo alongside Phonopy.
Both packages can be easily installed via the pip install command:
!pip install phonopy
!pip install sumo
In this example, we will use Silicon (Si) with a diamond structure.
To calculate phonons, structural optimization is the first required step. Phonon calculations are based on changes in forces when small atomic displacements are applied to a structure at its energy minimum. Therefore, it is recommended to use stricter convergence criteria for optimization.
Additionally, using a primitive unit cell for the input structure allows for results that correctly reflect the symmetry of the system.
Below, we define the get_primitive_atoms and optimize_atoms functions. We use the former to convert the input structure into a primitive unit cell and the latter to perform structural optimization (with a convergence criterion of fmax=0.0001).
from ase.build import bulk
from ase.optimize import FIRE, LBFGS
from ase.filters import FrechetCellFilter
from ase.constraints import FixSymmetry
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
from pymatgen.io.ase import AseAtomsAdaptor
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
def get_primitive_atoms(atoms_init):
atoms = atoms_init.copy()
adaptor = AseAtomsAdaptor()
struct = adaptor.get_structure(atoms)
sga = SpacegroupAnalyzer(struct, symprec=0.01)
struct_prim = sga.get_primitive_standard_structure()
atoms_prim = adaptor.get_atoms(struct_prim)
return atoms_prim
def optimize_atoms(atoms_init, calculator, opt_algo=LBFGS, fmax=0.0001, steps=10000):
atoms = atoms_init.copy()
atoms.calc = calculator
atoms.set_constraint(FixSymmetry(atoms))
cf = FrechetCellFilter(atoms) # Allows cell-variable optimization
Epot = atoms.get_potential_energy()
print(f"E_pot before opt: {Epot} eV\n")
opt = opt_algo(cf, logfile=None)
opt.run(fmax=fmax, steps=steps)
atoms.set_constraint()
cell_diff = (atoms.cell.cellpar() / atoms_init.cell.cellpar() - 1.0) * 100
print("Optimized Cell:\n", atoms.cell.cellpar(), "\n")
print("Optimized Cell diff (%):\n", cell_diff, "\n")
print("Scaled positions:\n", atoms.get_scaled_positions(), "\n")
Epot = atoms.get_potential_energy()
print(f"E_pot after opt: {Epot} eV")
return atoms
estimator = Estimator(calc_mode=EstimatorCalcMode.R2SCAN, model_version="v8.0.0")
calculator = ASECalculator(estimator)
atoms = bulk("Si")
atoms.calc = calculator
atoms_opt = optimize_atoms(atoms, calculator)
When performing phonon calculations with Phonopy, we use the PhonopyAtoms object. Since conversion between ASE's Atoms object and Phonopy's object is necessary, we define the following conversion functions:
from ase import Atoms
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
def phonopy_to_ase(phonopy_atoms: PhonopyAtoms):
atoms_ase = Atoms(symbols=phonopy_atoms.symbols,
positions=phonopy_atoms.positions,
cell=phonopy_atoms.cell,
pbc=True)
return atoms_ase
def ase_to_phonopy(atoms: Atoms):
atoms_phonopy = PhonopyAtoms(symbols=atoms.get_chemical_symbols(),
cell=atoms.cell,
positions=atoms.get_positions())
return atoms_phonopy
Since small atomic displacements do not necessarily maintain the periodicity of the primitive unit cell, we use a supercell for the calculation. Larger supercells yield more precise results for long-wavelength vibrations, but the computational cost increases significantly. Therefore, a reasonable cutoff is required. Here, we generate supercells with displacement patterns based on a 4x4x4 supercell of the primitive unit cell. The displacement distance is specified using the distance argument in the phonon.generate_displacements() method.
# Specify supercell size
supercell_matrix = [4,4,4]
# Convert ASE Atoms object to Phonopy PhonopyAtoms object
unitcell = ase_to_phonopy(atoms_opt)
# Create Phonopy object
phonon = Phonopy(unitcell, supercell_matrix)
# Generate atomic displacements
phonon.generate_displacements(distance=0.1)
supercells = phonon.supercells_with_displacements
Phonon frequencies are calculated by creating structures where atoms are slightly displaced from their equilibrium positions in various patterns. The force constants are then calculated using the finite difference method based on the calculated forces.
Calculate the forces for the generated displaced structures. The force constants are then calculated using the finite difference method based on the calculated forces.
# Calculate forces
forces = []
for sc in supercells:
sc_ase = phonopy_to_ase(sc)
sc_ase.calc = calculator
forces.append(sc_ase.get_forces())
# Calculate force constants
phonon.forces = forces
phonon.produce_force_constants()
Once the force constant matrix is obtained, we can perform a Fourier transform to find the eigenfrequencies at any arbitrary point in reciprocal space. Plotting these along a specific path in reciprocal space results in a phonon band structure.
Generally, band structures are plotted along paths connecting high-symmetry points that reflect the symmetry of the cell. Standard paths are determined for each cell symmetry. Here, we use sumo to automatically generate the path for plotting the bands.
from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections
from sumo.symmetry.kpoints import get_path_data
from pymatgen.io.ase import AseAtomsAdaptor
# Convert ASE Atoms object to pymatgen Structure object
structure = AseAtomsAdaptor.get_structure(atoms_opt * supercell_matrix)
kpath, band_paths, labels = get_path_data(structure, mode="pymatgen", symprec=0.0001, line_density=60, phonopy=True)
print("Generated path labels:", labels)
print(len(band_paths))
Now, we plot the band structure using the generated paths and high-symmetry point labels.
labels = ['$'+label+'$' for label in labels]
# Plot phonon band structure
phonon.run_band_structure(paths=band_paths, labels=labels)
phonon.plot_band_structure().show()
Unlike the band structure, the Density of States (DOS) is calculated by sampling phonon frequencies on a grid (mesh) within the unit structure in reciprocal space (the first Brillouin zone) and counting the number of vibrational states per frequency. Here, we calculate and plot the DOS using a 20x20x20 mesh.
q_mesh = [20,20,20]
phonon.run_mesh(q_mesh)
phonon.run_total_dos()
phonon.plot_total_dos().show()
It is also possible to plot the phonon band structure and DOS side-by-side. This helps in understanding which bands correspond to specific peaks in the DOS.
phonon.run_mesh(q_mesh)
phonon.run_total_dos()
phonon.plot_band_structure_and_dos().show()
From the phonon DOS, you can calculate the temperature dependence of thermodynamic properties such as Helmholtz free energy, entropy, and heat capacity. Note that these results represent contributions derived solely from lattice vibrations (phonons) and do not include electronic contributions. For semiconductors and insulators, electronic contributions are often negligible, and lattice vibrations account for most of the observed properties.
phonon.run_mesh(q_mesh)
phonon.run_thermal_properties(t_step=10, t_max=1000, t_min=0)
tp_dict = phonon.get_thermal_properties_dict()
temperatures = tp_dict['temperatures']
free_energy = tp_dict['free_energy']
entropy = tp_dict['entropy']
heat_capacity = tp_dict['heat_capacity']
for t, F, S, cv in zip(temperatures, free_energy, entropy, heat_capacity):
print(("%12.3f " + "%15.7f" * 3) % ( t, F, S, cv ))
phonon.plot_thermal_properties().show()
Note: The calculations above are based on the harmonic approximation. Physically, this assumes a scenario where there are no interactions between phonons, meaning effects like phonon-phonon scattering and phonon lifetime are not included. While this is a good approximation at low temperatures (specifically below the Debye temperature), discrepancies with reality increase at high temperatures. In such cases, methods that incorporate third-order or higher anharmonic effects are required.