日本語 | English

Chapter 1: First-Principles Phonon Calculations

DFPT, Frozen Phonon Method, and Modern Computational Tools

Learning Objectives

  • Understand the theoretical foundations of Density Functional Perturbation Theory (DFPT)
  • Master linear response theory and the 2n+1 theorem
  • Implement and compare frozen phonon and DFPT methods
  • Calculate force constants and perform Fourier interpolation
  • Use modern software tools (Quantum ESPRESSO, VASP, Phonopy) effectively
  • Understand convergence parameters (k-points, q-points, supercell size)
  • Apply non-analytic corrections for LO-TO splitting in polar materials
  • Develop complete Python workflows for phonon calculations

Introduction

First-principles phonon calculations have revolutionized materials science by enabling accurate prediction of vibrational properties without empirical parameters. Unlike classical force-field approaches, ab initio methods compute phonon properties directly from electronic structure calculations based on density functional theory (DFT).

This chapter covers the two main approaches to first-principles phonon calculations: Density Functional Perturbation Theory (DFPT) and the frozen phonon (finite displacement) method. We will explore their theoretical foundations, practical implementations, advantages and limitations, and integration with modern software packages. Understanding these methods is essential for computational materials research and phonon engineering.

1.1 Density Functional Perturbation Theory (DFPT)

1.1.1 Motivation and Overview

The dynamical matrix requires second derivatives of the total energy with respect to atomic displacements. A naive approach would use finite differences, requiring many self-consistent calculations for displaced configurations. DFPT provides an elegant alternative by computing the response of the electronic structure to perturbations analytically.

Definition: DFPT Philosophy

DFPT calculates the response of the electronic ground state to external perturbations (atomic displacements, electric fields) self-consistently, yielding second-order derivatives of the energy directly without computing finite differences.

1.1.2 Linear Response Theory

Consider a perturbation to the Hamiltonian:

$$\hat{H} = \hat{H}_0 + \lambda \hat{V}$$

where $\hat{H}_0$ is the unperturbed Hamiltonian and $\lambda$ is the perturbation strength. The electronic ground state energy can be expanded as:

$$E = E^{(0)} + \lambda E^{(1)} + \frac{1}{2}\lambda^2 E^{(2)} + \mathcal{O}(\lambda^3)$$

The density matrix also expands perturbatively:

$$\rho(\mathbf{r}) = \rho^{(0)}(\mathbf{r}) + \lambda \rho^{(1)}(\mathbf{r}) + \frac{1}{2}\lambda^2 \rho^{(2)}(\mathbf{r}) + \cdots$$

1.1.3 The 2n+1 Theorem

Theorem: 2n+1 Theorem (Wigner, 1935)

If the ground state wavefunction (or density) is known to order $n$ in the perturbation, the energy can be computed to order $2n+1$.

Corollary for Phonons:

Knowing only the unperturbed ground state ($n=0$), we can compute second-order energy derivatives exactly by solving for the first-order response of the wavefunction/density.

Proof Sketch:

The energy functional in DFT is:

$$E[\rho] = T[\rho] + \int V_\text{ext}(\mathbf{r})\rho(\mathbf{r})d\mathbf{r} + E_\text{H}[\rho] + E_\text{xc}[\rho]$$

The second-order energy change involves:

$$E^{(2)} = \int \frac{\delta^2 E}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} \rho^{(1)}(\mathbf{r})\rho^{(1)}(\mathbf{r}')d\mathbf{r}d\mathbf{r}' + \int \frac{\delta E}{\delta\rho(\mathbf{r})}\rho^{(2)}(\mathbf{r})d\mathbf{r}$$

However, the variational principle ensures $\delta E/\delta\rho = 0$ at the ground state, so the second term vanishes! Thus, $E^{(2)}$ depends only on $\rho^{(1)}$, not $\rho^{(2)}$.

1.1.4 Self-Consistent Phonon Equation

For an atomic displacement perturbation $\lambda u_{\kappa\alpha}$ at atom $\kappa$ in direction $\alpha$, the first-order change in the Kohn-Sham potential must be solved self-consistently:

$$\delta V_\text{KS}(\mathbf{r}) = \delta V_\text{ext}(\mathbf{r}) + \int \frac{\delta^2 E_\text{H}}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} \delta\rho(\mathbf{r}')d\mathbf{r}' + \int \frac{\delta^2 E_\text{xc}}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} \delta\rho(\mathbf{r}')d\mathbf{r}'$$

where the first-order density change is:

$$\delta\rho(\mathbf{r}) = \sum_{n,\mathbf{k}} f_{n\mathbf{k}} \left[\psi^*_{n\mathbf{k}}(\mathbf{r})\delta\psi_{n\mathbf{k}}(\mathbf{r}) + \text{c.c.}\right]$$

The first-order wavefunctions satisfy:

$$(\hat{H}_\text{KS}^{(0)} - \epsilon_{n\mathbf{k}})\delta\psi_{n\mathbf{k}} = -(\delta V_\text{KS} - \delta\epsilon_{n\mathbf{k}})\psi_{n\mathbf{k}}^{(0)}$$

This is solved iteratively until self-consistency, yielding $\delta\rho$ and ultimately the force constant matrix via:

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

Computational Advantage

DFPT requires solving the linear response equation for each phonon mode, which scales as $\mathcal{O}(N_\text{atoms} \times N_\text{q-points})$. For each $\mathbf{q}$-point, one self-consistent calculation yields the entire dynamical matrix, making it highly efficient for commensurate phonon wavevectors.

1.2 Frozen Phonon (Finite Displacement) Method

1.2.1 Conceptual Framework

The frozen phonon method is conceptually simpler: displace atoms in a supercell by small amounts, compute forces (negative energy gradients) via DFT, and extract force constants from the finite-difference approximation.

Definition: Frozen Phonon Approach

For a displacement $u_{\kappa'\alpha'}$ of atom $\kappa'$ in direction $\alpha'$, the force constant is approximated by:

$$\Phi_{\alpha\alpha'}(\kappa, \kappa') \approx -\frac{F_\alpha(\kappa; u_{\kappa'\alpha'}) - F_\alpha(\kappa; 0)}{u_{\kappa'\alpha'}}$$

where $F_\alpha(\kappa; u)$ is the force on atom $\kappa$ in direction $\alpha$ when displacement $u$ is applied.

1.2.2 Supercell Construction

To compute force constants in real space, we construct a supercell of size $N_1 \times N_2 \times N_3$ unit cells. The supercell must be large enough to:

  • Ensure force constants decay to zero at the supercell boundaries
  • Avoid artificial interactions between periodic images
  • Capture long-range interactions in polar materials

Example: Silicon (Diamond Structure)

For silicon with a diamond structure (2 atoms per primitive cell), a $3 \times 3 \times 3$ supercell contains 54 atoms. With 3 Cartesian directions, we have 162 degrees of freedom. However, symmetry reduces the number of independent displacements dramatically.

Typical procedure:

  1. Displace one atom in the supercell by $\pm\delta u$ in each Cartesian direction
  2. Compute forces on all atoms via DFT
  3. Apply central differences: $\Phi \approx [F(+\delta u) - F(-\delta u)]/(2\delta u)$
  4. Use symmetry to generate remaining force constants

1.2.3 Displacement Magnitude

The displacement magnitude $\delta u$ must be chosen carefully:

  • Too small: Numerical noise dominates, poor force accuracy
  • Too large: Anharmonic effects contaminate harmonic force constants
  • Optimal range: Typically 0.01 – 0.03 Å for most materials

Higher-order finite differences can improve accuracy:

$$\Phi \approx \frac{-F(+2\delta u) + 8F(+\delta u) - 8F(-\delta u) + F(-2\delta u)}{12\delta u} + \mathcal{O}(\delta u^4)$$

1.2.4 Workflow Diagram

graph TD A[Optimized Crystal Structure] --> B[Construct Supercell] B --> C[Apply Symmetry Analysis] C --> D[Generate Displacement Patterns] D --> E{For Each Displacement} E --> F[Run DFT SCF Calculation] F --> G[Extract Forces on All Atoms] G --> E E --> H[Collect Force Constants] H --> I[Apply Acoustic Sum Rule] I --> J[Fourier Transform to q-space] J --> K[Compute Phonon Dispersion]

1.2.5 Advantages and Limitations

Aspect Advantages Limitations
Conceptual Simplicity Intuitive, easy to understand and implement —
Software Compatibility Works with any DFT code (only needs forces) —
Computational Cost Can exploit symmetry to reduce calculations Requires many DFT runs for large supercells
q-point Sampling Provides real-space force constants for all q Limited q-resolution by supercell size
Numerical Accuracy Higher-order finite differences available Sensitive to force convergence and $\delta u$
Anharmonicity Can extend to third/fourth-order force constants Harmonic approximation may break for large $\delta u$

1.3 DFPT vs. Frozen Phonon: Detailed Comparison

1.3.1 Accuracy Considerations

Accuracy Comparison

DFPT:

  • Exact within DFT for infinitesimal perturbations
  • No finite-difference errors
  • Limited by DFT functional approximations only

Frozen Phonon:

  • Finite-difference truncation error $\mathcal{O}(\delta u^2)$ or higher
  • Requires tight force convergence (typically $< 10^{-5}$ eV/Å)
  • Anharmonic contamination for large displacements

In practice, with careful convergence testing, both methods yield nearly identical results for harmonic phonon frequencies (differences $< 1$ cm$^{-1}$).

1.3.2 Computational Cost Analysis

For a system with $N_\text{atoms}$ atoms per primitive cell and $N_\mathbf{q}$ q-points:

Method Number of SCF Calculations Cost per Calculation Total Scaling
DFPT $N_\mathbf{q}$ (one per q-point) $\sim 3 \times$ ground state $\mathcal{O}(N_\mathbf{q})$
Frozen Phonon $6N_\text{atoms}^{\text{super}}$ (symmetry-reduced) Standard ground state $\mathcal{O}(N_\text{atoms}^{\text{super}})$

Key Trade-offs:

  • Few q-points needed: DFPT is more efficient
  • Dense q-grid needed: Frozen phonon may be better (compute once, interpolate everywhere)
  • Large primitive cell: Frozen phonon often faster
  • Small primitive cell: DFPT typically faster

1.3.3 Special Cases and Best Practices

When to Use Each Method

Prefer DFPT when:

  • Working with simple structures (1–4 atoms per cell)
  • Need phonons at specific high-symmetry q-points
  • Computing dielectric responses and Born charges
  • Polar materials requiring non-analytic corrections

Prefer Frozen Phonon when:

  • Large primitive cells (> 10 atoms)
  • Software doesn't support DFPT (e.g., some DFT codes)
  • Need third/fourth-order anharmonic force constants
  • Studying disordered or defect systems

1.4 Force Constant Calculation and Fourier Interpolation

1.4.1 Real-Space Force Constants

Both DFPT and frozen phonon methods ultimately provide force constants in real space:

$$\Phi_{\alpha\beta}(\mathbf{R}_0 + \boldsymbol{\tau}_\kappa, \mathbf{R}_n + \boldsymbol{\tau}_{\kappa'})$$

where $\mathbf{R}_n$ is a lattice vector, $\boldsymbol{\tau}_\kappa$ is the position of atom $\kappa$ within the unit cell, and $\alpha, \beta$ are Cartesian indices.

Due to translational symmetry, force constants depend only on the relative lattice vector:

$$\Phi_{\alpha\beta}(\mathbf{R}_0 + \boldsymbol{\tau}_\kappa, \mathbf{R}_n + \boldsymbol{\tau}_{\kappa'}) = \Phi_{\alpha\beta}(\mathbf{R}_n; \kappa, \kappa')$$

1.4.2 Fourier Transform to q-Space

The dynamical matrix at wavevector $\mathbf{q}$ is obtained via Fourier transform:

$$D_{\kappa\alpha,\kappa'\beta}(\mathbf{q}) = \frac{1}{\sqrt{M_\kappa M_{\kappa'}}} \sum_{\mathbf{R}_n} \Phi_{\alpha\beta}(\mathbf{R}_n; \kappa, \kappa') e^{i\mathbf{q}\cdot(\mathbf{R}_n + \boldsymbol{\tau}_{\kappa'} - \boldsymbol{\tau}_\kappa)}$$

This allows computing phonons at any $\mathbf{q}$-point in the Brillouin zone from a finite set of real-space force constants.

1.4.3 Fourier Interpolation

The frozen phonon method with an $N_1 \times N_2 \times N_3$ supercell directly yields force constants on a commensurate grid:

$$\mathbf{q} = \frac{m_1}{N_1}\mathbf{b}_1 + \frac{m_2}{N_2}\mathbf{b}_2 + \frac{m_3}{N_3}\mathbf{b}_3$$

To obtain phonons at arbitrary $\mathbf{q}$, we use Fourier interpolation:

  1. Compute real-space force constants from supercell displacements
  2. Fourier transform to get $D(\mathbf{q})$ at any desired $\mathbf{q}$
  3. Diagonalize to obtain $\omega_j(\mathbf{q})$ and eigenvectors

Convergence with Supercell Size

The accuracy of Fourier interpolation depends on how well the supercell captures the range of force constants. If $\Phi(\mathbf{R})$ has not decayed to zero at the supercell boundary, phonon frequencies will show unphysical oscillations. Always test convergence by increasing supercell size until phonon frequencies change by less than your target accuracy (e.g., 1 cm$^{-1}$).

1.4.4 Acoustic Sum Rule Enforcement

Numerical errors can violate the acoustic sum rule:

$$\sum_{\mathbf{R}_n, \kappa'} \Phi_{\alpha\beta}(\mathbf{R}_n; \kappa, \kappa') = 0$$

This results in spurious imaginary frequencies at the $\Gamma$ point. Post-processing tools (e.g., Phonopy) can enforce the sum rule by symmetrizing force constants:

$$\Phi^\text{corrected}_{\alpha\beta}(0; \kappa, \kappa) = -\sum_{\mathbf{R}_n \neq 0, \kappa'} \Phi_{\alpha\beta}(\mathbf{R}_n; \kappa, \kappa')$$

1.5 Modern Software Tools

1.5.1 Quantum ESPRESSO (PHonon Package)

Quantum ESPRESSO provides the PHonon package implementing DFPT for phonon calculations. It uses plane-wave basis sets and pseudopotentials.

Quantum ESPRESSO: Basic Phonon Workflow
# Step 1: Ground state SCF calculation pw.x < scf.in > scf.out # scf.in structure: &control calculation = 'scf' prefix = 'silicon' outdir = './tmp/' / &system ibrav = 2, celldm(1) = 10.26 nat = 2, ntyp = 1 ecutwfc = 60.0 / &electrons conv_thr = 1.0d-10 / ATOMIC_SPECIES Si 28.086 Si.pbe-n-rrkjus_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) 12 12 12 0 0 0 # Step 2: Phonon calculation at Gamma point ph.x < ph.in > ph.out # ph.in: &inputph prefix = 'silicon' outdir = './tmp/' fildyn = 'silicon.dyn' tr2_ph = 1.0d-14 / 0.0 0.0 0.0 ! q-point (Gamma) # Step 3: Compute phonon dispersion along high-symmetry path # Use q2r.x and matdyn.x for Fourier interpolation

1.5.2 VASP with Phonopy

VASP does not include DFPT for phonons, so the frozen phonon method with Phonopy is standard.

VASP + Phonopy: Complete Workflow
#!/bin/bash # Phonon calculation workflow with VASP and Phonopy # Step 1: Structure optimization # Run VASP with ISIF=3 to optimize cell and atomic positions # (INCAR, POSCAR, POTCAR, KPOINTS prepared) # Step 2: Create SPOSCAR (supercell) with Phonopy phonopy -d --dim="3 3 3" -c POSCAR # This creates: # SPOSCAR: supercell structure # phonopy_disp.yaml: displacement information # POSCAR-{001,002,...}: displaced configurations # Step 3: Calculate forces for each displacement for i in POSCAR-{001..027}; do mkdir disp-${i##*-} cp $i disp-${i##*-}/POSCAR cd disp-${i##*-} # Copy INCAR, POTCAR, KPOINTS # INCAR must have: NSW=0, IBRION=-1, PREC=Accurate mpirun -np 16 vasp_std > log cd .. done # Step 4: Collect forces and create FORCE_SETS phonopy -f disp-*/vasprun.xml # Step 5: Compute phonon band structure # Create band.conf: cat > band.conf << EOF ATOM_NAME = Si DIM = 3 3 3 BAND = 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.5 BAND_POINTS = 101 EOF phonopy -p band.conf # Generates band.yaml and band.pdf

1.5.3 ABINIT DFPT

ABINIT has comprehensive DFPT implementation including electric field perturbations for Born effective charges and dielectric tensors.

