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$ denote a thermal average which has to be computed over different molecular dynamics 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 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-NVT simulation

#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 as follows

Equilibrate.sh

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

This bash script will create POSCAR_i where $i$ runs from 1 to 10. The number of time steps between successive starting structures should be determined according to the characteristics of the system under investigation. Usually equilibrating the system of interest for one or two times the self-correlation time is a suitable choice. 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 for each POSCAR-file. An INCAR file for NVE simulations can look as follows:

INCAR NVE simulations

#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 bash 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, a set of NVE trajectories is 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 vaspout.h5-file. Those are needed to compute the normalized velocity auto correlation functions.

Analyzing data

Step 3 to 4 Computing phonon density of states and checking for convergence

For the following time signal analysis and computation of the phonon density of states the following script are needed

The following python script can be used to compute normalized velocity auto correlation functions

Click to show ComputeCorrelation.py

import numpy as np
class AutoCorrelation:
    """
    A class to compute the velocity auto-correlation function for a given set of velocity data.

    Attributes:
    -----------
    delta : int, optional
        The step size for time intervals in the computation (default is 1).

    Methods:
    --------
    velocity_auto_correlation(velos):
        Computes the velocity auto-correlation function for the input velocity data.
    """
    def __init__( self, delta = 1 ):
        """
        Initializes the AutoCorrelation object with a specified time step size.

        Parameters:
        -----------
        delta : int, optional
            The step size for time intervals in the computation (default is 1).
        """
        self.delta = delta
    def velocity_auto_correlation( self, velos ):
        """
        Computes the velocity auto-correlation function for the given velocity data.

        Parameters:
        -----------
        velos : numpy.ndarray
            A 3D array of shape (Nt, Nx, Ndim) representing the velocity data, where:
            - Nt is the number of time steps,
            - Nx is the number of particles,
            - Ndim is the number of spatial dimensions.

        Returns:
        --------
        numpy.ndarray
            A 2D array of shape (Nt // 2, Nx) representing the velocity auto-correlation function
            for each particle over time.

        Notes:
        ------
        - The function normalizes the correlation values using the squared norm of the initial velocities.
        - The computation is performed for time intervals up to Nt // 2.
        """
        Nt, Nx, Ndim = velos.shape
        deltaT = self.delta
        corr_func = np.zeros( [ Nt // 2, Nx ] )
        counter   = np.zeros( [ Nt // 2, 1 ] )
        for dt in range( 0, Nt//2, deltaT ):
            v0   = velos[ dt, :, : ]
            norm = np.asarray( [ np.linalg.norm( v0[ i, : ] )**2 for i in range( Nx ) ] )
            for t in range( dt, Nt//2 ):
                vt = velos[ t, :, : ]
                value = np.asarray( [ np.dot( vt[i,:], v0[ i, : ] ) for i in range( Nx ) ] )
                corr_func[ t-dt, : ] += value / norm
                counter[ t-dt ] += 1
        return corr_func / counter

The following python script can be used to the phonon density of states by computing the power spectra of the normalized velocity auto correlation functions.

Click to show PhononDOS.py

import sys
import py4vasp
import numpy as np
import matplotlib.pyplot as plt


import ComputeCorrelation
    
class ComputePhonons:
    """
    @brief Class to compute phonon-related properties such as autocorrelation, power spectra, and averages.
    
    This class provides methods to compute velocity autocorrelation, power spectra, and averages for atomic systems 
    based on velocity data. It also includes functionality to write the computed data to files.
    
    @class ComputePhonons
    """
    def __init__( self, fname, dt = 1.0, timeShift=50 ):
        """
        @brief Constructor to initialize the ComputePhonons object.
        @param fname Path to the input file for the calculation.
        @param dt Time step in femtoseconds (default: 1.0).
        @param timeShift Time shift for autocorrelation computation (default: 50).
        """
        self.fname  =  fname
        self.calc   =  py4vasp.Calculation.from_path( self.fname )
        self.velos  =  self.calc.velocity[:].read()
        self.time_step =  dt /1000 # thz output
        self.timeShift = timeShift

    def compute_ac( self ):
        """
        @brief Compute the velocity autocorrelation function.
        This method calculates the velocity autocorrelation function using the provided velocity data.
        """
        dos     =  ComputeCorrelation.AutoCorrelation( self.timeShift )
        self.ac =  dos.velocity_auto_correlation( self.velos["velocities"] )

    def compute_averages( self ):
        """
        @brief Compute averages of the autocorrelation function for total and per-atom contributions.
        This method calculates the total autocorrelation and groups the autocorrelation by atomic species.
        """
        unique, counts = np.unique_counts( self.velos["structure"]["elements"] )
        self.total_ac = np.sum( self.ac, axis=1 )
        labels = self.velos["structure"]["elements"]
        unique_labels, inverse = np.unique(labels, return_inverse=True)
        result = np.zeros((self.ac.shape[0], len(unique_labels)), dtype=self.ac.dtype)
        np.add.at(result, (slice(None), inverse), self.ac )
        self.atom_ac = {label: result[:, i] for i, label in enumerate(unique_labels)}

    def compute_power_spectra( self ):
        """
        @brief Compute the power spectra for total and per-atom contributions.
        This method calculates the power spectra using the Fourier transform of the autocorrelation functions.
        """
        self.ps_total = np.abs( np.fft.fft( self.total_ac ) )**2
        self.ps_atom  = {}
        for key in self.atom_ac.keys():
            self.ps_atom[key] = np.abs( np.fft.fft( self.atom_ac[key] ) )**2
        
        freqs = np.fft.fftfreq( self.ps_total.shape[0], self.time_step )
        self.ps_total = np.vstack( [freqs, self.ps_total/np.max(self.ps_total)] ).T
        self.ps_total = self.ps_total[ :self.ps_total.shape[0]//2, : ]
        for key in self.ps_atom.keys():
            self.ps_atom[key] = np.vstack( [freqs, self.ps_atom[key]/np.max( self.ps_atom[key] )] ).T
            self.ps_atom[key] = self.ps_atom[key][ :self.ps_atom[key].shape[0]//2, : ] 

    def write_total_ps( self, fname="total_ps.dat" ):
        """
        @brief Write the total power spectrum to a file.
        @param fname Name of the output file (default: "total_ps.dat").
        """
        np.savetxt( fname, self.ps_total )

    def write_total_ac( self, fname="total_ac.dat" ):
        """
        @brief Write the total autocorrelation function to a file.
        @param fname Name of the output file (default: "total_ac.dat").
        """
        x = np.linspace( 0, self.time_step*self.total_ac.shape[0], self.total_ac.shape[0] )
        result = np.vstack( [x, self.total_ac] ).T
        np.savetxt( fname, result )
    
    def write_atom_ac( self ):
        """
        @brief Write the per-atom autocorrelation functions to files.
        Each atomic species' autocorrelation function is written to a separate file.
        """
        for key in self.ps_atom.keys():
            np.savetxt( f"{key}_ps.dat", self.ps_atom[key] )
 
     def write_atom_ps( self ):
         """
         @brief Write the per-atom power spectra to files.
         Each atomic species' power spectrum is written to a separate file.
         """
         for key in self.ps_atom.keys():
             np.savetxt( f"{key}_ps.dat", self.ps_atom[key] )
 
 
 if __name__=="__main__":
     x = ComputePhonons( sys.argv[1], float(sys.argv[2]) )
     x.compute_ac()
     x.compute_averages()
     x.compute_power_spectra()
     x.write_total_ps()
     x.write_total_ac()
     x.write_atom_ps()
     x.write_atom_ac()

The PhononDOS.py script can be used to compute the phonon density of states for a given NVE simulation folder containing an vaspout.h5 file created with the before mentioned INCAR NVE simulation. The script will create a file called total_ps.dat containing the total phonon density of states. The partial phonon density of states of the atomic species are written to files ElementKey_ps.dat. As input the script needs a folder name containing a vaspout.h5-file and the second input argument has to be the simulation time step of your simulation in fs. The written files will contain the frequency in THz as the first column. The second column will contain the phonon density of states computed from the power spectrum of the velocity auto correlation function.

Step 3 to 4 Computing phonon density of states and checking for convergence

When sampling phonons from molecular dynamics simulations a careful convergence analysis has to be done. To check results for convergence, a convergence plot can be made. To make a convergence analysis plot, the results obtained for the NVE trajectories can be averaged. To do so a single trajectory can be plotted, compared to an average over 2 trajectories and so on. The so obtained lines can be collected in a single plot as shown in Fig. 1. The yellow line shows an average over a single trajectory. The more red the lines are the more trajectories have been used for computing the average. The red line shows the average computed over all trajectories computed.

Fig. 1: $\mathbf{Left:}$ Shows convergence analysis of normalized velocity correlation function. $\mathbf{Right:}$Convergence analysis of phonon density of states.

From the plot it is possible to conclude that enough data was obtained to properly converge the phonon density of states. Additionally to the total phonon density of states the atom resolved normalized auto correlations and phonon density of states can be obtained. These plots are shown in Fig. 2.

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