Solving Negative Phonon Frequencies in Quantum ESPRESSO: A Comprehensive Troubleshooting Guide

Abstract

This guide provides a systematic approach for researchers and computational materials scientists to diagnose and resolve the common issue of negative (imaginary) phonon frequencies in Quantum ESPRESSO ph.x calculations. It covers foundational concepts explaining what negative frequencies signify, detailed methodological workflows for robust calculations, targeted troubleshooting strategies for numerical and physical instabilities, and validation techniques to confirm result reliability. The article synthesizes official documentation, community expertise, and practical examples to offer a complete solution framework applicable to various material systems, enabling users to distinguish between numerical artifacts and genuine physical instabilities in their simulations.

Understanding Negative Phonon Frequencies: From Numerical Artifacts to Physical Instabilities

Conceptual Foundations: What Are Negative Frequencies?

The Mathematical Meaning of "Imaginary" Frequencies

In computational physics, particularly in phonon calculations, a negative frequency does not imply a frequency value less than zero in the ordinary sense. Instead, it signifies an imaginary frequency. Mathematically, if the eigenvalue of the dynamical matrix is negative, the resulting frequency ( \omega ) is calculated as ( \omega = \sqrt{\lambda} ). When ( \lambda ) is negative, the frequency becomes an imaginary number ( \omega = i\sqrt{|\lambda|} ). In the output of codes like Quantum ESPRESSO, this is often reported as a "negative" value for convenience, but it physically represents an imaginary mode [1].

The Physical Interpretation: Instability

Imaginary frequencies reveal a dynamic instability in the system's structure [1]. A stable structure at a local energy minimum should have only real, positive vibrational frequencies. An imaginary frequency corresponds to a mode—a specific pattern of atomic displacement—along which the system's potential energy decreases, not increases. If the atoms are displaced along the direction of this mode, they will not oscillate back to their starting positions but will instead move towards a different, more stable configuration [1]. In the context of transition state theory, a single imaginary frequency is the hallmark of a saddle point on the potential energy surface, which is the unstable point separating reactants from products [1]. The presence of multiple imaginary frequencies can indicate a high-dimensional saddle point connecting multiple reaction pathways [1].

Troubleshooting Guide: Resolving Negative Frequencies in Quantum ESPRESSO

This guide addresses the most common causes and solutions for unphysical imaginary frequencies.

Inadequate Structural Relaxation

The most frequent cause of unexpected imaginary frequencies is an insufficiently relaxed structure [2] [3]. The phonon calculation must be performed on a structure that is at a true local minimum of the potential energy surface.

• Problem: The forces on atoms and/or the stress on the unit cell have not been reduced to a sufficiently small value during the initial vc-relax or relax calculation. The structure is still "feeling" a slope in the energy landscape, leading to instabilities [3].
• Symptoms: Multiple imaginary frequencies, often with large magnitudes, appearing across the Brillouin zone [4] [2].
• Solutions:
  • Tighten Convergence Criteria: In your pw.x relaxation input, use stricter convergence thresholds.

  • Verify Final Forces: Always check the output of the relaxation calculation to ensure the final forces on all atoms are very close to zero.

Inconsistent Computational Parameters

The phonon calculation (ph.x) must be performed on exactly the same potential energy surface as the self-consistent field (SCF) calculation that produced the ground state. Using different parameters creates an inconsistency [3].

• Problem: Parameters like the pseudo_dir, k-point grid, exchange-correlation functional (LDA/GGA), or vdw_corr differ between the pw.x (SCF) and ph.x calculations [3].
• Symptoms: Imaginary frequencies, especially at the gamma point, that should not be present.
• Solutions:
  • Parameter Consistency: Double-check that all parameters in the &SYSTEM and &ELECTRONS namelists of your SCF calculation are identical to those used to generate the charge density read by ph.x.
  • Functional Choice: Be aware that some structures are known to be unstable with certain functionals (e.g., rutile with GGA) and may require LDA or a hybrid functional [4].

Insufficient Phonon Calculation Convergence

The phonon calculation itself may not be fully self-consistent.

• Problem: The tr2_ph threshold in the ph.x input is too loose, or the k-point grid for the phonons (nk1, nk2, nk3) is too coarse [2] [3].
• Symptoms: Small, spurious imaginary frequencies, or inconsistencies between ph.x and matdyn.x results [4].
• Solutions:
  • Tighten Phonon Convergence: In your ph.x input, use a lower threshold.

  • Use a Denser q-Point Grid: Ensure the nq1, nq2, nq3 grid in ph.x is dense enough for your system [2].

Acoustic Sum Rule (ASR) Violations

At the gamma point (q=0), the three acoustic modes should have frequencies exactly equal to zero, corresponding to translational invariance. Finite-size effects and numerical approximations can cause these to be slightly non-zero, sometimes appearing as small imaginary frequencies [5].

• Problem: The calculated frequency of acoustic modes at Γ is not zero [5].
• Symptoms: Small imaginary frequencies (typically < 10 cm⁻¹, but sometimes higher) specifically at the gamma point [5].
• Solutions:
  • Impose the ASR: When running q2r.x and matdyn.x, use the zasr and asr flags to impose the Acoustic Sum Rule.

  • The Ultimate Test: As recommended in the Quantum ESPRESSO guide, the definitive test is to diagonalize the dynamical matrix using dynmat.x with the ASR imposed. If the acoustic mode then has a very small frequency (e.g., < 1 cm⁻¹) and other modes are unchanged, you can trust your results [5].

The relationship between these primary issues and their solutions can be visualized in the following troubleshooting workflow:

Frequently Asked Questions (FAQs)

My ph.x calculation had only positive frequencies, but matdyn.x shows imaginary ones. Why?

This discrepancy often arises from the interpolation process [4]. ph.x calculates the dynamical matrix exactly at a specific set of q-points. q2r.x then uses these to create a force constant matrix in real space, and matdyn.x interpolates this matrix to get frequencies at any q-point. If the original ph.x q-point grid (nq1, nq2, nq3) is too coarse, the interpolation can be inaccurate and introduce imaginary frequencies at q-points not explicitly calculated by ph.x [4]. The solution is to use a denser grid in your ph.x calculation.

I have one imaginary frequency at my saddle point. Is this correct?

Yes, this is expected and correct. In transition state theory, a saddle point on the potential energy surface is characterized by exactly one imaginary frequency [1]. The eigenvector of this mode points along the reaction coordinate—the path atoms take to convert from reactants to products.

What does it mean if I have two or more imaginary frequencies?

Multiple imaginary frequencies at a single point indicate a higher-order saddle point [1]. This means there is more than one direction in which the energy decreases. In the context of a reaction path, this might mean that the same saddle point connects to more than two potential energy minima (e.g., multiple reaction pathways) [1]. It may also suggest that an even lower-energy saddle point exists nearby, and you should follow one of the imaginary modes to find it [1].

How do I know if my imaginary frequencies are a physical instability or a numerical error?

The following table helps distinguish between physical and numerical instabilities:

Feature	Physical Instability	Numerical Instability
Magnitude	Often large (e.g., > 50 cm⁻¹) [3]	Often small (e.g., < 20 cm⁻¹)
Location	May occur at multiple q-points	Often only at gamma (Γ) or a few points [5]
Persistence	Remains after tightening convergence	Disappears or shrinks with better convergence [4]
ASR Test	Unaffected by imposing the ASR	Often fixed by imposing the ASR [5]

Essential Parameters for Stable Phonon Calculations

The following table summarizes key parameters in Quantum ESPRESSo that are critical for obtaining accurate, physical phonon spectra.

Table 1: Critical Parameters for Phonon Calculations

Calculation Step	Parameter	Description & Function	Recommended Value
pw.x (Relax/SCF)	forc_conv_thr	Force convergence threshold; crucial for finding a true minimum.	1.0d-5 Ry/au or tighter [2]
	etot_conv_thr	Total energy convergence threshold.	1.0d-6 Ry or tighter [2]
	conv_thr	SCF cycle convergence for electron density.	1.0d-12 or tighter [2]
	ecutwfc / ecutrho	Plane-wave kinetic energy cutoffs.	Must be high enough for convergence [5]
ph.x	tr2_ph	Phonon SCF convergence threshold.	1.0d-14 or tighter [4] [2]
	nq1, nq2, nq3	Grid of q-points for phonon calculation.	System-dependent; must be converged [2]
q2r.x	zasr	Applies Acoustic Sum Rule to force constants.	'crystal' [4] [2]
matdyn.x	asr	Applies Acoustic Sum Rule during interpolation.	'crystal' [4] [2]

A comprehensive guide to diagnosing the root causes of negative frequencies in your phonon calculations

Negative frequencies in Quantum ESPRESSO ph.x calculations represent imaginary phonon modes that can either signal a computational problem or reveal genuine physical instability in your material's structure. Correctly diagnosing the origin is crucial for producing reliable research. This guide provides a systematic approach to troubleshooting these issues.

Understanding Negative Frequencies: Artifact vs. Reality

Imaginary frequencies (reported as negative values in phonon dispersion plots) indicate that the atomic configuration is not at a true energy minimum. The critical distinction researchers must make is whether these represent:

• Computational artifacts from inadequate convergence parameters or methodological errors
• Genuine structural instabilities where the current structure is mechanically unstable

The acoustic sum rule (ASR) violation typically results in small negative frequencies (< 10 cm⁻¹) for acoustic modes at Γ-point, which are usually artifacts. In contrast, large imaginary frequencies (> 20-50 cm⁻¹) often suggest genuine instability, particularly if they persist after thorough convergence testing [5].

Diagnostic Framework: Systematic Troubleshooting Approach

Follow this logical workflow to identify the root cause of negative frequencies in your calculations:

Computational Artifacts: Identification and Solutions

Inadequate Convergence Parameters

The most common source of spurious negative frequencies stems from insufficiently tight convergence criteria across multiple parameters:

Table 1: Key Convergence Parameters and Recommended Values

Parameter	Typical Problematic Value	Recommended Value	Purpose
tr2_ph	1.0d-10	1.0d-14 to 1.0d-16 [2] [6]	Phonon SCF convergence threshold
conv_thr (pw.x)	1.0d-6	1.0d-8 to 1.0d-10 [2]	Electron SCF convergence threshold
forc_conv_thr	1.0d-3	1.0d-4 to 1.0d-5 [7] [6]	Force convergence in relaxation
etot_conv_thr	1.0d-4	1.0d-5 to 1.0d-6 [7]	Energy convergence in relaxation
k-point grid	Sparse (e.g., 4×4×1)	Denser grid (e.g., 12×12×1) [2] [6]	Brillouin zone sampling

A real-world example demonstrates how adjusting these parameters resolved issues: a user calculating phonons for bilayer MoS₂ reported that reducing tr2_ph and using a 4×4×1 k-point grid eliminated negative frequencies that appeared with looser settings [6].

Pseudopotential and Basis Set Issues

Inadequate plane-wave cutoffs or poor-quality pseudopotentials can artificially destabilize phonons:

• Wavefunction cutoff (ecutwfc): Insufficient values cause inaccurate force constants
• Charge density cutoff (ecutrho): Typically 4×ecutwfc for ultrasoft pseudopotentials, 8-12× for norm-conserving [8]
• Pseudopotential quality: Inconsistent behavior between different pseudopotential types (e.g., LDA vs. GGA, standard vs. stringent) may indicate pseudopotential issues [5] [6]

Acoustic Sum Rule Violations

The ASR requires that the frequency of acoustic modes at the Γ-point should be zero. Finite grid effects break translational invariance, typically resulting in small violations [5]:

• Diagnosis: Small imaginary frequencies (< 10 cm⁻¹) only for acoustic modes at Γ-point
• Solution: Apply ASR corrections using dynmat.x or matdyn.x with asr='crystal' [2]
• Verification: After applying ASR, acoustic modes should approach zero frequency (< 1 cm⁻¹) while other modes remain unchanged [5]

Genuine Structural Instability: Detection and Validation

Incomplete Structural Relaxation

A poorly relaxed structure is a common cause of genuine imaginary frequencies:

Table 2: Structural Relaxation Quality Checklist

Checkpoint	Indicator of Quality	Problematic Signature
Final forces	< 0.001 Ry/bohr [7]	Significant residual forces (> 0.01 Ry/bohr)
Stress tensor	Diagonal components < 0.5 GPa	High residual stress
Energy convergence	< 1.0d-5 Ry between steps [7]	Oscillating energy values
Atomic positions	Symmetry-preserved	Atoms drifted from symmetric positions

For 2D materials, ensure sufficient vacuum separation (typically 15-20 Å) to prevent spurious interactions between periodic images [6]. The assume_isolated parameter should be properly configured for low-dimensional systems [2].

Validating Genuine Instabilities

When computational artifacts have been ruled out, consider these valid causes of structural instability:

• True mechanical instability: The calculated structure may be in a saddle point, not a minimum
• Phase transitions: Imaginary modes may indicate a transition to a lower-energy structure
• Incorrect initial structure: Atomic positions or lattice parameters may be unphysical
• Exchange-correlation functional limitations: LDA/GGA approximations can sometimes incorrectly stabilize/destabilize certain phases [9]

Research on Y₂C MXenes demonstrates proper validation: after thorough convergence and multiple relaxation steps, the persistence of certain vibrational modes confirmed genuine dynamical stability [9].

The Researcher's Toolkit: Essential Workflow Solutions

Table 3: Research Reagent Solutions for Stable Phonon Calculations

Solution Component	Specific Implementation	Function
Structural Relaxation	vc-relax with cell_dofree='ibrav' [7]	Simultaneous optimization of atomic positions and lattice vectors
Convergence Testing	Systematic cutoff & k-point grids [2]	Determines parameters for numerically converged forces
Phonon Calculation	ph.x with ldisp=.true. and nq1,nq2,nq3 [10] [2]	Direct DFPT calculation on q-point grid
Post-processing	q2r.x + matdyn.x with asr='crystal' [2]	Fourier interpolation of force constants with ASR imposition
Validation	Phonon dispersion across entire Brillouin zone	Confirms absence of imaginary modes in stable materials

Step-by-Step Remediation Protocol

• Initial Diagnostic

  • Run matdyn.x with asr='crystal' to address acoustic mode issues [5] [2]
  • Check if imaginary frequencies persist only at specific q-points or throughout the Brillouin zone

• Systematic Convergence

  • Tighten tr2_ph to 1.0d-14 or lower [2] [6]
  • Increase k-point density, particularly for metallic systems [5]
  • Verify plane-wave cutoffs with convergence tests

• Structural Re-examination

  • Re-relax structure with tighter force convergence (forc_conv_thr=1d-4) [7]
  • For 2D systems, ensure adequate vacuum (~15 Å) and proper cell_dofree settings [6]
  • Verify atomic masses in ph.x input match those used in relaxation [5]

• Final Validation

  • Recaculate phonons with converged parameters
  • Apply ASR corrections consistently
  • For persistent imaginary modes, consider whether they indicate genuine instability

Key Recommendations for Reliable Phonon Calculations

• Always start with tighter convergence parameters than default, particularly for metallic systems or those with vacuum [5] [6]

• Validate your relaxed structure by checking residual forces and stresses before phonon calculations

• Apply ASR corrections systematically as part of your post-processing workflow [2]

• Use consistent parameters throughout your workflow - identical pseudopotentials, k-point grids, and convergence criteria between relaxation and phonon calculations

• When imaginary frequencies persist after thorough convergence testing, seriously consider they may reflect genuine physical instability requiring structural re-evaluation

By following this systematic approach, researchers can confidently distinguish between computational artifacts and genuine physical phenomena, ensuring the reliability of their phonon calculations and the scientific conclusions drawn from them.

The Role of the Acoustic Sum Rule (ASR) and Its Violations in Finite Calculations

Frequently Asked Questions (FAQs)

What is the Acoustic Sum Rule (ASR) and why is it important in phonon calculations?

