Merge pull request #371 from peastman/cpucustom

Created an optimized CPU implementation of CustomNonbondedForce

platforms/cpu/include/CpuCustomNonbondedForce.h

+
+/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef OPENMM_CPU_CUSTOM_NONBONDED_FORCE_H__
+#define OPENMM_CPU_CUSTOM_NONBONDED_FORCE_H__
+
+#include "AlignedArray.h"
+#include "CpuNeighborList.h"
+#include "openmm/internal/ThreadPool.h"
+#include "openmm/internal/vectorize.h"
+#include "lepton/CompiledExpression.h"
+#include <map>
+#include <set>
+#include <utility>
+#include <vector>
+
+namespace OpenMM {
+    
+class CpuCustomNonbondedForce {
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       CpuCustomNonbondedForce(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression,
+                                   const std::vector<std::string>& parameterNames, const std::vector<std::set<int> >& exclusions, ThreadPool& threads);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~CpuCustomNonbondedForce();
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use a cutoff.
+
+         @param distance            the cutoff distance
+         @param neighbors           the neighbor list to use
+
+         --------------------------------------------------------------------------------------- */
+
+      void setUseCutoff(RealOpenMM distance, const CpuNeighborList& neighbors);
+
+      /**---------------------------------------------------------------------------------------
+
+         Restrict the force to a list of interaction groups.
+
+         @param distance            the cutoff distance
+         @param neighbors           the neighbor list to use
+
+         --------------------------------------------------------------------------------------- */
+
+      void setInteractionGroups(const std::vector<std::pair<std::set<int>, std::set<int> > >& groups);
+
+      /**---------------------------------------------------------------------------------------
+      
+         Set the force to use a switching function.
+      
+         @param distance            the switching distance
+      
+         --------------------------------------------------------------------------------------- */
+      
+      void setUseSwitchingFunction(RealOpenMM distance);
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use periodic boundary conditions.  This requires that a cutoff has
+         already been set, and the smallest side of the periodic box is at least twice the cutoff
+         distance.
+
+         @param boxSize             the X, Y, and Z widths of the periodic box
+
+         --------------------------------------------------------------------------------------- */
+
+      void setPeriodic(OpenMM::RealVec& boxSize);
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom pair ixn
+
+         @param numberOfAtoms    number of atoms
+         @param posq             atom coordinates in float format
+         @param atomCoordinates  atom coordinates
+         @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
+         @param fixedParameters  non atom parameters (not currently used)
+         @param globalParameters the values of global parameters
+         @param forces           force array (forces added)
+         @param totalEnergy      total energy
+         @param threads          the thread pool to use
+
+         --------------------------------------------------------------------------------------- */
+
+    void calculatePairIxn(int numberOfAtoms, float* posq, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                          RealOpenMM* fixedParameters, const std::map<std::string, double>& globalParameters,
+                          std::vector<AlignedArray<float> >& threadForce, bool includeForce, bool includeEnergy, double& totalEnergy);
+private:
+    class ComputeForceTask;
+    class ThreadData;
+
+    bool cutoff;
+    bool useSwitch;
+    bool periodic;
+    const CpuNeighborList* neighborList;
+    RealOpenMM periodicBoxSize[3];
+    RealOpenMM cutoffDistance, switchingDistance;
+    ThreadPool& threads;
+    const std::vector<std::set<int> > exclusions;
+    std::vector<ThreadData*> threadData;
+    std::vector<std::string> paramNames;
+    std::vector<std::pair<int, int> > groupInteractions;
+    std::vector<double> threadEnergy;
+    // The following variables are used to make information accessible to the individual threads.
+    int numberOfAtoms;
+    float* posq;
+    RealVec const* atomCoordinates;
+    RealOpenMM** atomParameters;        
+    const std::map<std::string, double>* globalParameters;
+    std::vector<AlignedArray<float> >* threadForce;
+    bool includeForce, includeEnergy;
+    void* atomicCounter;
+
+    /**
+     * This routine contains the code executed by each thread.
+     */
+    void threadComputeForce(ThreadPool& threads, int threadIndex);
+
+    /**
+     * Calculate the interaction between two atoms.
+     * 
+     * @param atom1            the index of the first atom
+     * @param atom2            the index of the second atom
+     * @param data             workspace for the current thread
+     * @param forces           force array (forces added)
+     * @param totalEnergy      total energy
+     * @param boxSize          the size of the periodic box
+     * @param boxSize          the inverse size of the periodic box
+     */
+    void calculateOneIxn(int atom1, int atom2, ThreadData& data, float* forces, double& totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * Compute the displacement and squared distance between two points, optionally using
+     * periodic boundary conditions.
+     */
+    void getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, const fvec4& boxSize, const fvec4& invBoxSize) const;
+};
+
+class CpuCustomNonbondedForce::ThreadData {
+public:
+    ThreadData(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression, const std::vector<std::string>& parameterNames);
+    Lepton::CompiledExpression energyExpression;
+    Lepton::CompiledExpression forceExpression;
+    std::vector<double*> energyParticleParams;
+    std::vector<double*> forceParticleParams;
+    double* energyR;
+    double* forceR;
+};
+
+} // namespace OpenMM
+
+#endif // OPENMM_CPU_CUSTOM_NONBONDED_FORCE_H__

