Pheasy: First-Principles Phonon Package

• Pheasy is a Python package for lattice dynamics that uses a Taylor expansion framework to extract harmonic and anharmonic interatomic force constants with robust regression methods.
• It unifies workflows for phonon dispersions, temperature renormalization, and thermal transport by enforcing exact symmetry and invariance constraints within a single computational ecosystem.
• Pheasy integrates with DFT codes and external BTE solvers, employing LASSO/adaptive LASSO and multipolar electrostatics to accurately model complex anharmonic materials.

Searching arXiv for the specified paper and a few directly related keywords. Searching arXiv. Pheasy is a Python package for first-principles lattice dynamics that represents the Born–Oppenheimer lattice potential of crystalline solids as a Taylor expansion to arbitrarily high order, extracts interatomic force constants (IFCs) from force–displacement datasets using LASSO or adaptive LASSO regression, and uses those IFCs as a unified basis for harmonic and anharmonic phonon calculations, including phonon dispersions, temperature renormalization, thermodynamics, and thermal transport with three- and four-phonon scattering (Lin et al., 1 Aug 2025). It is designed for the regime in which conventional finite-difference or DFPT workflows become challenging beyond third-order anharmonicity because of the combinatorial growth of higher-order IFCs, and it is explicitly positioned as a phonon code ecosystem that connects multiple simulation platforms.

1. Definition, scope, and target problems

Pheasy targets parameter-free calculations of lattice dynamics from first principles in settings where high-order anharmonicity, structural instability, or thermal transport require a systematic treatment of IFCs beyond the harmonic and cubic levels (Lin et al., 1 Aug 2025). Its stated goals are to reconstruct the Born–Oppenheimer potential energy surface of crystals as a high-order Taylor expansion in atomic displacements, to extract harmonic and anharmonic IFCs of any order from force–displacement data with robust regression, and to compute a wide range of phonon-related properties.

The code is aimed at lattice dynamics, anharmonic renormalization and phase stability, thermodynamics, and lattice thermal conductivity. In this framework, second-order IFCs determine the harmonic dynamical matrix and therefore phonon frequencies and eigenvectors; third-order IFCs determine three-phonon interactions relevant to phonon lifetimes and thermal conductivity in weakly anharmonic solids; and fourth- and higher-order IFCs encode higher-order phonon–phonon interactions, including four-phonon scattering and the inputs required by non-perturbative renormalization methods such as SCP and SCHA.

A central feature is that all IFC orders are treated within one expansion framework. This unifies workflows that are often separated across codes: harmonic spectra, thermodynamic quantities, temperature-dependent effective harmonic IFCs, finite-temperature structural optimization, and IFC export to external Boltzmann transport solvers. Pheasy is also described as complementing Phonopy, Phono3py, thirdorder.py, CSLD, Alamode, hiPhive, and DFPT implementations by pushing toward systematic high order, rigorously constrained IFC fitting, long-range multipolar electrostatics including 2D corrections, and interoperability with external tools.

2. Lattice-potential expansion and cluster formalism

The theoretical basis is a Taylor expansion of the lattice potential in atomic displacements $u_{\mathbf a}$ from equilibrium positions $\mathcal R_{\mathbf a}$ (Lin et al., 1 Aug 2025):

$$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$

with bold indices $\mathbf a \equiv \{a,i\}$ combining atom index and Cartesian direction. The $n$-th order IFCs are defined by

$$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$

and at equilibrium $\Phi_{\mathbf a}=0$. Differentiation gives the force expansion,

$$F_{\mathbf a} = -\frac{\partial V}{\partial u_{\mathbf a}} = -\Phi_{\mathbf{ab}}u_{\mathbf b} -\frac{1}{2}\Phi_{\mathbf{abc}}u_{\mathbf b}u_{\mathbf c} -\cdots.$$

Pheasy adopts a cluster-expansion notation, following CSLD, in which

$$V=\frac{1}{\alpha!}\Phi_I(\alpha)u_I^\alpha, \qquad F_{\mathbf a}=-\frac{1}{\alpha!}\Phi_I(\alpha)\,\partial_{\mathbf a}u_I^\alpha.$$

Here a cluster $\alpha$ is a multiset of lattice sites, $\mathcal R_{\mathbf a}$0 is the set of Cartesian components, $\mathcal R_{\mathbf a}$1 is the corresponding displacement monomial, $\mathcal R_{\mathbf a}$2 is the cluster order, and $\mathcal R_{\mathbf a}$3 collects combinatorial factors. This notation is important because it allows symmetry reduction and real-space truncation to be handled in a compact way. The working assumption is that anharmonic IFCs are localized, so the expansion can be truncated at a chosen neighbor-shell radius.

