Trajectory Analysis¶

Description¶

A radial distribution function \(g(r)\), or pair correlation function, is a density of distances between particles, relative to the average distance density. The x-axis variable represents a distance \(r\), while the y-axis represents the relative density of that distance. For a complete homogeneous system of particles the \(g(r)\) values for the distances between all particles equals 1 everywhere.

Two sets of atoms \(\mathbb{S}_{\textrm{from}}\) and \(\mathbb{S}_{\textrm{to}}\), of length \(n_{\textrm{from}}\) and \(n_{\textrm{to}}\) respectively, are specified with the keywords AtomsFrom and AtomsTo in the RadialDistribution block. As a result the program computes \(n_{\textrm{from}}*n_{\textrm{to}}\) distances \(r_{ij}^s\) between atom \(i\) in \(\mathbb{S}_{\textrm{from}}\) and atom \(j\) in \(\mathbb{S}_{\textrm{to}}\) for each trajectory frame \(s\) out of a total of \(n_{\textrm{frames}}\) frames.

A normalized histogram is then computed from these distances, resulting in a function \(N(r)\).

\(N(r)=\frac{1}{n_{\textrm{frames}}} \sum_{s=1}^{n_{\textrm{frames}}} \sum_{i=1}^{n_{\textrm{from}}}\sum_{j=1}^{n_{\textrm{to}}} \delta(r_{ij}^s-r)\).

This histogram is converted to a density, by dividing all values \(N(r)\) with the volume \(V(r)= 4 \pi r^2 dr\) of a sphere-slice at radius \(r\) with thickness \(dr\).

The density is further converted to a relative density by dividing with the total density of the system \(\rho_{\textrm{tot}} = \frac{n_{\textrm{from}}*n_{\textrm{to}}}{V_{\textrm{tot}}}\), yielding the final radial distribution function \(g(r)\).

\(g(r) = \frac{N(r)}{V(r)*\rho_{\textrm{tot}}}\)

Options¶

Non-periodic systems The above equation assumes that the volume \(V_{\textrm{tot}}\) of the system is a well-defined quantity. This assumption is correct for systems with 3D periodicity, where the \(V_{\textrm{tot}}\) is defined as the volume of the periodic cell. In such a system the value of \(r\) can be no larger than \(r_{\textrm{max}}\), the radius of the largest sphere that can be placed inside the periodic cell.

If a system is non-periodic in one or more direction, then the program still computes a \(g(r)\), only if the radius \(r_{max}\) is supplied by the user with the Range keyword in the RadialDistribution block. The radius is the second value supplied.

RadialDistribution
   Range float_list
End

RadialDistribution

Type:	Block
Recurring:	True
Description:	All input related to radial distribution functions.

Range

Type:	Float List
Description:	Either one, two, or three real values. If one it is the stepsize. If two, it is the minimum value and the maximum value. If three, it is the minimum value, the maximum value, and the stepsize. The stepsize overrides NBins.

In this case the volume \(V_{\textrm{tot}}\) is assumed to be the volume of a sphere with radius \(r_{\textrm{max}}\).

NPT simulations The above equation further assumes that the volume \(V_{\textrm{tot}}\) is constant throughout the simulation. The \(g(r)\) of the trajectory from an NPT simulation can still be computed, and in this case \(V_{\textrm{tot}}\) is the average value of the volume of the periodic cell.

Simulations with varying numbers of atoms The above equation also assumes that \(n_{\textrm{from}}\) and \(n_{\textrm{to}}\) remain constant throughout the simulation. However, in a Molecular Gun simulation particles can be added to the system, and in a GCMC simulation particles can be both added and removed from the system. Nonetheless, the program still computes a \(g(r)\) in these situations.

If the AtomsFrom and AtomsTo blocks contain element names (supplied with the recurring Element keyword), then every time atoms are added to or removed from the system, the sets of atoms \(\mathbb{S}_{\textrm{from}}\) and \(\mathbb{S}_{\textrm{to}}\) are re-evaluated.

