Background Theory#

Mean square displacement#

Mean square displacement (MSD) is used to characterize moving behavior of particle or molecule of interest. Originally, it is used for random motion (Brownian motion) and it is extended to obtain self-diffusivity of the particle or molecule. MSD is measured by the sqaure of displacement of a particle from its reference positions. MSD is calculated by the equation

\[MSD_{\alpha}(\tau) = \Biggl\langle \frac{1}{N^{\alpha}} \sum_{i=1}^{N_{\alpha}} |\mathbf{r}_{\alpha} (\tau) - \mathbf{r}_{\alpha}(t_0)|^2 \Biggr\rangle _ {t_0}\]

where \(r(\tau)\) is the position of particle \(\alpha\) at time \(\tau\), \(r(t_0)\) is the position of particle \(\alpha\) at time origin \(\tau=t_0\), and \(N_{\alpha}\) is the number of equivalent particles in the summation. Here, the everage \(\langle \cdot \cdot \cdot \rangle _ {t_0}\) represents the average over different time origin \(\tau=t_0\).

For self-diffusivity, time derivative of the MSD give the diffusion coeficient as expressed in the equation below.

\[D = \frac{1}{2d}\frac{d}{d\tau}MSD(\tau)\]

where \(d\) is a dimension of the vectors in MSD.

Note

  • It is recommended to average over independent replicate trajectories. One can perform production run starting from different initial configurations or with different random seeds to initialize velocities.

  • For molecule cases, it is recommended using the center of mass (COM) instead of the individual atomic coordinates to avoid the short-time vibrational and rotational displacements of the atom group (molecule).

  • If the system have perioidicy, the unwraped coordinates must be used. This program will unwraped automatically.

Caution

  • Barostats in NPT emsemble will greatly affect the dynamics of a system by altering positions through the volume changes.

  • The interval between time origin for averaging should be longer than the correlation time to gaurantee the independence of the windows separated by an interval.

  • The length of the sampling windows (\(\tau=0 \sim \tau_{max}\)) must be long enough to capture the diffusive regime of the molecule or atom of interest.

  • For best practice, reference [1] is recommended.



Intermittent bond autocorrelation#

To evaluate lifetime of bond of interest, the intermittent bond autocorrelation function is used, following the hydrogen bond lifetime analysis proposed by Rapaport [2]. The presence of a bond is defined by a function \(h(t)\). When bond is formed within the pre-defined range, \(h(t)=1\). If bond breaks and bond length becomes longer than criteria, then \(h(t)=0\). Therefore, the time autocorrelation function of the step function \(h(t)\) can be represented by the following equation.

\[C_I(\tau)=\Biggl\langle \frac{\sum h_{ij}(t_0) h_{ij}(t_0+\tau) }{\sum h_{ij}(t_0)^2} \Biggr\rangle _ {t_0}\]

Here, the everage \(\langle \cdot \cdot \cdot \rangle _ {t_0}\) represents the average over different time origin \(\tau=t_0\). The denominator \(h_{ij}(t_0)^2\) is considered to normalize the value by the initial correlation, ensuring autocorrelation funciton to start decaying from 1.

Note

  • The intermittent scheme is independent of time step between neighboring MD snapshots as it is measuring the correlation between \(\tau=0\) and \(\tau=\tau\) without considering the event in between.

  • In intermittent scheme, fast vibration and oscillation of atoms that cause the temporary bond breaking and forming will not be captured. Instead, intermittent scheme measures the formation of the bond until the its last breaking.

Caution

  • For continuous scheme, accounting for every small the formation and breaking event, new implementation is required.



Angle autocorrelation#

