Sampling phonon spectra from molecular-dynamics simulations

Sampling phonon DOS from molecular dynamics simulation

The phonon density of states can be obtained as the power spectrum from the normalized velocity auto correlation function. The normalized velocity auto correlation function for a $N$-particle system is given by \begin{equation} f(t)=\sum_{s=1}^{types}f_{s}(t)=\frac{\langle \sum_{s=1}^{types}\sum_{i=1}^{N_{s}}\mathbf{v}_{i}(\Delta T)\mathbf{v}_{i}(\Delta T+t) \rangle}{\mathbf{v}_{i}(\Delta T)\mathbf{v}_{i}(\Delta T)}. \end{equation} The brackets $\langle ,\rangle$ denotes a thermal average which has to be computed over different trajectories and and starting times $\Delta T$ within each trajectory. The sum over $i$ runs over the atoms within each species and the sum $s$ is over all atomic species contained in the simulated system. From this the phonon density of states is obtained by computing the power spectrum of $f_{s}(t)$: \begin{equation} g(\omega)=g_{s}(\omega)=\sum_{s=1}^{types}g_{s}(\omega)=\left| \sum_{s=1}^{types}\int_{-\infty}^{\infty}f_{s}(t)e^{-i\omega t}\right|^{2}. \end{equation} To properly sample the phonon density of states from molecular dynamics simulations the following the steps have to be accomplished:

code	publication
step 1	generate starting structures in NVT simulation
step 2	generate a NVE simulation for every starting structure to obtain velocity fields
step 3	compute normalized velocity auto correlation function for every NVE simulation
step 4	compute power spectrum for every normalized velocity auto correlation function
step 5	compute averages and check for convergence

Creating trajectories

Step 1: Generating starting configurations from the NVT ensemble

The phonon density of states has to be sampled from NVE simulations. NVE simulations have to be used because otherwise the thermostat would add perturbations to the atomic velocities. To generate starting configurations for the NVE simulations NVT simulations with the Langevin thermostat are done. The Langevin thermostat is used for the equilibartion of our starting structures because it is a stochastic thermostat and therefore will be suitable to populate the available phonon modes of our system uniformly. The uniform population of the phonon modes originates from the Langevin thermostat's property to add white noise onto the velocity auto correlation function. A simple INCAR file which will perform an NVT simulation could look as follows

INCAR

#INCAR molecular-dynamics tags NVE ensemble 
IBRION = 0                   # choose molecular-dynamics 
MDALGO = 1                   # using Andersen thermostat
ISIF = 0                     # compute stress tensor but do not change box volume/shape 
TEBEG = 300                  # set temperature 
NSW = 10000                  # number of time steps 
POTIM = 2.0                  # time step in femto seconds 
LANGEVIN_GAMMA = 0.5 0.5 0.5 $ Langevine friction coefficient for 3 atomic species.

A bash script to produce 10 starting configurations in the from of POSCAR files could look ass follows

Equilibrate.sh

#Equi.sh script to generate POSCAR_1 to POSCAR_10 
for i in {1..20}; do
   cp POSCAR POSCAR_$i
   mpirun -np 32 vasp_std
   cp CONTCAR CONTCAR_$i
   cp CONTCAR POSCAR
done

Related tags and articles

Molecular-dynamics calculations, Computing the phonon dispersion and DOS, IBRION, MDALGO, ISIF, TEBEG, NSW, POTIM, ANDERSEN_PROB, QPOINTS, LPHON_DISPERSION, PHON_NWRITE, LPHON_POLAR, PHON_DIELECTRIC, PHON_BORN_CHARGES,PHON_G_CUTOFF