CHARMM c42b2 repdstr.doc



File: Repdstr -=- Node: Top
Up: (commands.doc) -=- Next: Syntax


                       The Parallel Distributed Replica

                      By Paul Maragakis and Milan Hodoscek, 2005

     Parallel ditributed replica allows independent replicated systems
over specified number of processors. It mainly works with CMPI
pref.dat keyword (YMMV). REPDSTR is still not the default pref.dat
keyword so the recommended way to compile CHARMM is the following:

install.com gnu M mpif90 +REPDSTR +ASYNC_PME +GENCOMM [+MSCALE]

for install.com em64t add +CMPI to the above list.

Also MSCALE is not really needed for pure REPDSTR runs, but it is
needed for triple parallel CHARMM.

For one of the examples of REPDSTR usage see this reference:

Jiang, W; Hodoscek, M; Roux, B; "Computation of Absolute Hydration and
Binding Free Energy with Free Energy Perturbation Distributed
Replica-Exchange Molecular Dynamics", J. Chem. Theo. and Comp., 2009,
Vol. 5, 2583-2588.

For the reservoir replica exchange code, please cite the following
references:

Boltzmann reservoir REX--

Okur A., Roe D., Cui G., Hornak V., Simmerling C. J. Chem Theo. Comput. 
   3, 557-568 (2007).

Non-boltzmann reservoir REX--

Roitberg A., Okur. A., Simmerling C. J. Phys. Chem. B. 111, 2415-2418
   (2007).

CHARMM implementation--

Okur A., Miller B. T, Joo K., Lee J., Brooks B. R. J. Chem Theo. Comput.
   9, 1115-1124 (2013).

* Menu:

* Syntax::              Syntax of the REPD command
* I/O::                 Input and output functionality
* FAST::                Usage of FAST replica exchange
* Examples::            Examples to show the possibilities
* Output::              Explanation of the replica exchange printout



File: Repdstr -=- Node: Syntax
Up: Top -=- Previous: Top -=- Next: I/O



REPDstr [ FAST ] NREP <int>  [ NATRep <int> ]


     Replicates the system <int> times. Optional NATRep limits the
number of atoms to be included in the path calculations (RPATH
commands). It also reduces the size of arrays that need to be
transfered between replicas in the RPATH calculations.


REPDstr NREP <int>  [REPEat <int> [ LOGLevel <int> ] { EXCHange replica-exchange-spec }
                                                     { PHREx ph-replica-exchange-spec }

replica-exchange-spec::= FREQuency <int> temperature-spec 
                         [UNIT <int>] [SUMP]
                         [NREP <int>]
                         [SGLD] sgld-replica-exchange-spec
                         [TIGEr] tiger-spec
                         [RSRV] reservoir-spec
                         [TWOD] twod-spec
                         [EWRU <int>]

temperature-spec::= [ [ TEMP <real> ] TEMP <real> ... ]
                    [STEMperature <real> DTEMperature <real> MTEMperature <real>]
                    
tiger-spec::= [ITER <int> ] [NEQU <int>] [NMIN <int>] [TOLG <real>]

ph-replica-exchange-spec ::= FREQuency <int> ph-spec [UNIT <int>] [ RSVR reservoir-spec ]

reservoir-spec ::=  RESHigh  { BOLTzmann   }  RHTEmp <real> RHUNit <int> RHSZ <int> [ RHEN <int> ]
                    RESLow   { NOBOltzmann }  RLTEmp <real> RLUNit <int> RLSZ <int> [ RLEN <int> ]
                    [ ECOR ] [ FHEN <int> ] [ FLEN <int> ]

ph-spec ::= [ [PHVAl <real] PHVAl <real> ... ]

sgld-replica-exchange-spec::=  [ [ SGTT <real> ] SGTT <real> ... ]
                               [SGTE <real> DSGT <real> MSGT <real>]
                               [SGFT <real> ] DSGF <real> ... ]


twod-spec::= DIM1 criteria-spec DIM2 criteria-spec D1FRequency <int> D2FRequency <int>
             temp-spec ph-spec 

criteria-spec::= { TEMPerature }
                 { HAM         }
                 { PH          }

     This is for the replica exchange method (see details in
c34test/rexc.inp, c36test/rexcpt.inp, c36test/rexsgld.inp) Currently
it works so that when exchange occurs all the coordinates and
velocities are exchanged, thus the lowest temperature is always on the
first replica. This also implies that the number of atoms in the
replicas have to be the same. The other method which exchange only the
temperature will be implemented later. With just temperature exchange
the replicas do not need to be the same anymore.

     Current implementation of replica exchange methods has its own
temperature control independent of the CHARMM's one. So in the case of
exchanging the coordinates and velocities also the appropriate
temperature scaling is perfomed. Perhaps it is best to turn CHARMM's
own temperature controls off, but one can also combine the two. To get
both temperature control mechanisms at the same time one need to define
different temperature for each replica. This can be accomplished by
the following commands in the CHARMM input script:

set st 300
set dt 10

repd nrep @nreps EXCHange FREQuency 50 STEMp @st DTEMp @dt sump unit 17

mult dt by ?myrep
incr st by @dt

dyna cpt start nstep 1000 timestep 0.001 -
....
    hoover reft @st tmass 2000.0 tbath @st -
....

     As of CHARMM version c37a2, replica exchange in pH space is also
supported. This uses the formalism described in the reference given in
consph.doc. The CONSPH key word must be in pref.dat along with REPDSTR
to activate this functionality.


   Details about each keyword:

   UNIT <int>  - Optional keyword for exchange output file. Default for
                 int is OUTU. Must be opened after repd: each replica
                 writes to its own file. If no open statement all
                 replicas write to the same file with the default
                 fortran file name for this unit. Open before the repd
                 command is not very useful. Can be also the same unit
                 as on the OUTU command so exchange info is written to
                 the same output files.

   SUMPrint    - Summary printout. On replica zero the summary from
                 all other replicas is printed to UNIT, and on the
                 rest of the repicas just their own data. This is
                 flaged since it requires extra communication just for
                 printouts.

   FREQuency   - when to exchange

   REPEat      - Number of times to repeat an exchange attempt every FREQ
                 steps.

   STEM        - Starting-temperature.
   DTEM        - temperature-increase
   MTEM        - Top temperature. When MTEM>0, DTEM is ignored and temperatures 
                 of replicas are expoenentially spaced.
   TEMPerature - Temperature of each replica if STEM is not set. It must be 
                 repeated NREP times.
   PHVAl       - pH value of each replica when replica exchange in pH space is
                 used.

   SGLD        - Flag to do RXSGLD with the self-guiding
                 temperature. It maybe used with the standard replica
                 exchange or one can specify all the temperature the
                 same, most convenient with the STEM <temp> DTEM 0.0.
   SGTE    0   - The self-guiding temperature for the first replica.
   DSGT    0   - Increment for the self-guiding temperature
   MSGT    0   - The top self-guiding temperature. If MSGT>0, DSGT is ignored and
                 the guiding temperatures of replicas are expoenentially spaced.
   SGTT        - Self-guiding temperature of each replica if SGTE is not set. 
                 It must be repeated NREP times.
   SGFT    0   - The guiding factor of the first replica. When the self-guiding
                 temperatures are set with SGTE..., SGFT will be adjusted 
                 automatically during simulation.
   DSGF    0   - Guiding factor increment.
   EWRU <int>  - Energy write out unit - this parameter is only active with
                 Hamiltonian replica exchange. It writes a log of the energies
                 for the replica and its partner at each exchange attempt. These
                 energies can be read into the FREN command to calculate overlap
                 between replicas (see fren.doc for details). If this keyword
                 is not specified, this information is written out to the unit
                 given by the UNIT keyword, if greater than 0.


   TIGEr       - Flag to start TIGER replica exchange.

   ITER        - number of iteration steps for minimization and
                 equlibration procedures before the exchange.
                 Default: 1

   NEQU        - number of steps in the equlibration process.
                 Default: 1000

   NMIN        - number of steps in the minimization process.
                 Default: 100

   TOLG        - gradient tolerance in the minimization step.
                 Default:0.0

   PHMD        - Flag to allow exchange of theta variables from CPHMD
                 along with spatial coordinates. Thus, replicas can be run 
                 at different pH (Hamiltonian replica exchange) or temperture. 

     It is also possible to couple the top or bottom replicas (or
both) to a reservoir of structures. To do so, the RSVR keyword is
used. When RSVR is used, at least one of the following keywords
must also be used with the corresponding unit numbers to tell CHARMM
which replica(s) should be coupled to the reservoir.

   RESH         - couple the top replica to a reservoir. RHUN must
                  be specified.

   RESL         - couple the bottom replica to a reservoir. RLUN
                  must be specified.

   RHSZ         - The number of elements in the top reservoir.

   RLSZ         - The number of elements in the bottom reservoir.

   RHEN         - A unit number pointing to a data file listing the
                  potential energies of each reservoir structure in
                  the top reservoir. The data file should be formatted
                  in order with one energy per line. This is ONLY 
                  required for Hamiltonian reservoir replica exchange.

   RLEN         - Identical to RHEN, but for the bottom reservoir.

