Merge pull request #1850 from peastman/customgb

Fixed bug in CustomGBForce

platforms/cpu/include/CpuKernels.h

@@ -373,7 +373,7 @@ private:
 class CpuCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
 public:
     CpuCalcCustomGBForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) :
-            CalcCustomGBForceKernel(name, platform), data(data) {
+            CalcCustomGBForceKernel(name, platform), data(data), ixn(NULL), neighborList(NULL) {
     }
     ~CpuCalcCustomGBForceKernel();
     /**
@@ -406,6 +406,7 @@ private:
     double **particleParamArray;
     double nonbondedCutoff;
     CpuCustomGBForce* ixn;
+    CpuNeighborList* neighborList;
     std::vector<std::set<int> > exclusions;
     std::vector<std::string> particleParameterNames, globalParameterNames, energyParamDerivNames, valueNames;
     std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;

platforms/cpu/src/CpuKernels.cpp

@@ -1018,6 +1018,8 @@ CpuCalcCustomGBForceKernel::~CpuCalcCustomGBForceKernel() {
     }
     if (ixn != NULL)
         delete ixn;
+    if (neighborList != NULL)
+        delete neighborList;
 }
 
 void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
@@ -1064,7 +1066,7 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
     nonbondedMethod = CalcCustomGBForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = force.getCutoffDistance();
     if (nonbondedMethod != NoCutoff)
-        data.requestNeighborList(nonbondedCutoff, 0.25*nonbondedCutoff, force.getNumExclusions() > 0, exclusions);
+        neighborList = new CpuNeighborList(4);
 
     // Create custom functions for the tabulated functions.
 
@@ -1171,7 +1173,8 @@ double CpuCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeFor
         ixn->setPeriodic(extractBoxSize(context));
     if (nonbondedMethod != NoCutoff) {
         vector<set<int> > noExclusions(numParticles);
-        ixn->setUseCutoff(nonbondedCutoff, *data.neighborList);
+        neighborList->computeNeighborList(numParticles, data.posq, noExclusions, boxVectors, data.isPeriodic, nonbondedCutoff, data.threads);
+        ixn->setUseCutoff(nonbondedCutoff, *neighborList);
     }
     map<string, double> globalParameters;
     for (auto& name : globalParameterNames)

platforms/reference/src/ReferenceKernels.cpp

@@ -1474,7 +1474,8 @@ double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool incl
     if (periodic)
         ixn.setPeriodic(extractBoxVectors(context));
     if (nonbondedMethod != NoCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
+        vector<set<int> > empty(context.getSystem().getNumParticles()); // Don't omit exclusions from the neighbor list
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, empty, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
         ixn.setUseCutoff(nonbondedCutoff, *neighborList);
     }
     map<string, double> globalParameters;

tests/TestCustomGBForce.h

@@ -118,6 +118,7 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
             params[1] = 0.1;
             custom->addParticle(params);
         }
+        custom->addExclusion(2*i, 2*i+1);
         positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
         positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
         velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));