Merge branch 'master' into hardwall

openmmapi/include/openmm/internal/vectorize_neon.h

@@ -40,6 +40,10 @@ typedef int int32_t;
 
 // This file defines classes and functions to simplify vectorizing code with NEON.
 
+// These two functions are defined in the vecmath library, which is linked into OpenMM.
+float32x4_t exp_ps(float32x4_t);
+float32x4_t log_ps(float32x4_t);
+
 /**
  * Determine whether ivec4 and fvec4 are supported on this processor.
  */
@@ -262,6 +266,14 @@ static inline fvec4 sqrt(const fvec4& v) {
     return rsqrt(v)*v;
 }
 
+static inline fvec4 exp(const fvec4& v) {
+    return fvec4(exp_ps(v.val));
+}
+
+static inline fvec4 log(const fvec4& v) {
+    return fvec4(log_ps(v.val));
+}
+
 static inline float dot3(const fvec4& v1, const fvec4& v2) {
     fvec4 result = v1*v2;
     return vgetq_lane_f32(result, 0) + vgetq_lane_f32(result, 1) + vgetq_lane_f32(result, 2);

openmmapi/include/openmm/internal/vectorize_pnacl.h

@@ -233,6 +233,14 @@ static inline fvec4 abs(const fvec4& v) {
     return v&(__m128) ivec4(0x7FFFFFFF);
 }
 
+static inline fvec4 exp(const fvec4& v) {
+    return fvec4(expf(v[0]), expf(v[1]), expf(v[2]), expf(v[3]));
+}
+
+static inline fvec4 log(const fvec4& v) {
+    return fvec4(logf(v[0]), logf(v[1]), logf(v[2]), logf(v[3]));
+}
+
 static inline float dot3(const fvec4& v1, const fvec4& v2) {
     fvec4 r = v1*v2;
     return r[0]+r[1]+r[2];

openmmapi/include/openmm/internal/vectorize_sse.h

@@ -37,6 +37,10 @@
 
 // This file defines classes and functions to simplify vectorizing code with SSE.
 
+// These two functions are defined in the vecmath library, which is linked into OpenMM.
+__m128 exp_ps(__m128);
+__m128 log_ps(__m128);
+
 /**
  * Determine whether ivec4 and fvec4 are supported on this processor.
  */
@@ -253,6 +257,14 @@ static inline fvec4 rsqrt(const fvec4& v) {
     return y;
 }
 
+static inline fvec4 exp(const fvec4& v) {
+    return fvec4(exp_ps(v.val));
+}
+
+static inline fvec4 log(const fvec4& v) {
+    return fvec4(log_ps(v.val));
+}
+
 static inline float dot3(const fvec4& v1, const fvec4& v2) {
     return _mm_cvtss_f32(_mm_dp_ps(v1, v2, 0x71));
 }

openmmapi/src/CustomCentroidBondForceImpl.cpp

@@ -150,24 +150,37 @@ ParsedExpression CustomCentroidBondForceImpl::prepareExpression(const CustomCent
     functions["dihedral"] = &dihedral;
     ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
     map<string, int> groups;
+    set<string> variables;
     for (int i = 0; i < force.getNumGroupsPerBond(); i++) {
-        stringstream name;
+        stringstream name, x, y, z;
         name << 'g' << (i+1);
+        x << 'x' << (i+1);
+        y << 'y' << (i+1);
+        z << 'z' << (i+1);
         groups[name.str()] = i;
+        variables.insert(x.str());
+        variables.insert(y.str());
+        variables.insert(z.str());
     }
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), groups, distances, angles, dihedrals)).optimize();
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    for (int i = 0; i < force.getNumPerBondParameters(); i++)
+        variables.insert(force.getPerBondParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), groups, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomCentroidBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> groups,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomCentroidBondForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], groups, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], groups, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

openmmapi/src/CustomCompoundBondForceImpl.cpp

