| |
Handy Routines for Ptraj (AmberTools 1.x and AMBER 9)
We (contact Hannes
Loeffler for questions) have extended the functionality of the ptraj
utility (AmberTools 1.x/AMBER9) to include some useful routines. A key feature
of the diffusion and the mean residence time code 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
reponses from the Amber mail reflector).
- Calculation of the density along a coordiante useful for
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.
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 orginial 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 orginial
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.
All files for Amber10/AmberTools 1.x. Older files for AMBER 9 still available
but unsupported:
old9_actions.c and
old9_actions.h.
Updates
- Feb 18, 2010: The order command now optionally computes tail end-to-end
distances.
- Feb 16, 2010: Order parameter calculations for lipids in planar membranes
implemented.
- Feb 9, 2010: The diffusion command has now 'convenience' switches to
calculate the MSD in only one, two, or all three dimensions. Data will be
written in the last column of the output file. A simple
Perl script is provided to compute window
averages.
- Feb 1, 2010: The density command now optionally reads a file with electron
numbers.
Code Modifications
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 funcions 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.
Usage
MRT
| 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 (center 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 intepreted as a single reference site. If it contains several
atoms the center of mass will be taken as reference point. |
Diffusion
| 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 computated. |
| 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 center of mass. (off) |
Density
The atom selections must be centered 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. (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. |
Order Parameter
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. |
Pair-Distance Distribution Function
| 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 Å) |
Distance to Images
| 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. |
Contact
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.
|