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-2025
===============================

* Ionic conductivity

* Mean square displacement (including ionic conductivity) can be
  computed for molecular centers of mass


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


Quickstart Guide
================

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

   #!/bin/sh

   "$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


Trajectory selection
====================

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


Atom Selection
==============

All tasks in the analysis program allow a subset of atoms to be
selected, to be used for the analysis. They can be selected using
element names ("Element"), region names ("Region") or atom indices
("Atom"). In case "Region" is used, the user can overwrite the regions
set in the original simulation in a new System block. The elements in
the system block do have to match the elements in the original
simulation.

Note:

  In the case of a MolecularGun simulation, the input regions will be
  overwritten with the regions from the trajectory file as soon as
  atoms are added to or removed from the system.

   #!/bin/sh

   "$AMSBIN/analysis" << eor

   Task RadialDistribution

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

   TrajectoryInfo
       Trajectory
           KFFilename ams.results/ams.rkf
           Range 1 1000 2
       End
   End

   RadialDistribution
       NBins 1000
       AtomsFrom
           Region Oxygen
       End
       AtomsTo
           Region Oxygen
       End
   End
   eor


Convergence Check
=================

The histogram and radial distribution tasks 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}}}


Atom Selection
--------------

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 directions, 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 or
MDHistory.  Examples are scalar quantities like "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.

It is also possible to select quantities that have scalar or vector
values for each atom, such as "Coords" or "Velocities". If this is
combined with multiple axes, then the number of values along all axes
must be the same. It is, for example, not possible to plot the atom
velocities along one axis, and the temperature per frame along the
other. A subset of atoms for which the atom-wise quantity should be
selected can be provided in the subblock "Atoms". The block can
contain element names (recurring keyword "Element"), region names
(recurring keyword "Region"), or individual atom numbers (recurring
keyword "Atom").

   Atoms
      Atom integer
      Element string
      Region string
   End

"Atoms"
   Type:
      Block

   Description:
      Relevant if variable has a value per atom (e.g. Coords,
      Velocities). This block specifies indices or elements for the
      set of atoms for which the variable is to be read. By default
      all atoms are used.

   "Atom"
      Type:
         Integer

      Recurring:
         True

      Description:
         Atom index.

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

   VecElements
      Index integer
   End

"VecElements"
   Type:
      Block

   Description:
      A set of indices referring to a subset of a vector. If the
      variable to be plotted has non-scalar values per step, then this
      block allows the selection of a subset of vector elements (e.g.
      1 and 2 for the x and y values). Can be used in combination with
      the Atoms block.

   "Index"
      Type:
         Integer

      Recurring:
         True

      Description:
         Element of the property vector.

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
            Atoms
               Atom integer
               Element string
               Region string
            End
            NBins integer
            Range float_list
            Variable string
            VecElements
               Index integer
            End
         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.

         "Atoms"
            Type:
               Block

            Description:
               Relevant if variable has a value per atom (e.g. Coords,
               Velocities). This block specifies indices or elements
               for the set of atoms for which the variable is to be
               read. By default all atoms are used.

            "Atom"
               Type:
                  Integer

               Recurring:
                  True

               Description:
                  Atom index.

            "Element"
               Type:
                  String

               Recurring:
                  True

               Description:
                  Element Symbol Atom.

            "Region"
               Type:
                  String

               Recurring:
                  True

               Description:
                  Region name

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

         "VecElements"
            Type:
               Block

            Description:
               A set of indices referring to a subset of a vector. If
               the variable to be plotted has non-scalar values per
               step, then this block allows the selection of a subset
               of vector elements (e.g. 1 and 2 for the x and y
               values). Can be used in combination with the Atoms
               block.

            "Index"
               Type:
                  Integer

               Recurring:
                  True

               Description:
                  Element of the property vector.


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


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.


Standard Deviation
~~~~~~~~~~~~~~~~~~

By setting the keyword "ComputeStandardDeviation" (by default this is
off), the standard deviation \sigma of the un-normalized
autocorrelation function is added to the plot. This value takes into
account the variation of all contributing values, both from the time
intervals and from the different contributing atoms. The computation
of the standard deviation increases the analysis time by approximately
a factor two.

   AutoCorrelation
      ComputeStandardDeviation Yes/No
   End