If the AtomsFrom and AtomsTo blocks contain atom numbers (supplied with the recurring Atom keyword), these numbers are updated in the sets \(\mathbb{S}_{from}\) and \(\mathbb{S}_{to}\) every time atoms are added to or removed from the system. If one of the atoms from the set disappears, the number of distances contributing to the \(g(r)\) decreases.

Note: Currently, the values of \(n_{from}\) and \(n_{to}\) in the normalization factor are taken from the last frame of the simulation.

Warning: If multiple trajectories are supplied, and the number of atoms changes between the end of one trajectory and the beginning of another, this may result in an error in the atom numbers used by the program internally.

Description¶

An autocorrelation function \(C(t)\) describes the average correlation (overlap) of a (vector) property \(\textbf{A}\) with itself as a function of time.

\(C(t) = \langle \textbf{A}(0) \cdot \textbf{A}(t)) \rangle\)

The average runs over all time-intervals \(\left( t_{0}, t_{0}+t \right),\left( t_{1}, t_{1}+t \right),...,\left( t_{N}, t_{N}+t \right)\), with \(t_{N} = t_{n} - t_{m}\). Here \(n\) is the total number of simulation steps in the trajectory, and \(m\) is the number of discrete \(t\) values for which \(C(t)\) is computed. The value \(m\) can be set with the keyword MaxStep, and defaults to half the total number of simulation steps. As stated above, the default average runs over the same number of time intervals for each value of \(t\). If UseAllValues is selected, the average runs over all available time intervals for each value of \(t\), which is large for small \(t\), and smallest (\(N\)) for \(t_{m}\).

The normalized autocorrelation function \(c(t)\) describes the decorrelation of the property with time, and always starts at 1.0 at \(t=0\).

\(c(t) = \frac{\langle \textbf{A}(0) \cdot \textbf{A}(t)) \rangle}{\langle \textbf{A}(0) \cdot \textbf{A}(0)) \rangle}\)

In most cases short timescale fluctuations are important, so frequent storage of the desired property is required (when preparing the molecular dynamics simulation, set the Frequency keyword in the Trajectory block of the MolecularDynanimcs settings low, preferably to 1).

A power spectrum is automatically computed by Fourier transform of the autocorrelation function, and provides information on the frequencies of the signal. When the selected property is the dipole moment or the dipole moment derivative, the power spectrum matches the IR spectrum.

Options¶

Autocorrelation functions can be computed for different simulation properties: 1) Dipole moments from coordinates and atomic charges 2) Dipole moment derivatives from velocities and atomic charges 3) Velocities 4) The pressure tensor 5) User provided values. Selecting 6) Diffusion coefficient is equivalent to selecting Velocities.

AutoCorrelation
   Property [Velocities | DipoleMomentFromCharges | InputValues | DiffusionCoefficient | 
             DipoleDerivativeFromCharges | PressureTensor | Viscosity]
End

AutoCorrelation

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

Property

Type:	Multiple Choice
Default value:	DipoleDerivativeFromCharges
Options:	[Velocities, DipoleMomentFromCharges, InputValues, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor, Viscosity]
Description:	Compute the ACF either from velocities (from rkf), the dipole moment (from coordinates and atomic charges in rkf), the dipole moment derivative (from velocities and atomic charges in rkf), from the pressure tensor (from rkf), or from values specified in input. Selecting DiffusionCoefficient is equivalent to selecting Velocities. The default, DipoleDerivativeFromCharges, results in the computation of an IR spectrum.

With the keyword MaxStep the number of values \(m\) in the autocorrelation function (\(t = [0,t_{1},t_{2},....,t_{m}]\)) can be set. The default value is half of the total number of simulation steps \(n\) used.

A subset of atoms for which the property \(\textbf{A}\) should be selected/computed can be provided in the block Atoms. The block can contain element names (recurring keyword Element), or individual atom numbers (recurring keyword Atom).

AutoCorrelation
   Atoms
      Atom integer
      Element string
   End
End

AutoCorrelation

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

Atoms

