Multi-body Dynamics: Overview
In multi-body dynamics, aggregates of atoms are gathered into
"bodies". For a dynamics run, the system comprises one or more bodies
and zero or more atoms which are not part of any body. By gathering
the atoms in this way, the total number of variables in the system is
considerably reduced which is expected to significantly improve the
computational performance. Furthermore because such a simulation aims
to reproduce the characteristic (i.e. low-frequency) motion of the
system, relatively long time steps are possible. The final advantage
of this scheme is that bond-lengths may be explicitly constrained
(between bodies and in the atomistic regions) in a computationally
efficient manner.
There are two steps needed before starting a dynamics run
Identifying the bodies and the atomistic regions and
Generating (or loading pre-calculated) modes for each of the
bodies.
All of the standard CHARMM output and analysis mechanisms work with
multi-body dynamics. One new file format is used (to store computed
mode shapes).
The MBOND command is used for setting up the system and controlling
its activation. It can be used either as a single line command (mainly for
control and status reports) or as an opening for a command block (for
setting up the system substructuring and mode assignments).
All the single line commands can be also used from within the command block.
The single line commands are:
MBOND { [ON ] }
{ [OFF ] }
{ [CLEAr ] }
{ [STATus { [SHORt ]} ]}
{ [LONG [BODY char*4] ]}
{ [VLONg [BODY char*4] ]}
{ [ELECtrostatics]}
{ [ATOM { [SCBOnd real] } ]}
{ [NOBOnd ] }
{ [STIFF ] }
The command block has the following syntax:
MBOND [NBODY integer]
SUBStructure
.
.
END
MODES
.
.
.
END
END
NBODY is used to specify the MAXIMUM number of bodies to be supported.
It only needs to be specified once per CHARMM run. To minimize memory
usage, it is useful to keep this number small. If you change the
value of NBODY within the same CHARMM run, then all previously defined
structures are erased. NBODY must be a positive number.
Within the MBOND command block all of CHARMM's miscellaneous commands
can be used. This allows for command-line control: STREam files, GOTO,
LABEl, IF etc.
* Menu:
* Substructuring:: How to define bodies
* Modes:: Description of the keywords and options
* Other:: One line MBOND commands
* Dynamics:: Recommended input options and values
* Output:: Output from a dynamics run
To identify collections of atoms to be grouped into bodies, use the
SUBSTRUCTURING sub-comands of MBOND. A body is identified using the
SELECTION commands of CHARMM. Generally bodies will be contiguous
atoms, though there is no restriction that they be so (e.g. you could
define a set of water molecules as a body).
Bodies can have any name of four characters, except 'ALL '.
By default, the bonds between bodies and all bonds in the atomistic
regions will be contrained. To change this behavior, use the HINGE
subcommand (currently only BOND and OFF are available).
It is possible to specify that all atoms in the system should be
modelled as separate bodies i.e. do an atomistic simulation using the
multi-body formulation. This is the function of the PARTICLE
sub-command.
The syntax for the SUBSTRUCTURE sub-command is:
MBOND [NBODY integer]
SUBStructure
BODY char*4 SELEct (...) END [SHOW]
HINGE { [BOND]} [ANGLe] [DIHEdral] [OFF]
[ PARTICLE ]
END
END
These options are documented on the next screen.
SUBStructure
Begins the substructuring commands. When you exit (END
subcommand), the "topology" information is passed to MBOND.
BODY char*4 SELEct (....) END [SHOW]
Define a BODY using the CHARMM regular SELECT facility. The body
does not have to be continuous in terms of atomic numbers.
The "body_name" can be any character*4 string.
SHOW flashes out information about the body just defined.
HINGE { [BOND ]} [ANGLe] [DIHEdral] [OFF]
Define the type of constraints between bodies. Currently only
BOND constraints are available.
options:
BOND = constrain all bonds
(i.e. except those within a body) (default)
ANGLe = constrain hinge angles (currently not available)
DIHEdral = constrain hinge dihedrals (currently not available)
OFF = NO constraints.
PARTICLE
Each atom in the system is to be considered as a separate
body. The MBOND formulation will be used for dynamics
simulations. This is sometimes useful for direct comparison
with CHARMM atomistic simulations.
If you simply select sets of atoms as bodies, then these will be
modelled as rigid bodies by default. To identify a body as flexible,
a set of modes can either be generated or loaded for it. A new file
format (see ) has been created with all necessary information for a
multi-body dynamics simulation. In addtion, modes can be sorted
(either by frequency or by delocalization).
The syntax for the MODE sub-command is:
MBOND [NBODY integer]
MODEs { [KEEP ] } [NMODe integer]
{ [DELEte] }
LOAD char*4 {[NMODe integer] } UNIT integer
{[MODE integer THRU integer]} [SHOW] [NOVEC]
UNLOad [char*4 ]
[ALL ]
USEM {[char*4]} {[NMODe integer] }
{[ALL]} {[MODE integer THRU integer]}
GENErate char*4 [NMODe integer] [DELO] [SHOW]
MINI minimiz_spec
DIAG {[VACUum]}
{FIXEd}
SORT {FREQ} [SHOW]
{DELOcalization}
WRITE {[NMODe integer] } UNIT integer
{[MODE integer THRU integer]}
{[ALL]}
END
SORT {[char*4]} {[FREQ]}
{[ALL]} {[DELOcalization]}
PED {[char*4]} [Mode-Spec] [Magnitude Spec] [TOL real]
END
These options are documented on the next screen.
MODES { [KEEP ] } [NMODe integer]
{ [DELEte] }
Invokes the flexible-body mode-definition module. All multi-atom
bodies (defined in the substructuring module) are assumed "rigid"
unless explicitly defined as "flexible" in the MODES module (by LOAD,
USEM or GENErate).
options:
KEEP (the default) keeps the mode information (eigenvectors,
frequencies, nominal coordinates) in CHARMM memory after
passing it to MBOND.
DELETE removes this information after interfacing with MBOND.
NMODE sets a default value for the NMODE keyword in the LOAD
command. Default is 20.
LOAD char*4 {NMODe integer } UNIT integer
{MODE integer THRU integer} [SHOW] [NOVEC]
Load the mode-information of the body with "body_name" character*4
from the formatted file opened in UNIT (UNFOrmatted option not
available).
The LOAD command reads the body-based nominal coordinates, the
"modal stiffness" matrix (F^-1 K F), and either the first NMODE
eigenvectors or the eigenvectors in the range "MODE ... THRU ...".
If this body is already defined as "flexible" a call to LOAD
will overwrite the previous mode-information (frees the
currently allocated space and reallocates a new one).
options:
SHOW echoes the information read from the file.
NOVEC a flag that allows for reading the modal-stiffness
matrix and properties without reading the eigenvectors.
[An appropriate unit has to be opened before using the LOAD
command.]
UNLOad {[char*4]}
{[ALL]}
Frees the space where the mode information for the body with
"body_name" character*4 is stored, effectively turning it to a
"rigid" body.
Can be used for changing a "flexible" body into a "rigid" body, and for
explicitly "cleaning up" before changing the character of an already
defined "flexible" body.
UNLOAD is identical to: "LOAD char*4 NMODE 0", and is not
necessary when LOADing or GENErating new modes for a "flexible" body.
UNLOAD ALL removes from memory all mode information for all bodies.
USEM {[char*4]} {NMODE no_of_modes }
{[ALL]} {MODE ifirst THRU ilast}
Use the first NMODE modes or MODE/THRU modes in the dynamics
calculation. The default is all the loaded modes. The number
of modes used cannot exceed the number of modes loaded by
LOAD (or GENErate).
The USEM command does not free memory so that the loaded modes can be
reached again if needed.
USEM ALL ... will use the specified mode numbers for all bodies.
GENErate char*4 {[NMODe integer]} [DELO] [SHOW]
{ALL}
MINI minimiz_spec
DIAG {[VACUum]}
{FIXE}
SORT {FREQ} [SHOW]
{DELOcalization}
WRITE {[NMODe integer] } UNIT integer
{[MODE integer THRU integer]}
END
Generate calculates modes for the specified body. At most,
NMODes (default 20) are generated. ALL modes for that body
may be determined if sufficient memory is available. Any
previously loaded (or generated) modes for that body are cleared.
The DELOcalization is a measure of how localized a particular
mode is. It is the ratio of the fourth moment of the
displacement along the mode vector and the second-moment
squared, i.e.,
Sum_over_atoms (r^4/(r^2)^2).
In very localized modes this ratio is close to 1,
in very delocalized modes it is close to 0.
[NOTE, at very high frequencies multiple occurrences of
localized modes (e.g., C-H stretching), will have a small
"DELOcalization" value.]
Selecting this keyword computes the DELOcalization of each
generated mode and stores it as the second property of that
mode (the frequency is the first).
All existing MINImization options are suppored (by the
existing CHARMM minimization module). All atoms which do not
belong to the specified body are fixed for the duration of the
GENErate command (all atoms are restored to their initial
coordinates on exit).
DIAG supports two different methods. In both, it computes
the Hessian for atoms in the selected body, and diagonalizes
it. The 6 translational/rotational modes are discarded.
FIXED environment includes the interactions between atoms in
this body and other atoms in the system when computing the
Hessian. VACUUM (the default) ignores these external
interactions .
By default, generated modes are SORTed by frequency (most
negative to most positive). You can change the order by
SORTing according to the delocalization factor (and re-sorting
by frequency).
WRITe saves the range of selected modes to the specified UNIT
in the format which LOAD uses. They are saved in the current
order (i.e. sorted by frequency or by delocalization).
All of CHARMMs control and command parsing options are
supported within the generate block. In particular, this
allows for different stream files to be created and used to
generate modes for each of the bodies in a system.
SORT {[char*4]} {[FREQ]}
{[ALL]} {[DELOcalization]}
Modes can be sorted at other times other than GENEration. For
example if modes are computed using some other program, then
loaded into CHARMM, they can be still be SORTed. The possible
keys are FREQuency and DELOcalization. If the DELOcalization
has not already been computed, it will be automatically.
PED {[char*4]} [Mode-Spec] [Magnitude Spec] [TOL real]
For a given magnitude specification, it computes the expectation
value for the energy contribution change for each internal
coordinate term (bond, angle, dihedral, and improper dihedral)
and prints that term if the fluctuation is greater than the
tolerance (default TOL 0.0001).
[Mode spec] is the usual choice of either NMODES N (i.e. the first N
modes) or MODE M THRU N. See *note Normal modes:(vibran.doc)
for a description of [Magnitude Spec].
Modes must have been loaded for the specified body (either by LOAD or
GENERATE). Right now, only VACUUM modes are supported (i.e.
contributions from atoms outside the body are NOT included).
MBOND { [ON ] }
{ [OFF ] }
{ [CLEAr ] }
{ [STATus { [SHORt ]} ]}
{ [LONG [BODY char*4] ]}
{ [VLONg [BODY char*4] ]}
{ [ELECtrostatics]}
{ [ATOM { [SCBOnd real ]} ]}
{ [NOBOnd ]}
{ [STIFfness]}
All these single line commands can be also used from within the MBOND
command block (but not from within the SUBS and MODE sub-blocks
ON Activate MBOND and take the substructuring into
account when calculating energy etc. This is the
default after setting up the substructuring (in the
MBOND command block).
OFF Un-activate MBOND (all data structures are kept intact).
CLEAR Un-activate MBOND and remove all related data
structures from memory. Also resets ELEC, ATOM etc.
STATus Printout MBOND status report (including number of
bodies, size, mode information, etc). For the LONG and
VLONG options one can specify a specific BODY name,
the default is all bodies. You can also specify ALL.
option:
SHORt - give global status
LONG (default) - adds body information
VLONg - adds mode information
ELEC Use MBOND's multipole expansion algorithm to calculate
electrostatic interactions (skip CHARMM's electrostatic
energy evaluations). The default is to use CHARMM's
electrostatics (!). CURRENTLY DISABLED.
ATOM Intra-body forces can be computed either explicitly or
by employing a linear force approximation (using the
STIFFNESS MATRIX). The default is the full atomistic
calculation specified by the ATOM option.
To compensate for Cartesian curvelinear effects the SCBOND
option may be invoked. This will scale all the intra-body
bond-energy terms by the factor <real>. The scaling factor
should be in the range [0.0,1.0] (where 0.0 is identical to
NOBOnd, and 1.0 is identical to the regulat MBOND ATOM option).
Inputs outside of this range are reset to 1.0).
The NOBOND keyword will skip intra-body bond-energy terms
altogether (equivalent to SCBOnd 0.0).
STIFf Use a linear force approximation to compute the
intra-body forces (using the STIFFNESS MATRIX).
This is the opposite of the ATOM option.
Multi body dynamics are invoked by using the option MBOND on the regular CHARMM DYNAmics command. This enables a number of keywords particular to MBOND and disables some which are not. In general, however, all the standard DYNAmics keywords are supported unless specifically mentioned here. Bodies must have been defined prior to invoking multi-body dynamics (see MBOND SUBStructure for details). If no modes have been loaded (or generated) for the bodies, they are assumed to be rigid. Control is passed from CHARMM to the MBO(N)D package until the end of the dynamics simulation. During the run, CHARMM is utilized to evaluate the forces. It also implements some of the controlling logic for heating, equilibration, thermostats etc. Only velocity scaling is supported for heating and equilibration protocols. However the initial velocities may be assigned using any of the CHARMM options (including using a RESTART file). The thermostat which is available is a simple Berendsen method. The general protocol is to heat a system using an atomistic simulation, then equilibrate it and generate modes. When the changeover is made to a multi-body simulation, the system needs to be equilibrated once more before production dynamics can begin. Multiple Time Scales are supported in MBOND. There is a special MTS keyword within the MBOND command to enable/control Multiple Time Scale dynamics. See *note MTS:(mbond.doc)Mbond for details. The following CHARMM features are NOT supported in multi-body simulations: SHAKE, CONSTANT PRESSURE, NOSE, LANGEVIN. In addition, the integrators the MBOND supports are Lobatto, Velocity Verlet, Velocity Central Difference and 4th Order Runge Kutta. There is no LEAPFROG or 4D Verlet. These MBOND dynamics options are documented on the next screen.
mbond-spec::= [[ {LOBAtto} ] ]
[[ VVER ] [NCYC integer] [VTOL real] ]
[[ VCD ] ]
[[ RK ] ]
[ MBPRLev integer ] [ VECFreq integer ] [ TFRQ integer ]
[ IUVEctor integer] [ ATOM ]
The following four integrator choices are available for multi-body dynamics
LOBATTO The default integrator.
VVER Iterated velocity-verlet integration.
VCD Verlet Central difference
RK 4th Order Runge-Kutta
NCYC The maximum number of iterations in any verlet style
integration step (default 1)
VTOL Convergence criterion for verlet style integrators
in multi-body dynamics (default 1E-10)
ATOM If this is not specified (or has not been
previously), the linear force approximation is used
within bodies. If specified, the full force
computation is performed.
MBPRlev 0 Print level for multi-body dynamics.
MBPRlev = -1 : No printout from within MBOND, only the
data passed to CHARMM will be printed out
(every NPRInt timesteps).
MBPRlev = 0 : Basic printout from within MBOND.
MBPRlev = 1 : Print regular (un-substructured) CHARMM energy
MBPRlev = 3 : Full printout from within MBOND.
MBPRlev = 9 : Extended (debugging) printout from within MBOND.
VECFreq NPRINT Printout frequency for MBOND vector information
(MBPRLEV at least 3).
IUVEctor 6 Unit to write MBOND debugging information.
TFRQ NPRINT Frequency for writing out the temperatures.
Multiple Time Scale dynamics are a natural complement to body based
dynamics. MBOND currently supports segregation of the computation of
non-bonded forces into different timescales. By default, the number
of bins and the relative ratios are computed automatically (although
these may also be specified exactly).
Multiple time step dynamics in MBOND is controlled by the MTS keyword
to the MBOND command. This must be invoked BEFORE dynamics is begun.
Once you have selected options in this way they remain in effect for
all subsequent dynamics runs.
MBOND
MTS
MASS [CALIbration f]
DIST
LINEar A B
CUTOn R
END
RATIos K L M..... ! 1 < K < L < M < ...
! K,L,M specify the different
! update frequencies
MAXStep t
CLEAr
END
The MTS keyword sigifies that subsequent dynamics runs are to utilize
multiple time steps. This keyword sets up the necessary data
structures to support multiple time steps. Dynamics is only initiated
with the DYNA MBOND command.
The number of stages for multiple time steps can either be specified
by the user (implicitly through the RATIos keyword) or computed
automatically. If it is not provided, it will always be computed
automatically. In addition, the optimal frequencies for updating the
various stages will be automatically computed if not provided.
The MASS keyword specifies that the "interaction mass" is used to segregate
the non-bonded force. The interaction mass depends on the
substructuring employed: for an atom in an atomistic region, it is the
mass of that atom; for an atom in a rigid body, it is the combined
mass of all atoms in that body; for an atom in a flexible body, it
is currently the mass of the atom.
When examining atom pairs in the construction of the non-bonded list,
the quantity SQRT((M1*M2)/M1+M2)) [where Mi is the interaction mass of
the i'th atom; which is the mass of the body if the atom belongs to a rigid
body, or is the mass of the atom itself if it belongs to a flexible body
or is a particle] is used to determine at what rate the non-bonded forces
for this pair will be computed. This ensures that interaction involving
atoms in bodies will, in general, be computed at a slower rate.
It is possible to calibrate the interpretation of this quantity (use
the CALIbrate keyword). The default value is 0.5 and typically values
are close to the fastest required update frequency (bearing in mind the
timestep which is specified later in the dynamics run).
If you use the DISTANCE keyword (recommended), the interaction mass
is multiplied by a DISTANCE factor C(r) (where r is the distance
between the two atoms). The form of the function C is specified by
the sub-options
CUTON r0 C(r) = 1 r < r0
LINEAR A B C(r) = A*(r - r0) + B r >= r0
Usually r0 ~ 4A. A must be positive and B should be 1.
The DISTance keyword implies the MASS keyword.
If you wish to expclicity define both the number of bins and the
different non-bond computation rates, use the RATIos keyword. Up to
four (integer) values may be provided. The fastest rate is always 1
and the values provided must be monotonically increasing and greater
than one. The number of values provided determines the number of bins
i.e. if one value is provided, the non-bonded interactions will be
separated into two stages; if 3 values are provided, 4 stages will be
used.
MAXStep is used to limit the biggest timesetep that is taken when
automatic ratios are computed. This is useful if you are using a
large value for the timestep in the DYNAMICS command and the automatic
computation may be producing huge net timesteps for the slowest update
rate.
CLEAR is used to disable body based multiple time steps. All ratios
(explicitly provided or automatically generated) are cleared.
NOTES:
1) A valid substructuring must have been defined before the MTS keyword
can be used.
2) Update frequencies in dynamics (non-bonded list generation, reporting
etc) are modified to be multiples of the lowest common multiple of the
MTS frequencies.
3) The MBOND Multiple Time Step only segregates non-bonded interactions.
Other force computations (e.g. dihedrals etc) are all computed at the
fastest rate. (We will be adding other keywords extend this).
4) There's not much benefit to using the MASS keyword by itself.
5) In theory, the interaction mass of atoms in a flexible body can be
larger than the mass of the individual atom. But in practice, the
benefits of doing this are not very great.
6) If you already have a large base timestep (15 or 20 fs), it is best to
just specify the exact ratios rather than computing them
automatically.
7) When specifying exact ratios, bear in mind that the LCM will determine
the non-bond update frequency. So
RATIo 2 3 4 6 LCM = 12
may be a better choice than
RATIo 2 3 4 5 LCM = 60
EXAMPLES
MBOND
SUB
(substructuring commands)
END
MTS ! Enable Multiple Time Step Dynamics
MASS CALIB .707 ! Use interactio mass with a factor of sqrt(2)
END ! The number of bins and the frequencies
END ! have not been specified and will be
! computed automatically
MBOND
SUB
(substructuring commands)
END
MTS ! Enable Multiple Time Step Dynamics
CLEAR ! Eliminate previously assigned parameters
DIST
LINEAR 2 1
CUTON 4
END
RATIo 2 4 6 ! 4 stages of force computation
END
The output of a multi-body dynamics simulation includes almost all of the same information generated during an atomistic simulation, together with additonal information pertinent to the body based approach. The standard dynamics files (trajectory, velocity, energy and restart files) are all supported during multi-body dynamics (there are no changes in format). Energy files have an additional line of BODY specific information per entry. Trajectories can be analyzed as usual and runs restarted if need be. NOTE that currently RESTART fits the current coordinates and velocities, the available modes/modal amplitudes and the coordinates at which the modes were generated. Care should be taken to check that this fitting process has not significantly modified the trjaectory. The energy reporting, controlled by NPRINT, has an additional line of information each step, which specifies the following terms: Deformation Energy Electrostatic energy as computed by MBOND's HBMA Momentum Aggregate Body Temperature Temperature in the atomistic regions The overall temperature and that of the body and atomistic regions can be separately controlled by the TFRQ option in the DYNAmics command. In addition, this will report the individual kinetic energy of each body. MBPRLEV controls the reporting of the following, MBPRLEV = 0 Modal and Kinetic Energy for each body Linear and Angular Momentum for each body MBPRLEV = 1 - The above plus Convergence rate during initial fitting MBPRLEV = 2 - The above plus All initial input data for each body and particle MBPRLEV = 9 - The above plus All sizes for the multibody dynamics arrays Complete description of the initial data Modal Amplitudes, body velocities and angular velocities
CHARMM Documentation / rvenable@deimos.cber.nih.gov