RHUN and RLUN must be units that point to open files in a simplified
trajectory format. This format is a standard CHARMM binary trajectory 
file, but has the header and crystal information stripped out. A utility,
simpletraj.py, is provided in the support/programs directory to convert a standard
CHARMM trajectory into a stripped down version suitable for use. These
files musrt be opened in DIREct mode with a record size specified (which
is four times the number of atoms except for pH reservoirs).

Two exchange schemes have been implemented to govern coupling of the
reservoir with its neighboring replica.

   BOLTzmann    - The standard Boltzmann temperature replica exchange
                  criterion is used. Use of this keyword implies that
                  the reservoir is a sample from a Boltzmann distribution.
                  If this option is used, RHTEmp and/or RLTEmp must
                  be used to specify the temperatures of the high and
                  low reservoirs, respectively.

   NOBOltzmann  - This option allows for a non-Boltzmann weighted
                  reservoir, using a slightly different exchange
                  criterion.

The RHUNit and RLUNit key-words tell CHARMM how many structures are in the 
high and low reservoirs. If FHEN or FLEN are used, then the energies of
the high and low structures are read from the given unit, one floating
point number per line. The number of lines in the file must match the number
of structures in the reservoir, and the order of the lines must correspond
to the order of the structures in the reservoirs. If these options are
omitted, the energy of each structure is calculated by CHARMM. If the ECOR
key word is specified, 0.5*kT of energy is added to the structure energy
for each degree of freedom in the system, which provides a quick and dirty
way of adjust the structure energy to the desired temperature.

    As of CHARMM c40a1, the reservoir replica exchange scheme has been
extended to Hamiltonian and discrete-state (CONSPH key-word) pH replica
exchange. It does NOT work with the continuous-state pH replica exchange
code (PHMD key word). Hamiltonian and pH reservoir replica exchange are
only implemented for BOLTzmann reservoir replica exchange; use of the
NOBOltzmann exchange criteria will throw an error.

REPDstr RESEt [ SYNC ] [ PONE ]

     Resets the run to a normal parallel run. SYNC does the global
sync before that. PONE is making for everybody NUMNOD=1. As of March
2010 RESET is still not fully supported.

REPDstr IORES

     Sometimes within the REPDstr run one wants to access the files
created by other replicas. After this command is executed the names in
the open command do not get _myrep appended!

REPDstr IOSET

     Sets the appending of the replica number back to original nameing
scheme in REPEDstr.


REPDstr NREP <int> EXLM [EXPT NRPT <int>] FREQuency <int>

     This is for Hamiltonian exchange method. Currently it works so that
