Phonon Calculations with DFPT in Quantum ESPRES…

Phonons — the quantized vibrations of a crystal lattice — govern thermal conductivity, heat capacity, thermal expansion, superconductivity, and dynamical stability. Density functional perturbation theory (DFPT) computes them directly from the ground-state electronic structure without the large supercells that finite-displacement methods require. This tutorial walks through the complete Quantum ESPRESSO phonon workflow using silicon.

Why DFPT Instead of Finite Displacements

The finite-displacement (frozen-phonon) approach displaces atoms one at a time in a large supercell and measures the restoring forces. It works, but the supercell size scales with the phonon wavelength you want to resolve, becoming expensive fast.

DFPT instead treats an atomic displacement as a small perturbation and solves a self-consistent linear-response equation for the change in the wavefunctions and density. Phonon frequencies $\omega_{\nu\mathbf{q}}$ are eigenvalues of the dynamical matrix

$$ D_{\alpha\beta}^{IJ}(\mathbf{q}) = \frac{1}{\sqrt{M_{I}M_{J}}}\sum_{\mathbf{R}} \frac{\partial^{2}E}{\partial u_{I\alpha}^{0}\,\partial u_{J\beta}^{\mathbf{R}}}\,e^{i\mathbf{q}\cdot\mathbf{R}} $$

Its decisive advantage: it computes the response at any $\mathbf{q}$-vector using only the primitive cell, because a phonon of wavevector $\mathbf{q}$ couples to electronic states shifted by $\mathbf{q}$ without needing a real-space supercell. This makes full phonon dispersions across the Brillouin zone tractable.

The Phonon Workflow at a Glance

The pipeline has four executables:

flowchart LR
  A["pw.x scf<br/>tight conv_thr"] --> B["ph.x<br/>q-grid DFPT"]
  B --> C[Dynamical matrices]
  C --> D["q2r.x<br/>IFCs + ASR"]
  D --> E["matdyn.x path<br/>dispersion"]
  D --> F["matdyn.x dos<br/>DOS + thermo"]

pw.x — a tightly converged SCF calculation on the relaxed structure.
ph.x — DFPT response on a uniform $\mathbf{q}$-point grid; produces dynamical matrices.
q2r.x — Fourier-transforms the dynamical matrices into real-space interatomic force constants (IFCs).
matdyn.x — interpolates the IFCs to any $\mathbf{q}$-point for dispersions, DOS, and thermodynamics.

A cardinal rule precedes all of it: relax the geometry first. Phonon frequencies are second derivatives of the energy, so residual forces on an unrelaxed structure produce spurious imaginary modes.

Step 1: A Well-Converged SCF Calculation

Start from a relaxed silicon primitive cell. Phonons demand tighter convergence than an ordinary SCF because you are differentiating the energy twice.

&CONTROL
    calculation = 'scf'
    prefix      = 'si'
    outdir      = './out'
    pseudo_dir  = './pseudo'
    tprnfor     = .true.
    tstress     = .true.
/
&SYSTEM
    ibrav     = 2
    celldm(1) = 10.20
    nat       = 2
    ntyp      = 1
    ecutwfc   = 45.0
    ecutrho   = 360.0
/
&ELECTRONS
    conv_thr    = 1.0d-12
    mixing_beta = 0.7
/
ATOMIC_SPECIES
 Si  28.0855  Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS (alat)
 Si  0.00  0.00  0.00
 Si  0.25  0.25  0.25
K_POINTS (automatic)
 8 8 8 0 0 0

Note conv_thr = 1.0d-12 — several orders tighter than a routine run. The linear-response step inherits any error in the ground state.

pw.x -in si.scf.in > si.scf.out

Step 2: Running ph.x on a q-Point Grid

The ph.x code computes dynamical matrices. Requesting a uniform grid with ldisp = .true. and nq1 nq2 nq3 produces the matrices at every inequivalent q needed for a full dispersion.

