Update behavior for alchemical bonded scaling 34/4934/2
authorradakb <brian.radak@gmail.com>
Fri, 1 Feb 2019 19:11:09 +0000 (14:11 -0500)
committerradakb <brian.radak@gmail.com>
Fri, 1 Feb 2019 20:03:00 +0000 (15:03 -0500)
Molecule::get_fep_bonded_type now explicitly catches bonded terms
involving alchemical atoms that should NOT be scaled and assigns
them to group 0 (as opposed to 1 or 2). This may have been a quiet
bug for the case alchBondLambdaEnd > 0.0 and alchBondDecouple off,
neither of which are default and are not often (or never?) used.

Change-Id: I404703eadd30951270ff2a453e8cdff31c835bdc

src/Molecule.h

index 7a7dd80..638ba33 100644 (file)
@@ -1378,14 +1378,14 @@ public:
     // This really only applies when order = 2.
     if ( simParams->alchBondDecouple ) order++;
 
-    if ( typeSum == 0 ) return 0; // Most interactions get caught here.
+    // The majority of interactions are type 0, so catch those first.
+    if ( typeSum == 0 || abs(typeSum) == order ) return 0;
     else if ( 0 < typeSum && typeSum < order ) return 1;
-    else if ( 0 > typeSum && typeSum > -order ) return 2;
+    else if ( -order < typeSum && typeSum < 0 ) return 2;
 
-    if ( simParams->alchBondDecouple ) {
-      // Alchemify should always keep this from bombing, but just in case...
-      NAMD_die("Unexpected alchemical bonded interaction!");
-    }
+    // Alchemify should always keep this from bombing, but just in case...
+    NAMD_die("Unexpected alchemical bonded interaction!");
+    return 0;
   }
 //fepe