"AutoCorrelation"
   Type:
      Block

   Recurring:
      True

   Description:
      All input related to auto correlation functions.

   "ComputeStandardDeviation"
      Type:
         Bool

      Default value:
         No

      Description:
         Compute the standard deviation of the autocorrelation
         function. This will slow down the calculation by
         approximately a factor 2. The standard deviation takes into
         account the variance introduced by both the simulation time
         and by the contributing atoms.


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
      ComputeStandardDeviation Yes/No
      DataReading [Auto | AtOnce | BlockWise]
      LimitMaxFrame Yes/No
      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.

   "ComputeStandardDeviation"
      Type:
         Bool

      Default value:
         No

      Description:
         Compute the standard deviation of the autocorrelation
         function. This will slow down the calculation by
         approximately a factor 2. The standard deviation takes into
         account the variance introduced by both the simulation time
         and by the contributing atoms.

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

   "LimitMaxFrame"
      Type:
         Bool

      Default value:
         Yes

      Description:
         If set to true (default), the program will exit if MaxFrame
         is selected larger than half the provided number of frames.
         This limit is only active is UseAllValues is set to false.

   "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


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


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


Standard Deviation
~~~~~~~~~~~~~~~~~~

By setting the keyword "ComputeStandardDeviation" (by default this is
off), the standard deviation \sigma of the mean square displacement is
added to the plot. This value takes into account the variation of all
contributing values, both from the time intervals and from the
different contributing atoms. The computation of the standard
deviation increases the analysis time by approximately a factor two.

   MeanSquareDisplacement
      ComputeStandardDeviation Yes/No
   End

"MeanSquareDisplacement"
   Type:
      Block

   Recurring:
      True

   Description:
      All input related to mean square displacement functions.

   "ComputeStandardDeviation"
      Type:
         Bool

      Default value:
         No

      Description:
         Compute the standard deviation of the mean square
         displacement. This will slow down the calculation by
         approximately a factor 2. The standard deviation takes into
         account the variance introduced by both the simulation time
         and by the contributing atoms.


Slope
~~~~~

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 \beta is automatically
computed using linear regression.

\beta = \frac{\sum_i(t_i - \overline{t})(MSD_i - \overline{MSD})}
{\sum_i(t_i - \overline{t})^2}

Here i runs over all contributions in the linear regime of the MSD
graph.

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. This is done by computing the slope from t_0 to t_m for
increasing t_0 (from 0 to \frac{1}{2}t_m). At small t_0 values we
expect the linear fit to keep improving, so the correlation
coefficient r_{xy} should increase. At the switching point, where the
MSD becomes linear, r_{xy} will no longer improve, so we select the
t_0 that results in a smaller r_{xy} than the previous one.

r_{xy} = \frac{\sum_i(t_i - \overline{t})(MSD_i - \overline{MSD})} {
\sqrt{\sum_i(t_i - \overline{t})^2 \sum_i(MSD_i - \overline{MSD})^2} }

To allow the user to determine if the linear regime has been
sufficiently sampled, the slope of MSD(t) as a function of en end time
t_e is computed as well. If the slope does not converge to a stable
value, the user should select a larger value of t_m or continue the
molecular dynamics simulation for a longer time.

The convergence plot of the slope takes the MSD data, and computes
slope values over intervals (t_0,t_e), where t_e increases from t_0 to
t_m. If the slope is to represent the diffusion coefficient (see
DiffusionFromMSD), then the slope \beta(t_e) of the MSD up to t_e is
divided by the number of dimensions (usually 3), and multiplied by a
factor \frac{1}{2}.

D(t_e) = \frac{1}{2d} \beta(t_e)


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

      "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 false, 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 Conductivity. The bonds from the
            MD input system are used.


Keywords
--------

The "MeanSquareDisplacement" block holds the following keywords.

   MeanSquareDisplacement
      Atoms
         Atom integer
         Element string
         Region string
      End
      ComputeStandardDeviation Yes/No
      DataReading [Auto | AtOnce | BlockWise]
      LimitMaxFrame Yes/No
      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.

   "ComputeStandardDeviation"
      Type:
         Bool

      Default value:
         No

      Description:
         Compute the standard deviation of the mean square
         displacement. This will slow down the calculation by
         approximately a factor 2. The standard deviation takes into
         account the variance introduced by both the simulation time
         and by the contributing atoms.

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

   "LimitMaxFrame"
      Type:
         Bool

      Default value:
         Yes

      Description:
         If set to true (default), the program will exit if MaxFrame
         is selected larger than half the provided number of frames.
         This limit is only active is UseAllValues is set to false.

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

      "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 false, 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 Conductivity. 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


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
      Atoms
         Atom integer
         Element string
         Region string
      End
      Nbins integer
      Property
         Axis float_list
         Name [FrictionCoefficient | Viscosity | Velocities | EngineGradients]
      End
      XProperty
         Name [Time | Coords]
         VecElements
            Index integer
         End
      End
   End