Phonons on a 4x4x4 q-grid for silicon
&INPUTPH
    prefix   = 'si'
    outdir   = './out'
    fildyn   = 'si.dyn'
    ldisp    = .true.
    nq1      = 4
    nq2      = 4
    nq3      = 4
    tr2_ph   = 1.0d-14
    fildvscf = 'dvscf'
/

ph.x -in si.ph.in > si.ph.out

Key parameters:

tr2_ph: the self-consistency threshold for the linear-response equation; 1.0d-14 to 1.0d-16 is typical.
ldisp = .true. with nq1/nq2/nq3: generates the regular q-grid; without it, you specify a single q-point manually.
fildyn: prefix for the dynamical-matrix output files, one per irreducible q.
epsil = .true.: add this for polar/insulating materials to compute the dielectric tensor and Born effective charges, needed for the LO-TO splitting correction.

This is by far the most expensive stage. It can be parallelized over q-points and irreducible representations (start_q/last_q, start_irr/last_irr) and run as independent jobs, then collected.

Step 3: Real-Space Force Constants with q2r.x

q2r.x inverse-Fourier-transforms the set of dynamical matrices on the coarse q-grid into short-ranged interatomic force constants in real space.

&INPUT
    fildyn = 'si.dyn'
    zasr   = 'simple'
    flfrc  = 'si.fc'
/

q2r.x -in si.q2r.in > si.q2r.out

The zasr keyword enforces the acoustic sum rule, which guarantees the three acoustic branches go exactly to zero frequency at Gamma — a consequence of translational invariance that numerical noise otherwise violates. The output si.fc file contains the IFCs used for interpolation.

Step 4: Dispersions and DOS with matdyn.x

With the IFCs in hand, matdyn.x cheaply evaluates phonon frequencies at any q. First, the dispersion along a high-symmetry path:

&INPUT
    flfrc  = 'si.fc'
    asr    = 'simple'
    flfrq  = 'si.freq'
    q_in_band_form = .true.
/
5
  0.500 0.500 0.500 20   ! L
  0.000 0.000 0.000 20   ! Gamma
  0.500 0.000 0.500 20   ! X
  0.750 0.250 0.750 20   ! W
  0.000 0.000 0.000 1    ! Gamma

matdyn.x -in si.matdyn.in > si.matdyn.out

For the phonon DOS, switch to a dense uniform grid instead of a path:

&INPUT
    flfrc = 'si.fc'
    asr   = 'simple'
    dos   = .true.
    nk1   = 24
    nk2   = 24
    nk3   = 24
    fldos = 'si.phdos'
    deltaE = 1.0
/

The si.freq.gp file (dispersion) and si.phdos (DOS) are ready for gnuplot. Silicon shows six branches — three acoustic starting from zero at Gamma and three optical near 500-520 cm^-1 — with the characteristic degeneracy at the zone center.

Step 5: Diagnosing Imaginary Frequencies

matdyn.x prints frequencies in cm^-1. A negative (imaginary) frequency signals that the structure is a saddle point or maximum along that vibrational mode, not a true minimum — the lattice is dynamically unstable and will spontaneously distort. This is a physically meaningful result: it is how DFPT predicts ferroelectric, charge-density-wave, and structural phase transitions.

Observation	Interpretation
All frequencies real, positive	Dynamically stable structure
Small imaginary near Gamma (< ~10 cm^-1)	Usually numerical: tighten `tr2_ph`, `conv_thr`, k-grid, ASR
Large imaginary at a specific q	Real instability; a supercell distortion at that q lowers the energy
Imaginary acoustic modes at Gamma	ASR not enforced or unrelaxed structure

When you see genuine imaginary modes, follow the eigenvector: displace the atoms along that mode, relax, and re-run to find the true ground-state phase.

Step 6: Thermal Properties

Within the harmonic approximation, the phonon DOS yields temperature-dependent thermodynamic functions — the Helmholtz free energy, internal energy, entropy, and constant-volume heat capacity. The vibrational free energy is

