# 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 also used under the hood in **ADFMovie** (**MD Properties** menu bar).

This is an example showing how to compute the oxygen-oxygen radial distribution function of a MD simulation using the analysis utility program:

```
$ADFBIN/analysis <<eor
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) simulaton.
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
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: Two or three values: start frame, end frame, step size.

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¶

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

keyword is set to RadialDistribution.

```
Task [RadialDistribution | Histogram]
```

`Task`

Type: Multiple Choice Options: [RadialDistribution, Histogram] 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 (analysis.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.

## Histogram¶

The Analysis program computes histograms if the `Task`

keyword is set to Histogram.

```
Task [RadialDistribution | Histogram]
```

`Task`

Type: Multiple Choice Options: [RadialDistribution, Histogram] 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 (analysis.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 [True | False]
End
```

`Histogram`

`Normalized`

Type: Bool Default value: False 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.