In this implementation, 5 coordinates are required. \(\vec{R}\) is a center where two planes cross. Plane \(A\) is defined by three points: \(\vec{R}\) , \(\vec{A}_1\) and \(\vec{A}_2\). Similarly, plane \(B\) is also defined by three points: \(\vec{R}\) , \(\vec{B}_1\) and \(\vec{B}_2\). To measure the degree of distortion at the metal node, the geometry index defined for four- and five-coordinate complexes were prepared. For the measuring of the angle between two planes centered at metal node, 2 schemes were implemented.

  • Geometry index \(\tau_5\) for five-coordinate complexes [3] [5]

    \[\tau_5 = \frac{\alpha - \beta}{60}\]

    where \(\alpha\) and \(\beta\) are the two greatest angles at the coordination center. Two limits connect between square pyramid \(\tau=0\) and trigonal bipyramid \(\tau=1\). Here, \(\alpha\) and \(\beta\) can be obtain by

    \[\alpha= \cos^{-1} \frac{ (\vec{A}_1 - \vec{R}) \cdot (\vec{B}_1 - \vec{R}) } { \left\lVert \vec{A}_1 - \vec{R} \right\rVert \left\lVert \vec{B}_1 - \vec{R} \right\rVert } ~~\mathrm{and}~~ \beta= \cos^{-1} \frac{ (\vec{A}_2- \vec{R}) \cdot (\vec{B}_2 - \vec{R}) } { \left\lVert \vec{A}_2 - \vec{R} \right\rVert \left\lVert \vec{B}_2 - \vec{R} \right\rVert }.\]

  • Geometry index \(\tau_{\delta}\) for four-coordinate complexes [4] [5]

    \[\tau_{\delta} = \frac{360 - (\alpha + \beta)}{141} \cdot \frac{\alpha}{\beta}\]

    where \(\alpha\) and \(\beta\) are the two greatest angles at the coordination center. Two limits connect between square planar \(\tau=0\) and tetahedral \(\tau=1\). Here, \(\alpha\) and \(\beta\) can be obtain by

    \[\alpha= \cos^{-1} \frac{ (\vec{A}_1 - \vec{R}) \cdot (\vec{B}_1 - \vec{R}) } { \left\lVert \vec{A}_1 - \vec{R} \right\rVert \left\lVert \vec{B}_1 - \vec{R} \right\rVert } ~~\mathrm{and}~~ \beta= \cos^{-1} \frac{ (\vec{A}_2- \vec{R}) \cdot (\vec{B}_2 - \vec{R}) } { \left\lVert \vec{A}_2 - \vec{R} \right\rVert \left\lVert \vec{B}_2 - \vec{R} \right\rVert }.\]

  • Angle between plane defined by normal vectors

    \[\theta = \cos^{-1} \frac{\vec{n}_A \cdot \vec{n}_B}{ \Vert \vec{n}_A \Vert \Vert \vec{n}_B \Vert }\]

    where normar vectors \(\vec{n}_1\) and \(\vec{n}_2\) are defined by

    \[\vec{n}_A = (\vec{A}_1 - \vec{R}) \times (\vec{A}_2 - \vec{R}) ~~\mathrm{and}~~ \vec{n}_B = (\vec{B}_1 - \vec{R}) \times (\vec{B}_2 - \vec{R})\]

  • Angle between plane defined by parallel vectors

    \[\varepsilon = \cos^{-1} \frac{\vec{p}_A \cdot \vec{p}_B}{ \Vert \vec{p}_A \Vert \Vert \vec{p}_B \Vert }\]

    where parallel vectors \(\vec{p}_1\) and \(\vec{p}_2\) are defined by

    \[\vec{p}_A = \frac{1}{2} \times \Bigl[(\vec{A}_1 - \vec{R}) + (\vec{A}_2 - \vec{R}) \Bigr] ~~\mathrm{and}~~ \vec{p}_B = \frac{1}{2} \times \Bigl[(\vec{B}_1 - \vec{R}) + (\vec{B}_2 - \vec{R}) \Bigr]\]

For the persistance of the geometry index \(\tau_5\) and \(\tau_{\delta}\) as time propagated, the autocorrelation function is defined by measuring the correlation between themself at time \(\tau=t_0\) and \(\tau=t_0+\xi\).

\[C(\xi)=\Biggl\langle \frac{\sum \tau_{}(t_0) \tau(t_0+\xi) }{\sum \tau(t_0)^2} \Biggr\rangle _ {t_0}\]

For the persistance of the angle \(\theta\) and \(\varepsilon\) as time propagated, the autocorrelation function is defined in a same manner with the dihedral autocorrelation function [6].

\[C(\xi)= \biggl\langle \cos \bigl[ \theta (t_0) - \theta (t_0 + \xi) \bigr] \biggl\rangle _ {t_0}\]

Although this is not a product of two functions separated by time as used in a general autocorrelation functions, the above equation can be written in a form of the sum of two products. \(C(\xi)\) gives a value of 1 when \(\theta (t_0) = \theta (t_0 + \xi)\), while gives a value of 0 when \(\theta (t_0)\) and \(\theta (t_0 + \xi)\) are crossed, essentially representing the measurement of autocorrelation function.

