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}}}\)

Atom Selecion¶

The variables \(n_{\textrm{from}}\) and \(n_{\textrm{to}}\) determining the atom pairs that contribute to the radial distribution can be selected in the AtomsFrom and AtomsTo blocks. Both can be selected using element names (Element), region names (Region) or atom indices (Atom).

RadialDistribution
   AtomsFrom
      Atom integer
      Element string
      Region string
   End
End

RadialDistribution

Type:

Block

Recurring:

True

Description:

All input related to radial distribution functions.

AtomsFrom

Type:

Block

Description:

Atom numbers or elements for the first set of atoms in the radial distribution.

Atom

Type:

Integer

Default value:

0

Recurring:

True

Description:

Atom number.

Element

Type:

String

Recurring:

True

Description:

Element Symbol Atom.

Region

Type:

String

Recurring:

True

Description:

Region name.

RadialDistribution

Type:

Block

Recurring:

True

Description:

All input related to radial distribution functions.

AtomsTo

Type:

Block

Description:

Atom numbers or elements for the second set of atoms in the radial distribution.

Atom

Type:

Integer

Default value:

0

Recurring:

True

Description:

Atom number.

Element

Type:

String

Recurring:

True

Description:

Element Symbol Atom.

Region

Type:

String

Recurring:

True

Description:

Region name.

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.

Inter- or intra-molecular atom-pairs By default, all \(n_{\textrm{from}}*n_{\textrm{to}}\) distances are included. Sometimes it can be convenient to view exclusively the distances between atoms within a molecule, or those between different molecules. This can be controlled with the keyword DistanceTypeSelection in the RadialDistribution block.

RadialDistribution
   DistanceTypeSelection [All | InterMolecular | IntraMolecular]
End

RadialDistribution

Type:

Block

Recurring:

True

Description:

All input related to radial distribution functions.

DistanceTypeSelection

Type:

Multiple Choice

Default value:

All

Options:

[All, InterMolecular, IntraMolecular]

Description:

Select only a certain type of interatomic distances.

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 frames read from the trajectory (as determined in the block TrajectoryInfo), and \(m\) is the number of discrete \(t\) values for which \(C(t)\) is computed. The value \(m\) can be set with the keywords MaxCorrelationTime or MaxFrame. This generally defaults to half the total number of frames read. Only when the pressure tensor is used (e.g. in the computation of the viscosity), the default value is smaller, at 10% of the total number of frames read. 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 SamplingFreq keyword in the Trajectory block of the MolecularDynanimcs settings low, preferably to 1).

Setting SamplingFreq to 1 results in the storage of a lot of superfluous information, which takes both time and disc space. For a few selected properties this can be circumvented with the MolecularDynamics%BinLog block. Properties requested in this block will be stored every step, even if SamplingFreq is set to a higher value. The trajectory analysis tool can take advantage of this for the computation of some autocorrelation functions, through the Property keyword.

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 time derivative, the power spectrum matches the IR spectrum.

Properties¶

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. Diffusion coefficient (equivalent to selecting Velocities).

6. Viscosity (equivalent to selecting PressureTensor, if only the off-diagonal elements are selected using VecElements)

7. DipoleMomentFromBinLog (approximately equivalent to DipoleMomentFromCharges if available)

8. ViscosityFromBinlog (equivalent to Viscosity if available)

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

AutoCorrelation

Type:

Block

Recurring:

True

Description:

All input related to auto correlation functions.

Property

Type:

Multiple Choice

Default value:

DipoleDerivativeFromCharges

Options:

[Velocities, DipoleMomentFromCharges, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor, Viscosity, DipoleMomentFromBinLog, ViscosityFromBinLog]

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.DipoleMomentFromBinLog and ViscosityFromBinLog allow the relevant properties (dipole moment and pressure tensor respectively) to be read from the BinLog section of the trajectory file. In the BinLog section requested properties are stored every step (even if SamplingFreq was set to a higher number than 1) but only if this was specifically requested at the start of the MD simulation.

When selecting options 7) or 8) the program attempts to read the dipole moment or the pressure tensor respectively from the BinLog section of the ams.rkf file. This will fail if the data is not stored there. The data will only be stored there if the user selected MoleclarDynamics%BinLog%DipoleMoment and MoleculatDynamics%BinLog%PressureTensor when setting up the MD simulation.