platforms/cpu/include/CpuKernels.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "CpuBondForce.h"
+#include "CpuCustomNonbondedForce.h"
 #include "CpuGBSAOBCForce.h"
 #include "CpuLangevinDynamics.h"
 #include "CpuNeighborList.h"
@@ -217,6 +218,52 @@ private:
     Kernel optimizedPme;
 };
 
+/**
+ * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
+ */
+class CpuCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
+public:
+    CpuCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data);
+    ~CpuCalcCustomNonbondedForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomNonbondedForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomNonbondedForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomNonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
+private:
+    CpuPlatform::PlatformData& data;
+    int numParticles;
+    double **particleParamArray;
+    double nonbondedCutoff, switchingDistance, periodicBoxSize[3], longRangeCoefficient;
+    bool useSwitchingFunction, hasInitializedLongRangeCorrection;
+    CustomNonbondedForce* forceCopy;
+    std::map<std::string, double> globalParamValues;
+    std::vector<std::set<int> > exclusions;
+    std::vector<std::string> parameterNames, globalParameterNames;
+    std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
+    NonbondedMethod nonbondedMethod;
+    CpuNeighborList* neighborList;
+    CpuCustomNonbondedForce* nonbonded;
+};
+
 /**
  * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
  */
@@ -261,7 +308,7 @@ private:
 class CpuIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
 public:
     CpuIntegrateLangevinStepKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : IntegrateLangevinStepKernel(name, platform),
-        data(data), dynamics(NULL) {
+            data(data), dynamics(NULL) {
     }
     ~CpuIntegrateLangevinStepKernel();
     /**

platforms/cpu/src/CpuCustomNonbondedForce.cpp

+
+/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+#include "ReferenceForce.h"
+#include "CpuCustomNonbondedForce.h"
+#include "gmx_atomic.h"
+
+using namespace OpenMM;
+using namespace std;
+
+class CpuCustomNonbondedForce::ComputeForceTask : public ThreadPool::Task {
+public:
+    ComputeForceTask(CpuCustomNonbondedForce& owner) : owner(owner) {
+    }
+    void execute(ThreadPool& threads, int threadIndex) {
+        owner.threadComputeForce(threads, threadIndex);
+    }
+    CpuCustomNonbondedForce& owner;
+};
+
+CpuCustomNonbondedForce::ThreadData::ThreadData(const Lepton::CompiledExpression& energyExpression, const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames) :
+            energyExpression(energyExpression), forceExpression(forceExpression) {
+    energyR = ReferenceForce::getVariablePointer(this->energyExpression, "r");
+    forceR = ReferenceForce::getVariablePointer(this->forceExpression, "r");
+    for (int i = 0; i < (int) parameterNames.size(); i++) {
+        for (int j = 1; j < 3; j++) {
+            stringstream name;
+            name << parameterNames[i] << j;
+            energyParticleParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, name.str()));
+            forceParticleParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, name.str()));
+        }
+    }
+}
+
+CpuCustomNonbondedForce::CpuCustomNonbondedForce(const Lepton::CompiledExpression& energyExpression,
+            const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, const vector<set<int> >& exclusions,ThreadPool& threads) :
+            cutoff(false), useSwitch(false), periodic(false), paramNames(parameterNames), exclusions(exclusions), threads(threads) {
+    for (int i = 0; i < threads.getNumThreads(); i++)
+        threadData.push_back(new ThreadData(energyExpression, forceExpression, parameterNames));
+}
+
+CpuCustomNonbondedForce::~CpuCustomNonbondedForce() {
+    for (int i = 0; i < (int) threadData.size(); i++)
+        delete threadData[i];
+}
+
+void CpuCustomNonbondedForce::setUseCutoff(RealOpenMM distance, const CpuNeighborList& neighbors) {
+    cutoff = true;
+    cutoffDistance = distance;
+    neighborList = &neighbors;
+  }
+
+void CpuCustomNonbondedForce::setInteractionGroups(const vector<pair<set<int>, set<int> > >& groups) {
+    for (int group = 0; group < (int) groups.size(); group++) {
+        const set<int>& set1 = groups[group].first;
+        const set<int>& set2 = groups[group].second;
+        for (set<int>::const_iterator atom1 = set1.begin(); atom1 != set1.end(); ++atom1) {
+            for (set<int>::const_iterator atom2 = set2.begin(); atom2 != set2.end(); ++atom2) {
+                if (*atom1 == *atom2 || exclusions[*atom1].find(*atom2) != exclusions[*atom1].end())
+                    continue; // This is an excluded interaction.
+                if (*atom1 > *atom2 && set1.find(*atom2) != set1.end() && set2.find(*atom1) != set2.end())
+                    continue; // Both atoms are in both sets, so skip duplicate interactions.
+                groupInteractions.push_back(make_pair(*atom1, *atom2));
+            }
+        }
+    }
+}
+
+void CpuCustomNonbondedForce::setUseSwitchingFunction(RealOpenMM distance) {
+    useSwitch = true;
+    switchingDistance = distance;
+}
+
+void CpuCustomNonbondedForce::setPeriodic(RealVec& boxSize) {
+    assert(cutoff);
+    assert(boxSize[0] >= 2.0*cutoffDistance);
+    assert(boxSize[1] >= 2.0*cutoffDistance);
+    assert(boxSize[2] >= 2.0*cutoffDistance);
+    periodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+}
+
+
+void CpuCustomNonbondedForce::calculatePairIxn(int numberOfAtoms, float* posq, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
+                                             RealOpenMM* fixedParameters, const map<string, double>& globalParameters,
+                                             vector<AlignedArray<float> >& threadForce, bool includeForce, bool includeEnergy, double& totalEnergy) {
+    // Record the parameters for the threads.
+    
+    this->numberOfAtoms = numberOfAtoms;
+    this->posq = posq;
+    this->atomCoordinates = &atomCoordinates[0];
+    this->atomParameters = atomParameters;
+    this->globalParameters = &globalParameters;
+    this->threadForce = &threadForce;
+    this->includeForce = includeForce;
+    this->includeEnergy = includeEnergy;
+    threadEnergy.resize(threads.getNumThreads());
+    gmx_atomic_t counter;
+    gmx_atomic_set(&counter, 0);
+    this->atomicCounter = &counter;
+    
+    // Signal the threads to start running and wait for them to finish.
+    
+    ComputeForceTask task(*this);
+    threads.execute(task);
+    threads.waitForThreads();
+    
+    // Combine the energies from all the threads.
+    
+    if (includeEnergy) {
+        int numThreads = threads.getNumThreads();
+        for (int i = 0; i < numThreads; i++)
+            totalEnergy += threadEnergy[i];
+    }
+}
+
+void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int threadIndex) {
+    // Compute this thread's subset of interactions.
+
+    int numThreads = threads.getNumThreads();
+    threadEnergy[threadIndex] = 0;
+    double& energy = threadEnergy[threadIndex];
+    float* forces = &(*threadForce)[threadIndex][0];
+    ThreadData& data = *threadData[threadIndex];
+    for (map<string, double>::const_iterator iter = globalParameters->begin(); iter != globalParameters->end(); ++iter) {
+        ReferenceForce::setVariable(ReferenceForce::getVariablePointer(data.energyExpression, iter->first), iter->second);
+        ReferenceForce::setVariable(ReferenceForce::getVariablePointer(data.forceExpression, iter->first), iter->second);
+    }
+    fvec4 boxSize(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2], 0);
+    fvec4 invBoxSize((1/periodicBoxSize[0]), (1/periodicBoxSize[1]), (1/periodicBoxSize[2]), 0);
+    if (groupInteractions.size() > 0) {
+        // The user has specified interaction groups, so compute only the requested interactions.
+        
+        while (true) {
+            int i = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (i >= groupInteractions.size())
+                break;
+            int atom1 = groupInteractions[i].first;
+            int atom2 = groupInteractions[i].second;
+            for (int j = 0; j < (int) paramNames.size(); j++) {
+                ReferenceForce::setVariable(data.energyParticleParams[j*2], atomParameters[atom1][j]);
+                ReferenceForce::setVariable(data.energyParticleParams[j*2+1], atomParameters[atom2][j]);
+                ReferenceForce::setVariable(data.forceParticleParams[j*2], atomParameters[atom1][j]);
+                ReferenceForce::setVariable(data.forceParticleParams[j*2+1], atomParameters[atom2][j]);
+            }
+            calculateOneIxn(atom1, atom2, data, forces, energy, boxSize, invBoxSize);
+        }
+    }
+    else if (cutoff) {
+        // We are using a cutoff, so get the interactions from the neighbor list.
+
+        while (true) {
+            int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (blockIndex >= neighborList->getNumBlocks())
+                break;
+            int blockAtom[4];
+            for (int i = 0; i < 4; i++)
+                blockAtom[i] = neighborList->getSortedAtoms()[4*blockIndex+i];
+            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
+            const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
+            for (int i = 0; i < (int) neighbors.size(); i++) {
+                int first = neighbors[i];
+                for (int j = 0; j < (int) paramNames.size(); j++) {
+                    ReferenceForce::setVariable(data.energyParticleParams[j*2], atomParameters[first][j]);
+                    ReferenceForce::setVariable(data.forceParticleParams[j*2], atomParameters[first][j]);
+                }
+                for (int k = 0; k < 4; k++) {
+                    if ((exclusions[i] & (1<<k)) == 0) {
+                        int second = blockAtom[k];
+                        for (int j = 0; j < (int) paramNames.size(); j++) {
+                            ReferenceForce::setVariable(data.energyParticleParams[j*2+1], atomParameters[second][j]);
+                            ReferenceForce::setVariable(data.forceParticleParams[j*2+1], atomParameters[second][j]);
+                        }
+                        calculateOneIxn(first, second, data, forces, energy, boxSize, invBoxSize);
+                    }
+                }
+            }
+        }
+    }
+    else {
+        // Every particle interacts with every other one.
+        
+        while (true) {
+            int ii = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
+            if (ii >= numberOfAtoms)
+                break;
+            for (int jj = ii+1; jj < numberOfAtoms; jj++) {
+                if (exclusions[jj].find(ii) == exclusions[jj].end()) {
+                    for (int j = 0; j < (int) paramNames.size(); j++) {
+                        ReferenceForce::setVariable(data.energyParticleParams[j*2], atomParameters[ii][j]);
+                        ReferenceForce::setVariable(data.energyParticleParams[j*2+1], atomParameters[jj][j]);
+                        ReferenceForce::setVariable(data.forceParticleParams[j*2], atomParameters[ii][j]);
+                        ReferenceForce::setVariable(data.forceParticleParams[j*2+1], atomParameters[jj][j]);
+                    }
+                    calculateOneIxn(ii, jj, data, forces, energy, boxSize, invBoxSize);
+                }
+            }
+        }
+    }
+}
+
+void CpuCustomNonbondedForce::calculateOneIxn(int ii, int jj, ThreadData& data, 
+        float* forces, double& totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
+    // Get deltaR, R2, and R between 2 atoms
+
+    fvec4 deltaR;
+    fvec4 posI(posq+4*ii);
+    fvec4 posJ(posq+4*jj);
+    float r2;
+    getDeltaR(posI, posJ, deltaR, r2, boxSize, invBoxSize);
+    if (cutoff && r2 >= cutoffDistance*cutoffDistance)
+        return;
+    float r = sqrtf(r2);
+
+    // accumulate forces
+
+    ReferenceForce::setVariable(data.energyR, r);
+    ReferenceForce::setVariable(data.forceR, r);
+    double dEdR = (includeForce ? data.forceExpression.evaluate()/r : 0.0);
+    double energy = (includeEnergy ? data.energyExpression.evaluate() : 0.0);
+    if (useSwitch) {
+        if (r > switchingDistance) {
+            RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
+            RealOpenMM switchValue = 1+t*t*t*(-10+t*(15-t*6));
+            RealOpenMM switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
+            dEdR = switchValue*dEdR + energy*switchDeriv/r;
+            energy *= switchValue;
+        }
+    }
+    fvec4 result = deltaR*dEdR;
+    (fvec4(forces+4*ii)+result).store(forces+4*ii);
+    (fvec4(forces+4*jj)-result).store(forces+4*jj);
+
+    // accumulate energies
+
+    totalEnergy += energy;
+}
+
+void CpuCustomNonbondedForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, const fvec4& boxSize, const fvec4& invBoxSize) const {
+    deltaR = posJ-posI;
+    if (periodic) {
+        fvec4 base = round(deltaR*invBoxSize)*boxSize;
+        deltaR = deltaR-base;
+    }
+    r2 = dot3(deltaR, deltaR);
+}

platforms/cpu/src/CpuKernelFactory.cpp

@@ -47,6 +47,8 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
         return new CpuCalcRBTorsionForceKernel(name, platform, data);
     if (name == CalcNonbondedForceKernel::Name())
         return new CpuCalcNonbondedForceKernel(name, platform, data);
+    if (name == CalcCustomNonbondedForceKernel::Name())
+        return new CpuCalcCustomNonbondedForceKernel(name, platform, data);
     if (name == CalcGBSAOBCForceKernel::Name())
         return new CpuCalcGBSAOBCForceKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/cpu/src/CpuKernels.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2013-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -37,12 +37,18 @@
 #include "ReferenceLJCoulomb14.h"
 #include "ReferenceProperDihedralBond.h"
 #include "ReferenceRbDihedralBond.h"
+#include "ReferenceTabulatedFunction.h"
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/CustomNonbondedForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "openmm/internal/vectorize.h"
 #include "RealVec.h"
+#include "lepton/CompiledExpression.h"
+#include "lepton/CustomFunction.h"
+#include "lepton/Parser.h"
+#include "lepton/ParsedExpression.h"
 
 using namespace OpenMM;
 using namespace std;
@@ -618,6 +624,168 @@ void CpuCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
 }
 
+CpuCalcCustomNonbondedForceKernel::CpuCalcCustomNonbondedForceKernel(string name, const Platform& platform, CpuPlatform::PlatformData& data) :
+            CalcCustomNonbondedForceKernel(name, platform), data(data), forceCopy(NULL), neighborList(NULL), nonbonded(NULL) {
+}
+
+CpuCalcCustomNonbondedForceKernel::~CpuCalcCustomNonbondedForceKernel() {
+    if (particleParamArray != NULL) {
+        for (int i = 0; i < numParticles; i++)
+            delete[] particleParamArray[i];
+        delete[] particleParamArray;
+    }
+    if (neighborList != NULL)
+        delete neighborList;
+    if (nonbonded != NULL)
+        delete nonbonded;
+    if (forceCopy != NULL)
+        delete forceCopy;
+}
+
+void CpuCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
+
+    // Record the exclusions.
+
+    numParticles = force.getNumParticles();
+    exclusions.resize(numParticles);
+    for (int i = 0; i < force.getNumExclusions(); i++) {
+        int particle1, particle2;
+        force.getExclusionParticles(i, particle1, particle2);
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+    }
+
+    // Build the arrays.
+
+    int numParameters = force.getNumPerParticleParameters();
+    particleParamArray = new double*[numParticles];
+    for (int i = 0; i < numParticles; i++)
+        particleParamArray[i] = new double[numParameters];
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = parameters[j];
+    }
+    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    nonbondedCutoff = force.getCutoffDistance();
+    if (nonbondedMethod == NoCutoff)
+        useSwitchingFunction = false;
+    else {
+        neighborList = new CpuNeighborList(4);
+        useSwitchingFunction = force.getUseSwitchingFunction();
+        switchingDistance = force.getSwitchingDistance();
+    }
+
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < force.getNumFunctions(); i++)
+        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
+
+    // Parse the various expressions used to calculate the force.
+
+    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
+    Lepton::CompiledExpression energyExpression = expression.createCompiledExpression();
+    Lepton::CompiledExpression forceExpression = expression.differentiate("r").createCompiledExpression();
+    for (int i = 0; i < numParameters; i++)
+        parameterNames.push_back(force.getPerParticleParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+        globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
+    }
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+    
+    // Record information for the long range correction.
+    
+    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic && force.getUseLongRangeCorrection()) {
+        forceCopy = new CustomNonbondedForce(force);
+        hasInitializedLongRangeCorrection = false;
+    }
+    else {
+        longRangeCoefficient = 0.0;
+        hasInitializedLongRangeCorrection = true;
+    }
+    
+    // Record the interaction groups.
+    
+    for (int i = 0; i < force.getNumInteractionGroups(); i++) {
+        set<int> set1, set2;
+        force.getInteractionGroupParameters(i, set1, set2);
+        interactionGroups.push_back(make_pair(set1, set2));
+    }
+    data.isPeriodic = (nonbondedMethod == CutoffPeriodic);
+    nonbonded = new CpuCustomNonbondedForce(energyExpression, forceExpression, parameterNames, exclusions, data.threads);
+    if (interactionGroups.size() > 0)
+        nonbonded->setInteractionGroups(interactionGroups);
+}
+
+double CpuCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& forceData = extractForces(context);
+    RealVec& box = extractBoxSize(context);
+    float floatBoxSize[3] = {(float) box[0], (float) box[1], (float) box[2]};
+    double energy = 0;
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (nonbondedMethod != NoCutoff) {
+        neighborList->computeNeighborList(numParticles, data.posq, exclusions, floatBoxSize, data.isPeriodic, nonbondedCutoff, data.threads);
+        nonbonded->setUseCutoff(nonbondedCutoff, *neighborList);
+    }
+    if (periodic) {
+        double minAllowedSize = 2*nonbondedCutoff;
+        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+        nonbonded->setPeriodic(box);
+    }
+    bool globalParamsChanged = false;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++) {
+        double value = context.getParameter(globalParameterNames[i]);
+        if (globalParamValues[globalParameterNames[i]] != value)
+            globalParamsChanged = true;
+        globalParamValues[globalParameterNames[i]] = value;
+    }
+    if (useSwitchingFunction)
+        nonbonded->setUseSwitchingFunction(switchingDistance);
+    nonbonded->calculatePairIxn(numParticles, &data.posq[0], posData, particleParamArray, 0, globalParamValues, data.threadForce, includeForces, includeEnergy, energy);
+    
+    // Add in the long range correction.
+    
+    if (!hasInitializedLongRangeCorrection || (globalParamsChanged && forceCopy != NULL)) {
+        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
+        hasInitializedLongRangeCorrection = true;
+    }
+    energy += longRangeCoefficient/(box[0]*box[1]*box[2]);
+    return energy;
+}
+
+void CpuCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
+    if (numParticles != force.getNumParticles())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerParticleParameters();
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = parameters[j];
+    }
+    
+    // If necessary, recompute the long range correction.
+    
+    if (forceCopy != NULL) {
+        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(force, context.getOwner());
+        hasInitializedLongRangeCorrection = true;
+        *forceCopy = force;
+    }
+}
+
 CpuCalcGBSAOBCForceKernel::~CpuCalcGBSAOBCForceKernel() {
 }
 

platforms/cpu/src/CpuPlatform.cpp

@@ -62,6 +62,7 @@ CpuPlatform::CpuPlatform() {
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     platformProperties.push_back(CpuThreads());

platforms/cuda/src/CudaExpressionUtilities.cpp

@@ -290,11 +290,19 @@ void CudaExpressionUtilities::processExpression(stringstream& out, const Express
         case Operation::DIVIDE:
         {
             bool haveReciprocal = false;
-            for (int i = 0; i < (int) temps.size(); i++)
-                if (temps[i].first.getOperation().getId() == Operation::RECIPROCAL && temps[i].first.getChildren()[0] == node.getChildren()[1]) {
-                    haveReciprocal = true;
-                    out << getTempName(node.getChildren()[0], temps) << "*" << temps[i].second;
-                }
+            if (node.getChildren()[1].getOperation().getId() == Operation::RECIPROCAL) {
+                for (int i = 0; i < (int) temps.size(); i++)
+                    if (temps[i].first == node.getChildren()[1].getChildren()[1]) {
+                        haveReciprocal = true;
+                        out << getTempName(node.getChildren()[0], temps) << "*" << temps[i].second;
+                    }
+            }
+            if (!haveReciprocal)
+                for (int i = 0; i < (int) temps.size(); i++)
+                    if (temps[i].first.getOperation().getId() == Operation::RECIPROCAL && temps[i].first.getChildren()[0] == node.getChildren()[1]) {
+                        haveReciprocal = true;
+                        out << getTempName(node.getChildren()[0], temps) << "*" << temps[i].second;
+                    }
             if (!haveReciprocal)
                 out << getTempName(node.getChildren()[0], temps) << "/" << getTempName(node.getChildren()[1], temps);
             break;