The same formalism underlies both direct property calculations and finite-temperature renormalization. In particular, the effective harmonic IFCs used in self-consistent phonon theory are written as

$\mathcal R_{\mathbf a}$4

When odd moments and structural relaxation are neglected, this reduces to the loop-only SCP equation,

$\mathcal R_{\mathbf a}$5

This establishes a direct connection between the fitted high-order IFCs and temperature-dependent phonon spectra.

3. Constraint structure and IFC extraction

IFC extraction is formulated as a linear regression problem from DFT forces or energies (Lin et al., 1 Aug 2025). After symmetry reduction, the force and energy equations are written as

$\mathcal R_{\mathbf a}$6

$\mathcal R_{\mathbf a}$7

In this notation, $\mathcal R_{\mathbf a}$8 is the concatenated force vector, $\mathcal R_{\mathbf a}$9 is the total-energy vector when energy fitting is used, $$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$0 contains the full IFC tensor components for representative clusters, $$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$1 is the null-space matrix that imposes symmetry and invariance constraints, and $$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$2 contains the independent parameters. The core inverse problem is therefore

$$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$3

A defining aspect of Pheasy is the exact enforcement of several constraints. Space-group symmetry imposes covariance of IFCs under crystal symmetry operations and groups clusters into orbits with one representative per orbit. Permutation symmetry expresses the interchange symmetry of mixed derivatives. Translational invariance imposes acoustic sum rules at each IFC order so that uniform translations do not generate forces. Rotational invariance and equilibrium stress conditions are enforced on harmonic IFCs, including Born–Huang rotational invariance; this is stated to be crucial for the correct quadratic dispersion of 2D flexural ZA modes, the correct twisting acoustic modes in 1D, and consistent coupling to stress in polar systems. For infrared-active materials, polar invariance conditions are also enforced so that long-range Coulomb terms are treated consistently.

All of these linear conditions are assembled into a null-space problem. Pheasy replaces earlier approximate schemes based on iterative null spaces with tolerances by an iterative maximal-pivot Gauss–Jordan elimination that is described as producing an exact null space without ad hoc tolerances and remaining stable in low-symmetry structures and in the presence of rotational or stress constraints. This design choice is methodologically significant because the reliability of high-order IFC extraction depends strongly on the correctness of the constrained parameter basis.

Regression is performed with ordinary least squares, LASSO, or adaptive LASSO. OLS minimizes

$$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$4

while LASSO adds an $$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$5 penalty,

$$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$6

and adaptive LASSO uses feature-dependent weights,

$$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$7

The rationale given in the paper is that high-order IFCs are physically short-ranged and sparse, so sparse regression is particularly appropriate in the underdetermined regime. Before regression, the columns of $$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$8 are standardized to zero mean and unit standard deviation.

4. Training ensembles, long-range interactions, and supported calculations

The quality of the fit depends on the generation of informative and weakly correlated training structures (Lin et al., 1 Aug 2025). Pheasy supports three principal strategies. For low-order IFCs, random fixed-magnitude displacements are used, with approximately $$V = V_0 + \Phi_{\mathbf a}u_{\mathbf a} + \frac{1}{2!}\Phi_{\mathbf{ab}}u_{\mathbf a}u_{\mathbf b} + \frac{1}{3!}\Phi_{\mathbf{abc}}u_{\mathbf a}u_{\mathbf b}u_{\mathbf c} + \cdots,$$9 Å for harmonic IFCs and approximately $\mathbf a \equiv \{a,i\}$0 Å for cubic IFCs. For high-order IFCs and temperature-dependent IFCs, the preferred method is quantum canonical sampling of normal modes, based on the harmonic covariance matrix

$\mathbf a \equiv \{a,i\}$1

with displacements drawn from the Gaussian distribution

$\mathbf a \equiv \{a,i\}$2

Ab initio molecular dynamics snapshots can also be used, but the paper generally recommends Gaussian sampling unless AIMD is specifically required.

Once IFCs are available, Pheasy computes harmonic phonon dispersions and eigenvectors, thermodynamic quantities in harmonic or quasi-harmonic approximations, anharmonic renormalization with real-space SCP, and exports data for thermal-transport calculations. Thermal conductivity is not solved internally; instead, Pheasy provides IFCs to external BTE solvers such as ShengBTE, Phoebe, and GPU_PBTE. The transport equations are written in the usual form,

$\mathbf a \equiv \{a,i\}$3

and

$\mathbf a \equiv \{a,i\}$4

with lifetimes determined by three- and four-phonon scattering rates.

A further technical component is the treatment of long-range Coulomb interactions. In bulk, the long-range contribution to the dynamical matrix is formulated in the long-wavelength limit with analytic multipolar expansions of the charge-density response and dielectric screening. At lowest order, the familiar dipole–dipole term is recovered,

