Merge pull request #801 from swails/strbnd

Generalize Amoeba Stretch-Bend term

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

@@ -71,11 +71,12 @@ public:
      * @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 k             the force constant for the stretch-bend
+     * @param k1            the force constant of the product of bond ab and angle a-b-c
+     * @param k2            the force constant of the product of bond bc and angle a-b-c (optional, default is the same as k1)
      * @return the index of the stretch-bend that was added
      */
     int addStretchBend(int particle1, int particle2, int particle3, double lengthAB,  double lengthCB, double angle,
-                       double k );
+                       double k1, double k2);
 
     /**
      * Get the force field parameters for a stretch-bend term.
@@ -87,10 +88,11 @@ public:
      * @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 k             the force constant for the stretch-bend
+     * @param k1            the force constant of the product of bond ab and angle a-b-c
+     * @param k2            the force constant of the product of bond bc and angle a-b-c
      */
-    void getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3,
-                                  double& lengthAB, double& lengthCB, double& angle, double& k ) const;
+    void getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3, double& lengthAB,
+                                  double& lengthCB, double& angle, double& k1, double& k2) const;
 
     /**
      * Set the force field parameters for a stretch-bend term.
@@ -102,10 +104,11 @@ public:
      * @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 k             the force constant for the stretch-bend
+     * @param k1            the force constant of the product of bond ab and angle a-b-c
+     * @param k2            the force constant of the product of bond bc and angle a-b-c (optional, default is the same as k1)
      */
     void setStretchBendParameters(int index, int particle1, int particle2, int particle3, 
-                                  double lengthAB,  double lengthCB, double angle, double k );
+                                  double lengthAB,  double lengthCB, double angle, double k1, double k2);
     /**
      * Update the per-stretch-bend term parameters in a Context to match those stored in this Force object.  This method provides
      * an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
@@ -139,16 +142,15 @@ private:
 class AmoebaStretchBendForce::StretchBendInfo {
 public:
     int particle1, particle2, particle3;
-    double lengthAB, lengthCB, angle, k;
+    double lengthAB, lengthCB, angle, k1, k2;
     StretchBendInfo() {
         particle1 = particle2  = particle3 = -1;
-        lengthAB  = lengthCB   = angle     = k   = 0.0;
+        lengthAB  = lengthCB   = angle     = k1 = k2 = 0.0;
     }
     StretchBendInfo(int particle1, int particle2, int particle3, 
-                    double lengthAB,  double lengthCB, double angle, double k ) :
+                    double lengthAB,  double lengthCB, double angle, double k1, double k2 ) :
                     particle1(particle1), particle2(particle2), particle3(particle3), 
-                    lengthAB(lengthAB), lengthCB(lengthCB), angle(angle), k(k) {
-     
+                    lengthAB(lengthAB), lengthCB(lengthCB), angle(angle), k1(k1), k2(k2) {
     }
 };
 

plugins/amoeba/openmmapi/src/AmoebaStretchBendForce.cpp

@@ -40,31 +40,33 @@ AmoebaStretchBendForce::AmoebaStretchBendForce() {
 }
 
 int AmoebaStretchBendForce::addStretchBend(int particle1, int particle2, int particle3,
-                                           double lengthAB,  double lengthCB, double angle, double k) {
-    stretchBends.push_back(StretchBendInfo(particle1, particle2, particle3, lengthAB,  lengthCB, angle, k));
+                                           double lengthAB,  double lengthCB, double angle, double k1, double k2) {
+    stretchBends.push_back(StretchBendInfo(particle1, particle2, particle3, lengthAB,  lengthCB, angle, k1, k2));
     return stretchBends.size()-1;
 }
 
 void AmoebaStretchBendForce::getStretchBendParameters(int index, int& particle1, int& particle2, int& particle3,
-                                                      double& lengthAB, double& lengthCB, double& angle, double& k ) const {
+                                                      double& lengthAB, double& lengthCB, double& angle, double& k1, double& k2 ) const {
     particle1       = stretchBends[index].particle1;
     particle2       = stretchBends[index].particle2;
     particle3       = stretchBends[index].particle3;
     lengthAB        = stretchBends[index].lengthAB;
     lengthCB        = stretchBends[index].lengthCB;
     angle           = stretchBends[index].angle;
-    k               = stretchBends[index].k;
+    k1              = stretchBends[index].k1;
+    k2              = stretchBends[index].k2;
 }
 
-void AmoebaStretchBendForce::setStretchBendParameters(int index, int particle1, int particle2, int particle3, 
-                                                      double lengthAB,  double lengthCB, double angle, double k) {
+void AmoebaStretchBendForce::setStretchBendParameters(int index, int particle1, int particle2, int particle3,
+                                                      double lengthAB,  double lengthCB, double angle, double k1, double k2) {
     stretchBends[index].particle1  = particle1;
     stretchBends[index].particle2  = particle2;
     stretchBends[index].particle3  = particle3;
     stretchBends[index].lengthAB   = lengthAB;
     stretchBends[index].lengthCB   = lengthCB;
     stretchBends[index].angle      = angle;
-    stretchBends[index].k          = k;
+    stretchBends[index].k1         = k1;
+    stretchBends[index].k2         = k2;
 }
 
 ForceImpl* AmoebaStretchBendForce::createImpl() const {

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -464,8 +464,8 @@ public:
     }
     void getParticlesInGroup(int index, std::vector<int>& particles) {
         int particle1, particle2, particle3;
-        double lengthAB, lengthCB, angle, k;
-        force.getStretchBendParameters(index, particle1, particle2, particle3, lengthAB, lengthCB, angle, k);
+        double lengthAB, lengthCB, angle, k1, k2;
+        force.getStretchBendParameters(index, particle1, particle2, particle3, lengthAB, lengthCB, angle, k1, k2);
         particles.resize(3);
         particles[0] = particle1;
         particles[1] = particle2;
@@ -473,23 +473,25 @@ public:
     }
     bool areGroupsIdentical(int group1, int group2) {
         int particle1, particle2, particle3;
-        double lengthAB1, lengthAB2, lengthCB1, lengthCB2, angle1, angle2, k1, k2;
-        force.getStretchBendParameters(group1, particle1, particle2, particle3, lengthAB1, lengthCB1, angle1, k1);
-        force.getStretchBendParameters(group2, particle1, particle2, particle3, lengthAB2, lengthCB2, angle2, k2);
-        return (lengthAB1 == lengthAB2 && lengthCB1 == lengthCB2 && angle1 == angle2 && k1 == k2);
+        double lengthAB1, lengthAB2, lengthCB1, lengthCB2, angle1, angle2, k11, k12, k21, k22;
+        force.getStretchBendParameters(group1, particle1, particle2, particle3, lengthAB1, lengthCB1, angle1, k11, k12);
+        force.getStretchBendParameters(group2, particle1, particle2, particle3, lengthAB2, lengthCB2, angle2, k21, k22);
+        return (lengthAB1 == lengthAB2 && lengthCB1 == lengthCB2 && angle1 == angle2 && k11 == k21 && k12 == k22);
     }
 private:
     const AmoebaStretchBendForce& force;
 };
 
 CudaCalcAmoebaStretchBendForceKernel::CudaCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) :
-                   CalcAmoebaStretchBendForceKernel(name, platform), cu(cu), system(system), params(NULL) {
+                   CalcAmoebaStretchBendForceKernel(name, platform), cu(cu), system(system), params1(NULL), params2(NULL) {
 }
 
 CudaCalcAmoebaStretchBendForceKernel::~CudaCalcAmoebaStretchBendForceKernel() {
     cu.setAsCurrent();
-    if (params != NULL)
-        delete params;
+    if (params1 != NULL)
+        delete params1;
+    if (params2 != NULL)
+        delete params2;
 }
 
 void CudaCalcAmoebaStretchBendForceKernel::initialize(const System& system, const AmoebaStretchBendForce& force) {
@@ -501,16 +503,21 @@ void CudaCalcAmoebaStretchBendForceKernel::initialize(const System& system, cons
     if (numStretchBends == 0)
         return;
     vector<vector<int> > atoms(numStretchBends, vector<int>(3));
-    params = CudaArray::create<float4>(cu, numStretchBends, "stretchBendParams");
-    vector<float4> paramVector(numStretchBends);
+    params1 = CudaArray::create<float3>(cu, numStretchBends, "stretchBendParams");
+    params2 = CudaArray::create<float2>(cu, numStretchBends, "stretchBendForceConstants");
+    vector<float3> paramVector(numStretchBends);
+    vector<float2> paramVectorK(numStretchBends);
     for (int i = 0; i < numStretchBends; i++) {
-        double lengthAB, lengthCB, angle, k;
-        force.getStretchBendParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], lengthAB, lengthCB, angle, k);
-        paramVector[i] = make_float4((float) lengthAB, (float) lengthCB, (float) angle, (float) k);
+        double lengthAB, lengthCB, angle, k1, k2;
+        force.getStretchBendParameters(startIndex+i, atoms[i][0], atoms[i][1], atoms[i][2], lengthAB, lengthCB, angle, k1, k2);
+        paramVector[i] = make_float3((float) lengthAB, (float) lengthCB, (float) angle);
+        paramVectorK[i] = make_float2((float) k1, (float) k2);
     }
-    params->upload(paramVector);
+    params1->upload(paramVector);
+    params2->upload(paramVectorK);
     map<string, string> replacements;
-    replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params->getDevicePointer(), "float4");
+    replacements["PARAMS"] = cu.getBondedUtilities().addArgument(params1->getDevicePointer(), "float3");
+    replacements["FORCE_CONSTANTS"] = cu.getBondedUtilities().addArgument(params2->getDevicePointer(), "float2");
     replacements["RAD_TO_DEG"] = cu.doubleToString(180/M_PI);
     cu.getBondedUtilities().addInteraction(atoms, cu.replaceStrings(CudaAmoebaKernelSources::amoebaStretchBendForce, replacements), force.getForceGroup());
     cu.addForce(new ForceInfo(force));
@@ -532,14 +539,17 @@ void CudaCalcAmoebaStretchBendForceKernel::copyParametersToContext(ContextImpl&
     
     // Record the per-stretch-bend parameters.
     
-    vector<float4> paramVector(numStretchBends);
+    vector<float3> paramVector(numStretchBends);
+    vector<float2> paramVector1(numStretchBends);
     for (int i = 0; i < numStretchBends; i++) {
         int atom1, atom2, atom3;
-        double lengthAB, lengthCB, angle, k;
-        force.getStretchBendParameters(startIndex+i, atom1, atom2, atom3, lengthAB, lengthCB, angle, k);
-        paramVector[i] = make_float4((float) lengthAB, (float) lengthCB, (float) angle, (float) k);
+        double lengthAB, lengthCB, angle, k1, k2;
+        force.getStretchBendParameters(startIndex+i, atom1, atom2, atom3, lengthAB, lengthCB, angle, k1, k2);
+        paramVector[i] = make_float3((float) lengthAB, (float) lengthCB, (float) angle);
+        paramVector1[i] = make_float2((float) k1, (float) k2);
     }
-    params->upload(paramVector);
+    params1->upload(paramVector);
+    params2->upload(paramVector1);
     
     // Mark that the current reordering may be invalid.
     

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h

@@ -229,7 +229,8 @@ private:
     int numStretchBends;
     CudaContext& cu;
     const System& system;
-    CudaArray* params;
+    CudaArray* params1; // Equilibrium values
+    CudaArray* params2; // force constants
 };
 
 /**

plugins/amoeba/platforms/cuda/src/kernels/amoebaStretchBendForce.cu

@@ -25,7 +25,8 @@ real angle = ACOS(cosine);
 
 // find chain rule terms for the bond angle deviation
 
-float4 parameters = PARAMS[index];
+float3 parameters = PARAMS[index];
+float2 force_constants = FORCE_CONSTANTS[index];
 
 real dt = RAD_TO_DEG*(angle - parameters.z);
 real terma = rab*rp != 0 ? (-RAD_TO_DEG/(rab*rab*rp)) : (real) 0;
@@ -41,11 +42,15 @@ real ddtdzic = termc * (xcb*yp-ycb*xp);
 
 // find chain rule terms for the bond length deviations
 
-real dr = (parameters.x > 0 ? (rab - parameters.x) : (real) 0);
+real dr1 = (parameters.x > 0 ? (rab - parameters.x) : (real) 0);
 terma = (parameters.x > 0 ? RECIP(rab) : (real) 0);
 
-dr += (parameters.y > 0 ? (rcb - parameters.y) : (real) 0);
+real dr2 = (parameters.y > 0 ? (rcb - parameters.y) : (real) 0);
 termc = (parameters.y > 0 ? RECIP(rcb) : (real) 0);
+real frc1 = ((rp != 0) ? force_constants.x : (real) 0);
+real frc2 = ((rp != 0) ? force_constants.y : (real) 0);
+
+real drkk = dr1*frc1 + dr2*frc2;
 
 real ddrdxia = terma * xab;
 real ddrdyia = terma * yab;
@@ -57,9 +62,8 @@ real ddrdzic = termc * zcb;
 
 // get the energy and master chain rule terms for derivatives
 
-real term = ((rp != 0) ? parameters.w : (real) 0);
-energy += term*dt*dr;
+energy += dt*drkk;
 
-real3 force1 = make_real3(-term*(dt*ddrdxia+ddtdxia*dr), -term*(dt*ddrdyia+ddtdyia*dr), -term*(dt*ddrdzia+ddtdzia*dr));
-real3 force3 = make_real3(-term*(dt*ddrdxic+ddtdxic*dr), -term*(dt*ddrdyic+ddtdyic*dr), -term*(dt*ddrdzic+ddtdzic*dr));
+real3 force1 = make_real3(-frc1*dt*ddrdxia-ddtdxia*drkk, -frc1*dt*ddrdyia-ddtdyia*drkk, -frc1*dt*ddrdzia-ddtdzia*drkk);
+real3 force3 = make_real3(-frc2*dt*ddrdxic-ddtdxic*drkk, -frc2*dt*ddrdyic-ddtdyic*drkk, -frc2*dt*ddrdzic-ddtdzic*drkk);
 real3 force2 = make_real3(-force1.x-force3.x, -force1.y-force3.y, -force1.z-force3.z);

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.cpp

@@ -82,14 +82,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
@@ -128,9 +128,11 @@ static void computeAmoebaStretchBendForce(int bondIndex,  std::vector<Vec3>& pos
     double angle;
     if( cosine >= 1.0 ){
        angle = 0.0;
-    } else if( cosine <= -1.0 ){
+    }
+    else if( cosine <= -1.0 ){
        angle = PI_M;
-    } else {
+    }
+    else {
        angle = RADIAN*acos(cosine);
     }
  
@@ -146,12 +148,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;
  
     // ---------------------------------------------------------------------------------------
  
@@ -163,8 +166,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] );
     }
  
@@ -184,8 +187,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 ); 
@@ -274,7 +276,7 @@ void testOneStretchBend( FILE* log ) {
     //double kStretchBend     = 0.750491578E-01;
     double kStretchBend     = 1.0;
 
-    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend );
+    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend );
 
     system.addForce(amoebaStretchBendForce);
     Context context(system, integrator, Platform::getPlatformByName( "CUDA"));
@@ -290,7 +292,7 @@ void testOneStretchBend( FILE* log ) {
     
     // Try changing the stretch-bend parameters and make sure it's still correct.
     
-    amoebaStretchBendForce->setStretchBendParameters(0, 0, 1, 2, 1.1*abLength, 1.2*cbLength, 1.3*angleStretchBend, 1.4*kStretchBend);
+    amoebaStretchBendForce->setStretchBendParameters(0, 0, 1, 2, 1.1*abLength, 1.2*cbLength, 1.3*angleStretchBend, 1.4*kStretchBend, 1.4*kStretchBend);
     bool exceptionThrown = false;
     try {
         // This should throw an exception.

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

@@ -307,15 +307,16 @@ void ReferenceCalcAmoebaStretchBendForceKernel::initialize(const System& system,
     numStretchBends = force.getNumStretchBends();
     for ( int ii = 0; ii < numStretchBends; ii++) {
         int particle1Index, particle2Index, particle3Index;
-        double lengthAB, lengthCB, angle, k;
-        force.getStretchBendParameters(ii, particle1Index, particle2Index, particle3Index, lengthAB, lengthCB, angle, k);
+        double lengthAB, lengthCB, angle, k1, k2;
+        force.getStretchBendParameters(ii, particle1Index, particle2Index, particle3Index, lengthAB, lengthCB, angle, k1, k2);
         particle1.push_back( particle1Index ); 
         particle2.push_back( particle2Index ); 
         particle3.push_back( particle3Index ); 
         lengthABParameters.push_back( static_cast<RealOpenMM>(lengthAB) );
         lengthCBParameters.push_back( static_cast<RealOpenMM>(lengthCB) );
         angleParameters.push_back(    static_cast<RealOpenMM>(angle) );
-        kParameters.push_back(        static_cast<RealOpenMM>(k) );
+        k1Parameters.push_back(       static_cast<RealOpenMM>(k1) );
+        k2Parameters.push_back(       static_cast<RealOpenMM>(k2) );
     }
 }
 
@@ -324,7 +325,8 @@ double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context,
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferenceStretchBendForce amoebaReferenceStretchBendForce;
     RealOpenMM energy      = amoebaReferenceStretchBendForce.calculateForceAndEnergy( numStretchBends, posData, particle1, particle2, particle3,
-                                                                                      lengthABParameters, lengthCBParameters, angleParameters, kParameters, forceData );
+                                                                                      lengthABParameters, lengthCBParameters, angleParameters, k1Parameters,
+                                                                                      k2Parameters, forceData );
     return static_cast<double>(energy);
 }
 
@@ -336,14 +338,15 @@ void ReferenceCalcAmoebaStretchBendForceKernel::copyParametersToContext(ContextI
 
     for (int i = 0; i < numStretchBends; ++i) {
         int particle1Index, particle2Index, particle3Index;
-        double lengthAB, lengthCB, angle, k;
-        force.getStretchBendParameters(i, particle1Index, particle2Index, particle3Index, lengthAB, lengthCB, angle, k);
+        double lengthAB, lengthCB, angle, k1, k2;
+        force.getStretchBendParameters(i, particle1Index, particle2Index, particle3Index, lengthAB, lengthCB, angle, k1, k2);
         if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i])
             throw OpenMMException("updateParametersInContext: The set of particles in a stretch-bend has changed");
         lengthABParameters[i] = (RealOpenMM) lengthAB;
         lengthCBParameters[i] = (RealOpenMM) lengthCB;
         angleParameters[i] = (RealOpenMM) angle;
-        kParameters[i] = (RealOpenMM) k;
+        k1Parameters[i] = (RealOpenMM) k1;
+        k2Parameters[i] = (RealOpenMM) k2;
     }
 }
 

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

@@ -248,7 +248,8 @@ private:
     std::vector<RealOpenMM> lengthABParameters;
     std::vector<RealOpenMM> lengthCBParameters;
     std::vector<RealOpenMM> angleParameters;
-    std::vector<RealOpenMM> kParameters;
+    std::vector<RealOpenMM> k1Parameters;
+    std::vector<RealOpenMM> k2Parameters;
     const System& system;
 };
 

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

@@ -52,8 +52,8 @@ using OpenMM::RealVec;
 RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn( const RealVec& positionAtomA, const RealVec& positionAtomB,
                                                                      const RealVec& positionAtomC,
                                                                      RealOpenMM lengthAB,      RealOpenMM lengthCB,
-                                                                     RealOpenMM idealAngle,    RealOpenMM kParameter,
-                                                                     RealVec* forces ) const {
+                                                                     RealOpenMM idealAngle,    RealOpenMM k1Parameter,
+                                                                     RealOpenMM k2Parameter, RealVec* forces ) const {
 
    // ---------------------------------------------------------------------------------------
 
@@ -112,7 +112,9 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn( const RealV
        deltaR[CBxP][ii] *= termC;
     }
  
-    RealOpenMM dr    = rAB - lengthAB + rCB - lengthCB;
+    RealOpenMM dr1   = rAB - lengthAB;
+    RealOpenMM dr2   = rCB - lengthCB;
+    RealOpenMM drkk  = dr1*k1Parameter + dr2*k2Parameter;
     termA            = one/rAB;
     termC            = one/rCB;
  
@@ -129,8 +131,8 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn( const RealV
     }
     RealOpenMM dt = angle - idealAngle*RADIAN;
     for( int jj = 0; jj < 3; jj++ ){
-       subForce[A][jj] = kParameter*(dt*termA*deltaR[AB][jj] + dr*deltaR[ABxP][jj] );
-       subForce[C][jj] = kParameter*(dt*termC*deltaR[CB][jj] + dr*deltaR[CBxP][jj] );
+       subForce[A][jj] = k1Parameter*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj];
+       subForce[C][jj] = k2Parameter*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj];
        subForce[B][jj] = -( subForce[A][jj] + subForce[C][jj] );
     }
  
@@ -144,7 +146,7 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn( const RealV
  
     // ---------------------------------------------------------------------------------------
  
-    return (kParameter*dt*dr);
+    return dt*drkk;
 }
 
 RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( int numStretchBends, vector<RealVec>& posData,
@@ -154,7 +156,8 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( int numStre
                                                                        const std::vector<RealOpenMM>& lengthABParameters,
                                                                        const std::vector<RealOpenMM>& lengthCBParameters,
                                                                        const std::vector<RealOpenMM>&  angle,
-                                                                       const std::vector<RealOpenMM>&  kQuadratic,
+                                                                       const std::vector<RealOpenMM>&  k1Quadratic,
+                                                                       const std::vector<RealOpenMM>&  k2Quadratic,
                                                                        vector<RealVec>& forceData) const {
     RealOpenMM energy      = 0.0; 
     for (unsigned int ii = 0; ii < static_cast<unsigned int>(numStretchBends); ii++) {
@@ -164,10 +167,11 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( int numStre
         RealOpenMM abLength     = lengthABParameters[ii];
         RealOpenMM cbLength     = lengthCBParameters[ii];
         RealOpenMM idealAngle   = angle[ii];
-        RealOpenMM angleK       = kQuadratic[ii];
+        RealOpenMM angleK1      = k1Quadratic[ii];
+        RealOpenMM angleK2      = k2Quadratic[ii];
         RealVec forces[3];
         energy                 += calculateStretchBendIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index],
-                                                           abLength, cbLength, idealAngle, angleK, forces );
+                                                           abLength, cbLength, idealAngle, angleK1, angleK2, forces );
         // accumulate forces
     
         for( int jj = 0; jj < 3; jj++ ){

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

@@ -77,7 +77,8 @@ public:
                                         const std::vector<RealOpenMM>& lengthABParameters,
                                         const std::vector<RealOpenMM>& lengthCBParameters,
                                         const std::vector<RealOpenMM>&  angle,
-                                        const std::vector<RealOpenMM>&  kQuadratic,
+                                        const std::vector<RealOpenMM>&  k1Quadratic,
+                                        const std::vector<RealOpenMM>&  k2Quadratic,
                                         std::vector<OpenMM::RealVec>& forceData ) const;
 
 
@@ -93,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
@@ -103,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 ); 
@@ -274,7 +275,7 @@ void testOneStretchBend( FILE* log ) {
     //double kStretchBend     = 0.750491578E-01;
     double kStretchBend     = 1.0;
 
-    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend );
+    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend );
 
     system.addForce(amoebaStretchBendForce);
     ASSERT(!amoebaStretchBendForce->usesPeriodicBoundaryConditions());
@@ -292,7 +293,7 @@ void testOneStretchBend( FILE* log ) {
     
     // Try changing the stretch-bend parameters and make sure it's still correct.
     
-    amoebaStretchBendForce->setStretchBendParameters(0, 0, 1, 2, 1.1*abLength, 1.2*cbLength, 1.3*angleStretchBend, 1.4*kStretchBend);
+    amoebaStretchBendForce->setStretchBendParameters(0, 0, 1, 2, 1.1*abLength, 1.2*cbLength, 1.3*angleStretchBend, 1.4*kStretchBend, 1.4*kStretchBend);
     bool exceptionThrown = false;
     try {
         // This should throw an exception.

plugins/amoeba/serialization/src/AmoebaStretchBendForceProxy.cpp

@@ -47,9 +47,9 @@ void AmoebaStretchBendForceProxy::serialize(const void* object, SerializationNod
     SerializationNode& bonds = node.createChildNode("StretchBendAngles");
     for (unsigned int ii = 0; ii < static_cast<unsigned int>(force.getNumStretchBends()); ii++) {
         int particle1, particle2, particle3;
-        double distanceAB, distanceCB, angle, k;
-        force.getStretchBendParameters(ii, particle1, particle2, particle3, distanceAB, distanceCB, angle, k);
-        bonds.createChildNode("StretchBendAngle").setIntProperty("p1", particle1).setIntProperty("p2", particle2).setIntProperty("p3", particle3).setDoubleProperty("dAB", distanceAB).setDoubleProperty("dCB", distanceCB).setDoubleProperty("angle", angle).setDoubleProperty("k", k);
+        double distanceAB, distanceCB, angle, k1, k2;
+        force.getStretchBendParameters(ii, particle1, particle2, particle3, distanceAB, distanceCB, angle, k1, k2);
+        bonds.createChildNode("StretchBendAngle").setIntProperty("p1", particle1).setIntProperty("p2", particle2).setIntProperty("p3", particle3).setDoubleProperty("dAB", distanceAB).setDoubleProperty("dCB", distanceCB).setDoubleProperty("angle", angle).setDoubleProperty("k1", k1).setDoubleProperty("k2", k2);
     }
 
 }
@@ -62,7 +62,7 @@ void* AmoebaStretchBendForceProxy::deserialize(const SerializationNode& node) co
         const SerializationNode& bonds = node.getChildNode("StretchBendAngles");
         for ( unsigned int ii = 0; ii < (int) bonds.getChildren().size(); ii++) {
             const SerializationNode& bond = bonds.getChildren()[ii];
-            force->addStretchBend(bond.getIntProperty("p1"), bond.getIntProperty("p2"), bond.getIntProperty("p3"),  bond.getDoubleProperty("dAB"), bond.getDoubleProperty("dCB"), bond.getDoubleProperty("angle"), bond.getDoubleProperty("k"));
+            force->addStretchBend(bond.getIntProperty("p1"), bond.getIntProperty("p2"), bond.getIntProperty("p3"),  bond.getDoubleProperty("dAB"), bond.getDoubleProperty("dCB"), bond.getDoubleProperty("angle"), bond.getDoubleProperty("k1"), bond.getDoubleProperty("k2"));
 
         }
     }

plugins/amoeba/serialization/tests/TestSerializeAmoebaStretchBendForce.cpp

@@ -45,9 +45,9 @@ void testSerialization() {
     // Create a Force.
 
     AmoebaStretchBendForce force1;
-    force1.addStretchBend(0, 1, 3, 1.0, 1.2, 150.1, 83.2);
-    force1.addStretchBend(2, 4, 4, 1.1, 2.2, 180.1, 89.2);
-    force1.addStretchBend(5, 0, 1, 3.1, 8.2, 140.1, 98.2);
+    force1.addStretchBend(0, 1, 3, 1.0, 1.2, 150.1, 83.2, 100.);
+    force1.addStretchBend(2, 4, 4, 1.1, 2.2, 180.1, 89.2, 100.);
+    force1.addStretchBend(5, 0, 1, 3.1, 8.2, 140.1, 98.2, 100.);
 
     // Serialize and then deserialize it.
 
@@ -64,10 +64,10 @@ void testSerialization() {
         double dAB1, dAB2;
         double dCB1, dCB2;
         double angle1, angle2;
-        double k1, k2;
+        double k11, k12, k21, k22;
 
-        force1.getStretchBendParameters(ii, p11, p12, p13, dAB1, dCB1, angle1, k1);
-        force2.getStretchBendParameters(ii, p21, p22, p23, dAB2, dCB2, angle2, k2);
+        force1.getStretchBendParameters(ii, p11, p12, p13, dAB1, dCB1, angle1, k11, k12);
+        force2.getStretchBendParameters(ii, p21, p22, p23, dAB2, dCB2, angle2, k21, k22);
 
         ASSERT_EQUAL(p11, p21);
         ASSERT_EQUAL(p12, p22);
@@ -75,7 +75,8 @@ void testSerialization() {
         ASSERT_EQUAL(dAB1, dAB2);
         ASSERT_EQUAL(dCB1, dCB2);
         ASSERT_EQUAL(angle1, angle2);
-        ASSERT_EQUAL(k1, k2);
+        ASSERT_EQUAL(k11, k21);
+        ASSERT_EQUAL(k12, k22);
     }
 }
 

wrappers/python/simtk/openmm/app/forcefield.py

@@ -3086,7 +3086,7 @@ class AmoebaStretchBendGenerator:
                        raise ValueError(outputString)
 
                     else:
-                        force.addStretchBend(angle[0], angle[1], angle[2], bondAB, bondCB, angleDict['idealAngle']/radian, self.k1[i])
+                        force.addStretchBend(angle[0], angle[1], angle[2], bondAB, bondCB, angleDict['idealAngle']/radian, self.k1[i], self.k2[i])
 
                     break
 

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -277,7 +277,7 @@ UNITS = {
 ("AmoebaPiTorsionForce",                  "getPiTorsionParameters")                        :  ( None, (None, None, None, None, None,  None, 'unit.kilojoule_per_mole')),
 
 ("AmoebaStretchBendForce",                "getNumStretchBends")                            :  ( None, ()),
-("AmoebaStretchBendForce",                "getStretchBendParameters")                      :  ( None, (None, None, None, 'unit.nanometer', 'unit.nanometer', 'unit.radian', 'unit.kilojoule_per_mole/unit.nanometer/unit.degree')),
+("AmoebaStretchBendForce",                "getStretchBendParameters")                      :  ( None, (None, None, None, 'unit.nanometer', 'unit.nanometer', 'unit.radian', 'unit.kilojoule_per_mole/unit.nanometer/unit.degree', 'unit.kilojoule_per_mole/unit.nanometer/unit.degree')),
 
 ("AmoebaTorsionTorsionForce",             "getNumTorsionTorsions")                         :  ( None, ()),
 ("AmoebaTorsionTorsionForce",             "getNumTorsionTorsionGrids")                     :  ( None, ()),