\[C(\xi)= \biggl\langle \cos \bigl(\theta (t_0)\bigr) \cos \bigl(\theta (t_0 + \xi)\bigr) + \sin \bigl(\theta (t_0)\bigr) \sin \bigl(\theta (t_0 + \xi)\bigr) \biggl\rangle _ {t_0}\]

Caution

  • Be careful when assign the atomic index in input file.

  • The definition of two greatest angles in \(\tau_5\) and \(\tau_{\delta}\) is not general. Please check again the source code if the target system is changed.



Dipole-derivative autocorrelation#

In order to calculate infrared (IR) absorption spectrum, the dipole moment autocorrelation function is required. [7]

\[\Gamma_M (\tau) = \Biggl\langle \vec{M}(t_0) \cdot \vec{M}(t_0 + \tau) \Biggr\rangle = \Biggl\langle \sum_{j=1}^n e_j \vec{r}_j(t_0) \cdot \sum_{j=1}^n e_j \vec{r}_j(t_0 + \tau) \Biggr\rangle _ {t_0}\]

Here \(\vec{M}(\tau)\) is the sum of all the dipole moments of the atoms in the simulation box at time \(\tau\).

\[d\vec{M}(\tau) = \sum _ {j=1} ^ N e_j \vec{r}_j(\tau)\]

Under periodic boundary conditions, the wrapped atomic coordinates can give discontinuities when atom moves beyond cell boundary. To avoid this problem, distances between atoms will be used instead of the coordinates. Thus, the dipole moment time derivative autocorrelation function can be computed using the following equation. [8]

\[\Gamma_M ' (\tau) = \Biggl\langle \frac{d\vec{M}(t_0) }{d\tau} \cdot \frac{d\vec{M}(t_0 + \tau) }{d\tau} \Biggr\rangle = \Biggl\langle \sum_{j=1}^n e_j\frac{d\vec{r}_j(t_0) }{d\tau} \cdot \sum_{j=1}^n e_j \frac{d\vec{r}_j(t_0 + \tau) }{d\tau} \Biggr\rangle _ {t_0}\]

Here \(\vec{M}(\tau)\) is the sum of all the dipole moments of the atoms in the simulation box at time \(\tau\). The derivaitve of the dipole moment at atoms \(j=1,N\) can be approximated assumming constant charges \(e_j\) using the following relation.

\[\frac{d\vec{M}(\tau) }{d\tau} = \sum _ {j=1} ^ N e_j \frac{d\vec{r}_j(\tau) }{d\tau} = \sum _ {j=1} ^ N e_j \vec{v}_j (\tau)\]

Then, the IR spectra can be obtain from Fourier transform of the \(\Gamma_M (\tau)\) and \(\Gamma_M ' (\tau)\).

\[I(\omega) \propto \int _0^{\infty} \Gamma_M (\tau) \cos (\omega \tau ) d \tau\]

or

\[I(\omega) \propto \int _0^{\infty} \Gamma_M ' (\tau) \cos (\omega \tau ) d \tau\]

The two scheme, \(\Gamma_M (\tau)\) based and \(\Gamma_M ' (\tau)\) based scheme, are different by the intensity. [7] [8]

Note

  • Time autocorrelation functions should decay to zero rapidly in strongly coupled molecular systems with a large number of atoms. Therefore, the selection of time step in MD sampling is crucial to capture the vibration periods. [8]

Caution

  • Mind that the constant charges are used. Therefore, one should use with care with system in which electronic polarizability and the electron transfer are non-negligible.

  • Autocorrelation function is sensitive to time step in MD simulation.



Velocity autocorrelation#

In order to calculate phonon density of states (DOS), the velocity autocorrelation function is used based on the trajectory from molecular dynamics.

\[C(\tau)=\Biggl\langle \sum v_{ij}(t_0) v_{ij}(t_0+\tau) \Biggr\rangle _ {t_0}\]

where \(v_{ij}(t_0)\) is the velocity at time \(t_0\). Then, the power spectra can be obtain from Fourier transform of the \(C (\tau)\)

\[I(\omega) \propto \int _0^{\infty} C (\tau) \cos (\omega \tau ) d \tau\]

Note

  • Time autocorrelation functions should decay to zero rapidly in strongly coupled molecular systems with a large number of atoms. Therefore, the selection of time step in MD sampling is crucial to capture the vibration periods. [8]

Caution

  • Autocorrelation function is sensitive to time step in MD simulation.



Reference#