Acoustic Sum Rule (ASR) Violation in Phonon Calculations: A Comprehensive Guide from Fundamentals to Troubleshooting

Abstract

This article provides a comprehensive examination of the Acoustic Sum Rule (ASR) in phonon frequency calculations, a critical aspect of lattice dynamics in materials science. It explores the foundational theory behind ASR and the root causes of its violation in density functional theory (DFT) and other computational methods. The content details practical methodologies for performing robust phonon calculations, systematic troubleshooting approaches for resolving non-zero acoustic frequencies and imaginary modes, and modern validation techniques, including the use of machine learning potentials. Aimed at researchers and computational scientists, this guide synthesizes current knowledge to enhance the accuracy and reliability of phonon property predictions for materials design and discovery.

Understanding the Acoustic Sum Rule: Why Translational Invariance is Key to Phonon Calculations

The Acoustic Sum Rule (ASR) represents a fundamental cornerstone in the lattice dynamics of crystalline solids. It is a direct mathematical consequence of the translational invariance of a crystal lattice, which dictates that uniformly displacing every atom in a crystal by the same infinitesimal amount should not change the system's potential energy or generate any net force on the atoms [1]. In the context of phonon calculations—the study of quantized vibrational modes in materials—the ASR ensures the existence of three zero-frequency acoustic modes at the Brillouin zone center (Γ-point). These modes correspond to the three rigid translational degrees of freedom of the entire crystal [1]. The violation of the ASR in practical calculations leads to unphysical, non-zero frequencies for these acoustic modes, compromising the predictive accuracy and physical meaningfulness of the results. This guide details the principles underlying the ASR, the consequences of its violation, and the methodologies for its enforcement, framed within ongoing research on robust phonon frequency calculations.

The formal derivation of the ASR stems from the Born-Huang invariance conditions [2]. For a crystal at equilibrium, the potential energy ( E ) is expanded to second order in atomic displacements ( u_{\varkappa \alpha} ):

[ E = E0 + \sum{\varkappa \alpha} \Phi{\varkappa \alpha} u{\varkappa \alpha} + \frac{1}{2} \sum{\varkappa \alpha, \varkappa' \beta} \Phi{\varkappa \alpha, \varkappa' \beta} u{\varkappa \alpha} u{\varkappa' \beta} ]

Here, ( \Phi{\varkappa \alpha} ) and ( \Phi{\varkappa \alpha, \varkappa' \beta} ) are the first and second-order interatomic force constants (IFCs). The condition of translational invariance requires that the energy is unchanged under a uniform translation ( \mathbf{u} ), leading to the ASR for the second-order IFCs [2]:

[ \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \beta} = 0 ]

This sum rule must be satisfied for the dynamical matrix to correctly yield three acoustic phonons with zero frequency at the Γ-point.

The Physical Basis and Theoretical Foundation

Connection to Translational and Rotational Invariance

The ASR is one part of a broader set of constraints known as the Born-Huang conditions, which enforce both translational and rotational invariance on the interatomic force constants [2]. These invariances are direct consequences of the conservation of total crystal momentum and angular momentum, as stated by Noether's theorem [2]. The complete set of constraints includes:

• Equation (2): ( \sum{\varkappa} \Phi{\varkappa \alpha} = 0 ) (Translational invariance for first-order IFCs)
• Equation (3): ( \sum{\varkappa} \Phi{\varkappa \alpha} \tau{\varkappa \beta} = \sum{\varkappa} \Phi{\varkappa \beta} \tau{\varkappa \alpha} ) (Rotational invariance for first-order IFCs)
• Equation (4): ( \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \beta} = 0 ) (The ASR: Translational invariance for second-order IFCs)
• Equation (5): ( \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \beta} \tau{\varkappa' \gamma} + \delta{\alpha \gamma} \Phi{\varkappa \beta} = \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \gamma} \tau{\varkappa' \beta} + \delta{\alpha \beta} \Phi{\varkappa \gamma} ) (Rotational invariance for second-order IFCs)

Here, ( \tau{\varkappa \alpha} ) is the equilibrium position of atom ϰ, and ( \delta{\alpha\beta} ) is the Kronecker delta. For a system in equilibrium, the first-order IFCs ( \Phi_{\varkappa \alpha} ) (the forces on atoms) are zero, simplifying the rotational condition for the second-order IFCs [2]. These relationships are critical as they link different orders of the expansion of the potential energy surface.

The Critical Role of the ASR in Low-Dimensional Materials

The enforcement of rotational invariance, closely linked to the ASR, becomes critically important in low-dimensional (LD) materials like 1D nanotubes and 2D monolayers [2]. In these systems, the long-wavelength behavior of flexural phonons (out-of-plane bending modes) is highly sensitive to the correct implementation of these invariances. When rotational invariance is broken in the calculation, an unphysical linear dispersion is often observed for the ZA (out-of-plane acoustic) mode in 2D materials, whereas a physically correct quadratic dispersion is obtained only when both translational and rotational invariances are enforced [2]. This is not merely a numerical artifact; it has significant implications for accurately predicting thermal transport, carrier mobility, and thermodynamic stability in low-dimensional systems. Research has shown that applying these corrections can dramatically improve results, transforming dynamically unstable phonon spectra into stable ones for a significant number of 2D materials [2].

ASR Violations: Causes and Consequences in Computational Practice

In an ideal and perfectly converged calculation, the IFCs would naturally satisfy the ASR. However, in practice, various numerical approximations lead to its violation. The primary causes identified in the literature include:

• Discreteness of the FFT Grid: The use of a discrete Fast Fourier Transform (FFT) grid to represent charge density and potentials is cited as the "main and most irreducible" source of ASR violation in plane-wave DFT calculations [3]. This discretization breaks the perfect translational symmetry of the continuous crystal.
• Insufficient Numerical Convergence: Inadequate convergence with respect to key parameters can cause violations. This includes:
  • Phonon Convergence (tr2_ph): The self-consistency threshold for phonon calculations [3].
  • Ground-State Convergence (conv_thr): The convergence threshold for the initial electronic ground state [3].
  • Plane-Wave Cutoff Energy (ecutwfc, ecutrho): An insufficiently high cutoff can lead to spurious forces [1].
  • k-point Sampling: An insufficiently dense k-point grid for the electronic Brillouin zone integration.
• Finite Supercell Size: In the finite-displacement method, the use of a finite supercell truncates the long-range interatomic interactions. If the supercell is too small, the force constants do not decay to zero, violating the ASR [4] [5].
• Symmetry Detection Failures: If the atomic positions are close to, but not exactly at, high-symmetry positions, the code may fail to correctly identify the full symmetry of the system. This prevents the code from using symmetry to automatically enforce the ASR for symmetry-related atoms [3] [1].
• System-Specific Challenges:
  • Molecules: For isolated molecules, the rotational modes (which should also have zero frequency) often exhibit non-zero, and sometimes imaginary, frequencies due to finite supercell effects in periodic boundary condition calculations or imperfect convergence [3] [1].
  • Polar Materials: In ionic materials, long-range dipole-dipole interactions require special treatment (e.g., using the Ewald sum). Incorrect or missing values for the Born effective charges and the dielectric tensor will lead to a violation of the "polar" ASR [4] [2].

Impact on Calculated Phonon Frequencies

The most direct consequence of ASR violation is the appearance of non-zero frequencies for the acoustic modes at the Γ-point [3]. Instead of three modes with zero frequency, calculations yield small but non-zero values. This error can propagate, affecting the entire phonon dispersion and density of states, particularly at low frequencies. In severe cases, it can even lead to imaginary frequencies (often reported as negative frequencies) that falsely indicate a structural instability when the material is actually stable [3]. These inaccuracies render subsequent calculations of thermal, optical, and transport properties unreliable.

The following workflow diagram illustrates the relationship between the causes of ASR violation and the corresponding corrective strategies.

Methodologies for Enforcing the ASR

Explicit Correction of the Dynamical Matrix

The most straightforward method is to explicitly impose the ASR as a constraint on the calculated force constants or the dynamical matrix after the calculation. This is often a post-processing step. For example, in Quantum ESPRESSO, this can be activated by setting the input flag asr = .true. in the ph.x input file, which enforces the acoustic sum rule on the dynamical matrix, typically with excellent results if the initial violation is small [1]. The CASTEP code offers a similar option with the phonon_sum_rule : TRUE keyword [6]. The specific mathematical technique used can vary, with common approaches including the projection of the force constant matrix onto a subspace that satisfies the sum rules or the use of Lagrange multipliers to minimize the changes to the IFCs while enforcing the constraints [2].

Ensuring Proper Convergence and Symmetry

A more fundamental approach is to address the root causes of the violation:

• Convergence Checks: Systematically converging key parameters is essential. The table below summarizes critical parameters and their roles.

• Symmetry-Aware Setup: To avoid symmetry-related errors, it is recommended to define the crystal structure using the correct internal Bravais lattice index (ibrav) rather than a generic setting (ibrav=0). Using standardized Wyckoff positions for atomic coordinates also helps the code correctly identify the system's symmetry, which in turn aids in automatically satisfying the ASR [3].

• Treatment of Polar Materials: For polar materials, the long-range Coulomb interactions must be handled correctly. This involves:

  • Performing a DFPT calculation in the unit cell with LEPSILON=.TRUE. (or equivalent) to compute the Born effective charges (( Z^* )) and the static dielectric tensor (( \epsilon_\infty )) [4].
  • Providing these tensors as input (PHON_BORN_CHARGES, PHON_DIELECTRIC) to the phonon calculation and setting LPHON_POLAR=.TRUE. (in VASP) to account for the long-range dipole-dipole interactions via Ewald summation [4]. This enforces the "polar" ASR [2].

Table 1: Key Computational Parameters for Converged Phonon Calculations

Parameter	Description	Function in ASR Enforcement
tr2_ph	Self-consistency threshold for phonon calculations [3] [1]	Tighter convergence reduces numerical noise in force constants.
conv_thr	Convergence threshold for the electronic ground state [3]	A well-converged starting density is crucial for accurate forces.
ecutwfc / cut_off_energy	Plane-wave kinetic energy cutoff [6] [1]	Higher cutoff improves basis set completeness, aiding symmetry.
kpoints_mp_spacing	K-point grid density for electrons [6]	Denser grids improve Brillouin zone sampling.
Supercell Size	Size of supercell for finite-displacement [4] [5]	Larger supercells capture long-range interactions fully.

Advanced and Machine Learning Approaches

Emerging methods leverage machine learning (ML) to address the computational cost of high-quality phonon calculations. One strategy uses Machine Learning Interatomic Potentials (MLIPs), such as MACE or M3GNet, which are trained on a diverse dataset of structures and forces from DFT [5]. These models learn the underlying potential energy surface and can be designed to inherently respect physical symmetries, including translational and rotational invariance. Once trained, they can predict forces for large supercells at a fraction of the computational cost of DFT, allowing for the computation of well-converged force constants that satisfy the ASR [5]. Another, more direct, strategy involves training graph neural networks like ALIGNN or VGNN to predict the entire phonon dispersion or density of states directly from the crystal structure, bypassing the calculation of the force constants altogether [5].

The Scientist's Toolkit: Essential Computational Reagents

Table 2: Key Software and Computational Tools for ASR-Compliant Phonon Research

Tool / Resource	Type	Primary Function in ASR Context
Quantum ESPRESSO	Software Suite	PWscf for SCF; ph.x for DFPT phonons with asr input flag [3] [1].
VASP	Software Suite	Finite-displacement (IBRION=5,6) or DFPT (IBRION=7,8) phonons; requires LPHON_POLAR for polar materials [4].
CASTEP	Software Suite	DFPT and finite-displacement phonons with phonon_sum_rule keyword [6].
Materials Project	Database	Contains pre-calculated phonon data for thousands of materials, though convergence should be verified [5].
Phonopy	Software Tool	A widely used tool for post-processing finite-displacement force sets, which includes ASR correction options.
MACE Potential	ML Interatomic Potential	A state-of-the-art ML model for accurately predicting forces and energies for phonon calculations [5].
ALIGNN	ML Model	Graph neural network for direct prediction of phonon density of states, bypassing force constants [5].

The Acoustic Sum Rule is not merely a technicality in computational materials science; it is the numerical embodiment of the fundamental physical principle of translational invariance. Its violation, stemming from the practical necessities of discrete and finite numerical computation, introduces unphysical artifacts that can compromise the integrity of phonon spectra and all derived properties. A comprehensive strategy for robust phonon research must therefore include: 1) a thorough understanding of the Born-Huang invariance conditions, 2) systematic convergence of key computational parameters, 3) careful attention to crystal symmetry and the treatment of long-range forces, and 4) the judicious application of ASR correction techniques. As the field progresses, machine learning approaches promise to alleviate the computational burden associated with achieving naturally ASR-compliant force constants, thereby enabling the high-throughput discovery and accurate characterization of next-generation materials.

The Critical Role of ASR in Ensuring Zero Frequency for Acoustic Modes at Gamma Point

The Acoustic Sum Rule (ASR) represents a fundamental cornerstone in the first-principles calculation of lattice dynamical properties within the harmonic approximation. This invariance principle dictates that the phonon frequencies of the three acoustic modes must approach zero as the wavevector q approaches the Brillouin zone center, known as the Gamma point (Γ). Physically, this phenomenon corresponds to the trivial translation of the entire crystal lattice, which should not experience any restoring force. The ASR is mathematically expressed as a constraint on the interatomic force constants (IFCs) in reciprocal space, requiring that the sum of the force constants between any atom and all other atoms in the crystal must vanish [7]. In practical computational terms, this translates to the condition that the dynamical matrix must possess three eigenvalues approaching zero at the Γ-point, corresponding to the acoustic modes.

Violations of the ASR, however, frequently occur in ab initio phonon calculations and present significant challenges for accurate material property prediction. These violations manifest as imaginary frequencies (often reported as negative frequencies in computational outputs) near the Γ-point, which lack physical meaning and indicate numerical instabilities in the calculation rather than genuine lattice instabilities. Such violations primarily stem from insufficient k-point sampling in the Brillouin zone, inadequate plane-wave energy cutoffs, or numerical precision limitations in the computation of force constants [7]. These inaccuracies propagate errors in derived thermodynamic properties and can severely impact the predictive capability of high-throughput computational materials screening. Consequently, enforcement of the ASR through post-processing corrections or improved computational parameters remains an active area of research in computational materials science, particularly within large-scale projects like the Materials Project dedicated to phonon property calculation [7] [8].

Theoretical Foundation of ASR and Phonon Formalism

The theoretical framework for phonon calculations in periodic solids derives from density functional perturbation theory (DFPT), which provides an efficient approach for computing second-order derivatives of the total energy. For a generic point q in the Brillouin zone, the phonon frequencies ωq,m and eigenvectors Um(qκ′β) are solutions to the generalized eigenvalue problem [7] [8]:

Here, κ labels atoms in the unit cell, α and β are Cartesian coordinates, C̃κα,κ′β(q) represents the interatomic force constants in reciprocal space, and Mκ denotes the atomic mass. The force constants C̃ are obtained through Fourier interpolation of values calculated on a regular grid of q-points using DFPT. Within this formalism, the ASR emerges as a direct consequence of translational invariance, requiring that the energy remains unchanged under infinitesimal rigid translations of the entire crystal lattice. This invariance imposes the following constraint on the IFCs at the Γ-point [7]:

This mathematical statement of the ASR guarantees that the dynamical matrix at q=0 has three eigenvalues equal to zero, corresponding to the acoustic modes. For polar materials, additional considerations involving long-range dipole-dipole interactions necessitate the inclusion of Born effective charges (BECs) and the dielectric tensor to correctly describe the splitting between longitudinal optical (LO) and transverse optical (TO) modes. The BEC tensor Zκ,βα, defined as the second derivative of energy with respect to atomic displacements and electric field, must simultaneously satisfy a charge neutrality sum rule (CNSR): ∑κZκ,βα = 0 [7]. Violations of either the ASR or CNSR serve as important indicators of numerical convergence issues in DFPT calculations, potentially arising from insufficient k-point grids, inadequate plane-wave cutoffs, or other computational parameters.

Table 1: Fundamental Sum Rules in DFPT Phonon Calculations

Sum Rule	Physical Principle	Mathematical Expression	Consequence
Acoustic Sum Rule (ASR)	Translational invariance	∑κC̃κα,κ′β(q=0) = 0	Three acoustic modes at Γ point have ω = 0
Charge Neutrality Sum Rule (CNSR)	Conservation of charge	∑κZ*κ,βα = 0	Correct description of LO-TO splitting in polar materials

Methodologies for ASR Implementation and Enforcement

The practical implementation of ASR in high-throughput computational frameworks involves specific methodologies and convergence criteria to ensure numerical accuracy. The Materials Project, which employs the ABINIT software package for density functional theory (DFT) and DFPT calculations, utilizes the PBEsol generalized gradient approximation exchange-correlation functional, which has demonstrated improved accuracy for phonon frequencies compared to experimental data [7] [8]. Calculations typically employ norm-conserving pseudopotentials from the PseudoDojo table (version 0.3), with plane-wave cutoffs selected based on the hardest element in each compound [7].

Brillouin zone sampling represents a critical parameter for ASR compliance. The Materials Project utilizes equivalent k-point and q-point grids that respect crystal symmetries with a density of approximately 1500 points per reciprocal atom, with the q-point grid always Γ-centered [7]. Structural relaxation precedes phonon calculations with strict convergence criteria: forces on atoms must be below 10⁻⁶ Ha/Bohr and stresses below 10⁻⁴ Ha/Bohr³ [7]. These stringent thresholds help minimize initial numerical errors that could propagate to phonon calculations.

Despite optimized parameters, raw DFPT outputs frequently exhibit minor ASR violations. To address this, the ASR and CNSR are explicitly imposed during the Fourier interpolation process of the force constants [7]. This enforcement procedure corrects for numerical inaccuracies that manifest as small imaginary frequencies near the Γ-point. However, the degree of deviation from these sum rules before imposition provides valuable diagnostics for assessing calculation quality. Significant deviations typically indicate insufficient convergence with respect to key computational parameters such as plane-wave cutoff energy or k-point/q-point grid density. Researchers therefore monitor these deviations as quality metrics, with large values triggering recalculation with refined parameters [7].

Quantitative Diagnostics for ASR Violation and Convergence

Systematic evaluation of ASR compliance requires specific quantitative diagnostics that serve as convergence metrics. High-throughput phonon databases, such as those generated by the Materials Project, implement several numerical flags to identify potentially problematic calculations [7]. These indicators help researchers quickly assess the reliability of phonon data without manual inspection of every phonon dispersion curve.

The primary indicator involves direct examination of the acoustic phonon frequencies at the Γ-point after the DFPT calculation but before the explicit imposition of the ASR. A well-converged calculation should exhibit acoustic frequencies very close to zero. The Materials Project sets a warning flag if the largest acoustic mode at Γ exceeds 30 cm⁻¹ when the ASR is not explicitly imposed [7]. Additionally, the presence of small negative frequencies in the region 0 < |q| < 0.05 (in fractional coordinates) along high-symmetry directions often indicates poor k-point or q-point sampling rather than a genuine physical instability [7].

For the charge neutrality sum rule, the relevant metric is the maximum value of the residual Born effective charge tensor: maxα,β |∑κZ*κ,βα|. Calculations where this value exceeds 0.2 are flagged for potential CNSR violation [7]. These quantitative thresholds, summarized in Table 2, provide an automated screening mechanism for identifying calculations requiring further convergence refinement. Monitoring these parameters across large datasets reveals that ASR violations are often correlated with specific material classes, such those with soft phonon modes or complex electronic structures, guiding more targeted convergence testing.

Table 2: Quantitative Warning Flags for Phonon Calculation Convergence

Warning Indicator	Threshold Value	Physical Interpretation	Recommended Action
Acoustic Mode Frequency at Γ	> 30 cm⁻¹	Severe ASR violation before imposition	Increase k-point/q-point grid density
CNSR Violation	max⎪∑Z*⎪ > 0.2	Significant breaking of charge neutrality	Check BEC convergence; increase cutoff
Imaginary Frequencies near Γ	Negative values for 0<⎪q⎪<0.05	Likely numerical instability, not physical	Improve q-point sampling; verify structure

Experimental Validation and Case Studies

Experimental validation of computationally predicted phonon spectra provides crucial verification of ASR enforcement methodologies. For instance, a comprehensive study on germanium monosulfide (GeS) employed DFPT-calculated phonon dispersions to interpret resonant Raman scattering data [9]. Bulk GeS crystallizes in an orthorhombic structure with 24 phonon branches, and the DFPT calculations successfully reproduced the measured phonon mode energies and symmetries when proper convergence criteria and sum rules were applied [9]. The experimental observation of Raman-forbidden phonon modes under resonant conditions further confirmed the accuracy of the computed phonon eigenvectors and the theoretical framework underlying the calculations.

Similarly, investigations of MgB₂ phonon modes using Raman and infrared spectroscopy highlighted discrepancies between experimental observations and simple symmetry predictions [10]. While point group analysis of the P6/mmm symmetry predicted only one Raman and two IR-active peaks, experimental spectra revealed additional peaks [10]. These discrepancies were resolved through DFPT calculations that considered super-lattice structures, demonstrating the importance of accurate force constant calculations and proper ASR treatment for interpreting complex experimental data. The case of MgB₂ particularly underscores how ASR-compliant phonon calculations can elucidate phenomena beyond harmonic approximations, including phonon anomalies and electron-phonon coupling relevant to superconductivity.

Table 3: Comparison of Selected Calculated and Experimental Phonon Modes in GeS

Mode Notation	Symmetry	Experimental ω (cm⁻¹)	Calculated ω (cm⁻¹)	Activity
Ag₂	Γ	116	111.8	Raman Active
B1g₂	Γ	218	215	Raman Active
B1u₂	Γ	231	232	IR Active (Raman Forbidden)
B3u₂	Γ	282	285	IR Active (Raman Forbidden)

Successful implementation of ASR-compliant phonon calculations relies on specific software tools, computational resources, and methodological protocols. The following toolkit outlines essential components for researchers undertaking these investigations.

Table 4: Research Reagent Solutions for ASR-Compliant Phonon Calculations

Tool Category	Specific Solution	Function/Role in ASR Compliance
DFPT Software	ABINIT [7] [8]	Primary engine for calculating second-order energy derivatives and force constants
Pseudopotential Library	PseudoDojo v0.3 [7]	Provides optimized norm-conserving pseudopotentials for accurate force calculations
Exchange-Correlation Functional	PBEsol [7] [8]	Generalized gradient approximation optimized for solids; improves phonon frequency accuracy
Sum Rule Enforcement	Internal post-processing codes	Explicitly imposes ASR and CNSR during Fourier interpolation of force constants
Phonon Visualization	Materials Project phonon website [8]	Validates acoustic modes approach zero at Γ point through band structure visualization
Convergence Diagnostics	ASR/CNSR breaking metrics [7]	Quantitative flags (e.g., Γ-point frequencies >30 cm⁻¹) identify unconverged calculations

