2.16. MedeA MD Phonon: Vibrational Properties from Molecular Dynamics Simulations

Contents

Introduction
Setting up a MedeA MD Phonon Simulation
Assessing the Results
Theoretical Background

download:	`pdf`

2.16.1. Introduction

By applying quantum statistics to the normal vibrational modes of a system thermodynamic properties such as entropy or Helmholtz free energy can be calculated accurately and efficiently. In the MedeA MD Phonon module this is done using velocity auto-correlation functions obtained from molecular dynamics simulations.

Key Benefits

Study systems with low symmetry. For example amorphous systems or crystalline structures with defects.
Automated plot creation facilitates analysis of results.
View the plots for velocity auto-correlation functions of each element in each subset
Study the vibration density of states, its partial contributions of all elements, and all subsets in automatically generated plots on the JobServer and within MedeA’s analysis tool.
Evaluate thermodynamic properties such as internal energy, entropy, Helmholtz free energy, and heat capacity.

Hint

The MedeA MD Phonon module works with MedeA LAMMPS classical molecular dynamics simulations. Ab initio MD trajectories are not supported.

2.16.2. Setting up a MedeA MD Phonon Simulation

The MedeA MD Phonon module can be accessed in a MedeA LAMMPS Flowchart by inserting the MD Phonon stage into any LAMMPS stage. The vibrational density of state is calculated from the velocity auto-correlation function.

Note

Ensure that the structure has been equilibrated prior to the MD Phonon stage. Use appropriately configured NPT and NVT stages for this.

The following illustrates the MedeA MD Phonon stage:

The parameters are:

Time: Duration of a simulation run to calculate the velocity auto-correlation function. The appropriate value for this parameter is, in most cases, between 10 - 100 ps.

Lower values will reduce the resolution of the calculated vibrational density of states to a too low value. Higher values will not further improve the results. The density of states frequency resolution is defined by the reciprocal of the simulation time. A 1 ps simulation yields a resolution of 1 THz. Increasing this time to 10 ps yields a resolution of 0.1 THz.
Number of runs to average over: The specified number of runs with a time length as defined by the variable time is performed. These runs are performed in parallel and the calculated velocity auto-correlation functions of each run are then averaged. Increasing this number will provide smoother results. Note, that the total simulation time is defined by time \(\times\) number of runs.
Time step: Time step size employed in solving the equations of motion in an NVE ensemble.
Sampling: Number of samples employed in performing averaging.
Trajectory: Number of trajectory frames saved during the molecular dynamics calculation.
Gaussian low pass filtering sigma (fs): Define the width of the Gaussian used for the low pass filter in fs. For best results match this value with the decorrelation times of the velocity auto-correlation functions.
Subset: One or more subsets of atoms for which partial vibrational density of states and thermodynamic properties are calculated.

2.16.3. Assessing the Results

The Job.out file provides a brief summary on the calculation setup of the MD Phonon stage as well as the calculated thermodynamic properties at the simulated temperature:

       Property         Value   +/- Uncertainty    Units      After Steps   % Run
------------------  ----------     ----------- ------------- ------------ -------

               Evib:     9589.0                 kJ/mol
            Evib_QC:     3656.7                 kJ/mol
               Avib:     7336.8                 kJ/mol
                ZPE:     8147.3                 kJ/mol
              SvibT:     2252.2                 kJ/mol
               Svib:     7272.8                 J/K/mol
                 Cv:    10751.8                 J/K/mol
            Etot_QC:  -406758.6                 kJ/mol
                       Atotal:  -407226.0                 kJ/mol

   The thermodynamic properties printed above are given for the whole structure
   Al384N384 containing in total 768 atoms.

   Thermodynamic properties have been calculated at a temperature of 309.7 K by
   integrating over the phonon density of states obtained from the
   molecular dynamics simulation.


             Definitions:
                  Evib(T)        : vibrational internal energy (including quantum contributions)
                 Evib_QC(T)      : only the quantum contributions to vibrational internal energy
                 Etot_QC(T)      : sum of MD total energy Etotal, and Evib_QC
                  Avib(T)        : vibrational Helmholtz free energy (including quantum contributions)
                  Atotal(T)      : Etot_QC(T) - Svib(T) T
                    ZPE          : zero point energy (from quantum mechanics)
                   SvibT         : vibrational entropy times temperature (including quantum contributions)
                  Svib(T)        : vibrational entropy (inculding qunatum contributions)
                   Cv(T)         : vibrational heat capacity at constant volume (including quantum contributions)

   If not specified differently the thermodynamic functions contain both
   the classical and quantum contributions.

In addition, the calculated thermodynamic properties are provided per formula unit, per subset (as fraction of formula units), and per element.