Type:	Block
Description:	Relevant if Property is set to Velocities, DipoleMomentFromCharges, DipoleDerivativeFromCharges, or DiffusionCoefficient. Atom numbers or elements for the set of atoms for which the property is read/computed. By default all atoms are used.

Atom

Type:	Integer
Recurring:	True
Description:	Atom number.

Element

Type:	String
Recurring:	True
Description:	Element Symbol Atom.

A subset of vector elements can be provided with the subblock VecElements. By default all vector elements found on the RKF file will be used.

AutoCorrelation
   VecElements
      Index integer
   End
End

AutoCorrelation

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

VecElements

Type:	Block
Description:	A set of indices referring to a subset of the property vector. Works in combination with the atoms block. For example, in combination with the property Velocities, the Atoms block allows the selection of a subset of atoms, while the VecElelements block allows the selection of a subset of vector elements (e.g. 1 and 2 for the elements x and y). Currently not implemented with InputValues.

Index

Type:	Integer
Recurring:	True
Description:	Element of the property vector.

Description¶

The mean square displacement \(MSD(t)\) describes the average change of a (vector) property \(\textbf{A}\) over time. This property is usually a set of atom coordinate vectors, but the implementation is entirely general.

\(MSD(t) = \langle [\textbf{A}(0) - \textbf{A}(t)]^2 \rangle\)

The average runs over all time-intervals \(\left( t_{0}, t_{0}+t \right),\left( t_{1}, t_{1}+t \right),...,\left( t_{N}, t_{N}+t \right)\), with \(t_{N} = t_{n} - t_{m}\). Here \(n\) is the total number of simulation steps in the trajectory, and \(m\) is the number of discrete \(t\) values for which \(MSD(t)\) is computed. The value \(m\) can be set with the keyword MaxStep, and defaults to half the total number of simulation steps. As stated above, the default average runs over the same number of time intervals for each value of \(t\). If UseAllValues is selected, the average runs over all available time intervals for each value of \(t\), which is large for small \(t\), and smallest (\(N\)) for \(t_{m}\).

The most common use of the mean square displacement is for the computation of the diffusion coefficient, in which case the selected property is the set of atom coordinates (Coords). The diffusion coefficient is proportionate to the slope of this mean square displacement function, and therefore this slope is automatically computed. The function \(MSD(t)\) becomes linear only after an initial time interval, and the user can set this initial time with the keyword StartTimeSlope. If not provided, this start time is automatically determined. To allow the user to determine if the linear regime has been sufficiently sampled, the slope of \(MSD(t)\) as a function of \(t\) is computed as well. If the slope has not converged to a stable value, the user should select a larger value of \(t_m\) or continue the molecular dynamics simulation for a longer time.

Options¶

Autocorrelation functions can be computed for different simulation properties: 1) Coordinates 2) User provided values. Option 3) Diffusion coefficient is equivalent to option 1).

MeanSquareDisplacement
   Property [Coords | InputValues | DiffusionCoefficient]
End

MeanSquareDisplacement

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

Property

Type:	Multiple Choice
Default value:	Coords
Options:	[Coords, InputValues, DiffusionCoefficient]
Description:	Compute the MSD from the property selected here (from rkf). Selecting DiffusionCoefficient is equivalent to selecting the property Coords.

A subset of atoms for which the property \(\textbf{A}\) should be selected/computed can be provided in the block Atoms. The block can contain element names (recurring keyword Element), or individual atom numbers (recurring keyword Atom).

MeanSquareDisplacement
   Atoms
      Atom integer
      Element string
   End
End

MeanSquareDisplacement

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

Atoms

Type:	Block
Description:	Relevant if Property is set to any quantity that is available per atom (Coords, DiffusionCoefficient). Atom numbers or elements for the set of atoms for which the property is read/computed are provided here. By default all atoms are used.

Atom

Type:	Integer
Recurring:	True
Description:	Atom number.

Element

Type:	String
Recurring:	True
Description:	Element Symbol Atom.

