Monte Carlo
The Monte Carlo commands in CHARMM have been designed to allow
construction and use of an almost arbitrary move set with only a few atom
selections. This goal is accomplished by providing a pre-defined set of
move types which can be used in combination with each other for an arbitrary
CHARMM molecule. Speed and flexibility are gained by separating the
bookkeeping associated with a move from the actual application of that move
to the molecule.
* Menu:
* Syntax:: Syntax of MOVE and MC commands
* Description:: Description of MOVE and MC commands
* Examples:: Examples of MOVE and MC commands
* Data Structures:: Data structures shared by the MOVE and MC commands
* Shortcomings:: Known problems and limitations
* References:: Some references of use
Syntax for MOVE and MC commands
[Syntax MOVE < ADD | DELEte | EDIT | READ | WRITe > ]
MOVE ADD 1{ MVTP move-type } nsele{ SELE...END } -
[ WEIGht 1.0 ] [ DMAX 1.0 ] -
[ ARMP -1.0 ] [ ARMA 0.0 ] [ ARMB 0.0 ] -
[ DOMCf -1.0 ] [ ANISotropic 0 ] -
[ FEWEr 0 ] [ NLIMit 1 ] [LABEL move-label ]
where nsele, the number of SELE...END statements,
depends on move-type
move-type (nsele)::= < RTRN rig-unit ( 1 ) | ! Rigid translations
RROT rig-unit ( 1+) | ! Rigid rotations
CART ( 1 ) | ! Single atom Cartesians
TORS ( 2 ) | ! Simple torsion rotations
CROT ( 1+) > ! Concerted torsion rotations
rig-unit ::= < BYREsidue | BYALl >
MOVE DELEte < MVINdex move-index | LABEL move-label > -
MOVE EDIT < MVINdex move-index | LABEL move-label > -
[ WEIGht prev ] [ DMAX prev ] -
[ ARMP prev ] [ ARMA prev ] [ ARMB prev ] -
[ DOMCf prev ] [ ANISotropic prev ] [ NLIMit prev ]
prev ::= previous value
MOVE WRITE [UNIT -1]
MOVE READ [UNIT -1] [APPEnd 1]
[Syntax MC]
MC [ NSTEps 0 ] [ ISEEd prev ] [ TEMPerature 300.0 ] -
[ INBFrq 0 ] [ IMGFrq 0 ] [ IECHeck 0 ] -
[ IUNCrd -1 ] [ NSAVc 0 ] [ IMULti -1 ] [ IACCept 0 ]
[ IARMfrq 0 ] [ IDOMcfrq 0 ]
The flexibility and speed associated with the CHARMM Monte Carlo
commands derive from the separation of move bookkeeping (MOVE subcommands)
from the actual application of the moves to the molecule (MC).
MOVE
The MOVE subcommands are associated with construction of the move set.
The primary MOVE subcommand is MOVE ADD, which determines all the
locations that a move type can be applied to a subset of atoms (each location
is a move instance) and, for each instance, the relevant rotation axes and
centers, the moving atoms, and the relevant bonded terms. Thus, each call of
MOVE ADD results in a group of move instances with the same move type.
In addition, MOVE ADD associates with that group of move instances a set of
parameters (set by the optional arguments and described below).
Since using several different types of moves in conjunction typically
yields the most efficient and complete sampling, the MOVE ADD command is
normally called repeatedly.
The available pre-defined move types are rigid translations of residues
(RTRN), rigid rotations of residues (RROT), single atom displacements (CART),
simple torsion rotations (TORS), and concerted rotation of 7 (or, in the case
of an end, 6) torsions (CROT) (Go and Scheraga, 1970; Dodd et al., 1993).
Each of these can be applied to an arbitrary subset of atoms using
standard CHARMM SELE...END statements:
MVTP rig-unit nsele Description
---- -------- ----- -----------
RTRN BYALl 1 The entire atom selection is rigidly translated.
RTRN BYREsidue 1 The residue containing each selected atom is
rigidly translated. If more than one atom in
a residue is selected, each counts as a separate
move instance.
RROT BYALl 2 The allowed rotation centers are each of the atoms
in the first atom selection and the rigid body of
atoms being rotated is the entire second atom
selection. The first selection need not be a subset
of the second, so rotations around atoms outside
the rigid body can occur.
RROT BYREsidue 1 There is only a single atom selection, and each
selected atom is a center of rotation (around a
randomly selected axis) for its residue. If more
than one atom in a residue is selected, each counts
as a separate move instance.
CART 1 Each selected atom is an instance of a simple move
in Cartesian coordinates.
TORS 2 The two selections define the middle atoms (JK in IJKL)
of the rotatable torsions. If the FEWEr keyword is
set to 1, the directionality of the selection will be
ignored and each rotatable bond will be included only
once in the move set. Otherwise, each rotatable bond
will be included either once or twice depending on
whether the atom selections match the bond in only one
direction (JK) or both (JK and KJ). Pseudo-torsions
are not enumerated.
CROT 1+ The first atom selection defines the "backbone"
along which the 7 (or in the case of a chain end, 6)
dihedrals lie. Each additional pair of selections
defines non-rotatable bonds. The first bond in a set
of 6 or 7 is the driver torsion. Non-rotatable bonds
are not allowed at the third or fifth bonds following
the driver (counting only rotatable ones). Note that
there is no checking for whether bonds selected to be
rotatable are indeed so. NLIMit is the number of
torsions in addition to the driver torsion that are
restricted by the maximum rotation (DMAX); only 0 and 1
are implemented at present. In general, NLIMit 1 is
recommended since it speeds up the root finding process
and moves with large changes to the torsions are likely
to be rejected
anyway.
The values of the following optional parameters are used in all MC
calls, regardless of the values of its optional arguments (discussed below).
WEIGht The relative weight of that group of move instances in
the complete move set. The probability of picking a
group of move instances with weight w_i is w_i/(sum_j w_j)
where (sum_j w_j) is the total of all the WEIGht values.
DMAX The initial maximum displacement of each instance in a
group. Translations are in angstroms and rotations are in
degrees. In cases where anisotropic automatic optimization
is to be performed (described below), the initial assignment
is isotropic.
The following optional parameters are associated with automatic
optimization of the volumes in space from which individual move instances
are chosen. The two available methods are the Acceptance Ratio Method (ARM)
and Dynamically Optimized Monte Carlo (DOMC); both are described in detail
by Bouzida et al. (1992). The latter has the advantage that it allows
optimization of anisotropic volumes.
ARMP ARM target probability of move instance acceptance.
ARMA, ARMB Parameters to avoid taking the logarithm of zero in ARM:
DMAX(new) = DMAX(old)*ln(ARMA*ARMP+ARMB)/ln(ARMA*obsP+ARMB)
where obsP is the observed probability of accepting that
move instance.
DOMCF The F factor in DOMC:
DMAX(new) = DOMCF*SQRT[(d2ave*TEMP)/Eave]
where d2ave is the observed average square of the
displacement and Eave is the observed average change in
energy (both averages are done over all moves, not just those
accepted). DOMCF is used for the anisotropic version of
this equation as well. In the event that the square
root of a negative number must be taken, the routine
branches to ARM optimization, so ARMA, ARMB, and ARMP
should be set even if one plans on using DOMC.
ANISotropic DOMC anisotropic optimization of the volume from which the
moves are chosen. If ANISotropic is 0, it is off (isotropic)
and, if ANISotropic is non-zero, it is on. At present,
only 3D Cartesian moves (RTRN and CART) allow anisotropic
optimization.
LABEL An optional tag for the group of move instances.
Only the first four characters are retained. All sets of
move instances are also given an integer index which can
be used instead.
MOVE DELEte allows one to delete a group of move instances. The
group to be deleted is the first which matches the 4 character tag specified
by LABEL move-label or that specified by MVINdex move-index; if there is a
conflict, the LABEL is used.
MOVE EDIT allows one to change the values of the parameters associated
with a group of move instances. The matching rules are the same as those for
MOVE DELEte (as a result, the LABEL parameter itself cannot be changed with
MOVE EDIT). Any parameter not specified retains its current value. If a
positive value is specified for DMAX, all move instances will be reset to
that value (isotropic); if a negative value (default) is specified the
maximum displacements retain their current values. If DMAX is not specified
and the ANISotropic flag changes such that anisotropy is no longer allowed
(when it was allowed previously) the maximum displacements are assigned the
geometric mean of the eigenvalues of the matrix used to calculate the allowed
ellipsoid from the unit sphere.
MOVE WRITe writes out the current move set to a formatted file opened
with OPEN WRITE CARD.
MOVE READ reads in a move set. If APPEnd is 0, existing moves
are removed; otherwise they are preserved and the new moves are appended.
MOVE ADD calls can follow to further expand the move set.
MC
The MC command performs the loop over the appropriate number of dynamics
steps. Each step consists of: randomly picking a group of move instances
(weighted), randomly picking an instance from that group (unweighted),
calculating the energetic contribution of the moving atoms and their images,
applying the move, calculating the new energetic contribution, applying
the acceptance criterion, updating automatic move size optimization counters
if necessary, and finally performing any desired I/O (at present, only
trajectory writing is enabled).
NSTEps The number of loop iterations. Each pick of a single move
instance and subsequent application of the acceptance
criterion counts.
ISEEd The seed for the random number generator. If it is not
specified, it is unchanged, so that a script can be seeded
once initially and then loop over an MC command and yield
different behavior with each call.
TEMPerature The absolute temperature in degrees Kelvin.
IACCept The acceptance criterion to be used.
If IACCept is 0, Boltzmann (Metropolis) weighting is used.
If IACCept is 1, multicanonical (constant entropy) weighting
is used (in which case TEMPerature is ignored).
INBFrq The non-bond list update frequency.
If INBFrq is 0, the list is not updated.
Note: a call to ENERgy or UPDAte must be made before MC to
initialize parameters for non-bond list generation.
IMGFrq The image list update frequency (must be a multiple of
INBFrq). If IMGFrq is 0, the list is not updated.
IECHeck The total energy check frequency (must be a multiple of
INBFrq). The difference between the MC running total and
the current total is printed in the Delta-E column of the
table. If IECHeck is 0, the energy is not checked.
IUNCrd The I/O unit for trajectory writing.
NSAVc The frequency of writing out the trajectory.
If NSAVc is 0, no coordinates are written.
IMULti The I/O unit for reading in the multicanonical weights.
The file format (subject to change) is:
CHARMM title
Emin Emax Nbin
i E_i ln[n(E_i)]
.
.
.
Nbin E_Nbin ln[n(E_Nbin)]
Note that MC closes this file, so that it must be reopened
before each MC call with multicanonical weighting.
IARMfrq The frequency of updating the move size by ARM. Note that
this counter runs separately for each move instance.
If IARMfrq is 0, the move size is not updated.
IDOMcfrq The frequency of updating the move size by DOMC. Note that
this counter runs separately for each move instance.
If IDOMcfrq is 0, the move size is not updated.
If both IARMfrq and IDOMcfrq are non-zero, IARMfrq takes
priority.
EXAMPLES
No special actions must be taken during PSF generation for an MC run.
Essentially, input files set up for dynamics can be turned into ones for MC by
replacing the DYNAmics call with a series of MOVE ADD calls (or a MOVE READ
call) and an MC call.
Example: A peptide in water.
.
.
.
! Standard PSF generation above
! Create the MC move set
! Allow waters to move by rigid translations and rotations.
! Allow anisotropic optimization of the volume from which the
! translation vectors are chosen.
MOVE ADD MVTP RTRN BYREsidue WEIGht 2.0 DMAX 0.10 SELE (TYPE OH2) END -
ARMP 0.2 ARMA 0.8 ARMB 0.1 DOMCF 2.0 ANISo 1
MOVE ADD MVTP RROT BYREsidue WEIGht 2.0 DMAX 30.0 SELE (TYPE OH2) END -
ARMP 0.2 ARMA 0.8 ARMB 0.1 DOMCF 2.0 ANISo 0
! Allow all torsions to move by simple rotations
MOVE ADD MVTP TORS WEIGht 0.1 DMAX 30.0 FEWEr 1 -
SELE ALL END SELE ALL END
! Allow the backbone to move by concerted rotations with non-rotatable
! peptide bonds and N-CA proline bonds. If disulfides are present, the
! cysteine phi and psi should be restricted too.
MOVE ADD MVTP CROT WEIGht 0.5 DMAX 10.0 NLIMit 1 LABEL PEPTide -
SELE ((TYPE N).OR.(TYPE CA).OR.(TYPE C)) END -
SELE (TYPE C) END SELE (TYPE N) END -
SELE (RESNAME PRO .AND. TYPE CA) END -
SELE (RESNAME PRO .AND. TYPE N) END
! Seed the generator (for this example, this could be done below)
MC ISEEd 391004
OPEN UNIT 32 WRITE UNFO NAME EXAMPLE.TRJ
! Do 20000 steps at 300 K, writing energy 100 steps.
! Update the non-bonded list every 200 and
! the maximum displacements every 5 picks of that move instance
MC NSTEP 20000 TEMP 300.00 IUNC 32 NSAVc 100 INBFrq 200 IDOMcfrq 5
In this example, there are four groups of move instances (one for
each MOVE ADD call).
It should be mentioned that, since the MOVE READ command does not
do any checking as it reads in the necessary move set information,
it is possible to use moves in MC which MOVE ADD does not generate.
For example, it is straightforward to make rigid rotations around a
pseudo-dihedral simply by changing the pivot and moving atom lists
of a dihedral rotation. An understanding of the following section
(Data Structures) will aid in manual move creation.
Data Structures
MOVE ADD establishes each of the following pointers for all move types.
Each is a pointer to a dynamically allocated array that is n-instance elements
long, where n-instance is equal to the number of move instances in that group.
In all cases, if the array does not apply to a particular move, it is not
allocated.
MDXP This array contains the information about the limits of the
move. For isotropic or one-dimensional moves, it is simply
an n-instance-long array of reals containing the maximum
displacement. If the displacements are to be drawn from an
anisotropic volume, the array is a list of pointers, each of
which points to an array of 9 reals that make up the matrix
that transforms the unit sphere into the appropriate ellipsoid.
IBLSTP A list of n-instance pointers, each of which points to
a list of bonded terms changing under that move instance.
For each element in the four element array QBND (bonds=1,
angles=2, dihedrals=3, impropers=4) that is true, there is
an element listing the index of the final element containing
indices of that bonded term type followed by the list of
terms themselves. This list is then followed by a similar
one for the next bonded term type with QBND(i) set to true.
For example, if bonds 3, 8, and 10 and angles 16 and 17
were changing, the QBND array would contain (T T F F) and the
list would contain (4 3 8 10 7 16 17).
Urey-Bradley terms are handled with the lists generated for
angle terms, so do not get their own entries.
IPIVTP This array keeps track of any pivot or special atoms.
If there is only one pivot atom, then it is stored in the
array. If there are multiple (e.g., 2 for a TORS move
and 14 for a CROT move), the list stores a pointer to
a list containing the pivot atoms.
IMVNGP This array contains a compact list of the moving atoms.
Each element contains a pointer to a list of the following
form. The first element in the list is 1 more than the
number of rigid groups (NG). Elements 2 to NG contain the
index of the last array element with information about that
rigid group. The atoms in a rigid group are stored as
the first and last atoms in a contiguous range of atom indices.
Shortcomings
Attempts to move cross-linked residues will break MOVE ADD if
MVTP is CROT. If there is a large drift in the bond energies when
bonds lengths and angles are fixed, it is probably because a non-rotatable
bond (for example, the N-CA bond in proline) is being rotated by CROT.
Someday, a flag might be provided to choose between automatic elimination
of such moves and CROT-type relaxation of the cross-link (correct Jacobian
weighting is necessary to meet the detailed balance condition in the latter),
but such a flag is not on the immediate agenda of the MC developer.
The considered energy terms are bonds, angles, Urey-Bradley, dihedrals,
impropers, vdw, electrostatic, image vdw, image electrostatic, and user.
All non-bonded calculations can be either atom- or group-based. Terms not
listed above (e.g., constraints or explicit H-bond terms) are not included in
the present implementation.
Group-based calculations scale poorly with the size of the system
in the present implementation due to the structure of the CHARMM exclusion
list and the group non-bonded routines.
There is no heuristic update for the non-bonded list. However,
the structure of the MC symmetric non-bonded list is set up for such a
heuristic already, so one might be included in future versions.
No attempt has been made to see if the image structure in MC works with
the CRYStal command.
No warnings are printed for attempts to move a bonded (or patched)
residue by rigid translation and rotation.
References:
Berg, B. A. and Neuhaus, T. (1992) Multicanonical ensemble: A new approach
to simulate first-order phase transitions. Phys. Rev. Lett. 68, 9-12.
Bouzida, D., Kumar, S. and Swendsen, R. H. (1992) Efficient Monte Carlo
Methods for the Computer Simulation of Biological Molecules.
Phys. Rev. A 45, 8894-8901.
Dodd, L. R., Boone, T. D. and Theodorou, D. N. (1993) A Concerted
Rotation Algorithm for Atomistic Monte Carlo Simulation of Polymer
Melts and Glasses. Mol. Phys. 78, 961-996.
Go, N. and Scheraga, H. A. (1970) Ring Closure and Local Conformational
Deformations of Chain Molecules. Macromolecules 3, 178-187.
Leontidis, E., de Pablo, J. J., Laso, M. and Suter, U. W. (1994)
A critical evaluation of novel algorithms for the off-lattice Monte Carlo
simulation of condensed polymer phases. Adv. Polymer Sci. 116, 285-318.
Lee, J. (1993) New Monte Carlo algorithm: Entropic sampling.
Phys. Rev. Lett. 71, 211-214.
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and
Teller, E. (1953) Equation of State Calculations by Fast Computing
Machines. J. Chem. Phys. 21, 1087-1092.
Okamoto, Y. and Hansmann, U. H. E. (1995) Thermodynamics of helix-coil
transitions studied by multicanonical algorithms. J. Phys. Chem. 99,
11276-11287.
CHARMM Documentation / Rick_Venable@nih.gov