Fix a handful of issues -- Reference Amoeba plugin now builds and finishes the

stretch-bend test correctly.

The potential energy calculation done in the ReferenceStretchBendTest file was
updated to mirror that of the Reference kernel, but first I made sure that the
test passed *before* updating to the new energy expression, just to make sure
that the energies and forces didn't *change* with the new expression.

plugins/amoeba/openmmapi/include/openmm/AmoebaStretchBendForce.h

@@ -94,24 +94,6 @@ public:
     void getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3, double& lengthAB,
                                   double& lengthCB, double& angle, double& k1, double& k2) const;
 
-    /**
-     * Get the force field parameters for a stretch-bend term, provided for
-     * backwards-compatibility. Since the two force constants may be different,
-     * you are recommended to call the other version of getStretchBendParameters
-     * with both a k1 and k2
-     * 
-     * @param index         the index of the stretch-bend for which to get parameters
-     * @param particle1     the index of the first particle connected by the stretch-bend
-     * @param particle2     the index of the second particle connected by the stretch-bend
-     * @param particle3     the index of the third particle connected by the stretch-bend
-     * @param lengthAB      the equilibrium length of the stretch-bend in bond ab [particle1, particle2], measured in nm
-     * @param lengthCB      the equilibrium length of the stretch-bend in bond cb [particle3, particle2], measured in nm
-     * @param angle         the equilibrium angle in radians
-     * @param k1            the force constant
-     */
-    void getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3, double& lengthAB,
-                                  double& lengthCB, double& angle, double& k1) const;
-
     /**
      * Set the force field parameters for a stretch-bend term.
      * 

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -346,7 +346,7 @@ void ReferenceCalcAmoebaStretchBendForceKernel::copyParametersToContext(ContextI
         lengthCBParameters[i] = (RealOpenMM) lengthCB;
         angleParameters[i] = (RealOpenMM) angle;
         k1Parameters[i] = (RealOpenMM) k1;
-        k1Parameters[i] = (RealOpenMM) k2;
+        k2Parameters[i] = (RealOpenMM) k2;
     }
 }
 

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.cpp

@@ -171,7 +171,7 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( int numStre
         RealOpenMM angleK2      = k2Quadratic[ii];
         RealVec forces[3];
         energy                 += calculateStretchBendIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index],
-                                                           abLength, cbLength, idealAngle, angleK1, anglek2, forces );
+                                                           abLength, cbLength, idealAngle, angleK1, angleK2, forces );
         // accumulate forces
     
         for( int jj = 0; jj < 3; jj++ ){

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.h

@@ -94,7 +94,8 @@ private:
        @param lengthAB                ideal AB bondlength
        @param lengthCB                ideal CB bondlength
        @param idealAngle              ideal angle
-       @param kParameter              k
+       @param k1Parameter             k for distance A-B * angle A-B-C
+       @param k2Parameter             k for distance B-C * angle A-B-C
        @param forces                  force vector
     
        @return energy
@@ -104,8 +105,8 @@ private:
     RealOpenMM calculateStretchBendIxn( const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
                                         const OpenMM::RealVec& positionAtomC,
                                         RealOpenMM lengthAB,      RealOpenMM lengthCB,
-                                        RealOpenMM idealAngle,    RealOpenMM kParameter,
-                                        OpenMM::RealVec* forces ) const;
+                                        RealOpenMM idealAngle,    RealOpenMM k1Parameter,
+                                        RealOpenMM k2Parameter,   OpenMM::RealVec* forces ) const;
  
 };
 

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaStretchBendForce.cpp

@@ -83,14 +83,14 @@ static void computeAmoebaStretchBendForce(int bondIndex,  std::vector<Vec3>& pos
                                           std::vector<Vec3>& forces, double* energy, FILE* log ) {
 
     int particle1, particle2, particle3;
-    double abBondLength, cbBondLength, angleStretchBend, kStretchBend;
+    double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend;
 
-    amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend);
+    amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
     angleStretchBend *= RADIAN;
 #ifdef AMOEBA_DEBUG
     if( log ){
-        (void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k=%10.3e\n", 
-                             bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend );
+        (void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k1=%10.3e k2=%10.3e\n",
+                             bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend );
         (void) fflush( log );
     }
 #endif
@@ -149,12 +149,13 @@ static void computeAmoebaStretchBendForce(int bondIndex,  std::vector<Vec3>& pos
        deltaR[CBxP][ii] *= termC;
     }
  
-    double dr    = rAB - abBondLength + rCB - cbBondLength;
+    double dr1   = rAB - abBondLength;
+    double dr2   = rCB - cbBondLength;
  
     termA        = 1.0/rAB;
     termC        = 1.0/rCB;
  
-    double term  = kStretchBend;
+    double drkk = dr1 * kStretchBend + dr2 * k2StretchBend;
  
     // ---------------------------------------------------------------------------------------
  
@@ -166,8 +167,8 @@ static void computeAmoebaStretchBendForce(int bondIndex,  std::vector<Vec3>& pos
     double subForce[LastAtomIndex][3];
     double dt = angle - angleStretchBend;
     for( int jj = 0; jj < 3; jj++ ){
-        subForce[A][jj] = term*(dt*termA*deltaR[AB][jj] + dr*deltaR[ABxP][jj] );
-        subForce[C][jj] = term*(dt*termC*deltaR[CB][jj] + dr*deltaR[CBxP][jj] );
+        subForce[A][jj] = kStretchBend*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj];
+        subForce[C][jj] = k2StretchBend*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj];
         subForce[B][jj] = -( subForce[A][jj] + subForce[C][jj] );
     }
  
@@ -187,7 +188,7 @@ static void computeAmoebaStretchBendForce(int bondIndex,  std::vector<Vec3>& pos
     forces[particle3][1]       -= subForce[2][1];
     forces[particle3][2]       -= subForce[2][2];
 
-    *energy                    += term*dt*dr;
+    *energy                    += dt*drkk;
 #ifdef AMOEBA_DEBUG
     if( log ){
         (void) fprintf( log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr );