Bug fixes to CustomGBForce

platforms/opencl/src/OpenCLKernels.cpp

@@ -2144,8 +2144,8 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
         if (needParameterGradient) {
             chainSource << "float4 grad1_0_1 = dV0dR1*delta*invR;\n";
             chainSource << "float4 grad1_0_2 = dV0dR2*delta*invR;\n";
-            chainSource << "float4 grad2_0_1 = -grad1_0_1;\n";
-            chainSource << "float4 grad2_0_2 = -grad1_0_2;\n";
+            chainSource << "float4 grad2_0_1 = grad1_0_1;\n";
+            chainSource << "float4 grad2_0_2 = grad1_0_2;\n";
             chainSource << "tempForce1 -= grad1_0_1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "1") << ";\n";
             chainSource << "tempForce1 -= grad1_0_2*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "2") << ";\n";
             chainSource << "tempForce2 -= grad2_0_1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "1") << ";\n";
@@ -2158,6 +2158,18 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
         variables = globalVariables;
         map<string, string> rename1;
         map<string, string> rename2;
+        variables["x1"] = "posq1.x";
+        variables["y1"] = "posq1.y";
+        variables["z1"] = "posq1.z";
+        variables["x2"] = "posq2.x";
+        variables["y2"] = "posq2.y";
+        variables["z2"] = "posq2.z";
+        rename1["x"] = "x1";
+        rename1["y"] = "y1";
+        rename1["z"] = "z1";
+        rename2["x"] = "x2";
+        rename2["y"] = "y2";
+        rename2["z"] = "z2";
         for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
             const string& name = force.getPerParticleParameterName(i);
             variables[name+"1"] = prefix+"params"+params->getParameterSuffix(i, "1");

platforms/opencl/tests/TestOpenCLCustomGBForce.cpp

@@ -210,9 +210,9 @@ void testPositionDependence() {
     VerletIntegrator integrator(0.01);
     CustomGBForce* force = new CustomGBForce();
     force->addComputedValue("a", "r", CustomGBForce::ParticlePair);
-    force->addComputedValue("b", "a+y", CustomGBForce::SingleParticle);
+    force->addComputedValue("b", "a+x*y", CustomGBForce::SingleParticle);
     force->addEnergyTerm("b*z", CustomGBForce::SingleParticle);
-    force->addEnergyTerm("b1+b2", CustomGBForce::ParticlePair); // = 2*r+y1+y2
+    force->addEnergyTerm("b1+b2", CustomGBForce::ParticlePair); // = 2*r+x1*y1+x2*y2
     force->addParticle(vector<double>());
     force->addParticle(vector<double>());
     system.addForce(force);
@@ -230,15 +230,15 @@ void testPositionDependence() {
         const vector<Vec3>& forces = state.getForces();
         Vec3 delta = positions[0]-positions[1];
         double r = sqrt(delta.dot(delta));
-        double energy = 2*r+positions[0][1]+positions[1][1];
+        double energy = 2*r+positions[0][0]*positions[0][1]+positions[1][0]*positions[1][1];
         for (int j = 0; j < 2; j++)
-            energy += positions[j][2]*(r+positions[j][1]);
-        Vec3 force1(-(1+positions[0][2])*delta[0]/r-(1+positions[1][2])*delta[0]/r,
-                    -(1+positions[0][2])*delta[1]/r-positions[0][2]-(1+positions[1][2])*delta[1]/r-1,
-                    -(1+positions[0][2])*delta[2]/r-(r+positions[0][1])-(1+positions[1][2])*delta[2]/r);
-        Vec3 force2((1+positions[0][2])*delta[0]/r+(1+positions[1][2])*delta[0]/r,
-                    (1+positions[0][2])*delta[1]/r+(1+positions[1][2])*delta[1]/r-positions[1][2]-1,
-                    (1+positions[0][2])*delta[2]/r+(1+positions[1][2])*delta[2]/r-(r+positions[1][1]));
+            energy += positions[j][2]*(r+positions[j][0]*positions[j][1]);
+        Vec3 force1(-(1+positions[0][2])*delta[0]/r-(1+positions[0][2])*positions[0][1]-(1+positions[1][2])*delta[0]/r,
+                    -(1+positions[0][2])*delta[1]/r-(1+positions[0][2])*positions[0][0]-(1+positions[1][2])*delta[1]/r,
+                    -(1+positions[0][2])*delta[2]/r-(r+positions[0][0]*positions[0][1])-(1+positions[1][2])*delta[2]/r);
+        Vec3 force2((1+positions[0][2])*delta[0]/r+(1+positions[1][2])*delta[0]/r-(1+positions[1][2])*positions[1][1],
+                    (1+positions[0][2])*delta[1]/r+(1+positions[1][2])*delta[1]/r-(1+positions[1][2])*positions[1][0],
+                    (1+positions[0][2])*delta[2]/r+(1+positions[1][2])*delta[2]/r-(r+positions[1][0]*positions[1][1]));
         ASSERT_EQUAL_VEC(force1, forces[0], 1e-4);
         ASSERT_EQUAL_VEC(force2, forces[1], 1e-4);
         ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);

platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp

@@ -397,6 +397,9 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, RealO
     variables[valueNames[0]] = values[0][atom1];
     for (int i = 1; i < (int) valueNames.size(); i++) {
         variables[valueNames[i]] = values[i][atom1];
+        variables["x"] = atomCoordinates[atom1][0];
+        variables["y"] = atomCoordinates[atom1][1];
+        variables["z"] = atomCoordinates[atom1][2];
         for (int j = 0; j < i; j++) {
             RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[i][j].evaluate(variables);
             for (int k = 0; k < 3; k++) {

platforms/reference/tests/TestReferenceCustomGBForce.cpp

@@ -212,9 +212,9 @@ void testPositionDependence() {
     VerletIntegrator integrator(0.01);
     CustomGBForce* force = new CustomGBForce();
     force->addComputedValue("a", "r", CustomGBForce::ParticlePair);
-    force->addComputedValue("b", "a+y", CustomGBForce::SingleParticle);
+    force->addComputedValue("b", "a+x*y", CustomGBForce::SingleParticle);
     force->addEnergyTerm("b*z", CustomGBForce::SingleParticle);
-    force->addEnergyTerm("b1+b2", CustomGBForce::ParticlePair); // = 2*r+y1+y2
+    force->addEnergyTerm("b1+b2", CustomGBForce::ParticlePair); // = 2*r+x1*y1+x2*y2
     force->addParticle(vector<double>());
     force->addParticle(vector<double>());
     system.addForce(force);
@@ -232,15 +232,15 @@ void testPositionDependence() {
         const vector<Vec3>& forces = state.getForces();
         Vec3 delta = positions[0]-positions[1];
         double r = sqrt(delta.dot(delta));
-        double energy = 2*r+positions[0][1]+positions[1][1];
+        double energy = 2*r+positions[0][0]*positions[0][1]+positions[1][0]*positions[1][1];
         for (int j = 0; j < 2; j++)
-            energy += positions[j][2]*(r+positions[j][1]);
-        Vec3 force1(-(1+positions[0][2])*delta[0]/r-(1+positions[1][2])*delta[0]/r,
-                    -(1+positions[0][2])*delta[1]/r-positions[0][2]-(1+positions[1][2])*delta[1]/r-1,
-                    -(1+positions[0][2])*delta[2]/r-(r+positions[0][1])-(1+positions[1][2])*delta[2]/r);
-        Vec3 force2((1+positions[0][2])*delta[0]/r+(1+positions[1][2])*delta[0]/r,
-                    (1+positions[0][2])*delta[1]/r+(1+positions[1][2])*delta[1]/r-positions[1][2]-1,
-                    (1+positions[0][2])*delta[2]/r+(1+positions[1][2])*delta[2]/r-(r+positions[1][1]));
+            energy += positions[j][2]*(r+positions[j][0]*positions[j][1]);
+        Vec3 force1(-(1+positions[0][2])*delta[0]/r-(1+positions[0][2])*positions[0][1]-(1+positions[1][2])*delta[0]/r,
+                    -(1+positions[0][2])*delta[1]/r-(1+positions[0][2])*positions[0][0]-(1+positions[1][2])*delta[1]/r,
+                    -(1+positions[0][2])*delta[2]/r-(r+positions[0][0]*positions[0][1])-(1+positions[1][2])*delta[2]/r);
+        Vec3 force2((1+positions[0][2])*delta[0]/r+(1+positions[1][2])*delta[0]/r-(1+positions[1][2])*positions[1][1],
+                    (1+positions[0][2])*delta[1]/r+(1+positions[1][2])*delta[1]/r-(1+positions[1][2])*positions[1][0],
+                    (1+positions[0][2])*delta[2]/r+(1+positions[1][2])*delta[2]/r-(r+positions[1][0]*positions[1][1]));
         ASSERT_EQUAL_VEC(force1, forces[0], 1e-4);
         ASSERT_EQUAL_VEC(force2, forces[1], 1e-4);
         ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);