@@ -136,24 +136,37 @@ ParsedExpression CustomCompoundBondForceImpl::prepareExpression(const CustomComp
     functions["dihedral"] = &dihedral;
     ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
     map<string, int> atoms;
+    set<string> variables;
     for (int i = 0; i < force.getNumParticlesPerBond(); i++) {
-        stringstream name;
+        stringstream name, x, y, z;
         name << 'p' << (i+1);
+        x << 'x' << (i+1);
+        y << 'y' << (i+1);
+        z << 'z' << (i+1);
         atoms[name.str()] = i;
+        variables.insert(x.str());
+        variables.insert(y.str());
+        variables.insert(z.str());
     }
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    for (int i = 0; i < force.getNumPerBondParameters(); i++)
+        variables.insert(force.getPerBondParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomCompoundBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomCompoundBondForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

openmmapi/src/CustomExternalForce.cpp

@@ -119,3 +119,7 @@ ForceImpl* CustomExternalForce::createImpl() const {
 void CustomExternalForce::updateParametersInContext(Context& context) {
     dynamic_cast<CustomExternalForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
 }
+
+bool CustomExternalForce::usesPeriodicBoundaryConditions() const {
+    return (energyExpression.find("periodicdistance") != string::npos);
+}

openmmapi/src/CustomHbondForceImpl.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -215,19 +215,28 @@ ParsedExpression CustomHbondForceImpl::prepareExpression(const CustomHbondForce&
     atoms["d1"] = 3;
     atoms["d2"] = 4;
     atoms["d3"] = 5;
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
+    set<string> variables;
+    for (int i = 0; i < force.getNumPerDonorParameters(); i++)
+        variables.insert(force.getPerDonorParameterName(i));
+    for (int i = 0; i < force.getNumPerAcceptorParameters(); i++)
+        variables.insert(force.getPerAcceptorParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomHbondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomHBondForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || op.getNumArguments() < 2)
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

openmmapi/src/CustomManyParticleForceImpl.cpp

@@ -165,24 +165,40 @@ ParsedExpression CustomManyParticleForceImpl::prepareExpression(const CustomMany
     functions["dihedral"] = &dihedral;
     ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
     map<string, int> atoms;
+    set<string> variables;
     for (int i = 0; i < force.getNumParticlesPerSet(); i++) {
-        stringstream name;
+        stringstream name, x, y, z;
         name << 'p' << (i+1);
+        x << 'x' << (i+1);
+        y << 'y' << (i+1);
+        z << 'z' << (i+1);
         atoms[name.str()] = i;
+        variables.insert(x.str());
+        variables.insert(y.str());
+        variables.insert(z.str());
+        for (int j = 0; j < force.getNumPerParticleParameters(); j++) {
+            stringstream param;
+            param << force.getPerParticleParameterName(j) << (i+1);
+            variables.insert(param.str());
         }
-    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        variables.insert(force.getGlobalParameterName(i));
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals, variables)).optimize();
 }
 
 ExpressionTreeNode CustomManyParticleForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
-        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals, set<string>& variables) {
     const Operation& op = node.getOperation();
+    if (op.getId() == Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("CustomManyParticleForce: Unknown variable '"+op.getName()+"'");
     if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
     {
         // This is not an angle or dihedral, so process its children.
 
         vector<ExpressionTreeNode> children;
         for (int i = 0; i < (int) node.getChildren().size(); i++)
-            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
+            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals, variables));
         return ExpressionTreeNode(op.clone(), children);
     }
     const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);

platforms/cpu/include/CpuKernels.h

@@ -92,6 +92,46 @@ private:
     Kernel referenceKernel;
 };
 
+/**
+ * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CpuCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
+public:
+    CpuCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) :
+            CalcHarmonicAngleForceKernel(name, platform), data(data), angleIndexArray(NULL), angleParamArray(NULL) {
+    }
+    ~CpuCalcHarmonicAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param force      the HarmonicAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const HarmonicAngleForce& 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 HarmonicAngleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force);
+private:
+    CpuPlatform::PlatformData& data;
+    int numAngles;
+    int **angleIndexArray;
+    RealOpenMM **angleParamArray;
+    CpuBondForce bondForce;
+};
+
 /**
  * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
  */
@@ -229,6 +269,7 @@ private:
     CpuNeighborList* neighborList;
     CpuNonbondedForce* nonbonded;
     Kernel optimizedPme;