$$ F_{\mathrm{vib}}(T) = \sum_{\nu\mathbf{q}}\Biggl[\frac{\hbar\omega_{\nu\mathbf{q}}}{2}

k_{\mathrm{B}}T\ln\bigl(1 - e^{-\hbar\omega_{\nu\mathbf{q}}/k_{\mathrm{B}}T}\bigr)\Biggr] $$

The matdyn.x DOS run can compute these directly, or you can post-process the DOS. The heat capacity follows the Debye $T^{3}$ law at low temperature and approaches the Dulong–Petit limit of $3Nk_{\mathrm{B}}$ at high temperature. Combining phonon free energies at several volumes gives the quasi-harmonic approximation for thermal expansion and the temperature dependence of the bulk modulus.

Quantum ESPRESSO ships a dedicated utility for this. Run matdyn.x with dos = .true. to produce the phonon DOS, then feed it forward: the harmonic thermodynamic quantities are tabulated as a function of temperature, letting you plot heat capacity and entropy curves directly against calorimetry data.

Property	Low-T limit	High-T limit
Heat capacity C_v	proportional to T^3 (Debye)	3NkB (Dulong-Petit)
Entropy S	approaches 0 (third law)	grows logarithmically
Free energy F	E_0 + zero-point energy	dominated by -TS term

The zero-point energy $E_{\mathrm{ZPE}} = \sum_{\nu\mathbf{q}}\frac{1}{2}\hbar\omega_{\nu\mathbf{q}}$ is a purely quantum contribution that shifts total energies and is essential for accurate reaction energetics involving light elements like hydrogen.

Parallelization and Managing Cost

DFPT is the most compute-intensive standard workflow in Quantum ESPRESSO, so understanding how to split it up is essential for anything beyond a toy cell. The ph.x step is embarrassingly parallel across two independent axes: q-points and irreducible representations of the dynamical matrix at each q. You can partition a large calculation into many small jobs and merge the results.

Phonons for a single q-block, images 1-2
&INPUTPH
    prefix   = 'si'
    outdir   = './out'
    fildyn   = 'si.dyn'
    ldisp    = .true.
    nq1      = 4
    nq2      = 4
    nq3      = 4
    tr2_ph   = 1.0d-14
    start_q  = 1
    last_q   = 3
/

Running separate jobs with different start_q/last_q ranges — or start_irr/last_irr within a q-point — lets you spread the workload across nodes and reconstruct the full set of dynamical matrices afterward. Combined with the plane-wave and pool parallelism of the underlying SCF, this hierarchical decomposition is what makes full-Brillouin-zone phonons feasible for realistic materials.

The recover mechanism (recover = .true.) resumes an interrupted ph.x run from its last completed representation, which is critical because a single q-point can take many hours. Always enable it for production jobs so a walltime kill does not discard partial work.

Common Pitfalls and Best Practices

Unrelaxed geometry: the number-one cause of spurious imaginary modes. Always relax with tight forc_conv_thr first.
Loose SCF or tr2_ph: linear response amplifies ground-state error; keep conv_thr near 1.0d-12 and tr2_ph near 1.0d-14.
Forgetting the ASR: without it the acoustic branches do not vanish at Gamma.
Ignoring LO-TO splitting: for polar insulators, omitting epsil = .true. misses the optical-mode splitting at Gamma.
q-grid too coarse: a sparse nq grid gives poorly interpolated dispersions; converge it as you would a k-grid.

For related electronic-structure guides see the blog and the platform’s computational engines; the definitive reference is the official Quantum ESPRESSO documentation.

Run it on Simatra

DFPT phonon runs are among the heaviest workloads in DFT — ph.x sweeps every irreducible q-point and representation, and the cost climbs steeply with cell size. Simatra runs Quantum ESPRESSO on GPU-accelerated clusters (instance GPU-Opt-V2) with up to 5x faster convergence and support for supercells up to ~2,000 atoms, so full phonon dispersions finish in a fraction of the time. Get started with a free trial and $100 in credits at app.simatra.io, and explore the computational engines — Quantum ESPRESSO alongside our native C++20 KRONOS engine with CUDA/HIP/Metal backends — that power the platform.