Some of the properties for which an autocorrelation function can be computed are simply read as is from the trajectory RKF file, but others are quite complex. For example, the dipole moment is a property obtained by reading the coordinates and the atomic charges, multiplying each atomic position vector with the corresponding charge, and then summing over all atoms. With the keyword WritePropertyToKF in the AutoCorrelation block, the user can choose to store not only the autocorrelation function, but also the property used to produce it.

AutoCorrelation
   WritePropertyToKF Yes/No
End

AutoCorrelation

Type:

Block

Recurring:

True

Description:

All input related to auto correlation functions.

WritePropertyToKF

Type:

Bool

Default value:

No

Description:

Write the selected property to the KF files for every requested frame

If a property involves atomic coordinates, then this generally means atomic coordinates wrapped into the periodic box at every time step. As a result, coordinates can make seemingly big jumps from one side of the box to another between two consecutive frames. The code is able to unwrap the coordinates, so that all atoms follow a continuous trajectory. For the DipoleMoment property, the wrapped coordinates will be used by default, but an experienced user can determine whether the wrapped or unwrapped coordinates are used with the keyword UnwrapCoordinates in the AutoCorrelation block.

AutoCorrelation
   UnwrapCoordinates [Auto | Yes | No]
End

AutoCorrelation

Type:

Block

Recurring:

True

Description:

All input related to auto correlation functions.

UnwrapCoordinates

Type:

Multiple Choice

Default value:

Auto

Options:

[Auto, Yes, No]

Description:

If the coordinates are involved in the requested property, those coordinates are wrapped into the box at each time step. If set to true, this keyword unwraps those coordinates so that the trajectory is continuous. If not provided the code uses automatic defaults.

Options¶

With the keywords MaxCorrelationTime (in fs) or MaxFrame the number of values \(m\) in the autocorrelation function (\(t = [0,t_{1},t_{2},....,t_{m}]\)) can be set. When both are set, MaxCorrelationTime takes precedent. In most cases the default value is half of the total number of frames \(n\) read from the trajectory. Only when the selected Property is ‘PressureTensor’, ‘Viscosity’ or ‘ViscosityFromBinLog’, the default is 10% of the trajectory.

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), region names (recurring keyword Region),or individual atom numbers (recurring keyword Atom).

AutoCorrelation
   Atoms
      Atom integer
      Element string
      Region 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.

Region

Type:

String

Recurring:

True

Description:

Region name.

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).

Index

Type:

Integer

Recurring:

True

Description:

Element of the property vector.

Keywords¶

The AutoCorrelation block holds the following keywords.

AutoCorrelation
   Atoms
      Atom integer
      Element string
      Region string
   End
   DataReading [Auto | AtOnce | BlockWise]
   MaxCorrelationTime float
   MaxFrame integer
   NPointsHighestFreq integer
   PerElement Yes/No
   Property [Velocities | DipoleMomentFromCharges | DiffusionCoefficient | DipoleDerivativeFromCharges | 
             PressureTensor | Viscosity | DipoleMomentFromBinLog | ViscosityFromBinLog]
   UnwrapCoordinates [Auto | Yes | No]
   UseAllValues Yes/No
   UseTimeDerivative
      Enabled Yes/No
   End
   VecElements
      Index integer
   End
   WritePropertyToKF Yes/No
End

AutoCorrelation

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.

Region

Type:

String

Recurring:

True

Description:

Region name.

DataReading

Type:

Multiple Choice

Default value:

Auto

Options:

[Auto, AtOnce, BlockWise]

Description:

The KF data can be read in and handled once, or blockwise. The former is memory intensive, but mostly faster. If Auto is selected, the data is read at once if it is less than 1 GB, and blockwise if it is more.

MaxCorrelationTime

Type:

Float

Description:

The maximum correlation time in fs. The default is half the simulation time, except when the PressureTensor is read from the ams.rkf, in which case it is 10 percent of the simulation time. The PressureTensor is read when the Properties PressureTensor, Viscosity, or ViscosityFromBinLog are selected.

MaxFrame

Type:

Integer

Description:

The maximum number of frames for which the autocorrelation function will be computed. The default is half of the number of provided frames. Determines the same settings as MaxCorrelationTime. If both are set, MaxCorrelationTime will take precedence.

NPointsHighestFreq

Type:

Integer

Default value:

4

Description:

The number of points (timesteps) used for the highest frequency displayed in spectrum. This determines up to which frequency the spectrum is displayed. If the spacing between time-steps used for the ACF is 1 fs, then by default the maximum frequency displayed is 0.25 fs-1 (or 8339 cm-1). This corresponds to a (default) value of NPointsHighestFreq of 4. A higher number selected here, will result in a lower maximum frequency returned by the program. The lowest possible value (spectrum up to highest possible frequency) is 2.

PerElement

Type:

Bool

Default value:

No

Description:

Compute ACF for all elements in the system. Any other settings in the block will be used.

Property

Type:

Multiple Choice

Default value:

DipoleDerivativeFromCharges

Options:

[Velocities, DipoleMomentFromCharges, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor, Viscosity, DipoleMomentFromBinLog, ViscosityFromBinLog]

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.DipoleMomentFromBinLog and ViscosityFromBinLog allow the relevant properties (dipole moment and pressure tensor respectively) to be read from the BinLog section of the trajectory file. In the BinLog section requested properties are stored every step (even if SamplingFreq was set to a higher number than 1) but only if this was specifically requested at the start of the MD simulation.

UnwrapCoordinates

Type:

Multiple Choice

Default value:

Auto

Options:

[Auto, Yes, No]

Description:

If the coordinates are involved in the requested property, those coordinates are wrapped into the box at each time step. If set to true, this keyword unwraps those coordinates so that the trajectory is continuous. If not provided the code uses automatic defaults.

UseAllValues

Type:

Bool

Default value:

No

Description:

By default the same number of values are used for each t-step in the ACF. This has the advantage that all values in the ACF are equally reliable, but it does mean that for the smaller timesteps much of the data is not used. To switch this off and use all data, UseAllValues can be set to true

UseTimeDerivative

Type:

Block

Description:

Possibly use the time derivative of the selected property (e.g. velocity or dipole moments).

Enabled

Type:

Bool

Default value:

No

Description:

Enable the use of the time derivative of the property.

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).

Index

Type:

Integer

Recurring:

True

Description:

Element of the property vector.

WritePropertyToKF

Type:

Bool

Default value:

No

Description:

Write the selected property to the KF files for every requested frame

Description¶

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

\(MSD(t) = \langle \frac{1}{k} \sum_i^k ||\textbf{A}_i(t') - \textbf{A}_i(t)||^2 \rangle _{t'}\)

Here \(\textbf{A}_i\) is usually the coordinate vector for atom \(i\) out of a set of \(k\) atoms. The average runs over all values of \(t'\), or the 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 frames read from 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 keywords MaxCorrelationTime or MaxFrame. It defaults to half the total number of frames read from the trajectory. 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}\).

Properties¶

The only property \(\textbf{A}\) of real relevance for the computation of the mean square displacement is a set of coordinates. Nonetheless, the Property keyword accepts four different options.

1. Coordinates

2. Diffusion coefficient

3. Ionic conductivity

Option 2) is equivalent to option 1). Option 3) uses charge weighted coordinates as the property \(\textbf{A}\).

MeanSquareDisplacement
   Property [Coords | DiffusionCoefficient | Conductivity]
End

MeanSquareDisplacement

Type:

Block

Recurring:

True

Description:

All input related to mean square displacement functions.

Property

Type:

Multiple Choice

Default value:

Coords

Options:

[Coords, DiffusionCoefficient, Conductivity]

Description:

Compute the MSD from the property selected here (from rkf). Selecting DiffusionCoefficient is equivalent to selecting the property Coords.

When read from file, the atomic coordinates will be wrapped inside the periodic box at every time step. As a result, coordinates can make seemingly big jumps from one side of the box to another between two consecutive frames. By default, the coordinates are unwrapped before the computation of the mean square displacement, so that all atoms follow a continuous trajectory. For the experienced users, the option exist to manually determine whether the wrapped or unwrapped coordinates are used with the keyword UnwrapCoordinates in the MeanSquareDisplacement block.

