From 315e67461312a4cf1cad5e635292617b0d6d4ab7 Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Mon, 5 Feb 2018 19:40:08 +0100
Subject: [PATCH] Fixed #337

Fixed partial_tempering with gromacs 5 convention. Old behavior can be
forced with --gromacs4 option.
---
 scripts/partial_tempering.sh | 67 ++++++++++++++++++++++++------------
 1 file changed, 45 insertions(+), 22 deletions(-)

diff --git a/scripts/partial_tempering.sh b/scripts/partial_tempering.sh
index e3dc15dcf..2e5f1832b 100755
--- a/scripts/partial_tempering.sh
+++ b/scripts/partial_tempering.sh
@@ -1,3 +1,4 @@
+# vim:ft=awk
 if [ "$1" = --description ] ; then
   echo "create a new collective variable from a template"
   exit 0
@@ -7,7 +8,7 @@ if [ "$1" = --help ] || [ "$1" = -h ] ; then
   cat <<EOF
 Usage:
 
-  plumed partial_tempering scale < processed.top
+  plumed partial_tempering [--gromacs4] scale < processed.top
 
 where scale is the Hamiltonian scaling factor and
 processed.top is a post-processed topology file (i.e. produced with grompp -pp)
@@ -55,7 +56,14 @@ EOF
   exit
 fi
 
-awk -v scale=$1 '
+gromacs5=1
+
+if [ $1 == --gromacs4 ] ; then
+  gromacs5=0
+  shift
+fi
+
+awk -v scale=$1 -v gromacs5=$gromacs5 '
 BEGIN{
   combrule=1;
 }
@@ -73,6 +81,39 @@ function warning(msg)
 {
      print "WARNING:",msg | "cat 1>&2"
 }
+function find_matching_torsion(params,atype, a1,a2,a3,a4,iswitch,progression,test,array,param,countX,bestCountX,bestMatch)
+{
+   progression=NR
+   bestCountX=5
+   for(iswitch=0;iswitch<32;iswitch++){
+     countX=0
+     if(iswitch%2==0){
+       a1=atype[1]; a2=atype[2]; a3=atype[3]; a4=atype[4];
+     } else {
+       a1=atype[4]; a2=atype[3]; a3=atype[2]; a4=atype[1];
+     }
+     if(int(iswitch/2)%2==1){ a1="X"; countX++; }
+     if(int(iswitch/4)%2==1){ a2="X"; countX++; }
+     if(int(iswitch/8)%2==1){ a3="X"; countX++; }
+     if(int(iswitch/16)%2==1){a4="X"; countX++; }
+     test=a1"-"a2"-"a3"-"a4"-"$5;
+     if(test in params){
+       split(params[test],array,":");
+       bestMatch=0;
+       if(gromacs5) {
+         if(countX<bestCountX || (countX==bestCountX && array[1]<progression)) bestMatch=1
+       } else {
+         if(array[1]<progression) bestMatch=1;
+       }
+       if(bestMatch){
+         progression=array[1];
+         bestCountX=countX
+         param=params[test];
+       }
+     }
+   }
+   return param;
+}
 {
 # This is the suffix for "hot" atoms:
   suffix="_";
@@ -185,27 +226,9 @@ function warning(msg)
        atype[2]=bondtype[ato[$2]]
        atype[3]=bondtype[ato[$3]]
        atype[4]=bondtype[ato[$4]]
+
+       param=find_matching_torsion(params,atype);
        
-       progression=NR
-       for(iswitch=0;iswitch<32;iswitch++){
-         if(iswitch%2==0){
-           a1=atype[1]; a2=atype[2]; a3=atype[3]; a4=atype[4];
-         } else {
-           a1=atype[4]; a2=atype[3]; a3=atype[2]; a4=atype[1];
-         }
-         if(int(iswitch/2)%2==1) a1="X";
-         if(int(iswitch/4)%2==1) a2="X";
-         if(int(iswitch/8)%2==1) a3="X";
-         if(int(iswitch/16)%2==1) a4="X";
-         test=a1"-"a2"-"a3"-"a4"-"$5;
-         if(test in params){
-           split(params[test],array,":");
-           if(array[1]<progression){
-             progression=array[1];
-             param=params[test];
-           }
-         }
-       }
        
     n=split(param,array,":");
     if(n<=1) error("params not found "$1" "$2" "$3" "$4" "$5" "atype[1]" "atype[2]" "atype[3]" "atype[4]);
-- 
GitLab