The following thermodynamic functions are provided
for the empirical formula AlN.

             Evib   Evib_QC      Avib       ZPE     SvibT     Svib       Cv
          kJ/mol    kJ/mol    kJ/mol    kJ/mol    kJ/mol  J/K/mol  J/K/mol
       ---------- -------------------  --------- -------- -------- ----------
Total:    24.9714    9.5228   19.1064   21.2170    5.8651  18.9396  27.9996
   Al:    11.0427    3.3184    7.0766    8.5695    3.9662  12.8076  16.7912
    N:    13.8358    6.1114   11.8703   12.5162    1.9655   6.3470  11.3883

Besides information contained in the Job.out file graphical output is provided by the files:

VACF<subset>_<element>_avg.png shows the velocity auto-correlation function of a given subset and element/isotope:

phonon_DOS.png contains the total vibrational density of states and, if subsets were considered, the partial vibrational density of states of each subset:

Phonon_DOS_Total.png contains the total vibrational density of states as well as the partial contributions of each element/isotope:

Phonon_DOS_subset<number>.png contains for each subset the partial vibrational density of states as well as the partial contributions of each element/isotope:

2.16.3.1. Viewing Vibrational Density of States in MedeA

By clicking on Analysis >> Phonon Density of States and selecting the results from the relevant completed job the vibrational density of states can be viewed interactively in MedeA. For more information on this refer to the section on viewing the phonon DOS in the MedeA Phonon documentation.

On Windows, the vibrational Density of states can also be exported to an Excel spreadsheet with Analysis >> Export Phonon Dispersions and DOS. The section on exporting DOS in the MedeA Phonon documentation provides more information on this functionality.

2.16.4. Theoretical Background

The vibrational density of state, \(DOS(\nu)\), is defined as the distribution of normal modes of a system [1]. Hence, the vibrational intensity at a given frequency \(\nu\) is calculated by the sum of all \(N\) atomic contributions [2] in the x, y, and z direction, as indicated by the sum over \(k\), at temperature \(T\) with

\[DOS(\nu) = \frac{1} {3N (2\pi)^3 k_B T} \sum_{j=1}^{N} \sum_{k=1}^{3} m_j s_j^k(\nu) ,\]

where \(m_j\) defines the mass of atom \(j\) and \(k_B\) the Boltzmann constant. The spectral density, \(s_j^k(\nu)\), can be described by the Fourier transform of the integrated velocity auto-correlation function, \(<v_j^k(t+\tau)\cdot v_j^k(t)>\). This allows for a rewrite of the equation above to

\[DOS(\nu) = \frac{1} {3N (2\pi)^3 k_B T} \sum_{j=1}^{N} \sum_{k=1}^{3} m_j \int_{-\infty}^{\infty} <v_j^k(t+\tau)\cdot v_j^k(t)> e^{-i \nu \tau} d\tau ,\]

putting the density of state into direct relation with the velocity auto-correlation function. The velocity auto-correlation function is calculated on-the-fly by LAMMPS during a NVE MD simulation.

Once \(DOS(\nu)\) has been established the thermodynamic properties, such as internal energy, entropy, Helmholtz free energy, and heat capacity, can be calculated by appropriate integrals.

Helmholtz free energy:

\[A_{vib} = k_BT \int_{\nu=0}^{\infty} DOS(\nu) \ln \left ( 2 \sinh \left ( \frac{\beta h \nu}{2} \right) \right) d \nu\]

Potential energy:

\[E_{vib} = \frac{1}{2} \int_{\nu=0}^{\infty} DOS(\nu) h \nu \left( \frac{\exp(\beta h \nu)+1}{\exp(\beta h \nu)-1} \right) d \nu\]

Entropy:

\[S_{vib} = \frac{E_{vib} - A_{vib}} {T}\]

Specific heat:

\[C_V = k_B \int_{\nu=0}^{\infty} DOS(\nu) (\beta h \nu)^2 \frac{\exp(\beta h\nu)}{(\exp(\beta h \nu)-1)^2} d \nu\]

In the above equations \(h\) stands for the Planck constant and \(\beta=1/(k_BT)\). These equations contain both the classical and quantum mechanical contributions. The classical potential energy can be calculated simply with \(E^{cl}_{vib}=3k_B T\) and the classical Helmholtz free energy with \(A^{cl}_{vib}=k_B T\int DOS(\nu) \ln\left(\beta h \nu\right ) d\nu\).

[1]	Peter H. Berens, Donald H. J. Mackay, Gary M. White, and Kent R. Wilson, J. Chem. Phys. 79, 2375 (1983)

[2]	Shiang-Tai Lin, Mario Blanco, and William A. Goddard III, J. Chem. Phys. 119, 11792 (2003)

download:	`pdf`