The Acoustic Sum Rule (ASR) embodies the principle of translational invariance in physical systems. In the context of phonon calculations, it requires that the frequency of the three acoustic modes at the gamma point (Γ, q=0) must be precisely zero. This represents the fact that translating the entire crystal as a rigid body should not cost any energy. The ASR is a fundamental check on the quality and correctness of a phonon calculation.

I am not getting zero acoustic mode frequencies. Why does this violation occur?

Violations of the ASR in practical calculations occur because translational invariance is broken in approximated computations. The primary, most irreducible cause is the discreteness of the Fast Fourier Transform (FFT) grid used in Plane-Wave (PW) calculations. Other contributing factors can include insufficient convergence parameters (such as tr2_ph for phonons and conv_thr for the ground state) or, for molecules, a less-than-perfect convergence of the molecular geometry. If the nonzero frequencies are small (typically less than 10 cm⁻¹), this is usually not a cause for major concern, and the ASR can be imposed in post-processing with excellent results [11] [5].

Why do I get negative phonon frequencies, and what should I do?

"Negative" frequencies are technically imaginary frequencies (ω² ≤ 0). Their appearance can signal different issues:

• If they occur for acoustic frequencies at the Gamma point or for rotational modes of a molecule, they are likely a fictitious effect of finite supercell size or poor convergence, similar to ASR violations [11].
• In all other cases, it can be either a convergence problem (bad k-point grid, insufficient energy cutoffs ecutwfc/ecutrho, or overly large convergence thresholds) or it may signal a real mechanical instability of the chosen crystal structure with the approximations you used [5] [12]. A systematic convergence check and verification that your initial structure is reasonable are crucial first steps.

What are mysterious symmetry-related errors, and how can I fix them?

Errors like symmetry operation is non orthogonal, wrong representation, or wrong degeneracy are almost always a consequence of atomic positions that are close to, but not exactly at, high-symmetry positions. The code's symmetry detection has a finite tolerance, and small deviations can break it.

Solution: To prevent this, always define the Bravais lattice using the correct ibrav value (avoid ibrav=0), and use Wyckoff positions if known. This must be done in the initial self-consistent field (SCF) calculation, not just the phonon calculation [11] [5].

Troubleshooting Guides

Guide 1: Diagnosing and Fixing Acoustic Sum Rule (ASR) Violations

A small violation of the ASR is normal, but large violations indicate a problem. Follow this diagnostic workflow to identify and correct the issue.

Quantitative Data on Convergence Parameters

The following table summarizes key parameters that influence ASR violations and their recommended adjustments for mitigation.

Table 1: Key Parameters Affecting ASR Violations and Convergence

Parameter	Description	Role in ASR/Convergence	Recommended Action for Issues
conv_thr	SCF convergence threshold for the ground-state calculation.	A threshold that is too large introduces noise in forces, hindering a proper geometry optimization and phonon calculation [12].	Tighten to 1.0e-10 or lower for high-quality phonons [12].
tr2_ph	Convergence threshold for the phonon self-consistency.	Directly controls the accuracy of the phonon calculation. A larger value can cause bad frequencies and ASR violations [5].	Reduce to 1.0d-14 or lower [13].
ecutrho	Charge density cutoff (for USPP, PAW).	A discrete FFT grid is the primary cause of ASR violation. Increasing ecutrho densifies this grid, reducing the problem [5].	Increase significantly, especially for GGA calculations or systems with vacuum [5].
ASR	Acoustic Sum Rule imposition in matdyn.x.	Corrects for the inherent violation by enforcing the sum rule on the force constants [14].	Apply asr='crystal' or asr='all' for infrared-active solids [14] [13].

Guide 2: Resolving Negative (Imaginary) Frequencies

Imaginary frequencies can be the most challenging problem to solve, as they can be either a numerical artifact or a physical phenomenon.

Table 2: Troubleshooting Negative Phonon Frequencies

Symptom	Possible Cause	Diagnostic Steps	Solution
Imaginary frequencies only at Γ for acoustic/rotational modes	Finite-size effects or imperfect geometry convergence.	Check if the structure is fully optimized.	Improve the ground-state geometry optimization [11] [12].
Isolated imaginary frequencies at specific q-points	Numerical: Insufficient convergence.	Systematically converge ecutwfc, ecutrho, k-point grid, and conv_thr/tr2_ph [5].	Follow the convergence protocol in the next section.
Multiple or large-magnitude imaginary frequencies	Physical: Mechanical instability of the crystal structure.	Verify the initial structure is correct and stable.	The calculation might be correct, revealing an instability. Consider studying a different phase.
Imaginary frequencies in 2D materials	Incorrect treatment of electrostatic interactions.	Check the input parameters for 2D systems.	Use assume_isolated = '2D' in the SYSTEM namelist of the initial pw.x calculation [13].

Experimental Protocols

Systematic Convergence Protocol for Stable Phonons

To minimize numerical artifacts like ASR violations and imaginary frequencies, a rigorous convergence study is essential. Here is a detailed methodology:

• Ground-State Energy Cutoffs:

  • Perform a scf calculation on a single k-point (e.g., Γ) and determine the kinetic energy convergence for the wavefunctions (ecutwfc). The charge density cutoff (ecutrho) is typically 8-12 times ecutwfc (or 4 times for norm-conserving pseudopotentials). Using a higher cutoff than for standard scf calculations is often necessary for accurate phonons [13].

• K-Point Grid Convergence:

  • With the converged cutoffs, perform scf calculations with increasingly dense k-point grids until the total energy changes by less than 1 meV/atom. For metallic systems, pay special attention to the smearing parameter, as convergence can be very slow [5].

• Geometry Optimization:

  • Optimize the atomic positions and lattice parameters using the converged cutoffs and k-point grid. Use a tight force threshold (e.g., 1.0e-5 Ry/au or lower) to ensure the structure is at a true minimum. Crucially, converge the scf conv_thr to a very tight value (e.g., 1.0e-10) during this step to reduce noise in the forces [12].

• Phonon-Specific Parameters:

  • Using the fully optimized structure, perform a phonon calculation with a tight self-consistency threshold (tr2_ph=1d-14 or similar) [13].
  • For a full dispersion, use ldisp=.true. with a nq1, nq2, nq3 grid. The density of this q-grid must also be checked for convergence [15] [13].

• Post-Processing and Validation:

  • Use q2r.x to obtain real-space force constants.
  • Run matdyn.x to get the phonon dispersion. At this stage, apply the ASR (e.g., asr='crystal') to correct the small violations of the acoustic modes [14] [13].
  • Validate your results by ensuring acoustic modes are near zero at Γ and that no unphysical imaginary frequencies remain.

The Scientist's Toolkit: Essential Input Parameters for PHonon Calculations

Table 3: Key Input Parameters for ph.x and matdyn.x

Category	Parameter	Program	Function and Usage
File & Directory	prefix, outdir	ph.x	Must match the initial pw.x scf calculation to read the correct data [10].
	fildyn	ph.x	Specifies the output file for dynamical matrices [10].
	flfrc	matdyn.x	Input file containing the force constants from q2r.x [14].
Convergence	tr2_ph	ph.x	Self-consistency threshold for phonons; crucial for accuracy [10].
	niter_ph	ph.x	Maximum number of scf iterations in the phonon code [10].
Q-point Control	ldisp	ph.x	If .true., calculates phonons on a uniform grid defined by nq1, nq2, nq3 [10] [13].
	nq1, nq2, nq3	ph.x	Defines the uniform q-grid for ldisp=.true. [13].
Sum Rule & Corrections	asr	matdyn.x	Critical for fixing ASR. Options: 'no', 'simple', 'crystal', 'all'. Use 'crystal' for 3D, 'one-dim' for 2D, 'all' for infrared-active solids [14].
	loto_2d	matdyn.x	Set to .true. for correct 2D treatment of LO-TO splitting [14].

Frequently Asked Questions

• My ph.x calculation stops with an error about reading a file or a recover file. What should I do? This typically indicates a problem with the input data file from the pw.x SCF calculation. The file could be bad, incomplete, or produced by an incompatible version of the code. If the error mentions a recover file (cannot recover or error reading recover file), you likely have a corrupted restart file from a previous failed run. To fix this, remove all files named recover* in your outdir [5].

• Why does ph.x warn that "occupation numbers probably wrong"? This warning appears for metallic or spin-polarized systems when the smearing method for occupations is not used. You should set occupations = 'smearing' in your pw.x input file for the initial SCF calculation [5].

• My phonon dispersion from matdyn.x shows negative frequencies, but the original ph.x results were fine. Why? This common issue can have several causes. It often indicates that your initial structure is not fully relaxed, so ensure your geometry optimization uses tight convergence thresholds (e.g., forc_conv_thr = 1.0d-5 and etot_conv_thr = 1.0d-6). Other causes include an insufficient k-point grid for metals, or a charge density cutoff (ecutrho) that is too low, especially when using ultrasoft pseudopotentials (USPP) or the PAW method [5] [4]. For 2D materials, a specific fix is to set assume_isolated = '2D' in the SYSTEM namelist of your pw.x input to correctly handle long-range interactions and avoid imaginary acoustic frequencies near the gamma point [13].

• The acoustic modes at Γ do not have zero frequency. Is this an error? Not necessarily. The Acoustic Sum Rule (ASR) is never exactly fulfilled due to the system's lack of perfect translational invariance on a discrete grid. The frequency of acoustic modes is typically below 10 cm⁻¹ but can be higher. You can test this by diagonalizing the dynamical matrix with dynmat.x while imposing the ASR. If the acoustic mode frequency becomes very small (<1 cm⁻¹) and other modes remain unchanged, your results are trustworthy [5].

• I get "wrong degeneracy" or other symmetry-related errors. How can I fix them? These errors are often caused by atomic positions that are very close to, but not exactly at, high-symmetry positions. To prevent this, define the Bravais lattice using the correct ibrav value (avoid ibrav=0) and use Wyckoff positions when known in your self-consistent calculation with pw.x [5].

Troubleshooting Guide: Negative Frequencies

The appearance of negative (imaginary) frequencies in phonon calculations often signals a mechanical instability in the chosen structure. The table below summarizes the primary causes and solutions, tailored for different system types.

System Type	Common Causes	Recommended Fixes
All Systems	- Insufficient structure relaxation- SCF convergence threshold (conv_thr) too large- Phonon convergence threshold (tr2_ph) too large- Wrong atomic masses in input	- Use tighter force/energy convergence in pw.x relaxation (forc_conv_thr = 1.0d-5, etot_conv_thr = 1.0d-6) [4]- Reduce conv_thr and tr2_ph (e.g., 1.0d-14) [5] [4]- Verify amass for each atomic type [5]
2D Materials	- Incorrect treatment of long-range interactions due to periodic images	- In pw.x input, set assume_isolated = '2D' [13]- Apply ASR appropriate for 2D systems (e.g., asr = 'crystal' in q2r.x/matdyn.x) [13]
Metals	- Insufficient k-point grid- Slow convergence with smearing, especially with semicore states- Incommensurate phonon wave-vectors	- Use a denser k-point grid [5]- Employ a smaller smearing parameter (degauss) and ensure occupations='smearing' [5]
Systems with USPP/PAW	- Charge density cutoff (ecutrho) is too low	- Increase ecutrho (typically ecutrho = 8 * ecutwfc) [5]
General Numerical Issues	- Wavefunction cutoff (ecutwfc) is too low- Violation of the Acoustic Sum Rule (ASR)	- Systematically increase ecutwfc [5]- Impose ASR in q2r.x (zasr) and matdyn.x (asr); 'crystal' is often a good choice [4] [13]

Experimental Protocol: Phonon Dispersion Calculation

This workflow outlines the standard procedure for calculating phonon dispersion, incorporating system-specific checks.

Workflow for Phonon Dispersion

• SCF Ground-State Calculation (pw.x)

  • Purpose: Calculate the self-consistent ground-state charge density.
  • Input Key Points:
    • Use a high plane-wave kinetic energy cutoff (ecutwfc) and, for USPP/PAW, a high ecutrho [5].
    • For metals, set occupations='smearing' and use an appropriately dense k-point grid [5].
    • For 2D materials and isolated molecules, set assume_isolated='2D' or assume_isolated='1D' in the SYSTEM namelist to correctly handle electrostatic interactions [13].
    • Use a high convergence threshold (e.g., conv_thr = 1.0d-8 or tighter) [13].

• Phonon Calculations on a Mesh (ph.x)

  • Purpose: Compute the dynamical matrix on a uniform grid of q-points.
  • Input Key Points:
    • Set ldisp = .true. and define the mesh with nq1, nq2, nq3.
    • Use a tight phonon convergence threshold (e.g., tr2_ph = 1.0d-14) [4].
    • The fildyn parameter specifies the output file for the dynamical matrices.

• Generate Interatomic Force Constants (q2r.x)

  • Purpose: Perform an inverse Fourier transform to obtain the interatomic force constants in real space.
  • Input Key Points:
    • Link to the dynamical matrix file via fildyn.
    • Apply an Acoustic Sum Rule (ASR) to fix inaccuracies for acoustic modes at Γ (e.g., zasr = 'crystal') [4] [13].

• Calculate Phonon Dispersion and DOS (matdyn.x)

  • Purpose: Fourier transform the force constants to get the dynamical matrix and frequencies at any q-point, generating the band structure and density of states.
  • Input Key Points:
    • Set flfrc to the force constants file from q2r.x.
    • Use asr = 'crystal' to consistently apply the ASR [4].
    • For a band structure, set q_in_band_form = .true. and provide the q-point path.
    • For DOS, set dos = .true. and define a mesh with nk1, nk2, nk3.

Research Reagent Solutions

This table lists the key software components and their roles in a typical Quantum ESPRESSO phonon calculation workflow.

Item	Function in the Experiment
pw.x	The core plane-wave self-consistent field (SCF) code for calculating the electronic ground state and performing structural relaxations [13].
ph.x	The density functional perturbation theory (DFPT) code for calculating phonon frequencies and eigenvectors on a uniform q-point mesh [13].
q2r.x	Computes the interatomic force constants in real space by Fourier transforming the dynamical matrices from ph.x [13].
matdyn.x	Calculates the phonon dispersion along any q-point path and the phonon density of states using the force constants from q2r.x [13].
dynmat.x	A utility to diagonalize the dynamical matrix, often used to test the effect of different Acoustic Sum Rules (ASR) on the phonon frequencies [5].
Pseudopotential Files	Define the effective interaction between ions and valence electrons. The choice and quality of pseudopotentials (e.g., NC, USPP, PAW) are critical for accuracy.

Robust Phonon Calculation Methodology: Best Practices from SCF to Dispersion

Why is a Proper Ground-State SCF Calculation Crucial for Phonons?

A self-consistent field (SCF) calculation aims to find the electronic ground state of a system by iteratively solving the Kohn-Sham equations until the energy and electron density no longer change significantly. This ground-state charge density is the fundamental starting point for subsequent phonon calculations using ph.x. Imaginary frequencies (often reported as "negative frequencies") in phonon calculations can often be traced back to an inaccurate or un-converged ground state. These frequencies signify an instability, which can either be a real physical property of the system or, more often, an artifact of a poor-quality SCF calculation or insufficient structural relaxation [5] [11] [2].

The flowchart below outlines the logical progression for diagnosing and fixing the root causes of negative phonon frequencies.

SCF Parameter Checklist for Stable Phonons

The following table summarizes the key parameters in your SCF input that require careful attention to ensure a high-quality ground state for phonon calculations.

