Phonons by finite displacement

The supercell / frozen-phonon (direct) route to the harmonic force constants $\Phi$, and the one I actually run for systems where the code has no analytic perturbation theory: displace a small set of symmetry-inequivalent atoms by a finite $\delta$ in a commensurate supercell, read off the Hellmann–Feynman forces on every atom, and finite-difference them into $\Phi$. From $\Phi$ the dynamical matrix $D(\mathbf{q})$ and the phonon spectrum follow by the standard mass-weighted Fourier transform. It is the method-of-record in my CRYSTAL23 + phonopy setup, the complementary route to DFPT linear response, and the only first-principles option when analytic response is unavailable in the code.

Force constants from a finite displacement

The harmonic energy is the second-order Taylor expansion of the Born–Oppenheimer PES in the atomic displacements $u_{l\kappa \alpha }$ ($l$ = lattice vector, $\kappa$ = basis atom, $\alpha$ = Cartesian direction) about equilibrium,

$$E=E_0+\frac{1}{2}\sum _{l\kappa \alpha }\sum _{l^{\prime }{\kappa }^{\prime }\beta }{\Phi }_{\alpha \beta }^{l\kappa ,l^{\prime }{\kappa }^{\prime }}{u}_{l\kappa \alpha }{u}_{l^{\prime }{\kappa }^{\prime }\beta }+\mathcal{O}(u^3),{\Phi }_{\alpha \beta }^{l\kappa ,l^{\prime }{\kappa }^{\prime }}=\frac{{\partial }^{2}E}{\partial u_{l\kappa \alpha }\partial u_{l^{\prime }{\kappa }^{\prime }\beta }},$$

the linear term vanishing because the forces are zero at a relaxed geometry. The whole idea of the direct method is to never differentiate the energy: the force constant is minus a first derivative of the force,

$${\Phi }_{\alpha \beta }^{l\kappa ,l^{\prime }{\kappa }^{\prime }}=-\frac{\partial F_{l^{\prime }{\kappa }^{\prime }\beta }}{\partial u_{l\kappa \alpha }}\approx -\frac{F_{l^{\prime }{\kappa }^{\prime }\beta }[+\delta {\hat{e}}_{l\kappa \alpha }]-F_{l^{\prime }{\kappa }^{\prime }\beta }[-\delta {\hat{e}}_{l\kappa \alpha }]}{2\delta },$$

so a single displacement of atom $(l\kappa )$ along $\alpha$ by $\delta$ yields the entire column ${\Phi }_{\alpha \beta }^{l\kappa ,l^{\prime }{\kappa }^{\prime }}$ for all $(l^{\prime }{\kappa }^{\prime }\beta )$ at once — every other atom’s force in the displaced cell is one finite-difference data point. This is the leverage of using forces rather than energies: the energy gives one scalar per calculation; the force field gives $3N$ numbers, the full set of partners of the displaced atom. Because the Hellmann–Feynman forces are already a by-product of the SCF (an expectation value, no extra electronic response), each displaced supercell is just one more single-point calculation.

The displacement amplitude $\delta$ is a bias–variance compromise. Too large and the cubic and quartic terms $\mathcal{O}(u^3)$ contaminate the second derivative; too small and the force difference drowns in the SCF force noise. The two-sided (central) difference above is the standard choice precisely because it cancels the leading cubic contamination, so the residual systematic error is $\mathcal{O}({\delta }^2)$ rather than $\mathcal{O}(\delta )$; phonopy’s default $\delta =0.01$ Å sits comfortably in the harmonic window for typical covalent and ionic solids. The hard requirement is the one I flag in imaginary modes: the SCF forces must be converged tightly enough that the finite-difference noise floor on $\Phi$ sits below the acoustic-mode energy scale, or spurious near-$\Gamma$ frequencies appear that are numerical, not physical. And with the atom-centred Gaussian basis of CRYSTAL the forces must include the Pulay term (see Hellmann–Feynman forces) — omitting it corrupts $\Phi$ directly and fabricates imaginary modes.

Symmetry reduction of the displacement set

The naive cost is $3N_{\text{sc}}$ supercell calculations (each atom, each Cartesian direction, times two for the central difference) where $N_{\text{sc}}$ is the number of atoms in the supercell — prohibitive. The space group of the crystal collapses this dramatically, and this is where crystal symmetry does the work:

• Site-symmetry equivalence. Only one atom from each Wyckoff orbit (each set of symmetry-equivalent atoms) need be displaced; the force constants of the others are obtained by applying the space-group operations that relate them, ${\Phi }^{\prime }=R\Phi R^{\top }$ for the point part $R$ of the relating operation.
• Site point group. For a chosen representative atom, its own site-symmetry group fixes which Cartesian displacement directions are independent. An atom on a high-symmetry site (a mirror, a rotation axis) needs fewer than three displacements; the remaining components of $\Phi$ are generated by the site operations.
• Permutation symmetry ${\Phi }_{\alpha \beta }^{l\kappa ,l^{\prime }{\kappa }^{\prime }}={\Phi }_{\beta \alpha }^{l^{\prime }{\kappa }^{\prime },l\kappa }$ and the translational/rotational invariance conditions further constrain the tensor.

phonopy enumerates this reduced, symmetry-inequivalent displacement set automatically from the input symmetry, writes one displaced supercell per required displacement, and after the forces come back symmetrizes the assembled force constants to restore the full space-group and permutation symmetry exactly (its --sym-fc / FC_SYMMETRY step). For a high-symmetry cell this can be the difference between a handful of displacements and many dozens. The reduction is exact in principle; in practice the symmetrization also launders a little of the finite-difference noise, because it projects $\Phi$ onto the symmetry-allowed subspace.

The supercell-size vs. interaction-range tradeoff

This is the central approximation of the direct method, and the one the user must own per system. $D(\mathbf{q})$ is the Fourier transform of $\Phi$ over lattice vectors,

$$D_{\alpha \beta }^{\kappa {\kappa }^{\prime }}(\mathbf{q})=\frac{1}{\sqrt{M_{\kappa }{M}_{{\kappa }^{\prime }}}}\sum _{l^{\prime }}{\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}{e}^{i\mathbf{q}\cdot ({\mathbf{R}}_{l^{\prime }}-{\mathbf{R}}_0)},$$

and a finite Born–von Kármán supercell can only sample $\Phi$ out to the partners it physically contains: interactions beyond the supercell boundary are silently truncated (or, equivalently, folded back under the periodic boundary conditions). The method is therefore exact only if $\Phi$ has decayed to within the SCF noise before reaching the supercell edge. Two regimes:

• Short-ranged $\Phi$ (covalent and most non-polar crystals): the force constants fall off within a few bond lengths, so a $2\times 2\times 2$–$4\times 4\times 4$ supercell already converges the dispersion. A practical convergence test is to enlarge the supercell and watch the zone-boundary and zone-interior frequencies stop moving. Note the asymmetry at $\Gamma$ itself: the three acoustic $\Gamma$ frequencies are pinned to zero by the ASR regardless of supercell size, so they are not a useful diagnostic of range — but the optical $\Gamma$ frequencies are genuine Fourier sums of $\Phi$ over the supercell and do shift with the interaction range, so a too-small supercell silently corrupts exactly those optical $\Gamma$ modes (along with the zone-interior dispersion) while the acoustic modes still look perfect. The range diagnostic must therefore watch the optical and zone-boundary frequencies, not the acoustic-$\Gamma$ values.
• Long-ranged $\Phi$ (polar crystals): the dipole–dipole interaction mediated by Born effective charges decays only as $1/r^3$ and is not captured by any affordable supercell. This is handled separately — the long-range dipole part is subtracted analytically and reintroduced as the non-analytic LO–TO correction at $\mathbf{q}\to 0$ (Born charges + ${\epsilon }^{\infty }$, the inputs supplied by DFPT or the Berry-phase route), so the supercell only has to converge the short-range remainder. This is exactly the NAC machinery the dynamical matrix & ASR note invokes.

The tradeoff is brutal because cost and memory scale with the supercell. The number of electrons (hence the SCF cost of each displaced single-point) grows linearly with $N_{\text{sc}}$, and the supercell volume — and with it the number of $k$-points needed for a fixed reciprocal-space density — counterbalances: a supercell $n$ times larger in each direction needs a $k$-mesh $n$ times coarser in each direction to hold the sampling fixed, which partly offsets the per-cell cost but does nothing for the memory of holding the larger system and its integrals. The art is choosing the smallest supercell whose $\Phi$ has genuinely decayed inside it. For the layered van der Waals systems I work on — the WSe₂ heterostructure campaign (CRYSTAL23 + phonopy on the MSU HPCC) — the in-plane range is short but the cells are atom-heavy, so the supercell/memory bookkeeping is the binding constraint, not the interaction range.

