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

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:

../../_images/image0018.png

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:    45132.2                 kJ/mol
           Evib_QC:    17621.5                 kJ/mol
              Avib:    34837.2                 kJ/mol
           Avib_QC:     9608.6                 kJ/mol
               ZPE:    38554.9                 kJ/mol
             SvibT:    10295.1                 kJ/mol
              Svib:    34499.7                 J/K/mol
           Svib_QC:    26852.0                 J/K/mol
                Cv:    50859.3                 J/K/mol

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

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

Definitions:
     Evib(T)        : vibrational internal energy
    Evib_QC(T)      : quantum contributions to vibrational internal energy
     Avib(T)        : vibrational Helmholtz free energy
    Avib_QC(T)      : quantum contributions to vibrational Helmholtz free energy
       ZPE          : zero point energy
      SvibT         : vibrational entropy times temperature
     Svib(T)        : vibrational entropy
    Svib_QC(T)      : quantum contributions to vibrational entropy
      Cv(T)         : vibrational heat capacity at constant volume

 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 Avib_QC     ZPE   SvibT    Svib Svib_QC      Cv
        kJ/mol  kJ/mol  kJ/mol  kJ/mol  kJ/mol  kJ/mol J/K/mol J/K/mol J/K/mol
      -------- ------- ------- ------- ------- ------- ------- ------- -------
Total: 24.4222  9.5354 18.8513  5.1994 20.8631  5.5709 18.6687 14.5303 27.5213
   Al: 10.8269  3.3835  7.1789  1.8176  8.5279  3.6480 12.2246  5.2476 16.4391
    N: 13.5930  6.1496 11.6685  3.3806 12.3319  1.9244  6.4489  9.2793 11.0867

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:
../../_images/VACF.png

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

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

  • 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:
../../_images/phonon_DOS_subset1.png

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.

../../_images/GUI_phonon_DOS.png

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.

../../_images/EXCEL_phonon_DOS.png

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