Parameter (Namelist)	Recommended Setting for Phonons	Function & Rationale
calculation (&CONTROL)	'scf'	Specifies a self-consistent field calculation to obtain the ground-state density [16].
conv_thr (&ELECTRONS)	1.0d-12 to 1.0d-16 [2] [12]	The convergence threshold for the SCF cycle. A tighter threshold (smaller number) is critical to reduce numerical noise in forces [12].
ecutwfc (&SYSTEM)	System-dependent (converged)	Plane-wave kinetic energy cutoff for wavefunctions. Must be high enough for total energy convergence [5] [16].
ecutrho (&SYSTEM)	4 * ecutwfc (for USPP/PAW) [5]	Cutoff for charge density. Too low a value can violate translational invariance, affecting acoustic modes [5].
occupations (&SYSTEM)	'smearing' for metals/spin-polarized; 'fixed' for insulators [5] [17]	Incorrect smearing for metallic systems leads to errors. Macroscopic field contributions require 'fixed' occupations [5] [11].
K_POINTS	Dense grid for metals	A converged k-point grid is essential. Metallic systems with semicore states require particularly careful convergence [5].
forc_conv_thr (&IONS)	1.0d-8 or tighter (in VC-RELAX)	Force convergence threshold in geometry optimization. Poor structural relaxation is a primary cause of imaginary frequencies [2] [12].

Advanced Troubleshooting and Parameter Adjustment

Diagnosing Common Error Messages

• "ph.x does not yield acoustic modes with zero frequency at Γ": This is often due to the Acoustic Sum Rule (ASR) violation because the finite FFT grid breaks perfect translational invariance [5] [11]. The frequency of acoustic modes at Gamma is typically small (<10 cm⁻¹). You can impose the ASR in dynmat.x or matdyn.x using the asr flag (e.g., asr='crystal') to fix this [5] [2].
• "ph.x yields really lousy phonons, with bad or negative frequencies": This error can have multiple origins [5]. Check the following:
  • Atomic masses: Verify that the masses in ATOMIC_SPECIES and ph.x input are correct.
  • SCF convergence: Tighten conv_thr in the SCF calculation and tr2_ph in the ph.x input [5] [2].
  • Structural optimization: Ensure the structure is fully relaxed to its minimum energy configuration with sufficiently tight force convergence [12] [18].
• "Mysterious symmetry-related errors": These errors (e.g., wrong degeneracy) usually occur when atomic positions are close to, but not exactly on, high-symmetry positions. To prevent this, use the correct ibrav value and Wyckoff positions in the SCF calculation instead of ibrav=0 [5] [11].

The Scientist's Toolkit: Essential SCF and Phonon Reagents

Item / Utility	Function in Quantum ESPRESSO
pw.x	The main PWscf code for performing SCF, structural relaxation, and MD calculations [16].
ph.x	The PHonon code for calculating linear response phonons and vibrational spectra [5] [15].
dynmat.x / matdyn.x	Post-processing codes to diagonalize the dynamical matrix and enforce the Acoustic Sum Rule (ASR) [5] [2].
Pseudopotential File (UPF)	File containing the atomic pseudopotential, defining the electron-ion interaction [16].
conv_thr & tr2_ph	Numerical parameters controlling the convergence of the SCF (pw.x) and phonon (ph.x) calculations, respectively [5] [2].
vc-relax	A calculation type that performs a variable-cell relaxation to optimize both atomic positions and unit cell vectors [2].

A Protocol for Systematic Convergence Testing

To methodically eliminate numerical artifacts, follow this workflow for parameter convergence testing.

• Converge ecutwfc and ecutrho: Start with a single-point SCF calculation for your material. Systematically increase ecutwfc until the change in total energy is smaller than your desired accuracy (e.g., 1 meV/atom). Set ecutrho = 4 * ecutwfc if using ultrasoft pseudopotentials or the PAW method [5] [16].
• Converge the k-point grid: With your converged cutoffs, perform SCF calculations with increasingly dense k-point grids (K_POINTS automatic) until the total energy converges [19].
• Perform rigorous geometry relaxation: Use the converged parameters from steps 1 and 2 in a calculation = 'vc-relax' run. Set forc_conv_thr to a very tight value, such as 1.0d-8 Ry/Bohr or lower, to ensure forces are minimized [2] [12].
• Final high-accuracy SCF: Using the relaxed geometry from step 3, run a final scf calculation with an exceptionally tight conv_thr (e.g., 1.0d-12 to 1.0d-16) [2] [12]. The output density of this calculation is your high-quality ground state for ph.x.
• Run ph.x with tight convergence: In your ph.x input, set a tight convergence threshold with tr2_ph (e.g., 1.0d-14 to 1.0d-16) to match the quality of your SCF calculation [2].

By meticulously following this protocol and understanding the role of each critical parameter, you can significantly increase the reliability of your phonon calculations and confidently interpret the results, distinguishing numerical artifacts from true physical instabilities.

Density Functional Perturbation Theory (DFPT) calculations using Quantum ESPRESSO's ph.x code represent a powerful approach for determining lattice dynamical properties, yet researchers frequently encounter the troubling issue of negative (imaginary) frequencies in their results. These non-physical frequencies indicate dynamical instabilities in the calculated structure and often stem from inaccuracies in either the initial structural model or the computational parameters used in the phonon calculation. Within the context of thesis research on fixing negative frequencies, this technical guide synthesizes practical troubleshooting knowledge and parameter optimization strategies to help computational materials scientists, chemists, and drug development researchers achieve stable, physically meaningful phonon results. The occurrence of negative frequencies typically signals that the system is mechanically unstable with the chosen approximations, requiring methodical investigation across the entire computational workflow from structural relaxation to phonon calculation convergence.

Troubleshooting Guide: Diagnosing Negative Frequencies

Systematic Problem Identification and Resolution

Table: Common Causes and Solutions for Negative Frequencies

Problem Category	Specific Symptoms	Recommended Solutions	Key ph.x/pw.x Parameters to Adjust
Structural Issues	Negative frequencies persist across multiple q-points; structure not at energy minimum	Tighten relaxation convergence; verify forces and stresses	forc_conv_thr = 1e-5 (tighter) [7], etot_conv_thr = 1e-8 [7], cell_dofree = 'ibrav' [7]
Pseudopotential & Basis Set	Negative frequencies especially in covalent systems; convergence with ecut unclear	Increase plane-wave cutoffs; verify pseudopotential quality	ecutwfc (increase) [5], ecutrho = 8*ecutwfc (USPP/PAW) [5]
k-point Sampling	Metallic systems with semicore states; slow convergence with k-points	Use denser k-point grids; adjust smearing	ngkpt (increase) [5], smearing = 'gaussian' [5], degauss (reduce) [5]
Phonon Calculation Convergence	Inconsistent results across q-points; SCF convergence warnings	Tighten phonon convergence thresholds; adjust mixing parameters	tr2_ph = 1e-14 (tighter) [4], niter_ph = 200 (increase) [10], alpha_mix(niter) [10], nmix_ph = 8 [10]
Acoustic Sum Rule (ASR)	Significant acoustic mode violations at Γ-point (>10 cm⁻¹)	Apply ASR corrections in post-processing	asr = 'crystal' (in q2r.x/matdyn.x) [4]

Advanced Diagnostic Workflow

Figure 1: Diagnostic workflow for systematic identification and resolution of negative frequency issues in ph.x calculations.

The diagnostic pathway begins with verifying the quality of the initial structure relaxation, as an improperly relaxed structure represents the most common source of imaginary frequencies [4]. The workflow proceeds through checking residual forces, k-point convergence for the specific system type (with particular attention to metallic systems with semicore states) [5], phonon SCF convergence, and finally application of acoustic sum rule corrections where appropriate. At each stage, failure to meet convergence criteria necessitates returning to earlier computational stages with adjusted parameters.

Essential Input Parameters for Stable Phonon Calculations

Critical ph.x Input Parameters

Table: Key ph.x Input Parameters for Stable DFPT Results

Parameter	Recommended Value	Purpose	Effect on Stability
tr2_ph	1e-12 to 1e-14	Phonon SCF convergence threshold	Tighter values prevent premature convergence; essential for metallic systems [4]
niter_ph	100-200	Maximum phonon SCF iterations	Prevents early termination; allows difficult perturbations to converge [10]
alpha_mix(niter)	0.3-0.7 (iteration-dependent)	Potential mixing factor	Improves SCF convergence; lower values stabilize oscillating solutions [10]
nmix_ph	4-8	Number of iterations for potential mixing	Higher values can speed convergence but use more memory [10]
verbosity	'high'	Detailed output	Enables monitoring of convergence per perturbation [4]
recover	.false. (initially)	Restart from interrupted run	Set to .true. only to continue failed calculations [10]
ldisp	.true.	Automatic q-point grid generation	For full phonon dispersion with nq1, nq2, nq3 [4]
asr	(in q2r.x/matdyn.x)	Acoustic sum rule imposition	Corrects acoustic mode violations; 'crystal' type recommended [4]

Pre-Phonon Calculation Checklist

Before executing ph.x, researchers should systematically verify these precursor calculations:

• Ground-State Convergence: Confirm the SCF calculation achieves at least 1e-10 Ry energy convergence, with particularly tight thresholds (1e-8 Ry or better) for metallic systems or those with delicate electronic structure [7].

• Structural Optimization Quality: Ensure geometry optimization (relax or vc-relax) reaches force thresholds below 0.001 Ry/Bohr, with stress tensor components below 0.5 GPa for solid systems [7].

• k-point Grid Validation: Verify the k-point sampling produces well-converged Fermi surface properties for metals or band gaps for insulators, with particularly dense grids required for systems with nested Fermi surfaces or van Hove singularities [5].

• Pseudopotential Verification: Confirm pseudopotentials accurately reproduce all-electron equations of state and elastic properties, with particular attention to semicore states and magnetic systems [5].

Frequently Asked Questions (FAQs)

Common Implementation Challenges

Q1: Why do I obtain negative frequencies even after thorough structure relaxation?

This persistent issue often stems from multiple contributing factors. First, verify that your relaxation utilized sufficiently tight convergence criteria (forc_conv_thr=1e-5 Ry/Bohr or tighter, etot_conv_thr=1e-8 Ry) [7]. Second, certain exchange-correlation functionals (particularly GGAs) are known to produce dynamical instabilities for specific materials where experimental evidence shows stability [5]. In such cases, hybrid functionals or van der Waals corrections may be necessary. Third, ensure your k-point grid is commensurate with the phonon q-point grid, as mismatches can introduce spurious instabilities, particularly for metallic systems [5].

Q2: How can I distinguish physical instabilities from numerical artifacts?

Genuine physical instabilities typically manifest as consistent negative frequencies across multiple q-points, particularly at zone boundary points where soft modes often occur. Numerical artifacts, conversely, tend to be isolated to specific q-points, especially those non-commensurate with the k-point grid [5]. To test this, perform two calculations: first with your standard parameters, then with increased ecutwfc (by 20-30%) and denser k-point sampling (2× in each direction). If the negative frequencies persist with similar magnitude across both calculations, they likely represent genuine physical instabilities [5].

Q3: What does "Wrong degeneracy error in star_q" indicate and how can it be resolved?

This error arises when the specified q-vector differs slightly from a high-symmetry point due to numerical precision issues [5]. The code compares the input q-vector with symmetry-equivalent points using a tolerance of 10⁻⁵, and discrepancies within this range can cause symmetry recognition failures. The solution is to ensure high-symmetry q-points are specified with exact fractional coordinates (e.g., 0.5 instead of 0.500001) or to use the automatic q-point generation via ldisp=.true. with nq1, nq2, nq3 [5] [4].

Q4: Why do acoustic modes at Γ-point show significant non-zero frequencies?

This Acoustic Sum Rule (ASR) violation occurs because the exchange-correlation energy computation on discrete real-space grids breaks exact translational invariance [5]. The problem is typically more severe for GGA functionals than LDA, and for systems with vacuum regions (surfaces, molecules) where reduced gradient diverges [5]. While frequencies below 10 cm⁻¹ are generally acceptable, larger violations can be corrected by applying the ASR in q2r.x and matdyn.x using asr='crystal' [4]. The ultimate test is to diagonalize the dynamical matrix with imposed ASR; if acoustic modes then show frequencies <1 cm⁻¹ with other modes unchanged, the corrected results are reliable [5].

Q5: How can I optimize computational efficiency while maintaining accuracy?

For multi-q-point calculations, avoid computing all perturbations in a single execution. Instead, run separate calculations for different q-points or irreducible perturbations, as each has different computational workload and symmetry constraints [20]. Use more MPI processes for computationally heavy perturbations (those at low-symmetry q-points) and fewer for lighter ones (high-symmetry points). This approach significantly improves overall computational efficiency compared to uniform parallelization [20].

Research Reagent Solutions: Essential Computational Tools

Table: Critical Computational Components for Stable Phonon Calculations

Component	Function	Implementation Notes
Optimized Structures	Provides energy-minimized starting configuration	Always use vc-relax with tight convergence before phonons [7]
Hard Pseudopotentials	Accurate electron-ion interaction representation	Prefer PAW or NC pseudopotentials with high ecutwfc when negative frequencies persist [5]
k-point Grids	Brillouin zone sampling for electronic structure	Metallic systems require denser grids; ensure commensurability with q-points [5]
Mixing Parameters	Stabilize SCF convergence in DFPT	Use alpha_mix and nmix_ph to handle difficult perturbations [10]
ASR Corrections	Imposes translational invariance	Apply in post-processing with q2r.x and matdyn.x [4]
Symmetry Tools	Identifies irreducible representations	Use ibrav (not 0) and Wyckoff positions to minimize symmetry errors [5]

Frequently Asked Questions (FAQs)

Q1: My ph.x calculation shows positive frequencies, but matdyn.x produces negative ones at the same q-points. Why? This common discrepancy arises because the dynamical matrices from ph.x are calculated on a coarse q-grid. The interpolation performed by q2r.x and matdyn.x can introduce errors if the initial grid is too sparse, force constants are imperfect, or the Acoustic Sum Rule (ASR) is misapplied [4]. The ultimate test is to diagonalize the dynamical matrix with program dynmat.x, imposing the ASR [5].

Q2: What is the most critical factor for eliminating imaginary frequencies near the Γ-point? Ensuring your crystal structure is fully and accurately relaxed is paramount. A poorly relaxed structure, even with minimal residual forces, is a primary cause of unphysical imaginary frequencies, as the system is not in a true energy minimum [4] [2]. For 2D materials, correctly activating the 2D treatment with assume_isolated = '2D' in the pw.x input and loto_2d = .true. in the matdyn.x input is also crucial [13] [2].

Q3: How do I choose between different Acoustic Sum Rule (ASR) options? The ASR corrects for broken translational invariance. The 'crystal' type is often a robust choice as it imposes translational invariance via an optimized correction [14] [13]. For infrared-active solids, using asr = 'all' in matdyn.x along with read_lr = .true. is required to include long-range force constants, as recommended in recent research [14].

Q4: My dispersion curve looks "messed up" compared to my ph.x results. What went wrong? This typically indicates a problem during the interpolation workflow. Key suspects are an insufficiently dense q-grid in the initial ph.x calculation or an inconsistency in atomic masses between the ph.x and matdyn.x inputs [5] [4]. Always verify that the masses specified with the amass keyword are identical and correct in both inputs.

Troubleshooting Guide: Fixing Negative Frequencies

Precursor Checks: Ensuring a Solid Foundation

Before delving into complex fixes, verify these foundational elements.

• Validate Pseudopotentials and Input Files: Ensure you are using appropriate, well-tested pseudopotentials and that the outdir directory in all input files points to the correct location containing the SCF charge density [5].
• Confirm Accurate Atomic Masses: Wrong atomic masses given in input will yield wrong frequencies [5]. Double-check the amass values for each atom type in both ph.x and matdyn.x inputs.
• Check File Permissions: Ensure the code has permission to read and write all necessary files (e.g., flfrc, fildyn). Remove old recover* files from outdir if a previous calculation failed [5].

Systematic Workflow Optimization

Follow this optimized workflow, paying close attention to the key parameters in the tables below.

• Step 1: Self-Consistent Field (SCF) Calculation with pw.x This initial step must provide a highly converged electronic ground state. Use tighter thresholds than for a single-point energy calculation [13].

  Table: Key pw.x &ELECTRONS Namelist Parameters for Phonon Accuracy

  Parameter	Recommended Value	Function
  conv_thr	1.0d-12 to 1.0d-16	Convergence threshold for SCF cycle [4] [2].
  mixing_beta	0.2 to 0.7	Mixing factor for self-consistency [13] [2].
  ecutwfc / ecutrho	Higher than normal	Plane-wave cutoffs; increase if negative frequencies persist [5] [13].

• Step 2: Dynamical Matrix Calculation with ph.x This is the most computationally demanding step. A denser q-grid is critical for accurate interpolation [21].

  Table: Key ph.x &INPUTPH Namelist Parameters

  Parameter	Recommended Value	Function
  tr2_ph	1.0d-14 to 1.0d-16	Convergence threshold for phonon self-consistency [13] [4] [2].
  nq1, nq2, nq3	System-dependent (e.g., 6 6 6)	Dimensions of the uniform q-point grid for DFPT calculations [13] [21].
  ldisp	.true.	Enables calculation on a uniform grid of q-points [13] [2].

• Step 3: Fourier Transform to Real Space with q2r.x This step converts dynamical matrices into real-space force constants.

  Table: Key q2r.x &INPUT Namelist Parameters

  Parameter	Recommended Value	Function
  zasr	'crystal' or 'all'	Applies Acoustic Sum Rule during the force constant generation [13] [2].

• Step 4: Interpolated Dispersion with matdyn.x This final step generates the phonon dispersion along high-symmetry paths.

  Table: Key matdyn.x &INPUT Namelist Parameters for Frequency Accuracy

  Parameter	Recommended Value	Function
  asr	'crystal' or 'all'	Type of Acoustic Sum Rule to apply; must be consistent with q2r.x [14] [13].
  loto_2d	.true. for 2D systems	Activates correct treatment of LO-TO splitting in 2D [14] [13].
  loto_disable	.false. (default)	Set to .true. only for testing to disable LO-TO splitting at q=0 [14].
  read_lr	.true. for asr='all'	Reads long-range force constants, crucial for infrared-active solids [14].

Advanced Diagnostics and Solutions

If problems persist, consider these advanced strategies.

• Grid Convergence Study: Systematically increase the nq1, nq2, nq3 grid in the ph.x calculation. A converged phonon dispersion should not change significantly with a denser grid.
• Structural Re-inspection: For low-symmetry crystals, ensure atomic positions are exactly at their high-symmetry points. Even slight deviations can cause "mysterious symmetry-related errors" in ph.x [5].
• Functional Considerations: Be aware that certain exchange-correlation functionals (particularly some GGAs) may predict a dynamical instability for your material that others (like LDA) do not [4]. This could be a physical prediction rather than a numerical error.

The Scientist's Toolkit: Essential Computational Reagents

Table: Key "Research Reagent" Solutions for Quantum ESPRESSO Phonon Calculations

Item	Function in Workflow
pw.x	Computes the electronic ground state and charge density, forming the foundation for all subsequent perturbative calculations [13] [22].
ph.x	Implements Density Functional Perturbation Theory (DFPT) to calculate the dynamical matrices on a uniform grid of q-points [13] [22].
q2r.x	Performs an inverse Fourier transform of the dynamical matrices to obtain the interatomic force constants in real space [13] [22].
matdyn.x	Interpolates the real-space force constants to calculate phonon frequencies at any q-point, enabling the plotting of smooth dispersion curves and DOS [14] [13].
Acoustic Sum Rule (ASR)	A constraint correction applied to force constants to enforce translational invariance, crucial for having acoustic modes go to zero at the Γ-point [14] [5].

This guide provides targeted troubleshooting for Quantum ESPRESSO ph.x calculations, focusing on resolving negative frequencies in complex systems like 2D materials and metallic compounds. These issues often point to structural instabilities or computational settings mismatches rather than simple code errors. The following FAQs and protocols will help you diagnose and fix these problems systematically.

Frequently Asked Questions (FAQs)

FAQ 1: My ph.x calculation finishes successfully, but I get negative frequencies when I run matdyn.x to plot the phonon dispersion. What is wrong?

This is a common issue that typically originates from the initial self-consistent field (SCF) calculation performed with pw.x, not the phonon code itself. The most probable cause is that your crystal structure was not fully relaxed to its ground state. The residual forces and stresses on the atoms mean that the dynamical matrix is calculated for a structure that is not at a minimum of its potential energy surface, leading to imaginary phonon modes. To fix this, ensure your structure relaxation uses tight convergence criteria for both energy and forces before starting phonon calculations [4].

FAQ 2: What does the error "occupation numbers probably wrong" mean, and how do I fix it?

This warning appears when you are calculating a metallic or spin-polarized system but have not set the occupations variable to 'smearing' in your pw.x input. Metals require smearing (e.g., Marzari-Vanderbilt, Methfessel-Paxton) to properly handle partial occupancy of electronic states near the Fermi level. Without it, the calculation of total energy and forces can be incorrect. To resolve this, in your pw.x input, set occupations = 'smearing' and choose an appropriate smearing parameter (e.g., 'mv') and a small degauss value [5].

FAQ 3: The frequencies of the acoustic modes at the gamma point (q=0) are not zero. Is this an error?

Not necessarily. A small deviation (e.g., less than 10 cm⁻¹) is expected due to the Acoustic Sum Rule (ASR) violation. This occurs because the system is never perfectly translationally invariant in practice, as the exchange-correlation energy is computed on a discrete real-space grid. The problem can be more severe in 2D systems or those with vacuum, and with GGA functionals compared to LDA. You can impose the ASR during post-processing using the dynmat.x tool. If the acoustic mode frequency drops to nearly zero (e.g., < 1 cm⁻¹) while other modes remain unchanged, your results are trustworthy [5].

FAQ 4: My system is a 2D material. Are there any special considerations to avoid negative frequencies?

Yes. 2D materials, often simulated with a vacuum layer, can be particularly sensitive to several factors:

• Pseudopotentials: The use of ultrasoft pseudopotentials (USPP) or projector-augmented waves (PAW) requires paying close attention to both the wavefunction cutoff (ecutwfc) and the charge density cutoff (ecutrho), which is typically ecutrho = 4 * ecutwfc [5].
• k-point Grid: A sufficiently dense k-point grid is crucial for convergence, especially for metallic 2D systems [5].
• Vacuum Size: Ensure the vacuum layer is thick enough to prevent spurious interactions between periodic images.

Troubleshooting Guide: Fixing Negative Frequencies

Protocol 1: Verifying Structural Relaxation

A poorly relaxed structure is the most common source of unphysical negative frequencies.

• Objective: Ensure the crystal structure is at a local energy minimum.
• Procedure:
  • Return to your pw.x structural relaxation (calculation = 'vc-relax' or 'relax') input.
  • Apply tighter convergence thresholds. It is recommended to set:
    • etot_conv_thr = 1.0d-10 or lower (e.g., 1.0d-12)
    • forc_conv_thr = 1.0d-8 or lower (e.g., 1.0d-10) [4]
  • For variable-cell relaxations, set tstress = .true. and ensure the stress tensor components are also converged.
  • Re-run the full relaxation and SCF calculation before restarting the ph.x calculation.

Protocol 2: Optimizing SCF and Phonon Calculation Parameters

Insufficient convergence in the SCF or phonon steps can cause numerical instabilities.

• Objective: Achieve high numerical accuracy in the energy, potential, and phonon calculations.
• Procedure:
  • In your pw.x SCF input, tighten the convergence threshold: conv_thr = 1.0d-12 or lower [4].
  • In your ph.x input, tighten the phonon convergence threshold: tr2_ph = 1.0d-14 or lower [4].
  • Increase the planewave cutoff energies. Systematically test higher values for ecutwfc and, if using USPP/PAW, ecutrho.
  • For metallic systems, ensure you are using smearing and test denser k-point grids, as convergence can be slow, especially with semicore states [5].

Protocol 3: Addressing Specifics for Metallic and Topological Systems

Metals and systems with topological properties require special attention to sampling and electron-phonon coupling.

• Objective: Correctly handle metallic states and dense phonon spectra.
• Procedure:
  • Always use occupations = 'smearing' in the pw.x input.
  • Use a very dense k-point grid. Convergence should be checked by observing if phonon frequencies stabilize as the grid density increases.
  • The nqc parameter in a ph.x metals calculation should be commensurate with the k-point grid to avoid errors from non-commensurate q-vectors [5].
  • Research on 2D superconducting borides (e.g., Nb₂B₂H₄) highlights the importance of high-frequency phonons from light elements (B, H) in achieving high-Tc superconductivity. Accurate phonon calculations are essential for predicting such properties [23].

Workflow for Robust Phonon Calculations

The following diagram outlines the logical workflow for diagnosing and fixing negative frequencies, integrating the protocols above.

Parameter Tables for Challenging Systems

Table 1: Recommended Convergence Parameters

These values are starting points. Always perform convergence tests for your specific system.

Parameter	Typical Standard Value	Recommended Tight Value	Function
forc_conv_thr	1.0d-3 Ry/Bohr	1.0d-8 to 1.0d-10 Ry/Bohr	Force convergence in relaxation [4]
etot_conv_thr	1.0d-4 Ry	1.0d-10 to 1.0d-12 Ry	Energy convergence in relaxation [4]
conv_thr	1.0d-6 Ry	1.0d-11 to 1.0d-12 Ry	SCF energy convergence [4]
tr2_ph	1.0d-12	1.0d-14	Phonon calculation convergence [4]
ecutwfc	System-dependent	20-30% higher	Plane-wave kinetic energy cutoff [5]

Table 2: System-Specific Considerations

System Type	Key Parameters	Potential Pitfalls
2D Materials (e.g., MXC₃, M₂B₂H₄) [24] [23]	Large vacuum layer; Dense k-grid in 2D; ASR = 'crystal' in q2r.x/matdyn.x.	Spurious interaction from small vacuum; Incorrect treatment of long-range interactions.
Metallic Compounds	occupations = 'smearing'; Very dense k-point grid; Commensurate q-point grid.	Slow k-point convergence; "Occupation numbers probably wrong" error [5].
Systems with Vanadium (e.g., V₂B₂H₄) [23]	Check for magnetic ordering (FM, AFM).	Magnetic ground state might be missed in non-spin-polarized calculations.

The Scientist's Toolkit: Key Research Reagents & Materials

The study of novel 2D materials often reveals new platforms for topological superconductivity. The following table lists material classes from recent research that are relevant for computational studies of phonons and superconductivity.

Material / System	Function / Role in Research	Key Characteristic (from first-principles)
2D MXC₃(M:X = In:As, Se:As, In:Te) [24]	Platform for 2D topological superconductivity.	Predicted to be thermodynamically/dynamically stable; Phonon-mediated superconductors with Tc up to 33.6 K.
Hydrogenated Transition-Metal Borides (e.g., Nb₂B₂H₄, V₂B₂H₄) [23]	High-Tc 2D BCS superconductor at ambient pressure.	Hydrogenation enhances high-frequency phonons; Calculated Tc of 69 K and 83 K, respectively.
Rutile Structure (e.g., TiO₂) [4]	Prototype material for testing computational methods.	Known to exhibit phonon instabilities (negative frequencies) when calculated with GGA functionals.

Systematic Troubleshooting: Diagnosing and Eliminating Imaginary Frequency Sources

A systematic guide to eliminating unphysical negative phonon frequencies in your Quantum ESPRESSO calculations

Understanding Negative Phonon Frequencies

What does a negative phonon frequency indicate? A negative phonon frequency (often reported as an imaginary frequency) signals a mechanical instability in your simulated crystal structure [5]. Fundamentally, it means the system can lower its energy by distorting along the corresponding vibrational mode. While this can sometimes indicate a genuine physical instability, more often it results from insufficient convergence of key computational parameters in the ph.x or preceding pw.x calculation [5] [2].

The Acoustic Sum Rule (ASR) and Acoustic Modes For non-metallic systems at the gamma point (q=0), the frequencies of the three acoustic modes should be precisely zero, fulfilling the Acoustic Sum Rule. In practice, finite grid calculations lead to small violations [5]. Frequencies below 10 cm⁻¹ are typically acceptable, but values approaching 100 cm⁻¹ suggest problems. The ultimate test is to re-diagonalize the dynamical matrix using dynmat.x with the ASR imposed. If this yields acoustic modes with frequencies below 1 cm⁻¹, you can trust your other results [5].

Core Convergence Parameters: A Troubleshooting Guide

The table below summarizes the primary parameters to test when troubleshooting negative frequencies.

Parameter	Description & Impact	Typical Symptoms of Poor Convergence	Recommended Verification Steps
ecutwfc / ecutrho	Plane-wave kinetic energy cutoffs for wavefunctions/charge density [25] [8]. Too low values are a primary cause of bad frequencies [5].	"Lousy" phonons, negative frequencies, gross ASR violations [5].	Perform a dedicated convergence test for the total energy and forces in the SCF calculation before any phonon run [26].
k-point Grid	Sampling of the Brillouin zone for the SCF calculation [25]. Convergence is slow for metals, especially with semicore states [5].	Negative frequencies, particularly for metallic systems [5].	Test increasingly dense grids until total energy and phonon frequencies stabilize.
tr2_ph	Convergence threshold for the self-consistency of the phonon calculation [10]. A threshold that is too large can yield bad frequencies [5].	Phonon frequencies fail to converge, leading to inaccuracies.	Systematically lower tr2_ph (e.g., to 1.0d-14) and monitor the convergence of the dynamical matrix [2].
conv_thr	Convergence threshold for the self-consistency of the initial SCF calculation [25].	A structure that is not fully relaxed propagates errors to the phonon calculation.	Use a tight threshold (e.g., 1.0d-12 to 1.0d-16) in the pw.x relaxation and SCF calculations [2] [13].
Atomic Masses (amass)	Masses of atomic species used to build the dynamical matrix [10]. Wrong masses yield wrong frequencies [5].	All phonon frequencies are incorrectly scaled, though the force constants are valid [5].	Verify the amass values in the ph.x input match those in the pseudopotential files.

Systematic Convergence Verification Protocol

Follow this workflow to methodically identify and resolve the source of negative frequencies.

Step 1: Verify the Ground-State Structure

A poorly relaxed structure is a common source of instability [2].

• Perform a Variable-Cell Relaxation: Use calculation = 'vc-relax' in pw.x to optimize both atomic positions and lattice vectors.
• Converge Force and Energy Thresholds: In the &IONS namelist, set forc_conv_thr to a stringent value (e.g., 1.0d-5 Ry/Bohr or ~0.00026 eV/Å). In the &ELECTRONS namelist, ensure conv_thr is tight (e.g., 1.0d-12 or lower) [2].
• Check Final Forces: Confirm that the maximum force on any atom in the relaxed structure is significantly below your convergence threshold.

Step 2: Converge Core SCF Parameters inpw.x

Before running ph.x, you must establish a fully converged electronic ground state.

• ecutwfc and ecutrho Convergence:

  • Perform a series of SCF calculations with increasing ecutwfc.
  • Plot the total energy versus ecutwfc. The value is sufficiently converged when the energy change is smaller than the conv_thr used in your production run.
  • Set ecutrho = 4 * ecutwfc for norm-conserving pseudopotentials, or 8 * ecutwfc for ultrasoft pseudopotentials [8]. For higher accuracy in phonon calculations, consider using a higher cutoff than the minimum sufficient for energy convergence [13].

• k-point Grid Convergence:

  • With a converged cutoff, perform SCF calculations with increasingly dense k-point grids (e.g., 8x8x8, 12x12x12, 16x16x16).
  • Plot the total energy against the k-point grid density. Choose a grid where the energy change is negligible.
  • For Metals: Convergence is slower. Use a finer k-point grid and appropriate smearing (e.g., smearing = 'gaussian' and degauss value) [5] [8].

Step 3: Converge Phonon-Specific Parameters inph.x

Once the SCF is robust, focus on the phonon calculation itself.

• tr2_ph Threshold: This is the convergence threshold for the phonon self-consistent cycle [10]. If your pw.x conv_thr was 1.0d-12, try setting tr2_ph = 1.0d-14 for the phonon calculation [2].
• Mixing Parameters: Slow convergence in the phonon calculation can be improved by adjusting alpha_mix(niter) (the mixing factor) and increasing nmix_ph (the number of iterations used in potential mixing) to values between 8 and 20, which can speed up convergence at the cost of more memory [10].

