Admin: Created page with "To enable logging immediately on lonestar, "setenv LOG_BUFS 0". Otherwise the log may be empty when your job dies. Put in brun. == gromacs instruction == http://www.gromacs.org/WIKI-import/REMD Server: http://folding.bmc.uu.se/remd/ paper: http://www.rsc.org/ej/CP/2008/b716554d.pdf '''GROMACS Procedure''' (http://www.gromacs.org/WIKI-import/REMD) The setup of a REMD simulation is actually quite straightforward. The following describes the steps that lead to hav..."

2022-05-10T15:57:59Z

Created page with "To enable logging immediately on lonestar, "setenv LOG_BUFS 0". Otherwise the log may be empty when your job dies. Put in brun. == gromacs instruction == http://www.gromacs.org/WIKI-import/REMD Server: http://folding.bmc.uu.se/remd/ paper: http://www.rsc.org/ej/CP/2008/b716554d.pdf '''GROMACS Procedure''' (http://www.gromacs.org/WIKI-import/REMD) The setup of a REMD simulation is actually quite straightforward. The following describes the steps that lead to hav..."

New page

To enable logging immediately on lonestar, "setenv LOG_BUFS 0". Otherwise the log may be empty when your job dies. Put in brun.

== gromacs instruction ==

http://www.gromacs.org/WIKI-import/REMD

Server: http://folding.bmc.uu.se/remd/

paper: http://www.rsc.org/ej/CP/2008/b716554d.pdf