$\mathbf a \equiv \{a,i\}$5

Pheasy also includes quadrupolar and octupolar contributions up to second order in $\mathbf a \equiv \{a,i\}$6, and for 2D materials it adopts the 2D Coulomb kernel of Sohier and Royo. The operational procedure is to subtract the long-range force contribution from DFT forces before fitting so that the fitted IFCs remain short-ranged, then to add the analytic long-range terms back at the dynamical-matrix stage. This division between short-range regression and analytic electrostatics is one of the code’s defining methodological features.

5. Benchmarks, validation, and methodological conclusions

Three benchmark systems are used to validate the workflow and to compare IFC extraction strategies (Lin et al., 1 Aug 2025).

System	Setup	Selected outcomes
Bulk Si	$\mathbf a \equiv \{a,i\}$7 supercell; expansion to 6th order; adaptive LASSO	RMSE $\mathbf a \equiv \{a,i\}$8 meV/Å; $\mathbf a \equiv \{a,i\}$9 W/mK for 3ph and $n$0 W/mK for 3+4ph
Monolayer WS$n$1	$n$2 supercell; expansion to 6th order; LASSO	RMSE $n$3 meV/Å; $n$4 W/mK for 3ph and $n$5 W/mK for 3+4ph
Cubic SrTiO$n$6	$n$7 cubic supercell; SCHA ensembles at 300 K; adaptive LASSO	RMSE $n$8 meV/Å; $n$9 W/mK for 3ph and $$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$0 W/mK for 3+4ph

For bulk silicon, the expansion includes harmonic IFCs from five random displacement structures of $$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$1 Å, high-order terms up to sixth order, and quantum canonical sampling at 300 K with 64 structures per ensemble. The paper reports frequency softening with increasing temperature in loop-only SCP, especially in the optical branches, but less softening than in TDEP including bubble diagrams. The thermal-conductivity calculations from 200 to 1000 K show an approximately $$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$2 reduction at 300 K when four-phonon scattering is added. All tested levels—HA+3ph, HA+3,4ph, IFCs@T+3ph, and IFCs@T+3,4ph—yield similar room-temperature conductivity, which is interpreted as evidence that silicon is weakly anharmonic and that temperature-dependent IFCs have only a minor effect.

For monolayer WS$$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$3, the system is chosen because the large mass contrast produces a large acoustic–optical gap and therefore strong four-phonon scattering is expected. Harmonic IFCs are from DFPT, the expansion again reaches sixth order, and the training structures are generated from canonical sampling at 300 K using DFPT phonons. The reported phonon dispersions between 0 and 1000 K show very weak temperature dependence, consistent with weak anharmonic renormalization. Four-phonon processes reduce the 300 K thermal conductivity by $$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$4, from approximately $$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$5 to approximately $$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$6 W/mK. The paper also notes that previous theory and experiment show a large spread, with experiments around $$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$7–$$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$8 W/mK, and attributes this mismatch to experimental challenges rather than using it to validate any single theoretical number.

For cubic SrTiO$$\Phi_{\mathbf a_1\cdots \mathbf a_n} \equiv \frac{\partial^n V}{\partial u_{\mathbf a_1}\cdots \partial u_{\mathbf a_n}},$$9, the benchmark addresses strong anharmonicity and a structural phase transition. The harmonic cubic phase is unstable in DFPT, showing imaginary frequencies, but loop-only SCP applied to the fitted fourth-order IFCs produces temperature-dependent dispersions in which all modes are real at $\Phi_{\mathbf a}=0$0, dynamically stabilizing the cubic phase. At $\Phi_{\mathbf a}=0$1, zero-point fluctuations stabilize the $\Phi_{\mathbf a}=0$2 and $\Phi_{\mathbf a}=0$3 soft modes but not the $\Phi_{\mathbf a}=0$4 antiferrodistortive mode. The paper explicitly states that loop-only SCP cannot generate a second-order phase transition by construction because the cubic bubble diagram is essential. It further notes that a transition around 220 K reported by Tadano et al. can be reproduced only under an inconsistent setup using a dense $\Phi_{\mathbf a}=0$5 $\Phi_{\mathbf a}=0$6-grid with a $\Phi_{\mathbf a}=0$7 supercell, and presents this as evidence of the subtlety of the numerical procedure rather than as a confirmation of a physical transition in loop-only SCP.