The enforcement of the Acoustic Sum Rule represents an indispensable component of robust phonon frequency calculations within density functional perturbation theory. As demonstrated through high-throughput frameworks like the Materials Project, explicit imposition of the ASR during the Fourier interpolation of force constants corrects for numerical inaccuracies arising from finite k-point sampling and other computational limitations, ensuring that acoustic modes properly approach zero at the Γ-point. The quantitative diagnostics developed for ASR and CNSR violation provide crucial quality metrics for automated screening of computational results across extensive materials databases. Future research directions will likely focus on developing more efficient algorithms for ASR enforcement in large-scale systems, improving convergence monitoring in high-throughput workflows, and extending these principles to more complex physical scenarios such as anharmonic lattice dynamics and electron-phonon coupling calculations. By maintaining rigorous adherence to these fundamental sum rules, computational materials scientists can ensure the reliability of predicted vibrational properties essential for understanding thermal conductivity, phase stability, and superconducting behavior in novel materials.

The acoustic sum rule (ASR) serves as a fundamental cornerstone in computational lattice dynamics, enforcing the physical requirement that the net force on any atom must vanish under uniform displacements. Violations of the ASR, frequently precipitated by the discrete nature of Fast Fourier Transform (FFT) grids and inadequate convergence of numerical parameters, introduce significant inaccuracies in predicted phonon spectra and related material properties. This technical guide provides an in-depth examination of how FFT grid discreteness undermines ASR compliance, presents detailed methodologies for diagnosing and mitigating these issues, and situates these solutions within the broader context of ensuring predictive accuracy in first-principles phonon calculations. The discussions and protocols are designed to equip researchers with the practical tools necessary to achieve robust and reproducible results in the study of lattice dynamics.

