[gmx-users] Analyse REMD
Christopher Neale
chris.neale at alum.utoronto.ca
Tue May 31 06:53:13 CEST 2016
Dear Ouyang:
There is a script in the gromacs package (at least up to 5.1.2) called demux.pl . You can use it to get more human friendly information out of your REMD runs ( see http://www.gromacs.org/Documentation/How-tos/REMD for some very sparse info on demux.pl ). Basically you feed it the log file and it provides you with some useful output files. You can then use the -demux option to trjcat to get demuxed xtc files. The demux.pl script is not perfect. The last time I looked it could not handle more than X frames (I forget how many, but I and others posted fixes to the gmx mailing list a few years ago). Also, it relies on a sensible .log file and if you (1) run REMD then (2) have some type of crash or walltime overrun that does not give you .cpt files then (3) you restart REMD and it runs fine past the point of the previous crash, then demux.pl will give you good data. However, if you go through steps 1 and 2 above and then you restart REMS and it crashed again before you get to the point of the first crash then even if you again restart and you pass this point the log file will never be OK for demux.pl . I have some scrips that can recover from such an error but they are generally not necessary unless you have a run on a cluster that is unstable over some period. I have pasted my fix script below for posterity, though you very likely don't need it.
$ cat checkforproblems.sh
# This script finds and fixes any problems in the log file prior to demuxing
# When this script is done:
# if the ./problems file is empty then everything was fixed (or no problems ex isted)
# if the ./problems file is NOT empty, then some problem remain
#
# Command line input parameters:
# $1 = the number of MD integration steps between exchange attempts
# $2 = the location of the log file to use
STEP=$1
LOG=$2
function find_problems {
# $1 = log file
grep "^Replica exchange at step" $1 |awk '{if($5!=last+'$STEP')print $5,last;l ast=$5}' > problems.init
if (($(cat problems.init|wc -l)!=0)); then
cat problems.init | awk '{printf("%d\t",$2);if(NR%2==0)printf("\n")}' > prob lems
else
rm -f problems
touch problems
fi
grep "^Replica exchange at step" $1 |awk '{print $5,last;last=$5}' > z
}
function fix_problems {
# $1 = log file input
# $2 = log file output
# $3 = low cut
# $4 = high cut
# $5 = step
cat $1 |awk '
BEGIN{
part=0;
step='$5'
foundlow=0
foundhigh=0
printjustonemore=0
printnextround=0
##print "Fixing with:","'$1'","'$2'",'$3','$4','$5' >> "./FIX_PROBLEMS_INPUT "
}
{
if(part==0){
if($1=="Replica"&&$2=="exchange"&&$3=="at"&&$4="step"){
part=1
n=$5
if(n=='$3'){
foundlow++
}
if(n=='$4'){
foundhigh++
}
if(foundlow==0){
# output if we did not yet find the low cutoff
printit=1
}else{
# found the low cutoff
if(foundhigh==0){
# want to print the low cutoff once, so use the printjustonemore var iable to do so, then reset printit variable after output
printjustonemore=1
}else{
# Now only ouput if we have seen the first occurrence of the high cu toff, but dont want this value, so use the printnextround variable
printnextround=1
}
}
}
}else{
if($1!="Repl"){
part=0
# Print a sinfle blank line after each replica segment
print ""
}
}
if(part){
if(printit){
print $0
if(printjustonemore){
printjustonemore=0
printit=0
}
}
if(printnextround){
printit=1
}
}
}' > $2
}
# Step 0: cleanup
rm -f _ERROR_ problems* z* fixed.log*
# Step 1: Find any problems
# this tr command is necessary in case there was corruption in the .log file
cat $LOG | tr -d '\0' > my.log
find_problems my.log
# Step 2: Fix any problems
n=0
cp problems problems.$n
cp problems.init problems.init.$n
cp z z.$n
cat problems.0 | while read line; do
# exit the loop if there are no more problems
if (($(cat problems|wc -l)==0)); then
break
fi
let "n++"
fn=($line)
fix_problems my.log fixed.log ${fn[0]} ${fn[1]} $STEP
cp fixed.log fixed.log.${n}
cp fixed.log my.log
find_problems fixed.log
cp problems problems.$n
cp problems.init problems.init.$n
cp z z.$n
done
# Want this file fixed.log to always exist so that I can use it in the DEMUX scripts
if [ ! -e fixed.log ]; then
mv my.log fixed.log
else
rm my.log
fi
# Find any problems
find_problems fixed.log
________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of YanhuaOuyang <15901283893 at 163.com>
Sent: 22 May 2016 04:48:53
To: gmx-users at gromacs.org
Subject: [gmx-users] Analyse REMD
Hi, I have run a REMD, which including 46 replicas. I searched how to analyze my REMD statistics, but I am still confused how to analyze my result; how to generate a continuous a energy distribution curve; how to deal with 46 .log files. It seems so complicated.
one of my md.log file is as below:
Replica exchange at step 998000 time 1996.00000
Repl 18 <-> 19 dE_term = 1.690e+01 (kT)
Repl ex 0 1 2 3 4 x 5 6 x 7 8 x 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
Repl pr .01 .17 .03 1.0 1.0 .00 .00 .01 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00 .00
Replica exchange statistics
Repl 499 attempts, 250 odd, 249 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 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
Repl .07 .08 .05 .05 .06 .06 .06 .06 .06 .05 .05 .04 .03 .05 .03 .03 .03 .01 .03 .02 .04 .02 .02 .02 .01 .02 .02 .01 .01 .01 .01 .01 .01 .01 .01 .02 .01 .01 .02 .00 .02 .00 .02 .02 .01
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 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
Repl 14 20 13 12 16 14 15 17 14 15 9 12 8 13 7 5 7 2 5 6 10 6 4 6 2 5 3 1 6 3 3 3 2 3 4 5 2 3 4 2 4 2 4 4 3
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 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
Repl .06 .08 .05 .05 .06 .06 .06 .07 .06 .06 .04 .05 .03 .05 .03 .02 .03 .01 .02 .02 .04 .02 .02 .02 .01 .02 .01 .00 .02 .01 .01 .01 .01 .01 .02 .02 .01 .01 .02 .01 .02 .01 .02 .02 .01
I wonder whether the result is ideal, do I need to optimize my temperature distribution.
Best regards
ouyang
--
Gromacs Users mailing list
* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users
mailing list