Step 4: Apply Post-Processing Checks and Corrections

• Impose the Acoustic Sum Rule (ASR): After generating the dynamical matrices with ph.x and the force constants with q2r.x, impose the ASR when generating the final dispersion with matdyn.x. In the matdyn.x input, use asr = 'crystal' or asr = 'simple' [2] [13].
• 2D Systems: For monolayers and other 2D materials, always set assume_isolated = '2D' in the &SYSTEM namelist of the initial pw.x calculation to correctly handle the long-range interactions and avoid imaginary acoustic frequencies near the gamma point [13].

The Scientist's Toolkit: Essential Input Parameters

This table provides a quick reference for the key parameters discussed in this guide.

Input File	Parameter	Function	Recommendation
pw.x	ecutwfc / ecutrho	Kinetic energy cutoff for wavefunctions/charge density [25].	Converge first with SCF calculations; use higher values for phonons [13].
	conv_thr	SCF convergence threshold for the electronic energy [25].	Use a tight threshold (e.g., 1.0d-12).
	forc_conv_thr	Convergence threshold for ionic forces during relaxation [25].	Set to a stringent value (e.g., 1.0d-5 Ry/Bohr).
	k-point grid	Defines sampling of the Brillouin zone [25].	Test grids of increasing density until energy stabilizes.
ph.x	tr2_ph	SCF convergence threshold for the phonon calculation [10].	Set lower than pw.x conv_thr (e.g., 1.0d-14).
	nq1, nq2, nq3	Defines the q-point grid for phonon calculations when ldisp=.true. [10].	Determines the mesh for calculating dynamical matrices.
	nmix_ph	Number of iterations used in potential mixing [10].	Increase (8-20) to speed up convergence.
matdyn.x	asr	Applies the Acoustic Sum Rule to the force constants [2].	Use 'crystal' or 'simple' to fix acoustic modes at Γ.

Frequently Asked Questions (FAQs)

My calculation stopped with "occupation numbers probably wrong." What should I do? This warning appears for metallic or spin-polarized systems when the occupation type is not set to smearing. In your pw.x input, inside the &SYSTEM namelist, set occupations = 'smearing' and select an appropriate smearing function (e.g., smearing = 'gaussian') and a small broadening value (degauss) [5] [8].

I get "wrong degeneracy" or other symmetry-related errors. How can I fix them? These errors often occur when atomic positions are almost, but not exactly, at high-symmetry positions. The solution is to define your system using the correct ibrav value (not 0) and, if known, use Wyckoff positions in the initial pw.x calculation. This ensures the code correctly identifies all symmetry operations [5].

After a calculation fails, ph.x cannot restart and gives errors about recover files. What is the solution? This indicates corrupted restart files from the previous failed run. Remove all files named recover* in your outdir directory, and then restart the calculation [5].

For a 2D material, I get negative frequencies near the gamma point even after convergence. What else can I try? Ensure you have set assume_isolated = '2D' in the &SYSTEM namelist of your pw.x input file. This is crucial for correctly handling the long-rangeelectrostatics in low-dimensional systems and often resolves issues with acoustic modes [13].

A technical guide for computational researchers tackling unstable phonon calculations

Why Your Phonon Calculations Have Negative Frequencies

Negative frequencies (or imaginary frequencies) in Quantum ESPRESSO phonon calculations are not merely numerical noise; they are a fundamental indicator of a mechanical instability in your simulated structure. This means the atomic configuration you provided is not in a stable energy minimum, causing the phonon code to detect modes that lower the system's energy, which manifest as imaginary frequencies [5]. For researchers, especially in drug development where accurate biomolecular simulations are critical, this instability can invalidate subsequent property predictions.

The primary culprit is almost invariably inadequate structural relaxation [4] [5]. If the initial structure is not sufficiently relaxed—meaning the interatomic forces and internal stresses are not minimized to a tight enough threshold—the resulting "equilibrium" structure is not truly at equilibrium. When the phonon code (ph.x) constructs the dynamical matrix based on this pseudo-equilibrium state, it correctly identifies these residual instabilities [5]. Other contributing factors include insufficient SCF convergence criteria (conv_thr, tr2_ph), an inadequate plane-wave energy cutoff (ecutwfc, ecutrho), or a poorly converged k-point grid [5].

---

Troubleshooting Guide: Resolving Structural Relaxation Issues

Diagnosing the Problem

Before proceeding, confirm that the negative frequencies are an artifact of poor relaxation and not a real property of your system. Cross-check the following:

• Check Specific q-points: Run ph.x for a few high-symmetry q-points individually. If these show only positive frequencies, but the band structure from matdyn.x shows negatives, the problem likely lies with the force constants generated by q2r.x, pointing to a global relaxation issue [4].
• Verify the Acoustic Sum Rule (ASR): For the Gamma point (q=0), the acoustic modes should have frequencies very close to zero. A significant deviation (e.g., > 10 cm⁻¹) is a strong sign of Acoustic Sum Rule violation due to inadequate relaxation or numerical approximations in the code [5]. You can test this by using the dynmat.x program to impose the ASR and see if the acoustic modes correct themselves.

Step-by-Step Remediation Protocol

1. Tighten Relaxation Convergence Criteria The most critical step is to ensure your initial structure is perfectly relaxed. For a vc-relax (variable cell relaxation) calculation, use the following pw.x input parameters as a starting point for a high-quality relaxation [4] [7]:

Parameter	Recommended Value	Purpose
forc_conv_thr	1.0e-5 to 1.0e-4 Ry/Bohr	Convergence threshold for forces [7].
etot_conv_thr	1.0e-6 to 1.0e-5 Ry	Convergence threshold for total energy [7].
conv_thr	1.0e-8 to 1.0e-10 Ry	SCF cycle convergence threshold [7].
cell_dofree	'ibrav'	Allows cell to change shape according to its Bravais symmetry [7].

2. Optimize Key SCF and Phonon Parameters After a successful, tight relaxation, ensure the subsequent phonon calculation is also well-converged.

Parameter	Recommended Value	Function
ecutwfc / ecutrho	Increase by 20-30%	Higher cutoff improves accuracy of energy and force calculations [5].
tr2_ph	1.0e-14	Stricter convergence for the phonon SCF calculation [4] [13].
k-point grid	Use denser grid	Essential for metals and systems with complex Fermi surfaces [5].

3. Apply the Acoustic Sum Rule (ASR) Always impose the ASR in both q2r.x and matdyn.x calculations to correct for small violations of translational invariance.

4. Address System-Specific Issues

• For 2D Systems: If simulating 2D materials like graphene, add assume_isolated = '2D' in the SYSTEM namelist of your pw.x input to correctly handle long-range interactions and prevent imaginary frequencies near the Gamma point [13].
• For Metastable or Specific Materials: Be aware that some materials may exhibit genuine soft modes or instabilities with certain exchange-correlation functionals (like GGA). In such cases, using a different functional (e.g., LDA) or exploring meta-GGA/hybrid functionals might be necessary [4].

The following workflow diagram illustrates the integrated process of structural relaxation and phonon calculation, highlighting the critical steps to prevent negative frequencies.

---

Frequently Asked Questions (FAQs)

Q1: My ph.x calculation had positive frequencies, but matdyn.x shows negative ones. Why? This is a classic sign that the problem is with the overall force field, not a single q-point. The ph.x calculation on a discrete grid might be stable, but the Fourier interpolation performed by q2r.x and matdyn.x uses force constants from the entire grid. If the initial structure relaxation was inadequate, these force constants are inconsistent, leading to unphysical imaginary frequencies in the full phonon band structure [4].

Q2: How tight should the convergence criteria for relaxation be? There is no one-size-fits-all answer, but for reliable phonons, aim for forc_conv_thr of 1.0e-5 Ry/Bohr or tighter, and an SCF conv_thr of 1.0e-8 to 1.0e-10 Ry [4] [7]. The system is considered relaxed when the maximum force on any atom is below forc_conv_thr.

Q3: I have a negative frequency only at the Gamma point. What should I do? This is often an Acoustic Sum Rule (ASR) violation. First, ensure you are using asr = 'crystal' in your matdyn.x input. If the problem persists, it underscores that your structure is not perfectly relaxed. The ultimate test is to run dynmat.x to impose the ASR; if the acoustic modes then drop to nearly zero (< 1 cm⁻¹), you can trust the other modes, but should still improve the initial relaxation [5].

Q4: Are negative frequencies always a sign of a problem? Not necessarily. They can indicate a genuine structural instability, meaning the crystal structure you are simulating is metastable or unstable at the level of theory you are using (e.g., with GGA). It can suggest a phase transition. However, you must first rule out technical issues by following the troubleshooting steps above [5].

---

The Scientist's Toolkit: Essential Computational Reagents

This table lists the key "research reagents"—the core computational tools and parameters—within the Quantum ESPRESSO ecosystem that are essential for successful structural relaxation and phonon analysis.

Tool / Parameter	Function & Purpose
pw.x	The core plane-wave self-consistent field code used for structural relaxation (relax, vc-relax) and generating the ground-state electron density [13].
ph.x	The phonon code using Density-Functional Perturbation Theory (DFPT) to calculate the dynamical matrix on a grid of q-points [13].
q2r.x	Computes the interatomic force constants in real space by inverse Fourier transform of the dynamical matrices [13].
matdyn.x	Calculates the phonon dispersion across any path in the Brillouin zone and the phonon density of states by Fourier transforming the force constants [13].
forc_conv_thr	The primary tolerance in pw.x controlling structural relaxation quality; defines the target maximum residual force on any atom [7].
etot_conv_thr	The tolerance for the total energy change between successive relaxation steps [7].
tr2_ph	The convergence threshold for the self-consistency loop within the ph.x calculation [4] [13].
asr / zasr	The Acoustic Sum Rule imposition in matdyn.x (asr) and q2r.x (zasr) to enforce zero frequency for acoustic modes at the Gamma point [4] [13].

By systematically applying these diagnostic and remediation strategies, you can transform the frustrating occurrence of negative frequencies from a dead-end error into a valuable guide for achieving truly high-quality, physically meaningful simulation results.

A technical guide for resolving pervasive symmetry and convergence challenges in Quantum ESPRESSO phonon calculations

Understanding the Core Issues: A Troubleshooting Guide

Q1: What does the "wrong degeneracy" error mean and what causes it?

This error occurs during phonon calculations when the code detects that the degeneracy of phonon modes does not match what is expected from the crystal's symmetry operations. The root cause is almost invariably atomic positions that are close to, but not sufficiently close to, their exact symmetry positions [11] [5].

Even minute deviations—smaller than typical geometry optimization thresholds—can break the symmetry enough to trigger this error. Quantum ESPRESSO compares the rotated q-vector with the original with an acceptance tolerance of approximately 10⁻⁵ (set in PW/src/eqvect.f90). If your q-vector differs from a high-symmetry point by an amount on this order, the symmetry identification fails [5]. Using ibrav=0 (the free lattice parameter) exacerbates this problem because it doesn't enforce the intrinsic symmetry of your crystal system.

Q2: Why are my acoustic modes not zero at the Γ-point, and what is the Acoustic Sum Rule (ASR)?

The Acoustic Sum Rule embodies the translational invariance of the system, which requires that the sum of all atomic forces in any direction must be zero. In practice, this manifests as three phonon modes at the Γ-point (q=0) that should have zero frequency—the acoustic modes. When these frequencies are not zero, it indicates ASR violation due to numerical approximations in the calculation [11].

The primary source of this violation comes from the discreteness of the Fast Fourier Transform (FFT) grid [11]. The calculated frequency of the acoustic mode is typically less than 10 cm⁻¹, but in problematic cases can reach up to 100 cm⁻¹ [5]. This problem is generally more severe for GGA functionals than LDA because GGA functionals have more strongly varying forms, particularly in systems with significant vacuum portions [5].

Q3: What is the relationship between "negative frequencies" and these issues?

"Negative" frequencies (actually imaginary frequencies where ω² ≤ 0) can stem from multiple sources [11]:

• For acoustic frequencies at Γ-point or rotational modes of molecules: These are fictitious effects of numerical inaccuracies (including ASR violation) or insufficient convergence [11].
• At other q-points or for optical modes: May signal a real mechanical instability of the chosen structure with the approximations used, but only after rigorously eliminating numerical causes [5].

Table: Troubleshooting Guide for Common Phonon Calculation Errors

Error Symptom	Primary Causes	Diagnostic Steps
"Wrong degeneracy" in star_q	Atomic positions imperfectly symmetric; q-vector slightly off high-symmetry point [5]	Verify atomic positions against Wyckoff positions; check q-vector values
Non-zero acoustic modes	ASR violation from discrete FFT grid; insufficient convergence [11] [5]	Check convergence of FFT grid (ecutwfc, ecutrho); enforce ASR in post-processing
"Negative" (imaginary) frequencies	Insufficient convergence; real mechanical instability; ASR violation [11] [5]	Methodically tighten convergence parameters; verify with ASR correction
Mysterious symmetry-related errors	Using ibrav=0 with near-symmetric atomic positions [11] [5]	Switch to correct ibrav value; use Wyckoff positions in SCF calculation

Practical Solutions and Implementation Protocols

Protocol: Resolving "Wrong Degeneracy" Errors

• Use the correct ibrav value: Avoid ibrav=0 in your self-consistent calculation. Instead, specify the correct Bravais lattice type for your crystal [11] [5].
• Employ Wyckoff positions: If known, use Wyckoff positions to define atomic coordinates in the input file rather than free coordinates [11] [5]. This ensures exact symmetry.
• Verify q-vectors: For calculations at specific q-points, ensure they exactly match high-symmetry points rather than being slightly offset [5].
• Symmetrize the dynamical matrix: In some cases, you can apply dyn.Symmetrize() to impose symmetry and ASR on a calculated dynamical matrix [27].

Protocol: Enforcing the Acoustic Sum Rule

The ASR can be imposed during the post-processing stages using q2r.x and matdyn.x:

• In q2r.x: Use the zasr variable with one of these options [28]:

  • 'simple': Previous implementation correcting diagonal elements
  • 'crystal': 3 translational ASR via optimized correction (projection)
  • 'one-dim': 3 translational + 1 rotational ASR
  • 'zero-dim': 3 translational + 3 rotational ASR

• In matdyn.x: Use the asr variable with similar options [14]:

  • 'crystal': 3 translational ASR via optimized correction
  • 'all': 3 translational + 3 rotational ASR + 15 Huang conditions (requires write_lr=.true. in q2r.x and read_lr=.true. in matdyn.x for infrared-active solids) [14]

• For infrared-active solids: When using asr='all', ensure long-range force constants are handled properly by setting write_lr=.true. in q2r.x and read_lr=.true. in matdyn.x [14].

Protocol: Comprehensive Convergence Checklist

To prevent both ASR violations and spurious imaginary frequencies:

• Tighten SCF convergence: Reduce conv_thr to at least 1.0e-10 (from the typical 1.0e-6) to lower noise in forces [12].
• Increase plane-wave cutoffs: Systematically converge both ecutwfc and ecutrho (typically ecutrho = 4-8 × ecutwfc for ultrasoft pseudopotentials) [5].
• Dense k-point grid: Especially crucial for metallic systems where convergence is slower [5].
• Tighten phonon convergence: Reduce tr2_ph to 1.0e-16 or lower for more accurate phonon calculations [12].
• Check FFT grid density: Higher charge density cutoffs increase grid density, reducing ASR violation [5].

Phonon Calculation Troubleshooting Workflow

Essential Computational Tools and Parameters

Table: Key Parameters for Resolving Symmetry and ASR Issues