A subset of vector elements can be provided with the subblock VecElements. By default all vector elements found on the RKF file will be used.

MeanSquareDisplacement
   VecElements
      Index integer
   End
End

MeanSquareDisplacement

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

VecElements

Type:	Block
Description:	A set of indices referring to a subset of the property vector. Works in combination with the atoms block. For example, in combination with the property Coords, the Atoms block allows the selection of a subset of atoms, while the VecElelements block allows the selection of a subset of vector elements (e.g. 1 and 2 for the elements x and y). Currently not implemented with InputValues.

Index

Type:	Integer
Recurring:	True
Description:	Element of the property vector.

From Velocity Autocorrelation¶

The diffusion coefficient can be defined as an integral over the velocity autocorrelation function.

\(D = \frac{1}{d} \int_{t=0}^{t=t_{max}} \langle \textbf{v}(0) \cdot \textbf{v}(t)) \rangle dt\)

The factor \(\frac{1}{d}\) corrects for the dimension of the system, which we assume to equal the length of the provided vector \(\textbf{v}\). The dimension \(d\) equals 3, unless specified otherwise in the subblock VecElements. In a system that is periodic in less than 3 dimensions, it may make sense to provide only the vector elements along the periodic dimensions. By default, however, all vector elements provided are used.

The diffusion coefficient is computed if the task AutoCorrelation is selected, and if in the AutoCorrelation block DiffusionCoefficient is selected as the Property. The result is completely equivalent to selecting the task AutoCorrelation with Velocities as the Property keyword.

$AMSBIN/analysis <<eor
   Task AutoCorrelation
   AutoCorrelation
      Property DiffusionCoefficient
   End
eor

AutoCorrelation

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

Property

Type:	Multiple Choice
Default value:	DipoleDerivativeFromCharges
Options:	[Velocities, DipoleMomentFromCharges, InputValues, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor, Viscosity]
Description:	Compute the ACF either from velocities (from rkf), the dipole moment (from coordinates and atomic charges in rkf), the dipole moment derivative (from velocities and atomic charges in rkf), from the pressure tensor (from rkf), or from values specified in input. Selecting DiffusionCoefficient is equivalent to selecting Velocities. The default, DipoleDerivativeFromCharges, results in the computation of an IR spectrum.

From Mean Square Displacement¶

The mean square displacement becomes linear with time, after an initial time interval. We can therefore define the linear part of the function as follows.

\(MSD(t) = {\langle \textbf{r}(0) - \textbf{r}(t)) \rangle} = at + b\),

with \(a\) as the slope of the function. The diffusion coefficient is proportional to the slope \(a\).

\(D = \frac{1}{2d} a\)

Here, \(d\) is the dimensionality of the system, or the length of the provided vector \(\textbf{r}\). The dimension \(d\) equals 3, unless specified otherwise in the subblock VecElements. In a system that is periodic in less than 3 dimensions, it may make sense to provide only the vector elements along the periodic dimensions. By default, however, all vector elements provided are used.

The diffusion coefficient is computed if the task MeanSquareDisplacement is selected, and if in the MeanSquareDisplacement block DiffusionCoefficient is selected as the Property. The result is completely equivalent to selecting the task MeanSquareDisplacement with Coords as the Property keyword.

$AMSBIN/analysis <<eor
   Task MeanSquareDisplacement
   MeanSquareDisplacement
      Property DiffusionCoefficient
   End
eor

MeanSquareDisplacement

Type:	Block
Recurring:	True
Description:	All input related to auto correlation functions.

Property

Type:	Multiple Choice
Default value:	Coords
Options:	[Coords, InputValues, DiffusionCoefficient]
Description:	Compute the MSD from the property selected here (from rkf). Selecting DiffusionCoefficient is equivalent to selecting the property Coords.

In both cases, a subset of atoms can be selected with the sublock Atoms, and a subset of vector elements (in this case elements 1=X, 2=Y, 3=Z for the Cartesian coordinates) can be selected with the subblock VecElements.

The value of the diffusion coefficient is written to the output, as well as to the KF file.