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:

step	task
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 equilibration 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 functions. 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                     # don't compute stress tensor. Box shape has to be fix
TEBEG = 500                  # 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

The number of time steps that has to be taken between taking starting structures has to be decided based on the system of interest. Also the time step has to be chosen according to the system. Ideally the time step is chosen such that the frequency of the fastest phonon mode of interest can still be resolved.

Step 2: Generating NVE trajectories to obtain atom velocities

After obtaining a bunch of starting structures stored in files POSCAR_1 to POSCAR_10 (including velocities) NVE simulations have to be performed. An INCAR file for NVE simulations can look as follows:

INCAR

#INCAR molecular-dynamics tags NVE ensemble 
IBRION = 0                   # choose molecular-dynamics 
MDALGO = 1                   # using Andersen thermostat
ISIF = 0                     # don't compute stress tensor. Box shape has to be fix 
TEBEG = 500                  # set temperature 
VELOCITY = T                 # make sure to write velocities to vaspout.h5
NSW = 10000                  # number of time steps 
POTIM = 2.0                  # time step in femto seconds 
ANDERSEN_PROB = 0.0          # setting Andersen collision probability to zero to get NVE enseble

Again it is advisable to use a script to generate NVE trajectories. The following script will assume a base folder containing POSCAR files named POSCAR_1 to POSCAR_10, an INCAR file, a KPOINTS file and an POTCAR file. The script will create folders Run1 to Run10. Each folder will contain a vaspout.h5 file after script execution. These vaspout.h5 files will be needed for the analysis scripts of the next section.

Production.sh

#Run NVE MD simulation for every starting configuration
for i in {1..10}; do
   mkdir Run$i
   cd Run$i
   cp ../INCAR .
   cp ../KPOINTS .
   cp ../POSCAR_${i} POSCAR
   vasp_std
   cd ..
done

After finishing Step 1 to Step 2 to a set of NVE trajectories was obtained. The trajectories are stored to the vaspout.h5 file. By setting VELOCITY=T it was assured that the atomic velocities were written to the file. Those are needed to compute the normalized velocity correlation functions.

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