Optimizations to updateParametersInContext() for OpenCL

platforms/opencl/include/OpenCLContext.h

@@ -609,6 +609,10 @@ public:
     OpenCLNonbondedUtilities& getNonbondedUtilities() {
         return *nonbonded;
     }
+    /**
+     * Set the particle charges.  These are packed into the fourth element of the posq array.
+     */
+    void setCharges(const std::vector<double>& charges);
     /**
      * Get the thread used by this context for executing parallel computations.
      */
@@ -692,6 +696,12 @@ public:
      * and order to be revalidated.
      */
     void invalidateMolecules();
+    /**
+     * Mark that the current molecule definitions from one particular force (and hence the atom order)
+     * may be invalid.  This should be called whenever force field parameters change.  It will cause the
+     * definitions and order to be revalidated.
+     */
+    bool invalidateMolecules(OpenCLForceInfo* force);
 private:
     struct Molecule;
     struct MoleculeGroup;
@@ -739,6 +749,7 @@ private:
     cl::Kernel clearSixBuffersKernel;
     cl::Kernel reduceReal4Kernel;
     cl::Kernel reduceForcesKernel;
+    cl::Kernel setChargesKernel;
     std::vector<OpenCLForceInfo*> forces;
     std::vector<Molecule> molecules;
     std::vector<MoleculeGroup> moleculeGroups;
@@ -754,6 +765,7 @@ private:
     OpenCLArray* energyBuffer;
     OpenCLArray* energyParamDerivBuffer;
     OpenCLArray* atomIndexDevice;
+    OpenCLArray* chargeBuffer;
     std::vector<std::string> energyParamDerivNames;
     std::map<std::string, double> energyParamDerivWorkspace;
     std::vector<int> atomIndex;

platforms/opencl/include/OpenCLKernels.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) 2008-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -270,9 +270,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force);
 private:
+    class ForceInfo;
     int numBonds;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLArray* params;
 };
@@ -310,9 +312,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomBondForce& force);
 private:
+    class ForceInfo;
     int numBonds;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLParameterSet* params;
     OpenCLArray* globals;
@@ -353,9 +357,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force);
 private:
+    class ForceInfo;
     int numAngles;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLArray* params;
 };
@@ -393,9 +399,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomAngleForce& force);
 private:
+    class ForceInfo;
     int numAngles;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLParameterSet* params;
     OpenCLArray* globals;
@@ -436,9 +444,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
 private:
+    class ForceInfo;
     int numTorsions;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLArray* params;
 };
@@ -476,9 +486,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
 private:
+    class ForceInfo;
     int numTorsions;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLArray* params;
 };
@@ -516,9 +528,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
 private:
+    class ForceInfo;
     int numTorsions;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     std::vector<mm_int2> mapPositionsVec;
     OpenCLArray* coefficients;
@@ -559,9 +573,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force);
 private:
+    class ForceInfo;
     int numTorsions;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLParameterSet* params;
     OpenCLArray* globals;
@@ -625,12 +641,14 @@ private:
         const char* getMaxValue() const {return "(int2) (INT_MAX, INT_MAX)";}
         const char* getSortKey() const {return "value.y";}
     };
+    class ForceInfo;
     class PmeIO;
     class PmePreComputation;
     class PmePostComputation;
     class SyncQueuePreComputation;
     class SyncQueuePostComputation;
     OpenCLContext& cl;
+    ForceInfo* info;
     bool hasInitializedKernel;
     OpenCLArray* sigmaEpsilon;
     OpenCLArray* exceptionParams;
@@ -704,8 +722,10 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
 private:
+    class ForceInfo;
     void initInteractionGroups(const CustomNonbondedForce& force, const std::string& interactionSource, const std::vector<std::string>& tableTypes);
     OpenCLContext& cl;
+    ForceInfo* info;
     OpenCLParameterSet* params;
     OpenCLArray* globals;
     OpenCLArray* interactionGroupData;
@@ -756,10 +776,12 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
 private:
+    class ForceInfo;
     double prefactor, surfaceAreaFactor, cutoff;
     bool hasCreatedKernels;
     int maxTiles;
     OpenCLContext& cl;
+    ForceInfo* info;
     OpenCLArray* params;
     OpenCLArray* bornSum;
     OpenCLArray* longBornSum;
@@ -807,10 +829,12 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
 private:
+    class ForceInfo;
     double cutoff;
     bool hasInitializedKernels, needParameterGradient, needEnergyParamDerivs;
     int maxTiles, numComputedValues;
     OpenCLContext& cl;
+    ForceInfo* info;
     OpenCLParameterSet* params;
     OpenCLParameterSet* computedValues;
     OpenCLParameterSet* energyDerivs;
@@ -864,9 +888,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
 private:
+    class ForceInfo;
     int numParticles;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     const System& system;
     OpenCLParameterSet* params;
     OpenCLArray* globals;
@@ -908,9 +934,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force);
 private:
+    class ForceInfo;
     int numDonors, numAcceptors;
     bool hasInitializedKernel;
     OpenCLContext& cl;
+    ForceInfo* info;
     OpenCLParameterSet* donorParams;
     OpenCLParameterSet* acceptorParams;
     OpenCLArray* globals;
@@ -961,9 +989,11 @@ public:
     void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force);
 
 private:
+    class ForceInfo;
     int numGroups, numBonds;
     bool needEnergyParamDerivs;
     OpenCLContext& cl;
+    ForceInfo* info;
     OpenCLParameterSet* params;
     OpenCLArray* globals;
     OpenCLArray* groupParticles;
@@ -1013,8 +1043,10 @@ public:
     void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
 
 private:
+    class ForceInfo;
     int numBonds;
     OpenCLContext& cl;
+    ForceInfo* info;
     OpenCLParameterSet* params;
     OpenCLArray* globals;
     std::vector<std::string> globalParamNames;
@@ -1059,7 +1091,9 @@ public:
     void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
 
 private:
+    class ForceInfo;
     OpenCLContext& cl;
+    ForceInfo* info;
     bool hasInitializedKernel;
     NonbondedMethod nonbondedMethod;
     int maxNeighborPairs, forceWorkgroupSize, findNeighborsWorkgroupSize;
@@ -1119,9 +1153,11 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const GayBerneForce& force);
 private:
+    class ForceInfo;
     class ReorderListener;
     void sortAtoms();
     OpenCLContext& cl;
+    ForceInfo* info;
     bool hasInitializedKernels;
     int numRealParticles, maxNeighborBlocks;
     GayBerneForce::NonbondedMethod nonbondedMethod;

platforms/opencl/src/OpenCLContext.cpp

@@ -69,8 +69,8 @@ static void CL_CALLBACK errorCallback(const char* errinfo, const void* private_i
 
 OpenCLContext::OpenCLContext(const System& system, int platformIndex, int deviceIndex, const string& precision, OpenCLPlatform::PlatformData& platformData) :
         system(system), time(0.0), platformData(platformData), stepCount(0), computeForceCount(0), stepsSinceReorder(99999), atomsWereReordered(false), posq(NULL),
-        posqCorrection(NULL), velm(NULL), forceBuffers(NULL), longForceBuffer(NULL), energyBuffer(NULL), energyParamDerivBuffer(NULL), atomIndexDevice(NULL), integration(NULL),
-        expression(NULL), bonded(NULL), nonbonded(NULL), thread(NULL) {
+        posqCorrection(NULL), velm(NULL), forceBuffers(NULL), longForceBuffer(NULL), energyBuffer(NULL), energyParamDerivBuffer(NULL), atomIndexDevice(NULL),
+        chargeBuffer(NULL), integration(NULL), expression(NULL), bonded(NULL), nonbonded(NULL), thread(NULL) {
     if (precision == "single") {
         useDoublePrecision = false;
         useMixedPrecision = false;
@@ -309,6 +309,7 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
     reduceReal4Kernel = cl::Kernel(utilities, "reduceReal4Buffer");
     if (supports64BitGlobalAtomics)
         reduceForcesKernel = cl::Kernel(utilities, "reduceForces");
+    setChargesKernel = cl::Kernel(utilities, "setCharges");
 
     // Decide whether native_sqrt(), native_rsqrt(), and native_recip() are sufficiently accurate to use.
 
@@ -439,6 +440,8 @@ OpenCLContext::~OpenCLContext() {
         delete energyParamDerivBuffer;
     if (atomIndexDevice != NULL)
         delete atomIndexDevice;
+    if (chargeBuffer != NULL)
+        delete chargeBuffer;
     if (integration != NULL)
         delete integration;
     if (expression != NULL)
@@ -747,6 +750,28 @@ void OpenCLContext::reduceBuffer(OpenCLArray& array, int numBuffers) {
     executeKernel(reduceReal4Kernel, bufferSize, 128);
 }
 
+void OpenCLContext::setCharges(const vector<double>& charges) {
+    if (chargeBuffer == NULL)
+        chargeBuffer = new OpenCLArray(*this, numAtoms, useDoublePrecision ? sizeof(double) : sizeof(float), "chargeBuffer");
+    if (getUseDoublePrecision()) {
+        double* c = (double*) getPinnedBuffer();
+        for (int i = 0; i < charges.size(); i++)
+            c[i] = charges[i];
+        chargeBuffer->upload(c);
+    }
+    else {
+        float* c = (float*) getPinnedBuffer();
+        for (int i = 0; i < charges.size(); i++)
+            c[i] = (float) charges[i];
+        chargeBuffer->upload(c);
+    }
+    setChargesKernel.setArg<cl::Buffer>(0, chargeBuffer->getDeviceBuffer());
+    setChargesKernel.setArg<cl::Buffer>(1, posq->getDeviceBuffer());
+    setChargesKernel.setArg<cl::Buffer>(2, atomIndexDevice->getDeviceBuffer());
+    setChargesKernel.setArg<cl_int>(3, numAtoms);
+    executeKernel(setChargesKernel, numAtoms);
+}
+
 /**
  * This class ensures that atom reordering doesn't break virtual sites.
  */
@@ -945,9 +970,19 @@ void OpenCLContext::findMoleculeGroups() {
 }
 
 void OpenCLContext::invalidateMolecules() {
+    for (int i = 0; i < forces.size(); i++)
+        if (invalidateMolecules(forces[i]))
+            return;
+}
+
+bool OpenCLContext::invalidateMolecules(OpenCLForceInfo* force) {
     if (numAtoms == 0 || nonbonded == NULL || !nonbonded->getUseCutoff())
-        return;
+        return false;
     bool valid = true;
+    int forceIndex = -1;
+    for (int i = 0; i < forces.size(); i++)
+        if (forces[i] == force)
+            forceIndex = i;
     for (int group = 0; valid && group < (int) moleculeGroups.size(); group++) {
         MoleculeGroup& mol = moleculeGroups[group];
         vector<int>& instances = mol.instances;
@@ -962,22 +997,21 @@ void OpenCLContext::invalidateMolecules() {
             Molecule& m2 = molecules[instances[j]];
             int offset2 = offsets[j];
             for (int i = 0; i < (int) atoms.size() && valid; i++) {
-                for (int k = 0; k < (int) forces.size(); k++)
-                    if (!forces[k]->areParticlesIdentical(atoms[i]+offset1, atoms[i]+offset2))
-                        valid = false;
+                if (!force->areParticlesIdentical(atoms[i]+offset1, atoms[i]+offset2))
+                    valid = false;
             }
 
             // See if the force groups are identical.
 
-            for (int i = 0; i < (int) forces.size() && valid; i++) {
-                for (int k = 0; k < (int) m1.groups[i].size() && valid; k++)
-                    if (!forces[i]->areGroupsIdentical(m1.groups[i][k], m2.groups[i][k]))
+            if (valid && forceIndex > -1) {
+                for (int k = 0; k < (int) m1.groups[forceIndex].size() && valid; k++)
+                    if (!force->areGroupsIdentical(m1.groups[forceIndex][k], m2.groups[forceIndex][k]))
                         valid = false;
             }
         }
     }
     if (valid)
-        return;
+        return false;
 
     // The list of which molecules are identical is no longer valid.  We need to restore the
     // atoms to their original order, rebuild the list of identical molecules, and sort them
@@ -1045,6 +1079,7 @@ void OpenCLContext::invalidateMolecules() {
     for (int i = 0; i < (int) reorderListeners.size(); i++)
         reorderListeners[i]->execute();
     reorderAtoms();
+    return true;
 }
 
 void OpenCLContext::reorderAtoms() {

platforms/opencl/src/OpenCLKernels.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) 2008-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -529,9 +529,9 @@ void OpenCLVirtualSitesKernel::computePositions(ContextImpl& context) {
     cl.getIntegrationUtilities().computeVirtualSites();
 }
 
-class OpenCLHarmonicBondForceInfo : public OpenCLForceInfo {
+class OpenCLCalcHarmonicBondForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLHarmonicBondForceInfo(const HarmonicBondForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const HarmonicBondForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumBonds();
@@ -581,7 +581,8 @@ void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const H
     replacements["COMPUTE_FORCE"] = OpenCLKernelSources::harmonicBondForce;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float2");
     cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::bondForce, replacements), force.getForceGroup());
-    cl.addForce(new OpenCLHarmonicBondForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 }
 
 double OpenCLCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -610,12 +611,12 @@ void OpenCLCalcHarmonicBondForceKernel::copyParametersToContext(ContextImpl& con
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomBondForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomBondForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomBondForceInfo(const CustomBondForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const CustomBondForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumBonds();
@@ -667,7 +668,8 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
             paramVector[i][j] = (cl_float) parameters[j];
     }
     params->setParameterValues(paramVector);
-    cl.addForce(new OpenCLCustomBondForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 
     // Record information for the expressions.
 
@@ -761,12 +763,12 @@ void OpenCLCalcCustomBondForceKernel::copyParametersToContext(ContextImpl& conte
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLHarmonicAngleForceInfo : public OpenCLForceInfo {
+class OpenCLCalcHarmonicAngleForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLHarmonicAngleForceInfo(const HarmonicAngleForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const HarmonicAngleForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumAngles();
@@ -818,7 +820,8 @@ void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const
     replacements["COMPUTE_FORCE"] = OpenCLKernelSources::harmonicAngleForce;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float2");
     cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::angleForce, replacements), force.getForceGroup());
-    cl.addForce(new OpenCLHarmonicAngleForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 }
 
 double OpenCLCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -847,12 +850,12 @@ void OpenCLCalcHarmonicAngleForceKernel::copyParametersToContext(ContextImpl& co
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomAngleForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomAngleForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomAngleForceInfo(const CustomAngleForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const CustomAngleForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumAngles();
@@ -905,7 +908,8 @@ void OpenCLCalcCustomAngleForceKernel::initialize(const System& system, const Cu
             paramVector[i][j] = (cl_float) parameters[j];
     }
     params->setParameterValues(paramVector);
-    cl.addForce(new OpenCLCustomAngleForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 
     // Record information for the expressions.
 
@@ -999,12 +1003,12 @@ void OpenCLCalcCustomAngleForceKernel::copyParametersToContext(ContextImpl& cont
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLPeriodicTorsionForceInfo : public OpenCLForceInfo {
+class OpenCLCalcPeriodicTorsionForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLPeriodicTorsionForceInfo(const PeriodicTorsionForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const PeriodicTorsionForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumTorsions();
@@ -1057,7 +1061,8 @@ void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, cons
     replacements["COMPUTE_FORCE"] = OpenCLKernelSources::periodicTorsionForce;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float4");
     cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::torsionForce, replacements), force.getForceGroup());
-    cl.addForce(new OpenCLPeriodicTorsionForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 }
 
 double OpenCLCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1086,12 +1091,12 @@ void OpenCLCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl&
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLRBTorsionForceInfo : public OpenCLForceInfo {
+class OpenCLCalcRBTorsionForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLRBTorsionForceInfo(const RBTorsionForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const RBTorsionForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumTorsions();
@@ -1144,7 +1149,8 @@ void OpenCLCalcRBTorsionForceKernel::initialize(const System& system, const RBTo
     replacements["COMPUTE_FORCE"] = OpenCLKernelSources::rbTorsionForce;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float8");
     cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::torsionForce, replacements), force.getForceGroup());
-    cl.addForce(new OpenCLRBTorsionForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 }
 
 double OpenCLCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1173,12 +1179,12 @@ void OpenCLCalcRBTorsionForceKernel::copyParametersToContext(ContextImpl& contex
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCMAPTorsionForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCMAPTorsionForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCMAPTorsionForceInfo(const CMAPTorsionForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const CMAPTorsionForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumTorsions();
@@ -1257,7 +1263,8 @@ void OpenCLCalcCMAPTorsionForceKernel::initialize(const System& system, const CM
     replacements["MAP_POS"] = cl.getBondedUtilities().addArgument(mapPositions->getDeviceBuffer(), "int2");
     replacements["MAPS"] = cl.getBondedUtilities().addArgument(torsionMaps->getDeviceBuffer(), "int");
     cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::cmapTorsionForce, replacements), force.getForceGroup());
-    cl.addForce(new OpenCLCMAPTorsionForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 }
 
 double OpenCLCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1307,9 +1314,9 @@ void OpenCLCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& cont
     torsionMaps->upload(torsionMapsVec);
 }
 
-class OpenCLCustomTorsionForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomTorsionForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomTorsionForceInfo(const CustomTorsionForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const CustomTorsionForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumTorsions();
@@ -1363,7 +1370,8 @@ void OpenCLCalcCustomTorsionForceKernel::initialize(const System& system, const
             paramVector[i][j] = (cl_float) parameters[j];
     }
     params->setParameterValues(paramVector);
-    cl.addForce(new OpenCLCustomTorsionForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 
     // Record information for the expressions.
 
@@ -1457,12 +1465,12 @@ void OpenCLCalcCustomTorsionForceKernel::copyParametersToContext(ContextImpl& co
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLNonbondedForceInfo : public OpenCLForceInfo {
+class OpenCLCalcNonbondedForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLNonbondedForceInfo(int requiredBuffers, const NonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    ForceInfo(int requiredBuffers, const NonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
         double charge1, charge2, sigma1, sigma2, epsilon1, epsilon2;
@@ -1910,7 +1918,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         replacements["PARAMS"] = cl.getBondedUtilities().addArgument(exceptionParams->getDeviceBuffer(), "float4");
         cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::nonbondedExceptions, replacements), force.getForceGroup());
     }
-    cl.addForce(new OpenCLNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
+    info = new ForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force);
+    cl.addForce(info);
 }
 
 double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
@@ -2149,25 +2158,17 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
     
     // Record the per-particle parameters.
     
-    OpenCLArray& posq = cl.getPosq();
-    posq.download(cl.getPinnedBuffer());
-    mm_float4* posqf = (mm_float4*) cl.getPinnedBuffer();
-    mm_double4* posqd = (mm_double4*) cl.getPinnedBuffer();
+    vector<double> chargeVector(cl.getNumAtoms());
     vector<mm_float2> sigmaEpsilonVector(cl.getPaddedNumAtoms(), mm_float2(0,0));
     double sumSquaredCharges = 0.0;
-    const vector<cl_int>& order = cl.getAtomIndex();
     for (int i = 0; i < force.getNumParticles(); i++) {
-        int index = order[i];
         double charge, sigma, epsilon;
-        force.getParticleParameters(index, charge, sigma, epsilon);
-        if (cl.getUseDoublePrecision())
-            posqd[i].w = charge;
-        else
-            posqf[i].w = (float) charge;
-        sigmaEpsilonVector[index] = mm_float2((float) (0.5*sigma), (float) (2.0*sqrt(epsilon)));
+        force.getParticleParameters(i, charge, sigma, epsilon);
+        chargeVector[i] = charge;
+        sigmaEpsilonVector[i] = mm_float2((float) (0.5*sigma), (float) (2.0*sqrt(epsilon)));
         sumSquaredCharges += charge*charge;
     }
-    posq.upload(cl.getPinnedBuffer());
+    cl.setCharges(chargeVector);
     sigmaEpsilon->upload(sigmaEpsilonVector);
     
     // Record the exceptions.
@@ -2189,7 +2190,7 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
         ewaldSelfEnergy = (cl.getContextIndex() == 0 ? -ONE_4PI_EPS0*alpha*sumSquaredCharges/sqrt(M_PI) : 0.0);
     if (force.getUseDispersionCorrection() && cl.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME))
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
 void OpenCLCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
@@ -2205,9 +2206,9 @@ void OpenCLCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, in
     }
 }
 
-class OpenCLCustomNonbondedForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomNonbondedForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomNonbondedForceInfo(int requiredBuffers, const CustomNonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    ForceInfo(int requiredBuffers, const CustomNonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
         if (force.getNumInteractionGroups() > 0) {
             groupsForParticle.resize(force.getNumParticles());
             for (int i = 0; i < force.getNumInteractionGroups(); i++) {
@@ -2384,7 +2385,8 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
             cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(cl_float), globals->getDeviceBuffer()));
         }
     }
-    cl.addForce(new OpenCLCustomNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
+    info = new ForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force);
+    cl.addForce(info);
     
     // Record information for the long range correction.
     
@@ -2723,12 +2725,12 @@ void OpenCLCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl&
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLGBSAOBCForceInfo : public OpenCLForceInfo {
+class OpenCLCalcGBSAOBCForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLGBSAOBCForceInfo(int requiredBuffers, const GBSAOBCForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    ForceInfo(int requiredBuffers, const GBSAOBCForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
         double charge1, charge2, radius1, radius2, scale1, scale2;
@@ -2806,7 +2808,8 @@ void OpenCLCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOB
     nb.addInteraction(useCutoff, usePeriodic, false, cutoff, vector<vector<int> >(), source, force.getForceGroup());
     nb.addParameter(OpenCLNonbondedUtilities::ParameterInfo("obcParams", "float", 2, sizeof(cl_float2), params->getDeviceBuffer()));;
     nb.addParameter(OpenCLNonbondedUtilities::ParameterInfo("bornForce", "real", 1, elementSize, bornForce->getDeviceBuffer()));;
-    cl.addForce(new OpenCLGBSAOBCForceInfo(nb.getNumForceBuffers(), force));
+    info = new ForceInfo(nb.getNumForceBuffers(), force);
+    cl.addForce(info);
 }
 
 double OpenCLCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -2938,33 +2941,27 @@ void OpenCLCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context,
     
     // Record the per-particle parameters.
     
-    OpenCLArray& posq = cl.getPosq();
-    mm_float4* posqf = (mm_float4*) cl.getPinnedBuffer();
-    mm_double4* posqd = (mm_double4*) cl.getPinnedBuffer();
-    posq.download(cl.getPinnedBuffer());
+    vector<double> chargeVector(cl.getNumAtoms());
     vector<mm_float2> paramsVector(cl.getPaddedNumAtoms(), mm_float2(1,1));
     const double dielectricOffset = 0.009;
     for (int i = 0; i < numParticles; i++) {
         double charge, radius, scalingFactor;
         force.getParticleParameters(i, charge, radius, scalingFactor);
+        chargeVector[i] = charge;
         radius -= dielectricOffset;
         paramsVector[i] = mm_float2((float) radius, (float) (scalingFactor*radius));
-        if (cl.getUseDoublePrecision())
-            posqd[i].w = charge;
-        else
-            posqf[i].w = (float) charge;
     }
-    posq.upload(cl.getPinnedBuffer());
+    cl.setCharges(chargeVector);
     params->upload(paramsVector);
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomGBForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomGBForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomGBForceInfo(int requiredBuffers, const CustomGBForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    ForceInfo(int requiredBuffers, const CustomGBForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
         vector<double> params1;
@@ -3812,7 +3809,8 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
         for (int i = 0; i < (int) arguments.size(); i++)
             cl.getNonbondedUtilities().addArgument(arguments[i]);
     }
-    cl.addForce(new OpenCLCustomGBForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
+    info = new ForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force);
+    cl.addForce(info);
     if (useLong)
         cl.addAutoclearBuffer(*longEnergyDerivs);
     else {
@@ -4073,12 +4071,12 @@ void OpenCLCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& context
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomExternalForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomExternalForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomExternalForceInfo(const CustomExternalForce& force, int numParticles) : OpenCLForceInfo(0), force(force), indices(numParticles, -1) {
+    ForceInfo(const CustomExternalForce& force, int numParticles) : OpenCLForceInfo(0), force(force), indices(numParticles, -1) {
         vector<double> params;
         for (int i = 0; i < force.getNumParticles(); i++) {
             int particle;
@@ -4133,7 +4131,8 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
             paramVector[i][j] = (cl_float) parameters[j];
     }
     params->setParameterValues(paramVector);
-    cl.addForce(new OpenCLCustomExternalForceInfo(force, system.getNumParticles()));
+    info = new ForceInfo(force, system.getNumParticles());
+    cl.addForce(info);
 
     // Record information for the expressions.
 
@@ -4228,12 +4227,12 @@ void OpenCLCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& c
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomHbondForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomHbondForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomHbondForceInfo(int requiredBuffers, const CustomHbondForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    ForceInfo(int requiredBuffers, const CustomHbondForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
         return true;
@@ -4406,7 +4405,8 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
         maxBuffers = max(maxBuffers, donorBufferCounter[i]);
     for (int i = 0; i < (int) acceptorBufferCounter.size(); i++)
         maxBuffers = max(maxBuffers, acceptorBufferCounter[i]);
-    cl.addForce(new OpenCLCustomHbondForceInfo(maxBuffers, force));
+    info = new ForceInfo(maxBuffers, force);
+    cl.addForce(info);
 
     // Record exclusions.
 
@@ -4767,12 +4767,12 @@ void OpenCLCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& cont
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomCentroidBondForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomCentroidBondForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomCentroidBondForceInfo(const CustomCentroidBondForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const CustomCentroidBondForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumBonds();
@@ -4840,7 +4840,8 @@ void OpenCLCalcCustomCentroidBondForceKernel::initialize(const System& system, c
         return;
     if (!cl.getSupports64BitGlobalAtomics())
         throw OpenMMException("CustomCentroidBondForce requires a device that supports 64 bit atomic operations");
-    cl.addForce(new OpenCLCustomCentroidBondForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
     
     // Record the groups.
     
@@ -5229,12 +5230,12 @@ void OpenCLCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImp
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomCompoundBondForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomCompoundBondForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomCompoundBondForceInfo(const CustomCompoundBondForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const CustomCompoundBondForce& force) : OpenCLForceInfo(0), force(force) {
     }
     int getNumParticleGroups() {
         return force.getNumBonds();
@@ -5285,7 +5286,8 @@ void OpenCLCalcCustomCompoundBondForceKernel::initialize(const System& system, c
             paramVector[i][j] = (cl_float) parameters[j];
     }
     params->setParameterValues(paramVector);
-    cl.addForce(new OpenCLCustomCompoundBondForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 
     // Record the tabulated functions.
 
@@ -5543,12 +5545,12 @@ void OpenCLCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImp
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLCustomManyParticleForceInfo : public OpenCLForceInfo {
+class OpenCLCalcCustomManyParticleForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLCustomManyParticleForceInfo(const CustomManyParticleForce& force) : OpenCLForceInfo(0), force(force) {
+    ForceInfo(const CustomManyParticleForce& force) : OpenCLForceInfo(0), force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
         vector<double> params1, params2;
@@ -5635,7 +5637,8 @@ void OpenCLCalcCustomManyParticleForceKernel::initialize(const System& system, c
             paramVector[i][j] = (float) parameters[j];
     }
     params->setParameterValues(paramVector);
-    cl.addForce(new OpenCLCustomManyParticleForceInfo(force));
+    info = new ForceInfo(force);
+    cl.addForce(info);
 
     // Record the tabulated functions.
 
@@ -6209,12 +6212,12 @@ void OpenCLCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImp
     
     // Mark that the current reordering may be invalid.
     
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
 }
 
-class OpenCLGayBerneForceInfo : public OpenCLForceInfo {
+class OpenCLCalcGayBerneForceKernel::ForceInfo : public OpenCLForceInfo {
 public:
-    OpenCLGayBerneForceInfo(int requiredBuffers, const GayBerneForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    ForceInfo(int requiredBuffers, const GayBerneForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
         int xparticle1, yparticle1;
@@ -6426,7 +6429,8 @@ void OpenCLCalcGayBerneForceKernel::initialize(const System& system, const GayBe
     neighborsKernel = cl::Kernel(program, "findNeighbors");
     forceKernel = cl::Kernel(program, "computeForce");
     torqueKernel = cl::Kernel(program, "applyTorques");
-    cl.addForce(new OpenCLGayBerneForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
+    info = new ForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force);
+    cl.addForce(info);
     cl.addReorderListener(new ReorderListener(*this));
 }
 
@@ -6577,7 +6581,7 @@ void OpenCLCalcGayBerneForceKernel::copyParametersToContext(ContextImpl& context
         }
         exceptionParams->upload(exceptionParamsVec);
     }
-    cl.invalidateMolecules();
+    cl.invalidateMolecules(info);
     sortAtoms();
 }
 

platforms/opencl/src/kernels/utilities.cl

@@ -107,3 +107,11 @@ __kernel void determineNativeAccuracy(__global float8* restrict values, int numV
         values[i] = (float8) (v, native_sqrt(v), native_rsqrt(v), native_recip(v), native_exp(v), native_log(v), 0.0f, 0.0f);
     }
 }
+
+/**
+ * Record the atomic charges into the posq array.
+ */
+__kernel void setCharges(__global real* restrict charges, __global real4* restrict posq, __global int* restrict atomOrder, int numAtoms) {
+    for (int i = get_global_id(0); i < numAtoms; i += get_global_size(0))
+        posq[i].w = charges[atomOrder[i]];
+}
\ No newline at end of file