# Trajectory Analysis¶

`analysis`

is a standalone program that performs analysis of molecular dynamics trajectories created with AMS. It can produce histograms and radial distribution functions. It is used under the hood in AMSmovie (**MD Properties** menu bar).

## New in Trajectory Analysis-2023¶

Intra- and inter-molecular radial distribution functions

Properties used to compute autocorrelation or mean square displacement can be stored to file

Manual selection of coordinate unwrapping for autocorrelation or mean square displacement functions

## General Input Options¶

This example computes the oxygen-oxygen radial distribution function of a MD simulation using the analysis utility program:

```
Task RadialDistribution
TrajectoryInfo
Trajectory
KFFilename ams.results/ams.rkf
Range 1 1000 2
End
End
RadialDistribution
NBins 1000
AtomsFrom
Element O
End
AtomsTo
Element O
End
End
eor
```

The analysis program reads one or more trajectory files (filename.rkf) from an AMS molecular dynamics (MD) or a Grand Canonical Monte Carlo (GCMC) simulation.
The file information is supplied in the `TrajectoryInfo`

input block.
In this block, a separate `Trajectory`

subblock needs to be supplied for each trajectory file.
The `Trajectory`

subblock contains a mandatory keyword `KFFilename`

, and an optional keyword `Range`

.
The latter contains the initial frame to be read, the final frame to be read, and optionally the stepsize.
By default all frames on the trajectory file are read.

```
TrajectoryInfo
NBlocksToCompare integer
Trajectory
KFFilename string
Range integer_list
StepSize integer
End
End
```

`TrajectoryInfo`

- Type
Block

- Description
All the info regarding the reading of the trajectory files.

`NBlocksToCompare`

- Type
Integer

- Default value
1

- Description
Get an error estimate by comparing histograms for NBLocks time blocks of the trajectory.

`Trajectory`

- Type
Block

- Recurring
True

- Description
All info regarding the reading of a single trajectory file.

`KFFilename`

- Type
String

- Default value
ams.rkf

- Description
The name of the AMS trajectory file.

`Range`

- Type
Integer List

- Description
One or two values: start frame, and optionally end frame. By default the first and last frame are read.

`StepSize`

- Type
Integer

- Default value
1

- Description
The step size at which frames are read from the RKF (default 1, every frame is read).

All tools in the analysis program provide an option to obtain information on the equilibration of the simulation.
If the optional keyword `NBlocksToCompare`

in the `TrajectoryInfo`

block
is set to a value \(N\) higher than 1, the trajectory is divided into \(N\) blocks, and the analysis results for each block are compared.
The variation in the analysis result is provided as a standard deviation.

## Radial Distribution Function (RDF)¶

The Analysis tool computes radial distribution functions \(g(r)\) if the `Task`

keyword is set to RadialDistribution.

```
Task [RadialDistribution | Histogram | AutoCorrelation | MeanSquareDisplacement | AverageBinPlot]
```

`Task`

- Type
Multiple Choice

- Options
[RadialDistribution, Histogram, AutoCorrelation, MeanSquareDisplacement, AverageBinPlot]

- Description
The analysis task.

Further details on the radial distribution functions are then set in the `RadialDistribution`

block.
If more than one `RadialDistribution`

block is present in the input, more than one radial distribution function will be computed.
The result is printed to output as text, as well as stored in a binary file (plot.kf).

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

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

## Histogram¶

The Analysis program computes histograms if the `Task`

keyword is set to Histogram.

```
Task [RadialDistribution | Histogram | AutoCorrelation | MeanSquareDisplacement | AverageBinPlot]
```

`Task`

- Type
Multiple Choice

- Options
[RadialDistribution, Histogram, AutoCorrelation, MeanSquareDisplacement, AverageBinPlot]

- Description
The analysis task.

Further details on the histogram need to be specified in the `Histogram`

block.
If more than one `Histogram`

block is present in the input, more than one histogram will be computed.
The result is printed to output as text, as well as stored in a binary file (plot.kf).
By default the histogram contains the number of occurrences of a certain value,
but the normalized occurrence is provided if the keyword
`Normalized`

in the `Histogram`

block is specified.

```
Histogram
Normalized Yes/No
End
```

`Histogram`

`Normalized`

- Type
Bool

- Default value
No

- Description
Give the normalized histogram.

Histograms can be computed for every quantity stored on the molecular dynamics
trajectory file (ams.rkf) in the section History. Example quantities are
`PotentialEnergy`

, `KineticEnergy`

, `TotalEnergy`

, `Temperature`

. In
the histogram block, this quantity is selected with the keyword `Variable`

in
the `Axis`

subblock. If more than one `Axis`

subblock is present, the
dimensionality of the histogram is increased: Three `Axis`

subblocks result
in a 3D histogram.

For each histogram axis, the number of bins can be selected with the `NBins`

keyword in the `Axis block`

,
in which case the range of values along each axis is automatically determined.
The default `NBins`

value is 100.

Alternatively, a range and a stepsize can be selected with the keyword `Range`

in the `Axis`

subblock.
The keyword `Range`

can contain one, two, or three values:
1: Only a stepsize.
2: A smallest value and a largest value.
3: A smallest value, a largest value, and the stepsize.

```
Histogram
Axes
Axis
NBins integer
Range float_list
Variable string
End
End
End
```

`Histogram`

- Type
Block

- Recurring
True

- Description
All input related to histograms.

`Axes`

- Type
Block

- Description
Specifications for the histogram axes.

`Axis`

- Type
Block

- Recurring
True

- Description
Specifications for a single histogram axis.

`NBins`

- Type
Integer

- Default value
100

- Description
The number of bins along the histogram axis.

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

`Variable`

- Type
String

- Description
The quantity along the histogram axis.

## Autocorrelation Functions¶

The Analysis program computes autocorrelation functions (ACF) if the `Task`

keyword is set to AutoCorrelation.

```
Task [RadialDistribution | Histogram | AutoCorrelation | MeanSquareDisplacement | AverageBinPlot]
```

`Task`

- Type
Multiple Choice

- Options
[RadialDistribution, Histogram, AutoCorrelation, MeanSquareDisplacement, AverageBinPlot]

- Description
The analysis task.

Further details need to be specified in the `AutoCorrelation`

block.
If more than one `AutoCorrelation`

block is present in the input, more than one ACF will be computed.
The result is printed to output as text, as well as stored in a binary file (plot.kf).

```
AutoCorrelation
Atoms
Atom integer
Element string
End
DataReading [Auto | AtOnce | BlockWise]
InputValues
Values float_list
End
MaxFrame integer
NPointsHighestFreq integer
PerElement Yes/No
Property [Velocities | DipoleMomentFromCharges | InputValues | DiffusionCoefficient |
DipoleDerivativeFromCharges | PressureTensor | Viscosity]
TimeStep float
UnwrapCoordinates [Auto | Yes | No]
UseAllValues Yes/No
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.

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

`InputValues`

- Type
Block

- Description
Relevant if Property is set to InputValues. All input values (a vector on each line) need to be provided in this block, using the keyword Values (possibly multiple times).

`Values`

- Type
Float List

- Recurring
True

- Description
The values at each step (on a single line)

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

`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, 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.

`TimeStep`

- Type
Float

- Description
Relevant if Property is set to InputValues. The time separating the entries (in fs). If Property is set to any of the other quantities, it can be read from an RKF file, and the timestep is read from the RKF file as well. The read value then overrides this keyword.