+    CpuBondForce bondForce;
 };
 
 /**

platforms/cpu/include/CpuLangevinDynamics.h

 
-/* Portions copyright (c) 2013 Stanford University and Simbios.
+/* Portions copyright (c) 2013-2015 Stanford University and Simbios.
  * Authors: Peter Eastman
  * Contributors: 
  *
@@ -37,6 +37,7 @@ class CpuLangevinDynamics : public ReferenceStochasticDynamics {
 public:
     class Update1Task;
     class Update2Task;
+    class Update3Task;
     /**
      * Constructor.
      *
@@ -80,9 +81,21 @@ public:
     void updatePart2(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
                      std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& inverseMasses, std::vector<OpenMM::RealVec>& xPrime);
 
+    /**  
+     * Third update
+     * 
+     * @param numberOfAtoms       number of atoms
+     * @param atomCoordinates     atom coordinates
+     * @param velocities          velocities
+     * @param inverseMasses       inverse atom masses
+     */
+    void updatePart3(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities,
+                     std::vector<RealOpenMM>& inverseMasses, std::vector<OpenMM::RealVec>& xPrime);
+
 private:
     void threadUpdate1(int threadIndex);
     void threadUpdate2(int threadIndex);
+    void threadUpdate3(int threadIndex);
     OpenMM::ThreadPool& threads;
     OpenMM::CpuRandom& random;
     std::vector<OpenMM_SFMT::SFMT> threadRandom;

platforms/cpu/include/CpuNeighborList.h

@@ -35,6 +35,7 @@
 #include "AlignedArray.h"
 #include "RealVec.h"
 #include "windowsExportCpu.h"
+#include "openmm/internal/gmx_atomic.h"
 #include "openmm/internal/ThreadPool.h"
 #include <set>
 #include <utility>
@@ -74,6 +75,7 @@ private:
     int numAtoms;
     bool usePeriodic;
     float maxDistance;
+    gmx_atomic_t atomicCounter;
 };
 
 } // namespace OpenMM

platforms/cpu/src/CpuCustomGBForce.cpp

@@ -28,7 +28,7 @@
 #include "SimTKOpenMMUtilities.h"
 #include "ReferenceForce.h"
 #include "CpuCustomGBForce.h"
-#include "gmx_atomic.h"
+#include "openmm/internal/gmx_atomic.h"
 
 using namespace OpenMM;
 using namespace std;

platforms/cpu/src/CpuCustomManyParticleForce.cpp

@@ -32,7 +32,7 @@
 #include "ReferenceTabulatedFunction.h"
 #include "openmm/internal/CustomManyParticleForceImpl.h"
 #include "lepton/CustomFunction.h"
-#include "gmx_atomic.h"
+#include "openmm/internal/gmx_atomic.h"
 
 using namespace OpenMM;
 using namespace std;

platforms/cpu/src/CpuCustomNonbondedForce.cpp

@@ -28,7 +28,7 @@
 #include "SimTKOpenMMUtilities.h"
 #include "ReferenceForce.h"
 #include "CpuCustomNonbondedForce.h"
-#include "gmx_atomic.h"
+#include "openmm/internal/gmx_atomic.h"
 
 using namespace OpenMM;
 using namespace std;

platforms/cpu/src/CpuGBSAOBCForce.cpp

@@ -25,7 +25,7 @@
 #include "CpuGBSAOBCForce.h"
 #include "SimTKOpenMMRealType.h"
 #include "openmm/internal/vectorize.h"
-#include "gmx_atomic.h"
+#include "openmm/internal/gmx_atomic.h"
 #include <algorithm>
 #include <cmath>
 #include <cstdlib>