MeanSquareDisplacement
   UnwrapCoordinates [Auto | Yes | No]
End

MeanSquareDisplacement

Type:

Block

Recurring:

True

Description:

All input related to mean square displacement functions.

UnwrapCoordinates

Type:

Multiple Choice

Default value:

Auto

Options:

[Auto, Yes, No]

Description:

If the coordinates are involved in the requested property, those coordinates are wrapped into the box at each time step. If set to true, this keyword unwraps those coordinates so that the trajectory is continuous. If not provided the code uses automatic defaults.

Since the coordinates used to compute the mean square displacement often differ from the coordinates as read from the trajectory (they will be unwrapped, so that the trajectory becomes continuous). With the keyword WritePropertyToKF in the MeanSquareDisplacement block, the user can choose to store not only the mean square displacement results, but also the property used to produce it.

MeanSquareDisplacement
   WritePropertyToKF Yes/No
End

MeanSquareDisplacement

Type:

Block

Recurring:

True

Description:

All input related to mean square displacement functions.

WritePropertyToKF

Type:

Bool

Default value:

No

Description:

Write the selected property to the KF files for every requested frame

Options¶

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), region names (recurring keyword Region),or individual atom numbers (recurring keyword Atom).

MeanSquareDisplacement
   Atoms
      Atom integer
      Element string
      Region string
   End
End

MeanSquareDisplacement

Type:

Block

Recurring:

True

Description:

All input related to mean square displacement 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.

Region

Type:

String

Recurring:

True

Description:

Region name.

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 mean square displacement 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).

Index

Type:

Integer

Recurring:

True

Description:

Element of the property vector.

The center of mass of each selected molecule can be used in the MSD computation. In that case, the selected atoms will be divided into groups belonging to different molecules, and the center of mass of each group is then computed. This will only work if the molecules remain stable throughout the simulation, which can only really be guaranteed if the ForceField or the GFNFF engines are used. If the bonding does change, the computation will abort by default. If, however, UseMolecularCentersOfMass%CheckForBondChanges is set to No, then any bond changes will be ignored, and the computation continues as if the molecules did not change.

Atomic masses used for the molecular centers of mass can be altered by providing a system block in the input. The elements in the input system block need to match the elements in the molecular dynamics simulation. If no system block is present in the input, the system block from the molecular dynamics simulation will be used.

#!bin/sh

"$AMSBIN/analysis" << eor

Task MeanSquareDisplacement

System
   Atoms
           O      -0.7113934019       1.6140994026      -0.7075521284 mass=15.999
           H      -0.0425480297       2.3360319223      -0.8302756333 mass=2.014
           H      -0.5643276324       1.0450643077      -1.4556911815 mass=2.014
           O       1.1006174021      -7.1200763645      -3.6017659954 mass=15.999
           H       1.4647560784      -6.1967314913      -3.5522086623 mass=2.014
           H       1.8229257860       7.2170281524      -3.7443300356 mass=2.014
   End
End

TrajectoryInfo
    Trajectory
        KFFilename ams.results/ams.rkf
    End
End

MeanSquareDisplacement
    Property Coords
    UseMolecularCentersOfMass
      Enabled Yes
    End
End
eor

MeanSquareDisplacement
   UseMolecularCentersOfMass
      CheckForBondChanges Yes/No
      Enabled Yes/No
   End
End

MeanSquareDisplacement

Type:

Block

Recurring:

True

Description:

All input related to mean square displacement functions.

UseMolecularCentersOfMass

Type:

Block

Description:

Specifications on using the center of mass coordinates per molecule instead of coordinates per atom. Only relevant for properties that rely on coordinates: Coords, DiffusionCoefficient, and Conductitvity.

CheckForBondChanges

Type:

Bool

Default value:

Yes

Description:

Check for each frame if the bonds still correspond to the input bonds. If they changed, an error is thrown. If this is set to true, it is assumed that the molecules stay in tact.

Enabled

Type:

Bool

Default value:

No

Description:

Use the center of mass coordinates per molecule. Only relevant for properties that rely on coordinates: Coords, DiffusionCoefficient, and Conductitvity. The bonds from the MD input system are used.

Keyword¶

The MeanSquareDisplacement block holds the following keywords.