Parameter/Utility	Function	Recommended Settings
ibrav	Defines Bravais lattice type	Use specific value rather than 0 [11] [5]
Wyckoff positions	Atomic coordinates respecting symmetry	Use when known instead of free coordinates [11] [5]
zasr (in q2r.x)	Enforces ASR during IFC generation	'crystal' for translations; 'one-dim'/'zero-dim' for rotations [28]
asr (in matdyn.x)	Enforces ASR during phonon interpolation	'crystal' or 'all' (with Huang conditions) [14]
conv_thr	SCF convergence threshold	1.0e-10 or tighter for forces [12]
tr2_ph	Phonon calculation convergence	1.0e-16 or tighter [12]
ecutwfc/ecutrho	Plane-wave cutoffs	Systematically converged; higher values reduce ASR violation [5]
dyn.Symmetrize()	Python function to impose symmetry/ASR	Use on dynamical matrices from output [27]

Advanced Considerations for Specific Cases

For metallic systems: Convergence with respect to k-point grid and smearing is notoriously slow when semicore states are present, particularly for phonon wave-vectors not commensurate with the k-point grid [5]. Additionally, the current implementation of electron-phonon coefficients works for metals only, not insulators [11].

For infrared-active solids: When enforcing the complete ASR with Huang conditions (asr='all'), you must enable the handling of long-range force constants by setting write_lr=.true. in q2r.x and read_lr=.true. in matdyn.x [14]. This ensures proper treatment of the non-analytical part that contributes to LO-TO splitting.

Low-symmetry crystals: To obtain q=0 results for different directions, specify twice q=0 in the q-point list [14]. The direction q̂ (as q→0) for the non-analytic part is extracted from the sequence of q-points as q̂ = q(n) - q(n-1) or q̂ = q(n) - q(n+1) [14].

Successfully addressing symmetry and precision issues requires methodical parameter convergence and appropriate use of symmetry enforcement tools. By implementing these protocols, researchers can distinguish genuine physical instabilities from numerical artifacts, laying the foundation for reliable lattice dynamics calculations in materials research.

A fundamental challenge in computational materials science, particularly within density-functional perturbation theory (DFPT) calculations using Quantum ESPRESSO's ph.x, is the occurrence of unphysical negative frequencies (imaginary phonon modes) in phonon dispersion spectra. These results signal a mechanical instability in the chosen structure under the employed computational parameters [5] [2]. For researchers investigating low-dimensional materials like 2D systems (e.g., monolayers, surfaces), this problem is especially prevalent near the gamma point (Γ) due to violations of the Acoustic Sum Rule (ASR) and inadequate treatment of long-range electrostatics in non-periodic dimensions [5] [29] [13]. This guide provides targeted troubleshooting and methodologies to systematically resolve these issues, ensuring physically meaningful phonon calculations essential for predicting material properties and stability.

Troubleshooting Guide: Q&A for Commonph.xIssues

General Phonon Calculation Problems

• Q1: My ph.x calculation stops with errors about reading files or recover files. What should I do?

  • A: This typically indicates issues with data files from the prior pw.x calculation or corrupted restart files.
    • Error reading file: Ensure the data file produced by pw.x is complete and was generated by a compatible version of the code. Verify that your outdir and prefix variables in the ph.x input correctly point to the SCF calculation results [5].
    • Cannot recover or error reading recover file: This suggests a corrupted restart file from a previously failed execution. The solution is to remove all files named recover* in your specified outdir directory and restart the calculation [5].

• Q2: ph.x reports "occupation numbers probably wrong" and continues. Is this a problem?

  • A: This warning appears for metallic or spin-polarized systems when the occupation type is not set to smearing. While the calculation may continue, it is not recommended to ignore this. For such systems, you should set occupations = 'smearing' in the pw.x input file to ensure proper treatment of partial occupancies [5].

• Q3: Why are my acoustic modes at Γ not exactly zero, and when should I worry?

  • A: A small deviation (typically < 10 cm⁻¹) is expected due to the numerical breaking of translational invariance on a discrete grid. However, frequencies much higher than this (up to 100 cm⁻¹) can occur, especially in systems with vacuum (like 2D materials) or when using GGA functionals. The definitive test is to diagonalize the dynamical matrix using dynmat.x while imposing the ASR. If this yields an acoustic mode with ω < 1 cm⁻¹ and leaves other modes unchanged, your results are trustworthy [5].

• Q4: My phonons are completely wrong, with severe negative frequencies, wrong symmetries, or gross ASR violations. What went wrong?

  • A: This severe problem can have multiple causes [5]:
    • Incorrect Structure: The initial structure may not be fully relaxed. Ensure your geometry optimization (e.g., vc-relax) is converged with tight thresholds for forces (forc_conv_thr) and total energy (etot_conv_thr) [2].
    • Poor Convergence: The convergence thresholds for the SCF (conv_thr in pw.x) or phonon (tr2_ph in ph.x) calculations might be too large. Try reducing them (e.g., tr2_ph=1.0d-14) [5] [2] [13].
    • Insufficient Cutoffs: The plane-wave cutoff for wavefunctions (ecutwfc) and, crucially, for charge density (ecutrho) might be too low. Increasing ecutrho (often 4-8 times ecutwfc) can help [5].
    • Sparse k-point Grid: The k-point grid used in the SCF calculation might be inadequate, particularly for metals. Using a denser grid can resolve this [5].
    • Wrong Masses: Incorrect atomic masses (amass) in the ph.x input will directly lead to wrong frequencies [5].

Specific Issues with ASR and 2D Systems

• Q5: What is the Acoustic Sum Rule (ASR), and which scheme should I use for my system?

  • A: The ASR enforces that the sum of the atomic forces in any periodic direction must be zero, which guarantees that the acoustic modes at Γ have zero frequency. The choice of scheme depends on your system's dimensionality and is critical for fixing spurious imaginary modes [5] [2] [13].

  The table below summarizes the common ASR schemes and their applications:

  Table: Acoustic Sum Rule (ASR) Schemes in Quantum ESPRESSO

  Scheme	Keyword	System Type	Description
  No Correction	'no'	Testing	Applies no correction. Not recommended for production.
  Simple 1D-3D	'simple'	Isotropic 3D crystals	Applies a universal correction. Good for standard 3D materials.
  Crystal	'crystal'	Anisotropic crystals	Respects the crystalline symmetry. Recommended for most bulk systems. [2] [13]
  2D Crystal	'2D'	2D materials (e.g., graphene, hBN)	Applies the sum rule only in the plane of the material. Essential for 2D systems. [13]

  These schemes are applied in q2r.x (via zasr) and/or matdyn.x (via asr).

• Q6: I am calculating phonons for a 2D material (like graphene or MoS₂). What special settings are required to avoid negative frequencies?

  • A: 2D materials require specific corrections to account for their anisotropic electrostatic environment [29] [13]:
    • In pw.x for SCF: You must set assume_isolated = '2D' in the SYSTEM namelist. This activates a different Poisson solver that correctly handles the non-periodic (Z) direction, which is crucial for avoiding unphysical interactions between periodic images. This is often the most critical step for fixing imaginary acoustic modes in 2D systems [13].
    • In ph.x: No special change is needed beyond ensuring a well-converged calculation.
    • In q2r.x and matdyn.x: Apply the appropriate 2D ASR. Set zasr = '2D' in q2r.x and asr = '2D' in matdyn.x [13].

• Q7: I get "Wrong degeneracy" or other symmetry-related errors. How can I fix them?

  • A: These errors almost always stem from atomic positions that are close to but not exactly at high-symmetry positions. The code's symmetry detection tolerance is about 10⁻⁵. To resolve this [5]:
    • Use the correct ibrav value in your pw.x input instead of ibrav=0.
    • If possible, use well-defined Wyckoff positions to generate atomic coordinates.
    • Ensure your initial structure is fully and accurately relaxed.

Frequently Asked Questions (FAQs)

• The phonon dispersion from matdyn.x still shows small imaginary frequencies after applying ASR. Is my structure unstable?

  • Small imaginary frequencies (< 10-20 cm⁻¹) near Γ after applying a 2D or crystal ASR are often a numerical artifact and can be considered "zero" for practical purposes. However, large imaginary frequencies at high-symmetry points in the Brillouin zone are a strong indicator of a genuine mechanical instability in the crystal structure [5].

• Can I use Grimme's DFT-D3 dispersion correction with phonon calculations?

  • Be cautious. As of the information available, an error stating "The phonon code with Grimme's DFT-D3 is not yet available" can occur. It is recommended to check the latest version's documentation for support. Alternatively, you can calculate the D3 Hessian matrix separately using d3hess.x [2].

• How important is convergence in the initial SCF calculation for accurate phonons?

  • It is critical. Phonon energies are sensitive to the ground-state charge density. Using a tighter SCF convergence threshold (conv_thr = 1.0d-12 or lower in pw.x) and a higher plane-wave cutoff (ecutwfc and ecutrho) is highly recommended for phonon calculations compared to a simple single-point energy calculation [5] [13].

• My calculation is running very slowly. What parameters can I adjust?

  • You can increase the mixing parameters in ph.x. Setting nmix_ph = 8 (or higher) can significantly speed up convergence at the cost of slightly higher memory usage. Adjusting alpha_mix can also help [10].

Experimental Protocols & Methodologies

Standard Protocol for Phonon Dispersion in 3D Bulk Materials

This workflow outlines the standard steps for a robust phonon calculation in bulk systems.

Diagram: Standard Phonon Calculation Workflow for 3D Bulk Materials

• Geometry Optimization (pw.x):

  • Perform a variable-cell relaxation (calculation = 'vc-relax').
  • Use tight convergence criteria, e.g., forc_conv_thr = 1.0d-4 Ry/bohr and etot_conv_thr = 1.0d-5 Ry.
  • Ensure the final structure is stress-free and has minimal forces on all atoms.

• SCF Calculation (pw.x):

  • Use the relaxed structure from step 1.
  • Set calculation = 'scf'.
  • Use a denser k-point grid and higher energy cutoffs (ecutwfc, ecutrho) than those used for a simple energy calculation [13].
  • Achieve high SCF convergence, e.g., conv_thr = 1.0d-12 or tighter.

• Phonon Calculation (ph.x):

  • Set ldisp = .true. and define a q-point mesh with nq1, nq2, nq3.
  • Use a tight phonon convergence threshold, e.g., tr2_ph = 1.0d-14 [2] [13].
  • Specify the output file for dynamical matrices with fildyn.

• Q to R Fourier Transform (q2r.x):

  • Read the dynamical matrices using fildyn.
  • Apply the Acoustic Sum Rule. For bulk materials, zasr = 'crystal' is typically the best choice [2] [13].
  • Output the force constants file (flfrc).

• Phonon Dispersion (matdyn.x):

  • Read the force constants from flfrc.
  • Set asr = 'crystal' to be consistent with q2r.x.
  • Define a high-symmetry path through the Brillouin Zone using q_in_band_form = .true..
  • Output the frequencies (flfrq) and, optionally, the eigenvectors (flvec).

Advanced Protocol for 2D Materials

This protocol modifies the standard workflow to correctly handle the unique electrostatics of 2D systems.

Diagram: Modified Workflow for 2D Materials Highlighting Critical Settings

The key differences from the 3D protocol are:

• In pw.x (SCF): The SYSTEM namelist must include assume_isolated = '2D' to employ a Poisson solver tailored for slab-like geometries [29] [13].
• In q2r.x: The ASR must be set to the 2D variant: zasr = '2D' [13].
• In matdyn.x: Similarly, the ASR is set to asr = '2D' [13].

The Scientist's Toolkit: Essential Input Parameters

Table: Key Input Parameters for Fixing Negative Frequencies in ph.x Calculations

Input Parameter	Program	Function	Recommended Value for Troubleshooting
assume_isolated	pw.x	Corrects long-range electrostatics for 0D, 1D, and 2D systems. Set to '2D' for monolayers and surfaces.	'2D'
asr / zasr	matdyn.x / q2r.x	Applies Acoustic Sum Rule to enforce zero acoustic modes at Γ. Choice depends on system dimensionality.	'crystal' (3D), '2D' (2D)
tr2_ph	ph.x	Convergence threshold for the phonon self-consistent cycle. Tighter thresholds improve accuracy.	1.0d-14 [2] [13]
conv_thr	pw.x	Convergence threshold for the electronic SCF cycle. A tighter threshold is crucial for accurate forces.	1.0d-12 or lower
ecutwfc / ecutrho	pw.x	Plane-wave kinetic energy cutoffs. Higher values improve basis set completeness, critical for phonons.	Increase by 20-30% from standard values [5]
nmix_ph	ph.x	Number of iterations for potential mixing. A higher value can significantly speed up convergence.	8 (default is 4) [10]
forc_conv_thr	pw.x	Force convergence threshold in geometry optimization. Ensures the initial structure is truly stable.	1.0d-4 Ry/bohr or tighter [2]

Result Validation and Interpretation: Confirming Physical Meaning and Accuracy

Frequently Asked Questions

What does a "negative frequency" in my phonon calculation mean? A negative frequency, more accurately an imaginary frequency, indicates a dynamical instability in the calculated structure. This means the atomic configuration you started with is not a true local minimum on the potential energy surface and may undergo a phase transition to a more stable structure [2].

My calculation was running, but now ph.x has stopped. What should I check? First, verify that the data file produced by your preceding pw.x calculation is good, complete, and was generated by a compatible version of the code. An incomplete or corrupted data file is a common cause for ph.x to stop with an error [5].

Can I perform a phonon calculation at just the Gamma-point? Yes, you can calculate phonons at the Gamma-point only using Density Functional Perturbation Theory (DFPT). This is perfectly reasonable for studying specific properties like Raman spectra. However, if you need phonon dispersions or properties at other q-points, you must converge your calculations with respect to a full q-point grid [15].

Why are my acoustic modes at the Gamma-point not exactly zero? A small, non-zero frequency (typically <10 cm⁻¹) for acoustic modes at Γ is normal and results from the Acoustic Sum Rule (ASR) not being exactly verified due to the discrete real-space grid used in the calculation. You can diagonalize the dynamical matrix with dynmat.x while imposing the ASR to obtain corrected, near-zero frequencies [5].

Troubleshooting Guide

Issue 1: Negative or "Lousy" Phonon Frequencies

Problem: The phonon calculation completes but yields negative frequencies, frequencies with wrong symmetries, or gross Acoustic Sum Rule violations [5] [2].

Troubleshooting Step	Action and Purpose
Verify Structural Relaxation	Ensure the atomic structure is fully relaxed to the ground state. A poor or incomplete relaxation (e.g., residual forces or stresses) is a primary cause of imaginary frequencies. Re-run vc-relax with tight convergence thresholds (forc_conv_thr, etot_conv_thr) [2].
Check Convergence Parameters	Systematically tighten key convergence parameters in both the SCF (pw.x) and phonon (ph.x) calculations. The table below provides a starting point for investigation [5] [2].
Inspect Pseudopotentials	Ensure all pseudopotentials are generated using the same flavor of Density Functional Theory (DFT). Using inconsistent pseudopotentials can cause the code to stop [17].
Apply Acoustic Sum Rule	Impose the ASR in post-processing (e.g., in q2r.x or matdyn.x using asr='crystal') to correct the acoustic modes at the Gamma point [2].

Issue 2:ph.xStops with "occupation numbers probably wrong"

Problem: The code issues a warning about occupation numbers and continues. This typically occurs in metallic or spin-polarized systems [5].

• Solution: In your pw.x input, set occupations='smearing' and choose an appropriate smearing function (e.g., 'gaussian' or 'mv'). Avoid using occupations='fixed' for metallic systems [17] [5].

Issue 3: "Wrong degeneracy" or Other Symmetry Errors

Problem: Errors related to symmetry operations occur [5].

• Solution: Atomic positions are likely close to, but not exactly at, high-symmetry positions. To fix this:
  • In your pw.x input, define the Bravais lattice using the correct ibrav value instead of ibrav=0.
  • Use Wyckoff positions to generate atomic coordinates if known.
  • Ensure your initial structure is precise.

Experimental Protocols and Parameters

Protocol 1: Systematic Convergence Checklist