ABINIT: Phonon Calculation Example
# Dataset 1: Ground state ndtset 2 jdtset 1 2 # Common parameters nband 8 ecut 30.0 nstep 100 toldfe 1.0d-10 # Dataset 1: Ground state SCF kptopt1 1 ngkpt1 8 8 8 # Dataset 2: DFPT phonon calculation kptopt2 3 ngkpt2 8 8 8 qpt2 0.25 0.0 0.0 # Q-point for phonon rfphon2 1 rfatpol2 1 2 rfdir2 1 1 1 tolvrs2 1.0d-10

1.5.4 Phonopy: Universal Post-Processing

Phonopy is a Python-based tool that interfaces with multiple DFT codes (VASP, Quantum ESPRESSO, ABINIT, CASTEP, etc.) for phonon analysis.

Python: Phonopy API Usage
#!/usr/bin/env python """ Phonopy API example: Load force constants and compute phonon DOS """ from phonopy import Phonopy from phonopy.interface.vasp import read_vasp from phonopy.file_IO import parse_FORCE_SETS, parse_BORN import numpy as np import matplotlib.pyplot as plt # Read crystal structure unitcell = read_vasp("POSCAR") # Create Phonopy object phonon = Phonopy(unitcell, supercell_matrix=[[3, 0, 0], [0, 3, 0], [0, 0, 3]]) # Load force constants force_sets = parse_FORCE_SETS() phonon.dataset = force_sets phonon.produce_force_constants() # For polar materials: load Born charges and dielectric tensor try: nac_params = parse_BORN(phonon.primitive, filename="BORN") phonon.nac_params = nac_params except FileNotFoundError: print("BORN file not found, skipping non-analytic corrections") # Set q-mesh for DOS calculation mesh = [20, 20, 20] phonon.run_mesh(mesh) phonon.run_total_dos() # Plot phonon DOS plt.figure(figsize=(8, 6)) plt.plot(phonon.total_dos.frequency_points, phonon.total_dos.dos, linewidth=2) plt.xlabel('Frequency (THz)', fontsize=12) plt.ylabel('Density of States', fontsize=12) plt.title('Phonon Density of States', fontsize=14) plt.grid(True, alpha=0.3) plt.xlim(0, None) plt.tight_layout() plt.savefig('phonon_dos.png', dpi=150) plt.show() # Compute thermal properties temperatures = np.arange(0, 1001, 10) phonon.run_thermal_properties(temperatures) tp = phonon.get_thermal_properties_dict() # Plot heat capacity plt.figure(figsize=(8, 6)) plt.plot(tp['temperatures'], tp['heat_capacity'], 'r-', linewidth=2) plt.xlabel('Temperature (K)', fontsize=12) plt.ylabel('Heat Capacity (J/K/mol)', fontsize=12) plt.title('Phonon Heat Capacity', fontsize=14) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('heat_capacity.png', dpi=150) plt.show()

1.6 Convergence Parameters

1.6.1 Electronic Structure Convergence

Phonon calculations require tighter convergence than standard DFT calculations:

Parameter Standard DFT Phonon Calculations Rationale
Energy cutoff Converged to ~1 meV Converged to ~0.1 meV Forces are energy derivatives
k-point grid 0.03 Å$^{-1}$ spacing 0.02 Å$^{-1}$ spacing or denser Force convergence requires dense k-mesh
SCF threshold $10^{-6}$ eV $10^{-8}$ – $10^{-10}$ eV Tight convergence for accurate forces
Force threshold $10^{-3}$ eV/Å $10^{-5}$ – $10^{-6}$ eV/Å Essential for frozen phonon method

1.6.2 q-Point Sampling (DFPT)

In DFPT, each q-point requires a separate calculation. The required q-mesh depends on:

  • Purpose: Band structure (high-symmetry path) vs. DOS (dense grid)
  • Force constant range: Long-range interactions need dense q-mesh
  • Material type: Metals require denser meshes than insulators

Typical q-meshes:

  • Band structure: 20–50 points along high-symmetry paths
  • DOS: $20 \times 20 \times 20$ to $40 \times 40 \times 40$ grids
  • Thermal properties: Same as DOS requirements

1.6.3 Supercell Size (Frozen Phonon)

Convergence testing procedure:

  1. Start with $2 \times 2 \times 2$ supercell
  2. Calculate phonon frequencies at several q-points
  3. Increase to $3 \times 3 \times 3$, then $4 \times 4 \times 4$
  4. Check if frequencies change by less than target accuracy (e.g., 1 cm$^{-1}$ or 0.03 THz)
  5. If not converged, continue increasing until convergence

Material-Dependent Convergence

  • Covalent materials (Si, C): $3 \times 3 \times 3$ often sufficient
  • Ionic materials (NaCl, MgO): $4 \times 4 \times 4$ or larger (long-range Coulomb)
  • Metals: $5 \times 5 \times 5$ or larger (screening length)
  • Layered materials (graphene, MoS$_2$): Anisotropic supercells (e.g., $5 \times 5 \times 2$)

1.6.4 Displacement Magnitude Convergence

For frozen phonon calculations, test displacement magnitude:

Python: Displacement Convergence Test
import numpy as np import matplotlib.pyplot as plt # Test different displacement magnitudes displacements = np.array([0.005, 0.01, 0.015, 0.02, 0.025, 0.03, 0.04, 0.05]) # Angstrom # Example: phonon frequency at Gamma point (TO mode) for different displacements # (In practice, these would come from actual DFT calculations) frequencies_example = np.array([15.2, 15.25, 15.28, 15.29, 15.30, 15.30, 15.35, 15.45]) # THz plt.figure(figsize=(8, 6)) plt.plot(displacements, frequencies_example, 'o-', linewidth=2, markersize=8) plt.xlabel('Displacement Magnitude (Å)', fontsize=12) plt.ylabel('Phonon Frequency (THz)', fontsize=12) plt.title('Convergence of Phonon Frequency with Displacement', fontsize=14) plt.grid(True, alpha=0.3) plt.axhline(y=15.30, color='r', linestyle='--', linewidth=1, label='Converged value') plt.legend(fontsize=11) plt.tight_layout() plt.savefig('displacement_convergence.png', dpi=150) plt.show() # Optimal displacement: plateau region (0.015-0.03 Å typical)

