[gmx-users] Analyse REMD
YanhuaOuyang
15901283893 at 163.com
Thu Jun 2 04:43:47 CEST 2016
Thank you for your detailed steps!
> 在 2016年5月31日，下午12:53，Christopher Neale <chris.neale at alum.utoronto.ca> 写道：
>
> 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.
> --
> 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