Calculate phonon spectra by DFT

Phonons are quasiparticles representing the collective excitations of a crystal structure, describing the quantized modes of lattice vibrations that follow Bose-Einstein statistics. They are closely related to the thermodynamic properties of crystals. Therefore, calculating the phonon dispersion relations, or phonon spectra, is of great significance. This process involves obtaining the necessary force constants from first-principles (FP) calculations, followed by further processing and computations using tools like Phonopy. Two common FP methods for this purpose are the finite displacement method (FDM) and density functional perturbation theory (DFPT).

flowchart LR
A(Structure optimization) --> B(Supercell)
B --> D[FDM by VASP]
B --> E[DFPT by VASP]
D --> F[Force constants]
E --> F[Force constants]
F --> G[Phonon spectra]

Here I’m gonna use monolayer molybdenum disulfide (MoS₂) as an example to elaborate on the process:

1. Structure optimization

Perform high-precision structure optimization of the initial unit cell, which can be from some tools like ASE , optimizing both lattice constants and atomic positions.

INCAR tag [ISIF]:
Keep lattice constants, and only optimize atomic positions;
The lattice constants, atomic positions, and volume are all optimized;
The lattice constants and atomic positions are optimized, but the volume does not.

WARNING

Using ISIF=3 for relaxation, which is unsuitable for 2D systems with vacuum as it may collapse the cell. Recommend using ISIF=4 to keep cell volume constant while relaxing lateral lattice parameters.

Global Parameters
ISTART =  0            # Read existing wavefunction, if there
ISPIN  =  1            # Non-Spin polarised DFT
ICHARG =  2            # The initial guess charge density is generated by superposition of atomic charge density

LREAL  =  Auto         # Projection operators: False, important for phonon dispersion

ENCUT  =  650          # Cut-off energy for plane wave basis set, in eV
PREC   =  Accurate     # Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation
LWAVE  = .FALSE.       # Write WAVECAR or not
LCHARG = .FALSE.       # Write CHGCAR or not
ADDGRID= .TRUE.        # Increase grid, helps GGA convergence

Electronic Relaxation
ISMEAR =  -5           # Gaussian smearing, metals:1
SIGMA  =  0.05         # SIGMA is ignored for the tetrahedron method when ISMEAR = -5
ALGO   =  N            # To specify the electronic minimisation algorithm
NELM   =  120          # Max electronic SCF steps
NELMIN =  6            # Min electronic SCF steps
EDIFF  =  1E-06        # SCF energy convergence, in eV
# GGA  =  PS           # PBEsol exchange-correlation

Ionic Relaxation
NSW    =  120          # Max ionic steps
IBRION =  2            # Algorithm: 0-MD, 1-Quasi-New, 2-CG
ISIF   =  4            # Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions
EDIFFG =  -0.02        # Ionic convergence, eV/AA
ISYM   =  -1           # Symmetry: -1=switched off, 0=none, 2=GGA, 3=hybrids

Global Parameters
ISTART =  0            # Read existing wavefunction, if there
ISPIN  =  1            # Non-Spin polarised DFT
ICHARG =  2            # The initial guess charge density is generated by superposition of atomic charge density

LREAL  =  Auto         # Projection operators: False, important for phonon dispersion

ENCUT  =  650          # Cut-off energy for plane wave basis set, in eV
PREC   =  Accurate     # Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation
LWAVE  = .FALSE.       # Write WAVECAR or not
LCHARG = .FALSE.       # Write CHGCAR or not
ADDGRID= .TRUE.        # Increase grid, helps GGA convergence

Electronic Relaxation
ISMEAR =  -5           # Gaussian smearing, metals:1
SIGMA  =  0.05         # SIGMA is ignored for the tetrahedron method when ISMEAR = -5
ALGO   =  N            # To specify the electronic minimisation algorithm
NELM   =  120          # Max electronic SCF steps
NELMIN =  6            # Min electronic SCF steps
EDIFF  =  1E-08        # SCF energy convergence, in eV
# GGA  =  PS           # PBEsol exchange-correlation

Ionic Relaxation
NSW    =  120          # Max ionic steps
IBRION =  2            # Algorithm: 0-MD, 1-Quasi-New, 2-CG
ISIF   =  2            # Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions
EDIFFG =  -1E-04       # Ionic convergence, eV/AA
ISYM   =  -1           # Symmetry: -1=switched off, 0=none, 2=GGA, 3=hybrids

Here is the structure file of MoS₂ with the POSCAR format:

Mo  S
 1.0000000000000000
     3.1600000000000001    0.0000000000000000    0.0000000000000000
    -1.5800000000000001    2.7366402759588260    0.0000000000000000
     0.0000000000000000   -0.0000000000000000   17.2500000000000000
 Mo  S
 1   2
Cartesian
  0.0000000000000000  0.0000000000000000  8.6250000000000000
  1.5799999999999998  0.9122134253196086 10.2500000000000000
  1.5799999999999998  0.9122134253196086  7.0000000000000000

