Merge pull request #282 from peastman/master

Fixed a bug in CustomGBForce

platforms/cuda/include/CudaKernels.h

@@ -738,7 +738,7 @@ private:
 class CudaCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
 public:
     CudaCalcCustomGBForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : CalcCustomGBForceKernel(name, platform),
-            hasInitializedKernels(false), cu(cu), params(NULL), computedValues(NULL), energyDerivs(NULL), longEnergyDerivs(NULL), globals(NULL),
+            hasInitializedKernels(false), cu(cu), params(NULL), computedValues(NULL), energyDerivs(NULL), energyDerivChain(NULL), longEnergyDerivs(NULL), globals(NULL),
             valueBuffers(NULL), tabulatedFunctionParams(NULL), system(system) {
     }
     ~CudaCalcCustomGBForceKernel();
@@ -772,6 +772,7 @@ private:
     CudaParameterSet* params;
     CudaParameterSet* computedValues;
     CudaParameterSet* energyDerivs;
+    CudaParameterSet* energyDerivChain;
     CudaArray* longEnergyDerivs;
     CudaArray* globals;
     CudaArray* valueBuffers;

platforms/cuda/src/CudaKernels.cpp

@@ -2602,6 +2602,8 @@ CudaCalcCustomGBForceKernel::~CudaCalcCustomGBForceKernel() {
         delete computedValues;
     if (energyDerivs != NULL)
         delete energyDerivs;
+    if (energyDerivChain != NULL)
+        delete energyDerivChain;
     if (longEnergyDerivs != NULL)
         delete longEnergyDerivs;
     if (globals != NULL)
@@ -2745,6 +2747,7 @@ void CudaCalcCustomGBForceKernel::initialize(const System& system, const CustomG
     }
     longEnergyDerivs = CudaArray::create<long long>(cu, force.getNumComputedValues()*cu.getPaddedNumAtoms(), "customGBLongEnergyDerivatives");
     energyDerivs = new CudaParameterSet(cu, force.getNumComputedValues(), cu.getPaddedNumAtoms(), "customGBEnergyDerivatives", true);
+    energyDerivChain = new CudaParameterSet(cu, force.getNumComputedValues(), cu.getPaddedNumAtoms(), "customGBEnergyDerivativeChain", true);
  
     // Create the kernels.
 
@@ -3011,6 +3014,11 @@ void CudaCalcCustomGBForceKernel::initialize(const System& system, const CustomG
             extraArgs << ", " << buffer.getType() << "* __restrict__ derivBuffers" << index;
             compute << buffer.getType() << " deriv" << index << " = derivBuffers" << index << "[index];\n";
         }
+        for (int i = 0; i < (int) energyDerivChain->getBuffers().size(); i++) {
+            CudaNonbondedUtilities::ParameterInfo& buffer = energyDerivChain->getBuffers()[i];
+            string index = cu.intToString(i+1);
+            extraArgs << ", " << buffer.getType() << "* __restrict__ derivChain" << index;
+        }
         extraArgs << ", const long long* __restrict__ derivBuffersIn";
         for (int i = 0; i < energyDerivs->getNumParameters(); ++i)
             load << "derivBuffers" << energyDerivs->getParameterSuffix(i, "[index]") <<
@@ -3056,6 +3064,10 @@ void CudaCalcCustomGBForceKernel::initialize(const System& system, const CustomG
         
         // Record values.
         
+        for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
+            string index = cu.intToString(i+1);
+            compute << "derivBuffers" << index << "[index] = deriv" << index << ";\n";
+        }
         compute << "forceBuffers[index] += (long long) (force.x*0x100000000);\n";
         compute << "forceBuffers[index+PADDED_NUM_ATOMS] += (long long) (force.y*0x100000000);\n";
         compute << "forceBuffers[index+PADDED_NUM_ATOMS*2] += (long long) (force.z*0x100000000);\n";
@@ -3068,7 +3080,7 @@ void CudaCalcCustomGBForceKernel::initialize(const System& system, const CustomG
         }
         for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
             string index = cu.intToString(i+1);
-            compute << "derivBuffers" << index << "[index] = deriv" << index << ";\n";
+            compute << "derivChain" << index << "[index] = deriv" << index << ";\n";
         }
         map<string, string> replacements;
         replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