Before concluding a structure is unstable, methodically check the convergence of these parameters. The following table summarizes key parameters to test, drawing from a real example of troubleshooting monolayer InSe [2].

Parameter	Description	Function in Calculation	Investigative Range
ecutwfc	Plane-wave kinetic energy cutoff for wavefunctions.	Determines the basis set size and computational cost.	Increase incrementally (e.g., from 60 Ry to 80 Ry).
ecutrho	Cutoff for charge density and potential (often 4x or 8x ecutwfc).	Critical for ultrasoft pseudopotentials (USPP) and PAW.	Ensure it is a sufficient multiple of ecutwfc (e.g., 240 Ry to 480 Ry) [2].
k-point Grid	Grid of points in the Brillouin zone for electronic integration.	Samples electron states; crucial for metals.	Use a denser grid (e.g., from 12x12x1 to 24x24x1 for 2D) [2].
q-point Grid	Grid for phonon calculation in the Brillouin zone.	Samples phonon wavevectors.	Test convergence with a denser grid (e.g., 8x8x1) [15].
conv_thr	Convergence threshold for the SCF cycle.	Stops the electronic minimization when energy change is below this value.	Tighten significantly (e.g., from 1.0e-6 to 1.0e-12 or 1.0e-16) [2].
tr2_ph	Convergence threshold for the phonon calculation.	Stops the phonon self-consistency cycle.	Tighten (e.g., to 1.0e-14 or lower) [2].
mixing_beta	Mixing parameter for SCF charge density.	Aids convergence of difficult SCF calculations.	Adjust if SCF is not converging (e.g., 0.2 to 0.7) [2].

Protocol 2: Sample Input Parameters from a Monolayer InSe Study

The following input blocks are adapted from a published case study on monolayer InSe, which encountered negative frequencies. Use them as a reference for structuring your own calculations [2].

Sample pw.x Input (SCF Calculation)

Atomic species and positions follow as in the original input file [2].

Sample ph.x Input

The Scientist's Toolkit: Research Reagent Solutions

Item	Function in Calculation
Pseudopotentials	Replace the core electrons of an atom with a potential, simplifying the calculation while retaining chemical accuracy. All pseudopotentials in a single calculation must be generated using the same DFT functional [17].
Acoustic Sum Rule (ASR)	A physical rule requiring the frequencies of acoustic modes at the Gamma point (q=0) to be zero. It can be imposed in post-processing (e.g., in q2r.x or matdyn.x) to correct small numerical errors [5] [2].
Density Functional Perturbation Theory (DFPT)	The mathematical framework used by ph.x to compute phonon frequencies and other response properties by considering the linear response of the electronic system to a displacement perturbation [15].
k-points and q-points	k-points: Sample the electronic Brillouin zone in the SCF calculation. q-points: Sample the phonon Brillouin zone. Both grids must be converged for accurate results [15].
Smearing Functions	(e.g., Gaussian, Marzari-Vanderbilt) A technique to assign partial occupations to electron states near the Fermi level, which is essential for achieving SCF convergence in metallic systems [17].

Workflow for Troubleshooting Negative Frequencies

The following diagram outlines a systematic workflow for diagnosing and fixing negative phonon frequencies.

A technical guide for Quantum ESPRESSO users troubleshooting phonon calculations

Why do my phonon calculations show acoustic modes not starting at zero, and how candynmat.xhelp?

Problem: After a ph.x calculation, the acoustic modes at the gamma point (q=0) do not exhibit zero frequency, which is expected due to the system's translational invariance. This indicates a violation of the Acoustic Sum Rule (ASR).

Solution: The dynmat.x tool can be used to re-diagonalize the dynamical matrix while imposing different types of ASR, which can force the acoustic modes to zero frequency. The ultimate test is to diagonalize the dynamical matrix with dynmat.x, imposing the ASR. If you obtain an acoustic mode with a much smaller ω (let us say < 1cm⁻¹) with all other modes virtually unchanged, you can trust your results [5].

How to Usedynmat.xfor ASR Verification

The dynmat.x program reads a dynamical matrix file produced by the phonon code, can apply a chosen Acoustic Sum Rule, and diagonalizes the dynamical matrix [30].

Input Parameters for ASR Correction

The key to fixing the acoustic modes lies in the &INPUT namelist of the dynmat.x input file. Critical variables for ASR include [30]:

Parameter	Type	Default	Description & Allowed Values
fildyn	Character	'matdyn'	Input file containing the dynamical matrix (e.g., 'prefix.dyn' from ph.x).
asr	Character	'no'	Type of Acoustic Sum Rule: 'no', 'simple' (3 translational ASR), 'crystal' (3 translational, optimized), 'one-dim' (3 trans. + 1 rotational), 'zero-dim' (3 trans. + 3 rotational).
axis	Integer	3	For 1D systems (asr='one-dim'), indicates the rotation axis (1=Ox, 2=Oy, 3=Oz).

Sampledynmat.xInput File

This example, from a user calculating a one-dimensional system (silicon nanowire), applies a one-dimensional ASR including rotational invariance around the x-axis [31].

Types of Acoustic Sum Rules (ASR) and Their Applications

The choice of ASR depends on the dimensionality of your system. The following table summarizes the options [30] [32]:

ASR Type (asr parameter)	Translational ASR	Rotational ASR	Recommended Use
'no'	No	No	Testing purposes only.
'simple'	Yes (3)	No	Bulk, 3D periodic systems. Corrects diagonal elements.
'crystal'	Yes (3)	No	Bulk crystals. Uses an optimized correction (projection).
'one-dim'	Yes (3)	Yes (1)	1D systems (nanowires, nanotubes). Axis set by axis parameter.
'zero-dim'	Yes (3)	Yes (3)	0D systems (molecules, clusters).

Experimental Workflow for ASR Verification

The following diagram illustrates the recommended workflow to diagnose and correct acoustic sum rule violations in your phonon calculations, positioning dynmat.x as a verification and correction tool.

The Scientist's Toolkit: Essentialdynmat.xInputs

Item	Function in Experiment	Technical Specification
Dynamical Matrix File	Primary input for dynmat.x. Contains the force constants from the ph.x DFPT calculation.	Specified via fildyn parameter (e.g., 'prefix.dyn' or 'prefix.dyn1' for Gamma-point) [30].
ASR Type	Determines the mathematical correction applied to enforce translational/rotational invariance.	Choose from 'simple', 'crystal', 'one-dim', or 'zero-dim' via the asr parameter [30].
Atomic Masses (amass)	Critical for correct frequency calculation. Masses are used when diagonalizing the dynamical matrix.	Must be specified in atomic mass units for each atom type. Should match values used in ph.x [30].
Q-direction (q)	Defines the direction for calculating non-analytical contributions (LO-TO splitting) along a specific Cartesian direction.	A real vector q(1), q(2), q(3). Required for proper treatment of longitudinal and transverse optical modes [30] [33].

Troubleshooting FAQ

Q1: My dynmat.x calculation fixed the acoustic modes, but matdyn.x still shows negative frequencies at gamma. Why? [31]

A: This is a known phenomenon, especially in low-dimensional systems. While dynmat.x diagonalizes the dynamical matrix from ph.x directly and can perfectly impose the ASR for that specific point, matdyn.x performs a Fourier interpolation of force constants from q2r.x to get the dynamical matrix at any q-point. Small numerical inaccuracies in the force constants can be amplified during this interpolation, leading to residual ASR violations and negative frequencies at interpolated q-points, including gamma. Using the same, appropriate asr keyword in both dynmat.x and matdyn.x/q2r.x inputs is crucial for consistency.

Q2: What if imposing the ASR with dynmat.x does not fix my negative frequencies, or they are large and not just near the gamma point? [5]

A: Large or persistent negative frequencies (imaginary frequencies) often signal a deeper problem beyond mere ASR violation. This can indicate a mechanical instability in your simulated structure. You should check:

• Structure Relaxation: Ensure your geometry is fully optimized (vc-relax in pw.x) with tight convergence thresholds for forces and stresses [2].
• SCF Convergence: Use stricter convergence thresholds in both your SCF (pw.x) and phonon (ph.x) calculations (e.g., conv_thr = 1.0d-12 or lower in pw.x, tr2_ph = 1.0d-14 or lower in ph.x) [2] [34].
• Computational Parameters: Verify your plane-wave cutoff (ecutwfc, ecutrho) and k-point grid are sufficiently dense [5].
• 2D Systems: For 2D materials, use assume_isolated = '2D' in the SYSTEM namelist of your pw.x input to avoid spurious interactions between periodic images [13].

Key Technical Takeaways

• Primary Function: The dynmat.x tool is essential for verifying and enforcing the Acoustic Sum Rule (ASR) in phonon calculations, particularly for correcting non-zero acoustic modes at the gamma point [30] [5].
• Critical Inputs: Successful application requires proper specification of the asr parameter tailored to system dimensionality, accurate atomic masses, and correct dynamical matrix file reference [30].
• Diagnostic Value: A significant reduction in acoustic mode frequencies (to < 1 cm⁻¹) after dynmat.x ASR application, with other modes remaining unchanged, validates the underlying phonon calculation [5].
• Limitation Awareness: While dynmat.x can correct ASR violations, persistent large negative frequencies often indicate fundamental issues like inadequate structure relaxation or poor convergence that require addressing at the pw.x or ph.x level [2] [5].

Frequently Asked Questions

What is the most common cause of negative frequencies in ph.x calculations? Negative frequencies often indicate a mechanical instability in the calculated structure. The most common causes are your structure not being sufficiently relaxed or an issue with the choice of exchange-correlation functional (LDA vs. GGA) [5] [4]. It is critical to first ensure your structure is fully optimized with tight convergence criteria for forces and energies before starting a phonon calculation.

My ph.x calculation finished with positive frequencies, but matdyn.x shows negative frequencies in the band structure. Why? This discrepancy can arise from how the force constants are interpolated between q-points. While the dynamical matrices calculated directly by ph.x at specific q-points may be fine, the process of generating the force constant file (q2r.x) and interpolating to new q-points (matdyn.x) can introduce instabilities if the initial supercell used for the q-grid was too small, or if the Acoustic Sum Rule (ASR) is not properly imposed in all steps [5] [4]. Ensure you use the same ASR (e.g., zasr='crystal') in both q2r.x and matdyn.x inputs.

How do LDA and GGA typically differ in predicting phonon frequencies? Generally, GGA functionals tend to produce lattice constants that are slightly larger and more accurate than LDA, which typically over-binds and underestimates them. This directly impacts the calculated phonon spectra. While both can sometimes underestimate phonon frequencies, GGA often leads to a systematic underestimation of bulk moduli and phonon frequencies compared to LDA [35]. In some cases, such as for rutile TiO₂, GGA is known to predict unstable (negative) phonons where LDA does not, indicating a genuine sensitivity of the predicted ground state to the functional [4].

Troubleshooting Guide: Resolving Negative Frequencies

Verify Structural Relaxation

The ground-state structure from your pw.x calculation is the foundation for phonons. If it is not fully relaxed, the forces on atoms will not be zero, leading to imaginary phonons.

• Symptom: Negative frequencies appear, especially at the gamma point (q=0).
• Solution: Re-run your structure relaxation (pw.x) with tighter convergence thresholds.
• Recommended Protocol:
  • forc_conv_thr = 1.0d-5 (Ry/bohr)
  • etot_conv_thr = 1.0d-8 (Ry)

Impose the Acoustic Sum Rule (ASR)

The ASR states that the sum of all atomic forces in a unit cell must be zero, which guarantees that acoustic phonons have zero energy at the gamma point. Numerical errors can break this invariance.

• Symptom: Acoustic modes at Γ point have non-zero (and sometimes significantly large) frequencies [5].
• Solution: Impose the ASR when post-processing your results.
• Recommended Protocol:
  • In your q2r.x input, add: zasr = 'crystal'
  • In your matdyn.x input, add: asr = 'crystal'
  • You can then use dynmat.x to diagonalize the dynamical matrix with the imposed ASR to verify the results [5].

Improve SCF and Phonon Convergence

Insufficient convergence in either the self-consistent field (SCF) cycle of the ground state or the iterative solution in ph.x can lead to inaccurate forces and dynamical matrices.

• Symptom: Phonon frequencies seem "lousy" or change significantly when you tighten convergence [5].
• Solution: Systematically reduce the convergence thresholds.
• Recommended Protocol:
  • SCF (pw.x): conv_thr = 1.0d-11 (Ry)
  • Phonon (ph.x): tr2_ph = 1.0d-14

Check K-point Grid and Smearing for Metals

For metallic systems, convergence is notoriously slow. Using a coarse k-point grid or inappropriate smearing can cause problems.

• Symptom: Severe negative frequencies or incorrect phonon dispersions in metals.
• Solution: Increase the k-point grid density and ensure you are using an appropriate smearing function and width [5] [36].
• Recommended Protocol:
  • Always perform a k-point convergence test for the ground-state energy.
  • In your pw.x input, use: occupations = 'smearing' and choose a smearing (e.g., smearing = 'gaussian') with a suitable degauss value.

Address Functional-Specific Issues (LDA vs. GGA)

The choice between LDA and GGA can be decisive. If your system is known to have instabilities with one functional, switching may be necessary.

• Symptom: All other checks pass, but significant negative frequencies persist.
• Solution: Consult literature for your specific material. If GGA gives negatives and LDA does not (or vice versa), it may indicate a real physical instability or a known limitation of the functional [4] [35].
• Recommended Protocol:
  • Re-run the entire calculation (relaxation + phonons) with a different functional (e.g., switch from GGA to LDA).
  • Consider using a more advanced functional beyond standard GGA if the problem persists.

Comparative Performance: LDA vs. GGA

The table below summarizes general trends and performance characteristics of LDA and GGA in predicting properties relevant to phonon calculations, based on documented studies [37] [35].

Property	LDA (Local Density Approximation)	GGA (Generalized Gradient Approximation)	Experimental Reference & Notes
Lattice Constant	Typically underestimated (over-binding) [35]	Generally more accurate or slightly overestimated [35]	For NaI, exp. ~6.47 Å; LDA ~6.39 Å; GGA ~6.59 Å [37]
Bulk Modulus	Typically overestimated [35]	Often underestimated [35]	Softer lattices in GGA contribute to lower phonon frequencies.
Cohesive Energy	Overestimated [35]	Closer to experimental values [35]	GGA total energies are generally more accurate.
Phonon Frequencies	Often overestimated (stiffer bonds)	Often underestimated (softer bonds) [35]	For Si, GGA (PW) underestimates frequencies vs. LDA and experiment [35].
Structural Stability	Can incorrectly stabilize some phases	Can predict instabilities (negative frequencies) where LDA does not [4]	E.g., Rutile TiO2 phonons are unstable with GGA but not LDA [4].

Experimental Protocol: Comparing LDA and GGA for Phonons

This protocol provides a step-by-step methodology for a systematic comparison of LDA and GGA performance in predicting phonon properties for a material like Sodium Iodide (NaI) [37].

1. Computational Setup

• Code: Quantum ESPRESSO (PWscf for SCF, ph.x for phonons) [37].
• Pseudopotentials: Use norm-conserving or ultrasoft pseudopotentials from a standard library (e.g., PSLibrary). Use the same pseudopotential type (e.g., SG15) for both LDA and GGA, only changing the functional designation.
• System: NaI in the rocksalt (B1) crystal structure. Wyckoff positions: Na at (0,0,0), I at (0.5,0.5,0.5) [37].

2. Ground-State Self-Consistent Field (SCF) Calculation

• Objective: Obtain the converged charge density and total energy.
• Input Parameters (pw.x):
  • calculation = 'scf'
  • ecutwfc = 50 (Ry) /* Convergence testing required /
  • ecutrho = 400 (Ry) / Often 8-10 x ecutwfc /
  • occupations = 'smearing' / Critical for metals /
  • smearing = 'gaussian'
  • degauss = 0.01 (Ry)
  • conv_thr = 1.0d-11 (Ry)
  • forc_conv_thr = 1.0d-5 (Ry/bohr) / For relaxation /
  • etot_conv_thr = 1.0d-8 (Ry) / For relaxation */

