Hi Gromacs users,
I seem to face some problems with my REMD dynamics.
I follow as much as possible the set up description of:
http://jcp.aip.org/resource/1/jcpsa6/v135/i14/p145102_s1?isAuthorized=no
Just to summarise:
Alanine dipeptide in explicit water: minimization, nvt, and npt equilibration
at each temperature of my 24 replicas between 300K and 500K exponentially
spaced. The box volume is the same for each replica, when I do an nvt
production run.
I have my inputfiles named md*.mdp, where * goes from 0->24 i produce my tpr
files successfully.
I run a nvt REMD simulation for 20ns where I swap every 1ps. The resulting
exchange statistics seem reasonable to me. Here is an example from one of the
log files:
Replica exchange statistics
Repl 19999 attempts, 10000 odd, 9999 even
Repl average probabilities:
Repl 0 1 2 3 4 5 6 7 8 9 10 11 12 13
14 15 16 17 18 19 20 21 22 23
Repl .32 .32 .33 .32 .34 .35 .35 .37 .36 .38 .37 .39 .39 .40
.40 .42 .42 .43 .44 .44 .45 .45 .46
Repl number of exchanges:
Repl 0 1 2 3 4 5 6 7 8 9 10 11 12 13
14 15 16 17 18 19 20 21 22 23
Repl 3157 3173 3245 3244 3353 3404 3492 3684 3651 3752 3769 3854 3849 4074
4049 4102 4187 4250 4322 4349 4475 4519 4562
Repl average number of exchanges:
Repl 0 1 2 3 4 5 6 7 8 9 10 11 12 13
14 15 16 17 18 19 20 21 22 23
Repl .32 .32 .32 .32 .34 .34 .35 .37 .37 .38 .38 .39 .38 .41
.40 .41 .42 .43 .43 .43 .45 .45 .46
all good so far. Now to the analysis where the problem arises:
Here is my post production script in order to extract constant temperature
trajectories:
rm -f total.log
for j in {0..23}
do
cat md${j}.log >> total.log
done
echo "demuxing replica"
demux total.log
trjcat -f *.xtc -demux replica_index.xvg
echo "now getting the ramachandran diagrams"
for j in {0..23}
do
g_rama -f ${j}_trajout.xtc -s ../../tprInput/md${j}.tpr -o
rama${j}.xvg
done
This also runs through without any problems...
Now I construct a histogram based on my 6 metstable states in my Ramachandran
digram, as defined in:
http://jcp-bcp.aip.org/resource/1/jcpbcp/v5/i6/p06B613_s1?isAuthorized=no
I literally just count how many point in my diagram belong to states 1->6. This
should give me an estimate of my equilibrium distribution of each state at each
temperature (if i do the counting for each Ramachandran diagram at each
temperature). Now here comes the problem:
The probability for each temperature is the same, also the transition rates
between states at different temperatures changes only very slightly. I have
taken averages over 10 * 20ns REMD trajectories and my state probabilities look
much like this:
http://www.nottingham.ac.uk/~ppxasjsm/StatDist.png
This behaviour is much too flat if i compare these results to the two papers
mentioned above. Also at 300K data points from a 40ns equilibrium run were
included as well as at 500K the same, same states have the same colours. The
values for 300K correspond to the findings of the two papers. This equilibrium
run was started from the same mdp file as used for the REMD simulation, but
this time just a single temperature was run.
Thus I am not able to reproduce my equilibrium distribution from my REMD
simulations. Below is my sample input file.
Just as a sanity check I also tried to create the the distribution from each
replica file (i.e not building the continuous temperature trajectory) and I get
the same result...
Something else I noticed is that if I do a g_energy -f md0.edr on my output
files the average temperature I get from these is always that of the input
trajectory (this seems strange too me)
; Run parameters
integrator = sd ; leap-frog integrator
nsteps = 10000000 ; 20 ns
dt = 0.002 ; 2 fs
; Output control
nstxout = 0 ; don't save in trr format
nstvout = 0 ;
nstxtcout = 50 ; save xtc every 0.1ps
nstenergy = 50 ; save energies every 0.1 ps
nstlog = 2500 ; update log file every 5 ps
; Bond parameters
continuation = yes ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 5 ; 10 fs
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range
electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = System ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in
K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; continuation
energygrps = protein water ;
xtc_grps = Protein
gen_temp = 300
gen_seed = -1
ld_seed = -1
(I also tried md as an integrator. The seed set on different processors always
seems to be the same somehow..(that should only give me correlations in the
data i guess) I am not sure if I have to set gen_temp, but some tutorials
suggested this, so I did)
Just one last note. This is pretty much the first time I am using Gromacs and I
am obviously overseeing something, which I cannot find. I have written my own
parallel tempering code before, so I am very familiar with the concepts, just
not with the Gromacs software.
If the mystery could be solved this would be greatly appreciated, of if more
information is required I am happy to provide it.
Best,
Antonia
--
Antonia Mey
Freie Universität Berlin
FB Mathematik + Informatik
Institut für Mathematik
Arnimallee 6
D-14195 Berlin
Email: [email protected]
Tel: +49-30-83875879This message and any attachment are intended solely for the
addressee and may contain confidential information. If you have received this
message in error, please send it back to me, and immediately delete it.
Please do not use, copy or disclose the information contained in this message
or in any attachment. Any views or opinions expressed by the author of this
email do not necessarily reflect the views of the University of Nottingham.
This message has been checked for viruses but the contents of an attachment
may still contain software viruses which could damage your computer system:
you are advised to perform your own checks. Email communications with the
University of Nottingham may be monitored as permitted by UK legislation.--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