'''GROMACS Procedure''' (http://www.gromacs.org/WIKI-import/REMD)

The setup of a REMD simulation is actually quite straightforward. The following describes the steps that lead to have a REMD simulation running on a given system. The "success" of the simulation will depend entirely on the problem being addressed and the criterion used to judge it. Although REMD simulations are helping in increasing sampling they do not provide the ultimate answer. This should be kept in mind.

Once the peptide and its surrounds have been determined, need to weigh up the range of temperature space to be sampled, the number of processors prepared to use and the length of time to simulate. With all versions of GROMACS 3.3 (i.e. prior to 4.0), only one processor per replica is allowed, as an additional constraint. The system, number of replicas, range of temperature space and the distribution of temperatures determine the average exchange probability between each replica. These should be constant across the replicas. To achieve this given the expected increase in system potential energy with increasing temperature, in the absence of known bottlenecks in the free energy space, the replica temperatures should be distributed exponentially.

There is good evidence to suggest that the period between exchange attempts for peptide/protein systems should not be under 1ps. This will determine the average length of time a replica will spend exploring conformation space between exchanges, and this is more important than the actual exchange probability. In the GROMACS implementation, an exchange attempt of a particular pair will only happen at every second attempt, as the odd and even pairs attempt exchange on alternating occasions. See the GROMACS Manual. General Steps

1. define the system eg: peptide + solvent (implicit or explicit).
2. depending on the number of processors available and the range of temperature to sample (they are actually extremely dependent on each other) choose a distribution of the temperatures. Use an exponential distribution: Ti=T0*ek*i where k and T0 can be tuned to obtain reasonable temperature intervals to allow exchanges. The exponential allows an increase of temperature interval as the temperature increases. This is necessary because the distribution of the total energy increase with the temperature and thus the exchange rate increases. Keep the exchange rate constant with the temperature.
3. once have the temperature distribution, equilibrate the systems at the N temperatures separately (using a separate .mdp file for each). Then construct a series of N run input files, .tpr from the equilibrated structures, again using different .mdp files to generate the different .tpr files. If have N temperatures should have N .tpr files, named prefix_0.tpr, ... prefix_N-1.tpr.
4. then run short simulations to have an estimate of the exchange rate (can get a good estimate within ~100 ps) and modify the temperatures if it does not correspond to the wished exchange. Typically 0.2 to 0.3 is good. What is actually most important than the value of the exchange rate itself is the time that a replica will spend at a give temperature. This is given by the exchange rate x frequency of exchange. eg: an exchange rate of 0.2 with exchange trials every 2 ps gives an average time at a temperature of 10 ps. One point to consider is also the way exchanged. Try to exchange all pairs at each trial or only one pair chosen randomly. In the latter case should also take this into account, the residence time of replica has then to be multiplied by N-1, the number of pairs exchanged.
5. the initial conformations at each temperature may be identical or different. This a choice and will depend on the reason for performing the REMD.

A web-server that implements the step of choosing temperatures for T-REMD can be [found here]. It works based on known energy distributions for solvated proteins. You input the number of protein atoms, the number of water atoms, the temperature range and the desired exchange probability and then the algorithm generates temperatures. Execution Steps

1. make a set of .mdp files, each specifying a different temperature being used. Number the output .tpr files according to index, from 0 to whatever (i.e., if have 10 different temperatures have prefix_0.tpr, prefix_1.tpr ... prefix_9.tpr). Prior to GROMACS 4.0, only a single processor may be used per replica, so either omit the -np flag to grompp or use -np 1. For GROMACS 4.0 no gromacs requires no -np flag.
2. the interval for replica exchange is given by the mdrun flag -replex. The number of cores (and the -np flag for mdrun for GROMACS prior 4.0) must agree with the number of .tpr files (i.e., 10 for the above general example using prefix_0.tpr through prefix_9.tpr). Nomenclature for the input files is critical to getting mdrun to work.
3. For a command line example, see below. Number of steps is the number of steps required to give you the exchange attempt period desired. Other command line options may be required, e.g. -o prefix.

GROMACS 3.0: mpirun -np 10 mdrun -s prefix_.tpr -np 10 -multi -replex (number of steps) (followed by output options) 

GROMACS 4.0: mpirun -np 10 mdrun -s prefix_.tpr -multi (integer) -replex (number of steps) (followed by output options)

*With -replex replica exchange is attempted every given number of steps.
*With -multi multiple systems are simulated in parallel. As many input files are required as the number of systems. The system number is appended to the run input and each output filename, for instance topol.tpr becomes topol0.tpr, topol1.tpr etc. The number of nodes per system is the total number of nodes divided by the number of systems. 

Understanding REMD-related Output

mdrun writes important information to your md.log file, you can extract it like this:

% grep Repl md0.log gives:

Initializing Replica Exchange Repl There are 6 replicas: Repl 0 1 2 3 4 5 Repl T 300.0 350.0 410.0 480.0 560.0 650.0 Repl Repl exchange interval: 1000 Repl random seed: 525106 Repl below: x=exchange, pr=probability Replica exchange at step 1000 time 2 Repl 0 <-> 1 dE = 2.492e+00 Repl ex 0 1 2 3 4 5 Repl pr .08 .00 .13 Replica exchange at step 2000 time 4 Repl ex 0 1 2 3 x 4 5 Repl pr .00 1.0

 This is indeed a very short explanation, but it should hopefully provide enough information to comprehend the results. Post-Processing

If you would like to make your trajectories continuous again you can use the script src/contrib/scripts/demux.pl that will read your md0.log file (you can concatenate several if necessary) and produce a few output files. One of these is a .xvg file (replica_ndx.xvg) that can be passed to trjcat, along with the original trajectory files, in order to produce continuous trajectories. The other file (replica_temp.xvg) contains the temperatures for each replica, starting at the original temperature. So if your replica of interest starts at, say, 300 K, you can follow its trajectory through temperature space. It would be interesting to add some functionality to make histograms of temperature distributions for each replica, which according to most authors, should be flat. The demuxed trajectories have been used in an application (implemented in g_kinetics - not available in GROMACS 3.3.1 or earlier) to obtain protein folding kinetics from REMD trajectories Phys. Rev. Lett. 96, 238102 (2006).

== Jenny: Gromacs replica exchange ==

1. run a 3 ns annealing to raise T from 300K to 800K. This is to generate a range of structures

2. use remdgen.pl to create 70 structures (remdxxx.gro) from the ann.trj from above. In GROMACS call on last line, "grp" input file contains only the number "0". This indicates "group 0: System" for the structure dump.

#!/usr/bin/perl -w
for ($i=0;$i<=69;$i++) {
$fr=$i*40;# one frame every 40 ps
$ogrom="remd"."$i.gro"; #output structure
`trjconv -s ann.tpr -f ann.trj -dump $fr -o $ogrom < grp`;
}

grp is a file with "0" in it

3. Use the following perl to prepare corresponding tpr files (remdxxx.tpr)

#!/usr/bin/perl -w
$tmp=280;
for ($i=0; $i<=72;$i++){
$tmp=1.0/(1/$tmp-0.00003*300**0.5/$tmp**0.5);
# $tmp=1.0/(1/$tmp-0.00003); # temperature gap between simulation following Gaussian
$f="remd"."$i.mdp";
`cp remd.mdp $f`;
`replace.pl TTT $tmp $f`;
$f2="remd"."$i.tpr";
$f3="remd"."$i.gro";
` grompp -v -f $f -c $f3 -o $f2 -p cmat.top `;
}

The perl script needs "replace.pl" in your path and a "remd.mdp" (below). Modify accordingly:

title = replica exchange
cpp = /usr/bin/cpp
constraints = all-bonds
;define = -DPOSRES
integrator = md
dt = 0.002  ; 2 fs time step !
nsteps = 2500000  ; total 4 ns.
;init_step = 2000
;disre = simple
;disre_fc = 5000
; freezing polymer matrix
;pbc = full
;freezegrps = pyu pyc CL
;freezedim = Y Y Y Y Y Y Y Y Y
nstcomm = 10
nstxout = 1000
nstvout = 1000
nstfout = 0
nstlog = 1000
nstenergy = 1000
nstlist = 10
ns_type = grid
rlist = 0.9
rcoulomb = 0.9
rvdw = 0.9
; Berendsen temperature coupling is on in two groups
;Tcoupl = nose-hoover
Tcoupl = Berendsen
tc-grps = System
tau_t = 0.5
ref_t = TTT
; Energy monitoring
energygrps = SYSTEM
;ergygrp_excl = pep pep
; Isotropic pressure coupling is now on
Pcoupl = no
;Pcoupl = Berendsen
;tau_p = 20
;pcoultype = isotropic
;ref_p = 100
;compressibility = 4.5E-5
;; Generate velocites is off at 300 K.
gen_vel = no
gen_temp = no
gen_seed = -1
andersen_seed = 371854
;pme
coulombtype = pme
fourierspacing = 0.10
pme_order = 4
ewald_rtol = 1e-05
epsilon_surface = 0
optimize_fft = yes

4. To submit job to tacc (or other cluster), "bsub < brun"

brun: Note remd matches the gro and tpr files prepared above.

#!/bin/sh
#BSUB -n 70
#BSUB -J remd
#BSUB -o 1.o%J
#BSUB -q normal
#BSUB -W 24:00
module load gromacs
#setenv GMXLIB /home/utexas/bio/prentacc/ANPtop
ibrun mdrun_mpi -replex 500 -s remd -o remd -c remd.gro -e -v -g remd -np 70

The interval should be at least 1 ps (500 steps with 2 fs step).

'''Xiaojia's script (7/2014)'''

My lonestar script 1(using gromacs4.5.3):
#!/bin/tcsh
#$ -V
#$ -cwd
#$ -N remd_9
#$ -q normal
#$ -j y
#$ -o remd9.out
#$ -e remd9.err
#$ -A TG-MCB100057
#$ -l h_rt=24:00:00
#$ -m be
#$ -pe 16way 1920
#$ -hold_jid 2456421 Note this line is for continuing different scripts
setenv BUF_LOGS 0
module load gromacs/4.5.3
ibrun mdrun_mpi -s -deffnm remd -multi 120 -replex 250 -cpi -append Note this line is for Gromacs remd simulation

5. To restart a REMD, use the following script (change 69 to however many runs you have):

#!/usr/bin/perl -w

for ($i=0; $i<=69;$i++){
$trr="remd"."$i.trr";
$tpr="remd"."$i.tpr";
$otpr="rc"."$i.tpr";
` tpbconv -s $tpr -f $trr -o $otpr`;
}

 Then submit the job with "remd" replaced by "rc".

6. '''At the bottom of the log file (remd0.log), the replica exchange statistical is printed. Ideally you want the probability of switching uniform across different Ts, and > 20% at least.'''

7. to join two trjactories in case you restarted:

cp remd6.trr 298-1.trr
cp rc6.trr 298-2.trr
module load gromacs
trjcat -f 298-1.trr 298-2.trr -o 298.trr
rm 298-1.trr 298-2.trr

== Additional REMD Scripts from Fonner ==

Using this website (also listed at the tope of the page) [http://folding.bmc.uu.se/remd/], you can get a list of comma separated temperatures for your REMD simulation based on your system size, temperature range, and the acceptance ratio you want. You can save that list to a file and use this perl script to generate a set of .gro files. This script is also located in /home/johnfonner/bin/remdExtract.pl. A couple of notes on usage: 1.) type remdExtract.pl -h for a reminder on the input formatting. 2.) Make sure your trajectory only uses a given reference value once. If there is more than one instance of a temperature (like you heated and then cooled your system), the script will take a snapshot of the first instance only (which is probably not the one you want in a full annealing trajectory. If this affects you, I would just clip the first part of your trr file with the command <code>trjconv_mpi -f traj.trr -s tprfile.tpr -o traj_out.trr -b 1000</code> or something like that.

#!/usr/bin/perl

# Input: Gromacs md.log, traj.trr file, and a comma separated list of temperatures
# Output: series of .gro file snapshots taken at the temperatures in the input file

# PROCESS INPUT

while ($ARGV = shift(@ARGV)) {
if ($ARGV eq '-h') {
die "Usage: remdExtract.pl -t temps.csv -r traj.trr -p tpr -g md.log \n";
}
elsif ($ARGV eq '-g') {
open(LOG, shift(@ARGV))||die "can't open log file\n";
@mdlog=<LOG>;
close(LOG);
}
elsif ($ARGV eq '-t') {
open(TEMP, shift(@ARGV))||die "can't open temperature file\n";
@temps=sort { $a <=> $b } split(',',<TEMP>);
close(TEMP);
}
elsif ($ARGV eq '-p') { $tpr = shift(@ARGV); }
elsif ($ARGV eq '-r') { $traj = shift(@ARGV); }
else {die "could not process $ARGV\n";}
}
die "No log file given. \n" unless (@mdlog);
die "No temperature file given. \n" unless (@temps);
$traj = "traj.trr" unless ($traj);
die "Can't find trajectory file. \n" unless (`ls $traj`);
die "Can't find tpr file. \n" unless (`ls $tpr`);

# trim leading whitespace from md.log file because it causes trouble
foreach (@mdlog) {s/^\s+//;}

# search process
for ($i = 0; $i <= $#mdlog; $i++) {
if ($mdlog[$i] =~ /Current ref\_t for group/) {
for ($j = $i+2; $j < $i + 8; $j++) {
@tokens2 = split ('\s{2,}',$mdlog[$j]);
for ($k = 0; $k <= $#tokens2; $k++) {
if ($tokens2[$k] =~ /Temperature/) {$thisTemp = (split(' ',$mdlog[$j+1]))[$k];}
}
}
$time = (split(' ',$mdlog[$i-2]))[1];
for ($l=0; $l <= $#temps; $l++) {
if (abs($thisTemp - $temps[$l]) < 0.2) {
$grofile = "remd"."$l".".gro";
$temps[$l]=0; #zero out this temperature so that it isn't used again but it still holds its place in the array.
# "the echo command tells trjconv to use the "0" group, which is "System". Output is suppressed by redirecting to /dev/null.
`echo "0\n" | trjconv_mpi -f $traj -s $tpr -o $grofile -b $time -e $time -sep 2> /dev/null`; #the "-sep option makes sure you see it if you get multiple frames
print "$grofile created with reference temperature $thisTemp at timestep $time\.\n";
}
}
}
}

Ammm:REMD jenny - Revision history