1.7 LO-TO Splitting and Non-Analytic Corrections

1.7.1 Physical Origin

In polar materials (ionic crystals, polar semiconductors), long-range Coulomb interactions between ions cause longitudinal optical (LO) and transverse optical (TO) phonons to split at the $\Gamma$ point ($\mathbf{q} \to 0$).

The dynamical matrix has a non-analytic term proportional to $\mathbf{q}\mathbf{q}^T/|\mathbf{q}|^2$:

$$D^{\text{NA}}_{\kappa\alpha,\kappa'\beta}(\mathbf{q}) = \frac{4\pi}{\Omega\epsilon_\infty} \frac{(\mathbf{q} \cdot \mathbf{Z}^*_\kappa)_\alpha (\mathbf{q} \cdot \mathbf{Z}^*_{\kappa'})_\beta}{|\mathbf{q}|^2}$$

where:

  • $\Omega$ is the unit cell volume
  • $\epsilon_\infty$ is the high-frequency dielectric tensor
  • $\mathbf{Z}^*_\kappa$ is the Born effective charge tensor of atom $\kappa$

1.7.2 Direction Dependence

The non-analytic term depends on the direction of $\mathbf{q}$ as $\mathbf{q} \to 0$, leading to different limiting frequencies for different directions:

$$\lim_{\mathbf{q} \to 0} \omega(\mathbf{q}) \neq \omega(\Gamma)$$

This causes discontinuities in phonon dispersion plots near $\Gamma$ unless properly handled.

Example: GaAs (Zinc Blende Structure)

GaAs has two atoms per primitive cell (Ga, As), both with significant Born charges. At the $\Gamma$ point:

  • TO mode: $\omega_{\text{TO}}(\Gamma) \approx 8.0$ THz (transverse motion)
  • LO mode: $\omega_{\text{LO}}(\Gamma) \approx 8.8$ THz (longitudinal motion)
  • Splitting: $\Delta\omega \approx 0.8$ THz due to macroscopic electric field

1.7.3 Computing Born Charges and Dielectric Tensor

DFPT can compute Born effective charges and dielectric tensors directly:

Quantum ESPRESSO: Born Charges Calculation
# ph.in for Born charges and dielectric tensor &inputph prefix = 'gaas' outdir = './tmp/' fildyn = 'gaas.dyn' epsil = .true. ! Compute dielectric tensor zeu = .true. ! Compute Born effective charges tr2_ph = 1.0d-14 / 0.0 0.0 0.0 # Output contains: # - High-frequency dielectric tensor epsilon_infinity # - Born effective charge tensors Z* for each atom

For VASP (frozen phonon approach), use DFPT for dielectric properties:

VASP: BORN File for Phonopy
# BORN file format for Phonopy # Line 1: Default (scaling factors, usually 1.0) default # Lines 2-4: High-frequency dielectric tensor (3x3) 12.40 0.00 0.00 0.00 12.40 0.00 0.00 0.00 12.40 # Subsequent lines: Born effective charge tensor for each atom (3x3) # Atom 1 (Ga): 2.12 0.00 0.00 0.00 2.12 0.00 0.00 0.00 2.12 # Atom 2 (As): -2.12 0.00 0.00 0.00 -2.12 0.00 0.00 0.00 -2.12

1.7.4 Implementing Non-Analytic Corrections

Python: Non-Analytic Correction in Dynamical Matrix
import numpy as np def add_nac_correction(D, q, born_charges, epsilon_inf, volume, masses): """ Add non-analytic correction to dynamical matrix for polar materials Parameters: ----------- D : ndarray (3N, 3N) Dynamical matrix without NAC q : ndarray (3,) Wavevector (Cartesian coordinates) born_charges : list of ndarray (N x 3 x 3) Born effective charge tensors for each atom epsilon_inf : ndarray (3, 3) High-frequency dielectric tensor volume : float Unit cell volume masses : ndarray (N,) Atomic masses Returns: -------- D_corrected : ndarray (3N, 3N) Dynamical matrix with NAC included """ N_atoms = len(masses) q_norm = np.linalg.norm(q) if q_norm < 1e-8: # At Gamma, use a small finite q in desired direction # Or handle separately for each approach direction return D # q-direction unit vector q_hat = q / q_norm # Inverse dielectric tensor epsilon_inv = np.linalg.inv(epsilon_inf) # NAC contribution (4π/Ω in atomic units) prefactor = 4 * np.pi / volume D_nac = np.zeros_like(D, dtype=complex) for kappa in range(N_atoms): for kappa_p in range(N_atoms): # (q · Z*_κ) qZ_kappa = np.dot(q, born_charges[kappa]) # (q · Z*_κ') qZ_kappa_p = np.dot(q, born_charges[kappa_p]) # Construct 3x3 block for alpha in range(3): for beta in range(3): i = 3 * kappa + alpha j = 3 * kappa_p + beta D_nac[i, j] = prefactor / (q_norm**2) * \ qZ_kappa[alpha] * qZ_kappa_p[beta] / \ np.sqrt(masses[kappa] * masses[kappa_p]) return D + D_nac # Example usage # D = construct_dynamical_matrix(q, force_constants, masses, ...) # D_with_nac = add_nac_correction(D, q, born_charges, epsilon_inf, volume, masses)

