An added note on the ad hoc REX package in the Lecture 14 notes; I have put together an improved version of this package, which fixes a few minor issues, and generally reduces the amount of file editing required.
Please contact me for the NEW! IMPROVED! package, rather than attempting to use the script and progs found in the lecture notes.
I've translated the the C and Python code into a single Fortran program, with the swap step range, cleanup frequency, number of T baths, and the T values all read from a simple config file. The .csh scripts were re-written to use the config file as well. The program should only need to be compiled once for a particular computer system. Anyone interested in this greatly improved package should contact me via email.
Attached is a modified package which fixes a major bug not found in my initial tests (a user found it and reported to me). The symptom is seen in the log file as an energy of 0.00 for one of the configs, and a swap probability of INF (infinity). There are 2 changes for this, one in the Trex.f program, the other in the rexstart.inp and rex.inp CHARMM scripts. Other than recompiling the program, the other change for users is the following addition (in red) to the input scripts--
! KEEP THIS NEXT BIT open unit 12 write card name @TACT/rex.ene format (F20.6) write title unit 12 * ?ENER *
With this change, the energy values are transmitted between programs w/o errors (at least in all my tests so far).
I made a few other minor changes, esp. for the log file, which is now named based on the last swap step it will perform. The frxn-swap.csh script was also changed to reflect the log file name change.
Please check this thread for new info and updated packages; I am currently investigating changes to support the use of the NPT ensemble. For informational purpose, the latest AdHocTrex package README.txt follows:
The TREX ad hoc interface to CHARMM for REMD November 2006 Update
--------------------------------------------------------------------- N.B.: The initial release had a bug in passing energy values, and incorrect swaps could occur; this release fixes that bug and a couple other minor issues. The bug symptom is an infinite probability (seen as INF with g77), and an energy of 0.00 for one of the baths. ---------------------------------------------------------------------
An interface to CHARMM for temperature based replica exchange is described, using a custom Fortran program, and taking advantage of the dynamics restart capabilities. Some editing of the scripts is needed as part of the setup; however, the interface is fairly flexible for this same reason. The programs and scripts are presented in the form used to perform some test REMD simulations on gel phase bilayers of the lipid DMPC (40 T baths, 250 MD steps, 4000 swaps; 1 ns), and a single lipid test (10 baths, more widely spaced; 2 ns). This is an original translation of the C and Python code from the authors of "Optimal estimates of free energies from multi-state nonequilibrium work data", Paul Maragakis, Martin Spichty, and Martin Karplus, PHYS REV LETT 96 (10): Art. No. 100602 MAR 17 2006.
The Fortran program controls the synchronization of the MD simulations and swaps, and calls various other programs (including CHARMM). After each MD run, the final potential energy is used to evaluate the Metropolis criterion; for a swap, the CHARMM restart files are exchanged, and the velocities are scaled based on the T change. From dynamc.doc, the feature used to scale the velocities is:
---------------------------------------------------------------- ISCALE 0 This option is to allow the user to scale the velocities by a factor SCALE at the beginning of a restart run. This may be useful in changing the desired temperature. .eq. 0 no scaling done (usual input value) .ne. 0 scale velocities by SCALE. WARNING: Please use this option only when you are changing the temperature of the run. SCALE 1. Scale factor for the previous option. ----------------------------------------------------------------
The value of SCALE is 1.0 unless a swap occurs; for a swap, the value computed from the square root of the T ratio is used. The values are communicated to CHARMM by a formatted file, which is read via the STREAM command. The final energy (?ENER) from each MD run is written to a file using the WRITE TITLE feature.
REMD is embarassingly parallel; it is ideally suited to Linux clusters, as there is no communication between processes except via files between the short MD runs, and each process is self contained. Parallel MD, on the other hand, must exchange a fair amount of data between all processes on every integration step, esp. if particle-mesh Ewald (PME) summation is being used. It's simplest to use one processor for each T bath, although it is certainly possible and reasonable (but more complicated) to use both processors in a dual-processor compute node. However, for a fixed number of nodes, it may often be better to have twice as many T baths than to double the speed.
In order to help manage the output, a cleanup script is called every N swap steps from the Trex program, performing: [1] deletion of unneeded restart files, [2] gzip compression and relocation of output files, and [3] merge of the trajectories into larger files and deletion of the originals (via merge.inp). Without these steps, one can easily become overwhelmed by the number of files produced and the disk space used.
The following listing outlines the steps required to use this collection of scripts and programs to perform REMD with CHARMM:
On a given cluster, the program should only need to be compiled once in a while, for new machine types or possibly software updates. The remaining steps apply to each REMD simulation. * Choose the number of T baths and the spacing between baths; some iterative tests may be needed to get the optimum spacing and to maximize the T range for a fixed number of baths; integer T values are assumed * Decide how many MD steps to perform between swaps, and how many swap steps for a particular run; again, some initial tests may be needed. * Based on these choices, create the config file (see below) * Create a link named "charmm" in the working directory that points to a single threaded (non-parallel) version of CHARMM; on NIH Biowulf use something like
* Edit rexstart.inp, rex.inp, and init.str based on the number of MD steps and protocol, and the actual files (RTF, PARAM, PSF, COOR, etc.) that will be used; note that at least one coord set must be written to the .trj file. * Edit the merge.inp script, to set the SKIP interval for MERGE * Edit the run_rex.csh script if needed, e.g. for some other queuing system besides PBS, or for the location of the MPICH installation. * Run the rex_init.csh script; it sets up the subdirs for each T bath, and requires the config filename as an argument * Submit (PBS) the run_rex.csh script to start the replica exchange; the example passes 2 args, the number of baths and the config filename
Monitoring the swap acceptance
Also included are a couple means of monitoring the swaps and the acceptance rate:
* frxn-swap.csh; computes the overall acceptance, and that for each bath from data in the rexswapN.log files; config file arg, prints a table * smove-temp.py; a Python program, which must be edited (T bath list) prior to use; this produces data suitable for the 'xmgrace' program via
% smove-temp.py > aswap.dat % xmgrace aswap.dat
Restart; extending a run, or after a failure
Continuation after a successful completion is straightforward; simply change the first and last swap step in the config file, as indicated in the detailed description of the config file below. Then submit the run_rex.csh script to the PBS queue as before. If a run stops early because of a node failure, network disruption, or some other problem, it can also be continued, but a couple of additional steps are needed. The steps are:
[1] determine the last two file numbers, which should represent the swap step MD in progress, and the last completed swap step; I use
% ls -str 293
[2] verify the status for all baths; if the last 2 file numbers are e.g. 694 and 695, the command
% ls -s ???/rex.trj.69[45]
will produce a listing of the last two .trj files for all baths, and their size (in blocks); for the running step which failed, one or more of the rex.trj.695 files will have a zero size
[3] in order to preserve the log data, change the name of the log file to reflect the last completed step, e.g.
% mv rexswap4000.log rexswap694.log
[4] edit the config file, and change the first step number to that of the step that failed in progress (695 for this example)
[5] resubmit run_rex.csh to a PBS batch queue
Additional notes on CHARMM input scripts
The filenames rexstart.inp and rex.inp should probably be preserved in order to avoid breaking the package; however, there is a great deal of latitude for the contents of the setup (init.str) and the short dynamics run (rex*inp) scripts. The example subdirs indicate a couple choices, vacuum sampling and condensed phase sampling; the latter is more difficult and requires a fairly close T bath spacing. They may not be ideal examples and are not test cases, but are meant to illustrate the usage of the package. The number of dynamics steps and the T bath choices may not be ideal; I'm still expolring these issues myself. The merge.inp script currently prunes the trajectory files to 1 ps intervals for the lowest T bath, and to 5 ps intervals for all other baths. Other than changing the SKIP value for saving coord sets, this file should not require any further editing. The example1 subdir has another script which further merges the files to 1 ns sizes, in order to faciliate analysis of the results.
Configuration file format
My convention has been to use .cfg as the extension for the config file, but any extension may be used. The first 4 lines of the file specify key parameters; the remaining lines are the integer values for the T baths. For the first 4 lines, the first 5 chars of each is a label; the integer values should start in col 6 or 7. The data is read based on order; the first line should be the number of the first swap step for this particular run, etc. A simple example:
FIRST 1 LAST 4000 CLEAN 80 NBATH 5 300 330 370 420 480
The labels indicate:
FIRST the first swap step; 1 to start, e.g. 4001 to continue LAST the last swap step; 1+(LAST-FIRST) CHARMM runs of rex*.inp CLEAN frequency for running cleantemp.csh; must be < 90 NBATH the number of T bath values on the following lines
Note that rextart.inp is only run for the very first swap step, i.e. step 1 when FIRST is 1. To continue the above run, the config file would be changed as indicated below, and run_rex.csh re-submitted to the PBS queue.
FIRST 4001 LAST 8000 CLEAN 80 NBATH 5 300 330 370 420 480
File manifest
Trex.f master Fortran prog
rex_init.csh set up the T bath subdirs run_rex.csh script to submit to PBS queue for REX cleantemp.csh periodically called to cleanup rexdone.csh optional cleanup after completion frxn-swap.csh monitor swap acceptance smove-temp.py build data for moves plot (xmgrace)
example1 subdir illustrating vacuum sampling of a lipid dppc10.cfg sample config file; 10 T baths init.str CHARMM initialization; read some files rexstart.inp the first MD run; write .res and .ene files rex.inp all other MD runs; reads .res merge.inp merges .trj files (cleantemp.csh) postmerg.inp secondary merge to 1 ns chunks for analysis
example2 subdir illustrating condensed phase sampling dmpc40.cfg sample config file; 40 T baths init.str CHARMM initialization; read some files rexstart.inp the first MD run; write .res and .ene files rex.inp all other MD runs; reads .res merge.inp merges .trj files (cleantemp.csh)
README.txt this description Bibliography.txt 2004-2005 publications using replica exchange