The benchmark most strongly supports the paper’s methodological conclusion about IFC fitting strategy. In SrTiO$\Phi_{\mathbf a}=0$8, a one-shot extraction that fits second-, third-, fourth-, and higher-order IFCs simultaneously yields SCP frequencies that agree essentially perfectly with direct SCHA effective harmonic IFCs, whereas a cocktail extraction that fixes harmonic IFCs from DFPT or finite displacements and fits only the anharmonic terms gives significantly different dispersions. The discrepancy is stated to be amplified when harmonic IFCs are generated with $\Phi_{\mathbf a}=0$9- or $$F_{\mathbf a} = -\frac{\partial V}{\partial u_{\mathbf a}} = -\Phi_{\mathbf{ab}}u_{\mathbf b} -\frac{1}{2}\Phi_{\mathbf{abc}}u_{\mathbf b}u_{\mathbf c} -\cdots.$$0-grids inconsistent with the supercell DFT used for the anharmonic data. This suggests that cocktail workflows can redistribute physics incorrectly between harmonic and anharmonic terms and thereby bias renormalization and thermal-transport predictions.

6. Software architecture, interoperability, limitations, and outlook

Pheasy is implemented in Python and uses ASE to interface with DFT codes, with Quantum ESPRESSO and VASP currently supported and extension to other ASE-compatible engines described as straightforward (Lin et al., 1 Aug 2025). It uses scikit-learn for LASSO implementations and standard machine-learning utilities. Internally, the code is organized into symmetry and cluster generation, null-space and constraint construction, sensing-matrix assembly and regression, phonon and thermodynamic post-processing, and export interfaces to third-party phonon and BTE solvers.

Interoperability is a major design objective. On the electronic-structure side, the ASE layer connects Pheasy to DFT engines. On the phonon and transport side, the code exports IFCs and dynamical matrices to ShengBTE and its four-phonon extension, Phono3py, Phoebe, GPU_PBTE, and EPW. It is also described as having been integrated into high-throughput frameworks such as atomate2 for large-scale harmonic phonon calculations. In this sense, Pheasy is intended not only as a standalone phonon package but as infrastructure linking multiple stages of a lattice-dynamical workflow.

The recommended usage guidelines are correspondingly specific. For harmonic IFCs only, small random displacements of approximately $$F_{\mathbf a} = -\frac{\partial V}{\partial u_{\mathbf a}} = -\Phi_{\mathbf{ab}}u_{\mathbf b} -\frac{1}{2}\Phi_{\mathbf{abc}}u_{\mathbf b}u_{\mathbf c} -\cdots.$$1 Å are recommended. For third-order IFCs without higher orders, random displacements of approximately $$F_{\mathbf a} = -\frac{\partial V}{\partial u_{\mathbf a}} = -\Phi_{\mathbf{ab}}u_{\mathbf b} -\frac{1}{2}\Phi_{\mathbf{abc}}u_{\mathbf b}u_{\mathbf c} -\cdots.$$2 Å are recommended, with the caution that if fourth-order terms are omitted then the fitted second- and third-order IFCs become effective, renormalized quantities. For fourth order and beyond, there is stated to be no universally safe fixed displacement magnitude; instead, the recommended strategy is a sixth-order expansion including at least onsite and two-body fifth- and sixth-order terms, combined with quantum canonical sampling. For temperature-dependent IFCs, effective harmonic IFCs should first be fitted separately, for example with SCP or SCHA, and then the third- and fourth-order IFCs should be fitted to residual forces after subtracting the effective harmonic contribution. The paper also warns that mixing TDEP- or AIMD-derived effective harmonic IFCs with explicit three-phonon scattering can lead to double counting of bubble contributions.

Several limitations are emphasized. High-order anharmonicity remains computationally demanding because the number of irreducible parameters grows rapidly and convergence with cutoff radius and body-order truncation must be checked carefully. The paper focuses on loop-only SCP and states that a fuller many-body treatment including bubble and tadpole diagrams, together with their consistent use in BTE transport, is only partially addressed and deferred to separate work. The role of long-range Coulomb forces in IFC regression is described as subtle: although Pheasy implements rigorous subtraction of dipolar and higher-order long-range contributions before fitting, tests indicate that for LASSO-fitted IFCs the impact on phonons and thermal conductivity can be small, which the paper notes is contrary to earlier claims and requires more systematic study. A plausible implication is that the numerical benefit of analytic electrostatics may depend strongly on the regularization strategy and the degree of real-space sparsity in the fitted model.

Overall, Pheasy defines a high-order, symmetry-constrained, regression-based framework for first-principles phonon physics in which the Taylor expansion of the lattice potential is elevated from a low-order perturbative device to a general computational representation of harmonic and anharmonic lattice dynamics. Its principal significance lies in the combination of arbitrarily high-order IFC extraction, exact constraint handling, multipolar electrostatics including 2D treatments, and explicit interoperability with the broader phonon and thermal-transport software ecosystem.