1.8 Complete Python Workflow Example

1.8.1 End-to-End Phonon Calculation

Here we present a comprehensive example combining all concepts: force constant loading, Fourier interpolation, NAC corrections, and band structure plotting.

Python: Complete Phonon Calculation Workflow
#!/usr/bin/env python """ Complete phonon calculation workflow for silicon Demonstrates: Force constants, Fourier interpolation, band structure """ import numpy as np import matplotlib.pyplot as plt from phonopy import Phonopy from phonopy.interface.vasp import read_vasp from phonopy.file_IO import parse_FORCE_SETS from phonopy.structure.atoms import PhonopyAtoms def create_silicon_structure(): """Create silicon crystal structure (diamond)""" lattice = [[0, 2.73, 2.73], [2.73, 0, 2.73], [2.73, 2.73, 0]] # FCC lattice vectors, Angstrom positions = [[0, 0, 0], [0.25, 0.25, 0.25]] # Diamond structure atoms = PhonopyAtoms(symbols=['Si', 'Si'], cell=lattice, scaled_positions=positions) return atoms def compute_band_structure(phonon, path_labels, num_points=51): """ Compute phonon band structure along high-symmetry path Parameters: ----------- phonon : Phonopy object path_labels : list of tuples [(label, q-point), ...] defining high-symmetry path num_points : int Number of points per segment Returns: -------- distances : ndarray Cumulative distance along path frequencies : ndarray (n_bands, n_points) Phonon frequencies in THz """ # Extract path labels = [label for label, _ in path_labels] qpoints = [q for _, q in path_labels] # Generate path bands = [] for i in range(len(qpoints) - 1): segment = np.linspace(qpoints[i], qpoints[i+1], num_points) if i > 0: segment = segment[1:] # Avoid double-counting bands.append(segment) qpoints_path = np.vstack(bands) # Compute phonon frequencies phonon.run_band_structure(qpoints_path, is_band_connection=False) band_dict = phonon.get_band_structure_dict() frequencies = band_dict['frequencies'] distances = band_dict['distances'] return distances, frequencies, labels, [sum(len(b)-1 for b in bands[:i]) for i in range(len(bands)+1)] def plot_band_structure(distances, frequencies, labels, label_positions): """Plot phonon band structure""" fig, ax = plt.subplots(figsize=(10, 7)) # Plot all branches for band in range(frequencies.shape[1]): ax.plot(distances, frequencies[:, band], 'b-', linewidth=1.5) # High-symmetry point labels ax.set_xlim(distances[0], distances[-1]) ax.set_xticks([distances[pos] for pos in label_positions]) ax.set_xticklabels(labels) # Vertical lines at high-symmetry points for pos in label_positions[1:-1]: ax.axvline(x=distances[pos], color='k', linewidth=0.5, linestyle='--') ax.set_ylabel('Frequency (THz)', fontsize=14) ax.set_title('Silicon Phonon Band Structure', fontsize=16, fontweight='bold') ax.grid(True, alpha=0.3, axis='y') ax.set_ylim(0, None) plt.tight_layout() plt.savefig('silicon_phonon_bands.png', dpi=200) plt.show() def main(): """Main workflow""" print("=== Silicon Phonon Calculation Workflow ===\n") # Step 1: Create structure print("Step 1: Creating silicon crystal structure...") unitcell = create_silicon_structure() # Step 2: Initialize Phonopy print("Step 2: Initializing Phonopy with 3x3x3 supercell...") phonon = Phonopy(unitcell, supercell_matrix=[[3, 0, 0], [0, 3, 0], [0, 0, 3]]) # Step 3: Load force constants # (In practice, this would come from VASP/QE calculations) print("Step 3: Loading force constants from FORCE_SETS...") try: force_sets = parse_FORCE_SETS() phonon.dataset = force_sets phonon.produce_force_constants() print(" ✓ Force constants loaded successfully") except FileNotFoundError: print(" ✗ FORCE_SETS not found. Using example for demonstration.") # For demonstration, use Phonopy's built-in example or skip return # Step 4: Verify acoustic sum rule print("\nStep 4: Checking acoustic sum rule at Gamma...") phonon.run_qpoints([[0, 0, 0]]) freqs_gamma = phonon.get_qpoints_dict()['frequencies'][0] acoustic_freqs = freqs_gamma[:3] print(f" Acoustic mode frequencies at Gamma: {acoustic_freqs} THz") if np.all(np.abs(acoustic_freqs) < 0.01): print(" ✓ Acoustic sum rule satisfied") else: print(" ⚠ Acoustic sum rule may be violated. Consider enforcement.") # Step 5: Compute band structure print("\nStep 5: Computing phonon band structure...") path_labels = [ ('L', [0.5, 0.5, 0.5]), ('Γ', [0.0, 0.0, 0.0]), ('X', [0.5, 0.0, 0.5]), ('K', [0.375, 0.375, 0.75]), ('Γ', [0.0, 0.0, 0.0]) ] distances, frequencies, labels, label_pos = compute_band_structure( phonon, path_labels, num_points=51) print(f" ✓ Computed {frequencies.shape[0]} q-points, {frequencies.shape[1]} bands") # Step 6: Plot print("\nStep 6: Plotting band structure...") plot_band_structure(distances, frequencies, labels, label_pos) print(" ✓ Plot saved to silicon_phonon_bands.png") # Step 7: Compute DOS print("\nStep 7: Computing phonon density of states...") phonon.run_mesh([20, 20, 20]) phonon.run_total_dos() dos_freq = phonon.total_dos.frequency_points dos = phonon.total_dos.dos plt.figure(figsize=(8, 6)) plt.plot(dos_freq, dos, 'r-', linewidth=2) plt.fill_between(dos_freq, dos, alpha=0.3, color='red') plt.xlabel('Frequency (THz)', fontsize=12) plt.ylabel('Density of States', fontsize=12) plt.title('Silicon Phonon DOS', fontsize=14) plt.grid(True, alpha=0.3) plt.xlim(0, None) plt.tight_layout() plt.savefig('silicon_phonon_dos.png', dpi=150) plt.show() print(" ✓ DOS saved to silicon_phonon_dos.png") # Step 8: Thermal properties print("\nStep 8: Computing thermal properties...") temperatures = np.arange(0, 1001, 10) phonon.run_thermal_properties(temperatures) tp = phonon.get_thermal_properties_dict() fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # Heat capacity ax1.plot(tp['temperatures'], tp['heat_capacity'], 'b-', linewidth=2) ax1.set_xlabel('Temperature (K)', fontsize=12) ax1.set_ylabel('Heat Capacity (J/K/mol)', fontsize=12) ax1.set_title('Phonon Heat Capacity', fontsize=13) ax1.grid(True, alpha=0.3) # Free energy ax2.plot(tp['temperatures'], tp['free_energy'], 'g-', linewidth=2) ax2.set_xlabel('Temperature (K)', fontsize=12) ax2.set_ylabel('Free Energy (kJ/mol)', fontsize=12) ax2.set_title('Phonon Free Energy', fontsize=13) ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('silicon_thermal_properties.png', dpi=150) plt.show() print(" ✓ Thermal properties saved") print("\n=== Workflow Complete ===") if __name__ == '__main__': main()