when exchange occurs all the coordinates are exchanged and new nonbond list 
are generated. To guarantee stable md run after exchange, velocities also
are exchanged once an exchange attempt is accepted. The present Hamiltonian-
exchange scheme works for all integrators, including VV2 integrator for
Drude oscillator model.

   EXLM - Keyword invoking Hamiltonian exchange. Currently it can be used
        in Free Energy perturbation and umbrella sampling.
   EXPT NRPT <int> - Optional keyword introducing parallel tempering into
        Hamiltonian exchange. With this keyword, the replica-exchange
        consists of two alternative stages: parallel tempering and Halmiltonian
        exchange. In parallel tempering stage, the number of replicas
        participating exchange is NREP, while in Hamiltonian exchange stage,
        the number of replica is NREP/NRPT. Currently it can be used to
        accelerate the relaxation of internal degrees of freedom, such as
        sidechain dynamics and backbone dynamics

REPD NREP <int> EXLM EX2D NRPX <int> FREQ <int> 

     This is a new 2 Dimensional Hamiltonian Replica exchange scheme. 
Hamiltonian-Exchange is extended to PBC systems for either NVT and NPT
simulation. This new feature is especially useful to enhance samplings of
umbrella sampling that involve multiple reaction coordinates.

     EX2D NRPX <int> - Optional keyword introducing 2D replica exchange.
With this keyword, the number of replicas along X (one reaction coordinate)
is NREPX, then the number of replicas along the other reaction coordinate
is NREP/NREPX.

REPD TWOD DIM1 <int> DIM2 <int> D1FR <int> D2FR <int> -
     D1CR <string> D2CR <string> -
     temp-spec
     ph-spec

    This is a alternate 2D replica exchange implementation that allows
each dimension to be temperature, hamiltonian, or pH replica exchange.
This is set by the two dimension criteria D1CR and D2CR, each of which
may be TEMP, HAM, or PH. Note that only one TEMP or PH dimension may be
used. However, two hamiltonian dimension can be used. The number of
replicas in each dimension is set by DIM1 and DIM2. The exchange
frequency for each dimension is set by D1FR and D2FR. D1FR should be
set smaller than D2FR. If D2FR is a multiple of D1FR, then no dimension
1 exchange will occur on the same step as an exchange in dimension 2.

    One VERY important note -- in the current implementation, although
temperature and pH must be set in the REPD command line, they are
not set automatically in the DYNAmics command! Rather they must be
set using the substitution variables below. Please see an example
in test/c40test/rex-2d.inp of how to do this. This will be corrected
in a future release.

    After this command has been executed, the following substitution
valriables become available:

    ?nrep    - number of replicas overall
    ?myrep   - global index of the current replica
    ?nrepd1  - number of replicas in dimension 1
    ?nrepd2  - number of replicas in dimension 2
    ?myrepd1 - current replica's position in dimension 1
    ?myrepd2 - current replica's position in dimension 2



File: Repdstr -=- Node: I/O
Up: Top -=- Previous: Syntax -=- Next: Fast



    Once REPDstr command is activated the I/O capabilities of CHARMM
are expanded. In standard parallel mode CHARMM deals with I/O only on
the first process. The rest of processes get their data through
network or memory communication. So all I/O statements that are in the
script before REPDstr command are valid only on first process. In
distributed replica mode each replica needs its own and independent
I/O which is enabled after the REPDstr keyword in the input script.
Two substitution parameters are defined after the REPD command is
specified in the input script: ?NREP (number of replicas) and ?MYREP
(current executing replica).

As of May 2009 the following is working:

1. OPEN

   The command open read|write unit 1 card name somefile will open
   somefile_0 for replica 0, somefile_1 for replica 1, etc

2. READ/WRITE

   writes to individual files one for each replica. It works for all
   I/O operations.

3. STRE stream

   This will open stream_0 for replica 0, stream_1 for replica 1, etc
   It allows CHARMM to run different input files for each replica (or
   group of processors)

4. OUTU unit

   Will stream output to individual files as specified in the open
   command for particular unit. This command should precede STRE
   command if one wants both input and output files for each group of
   processors

5. IF ?MYREP .EQ. n THEN ....

   Works, too. Output only for the processor zero, unless OUTU is
   specified.