ASR / sum-rule restoration

The continuous space-group symmetry includes uniform translation: a rigid shift of the entire crystal costs no energy, which forces the acoustic sum rule (ASR) on the force constants,

$$\sum _{l^{\prime }{\kappa }^{\prime }}{\Phi }_{\alpha \beta }^{0\kappa ,l^{\prime }{\kappa }^{\prime }}=0⟹{\omega }_{\Gamma }^{\text{acoustic}}=0.$$

Finite-difference forces never satisfy this exactly — incomplete-basis and grid errors leave a residual drift — so the acoustic branches do not go cleanly to zero at $\Gamma$; the deficit shows up as small spurious frequencies (a few cm${}^{-1}$ when the forces are good, alarmingly large when the basis or mesh is poor). It is restored numerically by enforcing the sum rule on the assembled $\Phi$, phonopy’s FC_SYMMETRY / acoustic-sum-rule step (the simplest version subtracts the row sum from the on-site block ${\Phi }^{0\kappa ,0\kappa }$ to drive ${\sum }_{l^{\prime }{\kappa }^{\prime }}\Phi =0$). I treat ASR restoration as mandatory before reading any frequency near $\Gamma$. The theory — why translational invariance implies the sum rule, the LO–TO non-analytic term for polar systems, and the further rotational and Huang (vanishing-stress) invariance conditions that FC_SYMMETRY does not impose (and which matter for the flexural ZA branch of 2D systems) — is in the dynamical matrix & ASR and pursued in imaginary modes. The rule of thumb: passing FC_SYMMETRY rules out the translational artifact and nothing more.

The phonopy + CRYSTAL23 workflow

Concretely, the pipeline I run is four stages, with phonopy orchestrating and CRYSTAL23 as the force engine. I describe the interface only at the documented level; CRYSTAL-side keyword and filename specifics I am not certain of are left as explicit placeholders rather than guessed.

1. Displacement generation. Relax the primitive cell tightly first (residual forces and all stress components small — the vacuum-direction stress matters for 2D, see imaginary modes). Hand phonopy the unit cell and supercell matrix; it determines the space group, builds the symmetry-reduced displacement set, and writes one displaced supercell per required displacement. It records the displacement metadata in phonopy_disp.yaml and emits a perfect (undisplaced) supercell plus the displaced supercells. With the CRYSTAL interface phonopy writes the structures in CRYSTAL’s input geometry format rather than VASP POSCAR/SPOSCAR files. TODO(user): confirm the exact CRYSTAL23 geometry-input filenames/format phonopy writes for the displaced supercells, and the CRYSTAL keyword block used for the single-point force evaluation (gradient-on), against the phonopy CRYSTAL-interface docs.
2. Force harvesting. Run a single-point SCF with analytic gradients on each displaced supercell in CRYSTAL23 (Pulay-corrected forces, tight SCF convergence). These are independent calculations — the embarrassingly parallel part of the campaign, one job per displacement. TODO(user): verify the exact CRYSTAL23 output-file name(s) phonopy parses the forces from, and any required print/keyword to make CRYSTAL emit the gradient in the parsed format.
3. Assemble FORCE_SETS. phonopy reads the forces back from each CRYSTAL output and collates them, together with the displacement pattern, into its FORCE_SETS file (the displacement-and-forces dataset). From FORCE_SETS it builds the supercell force constants $\Phi$ (FORCE_CONSTANTS / force_constants.hdf5), applying the symmetrization and the ASR (FC_SYMMETRY) at this step.
4. $D(\mathbf{q})$ and post-processing. With $\Phi$ in hand, phonopy assembles $D(\mathbf{q})$ at any $\mathbf{q}$ by the Fourier sum above (adding the non-analytic LO–TO correction for polar systems, given Born charges and ${\epsilon }^{\infty }$ in BORN), diagonalizes for ${\omega }_{\mathbf{q}\nu }$, and produces band structures along a $\mathbf{q}$-path, the phonon DOS, and the harmonic thermodynamic functions (free energy $F(T)$, entropy, heat capacity $C_v$) by Brillouin-zone integration.