1.9 Best Practices and Common Pitfalls

1.9.1 Best Practices Checklist

Essential Best Practices

  1. Always perform convergence tests
    • k-point mesh, energy cutoff, SCF threshold
    • Supercell size (frozen phonon) or q-mesh (DFPT)
    • Displacement magnitude (frozen phonon)
  2. Verify acoustic sum rule
    • Check that 3 acoustic modes have $\omega(\Gamma) \approx 0$
    • Use symmetry to reduce numerical errors
    • Apply post-processing corrections if needed
  3. Use appropriate pseudopotentials/PAW
    • Validated for phonon calculations (check literature)
    • Include semicore states if necessary (e.g., Ga 3d)
    • Test different pseudopotentials for consistency
  4. Handle polar materials carefully
    • Compute Born charges and dielectric tensor
    • Apply non-analytic corrections for LO-TO splitting
    • Use BORN file in Phonopy for correct dispersions
  5. Validate against experiments or literature
    • Compare key phonon frequencies to Raman/IR spectra
    • Check elastic constants from phonon dispersion slopes
    • Verify thermal properties against calorimetry data

1.9.2 Common Pitfalls to Avoid

Pitfall Symptom Solution
Imaginary frequencies everywhere Negative eigenvalues across entire BZ Structure not at equilibrium; re-optimize geometry with tight thresholds
Imaginary acoustic modes at Γ $\omega^2 < 0$ for acoustic modes Acoustic sum rule violated; enforce symmetry or use tighter convergence
Missing LO-TO splitting LO and TO degenerate at Γ in polar material Apply non-analytic corrections (BORN file)
Discontinuous bands Jumps in phonon dispersion Supercell too small (frozen phonon) or band connection issues
Inconsistent frequencies Different software/methods give different results Check unit conventions (THz vs. cm$^{-1}$), convergence, mass conventions
Huge computational cost Calculations take weeks Exploit symmetry, use coarser initial mesh, consider DFPT vs. frozen phonon

1.9.3 Debugging Imaginary Frequencies

Python: Diagnostic Tool for Imaginary Frequencies
import numpy as np from phonopy import Phonopy def diagnose_imaginary_frequencies(phonon, tolerance=1e-3): """ Diagnose imaginary phonon frequencies and suggest fixes Parameters: ----------- phonon : Phonopy object tolerance : float Threshold for identifying imaginary frequencies (THz) """ print("=== Imaginary Frequency Diagnostics ===\n") # Check Gamma point phonon.run_qpoints([[0, 0, 0]]) freqs_gamma = phonon.get_qpoints_dict()['frequencies'][0] imaginary_mask = freqs_gamma < -tolerance n_imaginary = np.sum(imaginary_mask) print(f"Gamma point frequencies: {freqs_gamma} THz") print(f"Number of imaginary modes: {n_imaginary}") if n_imaginary == 0: print("✓ No imaginary frequencies at Gamma\n") elif n_imaginary <= 3: print("⚠ Acoustic modes may be slightly imaginary (numerical error)") print(" → Check acoustic sum rule enforcement") print(f" → Imaginary frequencies: {freqs_gamma[imaginary_mask]}\n") else: print("✗ Significant imaginary frequencies detected!") print(" Possible causes:") print(" 1. Structure not fully optimized (forces > 1e-5 eV/Å)") print(" 2. Phonon instability (material unstable at 0K)") print(" 3. Numerical errors in force constants") print(" 4. Supercell too small (frozen phonon method)") print(" 5. k-point mesh too coarse in DFT calculations\n") # Check a few other q-points test_qpoints = [[0.1, 0.1, 0.1], [0.5, 0.5, 0.5], [0.5, 0.0, 0.5]] phonon.run_qpoints(test_qpoints) all_freqs = phonon.get_qpoints_dict()['frequencies'] for i, (q, freqs) in enumerate(zip(test_qpoints, all_freqs)): n_imag = np.sum(freqs < -tolerance) if n_imag > 0: print(f"q = {q}: {n_imag} imaginary modes") print(f" Min frequency: {freqs.min():.3f} THz") # Overall assessment all_imaginary = np.sum(all_freqs < -tolerance) total_modes = all_freqs.size print(f"\nOverall: {all_imaginary}/{total_modes} imaginary modes ({100*all_imaginary/total_modes:.1f}%)") if all_imaginary == 0: print("✓ All frequencies real and positive") elif all_imaginary < 0.05 * total_modes: print("⚠ Minor imaginary frequencies, likely numerical") else: print("✗ Significant instability or calculation error") # Usage: # diagnose_imaginary_frequencies(phonon)