@@ -3206,9 +3218,9 @@ void CudaCalcCustomGBForceKernel::initialize(const System& system, const CustomG
             if (chainStr.find(paramName+"1") != chainStr.npos || chainStr.find(paramName+"2") != chainStr.npos)
                 parameters.push_back(CudaNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
         }
-        for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
+        for (int i = 0; i < (int) energyDerivChain->getBuffers().size(); i++) {
             if (needChainForValue[i]) { 
-                CudaNonbondedUtilities::ParameterInfo& buffer = energyDerivs->getBuffers()[i];
+                CudaNonbondedUtilities::ParameterInfo& buffer = energyDerivChain->getBuffers()[i];
                 string paramName = prefix+"dEdV"+cu.intToString(i+1);
                 parameters.push_back(CudaNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
             }
@@ -3354,6 +3366,8 @@ double CudaCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeFo
             perParticleEnergyArgs.push_back(&computedValues->getBuffers()[i].getMemory());
         for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++)
             perParticleEnergyArgs.push_back(&energyDerivs->getBuffers()[i].getMemory());
+        for (int i = 0; i < (int) energyDerivChain->getBuffers().size(); i++)
+            perParticleEnergyArgs.push_back(&energyDerivChain->getBuffers()[i].getMemory());
         perParticleEnergyArgs.push_back(&longEnergyDerivs->getDevicePointer());
         if (tabulatedFunctionParams != NULL) {
             for (int i = 0; i < (int) tabulatedFunctions.size(); i++)

platforms/cuda/tests/TestCudaCustomGBForce.cpp

@@ -248,7 +248,7 @@ void testMembrane() {
     for (int i = 0; i < (int) forces.size(); ++i)
         norm += forces[i].dot(forces[i]);
     norm = std::sqrt(norm);
-    const double stepSize = 1e-3;
+    const double stepSize = 1e-2;
     double step = 0.5*stepSize/norm;
     vector<Vec3> positions2(numParticles), positions3(numParticles);
     for (int i = 0; i < (int) positions.size(); ++i) {

platforms/opencl/include/OpenCLKernels.h

@@ -741,7 +741,7 @@ private:
 class OpenCLCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
 public:
     OpenCLCalcCustomGBForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomGBForceKernel(name, platform),
-            hasInitializedKernels(false), cl(cl), params(NULL), computedValues(NULL), energyDerivs(NULL), longEnergyDerivs(NULL), globals(NULL),
+            hasInitializedKernels(false), cl(cl), params(NULL), computedValues(NULL), energyDerivs(NULL), energyDerivChain(NULL), longEnergyDerivs(NULL), globals(NULL),
             valueBuffers(NULL), longValueBuffers(NULL), tabulatedFunctionParams(NULL), system(system) {
     }
     ~OpenCLCalcCustomGBForceKernel();
@@ -775,6 +775,7 @@ private:
     OpenCLParameterSet* params;
     OpenCLParameterSet* computedValues;
     OpenCLParameterSet* energyDerivs;
+    OpenCLParameterSet* energyDerivChain;
     OpenCLArray* longEnergyDerivs;
     OpenCLArray* globals;
     OpenCLArray* valueBuffers;

platforms/opencl/src/OpenCLKernels.cpp

@@ -2654,6 +2654,8 @@ OpenCLCalcCustomGBForceKernel::~OpenCLCalcCustomGBForceKernel() {
         delete computedValues;
     if (energyDerivs != NULL)
         delete energyDerivs;
+    if (energyDerivChain != NULL)
+        delete energyDerivChain;
     if (longEnergyDerivs != NULL)
         delete longEnergyDerivs;
     if (globals != NULL)
@@ -2804,6 +2806,7 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
     }
     else
         energyDerivs = new OpenCLParameterSet(cl, force.getNumComputedValues(), cl.getPaddedNumAtoms()*cl.getNonbondedUtilities().getNumForceBuffers(), "customGBEnergyDerivatives", true);
+    energyDerivChain = new OpenCLParameterSet(cl, force.getNumComputedValues(), cl.getPaddedNumAtoms(), "customGBEnergyDerivativeChain", true);
 
     // Create the kernels.
 
@@ -3094,6 +3097,11 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
             extraArgs << ", __global " << buffer.getType() << "* restrict derivBuffers" << index;
             compute << buffer.getType() << " deriv" << index << " = derivBuffers" << index << "[index];\n";
         }
+        for (int i = 0; i < (int) energyDerivChain->getBuffers().size(); i++) {
+            const OpenCLNonbondedUtilities::ParameterInfo& buffer = energyDerivChain->getBuffers()[i];
+            string index = cl.intToString(i+1);
+            extraArgs << ", __global " << buffer.getType() << "* restrict derivChain" << index;
+        }
         if (useLong) {
             extraArgs << ", __global const long* restrict derivBuffersIn";
             for (int i = 0; i < energyDerivs->getNumParameters(); ++i)
@@ -3145,6 +3153,10 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
         
         // Record values.
         
+        for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
+            string index = cl.intToString(i+1);
+            compute << "derivBuffers" << index << "[index] = deriv" << index << ";\n";
+        }
         compute << "forceBuffers[index] = forceBuffers[index]+force;\n";
         for (int i = 1; i < force.getNumComputedValues(); i++) {
             compute << "real totalDeriv"<<i<<" = dV"<<i<<"dV0";
@@ -3155,7 +3167,7 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
         }
         for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
             string index = cl.intToString(i+1);
-            compute << "derivBuffers" << index << "[index] = deriv" << index << ";\n";
+            compute << "derivChain" << index << "[index] = deriv" << index << ";\n";
         }
         map<string, string> replacements;
         replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
@@ -3292,9 +3304,9 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
             if (chainStr.find(paramName+"1") != chainStr.npos || chainStr.find(paramName+"2") != chainStr.npos)
                 parameters.push_back(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
         }
-        for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
+        for (int i = 0; i < (int) energyDerivChain->getBuffers().size(); i++) {
             if (needChainForValue[i]) { 
-                const OpenCLNonbondedUtilities::ParameterInfo& buffer = energyDerivs->getBuffers()[i];
+                const OpenCLNonbondedUtilities::ParameterInfo& buffer = energyDerivChain->getBuffers()[i];
                 string paramName = prefix+"dEdV"+cl.intToString(i+1);
                 parameters.push_back(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
             }
@@ -3487,6 +3499,8 @@ double OpenCLCalcCustomGBForceKernel::execute(ContextImpl& context, bool include
             perParticleEnergyKernel.setArg<cl::Memory>(index++, computedValues->getBuffers()[i].getMemory());
         for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++)
             perParticleEnergyKernel.setArg<cl::Memory>(index++, energyDerivs->getBuffers()[i].getMemory());
+        for (int i = 0; i < (int) energyDerivChain->getBuffers().size(); i++)
+            perParticleEnergyKernel.setArg<cl::Memory>(index++, energyDerivChain->getBuffers()[i].getMemory());
         if (useLong)
             perParticleEnergyKernel.setArg<cl::Memory>(index++, longEnergyDerivs->getDeviceBuffer());
         if (tabulatedFunctionParams != NULL) {

platforms/opencl/tests/TestOpenCLCustomGBForce.cpp

@@ -248,7 +248,7 @@ void testMembrane() {
     for (int i = 0; i < (int) forces.size(); ++i)
         norm += forces[i].dot(forces[i]);
     norm = std::sqrt(norm);
-    const double stepSize = 1e-3;
+    const double stepSize = 1e-2;
     double step = 0.5*stepSize/norm;
     vector<Vec3> positions2(numParticles), positions3(numParticles);
     for (int i = 0; i < (int) positions.size(); ++i) {

platforms/reference/tests/TestReferenceCustomGBForce.cpp

@@ -251,7 +251,7 @@ void testMembrane() {
     for (int i = 0; i < (int) forces.size(); ++i)
         norm += forces[i].dot(forces[i]);
     norm = std::sqrt(norm);
-    const double stepSize = 1e-3;
+    const double stepSize = 1e-2;
     double step = 0.5*stepSize/norm;
     vector<Vec3> positions2(numParticles), positions3(numParticles);
     for (int i = 0; i < (int) positions.size(); ++i) {