Mesh / memory tradeoffs on the DOS and free energy. The dispersion along a path is cheap — it is a handful of $D(\mathbf{q})$ diagonalizations. The expense in post-processing is the $\mathbf{q}$-mesh used for the DOS and the thermodynamic integrals: the free energy $F(T)=\frac{1}{2}{\sum }_{\mathbf{q}\nu }\hslash {\omega }_{\mathbf{q}\nu }+k_{B}T{\sum }_{\mathbf{q}\nu }\ln (1-e^{-\hslash {\omega }_{\mathbf{q}\nu }/k_{B}T})$ is a sum over the mesh, and because the Fourier interpolation of $D(\mathbf{q})$ from a fixed $\Phi$ is analytic, this mesh is essentially free of SCF cost — the tradeoff is post-processing memory and time (diagonalizing $D$ on every mesh point and storing the eigendata), not electronic-structure cost. The practical discipline is a convergence sweep on the mesh: coarse meshes under-sample the low-frequency acoustic region and bias $F(T)$ and $C_v$ at low $T$ (where the soft modes dominate the partition function), so I converge the mesh on the integrated quantity I actually report, not just on the visual DOS. A tetrahedron-method integration converges faster than smearing for the same mesh. This is also the seam where an MLIP earns its keep: once a potential reproduces the DFT forces, the supercell can be enlarged and the $\mathbf{q}$-mesh refined almost without cost, lifting exactly the range/mesh constraints that bind the DFT calculation here.

Contrast with DFPT — the linear-response route

DFPT reaches the same force constants by analytic linear response instead of finite differences. The structural difference is decisive: a phonon of wavevector $\mathbf{q}$ is a perturbation $\propto e^{i\mathbf{q}\cdot \mathbf{R}}$ whose self-consistent response couples a Bloch state only to its $\mathbf{k}+\mathbf{q}$ partner — a monochromatic response — so DFPT delivers $D(\mathbf{q})$ at arbitrary $\mathbf{q}$ from a single primitive-cell calculation, with no supercell at all. Finite displacement, by contrast, ties the accessible wavevectors to the supercell: a zone-boundary mode at $\mathbf{q}=(\frac{1}{2},0,0)$ needs a doubled cell, a general commensurate $\mathbf{q}$ needs the matching supercell, and a small incommensurate $\mathbf{q}$ is simply inaccessible (one Fourier-interpolates from the commensurate grid the supercell does provide). The methods agree where they overlap — both are the harmonic $\Phi$ — and the choice between them is largely about what the code offers and what the system costs.

The reason this note exists as a standalone method, and not merely a fallback, is that DFPT is not always available. Analytic perturbation theory has to be implemented for each Hamiltonian, basis, and perturbation; many production codes — CRYSTAL among them for general phonon perturbations — do not expose a full phonon-DFPT, whereas any code that can compute forces can be driven by phonopy in finite displacement. So the practical division is: DFPT when the code has it and arbitrary-$\mathbf{q}$ / incommensurate response or cheap third-order ($2n+1$) anharmonicity is needed; finite displacement when the code lacks analytic response, when an interatomic potential or MLIP is the force engine, or when the short-range $\Phi$ of a non-polar solid converges in a modest supercell. For me the deciding factor is simply that finite displacement is what my CRYSTAL23 + phonopy stack supports.

Prerequisites

dynamical matrix & ASR · Hellmann–Feynman forces · crystal symmetry

Builds toward: imaginary modes · WSe₂ · DFPT linear response · MLIP phonons

Key references

• The direct / supercell method — K. Parlinski, Z. Q. Li & Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997); G. Kresse, J. Furthmüller & J. Hafner, Europhys. Lett. 32, 729 (1995).
• phonopy — A. Togo & I. Tanaka, Scr. Mater. 108, 1 (2015); A. Togo, J. Phys. Soc. Jpn. 92, 012001 (2023).
• Phonons in CRYSTAL (atom-centred / direct method) — F. Pascale, C. M. Zicovich-Wilson, F. López Gejo, B. Civalleri, R. Orlando & R. Dovesi, J. Comput. Chem. 25, 888 (2004). TODO(cite): CRYSTAL23 release paper — Dovesi et al., add exact journal/volume/page once confirmed.