In the framework of lattice dynamics, the dynamical matrix is central to determining the vibrational properties of a crystal. The ASR is a direct consequence of the invariance of the potential energy when the entire crystal is translated. Mathematically, for a crystal with (N) atoms in the unit cell, the ASR states that the dynamical matrix (\mathbf{D}(\mathbf{q})) at wave vector (\mathbf{q}) must satisfy: [ \sum{\kappa'} \mathbf{D}{\kappa\kappa'}(\mathbf{q} \rightarrow 0) = 0 ] for each atom (\kappa), where (\kappa') runs over all atoms in the unit cell. This ensures that the three acoustic phonon branches exhibit a linear dispersion relation as (\mathbf{q}) approaches the Brillouin zone center ((\Gamma)-point), correctly yielding zero frequency for these modes at (\Gamma).

Theoretical and Practical Significance: Adherence to the ASR is not merely a mathematical formality; it is physically crucial for accurately capturing long-wavelength phenomena. ASR violations manifest as unphysical imaginary frequencies in acoustic modes, which subsequently corrupt the prediction of thermodynamic properties (e.g., specific heat, free energy) and transport properties (e.g., lattice thermal conductivity) [11] [12]. Within the broader thesis of phonon calculation research, enforcing the ASR represents a fundamental step in bridging the gap between ab initio force constants and physically meaningful lattice dynamical models.

The Critical Role of FFT Grids in Force Constant Calculations

First-principles phonon calculations, often employing the density functional perturbation theory (DFPT) or the finite displacement method, fundamentally rely on the computation of interatomic force constants (IFCs). The IFCs describe the force on a given atom when another atom is displaced from its equilibrium position. The precision of these IFCs is paramount for constructing an accurate dynamical matrix.

The computation of IFCs in planewave-based density functional theory (DFT) codes necessitates the use of FFT grids to switch between real and reciprocal space representations of the electron density and potential. The discrete sampling of the FFT grid can lead to two primary sources of error:

• Aliasing Effects: Inadequate grid sampling fails to capture the high-frequency components of the electron density, leading to aliasing and inaccurate forces.
• Broken Translational Invariance: The ASR is a manifestation of translational invariance in the crystal. The discrete and finite nature of the FFT grid can break this exact invariance, as the force on an atom due to a uniform displacement of all atoms may not integrate precisely to zero numerically [13].

Consequently, the choice of the FFT grid is not just a technical detail but a potential source of systematic error that can trigger ASR violations. The following sections detail how this discreteness leads to violations and outline protocols for its mitigation.

FFT Grid Discreteness as a Primary Source of ASR Violation

Mechanism of Violation

The core of the problem lies in the numerical integration of the Hellmann-Feynman forces. A uniform displacement of the crystal should produce no net force on any atom. However, when the electron density response is computed on a discrete FFT grid, the numerical integration is no longer exact. The error originates from the fact that the FFT grid samples the continuous charge density at discrete points, and the quality of this sampling is determined by the grid's density. A coarse grid fails to resolve subtle changes in the charge density induced by atomic displacements, leading to forces that do not perfectly satisfy the sum rule. This error then propagates into the calculated IFCs.

Quantitative Impact on Phonon Frequencies

The most immediate and visible consequence of ASR violation is the appearance of unphysical imaginary frequencies in the acoustic phonon branches near the (\Gamma)-point. The table below summarizes the typical manifestations and their implications for material properties.

Table 1: Manifestations and Implications of ASR Violations from FFT Grid Discreteness

Aspect Affected	Manifestation of Error	Impact on Material Property Prediction
Acoustic Phonons	Imaginary frequencies at the (\Gamma)-point; non-linear dispersion near (\Gamma)	Invalidates predictions of sound velocity, thermodynamic stability.
Phonon Dispersion	Incorrect splittings and gaps throughout the Brillouin zone	Reduces accuracy for electron-phonon coupling and phase transition analysis.
Phonon DOS	Spurious low-energy peaks from imaginary modes	Corrupts calculation of specific heat, entropy, and free energy.
Thermal Conductivity	Unphysical scattering rates and relaxation times	Leads to significant over- or under-estimation of (\kappa_l) [12].

Experimental Protocols for Diagnosing and Mitigating Grid-Related Issues

Protocol 1: FFT Grid Convergence Testing

This protocol provides a systematic methodology for determining the sufficient FFT grid density to minimize ASR violations.

Objective: To establish the minimum FFT grid parameter that ensures ASR compliance and yields converged phonon frequencies.

Materials and Computational Setup:

• Software: A DFT code with phonon calculation capabilities (e.g., Quantum ESPRESSO, VASP, ABINIT).
• System: A target material (e.g., Silicon with a 2-atom unit cell).
• Initial Parameters: A converged plane-wave kinetic energy cutoff.

Methodology:

• Baseline Calculation: Perform a ground-state self-consistent calculation with a "coarse" FFT grid, typically defined by the code's default settings or a low multiplier of the plane-wave cutoff.
• Force Constant Calculation: Compute the dynamical matrix or the IFCs using the finite displacement method or DFPT.
• ASR Check: Explicitly calculate the sum (\sum{\kappa'} \mathbf{D}{\kappa\kappa'}(\mathbf{q} \rightarrow 0)) for a small, finite (\mathbf{q})-vector. A non-zero result indicates ASR violation. Visually inspect the phonon band structure for imaginary acoustic modes at (\Gamma).
• Iterative Refinement: Systematically increase the FFT grid density (e.g., by increasing the multiplier for the plane-wave cutoff) and repeat steps 1-3.
• Convergence Criterion: The calculation is considered converged when the frequency of the highest acoustic branch at a point very close to (\Gamma) (e.g., (|\mathbf{q}| = 0.001) in fractional coordinates) changes by less than 0.1 THz upon a further increase of the grid density, and the ASR sum is numerically zero.

The workflow for this protocol can be visualized as follows:

Protocol 2: ASR Correction and Post-Processing

When computational constraints prevent the use of a fully converged FFT grid, or for high-throughput screening where manual convergence is impractical, post-processing corrections can be applied to the IFCs.

Objective: To enforce the ASR on the computed IFCs after their initial calculation, thereby correcting for residual numerical noise.

Materials:

• A set of computed IFCs, (\Phi_{\alpha\beta}(\kappa, \kappa')), where (\alpha, \beta) are Cartesian indices.

Methodology:

• Extraction: Extract the real-space IFCs from the DFT phonon calculation.
• Correction Algorithm: Apply a constraint to the IFCs. A common and simple approach is the "projection method": [ \Phi{\alpha\beta}^{\text{corrected}}(\kappa, \kappa') = \Phi{\alpha\beta}(\kappa, \kappa') - \frac{\delta{\kappa\kappa'}}{N} \sum{\kappa''} \Phi_{\alpha\beta}(\kappa, \kappa'') ] This formula subtracts the net force on atom (\kappa) uniformly across all its connections, ensuring the sum over (\kappa') is zero.
• Validation: Reconstruct the dynamical matrix using the corrected IFCs and verify that the acoustic phonon frequencies at the (\Gamma)-point are now zero.

Table 2: Research Reagent Solutions for Phonon and ASR Studies

Item / Software	Primary Function	Relevance to ASR/FFT Issues
Phonopy [11]	Open-source package for phonon calculations.	Widely used; includes functionality to impose ASR on force constants post-calculation.
ShengBTE [12]	Calculates lattice thermal conductivity from BTE.	Relies on accurate, ASR-compliant phonon dispersions to compute scattering rates and (\kappa_l).
Phono3py [12]	Computes three-phonon scattering.	Input force constants must satisfy ASR for physically meaningful scattering phase space.
Quantum ESPRESSO	Integrated suite for DFT & DFPT.	Allows explicit control over FFT grid density via ecutrho parameter to mitigate aliasing.
High-Performance Computing (HPC) Cluster	Provides computational resources.	Essential for running convergence tests with dense FFT grids, which are computationally demanding [13].

The relationship between the core calculation, common issues, and mitigation strategies is summarized in the following workflow:

FFT grid discreteness presents a fundamental and pervasive challenge in first-principles phonon calculations, directly leading to violations of the acoustic sum rule and a cascade of errors in predicted material properties. This guide has outlined the mechanistic link between numerical sampling and broken translational invariance, providing two robust, detailed experimental protocols for addressing the issue: systematic FFT grid convergence testing and post-processing ASR correction. As the field of computational materials science progresses towards high-throughput screening and the study of complex systems with large unit cells [12], a rigorous and standardized approach to managing these numerical parameters is not just beneficial—it is essential for ensuring the reliability and predictive power of lattice dynamics simulations. Future research will likely focus on developing more efficient algorithms that intrinsically preserve sum rules and on integrating these validation protocols seamlessly into automated computational workflows.

Distinguishing Physical Instabilities from Computational Artifacts in Phonon Spectra

Phonon spectra, which describe the vibrational properties of crystals, are fundamental for understanding material stability, thermal conductivity, and phase transitions. In computational materials science, a central challenge arises when phonon calculations predict imaginary frequencies (negative frequencies squared)—these results can either indicate genuine physical instabilities or stem from numerical artifacts. Distinguishing between these possibilities is critical for accurate materials characterization. This guide examines the primary sources of computational artifacts in phonon spectra, with particular focus on Acoustic Sum Rule (ASR) violations and improper handling of long-range interactions, providing researchers with diagnostic methodologies and mitigation protocols.

Fundamental Concepts: Physical vs. Computational Origins

Physical Instabilities

Genuine imaginary frequencies in phonon spectra signify structural instabilities where the crystal structure is not in its ground state. These often correspond to:

• Structural phase transitions to lower-symmetry phases
• Charge density waves or other electronically driven distortions
• Soft modes in ferroelectric or magnetic materials
• Surface reconstructions or low-dimensional systems with inherent instability

Computational Artifacts

Spurious imaginary frequencies arise from approximations and numerical errors in calculations:

• ASR violations from incomplete force constant convergence
• Insufficient treatment of long-range dipolar interactions in polar materials
• Inadequate k-point sampling or supercell size
• Numerical precision limits in finite-difference methods

Acoustic Sum Rule Violations

The ASR requires the sum of interatomic force constants (IFCs) between an atom and all other atoms to vanish, ensuring zero energy cost for uniform translations. In practice, finite supercell sizes and numerical approximations during IFC calculation violate this rule, particularly affecting acoustic modes near the Brillouin zone center (Γ-point).

Improper Treatment of Long-Range Interactions

In polar materials, the coupling between atomic vibrations and electric fields creates long-range dipole-dipole interactions that decay slowly in real space. Neglecting these interactions causes significant errors, particularly for longitudinal optical (LO) modes, leading to unphysical imaginary frequencies and incorrect phonon dispersion [4] [14].

Table: Characteristic Differences Between Physical and Artificial Imaginary Frequencies

Feature	Physical Instability	Computational Artifact
Frequency Magnitude	Typically large (>10 cm⁻¹)	Often small (<10 cm⁻¹), especially near Γ-point
q-point Dependence	Appears at specific high-symmetry points	Dominant near Γ-point, disappears with better convergence
Supercell Size	Persists with increasing supercell size	Diminishes or disappears with larger supercells
LO-TO Splitting	Follows physical symmetry patterns	Shows unphysical discontinuities without proper corrections [4]
Symmetry Analysis	Corresponds to specific irreducible representations	Violates point group symmetry expectations

Diagnostic Methodologies and Protocols

ASR Correction Techniques

Several post-processing methods can correct ASR violations in calculated IFCs:

• Simple ASR: Applying constraints to force the IFCs to satisfy sum rules
• Diagonal ASR: Enforcing sum rules separately for each Cartesian component
• Eigenvalue-shift approaches: Modifying the dynamical matrix eigenvalues
• Iterative symmetry-preserving methods: Projecting IFCs onto physically allowed space

Table: Computational Parameters for Reliable Phonon Calculations

Parameter	Recommended Value	Effect on Phonon Spectrum
Supercell Size	≥ 4×4×4 for complex materials	Determines maximum IFC range; insufficient size causes ASR violations
k-point Mesh	Dense enough for electronic convergence	Affects force accuracy; impacts Born effective charges and dielectric tensor
Energy Cutoff (ENCUT)	High enough for force convergence	Inadequate values cause numerical noise in forces
Displacement Step Size	0.01-0.015 Å (finite differences)	Too large or small values introduce numerical errors
Q-point Path	Dense sampling for dispersion	Sparse sampling misses key features

Long-Range Correction Protocol for Polar Materials

For polar materials, the long-range dipole-dipole contributions must be explicitly included. The implementation protocol requires [4]:

Diagram: Workflow for Accurate Phonon Calculations in Polar Materials

• Calculate dielectric properties in the primitive cell:

  • Perform DFPT calculation with IBRION = 7 or 8
  • Set LEPSILON = True or LCALCEPS = True to compute Born effective charges (Z*) and static dielectric tensor (ε∞) [4]
  • Ensure high k-point density and energy cutoff convergence

• Extract the dielectric tensors from output:

  • Located in OUTCAR, vasprun.xml, or vaspout.h5 files
  • Born effective charges: results/born_charges/born_charges in vaspout.h5
  • Dielectric tensor: results/dielectric/dielectric_dft in vaspout.h5

• Perform supercell phonon calculation with proper settings:

  • Set LPHON_POLAR = True in INCAR
  • Provide PHON_BORN_CHARGES and PHON_DIELECTRIC tensors
  • Optionally set reciprocal space cutoff PHON_G_CUTOFF for Ewald summation [4]

Convergence Testing Protocol

Systematic convergence tests are essential for distinguishing artifacts:

• Supercell size convergence:

  • Calculate phonons with increasing supercell sizes (2×2×2, 3×3×3, 4×4×4)
  • Monitor the magnitude of imaginary frequencies near Γ-point
  • Genuine instabilities persist while artifacts diminish

• Numerical parameters verification:

  • Test k-point mesh density effect on forces
  • Verify energy cutoff convergence for stress and forces
  • Check displacement step size in finite-difference approaches [15]

The Scientist's Toolkit: Research Reagent Solutions

Table: Essential Computational Tools for Phonon Analysis

Tool/Code	Primary Function	Application Context
VASP [4]	DFPT and finite-difference phonons	First-principles phonon dispersion and DOS
phonopy [11]	Post-processing force constants	Phonon dispersion, DOS, and thermal properties
CRYSTAL [16]	Hybrid functional phonon calculations	Accurate Raman and IR spectra prediction
Multibinit [14]	Effective atomic potential simulations	Large-scale phonon calculations with anharmonicity
AMS [15]	Molecular vibrational analysis	IR frequencies and normal modes for molecules

Case Studies and Experimental Validation

MgO: LO-TO Splitting Artifacts

In magnesium oxide (MgO), a polar material with rock-salt structure, neglecting long-range corrections produces significant artifacts [4]:

• Without corrections: Unphysical oscillations and imaginary frequencies in optical branches
• With LPHON_POLAR = True: Smooth dispersion curves with proper LO-TO splitting
• The magnitude of LO-TO splitting is comparable to phonon frequencies themselves

BaTiO₃: Long-Range Interactions in Ferroelectrics

In barium titanate (BaTiO₃), a prototypical ferroelectric, dipole-dipole interactions are crucial [14]:

• Models neglecting long-range interactions can capture the ferroelectric instability
• However, they produce artifacts in specific phonon branches and spurious minima
• The analytical dipole-dipole correction is essential for accurate phonon dispersion

Numerical Hessian Calculation in Molecular Systems

For molecular systems, the numerical calculation of the Hessian matrix requires careful attention [15]:

• Full Hessian calculation requires 6N single-point calculations for N atoms
• Mobile Block Hessian (MBH) method allows partial frequency calculations for large systems
• Symmetric displacements utilize molecular symmetry to reduce computational cost
• Mode rescanning techniques identify spurious imaginary frequencies

Diagram: Diagnostic Pathway for Imaginary Frequency Analysis

Advanced Considerations: Electron-Phonon Coupling and Berry Curvature

Recent research reveals more complex relationships between electronic structure and phonon properties. In Dirac materials like BaMnSb₂ with time-reversal symmetry, phonons can inherit Berry curvature through coupling with topologically nontrivial electrons [17]. This effect:

• Generates phonon Berry curvature proportional to electronic valley Chern number
• Creates helical phonons with locked angular momentum and wave vector direction
• Emerges even when "bare" phonons have vanishing Berry curvature
• Could influence thermal transport properties through thermal Hall effect

Such phenomena illustrate that what may appear as artifacts could sometimes represent novel physics, emphasizing the need for careful diagnostic protocols.

Distinguishing physical instabilities from computational artifacts in phonon spectra requires systematic methodology and rigorous validation. Key principles include: (1) implementing appropriate ASR corrections for IFCs, (2) including long-range dipole-dipole interactions for polar materials, (3) performing comprehensive convergence tests, and (4) utilizing symmetry analysis and comparison with experimental data where available. Following the protocols outlined in this guide will enable researchers to accurately identify genuine material instabilities while avoiding misinterpretation of numerical artifacts, leading to more reliable computational materials characterization.

Computational Methods for Phonon Calculations: DFPT, Finite Displacement, and Machine Learning Potentials

The calculation of phonons—the quantized vibrational modes of a crystal lattice—is fundamental to understanding numerous material properties, from thermal conductivity and phase stability to spectroscopic signatures. Within the framework of Density Functional Theory (DFT), two primary methodologies have emerged for computing these vibrational spectra: the Finite Displacement Method (FDM) and Density-Functional Perturbation Theory (DFPT). The choice between them significantly impacts the computational workflow, the accessible physical properties, and the treatment of fundamental constraints like the Acoustic Sum Rule (ASR).

This guide provides an in-depth technical comparison of these two approaches, with a particular focus on the context of ASR violation and its implications for phonon frequency calculation research. It is structured to equip computational researchers and scientists with the knowledge to select the most appropriate method for their specific investigations.

Theoretical Foundations and Methodologies

Finite Displacement Method (FDM): A Direct Numerical Approach

The Finite Displacement Method, also known as the frozen-phonon approach, is a direct and conceptually straightforward technique. It operates by numerically evaluating the second derivatives of the total energy through atomic displacements.

• Core Principle: The method approximates the interatomic force constants, which form the dynamical matrix, by calculating the forces induced by small, finite displacements of atoms from their equilibrium positions. The fundamental equation for the force constant is: [ C{\text{mai}}^{\text{nbj}} = \frac{F{-}^{\text{nbj}} - F{+}^{\text{nbj}}}{2 \cdot \delta} ] Here, ( F{-}^{\text{nbj}} ) and ( F_{+}^{\text{nbj}} ) represent the force on atom ( n ) in direction ( j ) when atom ( m ) is displaced by a small distance ( \pm \delta ) in direction ( i ) [18].
• Implementation Variants: In practice, codes like VASP offer two main tags for FDM:
  • IBRION = 5: Displaces all atoms in all three Cartesian directions, leading to 3N displacements (where N is the number of atoms) [19].
  • IBRION = 6: Uses crystal symmetry to identify and compute only the symmetry-inequivalent displacements, significantly reducing the computational number of calculations required [19].
• Key Parameters: The step size POTIM is a critical parameter, with a default of 0.015 Å often providing a reasonable compromise between numerical accuracy and stability [19]. The NFREE tag controls the number of displacements per ion and direction, with NFREE=2 (central differences) being the standard.

Density-Functional Perturbation Theory (DFPT): An Analytical Linear Response Approach

DFPT is a more sophisticated method that leverages perturbation theory to compute the linear response of the electron density to a perturbation, such as an atomic displacement, directly from the self-consistent ground state.

• Core Principle: DFPT solves the Sternheimer equation to obtain the first-order change in the wavefunctions without the need for explicit atomic displacements or supercells for different phonon wavevectors [20] [21]. This allows for the direct calculation of the dynamical matrix at any wave vector q within the primitive cell.
• Mathematical Formulation: The phonon frequencies ( \omega{\mathbf{q},m} ) and eigenvectors ( Um(\mathbf{q} \kappa'\beta) ) are obtained by solving the eigenvalue problem: [ \sum{\kappa'\beta} \tilde{C}{\kappa\alpha, \kappa'\beta}(\mathbf{q}) Um(\mathbf{q} \kappa'\beta) = M\kappa \omega{\mathbf{q},m}^2 Um(\mathbf{q} \kappa\alpha) ] where ( \tilde{C} ) are the interatomic force constants in reciprocal space [7].
• Implementation: In VASP, DFPT is activated with IBRION = 7 (all displacements) or IBRION = 8 (symmetry-reduced displacements) [20]. In codes like ABINIT and CASTEP, specific keywords like task: PHONON or optdriver: 3 are used to initiate the DFPT workflow [6] [22].

Direct Comparison: DFPT vs. FDM

The choice between DFPT and FDM involves trade-offs across computational cost, functionality, and ease of use. The table below summarizes the key differentiating factors.

Table 1: A comparative overview of DFPT and the Finite Displacement Method

Feature	Density-Functional Perturbation Theory (DFPT)	Finite Displacement Method (FDM)
Theoretical Basis	Analytical linear response [23]	Numerical finite differences [18]
Computational Cost for Full Dispersion	Lower; force constants for any q-point obtained from primitive cell calculation [23]	Higher; requires a supercell large enough to fit the phonon wavelength, with 3N force calculations [23] [24]
Handling of Polar Materials	Superior. Born effective charges ((Z^*)) and dielectric tensor ((\epsilon_\infty)) are computed self-consistently for accurate LO-TO splitting [6] [7]	Possible but complex. Requires separate calculation of (Z^*) and (\epsilon_\infty) (e.g., via Berry phase) to apply a posteriori corrections [6]
Functional & Potential Compatibility	Limited. Typically restricted to LDA/GGA; not available for HF/hybrids, meta-GGA, or ultrasoft pseudopotentials in some codes [6]	Broad. Compatible with any XC functional, DFT+U, and pseudopotential type, including ultrasoft [6] [19]
Key Advantage	High efficiency for phonon dispersions; direct access to dielectric properties	Generality and simplicity; easier to implement and debug
Main Disadvantage	Limited compatibility with advanced Hamiltonians	High computational cost for low-symmetry or large-unit-cell materials

Workflow and Protocol Comparison

The practical application of these methods involves distinct workflows, from initial setup to final analysis. The diagram below illustrates the key steps in each protocol.

Diagram 1: A comparison of the DFPT and Finite Displacement Method workflows.

Detailed Experimental and Computational Protocols

Finite Displacement Method Protocol (e.g., using VASP)

• Ground-State Relaxation: Perform a full relaxation of the atomic positions and cell vectors (ISIF=3) using IBRION=2 to ensure the structure is at a true energy minimum. Forces should be converged to a tight threshold (e.g., < 1 meV/Å) [24].
• Symmetry Enforcement: Check the final structure (CONTCAR) and, if necessary, enforce the expected crystallographic symmetry by rounding lattice parameters and atomic coordinates. A subsequent relaxation with fixed cell (ISIF=2) and symmetry enabled (ISYM=2) is recommended [24].
• Supercell Construction: Build a supercell commensurate with the phonon wavevectors of interest. Convergence with respect to supercell size must be tested.
• Force Calculations:
  • Prepare the INCAR file with IBRION = 5 or 6, NSW = 1, NFREE = 2, and POTIM = 0.015.
  • Set PREC = Accurate and use a k-point mesh that provides a sampling density equivalent to that of the primitive cell [19].
  • Execute VASP. The code will generate the displaced configurations, run the calculations, and compute the force constants.

DFPT Protocol (e.g., using ABINIT)

• Highly Converged Ground State: Run a ground-state calculation to generate the Kohn-Sham wavefunctions. This step requires an extremely well-converged calculation with a stringent SCF stopping criterion (tolvrs ~1.0d-18 in ABINIT) [21].
• DFPT Calculation:
  • In the input file, set the driver option for DFPT (e.g., optdriver 3 in ABINIT for phonons) [22].
  • Specify the q-points where phonons are to be calculated using a block like PHONON_KPOINT_LIST [6].
  • To compute dielectric properties, set the appropriate flags (e.g., lrfopt in ABINIT).
• Generation of the Derivative Database (DDB): The calculation produces a DDB file containing the second-order derivatives. This file is then used by a post-processing tool (e.g., anaddb in ABINIT) to compute phonon dispersions, densities of states, and thermodynamic properties [22].

The Acoustic Sum Rule (ASR): Violations and Enforcement

Understanding the ASR and Its Violation

The Acoustic Sum Rule is a direct consequence of the translational invariance of the total energy. It states that the sum of the force constants for all atoms interacting with a given atom must be zero in each direction, ensuring that the frequencies of the three acoustic modes at the Γ-point are exactly zero [7].

[ \sum{\kappa} \tilde{C}{\kappa\alpha, \kappa'\beta}(\mathbf{q}=0) = 0 ]

In practice, numerical approximations lead to ASR violations. In FDM, these can arise from using a supercell that is too small, insufficient k-point sampling, or inadequate force convergence. In DFPT, the violation can stem from the numerical evaluation of the exchange-correlation functional or the use of a finite plane-wave basis set [20] [7]. The magnitude of the ASR breaking, often quantified by the maximum frequency of the acoustic modes at Γ, is a key indicator of the calculation's numerical quality. Flags have been proposed to identify problematic calculations, such as when the largest acoustic mode at Γ exceeds 30 cm⁻¹ before ASR imposition [7].

Enforcement Techniques and Their Implications

Both methodologies require explicit enforcement of the ASR to obtain physically meaningful results. The following diagram illustrates the pathways to ASR violation and the corresponding correction strategies.

Diagram 2: Pathways to Acoustic Sum Rule (ASR) violation and enforcement strategies.

• Enforcement in High-Throughput Studies: In large-scale computational studies, the imposed ASR and Charge Neutrality Sum Rule (CNSR) are critical for data quality. The Materials Project phonon database, for instance, imposes these sum rules during the Fourier interpolation process. They also use the deviation from the exact condition as a metric to flag potentially unconverged calculations [7].
• Practical Recommendations:
  • For FDM, ensure forces are highly converged and use the largest feasible supercell. Post-processing tools like phonopy often include options to apply the ASR to the calculated force constants.
  • For DFPT, increasing the plane-wave cutoff (ENCUT in VASP) can reduce the intrinsic violation. Post-processing codes like anaddb in ABINIT provide options to enforce the ASR and CNSR, which is particularly recommended for PAW calculations due to the large number of real-space integrals [22].

The Scientist's Toolkit: Essential Computational Reagents

Table 2: Key software tools and computational parameters for phonon calculations.

Tool / Parameter	Category	Function / Purpose	Example Usage
VASP [20] [19]	Software Package	A widely used DFT code implementing both FDM and DFPT for phonons.	IBRION=5/6 (FDM), IBRION=7/8 (DFPT)
ABINIT [21] [22]	Software Package	A DFT code with a strong focus on DFPT for various response properties.	optdriver 3, rfphon 1
CASTEP [6]	Software Package	A DFT code featuring DFPT for phonons and spectroscopic intensities.	task : PHONON
Phonopy [18] [24]	Post-Processor	A tool to post-process force constants from FDM to obtain phonon DOS, dispersion, and thermodynamics.	phonopy --fc vasprun.xml
POTIM [19]	FDM Parameter	Step size for atomic displacements in FDM. Critical for numerical accuracy.	POTIM = 0.015 (default in VASP)
ENCUT [19]	Basis Set	Plane-wave kinetic energy cutoff. Must be increased for stress convergence in elastic tensor calculations.	ENCUT = 1.3 * default (recommended for FDM)
Born Effective Charges (Z*)	DFPT Property	Tensor describing polarization change per atomic displacement; crucial for LO-TO splitting.	Computed with LEPSILON = .TRUE. in VASP DFPT [24]
Dielectric Tensor (ε∞)	DFPT Property	High-frequency dielectric constant from electronic polarization.	Computed with LEPSILON = .TRUE. in VASP DFPT [24]

The decision between DFPT and the Finite Displacement Method is not a matter of which is universally superior, but which is optimal for a specific research context. DFPT is the more efficient and elegant choice for high-throughput studies of phonon dispersions and dielectric properties within its domain of applicability (LDA/GGA with norm-conserving pseudopotentials). Conversely, the Finite Displacement Method offers unparalleled generality, making it indispensable for systems involving advanced functionals like DFT+U, hybrids, or specific pseudopotential types.

A critical consideration for both methods is the rigorous management of numerical accuracy, with Acoustic Sum Rule violation serving as a key metric. Researchers must systematically converge computational parameters and actively enforce the ASR to ensure the reliability of their predicted phonon frequencies. As computational materials science progresses towards more complex systems and higher levels of accuracy, the interplay between these two powerful methodologies will continue to be a cornerstone of predictive lattice dynamics.

Phonons, the quantized lattice vibrations in crystalline materials, are fundamental determinants of a wide range of material properties including thermal conductivity, mechanical behavior, and phase stability [5]. In computational materials science, predicting these phonon properties relies heavily on calculating interatomic force constants (IFCs), which represent the second-order derivatives of the potential energy with respect to atomic displacements. These IFCs form the foundation of the dynamical matrix from which phonon frequencies and eigenvectors are obtained.

A persistent challenge in these calculations is the violation of the acoustic sum rule (ASR), which stems from the requirement that the total energy should remain invariant under rigid translations of the crystal. When this rule is broken—often due to numerical noise, approximations in calculations, or discrete sampling grids—it manifests as non-zero frequencies for acoustic phonons at the Brillouin zone center (Γ-point) [3] [25] [26]. Such violations can lead to unphysical imaginary frequencies in phonon spectra and compromised predictive accuracy for thermodynamic properties.

This technical guide provides a comprehensive workflow for phonon calculations, from supercell creation to force constant determination, while explicitly addressing ASR violations throughout the process. By integrating traditional first-principles approaches with emerging machine learning methodologies, we present robust protocols for obtaining reliable phonon properties that obey the fundamental symmetries of crystalline materials.

Fundamental Concepts and Challenges

Force Constants and Symmetry Constraints

In the harmonic approximation, interatomic force constants (Φ) are defined as:

$$ \Phi{\alpha\beta}(\ell\kappa,\ell'\kappa') = \frac{\partial^2 E}{\partial u{\alpha}(\ell\kappa)\partial u_{\beta}(\ell'\kappa')} $$

where E is the potential energy, and $u{\alpha}(\ell\kappa)$ represents the displacement of atom κ in unit cell ℓ along the Cartesian direction α. These force constants must satisfy translational invariance, which requires that the sum of all force constants between one atom and all others must vanish: $\sum{\ell'\kappa'} \Phi_{\alpha\beta}(\ell\kappa,\ell'\kappa') = 0$ for all κ and α,β [27]. This fundamental constraint gives rise to the acoustic sum rule.

Acoustic Sum Rule Violations: Origins and Implications

ASR violations occur in practical calculations due to several factors:

• Numerical discretization: The finite size of Fast Fourier Transform (FFT) grids in density functional theory (DFT) calculations represents the most irreducible source of ASR violation [3]
• Insufficient convergence: Inadequate k-point sampling or force convergence criteria during geometry optimization [3]
• Supercell approximations: Finite supercell sizes that truncate long-range interactions [28]
• Approximations in force calculations: Numerical noise in finite-difference approaches or approximations in density-functional perturbation theory [25]

These violations have direct computational consequences: they result in non-zero frequencies for acoustic phonons at the Γ-point and can introduce unphysical imaginary frequencies in phonon dispersions, potentially misinterpreted as structural instabilities [3] [25].

Computational Workflow for Phonon Calculations

The following diagram illustrates the comprehensive workflow for phonon calculations, integrating both traditional DFT and machine learning approaches:

Figure 1: Comprehensive workflow for phonon calculations, showing both traditional finite-displacement and machine learning-accelerated approaches with integrated ASR handling.

Supercell Creation and Displacement Generation

Supercell Construction

The first critical step involves constructing an appropriately sized supercell from the primitive unit cell. The supercell must be sufficiently large to ensure that force constants vanish at the boundary, typically requiring a convergence study with respect to supercell size [4]. For materials with long-range interactions, particularly polar materials, larger supercells may be necessary to capture the correct physics.

Displacement Generation Strategies

Two primary approaches exist for generating atomic displacements:

Traditional Finite-Displacement Method: This approach involves creating multiple supercell configurations where individual atoms are displaced slightly (typically 0.01 Å) from their equilibrium positions [4]. For a supercell containing N atoms, this method requires 6N single-point DFT calculations (3 displacements × 2 directions per atom), making it computationally demanding for large systems [29] [5].

Random Displacement Approach: Emerging as a more efficient alternative, this method generates a smaller subset of supercells where all atoms are simultaneously randomly displaced within a sphere of radius 0.01-0.05 Å [29] [5]. This strategy significantly reduces the number of required supercell calculations while still providing sufficient information for force constant extraction, particularly when combined with machine learning interatomic potentials (MLIPs).

Table 1: Comparison of Displacement Generation Strategies

Method	Displacement Magnitude	Number of Supercells	Computational Cost	Best For
Finite-Displacement	0.01 Å	6N (N = number of atoms)	High	Small systems, high-precision
Random Displacement	0.01-0.05 Å	40-100 total	Low to Moderate	Large systems, high-throughput
MLIP-Accelerated	0.01-0.05 Å	40-100 total	Very Low (after training)	Large-scale screening

Force Calculation Methods

Density Functional Theory

Traditional DFT calculations provide the reference forces for each displaced supercell configuration. These calculations must be performed with high numerical accuracy, including:

• Appropriate k-point sampling for the supercell
• Sufficiently high plane-wave energy cutoff
• Tight force convergence criteria (typically 1-10 meV/Å) [29]

The resulting forces and displacements form the dataset for subsequent force constant extraction.

Machine Learning Interatomic Potentials

MLIPs offer a powerful alternative for rapid force evaluation once trained. The "one defect, one potential" strategy has proven particularly effective for defect phonon calculations, where a specialized MLIP is trained on a limited set of perturbed supercells specific to a defect system [29]. Modern architectures like Neural Equivariant Interatomic Potentials (NequIP) [29] and Multi-Atomic Cluster Expansion (MACE) [5] demonstrate remarkable data efficiency, achieving accurate phonon predictions with as few as 40 training structures regardless of supercell size.

Force Constant Extraction

Force constants are extracted from the displacement-force datasets by solving the system of equations relating atomic displacements to the resulting forces. In the harmonic approximation, this relationship is linear:

$$ F{i\alpha} = -\sum{j\beta} \Phi{\alpha\beta}(i,j) u{j\beta} $$

where $F{i\alpha}$ is the force on atom i in direction α, and $u{j\beta}$ is the displacement of atom j in direction β.

Specialized software packages like symfc implement symmetry-constrained extraction of force constants, explicitly enforcing the crystal symmetry relationships throughout the extraction process [27]. This approach yields force constants that inherently satisfy the ASR and other symmetry requirements, reducing the need for post-processing corrections.

Addressing Acoustic Sum Rule Violations

ASR Enforcement Methods

When force constants violate the ASR, several correction schemes can be applied:

On-site Correction Methods: These approaches modify the self-interaction force constants (on-site IFCs) to enforce the sum rule. The ANADDB code (part of the ABINIT package) provides two variants [26]:

• asr = 1 (asymmetric): Fully imposes ASR but may break the symmetry of on-site IFCs
• asr = 2 (symmetric): Maintains symmetry of on-site IFCs but may not exactly enforce ASR

Iterative Schemes: VASP implements an iterative correction method controlled by the IFC_ASR tag, where the number of iterations determines the degree of ASR enforcement [25]. The default setting (IFC_ASR = 1) typically provides sufficient correction for most systems.

Projection Methods: The symfc software employs projector-based estimation that inherently satisfies symmetry constraints, including ASR, during the initial force constant extraction rather than as a post-processing step [27].

Table 2: ASR Enforcement Methods in Popular Computational Codes

Software	Control Variable	Methods Available	Key Features
VASP	IFC_ASR	Iterative correction	Number of iterations adjustable (IFC_ASR > 0)
ABINIT/ANADDB	asr	On-site IFC modification	Choice between symmetric (asr=2) and asymmetric (asr=1)
symfc	Built-in	Projection-based estimation	Force constants inherently satisfy ASR during extraction
QuantumATK	ASR checkbox	Gamma-point correction	Simple enforcement for acoustic modes at Γ-point

Special Considerations for Polar Materials

Polar materials require additional handling of long-range dipole-dipole interactions, which contribute to the characteristic LO-TO splitting of optical phonon modes at the Γ-point [4]. Proper treatment requires:

• Calculation of Born effective charges and the dielectric tensor through DFPT (e.g., using LEPSILON = .TRUE. in VASP)
• Inclusion of these tensors in the phonon calculation via parameters like PHON_BORN_CHARGES and PHON_DIELECTRIC in VASP [4]
• Setting LPHON_POLAR = .TRUE. to activate the non-analytical term correction

Neglecting these corrections leads to unphysical phonon dispersions, particularly evident as oscillations and improper behavior near the Γ-point.

Machine Learning Approaches for Accelerated Phonon Calculations

Machine Learning Force Fields

MLIPs have emerged as powerful tools for accelerating phonon calculations while maintaining DFT-level accuracy. The key advantage lies in their ability to rapidly evaluate forces once trained, bypassing the need for explicit DFT calculations for each displacement configuration [29] [5] [30].

The training workflow for defect-specific MLIPs involves [29]:

• Generating reference data using DFT calculations on randomly perturbed supercells
• Training neural network potentials (e.g., using NequIP or Allegro) on energies and forces
• Validating against held-out DFT calculations
• Using the trained MLIP for force evaluations in phonon calculations

This approach has demonstrated remarkable efficiency, achieving accurate phonon predictions with as few as 40 training structures, regardless of supercell size [29].

Universal Potentials for High-Throughput Screening

For high-throughput materials screening, universal MLIPs trained on diverse datasets offer compelling advantages. Recent work has demonstrated training on 15,670 supercell structures spanning 77 elements, achieving a mean absolute error of 0.18 THz for vibrational frequencies across a held-out set of 384 materials [5]. These universal potentials enable rapid phonon calculations without system-specific training, though with potentially reduced accuracy compared to specialized potentials.

Practical Implementation Protocols

Step-by-Step Workflow for VASP Calculations

• Supercell Creation: Generate a sufficiently large supercell using PHON_SUPERCELL or external tools
• Displacement Generation: Use IBRION = 5 or 6 for finite-differences, or IBRION = 7 or 8 for DFPT [4]
• Force Constant Calculation: Run calculations with LPHON = .TRUE. to compute and store force constants
• ASR Enforcement: Set IFC_ASR = 1 (or higher) to apply iterative ASR correction [25]
• Phonon Property Calculation:
  • For dispersions: Set LPHON_DISPERSION = .TRUE. and provide a QPOINTS file with the desired path
  • For DOS: Set PHON_DOS > 0 and provide a uniform q-mesh in QPOINTS [4]
• Polar Materials: Include Born charges and dielectric tensor with LPHON_POLAR = .TRUE.

Step-by-Step Workflow for phono3py Calculations

For third-order force constants and thermal conductivity calculations [31]:

• Displacement Dataset Generation: Create supercells with displacements using phono3py
• Force Set Calculation: Compute forces on displaced structures using DFT or MLIP
• Force Constant Extraction: Calculate second- and third-order force constants from displacements and forces
• Thermal Conductivity Calculation: Compute lattice thermal conductivity using the obtained force constants

Validation and Troubleshooting

Validating Phonon Calculations:

• Check for small or zero acoustic frequencies at Γ-point
• Verify absence of unphysical imaginary frequencies (except for genuine instabilities)
• Compare with experimental data where available
• Ensure convergence with respect to supercell size and k-point sampling

Troubleshooting ASR Issues:

• For persistent ASR violations, increase the IFC_ASR iteration count in VASP [25]
• In ANADDB, try both asr = 1 and asr = 2 to determine which provides better acoustic modes [26]
• For severe violations, check structural symmetry and ensure proper geometry optimization
• Consider using symmetry-constrained methods like symfc for more robust enforcement [27]

Table 3: Research Reagent Solutions for Phonon Calculations

Tool/Software	Primary Function	Key Features	Applicability
VASP	DFT & Phonon Calculations	IFC_ASR for ASR enforcement, DFPT support	Comprehensive materials suite
ABINIT/ANADDB	DFT & Phonon Analysis	Multiple ASR schemes (asr=1,2), dielectric properties	First-principles phonons
Phono3py	Phonon-Phonon Interactions	Third-order force constants, thermal conductivity	Lattice thermal conductivity
symfc	Symmetry-Adapted Force Constants	Projector-based estimation, inherent ASR satisfaction	Symmetry-constrained IFCs
MLIPs (NequIP, MACE)	Machine Learning Force Fields	Data-efficient training, rapid force evaluation	Large systems, high-throughput

The practical workflow for phonon calculations, from supercell creation to force constant determination, requires careful attention to both computational methodology and fundamental physical constraints. The integration of traditional DFT approaches with emerging machine learning methods offers powerful pathways to accelerate these calculations while maintaining accuracy. Throughout this process, proper handling of acoustic sum rule violations remains essential for obtaining physically meaningful phonon properties.

By employing the protocols and strategies outlined in this guide—including appropriate displacement generation, symmetry-aware force constant extraction, and systematic ASR enforcement—researchers can reliably compute phonon dispersions, densities of states, and thermal properties across a broad range of material systems. The continued development of machine learning interatomic potentials and specialized phonon codes promises to further enhance the efficiency and scope of these calculations, enabling high-throughput screening and investigation of complex defective systems that were previously computationally prohibitive.

The Acoustic Sum Rule (ASR) is a fundamental physical constraint in phonon calculations, requiring that the sum of the interatomic force constants between any atom and all other atoms in the crystal must be zero. This condition ensures the correct physical behavior that a rigid translation of the entire crystal produces no net force on any atom, guaranteeing that the three acoustic phonon branches exhibit zero frequency at the Brillouin zone center (Γ point). In practical density functional theory (DFT) simulations, this rule is often violated due to the use of finite numerical parameters, such as a truncated supercell size that is too small for force constants to decay to zero, insufficient k-point sampling for Brillouin zone integration, or numerical precision limits in force calculations. Such violations manifest as imaginary frequencies in acoustic phonon branches near the Γ point, rendering the phonon spectrum unphysical and complicating the calculation of derived thermodynamic properties.

The implications of ASR violations are particularly pronounced in the study of polar materials like NaCl, where improper treatment of long-range interactions can lead to significant inaccuracies in the optical phonon branches. In such materials, the violation is twofold: the standard ASR related to translational invariance, and the need to account for long-range dipole-dipole interactions that cause LO-TO splitting. The failure to properly handle these effects can severely impact the predictive accuracy of lattice dynamics simulations, especially for properties like thermal conductivity, phase stability, and spectroscopic response. This technical guide examines the implementation of ASR corrections within three major computational codes—Quantum ESPRESSO, CASTEP, and Phonopy—framed within broader thesis research on phonon frequency calculations.

Theoretical Foundation of ASR and Long-Range Corrections

The formal definition of the ASR stems from translational invariance in the crystal potential. For the force constant matrix Φ, this invariance requires that the force on atom α when all atoms are displaced uniformly must be zero, leading to the mathematical expression: ∑{β} Φ{αβ} = 0, where the sum runs over all atoms β in the crystal, including α itself. In practical supercell calculations, this sum is truncated at the supercell boundary, introducing an error that grows as the supercell size decreases.

For polar materials, an additional complication arises from the long-range nature of the Coulomb interaction, which requires a separate treatment beyond the standard ASR. The dynamical matrix in this case splits into a short-range part and a non-analytical part that depends on the direction in which the wavevector q approaches zero: D(q) = Dshort-range(q) + DNA(q) The non-analytical term DNA(q) contains contributions from the Born effective charges (Z*) and the macroscopic dielectric tensor (ε∞), which account for the polarization response to atomic displacements and electric fields, respectively. This term is responsible for the LO-TO splitting observed in polar materials, where the longitudinal optical (LO) and transverse optical (TO) phonon frequencies differ at the Γ point.

Table 1: Key Physical Quantities for Proper Treatment of Long-Range Interactions in Polar Materials

Quantity	Symbol	Physical Meaning	Calculation Method
Born Effective Charges	Z*	Response of polarization to atomic displacement	DFPT or finite differences
Dielectric Tensor	ε_∞	Response to macroscopic electric field	DFPT with ε_∞ = .true.
Non-analytical Correction	D_NA(q)	Direction-dependent long-range interaction	Ewald sum with Z* and ε_∞
LO-TO Splitting	ωLO - ωTO	Frequency difference from polarization	Included via NAC

The theoretical framework for handling these effects was established in the classic work by Wang et al. on mixed-space approaches for polar materials, where the reciprocal-space treatment of the long-range interaction complements the real-space calculation of short-range force constants. Implementation of these corrections requires careful extraction of the dielectric properties and proper integration with the phonon calculation workflow, as detailed in the following sections.

ASR Implementation in Quantum ESPRESSO

Workflow and Code Integration

Quantum ESPRESSO provides two primary approaches for phonon calculations: the finite displacement method (through Phonopy) and the density functional perturbation theory (DFPT) method (through the ph.x code). For ASR correction in the finite displacement approach, the workflow involves several interconnected calculations:

• Initial Self-Consistent Field (SCF) Calculation: A highly converged electronic structure calculation using pw.x provides the ground state charge density.
• DFPT Calculation for Dielectric Properties: The ph.x code calculates the Born effective charges and dielectric tensor with specific input parameters tr2_ph = 1.0d-14, epsil = .true..
• Supercell Force Calculations: Multiple pw.x calculations on displaced supercells provide the force constants.
• Phonon Analysis with Non-analytical Correction: Phonopy processes the forces and incorporates the long-range correction using the BORN file.

Table 2: Key Input Parameters for DFPT Dielectric Property Calculation in Quantum ESPRESSO

Parameter	Value	Purpose	Location in Input
calculation	'scf'	Ground state calculation	&control
tprnfor	.true.	Print forces	&control
tstress	.true.	Print stress	&control
ecutwfc	70.0	Plane-wave cutoff	&system
k-point grid	8x8x8	Brillouin zone sampling	K_POINTS
tr2_ph	1.0d-14	Convergence threshold	&inputph
epsil	.true.	Calculate dielectric properties	&inputph

Creating the BORN File

The BORN file serves as the crucial link between the DFPT dielectric properties calculation and the phonon analysis. It contains the macroscopic dielectric tensor and the Born effective charge tensors for each atom in the primitive cell. After successfully running the ph.x calculation, the dielectric tensor and Born effective charges must be extracted from the output file (NaCl.ph.out in the NaCl example) and formatted into the BORN file following this specific structure [32]:

The first line is a header, the second line contains the dielectric tensor in 3×3 row-major format, and subsequent lines contain the Born effective charge tensors for each atom in the primitive cell. For convenience, the phonopy-qe-born utility can automate this process: phonopy-qe-born NaCl.in | tee BORN [33].

DFPT for Phonons and ASR

Quantum ESPRESSO's native DFPT capability (ph.x with ldisp = .true.) provides an alternative approach that inherently handles the ASR. When using DFPT, the ASR can be imposed during the post-processing step using q2r.x by setting the appropriate flags. However, caution is needed when combining DFPT with Phonopy, as the non-analytical correction treatment differs between the two approaches [33]. Specifically, the dielectric constant and Born effective charges embedded in the Gamma point dynamical matrix files must be manually removed when using q2r.x to generate force constants for Phonopy.

ASR Implementation in CASTEP

Workflow and Code Integration

CASTEP, a proprietary DFT code, integrates with Phonopy through a similar finite-displacement approach. The workflow for CASTEP shares conceptual similarities with Quantum ESPRESSO but requires specific file preparations and command options:

• Generate Supercells: Use phonopy -d --dim="2 2 2" --castep -c unitcell.cell to create displaced supercells [34].
• Prepare Calculation Directories: A shell script (make_displ_dirs.sh) organizes the supercell calculations by combining the displacement information with CASTEP parameter files.
• Run CASTEP Calculations: Execute CASTEP for each displacement directory to compute forces.
• Collect Forces and Apply NAC: Use phonopy --castep -f to gather forces into FORCE_SETS, then run phonopy with the BORN file for NAC.

For CASTEP, the unitcell.cell input file must contain all necessary structural information, including initial spin values for magnetic materials and LDA+U parameters where applicable [34]. The parameter file (supercell.param) contains standard CASTEP settings for convergence criteria, energy cutoffs, and exchange-correlation functionals.

Special Considerations for Magnetic Materials

CASTEP calculations for magnetic materials like chromium require additional parameters in the unitcell.cell file, specifically SPINS in the positions_frac block [34]. The tail file (tail.cell) supplements these with k-point mesh specifications, U parameters, symmetry operations, and pseudopotential assignments. This separation allows for efficient management of parameters common to all displacement calculations versus those specific to the supercell configuration.

ASR Implementation in Phonopy

Role as an Integrator Platform

Phonopy serves as a versatile integrator platform that can leverage force calculations from various DFT codes, including Quantum ESPRESSO, CASTEP, and VASP, while providing a unified interface for phonon analysis and ASR correction. Its primary function concerning ASR is the post-processing of real-space force constants to impose the sum rule and incorporate non-analytical corrections for polar materials.

The general Phonopy workflow consists of:

• Generating displaced supercells with phonopy -d --dim="N N N" -c input_file
• Calculating forces for each supercell using the preferred DFT code
• Collecting forces into the FORCE_SETS file with phonopy -f
• Performing phonon analysis with optional NAC via phonopy --nac

ASR Correction Methods in Phonopy

Phonopy implements several approaches to handle ASR violations:

• Simple ASR Symmetrization: The --asr option allows applying simple sum rules to the force constants [18].
• Non-analytical Correction (NAC): For polar materials, the --nac option activates the mixed-space approach for LO-TO splitting when a BORN file is present [32] [33].
• Force Constants Symmetrization: Setting FC_SYMMETRY = .TRUE. enforces the symmetry of force constants, which indirectly helps with ASR compliance [33].

Table 3: Comparison of ASR Handling Across Different Computational Workflows

Code/Workflow	ASR Method	NAC Implementation	Key Inputs for NAC
QE+Phonopy	--nac option with Phonopy	Through BORN file	Z*, ε∞ from ph.x
CASTEP+Phonopy	--nac option with Phonopy	Through BORN file	Z*, ε∞ from separate CASTEP calculation
VASP	LPHON_POLAR = True	Direct in DFPT	PHON_BORN_CHARGES, PHON_DIELECTRIC [4]
ASE	acoustic=True in ph.read()	born=True in band_structure()	Born charges and dielectric tensor [18]

Comparative Analysis and Best Practices

Code-Specific Considerations

Each code presents unique advantages and considerations for ASR correction:

• Quantum ESPRESSO offers the most comprehensive and documented approach for NAC, with robust DFPT capabilities for dielectric property calculation. The integration with Phonopy is well-established but requires careful file management between different components (pw.x, ph.x, and phonopy).
• CASTEP provides a streamlined workflow through shell scripts but requires more manual preparation of parameter files. The documentation for obtaining Born charges and dielectric constants specifically for NAC is less detailed than for Quantum ESPRESSO.
• VASP, while not the primary focus of this guide, implements NAC through the LPHON_POLAR tag, requiring explicit input of Born charges and dielectric tensors in the INCAR file [4].

Computational Parameters and Convergence

The accuracy of ASR correction depends critically on several computational parameters:

• Supercell Size: Larger supercells reduce ASR violations by allowing force constants to decay naturally. A 2×2×2 supercell represents a minimum, with 4×4×4 often needed for convergence [4].
• k-point Sampling: Denser k-point meshes improve the accuracy of Born effective charges and dielectric constants. For the initial SCF calculation, an 8×8×8 mesh or denser is recommended [32].
• Energy Cutoff: Higher plane-wave cutoffs (ecutwfc) improve force accuracy but increase computational cost.
• DFPT Convergence: Strict convergence thresholds (tr2_ph = 1.0d-14) are necessary for accurate dielectric properties [32].

Table 4: Recommended Computational Parameters for Accurate ASR Treatment

Parameter	Minimum Value	Recommended Value	Effect on ASR
Supercell Size	2×2×2	4×4×4	Reduces force constant truncation
k-point Mesh	4×4×4	8×8×8	Improves dielectric properties
SCF Convergence	1×10⁻⁸	1×10⁻¹²	More accurate forces
DFPT Convergence	1×10⁻¹²	1×10⁻¹⁴	Better Z* and ε∞

The Scientist's Toolkit: Essential Research Reagents

Table 5: Key Computational Tools for ASR Correction in Phonon Calculations

Tool/Component	Function	Implementation Examples
Born Effective Charges (Z*)	Quantify polarization response to atomic displacement	DFPT calculation in QE (ph.x), VASP (LEPSILON), or finite differences
Dielectric Tensor (ε∞)	Macroscopic response to electric field	DFPT with epsil=.true. in QE, LEPSILON in VASP
BORN File	Interface for NAC parameters between DFPT and phonon codes	Specific format with ε∞ tensor and Z* tensors
Force Constants	Second derivatives of energy with respect to atomic displacements	From finite differences or DFPT, symmetrized with ASR
Phonopy Post-processor	Applies ASR and NAC to force constants	--nac option with BORN file for polar materials
DFPT Codes	Calculate dielectric properties and phonons	QE ph.x, VASP with IBRION=7,8

Proper implementation of ASR correction, particularly through the non-analytical term correction for polar materials, is essential for obtaining physically meaningful phonon spectra in solid-state simulations. The integration between DFT codes like Quantum ESPRESSO and CASTEP with phonon post-processing tools like Phonopy has created robust workflows for handling these corrections, though each approach requires careful attention to parameter convergence and file management.

Future developments in this field are likely to focus on increased automation through materials informatics frameworks like atomate2, which provides standardized workflows for heterogeneous computational methods [35], and the integration of machine learning interatomic potentials for accelerated phonon calculations [36] [37]. These advances will make high-throughput phonon screening more accessible while maintaining the physical rigor necessary for accurate prediction of lattice dynamical properties.

Leveraging Machine Learning Potentials like MACE for High-Throughput Phonon Properties

The accurate calculation of phonon properties is fundamental to understanding material behavior, governing phenomena from thermal conductivity and mechanical stability to phase transitions. Traditional approaches based on Density Functional Theory (DFT) provide high accuracy but remain computationally prohibitive for large systems and high-throughput screening. The emergence of machine learning interatomic potentials (MLIPs), particularly the MACE (Multi-Atomic Cluster Expansion) architecture, now offers a pathway to accelerate these calculations while maintaining near-DFT accuracy. This technical guide examines the application of MACE models for phonon property prediction, with specific attention to persistent challenges including Acoustic Sum Rule (ASR) violation and imaginary frequency correction—critical considerations for determining dynamic stability within high-throughput computational workflows.

Theoretical Foundations: Phonons, ASR, and MLIPs

Phonon Calculations and the Acoustic Sum Rule

Phonons represent quantized lattice vibrations whose properties are derived from the dynamical matrix, constructed from interatomic force constants (IFCs). The Acoustic Sum Rule (ASR) embodies the physical principle of translational invariance, requiring that the energy of a system remains unchanged under rigid translation. Mathematically, this dictates that the sum of the IFCs between one atom and all others must be zero, guaranteeing three zero-frequency acoustic modes at the Γ-point [26].

Violations of the ASR occur in practical calculations due to numerical approximations, including the discreteness of the FFT grid or insufficient convergence parameters [3]. These violations manifest as non-zero frequencies for acoustic modes, which can be corrected by enforcing the ASR a posteriori on the dynamical matrix [3] [28]. Similarly, imaginary frequencies (ω² ≤ 0) arise from instabilities in the calculated potential energy surface, potentially indicating either a genuine structural instability or numerical issues from poor convergence or inadequate supercell size [3] [28].

Machine Learning Interatomic Potentials and MACE

Machine Learning Interatomic Potentials are trained to approximate the potential energy surface (PES) derived from quantum mechanical calculations. The MACE architecture employs multilayer equivariant message-passing neural networks built upon a tensor product representation, enabling it to capture complex, many-body atomic interactions with high fidelity [38].

The foundational MACE model, MACE-MP-0, provides broad coverage across the periodic table. However, its direct application to phonon properties in complex materials like metal-organic frameworks (MOFs) can yield spurious imaginary modes and inaccuracies in properties like thermal expansion [36] [38]. This limitation is addressed through specialized fine-tuning, resulting in derivatives like MACE-MP-MOF0, which is specifically optimized for accurate phonon calculations in MOFs by training on curated datasets of representative structures [36] [39].

Computational Workflow for Phonon Calculations with MACE

The standard workflow for computing phonon properties using MLIPs involves a structured sequence of steps to ensure reliability, particularly concerning structural stability and the mitigation of numerical artifacts. The overall process is depicted in the following workflow diagram.

Diagram: Workflow for phonon calculations using MACE models, highlighting the iterative step for addressing imaginary frequencies.

Key Methodological Steps

• Step 1: Full Cell Relaxation: The initial structure undergoes full, unconstrained relaxation using the fine-tuned MACE potential (e.g., MACE-MP-MOF0) until the maximum force component is below a stringent threshold (e.g., ≤ 10⁻⁶ eV/Å). This ensures the atomic configuration is at a local energy minimum, which is crucial for avoiding imaginary frequencies originating from residual forces [36].

• Step 2: Generation of Displaced Supercells: Two primary approaches are used to sample the PES for IFC calculation. The finite-displacement method creates supercells with small, systematic displacements of individual atoms [5]. An efficient alternative is the random displacement of all atoms within a supercell by 0.01–0.05 Å, which generates a diverse set of force data with fewer supercell configurations [5] [40].

• Step 3: Force Calculation and IFC Construction: The MACE model predicts atomic forces for all generated supercell configurations. These forces are used to fit the IFCs, which form the basis of the dynamical matrix [5].

• Step 4: ASR Enforcement and Post-Processing: The ASR is imposed on the dynamical matrix to correct for translational invariance violations. Common methods involve adjusting the on-site IFCs symmetrically (preferred for phonon band structure) or asymmetrically (strictly enforces ASR) [26].

Addressing Imaginary Frequencies

The workflow includes a critical check for imaginary frequencies. Their presence necessitates an iterative review of convergence parameters, supercell size, and structural relaxation. For MOFs, using a specialized model like MACE-MP-MOF0, which was fine-tuned to correct such artifacts, is often essential [36] [39].

Performance Benchmarking and Applications

Accuracy of Phonon Properties

Recent benchmarking studies have systematically evaluated the performance of universal MLIPs, including MACE, for predicting phonon properties across diverse materials.

Table 1: Benchmarking performance of universal machine learning potentials for phonon calculations.

Model	Training Data	Force MAE	Phonon Frequency MAE	Dynamical Stability Accuracy	Key Strengths
MACE	15,670 structures, 77 elements [5]	~45 meV/Å [38]	0.18 THz (full dispersion) [5]	86.2% [5]	Strong force accuracy, efficient message-passing architecture [38]
EquiformerV2	OMat24 dataset [41]	Low (comparable to MACE) [41]	-	-	Best overall for IFCs and lattice thermal conductivity [41]
MatterSim	Broad materials set [41]	Higher than MACE [41]	-	-	Intermediate IFC prediction via error cancellation [41]
CHGNet	Materials Project [41]	Low (comparable to MACE) [41]	-	-	Good force accuracy, but poor LTC prediction [41]

Notably, while MACE demonstrates excellent force prediction accuracy, its performance in deriving certain properties like lattice thermal conductivity (LTC) can show notable discrepancies compared to fine-tuned models like EquiformerV2, highlighting the complex relationship between force accuracy and final phonon properties [41].

Application to Metal-Organic Frameworks (MOFs)

The application of the fine-tuned MACE-MP-MOF0 model demonstrates the successful deployment of MLIPs for complex material systems. The model was trained on a curated dataset of 127 diverse MOFs, incorporating configurations from molecular dynamics, strained cells, and geometry optimization trajectories [36]. This approach enabled high-throughput prediction of phonon density of states, thermal expansion, and bulk moduli within the quasi-harmonic approximation, showing strong agreement with DFT and experimental data [36] [39]. The model successfully corrected imaginary phonon modes present in the foundational MACE-MP-0 predictions, establishing its utility for reliable high-throughput screening of MOF dynamical properties [36].

The Scientist's Toolkit: Essential Research Reagents

Table 2: Key computational tools and resources for high-throughput phonon calculations with machine learning potentials.

Tool/Resource	Type	Primary Function	Relevance to Phonon Calculations
MACE-MP-MOF0	Fine-tuned MLIP	Potential energy and force prediction for MOFs	Provides accurate IFCs; corrects imaginary modes in MOFs [36] [39]
ASE (Atomic Simulation Environment)	Python package	Structure manipulation, IO, and workflow control	Used for structure relaxation and analysis within the MACE phonon workflow [36]
Phonopy	Software package	Post-processing analysis for phonon calculations	Calculates phonon dispersion, DOS, and thermodynamic properties from IFCs
Curated MOF Dataset	Training data	127 representative MOF structures [36]	Enables fine-tuning for chemical accuracy in complex porous materials [36]
ASR Corrector (e.g., in Phonopy)	Algorithm	Enforces translational invariance on dynamical matrix	Corrects non-zero acoustic modes at Γ-point [3] [26]

Machine learning potentials, particularly those based on the MACE architecture, are fundamentally advancing high-throughput computational materials science by enabling rapid, accurate phonon calculations. The integration of robust workflows—incorporating rigorous structural relaxation, efficient sampling strategies, and post-calibration procedures like ASR enforcement—is critical for generating reliable data on dynamical stability and thermal properties. Future developments will likely focus on expanding the chemical diversity and system size covered by foundation models, improving the automated handling of imaginary modes, and tighter integration of these tools into multi-scale simulation platforms. The continued refinement of specialized, fine-tuned models like MACE-MP-MOF0 promises to unlock large-scale screening efforts for technologically relevant applications in thermal management, energy storage, and thermoelectrics.

Solving ASR Violations and Imaginary Frequencies: A Step-by-Step Troubleshooting Guide

Theoretical Foundation: The Acoustic Sum Rule and Its Violations

The Acoustic Sum Rule (ASR) is a fundamental physical principle arising from the translational invariance of a crystal lattice. It dictates that the sum of the force constants for any atom with all other atoms in the crystal must be zero. A direct consequence of this rule is that the phonon frequencies of the three acoustic modes at the Brillouin zone center (the Γ-point, where ( \mathbf{q} = 0 )) should be precisely zero. These modes correspond to uniform translations of the entire crystal. In computational practice, however, this ideal condition is frequently violated due to numerical approximations, leading to non-zero, often imaginary, frequencies for these acoustic modes. Such violations are not merely numerical artifacts but can significantly impact the accuracy of calculated thermodynamic properties and lattice dynamics.

The primary source of ASR violation in plane-wave Density Functional Theory (DFT) calculations, such as those performed with Quantum ESPRESSO, stems from the discreteness of the Fast Fourier Transform (FFT) grid used to represent charge density and potentials [42]. The total energy of the system becomes invariant only to translations by vectors that are multiples of the FFT grid spacing, breaking the continuous translational symmetry of the underlying physical system. This problem is often more pronounced in systems with significant vacuum regions, such as surfaces or molecules, and when using Generalized Gradient Approximation (GGA) functionals compared to Local Density Approximation (LDA), as GGA functionals have more strongly varying forms [42]. Other contributing factors include insufficient convergence of key computational parameters, an inadequately converged self-consistent field (SCF) calculation for the ground state, or a structure that has not been fully relaxed to its equilibrium geometry [3].

Table 1: Common Causes of Acoustic Sum Rule Violations in Phonon Calculations

Cause of Violation	Physical Reason	Typical Manifestation
Discrete FFT Grid [42]	Breakdown of continuous translational invariance; energy is only invariant to translations on the grid.	Non-zero acoustic modes at Γ, typically < 100 cm⁻¹.
Insufficient Convergence [3]	Incompletely calculated forces and force constants.	Unphysical imaginary frequencies & general inaccuracies.
Unrelaxed Structure [3]	Atomic forces in the ground state are not zero.	Fictitious negative frequencies for rotational/translational modes.
Use of GGA Functionals [42]	More complex functional form varies more strongly with position.	Generally slower convergence and more severe ASR violation than LDA.

Diagnostic Protocols: Assessing the Need for Manual ASR Correction

Before applying any corrective measures, it is crucial to diagnose the nature and severity of the ASR violation. A systematic diagnostic workflow, as illustrated in the diagram below, is essential for determining the correct course of action.

Diagram 1: Diagnostic workflow for ASR violations

The first and most critical diagnostic step is to quantify the magnitude of the error. A generally accepted rule of thumb is that frequencies for acoustic modes at the Γ-point below 10 cm⁻¹ are considered minor violations, often directly attributable to the FFT grid effect [42]. In such cases, manual application of the ASR is justified and expected to yield reliable results. However, if the frequencies are significantly higher (e.g., approaching 100 cm⁻¹) or are "loudly" imaginary (large negative values of ( \omega^2 )), this indicates a more profound problem [42]. In these scenarios, applying the ASR would be masking a fundamental issue with the calculation, such as a poor-quality ground-state wavefunction, an incorrect structure, or insufficient computational parameters. The ASR should be seen as a corrective measure for small, well-understood numerical errors, not a remedy for a failed calculation.

The next step is a thorough convergence and structure check. Key parameters that must be verified include the plane-wave kinetic energy cutoff (ecutwfc), the charge density cutoff (ecutrho), the k-point grid for sampling the Brillouin zone, and the convergence thresholds for both the SCF calculation (conv_thr) and the phonon calculation itself (tr2_ph) [42]. For metallic systems, convergence with respect to the k-point grid and smearing parameters can be particularly slow. Furthermore, one must ensure that the atomic structure is correctly defined and fully relaxed. "Mysterious symmetry-related errors" or wrong degeneracies often trace back to atomic positions that are close to, but not exactly at, high-symmetry positions [42] [3]. Using the correct Bravais lattice index (ibrav) and known Wyckoff positions in the initial self-consistent calculation, rather than ibrav=0, can prevent these issues.

Table 2: Key Parameters to Check Before Applying ASR

Parameter	Description	Impact on ASR
tr2_ph	Convergence threshold for phonon calculation [42].	A threshold that is too large is a common cause of gross ASR violations.
conv_thr	Convergence threshold for the SCF ground-state calculation [42].	Affects the quality of the ground-state charge density and forces.
ecutwfc & ecutrho	Plane-wave cutoffs for wavefunctions and charge density [42].	Too low a cutoff leads to inaccurate forces and force constants.
k-point grid	Sampling density of the Brillouin zone.	Insufficient sampling, especially in metals, causes inaccurate integrals.
Atomic Structure	Precision of atomic positions and lattice symmetry [42].	Atoms slightly off high-symmetry positions break symmetry and violate ASR.

The ultimate test for validating a phonon calculation with minor ASR violations is to use the dynmat.x post-processing tool available in Quantum ESPRESSO [42]. This program can read the dynamical matrix file and diagonalize it after explicitly imposing the Acoustic Sum Rule. If, after applying the ASR, the acoustic modes drop to a very small frequency (e.g., < 1 cm⁻¹) while all other optical modes remain virtually unchanged, the results of the original calculation can be trusted, and the ASR-corrected dynamical matrix should be used for further analysis.

Methodologies for Manual ASR Correction and Symmetrization

Once a diagnosis of a minor, numerically-induced ASR violation is confirmed, the researcher can proceed with manual correction. The most straightforward method within the Quantum ESPRESSO ecosystem is to use the dynmat.x tool. This program is specifically designed for manipulating dynamical matrices, including imposing the ASR and analyzing the resulting phonon modes. The typical workflow involves running the ph.x calculation to obtain the dynamical matrix at the Γ-point (or on a q-grid), and then processing the generated output file (often named matdyn or similar) with dynmat.x with the appropriate flag to enforce the sum rule.

The core operation of ASR correction involves a mathematical procedure on the force constant matrix. The force constant, ( C{mai}^{nbj} ), represents the force in direction *j* on atom *nb* when atom *ma* is displaced in direction *i*. The ASR requires that the sum of force constants over all atoms *nb* for a given displacement of atom *ma* must be zero: ( \sum{nb} C_{mai}^{nbj} = 0 ). In the finite displacement method, as implemented in codes like ASE, the force constants are calculated from central differences of the forces [43]. The ASR correction algorithm effectively redistributes the "missing" force, ensuring this sum rule is satisfied for every atom and direction. The acoustic method in the ASE Phonons class, for instance, performs this restoration explicitly [43].

For a comprehensive correction, especially when dealing with a full phonon dispersion, it is often necessary to impose the ASR on the entire set of interatomic force constants before performing the Fourier transform to obtain the dynamical matrix at an arbitrary q-vector. Some packages may do this automatically, while others require explicit user intervention. The subsequent symmetrization of the force constants based on the crystal's space group symmetry is also a critical step for ensuring the physical correctness of the results and improving computational efficiency, although it is noted that the small displacement method in ASE, for example, does not initially exploit space-group symmetries [43].

Diagram 2: ASR correction workflow for reliable results

The Scientist's Toolkit: Essential Research Reagents & Computational Tools

Table 3: Key Software and Computational "Reagents" for Phonon Calculations

Tool / Component	Function	Role in Diagnosing/Correcting ASR
Quantum ESPRESSO (pw.x, ph.x) [42]	Suite for DFT and DFPT ground-state and phonon calculations.	Generates the initial dynamical matrix and force constants where ASR violations may occur.
dynmat.x [42]	Post-processing tool for dynamical matrices.	The primary tool for manually imposing the ASoustic Sum Rule on the calculated dynamical matrix.
Symmetry Analysis Tools (e.g., ISOTROPY) [3]	Analyzes the symmetry of phonon modes and crystal structures.	Helps identify symmetry-related errors and validates the symmetry of calculated modes post-correction.
ASE (Atomic Simulation Environment) [43]	Python package for atomistic simulations.	Provides an alternative implementation of the finite-displacement method and ASR correction.
Converged Planewave Cutoff (ecutwfc)	Computational parameter determining basis set size.	A key "reagent" to reduce numerical noise and minimize the intrinsic FFT-grid ASR violation [42].

In the computational research of lattice dynamics, the appearance of non-zero acoustic modes serves as a critical indicator of the calculation's numerical health. A disciplined approach is required: manual ASR correction is warranted only for minor violations (< 10 cm⁻¹) after rigorously verifying the convergence of all parameters and the correctness of the crystal structure. For larger errors or significant imaginary frequencies, applying the ASR is a cosmetic fix that obscures underlying problems; the correct response is to revisit the convergence protocol and structural model. The dynmat.x tool provides the definitive test and solution for valid calculations suffering from the inevitable numerical approximations of the plane-wave method. By adhering to this diagnostic framework, researchers can ensure the reliability of their phonon spectra, leading to accurate predictions of thermodynamic properties and material stability, which are foundational to fields ranging from fundamental condensed matter physics to rational drug development through crystal structure prediction.

This technical guide provides an in-depth analysis of critical convergence parameters—tr2_ph, conv_thr, and k-point grid settings—within the context of phonon frequency calculations and Acoustic Sum Rule (ASR) violation research. Achieving accurate phonon calculations in density functional perturbation theory (DFPT) requires careful balancing of numerical precision and computational efficiency. This whitepaper synthesizes current methodologies for parameter optimization, troubleshooting common pathologies in phonon spectra, and implementing robust protocols for handling ASR violations. By establishing clear relationships between convergence parameters and numerical accuracy, we provide researchers with a structured framework for obtaining reliable phonon dispersions, particularly crucial for investigating material stability, phase transitions, and lattice dynamics in pharmaceutical development and materials science.

Phonon calculations form the cornerstone of lattice dynamics research, enabling predictions of material stability, thermal properties, and vibrational spectra. Within Quantum ESPRESSO's PHonon package (ph.x), three parameters predominantly govern convergence: tr2_ph (phonon calculation threshold), conv_thr (electronic convergence threshold from preceding pw.x calculations), and k-point sampling density. These parameters exhibit complex interdependencies that researchers must navigate to obtain physically meaningful results.

The Acoustic Sum Rule presents a particular challenge in phonon calculations. In an ideal crystal with perfect translational invariance, the frequencies of acoustic phonons should approach zero at the Brillouin zone center (Γ point). However, discrete numerical grids and finite basis sets break this exact invariance, leading to small but non-zero acoustic frequencies even in well-converged calculations [42]. The severity of ASR violations serves as a crucial diagnostic for the quality of convergence parameters—excessively large violations indicate inadequate parameter selection that can compromise subsequent analysis.

Parameter Specifications and Quantitative Guidelines

Core Parameter Definitions

• tr2_ph: This real-valued parameter sets the convergence threshold for the self-consistency of the DFPT calculation. The default value is 1e-12 in Quantum ESPRESSO. Lower values enforce stricter convergence criteria for the phonon calculation itself, directly impacting the accuracy of the dynamical matrix [44].

• conv_thr: Defined in the preceding PWscf (pw.x) calculation, this parameter controls the convergence threshold for the electronic ground state. Typical values range from 1e-8 to 1e-12, with tighter thresholds required for accurate force constants. The electronic structure must be well-converged before phonon calculations commence [42].

• k-point grid: The sampling density in reciprocal space critically influences both electronic and phonon calculations. For phonons, the grid must be commensurate with the phonon wavevectors (q-points), especially in metallic systems where convergence can be particularly slow [42].

Quantitative Parameter Recommendations

Table 1: Recommended Parameter Ranges for Different System Types

System Type	conv_thr (pw.x)	tr2_ph (ph.x)	k-point Density	Additional Notes
Insulators	1e-10 to 1e-12	1e-12 to 1e-14	8×8×8 or finer	Higher precision needed for dielectric properties
Metals	1e-10 to 1e-12	1e-12 to 1e-14	12×12×12 or finer	Denser k-points crucial with semicore states
Semiconductors	1e-10 to 1e-12	1e-12 to 1e-14	10×10×10 or finer	Intermediate values typically sufficient
Systems with Vacuum	1e-11 or tighter	1e-13 or tighter	8×8×1 or similar	Careful with spurious interactions

Table 2: Troubleshooting Guide for Parameter-Related Pathologies

Pathology	Primary Parameters to Adjust	Typical Adjustments	Expected Outcome
Large ASR violations (>10 cm⁻¹)	conv_thr, tr2_ph, k-grid	Tighten thresholds by 1-2 orders; densify k-grid	Acoustic modes <10 cm⁻¹
Negative frequencies	ecutwfc, ecutrho, k-grid	Increase cutoffs; improve k-sampling	Elimination of unphysical instabilities
Wrong degeneracies	k-grid, atomic positions	Densify k-grid; refine structure	Correct symmetry representation
Slow convergence	alpha_mix(niter), nmix_ph	Increase nmix_ph to 8-20	Faster SCF with memory tradeoff

Experimental Protocols for Parameter Optimization

Systematic Convergence Testing Methodology

A robust parameter optimization protocol requires sequential testing of each parameter while holding others constant:

• Begin with k-point convergence: First, establish a well-converged k-point grid for the electronic ground state. Monitor total energy changes as the k-mesh density increases; convergence within 1-5 meV/atom typically indicates sufficient sampling. For metals, pay particular attention to the density of states at the Fermi level.

• Optimize conv_thr in PWscf: With a fixed k-point grid, systematically tighten conv_thr from 1e-8 to 1e-12 while monitoring the total energy and forces. The electronic structure must be fully converged before proceeding to phonon calculations, as residual errors propagate into the force constants [42].

• Refine tr2_ph for phonon convergence: Using the converged electronic ground state, perform phonon calculations while varying tr2_ph. Monitor the magnitude of ASR violations—acoustic mode frequencies at Γ should ideally be below 10 cm⁻¹, though values up to 100 cm⁻¹ may occur in problematic cases [42].

• Iterative refinement: Return to step 1 if phonon convergence proves unsatisfactory, as parameters exhibit interdependencies. For complex systems, this cycle may require multiple iterations.

ASR Violation Diagnostic Protocol

When ASR violations exceed acceptable thresholds (>10 cm⁻¹), implement this diagnostic procedure:

• Initial assessment: Calculate phonons at Γ and measure the acoustic mode frequencies. Frequencies significantly above zero indicate potential convergence issues or numerical instabilities.

• Parameter adjustment: Systematically tighten conv_thr and tr2_ph while monitoring computational cost. As Stefano de Gironcoli notes, "Increasing the charge density cutoff increases the grid density thus making the integral more exact thus reducing the problem, unfortunately rather slowly... This problem is usually more severe for GGA than with LDA" [42].

• Symmetry verification: Use dynmat.x to diagonalize the dynamical matrix with imposed ASR. If this yields acoustic modes with ω < 1 cm⁻¹ while other modes remain unchanged, the original results can be trusted despite the ASR violation [42].

• Alternative approach: For severe cases, consider finite-displacement methods (e.g., phonopy) as an alternative to DFPT, which may exhibit different convergence characteristics [33].

The Scientist's Toolkit: Essential Research Reagents

Table 3: Critical Computational Tools for Phonon Convergence Research

Tool/Component	Function	Implementation Example
Quantum ESPRESSO Suite	First-principles DFT/DFPT calculations	pw.x for SCDFT, ph.x for phonons [44]
phonopy	Finite displacement method for phonons	Alternative approach to DFPT: phonopy --qe -d --dim="2 2 2" -c NaCl.in [33]
dynmat.x	ASR imposition and dynamical matrix analysis	Post-processing tool to enforce ASR and verify acoustic modes [42]
Born Effective Charges	Non-analytical term correction	phonopy-qe-born to generate BORN file for LO-TO splitting [33]
K-point Convergence Tools	Reciprocal space sampling optimization	Automatic k-mesh generation with convergence monitoring
Symmetry Analysis Tools	Space group verification	Detection of symmetry-breaking atomic positions

Advanced Considerations and Methodological Extensions

System-Specific Convergence Challenges

Different material classes present distinct convergence challenges that require specialized parameter strategies:

• Metallic systems: Present particular difficulties due to Fermi surface effects. The convergence with respect to k-point grid and smearing is very slow in the presence of semicore states, and for phonon wave-vectors that are not commensurate with the k-point grid [42]. Denser k-meshes (12×12×12 or finer) combined with appropriate smearing techniques (e.g., Methfessel-Paxton, Marzari-Vanderbilt) are essential.

• Systems with vacuum regions: Surfaces, interfaces, and 2D materials require careful treatment of vacuum spacing and increased real-space cutoffs. The convergence problems are exacerbated because "in the exponential tail of the charge density a) the finite cutoff induces oscillations in rho and b) the reduced gradient is diverging" [42].

• Mixed metallic-covalent systems: Complex bonding environments may require hybrid functional approaches (HSE) for accurate electronic structures, which subsequently influence phonon convergence. As demonstrated in CaFI monolayer research, advanced functionals can improve band gap predictions and subsequent lattice dynamical properties [45].

Methodological Cross-Validation

Robust phonon research necessitates cross-validation between different computational approaches:

• DFPT vs. finite displacement: Quantum ESPRESSO's native DFPT implementation can be cross-validated against finite-displacement methods using phonopy, which creates supercells with atomic displacements for force calculations [33]. Discrepancies between methods can reveal underlying convergence issues.

• ASR imposition techniques: Compare phonon spectra before and after ASR imposition using dynmat.x. Significant changes in optical modes indicate fundamental convergence problems, while changes restricted to acoustic modes suggest manageable numerical errors.

• Advanced electron-phonon coupling: For properties dependent on electron-phonon interactions (e.g., superconductivity, temperature-dependent optical properties), ensure convergence of both electronic and phononic components. The electron_phonon flag in ph.x provides various approaches including interpolation and tetrahedron methods [44].

Optimizing tr2_ph, conv_thr, and k-point grid settings represents a critical step in obtaining reliable phonon frequencies with controlled ASR violations. This guide establishes that systematic parameter testing, beginning with k-point convergence and progressively refining electronic and phononic thresholds, provides the most robust pathway to accurate lattice dynamical properties. The Acoustic Sum Rule serves as an essential diagnostic tool throughout this process, with violations below 10 cm⁻¹ representing a typical target for well-converged calculations.

Researchers should maintain awareness of system-specific challenges, particularly for metals and low-dimensional systems where convergence may be slower and require specialized techniques. The provided protocols, troubleshooting guidelines, and computational toolkit offer a comprehensive framework for navigating these complex parameter spaces. By implementing these methodologies, computational researchers can enhance the reliability of their phonon calculations, leading to more accurate predictions of material stability, thermal properties, and lattice dynamics across pharmaceutical and materials science applications.

Imaginary frequencies in phonon calculations represent a significant challenge in computational materials science, often indicating structural instabilities or numerical inaccuracies. This technical guide provides an in-depth examination of the origins and solutions for imaginary frequencies within the broader context of acoustic sum rule (ASR) violation and phonon frequency calculation research. We systematically outline diagnostic protocols and remediation strategies, supported by quantitative data from first-principles calculations. Specifically, we detail the critical roles of structural optimization and cutoff energy selection in achieving dynamically stable configurations. By integrating advanced computational techniques, including electron-phonon coupling considerations and machine-learning assisted approaches, this work establishes a comprehensive framework for addressing phonon instabilities across diverse material systems, from metal-organic frameworks to topological quantum materials.

Phonon calculations constitute an essential component of computational materials science, enabling the prediction of thermal properties, phase stability, and vibrational spectra. The emergence of imaginary frequencies (negative values in phonon dispersion relations) presents a fundamental challenge, indicating that the atomic configuration resides at a saddle point rather than a local minimum on the potential energy surface. Within the framework of acoustic sum rule (ASR) violation research, these imaginary modes reflect violations of translational invariance that necessitate careful numerical treatment [28].

The accurate calculation of phonon frequencies relies on the precise determination of interatomic force constants (IFCs), which can be computed using either the finite-displacement (frozen-phonon) method or density functional perturbation theory (DFPT) [46]. While these approaches have become standard in materials simulation packages, the persistence of imaginary frequencies remains a common obstacle requiring systematic investigation. This guide examines the interconnected roles of structural optimization and cutoff energy parameters in resolving these instabilities, providing researchers with validated protocols for achieving dynamically stable configurations.

Theoretical Foundations: Imaginary Frequencies and ASR Violation

Physical Interpretation of Imaginary Frequencies

Imaginary frequencies in phonon spectra signify vibrational modes with negative curvature on the potential energy surface, implying that the system can lower its energy through atomic displacements along these modes. Physically, this indicates one of three scenarios:

• Structural instability: The crystal structure is not at a local energy minimum
• Saddle point configuration: The system resides at a transition state between minima
• Numerical artifacts: Insufficient computational accuracy creates spurious instabilities [28]

The presence of small imaginary frequencies near the Brillouin zone center (Γ-point) particularly suggests ASR violations, where the dynamical matrix fails to satisfy the condition that acoustic phonons should have zero frequency at Γ-point due to translational invariance [28].

Electron-Phonon Coupling and Advanced Considerations

Recent investigations have revealed that electron-phonon coupling (EPC) can significantly influence phonon stability and behavior. In materials with topologically nontrivial electronic structures, such as BaMnSb₂, phonons can "inherit" a Berry curvature through their coupling to electrons, potentially affecting their stability and spectral properties [17]. Furthermore, in magnetic systems, EPC induces phonon magnetic moments and Zeeman splittings, adding complexity to phonon dispersion relations [47]. These advanced considerations highlight the multifaceted nature of phonon instabilities beyond simple structural optimization.

Methodological Framework: Computational Approaches

Fundamental Phonon Calculation Techniques

Table 1: Computational Methods for Phonon Calculations

Method	Principle	Advantages	Limitations
Frozen-phonon (finite-displacement)	Direct supercell approach with atomic displacements	Simple implementation, anharmonic capabilities	Combinatorial explosion for high-order IFCs [46]
Density Functional Perturbation Theory (DFPT)	Linear response in reciprocal space	Efficient for harmonic IFCs, no supercells required	Complex implementation, limited anharmonicity [46]
Compressive Sensing Lattice Dynamics (CSLD)	Sarse recovery of IFCs via machine learning	Handles high-order anharmonicity efficiently	Requires careful regularization [46]
Variational Polaron Equations (VarPEq)	Self-consistent optimization of polaron wavefunction	Avoids charged supercells, includes strong EPC	Computationally intensive, complex implementation [48]

Research Reagent Solutions: Computational Tools

Table 2: Essential Software Tools for Phonon Calculations

Tool Name	Function	Application Context
Phonopy	Frozen-phonon calculations	Harmonic and thermal properties [11]
Phonopy	Third-order anharmonic IFCs	Lattice thermal conductivity [46]
Pheasy	High-order IFC extraction	Anharmonic lattice dynamics [46]
Abinit	DFPT and electron-phonon coupling	Polaronic systems and response properties [48]
Alamode	Compressive sensing lattice dynamics	Efficient anharmonic IFC extraction [46]
ShengBTE	Boltzmann transport equation	Thermal conductivity from phonons [46]

A systematic approach to diagnosing imaginary frequencies requires sequential investigation of potential sources, as illustrated in the following workflow:

Acoustic Sum Rule Verification

For imaginary frequencies specifically at the Γ-point, first verify ASR compliance. The ASR requires that the sum of the force constants between one atom and all other atoms must be zero, ensuring translational invariance. Violations occur due to numerical approximations in force constant calculations. Modern computational packages like Phonopy and Abinit include ASR correction options that enforce this condition by adjusting the dynamical matrix [28].

Structural Optimization Assessment

Incomplete geometry optimization represents the most common source of imaginary frequencies. The optimization process must converge to a true local minimum where all Hellmann-Feynman forces fall below a stringent threshold (typically < 1 meV/Å) and stress tensor components are minimized. For complex systems like metal-organic frameworks (MOFs) with flexible building units, optimization must carefully balance bond stretches, angle bends, and dihedral torsions to achieve equilibrium [49].

Remediation Strategies: Structural Optimization Protocols

Comprehensive Geometry Optimization

Achieving a fully optimized structure requires sequential refinement:

• Initial Optimization: Begin with standard optimization criteria (forces < 10 meV/Å)
• Lattice Parameter Optimization: For crystalline systems, optimize both atomic positions and lattice parameters simultaneously
• High-Precision Optimization: Apply stringent convergence criteria (forces < 1 meV/Å) using increased cutoff energy and k-point density
• Symmetry Breaking: For persistent imaginary frequencies, apply small random displacements (0.01-0.05 Å) to break symmetry, followed by reoptimization to escape saddle points [28]

For MOFs and complex crystals, specialized flexibility parameterization protocols using force matching to ab initio molecular dynamics (AIMD) trajectories have proven effective for obtaining accurate force constants [49].

Electron-Phonon Coupling Considerations

In materials with strong EPC, standard harmonic approximations may fail to capture essential physics. The Variational Polaron Equations (VarPEq) approach provides an advanced framework for addressing these cases, reformulating the polaron formation problem in reciprocal space and avoiding the need for charged supercell relaxation [48]. This method explicitly treats the coupling between electronic states and lattice vibrations, which can stabilize otherwise imaginary modes.

Cutoff Energy Optimization: Quantitative Guidelines

Cutoff energy selection critically influences the accuracy of plane-wave DFT calculations, directly impacting force constant determination. The following protocol establishes systematic criteria for cutoff optimization:

Cutoff Energy Convergence Protocol

• Initial Scan: Perform total energy calculations across a range of cutoff energies (e.g., 200-800 eV in 50 eV increments)
• Energy Convergence: Identify the cutoff where total energy changes < 1 meV/atom
• Force Convergence: Verify that atomic forces converge to < 0.1 meV/Å at the selected cutoff
• Stress Tensor Validation: Ensure stress tensor components show oscillatory variations < 0.1 GPa
• Phonon Validation: Confirm that phonon frequencies, particularly soft modes, show variations < 0.1 THz with further cutoff increases

Table 3: Cutoff Energy Convergence Data for Representative Materials

Material Class	Typical Cutoff Range (eV)	Convergence Threshold (meV/atom)	k-point Density (Å⁻¹)	Force Tolerance (meV/Å)
Elemental Metals	300-500	0.5	0.03-0.05	0.5
Semiconductors	400-600	0.2	0.02-0.04	0.2
Oxides/Ceramics	500-700	0.1	0.03-0.05	0.1
MOFs/Porous Materials	550-750	0.05	0.01-0.02	0.05
2D Materials	500-700	0.1	0.02-0.03	0.1

Reciprocal Space Integration Parameters

Alongside cutoff energy, k-point sampling significantly influences computational accuracy:

• k-point Convergence: Perform systematic convergence tests for phonon properties, not just total energy
• Γ-centered Meshes: Prefer Γ-centered k-point grids for phonon calculations to ensure accurate treatment of long-wavelength modes
• q-point Independence: Phonon calculations require dense q-point meshes independent of electronic k-point sampling

For systems with persistent imaginary frequencies, increasing the k-point density often proves more effective than further increasing cutoff energy beyond the convergence threshold [50].

Advanced Considerations: Anharmonicity and Machine Learning Approaches

Anharmonic Lattice Dynamics

For strongly anharmonic systems, the standard harmonic approximation inherently predicts imaginary frequencies that disappear when anharmonic effects are properly included. The self-consistent harmonic approximation (SCHA) and related techniques effectively address this limitation by incorporating temperature-dependent phonon renormalization [46]. Implementation in codes like Pheasy enables the extraction of high-order interatomic force constants using machine-learning algorithms, providing a more physically accurate treatment of materials with intrinsic anharmonicity.

Machine-Learning Assisted Protocols

Recent advances integrate machine learning to accelerate phonon property prediction and instability diagnosis. High-throughput screening of lattice dynamics signatures, as demonstrated for sodium superionic conductors, enables rapid identification of stability descriptors across large materials databases [51]. These approaches leverage phonon-derived features such as mean squared displacement (MSD) and acoustic cutoff frequencies to predict dynamic stability without full first-principles calculations for every candidate material.

Addressing imaginary frequencies in phonon calculations requires a multifaceted approach combining rigorous structural optimization, careful parameter convergence, and consideration of advanced physical effects like electron-phonon coupling and anharmonicity. The protocols outlined in this work provide a systematic framework for resolving phonon instabilities across diverse material systems. Future research directions will likely focus on integrating machine-learning potentials for enhanced force constant accuracy and developing more efficient treatments of anharmonic effects in complex materials. As computational methodologies advance, the robust prediction of phonon stability will continue to enable the discovery and design of novel materials with tailored thermal and dynamic properties.

In first-principles density functional theory (DFT) calculations, particularly those investigating lattice dynamics and phonon frequencies, the precise definition of the crystal structure is paramount. Symmetry-related errors represent a frequent and significant obstacle in quantum materials simulation, often leading to incorrect force constants, violation of the acoustic sum rule (ASR), and the appearance of unphysical imaginary phonon frequencies. These errors fundamentally compromise the predictive capability of computational materials science, especially in advanced fields like drug development where accurate modeling of intermolecular vibrations is crucial for understanding molecular interactions and stability.

The core of this problem often lies in the incorrect specification of the Bravais lattice and atomic positions within input files for computational codes such as Quantum ESPRESSO. When atomic positions deviate even slightly from their exact symmetry locations, or when the lattice type is misidentified, the code's symmetry detection algorithms fail, causing a cascade of computational inaccuracies. This technical guide provides a comprehensive framework for researchers to correctly configure the ibrav parameter and utilize Wyckoff positions, ensuring robust phonon calculations that respect the fundamental physical constraints imposed by crystal symmetry.

Theoretical Foundation: Symmetry, ASR, and Phonon Calculations

The Acoustic Sum Rule and Its Violation

The Acoustic Sum Rule (ASR) embodies the principle of translational invariance in crystals, mandating that the total energy remains unchanged under an infinitesimal translation of the entire crystal system. This physical principle mathematically guarantees that the three phonon modes at the Brillouin zone center (Γ-point) exhibit precisely zero frequency, corresponding to acoustic waves. In practical computation, however, this rule is often violated due to numerical approximations.

The primary source of ASR violation stems from the discreteness of the Fast Fourier Transform (FFT) grid used in plane-wave calculations [3]. Other contributing factors include insufficient convergence parameters (such as tr2_ph for phonons and conv_thr for the electronic ground state) and the use of finite supercells for isolated molecules. When the ASR is violated, phonon calculations yield non-zero frequencies for acoustic modes at the Γ-point, undermining the physical validity of the entire simulation. As stated in the Quantum ESPRESSO documentation, "if the nonzero frequencies are small, you can impose the ASR to the dynamical matrix, usually with excellent results" [3].

Symmetry's Role in Computational Efficiency and Accuracy

Beyond physical correctness, symmetry plays a crucial computational role by significantly reducing the workload. Symmetry operations allow codes to calculate properties (like wavefunctions, forces, and phonon frequencies) only for symmetry-inequivalent atoms and k-points, subsequently generating the full set through transformation. This reduction can lead to dramatic computational savings. More importantly, correct symmetry identification ensures that degeneracies in the electronic band structure and phonon dispersions are properly captured, which is essential for predicting transport properties and structural instabilities.

Correct Lattice Definition: TheibravParameter in Quantum ESPRESSO

The ibrav parameter in Quantum ESPRESSO's &SYSTEM namelist defines the Bravais lattice of the crystal structure [52]. Its proper use is the first and most critical step in avoiding symmetry-related errors.

1ibrav≠ 0: Utilizing Built-in Lattice Symmetry

When ibrav is set to a non-zero value, the code leverages its internal database of lattice symmetries. The user provides fundamental lattice constants (e.g., celldm(1) or A, B, C and angles), and the code constructs the full unit cell based on the chosen lattice type [53].

Table 1: Common ibrav Settings and Corresponding Lattice Parameters in Quantum ESPRESSO

ibrav	Lattice Type	Required Parameters	Key Consideration
2	Face-Centered Cubic (FCC)	celldm(1) = lattice constant in Bohr	For FCC silicon, only two atomic positions needed: Si at (0,0,0) and (0.25,0.25,0.25) [53]
4	Hexagonal	celldm(1) = a, celldm(3) = c/a	Preferred over ibrav=-5 for true hexagonal structures [54]
-5	Trigonal	celldm(1) = a, celldm(4) = cos(α)	Requires a=b=c and α=β=γ; conflicting parameters cause symmetry errors [54]

For any ibrav ≠ 0, the user must provide only the coordinates of symmetry-inequivalent atoms in the primitive unit cell. The code automatically generates all equivalent atoms by applying the space group symmetry operations [53]. This approach is generally preferred because it inherently minimizes the risk of symmetry-breaking input.

2ibrav= 0: The Flexible but Risky Alternative

Setting ibrav = 0 allows complete flexibility in defining the unit cell. The user must explicitly provide the full set of CELL_PARAMETERS and all atomic coordinates in ATOMIC_POSITIONS [53] [52]. While this offers maximum control, it places the full burden of symmetry correctness on the researcher. The code still attempts to identify symmetries from the provided geometry, but this process is highly sensitive to numerical precision. Atoms positioned even slightly away from their exact symmetry locations can cause the symmetry detection to fail, leading to the infamous "some problem with symmetries" error [54]. When using ibrav=0, you can provide coordinates for either a primitive or conventional cell, but the former is computationally more efficient [53].

Practical Implementation: Wyckoff Positions and Workflow

The Wyckoff Position Approach

For crystal structures with known space groups, the most robust method is to define the structure using Wyckoff positions. This involves specifying the space group number (or Hermann-Mauguin symbol) and providing only the independent atomic coordinates within the asymmetric unit. Quantum ESPRESSO then applies the space group's symmetry operations to build the complete crystal structure [53]. This method guarantees that all generated atoms respect the crystal symmetry by construction, effectively eliminating one major source of symmetry-related errors.

Step-by-Step Input Workflow

The following diagram illustrates the logical decision process for configuring crystal structures in Quantum ESPRESSO to avoid symmetry errors:

Research Reagent Solutions: Computational Tools

Table 2: Essential Computational Tools for Symmetry-Correct Phonon Calculations

Tool/Parameter	Function/Role in Symmetry Handling	Implementation Notes
ibrav	Defines Bravais lattice type	Non-zero values activate internal symmetry database; "If such errors occur, set the Bravais lattice using the correct ibrav value (i.e., don't use ibrav=0)" [3]
space_group	Directly specifies crystal space group	Used with Wyckoff positions to build complete structure from asymmetric unit [53] [52]
uniqueb	Determines unique axis for specific space groups	Must be set consistently with space group requirements [52]
origin_choice	Selects between alternate origin choices	Critical for certain space groups with multiple conventional settings [52]
Acoustic Sum Rule (ASR)	Corrects for translational invariance violation	"If the nonzero frequencies are small, you can impose the ASR to the dynamical matrix" [3]
Symmetry Analysis Tools	Diagnose symmetry issues in output structures	ISOTROPY package: http://stokes.byu.edu/iso/isotropy.html [3]

Troubleshooting Common Symmetry Errors

Diagnostic Patterns and Solutions

Symmetry-related errors manifest in several characteristic ways in Quantum ESPRESSO calculations:

• "some problem with symmetries" Error: This fatal error during the initial setup typically indicates a fundamental mismatch between the defined ibrav, the provided lattice parameters, and the atomic positions [54]. For example, setting ibrav=-5 (trigonal) while providing lattice parameters where a=b≠c and β=γ≠α creates conflicting symmetry information.

• Non-zero Acoustic Phonon Frequencies: Small but non-zero frequencies for acoustic modes at Γ-point indicate ASR violation, which may be addressed by enforcing the sum rule in phonon calculations [3].

• Mysterious Symmetry-Related Errors: Errors mentioning "symmetry operation is non orthogonal," "wrong representation," or "wrong degeneracy" almost invariably stem from atomic positions that are close to, but not exactly at, their symmetry positions [3].

Case Study: Correcting a Trigonal/Hexagonal Confusion

A documented case on the PW Forum illustrates a classic symmetry error [54]. A user attempted to calculate a structure with ibrav=-5 (trigonal) but provided lattice parameters corresponding to a hexagonal system. The expert response clarified:

"ibrav=-5 corresponds to a trigonal lattice, for which a=b=c and the same angle alpha between ANY pair of lattice vectors (alpha=beta=gamma). Probably, the setting you give a=b/=c and beta=gamma/=alpha OVERRIDES any default and you end up with conflicting settings: ibrav=-5 tells QE to use a trigonal lattice and on this basis some symmetries are expected, but based on the input lattice parameters (A, B, C, cosAB, cosAC, cosBC) a lattice with different symmetries is set up."

The solution was to use ibrav=4 for the hexagonal lattice instead, which immediately resolved the symmetry error and produced a physically reasonable crystal structure [54].

Atomic Position Pitfalls: Avoiding Overlapping Atoms

A critical and often overlooked source of symmetry errors involves specifying equivalent atomic positions in the input file, particularly when using crystal coordinates. As highlighted in the forum discussion, providing coordinates like:

results in atoms that are translationally equivalent through periodic boundary conditions [54]. These "overlapping" atoms break the symmetry detection and lead to unphysical results. The solution is to provide only the symmetry-inequivalent atoms in the primitive unit cell, letting the code generate all periodic images automatically.

Integration with Phonon Frequency Calculations

From Correct Symmetry to Accurate Phonons

The integrity of the initial crystal structure symmetry directly propagates to the quality of phonon calculations. Correct symmetry ensures that:

• Force constants are properly symmetrized across all equivalent atom pairs
• Phonon degeneracies are correctly captured along high-symmetry paths
• Acoustic modes naturally satisfy the ASR at the Γ-point without aggressive numerical corrections
• Phonon dispersions display the correct physical behavior throughout the Brillouin Zone

When phonon calculations yield unexpected imaginary frequencies (negative ω²), the first checkpoint should be the structural symmetry, followed by convergence parameters like tr2_ph and the electronic conv_thr [3].

Advanced Considerations: Electron-Phonon Interactions

For researchers investigating charge transport properties or ultrafast carrier dynamics using packages like Perturbo that build upon Quantum ESPRESSO outputs, correct phonon spectra are indispensable [55] [56]. These advanced calculations require precise electron-phonon matrix elements that rely heavily on the underlying symmetry of both electronic and vibrational states. Symmetry errors in the initial DFT calculation propagate through the entire computational pipeline, compromising the predictive accuracy for properties like electrical conductivity, mobility, and Seebeck coefficient.

Ensuring symmetry correctness in Quantum ESPRESSO calculations requires a systematic approach from the initial structure definition. The following consolidated practices will significantly reduce symmetry-related errors:

• Prefer ibrav ≠ 0 whenever possible to leverage the code's internal symmetry database.
• Utilize Wyckoff positions for structures with known space groups to guarantee symmetry by construction.
• Validate atomic positions to prevent inclusion of translationally equivalent atoms.
• Verify lattice consistency by ensuring provided lattice parameters match the chosen ibrav symmetry.
• Apply ASR corrections in phonon calculations when small non-zero acoustic frequencies persist.
• Use high numerical precision in atomic coordinates, especially when employing ibrav=0.

By adhering to these guidelines and understanding the intimate connection between symmetry definition, ASR fulfillment, and phonon frequency accuracy, researchers can establish a robust foundation for reliable materials simulations across diverse applications, from fundamental condensed matter physics to pharmaceutical development.

Validating Phonon Results and Cross-Method Comparisons: Ensuring Computational Accuracy

Validating computational methods against known results from prototypical materials is a fundamental practice in phonon physics research. This process is crucial for verifying the accuracy of force constant calculations, assessing the proper handling of physical constraints like the acoustic sum rule (ASR), and ensuring the correct treatment of long-range interactions in polar materials. Within the broader context of acoustic sum rule violation and phonon frequency calculation research, benchmarking serves as an essential quality control measure, helping researchers identify computational artifacts, diagnose convergence issues, and build confidence in their methodologies before applying them to novel materials systems. This technical guide examines established benchmarking protocols using well-characterized prototypical materials, with particular emphasis on addressing ASR violations and other common computational challenges in lattice dynamics simulations.

Core Concepts in Phonon Calculation Validation

Acoustic Sum Rule Violations

The acoustic sum rule embodies the physical requirement of translational invariance in crystalline solids, mandating that the sum of interatomic forces acting on any atom when all atoms are displaced uniformly must be zero. Violations of this rule represent a significant challenge in phonon calculations. In practical implementations, approximated calculations frequently violate the ASR due to numerical discretization, with the discreteness of the FFT grid constituting the "most irreducible violation" in plane-wave calculations [3]. Additional factors include insufficient convergence parameters, particularly the tr2ph threshold for phonons and convthr for electronic ground state calculations [3].

These violations manifest computationally as non-zero acoustic mode frequencies at the Brillouin zone center (Γ-point) when they should theoretically be zero. When these frequencies are small, the ASR can be imposed directly on the dynamical matrix, typically with "excellent results" [3]. For molecular systems, non-zero frequencies for rotational modes represent a "fictitious effect of the finite supercell size" or imperfect convergence of the molecular geometry [3].

Imaginary Frequencies and Physical Instabilities

The appearance of imaginary frequencies (mathematically represented as ω² ≤ 0) in phonon spectra requires careful interpretation. While sometimes indicative of computational issues such as "bad convergence," they may also signal a "real instability" corresponding to a structural phase transition [3]. Proper benchmarking helps distinguish between these possibilities, enabling researchers to identify when imaginary frequencies reflect physical phenomena versus numerical artifacts.

Long-Range Interactions in Polar Materials

Polar materials present unique challenges due to the long-range dipole-dipole interactions that arise when two or more atoms in the unit cell carry non-zero Born effective charge tensors [4]. These interactions necessitate specialized treatment through Ewald summation techniques to properly account for the macroscopic electric fields that contribute to phonon frequencies [4]. Without proper handling, calculations fail to capture the characteristic LO-TO splitting observed experimentally in polar semiconductors and insulators.

Table 1: Common Computational Challenges in Phonon Calculations

Challenge	Physical Cause	Computational Manifestation	Remedial Approaches
ASR Violation	Broken translational invariance	Non-zero acoustic modes at Γ-point	Impose ASR on dynamical matrix
Imaginary Frequencies	Structural instability or poor convergence	Negative ω² in phonon spectra	Improve convergence or investigate instability
Improper LO-TO Splitting	Unaccounted dipole-dipole interactions	Discontinuous phonon dispersions near Γ	Enable Ewald summation with Born charges
Slow Convergence	Insufficient k-point sampling or energy cutoff	Unphysical oscillations in phonon bands	Increase mesh density and cutoff energy

Prototypical Materials for Benchmarking

Selection Criteria for Benchmark Materials

Effective benchmark materials possess well-characterized phonon properties with extensive experimental and theoretical data available for comparison. Ideal candidates exhibit high structural symmetry, precisely known lattice parameters, and comprehensive experimental phonon measurements. For validating treatments of long-range interactions, materials with significant polar character and measurable LO-TO splitting are essential. The materials should also represent different crystal classes and bonding types to ensure methodological transferability.

Magnesium Oxide (MgO) - A Cubic Ionic Crystal

Magnesium oxide crystallizes in the rock-salt structure (face-centered cubic) and serves as an excellent benchmark for ionic systems. Its high symmetry simplifies calculations while its strongly polar nature produces substantial LO-TO splitting. For MgO, typical Born effective charges are approximately ±1.97|e|, significantly enhanced over the nominal ionic values, while the static dielectric tensor is approximately 3.13 along all principal axes [4].

When LO-TO splitting is properly included using the Ewald summation technique with the correct dielectric tensor and Born effective charges, MgO's phonon dispersion shows smooth, physically realistic bands without unphysical oscillations [4]. The magnitude of the LO-TO splitting in MgO is "considerably large, on the same order of magnitude as the LO phonon frequencies" [4], making it a stringent test for computational methods.

Aluminum Nitride (AlN) - A Hexagonal Polar Semiconductor

Aluminum nitride adopts the wurtzite structure (hexagonal) with reduced symmetry compared to MgO, presenting additional benchmarking challenges. Its anisotropic structure results in direction-dependent Born effective charges and dielectric constants [4]. This anisotropy produces characteristic direction-dependent phonon frequencies and discontinuities around the Γ-point as q→0 [4].

The successful reproduction of these directional discontinuities provides critical validation for the proper implementation of long-range corrections in phonon codes. AlN's combination of anisotropy and polar character makes it an essential benchmark for verifying that computational methods can handle complex electrostatic interactions in low-symmetry systems.

Table 2: Prototypical Benchmark Materials and Their Key Properties

Material	Crystal Structure	Key Phonon Features	Validation Parameters
Magnesium Oxide (MgO)	Cubic (rock-salt)	Strong LO-TO splitting, acoustic modes	Born charges (~1.97	e	), dielectric tensor (~3.13)
Aluminum Nitride (AlN)	Hexagonal (wurtzite)	Anisotropic LO-TO splitting, directional discontinuities	Direction-dependent Born charges and dielectric constants
Silicon	Cubic (diamond)	Non-polar, well-known dispersion	No LO-TO splitting, three acoustic branches
Germanium	Cubic (diamond)	Similar to silicon but with softer modes	Test of transferability across isoelectronic systems

Computational Methodologies and Protocols

Force Constant Calculation Approaches

Two primary methodologies exist for computing interatomic force constants (IFCs) from first principles:

The finite displacement method (frozen phonon approach) employs real-space supercells where atoms are systematically displaced, and the resulting forces are calculated using density functional theory [4]. This approach is implemented in codes such as Phonopy and Phon [46] and is accessible through VASP with the IBRION = 5 or 6 tags [4]. While conceptually straightforward, this method requires careful convergence with respect to supercell size to ensure that force constants vanish appropriately at large distances [4].

The density functional perturbation theory (DFPT) approach computes force constants directly from the linear response of the electronic system to atomic displacements [46]. Implemented in packages like Quantum ESPRESSO and Abinit, and accessible in VASP via IBRION = 7 or 8 [4], DFPT operates in reciprocal space and generally offers superior computational efficiency for larger systems. Recent advances in compressive-sensing lattice dynamics (CSLD) and other machine-learning-inspired approaches, as implemented in the Pheasy code, enable efficient extraction of high-order IFCs by leveraging the physical reality that "IFCs of any system are generally sparse and decay rapidly with increasing interatomic distance" [46].

Treatment of Long-Range Electrostatic Interactions

For polar materials, the accurate calculation of phonon dispersions requires special treatment of long-range dipole-dipole interactions. The standard approach involves:

First, obtaining the static dielectric tensor and Born effective charge tensors through a separate DFPT calculation in the primitive unit cell with LEPSILON or LCALCEPS tags enabled [4]. These calculations must be properly converged with respect to k-point mesh and energy cutoff since "optical phonon frequencies depend strongly on the dielectric properties" [4].

Next, these properties are incorporated into the supercell calculation by setting LPHONPOLAR = True and providing the dielectric tensor (PHONDIELECTRIC) and Born charges (PHONBORNCHARGES) as inputs [4]. Optionally, a reciprocal space cutoff radius (PHONGCUTOFF) can be specified for the Ewald summation [4].

Recent methodological advances in codes like Pheasy incorporate "dynamical quadrupolar and octupolar effects that are crucial for accurately interpolating and converging the phonon dispersions of piezoelectric materials," going beyond the standard dipole-dipole approximation [46].

Phonon Dispersion and DOS Workflow

The following diagram illustrates the comprehensive workflow for computing phonon dispersions and density of states, incorporating validation against prototypical materials:

Experimental Protocols and Validation Metrics

Protocol for Phonon Dispersion Calculation

Step 1: Force Constant Computation

• For finite displacement method: Use IBRION = 5 or 6 in VASP to generate displaced configurations and compute forces [4]
• For DFPT approach: Use IBRION = 7 or 8 for efficient force constant calculation [4]
• Ensure supercell size is sufficiently large to converge force constants [4]

Step 2: Q-Path Specification for Dispersion

• Create a QPOINTS file containing a high-symmetry path through the Brillouin zone [4]
• Use established tools to determine appropriate high-symmetry paths and labels for the crystal structure [4]
• Set LPHON_DISPERSION = .TRUE. in the INCAR file [4]

Step 3: LO-TO Splitting Treatment for Polar Materials

• Obtain static dielectric tensor and Born effective charges from a separate DFPT calculation [4]
• Set LPHONPOLAR = .TRUE. and provide PHONDIELECTRIC and PHONBORNCHARGES tensors [4]
• Optionally specify PHONGCUTOFF for Ewald summation [4]

Protocol for Phonon Density of States Calculation

Step 1: Uniform Q-Mesh Generation

• Create a QPOINTS file specifying a sufficiently dense, uniform q-point mesh [4]
• Ensure mesh density provides converged DOS features

Step 2: DOS Computation

• Set PHON_DOS > 0 in the INCAR file [4]
• The DOS is automatically computed between [ωmin-5σ, ωmax+5σ] where σ is the broadening (PHON_SIGMA) [4]
• Specify the number of energy points with PHON_NEDOS [4]
• Choose Gaussian (PHONDOS = 1) or tetrahedron (PHONDOS = 2) smearing method [4]

Key Validation Metrics and Tolerances

Successful benchmarking requires quantitative assessment against established tolerances:

Acoustic Mode Compliance: Acoustic modes at Γ-point should approach zero within a reasonable tolerance. Frequencies below 1-2 cm⁻¹ typically indicate acceptable ASR compliance, while larger values may require imposition of the ASR or improved convergence [3].

LO-TO Splitting Verification: For polar materials like MgO and AlN, verify that the characteristic splitting between longitudinal and transverse optical modes matches established experimental and theoretical references [4].

Phonon Band Smoothness: Interpolated phonon bands should appear smooth without unphysical oscillations, particularly when including long-range corrections for polar materials [4].

Imaginary Frequency Assessment: Distinguish between physically meaningful imaginary frequencies indicating structural instabilities versus computational artifacts resulting from poor convergence [3].

The Scientist's Toolkit: Essential Research Reagents

Table 3: Essential Computational Tools for Phonon Calculations

Tool/Code	Primary Function	Key Features	Applicable Systems
VASP	DFT Force Calculation	IBRION parameters for phonons, DFPT support	Metals, insulators, surfaces
Phonopy	Force Constant Extraction	Finite displacement method, post-processing	Crystals with periodic boundary conditions
Pheasy	Advanced IFC Extraction	Machine-learning IFCs, high-order anharmonicity	Strongly anharmonic systems
Quantum ESPRESSO	DFPT Phonon Calculation	Open-source, comprehensive phonon capabilities	Insulators, metals (limited)
ALAMODE	Anharmonic Phonon Calculations	Compressive sensing lattice dynamics	Thermal transport, anharmonic systems
ISOTROPY	Symmetry Analysis	Phonon mode symmetry characterization	High-symmetry crystals

Advanced Considerations: Anharmonicity and High-Order Force Constants

Going beyond the harmonic approximation presents additional benchmarking challenges. The calculation of anharmonic IFCs using conventional finite-displacement or DFPT approaches becomes computationally prohibitive due to the "combinatorial explosion in the number of high-order IFCs" [46]. For third-order anharmonicity, specialized codes like thirdorder.py and Phono3py are available, while fourthorder.py extends this to fourth-order IFCs, albeit with significant computational constraints [46].

Advanced approaches like compressive-sensing lattice dynamics (CSLD) leverage the physical reality that "IFCs of any system are generally sparse and decay rapidly with increasing interatomic distance" [46]. These techniques employ ℓ₁ regularization to construct sparse representations of lattice-dynamical models, avoiding overfitting issues common in linear regression [46].

For infrared-active solids, special consideration must be given to the separation of long-range and short-range components of interatomic forces. A standard approach applies the Ewald summation technique to remove dipole-dipole interactions, ensuring that "the resulting short-range IFCs decay faster than 1/d³" [46].

The following diagram illustrates the advanced workflow for handling anharmonic lattice dynamics:

Benchmarking against prototypical materials remains an indispensable component of rigorous phonon research, particularly in the context of investigating ASR violations and refining phonon frequency calculation methodologies. Through systematic validation using well-characterized systems like MgO and AlN, researchers can verify their computational approaches, identify potential artifacts, and build confidence in their simulation protocols. The continuous development of advanced codes like Pheasy, with capabilities for extracting high-order force constants and treating complex long-range interactions, expands the frontiers of lattice dynamics simulations while simultaneously increasing the importance of robust benchmarking practices. As the field progresses toward increasingly complex materials systems and more sophisticated treatments of anharmonicity, the foundational practice of validation against prototypical materials will remain essential for ensuring the reliability and physical meaningfulness of computational predictions.

The accurate and efficient calculation of phonon properties, such as phonon frequencies and the lattice thermal conductivity, is a cornerstone of modern materials science and drug development, particularly in the design of novel pharmaceutical polymorphs and functional materials. These calculations are essential for understanding thermodynamic stability, phase transitions, and thermal transport. However, a critical challenge in this domain is the violation of the acoustic sum rule (ASR), a physical principle requiring the sum of force constants between an atom and all other atoms to be zero, which ensures the frequencies of acoustic phonons go to zero at the Brillouin zone center. ASR violations, often arising from numerical approximations or the use of finite cutoffs, can lead to unphysical imaginary frequencies and significant inaccuracies in predicted material properties.

Within this research context, this whitepaper provides a comprehensive technical comparison of three predominant computational methods for phonon calculations: Density Functional Perturbation Theory (DFPT), the Finite Displacement (FD) method, and Machine Learning Potentials (MLPs). We evaluate these methods based on their predictive accuracy, computational cost, and specific handling of physical constraints like the ASR, providing researchers with the necessary insights to select the appropriate tool for their investigations.

Theoretical Background and Key Concepts

The Acoustic Sum Rule (ASR) and Its Violation

The ASR is a direct consequence of translational invariance in a crystal lattice. Mathematically, it states that the force constants ( C{ai, bj}(\mathbf{R}m) )—which describe the force in direction i on atom a due to the displacement of atom b in direction j in the unit cell separated by lattice vector ( \mathbf{R}m)—must satisfy: [ \sum{m, b} C{ai, bj}(\mathbf{R}m) = 0 ] for all atoms a and Cartesian directions i, j [18]. This guarantees that the dynamical matrix at the wavevector ( \mathbf{q} = 0 ) has eigenvalues corresponding to three acoustic modes with zero frequency. In practice, approximations in computational methods can break this invariance. For instance, in the FD method, the use of a finite supercell or numerical noise in force calculations can lead to ASR violations. DFPT, while formally exact, can also suffer from ASR violations due to numerical quadrature or the use of a finite plane-wave basis set. Restoring the ASR often requires a post-processing symmetrization step on the force constants, which is crucial for obtaining physically meaningful phonon dispersions, especially for acoustic modes near the zone center.

• Density Functional Perturbation Theory (DFPT): DFPT is an analytical, quantum-mechanical approach that directly calculates the second-order derivative of the total energy with respect to atomic displacements through a self-consistent linear response formalism. It is considered a gold standard for accuracy as it avoids supercell approximations and can, in principle, provide exact force constants within the chosen density functional [57].
• Finite Displacement (FD) Method: Also known as the "small displacement method," the FD approach is a direct numerical technique. It works by displacing each atom in a supercell by a small amount (δ) and using finite differences to approximate the force constants from the resulting Hellmann–Feynman forces [18]. Its simplicity is a key advantage, but its accuracy depends critically on supercell size and the displacement parameter δ.
• Machine Learning Potentials (MLPs): MLPs are fast, data-driven surrogates for the potential energy surface (PES) derived from quantum mechanical methods like Density Functional Theory (DFT). Models such as MACE, EquiformerV2, and CHGNet are trained on large datasets of DFT calculations and can predict energies, forces, and stress tensors with near-DFT accuracy but at a fraction of the computational cost, enabling rapid computation of forces for phonon property prediction [58] [59].

Quantitative Performance Comparison

The table below summarizes the key performance characteristics of DFPT, FD, and MLPs, synthesizing data from recent assessments.

Table 1: Comparative analysis of DFPT, FD, and MLP methods for phonon calculations.

Feature / Metric	Density Functional Perturbation Theory (DFPT)	Finite Displacement (FD) Method	Machine Learning Potentials (MLPs)
Theoretical Foundation	Analytical linear response to atomic displacement [57].	Numerical finite difference of forces from supercell displacements [18].	Data-driven regression on DFT data; e.g., MACE, EquiformerV2 architectures [58] [59].
Typical Force Error	N/A (Formally exact for given functional)	Depends on supercell size and δ; typically higher than DFPT.	~20-50 meV/Å for state-of-the-art models (e.g., EquiformerV2) [58].
Phonon Frequency Error	Benchmark reference.	Increases for smaller supercells and near Γ-point due to ASR violation.	< 0.1 THz for stable modes in high-performing models [58].
ASR Handling	Can be violated by numerical approximations; may require correction.	Inherently violated; requires explicit imposition (e.g., acoustic=True in ASE) [18].	Accuracy depends on training data; ASR not inherently guaranteed, but can be learned.
Computational Cost	Very high; scales poorly with system size (O(N⁴) or worse) [57].	High; scales with the number of atoms (O(N²)) and required displacements.	Very low after training; ~10,000x faster than DFT for force evaluation [60].
System Size Limit	Typically up to ~100 atoms [57].	Limited by supercell DFT calculations (~100s of atoms).	Effectively unlimited; scales linearly with atom count.
Key Advantages	High accuracy; no supercell approximation; direct access to dynamical matrix.	Conceptually simple; easy to implement; works with any electronic structure code.	Extremely fast; enables high-throughput screening; access to anharmonic properties.
Key Limitations	Extreme computational expense; complex implementation; not all DFT codes support it.	ASR violation; computationally demanding for large unit cells; choice of δ is non-trivial.	Requires high-quality training data; potential poor extrapolation to unseen chemistries/structures [59].

A recent large-scale benchmark study of universal MLPs provides specific quantitative data on model performance for phonon properties. The assessment of six models (EquiformerV2, MatterSim, MACE, CHGNet, etc.) on 2,429 materials revealed that EquiformerV2 consistently achieved the highest accuracy in predicting atomic forces, interatomic force constants (IFCs), and ultimately, lattice thermal conductivity [58]. The study highlighted that while models like MACE and CHGNet could achieve comparable force accuracy, notable discrepancies in IFC fitting led to poorer predictions of thermal conductivity, underscoring the complex relationship between force accuracy and derived phonon properties.

Table 2: Performance of select universal MLPs on phonon property prediction (adapted from [58]).

Model	Force Prediction Accuracy	Phonon Frequency Accuracy	Thermal Conductivity (κL) Accuracy
EquiformerV2	Highest	Highest	Most reliable, best agreement with DFT/experiment
MACE	High	Intermediate	Lower than EquiformerV2 due to IFC fitting issues
CHGNet	High	Lower	Poorer prediction
MatterSim	Lower	Intermediate	Intermediate

Detailed Experimental and Computational Protocols

Finite Displacement Method with ASR Correction

The following protocol, based on the ASE package, details the steps for a standard FD phonon calculation, including explicit ASR restoration [18].

Workflow:

Title: Finite Displacement Method Workflow

Step-by-Step Protocol:

• System Initialization: Begin with a fully optimized crystal structure and select an appropriate calculator, such as an Effective Medium Theory (EMT) potential for simple metals or a DFT calculator for first-principles accuracy.
• Phonons Object Setup: Initialize the phonon calculation object. In ASE, this is done using the Phonons class. Critical parameters must be defined:
  • supercell: The size of the supercell used for displacements (e.g., (5,5,5)). Larger supercells minimize the effect of long-range interactions being truncated but increase computational cost.
  • delta: The magnitude of the atomic displacement (e.g., 0.05 Å). This must be small enough to remain in the harmonic regime but large enough to avoid numerical noise.
• Run Calculation: Execute ph.run(). This automatically generates displaced supercells, runs the single-point force calculations for each displacement, and collects the resulting forces.
• Force Constant Matrix Assembly: Read the calculated forces and construct the real-space force constant matrix using ph.read(acoustic=True). The force constant between atom a in the central cell and atom b in cell m is calculated as: ( C{mai, nbj} \approx \frac{F{-} - F{+}}{2 \times \delta} ) where ( F{+} ) and ( F_{-} ) are the forces in direction j on atom nb when atom mai is displaced in the +i and -i directions, respectively [18].
• Impose Acoustic Sum Rule: This is a critical step to correct for ASR violations introduced by the finite supercell. Invoke a method like ph.acoustic() to apply the correction, which typically involves adjusting the force constants to satisfy ( \sum{m, b} C{ai, bj}(\mathbf{R}_m) = 0 ) [18].
• Post-Processing and Analysis: Calculate desired properties.
  • Band Structure: Obtain the phonon dispersion along a high-symmetry path in the Brillouin zone (e.g., path = atoms.cell.bandpath('GXULGK', npoints=100)).
  • Density of States (DOS): Calculate the phonon DOS over a uniform k-point grid (e.g., ph.get_dos(kpts=(20, 20, 20))).

Machine Learning Potential Protocol

This protocol outlines the use of a pre-trained universal MLP for phonon calculations, leveraging the significant speed advantage over direct DFT.

Workflow:

Title: MLP Phonon Calculation Workflow

Step-by-Step Protocol:

• Model Selection: Choose a suitable pre-trained foundational MLP. Current benchmarks indicate EquiformerV2 provides state-of-the-art accuracy for phonon properties, while models like MACE-MP-0 also offer robust, quantitative performance out-of-the-box [58] [61].
• System Setup: Load the atomic structure of the material to be studied.
• Calculator Assignment: Set the chosen MLP as the calculator for the structure. This replaces the DFT calculator in a workflow like the one described in Section 4.1.
• Force and Phonon Calculation:
  • Option 1 (FD with MLP): Use the MLP-powered calculator within the ASE Phonons framework, following the exact same FD protocol in Section 4.1. The key difference is that forces for each displaced supercell are computed nearly instantaneously by the MLP instead of via expensive DFT.
  • Option 2 (Integrated Methods): Some MLP packages or subsequent scripts may offer integrated functionality to compute force constants and phonons directly, potentially using more advanced approaches than the standard FD method.
• Validation (Critical): For systems under unusual strain or with chemistries far from the training set, it is crucial to validate the MLP's predictions against a small number of explicit DFT calculations. This is especially important for properties sensitive to the ASR.

Table 3: Essential computational tools and datasets for phonon calculations.

Category	Item / Resource	Description and Function
Software & Codes	ASE (Atomic Simulation Environment)	A Python library that provides a universal interface for different atomistic codes and includes built-in tools for the Finite Displacement method [18].
	DFT Codes (VASP, Quantum ESPRESSO, ABINIT)	First-principles electronic structure codes that perform the underlying energy and force calculations for FD and DFPT.
	MLP Frameworks (MACE, Equiformer, CHGNet)	Software libraries providing pre-trained models and training code for Machine Learning Potentials [58] [59].
Benchmark Datasets	Materials Project (MP)	A open-access database containing over 154,000 materials with DFT-calculated properties, including dielectric tensors and structures, useful for training and validation [57] [62].
	Open Molecules 2025 (OMol25)	An unprecedented dataset of ~100 million molecular simulation snapshots, designed for training robust MLIPs on chemically diverse systems [60].
	Alexandria Database	A large dataset of DFT calculations extending the Materials Project, often used for training and benchmarking universal MLIPs [59].
Validation Metrics	ASR Compliance Check	A post-calculation check (or correction) ensuring the force constants sum to zero, validating the physical soundness of acoustic modes.
	Phonon Frequency Stability	Analysis of the phonon band structure to check for the absence of significant imaginary frequencies, indicating dynamic instability.

The choice between DFPT, the Finite Displacement method, and Machine Learning Potentials for phonon calculations involves a fundamental trade-off between computational cost, ease of use, and accuracy. DFPT remains the most accurate method for small systems but is computationally prohibitive for high-throughput studies or large complexes. The FD method offers a versatile and widely applicable approach, though its accuracy is contingent on supercell size and requires careful handling of the ASR. MLPs, particularly state-of-the-art models like EquiformerV2, have emerged as a transformative technology, offering near-DFT accuracy at a fraction of the cost, enabling the screening of thousands of materials for phonon-related properties.

For researchers focused on the nuanced challenge of ASR violation, the FD method requires explicit correction, while MLPs rely on the quality and diversity of their training data to learn physically correct force constants. As foundational MLPs continue to be trained on larger and more diverse datasets—such as OMol25 and expanded Materials Project data—their reliability for predicting phonon properties and inherently respecting physical laws like the ASR is expected to become increasingly robust, solidifying their role as an indispensable tool in computational materials science and drug development.

Phonons, the quantized lattice vibrations in crystalline solids, are fundamental determinants of numerous material properties, including thermal conductivity, mechanical behavior, and thermodynamic stability. The accurate computational prediction of these phonon-derived properties hinges on satisfying fundamental quantum-mechanical and symmetry-derived constraints, most notably the Acoustic Sum Rule (ASR). The ASR embodies the translational and rotational invariance of the crystal lattice, requiring that the energy of the system remains unchanged under rigid translations or infinitesimal rotations of the entire crystal [2]. In practical calculations, violations of the ASR frequently occur due to numerical approximations, insufficient basis sets, or inadequate Brillouin zone sampling, leading to unphysical imaginary phonon frequencies and erroneous predictions of dynamic instability. This guide provides an in-depth technical framework for analyzing phonon-derived properties within the critical context of ASR compliance, detailing methodologies for ensuring calculation fidelity and interpreting results for materials research and development.

Theoretical Foundations: Invariance Conditions and Lattice Dynamics

The Born-Huang Invariance Conditions

The formal foundation for the ASR is established by the Born-Huang invariance conditions, which provide the acoustic sum rules for both translational and rotational invariances of a harmonic lattice [2]. These constraints on the first- and second-order interatomic force constants (IFCs) are derived from the requirement of an unchanged lattice potential under global translation or infinitesimal rotation.

For a crystal potential energy expanded in atomic displacements as ( E = E0 + \sum{\varkappa \alpha} \Phi{\varkappa \alpha} u{\varkappa \alpha} + \frac{1}{2} \sum{\varkappa \alpha, \varkappa' \beta} \Phi{\varkappa \alpha, \varkappa' \beta} u{\varkappa \alpha} u{\varkappa' \beta} ), the Born-Huang conditions are [2]:

• Translational Invariance (1st order): ( \sum{\varkappa} \Phi{\varkappa \alpha} = 0 )
• Rotational Invariance (1st order): ( \sum{\varkappa} \Phi{\varkappa \alpha} \tau{\varkappa \beta} = \sum{\varkappa} \Phi{\varkappa \beta} \tau{\varkappa \alpha} )
• Translational Invariance (2nd order): ( \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \beta} = 0 )
• Rotational Invariance (2nd order): ( \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \beta} \tau{\varkappa' \gamma} + \delta{\alpha \gamma} \Phi{\varkappa \beta} = \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \gamma} \tau{\varkappa' \beta} + \delta{\alpha \beta} \Phi{\varkappa \gamma} )

Here, ( \tau{\varkappa \alpha} ) is the equilibrium position of atom ( \varkappa ), and ( \delta{\alpha \beta} ) is the Kronecker delta. When atoms are at their equilibrium positions, the first-order IFCs ( \Phi_{\varkappa \alpha} ) vanish, simplifying the rotational invariance condition for the second-order IFCs [2].

Consequences of ASR Violation in Low-Dimensional Materials

ASR violations have particularly severe consequences in low-dimensional (1D, 2D) materials. A primary manifestation is an incorrect linear dispersion of the flexural (out-of-plane, ZA) acoustic phonon mode in 2D materials near the Brillouin zone center (Γ-point), instead of the physically correct quadratic dispersion [2]. This error fundamentally misrepresents the vibrational density of states at low frequencies, leading to inaccurate predictions of thermal conductivity, thermodynamic stability, and carrier transport properties. Ensuring rotational invariance is, therefore, not merely a numerical formality but a critical prerequisite for obtaining physically meaningful phonon spectra in low-dimensional systems.

Table 1: Key Invariance Conditions and Their Physical Significance

Condition Type	Mathematical Expression	Physical Significance	Consequence of Violation
Translational Invariance	( \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \beta} = 0 )	Conservation of total crystal momentum	Non-zero forces from rigid translation; unphysical phonon gaps at Γ-point
Rotational Invariance	( \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \beta} \tau{\varkappa' \gamma} = \sum{\varkappa'} \Phi{\varkappa \alpha, \varkappa' \gamma} \tau{\varkappa' \beta} )	Conservation of total angular momentum	Incorrect ZA phonon dispersion in 2D materials; erroneous thermal properties

Computational Methodologies for Phonon Property Analysis

First-Principles Phonon Calculations: DFPT and Finite-Displacement Methods

Two primary first-principles methodologies are employed for phonon calculations: Density Functional Perturbation Theory (DFPT) and the finite-displacement method.

• Density Functional Perturbation Theory (DFPT): This approach directly calculates the second-order IFCs in reciprocal space by computing the linear response of the electron density to a lattice perturbation. It is highly efficient for polar materials as it naturally incorporates non-analytical term corrections for proper treatment of the longitudinal optical-transverse optical (LO-TO) splitting. In VASP, this is activated with IBRION = 7 or 8 [4].
• Finite-Displacement Method: This real-space method constructs the dynamical matrix by calculating Hellmann-Feynman forces resulting from finite displacements of atoms in a supercell. The force constants are then determined via ( \Phi{\varkappa \alpha, \varkappa' \beta} = - \frac{\partial F{\varkappa' \beta}}{\partial u_{\varkappa \alpha}} ). This method is more computationally expensive, especially for large unit cells, and requires careful convergence with supercell size. In VASP, the tag IBRION = 5 or 6 is used [4].

Protocol for Robust Phonon Dispersion and DOS Calculation

A robust workflow for phonon property calculation must explicitly address ASR compliance [4]:

• Force Constant Calculation: Perform a supercell calculation using either DFPT (IBRION = 7, 8) or the finite-displacement method (IBRION = 5, 6). The supercell must be large enough for force constants to decay at the boundaries.
• ASR Enforcement: Explicitly impose the translational and rotational invariance conditions on the calculated IFCs in post-processing. Most computational packages offer options to enforce these sum rules (e.g., LPHON_READ_FORCE_CONSTANTS = True in VASP to read and correct pre-calculated constants) [4].
• Phonon Dispersion: Compute the phonon dispersion along a high-symmetry path in the Brillouin zone. A QPOINTS file specifying the path is required, and the calculation is activated with LPHON_DISPERSION = .TRUE. [4].
• Phonon DOS: Compute the phonon density of states on a dense, uniform q-point mesh specified in the QPOINTS file, using PHON_DOS > 0 and setting the broadening with PHON_SIGMA [4].

For polar materials, the long-range dipole-dipole interactions must be treated via Ewald summation by setting LPHON_POLAR = True and providing the static dielectric tensor (PHON_DIELECTRIC) and Born effective charges (PHON_BORN_CHARGES) obtained from a separate DFPT calculation (LEPSILON = .TRUE.) [4].

Machine Learning Accelerated Approaches

Machine learning interatomic potentials (MLIPs) have emerged as a powerful tool to accelerate high-throughput phonon calculations. The strategy involves training models like MACE (Message Passing with Equivariant Embeddings) on a diverse dataset of supercell structures with atomic displacements and the corresponding DFT-calculated forces [5]. Once trained, these potentials can predict forces for new structures, bypassing the need for expensive DFT calculations in the finite-displacement method. This approach can reduce the number of required supercell calculations by using a subset of structures with all atoms randomly perturbed, significantly accelerating the screening of phonon properties for dynamic and thermodynamic stability across vast material spaces [5].

Figure 1: Computational workflow for first-principles phonon calculations, highlighting critical decision points for method selection and ASR correction.

Advanced Phonon Phenomena and Derived Properties

Phonon Thermodynamics and Stability

Phonon spectra are the direct input for calculating the vibrational contributions to thermodynamic properties. The Helmholtz free energy ( F(T) ), entropy ( S(T) ), and heat capacity ( C_V(T) ) are obtained from the phonon density of states ( g(\omega) ) as [5]:

• Helmholtz Free Energy: ( F{vib}(T) = kB T \int0^\infty \ln\left[2 \sinh\left(\frac{\hbar \omega}{2 kB T}\right)\right] g(\omega) d\omega )
• Entropy: ( S{vib}(T) = kB \int0^\infty \left[ \frac{\hbar \omega / kB T}{e^{\hbar \omega / kB T} - 1} - \ln\left(1 - e^{-\hbar \omega / kB T}\right) \right] g(\omega) d\omega )
• Heat Capacity: ( C{vib}(T) = kB \int0^\infty \left(\frac{\hbar \omega}{kB T}\right)^2 \frac{e^{\hbar \omega / kB T}}{(e^{\hbar \omega / kB T} - 1)^2} g(\omega) d\omega )

A fundamental application is assessing dynamic stability: a material is dynamically stable if all its phonon frequencies are real across the entire Brillouin zone. The presence of imaginary frequencies (often indicated as negative values in computational outputs) signifies an unstable lattice that would spontaneously distort.

Phonon Transport and Thermal Conductivity

Phonons are the primary heat carriers in semiconductors and insulators. The lattice thermal conductivity ( \kappa_p ) is a key property for thermoelectric applications and thermal management, calculated by solving the Boltzmann Transport Equation (BTE). Advanced approaches can even measure the entropy dynamics associated with nonequilibrium phonon transport in mesoscopic systems [63]. In the relaxation time approximation, the conductivity is given by:

( \kappap = \frac{1}{3} \sum{\lambda} \int C{ph}(\omega\lambda(\mathbf{q})) \, v{g,\lambda}^2(\mathbf{q}) \, \tau{\lambda}(\mathbf{q}) \, d\mathbf{q} )

where ( \lambda ) is the phonon mode index, ( C{ph} ) is the spectral heat capacity, ( vg ) is the group velocity, and ( \tau ) is the phonon lifetime. The group velocity ( vg = \nabla\mathbf{q} \omega(\mathbf{q}) ) is directly derived from the phonon dispersion, underscoring the critical need for an accurate, ASR-compliant spectrum. For thermoelectric materials, a common design strategy is hierarchical nanostructuring to introduce scattering centers across a wide range of phonon mean-free-paths, thereby significantly reducing ( \kappa_p ) without severely compromising electronic transport [64].

Topological Phononics and Berry Curvature

Beyond classical transport, phonons can exhibit non-trivial topological properties. In crystals where time-reversal or inversion symmetry is broken, phonons can acquire a Berry curvature in momentum space [17]. This curvature can influence phonon dynamics, leading to phenomena such as a phonon thermal Hall effect. A novel mechanism has been identified in Dirac materials like BaMnSb₂, where phonons can "inherit" a Berry curvature through their coupling to topologically non-trivial electrons, even when time-reversal symmetry is preserved [17]. This inherited curvature is proportional to the electronic valley Chern number and contributes to the thermal Hall conductivity, representing an active frontier in phonon research.

Figure 2: Relationship between the foundational phonon dispersion and key derived properties analyzed in this guide.

Table 2: Key Computational Tools and Resources for Phonon Research

Tool/Resource	Type	Primary Function	Relevance to ASR/Stability
VASP [4]	Software Package	First-principles DFT/DFPT calculations	Calculates force constants; requires post-processing for ASR enforcement.
Phonopy	Software Package	Post-processing force constants	Extracts phonon dispersion/DOS; can enforce ASR on IFCs.
ShengBTE [65]	Software Package	Solves BTE for lattice thermal conductivity	Requires accurate, stable phonon dispersions as input.
ALIGNN [5]	Machine Learning Model	Direct prediction of phonon density of states	Bypasses force constant calculation; trained on ASR-compliant data.
MACE Potentials [5]	Machine Learning Interatomic Potential	Accelerated force and energy prediction	Enables high-throughput phonon screening with DFT accuracy.
Materials Project [5]	Database	Repository of calculated material properties	Source for phonon data and structures for validation.

The reliable analysis of phonon-derived properties is intrinsically linked to the rigorous adherence to the Acoustic Sum Rule. ASR violations represent a critical source of error, particularly for low-dimensional materials where they corrupt the fundamental physics of flexural phonons. By implementing robust computational protocols that enforce translational and rotational invariance, researchers can ensure the predictive accuracy of phonon dispersions, which form the basis for evaluating thermodynamic stability, thermal transport, and emerging topological phenomena. The integration of machine learning potentials promises to further accelerate this field, enabling the high-throughput discovery and design of materials with tailored vibrational properties for applications ranging from thermoelectrics to nanoelectronics.

Assessing the Impact of ASR Enforcement on Predicted Material Properties

The accurate prediction of material properties is a cornerstone of computational materials science and drug development. In the specific domain of lattice dynamics, the Acoustic Sum Rule (ASR) is a fundamental physical principle arising from translational invariance, which dictates that the frequencies of the three acoustic phonons at the gamma point (Γ) must be zero [3]. In practical computations, however, approximations and numerical discretization often lead to the violation of this rule, a phenomenon known as ASR violation. This violation manifests as non-zero frequencies for acoustic modes, which can be a significant source of error, compromising the predictive fidelity of subsequent material properties such as thermodynamic stability, thermal conductivity, and vibrational spectra [3] [66]. This technical guide examines the profound impact of ASR enforcement on the prediction of material properties, situating the discussion within a broader research context focused on managing ASR violations and ensuring the robustness of phonon frequency calculations.

Understanding ASR Violations in Phonon Calculations

Physical Origins and Computational Manifestations

The Acoustic Sum Rule is a direct consequence of the invariance of the total energy of a crystal under rigid translations. In an ideal, perfectly converged calculation, the dynamical matrix would inherently satisfy this condition. However, in real-world simulations using packages like Quantum ESPRESSO, several factors can lead to its violation. The primary and most irreducible cause is the discreteness of the Fast Fourier Transform (FFT) grid used in Plane-Wave (PW) calculations [3]. This discretization breaks the perfect translational symmetry at a numerical level.

Other prevalent causes include:

• Insufficient Convergence Parameters: Inadequate settings for the self-consistent field (conv_thr) or phonon (tr2_ph) calculations can prevent the dynamical matrix from fully converging, exacerbating ASR violations [3] [66].
• Finite Supercell Effects: For molecules, non-zero frequencies for rotational modes and acoustic modes are a fictitious effect of the finite supercell size used in the calculation [3].
• Incomplete Structural Relaxation: A structure that has not been fully relaxed to its ground state, with residual forces or stresses, will inherently possess an unstable lattice, leading to significant ASR violations and often imaginary phonon frequencies [66]. Insufficiently tight convergence thresholds for total energy (etot_conv_thr) and forces (forc_conv_thr) during the initial structural relaxation are a common root cause.

A key manifestation of underlying problems in the calculation is the appearance of imaginary frequencies (often reported as "negative" frequencies, as they correspond to negative values of ω²). While small imaginary frequencies can sometimes be a convergence artifact, their presence, particularly away from the Gamma point, can also signal a genuine material instability [3].

Consequences for Predicted Material Properties

Violations of the ASR introduce systematic errors that propagate through the calculation, affecting a wide range of derived material properties:

• Thermodynamic Properties: Properties like entropy, free energy, and heat capacity, which are obtained by integrating over the entire phonon spectrum, become inaccurate. The erroneous low-frequency acoustic modes significantly skew the results in the low-temperature regime.
• Thermal Conductivity: Lattice thermal conductivity, often computed from phonon lifetimes and group velocities, is highly sensitive to the low-frequency part of the spectrum. ASR violations directly compromise the accuracy of such predictions.
• Phase Stability: The assessment of a material's dynamical stability relies on the absence of imaginary frequencies. Spurious imaginary frequencies caused by poor convergence and unresolved ASR violations can falsely indicate an unstable structure, hindering the accurate prediction of phase diagrams.

Methodologies for ASR Enforcement and Analysis

Protocols for ASR Enforcement in Quantum ESPRESSO

The following workflow provides a detailed protocol for addressing ASR violations in a standard phonon calculation using the Quantum ESPRESSO suite.

Diagram 1: ASR Enforcement Workflow

Step 1: Rigorous Structural Relaxation

• Objective: Ensure the initial crystal structure is in its true ground state with negligible residual forces.
• Protocol: Perform a vc-relax calculation with tightly converged parameters.
  • forc_conv_thr = 1.0d-10 or lower (e.g., 1.0d-12) [66].
  • etot_conv_thr = 1.0d-10 or lower.
  • Verify that the absolute values of all residual forces are below 0.0001 eV/Å.

Step 2: Self-Consistent Field (SCF) and Phonon Calculations

• Objective: Generate a highly accurate ground state and linear response.
• Protocol:
  • In the SCF calculation, use a high kinetic energy cutoff and dense k-point grid. Avoid smearing for insulators if macroscopic electric field contributions are to be calculated [3].
  • In the ph.x input, set a stringent convergence threshold: tr2_ph = 1.0d-14 [66].

Step 3: Post-Processing with ASR Enforcement

• Objective: Impose the ASR on the calculated force constants.
• Protocol: In the q2r.x input file, use the zasr variable to apply a correction.
  • zasr = 'crystal': Applies the sum rule based on the crystal symmetry. This is the recommended option for most crystalline solids [66].
  • zasr = 'simple': Applies a simple sum rule on all atoms.
  • Ensure the same asr = 'crystal' setting is used in the subsequent matdyn.x calculation for consistency.

Advanced and Machine Learning Approaches

Beyond direct enforcement in force constants, advanced computational frameworks are being developed to improve property prediction, which can indirectly benefit from robust phonon data.

Bilinear Transduction for Property Extrapolation: A transductive machine learning method has been proposed for out-of-distribution (OOD) property prediction. This approach reparameterizes the prediction problem by learning how property values change as a function of material differences in representation space, rather than predicting absolute values. This method has been shown to improve extrapolative precision for materials by 1.8× and boost the recall of high-performing candidates by up to 3×, demonstrating its utility in identifying materials with exceptional properties [67].

Dual-Stream Graph Neural Networks (GNNs): To overcome the limitations of models that only capture topological information, a dual-stream architecture fusing both topological and spatial atomic information has been developed. This model uses a message-passing GNN for topological information and a spatial encoder that processes 3D atomic coordinates, with features fused via an attention mechanism. This approach has demonstrated superior performance in predicting formation energies, a key property for assessing stability [68].

Quantitative Impact of ASR Enforcement

The enforcement of the ASR has a quantifiable and significant impact on the numerical outcomes of phonon calculations and the subsequent prediction of material properties. The tables below summarize key quantitative findings from the literature.

Table 1: Impact of ASR Enforcement on Phonon Frequency Calculations

Material System	Without ASR Enforcement (Acoustic Mode Freq., cm⁻¹)	With ASR Enforcement (Acoustic Mode Freq., cm⁻¹)	Computational Setup
Rutile (TiO₂)	Imaginary/negative frequencies in dispersion [66]	Stable, positive frequencies achieved	zasr='crystal' in q2r.x & matdyn.x [66]
General Crystals	Small, non-zero frequencies at Gamma [3]	Frequencies reduced to zero	ASR imposed on dynamical matrix [3]
Molecules in Supercell	Non-zero rotational/acoustic modes [3]	Frequencies reduced to zero	Corrected via ASR for finite-size effects [3]

Table 2: Performance of Advanced Prediction Models Leveraging Structural Data

Prediction Model	Key Methodology	Performance Gain	Application / Property Predicted
Bilinear Transduction [67]	Transductive learning for OOD extrapolation	1.8x extrapolative precision, 3x recall of top candidates	Multi-property prediction across AFLOW, Matbench
Dual-Stream GNN [68]	Fuses topological and spatial atomic information	Superior performance vs. CGCNN, MEGNet	Formation energy (Materials Project database)
ElaTBot-DFT [69]	Domain-specific Large Language Model (LLM)	33.1% error reduction vs. Darwin LLM	Elastic constant tensors, bulk modulus

Table 3: Key Software and Computational Resources for Phonon and Property Prediction

Item Name	Function / Role in Workflow	Relevant Context
Quantum ESPRESSO [3] [66]	Integrated suite for SCF, structural relaxation, and phonon calculations. The ph.x, q2r.x, and matdyn.x modules are central.	Primary engine for first-principles phonon computation and ASR enforcement.
phonopy [11]	An open-source package for performing first-principles phonon calculations, widely used in materials science.	An alternative to Quantum ESPRESSO for calculating and analyzing phonon properties.
ISOTROPY/ACKJ/ACMI [3]	External packages for symmetry analysis of phonons and normal modes.	Used for determining the symmetry of phonon modes, complementing the internal analyzer in Quantum ESPRESSO.
MatEx [67]	Open-source implementation (matex) of the Bilinear Transduction method for OOD property prediction.	For machine learning-based screening and discovery of materials with extreme properties.
TSGNN [68]	Code for the dual-stream GNN model that fuses topological and spatial information.	For enhanced property prediction where 3D atomic arrangement is critical.

The enforcement of the Acoustic Sum Rule is not merely a numerical formality but a critical step for ensuring the physical correctness and predictive accuracy of phonon calculations. As demonstrated, failure to properly address ASR violations leads to qualitatively and quantitatively erroneous predictions for fundamental material properties, from thermodynamic behavior to phase stability. The methodologies outlined—ranging from stringent convergence protocols and built-in correction tools in established software suites to emerging machine learning models—provide a comprehensive toolkit for researchers to mitigate these errors. The quantitative improvements afforded by ASR enforcement and advanced prediction models underscore their indispensable role in the robust computational design and discovery of new materials and molecules.

Conclusion

The accurate handling of the Acoustic Sum Rule is not merely a technical detail but a fundamental requirement for reliable phonon calculations. As demonstrated, ASR violations often stem from computational approximations rather than physical instabilities, and can be systematically addressed through proper convergence, symmetry handling, and the application of corrective schemes. The emergence of machine learning potentials offers a promising path for high-throughput phonon property screening, yet their validation against established DFT methods remains crucial. For researchers in biomaterials and drug development, mastering these concepts is increasingly important for predicting the stability and thermal properties of complex materials, from metal-organic frameworks to pharmaceutical crystals. Future directions will likely involve more automated ASR compliance in code bases and the integration of robust ML potentials into standard materials discovery workflows, ultimately accelerating the design of novel materials with tailored dynamic properties.

References

• 1. Lab 6: Vibrational Frequencies and Normal Modes [https://mse404.gitlab.io/labs/lab06]
• 2. General invariance and equilibrium conditions for lattice ... [https://www.nature.com/articles/s41524-022-00920-6]
• 3. phonons [https://www.quantum-espresso.org/faq/faq-phonons/]
• 4. Computing the phonon dispersion and DOS [https://www.vasp.at/wiki/Computing_the_phonon_dispersion_and_DOS]
• 5. Accelerating high-throughput phonon calculations via ... [https://www.sciencedirect.com/science/article/abs/pii/S2542529325000446]
• 6. 2 Running phonon calculations [https://www.tcm.phy.cam.ac.uk/castep/Phonons_Guide/2-sec:examples.html]
• 7. High-throughput density-functional perturbation theory ... [https://www.nature.com/articles/sdata201865]
• 8. Phonon Dispersion [https://docs.materialsproject.org/methodology/materials-methodology/phonon-dispersion]
• 9. Resonant Exciton Scattering Reveals Raman Forbidden ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10165653/]
• 10. Phonon modes of MgB 2 : super-lattice structures and ... [https://pubs.rsc.org/en/content/articlehtml/2014/cp/c4cp03449j]
• 11. First principles phonon calculations in materials science [https://www.sciencedirect.com/science/article/pii/S1359646215003127]
• 12. Fast and accurate machine learning prediction of phonon ... [https://www.nature.com/articles/s41524-023-01020-9]
• 13. Fast and accessible phonon calculations [https://www.mat3ra.com/news-and-blog-posts/fast-and-accessible-phonon-calculations]
• 14. Role of long-range dipolar interactions in the simulation ... [https://arxiv.org/html/2506.22337v1]
• 15. Vibrational Spectroscopy — AMS 2025.1 documentation - SCM [https://www.scm.com/doc/AMS/Vibrational_Spectroscopy.html]
• 16. A database of computed Raman spectra of inorganic ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC10803741/]
• 17. Inherited Berry curvature of phonons in Dirac materials with ... [https://arxiv.org/html/2502.13814v2]
• 18. Phonon calculations — ASE documentation [https://wiki.fysik.dtu.dk/ase/ase/phonons.html]
• 19. Phonons from finite differences [https://vasp.at/wiki/Phonons_from_finite_differences]
• 20. Phonons from density-functional-perturbation theory [https://vasp.at/wiki/Phonons_from_density-functional-perturbation_theory]
• 21. First tutorial on DFPT [https://docs.abinit.org/tutorial/rf1/]
• 22. DFPT [https://docs.abinit.org/topics/DFPT/]
• 23. density functional theory - DFPT vs Finite Difference Method [https://mattermodeling.stackexchange.com/questions/2000/dfpt-vs-finite-difference-method]
• 24. Phonon Calculations via VASP - Danny Rehn by rehnd [https://rehnd.github.io/tutorials/vasp/phonons]
• 25. IFC_ASR - VASP Wiki [https://vasp.at/wiki/IFC_ASR]
• 26. anaddb input variables [https://docs.abinit.org/variables/anaddb/]
• 27. symfc/symfc: Generate symmetrized force constants [https://github.com/symfc/symfc]
• 28. Imaginary frequencies/ ASR - QuantumATK Forum [https://forum.quantumatk.com/index.php?topic=3635.0]
• 29. “One defect, one potential” strategy for accurate machine ... [https://arxiv.org/html/2509.00498v1]
• 30. Machine learning force field based phonon dispersion ... [https://www.sciencedirect.com/science/article/pii/S1567173924001500]
• 31. Workflow — Phono3py v.3.19.4 [https://phonopy.github.io/phono3py/workflow.html]
• 32. Phonon calculation for NaCl using Phonopy and Quantum ... [https://ma.issp.u-tokyo.ac.jp/en/app-post/7148]
• 33. Quantum ESPRESSO (QE) & phonopy calculation [https://phonopy.github.io/phonopy/qe.html]
• 34. CASTEP & phonopy calculation [https://phonopy.github.io/phonopy/castep.html]
• 35. Atomate2: modular workflows for materials science [https://pubs.rsc.org/en/content/articlehtml/2025/dd/d5dd00019j]
• 36. Machine learned potential for high-throughput phonon ... [https://www.nature.com/articles/s41524-025-01611-8]
• 37. Accelerating point defect photo-emission calculations with ... [https://www.nature.com/articles/s41524-025-01820-1]
• 38. MACE-MP-0: ML Potential for Atomistic Simulations [https://www.emergentmind.com/topics/mace-mp-0]
• 39. Machine Learned Potential for High-Throughput Phonon ... [https://arxiv.org/abs/2412.02877]
• 40. "Accelerating High-Throughput Phonon Calculations via ... [https://pdxscholar.library.pdx.edu/mengin_fac/521/]
• 41. A Comprehensive Assessment and Benchmark Study of ... [https://arxiv.org/abs/2509.03401]
• 42. 6 Troubleshooting [https://www.quantum-espresso.org/Doc/ph_user_guide/node18.html]
• 43. Phonon calculations — ASE documentation [https://ase-lib.org/ase/phonons.html]
• 44. ph.x: input description [https://www.quantum-espresso.org/Doc/INPUT_PH.html]
• 45. Electronic, Optical, and Elastic Properties of CaFI ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9096925/]
• 46. First-principles phonon physics using the Pheasy code [https://arxiv.org/html/2508.01020v1]
• 47. Ab Initio Theory of Phonon Magnetic Moment Induced by ... [https://arxiv.org/html/2503.10160v4]
• 48. eph for variational polaron equations [https://docs.abinit.org/tutorial/eph4vpq/]
• 49. An automated protocol to construct flexibility parameters for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11258866/]
• 50. How to enhance the accuracy of the Phonon calculations? [https://matsci.org/t/how-to-enhance-the-accuracy-of-the-phonon-calculations/64476]
• 51. Machine-learning-assisted discovery of lattice dynamics ... [https://pubs.rsc.org/en/content/articlehtml/2025/mh/d5mh01176k]
• 52. pw.x: input description [https://www.quantum-espresso.org/Doc/INPUT_PW.html]
• 53. How to dive between different input geometries in Quantum ... [https://mattermodeling.stackexchange.com/questions/14234/how-to-dive-between-different-input-geometries-in-quantum-espresso]
• 54. [Pw_forum] Symmetry Related Error [https://pw-forum.pwscf.narkive.com/rgl4xFEB/symmetry-related-error]
• 55. Perturbo: A software package for ab initio electron–phonon ... [https://www.sciencedirect.com/science/article/abs/pii/S0010465521000837]
• 56. Perturbo: A software package for ab initio electron–phonon ... [https://authors.library.caltech.edu/records/w5ppg-n5n92/latest]
• 57. Discovery of highly anisotropic dielectric crystals with ... [https://pubs.rsc.org/en/content/articlehtml/2025/fd/d4fd00096j]
• 58. Researchers Accurately Predict Phonon Properties In ... [https://quantumzeitgeist.com/429-ml-accurately-predict-phonon-properties-materials-universal-potentials/]
• 59. Universal Machine Learning Potentials under Pressure [https://arxiv.org/html/2508.17792v1]
• 60. A Record-Breaking Dataset to Train AI Models has Launched [https://newscenter.lbl.gov/2025/05/14/computational-chemistry-unlocked-a-record-breaking-dataset-to-train-ai-models-has-launched/]
• 61. Data-efficient fine-tuning of foundational models for first ... [https://pubs.rsc.org/en/content/articlehtml/2025/fd/d4fd00107a]
• 62. Recent advances and applications of deep learning ... [https://www.nature.com/articles/s41524-022-00734-6]
• 63. Entropy measurement of nonequilibrium phonon heat ... [https://www.nature.com/articles/s42005-025-02187-x]
• 64. Phonon and electron transport engineering for enhanced ... [https://www.oaepublish.com/articles/energymater.2025.32]
• 65. Thermal Management and Phonon Transport in BAs/WSe ... [https://www.sciencedirect.com/science/article/abs/pii/S2542529325001993]
• 66. How to get rid of this negative frequencies? [https://mattermodeling.stackexchange.com/questions/13277/how-to-get-rid-of-this-negative-frequencies]
• 67. Out-of-Distribution Property Prediction in Materials and ... [https://www.nature.com/articles/s41524-025-01808-x]
• 68. Enhancing material property prediction via a dual stream ... [https://www.sciencedirect.com/science/article/abs/pii/S0927025625000369]
• 69. elastic constant tensor prediction and materials design [https://pubs.rsc.org/en/content/articlelanding/2025/dd/d5dd00061k]