Here is it after relaxation (CONTCAR file):

Mo  S
  1.000000000000000
     3.1857089607817932    0.0000000218282313    0.0000000024895515
    -1.5928544581578377    2.7589049024719241    0.0000000011539288
     0.0000000310083557    0.0000000002516627   16.9727051455909148
   Mo   S
   1    2
Direct
 -0.0000000007065635 -0.0000000020440497  0.4999999968674809
  0.6666666662068135  0.3333333335809116  0.5920061449599645
  0.6666666678330786  0.3333333351298096  0.4079938581725475

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

After completing the full optimization, the CONTCAR file is obtained, allowing us to proceed to the next step.

2. Build supercell

The Phonopy package provides a convenient way for effortlessly expanding the cell.

phonopy -d --dim='8 8 1’ -c CONTCAR

where the CONTCAR file is the structure after optimization.

After that, we will get the following files:

POSCAR-{num}: Displacement files representing structures with slight atomic shifts, numbered by symmetry modes, which will be used to calculate forces for force constants. (for FDM)
SPOSCAR: Abbreviation for “Supercell POSCAR”. (for DFPT)
phonopy_disp.yaml: Displacement information file, including direction, amplitude, etc.

The following command is batch processing for FDM calculation preparation:

for i in 001 002 003; do mkdir -p $i && mv POSCAR-$i $i/POSCAR; done

3. Run VASP

FDM: Perform single-point calculations in VASP using each displaced POSCAR-{num} file.

DFPT: Perform phonon calculation in VASP using SPOSCAR file.

The POTCAR and KPOINTS files can be generated by the VASPKIT tool.

WARNING

You may need to adjust some parameters based on your tests and experience! Note: DFPT requires a significant amount of memory.

# Global Parameters
ISTART =  0            # Not Read existing wavefunction; if there
ISPIN  =  1            # Non-Spin polarised DFT, Default=1

ICHARG =  2            # The initial guess charge density is generated by superposition of atomic charge density
LREAL  = .FALSE.       # Projection operators: automatic
ENCUT  =  650          # Cut-off energy for plane wave basis set, in eV

LWAVE  = .FALSE.       # Write WAVECAR or not
LCHARG = .FALSE.       # Write CHGCAR or not

PREC     =  Accurate   # Precision level
KGAMMA   = .TRUE.      # GAMMA point
KSPACING = 0.25        # Automatically calculate k-points
ALGO     = Normal

# Static Calculation
NSW    =   1           # Default: NSW = 0
IBRION =  -1           # Ions are not moved
ISMEAR =  -5           # gaussian smearing method
SIGMA  =  0.02         # SIGMA is ignored for the tetrahedron method when ISMEAR = -5

EDIFF  =  1E-08        # SCF energy convergence; in eV

NELM   =  150          # Max electronic SCF steps

KPAR   =  4            # make 4 groups, each group working on one set of k-points
NCORE  =  8            # one orbital handled by 4 cores

# For phonon calculation
PREC    =  Accurate   # Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation

ENCUT   =  750        # Cut-off energy for plane wave basis set, in eV
NSW     =  1

IBRION  =  8          # IBRION=8: second derivatives, Hessian matrix, and phonon frequencies (perturbation theory use symmetry)
EDIFF   =  1.0e-8     # SCF energy convergence, in eV
EDIFFG  =  -1E-03     # Ionic convergence; negetive -- force (eV/AA)

IALGO   =  38         # IALGO=38: Blocked-Davidson algorithm (ALGO=N)
ISMEAR  =  0          # Gaussian smearing, metals:1
SIGMA   =  0.02       # Smearing value in eV, metals:0.2)

LREAL   =  .FALSE.    # Projection operators are evaluated in (F:reciprocal space) (T:real space(not recommended)
ADDGRID =  .TRUE.     # Increase grid; helps GGA convergence
LWAVE   =  .FALSE.    # Write WAVECAR file or not
LCHARG  =  .FALSE.    # Write CHGCAR file or not

4. Get force constants info.

After VASP run, you can get the force constants file by :

For FDM, get FORCE_SETS, which contains force and displacement information for perturbed structures:

phonopy -f 001/vasprun.xml 002/vasprun.xml 003/vasprun.xml

For DFPT, get FORCE_CONSTANTS:

phonopy --fc vasprun.xml

5. Get phonon properties

Prepare band.conf, and get the phonon spectra, also phonon density of states (DOS).

ATOM_NAME = Mo  S
DIM = 8  8  1
BAND= 0.0 0.0 0.0  0.5 0.0 0.0  0.3333 0.3333 0.0  0.0 0.0 0.0
FORCE_CONSTANTS = READ

EIGENVECTORS = .TRUE.
GROUP_VELOCITY = .TRUE.
BAND_POINTS = 2001
BAND_LABELS = $\Gamma$ M K $\Gamma$

MP = 30 30 1
DOS = .TRUE.
SIGMA = 0.1

Run the following two commands: