CPU platform supports multiple NonbondedForces

platforms/cpu/include/CpuKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -268,7 +268,7 @@ public:
 private:
     class PmeIO;
     CpuPlatform::PlatformData& data;
-    int numParticles, num14;
+    int numParticles, num14, posqIndex;
     int **bonded14IndexArray;
     double **bonded14ParamArray;
     double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, ewaldSelfEnergy, dispersionCoefficient;
@@ -277,6 +277,7 @@ private:
     std::vector<std::set<int> > exclusions;
     std::vector<std::pair<float, float> > particleParams;
     std::vector<float> C6params;
+    std::vector<float> charges;
     NonbondedMethod nonbondedMethod;
     CpuNonbondedForce* nonbonded;
     Kernel optimizedPme, optimizedDispersionPme;
@@ -363,7 +364,9 @@ public:
     void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
 private:
     CpuPlatform::PlatformData& data;
+    int posqIndex;
     std::vector<std::pair<float, float> > particleParams;
+    std::vector<float> charges;
     CpuGBSAOBCForce obc;
 };
 

platforms/cpu/include/CpuPlatform.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -92,6 +92,7 @@ public:
     PlatformData(int numParticles, int numThreads, bool deterministicForces);
     ~PlatformData();
     void requestNeighborList(double cutoffDistance, double padding, bool useExclusions, const std::vector<std::set<int> >& exclusionList);
+    int requestPosqIndex();
     AlignedArray<float> posq;
     std::vector<AlignedArray<float> > threadForce;
     ThreadPool threads;
@@ -101,6 +102,7 @@ public:
     CpuNeighborList* neighborList;
     double cutoff, paddedCutoff;
     bool anyExclusions, deterministicForces;
+    int currentPosqIndex, nextPosqIndex;
     std::vector<std::set<int> > exclusions;
 };
 

platforms/cpu/src/CpuKernels.cpp

@@ -138,6 +138,19 @@ static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>&
     return 0.5*energy;
 }
 
+/**
+ * Copy particle charges into the fourth element of the posq array.
+ */
+static void copyChargesToPosq(ContextImpl& context, const vector<float>& charges, int index) {
+    CpuPlatform::PlatformData& data = CpuPlatform::getPlatformData(context);
+    if (index == data.currentPosqIndex)
+        return;
+    data.currentPosqIndex = index;
+    AlignedArray<float>& posq = data.posq;
+    for (int i = 0; i < charges.size(); i++)
+        posq[4*i+3] = charges[i];
+}
+
 CpuCalcForcesAndEnergyKernel::CpuCalcForcesAndEnergyKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data, ContextImpl& context) :
         CalcForcesAndEnergyKernel(name, platform), data(data) {
     // Create a Reference platform version of this kernel.
@@ -530,6 +543,7 @@ CpuCalcNonbondedForceKernel::~CpuCalcNonbondedForceKernel() {
 }
 
 void CpuCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
+    posqIndex = data.requestPosqIndex();
 
     // Identify which exceptions are 1-4 interactions.
 
@@ -556,12 +570,13 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
     for (int i = 0; i < num14; i++)
         bonded14ParamArray[i] = new double[3];
     particleParams.resize(numParticles);
+    charges.resize(numParticles);
     C6params.resize(numParticles);
     double sumSquaredCharges = 0.0;
     for (int i = 0; i < numParticles; ++i) {
         double charge, radius, depth;
         force.getParticleParameters(i, charge, radius, depth);
-        data.posq[4*i+3] = (float) charge;
+        charges[i] = (float) charge;
         particleParams[i] = make_pair((float) (0.5*radius), (float) (2.0*sqrt(depth)));
         C6params[i] = 8.0*pow(particleParams[i].first, 3.0) * particleParams[i].second;
         sumSquaredCharges += charge*charge;
@@ -660,6 +675,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
             }
         }
     }
+    copyChargesToPosq(context, charges, posqIndex);
     AlignedArray<float>& posq = data.posq;
     vector<Vec3>& posData = extractPositions(context);
     vector<Vec3>& forceData = extractForces(context);
@@ -726,11 +742,12 @@ void CpuCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
 
     // Record the values.
 
+    posqIndex = data.requestPosqIndex();
     double sumSquaredCharges = 0.0;
     for (int i = 0; i < numParticles; ++i) {
         double charge, radius, depth;
         force.getParticleParameters(i, charge, radius, depth);
-        data.posq[4*i+3] = (float) charge;
+        charges[i] = (float) charge;
         particleParams[i] = make_pair((float) (0.5*radius), (float) (2.0*sqrt(depth)));
         sumSquaredCharges += charge*charge;
     }
@@ -964,12 +981,14 @@ CpuCalcGBSAOBCForceKernel::~CpuCalcGBSAOBCForceKernel() {
 }
 
 void CpuCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
+    posqIndex = data.requestPosqIndex();
     int numParticles = system.getNumParticles();
     particleParams.resize(numParticles);
+    charges.resize(numParticles);
     for (int i = 0; i < numParticles; ++i) {
         double charge, radius, scalingFactor;
         force.getParticleParameters(i, charge, radius, scalingFactor);
-        data.posq[4*i+3] = (float) charge;
+        charges[i] = (float) charge;
         radius -= 0.009;
         particleParams[i] = make_pair((float) radius, (float) (scalingFactor*radius));
     }
@@ -983,6 +1002,7 @@ void CpuCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCFo
 }
 
 double CpuCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    copyChargesToPosq(context, charges, posqIndex);
     if (data.isPeriodic) {
         Vec3& boxSize = extractBoxSize(context);
         float floatBoxSize[3] = {(float) boxSize[0], (float) boxSize[1], (float) boxSize[2]};
@@ -1000,10 +1020,11 @@ void CpuCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, co
 
     // Record the values.
 
+    posqIndex = data.requestPosqIndex();
     for (int i = 0; i < numParticles; ++i) {
         double charge, radius, scalingFactor;
         force.getParticleParameters(i, charge, radius, scalingFactor);
-        data.posq[4*i+3] = (float) charge;
+        charges[i] = (float) charge;
         radius -= 0.009;
         particleParams[i] = make_pair((float) radius, (float) (scalingFactor*radius));
     }

platforms/cpu/src/CpuPlatform.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2018 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -147,7 +147,7 @@ const CpuPlatform::PlatformData& CpuPlatform::getPlatformData(const ContextImpl&
 }
 
 CpuPlatform::PlatformData::PlatformData(int numParticles, int numThreads, bool deterministicForces) : posq(4*numParticles), threads(numThreads),
-        deterministicForces(deterministicForces), neighborList(NULL), cutoff(0.0), paddedCutoff(0.0), anyExclusions(false) {
+        deterministicForces(deterministicForces), neighborList(NULL), cutoff(0.0), paddedCutoff(0.0), anyExclusions(false), currentPosqIndex(-1), nextPosqIndex(0) {
     numThreads = threads.getNumThreads();
     threadForce.resize(numThreads);
     for (int i = 0; i < numThreads; i++)
@@ -184,3 +184,7 @@ void CpuPlatform::PlatformData::requestNeighborList(double cutoffDistance, doubl
     else if (!anyExclusions)
         exclusions = exclusionList;
 }
+
+int CpuPlatform::PlatformData::requestPosqIndex() {
+    return nextPosqIndex++;
+}
\ No newline at end of file

platforms/cuda/src/CudaKernels.cpp

@@ -1716,7 +1716,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             map<string, string> replacements;
             replacements["CHARGE"] = (usePosqCharges ? "pos.w" : "charges[atom]");
             CUmodule module = cu.createModule(CudaKernelSources::vectorOps+cu.replaceStrings(CudaKernelSources::pme, replacements), pmeDefines);
-            if (cu.getPlatformData().useCpuPme && !doLJPME) {
+            if (cu.getPlatformData().useCpuPme && !doLJPME && usePosqCharges) {
                 // Create the CPU PME kernel.
 
                 try {

platforms/opencl/src/OpenCLKernels.cpp

@@ -1692,7 +1692,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
             bool deviceIsCpu = (cl.getDevice().getInfo<CL_DEVICE_TYPE>() == CL_DEVICE_TYPE_CPU);
             if (deviceIsCpu)
                 pmeDefines["DEVICE_IS_CPU"] = "1";
-            if (cl.getPlatformData().useCpuPme && !doLJPME) {
+            if (cl.getPlatformData().useCpuPme && !doLJPME && usePosqCharges) {
                 // Create the CPU PME kernel.
 
                 try {