Constant Pressure/Temperature (CPT) Dynamics
Two types of constant pressure/temperature dynamics are available
in CHARMM. The weak coupling method for temperature and pressure
control described in the paper by Berendsen et al. (JCP 81(8) p3684
1984) was the first constant pressure and temperature algortihm
implemented in CHARMM. Extended system constant pressure and
temperature algorithms have now been implemented based on the work
of Andersen (JCP 72(4) p2384 1980), Nose & Klein (Mol Physics 50(5)
p1055 1983), Hoover (Phys. Review A 31(3) p1695 1985). Additionally,
a variant on the extended system method which treats the control
variables by means of a Langevin equation is available (Feller, Zhang,
Pastor & Brooks, JCP, 103, 4613 (1995)).
Shape matrix propagation and coordinates scaling for triclinic
unit cell is done according to D. Brown and J.H.R. Clarke in
Computer Physics Comm. 62 (1991) 360-369.
A constant surface tension algorithm is included which is useful
for studying interfacial systems where one wishes to allow the area
to change dynamically during the simulation. The dynamical equations
and statistical ensemble are discussed in (Zhang, Feller, Brooks &
Pastor, JCP, 103, 10252 (1995)).
* Syntax / Syntax of the CPT dynamics command
* Description / Description of the keywords and options
* Notes / Other points to be noted
* Examples / Isotropic and interfacial systems; constant tensor
* Pressure / The pressure command
Top
Description of CPT Dynamics Commands
Only a few changes are needed to a standard CHARMM dynamics input
file to run a CPT MD simulation. There are a few things to note :
a. The CPT algorithm is invoked by the CPT keyword.
b. It's not possible to use LANGEVIN dynamics with the constant
pressure or temperature algorithm.
c. All the non-Langevin dynamics keywords have the same meaning as
the in a standard dynamics input file. This includes the
keywords STRT and REST.
e. The CPT specific keywords (apart from CPT itself) are :
1. PCONstant - do a constant pressure calculation. Extended system
algorithm is the default, weak-coupling is available
with BEREndsen keyword.
2. TCONstant - do a constant temperature calculation with the weak-
coupling algorithm. HOOVer constant temperature is
only available with PCONstant simulations.
3a. To be used with Berendsen algorithm:
COMPressibility <real> - the isothermal compressibility
(atmospheres**-1).
PCOUple <real> - the pressure coupling constant
(picoseconds).
3b. To be used with extended system algorithm:
PMASs <real> - the mass of the pressure piston (amu)
PGAMma <real> - Langevin piston collision frequency (1/ps)
TBATh <real> - Langevin piston bath temperature
TENSion <real> - reference surface tension (dyne/cm)
IUPTEN <int> - unit number, P tensor at every step
P21XCEN - flag; enables dynamic XCEN for P21 TETRagonal
To be used with either constant pressure algorithm:
PREFerence <real> - the reference pressure (atmospheres).
(for isotropic pressure)
PREFInitial <real> - initial reference pressure tensor (atmospheres).
(for isotropic pressure)
PREFFnitial <real> - final reference pressure tensor (atmospheres).
(for isotropic pressure)
PRXX,PRYY,PRZZ <real> - the reference pressure tensor (atmospheres).
PRXY,PRXZ,PRYZ (for anisotropic pressure)
PIXX,PIYY,PIZZ <real> - initial reference pressure tensor (atmospheres).
PIXY,PIXZ,PIYZ (for anisotropic pressure)
PFXX,PFYY,PFZZ <real> - final reference pressure tensor (atmospheres).
PFXY,PFXZ,PFYZ (for anisotropic pressure)
PREFI,PREFF,PIXX...,PFXX... - are used for linear pressure ramping
4a. To be used with Berendsen algorithm
TCOUple <real> - the temperature coupling constant
(picoseconds).
TREFerence <real> - the berendsen reference temperature (K).
4b. To be used with extended system (HOOVer) algorithm
TMASs <real> - the mass of the thermal piston (kcal*mol^-1*ps^2).
REFT <real> - the hoover reference temperature (K).
Note : for full descriptions of these parameters and the suggested
values to use see the reference given above.
f. The CPT module is only available for use with the leap-frog integrator
Top
Other Points
Suggested values for solvated systems:
COMPressibility (beta) = 4.63e-5 /atm (for proteins)
PCOUple = 5.0 ps (or more)
PREF = 1.0 atm (default)
PMASs = 500 amus (default is infinity)
TCOUPle = 5.0 ps (or more)
TREF = 298.0 K (default)
REFT = 298.0 K (default)
TBATh = 298.0 K (default)
TMASs = 1000.0 kcal ps^2 (default is infinity)
TENSion = 0.0 (default, results in regular constant pressure)
Other Points
------------
1. Although the heating and equilibration commands are the same as
for standard dynamics it is possible to use the CPT algorithm
to perform both without velocity modification (c.f. Langevin
dynamics).
2. The algorithm requires the use of the CHARMM CRYSTAL facility
for constant pressure dynamics. If only Berendsen constant temperature
is requested, then the crystal code need not be used.
3. For the Berendsen algorithm, when a reference pressure term
( PRXX,PRYY,PRZZ,PRXY,PRXZ,PRYZ is
set to a very large negative number (less than -9999.0), then this
component of the pressure will not be considered. For orthorhombic
simulations, the particular box length will not change (for example,
if the command says: PRXX 1.0 PRYY -10000. PRZZ -10000. then
only the box length in the x direction will change during dynamics).
If a cubic simulation is performed, then the pressure terms corresponding
to the large negative values are not considered in the calculation of
the instentaneous pressure.
4. For the extended system pressure algorithm, setting any component of the
piston mass array to zero will cause the corresponding simulation cell
length to remain constant. For example when using the orthormobic
cell, setting pmxx=pmyy=0 results in only the zdirection changing during
the dynamics (and it changes according to the z component of the pressure
tensor). This is the standard method for interfacial NPAT systems.
5. A discussion of Hoover temperature control can be found in the documentation
file nose.doc. The temperature control implemented in the velocity verlet
integrator is very similar to the one used in the leapfrog integrator.
NOTE: Hoover temp control only works in conjunction with constant pressure
6. The Berendsen pressure/temperature control scheme may not be appropriate
for inhomegenous systems (protein in water, aqueous membrane, interfacial
systems). This is especially true if SHAKE constraints are used on one
component. A full discussion is given in the paper by Feller, Zhang, Pastor
and Brooks (JCP, 9/15/95).
7. The extended system pressure algorithm can be run with temperature control
(resulting in isothermal-isobaric ensemble) or without (resulting in
isoenthalpic ensemble). The Berendsen pressure method must be run with the
constant temperature control (the ensemble for these methods is unknown).
8. If PINTernal is used (default), then the pressure is determined by the
internal virial and the atoms in the box are instantaneously scaled in a
homogeneous response to the altered box dimensions. If PEXTernal is used,
then the external virial (related to the force required to maintain the
symmetry constraint) is used and the atom positions are not instantaneously
scaled by box size changes. The PEXTernal option is normally used for
minimization, but may be used with molecular dynamics. It is not recommended
for large systems.
9. The pressure tensor (extended system) on every integration step is saved
to the formatted file indicated by unit number specified by IUPTEN; this
allows calculation of viscosity using the Green-Kubo relationship. The
viscosity is computed from the integral of the autocorrelation function
of the off-diagonal elements of the P tensor, scaled by V/kT; see
J. Phys. Chem. 100:17011-17020 (1996). The column order of the data is:
Time PIXX PIXY PIXZ PIYX PIYY PIYZ PIZX PIZY PIZZ
10. For simulations using P21 symmetry with a tetragonal prism lattice, the
P21XCEN keyword enables dynamic updating of the XCEN value used for image
centering to be consistent with 0.25*A; for details on P21 symmetry, see
the example in (chmdoc/crystl.doc)
Top
Examples of Constant Pressure Usage
1. Basic constant pressure, appropriate for a cubic simulation cell, box
length is coupled to the trace of the pressure tensor, using langevin on
pressure piston degree of freedom. This also works for a non cubic
cell, but in that case each length moves independently to maintain a
constant pressure tensor. Constant volume is the limit pmass -> infinity
(implementation in CHARMM: Set pmass = 0 for constant V).
dynamics cpt leap restart time 0.001 nstep 10000 iseed 314159 -
pconstant pmass 400.0 pref 1.0 pgamma 20.0 -
tbath 300.0
2. Constant normal pressure, constant area. Appropriate for orthorombic
cell where only the z direction is allowed to change. The box
length in z direction is coupled to the z component of the pressure tensor.
dynamics cpt leap restart time 0.001 nstep 10000 iseed 314159 -
pconstant pmzz 225.0 pmxx 0.0 pmyy 0.0 pref 1.0
3. Constant pressure (stress) tensor. Each box length moves independently
to maintain the desired pressure tensor.
dynamics cpt leap restart time 0.001 nstep 10000 iseed 314159 -
pconstant pmass 225.0 przz 1.0 prxx 2.0 pryy 3.0
4. Constant pressure, constant surface tension. Z direction moves independently
of x and y and is coupled to the bulk pressure (z component of pressure
tensor). X and y box lengths move to maintain constant surface tension.
Note: this is only appropriate for interfacial systems where the interface
is perpendicular to the z axis.
dynamics cpt leap restart time 0.001 nstep 10000 iseed 314159 -
pconstant pmass 225.0 pref 1.0 surface tension 50.0
5. Constant pressure and temperature (NPT)
dynamics cpt leap restart time 0.001 nstep 10000 iseed 314159 -
pconstant pmass 400.0 pref 1.0 pgamma 20.0 -
tbath 300.0 tcons hoover reft 300. tmass 1000.
Examples of Constant Temperature Usage
1. Basic constant temperature using the Hoover method
dynamics cpt leap restart time 0.001 nstep 10000 iseed 314159 -
HOOVer TMASs 1000.0 REFT 298.0
2. Constant T with calculation of pressure data; this will also print
the surface tension in the output log, useful for studying
interfacial systems. The optional IUPTEN keyword will store the
pressure tensor data for every timestep.
open unit 29 card write name dyn.ptn
dynamics cpt leap restart time 0.001 nstep 10000 iseed 314159 -
pcons pmass 0.0 pint pref 1. iupten 29 -
hoover tmass 1000.0 reft 293.0
Top
The PRESsure command
Process the pressure commands for the system. There are three
modes :
Mode 1 : Initialise all pressure arrays.
Syntax:
PRESsure INITialise
Mode 2 : Calculate and print the instantaneous pressures for
a system.
Syntax:
PRESsure INSTantaneous TEMPerature <Real> VOLUme <Real> -
NDEGf <Integer> NOPRint
The external isotropic pressure and tensor are calculated
if a volume is present. The isotropic internal pressure is
calculated if a volume is present and a temperature has been
given. If no degrees of freedom are specified then a value
of 3*NATOM is taken be default. The virials are always printed.
NOPRint will suppress all printing.
Note: a previous call to energy is required so that the
virials (and volume if not specified) are available
in ENERGY.FCM. The command also accumulates the
average and the square of the pressure variables.
Mode 3 : Print the averages and fluctuations.
PRESsure STATistics