Summary

This chapter provided a comprehensive introduction to first-principles phonon calculations, covering both theoretical foundations and practical implementations:

  • DFPT: Theoretical basis via linear response theory and the 2n+1 theorem, enabling direct computation of second-order energy derivatives
  • Frozen Phonon Method: Finite-displacement approach using supercells and force calculations, more flexible but computationally intensive
  • Comparison: Accuracy, computational cost, and appropriate use cases for each method
  • Force Constants: Real-space calculation, Fourier interpolation to arbitrary q-points, and acoustic sum rule enforcement
  • Software Tools: Practical workflows with Quantum ESPRESSO, VASP+Phonopy, ABINIT, and Python APIs
  • Convergence: Systematic testing of k-points, q-points, supercell size, and displacement magnitude
  • LO-TO Splitting: Non-analytic corrections for polar materials using Born charges and dielectric tensors
  • Best Practices: Validation procedures, common pitfalls, and debugging strategies

These methods form the foundation for all advanced phonon physics topics covered in subsequent chapters, including anharmonicity, thermal transport, and phonon engineering.

Exercises

Problem 1: DFPT Theory (Fundamental)

Prove that the 2n+1 theorem allows second-order energy derivatives to be computed from first-order wavefunctions. Start from the variational principle $\delta E[\psi]/\delta\psi = 0$ at the ground state and show that the second-order energy change $E^{(2)}$ does not depend on $\psi^{(2)}$.

Problem 2: Frozen Phonon Convergence (Intermediate)

For a 1D diatomic chain, numerically test the convergence of force constants with respect to displacement magnitude. Use displacements from 0.001 to 0.1 Å and implement both central differences and 5-point stencil methods. Compare accuracy and identify the optimal displacement range.

Problem 3: Supercell Size Effects (Intermediate)

Using Phonopy, calculate the phonon dispersion of silicon with supercell sizes $2 \times 2 \times 2$, $3 \times 3 \times 3$, and $4 \times 4 \times 4$. Compare the phonon frequencies at the X and L points. Determine the minimum supercell size needed for 1 cm$^{-1}$ accuracy. Explain why convergence differs for different q-points.

Problem 4: LO-TO Splitting Calculation (Advanced)

For GaAs:

  1. Compute the Born effective charges and high-frequency dielectric tensor using DFPT
  2. Implement the non-analytic correction formula in Python
  3. Plot the phonon dispersion near the $\Gamma$ point along the [100], [110], and [111] directions, showing the LO-TO splitting explicitly
  4. Verify the Lyddane-Sachs-Teller relation: $\epsilon_0/\epsilon_\infty = (\omega_{\text{LO}}/\omega_{\text{TO}})^2$

Problem 5: Acoustic Sum Rule Enforcement (Advanced)

Given a set of force constants that violate the acoustic sum rule:

  1. Quantify the violation by computing $\sum_{n', \kappa'} \Phi_{\alpha\beta}(0\kappa, n'\kappa')$ for each atom and direction
  2. Implement a symmetry-preserving correction algorithm that enforces the sum rule while maintaining crystallographic symmetry
  3. Verify that the corrected force constants yield zero acoustic frequencies at $\Gamma$
  4. Compare the phonon DOS before and after correction

Problem 6: Complete Workflow Implementation (Research-Level)

Choose a material system (e.g., diamond, wurtzite GaN, or cubic perovskite) and perform a complete phonon calculation:

  1. Structure optimization with tight convergence criteria
  2. Convergence tests for all relevant parameters
  3. Phonon calculations using both DFPT and frozen phonon methods
  4. Comparison of results and computational cost
  5. Application of NAC corrections if applicable
  6. Validation against experimental Raman/IR spectra or literature
  7. Computation of thermodynamic properties and comparison with experimental heat capacity

References

  1. S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, "Phonons and related crystal properties from density-functional perturbation theory", Rev. Mod. Phys. 73, 515 (2001). — The definitive DFPT review
  2. X. Gonze and C. Lee, "Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory", Phys. Rev. B 55, 10355 (1997). — Theoretical foundations of DFPT
  3. A. Togo and I. Tanaka, "First principles phonon calculations in materials science", Scr. Mater. 108, 1 (2015). — Phonopy methodology and applications
  4. K. Parlinski, Z. Q. Li, and Y. Kawazoe, "First-Principles Determination of the Soft Mode in Cubic ZrO$_2$", Phys. Rev. Lett. 78, 4063 (1997). — Early frozen phonon methodology
  5. W. Kohn, "Nobel Lecture: Electronic structure of matter—wave functions and density functionals", Rev. Mod. Phys. 71, 1253 (1999). — DFT foundations
  6. D. C. Wallace, Thermodynamics of Crystals, Dover (1998). — Statistical mechanics and phonon thermodynamics
  7. M. T. Dove, Introduction to Lattice Dynamics, Cambridge University Press (1993). — Comprehensive lattice dynamics textbook
  8. G. Grimvall, Thermophysical Properties of Materials, North-Holland (1999). — Connection to experimental thermodynamics
← Series Index Chapter 2: Anharmonic Phonons →

Disclaimer

This educational content was generated with AI assistance for the Hashimoto Lab knowledge base. While efforts have been made to ensure accuracy, readers should verify critical information with primary sources and peer-reviewed literature.