"AverageBinPlot"
   "Atoms"
      Type:
         Block

      Description:
         Relevant if Properties are atom dependent. 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.

   "Nbins"
      Type:
         Integer

      Default value:
         10

      Description:
         Number of bins that are plotted

   "Property"
      Type:
         Block

      Description:
         Property to be plotted along the Y-axis

      "Axis"
         Type:
            Float List

         Description:
            If defined the dot_product along this axis will be taken.
            Otherwise, the length of the property vector will be used.

      "Name"
         Type:
            Multiple Choice

         Options:
            [FrictionCoefficient, Viscosity, Velocities,
            EngineGradients]

         Description:
            Name of the property

   "XProperty"
      Type:
         Block

      Description:
         Property to be plotted along the Y-axis

      "Name"
         Type:
            Multiple Choice

         Default value:
            Time

         Options:
            [Time, Coords]

         Description:
            Timestep used for the plotting

      "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

            Default value:
               3

            Description:
               Element of the x_property, in case it is a vector (For
               Coords: 1 for X, 2 for Y, 3 for Z).


Viscosity
=========

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

\eta = \frac{V}{k_BT} \int_{t=0}^{t=t_{max}} \langle P_{ij}(0) \cdot
P_{ij}(t)) \rangle_{i \not= j} 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, 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.

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.

If the pressure tensor was written to the BinLog section during the MD
simulation, then the viscosity can be computed with the keyword
"AutoCorrelation%Property" set to *ViscosityFromBinLog*. This is
equivalent to setting "AutoCorrelation%Property" to *Viscosity*, but
since values are written to the BinLog section at every step, there
will be more data available.


Diffusion Coefficient
=====================

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

1. As the integral over the velocity autocorrelation function.

2. 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:
         [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.


Ionic Conductivity
==================

Ionic conductivity can be computed from atomic charge weighted
diffusion coefficients. In this case, we compute the diffusion
coefficients from the mean square displacement. It can be computed by
selecting the "Property" *Conductivity* in the
"MeanSquareDisplacement" block.

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

The set of atoms for which the conductivity is to be computed can be
selected as described in the Mean Square Displacement section. If
"UseMolecularCentersOfMass" is enabled, then the ionic conductivity is
computed per molecule.

Atomic charges can be provided in 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. If no atomic charges are present there either, then the code
will attempt to read atomic charges from the ams.rkf file.

Note:

  The ionic conductivity should be computed only for simulations with
  constant volume and using constant atomic charges. By default, only
  the ForceField engine keeps the atomic charges constant. Any other
  engine (e.g. REAXFF) varies the atomic charges throughout the
  simulation. In those cases it is highly recommended to provide
  atomic charges in an input system block.

   #!/bin/sh

   "$AMSBIN/analysis" << eor

   Task MeanSquareDisplacement

   System
      Atoms
              O      -0.7113934019       1.6140994026      -0.7075521284 analysis.charge=-0.8
              H      -0.0425480297       2.3360319223      -0.8302756333 analysis.charge=0.4
              H      -0.5643276324       1.0450643077      -1.4556911815 analysis.charge=0.4
              O       1.1006174021      -7.1200763645      -3.6017659954 analysis.charge=-0.8
              H       1.4647560784      -6.1967314913      -3.5522086623 analysis.charge=0.4
              H       1.8229257860       7.2170281524      -3.7443300356 analysis.charge=0.4
      End
   End

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

   MeanSquareDisplacement
       Property Conductivity
       Atoms
           Element O
       End
   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.

The ionic conductivity is computed from the mean square displacement,
using coordinates multiplied by the charges (atomic or molecular).

MSD(t) = \frac{1}{N} \sum_i^N \langle (q_i(0){\bf{r_i}} (0) -
q_i(t){\bf{r_i}}(t))^2 \rangle

Here N is the number of atoms/molecules for which the conductivity is
computed, and i runs over those atoms/molecules. The average runs over
all time-intervals of length t throughout the simulation. If the
charges q_i remain constant throughout the simulation, then they can
be moved outside the average.

MSD(t) = \frac{1}{N} \sum_i^N q_i^2 \langle ({\bf{r_i}}(0) -
{\bf{r_i}}(t))^2 \rangle

The charge weighted diffusion coefficient is defined as the time
derivative of the MSD, which is assumed constant.

D_q = \lim_{t\to\infty} \frac{1}{2dN} \frac{d\sum_i^N q_i^2 \langle
({\bf{r_i}}(0) - {\bf{r_i}}(t))^2 \rangle}{dt}

Here d is the dimensionality of the systems (usually d=3).

The ionic conductivity can now be defined as

\sigma = \frac{N}{Vk_BT}D_q

Here, V is the volume of the system, k_B is the Boltzmann constant,
and T is the average temperature.