6. All the above works in parallel/parallel mode, ie each replica can
   be a parallel job in itself. The numeration of input and output
   files follows the replica numbers.

   The output is written only on a local process 0 for each replica,
   and similar is true also for stream command.  The limitation is
   that the number of replicas must divide the number of processes
   allocated for parallel. Otherwise it bombs out with the level -5.



File: Repdstr -=- Node: Fast
Up: Top -=- Previous: I/O -=- Next: Examples


The FAST keyword turns on fast replica exchange. When this option is
activated, all exchange decisions are made on processor 0, and new
temperatures are sent to each individual replica, as opposed to sending
coordinates and velocities between replicas. This method requires
substantially less communication, especially when the REPEat keyword
is used. The drawback is that the outputs are per-replica rather than
per-temperature. Additional functionality has been added to the
MERGE command (see dynamc.doc) to convert per-replica trajectories
to per-temperature.

Additionally, since all decisions are made on processor 0, only a single
exchange file is written out, showing the results of ALL exchanges at
a particular timestep. If REPEat is used the LOGLevel may be specified
to limit the number of exchanges written to this file. If LOGLevel is set
to N, every Nth exchange at a given step is written, however, the first
and last exchange is always written. NB, this only applies when the
REPEat key word is used, and at least one exchange at each time step
is always written.

Currently, FAST is only compatible with temperature replica exchange
without the use of a reservoir.



File: Repdstr -=- Node: Examples
Up: Top -=- Previous: Fast -=- Next: Output


NOTE: If you are using mpich-1.2.X then you need to use -p4wd with the
      absolute path or -p4wd `pwd`


Example 1:
==========

read psf
read coor
repd nrep 4

This will replicate PSF and coordinates, so after nrep 4 there are
four independent runs with the same coordinates

Example 2:
==========

read psf
repd nrep 4
read coor name system.crd

This will replicate PSF but the coordinates will be read from 4
separate files: system.crd_0, system.crd_1, etc


Example 3:
==========

repd nrep 4
stre inp

This will run for independent CHARMM jobs. Each inp_0, inp_1, inp_3,
and inp_4 can be different input files, with different PSFs,
parameters, etc

Example 4:
==========

open write unit 1 card name out

repd nrep 4
outu 1
stre inp

The same as example 3 but now also output files out_0, out_1, ... will
be written. Note that OUTU must precede STREam command.

Example 5: RXSGLD
==========

read psf
read coor name system.crd

!All stages have the same temperature of 300 K but have TEMPSG from 300 K to 500 K.
repd nrep 8 EXCHange FREQuency 1000 STEMp 300 DTEMp 0  -
  SGLD  SGTE 300 MSGT 500  DSGF 0.2

SCAL FBETA SET 1.0 SELE ALL END

!Perform SGLD with SGFT set to 0 to allow above RXSGLD setting in control
DYNA LANG SGLD SGFT 0






File: Repdstr -=- Node: Output
Up: Top -=- Previous: Examples -=- Next: Top



     The replica exchange printout is written to the unit specified in
the command after the UNIT keyword. The output of the current results 
is labeled by either REX> for temperature based replica exchange or 
RXSG> for self-guding replica exchange (RXSGLD). The RXSG> line contains the
following fields:
RXSG> Exchanges DynSteps StagID NeighborID ReplicaID Ep EpNeighbor 
TempScale TSGScale AcceptRatio Acceptance

A summary from all other 
replicas is labeled by REXSUM>. The labels in the output are
shortened, end the meaning of some of them is as the following:

Epot     - potential energy (current)
Tscale   - temperature scaling for exchange [Tscale=sqrt(Temp/NewTemp)]
Sratio   - success ratio [Srate=#-of-successful-exchanges/#-of-tried-exchanges]
NewTemp  - new temperature after the exchange
CurrTemp - current temperature
PROB     - probability to perform exchange P=exp(-Delta(1/kT)*Delta(Epot))
Rand     - random number used for exchange condition PROB>Rand => Success=T
NEIGHBOR - current neighbor with which the exchange occurs (or not)


NIH Helix/Biowulf Systems
charmm.org Homepage

CHARMM Documentation / Rick_Venable@nih.gov