Temperature Dependent Effective Potential

• TDEP is a method that replaces static force constants with temperature-renormalized parameters derived from MD or quantum simulations to capture anharmonic behavior.
• It employs linear least-squares fitting of force-displacement data while enforcing crystal symmetries to efficiently extract higher-order interatomic force constants.
• The framework enables accurate predictions of phonon frequencies, elastic constants, and thermal transport properties in strongly anharmonic and quantum-dominated materials.

The temperature dependent effective potential (TDEP) framework is a first-principles approach for capturing thermodynamic, vibrational, and transport properties of solids with strong anharmonicity by extracting effective force constants directly from finite-temperature atomistic trajectories. Rather than relying on the static, zero-temperature potential energy surface, TDEP replaces the bare Taylor expansion with a temperature-renormalized model Hamiltonian, whose parameters are fitted via statistical inference on molecular dynamics (MD) or quantum path-integral MD data. This variational, data-driven pipeline enables a self-consistent, quantitatively accurate description of phonons, thermodynamics, and elastic and transport coefficients across the full temperature range, including dynamically stabilized and quantum-fluctuation-dominated regimes. The method has been generalized to extract higher-order interatomic force constants, extended for quantum nuclei using path-integral sampling, and incorporated into multiple implementations for materials ranging from elemental solids to complex oxides and paramagnetic compounds (Hellman et al., 2013, Hellman et al., 2013, Mozafari et al., 2016, Geng, 2022, Folkner et al., 2024).

1. Mathematical Formalism and Model Hamiltonian

TDEP begins from a general expansion of the Born–Oppenheimer or Born–von Kármán potential energy $U(\{ \mathbf{u} \})$ in atomic displacements about temperature-dependent equilibrium positions: $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$ where $u_i$ is the displacement of atom $i$, and $\Phi^{(n)}$ are the $n$th-order interatomic force constants (IFCs), with $U_0$ an energy reference (Hellman et al., 2013, Folkner et al., 2024).

In TDEP, these IFCs are replaced by effective, temperature-dependent parameters $\Phi^{(n)}(T)$, defining an optimized model Hamiltonian

$$\hat H_{TDEP}(T) = U_0 + \sum_i \frac{p_i^2}{2m_i} + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}(T)u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}(T)u_iu_ju_k + \dots$$

which is fit to force-displacement data sampled from MD or quantum MD at temperature $T$. By construction, $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$0 incorporates anharmonic and quantum fluctuations renormalized at that temperature (Hellman et al., 2013, Geng, 2022, Folkner et al., 2024).

2. Fitting Procedure and Statistical Extraction from Simulation

To extract the $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$1 effective force constants, TDEP employs a linear least-squares minimization of the difference between the true forces (from ab initio MD or PIMD) and the model forces, for a set of $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$2 uncorrelated snapshots $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$3 (Hellman et al., 2013, Hellman et al., 2013, Folkner et al., 2024): $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$4 where the model force for atom $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$5 is

$$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$6

The covariance structure and force constants are reduced by imposing crystal translation, point group, and permutation symmetries, as well as acoustic sum rules. For higher-order IFCs, the least-squares system is enlarged to include cubic and (optionally) quartic terms (Hellman et al., 2013, Han et al., 2023). In quantum solids or at low temperature, path-integral MD sampling is essential to capture nuclear quantum effects without artifacts (Geng, 2022).

Variational, self-consistent algorithms (sometimes with re-weighting for all prior samples) are used to ensure convergence and statistical robustness (Bichelmaier et al., 2021). Symmetry-imposed variants (e.g., SIFC-TDEP for paramagnetic or disordered systems) enforce ideal lattice point-group averaging for force constants, increasing efficiency and statistical convergence (Mozafari et al., 2016).

3. Phonons, Elastic Constants, and Free Energy from Effective Force Constants

The fitted effective harmonic force constants $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$7 define the temperature-dependent dynamical matrix: $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$8 Diagonalization yields temperature-renormalized phonon frequencies $$U(\{ \mathbf{u} \}) = U_0 + \sum_i \Phi_i u_i + \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)}u_iu_j + \frac{1}{3!}\sum_{ijk}\Phi_{ijk}^{(3)}u_iu_ju_k + \dots$$9 and eigenvectors. The Helmholtz free energy is then constructed from harmonic and anharmonic (residual) terms: $u_i$0

$u_i$1

with $u_i$2 the phonon density of states. Optionally, further anharmonic corrections are computed via direct averaging or perturbative expansions in higher-order IFCs (Hellman et al., 2013, Bichelmaier et al., 2021).

The elastic constants at finite temperature are extracted using the long-wavelength limit (Born–Huang relations) of the dynamical matrix or via direct sums over IFCs: $u_i$3 Voigt–Reuss–Hill averaging provides polycrystalline moduli (Mozafari et al., 2016, Folkner et al., 2024).

4. Incorporation of Anharmonicity, Quantum Effects, and Sampling Schemes

TDEP, by design, captures the leading effect of strong anharmonicity via the thermal renormalization of force constants. Extensions have been developed to explicitly fit third- and fourth-order IFCs, enabling direct calculation of Grüneisen parameters, three- and four-phonon scattering rates, and phonon linewidths/lifetimes relevant for transport (Hellman et al., 2013, Han et al., 2023, Folkner et al., 2024). The quantum-classical boundary is addressed by full temperature-dependent potential (FTDP) methods, which use ab initio path-integral MD sampling to extract $u_i$4 unbiased with respect to the quantum distribution, ensuring physicality below the Debye temperature and for light-element compounds (Geng, 2022, Folkner et al., 2024).

The workflow for quantum TDEP (FTDP) is as follows:

1. Generate PIMD trajectories with a machine-learned or ab initio potential.
2. For each snapshot and “bead”, record atomic displacements and forces (excluding spring terms).
3. Fit the force constants as in the classical case, discarding bead–bead spring contributions (Geng, 2022, Folkner et al., 2024).
4. Check convergence with supercell and bead number; statistical averaging and error bars are recommended.

5. Applications and Benchmark Results

TDEP and its variants have enabled quantitatively accurate predictions for multiple classes of materials:

• Anharmonic metals and dynamically unstable phases: Dynamical stabilization of bcc Zr, which is harmonically unstable at $u_i$5 but stabilized by TDEP-renormalized force constants at $u_i$6 K, producing real phonon frequencies, correct free energy surfaces, and finite-$u_i$7 EOS (Hellman et al., 2013).
• Paramagnetic and disordered magnets: Assessment of elastic constants in high-$u_i$8 paramagnetic CrN via SIFC-TDEP + DLM-MD, revealing a 14% reduction in Young’s modulus from 300 to 1200 K and increasing isotropy with $u_i$9 (Mozafari et al., 2016).
• Quantum and ultra-anhamonic solids: In metallic hydrogen, FTDP enabled detection of colossal nuclear quantum effects and the breakdown of the quasi-phonon picture (Geng, 2022).
• Thermal transport: Calculation of lattice thermal conductivity in ceramics (e.g., CeO$i$0, MgO) requires temperature-renormalized phonons and explicit inclusion of both three- and four-phonon processes, with TDEP dramatically improving agreement with experiment (Han et al., 2023, Folkner et al., 2024).
• Critical and soft-mode systems: Modified TDEP formalism remedies unphysical limitations (e.g., absence of diverging susceptibility at second-order transitions), enabling correct dynamical and static response functions near ferroelectric or structural instabilities (Monacelli, 2024).

A summary table of selected applications is provided:

System	Extension/Variant	Remarkable Feature/Result
bcc Zr	TDEP, classical	Dynamical stabilization at high T, force localization with T (Hellman et al., 2013)
CrN paramagnet	SIFC-TDEP + DLM-MD	Single/polycrystal $i$1, enhanced isotropy, ~14% $i$2 reduction (Mozafari et al., 2016)
Metallic H	FTDP (quantum)	Colossal NQE, phonon nonconvergence, lattice quantum fluctuations (Geng, 2022)
Si, FeSi	TDEP + cubic	Mode Grüneisen parameters, $i$3-dependent lifetimes (Hellman et al., 2013)
CeO$i$4, MgO	TDEP + quartic	Four-phonon rates, $i$5-renormalized $i$6 matching experiment (Han et al., 2023)
PbTe, CsSnI$i$7	Modified TDEP	Critical-mode stabilization, correct spectral and static response (Monacelli, 2024)

6. Advantages, Limitations, and Recent Developments

TDEP provides direct inclusion of anharmonic effects, captures thermodynamic stabilization of otherwise unstable phases, and is extensible to arbitrary order with symmetry-imposed reduction in parameter space (Hellman et al., 2013, Mozafari et al., 2016, Folkner et al., 2024). The method is robust with respect to statistical noise (by leveraging symmetry and block-averaging), highly parallelizable, and compatible with machine-learned force fields (Folkner et al., 2024).

Limitations include: (i) standard TDEP models are quadratic (or, with extensions, cubic/quartic) and so only implicitly capture higher-order correlations unless higher IFCs are explicitly fit (Hellman et al., 2013, Han et al., 2023); (ii) classical TDEP is invalid below the Debye temperature without quantum sampling; (iii) conventional TDEP does not produce divergent susceptibilities at critical points—a new correction (subtracting the static bubble from the self-energy) is required to restore physicality in the dynamical spectra (Monacelli, 2024).

Recent work incorporates PGIMD for quantum sampling, machine-learning potentials for efficiency in large systems, self-consistent ensemble reweighting for improved convergence, and hybrid approaches to join TDEP with self-consistent phonon or renormalization group frameworks (Bichelmaier et al., 2021, Geng, 2022, Folkner et al., 2024, Monacelli, 2024).

7. Extensions and Variants: SIFC-TDEP, FTDP, and Beyond

• SIFC-TDEP: Enforces full point-group symmetry of the underlying lattice for efficient force-constant reduction and faster convergence, especially in large/disordered or magnetic supercells (Mozafari et al., 2016).
• FTDP (Full Temperature-Dependent Potential): Generalizes TDEP by fitting to path-integral ab initio trajectories, ensuring validity in the quantum limit and for systems with significant nuclear quantum effects (Geng, 2022).
• Step-by-step higher-order fitting: Successively fits cubic and quartic terms to force residuals for improved accuracy in scattering and thermal transport properties (Hellman et al., 2013, Han et al., 2023).
• Modified self-energies for criticality: Corrects the static limit of bubble diagrams to ensure divergence of susceptibility at phase transitions and proper frequency dependences of dynamical response (Monacelli, 2024).
• Integration with machine-learned potentials and compressive sensing: Enables tractable extraction of effective force constants from long, large-scale quantum MD trajectories (Folkner et al., 2024).

TDEP thus remains a leading approach for first-principles investigation of strongly anharmonic, thermally activated, and quantum-fluctuating solids across the material spectrum. Its implementation, validation, and further development continue to be active research areas (Hellman et al., 2013, Hellman et al., 2013, Mozafari et al., 2016, Bichelmaier et al., 2021, Geng, 2022, Han et al., 2023, Monacelli, 2024, Folkner et al., 2024).