rotdif

Calculate rotational diffusion using given rotation matrices (from e.g. rms).

rotdif [outfile <outfilename>] [usefft]

Options for generating random vectors:

[nvecs <nvecs>] [rvecin <randvecIn>] [rseed <random seed>] [rvecout <randvecOut>] [rmatrix <set name> [rmout <rmOut>]]

Options for calculating vector time correlation functions:

[order <olegendre>] [ncorr <ncorr>] [corrout <corrOut>]

*** The options below only apply if ’usefft’ IS NOT specified. ***

Options for calculating local effective D, small anisotropy:
[deffout <deffOut>] [itmax <itmax>] [tol <tolerance>] [d0 <d0>][nmesh <NmeshPoints>] dt <tfac> [ti <ti>] tf <tf>

Options for calculating D with full anisotropy:
[amoeba_tol <tolerance>] [amoeba_itmax <iterations>][amoeba_nsearch <n>] [scalesimplex <scale>] [gridsearch]

*** The options below only apply if ’usefft’ IS specified. ***

Options for curve-fitting:
[fit_tol <tolerance>] [fit_itmax <max # iterations>]

outfile <outfilename> File to write all output from rotdif command to.

Options for generating random vectors:
nvecs <nvecs> Number of random vectors to generate (default 1000).
rvecin <randvecIn> File to read random vectors from (format is 1 per line, 4 columns, <#> <VX> <VY> <VZ>).
rseed <random seed> Seed for random number generator (default 80531). Specify -1 to use wallclock time.
rvecout <randvecOut> File to write random vectors to (format is 1 per line, 4 columns, <#> <VX> <VY> <VZ>).
rmatrix <set name> Data set to read rotation matrices from. Rotation matrices will be used to rotate random vectors.
rmout <rmOut> Write rotation matrices to file, 1 per line, frame # followed by matrix in row-major order.

Options for calculating vector time correlation functions:
order <olegendre> The order of Legendre polynomials to use when calculating vector time correlation functions (default 2).
ncorr <ncorr> Maximum length of time correlation functions in frames. If this is not specified it will be set to (tf – ti) / dt (recommended).
corrout <corrOut> If specified write vector time correlation functions to <corrOut>.X with format: <Time> <Px>

Options for calculating local effective D, small anisotropy:
deffout <deffOut> File to write out local effective diffusion constants determined in the limit of small anisotropy.
itmax <itmax> Maximum number of iterations to determine each local effective diffusion constant (small anisotropy) assuming fit to single exponential form (default 500).
tol <tolerance> Tolerance for determining local effective diffusion constant (small anisotropy) assuming fit to single exponential form (default 1E-6).
d0 <d0> Initial guess for small anisotropy diffusion constant in radians^2/ns (default 0.03).
nmesh <NmeshPoints> Number of points per frame to use when creating cubic-splined-smoothed forms of vector time correlation curves (default 2).
dt <tfac> Time interval between frames (used in integrating vector time correlation curves) in ns.
ti <ti> Initial time value in ns for integrating the time correlation functions (default 0.0).
tf <tf> Final time value in ns for integrating the time correlation functions. It is recommended this be less than the maximum simulation time since the tails of time correlation functions tend to be noisy.

Options for calculating D with full anisotropy:
amoeba_tol <tolerance> Tolerance for downhill-simplex minimizer (default 1E-7).
amoeba_itmax <iterations> Number of iterations to run downhill-simplex minimizer (default 10000).
amoeba_nsearch <n> Number of searches to perform with downhill-simplex minimizer (default 1).
scalesimplex <scale> Factor to use when scaling simplexes (default 0.5).
gridsearch If specified, perform a brute-force grid search to attempt to find a better solution for diffusion tensor with full anisotropy (may be
expensive).

Evaluate rotational diffusion properties of a molecule over a trajectory according to an expanded version of the procedure laid out by Wong & Case. Briefly, random vectors (representing the orientation of the molecule) are rotated according to rotation matrices obtained from an RMS fit to a reference structure (typically an averaged structure). For each random vector the time correlation function of the rotated vector is calculated using Legendre polynomials of the specified order. The integral over this time correlation function (which may be smoothed using cubic splines to improve the integration) is then used to find the effective diffusion constant (D) in the limit of small anisotropy. Then, using each calculated D, the diffusion tensor is determined with full anisotropy. Finally, a downhill simplex minimizer is used to optimize D with full anisotropy; (this last step is not described in the original paper).

Rotation matrices are generated via an RMS fit to a reference structure. It is recommended that the RMS fit be done to an average structure. These rotation matrices are used to rotate each random vector M times (where M is the total number of frames), which creates a time series for each random vector. The time correlation functions are calculated for each random vector time series using Legendre polynomials of the specified order (default 2). The maximum length of the correlation function (or lag) can be specified by ncorr (in frames). If ncorr is not specified it will be set internally based on the specified values of ti, tf, and dt; this is recommended. Note that if ncorr is specified it should be set to a number less than the total number of frames since noise in time correlation functions increases as ncorr approaches the # of frames.

The integration over the correlation function is from ti (in whatever units are used of dt, generally ns; 0.0 ns if not specified) to tf (same units as ti), with the time between frames specified by dt; the final time should be less than the total simulation time (see example below).

The relative size of the mesh used with cubic spline interpolation for integration is controlled by nmesh (size of the mesh is ncorr points * nmesh); nmesh = 1 means no interpolation, default is 2. The iterative solver for effective value of the diffusion constant from the correlation functions is controlled by itmax, tol, and d0, where itmax specifies the number of iterations to perform (default 500), tol specifies the tolerance (default 1E-6), and d0
specifies the initial guess for the diffusion constant in radians^2 / ns (default 0.03). Effective diffusion constants for each random vector can be written out to a file specified by deffout. Results are printed to the file specified by outfile. Details on the Q and D tensors are given, as well as observed and calculated tau for each random vector.

First, results are printed for analysis in the limit of small anisotropy. Next, results are printed for analysis with full anisotropy. The results of the full anisotropic calculation are first given using results from the small anisotropic analysis as an initial guess, followed by the final results after minimization using the downhill simplex (amoeba) minimizer.

Example
There are two important things to keep in mind when using rotdif analysis:

1. When calculating any kind of diffusive property it is best to simulate in the microcanonical (NVE) ensemble with a shorter time step and increased SHAKE tolerance; thermostats and barostats will effect diffusion calculations.
2. Time correlation functions become noisier as the length of the function approaches the maximum. Therefore in general one should choose parameters for the time correlation function that are much shorter than the total simulation length.

For example, given a trajectory ’mdcrd.nc’ containing 10000 frames with a total simulation time of 200 ns (so the time between frames is 0.02 ns), to calculate rotational diffusion using 100 vectors using rotation matrices generated via an RMS fit to ’avgstruct.pdb’, computing and integrating the time correlation function for each vector from 0 to 5 ns (1/40th of the simulation), and writing out the effective diffusion constants and results to
’deffs.dat’ and ’rotdif.out’ respectively:

reference avgstruct.pdb [avg]
rms R0 @CA,C,N,O ref [avg] savematrices
trajin mdcrd.nc
rotdif nvecs 100 rmatrix R0[RM] \
ti 0.0 tf 5.0 dt 0.02 deffout deffs.dat \
outfile rotdif.out