MeanSquareDisplacement
   Atoms
      Atom integer
      Element string
      Region string
   End
   DataReading [Auto | AtOnce | BlockWise]
   MaxCorrelationTime float
   MaxFrame integer
   PerElement Yes/No
   Property [Coords | DiffusionCoefficient | Conductivity]
   StartTimeSlope float
   UnwrapCoordinates [Auto | Yes | No]
   UseAllValues Yes/No
   UseMolecularCentersOfMass
      CheckForBondChanges Yes/No
      Enabled Yes/No
   End
   VecElements
      Index integer
   End
   WritePropertyToKF Yes/No
End

MeanSquareDisplacement

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.

Region

Type:

String

Recurring:

True

Description:

Region name.

DataReading

Type:

Multiple Choice

Default value:

Auto

Options:

[Auto, AtOnce, BlockWise]

Description:

The KF data can be read in and handled once, or blockwise. The former is memory intensive, but mostly faster. If Auto is selected, the data is read at once if it is less than 1 GB, and blockwise if it is more.

MaxCorrelationTime

Type:

Float

Description:

The maximum correlation time in fs. The default is half the simulation time.

MaxFrame

Type:

Integer

Description:

The maximum number of frames for which the mean square displacement function will be computed. The default is half of the number of provided frames. Determines the same settings as MaxCorrelationTime. If both are set, MaxCorrelationTime will take precedence.

PerElement

Type:

Bool

Default value:

No

Description:

Compute MSD for all elements in the system. Any other settings in thie block will be used.

Property

Type:

Multiple Choice

Default value:

Coords

Options:

[Coords, DiffusionCoefficient, Conductivity]

Description:

Compute the MSD from the property selected here (from rkf). Selecting DiffusionCoefficient is equivalent to selecting the property Coords.

StartTimeSlope

Type:

Float

Default value:

0.0

Description:

The MSD has a nonlinear regime at short timescales, and a linear regime at long timescales. To determine the slope, the starting point for the linear regime has to be determined. This keyword sets the starting time in fs. If set to zero, the starttime will be automatically determined.

UnwrapCoordinates

Type:

Multiple Choice

Default value:

Auto

Options:

[Auto, Yes, No]

Description:

If the coordinates are involved in the requested property, those coordinates are wrapped into the box at each time step. If set to true, this keyword unwraps those coordinates so that the trajectory is continuous. If not provided the code uses automatic defaults.

UseAllValues

Type:

Bool

Default value:

No

Description:

By default the same number of values are used for each t-step in the MSD. This has the advantage that all values in the MSD are equally reliable, but it does mean that for the smaller timesteps much of the data is not used. To switch this off and use all data, UseAllValues can be set to true

UseMolecularCentersOfMass

Type:

Block

Description:

Specifications on using the center of mass coordinates per molecule instead of coordinates per atom. Only relevant for properties that rely on coordinates: Coords, DiffusionCoefficient, and Conductitvity.

CheckForBondChanges

Type:

Bool

Default value:

Yes

Description:

Check for each frame if the bonds still correspond to the input bonds. If they changed, an error is thrown. If this is set to true, it is assumed that the molecules stay in tact.

Enabled

Type:

Bool

Default value:

No

Description:

Use the center of mass coordinates per molecule. Only relevant for properties that rely on coordinates: Coords, DiffusionCoefficient, and Conductitvity. The bonds from the MD input system are used.

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).

Index

Type:

Integer

Recurring:

True

Description:

Element of the property vector.

WritePropertyToKF

Type:

Bool

Default value:

No

Description:

Write the selected property to the KF files for every requested frame

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, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor, Viscosity, DipoleMomentFromBinLog, ViscosityFromBinLog]

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.DipoleMomentFromBinLog and ViscosityFromBinLog allow the relevant properties (dipole moment and pressure tensor respectively) to be read from the BinLog section of the trajectory file. In the BinLog section requested properties are stored every step (even if SamplingFreq was set to a higher number than 1) but only if this was specifically requested at the start of the MD simulation.

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 ({\bf{r}}(0) - {\bf{r}}(t))^2 \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 mean square displacement functions.

Property

Type:

Multiple Choice

Default value:

Coords

Options:

[Coords, DiffusionCoefficient, Conductivity]

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.