3. Structural Relaxation

• Objective: Find the theoretical equilibrium lattice constant and atomic positions for each functional.
• Input Parameters (pw.x):
  • calculation = 'vc-relax' /* For full cell and ionic relaxation */
  • Use the same convergence criteria as above.
• Output Analysis: Plot energy vs. volume for a range of volumes and fit to an equation of state (e.g., Murnaghan) to extract the equilibrium lattice constant (a₀) and bulk modulus (B) for both LDA and GGA.

4. Phonon Calculation via DFPT

• Objective: Calculate the dynamical matrices on a coarse q-grid.
• Input Parameters (ph.x):
  • tr2_ph = 1.0d-14
  • ldisp = .true.
  • nq1 = 4, nq2 = 4, nq3 = 4 /* Set a suitable q-grid */
  • outdir = './'
  • fildyn = 'prefix.dyn'

5. Post-Processing for Phonon Dispersion

• Objective: Obtain the full phonon dispersion curve along high-symmetry paths.
  • a. Generate Force Constants: Use q2r.x to create the force constant file from the dynamical matrices. Input: &INPUT fildyn='prefix.dyn' flfrc='prefix.fc' zasr='crystal' /
  • b. Interpolate Phonons: Use matdyn.x to calculate frequencies along a path. Input: &INPUT asr='crystal' flfrc='prefix.fc' flfrq='prefix.freq' flvec='prefix.modes' /

6. Data Analysis and Comparison

• Compare the calculated LDA and GGA lattice constants and bulk moduli against experimental values.
• Plot the LDA and GGA phonon dispersion curves on the same graph for visual comparison of frequencies and any soft modes.
• Calculate the phonon density of states (PhDOS) and compare specific mode frequencies (e.g., transverse optic - TO - mode at Γ).

The Scientist's Toolkit: Research Reagent Solutions

Tool / Reagent	Function in Experiment
Quantum ESPRESSo Suite	Integrated suite of open-source computer codes for electronic-structure calculations and materials modeling at the nanoscale. Core components are pw.x (SCF/relax) and ph.x (phonons) [37].
Norm-Conserving Pseudopotentials	A type of pseudopotential that preserves the charge density outside the core region. Used to represent the core electrons, making plane-wave DFT calculations computationally feasible [37].
PBE-GGA Functional	A specific and very popular parameterization of the GGA functional. Often used as a standard for comparison against LDA [37].
XCrysDen	A program for visualizing crystal structures and densities. Useful for visualizing your input structure and verifying atomic positions [37].
dynmat.x Code	A tool in the Quantum ESPRESSO distribution used to diagonalize the dynamical matrix and analyze phonon modes, with the option to impose the Acoustic Sum Rule (ASR) post-calculation [5].

Experimental Workflow for Phonon Calculation

The diagram below outlines the critical steps for performing and troubleshooting a phonon calculation, highlighting key decision points for dealing with negative frequencies.

Frequently Asked Questions

What does an imaginary frequency in my phonon calculation mean?

An imaginary frequency (often reported as a negative value in output files) indicates that your system is dynamically unstable at the calculated level of theory [5] [38]. Mathematically, it arises from a negative eigenvalue of the dynamical matrix, which corresponds to a direction in the potential energy surface where the energy decreases, not increases [38]. This can be either a computational artifact or the signature of a genuine soft mode.

My ph.x calculation had positive frequencies, but matdyn.x shows imaginary ones. Why?

This discrepancy often points to issues with the interpolation process [4]. The ph.x program calculates the dynamical matrix on a discrete grid of q-points. The q2r.x and matdyn.x programs then use these to build a real-space force constant matrix and interpolate it to any q-point for the band structure. If the original q-point grid in the ph.x calculation is too coarse, the interpolation can introduce inaccuracies that manifest as imaginary frequencies, even if the frequencies at the calculated points were real [4] [13].

How can I distinguish a computational error from a real soft mode?

A systematic approach is needed to distinguish between them. The table below outlines common causes and their solutions.

Cause	Description	Solution
Insufficient Structural Relaxation	The atomic positions, cell volume, and shape are not at a true energy minimum, so forces or stresses remain [4].	Tighten convergence criteria in the pw.x relaxation (forc_conv_thr, etot_conv_thr) [4].
Poor SCF Convergence	The self-consistent field (SCF) calculation for the ground-state electron density is not fully converged [5].	Reduce conv_thr in the pw.x input and increase tr2_ph in the ph.x input [5].
Insufficient k-points or q-points	The sampling of the Brillouin zone is too sparse, failing to capture the electronic or vibrational structure accurately [5].	Use a denser k-point grid in the pw.x SCF calculation and a denser q-point grid (nq1, nq2, nq3) in the ph.x calculation [4] [13].
Inadequate Plane-Wave Cutoff	The energy cutoffs (ecutwfc, ecutrho) are too low, leading to an inaccurate description of the wavefunctions and charge density [5].	Systematically increase ecutwfc and ecutrho until total energy and phonon frequencies converge.

What does "following the imaginary mode" mean in practice?

It is a systematic procedure to find a lower-energy structure. The eigenvector of the dynamical matrix associated with the imaginary phonon describes the pattern of atomic displacements that lowers the energy. You "follow" this mode by:

• Identifying the wavevector (q-point) and branch (ν) of the imaginary mode.
• Extracting its eigenvector, which defines the displacement pattern.
• Generating new structures by displacing the atoms from their high-symmetry positions along this eigenvector with a scaling factor α: S_α = S_reference + α * U_mode [38].
• Calculating the total energy of each displaced structure S_α using pw.x. Plotting energy versus α will reveal a potential well, and the minimum of this well corresponds to the new, lower-symmetry, stable structure [38].

Troubleshooting Guide: Fixing Negative Frequencies

Step 1: Verify and Improve the Initial Structure

A structurally sound input is the most critical factor for obtaining physically meaningful phonons.

• Action: Re-run your geometry optimization (pw.x vc-relax or relax) with tighter convergence thresholds. For example, set forc_conv_thr to 1.0d-5 Ry/Bohr or lower and etot_conv_thr to 1.0d-6 Ry or lower [4].
• Action: For systems with vacuum (like surfaces or 2D materials), specify assume_isolated = '2D' in the SYSTEM namelist of your pw.x input to correctly handle the long-range Coulomb interactions and prevent spurious imaginary acoustic modes [13].

Step 2: Ensure Computational Parameters are Converged

Inaccurate computational parameters are a common source of unphysical instabilities.

• Action: Perform a systematic convergence test for the plane-wave energy cutoffs (ecutwfc and ecutrho) and the k-point grid. Phonon properties often require higher cutoffs than those needed for ground-state energy convergence [5] [13].
• Action: In your ph.x input, use a denser q-point grid by increasing nq1, nq2, and nq3. A 4x4x4 grid is often a minimum starting point for simple systems, with 6x6x6 or 8x8x8 being more reliable [4] [13].
• Action: Tighten the convergence thresholds for both the electronic (pw.x) and phonon (ph.x) calculations. In your pw.x input, set conv_thr = 1.0d-10 or lower. In your ph.x input, set tr2_ph = 1.0d-14 or lower [4] [5].

Step 3: Apply the Acoustic Sum Rule (ASR)

The ASR ensures that the frequencies of the three acoustic modes at the Brillouin zone center (Γ-point) are exactly zero, a consequence of translational invariance.

• Action: Both q2r.x and matdyn.x have an asr (or zasr) input flag. Apply the 'crystal' ASR to clean up small imaginary frequencies in the acoustic modes [4] [13]. If you have imaginary frequencies only for the acoustic modes at Gamma and they are small (e.g., < 20 cm⁻¹), applying the ASR is usually sufficient [5].

Case Study: Rutile and the High-Pressure Ice Phase Transition

The Rutile Computational Problem

A user on the Matter Modeling Stack Exchange reported a classic problem: their ph.x calculation for rutile showed positive frequencies at computed q-points, but the interpolated band structure from matdyn.x contained negative frequencies [4]. This is a strong indicator of an interpolation error due to an insufficiently coarse q-point grid in the initial DFPT calculation. The solution was to re-run the ph.x calculation with a denser grid (nq1=6, nq2=6, nq3=6 or higher) to capture the force constants more accurately before interpolation [4] [13].

Genuine Soft Mode in High-Pressure Ice

A study of high-pressure ice phases (VII, VIII, and X) provides a clear example of a physically real soft mode [39]. Theoretical calculations using classical molecular dynamics failed to reproduce the experimentally observed phase transition from proton-ordered ice VIII to proton-disordered ice VII at ~60 GPa and 200 K. However, when nuclear quantum effects (proton tunneling and zero-point motion) were included via ab initio path-integral molecular dynamics (PIMD), the calculations successfully showed the transition [39]. The "softening" of the potential energy barrier, driven by quantum tunneling, was the physical origin of the instability. In this case, the imaginary mode (if calculated in the higher-symmetry phase) pointed directly to a real, physically meaningful phase transition.

The Scientist's Toolkit

Research Reagent / Computational Tool	Function / Explanation
ph.x	The main Quantum ESPRESSO program for calculating phonons using Density Functional Perturbation Theory (DFPT) [13].
q2r.x	Computes the interatomic force constants in real space by performing an inverse Fourier transform of the dynamical matrices generated by ph.x [4] [13].
matdyn.x	Calculates the phonon dispersion across the entire Brillouin zone by Fourier interpolating the real-space force constants [4] [13].
Acoustic Sum Rule (ASR)	A constraint applied in q2r.x (zasr) or matdyn.x (asr) to enforce zero frequency for acoustic modes at Γ, fixing small errors from finite grid effects [4] [5] [13].
Eigenvector of Imaginary Mode	The atomic displacement pattern associated with a soft mode; used to "follow" the mode to a lower-energy structure [38].
Path-Integral Molecular Dynamics (PIMD)	A simulation method that incorporates nuclear quantum effects, essential for correctly modeling phase transitions in systems with light atoms like hydrogen [39].

Conclusion

Successfully resolving negative phonon frequencies in Quantum ESPRESSO requires a multifaceted approach that combines rigorous computational methodology with careful physical interpretation. By systematically addressing convergence parameters, structural relaxation quality, and appropriate application of sum rules, researchers can distinguish numerical artifacts from genuine physical instabilities. The implementation of robust validation protocols ensures reliable phonon spectra that accurately represent material properties. Mastering these techniques enables confident prediction of dynamical stability and lattice vibrational properties, providing crucial insights for materials design and drug development applications where understanding molecular interactions and stability is paramount. Future directions include leveraging improved exchange-correlation functionals, advanced solvation models for biological systems, and machine learning approaches to accelerate convergence in complex molecular simulations.

References

• 1. [Pw_forum] Physical meaning of two imaginary phonon ... [https://lists.quantum-espresso.org/pipermail/users/2017-November/039744.html]
• 2. Quantum Espresso - Obtaining Negative Phonon ... [https://mattermodeling.stackexchange.com/questions/13031/quantum-espresso-obtaining-negative-phonon-frequencies]
• 3. Re: [QE-users] Phonon: negative frequencies for adsorbate ... [https://www.mail-archive.com/users@lists.quantum-espresso.org/msg39359.html]
• 4. How to get rid of this negative frequencies? [https://mattermodeling.stackexchange.com/questions/13277/how-to-get-rid-of-this-negative-frequencies]
• 5. 6 Troubleshooting [https://www.quantum-espresso.org/Doc/ph_user_guide/node18.html]
• 6. [QE-users] Negative frequencies for phonon dispersion in ... [https://lists.quantum-espresso.org/pipermail/users/2023-August/050750.html]
• 7. Structure optimization • Quantum Espresso Tutorial [https://pranabdas.github.io/espresso/hands-on/structure-optimization/]
• 8. Quantum Espresso [https://www.materialssquare.com/docs/en/appendix]
• 9. First-principles investigation of the structural and dynamical ... [https://www.sciencedirect.com/science/article/abs/pii/S2452262721001070]
• 10. ph.x: input description [https://www.quantum-espresso.org/Doc/INPUT_PH.html]
• 11. phonons [https://www.quantum-espresso.org/faq/faq-phonons/]
• 12. quantum espresso - Negative phonon calculation [https://mattermodeling.stackexchange.com/questions/13046/negative-phonon-calculation]
• 13. Phonon dispersion • Quantum Espresso Tutorial - of Pranab Das [https://pranabdas.github.io/espresso/hands-on/phonon]
• 14. matdyn.x: input description [https://www.quantum-espresso.org/Doc/INPUT_MATDYN.html]
• 15. quantum espresso - Phonon calculations related problems [https://mattermodeling.stackexchange.com/questions/12940/phonon-calculations-related-problems]
• 16. Self consistent field calculation for silicon • Quantum Espresso ... [https://pranabdas.github.io/espresso/hands-on/scf/]
• 17. 5 Troubleshooting [https://www.quantum-espresso.org/Doc/pw_user_guide/node21.html]
• 18. Negative frequencies in linewidth claculation [https://forum.epw-code.org/viewtopic.php?t=233]
• 19. 1. Running an SCF job and calculating the band structure ... [https://www.cp2k.org/exercises:2014_ethz_mmm:bs]
• 20. DFPT FAQs [https://docs.abinit.org/faq/dfpt/]
• 21. Phonon Dispersions and Density of States Calculation [https://docs.mat3ra.com/tutorials/dft/vibrational/phonon-dispersion-dos/]
• 22. Phonon Band Structure, DOS, and Thermodynamics - SCM [https://www.scm.com/doc/Tutorials/ExternalPrograms/QEPhonons.html]
• 23. High-Tc 2D ambient BCS superconductors in ... [https://www.nature.com/articles/s41699-025-00590-0]
• 24. Topological and superconducting properties of two- ... [https://pubs.rsc.org/en/content/articlehtml/2025/cp/d5cp01206f]
• 25. pw.x: input description [https://www.quantum-espresso.org/Doc/INPUT_PW.html]
• 26. Hands-on • Quantum Espresso Tutorial - of Pranab Das [https://pranabdas.github.io/espresso/category/hands-on/]
• 27. Setting up automatic SSCHA simulations on clusters [http://sscha.eu/Tutorials/tutorial_02_advanced_submission/]
• 28. q2r.x: input description [https://www.quantum-espresso.org/Doc/INPUT_Q2R.html]
• 29. Defects in 2D materials [http://jdftx.org/Defects2D.html]
• 30. dynmat.x: input description [https://www.quantum-espresso.org/Doc/INPUT_DYNMAT.html]
• 31. [Pw_forum] question of matdyn.x and dynmat.x [https://lists.quantum-espresso.org/pipermail/users/2013-July/027668.html]
• 32. matdyn.x: input description [https://www.densityflow.com/about/INPUT/INPUT_MATDYN.html]
• 33. Error in calculating FTIR and RAMAM using Quantum ... [https://mattermodeling.stackexchange.com/questions/12689/error-in-calculating-ftir-and-ramam-using-quantum-espresso]
• 34. Tutorial: Phonons, EELS and magnons for HPC and GPUs [https://enccs.github.io/efficient-materials-modelling-on-hpc/espresso-tutorial-phonons/]
• 35. GGA Vs LDA | PDF | Density Functional Theory | Phonon [https://www.scribd.com/document/851655398/GGA-vs-LDA]
• 36. Quantum ESPRESSO: Magnetism, Band Structure and pDOS [https://www.scm.com/doc/Tutorials/ExternalPrograms/QEMagnetismBandsPDOS.html]
• 37. First-principles calculations to investigate structural ... [https://www.sciencedirect.com/science/article/abs/pii/S2352214322000405]
• 38. How to deduce phase transitions from a phonon calculation? [https://mattermodeling.stackexchange.com/questions/1635/how-to-deduce-phase-transitions-from-a-phonon-calculation]
• 39. Quantum simulation of thermally-driven phase transition ... [https://www.nature.com/articles/srep03272]