First Principles Phonon Calculations with Phonopy: A Comprehensive Guide for Materials Research

Abstract

This article provides a comprehensive guide to performing first-principles phonon calculations using the open-source Phonopy code, a cornerstone tool in computational materials science. It covers foundational concepts of lattice dynamics and the essential roles phonons play in material properties, from thermodynamic behavior to phase stability. A detailed, step-by-step methodology for setting up and running Phonopy calculations is presented, alongside advanced applications that connect phonon properties to functional material behavior. The guide also addresses common troubleshooting scenarios, performance optimization techniques, and provides a critical evaluation of Phonopy's capabilities, including its integration with modern machine-learning interatomic potentials and its position within the broader ecosystem of phonon simulation software.

Understanding Phonons and Why Phonopy is Indispensable in Materials Science

Phonons, the quantized lattice vibrations in crystalline materials, play an essential role in determining a wide array of material properties and behaviors. These collective atomic oscillations govern fundamental aspects such as thermal conductivity, phase stability, and ionic diffusion mechanisms. The importance of first principles phonon calculations in materials science cannot be overly emphasized, as they provide crucial insights into dynamical behaviors and thermal properties that are central topics in fundamental issues of materials science [1]. With the advent of sophisticated computational tools like Phonopy, an open-source code world-widely used for such calculations, researchers can now accurately predict and analyze phonon-related phenomena across diverse material systems [1].

The investigation of phonon behaviors extends from nearly instable systems to superionic conductors exhibiting sublattice melting phenomena. In materials such as Mg₃Sb₂, AgCrSe₂, and Li₆PS₅Cl, phonon properties directly influence transport mechanisms through increasingly anharmonic interactions [2]. Understanding these phonon dynamics through computational methods enables researchers to design novel materials with tailored thermal and ionic conduction properties for specific applications in energy storage and conversion systems.

Theoretical Framework: Phonon Fundamentals

Basic Principles of Lattice Dynamics

Phonons represent the normal modes of vibration in a crystalline lattice, described mathematically through the dynamical matrix derived from interatomic forces. The harmonic approximation serves as the foundational starting point, where atoms oscillate with small displacements around their equilibrium positions. However, real materials invariably exhibit anharmonic effects that become particularly significant at elevated temperatures and in materials with complex bonding environments [2]. These anharmonic interactions are responsible for crucial phenomena including thermal expansion, phonon-phonon scattering, and the temperature dependence of thermal conductivity.

The breakdown of the harmonic approximation becomes evident in superionic conductors where anomalous phonon behaviors emerge. In such materials, the traditional phonon-liquid picture fails to accurately describe the observed dynamics, as demonstrated by the persistence of long wavelength transverse acoustic (TA) phonons in the superionic state of AgCrSe₂ [2]. This deviation from classical models underscores the necessity of employing advanced computational methods that can capture the full complexity of lattice dynamics in functional materials.

Phonon DOS and Thermal Properties

The Phonon Density of States (DOS) provides critical insights into the vibrational characteristics of materials, representing the distribution of phonon modes across different frequencies. Experimental techniques such as inelastic neutron scattering (INS) and inelastic x-ray scattering (IXS) coupled with computational approaches like ab initio molecular dynamics (AIMD) enable comprehensive characterization of phonon DOS [2]. Abnormal features in the phonon DOS, such as low-energy shoulders observed in Mg₃X₂ (X = Sb, Bi), often indicate lattice instabilities arising from ionic size mismatches and weakened atomic bonds [2].

The connection between phonon DOS and macroscopic material properties is particularly evident in thermal transport behavior. Materials exhibiting low thermal conductivity typically display phonon DOS characteristics including soft phonon modes, low-energy anomalies, and enhanced scattering phase space [2]. In superionic conductors like Li₆PS₅Cl, the strong coupling between low-energy Li partial DOS below 10meV and Li-ion diffusion pathways directly links vibrational properties to ionic transport mechanisms [2].

Computational Methods for Phonon Analysis

First-Principles Phonon Calculation Approaches

Computational methods for phonon analysis primarily employ two dominant approaches: Density Functional Perturbation Theory (DFPT) and the finite displacement method. Each technique offers distinct advantages and limitations, making them suitable for different material systems and research objectives. DFPT represents the most efficient method for calculating phonon properties within semilocal DFT frameworks (LDA, PBE), particularly as it enables computation of infrared and Raman intensities for spectroscopic modeling [3]. However, DFPT implementation faces restrictions with ultrasoft pseudopotentials, Hubbard U corrections, and certain dispersion-corrected DFT methods, necessitating alternative approaches in these cases [3].

The finite displacement method, implemented in codes like Phonopy, provides a more generally applicable approach across diverse material systems and Hamiltonian types. This method involves creating supercell structures with atomic displacements to compute force constants, which subsequently yield phonon frequencies and eigenvectors through construction of the dynamical matrix [4]. While computationally more demanding than DFPT for some applications, the finite displacement method offers robust performance across various material classes, including those with complex anharmonicities [3].

Table 1: Comparison of Phonon Calculation Methods

Method	Applicable Systems	Key Advantages	Limitations
DFPT	Semilocal DFT (LDA, PBE) with NCP potentials	Most efficient; Enables IR/Raman intensity calculation	Not implemented for USP, DFT+U, MGGA, hybrid functionals [3]
Finite Displacement	All systems (NCP, USP, DFT+U, MGGA, hybrid)	Universal applicability; Robust for anharmonic systems	Computationally intensive for large systems [3]

Phonopy Protocol for VASP

The Phonopy code provides a comprehensive workflow for calculating phonon properties using density functional theory calculations, particularly with the VASP (Vienna Ab initio Simulation Package) software. The protocol consists of three main stages: pre-processing, force calculations, and post-processing, each with specific requirements and procedures [4].

Pre-process Stage: The initial stage involves generating supercell structures with atomic displacements from the original unit cell. Using the command phonopy -d --dim 2 2 3 --pa auto creates a 2×2×3 supercell with displacements fully considering crystal symmetry [4]. This step produces several critical files: SPOSCAR (perfect supercell structure), phonopy_disp.yaml (displacement information), and multiple POSCAR-{number} files (supercells with specific atomic displacements) [4].

Force Calculation Stage: Each POSCAR-{number} file serves as input for independent VASP calculations to determine the forces on all atoms. The VASP INCAR file must be carefully configured with specific parameters to ensure accurate force calculations without structural relaxation [4]:

After completing all VASP calculations, the FORCE_SETS file is created using the Phonopy VASP interface: phonopy -f disp-{001..003}/vasprun.xml [4].

Post-process Stage: The final stage involves computing various phonon properties from the force constants. Key analyses include:

• Phonon DOS: phonopy --mesh 20 20 20 -p [4]
• Thermal properties: phonopy --mesh 20 20 20 -t [4]
• Projected DOS: phonopy --mesh 20 20 20 --pdos "1 2, 3 4 5 6" -p [4]
• Band structure: phonopy --band "0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0" -p [4]

Non-Analytical Term Correction

For accurate phonon calculations in polar materials, particularly for optical modes near the Brillouin zone center, non-analytical term correction (NAC) must be applied to account for long-range electrostatic interactions. This correction requires the Born effective charge tensor and dielectric constant matrix, which can be computed using VASP with specific parameter settings [4]:

The LEPSILON = .TRUE. setting enables calculation of the dielectric tensor and Born effective charges. After the VASP calculation, the phonopy-vasp-born utility processes the vasprun.xml or OUTCAR files to generate the required BORN file, optionally applying symmetry constraints with the --st parameter for improved numerical stability [4]. Once the BORN file is placed in the working directory, Phonopy automatically activates the non-analytical term correction for subsequent phonon calculations.

Advanced Phonon Characterization Techniques

Experimental- Computational Integration

The integration of experimental measurements with computational approaches provides the most comprehensive framework for understanding phonon behaviors in complex materials. Neutron scattering and x-ray scattering spectroscopy offer direct experimental probes of phonon dynamics, while first-principles calculations provide the theoretical foundation for interpreting these observations [2]. This combined approach is particularly powerful for investigating anharmonic lattice dynamics, where traditional harmonic approximations prove inadequate.

In superionic conductors like AgCrSe₂, combined INS, IXS measurements, and theoretical calculations have revealed complex atomic dynamics that defy conventional models. Contrary to the traditional phonon-liquid picture, long wavelength transverse acoustic (TA) phonons persist in the superionic state, with the ultralow thermal conductivity originating from intrinsic anharmonicity present even at low temperatures [2]. Similarly, in solid-state electrolyte Li₆PS₅Cl, neutron scattering measurements combined with machine-learned molecular dynamics (MLMD) calculations have elucidated the coupling between vibrational dynamics and Li-ion diffusion mechanisms, revealing strong softening of low-energy modes that facilitate ionic conduction [2].

CASTEP Phonon Calculations

The CASTEP software offers robust capabilities for phonon calculations, particularly through its implementation of density-functional perturbation theory (DFPT). A typical CASTEP phonon calculation requires specific keywords in the seedname.cell file, including phonon_kpoint_list to specify phonon wavevectors and appropriate pseudopotential definitions [3]. The following example illustrates a complete input for boron nitride in the wurtzite structure:

CASTEP phonon calculations yield comprehensive output including vibrational frequencies, irreducible representations, infrared intensities, and Raman activities, providing complete spectroscopic characterization for comparison with experimental observations [3].

Table 2: Phonon Calculation Software and Capabilities

Software	Calculation Methods	Key Features	Typical Applications
Phonopy	Finite displacement	Open-source, works with multiple DFT codes	Phonon DOS, band structure, thermal properties [4]
CASTEP	DFPT, Finite displacement	Integrated with CASTEP DFT code	IR/Raman spectra, full phonon dispersion [3]
VASP	DFPT, Finite displacement (via Phonopy)	High accuracy, extensive functionality	Complex materials, anharmonic systems [4]

Case Studies: Phonons in Functional Materials

Thermoelectric Materials: Mg₃Sb₂

In thermoelectric materials like Mg₃Sb₂, phonon engineering plays a crucial role in optimizing the thermoelectric figure of merit by reducing lattice thermal conductivity while maintaining electronic properties. The anomalously low thermal conductivity in Mg₃X₂ (X = Sb, Bi) compared to isostructural ternaries CaMg₂X₂ and YbMg₂X₂ with heavier elements originates from abnormally soft phonons and extra low-energy shoulders in the phonon density of states [2]. These phonon anomalies reflect a near-instability arising from ionic size mismatch that weakens the Mg-X bonds, significantly enhancing phonon scattering phase space and disrupting phonon propagation [2].

This case study demonstrates how first-principles phonon calculations can identify the fundamental origins of unusual thermal transport properties, guiding the design of improved thermoelectric materials through targeted phonon spectrum engineering. The combination of time-of-flight inelastic neutron scattering (INS), inelastic x-ray scattering (IXS), and ab initio molecular dynamics (AIMD) provides a powerful toolkit for correlating microscopic phonon behaviors with macroscopic thermal properties [2].

Superionic Conductors: AgCrSe₂

Superionic conductors represent an extreme case of anharmonic lattice dynamics where phonon behaviors transition toward liquid-like characteristics. In AgCrSe₂, combined experimental and computational investigations have revealed complex atomic dynamics that challenge traditional models [2]. Transverse acoustic (TA) phonons breakdown selectively depending on the involvement of mobile ions, with spectral weight transferring from overdamped vibrational modes to diffusive dynamics in the superionic regime [2].

Contrary to the traditional phonon-liquid picture that assumes complete disappearance of acoustic phonons, IXS measurements show the persistence of long wavelength TA phonons in the superionic state of AgCrSe₂ [2]. This observation, coupled with theoretical calculations, demonstrates that the ultralow thermal conductivity originates from intrinsic anharmonicity present even at low temperature, rather than complete lattice melting. The understanding of these phonon dynamics provides crucial insights for designing improved superionic materials for energy applications.

Solid-State Electrolytes: Li₆PS₅Cl

Solid-state electrolytes for battery applications represent another material class where phonon behaviors directly influence technological performance. In Li₆PS₅Cl, neutron scattering measurements combined with machine-learned molecular dynamics (MLMD) calculations have elucidated the coupling between vibrational dynamics and Li-ion diffusion [2]. INS measurements reveal strong softening of low-energy modes, with the finite value in DOS at energy zero directly linked to Li-ion diffusion processes [2].

Constrained MLMD simulations further demonstrate stronger effects from translation of PS₄ units than rotation on Li diffusion, highlighting the complex interplay between different vibrational modes and ionic transport [2]. The low-energy Li partial DOS below 10meV shows strong coupling with diffusion, providing a quantitative metric for assessing and optimizing ionic conductivity in solid electrolyte materials.

Research Reagent Solutions

Table 3: Essential Computational Tools for Phonon Research

Research Reagent	Function/Purpose	Application Context
Phonopy Code	Open-source package for phonon calculations	Finite displacement method calculations; Post-processing force constants [4]
VASP	First-principles DFT package	Force calculations for displaced structures; Born effective charges [4]
CASTEP	DFT simulation package with phonon capabilities	DFPT phonon calculations; IR/Raman spectrum modeling [3]
INELASTIC Neutron Scattering	Experimental phonon density of states measurement	Validating computational results; Probing low-energy modes [2]
IXS (Inelastic X-ray Scattering)	High-resolution phonon dispersion measurement	Investigating phonon dynamics in small samples [2]
Machine-Learned Molecular Dynamics	Accelerated molecular dynamics with DFT accuracy	Modeling diffusive dynamics in superionic conductors [2]

Phonons play a central role in determining fundamental material properties across diverse applications, from thermoelectrics to solid-state batteries. The continuing development of computational methods, particularly first-principles approaches implemented in codes like Phonopy and CASTEP, has dramatically enhanced our understanding of phonon behaviors in complex materials [4] [3]. The integration of these computational tools with advanced experimental techniques such as inelastic neutron and x-ray scattering provides a powerful framework for investigating anharmonic lattice dynamics and their relationship to macroscopic material properties [2].

Future research directions will likely focus on extending phonon calculations to increasingly complex and disordered materials, developing more efficient methods for treating strong anharmonicity, and integrating machine learning approaches to accelerate phonon calculations and analysis. As these methodologies continue to evolve, they will enable the rational design of materials with optimized thermal and ionic transport properties for next-generation energy technologies.

Theoretical Foundations of Phonon Calculations

Fundamental Theory of Atomic Vibrations

Atomic vibrations in molecules and solids represent fundamental processes that underlie many material properties. According to quantum mechanics, atoms never cease moving, even at absolute zero temperature, and their vibrational behavior changes with internal structure, surrounding environment, and external stimuli. In crystalline solids, these discrete normal modes spread throughout the solid, with coupled motions forming continuous bands in the energy domain and dispersion curves in reciprocal space. The quanta of these vibrational excitations are known as phonons [5].

The potential energy of an atomistic system can be described through a Taylor expansion:

[ V = V0 + \sum{i}\frac{\partial V}{\partial \mathbf{r}i} \cdot \delta\mathbf{r}i + \frac{1}{2}\sum{i,j}\frac{\partial^2 V}{\partial \mathbf{r}i\partial \mathbf{r}j} : \delta\mathbf{r}i\delta\mathbf{r}_j + \cdots ]

Where the negatives of the first derivatives represent interatomic forces, and the second derivatives represent force constants. Within the harmonic approximation, higher-order derivatives are neglected, assuming a parabolic potential energy profile near equilibrium [5].

A normal mode with angular frequency ω and wavevector q corresponds to the displacement of atom α in unit cell k:

[ \mathbf{u}\alpha(k) = \frac{1}{\sqrt{m\alpha}}\mathbf{w}_\alpha(\mathbf{q})e^{i(\mathbf{q}\cdot\mathbf{r}(k) - \omega t)} ]

Here, (\mathbf{w}_\alpha(\mathbf{q})) is the mode's polarization vector of atom α. Applying Newton's second law to each atom yields the equation describing the vibrational dynamics of the solid:

[ \omega^2 w\alpha(\mathbf{q}) = \sum{\alpha'} D{\alpha\alpha'}(\mathbf{q}) w{\alpha'}(\mathbf{q}) ]

Where (D_{\alpha\alpha'}(\mathbf{q})) is a block of the dynamical matrix. Solving for eigenvalues and eigenvectors provides phonon frequencies and polarization vectors at each wavevector [5].

Macroscopic Thermal Properties from Phonons

On a macroscopic level, phonons directly determine several key thermodynamic properties of materials. The vibrational free energy (F), entropy (S), and constant-volume specific heat (C_v) can be calculated using the following expressions within the harmonic approximation [5]:

[ F(T) = kB T \int0^\infty \ln\left[2\sinh\left(\frac{\hbar\omega}{2kB T}\right)\right] g(\omega)d\omega ] [ S(T) = kB \int0^\infty \left[\frac{\hbar\omega}{2kB T}\coth\left(\frac{\hbar\omega}{2kB T}\right) - \ln\left(2\sinh\left(\frac{\hbar\omega}{2kB T}\right)\right)\right] g(\omega)d\omega ] [ Cv(T) = kB \int0^\infty \left(\frac{\hbar\omega}{2kB T}\right)^2 \mathrm{csch}^2\left(\frac{\hbar\omega}{2k_B T}\right) g(\omega)d\omega ]

Where (g(ω)) represents the phonon density of states, (k_B) is Boltzmann's constant, ħ is the reduced Planck constant, and T is temperature [5].

Table 1: Key Phonon-Related Thermal Properties and Their Physical Significance

Property	Mathematical Expression	Physical Significance
Vibrational Free Energy	(F(T) = kB T \int \ln\left[2\sinh\left(\frac{\hbar\omega}{2kB T}\right)\right] g(\omega)d\omega)	Contribution of atomic vibrations to the system's total free energy
Vibrational Entropy	(S(T) = kB \int \left[\frac{\hbar\omega}{2kB T}\coth\left(\frac{\hbar\omega}{2kB T}\right) - \ln\left(2\sinh\left(\frac{\hbar\omega}{2kB T}\right)\right)\right] g(\omega)d\omega)	Measure of disorder arising from atomic vibrations
Heat Capacity	(Cv(T) = kB \int \left(\frac{\hbar\omega}{2kB T}\right)^2 \mathrm{csch}^2\left(\frac{\hbar\omega}{2kB T}\right) g(\omega)d\omega)	Amount of heat required to raise the temperature of the material

In anharmonic systems where higher-order derivatives are non-zero, phonon-phonon scattering occurs, leading to finite phonon lifetime (τ_ω) and thermal conductivity. The intrinsic phonon-phonon scattering rate due to anharmonic three-phonon processes can be expressed as:

[ \tau\omega^{-1} = \frac{\hbar\pi}{16} \sum{\omega',\omega''} |V{\omega,\omega',\omega''}|^2 \left[ (n{\omega'} + n{\omega''} + 1)\delta(\omega - \omega' - \omega'') + 2(n{\omega'} - n_{\omega''})\delta(\omega + \omega' - \omega'') \right] ]

Where (W^{\pm}_{\omega,\omega',\omega''}) represents the three-phonon scattering rates [5]. This anharmonic phonon-phonon scattering causes heat dissipation and finite thermal conductivity, which can be calculated using:

[ \kappa\alpha = \frac{1}{NV0} \sum\omega C{v,\omega} v{\alpha,\omega}^2 \tau\omega ]

Where (κα) denotes the lattice thermal conductivity in the αth direction, (v{α,ω}) is the phonon group velocity, (τω) is the phonon lifetime, and (C{v,ω}) refers to the phonon volumetric specific heat [5].

Phonopy Code Structure and Architecture

Phonopy is an open-source package for phonon calculations at harmonic and quasi-harmonic levels, while Phono3py is a complementary package for phonon-phonon interaction and lattice thermal conductivity calculations [6]. The code has become an essential tool in computational materials science, with its importance "cannot be overly emphasized" according to the developers [1].

Phonopy supports various external force calculators, primarily ab initio codes, with VASP being the default. The code structure employs three primary unit cells [7]:

• Unitcell: The initial input crystal structure
• Supercell: Created from unitcell by the supercell matrix
• Primitive: Created from supercell using primitive matrix

The transformation between these cells follows specific mathematical relationships. Having the supercell matrix (Ms) and primitive matrix (Mp), the basis vectors of supercell are built from those of unitcell by:

[ ( \mathbf{a}\mathrm{s} \; \mathbf{b}\mathrm{s} \; \mathbf{c}\mathrm{s} ) = ( \mathbf{a}\mathrm{u} \; \mathbf{b}\mathrm{u} \; \mathbf{c}\mathrm{u} ) \mathrm{M}_\mathrm{s} ]

The basis vectors of primitive are built from those of supercell by:

[ ( \mathbf{a}\mathrm{p} \; \mathbf{b}\mathrm{p} \; \mathbf{c}\mathrm{p} ) = ( \mathbf{a}\mathrm{s} \; \mathbf{b}\mathrm{s} \; \mathbf{c}\mathrm{s} ) \mathrm{M}\mathrm{s}^{-1} \mathrm{M}\mathrm{p} ]

Once supercell and primitive cells are established, the original unitcell is generally no longer used in subsequent calculations [7].

Computational Workflow

The following diagram illustrates the complete phonon calculation workflow using Phonopy:

Computational Protocols and Methodologies

Installation and Environment Setup

Phonopy can be installed via conda package manager, which provides the most straightforward installation method [8]:

Alternative compilation from source requires several prerequisites including lxml, matplotlib, and PyYAML, as detailed in the search results [9].

Pre-Processing: Generating Displacements

The initial step involves creating a Phonopy object with the unit cell and supercell matrix, followed by displacement generation [7]:

For command-line interface, the displacement generation can be performed using [8]:

This command writes a series of POSCAR files, including SPOSCAR (the supercell with no distortions) and multiple POSCAR files with atomic displacements.

Force Calculator Configuration

The search results provide specific VASP INCAR parameters recommended for phonon calculations [8] [9]:

Table 2: Recommended VASP INCAR Parameters for Phonon Calculations

Parameter	Recommended Value	Purpose
PREC	Accurate	High precision calculation
ENCUT	480 (system dependent)	Plane-wave cutoff energy
EDIFF	1e-8	Strict electronic convergence
IBRION	-1	No ionic relaxation (forces only)
NSW	0	No ionic steps
ISMEAR	0	Gaussian smearing
SIGMA	0.01	Smearing width
LREAL	.FALSE.	Exact projection operators
LCHARG	.FALSE.	No CHGCAR output
LWAVE	.FALSE.	No WAVECAR output
NELM	200	Maximum electronic steps

The ENCUT parameter should be carefully determined based on the POTCAR files, typically 130 eV larger than the largest ENMAX value found using grep 'ENMAX' POTCAR [8].

For other calculators such as Quantum ESPRESSO, the calculator keyword must be specified and set_factor_by_calculator should be set to True for proper unit conversion [7].

Post-Processing and Force Constants

After completing force calculations, the force sets are collected and force constants are produced [7]:

From the command line, force sets can be generated using [10]:

This command processes VASP output files from directories 001 through 010 and creates the FORCE_SETS file needed for subsequent phonon calculations.

Phonon Property Calculations and Analysis

Phonon Band Structure and DOS

Phonon band structures and density of states (DOS) represent fundamental phonon properties. The following diagram illustrates the relationships between different phonon properties and their computational pathways:

For band structure calculation, a configuration file (band.conf) specifies the k-point path [9]:

The band structure is then calculated and plotted using:

For DOS calculations, a mesh sampling is specified [10]:

Thermal Property Calculations

Thermal properties including free energy, entropy, and heat capacity are calculated using the -t flag [8] [10]:

Where mesh.conf contains:

The output provides temperature-dependent thermal properties in tabular form [10]:

Projected DOS and Band Structure with Specified Paths

For more detailed analysis, projected DOS (PDoS) and band structures with specific high-symmetry paths can be calculated. The following configuration file (band-pdos.conf) demonstrates these capabilities [8]:

The corresponding phonopy command generates both band structure and projected DOS:

Application Examples and Case Studies

Silicon Phonon Calculations

Silicon serves as a canonical test case for phonon calculations. The phonon spectrum of silicon exhibits characteristic features including acoustic branches, optical branches, and the absence of imaginary frequencies indicating dynamic stability [10]. The DOS calculation for silicon using a 31×31×31 mesh samples 816 irreducible q-points from 29791 total points on the sampling mesh [10].

NaCl with Non-Analytical Term Correction

For polar materials like NaCl, the non-analytical term correction (NAC) is essential for accurate LO-TO splitting at the Γ-point. This requires preparing a BORN file containing Born effective charges and the dielectric tensor [10]. The band structure calculation for NaCl demonstrates this capability:

Fluorite-Type Metallic Fluorides

First-principles investigations of fluorite-type metallic fluorides (MF₂, where M = Ca, Sr, Ba, and Pb) demonstrate phonopy's application in predicting dynamical stability and thermal properties. Phonon dispersion analyses confirm that all these compounds are dynamically stable, while subsequent calculations provide insights into their thermal conductivity and thermodynamic behavior [11].

The Scientist's Toolkit: Essential Research Reagents

Table 3: Essential Computational Tools and Files for Phonon Calculations

Tool/File	Function	Format/Example
Phonopy Code	Main program for phonon calculations	Python package
Force Calculator	External code for force computations	VASP, Quantum ESPRESSO, ABINIT
POSCAR	Input crystal structure	VASP format
SPOSCAR	Supercell structure with no displacements	VASP format
FORCE_SETS	Collected forces from displacement calculations	YAML format
FORCE_CONSTANTS	Force constants matrix	HDF5 or dat format
BORN File	Born effective charges and dielectric tensor	Required for NAC in polar materials
band.conf	Band structure configuration	Text file with k-path
mesh.conf	Mesh sampling configuration	Text file with mesh parameters

Advanced Topics and Future Directions

AI-Enhanced Phonon Calculations

Recent advancements in data-driven artificial intelligence methodologies have substantially enhanced the efficiency of phonon studies. AI-driven approaches provide orders-of-magnitude improvement in computation efficiency with comparable accuracy, enabling real-time experimental data simulation and interpretation [5]. Key developments include advancements in databases, structural representations, machine-learning interatomic potentials, and graph neural networks.

Compared to traditional techniques, AI methods exhibit transformative potential, dramatically improving the efficiency and scope of research in materials science. These approaches are particularly valuable for high-throughput screening and inverse design of materials with superior thermal properties [5].

Quasi-Harmonic Approximation and Thermal Conductivity

Beyond the harmonic approximation, phonopy supports quasi-harmonic approximation (QHA) for calculating thermal expansion and temperature-dependent phonon properties. Additionally, Phono3py extends capabilities to compute phonon-phonon interactions and lattice thermal conductivity, which are crucial for thermoelectric materials and thermal management applications [6].

The calculation of lattice thermal conductivity involves solving the Boltzmann transport equation with three-phonon scattering processes, requiring second- and third-order force constants. These advanced calculations provide insights for developing materials with tailored thermal transport properties.

The study of atomic vibrations in solids, or lattice dynamics, is a cornerstone of modern condensed matter physics, essential for understanding phenomena such as superconductivity, optical processes, phase transitions, and thermal transport properties [12]. The formal microscopic treatment of lattice dynamics in crystals originates from the work of Born and Von Kármán, who modeled crystals as atoms interacting via spring-like forces within the harmonic approximation [12]. This approach, coupled with the Born-Oppenheimer approximation that decouples electronic and ionic motions, provides the foundational framework for the modern theory of lattice dynamics. Within this harmonic theory, the quantization of atomic eigenvibrations yields quasiparticles known as phonons, which possess well-defined energy and an infinite lifetime, as they do not interact with one another [12]. The harmonic approximation has proven remarkably successful in explaining many physical properties of solids, including thermodynamic functions and phonon-mediated carrier dynamics [12].

The central quantity in lattice dynamics is the Interatomic Force Constant (IFC), which is formally defined as the derivative of the Born-Oppenheimer potential energy surface with respect to atomic displacements [12]. In practice, the potential energy of interacting atoms is approximated by a Taylor expansion with respect to atomic displacements from their equilibrium positions [13]:

[ U - U{0} = \sum{n=1}^{N} U{n} = U{1} + U{2} + U{3} + \cdots, ]

where

[ U{n} = \frac{1}{n!} \sum{\substack{\ell{1}\kappa{1}, \dots, \ell{n}\kappa{n} \ \mu{1},\dots, \mu{n}}} \Phi{\mu{1}\dots\mu{n}}(\ell{1}\kappa{1};\dots;\ell{n}\kappa{n}) \; u{\mu{1}}(\ell{1}\kappa{1})\cdots u{\mu{n}}(\ell{n}\kappa_{n}). ]

Here, (u{\mu}(\ell\kappa)) represents the atomic displacement of the (\kappa)th atom in the (\ell)th unit cell along the (\mu)th direction, and (\Phi{\mu{1}\dots\mu{n}}(\ell{1}\kappa{1};\dots;\ell{n}\kappa{n})) is the (n)th-order interatomic force constant [13]. The second-order IFCs define the harmonic phonon dispersion relations, while the third- and higher-order IFCs govern anharmonic phonon-phonon interactions [12].

Table: Hierarchy of Interatomic Force Constants and Their Roles

IFC Order	Mathematical Representation	Physical Role	Computational Method
Second-order	(\Phi{\mu{1}\mu{2}}(\ell{1}\kappa{1};\ell{2}\kappa_{2}))	Defines harmonic phonon dispersion relations [12]	Finite-displacement [12], DFPT [14]
Third-order	(\Phi{\mu{1}\mu{2}\mu{3}}(\ell{1}\kappa{1};\ell{2}\kappa{2};\ell{3}\kappa{3}))	Governs three-phonon scattering processes [15]	Finite-displacement (Thirdorder.py, Phono3py) [12], CSLD [12]
Fourth-order	(\Phi{\mu{1}\mu{2}\mu{3}\mu{4}}(\ell{1}\kappa{1};\ell{2}\kappa{2};\ell{3}\kappa{3};\ell{4}\kappa_{4}))	Contributes to higher-order phonon interactions and thermal properties [16]	Fourthorder.py [12], CSLD [12], hiPhive [16]

Fundamental Theory of Interatomic Force Constants

Symmetry Properties and Constraints

IFCs possess several fundamental relationships that reduce the number of independent constants and ensure physical realism. These symmetry relationships are crucial for efficient computation and accurate physical modeling [13].

• Permutation Invariance: IFCs are invariant under the exchange of the triplet ((\ell,\kappa,\mu)). For example, a third-order IFC satisfies: [\Phi{\mu{1}\mu{2}\mu{3}}(\ell{1}\kappa{1};\ell{2}\kappa{2};\ell{3}\kappa{3})=\Phi{\mu{1}\mu{3}\mu{2}}(\ell{1}\kappa{1};\ell{3}\kappa{3};\ell{2}\kappa{2})=\dots.] This permutation symmetry significantly reduces the number of unique IFCs that need to be calculated [13].

• Periodicity: Since IFCs depend on interatomic distances, they are invariant under translation by lattice vectors: [\Phi{\mu{1}\mu{2}\dots\mu{n}}(\ell{1}\kappa{1};\ell{2}\kappa{2};\dots;\ell{n}\kappa{n})=\Phi{\mu{1}\mu{2}\dots\mu{n}}(0\kappa{1};\ell{2}-\ell{1}\kappa{2};\dots;\ell{n}-\ell{1}\kappa_{n}).] This property allows the problem to be reduced to interactions within a reference unit cell [13].

• Crystal Symmetry: IFCs must transform under crystal symmetry operations (rotation and translation) according to: [\sum{\nu{1},\dots,\nu{n}}\Phi{\nu{1}\dots\nu{n}}(L{1}K{1};\dots;L{n}K{n}) O{\nu{1}\mu{1}}\cdots O{\nu{n}\mu{n}} = \Phi{\mu{1}\dots\mu{n}}(\ell{1}\kappa{1};\dots;\ell{n}\kappa{n}),] where (O) is the rotational matrix of the symmetry operation. For non-hexagonal systems, this often simplifies to (\Phi{\nu{1}^{\prime}\dots\nu{n}^{\prime}}(L{1}K{1};\dots;L{n}K{n}) = s \Phi{\mu{1}\dots\mu{n}}(\ell{1}\kappa{1};\dots;\ell{n}\kappa_{n})) with (s=\pm 1) [13].

• Translational Invariance: The potential energy must be invariant under rigid translations, imposing the constraint: [\sum{\ell{1}\kappa{1}}\Phi{\mu{1}\mu{2}\dots\mu{n}}(\ell{1}\kappa{1};\ell{2}\kappa{2};\dots;\ell{n}\kappa{n}) = 0,] which must hold for any set of (\ell{2}\kappa{2},\dots,\ell{n}\kappa{n}) and (\mu{1},\dots,\mu_{n}) [13].

• Rotational Invariance: Additionally, the potential energy must be invariant under rigid rotations, leading to more complex constraints that relate (n)th-order IFCs to ((n+1))th-order IFCs. For example, the rotational invariance constraint for harmonic IFCs is given by: [ \sum{\ell{2}\kappa{2}} (\Phi{\mu{1}\nu}(\ell{1}\kappa{1};\ell{2}\kappa{2})r{\mu}(\ell{2}\kappa{2})-\Phi{\mu{1}\mu}(\ell{1}\kappa{1};\ell{2}\kappa{2})r{\nu}(\ell{2}\kappa{2})) + \Phi{\nu}(\ell{1}\kappa{1})\delta{\mu,\mu{1}} - \Phi{\mu}(\ell{1}\kappa{1})\delta{\nu,\mu{1}} = 0, ] where (r{\mu}(\ell\kappa)) represents atomic positions [13].

Calculation Methods for IFCs

Two principal first-principles methods have been developed for computing IFCs:

• Finite-Displacement (Frozen-Phonon) Method: This real-space approach involves calculating the forces acting on all atoms in a supercell after displacing one or more atoms from their equilibrium positions. The second-order IFCs are obtained from the force responses to single atomic displacements, while higher-order IFCs require multiple simultaneous displacements [12]. This method is implemented in codes such as Phonopy [12] [1], Phon [12], thirdorder.py [12], and fourthorder.py [12]. While formally straightforward, its application beyond third-order IFCs becomes computationally prohibitive due to the combinatorial explosion in the number of required displaced configurations [12] [16].

• Density Functional Perturbation Theory (DFPT): This reciprocal-space method employs linear response theory to compute the derivatives of the total energy with respect to atomic displacements directly. A key advantage of DFPT is its use of a variational principle to determine the changes in wavefunctions due to atomic displacements and to compute the corresponding energy change accurately [14]. DFPT is particularly efficient for calculating second-order IFCs and is implemented in packages like Quantum ESPRESSO and Abinit [12]. Its extension to third-order IFCs is available in the D3Q code [12], but implementations beyond third order are hindered by computational complexity and pseudopotential dependencies [12].

Practical Implementation and Protocols

High-Throughput Workflow for Lattice Dynamics

Recent advancements have enabled automated high-throughput frameworks for calculating lattice dynamical properties, including anharmonic effects. These workflows systematically bridge atomic-scale quantum simulations with macroscopic thermal properties [16]. A representative pipeline, as implemented in the atomate package, comprises four key steps [16]:

Table: Key Steps in a High-Throughput Lattice Dynamics Workflow

Step	Primary Task	Tools/Packages	Critical Parameters
1. Structure Preparation & DFT	Stringent optimization of primitive cell; SCF force calculations in displaced supercells [16]	VASP [16]	PBEsol functional (for accurate lattice parameters) [16]
2. IFC Fitting	Extraction of harmonic and anharmonic force constants from force-displacement data [16]	hiPhive [16]	Supercell size (~20 Å) [16], Cutoff radius [16]
3. Phonon Renormalization	Obtaining stable effective phonons at finite temperatures for dynamically unstable compounds [16]	Phonopy [16], Phono3py [16]	Temperature grid, Convergence threshold
4. Property Calculation	Computation of thermal conductivity, expansion coefficient, and free energy [16]	ShengBTE [16], Phono3py [16]	q-point mesh, Scattering processes

• Structure Optimization and Force Calculations: The initial primitive cell undergoes stringent structural optimization. Subsequently, self-consistent field (SCF) force calculations are performed in systematically displaced supercells. The PBEsol exchange-correlation functional is often preferred over PBE due to its more accurate prediction of lattice parameters and phonon frequencies [16].

• Force Constants Fitting: Harmonic and anharmonic IFCs are extracted from the force-displacement datasets using advanced fitting techniques. The hiPhive package, which leverages machine-learning algorithms like compressive sensing or ordinary least-squares regression, is particularly effective for this purpose. The fitting process minimizes ( \parallel {\bf{F}}-{\mathbb{A}}{\mathbf{\Phi }}{\parallel }_{2} ), where ( {\mathbb{A}} ) is the sensing matrix of atomic displacements and ( \mathbf{F} ) is the vector of forces [16]. Key parameters requiring convergence include supercell size (e.g., ~20 Å) and the cutoff radius for IFC interactions [16].

• Renormalization: For dynamically unstable compounds (which exhibit imaginary phonon frequencies within the harmonic approximation), a renormalization step is performed to obtain real effective phonon spectra at finite temperatures. This step is crucial for accurately modeling phase transitions and temperature-dependent stability [16].

• Thermal Property Calculation: The final step involves computing macroscopic thermal properties. The lattice thermal conductivity is typically obtained by solving the Boltzmann transport equation using packages like ShengBTE or Phono3py. The coefficient of thermal expansion and vibrational free energy are also calculated at this stage, incorporating anharmonic corrections [16].

Advanced Fitting Techniques for High-Order IFCs

Conventional finite-displacement methods become prohibitively expensive for high-order IFCs due to combinatorial explosion. Linear-regression-based supercell approaches have emerged as powerful alternatives [12]:

• Compressive Sensing Lattice Dynamics (CSLD): This technique leverages L1 regularization to construct a sparse representation of lattice-dynamical models. CSLD capitalizes on the physical reality that IFCs are generally sparse and decay rapidly with increasing interatomic distance. It effectively avoids overfitting and is implemented in codes like CSLD, Alamode, and hiPhive [12].

• Temperature-Dependent Effective Potential (TDEP): This method extracts effective IFCs from molecular dynamics trajectories, capturing temperature-dependent renormalization effects naturally [16].

• hiPhive Implementation: The hiPhive package provides a flexible framework for extracting IFCs using various fitting methods, including ordinary least-squares and compressive sensing. Its Python integration makes it suitable for automated high-throughput workflows [16].

The computational phonon community has developed a rich ecosystem of open-source software packages that enable first-principles lattice dynamics calculations. The following table details key tools and their primary functions.

Table: Essential Software Tools for Lattice Dynamics Calculations

Tool Name	Primary Function	Key Features	Reference
Phonopy	Harmonic phonon calculations	Finite-displacement method, Thermal properties, QHA, Interface with major DFT codes [1] [17]	[1]
Phono3py	Third-order anharmonic IFCs and thermal conductivity	Three-phonon interactions, LTC from BTE [12]	[12]
ShengBTE	Lattice thermal conductivity	Solves BTE using 2nd and 3rd-order IFCs [16]	[16]
ALAMODE	Anharmonic phonon calculations	Compressive sensing for high-order IFCs, LTC, spectral energy density [13]	[13]
hiPhive	High-order IFC extraction	Machine-learning force constant potentials, Multiple fitting methods [16]	[16]
Pheasy	Comprehensive phonon ecosystem	High-order IFC reconstruction, SCHA/SCP calculations, Interfaces with transport codes [12]	[12]
Quantum ESPRESSO	Ab-initio DFT calculations	DFPT for harmonic IFCs, Plane-wave basis set [12]	[12]
VASP	Ab-initio DFT calculations	Finite-displacement forces, PAW pseudopotentials [16]	[16]
κALDo	Thermal transport calculations	BTE solvers, QHGK method, Crystalline and amorphous solids [15]	[15]

Thermal Properties within the Harmonic Approximation

Theoretical Framework

Within the harmonic approximation, various thermal properties can be derived from the phonon spectrum. The fundamental quantity is the phonon density of states, from which thermodynamic properties are calculated [18]:

• Vibrational Free Energy: [ F(T) = \int{0}^{\infty} \frac{d\omega}{2\pi} g(\omega) \left[ \frac{\hbar\omega}{2} + kB T \ln\left(1 - e^{-\hbar\omega/k_B T}\right) \right] ] where (g(\omega)) is the phonon density of states.

• Entropy: [ S(T) = -kB \int{0}^{\infty} \frac{d\omega}{2\pi} g(\omega) \left[ \ln\left(1 - e^{-\hbar\omega/kB T}\right) - \frac{\hbar\omega}{kB T} \frac{1}{e^{\hbar\omega/k_B T} - 1} \right] ]

• Constant-Volume Heat Capacity: [ CV(T) = kB \int{0}^{\infty} \frac{d\omega}{2\pi} g(\omega) \left(\frac{\hbar\omega}{kB T}\right)^2 \frac{e^{\hbar\omega/kB T}}{\left(e^{\hbar\omega/kB T} - 1\right)^2} ]

• Mean-Square Displacements (MSD): The MSD quantifies thermal vibration amplitudes and is calculated as [18]: [ \text{MSD}(T) = \frac{1}{N} \sum{\mu} \frac{\hbar}{2\omega{\mu}} \coth\left(\frac{\hbar\omega{\mu}}{2kB T}\right) ] where the sum is over all phonon modes (\mu).

Quasi-Harmonic Approximation (QHA)

The harmonic approximation strictly predicts no thermal expansion, as phonon frequencies are volume-independent. The quasi-harmonic approximation (QHA) extends this by introducing volume dependence to phonon frequencies, thereby enabling the calculation of thermal expansion and other temperature-dependent properties at constant pressure [17].

The key implementation in phonopy-qha involves [17]:

• Volume-Energy Data: Electronic total energies at 0K are computed for a series of volumes, typically at least 5 volume points.

• Phonon Calculations at Multiple Volumes: Harmonic phonon calculations are performed at each volume to obtain the vibrational free energy (F_\text{phonon}(T;V)).

• Gibbs Free Energy Minimization: At each temperature, the Gibbs free energy is constructed and minimized with respect to volume: [ G(T, p) = \minV \left[ U(V) + F\mathrm{phonon}(T;\,V) + pV \right] ] where (U(V)) is the static electronic energy, (F_\mathrm{phonon}(T;V)) is the phonon free energy, and (pV) is the pressure-volume work term [17].

• Thermal Property Extraction: From the equilibrium volume at each temperature, the code computes:

  • Volume-temperature relationship
  • Thermal expansion coefficient
  • Constant-pressure heat capacity ((C_p))
  • Bulk modulus
  • Grüneisen parameter

Protocol: Calculating Harmonic Thermal Properties

This protocol outlines the steps for calculating thermal properties within the harmonic approximation using a force constant potential generated by hiPhive in combination with phonopy [18].

Preparation and Setup

• Define Calculation Parameters:

  • dim: Supercell dimension for the phonopy calculation (e.g., 5 for a 5x5x5 supercell)
  • mesh: q-point mesh for Brillouin zone integration (e.g., [32, 32, 32])
  • temperatures: List of temperatures for evaluating thermal properties (e.g., [300, 600, 900, 1200] K)

• Initialize Phonopy Object:

  • Create a primitive structure (e.g., using ase.build.bulk)
  • Generate a PhonopyAtoms object from the primitive structure
  • Initialize the Phonopy object with the supercell matrix

Force Constant Retrieval

• Load ForceConstantPotential:

  • Read the precomputed force constant potential from file: fcp = ForceConstantPotential.read('filename.fcp')
  • Retrieve the supercell from the phonopy object

• Generate Force Constants:

  • Extract the force constants for the supercell: fcs = fcp.get_force_constants(supercell)
  • The resulting fcs object contains the force constant matrices for all specified orders

Thermal Property Calculation

• Configure Phonopy for Harmonic Analysis:

  • Set the second-order force constants: phonopy.set_force_constants(fcs.get_fc_array(order=2))
  • Set the q-point mesh for Brillouin zone integration

• Compute Mean-Square Displacements (MSD):

  • Configure thermal displacement calculation: phonopy.set_thermal_displacements(temperatures=temperatures)
  • Retrieve and process MSD results: sum over Cartesian directions to get total MSD per atom

• Calculate Phonon Dispersion:

  • Define a path through the Brillouin zone (e.g., Γ-X-K-Γ-L)
  • Compute phonon frequencies along the path: phonopy.set_band_structure(bands)
  • Extract frequencies and q-point vectors for plotting

• Compute Thermodynamic Properties:

  • Use phonopy's built-in methods to calculate Helmholtz free energy, entropy, and heat capacity at constant volume

This protocol provides a robust framework for evaluating harmonic thermal properties, with the entire process typically completing in under a minute for systems of moderate size [18].

Phonopy is an open-source package for calculating phonon properties in crystalline solids at the harmonic and quasi-harmonic levels, playing a critical role in first-principles materials research [6]. By implementing the finite displacement method, it calculates interatomic force constants (IFCs), which are the second-order derivatives of the Born-Oppenheimer potential energy surface with respect to atomic displacements [19] [12]. These IFCs form the foundation of the dynamical matrix, from which key phonon properties such as dispersion relations, density of states (DOS), and various thermodynamic properties are derived [19]. This application note details the protocols for accessing these essential physical properties.

Accessible Physical Properties and Outputs

Phonopy calculates a wide range of phonon-related properties, which can be categorized into structural, dispersive, and thermodynamic properties.

Phonon Dispersions and Band Structure

Phonon dispersion relations describe the vibrational frequency as a function of the wave vector ( \mathbf{q} ) in the Brillouin zone. The dynamical matrix is constructed from the force constants, and its diagonalization yields the squares of the phonon frequencies ( \omega ) and the eigenvectors for each q-point [19]:

[ \sum{j'\beta} D{\alpha\beta}(jj',\mathbf{q}) e\beta(j', \mathbf{q}\nu) = [ \omega(\mathbf{q}\nu) ]^2 e\alpha(j, \mathbf{q}\nu) ]

• Output File: band.yaml [20].
• Contents: This file stores the calculated phonon frequencies along specified band paths in the Brillouin zone. For each q-point, it contains the q-vector in reduced coordinates, the distance from the origin, and for each phonon mode, the frequency and eigenvector [20].
• Visualization: The phonopy-bandplot tool can be used to plot the data stored in band.yaml [20].

Phonon Density of States (DOS) and Projected DOS

The phonon DOS quantifies the number of vibrational states at a specific frequency. The projected DOS (PDOS) further decomposes this quantity per atom, providing insights into the contribution of specific atomic species to the lattice dynamics.

• Output Files:
  • total_dos.dat: Contains the total density of states [20].
  • projected_dos.dat: Contains the projected density of states for atoms in the primitive cell. The first column is the frequency, and the subsequent columns are the PDOS for each atom [20].
• Usage: The DOS is calculated on a uniform mesh of q-points specified by the MESH tag. The --pdos option enables the calculation of projected DOS [4].
• Visualization: The phonopy-pdosplot tool can be used to plot the contents of these files [20].

Thermodynamic Properties

Within the harmonic approximation, phonopy can compute several thermodynamic properties by summing contributions from all phonon modes in the Brillouin zone [19]. The key properties and their statistical thermodynamic expressions are:

• Helmholtz Free Energy (F): [ F = \varphi + \frac{1}{2} \sum{\mathbf{q}\nu} \hbar\omega(\mathbf{q}\nu) + k\mathrm{B} T \sum{\mathbf{q}\nu} \ln \left[1 -\exp(-\hbar\omega(\mathbf{q}\nu)/k\mathrm{B} T) \right] ]
• Entropy (S): [ S = \frac{1}{2T} \sum{\mathbf{q}\nu} \hbar\omega(\mathbf{q}\nu) \coth\left(\frac{\hbar\omega(\mathbf{q}\nu)}{2k\mathrm{B}T}\right) - k\mathrm{B} \sum{\mathbf{q}\nu} \ln\left[2\sinh\left(\frac{\hbar\omega(\mathbf{q}\nu)}{2k_\mathrm{B}T}\right)\right] ]

• Constant-Volume Heat Capacity (Cv): [ CV = \sum{\mathbf{q}\nu} k\mathrm{B} \left( \frac{\hbar\omega(\mathbf{q}\nu)}{k\mathrm{B} T} \right)^2 \frac{\exp(\hbar\omega(\mathbf{q}\nu)/k\mathrm{B} T)}{[\exp(\hbar\omega(\mathbf{q}\nu)/k\mathrm{B} T)-1]^2} ]

• Output File: thermal_properties.yaml [20].

• Units: The default output units are kJ/mol for free energy, and J/K/mol for heat capacity and entropy, where "mol" refers to a mole of the input unit cell [21]. To compare with experimental data per formula unit, these values must be divided by the number of formula units in the calculated unit cell [21].
• Visualization: The phonopy-propplot tool can plot the contents of this file [20].

Table 1: Key Phonopy Output Files and Their Contents

Output File	Property	Content Description	Visualization Tool
band.yaml	Band Structure	Phonon frequencies & eigenvectors on band paths [20]	phonopy-bandplot [20]
mesh.yaml	Mesh Sampling	Phonon frequencies on irreducible q-points of a mesh [20]
thermal_properties.yaml	Thermodynamics	Free energy, entropy, and heat capacity vs. temperature [20]	phonopy-propplot [20]
total_dos.dat	Total DOS	Total phonon density of states [20]	phonopy-pdosplot [20]
projected_dos.dat	Projected DOS	Phonon density of states projected onto atoms [20]	phonopy-pdosplot [20]
phonopy_disp.yaml	Displacements	Information for creating supercells with displacements [20]

Experimental Protocols

Workflow for Basic Phonon Property Calculation

The following protocol outlines the standard procedure for calculating phonon dispersions, DOS, and thermal properties using Phonopy with VASP as the force calculator [4].

Figure 1: A high-level workflow for phonon calculations with Phonopy and VASP, showing the main stages from initial structure to final analysis [4].

Pre-Process: Generate Displacements

• Prepare Input: Start with a fully relaxed unit cell structure file, typically POSCAR.
• Run Phonopy: Execute the following command to create supercells with finite displacements:
  • --dim "2 2 2": Defines the supercell dimension (e.g., 2x2x2) [22] [4].
  • --pa auto: Automatically determines the transformation matrix to a primitive cell, which can reduce computational cost [22].
• Output Files: This step generates:
  • SPOSCAR: The perfect supercell structure.
  • phonopy_disp.yaml: Information about the displacements.
  • POSCAR-001, POSCAR-002, ...: Supercell structure files with individual atomic displacements [4].

Force Calculation

• VASP Setup: Use the generated POSCAR-{number} files as the POSCAR for a series of VASP calculations. The INCAR file must be set for a single-point energy and force calculation without ionic relaxation (e.g., IBRION = -1). Key settings often include [4]:

• Run Calculations: Execute VASP for each displacement file (POSCAR-001, POSCAR-002, ...).

Post-Process: Calculate Properties

• Collect Forces: Create the FORCE_SETS file, which contains the force-displacement data, by parsing the VASP output files (vasprun.xml):
  or using a wildcard pattern [4].
• Run Phonopy for Properties:
  • Band Structure:

  • Density of States and Thermal Properties:
    • --mesh "20 20 20": Defines a 20x20x20 q-point mesh for sampling the Brillouin zone [4].
    • --dos: Calculates the density of states.
    • -t: Calculates thermal properties.
    • -p: Plots the results on the screen. Use -s to save the plot as a PDF [4].
  • Projected DOS:
    This example projects DOS onto atom 1, atom 2 (as one group), and atoms 3,4,5,6 (as another group) [4].

Advanced Application: Quasi-Harmonic Approximation (QHA)

The quasi-harmonic approximation allows for the calculation of thermal expansion and heat capacity at constant pressure by considering the volume dependence of phonon frequencies [17] [12].

Protocol:

• Volume Sampling: Perform full phonon calculations (as in Section 3.1) for a set of different volumes (at least 5 points are recommended) [17].
• Input Files:
  • e-v.dat: A file containing the unit cell volume (in ų) and the corresponding electronic total energy (in eV) at 0K for each volume [17].
  • thermal_properties.yaml-{volume-id}: A series of thermal_properties.yaml files calculated at each volume. These files must be calculated over the same temperature range and step [17].
• Run phonopy-qha: Execute the QHA calculation. To plot the results directly, use:
  Without the -p option, the results are saved to text files without plotting [17].
• Output: The phonopy-qha script outputs several properties, including:
  • Volume vs. Temperature (volume-temperature.*)
  • Thermal expansion coefficient (thermal_expansion.*)
  • Gibbs free energy vs. Temperature (gibbs-temperature.*)
  • Heat capacity at constant pressure (Cp-temperature.*) [17].

The Scientist's Toolkit

Table 2: Essential "Reagents" for a Standard Phonopy-VASP Calculation

Research Reagent / Component	Function / Purpose
Unit Cell Structure (POSCAR)	The initial, fully optimized crystal structure. Serves as the input for generating displaced supercells [4].
Supercell Matrix (DIM tag)	Defines the size of the supercell used for the finite displacement calculation. A larger supercell captures longer-range interactions but increases computational cost [22].
Displacement Distance (DISPLACEMENT_DISTANCE tag)	Sets the amplitude of the finite atomic displacement. The default is 0.01 Å, which is generally suitable for VASP calculations [22].
Force Sets (FORCE_SETS file)	The compiled data file containing the relationship between each atomic displacement and the resulting forces on all atoms, used to compute the force constants [4].
Born Effective Charge & Dielectric Tensor (BORN file)	Required for non-analytical term correction (NAC), which correctly describes the splitting between longitudinal optical (LO) and transverse optical (TO) modes near the Γ-point for polar materials [4].

A Step-by-Step Workflow for Phonopy Calculations: From Setup to Analysis

The finite displacement method (FDM), also known as the frozen phonon approach, is a foundational technique in first-principles lattice dynamics for calculating phonon spectra and related properties of materials [12] [23]. This method relies on constructing a supercell of the original crystal structure and systematically displacing atoms from their equilibrium positions to compute the resulting forces [24]. The core principle involves using these force responses to construct the interatomic force constants (IFCs), which are the derivatives of the Born-Oppenheimer potential energy surface with respect to atomic displacements [12]. The second-order IFCs define the harmonic phonon dispersion relations, while higher-order IFCs govern anharmonic phonon-phonon interactions [12].

The importance of supercells in this context cannot be overstated [24]. A supercell is defined as a crystallographic cell that describes the same crystal as the primitive or conventional unit cell but has a larger volume [24]. For phonon calculations, supercells are essential because they allow for the sampling of force constant matrices that capture the interactions between atoms and their periodic images. The finite displacement method requires supercells large enough to ensure that the force constants between distant atoms decay to zero, thus minimizing the effects of periodic boundary conditions inherent in Density Functional Theory (DFT) calculations [23].

Supercell Generation Methodology

Mathematical Definition and Transformation

The generation of a supercell from an input unit cell is defined by a linear transformation. Given the basis vectors of the unit cell, (( \vec{a}, \vec{b}, \vec{c} )), the basis vectors of the supercell, (( \vec{a}', \vec{b}', \vec{c}' )), are obtained using a transformation matrix (\hat{P}) [24]:

[ ( \vec{a}' \; \vec{b}' \; \vec{c}' ) = ( \vec{a} \; \vec{b} \; \vec{c} ) \hat{P} = ( \vec{a} \; \vec{b} \; \vec{c} ) \begin{pmatrix} P{11} & P{12} & P{13} \ P{21} & P{22} & P{23} \ P{31} & P{32} & P_{33} \end{pmatrix} ]

All elements (P{ij}) of this matrix must be integers, and the determinant of (\hat{P}) must be greater than 1 to create a larger cell [24]. The most straightforward and commonly used transformation is the diagonal supercell expansion, where the matrix (\hat{P}) is diagonal (i.e., (P{i \neq j} = 0)). This represents a simple repetition of the initial cell along its crystallographic axes [24]. For example, a DIM = 2 2 2 tag in Phonopy creates a 2x2x2 supercell using the diagonal matrix [22]:

[ \hat{P} = \begin{pmatrix} 2 & 0 & 0 \ 0 & 2 & 0 \ 0 & 0 & 2 \end{pmatrix} ]

However, non-diagonal matrices are also possible and can be used to generate more complex or specific supercell types, such as transforming a primitive cell into its conventional counterpart [24] [22].

Practical Implementation in Phonopy

The open-source code Phonopy is widely used to automate the supercell generation and displacement creation process [1] [6]. The following provides a detailed protocol for the initial setup.

Table 1: Key Phonopy Configuration Tags for Supercell and Displacement Generation

Tag	Function	Example Usage
DIM	Defines the supercell matrix.	DIM = 2 2 2 (2x2x2 diagonal expansion) [22].
DISPLACEMENT_DISTANCE	Sets the finite displacement amplitude.	DISPLACEMENT_DISTANCE = 0.01 (0.01 Å default) [22].
CREATE_DISPLACEMENTS	Instructs Phonopy to generate supercells with displacements.	CREATE_DISPLACEMENTS = .TRUE. [22].
PRIMITIVE_AXES	(Optional) Defines a transformation to a primitive cell for phonon calculations.	PRIMITIVE_AXES = AUTO [22].

The primary command to generate displaced supercells with Phonopy is [25]:

• The -d argument triggers the creation of displacements.
• The --dim="2 2 2" specifies the supercell dimensions.
• The -c FeSe.struct provides the input unit cell file.
• The --wien2k flag indicates the file format (other codes like VASP use POSCAR by default).

This command generates several files, including a set of supercell structures with displacements (e.g., FeSe.structS-001, FeSe.structS-002, ...) and a phonopy_disp.yaml file that documents the displacements [25]. Each of these generated structure files must then be used as input for independent first-principles force calculations.

Workflow Visualization

The following diagram illustrates the complete workflow from the initial unit cell to the generation of displaced supercells and the subsequent force calculations, which form the foundation for extracting phonon properties.

Diagram Title: Finite Displacement Method Workflow

The Scientist's Toolkit: Essential Research Reagents and Solutions

Table 2: Essential Computational Tools for Finite Displacement Phonon Calculations

Tool / Resource	Function / Purpose	Role in the Protocol
Phonopy	Open-source package for phonon calculations [6].	Automates supercell generation, symmetry-aware displacement creation, and post-processing of forces into phonon properties [22].
DFT Code (e.g., VASP, WIEN2k, Abinit)	First-principles force calculator [26].	Computes the Hellmann-Feynman forces acting on all atoms in each displaced supercell configuration [23].
Optimized Unit Cell	The equilibrium crystal structure.	The essential starting input. Must be fully relaxed (forces and stresses minimized) to ensure accurate force constants [23].
Supercell Matrix ((\hat{P}))	Defines the size and shape of the supercell [24].	Determines the range of interatomic interactions captured and the computational cost. A 2x2x2 matrix is a common starting point [25].

Critical Protocol Considerations

• Structure Relaxation is Prerequisite: The initial unit cell must be fully geometrically relaxed to its equilibrium state (minimized forces and stresses) before generating supercells. Calculations on an unrelaxed structure will lead to inaccurate force constants and phonon spectra [23].
• Supercell Size Convergence: The choice of supercell size is critical. It must be large enough so that the force constants between an atom and its periodic images effectively decay to zero. Failure to use a sufficiently large supercell can introduce unphysical "image" interactions [23]. The required size should be tested for convergence.
• Symmetry and Computational Efficiency: Phonopy leverages crystal symmetry to minimize the number of unique displacements required. This significantly reduces the total number of independent DFT calculations needed, as only symmetrically inequivalent displacements are computed [22].

Running DFT Calculations to Compute Forces on Displaced Atoms

In the broader context of first-principles phonon calculations, computing interatomic force constants (IFCs) is a fundamental step for determining lattice dynamical properties such as phonon band structures, thermodynamic properties, and lattice thermal conductivity [27] [12]. These IFCs are directly related to the forces induced on all atoms in a supercell when one or more atoms are displaced from their equilibrium positions [12]. Within the Born-Oppenheimer approximation, the central quantity is the Born-Oppenheimer potential energy surface, and the IFCs are its derivatives with respect to atomic displacements [12].

The Phonopy code is a powerful and widely used tool that facilitates these calculations by automating the generation of systematically displaced supercells, orchestrating the required force calculations using an external Density Functional Theory (DFT) code, and post-processing the results to produce phonon properties [28] [27]. This application note provides a detailed protocol for performing these critical force computations, which serve as the foundation for reliable and accurate phonon simulations within a thesis research framework.

The general process for calculating forces on displaced atoms to obtain phonon properties with Phonopy involves a sequence of well-defined steps, connecting the initial crystal structure to the final phonon band structure. The workflow, integrating DFT calculations with Phonopy, is summarized in the following diagram:

Diagram Title: Workflow for Phonon Calculations with Phonopy and DFT.

Experimental Protocol

The following table outlines the detailed, step-by-step methodology for computing forces on displaced atoms, tailored for a VASP-Phonopy workflow [28]. The process can be adapted to other DFT codes (e.g., Quantum ESPRESSO [27]) or semi-empirical methods (e.g., DFTB+ [29]) with appropriate modifications to the file formats and commands.

Table: Detailed Protocol for Force Calculations with VASP and Phonopy

Step	Action	Key Commands / Inputs	Outputs & Purpose
1. Structure Preparation	Prepare a POSCAR file for the fully relaxed unit cell of the material under investigation.	POSCAR-unitcell	A POSCAR file containing the optimized ground-state structure.
2. Supercell Generation	Use Phonopy to create a set of supercells with finite atomic displacements.	phonopy -d --dim "2 2 2" -c POSCAR-unitcell	SPOSCAR (undisplaced supercell) and multiple POSCAR-001... files (displaced supercells). The --dim flag defines the supercell expansion.
3. DFT Calculation Setup	Prepare INCAR, KPOINTS, and POTCAR files for VASP. Key settings must ensure accurate force evaluation.	INCAR:PREC = AccurateENCUT = [High Enough]EDIFF = 1E-8IBRION = -1 (or 8 for DFPT)ISIF = 2NSW = 0LREAL = .FALSE.ISMEAR = 0; SIGMA = 0.1LWAVE = .FALSE.LCHARG = .FALSE. [28]	A set of VASP input files configured for single-point energy and force calculations on the displaced structures.
4. Force Computation	Run a VASP calculation for each displaced POSCAR file.	For each i in [001, 002, ...]:cp POSCAR-$i POSCARmpirun -np [N] vasp_std	vasprun.xml files containing the forces on every atom for each displacement.
5. Force Set Creation	Collect all the force results and compile them into Phonopy's FORCE_SETS file.	phonopy -f disp-*/vasprun.xml	FORCE_SETS file. This file is the critical input for Phonopy that maps displacements to calculated forces.
6. Phonon Property Calculation	Run Phonopy to compute harmonic properties like the phonon band structure.	phonopy --dim="2 2 2" -c POSCAR-unitcell --band="[Q-Path]"phonopy-bandplot --gnuplot band.yaml > band.dat	band.yaml and band.dat files containing the phonon band structure data for visualization.

Convergence and Best Practices

The accuracy of the computed forces, and consequently the phonon properties, is highly sensitive to the convergence parameters of the underlying DFT calculations [30]. Inaccurate forces can lead to significant errors, a problem notably identified in several large datasets used for machine learning interatomic potentials [31]. The following table summarizes key convergence parameters and their impact on the calculated forces.

Table: Essential DFT Convergence Parameters for Accurate Forces

Parameter	Description	Convergence Strategy & Impact on Forces
Plane-Wave Cutoff (ENCUT)	Kinetic energy cutoff for the plane-wave basis set.	Systematically increase ENCUT until the total energy and forces change by less than a target error (e.g., 1 meV/atom). An insufficient cutoff introduces a systematic error [30].
k-Point Mesh (KPOINTS)	Sampling of the Brillouin zone.	Use a dense enough k-mesh. Convergence should be checked for the supercell, which often requires a less dense k-mesh than the primitive cell due to its larger real-space size.
Force Convergence (EDIFF)	Stopping criterion for the electronic SCF cycle.	Use a tight criterion (e.g., EDIFF = 1E-7 or 1E-8) to ensure the electronic wavefunction is fully converged, which is necessary for accurate Hellmann-Feynman forces [32] [28].
Numerical Settings	Includes integration grids and exchange-correlation quadrature.	In codes like ORCA, using the DEFGRID3 keyword and disabling approximations like RIJCOSX can be critical for reducing force errors to below 1 meV/Å [31].
SCF Algorithm (ALGO)	Algorithm for wavefunction optimization.	For difficult-to-converge systems (e.g., metals), switching from the default DIIS algorithm to blocked-Davidson or conjugate gradient (ALGO = All or Damped in VASP) can improve stability [32].

Troubleshooting Common Problems

• SCF Non-Convergence: If the self-consistent field (SCF) cycle fails to converge, consider increasing the number of SCF steps (NELM in VASP), using a different mixing parameter, employing smearing for metallic systems (ISMEAR), or using a pre-converged wavefunction from a previous calculation as an initial guess [32].
• Geometry Optimization Problems: Before generating displacements, ensure the initial unit cell is fully relaxed. A poorly converged geometry can lead to imaginary phonon frequencies. Monitor the forces in the optimization to ensure they are sufficiently small (e.g., < 0.001 eV/Å) [32] [33].

The Scientist's Toolkit

Table: Key Research Reagent Solutions for DFT Phonon Calculations

Tool / Resource	Function in the Workflow
Phonopy Code	The primary orchestrator: generates atomic displacements and post-processes DFT forces to calculate phonon frequencies and dispersion [28].
DFT Software (VASP, Quantum ESPRESSO)	The "energy calculator": computes the total energy and, crucially, the Hellmann-Feynman forces for each displaced atomic configuration generated by Phonopy [28] [27].
Pseudopotential Library	Provides the atomic potential data (POTCAR for VASP). High-quality, consistent pseudopotentials are essential for transferable and accurate results [30].
High-Performance Computing (HPC) Cluster	Provides the computational power needed to run dozens to hundreds of independent DFT force calculations in a parallel high-throughput manner.
Tight Force Convergence Settings	A specific set of computational parameters (e.g., EDIFF=1E-8, high ENCUT) that acts as a "reagent" to ensure the output forces have minimal numerical error, which is critical for dataset quality [31].
Phono3py Code	An extension for calculating third-order force constants and anharmonic properties, such as lattice thermal conductivity, following a workflow very similar to Phonopy but requiring a larger number of displacements [27].

Phonopy is an open-source package for phonon calculations at harmonic and quasi-harmonic levels, serving as an essential tool in computational materials science [6]. The post-processing stage in Phonopy begins after the successful completion of force calculations on systematically displaced supercells. This phase involves extracting the interatomic force constants (IFCs) and computing fundamental phonon properties, transforming raw force data into physically meaningful quantities. The IFCs represent the second-order derivatives of the Born-Oppenheimer potential energy surface with respect to atomic displacements, forming the foundation of the harmonic lattice dynamical model [12]. Within the context of a broader thesis on first-principles phonon calculations, mastering Phonopy's post-processing workflow is crucial for investigating material properties such as dynamical stability, thermodynamic behavior, and vibrational spectra.

The post-processing workflow encompasses several key stages: generating force sets from calculated forces, constructing the dynamical matrix from the IFCs, and finally calculating and visualizing phonon properties. This process enables researchers to bridge the gap between atomic-scale interactions and macroscopic observable properties. The accuracy of these predictions hinges on the careful execution of each post-processing step and the correct interpretation of the resulting data.

Theoretical Foundations: From IFCs to Phonon Properties

Interatomic Force Constants and the Dynamical Matrix

Within the Born-Oppenheimer and harmonic approximations, the potential energy of a crystal can be expressed as a Taylor expansion in terms of atomic displacements. The IFCs are defined as the second-order force constants:

[ \Phi{\alpha\beta}(\kappa\kappa',l'l) = \frac{\partial^2 E}{\partial u{\alpha}(l\kappa)\partial u_{\beta}(l'\kappa')} ]

where ( E ) is the potential energy, ( u{\alpha}(l\kappa) ) is the displacement of atom ( \kappa ) in unit cell ( l ) along the Cartesian direction ( \alpha ), and ( \Phi{\alpha\beta}(\kappa\kappa',l'l) ) represents the force constant matrix [12]. These IFCs are fundamental building blocks that encode information about the restoring forces when atoms are displaced from their equilibrium positions.

The dynamical matrix ( D(\mathbf{q}) ) is then constructed via Fourier transformation of the IFCs:

[ D{\alpha\beta}(\kappa\kappa'|\mathbf{q}) = \frac{1}{\sqrt{m\kappa m{\kappa'}}} \sum{l'} \Phi_{\alpha\beta}(\kappa\kappa',0l') e^{i\mathbf{q}\cdot(\mathbf{r}(l')+\mathbf{r}(\kappa')-\mathbf{r}(\kappa))} ]

where ( m\kappa ) is the mass of atom ( \kappa ), ( \mathbf{q} ) is the wave vector in the Brillouin zone, and ( \mathbf{r} ) denotes atomic positions. Diagonalization of the dynamical matrix for each wave vector ( \mathbf{q} ) yields the square of the phonon frequencies ( \omega{\mathbf{q}j} ) and the corresponding polarization vectors, where ( j ) indexes the phonon branches [12].

Phonon Properties and Their Physical Significance

The computed phonon frequencies enable the calculation of various important physical properties:

• Phonon Dispersion Relations: The dependence of phonon frequencies on the wave vector ( \mathbf{q} ) along high-symmetry paths in the Brillouin zone provides crucial information about lattice dynamics, including acoustic and optical branches, and can reveal soft modes indicative of structural instabilities.
• Phonon Density of States (DOS): The phonon DOS, ( g(\omega) ), represents the number of vibrational states per unit frequency range at a given frequency ( \omega ), offering insights into the distribution of vibrational modes and their contributions from different atomic species.
• Thermodynamic Properties: Within the harmonic approximation, phonons contribute significantly to various thermodynamic properties. The Helmholtz free energy ( F ), entropy ( S ), and constant-volume heat capacity ( C_v ) can be derived from the phonon frequencies using statistical mechanics relations [21].

Table 1: Fundamental phonon properties calculable in Phonopy post-processing

Property	Physical Interpretation	Computational Utility
Phonon Dispersion	Relationship between phonon frequency and wave vector in reciprocal space	Identifies lattice dynamical stability, acoustic branches, and soft modes
Phonon DOS	Distribution of vibrational states as a function of frequency	Reveals contributions from different atomic species; identifies gap modes
Partial DOS	Projection of DOS onto specific atoms or groups	Quantifies contribution of specific atomic species to vibrational spectra
Thermal Properties	Temperature-dependent free energy, entropy, heat capacity	Predicts thermodynamic stability and behavior at finite temperatures

Computational Methodology: Phonopy Post-Processing Workflow

Generating the Force Sets

The post-processing workflow begins after performing VASP calculations on the displaced supercells generated by phonopy -d. To create the FORCE_SETS file, which contains the displacement-force data, navigate to the directory containing the VASP output files (vasprun.xml) and run:

Alternatively, use Bash expansion for efficiency:

This command reads the displacement information from phonopy_disp.yaml and the computed forces from each vasprun.xml file, then generates the FORCE_SETS file [4]. This file contains the essential data mapping atomic displacements to the resulting forces, from which the IFCs will be extracted.

Calculating Phonon Properties

With the FORCE_SETS file prepared, various phonon properties can be computed:

Phonon Density of States:

The --mesh flag specifies a 20×20×20 q-point mesh for sampling the Brillouin zone. For plotting the DOS:

Projected DOS: To decompose the DOS by atomic species:

This projects the DOS onto specified atoms (e.g., atoms 1-2 for one species, atoms 3-6 for another) [4].

Phonon Band Structure:

This command calculates the phonon dispersion along the specified reciprocal space path, with the example showing a path through high-symmetry points [4].

Thermal Properties:

This computes temperature-dependent thermal properties including free energy, entropy, and heat capacity [21]. The -p option can be added to plot the results.

Figure 1: Workflow for Phonopy post-processing showing the sequence from force calculations to phonon properties.

Advanced Post-Processing: Special Considerations

Non-Analytical Term Correction for Polar Materials

For polar materials, the long-range dipole-dipole interactions induce LO-TO splitting at the Brillouin zone center, which must be properly accounted for to obtain accurate phonon frequencies. This requires calculating the Born effective charges and dielectric tensor [34] [4].

First, perform a DFPT calculation in the unit cell with the following VASP INCAR settings:

After obtaining the dielectric tensor and Born effective charges from this calculation, create the BORN file using the phonopy-vasp-born utility:

For symmetry-constrained values, use:

The resulting BORN file should be placed in the working directory, and Phonopy will automatically include non-analytical term correction in subsequent calculations [4].

Convergence Testing

Ensuring converged results is crucial for reliable phonon calculations. Key parameters requiring convergence tests include:

• q-point mesh density for DOS and thermal properties calculations
• Supercell size for force calculations
• Energy cutoffs and k-point meshes in DFT calculations

Systematically increase the mesh density until the phonon frequencies and thermodynamic properties show negligible changes. The convergence of the supercell size can be monitored by checking if the IFCs have decayed to negligible values at the supercell boundaries [34].

Table 2: Key parameters for convergence testing in phonon calculations

Parameter	Physical Significance	Convergence Strategy
q-mesh density	Sampling quality in reciprocal space	Increase until DOS/thermal properties stabilize (~20×20×20 for complex structures)
Supercell size	Range of captured interatomic interactions	Increase until IFCs decay to zero at boundaries (typically 3×3×3 to 5×5×5)
ENCUT	Plane-wave basis set completeness	Increase until total energy changes by <1 meV/atom
k-point mesh	Brillouin zone sampling for electronic structure	Increase until phonon frequencies change by <0.1 THz

Data Interpretation and Analysis

Analyzing Phonon Dispersion and DOS

The phonon dispersion relations provide immediate insights into lattice dynamical stability. Imaginary frequencies (negative values in plots) indicate structural instabilities, often preceding phase transitions. The slope of acoustic branches near the Γ-point relates to elastic constants, with longitudinal acoustic (LA) modes typically having higher velocities than transverse acoustic (TA) modes.

The phonon DOS reveals the distribution of vibrational states, with low-frequency regions dominated by heavier atoms and high-frequency regions by lighter atoms. Gaps in the DOS indicate frequency ranges with no vibrational states, often reflecting large mass differences or specific bonding patterns. Projected DOS (PDOS) analysis enables identification of specific atomic contributions to particular vibrational features, which is particularly valuable for complex materials with multiple atomic species.

Interpreting Thermal Properties

Phonopy computes thermal properties using standard statistical mechanical relations for a harmonic crystal:

• Helmholtz free energy: ( F(T) = kB T \int0^\infty \ln\left[2\sinh\left(\frac{\hbar\omega}{2k_B T}\right)\right] g(\omega)d\omega )
• Entropy: ( S(T) = kB \int0^\infty \left[\frac{\hbar\omega}{2kB T} \coth\left(\frac{\hbar\omega}{2kB T}\right) - \ln\left(2\sinh\left(\frac{\hbar\omega}{2k_B T}\right)\right)\right] g(\omega)d\omega )
• Heat capacity: ( Cv(T) = kB \int0^\infty \left(\frac{\hbar\omega}{2kB T}\right)^2 \mathrm{csch}^2\left(\frac{\hbar\omega}{2k_B T}\right) g(\omega)d\omega )

where ( g(\omega) ) is the phonon density of states [21].

Important note on units: Phonopy reports thermal properties per unit cell. For comparison with experimental data, these values must be normalized appropriately. As noted in the Phonopy documentation, "The unit systems of free energy, heat capacity, and entropy are kJ/mol, J/K/mol, and J/K/mol, respectively, where 1 mol means ( N_A ) × your input unit cell (not formula unit)" [21]. For example, in a conventional MgO unit cell with 4 formula units, divide Phonopy outputs by 4 to obtain values per formula unit.

Figure 2: Relationship between phonon frequencies and thermal properties in Phonopy, showing the normalization requirement.

Research Applications and Protocol Implementation

Case Study: Lattice Dynamics in Energy Materials

Recent research has demonstrated the powerful combination of Phonopy with machine learning approaches for discovering materials with tailored phonon properties. In a groundbreaking study on sodium superionic conductors, researchers analyzed 3903 Na-containing structures and established a strong correlation between phonon mean squared displacement (MSD) of Na+ ions and their diffusion coefficients [35].

The key lattice dynamics signatures identified for superionic conductivity include:

• Low acoustic cutoff phonon frequencies
• Suppressed Na+ vibrational density of states slightly above the acoustic cutoff
• Enhanced low-frequency vibrational coupling between Na+ ions and the host sublattice

This research incorporated phonon-derived descriptors into machine learning frameworks to rapidly predict ionic transport properties, accelerating the discovery of high-performance solid electrolytes for all-solid-state sodium batteries [35].

Integration with Advanced Phonon Codes

The phonon code ecosystem has expanded significantly, with specialized packages now available for specific applications:

• Phono3py: For third-order IFCs and lattice thermal conductivity calculations [36]
• Pheasy: Employs machine-learning algorithms for efficient extraction of high-order IFCs [12]
• ALAMODE: Implements compressive sensing lattice dynamics for anharmonic calculations

These tools can complement Phonopy for investigating anharmonic properties beyond the harmonic approximation, enabling studies of thermal conductivity, phonon-phonon interactions, and temperature-dependent phonon renormalization.

Research Reagent Solutions: Computational Tools

Table 3: Essential computational tools for phonon calculations and their functions

Tool/Software	Primary Function	Application in Phonon Research
Phonopy	Harmonic phonon calculations	Calculating phonon dispersion, DOS, and thermodynamic properties
VASP	First-principles electronic structure calculations	Generating force constants via finite-displacement or DFPT methods
Phono3py	Third-order force constants and phonon-phonon interactions	Computing lattice thermal conductivity and anharmonic properties
Pheasy	High-order IFC extraction using machine learning	Studying strongly anharmonic systems and high-temperature properties
ShengBTE	Solving Boltzmann transport equation	Calculating thermal conductivity from second- and third-order IFCs

Troubleshooting and Best Practices

Common Issues and Solutions

Imaginary Frequencies: The appearance of significant imaginary frequencies in phonon spectra may indicate:

• Structural instabilities or unresolved symmetry issues
• Inadequate supercell size or insufficient k-point sampling
• Need for higher accuracy in force calculations

Memory Issues with Large Supercells: For large systems, Phonopy calculations may require substantial memory. Consider:

• Using a coarser q-mesh for initial calculations
• Increasing system memory or using computational clusters
• Utilizing Phonopy's HDF5 support for efficient data handling

LO-TO Splitting Artifacts: Improper treatment of non-analytical term correction may lead to unphysical phonon dispersions in polar materials. Ensure:

• Accurate calculation of Born effective charges and dielectric tensor
• Proper convergence of the DFPT calculation for dielectric properties
• Correct formatting of the BORN file

Validation Protocols

To ensure the reliability of phonon calculations:

• Compare with experimental data when available (Raman spectra, infrared spectroscopy, inelastic neutron scattering)
• Verify IFC convergence with respect to supercell size
• Check acoustic sum rules - acoustic modes should approach zero at the Γ-point
• Validate against known results for standard materials before studying new systems

Following these post-processing protocols systematically will ensure robust and reproducible phonon calculations, contributing valuable insights to materials research and facilitating the discovery of materials with tailored thermal and vibrational properties.

The Quasi-Harmonic Approximation (QHA) serves as a pivotal extension to the standard harmonic phonon model in solid-state physics, enabling the calculation of volume-dependent thermal properties that are fundamental to materials science and drug development. While the harmonic approximation assumes phonon frequencies are independent of volume and thus cannot describe thermal expansion, QHA introduces volume dependence to phonon frequencies, thereby incorporating a crucial anharmonic effect while maintaining much of the harmonic formalism's computational tractability [37]. This approach has become the workhorse method for predicting thermodynamic properties under arbitrary temperature and pressure conditions using density-functional theory (DFT) [38].

Within the context of first-principles phonon calculations using the Phonopy code, QHA provides the theoretical foundation for determining essential material characteristics including thermal expansion coefficients, heat capacities at constant pressure, and temperature-dependent free energies [17]. These properties are indispensable for researchers investigating material stability, phase transitions, and thermodynamic behavior across temperature ranges—particularly relevant for pharmaceutical development where crystal structure stability affects drug efficacy and shelf life.

Theoretical Foundation

Basic Principles of QHA

The quasi-harmonic approximation builds upon the harmonic phonon model by maintaining the mathematical framework of harmonic lattice dynamics while introducing volume dependence to phonon frequencies. In essence, QHA assumes that the harmonic approximation remains valid for every value of the lattice constant, which is treated as an adjustable parameter [37]. This fundamental extension allows the model to capture anharmonic effects responsible for thermal expansion, which pure harmonic models cannot describe since their equilibrium atomic distances remain temperature-independent.

The core physical insight of QHA is that phonon frequencies (ω) become volume-dependent (ω(V)), creating an implicit temperature dependence through the volume expansion (ω(V(T))) [38]. This relationship enables the approximation to partially incorporate anharmonicity while avoiding the computational complexity of fully anharmonic methods. The volume dependence is typically determined by calculating phonon spectra at multiple volumes, then fitting the results to establish the ω(V) relationship.

Thermodynamic Formalism

In the quasi-harmonic approximation, the Helmholtz free energy F of a crystal lattice is expressed as a function of temperature T and volume V:

$$ F(T,V) = E{\text{lat}}(V) + U{\text{vib}}(T,V) - TS(T,V) $$

where $E{\text{lat}}(V)$ represents the static internal lattice energy, $U{\text{vib}}(T,V)$ is the internal vibrational energy of the phonon system, and S is the vibrational entropy [37].

The vibrational energy component is given by:

$$ U{\text{vib}}(T,V) = \frac{1}{N} \sum{\mathbf{k},i} \left[ \frac{1}{2} + n{\mathbf{k},i}(T,V) \right] \hbar \omega{\mathbf{k},i}(V) $$

where $\omega{\mathbf{k},i}(V)$ are the volume-dependent phonon frequencies, and $n{\mathbf{k},i}(T,V)$ represents the Bose-Einstein distribution function for phonon occupation numbers [37].

For practical applications at constant pressure, a Legendre transformation converts the Helmholtz free energy to the Gibbs free energy:

$$ G(T,p) = \minV \left[ E{\text{lat}}(V) + F_{\text{phonon}}(T;V) + pV \right] $$

where the minimization over volume is performed at each temperature and pressure point to find the equilibrium volume $V(T,p)$ [17].

The Grüneisen Parameter

A central concept in QHA is the Grüneisen parameter, which quantifies the volume dependence of phonon frequencies and provides a direct measure of anharmonicity:

$$ \gammai = -\frac{\partial \ln \omegai}{\partial \ln V} $$

The Grüneisen parameter directly links phonon frequency shifts to volumetric strain and plays a crucial role in determining thermal expansion behavior. The total thermodynamic Grüneisen parameter is calculated as a weighted sum over all phonon modes [37].

Table: Key Thermodynamic Relationships in QHA

Property	Mathematical Expression	Physical Significance
Helmholtz Free Energy	$F(T,V) = E{\text{lat}}(V) + \frac{1}{N}\sum{\mathbf{k},i} \frac{1}{2}\hbar \omega{\mathbf{k},i}(V) + \frac{1}{N}\sum{\mathbf{k},i} kB T \ln \left[1 - \exp(-\Theta{\mathbf{k},i}(V)/T)\right]$	Fundamental thermodynamic potential at constant volume
Gibbs Free Energy	$G(T,p) = \min_V \left[ F(T,V) + pV \right]$	Fundamental thermodynamic potential at constant pressure
Thermal Expansion Coefficient	$\alphaV = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)p$	Measures volume change with temperature at constant pressure
Heat Capacity at Constant Volume	$CV = \frac{1}{N}\sum{\mathbf{k},i} kB \frac{(\Theta{\mathbf{k},i}/T)^2 \exp(\Theta{\mathbf{k},i}/T)}{[\exp(\Theta{\mathbf{k},i}/T) - 1]^2}$	Temperature derivative of internal energy at constant volume

Computational Protocols

Implementing QHA calculations requires a systematic approach spanning multiple computational stages. The following diagram illustrates the comprehensive workflow from initial structure preparation to final thermodynamic analysis:

Volume Sampling and Electronic Structure Calculations

The initial stage involves calculating the total energy and phonon properties across a range of volumes:

• Volume Grid Selection: Choose at least 5-7 volume points centered around the equilibrium volume, typically spanning ±5-10% of V₀. More volumes may be needed for systems with complex energy landscapes [17].

• Electronic Structure Calculations: For each volume, perform the following:

  • Structural relaxation while maintaining fixed volume
  • Single-point energy calculation to obtain Eₗₐₜ(V)
  • DFPT or finite-displacement calculation to determine phonon frequencies

  Example VASP INCAR parameters for force calculations:

  Critical: Ensure structural relaxation is disabled (IBRION = -1) during force calculations to maintain prescribed atomic displacements [4].

Phonon Calculations with Phonopy

The phonon properties are calculated using Phonopy through these specific steps:

• Pre-processing and Supercell Generation:

  This generates supercells with displacements (POSCAR-{number}), the perfect supercell (SPOSCAR), and displacement information (phonopy_disp.yaml) [4].

• Force Calculations and Force Sets:

  • Run VASP calculations for each POSCAR-{number} to compute atomic forces
  • Create FORCE_SETS file by parsing the vasprun.xml files:

• Thermal Property Calculation:

  This computes harmonic thermal properties at constant volume across a 20×20×20 q-point mesh, generating thermal_properties.yaml files containing Helmholtz free energy, entropy, and heat capacity at constant volume across a temperature range [4].

QHA Analysis with phonopy-qha

The final stage combines electronic and vibrational contributions to determine constant-pressure properties:

• Prepare Input Files:

  • e-v.dat: Volume-energy data (electronic total energies at 0K)
  • Multiple thermal_properties.yaml files: One for each volume point

• Run QHA Analysis:

  The -p flag generates plots of the results [17].

• Key Analysis Steps:

  • At each temperature, fit Helmholtz free energy F(T,V) to an equation of state (default: Vinet)
  • Minimize Gibbs energy G(T,p) = F(T,V) + pV with respect to volume
  • Extract equilibrium volume V(T,p) and compute derived properties

Data Analysis and Interpretation

Output Files and Properties

The phonopy-qha script generates multiple output files containing calculated thermodynamic properties:

Table: Key Output Files from phonopy-qha

Output File	Content	Physical Units
volume-temperature.*	Equilibrium volume as function of temperature	ų
thermal_expansion.*	Volumetric thermal expansion coefficient β vs. T	K⁻¹
bulk_modulus-temperature.*	Bulk modulus B_T vs. temperature	GPa
gibbs-temperature.*	Gibbs free energy vs. temperature	eV
Cp-temperature.*	Heat capacity at constant pressure vs. T	J/K/mol
gruneisen-temperature.dat	Thermodynamic Grüneisen parameter vs. T	Unitless

Practical Example: Silicon QHA Calculation

A specific example for silicon QHA calculations is provided in the example/Si-QHA directory of Phonopy. The step-by-step protocol includes:

• Volume-Energy Data Preparation (e-v.dat):

• Thermal Properties Calculation: Phonopy calculations at each volume generate thermal_properties.yaml files containing temperature-dependent Helmholtz free energies.

• QHA Execution:

• Results Interpretation: The script calculates and plots the thermal expansion coefficient showing positive expansion for silicon, heat capacity approaching the Dulong-Petit limit at high temperature, and smoothly varying bulk modulus with temperature.

Advanced Options and Considerations

The phonopy-qha implementation includes several important options for specialized applications:

• Equation of State Selection: --eos option allows choosing between vinet (default), birch_murnaghan, and murnaghan equations of state [17]
• Pressure Effects: --pressure option includes PV terms for non-ambient pressure conditions
• Electronic Free Energy: --efe option incorporates temperature-dependent electronic free energies beyond the phonon contribution, requiring fe-v.dat with electronic free energies at multiple volumes and temperatures

Research Reagent Solutions

Table: Essential Computational Tools for QHA Calculations

Tool/Code	Function	Application Context
Phonopy	Primary code for phonon calculations	Harmonic and quasi-harmonic lattice dynamics
VASP	First-principles electronic structure	Force calculations for force constants
phonopy-qha	QHA analysis script	Thermodynamic property calculation at constant pressure
Phonopy-VASP-Born	Auxiliary tool for non-analytical term correction	LO-TO splitting in polar materials
Phonopy-VASP-efe	Electronic free energy calculation	Temperature-dependent electronic contributions

Limitations and Advanced Methodologies

QHA Limitations

While QHA represents a significant advancement over the harmonic approximation, researchers should recognize its inherent limitations:

• High-Temperature Breakdown: QHA produces spurious predictions for volume and other properties at sufficiently high temperatures due to its approximate treatment of anharmonicity [38]. The method's temperature validity limit depends on the specific system and pressure, with systems containing soft phonons or light atoms being particularly affected.

• Dynamic Instability: QHA cannot model dynamically stabilized phases that exhibit imaginary frequencies at zero temperature but become stable at elevated temperatures due to anharmonic effects [38].

• Implicit Temperature Dependence: The approximation includes temperature dependence only through volume expansion (ω(V(T))), which is correct only to first order in anharmonicity [38].

Beyond QHA: Quasiparticle Theory

Recent advancements address QHA limitations through approaches combining quasiparticle (QP) theory with the self-consistent harmonic approximation (SCHA). This methodology:

• Calculates temperature-dependent effective harmonic frequencies ω(V,T) using SCHA
• Employs Allen's quasiparticle theory to compute anharmonic entropy from effective frequencies
• Uses a Debye-like numerical model to compute all thermodynamic properties from QP entropies
• Coincides with QHA at low temperature while eliminating spurious high-temperature blowout [38]

This approach maintains computational complexity similar to QHA but requires additional supercell calculations for higher-order force constants. Implementation examples for materials like MgO and CaO demonstrate successful prediction of thermodynamic properties up to melting temperatures, substantially improving upon QHA limitations [38].

The Quasi-Harmonic Approximation implemented through the Phonopy code provides researchers with a powerful, computationally efficient method for predicting thermodynamic properties of materials across temperature and pressure ranges. While the approach has limitations at high temperatures and for dynamically stabilized phases, its systematic protocol for calculating thermal expansion, heat capacities, and free energies makes it invaluable for materials research and pharmaceutical development. Emerging methodologies that incorporate more sophisticated treatments of anharmonicity promise to extend these capabilities while maintaining the foundational framework established by QHA.

The discovery and development of high-performance superionic conductors are crucial for advancing next-generation energy storage technologies, particularly all-solid-state batteries. These materials exhibit remarkably high ionic conductivity (σ ≈ 1 mS cm⁻¹ or above at room temperature), comparable to liquids, while maintaining solid-state stability. Traditional discovery approaches have relied heavily on insights from material defect chemistry and hopping theory, often overlooking the fundamental role of lattice vibrations (phonons) in ion transport mechanisms. Recent computational studies have revealed that phonon mean squared displacement (MSD) provides a quantitative correlation between lattice dynamics and ionic transport properties, serving as a powerful descriptor for identifying promising superionic conductors.

Within the framework of first-principles phonon calculations using the Phonopy code, researchers can now efficiently compute phonon MSD and related lattice dynamics signatures to screen for materials with inherent superionic potential. This approach has demonstrated particular value for sodium superionic conductors (NASICONs), which are gaining significant interest as solid-state electrolytes for all-solid-state sodium batteries due to their greater abundance, lower cost, and improved safety compared to lithium-based systems. The phonon MSD descriptor effectively captures the anharmonic effects and lattice softness that promote ionic diffusion, enabling computational screening of large material databases without resource-intensive molecular dynamics simulations for every candidate.

Theoretical Foundation: Linking Lattice Dynamics to Ionic Transport

Phonon Mean Squared Displacement Fundamentals

The phonon mean squared displacement represents the average square displacement of atoms from their equilibrium positions due to thermal vibrations within the lattice. In the context of superionic conductors, the MSD of the mobile ion species (e.g., Na⁺, Li⁺, Cu⁺) provides crucial insights into their potential migration behavior. Mathematically, the MSD can be derived from the phonon mode eigenvectors and frequencies, incorporating both harmonic and anharmonic contributions to atomic vibrations. For superionic conduction, the low-frequency acoustic and optic modes have been identified as particularly significant, as these modes dominate the large-amplitude vibrations that enable ion hopping between lattice sites.

First-principles calculations within the density functional theory (DFT) framework enable the computation of MSD through the interatomic force constants (IFCs), which are the derivatives of the Born-Oppenheimer potential energy surface with respect to atomic displacements. The Phonopy code efficiently calculates these IFCs using the finite displacement method, providing the foundation for MSD analysis. Recent advances have extended these calculations beyond the harmonic approximation through methods like the anharmonic special displacement method (ASDM), which better captures the strong anharmonicity and temperature-dependent effects characteristic of superionic conductors.

Correlation Between MSD and Ionic Diffusion

Substantial research has established a direct relationship between phonon MSD and ionic diffusion coefficients in superionic materials. A comprehensive study screening 3903 Na-containing structures demonstrated a strong positive correlation between phonon MSD of Na⁺ ions and their diffusion coefficients [35]. This correlation emerges because large MSD values indicate significant vibrational amplitudes that lower the effective energy barrier for ion migration between adjacent lattice sites. The MSD effectively captures the preferential vibrational pathways that precede actual ion jumps, serving as a proxy for ionic mobility without requiring explicit modeling of the diffusion process.

The connection between MSD and ionic conductivity can be understood through the Einstein relationship, which relates diffusivity (D) to mean square displacement: D = (1/6t)⟨ΔR²(t)⟩, where ⟨ΔR²(t)⟩ represents the ensemble-averaged mean square displacement over time t. This relationship forms the basis for computing diffusion coefficients from molecular dynamics simulations and establishes the theoretical foundation for using MSD as a predictive descriptor for superionic behavior. In superionic conductors with strong anharmonicity, the phonon lifetime and frequency renormalization effects further enhance this correlation, as these materials often exhibit overdamped vibrations that facilitate ionic motion.

Computational Workflow for Phonon MSD Screening

High-Throughput Screening Methodology

The integration of phonon MSD analysis into high-throughput computational screening pipelines enables efficient identification of promising superionic conductors from large material databases. The established workflow begins with structure selection from comprehensive databases such as the Open Quantum Materials Database (OQMD) and Inorganic Crystal Structure Database (ICSD), applying initial filters for sodium-containing structures with non-zero band gaps and negative formation energies to ensure thermodynamic stability. For the filtered structures, DFT-based structural optimization is performed with strict computational parameters to obtain accurate ground-state configurations.

Following structural optimization, lattice dynamics properties are calculated using the Phonopy code, which computes harmonic and anharmonic IFCs through the finite displacement method. For dynamically stable structures (those free of imaginary frequencies across the Brillouin zone), phonon MSD values are calculated specifically for the mobile ion species. These values are then correlated with diffusion coefficients to establish a quantitative relationship between MSD and ionic conductivity. Machine learning models can be incorporated at this stage, using phonon-derived descriptors to rapidly predict ionic transport properties across broad structural spaces, significantly accelerating the discovery process [35].

Workflow Diagram

The following diagram illustrates the integrated computational workflow for high-throughput screening of superionic conductors using phonon MSD:

Quantitative Relationships: Key Lattice Dynamics Descriptors

Correlation Between Phonon Properties and Ionic Conductivity

Extensive screening of sodium-containing structures has revealed several key lattice dynamics descriptors that strongly correlate with superionic conductivity. The table below summarizes these relationships based on analysis of 921 dynamically stable structures from a initial set of 3903 Na-containing compounds:

Table 1: Key Lattice Dynamics Descriptors Correlated with Superionic Conductivity

Descriptor	Correlation with Ionic Conductivity	Physical Significance	Optimal Range
Phonon MSD of Na⁺ ions	Strong positive correlation [35]	Measures vibrational amplitude preceding ion migration	Higher values preferred
Acoustic cutoff frequency	Negative correlation [35]	Indicates lattice softness	Lower values preferred
Center of Na⁺ VDOS	Moderate negative correlation [35]	Reflects vibrational energy distribution	Slightly above acoustic cutoff
Low-frequency phonon coupling	Positive correlation [35]	Enhances host-mobile ion energy transfer	Enhanced coupling preferred

The phonon MSD emerges as the most significant predictor, with studies demonstrating that only a small subset of low-frequency acoustic and optic modes contribute dominantly to large phonon MSDs and Na⁺ ion migration, while higher-energy modes contribute negligibly. This selective contribution enables efficient screening by focusing computational resources on characterizing the most relevant portions of the phonon spectrum.

Material-Specific Phonon MSD Signatures

Different classes of superionic conductors exhibit characteristic phonon MSD signatures that reflect their underlying conduction mechanisms. In fluorite-structured compounds, strong phonon-ion interactions triggered by geometrically frustrated configurations promote the activation of mobile species and enhance superionic conduction [39]. For copper-based superionic conductors like Cu₂Se, polymorphous structures that account for local disorder yield overdamped vibrations across the entire phonon spectrum while preserving transverse acoustic phonons, consistent with experimental measurements [40].

The temperature dependence of phonon MSD provides additional insights into the transition to superionic behavior. In many materials, the MSD of mobile ions shows anomalous enhancement at characteristic transition temperatures, reflecting the onset of liquid-like ionic mobility. This MSD enhancement often correlates with the emergence of low-energy vibrational modes and a broadening of the vibrational density of states, as observed in both fluorite and antifluorite structured superionic compounds [39].

Experimental Protocols for Phonon MSD Analysis

First-Principles Calculation Procedure

The computational protocol for phonon MSD analysis begins with density functional theory calculations to determine the optimized crystal structure and interatomic forces. For accurate results, we recommend using the following parameters:

• Software: Quantum ESPRESSO or VASP for DFT calculations
• Exchange-correlation functional: PBE for initial screening, HSE06 or PBE0 for accurate band gaps
• Pseudopotentials: Projector augmented-wave (PAW) or optimized norm-conserving Vanderbilt (ONCV) pseudopotentials
• Energy cutoff: 100 Ry for plane-wave basis sets
• k-point mesh: Γ-centered grid with minimum spacing of 0.02 Å⁻¹
• Convergence criteria: Energy threshold of 10⁻⁸ Ry for electronic minimization

Following structural optimization, phonon calculations are performed using Phonopy with the following specifications:

• Supercell size: 2×2×2 or 3×3×3 expansion of the primitive cell
• Displacement distance: 0.01-0.03 Å for finite differences
• Symmetry tolerance: 0.01 Å for detecting equivalent atoms
• Post-processing: Phonon band structure, density of states, and thermal property calculations

For superionic conductors with significant anharmonicity, we recommend employing the anharmonic special displacement method (ASDM) as implemented in complementary codes to account for temperature-dependent phonon renormalization and lifetime effects [40].

Phonon MSD Calculation and Analysis

The calculation of phonon mean squared displacement from Phonopy output involves the following steps:

• Extract phonon eigenvectors and eigenvalues from the Phonopy output files (band.yaml or mesh.yaml)
• Compute the MSD for each atom type using the formula: MSD(T) = (ℏ/2N)ΣₖΣ_j (1/ωₖⱼ)⟨nₖⱼ(T)+1/2⟩|eₖⱼ(i)|² where the summation is over wavevectors k and phonon branches j, ωₖⱼ is the phonon frequency, eₖⱼ(i) is the eigenvector component for atom i, and ⟨nₖⱼ(T)⟩ is the Bose-Einstein occupation number at temperature T
• Separate contributions from different phonon branches (acoustic vs. optical) and frequency ranges
• Calculate anisotropic MSD components to identify directional preferences for ion migration
• Correlate MSD values with ionic conductivity measurements or diffusion coefficients from molecular dynamics simulations

For high-throughput screening, we recommend focusing on the MSD of the mobile ion species (e.g., Na⁺) in the temperature range of intended application (typically 300-500 K for battery electrolytes). Structures exhibiting MSD values in the top quartile of the screening set should be prioritized for further investigation.

Research Toolkit for Phonon MSD Studies

Essential Software and Computational Tools

Table 2: Research Reagent Solutions for Phonon MSD Studies

Tool/Software	Function	Application in MSD Studies
Phonopy [6]	Phonon calculations	Harmonic force constants, phonon dispersion, and density of states
Phono3py	Anharmonic phonon calculations	Three-phonon interactions and thermal conductivity
Quantum ESPRESSO [40]	DFT calculations	Electronic structure and interatomic forces
VASP [40]	DFT calculations	Alternative platform for force calculations
ALAMODE	Anharmonic lattice dynamics	Higher-order force constants and anharmonic properties
Pheasy [12]	High-order IFC extraction	Potential energy surface reconstruction
Phononwebsite [41]	Visualization	Phonon mode animation and analysis

Machine Learning Integration

The integration of machine learning methods significantly accelerates the screening process by using phonon-derived descriptors to predict ionic transport properties. The recommended approach involves:

• Feature selection: Using phonon MSD, acoustic cutoff frequency, Na⁺ VDOS center, and low-frequency phonon coupling as primary features
• Model training: Employing neural networks or ensemble methods trained on DFT-calculated properties
• Validation: Comparing ML-predicted diffusion coefficients with explicit molecular dynamics simulations
• Active learning: Iteratively improving models by incorporating new calculations for promising candidates

This ML-integrated workflow has demonstrated robust performance in identifying sodium superionic conductors, with feature importance analysis confirming the dominant role of phonon MSD in predicting ionic conductivity [35].

Validation and Experimental Correlation

Experimental Verification Methods

The computational predictions based on phonon MSD analysis require validation through experimental techniques that probe ionic transport and lattice dynamics. Key validation methods include:

• Impedance spectroscopy: For direct measurement of ionic conductivity across relevant temperature ranges
• Quasielastic neutron scattering: To probe localized diffusive motions of mobile ions
• Inelastic neutron scattering: For experimental determination of phonon dispersion and density of states
• Nuclear magnetic resonance (NMR): To characterize local environments and hopping rates of mobile ions
• X-ray and neutron diffraction: For structural characterization and identification of anion sublattice disorder

For the copper-based superionic conductor Cu₂Se, polymorphous structures accounting for local disorder yield phonon density of states in excellent agreement with experimental measurements across the 0-35 meV energy range, validating the computational approach [40]. Similarly, the temperature-dependent band gap reduction calculated through anharmonic electron-phonon coupling (35 meV zero-point motion and 212 meV reduction from 0 to 780 K) matches well with experimental observations [40].

Case Study: Sodium Superionic Conductors

In a comprehensive study of sodium superionic conductors, screening of 3903 Na-containing structures identified specific lattice dynamics signatures associated with high ionic conductivity [35]. The successful candidates exhibited:

• Low acoustic cutoff frequencies indicating soft lattices that facilitate ion migration
• Suppressed Na⁺ vibrational density near the acoustic limit, enabling easier hopping between sites
• Enhanced low-frequency coupling between Na⁺ ions and the host lattice, promoting energy transfer
• Selective contribution from only a small subset of low-frequency acoustic and optic modes to Na⁺ migration

These features collectively create an environment where the energy barriers for ion hopping are effectively reduced through dynamic lattice fluctuations, as captured by the phonon MSD descriptor.

The integration of phonon mean squared displacement analysis into high-throughput screening frameworks represents a significant advancement in the discovery of superionic conductors. By establishing a quantitative correlation between lattice dynamics descriptors and ionic transport properties, this approach enables efficient identification of promising materials from extensive databases. The phonon MSD descriptor successfully captures the essential anharmonic effects and lattice softness that promote superionic behavior, providing insights beyond traditional static structure-based predictions.

Future developments in this field will likely focus on enhancing the accuracy of anharmonic phonon calculations through machine-learning potentials and extending the screening to more complex material systems including disordered crystals and interface structures. The integration of phonon MSD analysis with experimental characterization techniques will further validate and refine these computational predictions, accelerating the development of next-generation solid electrolytes for all-solid-state batteries and other energy storage technologies.

Solving Common Phonopy Challenges and Optimizing Calculation Performance

Imaginary frequencies, indicated by negative values in phonon dispersion calculations, represent a significant challenge in first-principles lattice dynamics using the Phonopy code. Within the context of density functional theory (DFT), these frequencies signify dynamical instabilities, suggesting that the atomic configuration is not at its minimum energy or that the calculation parameters are improperly converged. For researchers relying on phonon spectra to predict thermodynamic properties and material stability, distinguishing physical instabilities from numerical artifacts is paramount. This document details the common causes of these unphysical results and provides systematic protocols for their resolution, serving as a critical resource for ensuring the accuracy and reliability of computational materials science research.

Understanding Imaginary Frequencies

Fundamental Causes

Imaginary frequencies manifest when the calculated force constant matrix has negative eigenvalues. This can stem from two primary categories of issues:

• Physical Instabilities: The crystal structure under investigation may be mechanically or dynamically unstable at the calculated level of theory. This could indicate a genuine propensity for the structure to distort to a lower-energy phase.
• Numerical Artifacts: These are non-physical results caused by inadequacies in the computational setup. They often arise from insufficient convergence of key parameters, leading to inaccuracies in the calculated forces, which are the fundamental input for phonopy.

A Practical Case Study

A clear example from a VASP/DFPT calculation illustrates a common numerical artifact. A researcher observed that a structure showed a positive (real) frequency of 0.000928 THz with a k-point mesh of 7×4×2. However, upon increasing the k-point sampling to a denser 8×14×4 mesh, a small imaginary frequency of -0.000908 THz appeared [42].

The root cause was identified as a structure that was not re-optimized for the new, more accurate k-point set. The atomic positions optimized with a coarser k-point mesh were not the true minimum for the finer mesh. Consequently, residual forces made the structure appear slightly unstable. Re-optimizing the geometry with the improved k-point sampling eliminated the imaginary frequency, confirming it was a numerical artifact [42]. This underscores that any change in the electronic structure convergence parameters (like k-points or plane-wave cutoff) necessitates a fresh geometry optimization to maintain self-consistency.

Systematic Resolution Protocol

The following workflow provides a structured approach to diagnosing and resolving the issue of imaginary frequencies.

Diagnostic and Resolution Workflow

The diagram below outlines a step-by-step protocol for addressing imaginary frequencies.

Key Parameters and Convergence

Critical Computational Parameters

Achieving reliable phonon spectra requires careful convergence of several computational parameters. The table below summarizes the key parameters, their common issues, and recommended solutions.

Table 1: Key Parameters Affecting Phonon Calculation Stability

Parameter	Common Issue	Recommended Solution	Phonopy Tag/Setting
Geometry Optimization	Structure not at minimum for current DFT settings.	Re-optimize ionic positions (and cell) with the same KPOINTS, ENCUT, and functional used for the force calculation [42].	CELL_FILENAME (to provide optimized unit cell)
k-point Sampling	Insufficient sampling leads to inaccurate forces and energies.	Systematically increase k-point density until phonon frequencies converge.	(Set in DFT code input, e.g., VASP's KPOINTS)
Supercell Size (DIM)	Too small supercell fails to capture long-range interactions.	Increase supercell size until phonon dispersions, especially at low frequencies, are stable [43].	DIM = 2 2 2 (or larger)
Displacement Parameters	Large displacement distance can introduce anharmonicity.	Use the default displacement (0.01 Å) or a smaller value if numerical noise is low.	DISPLACEMENT_DISTANCE = 0.01
Force Constants	Inaccurate force constant calculation from poor forces.	Ensure high precision in the underlying DFT force calculations (tight EDIFF, PREC=Accurate in VASP).	FC_CALCULATOR

Advanced Configuration in Phonopy

Phonopy offers several tags that provide finer control over the calculation, which can be crucial for complex systems.

• PRIMITIVE_AXES: This tag defines the transformation from the input unit cell to the primitive cell used to build the dynamical matrix. Using PRIMITIVE_AXES = AUTO allows phonopy to automatically find a standard primitive cell, which is recommended to check for possible transformations that might simplify the calculation [22].
• CREATE_DISPLACEMENTS and PM: These tags control the generation of displaced supercells. The PM = AUTO setting is often optimal, as it generates "plus-minus" displacement pairs unless they are symmetrically equivalent, offering a good balance between computational cost and accuracy [22].

The Scientist's Toolkit

Essential Research Reagents and Computational Tools

Successful phonon calculations rely on a well-defined set of computational "reagents." The following table itemizes these essential components and their functions.

Table 2: Essential Computational Tools for Phonon Calculations

Tool/Solution	Function	Implementation Example
DFT Code	Calculates the electronic ground state and atomic forces.	VASP, Quantum ESPRESSO, ABINIT, CASTEP.
Phonopy Code	Post-processes forces from DFT to calculate phonons.	phonopy --dim="2 2 2" -c POSCAR-unitcell
Geometry Optimizer	Finds the atomic configuration with minimal forces.	VASP's IBRION=2, ISIF=3 (relax ions and cell).
Force Convergence Script	Automates convergence tests for k-points and cutoff.	Bash/Python script to run series of single-point calculations.
Symmetry Analysis Tool	Determines the space group and primitive cell.	phonopy --symmetry -c POSCAR (uses spglib)
Visualization Software	Plots phonon dispersion curves and density of states.	Phonopy's built-in plotting, VESTA, matplotlib.

Imaginary frequencies in phonon calculations are a common hurdle, but a systematic approach to diagnosing and resolving them is highly effective. The key lies in a rigorous pre-calculation workflow: ensuring the structure is fully optimized and that all parameters, especially k-point sampling and supercell size, are thoroughly converged before proceeding with the computationally expensive force calculations for phonopy. By adhering to the protocols outlined in this document, researchers can confidently distinguish numerical artifacts from genuine physical instabilities, thereby enhancing the predictive power and reliability of their first-principles materials research.

In the field of computational materials science, first-principles phonon calculations using codes like Phonopy have become indispensable for predicting a wide range of material properties, including thermodynamic behavior, phase transitions, and transport phenomena [12]. The reliability of these predictions, however, critically depends on the convergence of key computational parameters. Insufficient convergence can lead to unphysical results, such as imaginary frequencies in phonon dispersion relations, and quantitatively incorrect thermodynamic properties [44]. This application note provides detailed protocols for achieving proper convergence in supercell size, k-point mesh, and DFT parameters within the context of Phonopy-based research, ensuring high-fidelity lattice-dynamical simulations.

Core Concepts and Convergence Criteria

The Supercell Size Convergence

The finite displacement method in Phonopy requires a supercell large enough to capture all relevant interatomic interactions. The fundamental quantity being converged is the matrix of force constants, which describes how the displacement of one atom affects the forces on other atoms [45]. This matrix must be sufficiently sampled until the elements decay to approximately zero with increasing interatomic distance.

Table 1: Supercell Size Convergence Guidelines Based on System Characteristics

System Type	Minimum Suggested Dimension	Convergence Criteria	Special Considerations
Simple crystals (e.g., Diamond)	7-15 Å in each direction [45]	Force constants decay to near-zero	Small primitive cells require larger supercells
Complex unit cells (e.g., 40+ atoms)	2×2×2 may suffice [45]	Phonon frequencies stabilize	Larger primitive cells naturally sample more interactions
Anisotropic crystals	Non-uniform grid matching cell geometry	Direction-dependent convergence	Elongated cells need different sampling along each axis
Molecular crystals	Explicit convergence test essential	Properties (e.g., free energy) stabilize	Required for thermodynamic accuracy [44]

For systems with long-range interactions, such as ionic compounds, the convergence distance may extend further due to slower decay of electrostatic forces. A general guideline suggests using a supercell with at least 7 Å in each direction as a starting point, though rigorous convergence testing remains essential [45].

K-point Mesh Sampling

The k-point mesh for the electronic structure calculation must be dense enough to accurately represent the electronic wavefunctions, which directly impact the calculated forces. An insufficient k-point mesh can lead to errors in force constants and consequently inaccurate phonon frequencies.

Table 2: K-point Mesh Convergence Protocol

System Type	Initial Mesh	Convergence Test	Energy Cutoff Consideration
Semiconductors/Insulators	Γ-centered 2×3×2 [44]	Total energy (meV/atom)	Must be consistent with k-points
Metals	Denser mesh required	Force constants	Smearing parameter critical
Supercell Calculations	Often Γ-point only [44]	Phonon frequencies	Increased cutoff may be necessary

For supercell calculations, which are inherently more computationally expensive, using only the Γ-point is often sufficient, as demonstrated in studies of organic semiconductor materials [44]. However, this must be verified through explicit testing.

DFT Parameter Settings

The choice of DFT parameters significantly impacts the accuracy and numerical stability of phonon calculations. Based on VASP calculations for phonons, the following parameters have been recommended for high-accuracy results [4] [44]:

• PREC = Accurate: Ensures precise numerical settings
• ENCUT = 500 eV (or higher): Plane-wave energy cutoff
• EDIFF = 1.0e-08: Strict SCF energy convergence
• ISMEAR = 0 (Gaussian smearing) and SIGMA = 0.01: Appropriate for semiconductors/insulators
• LREAL = .FALSE.: Avoids real-space projection approximations
• IBRION = -1: Prevents ionic relaxation during force calculations
• LWAVE = .FALSE. and LCHARG = .FALSE.: Reduces I/O for displacement calculations

For systems with van der Waals interactions (e.g., molecular crystals), including an appropriate dispersion correction (such as DFT-D3 with Becke-Johnson damping) is essential for accurately capturing weak intermolecular forces [44].

Experimental Protocols

Workflow for Converged Phonon Calculations

The following diagram illustrates the complete workflow for performing converged phonon calculations using Phonopy:

Step-by-Step Protocol for Supercell Convergence Test

• Initial Structure Preparation

  • Begin with a fully optimized primitive cell structure
  • Confirm the absence of residual forces (< 1 meV/Å) and stresses

• Supercell Size Screening

  • Test progressively larger supercells (e.g., 2×2×2, 3×3×3, 4×4×4)
  • For anisotropic systems, use non-uniform expansions (e.g., 2×2×3, 3×3×2)
  • Consider employing nondiagonal supercells for more efficient sampling [45]

• Force Calculations

  • Generate displaced supercells using: phonopy -d --dim N M O
  • Perform single-point force calculations for each displacement
  • Use consistent DFT parameters across all calculations

• Convergence Metrics

  • Calculate phonon frequencies at high-symmetry points
  • Monitor key properties (e.g., acoustic modes near Γ-point, optical modes)
  • Check for the absence of imaginary frequencies (except known soft modes)
  • Compute thermodynamic properties (free energy, heat capacity) for stability

• Convergence Criteria

  • Phonon frequency changes < 0.1 THz between successive supercell sizes
  • Free energy differences < 1 meV/atom
  • Force constants decay to < 1% of maximum value at supercell boundary

Advanced Technique: Nondiagonal Supercells

Traditional "diagonal" supercells scale uniformly, leading to a rapid increase in computational cost (N³ for an N×N×N grid). Nondiagonal supercells provide a mathematical equivalent with significantly reduced computational requirements by using linear combinations of primitive lattice vectors [45]. For example, a 48×48×48 q-point grid requirement would need a supercell of 110,592 atoms with diagonal supercells but only 96 atoms with nondiagonal supercells.

The Scientist's Toolkit

Table 3: Essential Computational Tools for Phonon Research

Tool/Resource	Function	Application Context
Phonopy [4]	Supercell generation, force constants, phonon properties	Primary phonon calculation workflow
VASP [4] [44]	DFT force and energy calculations	First-principles force generation
ALAMODE/hiPhive [12]	High-order IFC extraction	Anharmonic lattice dynamics
Compressive Sensing [12]	Sparse IFC reconstruction	Efficient high-order force constants
Nondiagonal Supercells [45]	Efficient Brillouin zone sampling	Reduced computational cost for large q-grids
ShengBTE/Phono3py [12]	Thermal conductivity calculation	Anharmonic property calculation

Emerging Methods and Future Outlook

The field of computational lattice dynamics is rapidly evolving with the integration of machine learning approaches. Recent advancements include:

Machine-Learning Assisted Discovery: High-throughput screening combined with machine learning has enabled the identification of lattice dynamics signatures correlated with material properties, such as ionic conductivity in sodium superionic conductors [35]. These approaches use phonon-derived descriptors (e.g., mean squared displacements, vibrational density of states) to predict material behavior across broad structural spaces.

Advanced IFC Extraction: New codes like Pheasy utilize machine-learning algorithms to accurately reconstruct potential energy surfaces with arbitrarily high-order interatomic force constants, addressing combinatorial challenges in traditional approaches [12].

Beyond Harmonic Approximations: For strongly anharmonic systems, non-perturbative approaches like the self-consistent harmonic approximation (SCHA) and temperature-dependent effective potentials (TDEP) are becoming more accessible, enabling accurate treatment of temperature-dependent phonon renormalization [12].

The integration of these advanced methodologies with established convergence protocols will further enhance the predictive power of first-principles phonon calculations, opening new avenues for materials design and discovery.

In the realm of computational materials science, first-principles phonon calculations using codes like Phonopy and Phono3py are indispensable for predicting vibrational, thermal, and transport properties of materials [46]. However, these calculations, particularly for complex systems or those requiring high accuracy, are notoriously computationally intensive [47]. This application note provides a structured guide to managing computational costs and implementing workflow automation, framed within a broader thesis on first-principles phonon research. We summarize key performance optimization strategies, detail automated protocols, and provide visual workflows to enhance research efficiency and reproducibility.

Computational Cost Management Strategies

Effective management of computational resources is paramount. The strategies below are categorized to help researchers make informed decisions.

Table 1: Strategies for Computational Cost Management in Phonon Calculations

Strategy Category	Specific Method/Parameter	Key Consideration / Impact on Cost
Core Method Selection	Finite Displacement Method [47]	Requires multiple supercell calculations; cost scales with number of atoms and symmetry.
	Density Functional Perturbation Theory (DFPT) [48]	Can be more efficient for certain properties (e.g., electron-phonon coupling) but is mathematically complex.
Supercell & Sampling	Supercell Size (DIM tag) [22]	Larger supercells capture long-range interactions but increase cost dramatically. A 2x2x2 supercell multiplies atoms by 8.
	PRIMITIVE_AXES Tag [22]	Using the smallest possible primitive cell reduces the dimension of the dynamical matrix, lowering computational load.
	K-Points Mesh [49]	Denser k-points improve accuracy but increase cost. Convergence tests are essential.
Displacement Generation	DIAG = .FALSE. [22]	Reduces the number of displacements by restricting them to axial directions, recommended for anisotropic cells.
	PM = AUTO or .FALSE. [22]	Controls the generation of plus-minus displacement pairs. .FALSE. yields the minimum number of displacements.
	DISPLACEMENT_DISTANCE [22]	Default is 0.01 Å. A larger value (e.g., 0.02 Å) may be used with certain codes (WIEN2k, Abinit) but can affect accuracy.
Advanced & ML Methods	Machine Learning Interatomic Potentials (MLIPs) [47]	MLIPs like MACE can predict forces, drastically reducing the need for expensive DFT supercell calculations.
	Random Displacements [22]	Using RANDOM_DISPLACEMENTS with a subset of supercells can be a data-efficient way to train MLIPs or fit force constants.

Detailed Experimental Protocols

Protocol 1: Standard Harmonic Phonon Calculation with Finite Displacement

This protocol outlines the standard procedure for calculating harmonic phonon properties using the finite displacement method in Phonopy [50] [22].

• Geometry Optimization: Begin with a full relaxation of the unit cell's atomic positions and lattice parameters using DFT. Use a high k-point mesh and enforce symmetry (ISYM=2) if the structure is known to be correct [49].
• Supercell Generation: Generate a supercell using the DIM tag. A 2x2x2 or 3x3x3 supercell is common. Always check convergence with supercell size.
  This creates displacement supercells (e.g., SPOSCAR-001).
• Force Calculations: Perform single-point DFT calculations for each displaced supercell to compute the forces on all atoms. The INCAR should typically set NSW = 0 (or 1) and IBRION = -1.
• Force Sets Collection: Collect all the calculated forces into a FORCE_SETS file.

• Phonon Analysis: Calculate phonon band structure, density of states (DOS), and thermal properties.
  • Band Structure: Create a band.conf file specifying high-symmetry paths, then run:

  • DOS & Thermal Properties: Use a dense q-point mesh.

Protocol 2: Quasi-Harmonic Approximation (QHA) for Thermal Properties

The QHA approximates anharmonic effects by calculating phonons at multiple volumes to obtain volume-dependent vibrational free energy [50].

• Volume Sampling: Calculate the total energy and lattice parameter at equilibrium (E0, V0). Create several strained cells (e.g., ±2% of V0 in 0.5% increments).
• Phonon at Multiple Volumes: For each strained structure, perform a full harmonic phonon calculation following Protocol 1.
• Thermodynamic Integration: For each volume, compute the vibrational Helmholtz free energy ( F(T; V) ) using Phonopy.
• Equation of State Fitting: At each temperature, fit the ( F(T; V) ) vs. ( V ) data to an equation of state (e.g., Birch-Murnaghan). The minimum of this curve gives the equilibrium volume at that temperature, ( V(T) ).
• Property Extraction: Derive temperature-dependent properties like the thermal expansion coefficient and constant-pressure heat capacity from the fitted ( F(T, V(T)) ).

Protocol 3: Machine Learning-Accelerated Phonon Calculations

This protocol leverages MLIPs to drastically reduce the number of required DFT calculations [47].

• Training Set Generation:
  • For a set of materials, generate a small number of supercells (e.g., ~6 per material).
  • In each supercell, apply a collective random displacement to all atoms, with displacement distances ranging from 0.01 to 0.05 Å.
  • Perform DFT calculations on these randomly perturbed supercells to obtain the reference forces.
• MLIP Training: Train a universal MLIP model (e.g., MACE) on the aggregated dataset of structures and forces. The model learns the underlying potential energy surface across different chemistries and structures.
• Force Prediction: For a new material of interest, use the trained MLIP to predict the forces for the displacement patterns required by Phonopy, bypassing the need for explicit DFT force calculations.
• Phonon Property Calculation: Proceed with the standard Phonopy workflow using the MLIP-predicted forces to compute phonon dispersion and other properties.

Workflow Automation and Visualization

Automating the workflows described above is key to high-throughput materials screening. The following diagrams illustrate the logical flow for both standard and ML-accelerated protocols.

Standard Phonopy Calculation Workflow

ML-Accelerated Phonon Calculation Workflow

The Scientist's Toolkit: Research Reagent Solutions

Table 2: Essential Software and Computational "Reagents" for Phonon Research

Item Name	Function / Purpose	Key Notes
Phonopy	Primary code for harmonic and quasi-harmonic phonon calculations [6].	Calculates band structure, DOS, thermal properties; interfaces with many DFT codes.
Phono3py	Code for calculating phonon-phonon interactions and lattice thermal conductivity [46] [51].	Extends Phonopy to compute 3rd-order force constants; essential for κL.
VASP	Widely-used DFT code for computing energies and forces [49] [52].	Often used as the "force calculator" in Phonopy workflows.
DFPT	An alternative method to finite displacement for calculating force constants [49] [48].	Implemented in VASP (IBRION=7/8); can be more efficient for specific properties.
MACE MLIP	A state-of-the-art Machine Learning Interatomic Potential model [47].	Used to predict forces accurately, dramatically reducing need for DFT calculations.
EPIC STAR	A method for efficient charge transport calculation [48].	Uses DFPT results to compute carrier mobility and thermoelectric properties.

Troubleshooting Force Set Inconsistencies and Symmetry Errors

Within the framework of first-principles phonon calculations, the accuracy and reliability of computed lattice dynamics are paramount. The Phonopy code has emerged as a cornerstone tool in this domain, enabling researchers to predict fundamental materials properties such as thermal conductivity, phase stability, and spectroscopic characteristics [1]. However, a frequently encountered obstacle in these calculations is the appearance of errors related to force set inconsistencies and symmetry detection failures. These errors, often manifesting as "Input forces are not enough to calculate force constants" or issues in the compute_permutation function, can jeopardize the entire computational workflow [53] [54]. This document provides a systematic protocol for diagnosing and resolving these challenges, ensuring robust phonon calculations for materials science and drug development research.

Theoretical Background: Force Constants and Symmetry

The computation of phonons relies on the construction of the force constant matrix, which represents the second-order derivative of the system's total energy with respect to atomic displacements. For a supercell containing N primitive cells, each with M atoms, there are 3MN potential degrees of freedom. Fortunately, crystal symmetry drastically reduces the number of unique force constants that need to be calculated explicitly [53].

Phonopy leverages the space group symmetry of the crystal to minimize the number of necessary displacement calculations. It generates a minimal set of displaced supercells, computes the atomic forces for each displacement via a first-principles code (e.g., VASP), and then uses this force information to build the complete force constant matrix by applying symmetry operations [4]. The error "Input forces are not enough to calculate force constants" arises when the symmetry analysis fails or becomes inconsistent with the actual force data, preventing the successful reconstruction of this matrix [53] [54].

Diagnostic Protocol: Identifying the Root Cause

A systematic approach to diagnosing these errors is crucial. The following workflow outlines a step-by-step diagnostic procedure, with corresponding logic visualized in the diagram below.

Symptom Analysis

The initial error trace typically occurs during the post-processing phase, when Phonopy attempts to compute force constants from the FORCE_SETS file. The key error message is: "ValueError: Input forces are not enough to calculate force constants, or something wrong (e.g. crystal structure does not match)" [53] [54].

Diagnostic Steps

• Crystal Structure Verification: Confirm that the initial POSCAR (unit cell) and the generated SPOSCAR (supercell) are consistent. The atomic coordinates in all POSCAR-{number} files must be derivatives of the same base structure.
• Force Sets Completeness Check: Ensure the FORCE_SETS file contains force information for every displaced supercell generated in the pre-process. The number of force sets must match the number of POSCAR-{number} files.
• Symmetry Tolerance Assessment: The default symmetry tolerance (symprec) used by Phonopy might be too strict for your system, leading to an incorrect identification of symmetry operations. This is a common root cause [53] [54].

Resolution Strategies and Experimental Protocols

Primary Resolution: Symmetry Tolerance Adjustment

The most common and effective solution is to adjust the symmetry tolerance parameter. This tells the symmetry finding library (Spglib) within Phonopy to be more lenient when identifying symmetry operations.

Protocol:

• Execute Phonopy with an explicitly defined --tolerance flag.
• Start with a value one order of magnitude larger than the default (e.g., 1e-4 or 1e-5).

• If the error persists, gradually increase the tolerance (e.g., to 2e-2 as used in a reported case [54]) until the calculation succeeds. Be aware that an excessively high tolerance might incorrectly assign symmetry.

Comprehensive Computational Workflow

A robust end-to-end workflow for Phonopy calculations, integrating troubleshooting steps, is essential for reliable results. The following diagram and protocol detail this process.

Pre-Processing and Force Calculation:

• Generate Displacements: Use phonopy -d --dim "2 2 1" to create the supercell (SPOSCAR) and a set of displaced supercells (POSCAR-001, POSCAR-002, ...) [4].
• First-Principles Calculations: Run force calculations for each POSCAR-{number} file using your chosen DFT code (e.g., VASP). The critical requirement is that the atomic positions must not be relaxed (IBRION = -1 in VASP's INCAR) [4].
• Compile Force Sets: Use the Phonopy interface to collect all the force data into a single FORCE_SETS file. For VASP, this is done with: phonopy -f disp-{001..003}/vasprun.xml [4].

Post-Processing with Error Handling:

• Initial Analysis: Run the post-process command (e.g., phonopy -p band.conf).
• Error Diagnosis and Resolution: If the symmetry error occurs, implement the tolerance adjustment protocol from Section 4.1.
• Result Extraction: Once successful, proceed to calculate phonon band structures, density of states (DOS), and thermal properties.

Advanced Consideration: Non-Analytical Term Correction

For accurate phonons in polar materials, especially near the Brillouin zone center, including the non-analytical term correction (NAC) is essential.

Protocol:

• Perform a DFPT linear response calculation with your DFT code to obtain the Born effective charges and the dielectric tensor (e.g., in VASP, set LEPSILON = .TRUE.) [4].
• Use the Phonopy auxiliary tool to generate the BORN file. For VASP, the command is phonopy-vasp-born. The --st option can apply further symmetry constraints to the values [4].
• Place the BORN file in your working directory. Phonopy will automatically use it to apply the NAC during post-processing.

The Scientist's Toolkit: Essential Research Reagents

Table 1: Key Computational "Reagents" and Their Functions in Phonopy Calculations.

Research Reagent	Function	Critical Notes
POSCAR (Unit Cell)	The initial, optimized crystal structure input.	Must be a fully relaxed, stable structure. The foundation of all subsequent steps.
SPOSCAR (Supercell)	The perfect supercell structure generated by Phonopy.	Used as a reference for applying displacements and symmetry operations.
POSCAR-{number}	Supercell structures with finite atomic displacements.	Each file is the input for an independent force calculation.
FORCE_SETS	Data file collating forces from all displacement calculations.	The number of sets must match the number of displaced supercells. Inconsistency causes errors [53].
phonopy_disp.yaml	Metadata file recording displacement magnitudes and directions.	Useful for verifying the completeness of the force calculations.
BORN File	Contains Born effective charges and dielectric tensor for NAC.	Crucial for correcting the long-range Coulomb interaction in polar materials [4].
Symmetry Tolerance (--tolerance)	Numerical threshold for identifying crystal symmetry.	A primary parameter for troubleshooting symmetry-related errors. Adjust if the default fails [53] [54].

Data Presentation and Analysis

Table 2: Summary of Common Phonopy Errors and Corresponding Solutions.

Error Symptom	Potential Root Cause	Recommended Solution	Reference
"Input forces are not enough to calculate force constants"	Symmetry detection failure due to tight tolerance.	Rerun phonopy with increased --tolerance (e.g., 1e-4).	[53]
"compute_permutation" failure	Atomic mapping error after symmetry rotation.	Check structure consistency; increase --tolerance.	[54]
Incomplete phonon spectrum	Missing BORN file for polar material.	Calculate Born charges and dielectric constant to create BORN file.	[4]
Force constants calculation crash	Mismatch between expected and provided force sets.	Verify all POSCAR-{number} files have corresponding force data in FORCE_SETS.	[53] [4]

Benchmarking Phonopy and Its Role in the Modern Computational Ecosystem

The accuracy of phonon calculations is paramount in computational materials science, influencing the prediction of thermodynamic properties, phase stability, and electron-phonon interactions. For researchers relying on the Phonopy code, which implements the finite displacement method, validation of results against established experimental data and alternative computational methods like Density Functional Perturbation Theory (DFPT) represents a critical step in ensuring research reliability. This application note details comprehensive protocols for validating phonon calculations, providing structured comparisons and practical workflows tailored for the Phonopy environment. By framing these procedures within a broader thesis on first-principles phonon calculations, we aim to equip researchers with robust methodologies for verifying computational predictions against experimental measurements and other first-principles approaches.

Theoretical Background and Computational Methods

Finite Displacement Method in Phonopy

Phonopy determines second-order force constants by performing finite displacements of atoms in a supercell and calculating the resulting forces [55]. The dynamical matrix is then constructed from these force constants and diagonalized to obtain phonon frequencies and eigenvectors. The key parameters controlling this process include:

• Supercell size (DIM): Determines the size of the supercell used for displacements, affecting q-point sampling [22]
• Displacement distance (DISPLACEMENT_DISTANCE): Default is 0.01 Å, but can be adjusted [22]
• Displacement pattern (PM): Controls whether only least displacements or plus-minus displacements are generated [22]

Density Functional Perturbation Theory (DFPT)

In contrast to the finite displacement approach, DFPT calculates phonon properties by computing the linear response of the electron system to atomic displacements, typically providing phonon spectra at any q-point without requiring large supercells. While DFPT is implemented in codes like VASP (IBRION=7 or 8), Phonopy utilizes the finite displacement approach (IBRION=5 or 6 in VASP) [55], making methodological comparisons essential for validation.

Validation Protocols

Comparison with Experimental Data

Table 1: Primary Experimental Validation Methods for Phonon Calculations

Experimental Technique	Measurable Phonon Properties	Comparison Methodology	Key Considerations
Inelastic Neutron Scattering	Phonon dispersion curves throughout Brillouin Zone	Direct comparison of calculated and measured frequencies along high-symmetry paths	Sample quality, temperature effects, instrumental resolution
Raman Spectroscopy	Zone-center optical phonon frequencies	Compare measured and calculated Γ-point frequencies	Selection rules, resonant effects, polarization dependence
Infrared Spectroscopy	Zone-center optical phonons with LO-TO splitting	Compare frequencies when NAC included	Dielectric function accuracy, TO-LO splitting magnitude
Specific Heat Measurements	Phonon density of states	Compare calculated and measured Cv(T)	Anharmonic effects at elevated temperatures, experimental uncertainty

Protocol: Validating Against Experimental Dispersion Curves

• Calculation Setup: Perform Phonopy calculation with sufficiently large supercell (e.g., 2×2×2 or larger) to enable accurate Fourier interpolation of force constants [4].

• Band Structure Generation: Calculate phonon frequencies along the same high-symmetry paths used in experimental studies using the --band option [4].

• Data Extraction: From experimental papers, extract numerical data for phonon frequencies at specific q-points using digitization tools if necessary.

• Quantitative Comparison: Calculate root-mean-square deviation (RMSD) between calculated and experimental frequencies: RMSD = √[Σ(f_calc - f_exp)²/N]

• Systematic Error Analysis: Identify any systematic shifts (e.g., uniform overestimation) that might indicate exchange-correlation functional limitations.

Comparison with DFPT Methods

Table 2: Finite Displacement vs. DFPT Comparison Framework

Comparison Aspect	Finite Displacement (Phonopy)	DFPT	Validation Protocol
Computational Cost	Scales with supercell size (3N calculations)	Scales with number of q-points	Compare computational time for similar system sizes
Γ-Point Accuracy	Direct calculation in large supercell	Direct calculation	Compare frequencies at Γ-point with identical computational parameters
Brillouin Zone Sampling	Dependent on supercell size	Independent q-point sampling	Compare dispersion along high-symmetry directions
LO-TO Splitting	Requires separate NAC calculation with BORN file [4]	Can be included directly	Compare LO and TO frequencies in polar materials
Code Implementation	IBRION=5/6 in VASP [55] with Phonopy post-processing [4]	IBRION=7/8 in VASP	Use identical POTCAR, INCAR parameters except IBRION

Protocol: Direct Phonopy-DFPT Benchmarking

• Consistent Computational Parameters: Use identical pseudopotentials, plane-wave cutoff (ENCUT), k-point mesh, and exchange-correlation functional for both calculations.

• DFPT Reference Calculation: Perform DFPT phonon calculation (IBRION=8 in VASP) for the primitive cell, calculating the full dynamical matrix.

• Phonopy Calculation: Perform finite displacement calculation with Phonopy using a supercell that provides equivalent q-point sampling.

• Frequency Comparison: Extract and compare frequencies at high-symmetry points, particularly focusing on:

  • Γ-point optical modes
  • Acoustic modes near Γ-point
  • Anomalies or soft modes, if present

• Convergence Testing: Verify that both methods converge to the same results with increasing computational parameters (k-points, supercell size).

Workflow for Comprehensive Validation

The following diagram illustrates the integrated validation workflow combining both experimental and computational benchmarking:

Practical Implementation Guide

Phonopy Calculation with VASP: Complete Workflow

The diagram below details the specific steps for generating phonon properties with Phonopy and VASP, highlighting critical validation checkpoints:

Step-by-Step Protocol

• Pre-processing:

  This generates displacement supercells (POSCAR-001, POSCAR-002, ...) and a perfect supercell (SPOSCAR) [4].

• VASP Force Calculations: Use consistent INCAR parameters with PREC = Accurate, IBRION = -1 (no relaxation), and appropriate ENCUT [4]. For NAC, include LEPSILON = .TRUE. in a separate calculation to generate BORN file using phonopy-vasp-born [4].

• Force Sets Generation:

  Creates FORCE_SETS file from VASP outputs [4].

• Post-processing and Validation:

  • Band structure: phonopy --band "q-points" -p
  • DOS: phonopy --mesh "20 20 20" --dos -p
  • Thermal properties: phonopy --mesh "20 20 20" -t -p [4]

Troubleshooting Common Discrepancies

Systematic Overestimation/Underestimation of Frequencies

• Cause: Exchange-correlation functional limitations
• Solution: Compare multiple functionals (PBE, PBEsol, SCAN) or apply scaling factors based on known benchmarks

Imaginary Frequencies at Γ-point

• Cause: Structural instability or insufficient structural relaxation
• Solution: Re-relax structure with stricter criteria; verify symmetry; check for actual structural instability

Poor Agreement with DFPT at Specific q-points

• Cause: Insufficient supercell size leading to q-point sampling issues
• Solution: Increase supercell dimensions and repeat calculation

Inconsistent LO-TO Splitting

• Cause: Incorrect or missing BORN file
• Solution: Recalculate dielectric tensor and Born effective charges with high accuracy parameters [4]

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools for Phonon Validation

Tool/Resource	Type	Function in Validation	Implementation Notes
Phonopy Code	Software Package	Finite displacement phonon calculations	Use latest version; configure with proper symmetry settings [4]
VASP	DFT Code	Force calculator for Phonopy	Employ IBRION=5/6 for finite differences or IBRION=7/8 for DFPT [55]
BORN File	Input File	Enables non-analytical term correction for LO-TO splitting	Generate using phonopy-vasp-born with LEPSILON=.TRUE. calculation [4]
Experimental Data Repository	Data Resource	Reference data for validation	ICSD for structures; experimental phonon databases for dispersion
Phonopy-Dfpt Toolkit	Validation Scripts	Automated comparison between methods	Custom Python scripts for RMSD calculation and visualization
High-Performance Computing	Infrastructure	Enables large supercell calculations	Scale to appropriate core counts for force calculations

Robust validation of Phonopy results against both experimental data and DFPT calculations establishes essential confidence in computational predictions of lattice dynamical properties. The protocols detailed in this application note provide a systematic framework for such validation, emphasizing quantitative comparison metrics and troubleshooting strategies. By implementing these procedures within a broader research workflow, scientists can ensure the reliability of their phonon calculations while developing intuition for the limitations and strengths of the finite displacement method as implemented in Phonopy. As machine learning approaches continue to emerge in lattice dynamics [35] [56], these validation methodologies will become increasingly important for verifying new computational paradigms.

The calculation of phonons—the quantized lattice vibrations in materials—is a cornerstone of computational materials science, essential for understanding thermodynamic stability, thermal conductivity, and spectroscopic properties. For years, density functional theory (DFT) has served as the primary computational tool for these calculations, with the Phonopy code emerging as a standard for calculating harmonic phonon properties via the finite displacement method. However, the computational expense of DFT, which scales cubically with system size, has remained a significant bottleneck, particularly for large supercells or high-throughput screening.

The emergence of universal machine learning interatomic potentials (uMLIPs) represents a paradigm shift. These potentials, trained on vast datasets of DFT calculations, can accurately predict energies and atomic forces at a fraction of the computational cost of full DFT calculations. This advancement creates a powerful synergy with Phonopy, enabling rapid, high-throughput phonon calculations while maintaining near-DFT accuracy. This integration is particularly crucial for properties dependent on higher-order derivatives of the potential energy surface, such as phonon dispersion and thermal conductivity, where accuracy in force predictions is paramount.

Universal MLIPs: A New Tool for Force Evaluation

Core Principles of uMLIPs

Universal MLIPs are deep learning models designed to approximate the potential energy surface (PES) of diverse atomic systems. Unlike traditional empirical potentials or system-specific MLIPs, uMLIPs are trained on extensive datasets encompassing a wide range of elements and crystal structures, granting them remarkable transferability. The core task of any MLIP is to map atomic configurations (positions, species) to the total potential energy and, crucially, the atomic forces, which are the negative gradients of the energy [57] [58].

The key advancement enabling uMLIPs is the integration of geometric equivariance into model architectures. Specifically, E(3)-equivariance ensures that model predictions transform correctly under Euclidean symmetry operations (rotations, translations, and reflections), which is essential for generating physically consistent forces and vibrational modes [57]. Architectures such as equivariant graph neural networks (GNNs) and transformers have proven particularly successful in this domain.

Benchmark Performance for Phonon Properties

Recent comprehensive benchmarks have evaluated the performance of various uMLIPs specifically for phonon calculations. One study tested seven uMLIPs on a dataset of ~10,000 non-magnetic semiconductors, comparing their predictions of harmonic phonon properties against direct DFT results [59]. The findings reveal that while some uMLIPs achieve high accuracy, others show substantial inaccuracies, even if they perform well on energy and force predictions for equilibrium structures.

Table 1: Performance of Selected uMLIPs in Geometry Relaxation and Phonon Calculations

Model Name	Key Architectural Features	Force Conservation	Typical Energy MAE (eV/atom)	Performance on Phonons
M3GNet [59]	Three-body interactions, GNN	Conservative (Energy-derived)	~0.035 (on equilib. structures)	Foundational model, reasonable performance
CHGNet [59]	Incorporates magnetic moments	Conservative	Higher than others (no correction applied)	Reliable geometry convergence
MACE-MP-0 [59] [60]	Atomic cluster expansion	Conservative	Low	Good accuracy
MatterSim-v1 [60]	Active learning on diverse data	Conservative	Low	Identified as top-performer for defect PL spectra
eqV2-M [59]	Equivariant transformers	Non-conservative (Direct force prediction)	Low	High force errors can hinder relaxation
ORB [59]	Graph network simulator	Non-conservative	Low	High failure rate in geometry convergence

The benchmark highlights a critical point: excellent performance on energy and force prediction for structures near equilibrium does not automatically guarantee accuracy for phonon properties, which probe the curvature of the PES [59]. Furthermore, conservative potentials—those where forces are derived as the explicit gradient of the energy—generally demonstrate superior performance and reliability in geometry optimization and phonon calculations compared to non-conservative models that predict forces directly [61]. The failure of non-conservative models like ORB and eqV2-M in some geometry relaxations has been attributed to high-frequency errors in their predicted forces [59].

Integrated Protocol: Phonopy with uMLIPs

The following protocol details the steps for performing phonon calculations using Phonopy driven by a universal MLIP, replacing the traditional DFT force calculator.

The diagram below illustrates the integrated workflow, highlighting the replacement of the DFT force calculator with a uMLIP.

Step-by-Step Protocol

Step 1: Pre-processing and Generation of Displaced Supercells

This step generates the structurally perturbed supercells for which forces must be calculated.

• Prepare Input Files: You will need a POSCAR file containing the fully optimized structure of your system's primitive or conventional unit cell.
• Run Phonopy Displacement Generator:
  • --dim "2 2 2" defines the size of the supercell (e.g., 2x2x2 in this case). The supercell must be large enough to capture the relevant interatomic interactions.
  • --pa auto automatically determines the primitive cell, which is more efficient for phonon calculations.
• Output: This command generates:
  • SPOSCAR: The perfect supercell without displacements.
  • POSCAR-001, POSCAR-002, ...: A set of supercell files, each containing a single atomic displacement.
  • phonopy_disp.yaml: A file cataloging all the displacements [4].

Step 2: Force Calculation Using a Universal MLIP

This is the core step where the uMLIP replaces DFT. The following is a generic Python pseudocode illustrating how to compute forces for the displaced structures.

Key Considerations for Force Calculation:

• Model Selection: Refer to Table 1. For phonon calculations, prefer a conservative model (e.g., MatterSim, MACE-MP-0) as they generally provide more stable and reliable results for evaluating the curvature of the PES [59] [61].
• Consistency: Ensure the uMLIP was trained on data compatible with your system (e.g., chemical elements, level of theory). Most current uMLIPs are trained on PBE-DFT data [59] [62].

Step 3: Post-processing with Phonopy

With the FORCE_SETS file created, the remaining workflow is standard for Phonopy.

• Compute Force Constants and Phonons:

  This command reads the FORCE_SETS file and calculates the force constants.

• Calculate Phonon Properties:

  • Phonon Band Structure:

  • Phonon Density of States (DOS):

  • Thermal Properties:
    The -p flag plots the results, and -s saves the plot to a PDF file [4].

Application Notes and Case Studies

High-Throughput Screening of Superionic Conductors

A 2025 study successfully leveraged this integrated approach to discover lattice dynamics signatures of sodium superionic conductors (NASICONs). The workflow is a prime example of the power of uMLIPs and Phonopy [35].

• Objective: Identify lattice dynamics features correlated with high Na+ ionic conductivity from 3,903 Na-containing structures.
• uMLIP Selection: After evaluating several models, the EquiformerV2 model fine-tuned on the OMAT and MPtraj datasets was selected for its accuracy and reliability.
• Screening Workflow:
  • The candidate structures were first filtered for dynamic stability. This involved using the uMLIP with Phonopy to compute phonon band structures and discard structures with imaginary frequencies.
  • For the dynamically stable structures, the uMLIP was used to compute key phonon-derived descriptors, including phonon mean squared displacement (MSD) of Na+ ions and acoustic cutoff frequencies.
• Outcome: The high-throughput analysis revealed a strong positive correlation between phonon MSD and diffusion coefficients. It identified that low acoustic cutoff frequencies and enhanced low-frequency vibrational coupling promote superionic conductivity. This discovery was accelerated by orders of magnitude using the uMLIP-Phonopy workflow compared to standard DFT [35].

Accelerating Photoemission Spectra of Point Defects

Another frontier application is the calculation of photoluminescence (PL) spectra for point defects (e.g., quantum color centers). The phonon bottleneck is a major hurdle in these calculations, as they require harmonic and anharmonic phonon computations in large, low-symmetry defect supercells [60].

• Objective: Calculate the Huang-Rhys factor and PL lineshapes for 791 point defects in 2D materials.
• Method: The standard workflow requires calculating the excited-state relaxation and the phonon modes of the ground-state supercell. The researchers replaced the DFT phonon calculation with a uMLIP-based one.
• Benchmarking: Seven uMLIPs were benchmarked against ab initio results. The study found that MatterSim-v1 delivered the best performance, accurately predicting HR factors and PL lineshapes with minimal precision loss.
• Impact: This ML-accelerated framework achieved speed improvements exceeding an order of magnitude, making high-throughput screening of defect optical properties computationally tractable. Notably, the uMLIPs performed well regardless of the defect's charge or magnetic state, despite having no explicit knowledge of these properties [60].

The Scientist's Toolkit

Table 2: Essential Software and Computational Resources

Tool / Resource	Type	Primary Function	Reference/URL
Phonopy	Software Code	Calculates harmonic phonon properties from force sets.	https://phonopy.github.io/ [63] [4]
VASP	DFT Calculator	(Reference) Generates high-fidelity training data and benchmark results.	N/A
MatterSim-v1	uMLIP	A top-performing conservative potential for phonon and defect property prediction.	[59] [60]
MACE-MP-0	uMLIP	A high-performance model using atomic cluster expansion.	[59] [60]
CHGNet	uMLIP	A model incorporating magnetic moments, good for geometry convergence.	[59]
Matbench Discovery	Benchmark Platform	Leaderboard for comparing uMLIP performance on materials science tasks.	[59] [62]
Materials Project	Database	Source of crystal structures and data for training and testing.	[59]

In the field of first-principles materials science, computational modeling of lattice dynamics is indispensable for predicting thermal properties. The calculation of lattice thermal conductivity (κL) is a central challenge, with implications for thermoelectrics, thermal management, and electronics. This analysis examines the distinct roles and capabilities of three computational codes: Phonopy, a cornerstone for harmonic lattice dynamics, and the anharmonic codes Phono3py and Pheasy, which are specialized for calculating phonon-phonon interactions and thermal transport [12] [6]. Phonopy operates within the harmonic or quasi-harmonic approximation, providing foundational phonon properties like dispersion relations and density of states. In contrast, Phono3py and Pheasy address the critical limitation of the harmonic approximation by explicitly treating anharmonicity—the phonon-phonon scattering that governs finite phonon lifetimes and thus, thermal conductivity [12] [64]. This document provides a comparative analysis and detailed application protocols to guide researchers in selecting and employing the appropriate tool for their investigations into phonon-mediated thermal transport.

Code Comparison: Scope, Theory, and Application

• Phonopy: This open-source package is the workhorse for calculating harmonic and quasi-harmonic phonon properties. It employs the finite displacement method to compute second-order interatomic force constants (IFCs), which define the force-free phonon band structure and density of states [1] [6]. Its predictions are reliable for thermodynamic properties like the free energy and heat capacity within the harmonic approximation. However, as it does not compute third-order anharmonic IFCs, it cannot directly predict phonon lifetimes or intrinsic lattice thermal conductivity.

• Phono3py: An extension of the Phonopy ecosystem, Phono3py is specifically designed to compute phonon-phonon interactions and related properties, most notably κL [63] [65]. Its core capability is the efficient calculation of second-order (Φ) and third-order (Ψ) IFCs from supercell force calculations [64]. Using these IFCs, it solves the linearized phonon Boltzmann transport equation (LBTE) to obtain κL, typically within the single-mode relaxation time approximation (RTA) [66] [63]. It has been widely validated for numerous solids but faces a combinatorial explosion in computational cost when attempting to include fourth-order or higher anharmonicity [12].

• Pheasy: Positioned as a next-generation tool, Pheasy uses advanced machine-learning algorithms (e.g., compressive-sensing lattice dynamics) to accurately reconstruct the potential energy surface and extract IFCs to arbitrarily high orders from force-displacement datasets [12]. This approach mitigates the computational bottleneck of high-order IFCs, enabling robust non-perturbative treatments of strong anharmonicity via techniques like the self-consistent harmonic approximation (SCHA) [12]. It also aims to create a "phonon code ecosystem" by interfacing with other codes like ShengBTE and Phono3py for specific property calculations [12].

Comparative Analysis: Capabilities and Limitations

Table 1: Comparative analysis of phonon computational codes.

Feature	Phonopy	Phono3py	Pheasy
Primary Scope	Harmonic & Quasi-harmonic properties [6]	3rd-order anharmonicity & thermal conductivity [63]	Arbitrarily high-order anharmonicity [12]
Key Outputs	Phonon dispersion, DOS, Thermodynamics [1]	Phonon lifetimes, ΓL (RTA & beyond) [63] [67]	Anharmonic phonon spectra, ΓL, Thermal transport [12]
Anharmonic Treatment	Not applicable	Perturbative (3-phonon scattering) [12]	Non-perturbative (SCHA, SCP) [12]
IFC Extraction	Finite displacement method [1]	Finite displacement method [65]	Machine-learning (Compressive sensing) [12]
High-Order IFCs	Not supported	Up to 3rd-order [12]	Arbitrary order (Sparse representation) [12]
Code Interoperability	Phonopy ecosystem [6]	Interfaces with VASP, QE [65]	Interfaces with ShengBTE, Phono3py, EPW [12]

Protocols for Lattice Thermal Conductivity Calculation

Protocol 1: Standard Workflow with Phono3py and VASP

This protocol details the steps for a standard κL calculation, which is computationally intensive and requires careful convergence testing [66] [65].

Table 2: Key research reagents and computational tools for Phono3py calculations.

Reagent / Tool	Function / Description
POSCAR-unitcell	Input file defining the crystal structure of the primitive unit cell.
VASP	First-principles code for calculating electronic structure and atomic forces [66] [65].
Supercell Model	A larger cell built from the unit cell, used to capture the relevant interatomic interactions [65].
Displacement Dataset	Sets of atomic displacements in supercells, generated by Phono3py, for which forces are computed [63].
Force Sets	Collected results of force calculations for all displaced supercell configurations [63].
Force Constants (fc2, fc3)	Harmonic (fc2) and anharmonic (fc3) coefficients of the potential energy surface, extracted from force sets [63] [65].

Step-by-Step Procedure:

• Preparation of Input Files:

  • Begin with a fully relaxed primitive unit cell (POSCAR-unitcell).
  • Prepare a high-precision INCAR file for VASP force calculations. Critical settings include:
    • PREC = Accurate
    • IBRION = -1 (ions are not moved)
    • ISMEAR = 0; SIGMA = 0.01
    • LREAL = .FALSE.
    • ADDGRID = .TRUE. (crucial for accurate forces) [66].
  • Generate a WAVECAR file from a prior electronic convergence calculation to accelerate force computations [66].

• Generate Displacement Supercells:

  • Run Phono3py to create the displacement dataset. For a 2x2x2 supercell:

  • This generates many POSCAR-XXXXX files, each representing a supercell with a unique set of atomic displacements [65].

• Calculate Forces on Displaced Supercells:

  • For each POSCAR-XXXXX, run a VASP calculation to compute the forces on all atoms. This step is computationally intensive and can be automated using scripts to manage job submission on clusters [66] [65].
  • The output for each calculation is a vasprun.xml file.

• Collect Forces and Calculate Force Constants:

  • Collect all the vasprun.xml files to build the force sets:

  • Symmetrize the second- and third-order force constants and save them to fc2.hdf5 and fc3.hdf5 [65]:

• Compute Lattice Thermal Conductivity:

  • Perform the κL calculation by solving the LBTE. A dense q-point mesh is required for convergence [66] [63].

  • The --mesh parameter must be tested for convergence [66]. The calculation can be distributed across many computer nodes by calculating phonon lifetimes at irreducible q-points separately using the --wgp and --write-gamma options, then combined with --read-gamma [65].

• Visualization and Analysis:

  • Results are stored in HDF5 format (e.g., kappa-m111111.hdf5). Use tools like hdfview or the phono3py-kdeplot script to visualize results, such as the distribution of phonon lifetimes versus frequency [67].

Below is the workflow for calculating lattice thermal conductivity using Phono3py.

Protocol 2: High-Order Anharmonicity Workflow with Pheasy

Pheasy offers a different approach, leveraging machine learning to access high-order anharmonic effects that are prohibitively expensive for conventional methods [12].

Step-by-Step Procedure:

• Data Generation:

  • Perform first-principles calculations (using VASP, Quantum ESPRESSO, etc.) on a set of supercell configurations to generate a force-displacement dataset. This dataset does not necessarily require the systematic multiple displacements of the finite displacement method; it can include configurations from, for example, molecular dynamics snapshots.

• Machine-Learning Assisted IFC Extraction:

  • Use Pheasy's algorithms (e.g., based on compressive sensing) to fit a high-order (beyond 3rd) Taylor expansion of the potential energy surface to the force-displacement dataset [12]. This step efficiently identifies a sparse representation of the physically relevant IFCs.

• Anharmonic Lattice Dynamics:

  • Use the extracted high-order IFCs to perform non-perturbative calculations. Pheasy implements methods like the Self-Consistent Harmonic Approximation (SCHA) to compute anharmonically renormalized phonon frequencies and lifetimes, which include effects of multi-phonon scattering [12].

• Property Calculation and Interfacing:

  • Calculate thermal conductivity within Pheasy or by writing the obtained IFCs in a format compatible with other thermal transport codes like ShengBTE or Phono3py, leveraging Pheasy's ecosystem approach [12].

The logical structure of a Pheasy computation is outlined below.

The selection of a computational code for lattice dynamics should be guided by the specific scientific question. Phonopy remains an essential, efficient tool for harmonic property analysis and preliminary structure stability screening. For direct calculations of lattice thermal conductivity in materials where third-order anharmonicity dominates, Phono3py is the established, robust choice with well-documented protocols and extensive community use. For systems exhibiting strong anharmonicity—such as those with phase transitions, rattlers, or soft modes—where fourth and higher-order phonon scattering processes are significant, Pheasy represents the cutting edge. Its machine-learning-driven approach to high-order IFC extraction provides a path to more accurate and predictive simulations of thermal transport in complex materials, marking a significant advancement beyond conventional perturbative approaches.

Phonopy is an open-source package that has become a cornerstone in computational materials science for performing phonon calculations at the harmonic and quasi-harmonic levels. As a robust and versatile tool, it enables researchers to investigate lattice dynamics properties from first principles, providing critical insights into material behavior, stability, and thermal properties. The package serves as a central hub that interfaces with numerous electronic structure codes, creating an extensive ecosystem for property prediction. This interoperability allows scientists to leverage the strengths of different computational approaches while maintaining a standardized workflow for phonon analysis. Within the broader context of first-principles phonon calculations, Phonopy bridges the gap between electronic structure calculations and thermodynamic properties, enabling the prediction of experimentally relevant quantities such as thermal conductivity, phase diagrams, and vibrational spectra.

The Phonopy ecosystem extends beyond basic harmonic approximations through its companion code, Phono3py, which specializes in phonon-phonon interaction and lattice thermal conductivity calculations. This integrated approach provides researchers with a comprehensive toolkit for investigating lattice dynamics across different temperature regimes and material systems. The continued development of these tools, supported by the National Institute for Materials Science, ensures they remain at the forefront of computational materials research, adapting to new challenges and incorporating advances in both methodology and high-performance computing architectures [6].

Core Components of the Phonopy Ecosystem

The Phonopy architecture is built around a modular design that separates the core phonon calculation routines from the interfaces to external force calculators. This design principle enables its remarkable flexibility in working with diverse computational chemistry and materials physics packages. At its core, Phonopy implements the finite displacement method for constructing dynamical matrices through systematic atomic displacements in supercells. The key components include displacement pattern generators, symmetry analysis tools, force constants calculators, and post-processors for spectral and thermodynamic properties.

The software's notable features include the calculation of phonon band structures, density of states (DOS), projected DOS (PDOS), thermodynamic properties, and mean-square displacements. Phonopy also incorporates support for non-analytical term correction (NAC) for accurately modeling longitudinal-transverse optical (LO-TO) splitting in polar materials, which requires Born effective charges and dielectric tensors typically obtained from density functional perturbation theory (DFPT) calculations [6] [22]. The recent introduction of PRIMITIVE_AXES = AUTO functionality exemplifies the ongoing development to automate and simplify workflows while maintaining physical accuracy [22].

Phono3py: Extending to Three-Phonon Interactions

Phono3py serves as a specialized extension to the Phonopy ecosystem, focusing specifically on three-phonon interactions and their role in thermal transport properties. The calculation of lattice thermal conductivity using Phono3py follows a well-defined workflow consisting of three primary steps: (1) calculation of supercell-force sets, (2) computation of second- and third-order force constants, and (3) calculation of lattice thermal conductivity itself. This structured approach enables the efficient handling of computationally demanding anharmonic lattice dynamics problems [63].

When the input unit cell is not a primitive cell, Phono3py requires the primitive cell matrix to be specified for proper calculation setup. Like Phonopy, it also supports the inclusion of long-range dipole-dipole interactions through non-analytical term correction when the appropriate parameters are provided. This capability is particularly important for accurately modeling thermal transport in polar semiconductors and insulators, where LO-TO splitting significantly influences phonon scattering rates [63].

Interfacing Phonopy with Electronic Structure Codes

Standardized Workflows for Property Prediction

The interoperability of Phonopy with various computational chemistry packages follows standardized workflows that maintain consistency across different force calculators. Table 1 summarizes the primary interface capabilities with selected electronic structure codes, highlighting the specific requirements and adaptations for each.

Table 1: Phonopy Interfaces with Selected Electronic Structure Codes

Code	Displacement Distance (Default)	Force Output Files	Special Requirements
VASP	0.01 Å	vasprun.xml or OUTCAR	INCAR must prevent ionic relaxation (IBRION = -1)
WIEN2k	0.02 Bohr	case.scf	Lattice vectors in case.struct must be row vectors
Abinit	0.02 Bohr	netcdf or text output	--
Fleur	0.02 Bohr	force.txt	--
Quantum ESPRESSO	0.01 Å	pwscf output	--

A generalized workflow for phonon property prediction consists of three main stages: pre-processing, force calculation, and post-processing. In the pre-processing stage, Phonopy generates displacement supercells based on user-specified dimensions and displacement parameters. The force calculation stage involves running separate electronic structure calculations for each displaced supercell using the preferred computational package. Finally, in the post-processing stage, Phonopy reads the calculated forces, constructs the force constants, and computes the desired phonon properties [4].

The configuration of these workflows can be managed through either command-line options or a configuration file, with command-line options taking precedence when both are specified. The configuration file approach is recommended for complex calculations as it provides better reproducibility and documentation of calculation parameters [22].

Detailed VASP Integration Protocol

The integration between Phonopy and VASP represents one of the most refined and widely used interfaces in the ecosystem. A detailed protocol for this integration exemplifies the general principles that apply to other calculators as well:

• Pre-processing: Generate displacement supercells using the command: phonopy -d --dim 2 2 3 --pa auto. This creates SPOSCAR (the perfect supercell), phonopy_disp.yaml (displacement information), and POSCAR-{number} files (individual displacement supercells) [4].

• VASP Calculation Setup: For each POSCAR-{number} file, create a corresponding VASP calculation directory (disp-00{number}) with the following key INCAR settings:

  The critical setting is IBRION = -1 which prevents ionic relaxation, preserving the displaced structures exactly as generated by Phonopy [4].

• Force Sets Creation: After completing all VASP calculations, collect the force information using: phonopy -f disp-{001..003}/vasprun.xml. This command generates the FORCE_SETS file containing the displacement-force relationships needed for force constant calculation [4].

• Post-processing: With the FORCE_SETS file created, various phonon properties can be calculated:

  • Phonon DOS: phonopy-load --mesh 20 20 20 -p
  • Thermal properties: phonopy-load --mesh 20 20 20 -t
  • Projected DOS: phonopy-load --mesh 20 20 20 --pdos "1 2, 3 4 5 6" -p
  • Band structure: phonopy-load --band "0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0" -p [4]

The following diagram illustrates this integrated workflow:

Essential Research Reagents: Computational Tools and Parameters

Key Configuration Parameters

The accuracy and efficiency of phonon calculations within the Phonopy ecosystem depend critically on the appropriate selection of computational parameters. Table 2 summarizes the essential "research reagents" – the key parameters and their roles in property prediction.

Table 2: Essential Phonopy Configuration Parameters and Their Functions

Parameter/Tag	Function	Default Value	Impact on Results
DIM	Defines supercell dimensions for finite displacement method	--	Larger supercells improve q-space sampling but increase computational cost
DISPLACEMENT_DISTANCE	Controls finite displacement amplitude	0.01 Å (0.02 Bohr for some codes)	Too small: numerical errors; Too large: anharmonic effects
PRIMITIVE_AXES	Transformation matrix from conventional to primitive cell	Identity matrix	Reduces computational cost by working in primitive basis
BORN	Contains Born charges and dielectric tensor for NAC	--	Essential for accurate phonons in polar materials
MESH	q-point sampling mesh for DOS/thermal properties	--	Determines convergence of thermodynamic properties
BAND	Defines q-point path for band structure	--	Visualizes phonon dispersion along high-symmetry paths

The DIM tag is particularly important as it controls the supercell size used for finite displacements. This can be specified as three integers for a diagonal supercell matrix (e.g., DIM = 2 2 2) or nine integers for a full supercell matrix that can describe more complex supercell shapes. The choice of supercell dimensions involves a trade-off between computational cost and the ability to capture long-range interactions, with common practice suggesting dimensions that ensure at least 10 Å in each direction [22].

The PRIMITIVE_AXES tag allows transformation from the input unit cell to the primitive cell, which reduces the computational cost by minimizing the number of atoms in the primitive basis. The PRIMITIVE_AXES = AUTO option, introduced in v1.14.0, automatically selects an appropriate transformation matrix, simplifying setup and ensuring optimal performance [22].

Advanced Displacement Schemes

Beyond basic finite-difference approaches, Phonopy implements sophisticated displacement generation algorithms that optimize computational efficiency:

• CREATE_DISPLACEMENTS: When set to .TRUE., this tag triggers the generation of supercells with displacements as a post-process of the initial setup. The number and pattern of displacements are controlled by additional parameters including DISPLACEMENT_DISTANCE, DIAG, and PM [22].

• DIAG: When set to .FALSE., this tag restricts displacements to directions along the lattice vectors only, avoiding diagonal displacements. This approach is recommended when one lattice parameter is significantly different from the others, as it generates more numerically stable force constants [22].

• PM: This parameter controls the generation of plus-minus displacement pairs. The default PM = AUTO provides a balanced approach that includes symmetrically distinct opposite displacements where needed, improving force constant accuracy without unnecessary redundancy [22].

• RANDOM_DISPLACEMENTS: This advanced feature generates a specified number of supercells with random atomic displacements, useful for sampling configuration space more broadly. When combined with DISPLACEMENT_DISTANCE_MAX, it creates displacements with varying amplitudes, providing more diverse sampling of the potential energy surface [22].

Workflow Visualization: From Structure to Properties

The complete pathway from initial crystal structure to final phonon properties involves multiple stages with specific inputs and outputs at each step. The following diagram maps this comprehensive workflow, highlighting the role of Phonopy as the coordinating engine:

For thermal conductivity calculations requiring three-phonon interactions, Phono3py implements a specialized workflow that extends the basic Phonopy approach:

Emerging Integration with Machine Learning Interatomic Potentials

The Phonopy ecosystem is increasingly interfacing with machine learning interatomic potentials (MLIPs) that offer significant computational advantages for complex systems while maintaining quantum-mechanical accuracy. Recent benchmarking studies have evaluated universal MLIPs (uMLIPs) for phonon property prediction, revealing that some models achieve high accuracy in predicting harmonic phonon properties, though substantial variations exist between different architectures [59].

The integration workflow with MLIPs follows similar patterns as with traditional force calculators but with important distinctions:

• Model Selection: Choosing appropriate MLIPs based on their proven performance for phonon properties. Current evidence suggests that models like M3GNet, CHGNet, MACE-MP-0, and fine-tuned versions of EquiformerV2 show promising results for lattice dynamics calculations [59] [35].

• Force Evaluation: Using MLIPs to compute forces on displacement supercells instead of traditional electronic structure codes. This approach can reduce computational cost by several orders of magnitude while maintaining acceptable accuracy for many applications.

• Validation: Comparing MLIP-predicted phonon properties with reference DFT calculations for representative materials to establish error bounds and applicability domains.

A recent study demonstrated this approach successfully by employing a fine-tuned EquiformerV2 model to screen 921 dynamically stable Na-containing structures from an initial set of 3903 candidates, calculating their lattice dynamics features and molecular dynamics simulations to establish correlations between phonon properties and ionic conductivity [35].

The critical consideration when using MLIPs with Phonopy is that many universal machine learning potentials are trained primarily on equilibrium or near-equilibrium geometries, which may limit their accuracy for large-displacement phonon calculations. Models that incorporate off-equilibrium data through active learning or augmented datasets generally show superior performance for force constant calculations [59].

Best Practices and Validation Protocols

Convergence Testing and Quality Assurance

Robust phonon calculations require systematic convergence testing to ensure results are physically meaningful. Key parameters to test for convergence include:

• Supercell size: Test increasing supercell dimensions until phonon frequencies at high-symmetry points remain unchanged within acceptable tolerances (typically 0.1-0.5 THz).

• q-point mesh: For DOS and thermal properties, increase the mesh density until properties stabilize. A common practice is to use a mesh density that yields at least 1000 q-points per reciprocal atom.

• Displacement distance: Verify that phonon frequencies are insensitive to small variations in displacement distance (e.g., between 0.01-0.03 Å).

• Force convergence criteria: Ensure electronic structure calculations are fully converged with respect to energy, forces, and stress, as incomplete convergence introduces noise in force constants.

Validation against experimental data or higher-level calculations is essential, particularly when exploring new materials systems. Common validation approaches include comparing calculated phonon frequencies with experimental Raman or infrared spectra, checking for imaginary frequencies that may indicate dynamical instability, and verifying that thermal expansion coefficients match experimental values in quasi-harmonic approximations.

Data Management and Reproducibility

With the complexity of workflows involving multiple software packages and parameter sets, maintaining reproducibility requires careful data management. Recent initiatives have proposed metadata schemas for lattice thermal conductivity calculations that align with FAIR (Findable, Accessible, Interoperable, and Reusable) principles [68].

The recommended practice includes documenting:

• Software versions: Phonopy/Phono3py versions and external code versions
• Input parameters: Complete set of calculation parameters including those passed to external codes
• Workflow provenance: Relationships between generated files and calculation steps
• Convergence tests: Results of parameter convergence studies

This documentation facilitates both reproducibility and data sharing, addressing a critical challenge in computational materials science where published results often lack sufficient detail for exact replication [68].

The Phonopy ecosystem provides a sophisticated yet adaptable framework for phonon property prediction through its extensive interfaces with electronic structure codes and emerging machine learning potentials. Its modular architecture separates the core phonon calculation routines from force calculator specifics, creating a versatile platform that continues to incorporate new computational methodologies while maintaining consistent workflows.

The protocols and application notes detailed in this work highlight the importance of appropriate parameter selection, systematic validation, and comprehensive documentation for obtaining reliable results. As computational materials science increasingly emphasizes high-throughput screening and machine learning acceleration, the robust interoperability standards established by Phonopy position it as a critical infrastructure component for next-generation materials discovery.

The ongoing development of both Phonopy and Phono3py, coupled with their expanding integration with machine learning approaches, ensures that these tools will continue to enable researchers to unravel the complex relationships between atomic structure, lattice dynamics, and macroscopic material properties across diverse chemical and materials systems.

Conclusion

Phonopy has established itself as a critical, robust, and accessible tool for performing first-principles phonon calculations, enabling researchers to accurately predict vibrational, thermodynamic, and dynamic stability properties of materials. As demonstrated, its application ranges from foundational studies to the high-throughput screening of advanced functional materials like superionic conductors. The future of phonon computations is being shaped by a tighter integration with machine learning, as evidenced by the development of universal ML interatomic potentials that show promising accuracy for harmonic properties. Furthermore, the growing ecosystem around Phonopy, including its interfaces with codes for calculating anharmonic properties like lattice thermal conductivity, ensures its continued relevance. For the materials science community, mastering Phonopy is not just about computing phonons—it is about unlocking a deeper, quantum-mechanical understanding of the lattice dynamics that govern material behavior and performance from the ground up.

References

• 1. First principles phonon calculations in materials science [https://www.sciencedirect.com/science/article/pii/S1359646215003127]
• 2. anharmonicity and Phonon diffusion in emergent energy ionic materials [https://dukespace.lib.duke.edu/items/b787a6a7-6f49-4c77-8ee0-9fc17cdffa1b]
• 3. 2 Running phonon calculations [https://www.tcm.phy.cam.ac.uk/castep/Phonons_Guide/2-sec:examples.html]
• 4. VASP & phonopy calculation [https://phonopy.github.io/phonopy/vasp.html]
• 5. AI-powered exploration of molecular vibrations, phonons, and ... [https://pubs.rsc.org/en/content/articlehtml/2025/dd/d4dd00353e]
• 6. Welcome to phonopy [https://phonopy.github.io/phonopy/]
• 7. Phonopy API for Python [https://phonopy.github.io/phonopy/phonopy-module.html]
• 8. A quick guide for phonopy [https://materials-lab.io/blog/2021/11/19/phonopy-guide/]
• 9. phonopy — Marianetti Group Documentation [http://grandcentral.apam.columbia.edu:5555/documentation/phonopy/phonopy.html]
• 10. Examples — Phonopy v.2.44.0 [https://phonopy.github.io/phonopy/examples.html]
• 11. First principles investigation of the structural, anisotropic ... [https://www.sciencedirect.com/science/article/pii/S2949822825004113]
• 12. First-principles phonon physics using the Pheasy code [https://arxiv.org/html/2508.01020v1]
• 13. 5.3. ALM: Theoretical background - ALAMODE - Read the Docs [https://alamode.readthedocs.io/en/latest/almdir/formalism_alm.html]
• 14. Interatomic Force Constants in Periodic Solids from Density ... [https://www.sciencedirect.com/science/article/abs/pii/S0065327608604384]
• 15. Efficient Anharmonic Lattice Dynamics Calculations of ... [https://education.molssi.org/2020-software-fellow-posters/giuseppe-barbalinardo/]
• 16. A high-throughput framework for lattice dynamics [https://www.nature.com/articles/s41524-024-01437-w]
• 17. Quasi harmonic approximation — Phonopy v.2.44.0 [https://phonopy.github.io/phonopy/qha.html]
• 18. Thermal properties in the harmonic approximation - hiPhive [https://hiphive.materialsmodeling.org/get_started/compute_harmonic_thermal_properties.html]
• 19. Formulations — Phonopy v.2.44.0 - GitHub Pages [https://phonopy.github.io/phonopy/formulation.html]
• 20. Output files — Phonopy v.2.44.0 [https://phonopy.github.io/phonopy/output-files.html]
• 21. Thermal Properties Calculated with Phonopy [https://mattermodeling.stackexchange.com/questions/12234/thermal-properties-calculated-with-phonopy]
• 22. Setting tags — Phonopy v.2.44.0 [https://phonopy.github.io/phonopy/setting-tags.html]
• 23. A Practical Guide to Frozen Phonon Calculations [https://sites.google.com/site/wolvertonresearchgroup/resources/a-practical-guide-to-frozen-phonon-calculations]
• 24. Advanced > Supercell [https://docs.mat3ra.com/materials-designer/header-menu/advanced/supercell/]
• 25. DMFT_W2K Tutorial 8: FeSe Phonons [http://hauleweb.rutgers.edu/tutorials/Tutorial_phonons.html]
• 26. PH-NODE: A DFPT and finite displacement supercell ... [https://arxiv.org/html/2403.12643v1]
• 27. Phonon Transport Tutorial — Phoebe 1.1.0 documentation [https://phoebe.readthedocs.io/en/develop/tutorials/phononTransport.html]
• 28. VASP-DFPT & phonopy calculation [https://phonopy.github.io/phonopy/vasp-dfpt.html]
• 29. DFTB+ & phonopy calculation [https://phonopy.github.io/phonopy/dftb%2B.html]
• 30. Automated optimization and uncertainty quantification of ... [https://www.nature.com/articles/s41524-024-01388-2]
• 31. How Accurate Are DFT Forces? Unexpectedly Large ... [https://arxiv.org/html/2510.19774v1]
• 32. Convergence problems in DFT [https://storion.ru/all/dft/dft-conv]
• 33. Geometry optimization — AMS 2025.1 documentation - SCM [https://www.scm.com/doc/AMS/Tasks/Geometry_Optimization.html]
• 34. Computing the phonon dispersion and DOS [https://www.vasp.at/wiki/Computing_the_phonon_dispersion_and_DOS]
• 35. Machine-learning-assisted discovery of lattice dynamics ... [https://pubs.rsc.org/en/content/articlehtml/2025/mh/d5mh01176k]
• 36. Change Log — Phono3py v.3.19.4 [https://phonopy.github.io/phono3py/changelog.html]
• 37. Quasi-harmonic approximation [https://en.wikipedia.org/wiki/Quasi-harmonic_approximation]
• 38. Thermodynamics of solids including anharmonicity through ... [https://www.nature.com/articles/s41524-024-01447-8]
• 39. Article Phonon-promoted superionic conduction in fluorite ... [https://www.sciencedirect.com/science/article/pii/S2451929423004047]
• 40. Efficient First - Principles Framework for Overdamped Phonon ... [https://arxiv.org/html/2502.07217v2]
• 41. Phonon website - Henrique Miranda [https://henriquemiranda.github.io/phononwebsite/]
• 42. vasp - The Appearance of Imaginary Frequencies in Lattice ... [https://mattermodeling.stackexchange.com/questions/11991/the-appearance-of-imaginary-frequencies-in-lattice-vibration-calculations-with-i]
• 43. vasp - Inconsistencies in Phonon Dispersion Results Using ... [https://mattermodeling.stackexchange.com/questions/12585/inconsistencies-in-phonon-dispersion-results-using-phonopy-dfpt-with-different-s]
• 44. Evaluating Computational Shortcuts in Supercell-Based ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC7205391/]
• 45. How big should a supercell be in phonon calculations? [https://mattermodeling.stackexchange.com/questions/1480/how-big-should-a-supercell-be-in-phonon-calculations]
• 46. Implementation strategies in phonopy and phono3py | CoLab [https://colab.ws/articles/10.1088%2F1361-648x%2Facd831]
• 47. Accelerating high-throughput phonon calculations via ... [https://www.sciencedirect.com/science/article/abs/pii/S2542529325000446]
• 48. EPIC STAR: a reliable and efficient approach for phonon [https://www.nature.com/articles/s41524-020-0316-7]
• 49. vasp - Phonon calculations using DFPT - Matter Modeling Stack... [https://mattermodeling.stackexchange.com/questions/12494/phonon-calculations-using-dfpt]
• 50. Work flow — Phonopy v.2.44.0 [https://phonopy.github.io/phonopy/workflow.html]
• 51. Ph3pyWF: An automated workflow software package for ... [https://www.sciencedirect.com/science/article/pii/S0010465524003643]
• 52. First-Principles Calculations of Phonon Dispersion and ... [https://bio-protocol.org/exchange/minidetail?id=3800856&type=30]
• 53. "Input forces are not enough to calculate force constants ... [https://mattermodeling.stackexchange.com/questions/1922/input-forces-are-not-enough-to-calculate-force-constants-phonopy-error]
• 54. get_positions_sent_by_rot_inv is occasionally given inputs ... [https://github.com/atztogo/phonopy/issues/65]
• 55. Phonons from finite differences [https://vasp.at/wiki/Phonons_from_finite_differences]
• 56. “One defect, one potential” strategy for accurate machine ... [https://arxiv.org/html/2509.00498v1]
• 57. A critical review of machine learning interatomic potentials ... [https://www.oaepublish.com/articles/jmi.2025.17]
• 58. Machine learning interatomic potential: Bridge the gap ... [https://www.sciencedirect.com/science/article/pii/S2589004224008952]
• 59. Universal machine learning interatomic potentials are ... [https://www.nature.com/articles/s41524-025-01650-1]
• 60. Accelerating point defect photo-emission calculations with ... [https://www.nature.com/articles/s41524-025-01820-1]
• 61. Learning Smooth and Expressive Interatomic Potentials for ... [https://arxiv.org/html/2502.12147v1]
• 62. Universal Machine Learning Potential for Systems with ... [https://arxiv.org/html/2508.15614v1]
• 63. Workflow — Phono3py v.3.19.4 [https://phonopy.github.io/phono3py/workflow.html]
• 64. THERMACOND, a code to compute lattice thermal ... [https://www.nature.com/articles/s41524-025-01673-8]
• 65. VASP & phono3py calculation [https://phonopy.github.io/phono3py/vasp.html]
• 66. Lattice thermal conductivity using Phono3py and VASP [http://s-utkarsh.github.io/k_lattice/]
• 67. Auxiliary tools — Phono3py v.3.19.4 [https://phonopy.github.io/phono3py/auxiliary-tools.html]
• 68. A metadata schema for lattice thermal conductivity from first ... [https://www.oaepublish.com/articles/jmi.2022.20]