Crystal Phonons and Lattice Dynamics
On this page
- Harmonic energy of a periodic solid
- Bloch form and the dynamical matrix
- Acoustic and optical branches
- Finite-displacement direct method
- Supercell size and real-space truncation
- Displacement amplitude
- Primitive-cell standardization
- Reciprocal-space paths
- What the current DOS means
- Polar crystals and the non-analytic limit
- Constraints, isotope masses, and primitive reconstruction
- Reading stability evidence
- Execution sequence and partial results
- Cost and memory
- Convergence design
- Limitations and reporting
- References
A molecular vibration describes motion around one isolated equilibrium geometry. A crystal phonon describes a collective displacement pattern that repeats through an infinite lattice with a wave vector . Periodicity turns one finite Hessian into a real-space family of intercell force constants, and a Fourier transform of that family produces a different dynamical matrix at every point in the Brillouin zone. Stability at alone therefore does not establish dynamical stability throughout the crystal.
The visible command path belongs to Calculate Crystal Phonons. Script construction and result handling belong to the separate Crystal Phonons with Tako Script. This chapter defines the physical and numerical objects those procedures manipulate.
Harmonic energy of a periodic solid
Let be the lattice vector of cell , let locate basis atom inside a reference cell, and let be its displacement. Expanding the Born–Oppenheimer energy around an equilibrium crystal gives
where are Cartesian directions and
is a harmonic interatomic force constant. Translational symmetry makes the coupling depend only on the cell difference. At an exact equilibrium, the first derivatives vanish. In practice, residual forces and stress matter because the force constants are evaluated about the supplied geometry, not about an automatically repaired structure.
The harmonic approximation discards cubic and higher derivatives. It describes infinitesimal motion on one potential-energy surface. Thermal expansion, finite-temperature phonon renormalization, phonon–phonon lifetimes, strong anharmonic stabilization, diffusion, disorder, and phase hopping require additional theory. A harmonic imaginary mode can be valuable evidence of a zero-temperature instability, but it is not by itself a prediction that a finite-temperature phase cannot exist.
Bloch form and the dynamical matrix
The classical equation of motion is
A normal mode of a periodic lattice has Bloch form
where is expressed in fractional reciprocal coordinates and labels the branch. Substitution produces the Hermitian eigenproblem
with
Tako uses this fractional-coordinate phase convention. Every compact force-constant block is associated with a wrapped integer lattice vector , and its real and imaginary phase contributions are accumulated before the Hermitian matrix is diagonalized. Positive eigenvalues give real frequencies. Negative eigenvalues give imaginary angular frequencies; Tako stores their signed square roots as negative frequencies rather than appending the symbol .
For active atoms in the primitive cell, has dimension and produces branches. The eigenvector phase is arbitrary, and at a degeneracy any orthonormal rotation within the degenerate subspace is equally valid. Compare a degenerate subspace or its projector, not individual eigenvector signs.
Acoustic and optical branches
Uniform translation of the complete crystal costs no energy. Consequently, three acoustic branches approach zero at . Their long-wavelength slopes determine sound velocities and relate to elastic response. A primitive cell containing more than one atom also has optical branches at , corresponding to relative motion of basis atoms.
The labels “acoustic” and “optical” describe branch character near , not a universal frequency ordering over the whole Brillouin zone. Branches can cross, avoid one another when symmetry allows coupling, exchange eigenvector character, or become degenerate along high-symmetry lines. A plot line connecting sorted frequencies does not by itself track modal identity. Assignment requires eigenvector overlap or symmetry information.
The acoustic condition is more than a visual expectation. Translational invariance requires the acoustic sum rule
for every source atom and pair of Cartesian components. Finite force noise, incomplete supercell convergence, and numerical differentiation violate this identity. Tako restores it during its default force-constant postprocessing, but that algebraic correction cannot turn an inaccurate potential or truncated interaction range into accurate physics.
Finite-displacement direct method
Tako obtains force constants from central differences of supercell forces. Displacing active atom in the reference primitive cell by along direction gives
Unlike molecular vibrationNfree, the public phonon operation exposes only this two-sided central stencil. For active atoms in the standardized primitive cell, the force-constant calculation requires displaced force evaluations. Each evaluation contains atoms only when every primitive atom is active; more generally it contains the complete primitive basis repeated by the supercell volume. Computational cost is therefore controlled separately by the number of displaced primitive atoms and the number of atoms evaluated in each repeated cell.
Tako’s library also evaluates an undisplaced supercell in cache-oriented lower-level paths, but the public browser operation currently has no restartable displacement cache. Cancellation or a failed force evaluation ends the operation; a new run repeats the required displaced calculations.
Frederiksen force correction
The raw calculator forces of a displaced periodic structure can have a nonzero total residual. Tako follows the ASE phonon default called the Frederiksen method: for each displaced structure it sums the force over all supercell atoms and subtracts that total from the displaced atom before differencing. This enforces a consistent net-force reference at the level of each displacement.
The correction is numerical conditioning, not a replacement for structural relaxation. If the undisplaced crystal has substantial forces or stress, the local harmonic expansion still describes the wrong reference state. Report the relaxation threshold and residual stress together with the phonon settings.
Symmetry and acoustic-sum-rule iterations
After forming compact image-resolved force constants, Tako performs three ASE-style postprocessing iterations. Each iteration first enforces the transpose relation between a lattice image and its opposite image,
by averaging the two blocks. It then adjusts the on-site reference-cell block so the acoustic sum rule holds. Repeating the two operations reduces the way one correction perturbs the other.
The public UI does not expose the iteration count, force method, a real-space cutoff, or symmetry-reduced displacement generation. The lower-level lattice-dynamics library supports those variants, but they are not settings of tako.phonon. Documentation and reports must describe what the public operation actually ran: central differences, Frederiksen correction, three symmetrization/ASR iterations, and no explicit real-space cutoff.
Supercell size and real-space truncation
A diagonal supercell [n1,n2,n3] samples force responses in a finite range of lattice images. Periodic repetition means that interactions beyond the effective half-supercell can fold onto nearer images. The direct-method force constants are converged only when increasing the supercell no longer changes the frequencies and branch shapes relevant to the question.
Supercell convergence is anisotropic. A layered material may need more repeats normal to or within the layers depending on the force range; a one-dimensional chain embedded in vacuum should not be enlarged blindly in all three directions. Compare physically nested cells and preserve the primitive cell, calculator, displacement, path, and relaxation protocol while changing the repeats.
The current application normalizes each repeat to a positive integer and caps it at 6. The default is [3,3,3] for nequix-phono and [2,2,2] for every other level. These are runtime defaults, not convergence claims. A 3×3×3 default can be less adequate than an anisotropic 4×4×1 cell for one material and unnecessarily expensive for another.
Larger band-path interpolation does not repair a small supercell. The path merely evaluates the Fourier interpolation of already-computed force constants at more points. If the real-space tensor is truncated or noisy, a smoother line is a smoother representation of the same error.
Displacement amplitude
The finite difference balances truncation and force noise. With displacement , the central derivative has leading truncation error of order , while cancellation and calculator noise are amplified approximately as . Too large a displacement samples anharmonic forces; too small a displacement subtracts nearly equal noisy forces.
Tako’s default is 0.01 Å. Establish numerical stability by repeating representative calculations at another physically reasonable amplitude while keeping the geometry, supercell, calculator, and path fixed. Examine both frequencies and eigenvector character. A higher-resolution band plot is not an amplitude test.
MLIP forces can be smooth even outside the training distribution, so numerical convergence of does not validate the learned potential. Conversely, an SCF-based force model can show displacement sensitivity because electronic convergence changes between displaced cells. Numerical and model validation are separate layers.
Primitive-cell standardization
For a periodic input with a valid cell, Tako constructs a SeekPath object and converts the input to its standardized primitive cell before making the displacement supercell. The result records phonon_cell: "primitive", and the console reports an atom-count change when standardization reduces the structure.
This behavior has several consequences:
- the number of branches is determined by the standardized primitive basis, not necessarily by the input conventional-cell atom count;
- supercell repeats apply to the standardized primitive lattice vectors;
- input constraint and custom-mass metadata are not transferred by the current primitive-cell reconstruction;
- comparisons with another code require the same primitive setting or an explicit mapping;
- a user-supplied reciprocal path must be expressed for the standardized primitive reciprocal basis.
Standardization does not remove residual forces, relax cell shape, select magnetic order, or repair occupancies. A slightly distorted structure can also produce a lower-symmetry primitive cell and a different recommended path than an idealized reference.
Reciprocal-space paths
A phonon band diagram samples selected lines between special points. It does not scan the complete Brillouin zone. Tako uses SeekPath to determine a crystallographic Bravais-lattice description, special points, and a recommended path when phononBandPath is empty. In the interface, Generate SeekPath writes the path string into the control so it can be inspected or edited before the run.
An explicit path such as G X W L G is parsed against the special points of the standardized primitive cell. G or a Gamma spelling denotes . phononBandPointsPerSegment is clamped to 1–64. The returned band object distinguishes path_source: "custom" from "seekpath", includes the fractional q-points, and—when path metadata is available—includes a linear plot coordinate, label positions, labels, special-point coordinates, and segment endpoints.
A path can miss an instability at a general wave vector. For a consequential stability claim, supplement the conventional path with an appropriate reciprocal-space mesh or independent calculation. Tako’s current top-level operation does not expose a q-mesh stability scan.
Plot coordinate, units, and group velocity
The fractional q-point coordinates are dimensionless coefficients of the primitive reciprocal basis. The plotted horizontal axis is a cumulative linear coordinate constructed along each selected segment; it is not itself a Cartesian reciprocal-vector component and should not be compared between unrelated paths without the segment metadata. Repeated labels can mark a discontinuity between nonconnected path segments. A line drawn across such a discontinuity would imply phonon information that was never evaluated, so the returned segments, label positions, and special points are part of the scientific result rather than decorative plot metadata.
The runtime returns band frequencies in cm⁻¹. Tako’s band viewer converts them to THz using 1 THz = 33.3564095 cm⁻¹ and labels the vertical axis Frequency (THz). Negative values remain negative after conversion and represent imaginary signed square roots. A screenshot of the viewer and the underlying JSON therefore use different numerical units; a report must name which representation supplied each value.
For a real branch, the group velocity is
where is the physical reciprocal-space wave vector rather than the dimensionless plot coordinate. Tako does not return group velocities. Estimating them from a screenshot slope is invalid because the horizontal coordinate changes direction and scale between path segments. A quantitative velocity calculation needs reciprocal-lattice vectors, consistent branch tracking, sufficiently dense q sampling, and a derivative with respect to physical .
The number of band points controls visual and numerical interpolation along the chosen lines, but not the force-constant fit. Too few points can hide a narrow dip between special points; more points can reveal it on that line. Neither choice adds off-path coverage. Preserve the exact path string and points per segment whenever plots are compared.
Branch continuity and degeneracy
At each q point, the runtime diagonalizes the dynamical matrix and reports eigenvalues in sorted order. Sorted index is not guaranteed to follow one physical branch through a crossing or near-degeneracy. Colorless line connection by array index can swap mode character even while the eigenvalues remain smooth. For assignments such as longitudinal/transverse character or atom participation, compare the complex eigenvector subspaces at adjacent q points and account for their arbitrary phase.
Crystal symmetry protects some degeneracies and constrains polarization directions. A small numerical splitting can reflect force-constant noise or a slight structural distortion, while a real splitting can signal lowered symmetry. Decide which by checking the standardized cell, symmetry tolerance used by SeekPath, relaxation, and supercell/displacement convergence; do not restore degeneracy by manually averaging plotted frequencies unless that operation is part of a documented symmetry model.
What the current DOS means
The public phonon_dos output is not a Brillouin-zone-integrated phonon density of states. It takes the positive -point harmonic energies from the standardized primitive cell and bins them into 80-scale histogram spacing, with a minimum bin width equivalent to 1 cm⁻¹. It returns bin-center energies in Hartree and cm⁻¹ plus integer-like weights.
A true crystal phonon DOS approximates
over a reciprocal-space mesh. The underlying lattice-dynamics library can perform q-point or Monkhorst–Pack DOS calculations, but the current public tako.phonon workflow does not call those paths. Therefore the artifact should be described as a histogram or quick fingerprint, not used for vibrational thermodynamics, heat capacity, entropy, or comparison with a measured neutron-weighted DOS.
Polar crystals and the non-analytic limit
Long-range dipole–dipole interactions cause longitudinal-optical/transverse-optical splitting near in polar crystals. The non-analytic correction depends on Born effective-charge tensors, the high-frequency dielectric tensor, cell volume, and the direction in which approaches zero.
The underlying lattice-dynamics library can attach Born charges and a dielectric tensor for non-analytic corrections. The current browser phonon operation does not request or return those data and builds its band without NAC. Consequently, the calculated limit of a polar material lacks explicit LO–TO splitting. This is a method boundary, not a plotting defect. State it when interpreting polar optical modes.
Constraints, isotope masses, and primitive reconstruction
The lower-level lattice-dynamics library can exclude completely fixed atoms and can use custom atomic masses. The current public periodic workflow does not preserve those metadata. SeekPath.primitive_atoms() reconstructs a new primitive Atoms object from atomic numbers, standardized positions, and the primitive lattice. It does not map input constraints or custom isotope masses onto the representatives. The public periodic calculation therefore displaces every reconstructed primitive atom and uses the element-default masses.
This is a consequential boundary. Applying a fixed-substrate mask before tako.phonon does not currently produce a periodic partial-Hessian adsorbate calculation, and isotope substitution on the input does not propagate through primitive standardization. Do not report either effect as modeled. Use a workflow that explicitly preserves the intended active selection or masses when those questions matter.
The nonperiodic compatibility branch does not call primitive_atoms() and retains its input Atoms object, so its constraint behavior follows molecular vibrations. That difference is another reason not to treat the periodic and fallback branches as interchangeable. A periodic result should report that all standardized primitive atoms were active and that default elemental masses were used.
Reading stability evidence
A negative signed frequency means the harmonic dynamical matrix has a negative eigenvalue at that sampled . The corresponding eigenvector indicates a direction in which the model energy falls to second order. Several cases must be separated:
- Structural nonstationarity. Residual atomic force or cell stress tilts the local surface. Re-relax with the same model and tighter criteria.
- Finite-displacement error. Force noise or an unsuitable amplitude corrupts the derivative. Test .
- Supercell truncation. Long-range force constants wrap into the finite repeat. Increase the physically relevant repeats.
- Potential error or extrapolation. The selected calculator gives an inaccurate curvature. Benchmark forces and frequencies for representative structures.
- Missing long-range electrostatics. A polar material needs non-analytic treatment near .
- Genuine harmonic instability. The mode persists in magnitude and eigenvector character after numerical, structural, and model checks.
The location of the instability is physical information. A mode preserves the primitive translational periodicity; a zone-boundary mode implies a larger real-space modulation. Following a soft-mode eigenvector into a lower-symmetry supercell can reveal a distorted phase, but that structure must be relaxed and its stability recalculated.
Three near-zero acoustic branches at are expected only for a complete crystal calculation. Because Tako applies the acoustic sum rule, exact or small acoustic offsets cannot independently validate the raw forces. Check optical modes, branch smoothness, supercell convergence, and independent force evidence as well.
Execution sequence and partial results
For a periodic calculation, Tako performs the following sequence:
- Parse the structure and require a valid periodic cell.
- Standardize it to a SeekPath primitive cell.
- Normalize repeats to 1–6, displacement to a positive value, and band points to 1–64.
- Construct the calculator for the standardized primitive/supercell representation.
- Evaluate central finite displacements and condition the force constants.
- Sum force-constant images to form the Hessian and emit the first partial result containing
phonon, top-levelfrequencies_cm_inv, andmodes. - Build and emit the positive-mode histogram as
phonon_dos. - Construct the custom or SeekPath band. If path construction or band evaluation fails, log the reason and finish without
phonon_band; otherwise emit a final partial result with the band.
This staged behavior means a calculation can possess valid data and DOS while lacking a band object. Consumers must test field availability rather than assume every successful return has all three. The public result does not contain the image-resolved force constants, displaced structures, or a restart cache, even though the lower-level library represents those objects internally.
For a nonperiodic input, the operation deliberately falls back to molecular finite-displacement vibrations with a two-point stencil. It returns kind: "phonon", phonon, top-level frequencies/modes, and a molecular-style DOS, but no phonon_cell and no band. Use the explicit vibration operation for molecules; a fallback is compatibility behavior, not a route to crystallographic bands.
Cost and memory
If the standardized primitive cell has atoms, of them active, and supercell volume , the central-difference work is approximately
The calculator cost per evaluation can scale roughly linearly with atoms for a local MLIP at fixed cutoff, but more steeply for self-consistent electronic methods. Memory also includes the repeated coordinates, calculator workspace, forces, compact force-constant blocks, and repeated partial-result serialization. Browser feasibility can fail before a formal repeat cap is reached.
Primitive-cell reduction can drastically lower , while a larger supercell raises the per-evaluation atom count. Band points are cheap relative to force evaluation because they reuse the same force constants, although dense diagonalization still scales with the cube of the active dynamical-matrix dimension.
Convergence design
A defensible study varies one numerical layer at a time:
| Layer | Controlled comparison | Evidence |
|---|---|---|
| Geometry | tighter atomic and cell relaxation | important branches stabilize; residual-force artifacts diminish |
| Displacement | at least two amplitudes | frequencies and eigenvector subspaces remain stable |
| Supercell | nested, physically motivated repeat sizes | branch frequencies cease changing within the required tolerance |
| Path | same standardized cell and labels | visual comparison samples identical reciprocal lines |
| Reciprocal coverage | path plus an appropriate independent mesh | no missed general- instability for a stability claim |
| Calculator | matched reference force/frequency tests | potential is credible for the material and distorted configurations |
| Polar correction | independent NAC-capable reference where needed | LO–TO behavior is not misassigned |
Changing the cell, calculator, supercell, displacement, and path together does not establish which error was reduced. Preserve raw result objects and state the comparison tolerance. Near-degenerate branches can exchange ordering; compare subspaces and displacement character rather than only array indices.
Limitations and reporting
Tako’s public phonon operation is a browser-local harmonic direct method. It does not expose third- or fourth-order force constants, phonon lifetimes, thermal conductivity, self-consistent phonons, quasi-harmonic volume sampling, isotope disorder, electron–phonon coupling, neutron or Raman cross sections, q-mesh thermodynamics, Born-charge NAC, or restartable displacement caches. It does not automatically search for a lower-symmetry phase after finding an instability.
Report the input structure and primitive-cell transformation, periodicity, calculator/model and asset provenance, charge/spin settings where relevant, atomic and cell relaxation criteria, residual force and stress, constraints/active atoms, supercell repeats, displacement, force-constant conditioning, path source and string, points per segment, frequency units, whether the DOS is the current histogram, missing NAC, imaginary-mode convention, convergence comparisons, and every omitted result field.
For a stability conclusion, include the minimum sampled frequency and its location, eigenvector character, supercell/displacement sensitivity, and reciprocal-space coverage. “No imaginary modes on the plotted path” is narrower than “dynamically stable throughout the Brillouin zone” and must be written that way.
References
- M. T. Dove, Introduction to Lattice Dynamics, Cambridge University Press, Cambridge (1993).
- M. Born, K. Huang, Dynamical Theory of Crystal Lattices, Oxford University Press, Oxford (1954).
- K. Parlinski, Z. Q. Li, Y. Kawazoe, “First-Principles Determination of the Soft Mode in Cubic ZrO₂,” Phys. Rev. Lett. 78, 4063–4066 (1997).
- G. Kresse, J. Furthmüller, J. Hafner, “Ab initio Force Constant Approach to Phonon Dispersion Relations of Diamond and Graphite,” Europhys. Lett. 32, 729–734 (1995).
- A. Togo, I. Tanaka, “First Principles Phonon Calculations in Materials Science,” Scr. Mater. 108, 1–5 (2015).
- Y. Hinuma, G. Pizzi, Y. Kumagai, F. Oba, I. Tanaka, “Band Structure Diagram Paths Based on Crystallography,” Comput. Mater. Sci. 128, 140–184 (2017).
- G. Pizzi et al., “Seekpath: A Python Tool for Uniform k-Point Path Selection and Crystal Structure Standardization,” Comput. Mater. Sci. 128, 436–446 (2017).
- A. H. Larsen et al., “The Atomic Simulation Environment—A Python Library for Working with Atoms,” J. Phys.: Condens. Matter 29, 273002 (2017).