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.
For detail description of MBO(N)D method, refer to the paper of
"MBO(N)D: A multibody method for long-time molecular dynamics simulations"
in Journal of Computational Chemistry, Vol 21, 1 (2000).
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
* Harmonic:: Body-based modes from Harmonic Poteintial
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]
MAKEBODY SELE (...) 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.
MAKEBODY SELEct (....) END [SHOW]
Define bodies using the CHARMM regular SELECT facility.
i.e.
Makebody sele resname tip3 end show
- put each residue name "TIP3" to be in each separated body.
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] }
{[HARMonic] [FCON] [RCUT] [VDWR]} Optional-selection (see below)
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
***Optional-selection { [MULTiple] [BOND] [ANGLE] [ALL ] [VDWF] }
These options are documented on the next screen.
MODES { [KEEP ] } [NMODe integer]
{ [DELEte] }
{ [HARMomic] [FCON] [RCUT] [VDWR]} Optional-selection (see below)
***Optional-selection { [MULTiple] [BOND] [ANGLE] [ALL ] [VDWF]
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).
By applying "HARMonic" command after "MODES", the
simple harmonic potential, instead of the regular CHARMM potential,
will be used to generate modes for the body. (For more detail of the
harmonic potential mode, see *note Harmonic:(chmdoc/harmonic.doc).
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.
HARMonic invokes the harmonic potential mode generation for the
flexible-body.
----- Following commands works only if "HARM" command is defined.
FCON defines a single force constant of the harmonic potential
to be used for mode-generation. (No need to be defined when
MULT command is on.)
RCUT defines the cut-off distance for pairing atoms. Only distance
criteria is used to define a pair of atoms.
VDWR sets an initial cut-off distance to be a sum of van der Waals
radius of interacting atoms.
Final Cut-off distance = RCUT + VDWR-atom1 + VDWR-atom2
### Optional-selection Commands for Harmonic potential modes ###
BOND
ANGLE
ALL These options overide the distance criteria for defineing
pair of atoms. BOND command forces to make bonded atoms as
a interaction pair. ANGLE command forces to connect
atoms which have a common atom center. ALL command will involke
both BOND and ANGLE command. The remaining of atoms are
paired by the distance criteria.
MULTIple forces to use the different force constant based on the
nature of interactions. Each harmonic force constant is
obtained from the CHARMM parameters, such as bond, angle,
Lennard-Jones (VDW) parameters.
VDWF In MULT selection, the force constant (K) between paired atoms,
that are connected by the distance criteria, is defined by the
second derivative of the Lennard-Jones potential.
K = d^2 V /dr^2 at where r is the given distance.
IF VDWF is on,
K = d^2 V/ d r^2 at where r is the potential minima
#######################################################################
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. 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
SLFG BIN #1 RSCUT #2 BUFF #3
! #1 - bin number
! #2 - Cut off distance in Angstrom
! #3 - Buffer region in Angstrom
END
RATIos K L M..... ! 1 < K < L < M < ...
! K,L,M specify the different
! update frequencies
MAXStep t
ADD ## SELE (...) END ! ## - bin number
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.
In addition, you can use the capabilitiy to select distance-dependent
binning (Similar to the CHARMM MTS - refer to mts.doc) by the keyword, SLFG.
User can select only a distance-depenent binning, instead of mass-weighted
one described previously. By using this method, nonbond forces are
calculated by using switching functions in order to have smooth transitions
of forces. This method currently works with atom-based cut-off only.
Also by using the keyword, SLFG, you can invoke the capability to distribute
only interanal motions (bond, angle, dihedral, etc.) into the fastest bin.
If user uses "SLFG BIN 1 RSCUT 0.0 BUFF 0.0", user distribute all interanl
force calculations to the 1st bin. This can help to get the better stability
for hinge-off simulation.
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.
The ADD keyword is used to specify interaction, which involve selected atoms,
into the user specified bins in MTS. i.e. ADD 1 SELEct RESNAME TIP3 END -
This means that all interaction involving "tip3" residues are place into 1st
bin.
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
CUTON 4
END
RATIo 2 4 6 ! 4 stages of force computation
END
MBOND
MTS
RATIO 2 5
DIST
SLFG BIN 1 RSCUT 5.0 BUFF 1.0 ! Bin #1 has the interactions
! between 0 and 5A
SLFG BIN 2 RSCUT 8.0 BUFF 1.0 ! Bin #2 has the all interactions
! between 5 and 8.0A
! Rest of interactions are in Bin #3
END
END
END
MBOND
MTS
MASS CALI 0.5
RATIO 3 6
DIST ! 3 stages of force computation
Linear 1.0 ! but all interactions involving TIP3
CUTON 4 ! residues are place into 1st bin.
END
ADD 1 SELEct RESNAME TIP3 END
END
END
Langevin Dynamics is invoked for multibody dynamics by specifying the
LANGEVIN keyword in the DYNAMICS MBOND command. The usual CHARMM
control parameters, ISEED and TBATH, are read (the defaults are 314159
and FINALT (which itself defaults to 298.0)).
Friction coefficients for the forces are specified separately using
the MBOND LANGEVIN keyword. The available options are
MBOND
LANGEVIN
BODY {Name} COEF real [real, real,... real]
BODY ALL COEF real
ATOM ALL COEF real [real, real]
ATOM SELEct (standard selection) END COEF real [real, real]
LIST
CLEAR
END
END
Currently, a different coefficient can be specified for EVERY degree
of freedom that a body or particle has. So for each body, there are 3
translational, 3 rotational and 1 for each mode USEd for that body.
To simplify things, a single coefficient may be specified
(implementation detail - this is replicated 6+nModes times) for a
body. If you specify more than one, you MUST specify exactly
6+nModes.
Likewise, each particle has 3 degrees of freedom and either 1 or 3
coefficients may be specified. You can specify ALL particles or
select them using any CHARMM selection sequence. Likewise ALL bodies
can be done at once or you can deal with one at a time.
1) It is OK to specify values for everything using ALL and then
overwrite the specific values for a specific set of particles or a
specific body.
2) As bodies may have different numbers of modes loaded, currently you
can only specify a single coefficient when using ALL bodies.
Specification of the coefficients obviously requires that
substructuring has been defined AND that the appropriate mode
specifications made. If you subsequently change any of this, then
there is no way to track how the coefficients map to the new scheme
and they are ALL flagged as invalid.
Some examples
MBOND
SUBST
STREAM blah.str ! Defines substructuring
END
MODES
OPEN READ UNIT 16 modes.mod
LOAD bdy1 NMODe 4 UNIT 16
END
LANGEVIN
BODY ALL COEF 0.8
BODY bdy1 COEF 0.5 0.5 0.5 0.6 0.6 0.6 1.1 1.2 1.1 0.9
ATOM ALL COEF 1.0
END
END
DYNAMICS MBOND -
LANGEVI ISEED 3250 TBATH 200 -
LOBATTO NCYC 3 -
etc
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 ^_
Body-Based Mode Generation by Harmonic Potential
=============================================
1) Introduction
Generating reasonable body-based modes is a crucial part of MBO(N)D. Those
modes are used for modeling flexible bodies. In recent years, there have
been several research papers in which people described a simple approach
to generate reasonable thermal fluctuations of proteins. In the approach,
a single parameter harmonic potential was used instead of the more complicated
force field. Based on their observations, low-frequency mode-vectors
generated by this simple potential was in good agreement with those calculated
from the full potential.
From the MBO(N)D point of view, this approach for generating body-based
modes is attractive. The instantaneous structure can be treated as an
equilibrium structure in this approach so that no artificial energy minimization
is necessary. Other approach, using the full CHARMM poteintail, requires
small steps of energy minimization before generating modes because
inappropriate modes produced from slight distorted peptide-plains. In the same
vein, no negative frequency eigenvectors are generated. In the regular
MBO(N)D mode generation approach, several of the body-based modes are usually
of negative-frequency. However, harmonic potential approach will eliminate
those problems altogether.
References:
1. M.M. Tirion, Phys. Rev. Letts., 77, 1905 (1996)
2. I. Bahar, A. R. Atilgan, and B. Erman, Folding & Design, 2, 173 (1997)
2) Command Description
- Generation of body-based modes by using the harmonic potential
instead of the regular CHARMM potential.
By applying "HARMonic" command after "MODES", the simple harmonic
potential, instead of the regular CHARMM potential, will be used to
generate modes for the body. This "HARM" command works only when
the body-based modes are generating. If you want to read the existed
modes, this command will not do any effects.
The syntax for the MODE sub-command for using the harmonic potential, is:
MBOND [NBODY integer]
MODEs { [HARMomic] [FCON] [RCUT] [VDWR]} Optional-selection (see below)
- Regular mode generation command, such as USE, SORT, WRITE, etc.
(Check interface.note)
END
END
***Optional-selection { [MULTiple] [BOND] [ANGLE] [ALL ] [VDWF]}
----------------------------------------------------------
MODES { [HARMomic] [FCON] [RCUT] [VDWR] } Optional-selection
options:
HARMonic invokes the harmonic potential mode generation for the
flexible-body.
FCON defines a single force constant of the harmonic potential
to be used for mode-generation. (No need to be defined when
MULT command is on.)
RCUT defines the cut-off distance for pairing atoms. Only distance
criteria is used to define a pair of atoms.
VDWR sets an initial cut-off distance to be a sum of van der Waals
radius of interacting atoms.
Final Cut-off distance = RCUT + VDWR-atom1 + VDWR-atom2
### Optional-selection Commands for Harmonic potential modes ###
BOND
ANGLE
ALL These options overide the distance criteria for defineing
pair of atoms. BOND command forces to make bonded atoms as
a interaction pair. ANGLE command forces to connect
atoms which have a common atom center. ALL command will involke
both BOND and ANGLE command. The remaining of atoms are
paired by the distance criteria.
MULTIple forces to use the different force constant based on the
nature of interactions. Each harmonic force constant is
obtained from the CHARMM parameters, such as bond, angle,
Lennard-Jones (VDW) parameters.
VDWF In MULT selection, the force constant (K) between paired atoms,
that are connected by the distance criteria, is defined by the
second derivative of the Lennard-Jones potential.
K = d^2 V /dr^2 at where r is the given distance. (1)
IF VDWF is on,
K = d^2 V/ d r^2 at where r is the potential minima (2)
3) Examples
----- A simple single parameter approach --------
(1)
MODEs HARM VDWR RCUT 2.0 FCON 100.0
Open write ...
generate ....
Diag ...
write ...
END
- A single parameter Harmonic potential approach.
- Connecting pairs are defined by the distance cutoff.
When The distance between two atoms are less than
2.0 + vdw radius of atom 1 + vdw radius of atom2,
those atoms are connected by a single spring.
(2)
MODEs HARM BOND RCUT 5.0 FCON 100.0
.....
END
- A single parameter harmonic potential approach
- First, atoms which are covalent bonded to each other
are connected by a spring with force constant = 100kcal/mole*A^2
Then, rest of atoms are examined by the distance criteria.
Here, two atoms, not covalent bonded, are connected by a spring
if the distance between those atoms are less than 5A.
(3)
MODEs HARM ALL RCUT 2.0 FCON 50.0
....
END
- A single parameter harmonic potential approach
- First, atoms which are covalent bonded to each other or are connected
by a common atoms (ANGLE) are connect by a spring with
force constant = 50Kcal/mole*A^2. Then, rest of atoms are checked
by the distance. If the distance between two atoms are less than
2A, those atoms, which are not already paired by the first check
(BOND and ANGLE) are connected by a spring.
-------- Multipe harmonic potential approach -------
(4)
MODEs HARM MULTI VDWF RCUT 5.0
.....
END
- Multiple parameter harmonic potential approach
- Atom pairs are solely created by checking the distance between
atoms. Then, each pair's force constant is calculated by the
second derivative of Lennard-Jones potential. "VDWF" command
involk to use Eq. (2), defined above, for the calculations of
the force constants.
(5)
MODEs HARM MULTI ALL VDWF RCUT 5.0
....
END
- Multiple parameter harmonic potential approach
- First, atoms which are covalent bonded to each other or
are connected by a common atoms (ANGLE) are connect by a spring
with force constants which are calculated from the CHARMM
parameters (BOND and angle parameters) Then, the rest of atoms
are checked by the distance criteria (<5.0A). IF the distance
is less than 5A, the atoms are connected by a spring with a
force constant calcualted by Eq. (2).
(6)
MODEs HARAM
.....
END
- This simple command will be treated as the same command as
(5).
If you have any question, please contact Masa (watanabe@moldyn.com)
NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute