| |
Handy Routines for Ptraj (AmberTools 1.x)
We have extended the functionality of the ptraj
utility (part of the AmberTools
which are freely downloadable from the Amber web
page.) to include some useful routines. One key feature of the new
diffusion and mean residence time codes is that it allows dynamic selection of
regions (at the moment based on fixed distances between sets of atoms). The
functions in detail are:
- A mostly rewritten diffusion code.
- Mean residence time calculations based on a) the survival function and
b) the "standard" autocorrelation function [1-3].
(Alternative
method described in
selected
responses from the Amber mail reflector).
- Calculation of the density along a coordinate useful
e.g. for membrane simulations. (Alternatively, the output of the grid
command can be post-processed to get similar results.)
- Calculation of order parameters for lipids in planar membranes.
- Calculation of the pair-distance distribution function P(r).
- A simple routine to calculate distances to nearest neighbour
boxes.
- Dec 7, 2011: updated to version 1.5, see new
archive file (with
bugfixes 1–18).
For convenience we provide the modifications in a single archive file (currently contains bugfixes
1–18). Extract the files in the
AmberTools/src directory and the files will be written to a new subdirectory
ptraj.mod. Configure the build of AmberTools as usual, change into this
subdirectory and do a make. Alternatively, carry out the modifications
by hand as described below.
The following three files contain the routines and header information.
- actions.c contains the main routines
- transformDiffusion (modified diffusion code)
- transformMRT (mean residence time correlation functions)
- transformDensity (density along coordinate)
- transformOrderParameter (order parameters)
- transformPairDist (pair-distance distribution function P(r))
- transformImgDist (distances to neighbour images)
and the helper functions calculateMSD, inout_MRT,
increaseMem plus some basic vector routines.
These functions need to be copied and pasted to the original actions.c from the
distribution (simply appending the contents of the file provided here should
be fine). The original function transformDiffusion should be commented out
if the new one is to be used.
- actions.h contains the header information
for the routines above and needs to be copied and pasted to the original
actions.h. Please note that the contents of the file provided here goes
before the final #endif /* ACTION_MODULE */ (delete this final line from the
original action.h and append).
- version.h contains version information.
(Older files for AMBER 9 still available
but unsupported:
old9_actions.c and
old9_actions.h.)
The following modifications to the code are also required in addition to
copying and pasting the contents of the files as described above to compile the
code.
- Prototypes for the new functions and new enumeration types have to be
added to action.h.
typedef enum _actionType {
.
.
TRANSFORM_MRT,
TRANSFORM_DENSITY,
TRANSFORM_PAIRDIST,
TRANSFORM_IMGDIST,
TRANSFORM_ORDERPARAMETER
} actionType;
.
.
extern int transform2dRMS(actionInformation *,
double *, double *, double *, double *, int);
/* the next four prototypes are new */
extern int transformMRT(actionInformation *, double *, double *, double *,
double *, int);
extern int transformDensity(actionInformation *,
double *, double *, double *, double *, int);
extern int transformPairDist(actionInformation *, double *, double *, double *,
double *, int);
extern int transformImgDist(actionInformation *, double *, double *, double *,
double *, int);
extern int transformOrderParameter(actionInformation *, double *, double *, double *,
double *, int);
extern actionInformation* ptrajCopyAction(actionInformation**);
.
.
- The new commands have to be added to dispatch.c.
Token ptrajTokenlist[] = {
.
.
{"2drms", 5, -1, TRANSFORM_2DRMS, ptrajSetup},
{"mrt", 3, -1, TRANSFORM_MRT, ptrajSetup},
{"density", 4, -1, TRANSFORM_DENSITY, ptrajSetup},
{"order", 5, -1, TRANSFORM_ORDERPARAMETER, ptrajSetup},
{"pairdist", 5, -1, TRANSFORM_PAIRDIST, ptrajSetup},
{"imgdist", 4, -1, TRANSFORM_IMGDIST, ptrajSetup},
{"go", 2, -1, TRANSFORM_TRANSFORM, ptrajSetup},
.
.
- Also ptraj.c must be updated.
switch ( type ) {
.
.
case TRANSFORM_2DRMS:
action->type = TRANSFORM_2DRMS;
action->fxn = (actionFunction) transform2dRMS;
break;
case TRANSFORM_MRT: /* new */
action->type = TRANSFORM_MRT;
action->fxn = (actionFunction) transformMRT;
break;
case TRANSFORM_DENSITY: /* new */
action->type = TRANSFORM_DENSITY;
action->fxn = (actionFunction) transformDensity;
break;
case TRANSFORM_ORDERPARAMETER: /* new */
action->type = TRANSFORM_ORDERPARAMETER;
action->fxn = (actionFunction) transformOrderParameter;
break;
case TRANSFORM_PAIRDIST: /* new */
action->type = TRANSFORM_PAIRDIST;
action->fxn = (actionFunction) transformPairDist;
break;
case TRANSFORM_IMGDIST: /* new */
action->type = TRANSFORM_IMGDIST;
action->fxn = (actionFunction) transformImgDist;
break;
case TRANSFORM_TRANSFORM:
case TRANSFORM_NOOP:
.
.
- Finally, compile as usual.
| mrt |
out <filename>[autocorr <filename>tcorr <time>toffset
<time>] |
|
[lower <dist>] [upper <dist>] [time <t>] [tstar
<t>] [noimage] |
|
[solvent <mask>| solute <mask>]
(siteatoms <mask>| [onemol]
<sitemask>| |
|
<sitemask1> ... <sitemaskN>) |
Square brackets mean optional parameters. In the following defaults are
given in parenthesis.
| out |
Output file: time vs. MRT correlation function for each mask. |
| autocorr |
Optional autocorrelation code. (off) |
| tcorr |
(autocorr) Maximum correlation time = size of windows in ps. (400.0 ps) |
| toffset |
(autocorr) Time step/gap between windows (code computes parallel
windows to improve statistics). (10.0 ps) |
| solvent |
Optional solvent mask which overrides the internal definition of
solvent. May also be changed with the solvent command (see ptraj
manual). |
| solute |
Optional mask representing a single site (centre of mass). |
| lower |
Smaller distance from reference point(s). (0.01 Å) |
| upper |
Larger distance from reference point(s). (3.5 Å) |
| time |
Time step in the trajectory. (1.0 ps) |
| tstar |
t* to account for short term excursions (see Impey paper[1]). (0.0 ps) |
The reference point(s)/surface can be set in the following ways at the
moment.
| siteatoms <mask> |
Atoms in the mask will be considered as individual reference points,
i.e. if the mask expands to 100 atoms then 100 MRTs will be computed,
useful for "quick scanning" of many sites. |
| onemole <sitemask> |
All atoms in the mask constitute a single reference site, i.e. if the
mask contains all atoms of a protein then a single MRT around that
protein will be calculated. IMPORTANT: may be very slow!!! |
| <sitemask1>... |
After all named parameters have been processed each remaining parameter
will be interpreted as a single reference site. If it contains several
atoms the centre of mass will be taken as reference point. |
| diffusion |
mask <mask>[out <file>] [time <time per frame>] |
|
[mask2 <mask> lower <distance> upper <distance> nwout <file>] |
|
[avout <file>] [x|y|z|xy|xz|yz|xyz] [distances] [com] |
Square brackets mean optional parameters. In the following defaults are
given in parenthesis. A simple Perl script is provided
to compute window averages.
| mask |
Atoms for which MSDs will be computed. |
| out |
Output file: time vs. MSD. |
| time |
Time step in the trajectory. (1.0 ps) |
| mask2 |
Compute MSDs only within the lower and upper limit of mask2.
IMPORTANT: may be very slow!!! |
| lower |
Smaller distance from reference point(s). (0.01 Å) |
| upper |
Larger distance from reference point(s). (3.5 Å) |
| nwout |
Output file containing number of water molecules in the chosen
region, see mask2. (off) |
| avout |
Output file containing average distances. (off) |
| x|y|z|xy|xz|yz|xyz |
Computation of the mean square displacement in the chosen dimension. (xyz) |
| distances |
Dump un-imaged distances. By default only averages are output. (off) |
| com |
Calculate MSD for centre of mass. (off) |
The atom selections must be centred at the origin before calculating the density. See also example input file.
| density |
out <filename> [delta <resolution>] [x|y|z] |
|
[number|mass|charge|electron] [efile <filename>] <mask1> ... <maskN> |
Square brackets mean optional parameters. In the following defaults are
given in parenthesis.
| out |
Output file for histogram: relative distances vs. densities for each mask. |
| delta |
Resolution, i.e. determines number of slices. (0.25 Å) |
| x|y|z |
Coordinate for density calculation. (z) |
| number|mass|charge|electron |
Number, mass, partial charge (q) or electron
(Ne - q) density. To convert the electron density to
e-/Å3 divide by the (average) area spanned
by the other two dimensions. (number) |
| efile |
Optional file with electron numbers Ne for each atom
type. If this file is not supplied the command assumes that atom types
starting with 'H', 'C', 'N', 'O', 'P' and 'S' are the actual elements.
The format of the file is as follows. Comments are lines starting with
'#' or empty lines. All other lines must contain the atom type followed
by an integer number for the electron number. Entries must be separated
by spaces or '='. See example input file. |
| mask1 ... maskN |
Arbitrary number of masks for atom selection, the output will contain a
column for each mask. |
The order parameters Sx, Sy, Sz and
|SCD| are calculated. Sz is the vector joining carbons
Cn-1 and Cn+1, Sx the vector normal to the
Cn-1–Cn and Cn–Cn+1
plane and Sy is the third axis in the molecular coordinate system.
The order parameter is then calculated from Sc =
0.5<3cos2θ> - 1, where θ is the angle to the
chosen reference axis. See example input file.
| order |
out <filename> [x|y|z] [unsat <mask>] |
[taildist <filename> [delta <resolution>] tailstart <mask> tailend <mask>] |
<mask0> ... <maskN> |
| out |
Output file for order parameters: Sx, Sy,
Sz (each succeeded by the standard deviation), and two
estimates for the deuterium-order parameter |SCD| =
0.5Sz and |SCD| = -(2Sx +
Sy)/3. If scd is set then the order parameter directly
computed from the C-H vectors is output. |
| x|y|z |
Reference axis. (z) |
| unsat |
Mask for unsaturated bonds. Sz is calculated for vector
Cn–Cn+1. |
| scd |
Calculate the deuterium-order parameter |SCD| directly from
the C-H vectors. |
| taildist |
Optional output file for end-to-end distances. |
| delta |
Optional resolution for taildist. (0.1) |
| tailstart |
Mask for the start of the tail. Must be given if taildist. |
| tailend |
Mask for the end of the tail. Must be given if taildist. |
| mask0 ... maskN |
Masks for each group in the lipid chain. Carbons must be
given in bonding order. If scd the masks must be made up of C-H-H
triples, hence hydrogens to double bonds must be enumerated twice while
methyl groups require an additional mask which will also create two
entries in the output. See example input file. |
| pairdist |
out <filename> mask <mask> [delta <resolution>] |
Square brackets mean optional parameters. In the following defaults are
given in parenthesis.
| out |
Output file for histogram: distance, P(r), s(P(r)). |
| mask |
Atoms for which distances should be computed. |
| delta |
Resolution. (0.1 Å) |
| imgdist |
out <filename> mask <mask> threshold <distance> |
| out |
Output file listing distances shorter than the threshold. |
| mask |
Atoms for which distances should be computed. |
| threshold |
Report all distances smaller than this value. |
If you have any questions regarding the code, installation/compilation, ideas for further
development, contributions, etc. feel free to contact Hannes
Loeffler.
- R. W. Impey, P. A. Madden and I. R. McDonald, J. Phys. Chem. 87(1983),
5071.
- A. E. Garcia and L. Stiller, J. Comput. Chem. 14(1993), 1396.
- C. Schröder, T. Rudas, S. Boresch, O. Steinhauser,
J. Chem. Phys. 124(2006), 234907.
|