`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

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

`WritePropertyToKF`

- Type
Bool

- Default value
No

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

### 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 keyword `MaxFrame`

, and defaults to half 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 `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.

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

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 keyword `MaxFrame`

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 frames \(n\) read from 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`

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

## Mean Square Displacement¶

The Analysis program computes mean square displacements (MSD) if the `Task`

keyword is set to MeansSquareDisplacement.

```
Task [RadialDistribution | Histogram | AutoCorrelation | MeanSquareDisplacement | AverageBinPlot]
```

`Task`

- Type
Multiple Choice

- Options
[RadialDistribution, Histogram, AutoCorrelation, MeanSquareDisplacement, AverageBinPlot]

- Description
The analysis task.

Further details need to be specified in the `MeanSquareDisplacement`

block.
If more than one `MeanSquareDisplacement`

block is present in the input, more than one MSD will be computed.
The result is printed to output as text, as well as stored in a binary file (plot.kf).

```
MeanSquareDisplacement
Atoms
Atom integer
Element string
End
DataReading [Auto | AtOnce | BlockWise]
InputValues
Values float_list
End
MaxFrame integer
PerElement Yes/No
Property [Coords | InputValues | DiffusionCoefficient]
StartTimeSlope float
TimeStep float
UnwrapCoordinates [Auto | Yes | No]
UseAllValues Yes/No
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.

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

`InputValues`

- Type
Block

- Description
Relevant if Property is set to InputValues. All input values (a vector on each line) need to be provided in this block, using the keyword Values (possibly multiple times).

`Values`

- Type
Float List

- Recurring
True

- Description
The values at each step (on a single line)

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

`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, InputValues, DiffusionCoefficient]

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

`TimeStep`

- Type
Float

- Description
Relevant if Property is set to InputValues. The time separating the entries (in fs). If Property is set to any of the other quantities, it can be read from an RKF file, and the timestep is read from the RKF file as well. The read value then overrides this keyword.

`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 seti 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

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

`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 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 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 keyword `MaxFrame`

, and 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}\).

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.

### Properties¶

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.

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 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 seti 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 auto correlation 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`

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

## AverageBinPlot¶

The Analysis program can plot two arbitrary properties, present on the RKF file, against each other averaged over each bin if the `Task`

keyword is set to AverageBinPlot.

```
Task [RadialDistribution | Histogram | AutoCorrelation | MeanSquareDisplacement | AverageBinPlot]
```

`Task`

- Type
Multiple Choice

- Options
[RadialDistribution, Histogram, AutoCorrelation, MeanSquareDisplacement, AverageBinPlot]

- Description
The analysis task.

Further details need to be specified in the `AverageBinPlot`

block.
If more than one `AverageBinPlot`

block is present in the input, more than one AverageBinPlot will be computed.
The result is printed to output as text, as well as stored in a binary file (analysis.kf).

```
AverageBinPlot
Nbins integer
Property
Axis float_list
Name string
Region string
Vector string
End
Timestep integer
End
```

`AverageBinPlot`

`Nbins`

- Type
Integer

- Default value
10

- Description
Number of bins that are plotted

`Property`

- Type
Block

- Recurring
True

- Description
Property 1

`Axis`

- Type
Float List

- Description
If defined the dot_product along this axis will be taken

`Name`

- Type
String

- Description
Name of the property

`Region`

- Type
String

- Description
Name of the atom region among which the forces will be averaged

`Vector`

- Type
String

- Description
is it a vector? yes/no

`Timestep`

- Type
Integer

- Default value
1

- Description
Timestep used for the plotting

## Viscosity¶

The viscosity can be computed from an equilibrated Molecular Dynamics run as the integral over the pressure tensor autocorrelation function.

\(\eta = \frac{1}{10T} \int_{t=0}^{t=t_{max}} \langle \P_ij(0) \cdot \P_ij{v}(t)) \rangle dt\)

The viscosity is computed if the task `AutoCorrelation`

is selected,
and if in the `AutoCorrelation`

block Viscosity is selected as the `Property`

.

```
$AMSBIN/analysis <<eor
Task AutoCorrelation
AutoCorrelation
Property Viscosity
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.

Again, a subset of atoms can be selected with the sublock `Atoms`

.

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

## Diffusion Coefficient¶

The diffusion coefficient can be computed from a molecular dynamics trajectory in two ways.

As the integral over the velocity autocorrelation function.

As the slope of the mean square displacement of the atomic coordinates.

The latter is more commonly used, as the former requires trajectory information to be stored at very short time intervals. Note that the obtained values for the diffusion coefficients correspond to the temperature of the molecular dynamics simulation.

### 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
- Description

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