@@ -279,7 +279,7 @@ void CpuGBSAOBCForce::threadComputeForce(ThreadPool& threads, int threadIndex) {
             fvec4 r = sqrt(r2);
             fvec4 alpha2_ij = radii*bornRadii[atomJ];
             fvec4 D_ij = r2/(4.0f*alpha2_ij);
-            fvec4 expTerm(expf(-D_ij[0]), expf(-D_ij[1]), expf(-D_ij[2]), expf(-D_ij[3]));
+            fvec4 expTerm = exp(-D_ij);
             fvec4 denominator2 = r2 + alpha2_ij*expTerm;
             fvec4 denominator = sqrt(denominator2);
             fvec4 Gpol = (partialChargeI*posJ[3])/denominator; 

platforms/cpu/src/CpuKernelFactory.cpp

@@ -41,6 +41,8 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
     CpuPlatform::PlatformData& data = CpuPlatform::getPlatformData(context);
     if (name == CalcForcesAndEnergyKernel::Name())
         return new CpuCalcForcesAndEnergyKernel(name, platform, data, context);
+    if (name == CalcHarmonicAngleForceKernel::Name())
+        return new CpuCalcHarmonicAngleForceKernel(name, platform, data);
     if (name == CalcPeriodicTorsionForceKernel::Name())
         return new CpuCalcPeriodicTorsionForceKernel(name, platform, data);
     if (name == CalcRBTorsionForceKernel::Name())

platforms/cpu/src/CpuKernels.cpp

@@ -30,6 +30,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "CpuKernels.h"
+#include "ReferenceAngleBondIxn.h"
 #include "ReferenceBondForce.h"
 #include "ReferenceConstraints.h"
 #include "ReferenceKernelFactory.h"
@@ -47,6 +48,7 @@
 #include "RealVec.h"
 #include "lepton/CompiledExpression.h"
 #include "lepton/CustomFunction.h"
+#include "lepton/Operation.h"
 #include "lepton/Parser.h"
 #include "lepton/ParsedExpression.h"
 
@@ -83,6 +85,17 @@ static ReferenceConstraints& extractConstraints(ContextImpl& context) {
     return *(ReferenceConstraints*) data->constraints;
 }
 
+/**
+ * Make sure an expression doesn't use any undefined variables.
+ */
+static void validateVariables(const Lepton::ExpressionTreeNode& node, const set<string>& variables) {
+    const Lepton::Operation& op = node.getOperation();
+    if (op.getId() == Lepton::Operation::VARIABLE && variables.find(op.getName()) == variables.end())
+        throw OpenMMException("Unknown variable in expression: "+op.getName());
+    for (int i = 0; i < (int) node.getChildren().size(); i++)
+        validateVariables(node.getChildren()[i], variables);
+}
+
 /**
  * Compute the kinetic energy of the system, possibly shifting the velocities in time to account
  * for a leapfrog integrator.
@@ -240,6 +253,64 @@ double CpuCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, boo
     return referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().finishComputation(context, includeForce, includeEnergy, groups, valid);
 }
 
+CpuCalcHarmonicAngleForceKernel::~CpuCalcHarmonicAngleForceKernel() {
+    if (angleIndexArray != NULL) {
+        for (int i = 0; i < numAngles; i++) {
+            delete[] angleIndexArray[i];
+            delete[] angleParamArray[i];
+        }
+        delete[] angleIndexArray;
+        delete[] angleParamArray;
+    }
+}
+
+void CpuCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
+    numAngles = force.getNumAngles();
+    angleIndexArray = new int*[numAngles];
+    for (int i = 0; i < numAngles; i++)
+        angleIndexArray[i] = new int[3];
+    angleParamArray = new RealOpenMM*[numAngles];
+    for (int i = 0; i < numAngles; i++)
+        angleParamArray[i] = new RealOpenMM[2];
+    for (int i = 0; i < numAngles; ++i) {
+        int particle1, particle2, particle3;
+        double angle, k;
+        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
+        angleIndexArray[i][0] = particle1;
+        angleIndexArray[i][1] = particle2;
+        angleIndexArray[i][2] = particle3;
+        angleParamArray[i][0] = (RealOpenMM) angle;
+        angleParamArray[i][1] = (RealOpenMM) k;
+    }
+    bondForce.initialize(system.getNumParticles(), numAngles, 3, angleIndexArray, data.threads);
+}
+
+double CpuCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& forceData = extractForces(context);
+    RealOpenMM energy = 0;
+    ReferenceAngleBondIxn angleBond;
+    bondForce.calculateForce(posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, angleBond);
+    return energy;
+}
+
+void CpuCalcHarmonicAngleForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force) {
+    if (numAngles != force.getNumAngles())
+        throw OpenMMException("updateParametersInContext: The number of angles has changed");
+
+    // Record the values.
+
+    for (int i = 0; i < numAngles; ++i) {
+        int particle1, particle2, particle3;
+        double angle, k;
+        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
+        if (particle1 != angleIndexArray[i][0] || particle2 != angleIndexArray[i][1] || particle3 != angleIndexArray[i][2])
+            throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
+        angleParamArray[i][0] = (RealOpenMM) angle;
+        angleParamArray[i][1] = (RealOpenMM) k;
+    }
+}
+
 CpuCalcPeriodicTorsionForceKernel::~CpuCalcPeriodicTorsionForceKernel() {
     if (torsionIndexArray != NULL) {
         for (int i = 0; i < numTorsions; i++) {
@@ -467,6 +538,7 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
         bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
         bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
     }
+    bondForce.initialize(system.getNumParticles(), num14, 2, bonded14IndexArray, data.threads);
     
     // Record other parameters.
     
@@ -527,7 +599,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
     if (nonbondedMethod != NoCutoff) {
         // Determine whether we need to recompute the neighbor list.
         
-        double padding = 0.15*nonbondedCutoff;
+        double padding = 0.25*nonbondedCutoff;
         bool needRecompute = false;
         double closeCutoff2 = 0.25*padding*padding;
         double farCutoff2 = 0.5*padding*padding;
@@ -599,9 +671,8 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
     }
     energy += nonbondedEnergy;
     if (includeDirect) {
-        ReferenceBondForce refBondForce;
         ReferenceLJCoulomb14 nonbonded14;
-        refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
+        bondForce.calculateForce(posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
         if (data.isPeriodic)
             energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
     }
@@ -737,6 +808,14 @@ void CpuCalcCustomNonbondedForceKernel::initialize(const System& system, const C
         globalParameterNames.push_back(force.getGlobalParameterName(i));
         globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
     }
+    set<string> variables;
+    variables.insert("r");
+    for (int i = 0; i < numParameters; i++) {
+        variables.insert(parameterNames[i]+"1");
+        variables.insert(parameterNames[i]+"2");
+    }
+    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    validateVariables(expression.getRootNode(), variables);
 
     // Delete the custom functions.
 
@@ -950,6 +1029,18 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
     vector<vector<Lepton::CompiledExpression> > valueGradientExpressions(force.getNumComputedValues());
     vector<Lepton::CompiledExpression> valueExpressions;
     vector<Lepton::CompiledExpression> energyExpressions;
+    set<string> particleVariables, pairVariables;
+    pairVariables.insert("r");
+    particleVariables.insert("x");
+    particleVariables.insert("y");
+    particleVariables.insert("z");
+    for (int i = 0; i < numPerParticleParameters; i++) {
+        particleVariables.insert(particleParameterNames[i]);
+        pairVariables.insert(particleParameterNames[i]+"1");
+        pairVariables.insert(particleParameterNames[i]+"2");
+    }
+    particleVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
+    pairVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
     for (int i = 0; i < force.getNumComputedValues(); i++) {
         string name, expression;
         CustomGBForce::ComputationType type;
@@ -958,15 +1049,21 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
         valueExpressions.push_back(ex.createCompiledExpression());
         valueTypes.push_back(type);
         valueNames.push_back(name);
-        if (i == 0)
+        if (i == 0) {
             valueDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
+            validateVariables(ex.getRootNode(), pairVariables);
+        }
         else {
             valueGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
             valueGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
             valueGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
             for (int j = 0; j < i; j++)
                 valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
+            validateVariables(ex.getRootNode(), particleVariables);
         }
+        particleVariables.insert(name);
+        pairVariables.insert(name+"1");
+        pairVariables.insert(name+"2");
     }
 
     // Parse the expressions for energy terms.
@@ -988,10 +1085,12 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
                 energyGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
                 energyGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
                 energyGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
+                validateVariables(ex.getRootNode(), particleVariables);
             }
             else {
                 energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").createCompiledExpression());
                 energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").createCompiledExpression());
+                validateVariables(ex.getRootNode(), pairVariables);
             }
         }
     }

platforms/cpu/src/CpuLangevinDynamics.cpp

 
-/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2015 Stanford University and Simbios.
  * Authors: Peter Eastman
  * Contributors: 
  *
@@ -49,6 +49,16 @@ public:
     CpuLangevinDynamics& owner;
 };
 
+class CpuLangevinDynamics::Update3Task : public ThreadPool::Task {
+public:
+    Update3Task(CpuLangevinDynamics& owner) : owner(owner) {
+    }
+    void execute(ThreadPool& threads, int threadIndex) {
+        owner.threadUpdate3(threadIndex);
+    }
+    CpuLangevinDynamics& owner;
+};
+
 CpuLangevinDynamics::CpuLangevinDynamics(int numberOfAtoms, RealOpenMM deltaT, RealOpenMM tau, RealOpenMM temperature, ThreadPool& threads, CpuRandom& random) : 
            ReferenceStochasticDynamics(numberOfAtoms, deltaT, tau, temperature), threads(threads), random(random) {
 }
@@ -92,6 +102,23 @@ void CpuLangevinDynamics::updatePart2(int numberOfAtoms, vector<RealVec>& atomCo
     threads.waitForThreads();
 }
 
+void CpuLangevinDynamics::updatePart3(int numberOfAtoms, vector<RealVec>& atomCoordinates, vector<RealVec>& velocities,
+                                       vector<RealOpenMM>& inverseMasses, vector<RealVec>& xPrime) {
+    // Record the parameters for the threads.
+    
+    this->numberOfAtoms = numberOfAtoms;
+    this->atomCoordinates = &atomCoordinates[0];
+    this->velocities = &velocities[0];
+    this->inverseMasses = &inverseMasses[0];
+    this->xPrime = &xPrime[0];
+    
+    // Signal the threads to start running and wait for them to finish.
+    
+    Update3Task task(*this);
+    threads.execute(task);
+    threads.waitForThreads();
+}
+
 void CpuLangevinDynamics::threadUpdate1(int threadIndex) {
     const RealOpenMM tau = getTau();
     const RealOpenMM vscale = EXP(-getDeltaT()/tau);
@@ -122,3 +149,16 @@ void CpuLangevinDynamics::threadUpdate2(int threadIndex) {
         }
    }
 }
+
+void CpuLangevinDynamics::threadUpdate3(int threadIndex) {
+   const RealOpenMM invStepSize = 1.0/getDeltaT();
+    int start = threadIndex*numberOfAtoms/threads.getNumThreads();
+    int end = (threadIndex+1)*numberOfAtoms/threads.getNumThreads();
+
+   for (int i = start; i < end; ++i)
+       if (inverseMasses[i] != 0.0) {
+            velocities[i] = (xPrime[i]-atomCoordinates[i])*invStepSize;
+            atomCoordinates[i] = xPrime[i];
+       }
+}
+

platforms/cpu/src/CpuNeighborList.cpp

@@ -59,22 +59,25 @@ public:
  */
 class CpuNeighborList::Voxels {
 public:
-    Voxels(int blockSize, float vsy, float vsz, float miny, float maxy, float minz, float maxz, const RealVec* periodicBoxVectors, bool usePeriodic) :
-            blockSize(blockSize), voxelSizeY(vsy), voxelSizeZ(vsz), miny(miny), maxy(maxy), minz(minz), maxz(maxz), periodicBoxVectors(periodicBoxVectors), usePeriodic(usePeriodic) {
-        periodicBoxSize[0] = (float) periodicBoxVectors[0][0];
-        periodicBoxSize[1] = (float) periodicBoxVectors[1][1];
-        periodicBoxSize[2] = (float) periodicBoxVectors[2][2];
-        recipBoxSize[0] = (float) (1/periodicBoxVectors[0][0]);
-        recipBoxSize[1] = (float) (1/periodicBoxVectors[1][1]);
-        recipBoxSize[2] = (float) (1/periodicBoxVectors[2][2]);
-        triclinic = (periodicBoxVectors[0][1] != 0.0 || periodicBoxVectors[0][2] != 0.0 ||
-                     periodicBoxVectors[1][0] != 0.0 || periodicBoxVectors[1][2] != 0.0 ||
-                     periodicBoxVectors[2][0] != 0.0 || periodicBoxVectors[2][1] != 0.0);
+    Voxels(int blockSize, float vsy, float vsz, float miny, float maxy, float minz, float maxz, const RealVec* boxVectors, bool usePeriodic) :
+            blockSize(blockSize), voxelSizeY(vsy), voxelSizeZ(vsz), miny(miny), maxy(maxy), minz(minz), maxz(maxz), usePeriodic(usePeriodic) {
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                periodicBoxVectors[i][j] = (float) boxVectors[i][j];
+        periodicBoxSize[0] = (float) boxVectors[0][0];
+        periodicBoxSize[1] = (float) boxVectors[1][1];
+        periodicBoxSize[2] = (float) boxVectors[2][2];
+        recipBoxSize[0] = (float) (1/boxVectors[0][0]);
+        recipBoxSize[1] = (float) (1/boxVectors[1][1]);
+        recipBoxSize[2] = (float) (1/boxVectors[2][2]);
+        triclinic = (boxVectors[0][1] != 0.0 || boxVectors[0][2] != 0.0 ||
+                     boxVectors[1][0] != 0.0 || boxVectors[1][2] != 0.0 ||
+                     boxVectors[2][0] != 0.0 || boxVectors[2][1] != 0.0);
         if (usePeriodic) {
-            ny = (int) floorf(periodicBoxVectors[1][1]/voxelSizeY+0.5f);
-            nz = (int) floorf(periodicBoxVectors[2][2]/voxelSizeZ+0.5f);
-            voxelSizeY = periodicBoxVectors[1][1]/ny;
-            voxelSizeZ = periodicBoxVectors[2][2]/nz;
+            ny = (int) floorf(boxVectors[1][1]/voxelSizeY+0.5f);
+            nz = (int) floorf(boxVectors[2][2]/voxelSizeZ+0.5f);
+            voxelSizeY = boxVectors[1][1]/ny;
+            voxelSizeZ = boxVectors[2][2]/nz;
         }
         else {
             ny = max(1, (int) floorf((maxy-miny)/voxelSizeY+0.5f));
@@ -282,9 +285,10 @@ public:
                 
                 // Loop over atoms and check to see if they are neighbors of this block.
                 
+                const vector<pair<float, int> >& voxelBins = bins[voxelIndex.y][voxelIndex.z];
                 for (int range = 0; range < numRanges; range++) {
                     for (int item = rangeStart[range]; item < rangeEnd[range]; item++) {
-                        const int sortedIndex = bins[voxelIndex.y][voxelIndex.z][item].second;
+                        const int sortedIndex = voxelBins[item].second;
 
                         // Avoid duplicate entries.
                         if (sortedIndex >= lastSortedIndex)
@@ -361,7 +365,7 @@ private:
     int ny, nz;
     float periodicBoxSize[3], recipBoxSize[3];
     bool triclinic;
-    const RealVec* periodicBoxVectors;
+    float periodicBoxVectors[3][3];
     const bool usePeriodic;
     vector<vector<vector<pair<float, int> > > > bins;
 };
@@ -444,6 +448,7 @@ void CpuNeighborList::computeNeighborList(int numAtoms, const AlignedArray<float
 
     // Signal the threads to start running and wait for them to finish.
     
+    gmx_atomic_set(&atomicCounter, 0);
     threads.resumeThreads();
     threads.waitForThreads();
     
@@ -500,7 +505,11 @@ void CpuNeighborList::threadComputeNeighborList(ThreadPool& threads, int threadI
     vector<int> blockAtoms;
     vector<float> blockAtomX(blockSize), blockAtomY(blockSize), blockAtomZ(blockSize);
     vector<VoxelIndex> atomVoxelIndex;
-    for (int i = threadIndex; i < numBlocks; i += numThreads) {
+    while (true) {
+        int i = gmx_atomic_fetch_add(&atomicCounter, 1);
+        if (i >= numBlocks)
+            break;
+
         // Find the atoms in this block and compute their bounding box.
         
         int firstIndex = blockSize*i;
@@ -532,14 +541,24 @@ void CpuNeighborList::threadComputeNeighborList(ThreadPool& threads, int threadI
 
         // Record the exclusions for this block.
 
+        map<int, char> atomFlags;
         for (int j = 0; j < atomsInBlock; j++) {
             const set<int>& atomExclusions = (*exclusions)[sortedAtoms[firstIndex+j]];
             char mask = 1<<j;
-            for (int k = 0; k < (int) blockNeighbors[i].size(); k++) {
-                int atomIndex = blockNeighbors[i][k];
-                if (atomExclusions.find(atomIndex) != atomExclusions.end())
-                    blockExclusions[i][k] |= mask;
+            for (set<int>::const_iterator iter = atomExclusions.begin(); iter != atomExclusions.end(); ++iter) {
+                map<int, char>::iterator thisAtomFlags = atomFlags.find(*iter);
+                if (thisAtomFlags == atomFlags.end())
+                    atomFlags[*iter] = mask;
+                else
+                    thisAtomFlags->second |= mask;
+            }
         }
+        int numNeighbors = blockNeighbors[i].size();
+        for (int k = 0; k < numNeighbors; k++) {
+            int atomIndex = blockNeighbors[i][k];
+            map<int, char>::iterator thisAtomFlags = atomFlags.find(atomIndex);
+            if (thisAtomFlags != atomFlags.end())
+                blockExclusions[i][k] |= thisAtomFlags->second;
         }
     }
 }