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 AMSmovie (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:

$AMSBIN/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) 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
   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 Function (RDF)

The Analysis tool computes radial distribution functions \(g(r)\) if the Task keyword is set to RadialDistribution.

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

AutoCorrelation
   Atoms
      Atom integer
      Element string
   End
   DataReading [Auto | AtOnce | BlockWise]
   InputValues
      Values float_list
   End
   MaxStep integer
   NPointsHighestFreq integer
   Normalized Yes/No
   Property [Velocities | DipoleMomentFromCharges | InputValues | DiffusionCoefficient | DipoleDerivativeFromCharges | PressureTensor]
   TimeStep float
   UseAllValues Yes/No
   UseTimeDerivative
      Enabled Yes/No
      ProjectOutRotations Yes/No
   End
   VecElements
      Index integer
   End
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 handledt 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 is Property is set to InputValues. All input values (a vector on each line).
Values
Type:Float List
Recurring:True
Description:The values at each step (on a single line)
MaxStep
Type:Integer
Description:The maximum interval of the autocorrelation. 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). A higher number selected here, will result in a lower maximum frequency returned by the program. The default value is 4. and the lowest possible value (spectrum up to highest possible frequency) is 2.
Normalized
Type:Bool
Default value:Yes
Description:Determines if the ACF is normalized. Keyword is overruled (set to False) if Property is set to DiffusionCoefficient.
Property
Type:Multiple Choice
Default value:DipoleMomentFromCharges
Options:[Velocities, DipoleMomentFromCharges, InputValues, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor]
Description:Compute the ACF either from velocities (from rkf), the dipole moment (from atomic charges in rkf), or from values specified in input. If DiffusionCoefficient is selected the unnormalized velocity autocorrelation function is computed and integrated.
TimeStep
Type:Float
Description:Relevant if Property is set to InputValues. The time separating the entries (in fs). If Property is set to Velocities, DipoleMomentFromCharges, or DiffusionCoefficient, then the property can be obtained from an RKF file, and the timestep is read from the RKF file as well. The read value then overrides this keyword.
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.
ProjectOutRotations
Type:Bool
Default value:No
Description:Take the rotations out of the time derivative.
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

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. If applicable, the average also runs over all possible contributions to \(\textbf{A}\) at each simulation timestep. 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, the power spectrum matches the IR spectrum.

Options

Autocorrelation functions can be computed for different simulation properties: 1) Dipole moments from atomic charges 2) Velocities 3) User provided values.

AutoCorrelation
   Property [Velocities | DipoleMomentFromCharges | InputValues | DiffusionCoefficient | DipoleDerivativeFromCharges | PressureTensor]
End
AutoCorrelation
Type:Block
Recurring:True
Description:All input related to auto correlation functions.
Property
Type:Multiple Choice
Default value:DipoleMomentFromCharges
Options:[Velocities, DipoleMomentFromCharges, InputValues, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor]
Description:Compute the ACF either from velocities (from rkf), the dipole moment (from atomic charges in rkf), or from values specified in input. If DiffusionCoefficient is selected the unnormalized velocity autocorrelation function is computed and integrated.

With the keyword Normalized a normalized ACF is computed, and with the keyword MaxStep the number of values \(n\) in the autocorrelation function (\(t = [0,t_{1},t_{2},....,t_{n}]\)) can be set. The default value is half of the total number of simulation steps 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.

Diffusion Coefficient

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

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

The factor \(\frac{1}{3}\) corrects for the dimension of the system, which we assume to be always 3.

The diffusion coefficient is computed if the task AutoCorrelation is selected, and if in the AutoCorrelation block DiffusionCoefficient is selected as the Property.

$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:DipoleMomentFromCharges
Options:[Velocities, DipoleMomentFromCharges, InputValues, DiffusionCoefficient, DipoleDerivativeFromCharges, PressureTensor]
Description:Compute the ACF either from velocities (from rkf), the dipole moment (from atomic charges in rkf), or from values specified in input. If DiffusionCoefficient is selected the unnormalized velocity autocorrelation function is computed and integrated.

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

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