Merge branch 'v2.4'

CHANGES/v2.3.md

@@ -218,6 +218,10 @@ For developers:
 
 ## Version 2.3.5
 
+For users:
+- Fixed `plumed partial_tempering` to agree with GROMACS conventions for choice of dihedrals (see \issue{337}).
+  Should be irrelevant for the vast majority of cases.
+
 For developers:
 - Doxygen on travis-ci has been updated to 1.8.14
 - `make clean` now correctly removes the `src/lib/plumed` executable.

regtest/basic/rt-hrex-0/out.reference

@@ -578,8 +578,8 @@ N   C   N            1   120.000    585.760  ; Added for Urea (same as N2-CA-N2)
 ; LINE(512) X   X   CM  HA       4      180.00     4.60240     2	 ; 
 ; LINE(513) X   X   CA  H4       4      180.00     4.60240     2	 ; bsd.on C6H6 nmodes 
 ; LINE(514) X   X   CA  H5       4      180.00     4.60240     2	 ; bsd.on C6H6 nmodes
-
-
+[ dihedraltypes ] 
+; LINE(516) X   X   X   X     9       0.0      1.60247     3   ; fake choice, just here to test gromacs 5.1.4 behavior
 [ dihedraltypes ] 
 ;i   j   k   l	   func	
 ; LINE(519)  CT  CT  OS  CT    9       0.0      1.60247     3   ; 

regtest/basic/rt-hrex-0/processed.top

@@ -512,8 +512,8 @@ X   X   CM  H4       4      180.00     4.60240     2	;
 X   X   CM  HA       4      180.00     4.60240     2	; 
 X   X   CA  H4       4      180.00     4.60240     2	; bsd.on C6H6 nmodes 
 X   X   CA  H5       4      180.00     4.60240     2	; bsd.on C6H6 nmodes
-
-
+[ dihedraltypes ]
+X   X   X   X     9       0.0      1.60247     3  ; fake choice, just here to test gromacs 5.1.4 behavior
 [ dihedraltypes ]
 ;i   j   k   l	   func	
  CT  CT  OS  CT    9       0.0      1.60247     3  ; 

scripts/partial_tempering.sh

+# 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]);