First stab at adjusting AmoebaStretchBendForce to take 2 force constants instead

of just 1.  Backwards compatibility is provided by making the 2nd force constant
default to -1, which is reinterpreted as "copy the first force constant".

Updates both the reference and CUDA kernels.

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=-1.0);
 
     /**
      * 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=-1.0 );
     /**
      * 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,37 @@ 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) {
+    if (k2 == -1.0) k2 = k1
+    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;
+    if (k2 == -1.0)
+        stretchBends[index].k2     = k1;
+    else
+        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,10 +473,10 @@ 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;
@@ -488,8 +488,10 @@ CudaCalcAmoebaStretchBendForceKernel::CudaCalcAmoebaStretchBendForceKernel(std::
 
 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));

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/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;
+        k1Parameters[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

@@ -53,7 +53,7 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn( const RealV
                                                                      const RealVec& positionAtomC,
                                                                      RealOpenMM lengthAB,      RealOpenMM lengthCB,
                                                                      RealOpenMM idealAngle,    RealOpenMM kParameter,
-                                                                     RealVec* forces ) const {
+                                                                     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,

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

@@ -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);
     }
 }