Merge branch 'master' of https://github.com/SimTk/openmm

platforms/cpu/include/CpuNeighborList.h

+#ifndef OPENMM_CPU_NEIGHBORLIST_H_
+#define OPENMM_CPU_NEIGHBORLIST_H_
+
+#include "windowsExportCpu.h"
+#include <pthread.h>
+#include <set>
+#include <utility>
+#include <vector>
+
+namespace OpenMM {
+    
+class OPENMM_EXPORT_CPU CpuNeighborList {
+public:
+    class ThreadData;
+    class VoxelHash;
+    CpuNeighborList();
+    ~CpuNeighborList();
+    void computeNeighborList(int numAtoms, const std::vector<float>& atomLocations, const std::vector<std::set<int> >& exclusions,
+            const float* periodicBoxSize, bool usePeriodic, float maxDistance);
+    const std::vector<std::pair<int, int> >& getNeighbors();
+    /**
+     * This routine contains the code executed by each thread.
+     */
+    void runThread(int index, std::vector<std::pair<int, int> >& threadNeighbors);
+private:
+    bool isDeleted;
+    int numThreads, waitCount;
+    std::vector<std::pair<int, int> > neighbors;
+    std::vector<pthread_t> thread;
+    std::vector<ThreadData*> threadData;
+    pthread_cond_t startCondition, endCondition;
+    pthread_mutex_t lock;
+    // The following variables are used to make information accessible to the individual threads.
+    VoxelHash* voxelHash;
+    const std::vector<std::set<int> >* exclusions;
+    const float* atomLocations;
+    const float* periodicBoxSize;
+    int numAtoms;
+    bool usePeriodic;
+    float maxDistance;
+};
+
+} // namespace OpenMM
+
+#endif // OPENMM_REFERENCE_NEIGHBORLIST_H_

platforms/reference/include/ReferenceShakeAlgorithm.h → platforms/cpu/include/CpuNonbondedForce.h

 
-/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -22,44 +22,28 @@
  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  */
 
-#ifndef __ReferenceShakeAlgorithm_H__
-#define __ReferenceShakeAlgorithm_H__
-
-#include "ReferenceConstraintAlgorithm.h"
+#ifndef OPENMM_CPU_NONBONDED_FORCE_H__
+#define OPENMM_CPU_NONBONDED_FORCE_H__
 
+#include "ReferencePairIxn.h"
+#include <pthread.h>
+#include <set>
+#include <utility>
+#include <vector>
+#include <smmintrin.h>
 // ---------------------------------------------------------------------------------------
 
-class ReferenceShakeAlgorithm : public ReferenceConstraintAlgorithm {
-
-   protected:
-
-      int _maximumNumberOfIterations;
-      RealOpenMM _tolerance;
-
-      int _numberOfConstraints;
-      int** _atomIndices;
-      RealOpenMM* _distance;
-
-      RealOpenMM** _r_ij;
-      RealOpenMM* _d_ij2;
-      RealOpenMM* _distanceTolerance;
-      RealOpenMM* _reducedMasses;
-      bool _hasInitializedMasses;
-      
+class CpuNonbondedForce {
     public:
+        class ThreadData;
 
       /**---------------------------------------------------------------------------------------
       
-         ReferenceShakeAlgorithm constructor
-      
-         @param numberOfConstraints      number of constraints 
-         @param atomIndices              atom indices for contraints
-         @param distance                 distances for constraints
-         @param tolerance                constraint tolerance
+         Constructor
       
          --------------------------------------------------------------------------------------- */
 
-      ReferenceShakeAlgorithm( int numberOfConstraints, int** atomIndices, RealOpenMM* distance, RealOpenMM tolerance );
+       CpuNonbondedForce();
 
       /**---------------------------------------------------------------------------------------
       
@@ -67,107 +51,187 @@ class ReferenceShakeAlgorithm : public ReferenceConstraintAlgorithm {
       
          --------------------------------------------------------------------------------------- */
 
-       ~ReferenceShakeAlgorithm( );
+       ~CpuNonbondedForce();
 
       /**---------------------------------------------------------------------------------------
       
-         Get number of constraints
+         Set the force to use a cutoff.
       
-         @return number of constraints
+         @param distance            the cutoff distance
+         @param neighbors           the neighbor list to use
+         @param solventDielectric   the dielectric constant of the bulk solvent
       
          --------------------------------------------------------------------------------------- */
       
-      int getNumberOfConstraints( void ) const;
+      void setUseCutoff(float distance, const std::vector<std::pair<int, int> >& neighbors, float solventDielectric);
 
       /**---------------------------------------------------------------------------------------
       
-         Get maximum number of iterations
+         Set the force to use a switching function on the Lennard-Jones interaction.
       
-         @return maximum number of iterations
+         @param distance            the switching distance
       
          --------------------------------------------------------------------------------------- */
       
-      int getMaximumNumberOfIterations( void ) const;
+      void setUseSwitchingFunction(float distance);
       
       /**---------------------------------------------------------------------------------------
       
-         Set maximum number of iterations
+         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 maximumNumberOfIterations   new maximum number of iterations
+         @param boxSize             the X, Y, and Z widths of the periodic box
       
          --------------------------------------------------------------------------------------- */
       
-      void setMaximumNumberOfIterations( int maximumNumberOfIterations );
+      void setPeriodic(float* periodicBoxSize);
        
       /**---------------------------------------------------------------------------------------
       
-         Get tolerance
+         Set the force to use Ewald summation.
       
-         @return tolerance
+         @param alpha  the Ewald separation parameter
+         @param kmaxx  the largest wave vector in the x direction
+         @param kmaxy  the largest wave vector in the y direction
+         @param kmaxz  the largest wave vector in the z direction
       
          --------------------------------------------------------------------------------------- */
       
-      RealOpenMM getTolerance( void ) const;
+      void setUseEwald(float alpha, int kmaxx, int kmaxy, int kmaxz);
+
      
       /**---------------------------------------------------------------------------------------
       
-         Set tolerance
+         Set the force to use Particle-Mesh Ewald (PME) summation.
       
-         @param tolerance new tolerance
+         @param alpha    the Ewald separation parameter
+         @param gridSize the dimensions of the mesh
       
          --------------------------------------------------------------------------------------- */
       
-      void setTolerance( RealOpenMM tolerance );
+      void setUsePME(float alpha, int meshSize[3]);
 
       /**---------------------------------------------------------------------------------------
       
-         Apply Shake algorithm
+         Calculate Ewald ixn
       
          @param numberOfAtoms    number of atoms
-         @param atomCoordinates  atom coordinates
-         @param atomCoordinatesP atom coordinates prime
-         @param inverseMasses    1/mass
-      
-         @return SimTKOpenMMCommon::DefaultReturn if converge; else
-          return SimTKOpenMMCommon::ErrorReturn
+         @param posq             atom coordinates and charges
+         @param atomCoordinates  atom coordinates (in format needed by PME)
+         @param atomParameters   atom parameters (sigma/2, 2*sqrt(epsilon))
+         @param exclusions       atom exclusion indices
+                                 exclusions[atomIndex] contains the list of exclusions for that atom
+         @param forces           force array (forces added)
+         @param totalEnergy      total energy
             
          --------------------------------------------------------------------------------------- */
           
-      int apply( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
-                       std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses );
+      void calculateReciprocalIxn(int numberOfAtoms, float* posq, std::vector<OpenMM::RealVec>& atomCoordinates,
+                            const std::vector<std::pair<float, float> >& atomParameters, const std::vector<std::set<int> >& exclusions,
+                            std::vector<OpenMM::RealVec>& forces, float* totalEnergy) const;
       
       /**---------------------------------------------------------------------------------------
       
-         Apply constraint algorithm to velocities.
+         Calculate LJ Coulomb pair ixn
       
          @param numberOfAtoms    number of atoms
-         @param atomCoordinates  atom coordinates
-         @param velocities       atom velocities
-         @param inverseMasses    1/mass
-
-         @return SimTKOpenMMCommon::DefaultReturn if converge; else
-          return SimTKOpenMMCommon::ErrorReturn
+         @param posq             atom coordinates and charges
+         @param atomParameters   atom parameters (sigma/2, 2*sqrt(epsilon))
+         @param exclusions       atom exclusion indices
+                                 exclusions[atomIndex] contains the list of exclusions for that atom
+         @param forces           force array (forces added)
+         @param totalEnergy      total energy
       
          --------------------------------------------------------------------------------------- */
           
-      int applyToVelocities(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
-                     std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses);
+      void calculateDirectIxn(int numberOfAtoms, float* posq, const std::vector<std::pair<float, float> >& atomParameters,
+            const std::vector<std::set<int> >& exclusions, float* forces, float* totalEnergy);
+
+    /**
+     * This routine contains the code executed by each thread.
+     */
+    void runThread(int index, std::vector<float>& threadForce, double& threadEnergy);
+
+private:
+        bool cutoff;
+        bool useSwitch;
+        bool periodic;
+        bool ewald;
+        bool pme;
+        const std::vector<std::pair<int, int> >* neighborList;
+        float periodicBoxSize[3];
+        float cutoffDistance, switchingDistance;
+        float krf, crf;
+        float alphaEwald;
+        int numRx, numRy, numRz;
+        int meshDim[3];
+        std::vector<float> ewaldScaleX, ewaldScaleY, ewaldScaleDeriv;
+        float ewaldDX, ewaldDXInv;
+        bool isDeleted;
+        int numThreads, waitCount;
+        std::vector<pthread_t> thread;
+        std::vector<ThreadData*> threadData;
+        pthread_cond_t startCondition, endCondition;
+        pthread_mutex_t lock;
+        // The following variables are used to make information accessible to the individual threads.
+        int numberOfAtoms;
+        float* posq;
+        std::pair<float, float> const* atomParameters;        
+        std::set<int> const* exclusions;
+        bool includeEnergy;
+
+        static const float TWO_OVER_SQRT_PI;
+        static const int NUM_TABLE_POINTS;
             
       /**---------------------------------------------------------------------------------------
       
-         Report any violated constraints
+         Calculate LJ Coulomb pair ixn between two atoms
       
-         @param numberOfAtoms    number of atoms
-         @param atomCoordinates  atom coordinates
-         @param message          report
+         @param atom1            the index of the first atom
+         @param atom2            the index of the second atom
+         @param forces           force array (forces added)
+         @param totalEnergy      total energy
             
-         @return number of violated constraints
+         --------------------------------------------------------------------------------------- */
+          
+      void calculateOneIxn(int atom1, int atom2, float* forces, double* totalEnergy, const __m128& boxSize, const __m128& invBoxSize);
+            
+      /**---------------------------------------------------------------------------------------
+      
+         Calculate LJ Coulomb pair ixn between two atoms
+      
+         @param atom1            the index of the first atom
+         @param atom2            the index of the second atom
+         @param forces           force array (forces added)
+         @param totalEnergy      total energy
             
          --------------------------------------------------------------------------------------- */
           
-      int reportShake( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::stringstream& message );
+      void calculateOneEwaldIxn(int atom1, int atom2, float* forces, double* totalEnergy, const __m128& boxSize, const __m128& invBoxSize);
+
+      /**
+       * Compute the displacement and squared distance between two points, optionally using
+       * periodic boundary conditions.
+       */
+      void getDeltaR(const __m128& posI, const __m128& posJ, __m128& deltaR, float& r2, bool periodic, const __m128& boxSize, const __m128& invBoxSize) const;
+
+      /**
+       * Compute a fast approximation to erfc(x).
+       */
+      static float erfcApprox(float x);
+
+      /**
+       * Create a lookup table for the scale factor used with Ewald and PME.
+       */
+      void tabulateEwaldScaleFactor();
+      
+      /**
+       * Evaluate the scale factor used with Ewald and PME: erfc(alpha*r) + 2*alpha*r*exp(-alpha*alpha*r*r)/sqrt(PI)
+       */
+      float ewaldScaleFunction(float x);
 };
 
 // ---------------------------------------------------------------------------------------
 
-#endif // __ReferenceShakeAlgorithm_H__
+#endif // OPENMM_CPU_NONBONDED_FORCE_H__

platforms/cpu/include/CpuPlatform.h

+#ifndef OPENMM_CPUPLATFORM_H_
+#define OPENMM_CPUPLATFORM_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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 "ReferencePlatform.h"
+#include "windowsExportCpu.h"
+
+namespace OpenMM {
+    
+/**
+ * This Platform subclass uses CPU implementations of the OpenMM kernels.
+ */
+
+class OPENMM_EXPORT_CPU CpuPlatform : public ReferencePlatform {
+public:
+    CpuPlatform();
+    const std::string& getName() const {
+        static const std::string name = "CPU";
+        return name;
+    }
+    double getSpeed() const;
+    bool supportsDoublePrecision() const;
+    static bool isProcessorSupported();
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CPUPLATFORM_H_*/

platforms/cpu/include/windowsExportCpu.h

+#ifndef OPENMM_WINDOWSEXPORTCPU_H_
+#define OPENMM_WINDOWSEXPORTCPU_H_
+
+/*
+ * Shared libraries are messy in Visual Studio. We have to distinguish three
+ * cases:
+ *   (1) this header is being used to build the OpenMM shared library
+ *       (dllexport)
+ *   (2) this header is being used by a *client* of the OpenMM shared
+ *       library (dllimport)
+ *   (3) we are building the OpenMM static library, or the client is
+ *       being compiled with the expectation of linking with the
+ *       OpenMM static library (nothing special needed)
+ * In the CMake script for building this library, we define one of the symbols
+ *     OPENMM_CPU_BUILDING_{SHARED|STATIC}_LIBRARY
+ * Client code normally has no special symbol defined, in which case we'll
+ * assume it wants to use the shared library. However, if the client defines
+ * the symbol OPENMM_USE_STATIC_LIBRARIES we'll suppress the dllimport so
+ * that the client code can be linked with static libraries. Note that
+ * the client symbol is not library dependent, while the library symbols
+ * affect only the OpenMM library, meaning that other libraries can
+ * be clients of this one. However, we are assuming all-static or all-shared.
+ */
+
+#ifdef _MSC_VER
+    // We don't want to hear about how sprintf is "unsafe".
+    #pragma warning(disable:4996)
+    // Keep MS VC++ quiet about lack of dll export of private members.
+    #pragma warning(disable:4251)
+    #if defined(OPENMM_CPU_BUILDING_SHARED_LIBRARY)
+        #define OPENMM_EXPORT_CPU __declspec(dllexport)
+    #elif defined(OPENMM_CPU_BUILDING_STATIC_LIBRARY) || defined(OPENMM_CPU_USE_STATIC_LIBRARIES)
+        #define OPENMM_EXPORT_CPU
+    #else
+        #define OPENMM_EXPORT_CPU __declspec(dllimport)   // i.e., a client of a shared library
+    #endif
+#else
+    #define OPENMM_EXPORT_CPU // Linux, Mac
+#endif
+
+#endif // OPENMM_WINDOWSEXPORTCPU_H_

platforms/cpu/sharedTarget/CMakeLists.txt

+SET_SOURCE_FILES_PROPERTIES(${SOURCE_FILES} PROPERTIES COMPILE_FLAGS "-msse4.1")
+ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
+
+IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME}_d)
+ELSE (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+    SET(MAIN_OPENMM_LIB ${OPENMM_LIBRARY_NAME})
+ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
+TARGET_LINK_LIBRARIES(${SHARED_TARGET} ${MAIN_OPENMM_LIB} ${PTHREADS_LIB})
+SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_CPU_BUILDING_SHARED_LIBRARY")
+
+INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${SHARED_TARGET})

platforms/cpu/src/CpuKernelFactory.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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 "CpuKernelFactory.h"
+#include "CpuKernels.h"
+#include "CpuPlatform.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/OpenMMException.h"
+
+using namespace OpenMM;
+
+KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
+    ReferencePlatform::PlatformData& data = *static_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    if (name == CalcNonbondedForceKernel::Name())
+        return new CpuCalcNonbondedForceKernel(name, platform);
+    throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str());
+}

platforms/cpu/src/CpuKernels.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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 "CpuKernels.h"
+#include "ReferenceBondForce.h"
+#include "ReferenceLJCoulomb14.h"
+#include "openmm/Context.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/NonbondedForceImpl.h"
+#include "RealVec.h"
+
+using namespace OpenMM;
+using namespace std;
+
+static vector<RealVec>& extractPositions(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *((vector<RealVec>*) data->positions);
+}
+
+static vector<RealVec>& extractVelocities(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *((vector<RealVec>*) data->velocities);
+}
+
+static vector<RealVec>& extractForces(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *((vector<RealVec>*) data->forces);
+}
+
+static RealVec& extractBoxSize(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *(RealVec*) data->periodicBoxSize;
+}
+
+class CpuCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO {
+public:
+    PmeIO(float* posq, float* force, int numParticles) : posq(posq), force(force), numParticles(numParticles) {
+    }
+    float* getPosq() {
+        return posq;
+    }
+    void setForce(float* f) {
+        for (int i = 0; i < numParticles; i++) {
+            force[4*i] += f[4*i];
+            force[4*i+1] += f[4*i+1];
+            force[4*i+2] += f[4*i+2];
+        }
+    }
+private:
+    float* posq;
+    float* force;
+    int numParticles;
+};
+
+CpuCalcNonbondedForceKernel::~CpuCalcNonbondedForceKernel() {
+    if (bonded14ParamArray != NULL) {
+        for (int i = 0; i < num14; i++) {
+            delete[] bonded14IndexArray[i];
+            delete[] bonded14ParamArray[i];
+        }
+        delete bonded14IndexArray;
+        delete bonded14ParamArray;
+    }
+}
+
+void CpuCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
+
+    // Identify which exceptions are 1-4 interactions.
+
+    numParticles = force.getNumParticles();
+    posq.resize(4*numParticles, 0);
+    forces.resize(4*numParticles, 0);
+    exclusions.resize(numParticles);
+    vector<int> nb14s;
+    for (int i = 0; i < force.getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+        if (chargeProd != 0.0 || epsilon != 0.0)
+            nb14s.push_back(i);
+    }
+
+    // Record the particle parameters.
+
+    num14 = nb14s.size();
+    bonded14IndexArray = new int*[num14];
+    for (int i = 0; i < num14; i++)
+        bonded14IndexArray[i] = new int[2];
+    bonded14ParamArray = new double*[num14];
+    for (int i = 0; i < num14; i++)
+        bonded14ParamArray[i] = new double[3];
+    particleParams.resize(numParticles);
+    double sumSquaredCharges = 0.0;
+    for (int i = 0; i < numParticles; ++i) {
+        double charge, radius, depth;
+        force.getParticleParameters(i, charge, radius, depth);
+        posq[4*i+3] = (float) charge;
+        particleParams[i] = make_pair((float) (0.5*radius), (float) (2.0*sqrt(depth)));
+        sumSquaredCharges += charge*charge;
+    }
+    
+    // Recorded exception parameters.
+    
+    for (int i = 0; i < num14; ++i) {
+        int particle1, particle2;
+        double charge, radius, depth;
+        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
+        bonded14IndexArray[i][0] = particle1;
+        bonded14IndexArray[i][1] = particle2;
+        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
+        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
+        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
+    }
+    
+    // Record other parameters.
+    
+    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    nonbondedCutoff = force.getCutoffDistance();
+    if (nonbondedMethod == NoCutoff)
+        useSwitchingFunction = false;
+    else {
+        useSwitchingFunction = force.getUseSwitchingFunction();
+        switchingDistance = force.getSwitchingDistance();
+    }
+    if (nonbondedMethod == Ewald) {
+        double alpha;
+        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
+        ewaldAlpha = alpha;
+    }
+    else if (nonbondedMethod == PME) {
+        double alpha;
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
+        ewaldAlpha = alpha;
+    }
+    if (nonbondedMethod == Ewald || nonbondedMethod == PME)
+        ewaldSelfEnergy = -ONE_4PI_EPS0*ewaldAlpha*sumSquaredCharges/sqrt(M_PI);
+    else
+        ewaldSelfEnergy = 0.0;
+    rfDielectric = force.getReactionFieldDielectric();
+    if (force.getUseDispersionCorrection())
+        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
+    else
+        dispersionCoefficient = 0.0;
+    lastPositions.resize(numParticles, Vec3(1e10, 1e10, 1e10));
+}
+
+double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
+    if (!hasInitializedPme) {
+        hasInitializedPme = true;
+        useOptimizedPme = false;
+        if (nonbondedMethod == PME) {
+            // If available, use the optimized PME implementation.
+
+            try {
+                optimizedPme = getPlatform().createKernel(CalcPmeReciprocalForceKernel::Name(), context);
+                optimizedPme.getAs<CalcPmeReciprocalForceKernel>().initialize(gridSize[0], gridSize[1], gridSize[2], numParticles, ewaldAlpha);
+                useOptimizedPme = true;
+            }
+            catch (OpenMMException& ex) {
+                // The CPU PME plugin isn't available.
+            }
+        }
+    }
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& forceData = extractForces(context);
+    RealVec boxSize = extractBoxSize(context);
+    float floatBoxSize[3] = {(float) boxSize[0], (float) boxSize[1], (float) boxSize[2]};
+    double energy = ewaldSelfEnergy;
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    bool ewald  = (nonbondedMethod == Ewald);
+    bool pme  = (nonbondedMethod == PME);
+    
+    // Convert the positions to single precision.
+    
+    if (periodic)
+        for (int i = 0; i < numParticles; i++)
+            for (int j = 0; j < 3; j++) {
+                RealOpenMM x = posData[i][j];
+                double base = floor(x/boxSize[j]+0.5)*boxSize[j];
+                posq[4*i+j] = (float) (x-base);
+            }
+    else
+        for (int i = 0; i < numParticles; i++) {
+            posq[4*i] = (float) posData[i][0];
+            posq[4*i+1] = (float) posData[i][1];
+            posq[4*i+2] = (float) posData[i][2];
+        }
+    for (int i = 0; i < 4*numParticles; i++)
+        forces[i] = 0.0f;
+    if (nonbondedMethod != NoCutoff) {
+        // Determine whether we need to recompute the neighbor list.
+        
+        double padding = 0.1*nonbondedCutoff;
+        bool needRecompute = false;
+        for (int i = 0; i < numParticles; i++) {
+            RealVec delta = posData[i]-lastPositions[i];
+            if (delta.dot(delta) > 0.25*padding*padding) {
+                needRecompute = true;
+                break;
+            }
+        }
+        if (needRecompute) {
+            neighborList.computeNeighborList(numParticles, posq, exclusions, floatBoxSize, periodic || ewald || pme, nonbondedCutoff+padding);
+            lastPositions = posData;
+        }
+        nonbonded.setUseCutoff(nonbondedCutoff, neighborList.getNeighbors(), rfDielectric);
+    }
+    if (periodic || ewald || pme) {
+        double minAllowedSize = 1.999999*nonbondedCutoff;
+        if (boxSize[0] < minAllowedSize || boxSize[1] < minAllowedSize || boxSize[2] < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+        nonbonded.setPeriodic(floatBoxSize);
+    }
+    if (ewald)
+        nonbonded.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
+    if (pme)
+        nonbonded.setUsePME(ewaldAlpha, gridSize);
+    if (useSwitchingFunction)
+        nonbonded.setUseSwitchingFunction(switchingDistance);
+    float nonbondedEnergy = 0;
+    if (includeDirect)
+        nonbonded.calculateDirectIxn(numParticles, &posq[0], particleParams, exclusions, &forces[0], includeEnergy ? &nonbondedEnergy : NULL);
+    if (includeReciprocal) {
+        if (useOptimizedPme) {
+            PmeIO io(&posq[0], &forces[0], numParticles);
+            Vec3 periodicBoxSize(boxSize[0], boxSize[1], boxSize[2]);
+            optimizedPme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, periodicBoxSize, includeEnergy);
+            optimizedPme.getAs<CalcPmeReciprocalForceKernel>().finishComputation(io);
+        }
+        else
+            nonbonded.calculateReciprocalIxn(numParticles, &posq[0], posData, particleParams, exclusions, forceData, includeEnergy ? &nonbondedEnergy : NULL);
+    }
+    energy += nonbondedEnergy;
+    for (int i = 0; i < numParticles; i++) {
+        forceData[i][0] += forces[4*i];
+        forceData[i][1] += forces[4*i+1];
+        forceData[i][2] += forces[4*i+2];
+    }
+    if (includeDirect) {
+        ReferenceBondForce refBondForce;
+        ReferenceLJCoulomb14 nonbonded14;
+        refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
+        if (periodic || ewald || pme)
+            energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);
+    }
+    return energy;
+}
+
+void CpuCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const NonbondedForce& force) {
+    if (force.getNumParticles() != numParticles)
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    vector<int> nb14s;
+    for (int i = 0; i < force.getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        if (chargeProd != 0.0 || epsilon != 0.0)
+            nb14s.push_back(i);
+    }
+    if (nb14s.size() != num14)
+        throw OpenMMException("updateParametersInContext: The number of non-excluded exceptions has changed");
+
+    // Record the values.
+
+    double sumSquaredCharges = 0.0;
+    for (int i = 0; i < numParticles; ++i) {
+        double charge, radius, depth;
+        force.getParticleParameters(i, charge, radius, depth);
+        posq[4*i+3] = (float) charge;
+        particleParams[i] = make_pair((float) (0.5*radius), (float) (2.0*sqrt(depth)));
+        sumSquaredCharges += charge*charge;
+    }
+    if (nonbondedMethod == Ewald || nonbondedMethod == PME)
+        ewaldSelfEnergy = -ONE_4PI_EPS0*ewaldAlpha*sumSquaredCharges/sqrt(M_PI);
+    else
+        ewaldSelfEnergy = 0.0;
+    for (int i = 0; i < num14; ++i) {
+        int particle1, particle2;
+        double charge, radius, depth;
+        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
+        bonded14IndexArray[i][0] = particle1;
+        bonded14IndexArray[i][1] = particle2;
+        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
+        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
+        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
+    }
+    
+    // Recompute the coefficient for the dispersion correction.
+
+    NonbondedForce::NonbondedMethod method = force.getNonbondedMethod();
+    if (force.getUseDispersionCorrection() && (method == NonbondedForce::CutoffPeriodic || method == NonbondedForce::Ewald || method == NonbondedForce::PME))
+        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
+}

platforms/cpu/src/CpuNeighborList.cpp

+#include "CpuNeighborList.h"
+#include "openmm/internal/hardware.h"
+#include <set>
+#include <map>
+#include <cmath>
+#include <smmintrin.h>
+
+using namespace std;
+
+namespace OpenMM {
+
+class VoxelIndex 
+{
+public:
+    VoxelIndex(int xx, int yy, int zz) : x(xx), y(yy), z(zz) {
+    }
+
+    // operator<() needed for map
+    bool operator<(const VoxelIndex& other) const {
+        if      (x < other.x) return true;
+        else if (x > other.x) return false;
+        else if (y < other.y) return true;
+        else if (y > other.y) return false;
+        else if (z < other.z) return true;
+        else return false;
+    }
+    
+    int x;
+    int y;
+    int z;
+};
+
+typedef pair<const float*, int> VoxelItem;
+typedef vector< VoxelItem > Voxel;
+
+class CpuNeighborList::VoxelHash {
+public:
+    VoxelHash(float vsx, float vsy, float vsz, const float* periodicBoxSize, bool usePeriodic) :
+            voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxSize(periodicBoxSize), usePeriodic(usePeriodic) {
+        if (usePeriodic) {
+            nx = (int) floorf(periodicBoxSize[0]/voxelSizeX+0.5f);
+            ny = (int) floorf(periodicBoxSize[1]/voxelSizeY+0.5f);
+            nz = (int) floorf(periodicBoxSize[2]/voxelSizeZ+0.5f);
+            voxelSizeX = periodicBoxSize[0]/nx;
+            voxelSizeY = periodicBoxSize[1]/ny;
+            voxelSizeZ = periodicBoxSize[2]/nz;
+        }
+    }
+
+    void insert(const int& item, const float* location) {
+        VoxelIndex voxelIndex = getVoxelIndex(location);
+        if (voxelMap.find(voxelIndex) == voxelMap.end())
+            voxelMap[voxelIndex] = Voxel(); 
+        Voxel& voxel = voxelMap.find(voxelIndex)->second;
+        voxel.push_back(VoxelItem(location, item));
+    }
+    
+
+    VoxelIndex getVoxelIndex(const float* location) const {
+        float xperiodic, yperiodic, zperiodic;
+        if (!usePeriodic) {
+            xperiodic = location[0];
+            yperiodic = location[1];
+            zperiodic = location[2];
+        }
+        else {
+            xperiodic = location[0]-periodicBoxSize[0]*floorf(location[0]/periodicBoxSize[0]);
+            yperiodic = location[1]-periodicBoxSize[1]*floorf(location[1]/periodicBoxSize[1]);
+            zperiodic = location[2]-periodicBoxSize[2]*floorf(location[2]/periodicBoxSize[2]);
+        }
+        int x = int(floorf(xperiodic / voxelSizeX));
+        int y = int(floorf(yperiodic / voxelSizeY));
+        int z = int(floorf(zperiodic / voxelSizeZ));
+        
+        return VoxelIndex(x, y, z);
+    }
+
+    void getNeighbors(vector<pair<int, int> >& neighbors, const VoxelItem& referencePoint, const vector<set<int> >& exclusions, float maxDistance) const {
+
+        // Loop over neighboring voxels
+        // TODO use more clever selection of neighboring voxels
+
+        const int atomI = referencePoint.second;
+        const float* locationI = referencePoint.first;
+        __m128 posI = _mm_loadu_ps(locationI);
+        __m128 boxSize = _mm_set_ps(0, periodicBoxSize[2], periodicBoxSize[1], periodicBoxSize[0]);
+        __m128 invBoxSize = _mm_set_ps(0, (1/periodicBoxSize[2]), (1/periodicBoxSize[1]), (1/periodicBoxSize[0]));
+        __m128 half = _mm_set1_ps(0.5);
+        
+        float maxDistanceSquared = maxDistance * maxDistance;
+
+        int dIndexX = int(maxDistance / voxelSizeX) + 1; // How may voxels away do we have to look?
+        int dIndexY = int(maxDistance / voxelSizeY) + 1;
+        int dIndexZ = int(maxDistance / voxelSizeZ) + 1;
+        VoxelIndex centerVoxelIndex = getVoxelIndex(locationI);
+        int lastx = centerVoxelIndex.x+dIndexX;
+        int lasty = centerVoxelIndex.y+dIndexY;
+        int lastz = centerVoxelIndex.z+dIndexZ;
+        if (usePeriodic) {
+            lastx = min(lastx, centerVoxelIndex.x-dIndexX+nx-1);
+            lasty = min(lasty, centerVoxelIndex.y-dIndexY+ny-1);
+            lastz = min(lastz, centerVoxelIndex.z-dIndexZ+nz-1);
+        }
+        VoxelIndex voxelIndex(0, 0, 0);
+        for (int x = centerVoxelIndex.x - dIndexX; x <= lastx; ++x) {
+            voxelIndex.x = x;
+            if (usePeriodic)
+                voxelIndex.x = (x < 0 ? x+nx : (x >= nx ? x-nx : x));
+            for (int y = centerVoxelIndex.y - dIndexY; y <= lasty; ++y) {
+                voxelIndex.y = y;
+                if (usePeriodic)
+                    voxelIndex.y = (y < 0 ? y+ny : (y >= ny ? y-ny : y));
+                for (int z = centerVoxelIndex.z - dIndexZ; z <= lastz; ++z) {
+                    voxelIndex.z = z;
+                    if (usePeriodic)
+                        voxelIndex.z = (z < 0 ? z+nz : (z >= nz ? z-nz : z));
+                    const map<VoxelIndex, Voxel>::const_iterator voxelEntry = voxelMap.find(voxelIndex);
+                    if (voxelEntry == voxelMap.end())
+                        continue; // no such voxel; skip
+                    const Voxel& voxel = voxelEntry->second;
+                    for (Voxel::const_iterator itemIter = voxel.begin(); itemIter != voxel.end(); ++itemIter) {
+                        const int atomJ = itemIter->second;
+
+                        // Avoid duplicate entries.
+                        if (atomJ >= atomI)
+                            break;
+                        
+                        __m128 posJ = _mm_loadu_ps(itemIter->first);
+                        __m128 delta = _mm_sub_ps(posJ, posI);
+                        if (usePeriodic) {
+                            __m128 base = _mm_mul_ps(_mm_floor_ps(_mm_add_ps(_mm_mul_ps(delta, invBoxSize), half)), boxSize);
+                            delta = _mm_sub_ps(delta, base);
+                        }
+                        float dSquared = _mm_cvtss_f32(_mm_dp_ps(delta, delta, 0x71));
+                        if (dSquared > maxDistanceSquared)
+                            continue;
+                        
+                        // Ignore exclusions.
+                        if (exclusions[atomI].find(atomJ) != exclusions[atomI].end())
+                            continue;
+
+                        neighbors.push_back(make_pair(atomI, atomJ));
+                    }
+                }
+            }
+        }
+    }
+
+private:
+    float voxelSizeX, voxelSizeY, voxelSizeZ;
+    int nx, ny, nz;
+    const float* periodicBoxSize;
+    const bool usePeriodic;
+    map<VoxelIndex, Voxel> voxelMap;
+};
+
+class CpuNeighborList::ThreadData {
+public:
+    ThreadData(int index, CpuNeighborList& owner) : index(index), owner(owner) {
+    }
+    int index;
+    CpuNeighborList& owner;
+    vector<pair<int, int> > threadNeighbors;
+};
+
+static void* threadBody(void* args) {
+    CpuNeighborList::ThreadData& data = *reinterpret_cast<CpuNeighborList::ThreadData*>(args);
+    data.owner.runThread(data.index, data.threadNeighbors);
+    delete &data;
+    return 0;
+}
+
+CpuNeighborList::CpuNeighborList() {
+    isDeleted = false;
+    numThreads = getNumProcessors();
+    pthread_cond_init(&startCondition, NULL);
+    pthread_cond_init(&endCondition, NULL);
+    pthread_mutex_init(&lock, NULL);
+    thread.resize(numThreads);
+    pthread_mutex_lock(&lock);
+    waitCount = 0;
+    for (int i = 0; i < numThreads; i++) {
+        ThreadData* data = new ThreadData(i, *this);
+        threadData.push_back(data);
+        pthread_create(&thread[i], NULL, threadBody, data);
+    }
+    while (waitCount < numThreads)
+        pthread_cond_wait(&endCondition, &lock);
+    pthread_mutex_unlock(&lock);
+}
+
+CpuNeighborList::~CpuNeighborList() {
+    isDeleted = true;
+    pthread_mutex_lock(&lock);
+    pthread_cond_broadcast(&startCondition);
+    pthread_mutex_unlock(&lock);
+    for (int i = 0; i < (int) thread.size(); i++)
+        pthread_join(thread[i], NULL);
+    pthread_mutex_destroy(&lock);
+    pthread_cond_destroy(&startCondition);
+    pthread_cond_destroy(&endCondition);
+}
+
+void CpuNeighborList::computeNeighborList(int numAtoms, const vector<float>& atomLocations, const vector<set<int> >& exclusions,
+            const float* periodicBoxSize, bool usePeriodic, float maxDistance) {
+    // Build the voxel hash.
+    
+    float edgeSizeX, edgeSizeY, edgeSizeZ;
+    if (!usePeriodic)
+        edgeSizeX = edgeSizeY = edgeSizeZ = maxDistance; // TODO - adjust this as needed
+    else {
+        edgeSizeX = 0.5f*periodicBoxSize[0]/floorf(periodicBoxSize[0]/maxDistance);
+        edgeSizeY = 0.5f*periodicBoxSize[1]/floorf(periodicBoxSize[1]/maxDistance);
+        edgeSizeZ = 0.5f*periodicBoxSize[2]/floorf(periodicBoxSize[2]/maxDistance);
+    }
+    VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxSize, usePeriodic);
+    for (int i = 0; i < numAtoms; i++)
+        voxelHash.insert(i, &atomLocations[4*i]);
+    
+    // Record the parameters for the threads.
+    
+    this->voxelHash = &voxelHash;
+    this->exclusions = &exclusions;
+    this->atomLocations = &atomLocations[0];
+    this->periodicBoxSize = periodicBoxSize;
+    this->numAtoms = numAtoms;
+    this->usePeriodic = usePeriodic;
+    this->maxDistance = maxDistance;
+    
+    // Signal the threads to start running and wait for them to finish.
+    
+    pthread_mutex_lock(&lock);
+    waitCount = 0;
+    pthread_cond_broadcast(&startCondition);
+    while (waitCount < numThreads)
+        pthread_cond_wait(&endCondition, &lock);
+    pthread_mutex_unlock(&lock);
+    
+    // Combine the results from all the threads.
+    
+    neighbors.clear();
+    for (int i = 0; i < numThreads; i++)
+        neighbors.insert(neighbors.end(), threadData[i]->threadNeighbors.begin(), threadData[i]->threadNeighbors.end());
+}
+
+const vector<pair<int, int> >& CpuNeighborList::getNeighbors() {
+    return neighbors;
+}
+
+void CpuNeighborList::runThread(int index, vector<pair<int, int> >& threadNeighbors) {
+    while (true) {
+        // Wait for the signal to start running.
+        
+        pthread_mutex_lock(&lock);
+        waitCount++;
+        pthread_cond_signal(&endCondition);
+        pthread_cond_wait(&startCondition, &lock);
+        pthread_mutex_unlock(&lock);
+        if (isDeleted)
+            break;
+        
+        // Compute this thread's subset of neighbors.
+        
+        threadNeighbors.clear();
+        for (int i = index; i < numAtoms; i += numThreads)
+            voxelHash->getNeighbors(threadNeighbors, VoxelItem(&atomLocations[4*i], i), *exclusions, maxDistance);
+    }
+}
+
+} // namespace OpenMM

platforms/cpu/src/CpuNonbondedForce.cpp

+
+/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * 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 <complex>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMUtilities.h"
+#include "CpuNonbondedForce.h"
+#include "ReferenceForce.h"
+#include "ReferencePME.h"
+#include "openmm/internal/hardware.h"
+#include "openmm/internal/SplineFitter.h"
+
+// In case we're using some primitive version of Visual Studio this will
+// make sure that erf() and erfc() are defined.
+#include "openmm/internal/MSVC_erfc.h"
+
+using namespace std;
+using namespace OpenMM;
+
+const float CpuNonbondedForce::TWO_OVER_SQRT_PI = (float) (2/sqrt(PI_M));
+const int CpuNonbondedForce::NUM_TABLE_POINTS = 1025;
+
+class CpuNonbondedForce::ThreadData {
+public:
+    ThreadData(int index, CpuNonbondedForce& owner) : index(index), owner(owner) {
+    }
+    int index;
+    CpuNonbondedForce& owner;
+    vector<float> threadForce;
+    double threadEnergy;
+};
+
+static void* threadBody(void* args) {
+    CpuNonbondedForce::ThreadData& data = *reinterpret_cast<CpuNonbondedForce::ThreadData*>(args);
+    data.owner.runThread(data.index, data.threadForce, data.threadEnergy);
+    delete &data;
+    return 0;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   CpuNonbondedForce constructor
+
+   --------------------------------------------------------------------------------------- */
+
+CpuNonbondedForce::CpuNonbondedForce() : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false) {
+    isDeleted = false;
+    numThreads = getNumProcessors();
+    pthread_cond_init(&startCondition, NULL);
+    pthread_cond_init(&endCondition, NULL);
+    pthread_mutex_init(&lock, NULL);
+    thread.resize(numThreads);
+    pthread_mutex_lock(&lock);
+    waitCount = 0;
+    for (int i = 0; i < numThreads; i++) {
+        ThreadData* data = new ThreadData(i, *this);
+        threadData.push_back(data);
+        pthread_create(&thread[i], NULL, threadBody, data);
+    }
+    while (waitCount < numThreads)
+        pthread_cond_wait(&endCondition, &lock);
+    pthread_mutex_unlock(&lock);
+}
+
+/**---------------------------------------------------------------------------------------
+
+   CpuNonbondedForce destructor
+
+   --------------------------------------------------------------------------------------- */
+
+CpuNonbondedForce::~CpuNonbondedForce(){
+    isDeleted = true;
+    pthread_mutex_lock(&lock);
+    pthread_cond_broadcast(&startCondition);
+    pthread_mutex_unlock(&lock);
+    for (int i = 0; i < (int) thread.size(); i++)
+        pthread_join(thread[i], NULL);
+    pthread_mutex_destroy(&lock);
+    pthread_cond_destroy(&startCondition);
+    pthread_cond_destroy(&endCondition);
+}
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use a cutoff.
+
+     @param distance            the cutoff distance
+     @param neighbors           the neighbor list to use
+     @param solventDielectric   the dielectric constant of the bulk solvent
+
+     --------------------------------------------------------------------------------------- */
+
+  void CpuNonbondedForce::setUseCutoff(float distance, const vector<pair<int, int> >& neighbors, float solventDielectric) {
+
+    cutoff = true;
+    cutoffDistance = distance;
+    neighborList = &neighbors;
+    krf = pow(cutoffDistance, -3.0f)*(solventDielectric-1.0)/(2.0*solventDielectric+1.0);
+    crf = (1.0/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
+  }
+
+/**---------------------------------------------------------------------------------------
+
+   Set the force to use a switching function on the Lennard-Jones interaction.
+
+   @param distance            the switching distance
+
+   --------------------------------------------------------------------------------------- */
+
+void CpuNonbondedForce::setUseSwitchingFunction(float distance) {
+    useSwitch = true;
+    switchingDistance = distance;
+}
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use periodic boundary conditions.  This requires that a cutoff has
+     also 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 CpuNonbondedForce::setPeriodic(float* periodicBoxSize) {
+
+    assert(cutoff);
+    assert(periodicBoxSize[0] >= 2*cutoffDistance);
+    assert(periodicBoxSize[1] >= 2*cutoffDistance);
+    assert(periodicBoxSize[2] >= 2*cutoffDistance);
+    periodic = true;
+    this->periodicBoxSize[0] = periodicBoxSize[0];
+    this->periodicBoxSize[1] = periodicBoxSize[1];
+    this->periodicBoxSize[2] = periodicBoxSize[2];
+  }
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use Ewald summation.
+
+     @param alpha  the Ewald separation parameter
+     @param kmaxx  the largest wave vector in the x direction
+     @param kmaxy  the largest wave vector in the y direction
+     @param kmaxz  the largest wave vector in the z direction
+
+     --------------------------------------------------------------------------------------- */
+
+  void CpuNonbondedForce::setUseEwald(float alpha, int kmaxx, int kmaxy, int kmaxz) {
+      alphaEwald = alpha;
+      numRx = kmaxx;
+      numRy = kmaxy;
+      numRz = kmaxz;
+      ewald = true;
+      tabulateEwaldScaleFactor();
+  }
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use Particle-Mesh Ewald (PME) summation.
+
+     @param alpha  the Ewald separation parameter
+     @param gridSize the dimensions of the mesh
+
+     --------------------------------------------------------------------------------------- */
+
+  void CpuNonbondedForce::setUsePME(float alpha, int meshSize[3]) {
+      alphaEwald = alpha;
+      meshDim[0] = meshSize[0];
+      meshDim[1] = meshSize[1];
+      meshDim[2] = meshSize[2];
+      pme = true;
+      tabulateEwaldScaleFactor();
+  }
+
+  
+void CpuNonbondedForce::tabulateEwaldScaleFactor() {
+    ewaldDX = cutoffDistance/(NUM_TABLE_POINTS-2);
+    ewaldDXInv = 1.0f/ewaldDX;
+    vector<double> x(NUM_TABLE_POINTS);
+    vector<double> y(NUM_TABLE_POINTS);
+    vector<double> deriv;
+    for (int i = 0; i < NUM_TABLE_POINTS; i++) {
+        double r = i*cutoffDistance/(NUM_TABLE_POINTS-2);
+        double alphaR = alphaEwald*r;
+        x[i] = r;
+        y[i] = erfc(alphaR) + TWO_OVER_SQRT_PI*alphaR*exp(-alphaR*alphaR);
+    }
+    SplineFitter::createNaturalSpline(x, y, deriv);
+    ewaldScaleX.resize(NUM_TABLE_POINTS);
+    ewaldScaleY.resize(NUM_TABLE_POINTS);
+    ewaldScaleDeriv.resize(NUM_TABLE_POINTS);
+    for (int i = 0; i < NUM_TABLE_POINTS; i++) {
+        ewaldScaleX[i] = (float) x[i];
+        ewaldScaleY[i] = (float) y[i];
+        ewaldScaleDeriv[i] = (float) (deriv[i]*ewaldDX*ewaldDX/6);
+    }
+}
+  
+void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, vector<RealVec>& atomCoordinates,
+                                             const vector<pair<float, float> >& atomParameters, const vector<set<int> >& exclusions,
+                                             vector<RealVec>& forces, float* totalEnergy) const {
+    typedef std::complex<float> d_complex;
+
+    static const float epsilon     =  1.0;
+
+    int kmax                            = (ewald ? std::max(numRx, std::max(numRy,numRz)) : 0);
+    float factorEwald              = -1 / (4*alphaEwald*alphaEwald);
+    float TWO_PI                   = 2.0 * PI_M;
+    float recipCoeff               = (float)(ONE_4PI_EPS0*4*PI_M/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
+
+    if (pme) {
+        pme_t pmedata;
+        RealOpenMM virial[3][3];
+        pme_init(&pmedata, alphaEwald, numberOfAtoms, meshDim, 5, 1);
+        vector<RealOpenMM> charges(numberOfAtoms);
+        for (int i = 0; i < numberOfAtoms; i++)
+            charges[i] = posq[4*i+3];
+        RealOpenMM boxSize[3] = {periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2]};
+        RealOpenMM recipEnergy = 0.0;
+        pme_exec(pmedata, atomCoordinates, forces, charges, boxSize, &recipEnergy, virial);
+        if (totalEnergy)
+            *totalEnergy += recipEnergy;
+        pme_destroy(pmedata);
+    }
+
+    // Ewald method
+
+    else if (ewald) {
+
+        // setup reciprocal box
+
+        float recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
+
+
+        // setup K-vectors
+
+        #define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
+        vector<d_complex> eir(kmax*numberOfAtoms*3);
+        vector<d_complex> tab_xy(numberOfAtoms);
+        vector<d_complex> tab_qxyz(numberOfAtoms);
+
+        for (int i = 0; (i < numberOfAtoms); i++) {
+            float* pos = posq+4*i;
+            for (int m = 0; (m < 3); m++)
+              EIR(0, i, m) = d_complex(1,0);
+
+            for (int m=0; (m<3); m++)
+              EIR(1, i, m) = d_complex(cos(pos[m]*recipBoxSize[m]),
+                                       sin(pos[m]*recipBoxSize[m]));
+
+            for (int j=2; (j<kmax); j++)
+              for (int m=0; (m<3); m++)
+                EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
+        }
+
+        // calculate reciprocal space energy and forces
+
+        int lowry = 0;
+        int lowrz = 1;
+
+        for (int rx = 0; rx < numRx; rx++) {
+            float kx = rx * recipBoxSize[0];
+            for (int ry = lowry; ry < numRy; ry++) {
+                float ky = ry * recipBoxSize[1];
+                if (ry >= 0) {
+                    for (int n = 0; n < numberOfAtoms; n++)
+                      tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
+                }
+                else {
+                    for (int n = 0; n < numberOfAtoms; n++)
+                      tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
+                }
+                for (int rz = lowrz; rz < numRz; rz++) {
+                    if (rz >= 0) {
+                        for (int n = 0; n < numberOfAtoms; n++)
+                            tab_qxyz[n] = posq[4*n+3] * (tab_xy[n] * EIR(rz, n, 2));
+                    }
+                    else {
+                        for (int n = 0; n < numberOfAtoms; n++)
+                            tab_qxyz[n] = posq[4*n+3] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
+                    }
+                    float cs = 0.0f;
+                    float ss = 0.0f;
+
+                    for (int n = 0; n < numberOfAtoms; n++) {
+                        cs += tab_qxyz[n].real();
+                        ss += tab_qxyz[n].imag();
+                    }
+
+                    float kz = rz * recipBoxSize[2];
+                    float k2 = kx * kx + ky * ky + kz * kz;
+                    float ak = exp(k2*factorEwald) / k2;
+
+                    for (int n = 0; n < numberOfAtoms; n++) {
+                        float force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
+                        forces[n][0] += 2 * recipCoeff * force * kx;
+                        forces[n][1] += 2 * recipCoeff * force * ky;
+                        forces[n][2] += 2 * recipCoeff * force * kz;
+                    }
+
+                    if (totalEnergy)
+                        *totalEnergy += recipCoeff * ak * (cs * cs + ss * ss);
+
+                    lowrz = 1 - numRz;
+                }
+                lowry = 1 - numRy;
+            }
+        }
+    }
+}
+
+
+void CpuNonbondedForce::calculateDirectIxn(int numberOfAtoms, float* posq, const vector<pair<float, float> >& atomParameters,
+                const vector<set<int> >& exclusions, float* forces, float* totalEnergy) {
+    // Record the parameters for the threads.
+    
+    this->numberOfAtoms = numberOfAtoms;
+    this->posq = posq;
+    this->atomParameters = &atomParameters[0];
+    this->exclusions = &exclusions[0];
+    includeEnergy = (totalEnergy != NULL);
+    
+    // Signal the threads to start running and wait for them to finish.
+    
+    pthread_mutex_lock(&lock);
+    waitCount = 0;
+    pthread_cond_broadcast(&startCondition);
+    while (waitCount < numThreads)
+        pthread_cond_wait(&endCondition, &lock);
+    pthread_mutex_unlock(&lock);
+    
+    // Combine the results from all the threads.
+    
+    double directEnergy = 0;
+    for (int i = 0; i < numThreads; i++)
+        directEnergy += threadData[i]->threadEnergy;
+    for (int i = 0; i < numberOfAtoms; i++) {
+        __m128 f = _mm_loadu_ps(forces+4*i);
+        for (int j = 0; j < numThreads; j++)
+            f = _mm_add_ps(f, _mm_loadu_ps(&threadData[j]->threadForce[4*i]));
+        _mm_storeu_ps(forces+4*i, f);
+    }
+
+    if (ewald || pme) {
+        // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
+
+        __m128 boxSize = _mm_set_ps(0, periodicBoxSize[2], periodicBoxSize[1], periodicBoxSize[0]);
+        __m128 invBoxSize = _mm_set_ps(0, (1/periodicBoxSize[2]), (1/periodicBoxSize[1]), (1/periodicBoxSize[0]));
+        for (int i = 0; i < numberOfAtoms; i++)
+            for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
+                if (*iter > i) {
+                    int ii = i;
+                    int jj = *iter;
+                    __m128 deltaR;
+                    __m128 posI = _mm_loadu_ps(posq+4*ii);
+                    __m128 posJ = _mm_loadu_ps(posq+4*jj);
+                    float r2;
+                    getDeltaR(posJ, posI, deltaR, r2, false, boxSize, invBoxSize);
+                    float r         = sqrtf(r2);
+                    float inverseR  = 1/r;
+                    float chargeProd = ONE_4PI_EPS0*posq[4*ii+3]*posq[4*jj+3];
+                    float alphaR = alphaEwald*r;
+                    float erfcAlphaR = erfcApprox(alphaR);
+                    float dEdR      = (float) (chargeProd * inverseR * inverseR * inverseR);
+                          dEdR      = (float) (dEdR * (1.0f-erfcAlphaR-TWO_OVER_SQRT_PI*alphaR*exp(-alphaR*alphaR)));
+                    __m128 result = _mm_mul_ps(deltaR, _mm_set1_ps(dEdR));
+                    _mm_storeu_ps(forces+4*ii, _mm_sub_ps(_mm_loadu_ps(forces+4*ii), result));
+                    _mm_storeu_ps(forces+4*jj, _mm_add_ps(_mm_loadu_ps(forces+4*jj), result));
+                    if (includeEnergy)
+                        directEnergy -= chargeProd*inverseR*(1.0f-erfcAlphaR);
+                }
+            }
+    }
+    if (totalEnergy != NULL)
+        *totalEnergy += (float) directEnergy;
+}
+
+
+void CpuNonbondedForce::runThread(int index, vector<float>& threadForce, double& threadEnergy) {
+    while (true) {
+        // Wait for the signal to start running.
+        
+        pthread_mutex_lock(&lock);
+        waitCount++;
+        pthread_cond_signal(&endCondition);
+        pthread_cond_wait(&startCondition, &lock);
+        pthread_mutex_unlock(&lock);
+        if (isDeleted)
+            break;
+        
+        // Compute this thread's subset of interactions.
+        
+        threadEnergy = 0;
+        double* energyPtr = (includeEnergy ? &threadEnergy : NULL);
+        threadForce.resize(4*numberOfAtoms, 0.0f);
+        for (int i = 0; i < 4*numberOfAtoms; i++)
+            threadForce[i] = 0.0f;
+        __m128 boxSize = _mm_set_ps(0, periodicBoxSize[2], periodicBoxSize[1], periodicBoxSize[0]);
+        __m128 invBoxSize = _mm_set_ps(0, (1/periodicBoxSize[2]), (1/periodicBoxSize[1]), (1/periodicBoxSize[0]));
+        if (ewald || pme) {
+            // Compute the interactions from the neighbor list.
+
+            for (int i = index; i < (int) neighborList->size(); i += numThreads) {
+                pair<int, int> pair = (*neighborList)[i];
+                calculateOneEwaldIxn(pair.first, pair.second, &threadForce[0], energyPtr, boxSize, invBoxSize);
+            }
+        }
+        else if (cutoff) {
+            // Compute the interactions from the neighbor list.
+
+            for (int i = index; i < (int) neighborList->size(); i += numThreads) {
+                pair<int, int> pair = (*neighborList)[i];
+                calculateOneIxn(pair.first, pair.second, &threadForce[0], energyPtr, boxSize, invBoxSize);
+            }
+        }
+        else {
+            // Loop over all atom pairs
+
+            for (int i = index; i < numberOfAtoms; i += numThreads){
+                for (int j = i+1; j < numberOfAtoms; j++)
+                    if (exclusions[j].find(i) == exclusions[j].end())
+                        calculateOneIxn(i, j, &threadForce[0], energyPtr, boxSize, invBoxSize);
+            }
+        }
+    }
+}
+
+void CpuNonbondedForce::calculateOneIxn(int ii, int jj, float* forces, double* totalEnergy, const __m128& boxSize, const __m128& invBoxSize) {
+    // get deltaR, R2, and R between 2 atoms
+
+    __m128 deltaR;
+    __m128 posI = _mm_loadu_ps(posq+4*ii);
+    __m128 posJ = _mm_loadu_ps(posq+4*jj);
+    float r2;
+    getDeltaR(posJ, posI, deltaR, r2, periodic, boxSize, invBoxSize);
+    if (cutoff && r2 >= cutoffDistance*cutoffDistance)
+        return;
+    float r = sqrtf(r2);
+    float inverseR = 1/r;
+    float switchValue = 1, switchDeriv = 0;
+    if (useSwitch && r > switchingDistance) {
+        float t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
+        switchValue = 1+t*t*t*(-10+t*(15-t*6));
+        switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
+    }
+    float sig       = atomParameters[ii].first + atomParameters[jj].first;
+    float sig2      = inverseR*sig;
+          sig2     *= sig2;
+    float sig6      = sig2*sig2*sig2;
+
+    float eps       = atomParameters[ii].second*atomParameters[jj].second;
+    float dEdR      = switchValue*eps*(12.0f*sig6 - 6.0f)*sig6;
+    float chargeProd = ONE_4PI_EPS0*posq[4*ii+3]*posq[4*jj+3];
+    if (cutoff)
+        dEdR += (float) (chargeProd*(inverseR-2.0f*krf*r2));
+    else
+        dEdR += (float) (chargeProd*inverseR);
+    dEdR *= inverseR*inverseR;
+    float energy = eps*(sig6-1.0f)*sig6;
+    if (useSwitch) {
+        dEdR -= energy*switchDeriv*inverseR;
+        energy *= switchValue;
+    }
+
+    // accumulate energies
+
+    if (totalEnergy) {
+        if (cutoff)
+            energy += (float) (chargeProd*(inverseR+krf*r2-crf));
+        else
+            energy += (float) (chargeProd*inverseR);
+        *totalEnergy += energy;
+    }
+
+    // accumulate forces
+
+    __m128 result = _mm_mul_ps(deltaR, _mm_set1_ps(dEdR));
+    _mm_storeu_ps(forces+4*ii, _mm_add_ps(_mm_loadu_ps(forces+4*ii), result));
+    _mm_storeu_ps(forces+4*jj, _mm_sub_ps(_mm_loadu_ps(forces+4*jj), result));
+  }
+
+void CpuNonbondedForce::calculateOneEwaldIxn(int ii, int jj, float* forces, double* totalEnergy, const __m128& boxSize, const __m128& invBoxSize) {
+    __m128 deltaR;
+    __m128 posI = _mm_loadu_ps(posq+4*ii);
+    __m128 posJ = _mm_loadu_ps(posq+4*jj);
+    float r2;
+    getDeltaR(posJ, posI, deltaR, r2, true, boxSize, invBoxSize);
+    if (r2 < cutoffDistance*cutoffDistance) {
+        float r         = sqrtf(r2);
+        float inverseR  = 1/r;
+        float switchValue = 1, switchDeriv = 0;
+        if (useSwitch && r > switchingDistance) {
+            float t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
+            switchValue = 1+t*t*t*(-10+t*(15-t*6));
+            switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
+        }
+        float chargeProd = ONE_4PI_EPS0*posq[4*ii+3]*posq[4*jj+3];
+        float dEdR      = chargeProd*inverseR*ewaldScaleFunction(r);
+        float sig       = atomParameters[ii].first +  atomParameters[jj].first;
+        float sig2      = inverseR*sig;
+              sig2     *= sig2;
+        float sig6      = sig2*sig2*sig2;
+        float eps       = atomParameters[ii].second*atomParameters[jj].second;
+              dEdR     += switchValue*eps*(12.0f*sig6 - 6.0f)*sig6;
+        dEdR *= inverseR*inverseR;
+        float energy = eps*(sig6-1.0f)*sig6;
+        if (useSwitch) {
+            dEdR -= energy*switchDeriv*inverseR;
+            energy *= switchValue;
+        }
+
+        // accumulate forces
+
+        __m128 result = _mm_mul_ps(deltaR, _mm_set1_ps(dEdR));
+        _mm_storeu_ps(forces+4*ii, _mm_add_ps(_mm_loadu_ps(forces+4*ii), result));
+        _mm_storeu_ps(forces+4*jj, _mm_sub_ps(_mm_loadu_ps(forces+4*jj), result));
+
+        // accumulate energies
+
+        if (totalEnergy) {
+            energy += (float) (chargeProd*inverseR*erfcApprox(alphaEwald*r));
+            *totalEnergy += energy;
+        }
+    }
+}
+
+void CpuNonbondedForce::getDeltaR(const __m128& posI, const __m128& posJ, __m128& deltaR, float& r2, bool periodic, const __m128& boxSize, const __m128& invBoxSize) const {
+    deltaR = _mm_sub_ps(posJ, posI);
+    if (periodic) {
+        __m128 base = _mm_mul_ps(_mm_floor_ps(_mm_add_ps(_mm_mul_ps(deltaR, invBoxSize), _mm_set1_ps(0.5))), boxSize);
+        deltaR = _mm_sub_ps(deltaR, base);
+    }
+    r2 = _mm_cvtss_f32(_mm_dp_ps(deltaR, deltaR, 0x71));
+}
+
+float CpuNonbondedForce::erfcApprox(float x) {
+    // This approximation for erfc is from Abramowitz and Stegun (1964) p. 299.  They cite the following as
+    // the original source: C. Hastings, Jr., Approximations for Digital Computers (1955).  It has a maximum
+    // error of 3e-7.
+
+    float t = 1.0f+(0.0705230784f+(0.0422820123f+(0.0092705272f+(0.0001520143f+(0.0002765672f+0.0000430638f*x)*x)*x)*x)*x)*x;
+    t *= t;
+    t *= t;
+    t *= t;
+    return 1.0f/(t*t);
+}
+
+float CpuNonbondedForce::ewaldScaleFunction(float x) {
+    // Compute the tabulated Ewald scale factor: erfc(alpha*r) + 2*alpha*r*exp(-alpha*alpha*r*r)/sqrt(PI)
+
+    int lower = (int) (x*ewaldDXInv);
+    int upper = lower+1;
+    float a = (ewaldScaleX[upper]-x)*ewaldDXInv;
+    float b = 1.0f-a;
+    return a*ewaldScaleY[lower]+b*ewaldScaleY[upper]+((a*a*a-a)*ewaldScaleDeriv[lower] + (b*b*b-b)*ewaldScaleDeriv[upper]);
+}

platforms/cpu/src/CpuPlatform.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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 "CpuPlatform.h"
+#include "CpuKernelFactory.h"
+#include "CpuKernels.h"
+#include "openmm/internal/hardware.h"
+
+using namespace OpenMM;
+
+extern "C" OPENMM_EXPORT_CPU void registerPlatforms() {
+    // Only register this platform if the CPU supports SSE 4.1.
+
+    if (CpuPlatform::isProcessorSupported())
+        Platform::registerPlatform(new CpuPlatform());
+}
+
+CpuPlatform::CpuPlatform() {
+    CpuKernelFactory* factory = new CpuKernelFactory();
+    registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
+}
+
+double CpuPlatform::getSpeed() const {
+    return 10;
+}
+
+bool CpuPlatform::supportsDoublePrecision() const {
+    return false;
+}
+
+bool CpuPlatform::isProcessorSupported() {
+    // Make sure the CPU supports SSE 4.1.
+    
+    int cpuInfo[4];
+    cpuid(cpuInfo, 0);
+    if (cpuInfo[0] >= 1) {
+        cpuid(cpuInfo, 1);
+        return ((cpuInfo[2] & ((int) 1 << 19)) != 0);
+    }
+    return false;
+}

platforms/cpu/tests/CMakeLists.txt

+#
+# Testing
+#
+
+ENABLE_TESTING()
+
+SET( INCLUDE_SERIALIZATION FALSE )
+#SET( INCLUDE_SERIALIZATION TRUE )
+
+IF( INCLUDE_SERIALIZATION )
+    INCLUDE_DIRECTORIES(${OPENMM_DIR}/serialization/include)
+    SET( SHARED_OPENMM_SERIALIZATION "OpenMMSerialization" )
+ENDIF( INCLUDE_SERIALIZATION )
+
+# Automatically create tests using files named "Test*.cpp"
+FILE(GLOB TEST_PROGS "*Test*.cpp")
+FOREACH(TEST_PROG ${TEST_PROGS})
+    GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
+
+    # Link with shared library
+    ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
+    TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET})
+    ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT} single)
+
+ENDFOREACH(TEST_PROG ${TEST_PROGS})

platforms/cpu/tests/TestCpuEwald.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Ewald summation method CPU implementation of NonbondedForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CpuPlatform.h"
+#include "ReferencePlatform.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/internal/ContextImpl.h"
+#include "SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+CpuPlatform platform;
+
+const double TOL = 1e-5;
+
+void testEwaldPME(bool includeExceptions) {
+
+//      Use amorphous NaCl system for the tests
+
+    const int numParticles = 894;
+    const double cutoff = 1.2;
+    const double boxSize = 3.00646;
+    double tol = 1e-5;
+
+    ReferencePlatform reference;
+    System system;
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setEwaldErrorTolerance(tol);
+
+    for (int i = 0; i < numParticles/2; i++)
+        system.addParticle(22.99);
+    for (int i = 0; i < numParticles/2; i++)
+        system.addParticle(35.45);
+    for (int i = 0; i < numParticles/2; i++)
+        nonbonded->addParticle(1.0, 1.0,0.0);
+    for (int i = 0; i < numParticles/2; i++)
+        nonbonded->addParticle(-1.0, 1.0,0.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+
+    vector<Vec3> positions(numParticles);
+    #include "nacl_amorph.dat"
+    if (includeExceptions) {
+        // Add some exclusions.
+
+        for (int i = 0; i < numParticles-1; i++) {
+            Vec3 delta = positions[i]-positions[i+1];
+            if (sqrt(delta.dot(delta)) < 0.5*cutoff)
+                nonbonded->addException(i, i+1, i%2 == 0 ? 0.0 : 0.5, 1.0, 0.0);
+        }
+    }
+
+//    (1)  Check whether the Reference and CPU platforms agree when using Ewald Method
+
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context cpuContext(system, integrator1, platform);
+    Context referenceContext(system, integrator2, reference);
+    cpuContext.setPositions(positions);
+    referenceContext.setPositions(positions);
+    State cpuState = cpuContext.getState(State::Forces | State::Energy);
+    State referenceState = referenceContext.getState(State::Forces | State::Energy);
+    tol = 1e-2;
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(referenceState.getForces()[i], cpuState.getForces()[i], tol);
+    }
+    tol = 1e-5;
+    ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cpuState.getPotentialEnergy(), tol);
+
+//    (2) Check whether Ewald method in CPU is self-consistent
+
+    double norm = 0.0;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 f = cpuState.getForces()[i];
+        norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+    }
+
+    norm = std::sqrt(norm);
+    const double delta = 5e-3;
+    double step = delta/norm;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = cpuState.getForces()[i];
+        positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+    }
+    VerletIntegrator integrator3(0.01);
+    Context cpuContext2(system, integrator3, platform);
+    cpuContext2.setPositions(positions);
+
+    tol = 1e-2;
+    State cpuState2 = cpuContext2.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (cpuState2.getPotentialEnergy()-cpuState.getPotentialEnergy())/delta, tol)
+
+//    (3)  Check whether the Reference and CPU platforms agree when using PME
+
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    cpuContext.reinitialize();
+    referenceContext.reinitialize();
+    cpuContext.setPositions(positions);
+    referenceContext.setPositions(positions);
+    cpuState = cpuContext.getState(State::Forces | State::Energy);
+    referenceState = referenceContext.getState(State::Forces | State::Energy);
+    tol = 1e-2;
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(referenceState.getForces()[i], cpuState.getForces()[i], tol);
+    }
+    tol = 1e-5;
+    ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cpuState.getPotentialEnergy(), tol);
+
+//    (4) Check whether PME method in CPU is self-consistent
+
+    norm = 0.0;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 f = cpuState.getForces()[i];
+        norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+    }
+
+    norm = std::sqrt(norm);
+    step = delta/norm;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = cpuState.getForces()[i];
+        positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+    }
+    VerletIntegrator integrator4(0.01);
+    Context cpuContext3(system, integrator4, platform);
+    cpuContext3.setPositions(positions);
+
+    tol = 1e-2;
+    State cpuState3 = cpuContext3.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (cpuState3.getPotentialEnergy()-cpuState.getPotentialEnergy())/delta, tol)
+}
+
+void testEwald2Ions() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(-1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setEwaldErrorTolerance(TOL);
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(3.048000,2.764000,3.156000);
+    positions[1] = Vec3(2.809000,2.888000,2.571000);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+
+    ASSERT_EQUAL_VEC(Vec3(-123.711,  64.1877, -302.716), forces[0], 10*TOL);
+    ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877,  302.716), forces[1], 10*TOL);
+    ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
+}
+
+void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 5.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
+        positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+    }
+    force->setNonbondedMethod(method);
+
+    // For various values of the cutoff and error tolerance, see if the actual error is reasonable.
+
+    for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff *= 1.2) {
+        force->setCutoffDistance(cutoff);
+        vector<Vec3> refForces;
+        double norm = 0.0;
+        for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
+            force->setEwaldErrorTolerance(tol);
+            VerletIntegrator integrator(0.01);
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces);
+            if (refForces.size() == 0) {
+                refForces = state.getForces();
+                for (int i = 0; i < numParticles; i++)
+                    norm += refForces[i].dot(refForces[i]);
+                norm = sqrt(norm);
+            }
+            else {
+                double diff = 0.0;
+                for (int i = 0; i < numParticles; i++) {
+                    Vec3 delta = refForces[i]-state.getForces()[i];
+                    diff += delta.dot(delta);
+                }
+                diff = sqrt(diff)/norm;
+                ASSERT(diff < 2*tol);
+            }
+        }
+    }
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        if (!CpuPlatform::isProcessorSupported()) {
+            cout << "CPU is not supported.  Exiting." << endl;
+            return 0;
+        }
+        testEwaldPME(false);
+        testEwaldPME(true);
+//        testEwald2Ions();
+        testErrorTolerance(NonbondedForce::Ewald);
+        testErrorTolerance(NonbondedForce::PME);
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cpu/tests/TestCpuNeighborList.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests all the CPU implementation of neighbor list construction.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "CpuNeighborList.h"
+#include "CpuPlatform.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <set>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testNeighborList(bool periodic) {
+    const int numParticles = 500;
+    const float cutoff = 2.0f;
+    const float boxSize[3] = {20.0f, 15.0f, 22.0f};
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<float> positions(4*numParticles);
+    for (int i = 0; i < 4*numParticles; i++)
+        if (i%4 < 3)
+            positions[i] = boxSize[i%4]*genrand_real2(sfmt);
+    vector<set<int> > exclusions(numParticles);
+    for (int i = 0; i < numParticles; i++) {
+        int num = min(i+1, 10);
+        for (int j = 0; j < num; j++) {
+            exclusions[i].insert(i-j);
+            exclusions[i-j].insert(i);
+        }
+    }
+    CpuNeighborList neighborList;
+    neighborList.computeNeighborList(numParticles, positions, exclusions, boxSize, periodic, cutoff);
+    
+    // Convert the neighbor list to a set for faster lookup.
+    
+    set<pair<int, int> > neighbors;
+    for (int i = 0; i < (int) neighborList.getNeighbors().size(); i++) {
+        pair<int, int> entry = neighborList.getNeighbors()[i];
+        ASSERT(neighbors.find(entry) == neighbors.end() && neighbors.find(make_pair(entry.second, entry.first)) == neighbors.end()); // No duplicates
+        neighbors.insert(entry);
+    }
+    
+    // Check each particle pair and figure out whether they should be in the neighbor list.
+    
+    for (int i = 0; i < numParticles; i++)
+        for (int j = 0; j <= i; j++) {
+            bool shouldInclude = (exclusions[i].find(j) == exclusions[i].end());
+            float dx = positions[4*i]-positions[4*j];
+            float dy = positions[4*i+1]-positions[4*j+1];
+            float dz = positions[4*i+2]-positions[4*j+2];
+            if (periodic) {
+                dx -= floor(dx/boxSize[0]+0.5f)*boxSize[0];
+                dy -= floor(dy/boxSize[1]+0.5f)*boxSize[1];
+                dz -= floor(dz/boxSize[2]+0.5f)*boxSize[2];
+            }
+            if (dx*dx + dy*dy + dz*dz > cutoff*cutoff)
+                shouldInclude = false;
+            bool isIncluded = (neighbors.find(make_pair(i, j)) != neighbors.end() || neighbors.find(make_pair(j, i)) != neighbors.end());
+            ASSERT_EQUAL(shouldInclude, isIncluded);
+        }
+}
+
+int main() {
+    try {
+        if (!CpuPlatform::isProcessorSupported()) {
+            cout << "CPU is not supported.  Exiting." << endl;
+            return 0;
+        }
+        testNeighborList(false);
+        testNeighborList(true);
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cpu/tests/TestCpuNonbondedForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests all the different force terms in the CUDA implementation of NonbondedForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CpuPlatform.h"
+#include "ReferencePlatform.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/internal/ContextImpl.h"
+#include "SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+CpuPlatform platform;
+
+const double TOL = 1e-5;
+
+void testCoulomb() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->addParticle(0.5, 1, 0);
+    forceField->addParticle(-1.5, 1, 0);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double force = ONE_4PI_EPS0*(-0.75)/4.0;
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
+}
+
+void testLJ() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->addParticle(0, 1.2, 1);
+    forceField->addParticle(0, 1.4, 2);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double x = 1.3/2.0;
+    double eps = SQRT_TWO;
+    double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/2.0;
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
+}
+
+void testExclusionsAnd14() {
+    System system;
+    NonbondedForce* nonbonded = new NonbondedForce();
+    for (int i = 0; i < 5; ++i) {
+        system.addParticle(1.0);
+        nonbonded->addParticle(0, 1.5, 0);
+    }
+    vector<pair<int, int> > bonds;
+    bonds.push_back(pair<int, int>(0, 1));
+    bonds.push_back(pair<int, int>(1, 2));
+    bonds.push_back(pair<int, int>(2, 3));
+    bonds.push_back(pair<int, int>(3, 4));
+    nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
+    int first14, second14;
+    for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
+            first14 = i;
+        if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
+            second14 = i;
+    }
+    system.addForce(nonbonded);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    for (int i = 1; i < 5; ++i) {
+
+        // Test LJ forces
+
+        vector<Vec3> positions(5);
+        const double r = 1.0;
+        for (int j = 0; j < 5; ++j) {
+            nonbonded->setParticleParameters(j, 0, 1.5, 0);
+            positions[j] = Vec3(0, j, 0);
+        }
+        nonbonded->setParticleParameters(0, 0, 1.5, 1);
+        nonbonded->setParticleParameters(i, 0, 1.5, 1);
+        nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
+        positions[i] = Vec3(r, 0, 0);
+        context.reinitialize();
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double x = 1.5/r;
+        double eps = 1.0;
+        double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
+        double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        if (i == 3) {
+            force *= 0.5;
+            energy *= 0.5;
+        }
+        if (i < 3) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+
+        // Test Coulomb forces
+
+        nonbonded->setParticleParameters(0, 2, 1.5, 0);
+        nonbonded->setParticleParameters(i, 2, 1.5, 0);
+        nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
+        context.reinitialize();
+        context.setPositions(positions);
+        state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces2 = state.getForces();
+        force = ONE_4PI_EPS0*4/(r*r);
+        energy = ONE_4PI_EPS0*4/r;
+        if (i == 3) {
+            force /= 1.2;
+            energy /= 1.2;
+        }
+        if (i < 3) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+    }
+}
+
+void testCutoff() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->addParticle(1.0, 1, 0);
+    forceField->addParticle(1.0, 1, 0);
+    forceField->addParticle(1.0, 1, 0);
+    forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    const double cutoff = 2.9;
+    forceField->setCutoffDistance(cutoff);
+    const double eps = 50.0;
+    forceField->setReactionFieldDielectric(eps);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(0, 3, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    const double force1 = ONE_4PI_EPS0*(1.0)*(0.25-2.0*krf*2.0);
+    const double force2 = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
+    ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
+    const double energy1 = ONE_4PI_EPS0*(1.0)*(0.5+krf*4.0-crf);
+    const double energy2 = ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf);
+    ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
+}
+
+void testCutoff14() {
+    System system;
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    for (int i = 0; i < 5; ++i) {
+        system.addParticle(1.0);
+        nonbonded->addParticle(0, 1.5, 0);
+    }
+    const double cutoff = 3.5;
+    nonbonded->setCutoffDistance(cutoff);
+    const double eps = 30.0;
+    nonbonded->setReactionFieldDielectric(eps);
+    vector<pair<int, int> > bonds;
+    bonds.push_back(pair<int, int>(0, 1));
+    bonds.push_back(pair<int, int>(1, 2));
+    bonds.push_back(pair<int, int>(2, 3));
+    bonds.push_back(pair<int, int>(3, 4));
+    nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
+    int first14, second14;
+    for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
+            first14 = i;
+        if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
+            second14 = i;
+    }
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(2, 0, 0);
+    positions[3] = Vec3(3, 0, 0);
+    positions[4] = Vec3(4, 0, 0);
+    for (int i = 1; i < 5; ++i) {
+
+        // Test LJ forces
+
+        nonbonded->setParticleParameters(0, 0, 1.5, 1);
+        for (int j = 1; j < 5; ++j)
+            nonbonded->setParticleParameters(j, 0, 1.5, 0);
+        nonbonded->setParticleParameters(i, 0, 1.5, 1);
+        nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
+        context.reinitialize();
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double r = positions[i][0];
+        double x = 1.5/r;
+        double e = 1.0;
+        double force = 4.0*e*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
+        double energy = 4.0*e*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        if (i == 3) {
+            force *= 0.5;
+            energy *= 0.5;
+        }
+        if (i < 3 || r > cutoff) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+
+        // Test Coulomb forces
+
+        const double q = 0.7;
+        nonbonded->setParticleParameters(0, q, 1.5, 0);
+        nonbonded->setParticleParameters(i, q, 1.5, 0);
+        nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
+        nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
+        context.reinitialize();
+        context.setPositions(positions);
+        state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces2 = state.getForces();
+        force = ONE_4PI_EPS0*q*q/(r*r);
+        energy = ONE_4PI_EPS0*q*q/r;
+        if (i == 3) {
+            force /= 1.2;
+            energy /= 1.2;
+        }
+        if (i < 3 || r > cutoff) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+    }
+}
+
+void testPeriodic() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addException(0, 1, 0.0, 1.0, 0.0);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    positions[2] = Vec3(3, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    const double eps = 78.3;
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    const double force = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
+}
+
+
+void testLargeSystem() {
+    const int numMolecules = 600;
+    const int numParticles = numMolecules*2;
+    const double cutoff = 2.0;
+    const double boxSize = 20.0;
+    const double tol = 2e-3;
+    ReferencePlatform reference;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            nonbonded->addParticle(-1.0, 0.2, 0.1);
+            nonbonded->addParticle(1.0, 0.1, 0.1);
+        }
+        else {
+            nonbonded->addParticle(-1.0, 0.2, 0.2);
+            nonbonded->addParticle(1.0, 0.1, 0.2);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
+        bonds->addBond(2*i, 2*i+1, 1.0, 0.1);
+        nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
+    }
+
+    // Try with cutoffs but not periodic boundary conditions, and make sure the cl and Reference
+    // platforms agree.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.addForce(nonbonded);
+    system.addForce(bonds);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context cuContext(system, integrator1, platform);
+    Context referenceContext(system, integrator2, reference);
+    cuContext.setPositions(positions);
+    cuContext.setVelocities(velocities);
+    referenceContext.setPositions(positions);
+    referenceContext.setVelocities(velocities);
+    State cuState = cuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(cuState.getPositions()[i], referenceState.getPositions()[i], tol);
+        ASSERT_EQUAL_VEC(cuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    }
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+
+    // Now do the same thing with periodic boundary conditions.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    cuContext.reinitialize();
+    referenceContext.reinitialize();
+    cuContext.setPositions(positions);
+    cuContext.setVelocities(velocities);
+    referenceContext.setPositions(positions);
+    referenceContext.setVelocities(velocities);
+    cuState = cuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++) {
+        double dx = cuState.getPositions()[i][0]-referenceState.getPositions()[i][0];
+        double dy = cuState.getPositions()[i][1]-referenceState.getPositions()[i][1];
+        double dz = cuState.getPositions()[i][2]-referenceState.getPositions()[i][2];
+        ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][0]-referenceState.getPositions()[i][0], boxSize), 0, tol);
+        ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][1]-referenceState.getPositions()[i][1], boxSize), 0, tol);
+        ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][2]-referenceState.getPositions()[i][2], boxSize), 0, tol);
+        ASSERT_EQUAL_VEC(cuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    }
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+}
+
+void testDispersionCorrection() {
+    // Create a box full of identical particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.7;
+    double cutoff = boxSize/3;
+    System system;
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    vector<Vec3> positions(numParticles);
+    int index = 0;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                system.addParticle(1.0);
+                nonbonded->addParticle(0, 1.1, 0.5);
+                positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
+                index++;
+            }
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+
+    // See if the correction has the correct value.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(false);
+    context.reinitialize();
+    context.setPositions(positions);
+    double energy2 = context.getState(State::Energy).getPotentialEnergy();
+    double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
+
+    // Now modify half the particles to be different, and see if it is still correct.
+
+    int numType2 = 0;
+    for (int i = 0; i < numParticles; i += 2) {
+        nonbonded->setParticleParameters(i, 0, 1, 1);
+        numType2++;
+    }
+    int numType1 = numParticles-numType2;
+    nonbonded->updateParametersInContext(context);
+    energy2 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy1 = context.getState(State::Energy).getPotentialEnergy();
+    term1 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    term2 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    term1 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 12)/pow(cutoff, 9))/9;
+    term2 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 6)/pow(cutoff, 3))/3;
+    double combinedSigma = 0.5*(1+1.1);
+    double combinedEpsilon = sqrt(1*0.5);
+    term1 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 12)/pow(cutoff, 9))/9;
+    term2 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 6)/pow(cutoff, 3))/3;
+    term1 /= (numParticles*(numParticles+1))/2;
+    term2 /= (numParticles*(numParticles+1))/2;
+    expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
+}
+
+void testChangingParameters() {
+    const int numMolecules = 600;
+    const int numParticles = numMolecules*2;
+    const double cutoff = 2.0;
+    const double boxSize = 20.0;
+    const double tol = 2e-3;
+    ReferencePlatform reference;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            nonbonded->addParticle(-1.0, 0.2, 0.1);
+            nonbonded->addParticle(1.0, 0.1, 0.1);
+        }
+        else {
+            nonbonded->addParticle(-1.0, 0.2, 0.2);
+            nonbonded->addParticle(1.0, 0.1, 0.2);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        system.addConstraint(2*i, 2*i+1, 1.0);
+        nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
+    }
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setCutoffDistance(cutoff);
+    system.addForce(nonbonded);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    
+    // See if Reference and CPU give the same forces and energies.
+    
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context cuContext(system, integrator1, platform);
+    Context referenceContext(system, integrator2, reference);
+    cuContext.setPositions(positions);
+    referenceContext.setPositions(positions);
+    State cuState = cuContext.getState(State::Forces | State::Energy);
+    State referenceState = referenceContext.getState(State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+    
+    // Now modify parameters and see if they still agree.
+
+    for (int i = 0; i < numParticles; i += 5) {
+        double charge, sigma, epsilon;
+        nonbonded->getParticleParameters(i, charge, sigma, epsilon);
+        nonbonded->setParticleParameters(i, 1.5*charge, 1.1*sigma, 1.7*epsilon);
+    }
+    nonbonded->updateParametersInContext(cuContext);
+    nonbonded->updateParametersInContext(referenceContext);
+    cuState = cuContext.getState(State::Forces | State::Energy);
+    referenceState = referenceContext.getState(State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
+    ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
+}
+
+void testSwitchingFunction(NonbondedForce::NonbondedMethod method) {
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(0, 1.2, 1);
+    nonbonded->addParticle(0, 1.4, 2);
+    nonbonded->setNonbondedMethod(method);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    nonbonded->setUseDispersionCorrection(false);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    double eps = SQRT_TWO;
+    
+    // Compute the interaction at various distances.
+    
+    for (double r = 1.0; r < 2.5; r += 0.1) {
+        positions[1] = Vec3(r, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        
+        // See if the energy is correct.
+        
+        double x = 1.3/r;
+        double expectedEnergy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        double switchValue;
+        if (r <= 1.5)
+            switchValue = 1;
+        else if (r >= 2.0)
+            switchValue = 0;
+        else {
+            double t = (r-1.5)/0.5;
+            switchValue = 1+t*t*t*(-10+t*(15-t*6));
+        }
+        ASSERT_EQUAL_TOL(switchValue*expectedEnergy, state.getPotentialEnergy(), TOL);
+        
+        // See if the force is the gradient of the energy.
+        
+        double delta = 1e-3;
+        positions[1] = Vec3(r-delta, 0, 0);
+        context.setPositions(positions);
+        double e1 = context.getState(State::Energy).getPotentialEnergy();
+        positions[1] = Vec3(r+delta, 0, 0);
+        context.setPositions(positions);
+        double e2 = context.getState(State::Energy).getPotentialEnergy();
+        ASSERT_EQUAL_TOL((e2-e1)/(2*delta), state.getForces()[0][0], 1e-3);
+    }
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        if (!CpuPlatform::isProcessorSupported()) {
+            cout << "CPU is not supported.  Exiting." << endl;
+            return 0;
+        }
+        testCoulomb();
+        testLJ();
+        testExclusionsAnd14();
+        testCutoff();
+        testCutoff14();
+        testPeriodic();
+        testLargeSystem();
+        testDispersionCorrection();
+        testChangingParameters();
+        testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
+        testSwitchingFunction(NonbondedForce::PME);
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cpu/tests/nacl_amorph.dat

+positions[0] = Vec3(1.066000,1.628000,0.835000);
+positions[1] = Vec3(1.072000,0.428000,0.190000);
+positions[2] = Vec3(0.524000,1.442000,1.160000);
+positions[3] = Vec3(2.383000,1.524000,1.119000);
+positions[4] = Vec3(0.390000,1.441000,0.575000);
+positions[5] = Vec3(0.618000,0.399000,0.819000);
+positions[6] = Vec3(1.003000,1.257000,1.543000);
+positions[7] = Vec3(2.933000,1.569000,0.642000);
+positions[8] = Vec3(0.849000,0.739000,0.089000);
+positions[9] = Vec3(0.060000,0.794000,0.766000);
+positions[10] = Vec3(1.652000,1.405000,1.010000);
+positions[11] = Vec3(2.843000,1.533000,1.781000);
+positions[12] = Vec3(0.952000,1.309000,0.996000);
+positions[13] = Vec3(1.847000,1.402000,0.313000);
+positions[14] = Vec3(2.674000,0.083000,1.691000);
+positions[15] = Vec3(1.763000,2.104000,0.728000);
+positions[16] = Vec3(0.914000,0.574000,0.982000);
+positions[17] = Vec3(0.514000,0.078000,0.891000);
+positions[18] = Vec3(0.538000,0.766000,1.110000);
+positions[19] = Vec3(0.808000,0.676000,0.570000);
+positions[20] = Vec3(0.178000,0.014000,0.628000);
+positions[21] = Vec3(1.329000,1.333000,0.339000);
+positions[22] = Vec3(1.029000,1.678000,0.503000);
+positions[23] = Vec3(1.423000,1.767000,1.104000);
+positions[24] = Vec3(1.966000,1.051000,0.282000);
+positions[25] = Vec3(1.596000,1.971000,0.194000);
+positions[26] = Vec3(1.025000,1.043000,2.809000);
+positions[27] = Vec3(1.628000,2.614000,0.088000);
+positions[28] = Vec3(0.440000,0.606000,0.141000);
+positions[29] = Vec3(1.050000,2.821000,2.517000);
+positions[30] = Vec3(0.644000,1.604000,0.770000);
+positions[31] = Vec3(0.637000,0.917000,0.392000);
+positions[32] = Vec3(0.611000,2.768000,0.013000);
+positions[33] = Vec3(1.892000,0.660000,0.473000);
+positions[34] = Vec3(1.052000,2.081000,0.982000);
+positions[35] = Vec3(1.508000,2.300000,0.439000);
+positions[36] = Vec3(2.617000,0.328000,1.099000);
+positions[37] = Vec3(0.910000,0.040000,0.259000);
+positions[38] = Vec3(1.195000,1.494000,1.202000);
+positions[39] = Vec3(2.657000,0.997000,0.564000);
+positions[40] = Vec3(1.465000,1.580000,0.648000);
+positions[41] = Vec3(0.154000,2.538000,1.331000);
+positions[42] = Vec3(0.849000,1.476000,1.365000);
+positions[43] = Vec3(0.898000,0.987000,1.178000);
+positions[44] = Vec3(0.958000,0.656000,1.358000);
+positions[45] = Vec3(1.067000,0.934000,0.211000);
+positions[46] = Vec3(1.030000,0.319000,1.281000);
+positions[47] = Vec3(2.709000,0.807000,0.240000);
+positions[48] = Vec3(0.837000,1.362000,0.588000);
+positions[49] = Vec3(2.080000,0.791000,2.947000);
+positions[50] = Vec3(0.200000,0.266000,1.474000);
+positions[51] = Vec3(0.848000,0.379000,1.625000);
+positions[52] = Vec3(0.637000,1.071000,0.821000);
+positions[53] = Vec3(1.324000,0.757000,2.951000);
+positions[54] = Vec3(2.666000,0.935000,1.373000);
+positions[55] = Vec3(1.584000,1.025000,1.703000);
+positions[56] = Vec3(1.699000,0.636000,0.038000);
+positions[57] = Vec3(1.099000,1.644000,1.879000);
+positions[58] = Vec3(2.897000,1.302000,1.522000);
+positions[59] = Vec3(1.753000,0.949000,2.885000);
+positions[60] = Vec3(2.502000,1.321000,0.752000);
+positions[61] = Vec3(0.545000,0.193000,1.959000);
+positions[62] = Vec3(1.098000,2.646000,1.706000);
+positions[63] = Vec3(0.001000,1.205000,0.670000);
+positions[64] = Vec3(2.997000,0.061000,1.040000);
+positions[65] = Vec3(0.662000,0.828000,1.535000);
+positions[66] = Vec3(1.252000,1.246000,0.780000);
+positions[67] = Vec3(1.173000,0.472000,0.810000);
+positions[68] = Vec3(0.124000,0.622000,2.992000);
+positions[69] = Vec3(1.036000,0.883000,0.848000);
+positions[70] = Vec3(1.423000,2.146000,1.340000);
+positions[71] = Vec3(2.391000,1.136000,1.165000);
+positions[72] = Vec3(1.189000,2.961000,0.425000);
+positions[73] = Vec3(1.584000,2.500000,0.782000);
+positions[74] = Vec3(0.565000,1.122000,1.240000);
+positions[75] = Vec3(1.733000,1.716000,1.763000);
+positions[76] = Vec3(1.548000,1.522000,0.041000);
+positions[77] = Vec3(1.485000,0.561000,0.369000);
+positions[78] = Vec3(0.350000,1.661000,0.928000);
+positions[79] = Vec3(1.653000,1.223000,0.578000);
+positions[80] = Vec3(0.648000,1.349000,0.253000);
+positions[81] = Vec3(0.340000,1.820000,0.483000);
+positions[82] = Vec3(2.926000,0.119000,1.421000);
+positions[83] = Vec3(1.512000,1.084000,0.156000);
+positions[84] = Vec3(1.600000,2.115000,1.792000);
+positions[85] = Vec3(1.089000,0.934000,1.584000);
+positions[86] = Vec3(1.276000,1.104000,1.230000);
+positions[87] = Vec3(0.485000,0.305000,0.428000);
+positions[88] = Vec3(1.317000,1.261000,1.795000);
+positions[89] = Vec3(0.039000,1.413000,1.085000);
+positions[90] = Vec3(0.453000,0.701000,0.605000);
+positions[91] = Vec3(1.283000,1.937000,0.752000);
+positions[92] = Vec3(0.212000,1.416000,1.447000);
+positions[93] = Vec3(0.203000,0.358000,0.723000);
+positions[94] = Vec3(0.556000,0.445000,1.364000);
+positions[95] = Vec3(1.436000,0.861000,0.911000);
+positions[96] = Vec3(0.358000,0.966000,0.176000);
+positions[97] = Vec3(1.478000,2.715000,0.427000);
+positions[98] = Vec3(1.581000,0.575000,0.809000);
+positions[99] = Vec3(1.007000,2.153000,2.887000);
+positions[100] = Vec3(2.343000,0.663000,2.513000);
+positions[101] = Vec3(2.105000,0.649000,1.635000);
+positions[102] = Vec3(0.875000,0.743000,2.459000);
+positions[103] = Vec3(0.229000,1.315000,1.879000);
+positions[104] = Vec3(0.285000,0.935000,1.700000);
+positions[105] = Vec3(2.269000,1.284000,2.234000);
+positions[106] = Vec3(1.406000,1.149000,2.767000);
+positions[107] = Vec3(1.076000,0.220000,1.849000);
+positions[108] = Vec3(2.001000,1.532000,2.881000);
+positions[109] = Vec3(2.893000,0.485000,1.860000);
+positions[110] = Vec3(1.621000,1.786000,2.624000);
+positions[111] = Vec3(0.500000,0.616000,1.818000);
+positions[112] = Vec3(0.938000,2.978000,2.104000);
+positions[113] = Vec3(0.550000,2.081000,0.454000);
+positions[114] = Vec3(1.121000,0.685000,2.196000);
+positions[115] = Vec3(1.088000,1.385000,2.184000);
+positions[116] = Vec3(1.122000,2.705000,2.080000);
+positions[117] = Vec3(0.918000,1.767000,2.861000);
+positions[118] = Vec3(2.748000,0.232000,2.126000);
+positions[119] = Vec3(1.238000,2.766000,0.109000);
+positions[120] = Vec3(1.380000,0.785000,1.961000);
+positions[121] = Vec3(1.236000,1.757000,0.150000);
+positions[122] = Vec3(1.339000,2.187000,2.592000);
+positions[123] = Vec3(1.414000,0.342000,2.714000);
+positions[124] = Vec3(1.310000,0.770000,2.589000);
+positions[125] = Vec3(1.686000,0.765000,2.321000);
+positions[126] = Vec3(1.659000,1.367000,2.780000);
+positions[127] = Vec3(0.141000,0.095000,1.903000);
+positions[128] = Vec3(2.084000,1.002000,2.520000);
+positions[129] = Vec3(2.819000,1.286000,2.626000);
+positions[130] = Vec3(1.257000,1.044000,2.401000);
+positions[131] = Vec3(1.064000,0.546000,2.839000);
+positions[132] = Vec3(0.078000,1.246000,0.010000);
+positions[133] = Vec3(1.506000,0.420000,2.223000);
+positions[134] = Vec3(1.778000,0.699000,1.920000);
+positions[135] = Vec3(1.315000,1.721000,2.733000);
+positions[136] = Vec3(0.114000,0.281000,0.279000);
+positions[137] = Vec3(1.082000,1.421000,2.596000);
+positions[138] = Vec3(3.001000,0.592000,2.276000);
+positions[139] = Vec3(1.384000,2.227000,2.992000);
+positions[140] = Vec3(1.353000,0.044000,1.985000);
+positions[141] = Vec3(1.367000,1.832000,2.383000);
+positions[142] = Vec3(0.853000,1.119000,2.230000);
+positions[143] = Vec3(1.675000,1.482000,2.295000);
+positions[144] = Vec3(1.334000,1.890000,1.904000);
+positions[145] = Vec3(1.630000,0.140000,2.939000);
+positions[146] = Vec3(0.195000,1.290000,2.300000);
+positions[147] = Vec3(2.178000,1.173000,3.001000);
+positions[148] = Vec3(0.637000,0.655000,2.126000);
+positions[149] = Vec3(0.993000,1.796000,2.529000);
+positions[150] = Vec3(0.910000,0.701000,1.845000);
+positions[151] = Vec3(0.191000,2.128000,0.355000);
+positions[152] = Vec3(0.692000,1.163000,2.578000);
+positions[153] = Vec3(1.154000,1.052000,1.974000);
+positions[154] = Vec3(1.682000,0.335000,2.509000);
+positions[155] = Vec3(0.569000,1.032000,1.895000);
+positions[156] = Vec3(1.800000,2.796000,1.295000);
+positions[157] = Vec3(2.517000,2.347000,2.878000);
+positions[158] = Vec3(0.639000,2.470000,1.678000);
+positions[159] = Vec3(0.634000,2.006000,1.829000);
+positions[160] = Vec3(0.892000,0.215000,0.566000);
+positions[161] = Vec3(1.800000,2.143000,1.491000);
+positions[162] = Vec3(1.898000,0.226000,2.765000);
+positions[163] = Vec3(0.791000,1.738000,0.260000);
+positions[164] = Vec3(0.437000,1.740000,2.048000);
+positions[165] = Vec3(1.687000,2.438000,1.166000);
+positions[166] = Vec3(1.337000,2.304000,1.643000);
+positions[167] = Vec3(1.270000,2.397000,1.033000);
+positions[168] = Vec3(0.702000,2.429000,0.591000);
+positions[169] = Vec3(0.842000,1.976000,0.724000);
+positions[170] = Vec3(1.965000,0.095000,1.206000);
+positions[171] = Vec3(0.355000,2.710000,0.618000);
+positions[172] = Vec3(0.745000,1.434000,2.781000);
+positions[173] = Vec3(0.707000,2.171000,1.502000);
+positions[174] = Vec3(1.294000,2.696000,0.847000);
+positions[175] = Vec3(1.143000,2.075000,0.276000);
+positions[176] = Vec3(1.111000,2.474000,0.312000);
+positions[177] = Vec3(1.144000,2.316000,0.633000);
+positions[178] = Vec3(1.335000,0.292000,0.515000);
+positions[179] = Vec3(1.926000,2.813000,2.703000);
+positions[180] = Vec3(0.559000,2.314000,2.904000);
+positions[181] = Vec3(1.308000,1.605000,1.534000);
+positions[182] = Vec3(0.773000,2.913000,1.217000);
+positions[183] = Vec3(1.612000,0.082000,1.027000);
+positions[184] = Vec3(1.510000,0.287000,1.787000);
+positions[185] = Vec3(0.716000,1.424000,1.843000);
+positions[186] = Vec3(1.267000,0.685000,1.160000);
+positions[187] = Vec3(0.306000,1.653000,1.717000);
+positions[188] = Vec3(0.349000,0.020000,1.275000);
+positions[189] = Vec3(0.166000,1.979000,0.804000);
+positions[190] = Vec3(1.523000,2.992000,0.711000);
+positions[191] = Vec3(1.998000,2.146000,0.088000);
+positions[192] = Vec3(0.047000,2.513000,0.642000);
+positions[193] = Vec3(0.501000,1.793000,1.438000);
+positions[194] = Vec3(1.099000,2.010000,1.626000);
+positions[195] = Vec3(2.580000,2.854000,1.328000);
+positions[196] = Vec3(1.080000,2.779000,1.190000);
+positions[197] = Vec3(0.901000,2.561000,0.948000);
+positions[198] = Vec3(0.920000,2.990000,0.844000);
+positions[199] = Vec3(0.819000,2.924000,1.711000);
+positions[200] = Vec3(0.434000,1.516000,0.063000);
+positions[201] = Vec3(1.470000,0.058000,0.231000);
+positions[202] = Vec3(0.530000,3.005000,1.550000);
+positions[203] = Vec3(0.447000,2.330000,1.277000);
+positions[204] = Vec3(1.632000,2.683000,1.593000);
+positions[205] = Vec3(0.885000,1.835000,2.072000);
+positions[206] = Vec3(0.868000,2.601000,1.425000);
+positions[207] = Vec3(0.720000,2.242000,0.907000);
+positions[208] = Vec3(1.194000,0.144000,1.065000);
+positions[209] = Vec3(0.448000,2.485000,0.959000);
+positions[210] = Vec3(1.377000,2.694000,1.352000);
+positions[211] = Vec3(1.305000,2.928000,2.713000);
+positions[212] = Vec3(1.784000,2.456000,1.981000);
+positions[213] = Vec3(0.354000,2.136000,1.563000);
+positions[214] = Vec3(0.489000,2.000000,1.108000);
+positions[215] = Vec3(1.884000,2.221000,0.461000);
+positions[216] = Vec3(1.860000,2.540000,0.306000);
+positions[217] = Vec3(1.753000,2.335000,2.768000);
+positions[218] = Vec3(1.536000,2.441000,2.344000);
+positions[219] = Vec3(0.531000,0.025000,2.235000);
+positions[220] = Vec3(0.809000,0.011000,2.834000);
+positions[221] = Vec3(0.289000,2.614000,2.879000);
+positions[222] = Vec3(0.613000,1.891000,2.337000);
+positions[223] = Vec3(0.507000,0.037000,2.694000);
+positions[224] = Vec3(0.882000,2.185000,2.583000);
+positions[225] = Vec3(0.503000,2.051000,2.615000);
+positions[226] = Vec3(1.907000,1.956000,2.831000);
+positions[227] = Vec3(2.833000,2.769000,1.644000);
+positions[228] = Vec3(1.141000,0.113000,2.945000);
+positions[229] = Vec3(0.600000,1.338000,2.200000);
+positions[230] = Vec3(0.904000,2.360000,1.952000);
+positions[231] = Vec3(0.738000,1.568000,2.437000);
+positions[232] = Vec3(1.136000,2.535000,2.805000);
+positions[233] = Vec3(1.430000,2.767000,2.321000);
+positions[234] = Vec3(1.078000,2.470000,2.385000);
+positions[235] = Vec3(0.296000,2.376000,2.560000);
+positions[236] = Vec3(0.719000,0.300000,0.075000);
+positions[237] = Vec3(0.518000,1.911000,0.080000);
+positions[238] = Vec3(0.381000,1.570000,2.450000);
+positions[239] = Vec3(0.716000,2.581000,2.697000);
+positions[240] = Vec3(1.473000,2.617000,1.936000);
+positions[241] = Vec3(0.421000,2.449000,0.229000);
+positions[242] = Vec3(0.425000,2.817000,1.910000);
+positions[243] = Vec3(1.312000,2.294000,2.057000);
+positions[244] = Vec3(1.239000,0.007000,1.539000);
+positions[245] = Vec3(0.822000,0.379000,2.086000);
+positions[246] = Vec3(0.560000,2.562000,2.227000);
+positions[247] = Vec3(0.863000,2.417000,0.050000);
+positions[248] = Vec3(1.263000,0.151000,2.332000);
+positions[249] = Vec3(0.237000,0.208000,2.336000);
+positions[250] = Vec3(0.437000,2.370000,1.910000);
+positions[251] = Vec3(1.119000,2.058000,2.207000);
+positions[252] = Vec3(1.960000,1.749000,0.118000);
+positions[253] = Vec3(2.415000,0.870000,2.757000);
+positions[254] = Vec3(1.781000,0.342000,0.366000);
+positions[255] = Vec3(2.172000,1.279000,1.421000);
+positions[256] = Vec3(1.986000,0.715000,1.301000);
+positions[257] = Vec3(1.657000,1.804000,0.810000);
+positions[258] = Vec3(2.405000,1.202000,0.416000);
+positions[259] = Vec3(1.932000,1.457000,0.786000);
+positions[260] = Vec3(1.901000,1.271000,1.207000);
+positions[261] = Vec3(1.864000,0.301000,0.810000);
+positions[262] = Vec3(1.658000,0.673000,1.558000);
+positions[263] = Vec3(2.637000,2.247000,0.396000);
+positions[264] = Vec3(1.353000,0.369000,1.438000);
+positions[265] = Vec3(0.530000,2.688000,1.346000);
+positions[266] = Vec3(0.237000,0.485000,1.047000);
+positions[267] = Vec3(2.806000,0.601000,0.822000);
+positions[268] = Vec3(1.617000,2.018000,2.136000);
+positions[269] = Vec3(2.000000,2.898000,0.022000);
+positions[270] = Vec3(2.049000,1.883000,1.001000);
+positions[271] = Vec3(2.477000,0.355000,1.786000);
+positions[272] = Vec3(1.646000,0.983000,1.266000);
+positions[273] = Vec3(1.683000,2.097000,1.114000);
+positions[274] = Vec3(2.161000,0.921000,1.065000);
+positions[275] = Vec3(2.099000,0.463000,1.942000);
+positions[276] = Vec3(2.561000,1.638000,0.572000);
+positions[277] = Vec3(2.205000,0.395000,1.005000);
+positions[278] = Vec3(2.836000,0.203000,0.698000);
+positions[279] = Vec3(2.662000,0.909000,0.966000);
+positions[280] = Vec3(0.334000,0.350000,2.767000);
+positions[281] = Vec3(2.241000,2.934000,1.248000);
+positions[282] = Vec3(2.599000,2.953000,0.921000);
+positions[283] = Vec3(2.219000,0.262000,0.058000);
+positions[284] = Vec3(0.274000,0.656000,1.456000);
+positions[285] = Vec3(1.814000,1.008000,0.882000);
+positions[286] = Vec3(2.793000,1.395000,0.316000);
+positions[287] = Vec3(0.773000,1.753000,1.639000);
+positions[288] = Vec3(2.321000,0.994000,1.591000);
+positions[289] = Vec3(2.243000,2.255000,1.690000);
+positions[290] = Vec3(0.178000,1.342000,0.327000);
+positions[291] = Vec3(1.623000,1.756000,1.426000);
+positions[292] = Vec3(2.252000,0.109000,0.375000);
+positions[293] = Vec3(3.003000,1.895000,1.895000);
+positions[294] = Vec3(0.407000,0.831000,2.756000);
+positions[295] = Vec3(2.193000,0.956000,0.632000);
+positions[296] = Vec3(2.405000,0.641000,1.107000);
+positions[297] = Vec3(2.361000,0.958000,0.162000);
+positions[298] = Vec3(2.173000,1.544000,0.528000);
+positions[299] = Vec3(1.565000,1.380000,1.428000);
+positions[300] = Vec3(2.342000,0.538000,0.253000);
+positions[301] = Vec3(1.910000,0.701000,0.954000);
+positions[302] = Vec3(2.910000,0.288000,2.938000);
+positions[303] = Vec3(0.257000,1.189000,0.958000);
+positions[304] = Vec3(0.134000,1.775000,1.243000);
+positions[305] = Vec3(2.476000,1.583000,1.956000);
+positions[306] = Vec3(1.838000,1.791000,2.354000);
+positions[307] = Vec3(1.906000,1.338000,1.696000);
+positions[308] = Vec3(2.413000,2.869000,0.166000);
+positions[309] = Vec3(3.006000,1.038000,1.322000);
+positions[310] = Vec3(1.961000,0.962000,1.605000);
+positions[311] = Vec3(0.082000,2.857000,0.020000);
+positions[312] = Vec3(2.408000,1.499000,0.062000);
+positions[313] = Vec3(2.349000,0.267000,1.415000);
+positions[314] = Vec3(2.327000,1.717000,2.350000);
+positions[315] = Vec3(2.928000,0.810000,1.582000);
+positions[316] = Vec3(2.150000,0.336000,0.576000);
+positions[317] = Vec3(2.664000,1.085000,2.962000);
+positions[318] = Vec3(2.851000,0.670000,1.174000);
+positions[319] = Vec3(1.954000,1.013000,1.975000);
+positions[320] = Vec3(2.474000,1.542000,1.545000);
+positions[321] = Vec3(2.826000,0.455000,1.490000);
+positions[322] = Vec3(2.140000,2.826000,0.558000);
+positions[323] = Vec3(2.151000,1.684000,1.780000);
+positions[324] = Vec3(0.174000,0.673000,0.397000);
+positions[325] = Vec3(0.066000,1.708000,0.160000);
+positions[326] = Vec3(2.158000,0.303000,2.582000);
+positions[327] = Vec3(2.602000,1.611000,2.632000);
+positions[328] = Vec3(2.566000,1.138000,2.465000);
+positions[329] = Vec3(2.026000,1.443000,2.477000);
+positions[330] = Vec3(2.365000,0.309000,2.255000);
+positions[331] = Vec3(1.636000,1.107000,2.058000);
+positions[332] = Vec3(2.522000,2.584000,2.399000);
+positions[333] = Vec3(2.537000,2.900000,2.158000);
+positions[334] = Vec3(2.660000,0.537000,2.577000);
+positions[335] = Vec3(2.679000,1.158000,1.724000);
+positions[336] = Vec3(0.220000,1.894000,2.498000);
+positions[337] = Vec3(2.266000,1.248000,1.837000);
+positions[338] = Vec3(0.055000,1.656000,2.128000);
+positions[339] = Vec3(2.899000,1.902000,2.823000);
+positions[340] = Vec3(0.085000,2.994000,2.720000);
+positions[341] = Vec3(0.013000,0.889000,2.468000);
+positions[342] = Vec3(1.804000,0.372000,1.636000);
+positions[343] = Vec3(0.201000,1.616000,2.824000);
+positions[344] = Vec3(0.369000,1.273000,2.699000);
+positions[345] = Vec3(2.996000,0.355000,2.596000);
+positions[346] = Vec3(2.867000,1.314000,2.107000);
+positions[347] = Vec3(2.611000,0.563000,2.140000);
+positions[348] = Vec3(2.676000,2.954000,2.955000);
+positions[349] = Vec3(0.256000,0.848000,2.062000);
+positions[350] = Vec3(2.530000,0.028000,2.528000);
+positions[351] = Vec3(0.537000,1.273000,1.596000);
+positions[352] = Vec3(0.004000,1.004000,0.401000);
+positions[353] = Vec3(1.676000,1.060000,2.463000);
+positions[354] = Vec3(2.622000,1.473000,2.257000);
+positions[355] = Vec3(2.917000,2.991000,2.316000);
+positions[356] = Vec3(0.672000,1.123000,2.984000);
+positions[357] = Vec3(2.229000,1.806000,2.673000);
+positions[358] = Vec3(0.463000,0.951000,2.383000);
+positions[359] = Vec3(2.126000,0.049000,2.037000);
+positions[360] = Vec3(2.868000,0.876000,2.015000);
+positions[361] = Vec3(2.720000,2.582000,0.079000);
+positions[362] = Vec3(1.966000,0.693000,2.624000);
+positions[363] = Vec3(1.971000,0.398000,2.318000);
+positions[364] = Vec3(0.337000,0.630000,2.458000);
+positions[365] = Vec3(2.562000,1.044000,2.040000);
+positions[366] = Vec3(2.817000,1.485000,2.963000);
+positions[367] = Vec3(2.514000,0.621000,2.992000);
+positions[368] = Vec3(3.000000,1.551000,2.496000);
+positions[369] = Vec3(0.698000,2.167000,2.180000);
+positions[370] = Vec3(2.693000,0.849000,2.389000);
+positions[371] = Vec3(2.092000,2.565000,2.986000);
+positions[372] = Vec3(2.010000,3.001000,0.819000);
+positions[373] = Vec3(2.392000,2.622000,1.636000);
+positions[374] = Vec3(2.086000,2.325000,1.340000);
+positions[375] = Vec3(2.578000,2.971000,0.502000);
+positions[376] = Vec3(1.871000,2.789000,2.225000);
+positions[377] = Vec3(2.230000,2.985000,1.594000);
+positions[378] = Vec3(2.860000,2.788000,0.729000);
+positions[379] = Vec3(2.051000,1.928000,1.472000);
+positions[380] = Vec3(2.307000,2.219000,1.067000);
+positions[381] = Vec3(2.369000,2.572000,1.289000);
+positions[382] = Vec3(2.206000,1.924000,0.693000);
+positions[383] = Vec3(1.984000,2.058000,2.005000);
+positions[384] = Vec3(2.287000,1.854000,0.317000);
+positions[385] = Vec3(2.525000,0.345000,0.686000);
+positions[386] = Vec3(2.933000,1.920000,1.053000);
+positions[387] = Vec3(0.324000,2.324000,0.601000);
+positions[388] = Vec3(2.042000,1.576000,1.277000);
+positions[389] = Vec3(0.031000,2.376000,0.949000);
+positions[390] = Vec3(2.519000,2.250000,1.465000);
+positions[391] = Vec3(0.221000,2.722000,1.652000);
+positions[392] = Vec3(2.409000,2.361000,2.051000);
+positions[393] = Vec3(2.472000,1.917000,1.673000);
+positions[394] = Vec3(0.999000,2.715000,0.562000);
+positions[395] = Vec3(1.669000,0.017000,1.508000);
+positions[396] = Vec3(1.924000,1.777000,0.542000);
+positions[397] = Vec3(2.635000,2.634000,1.905000);
+positions[398] = Vec3(2.042000,2.628000,1.025000);
+positions[399] = Vec3(2.694000,1.974000,2.009000);
+positions[400] = Vec3(2.988000,2.221000,1.333000);
+positions[401] = Vec3(1.772000,0.196000,1.978000);
+positions[402] = Vec3(2.893000,2.961000,0.283000);
+positions[403] = Vec3(2.615000,0.261000,0.245000);
+positions[404] = Vec3(2.797000,2.521000,1.412000);
+positions[405] = Vec3(0.013000,2.497000,0.246000);
+positions[406] = Vec3(1.875000,2.861000,1.801000);
+positions[407] = Vec3(2.800000,2.617000,1.049000);
+positions[408] = Vec3(2.824000,1.858000,1.487000);
+positions[409] = Vec3(2.434000,1.868000,1.275000);
+positions[410] = Vec3(2.814000,0.526000,0.384000);
+positions[411] = Vec3(2.844000,2.545000,2.246000);
+positions[412] = Vec3(1.896000,2.587000,0.719000);
+positions[413] = Vec3(0.350000,0.055000,0.076000);
+positions[414] = Vec3(2.686000,1.784000,0.222000);
+positions[415] = Vec3(2.724000,1.604000,0.989000);
+positions[416] = Vec3(0.807000,1.761000,1.122000);
+positions[417] = Vec3(2.120000,2.382000,2.226000);
+positions[418] = Vec3(2.058000,1.587000,2.067000);
+positions[419] = Vec3(2.904000,2.571000,2.686000);
+positions[420] = Vec3(2.228000,2.910000,2.410000);
+positions[421] = Vec3(2.797000,2.142000,0.114000);
+positions[422] = Vec3(2.905000,1.875000,0.480000);
+positions[423] = Vec3(1.881000,2.565000,2.469000);
+positions[424] = Vec3(2.404000,1.929000,2.999000);
+positions[425] = Vec3(2.389000,2.814000,2.782000);
+positions[426] = Vec3(2.520000,0.301000,2.815000);
+positions[427] = Vec3(2.726000,1.907000,2.339000);
+positions[428] = Vec3(2.880000,2.273000,2.500000);
+positions[429] = Vec3(2.574000,2.045000,2.716000);
+positions[430] = Vec3(2.988000,2.288000,2.001000);
+positions[431] = Vec3(0.011000,2.341000,2.904000);
+positions[432] = Vec3(0.215000,2.265000,2.257000);
+positions[433] = Vec3(2.268000,2.311000,0.234000);
+positions[434] = Vec3(2.462000,2.621000,0.550000);
+positions[435] = Vec3(1.530000,2.540000,2.728000);
+positions[436] = Vec3(2.162000,2.306000,2.687000);
+positions[437] = Vec3(2.748000,2.301000,1.734000);
+positions[438] = Vec3(2.334000,1.976000,2.041000);
+positions[439] = Vec3(1.981000,2.076000,2.443000);
+positions[440] = Vec3(2.301000,1.367000,2.665000);
+positions[441] = Vec3(2.399000,2.164000,2.403000);
+positions[442] = Vec3(0.244000,2.713000,2.257000);
+positions[443] = Vec3(0.683000,0.488000,2.781000);
+positions[444] = Vec3(2.194000,2.711000,1.993000);
+positions[445] = Vec3(2.947000,2.848000,2.001000);
+positions[446] = Vec3(0.223000,1.981000,2.913000);
+positions[447] = Vec3(0.010000,1.226000,0.917000);
+positions[448] = Vec3(1.911000,0.426000,0.582000);
+positions[449] = Vec3(2.204000,0.015000,0.136000);
+positions[450] = Vec3(0.927000,0.138000,1.645000);
+positions[451] = Vec3(0.155000,0.885000,1.479000);
+positions[452] = Vec3(1.550000,1.933000,1.261000);
+positions[453] = Vec3(1.304000,0.407000,0.287000);
+positions[454] = Vec3(0.270000,1.384000,2.910000);
+positions[455] = Vec3(0.516000,1.817000,1.695000);
+positions[456] = Vec3(1.458000,2.879000,1.523000);
+positions[457] = Vec3(1.702000,1.670000,0.593000);
+positions[458] = Vec3(1.974000,1.380000,0.534000);
+positions[459] = Vec3(2.835000,1.185000,0.479000);
+positions[460] = Vec3(0.548000,2.979000,1.126000);
+positions[461] = Vec3(1.202000,2.174000,1.466000);
+positions[462] = Vec3(1.237000,1.701000,0.653000);
+positions[463] = Vec3(2.939000,0.761000,0.349000);
+positions[464] = Vec3(1.667000,2.119000,0.377000);
+positions[465] = Vec3(1.196000,0.552000,1.372000);
+positions[466] = Vec3(1.416000,0.901000,1.178000);
+positions[467] = Vec3(0.519000,1.577000,2.227000);
+positions[468] = Vec3(1.214000,1.281000,1.063000);
+positions[469] = Vec3(0.822000,0.433000,1.375000);
+positions[470] = Vec3(0.095000,2.760000,0.604000);
+positions[471] = Vec3(1.325000,2.144000,1.848000);
+positions[472] = Vec3(0.681000,0.896000,1.285000);
+positions[473] = Vec3(0.406000,2.936000,0.717000);
+positions[474] = Vec3(0.565000,1.852000,0.349000);
+positions[475] = Vec3(0.597000,1.651000,1.020000);
+positions[476] = Vec3(1.236000,0.170000,1.335000);
+positions[477] = Vec3(0.586000,0.441000,1.980000);
+positions[478] = Vec3(1.443000,1.208000,1.575000);
+positions[479] = Vec3(0.247000,0.243000,0.502000);
+positions[480] = Vec3(1.386000,1.564000,0.236000);
+positions[481] = Vec3(0.871000,1.063000,0.930000);
+positions[482] = Vec3(0.136000,0.992000,0.621000);
+positions[483] = Vec3(0.889000,0.986000,0.010000);
+positions[484] = Vec3(1.107000,2.731000,1.452000);
+positions[485] = Vec3(0.942000,2.471000,0.517000);
+positions[486] = Vec3(0.989000,0.652000,0.747000);
+positions[487] = Vec3(0.899000,1.235000,2.707000);
+positions[488] = Vec3(1.105000,0.684000,0.068000);
+positions[489] = Vec3(1.660000,1.235000,2.276000);
+positions[490] = Vec3(1.593000,1.883000,1.915000);
+positions[491] = Vec3(1.528000,1.613000,0.920000);
+positions[492] = Vec3(0.459000,1.046000,1.011000);
+positions[493] = Vec3(0.213000,0.612000,0.644000);
+positions[494] = Vec3(0.078000,1.392000,1.676000);
+positions[495] = Vec3(0.605000,0.491000,0.574000);
+positions[496] = Vec3(0.990000,1.586000,1.076000);
+positions[497] = Vec3(0.297000,1.434000,1.028000);
+positions[498] = Vec3(1.101000,1.471000,1.443000);
+positions[499] = Vec3(0.072000,0.139000,1.653000);
+positions[500] = Vec3(0.633000,0.884000,0.645000);
+positions[501] = Vec3(0.352000,2.841000,1.463000);
+positions[502] = Vec3(0.418000,0.774000,0.350000);
+positions[503] = Vec3(2.641000,0.198000,0.869000);
+positions[504] = Vec3(0.608000,1.341000,0.695000);
+positions[505] = Vec3(1.778000,1.151000,1.830000);
+positions[506] = Vec3(1.669000,0.342000,2.768000);
+positions[507] = Vec3(1.256000,0.994000,0.798000);
+positions[508] = Vec3(1.068000,0.375000,1.036000);
+positions[509] = Vec3(0.910000,0.758000,1.589000);
+positions[510] = Vec3(0.243000,2.452000,0.805000);
+positions[511] = Vec3(1.018000,0.764000,1.122000);
+positions[512] = Vec3(2.464000,1.089000,1.404000);
+positions[513] = Vec3(0.670000,0.564000,0.034000);
+positions[514] = Vec3(0.030000,1.296000,1.310000);
+positions[515] = Vec3(1.210000,1.785000,1.691000);
+positions[516] = Vec3(0.022000,0.620000,0.974000);
+positions[517] = Vec3(1.499000,1.277000,2.986000);
+positions[518] = Vec3(1.227000,1.896000,1.006000);
+positions[519] = Vec3(0.528000,1.022000,1.635000);
+positions[520] = Vec3(1.887000,2.670000,0.089000);
+positions[521] = Vec3(1.661000,0.825000,0.793000);
+positions[522] = Vec3(0.831000,1.494000,0.374000);
+positions[523] = Vec3(1.356000,0.613000,0.930000);
+positions[524] = Vec3(0.667000,0.600000,0.968000);
+positions[525] = Vec3(1.154000,1.702000,2.925000);
+positions[526] = Vec3(1.420000,1.581000,1.289000);
+positions[527] = Vec3(1.383000,0.041000,0.932000);
+positions[528] = Vec3(1.727000,0.140000,1.725000);
+positions[529] = Vec3(0.711000,1.215000,2.004000);
+positions[530] = Vec3(1.061000,1.067000,1.366000);
+positions[531] = Vec3(0.377000,0.597000,1.224000);
+positions[532] = Vec3(0.274000,0.719000,1.842000);
+positions[533] = Vec3(0.840000,1.658000,1.874000);
+positions[534] = Vec3(0.877000,0.290000,0.311000);
+positions[535] = Vec3(2.130000,1.153000,1.196000);
+positions[536] = Vec3(1.028000,1.379000,0.747000);
+positions[537] = Vec3(1.107000,2.450000,2.079000);
+positions[538] = Vec3(1.419000,1.333000,0.585000);
+positions[539] = Vec3(0.430000,1.305000,1.369000);
+positions[540] = Vec3(0.775000,1.363000,1.596000);
+positions[541] = Vec3(1.522000,2.009000,0.736000);
+positions[542] = Vec3(0.857000,1.722000,0.696000);
+positions[543] = Vec3(0.722000,2.831000,1.478000);
+positions[544] = Vec3(0.411000,1.673000,0.681000);
+positions[545] = Vec3(1.511000,0.456000,0.597000);
+positions[546] = Vec3(2.684000,0.820000,2.996000);
+positions[547] = Vec3(1.593000,1.713000,2.369000);
+positions[548] = Vec3(1.113000,0.803000,1.958000);
+positions[549] = Vec3(1.267000,1.095000,0.254000);
+positions[550] = Vec3(2.120000,0.540000,2.477000);
+positions[551] = Vec3(0.566000,1.409000,2.588000);
+positions[552] = Vec3(0.261000,0.872000,2.546000);
+positions[553] = Vec3(1.878000,1.446000,2.680000);
+positions[554] = Vec3(0.878000,1.606000,2.658000);
+positions[555] = Vec3(1.564000,0.749000,1.786000);
+positions[556] = Vec3(1.412000,1.942000,2.625000);
+positions[557] = Vec3(1.660000,1.114000,2.710000);
+positions[558] = Vec3(1.118000,0.813000,2.424000);
+positions[559] = Vec3(1.482000,0.893000,2.434000);
+positions[560] = Vec3(1.093000,1.129000,1.740000);
+positions[561] = Vec3(2.163000,0.849000,2.709000);
+positions[562] = Vec3(1.201000,1.429000,1.957000);
+positions[563] = Vec3(0.235000,2.258000,2.002000);
+positions[564] = Vec3(0.413000,1.444000,0.314000);
+positions[565] = Vec3(0.164000,0.450000,2.408000);
+positions[566] = Vec3(1.551000,0.851000,0.033000);
+positions[567] = Vec3(0.659000,0.228000,2.807000);
+positions[568] = Vec3(0.741000,0.131000,2.124000);
+positions[569] = Vec3(0.455000,0.567000,2.682000);
+positions[570] = Vec3(0.729000,0.971000,2.408000);
+positions[571] = Vec3(1.487000,2.820000,0.162000);
+positions[572] = Vec3(1.855000,0.700000,2.858000);
+positions[573] = Vec3(0.305000,1.074000,1.926000);
+positions[574] = Vec3(1.300000,0.153000,1.747000);
+positions[575] = Vec3(1.272000,1.249000,2.568000);
+positions[576] = Vec3(0.431000,0.743000,2.238000);
+positions[577] = Vec3(0.493000,0.240000,0.184000);
+positions[578] = Vec3(1.734000,0.506000,2.317000);
+positions[579] = Vec3(0.874000,0.631000,2.692000);
+positions[580] = Vec3(0.473000,2.790000,2.161000);
+positions[581] = Vec3(1.310000,0.571000,2.759000);
+positions[582] = Vec3(0.677000,0.798000,1.916000);
+positions[583] = Vec3(0.944000,0.442000,1.858000);
+positions[584] = Vec3(3.006000,2.098000,2.976000);
+positions[585] = Vec3(0.864000,0.592000,2.231000);
+positions[586] = Vec3(1.366000,0.611000,2.147000);
+positions[587] = Vec3(2.871000,1.217000,2.880000);
+positions[588] = Vec3(1.674000,2.664000,2.336000);
+positions[589] = Vec3(1.757000,0.879000,2.101000);
+positions[590] = Vec3(1.293000,2.939000,2.457000);
+positions[591] = Vec3(1.108000,1.131000,2.206000);
+positions[592] = Vec3(1.207000,1.658000,2.498000);
+positions[593] = Vec3(0.850000,1.373000,2.312000);
+positions[594] = Vec3(1.413000,1.060000,1.939000);
+positions[595] = Vec3(1.138000,0.140000,2.102000);
+positions[596] = Vec3(0.752000,1.307000,1.190000);
+positions[597] = Vec3(1.254000,0.942000,2.790000);
+positions[598] = Vec3(1.544000,1.614000,2.800000);
+positions[599] = Vec3(2.128000,0.302000,2.833000);
+positions[600] = Vec3(0.300000,1.744000,0.027000);
+positions[601] = Vec3(1.878000,2.986000,2.060000);
+positions[602] = Vec3(1.528000,0.233000,2.045000);
+positions[603] = Vec3(1.146000,1.817000,2.067000);
+positions[604] = Vec3(1.037000,2.746000,0.813000);
+positions[605] = Vec3(0.524000,0.610000,1.566000);
+positions[606] = Vec3(0.945000,2.964000,0.503000);
+positions[607] = Vec3(1.788000,2.565000,0.965000);
+positions[608] = Vec3(0.471000,2.510000,0.491000);
+positions[609] = Vec3(0.512000,2.043000,1.371000);
+positions[610] = Vec3(2.316000,2.423000,1.494000);
+positions[611] = Vec3(1.575000,2.394000,2.953000);
+positions[612] = Vec3(2.845000,2.869000,0.985000);
+positions[613] = Vec3(1.016000,2.335000,1.003000);
+positions[614] = Vec3(0.998000,2.830000,1.879000);
+positions[615] = Vec3(0.624000,2.508000,0.075000);
+positions[616] = Vec3(1.362000,2.808000,2.069000);
+positions[617] = Vec3(1.747000,0.068000,0.810000);
+positions[618] = Vec3(1.768000,2.355000,0.661000);
+positions[619] = Vec3(1.535000,2.410000,2.085000);
+positions[620] = Vec3(0.844000,2.004000,1.646000);
+positions[621] = Vec3(1.124000,0.280000,0.649000);
+positions[622] = Vec3(0.689000,2.170000,0.648000);
+positions[623] = Vec3(0.849000,2.666000,1.175000);
+positions[624] = Vec3(2.975000,1.963000,1.308000);
+positions[625] = Vec3(1.074000,2.082000,0.714000);
+positions[626] = Vec3(1.284000,2.651000,1.110000);
+positions[627] = Vec3(1.669000,0.205000,0.180000);
+positions[628] = Vec3(1.716000,0.047000,1.253000);
+positions[629] = Vec3(0.501000,2.241000,1.043000);
+positions[630] = Vec3(1.038000,1.833000,0.305000);
+positions[631] = Vec3(0.646000,2.431000,1.424000);
+positions[632] = Vec3(1.383000,2.059000,2.230000);
+positions[633] = Vec3(0.370000,2.566000,1.192000);
+positions[634] = Vec3(1.355000,2.006000,0.120000);
+positions[635] = Vec3(2.113000,0.075000,0.589000);
+positions[636] = Vec3(1.850000,0.448000,1.890000);
+positions[637] = Vec3(1.215000,2.704000,0.405000);
+positions[638] = Vec3(0.575000,2.997000,1.798000);
+positions[639] = Vec3(0.967000,2.586000,2.603000);
+positions[640] = Vec3(0.276000,1.669000,1.430000);
+positions[641] = Vec3(1.483000,2.284000,1.128000);
+positions[642] = Vec3(0.983000,3.003000,1.099000);
+positions[643] = Vec3(0.539000,2.222000,1.720000);
+positions[644] = Vec3(0.648000,2.826000,2.751000);
+positions[645] = Vec3(0.803000,1.994000,0.993000);
+positions[646] = Vec3(0.451000,0.216000,1.438000);
+positions[647] = Vec3(1.604000,2.512000,0.334000);
+positions[648] = Vec3(1.980000,2.022000,0.588000);
+positions[649] = Vec3(1.843000,2.834000,1.544000);
+positions[650] = Vec3(1.835000,3.005000,2.858000);
+positions[651] = Vec3(0.679000,2.499000,0.838000);
+positions[652] = Vec3(0.012000,2.637000,1.524000);
+positions[653] = Vec3(0.314000,2.065000,0.602000);
+positions[654] = Vec3(1.157000,0.004000,0.173000);
+positions[655] = Vec3(0.736000,1.705000,1.382000);
+positions[656] = Vec3(1.511000,2.736000,0.690000);
+positions[657] = Vec3(1.330000,2.541000,1.735000);
+positions[658] = Vec3(0.744000,0.170000,0.785000);
+positions[659] = Vec3(2.593000,2.794000,0.703000);
+positions[660] = Vec3(0.275000,1.872000,1.043000);
+positions[661] = Vec3(1.624000,2.608000,1.341000);
+positions[662] = Vec3(1.382000,0.122000,2.855000);
+positions[663] = Vec3(1.326000,2.434000,0.783000);
+positions[664] = Vec3(0.117000,0.116000,1.254000);
+positions[665] = Vec3(1.045000,2.970000,2.748000);
+positions[666] = Vec3(1.341000,2.692000,2.799000);
+positions[667] = Vec3(1.797000,2.586000,2.709000);
+positions[668] = Vec3(0.890000,2.484000,1.716000);
+positions[669] = Vec3(2.373000,2.558000,1.889000);
+positions[670] = Vec3(1.566000,2.323000,2.574000);
+positions[671] = Vec3(1.257000,2.280000,0.399000);
+positions[672] = Vec3(0.679000,2.130000,2.434000);
+positions[673] = Vec3(2.016000,2.334000,2.462000);
+positions[674] = Vec3(1.077000,2.213000,2.416000);
+positions[675] = Vec3(0.581000,1.950000,2.081000);
+positions[676] = Vec3(0.805000,2.315000,2.810000);
+positions[677] = Vec3(0.844000,1.787000,2.322000);
+positions[678] = Vec3(0.980000,2.205000,0.129000);
+positions[679] = Vec3(2.468000,0.603000,2.740000);
+positions[680] = Vec3(2.366000,2.403000,2.299000);
+positions[681] = Vec3(0.337000,2.487000,2.329000);
+positions[682] = Vec3(2.007000,2.793000,2.452000);
+positions[683] = Vec3(1.072000,2.571000,0.063000);
+positions[684] = Vec3(1.217000,2.283000,2.806000);
+positions[685] = Vec3(0.459000,2.477000,2.728000);
+positions[686] = Vec3(0.958000,1.975000,2.710000);
+positions[687] = Vec3(0.914000,2.111000,2.052000);
+positions[688] = Vec3(0.768000,2.958000,0.075000);
+positions[689] = Vec3(0.474000,1.805000,2.533000);
+positions[690] = Vec3(1.313000,2.552000,2.395000);
+positions[691] = Vec3(1.853000,2.014000,2.229000);
+positions[692] = Vec3(2.405000,2.230000,2.658000);
+positions[693] = Vec3(0.727000,1.781000,0.016000);
+positions[694] = Vec3(0.974000,2.791000,2.271000);
+positions[695] = Vec3(0.438000,0.096000,2.457000);
+positions[696] = Vec3(0.652000,2.392000,2.064000);
+positions[697] = Vec3(1.972000,2.209000,2.834000);
+positions[698] = Vec3(0.333000,0.141000,2.088000);
+positions[699] = Vec3(1.813000,1.952000,0.063000);
+positions[700] = Vec3(0.166000,2.838000,1.877000);
+positions[701] = Vec3(1.772000,0.487000,0.951000);
+positions[702] = Vec3(1.924000,1.404000,1.434000);
+positions[703] = Vec3(2.734000,0.348000,1.712000);
+positions[704] = Vec3(2.874000,0.729000,1.811000);
+positions[705] = Vec3(1.841000,0.877000,1.137000);
+positions[706] = Vec3(2.327000,1.491000,1.768000);
+positions[707] = Vec3(1.916000,1.483000,1.057000);
+positions[708] = Vec3(2.783000,0.850000,0.745000);
+positions[709] = Vec3(1.829000,1.526000,0.085000);
+positions[710] = Vec3(2.426000,1.082000,0.652000);
+positions[711] = Vec3(1.645000,1.241000,1.217000);
+positions[712] = Vec3(2.286000,0.725000,0.084000);
+positions[713] = Vec3(2.755000,0.691000,1.421000);
+positions[714] = Vec3(2.651000,0.591000,1.023000);
+positions[715] = Vec3(2.040000,0.863000,0.442000);
+positions[716] = Vec3(0.035000,0.109000,2.497000);
+positions[717] = Vec3(0.127000,1.410000,0.572000);
+positions[718] = Vec3(2.174000,0.357000,0.307000);
+positions[719] = Vec3(2.705000,1.508000,0.758000);
+positions[720] = Vec3(2.223000,1.407000,2.913000);
+positions[721] = Vec3(2.528000,1.722000,1.088000);
+positions[722] = Vec3(2.860000,0.345000,0.198000);
+positions[723] = Vec3(2.580000,1.789000,1.479000);
+positions[724] = Vec3(2.779000,0.295000,1.295000);
+positions[725] = Vec3(0.097000,0.434000,2.826000);
+positions[726] = Vec3(2.952000,1.654000,1.091000);
+positions[727] = Vec3(0.119000,1.878000,0.343000);
+positions[728] = Vec3(1.718000,1.173000,0.327000);
+positions[729] = Vec3(2.833000,0.016000,0.527000);
+positions[730] = Vec3(2.085000,1.779000,2.888000);
+positions[731] = Vec3(2.754000,2.952000,1.485000);
+positions[732] = Vec3(2.826000,0.935000,1.162000);
+positions[733] = Vec3(1.564000,1.585000,1.615000);
+positions[734] = Vec3(2.132000,0.645000,1.093000);
+positions[735] = Vec3(2.294000,1.490000,1.350000);
+positions[736] = Vec3(0.081000,0.490000,1.479000);
+positions[737] = Vec3(2.118000,1.165000,1.642000);
+positions[738] = Vec3(2.141000,0.121000,1.390000);
+positions[739] = Vec3(2.385000,0.389000,1.196000);
+positions[740] = Vec3(0.049000,0.166000,0.817000);
+positions[741] = Vec3(1.993000,0.806000,1.814000);
+positions[742] = Vec3(0.006000,1.450000,0.171000);
+positions[743] = Vec3(2.297000,0.428000,0.764000);
+positions[744] = Vec3(2.851000,0.469000,2.114000);
+positions[745] = Vec3(1.814000,1.957000,0.945000);
+positions[746] = Vec3(0.386000,0.327000,0.902000);
+positions[747] = Vec3(2.452000,1.070000,1.807000);
+positions[748] = Vec3(2.309000,1.537000,2.159000);
+positions[749] = Vec3(2.712000,1.497000,2.007000);
+positions[750] = Vec3(1.727000,0.924000,1.503000);
+positions[751] = Vec3(0.861000,0.801000,0.344000);
+positions[752] = Vec3(1.740000,1.245000,0.819000);
+positions[753] = Vec3(0.117000,0.042000,0.197000);
+positions[754] = Vec3(2.557000,0.996000,0.317000);
+positions[755] = Vec3(2.228000,1.588000,2.548000);
+positions[756] = Vec3(2.849000,1.557000,2.708000);
+positions[757] = Vec3(0.152000,1.107000,0.219000);
+positions[758] = Vec3(2.460000,1.318000,1.002000);
+positions[759] = Vec3(2.405000,1.436000,0.528000);
+positions[760] = Vec3(2.135000,1.179000,2.046000);
+positions[761] = Vec3(1.726000,0.588000,0.286000);
+positions[762] = Vec3(2.831000,1.053000,1.538000);
+positions[763] = Vec3(1.932000,1.556000,1.833000);
+positions[764] = Vec3(2.423000,0.900000,1.064000);
+positions[765] = Vec3(3.001000,1.807000,0.709000);
+positions[766] = Vec3(0.578000,1.095000,0.223000);
+positions[767] = Vec3(2.215000,1.160000,0.252000);
+positions[768] = Vec3(2.050000,0.921000,0.835000);
+positions[769] = Vec3(2.080000,1.682000,0.738000);
+positions[770] = Vec3(2.851000,1.753000,0.027000);
+positions[771] = Vec3(0.203000,0.509000,0.202000);
+positions[772] = Vec3(1.967000,1.018000,0.018000);
+positions[773] = Vec3(1.869000,0.878000,2.472000);
+positions[774] = Vec3(1.917000,0.228000,2.507000);
+positions[775] = Vec3(0.316000,0.795000,2.991000);
+positions[776] = Vec3(2.175000,1.229000,2.472000);
+positions[777] = Vec3(2.405000,1.062000,2.931000);
+positions[778] = Vec3(2.501000,0.511000,2.369000);
+positions[779] = Vec3(2.641000,0.819000,2.141000);
+positions[780] = Vec3(0.649000,1.384000,3.006000);
+positions[781] = Vec3(1.012000,0.329000,2.963000);
+positions[782] = Vec3(2.755000,0.350000,2.718000);
+positions[783] = Vec3(2.315000,0.153000,2.454000);
+positions[784] = Vec3(2.583000,1.696000,2.389000);
+positions[785] = Vec3(0.439000,2.593000,1.776000);
+positions[786] = Vec3(2.630000,1.390000,0.116000);
+positions[787] = Vec3(2.854000,0.669000,2.478000);
+positions[788] = Vec3(2.551000,1.342000,2.621000);
+positions[789] = Vec3(2.533000,2.734000,2.987000);
+positions[790] = Vec3(2.772000,2.446000,2.875000);
+positions[791] = Vec3(2.817000,1.051000,2.498000);
+positions[792] = Vec3(2.688000,1.404000,1.621000);
+positions[793] = Vec3(0.083000,2.737000,2.775000);
+positions[794] = Vec3(2.514000,0.322000,2.041000);
+positions[795] = Vec3(2.470000,0.900000,2.504000);
+positions[796] = Vec3(2.790000,0.444000,0.624000);
+positions[797] = Vec3(0.040000,0.840000,2.202000);
+positions[798] = Vec3(0.530000,1.067000,2.764000);
+positions[799] = Vec3(0.191000,1.385000,2.541000);
+positions[800] = Vec3(2.465000,0.363000,0.051000);
+positions[801] = Vec3(1.850000,1.902000,2.592000);
+positions[802] = Vec3(1.432000,0.306000,2.449000);
+positions[803] = Vec3(2.259000,0.489000,1.753000);
+positions[804] = Vec3(2.803000,1.118000,1.956000);
+positions[805] = Vec3(2.426000,0.147000,1.636000);
+positions[806] = Vec3(2.880000,1.846000,2.133000);
+positions[807] = Vec3(2.862000,2.110000,1.867000);
+positions[808] = Vec3(0.424000,1.184000,2.299000);
+positions[809] = Vec3(2.518000,1.218000,2.228000);
+positions[810] = Vec3(2.153000,0.834000,1.468000);
+positions[811] = Vec3(0.105000,1.397000,2.088000);
+positions[812] = Vec3(2.579000,0.601000,0.316000);
+positions[813] = Vec3(2.594000,2.106000,2.968000);
+positions[814] = Vec3(0.448000,1.435000,1.783000);
+positions[815] = Vec3(2.125000,0.299000,2.132000);
+positions[816] = Vec3(2.849000,1.402000,2.356000);
+positions[817] = Vec3(2.956000,0.091000,2.078000);
+positions[818] = Vec3(0.156000,1.696000,2.357000);
+positions[819] = Vec3(1.566000,2.211000,1.557000);
+positions[820] = Vec3(2.047000,0.194000,0.985000);
+positions[821] = Vec3(1.947000,2.680000,0.488000);
+positions[822] = Vec3(2.343000,2.796000,1.447000);
+positions[823] = Vec3(2.006000,2.332000,0.265000);
+positions[824] = Vec3(2.396000,1.834000,0.546000);
+positions[825] = Vec3(2.538000,2.059000,2.207000);
+positions[826] = Vec3(0.110000,2.360000,0.447000);
+positions[827] = Vec3(2.198000,2.448000,1.136000);
+positions[828] = Vec3(2.420000,2.121000,1.271000);
+positions[829] = Vec3(0.422000,2.192000,0.260000);
+positions[830] = Vec3(2.145000,2.767000,2.839000);
+positions[831] = Vec3(2.434000,2.398000,0.421000);
+positions[832] = Vec3(2.489000,2.175000,1.718000);
+positions[833] = Vec3(2.870000,2.527000,0.814000);
+positions[834] = Vec3(2.741000,2.016000,0.337000);
+positions[835] = Vec3(1.997000,2.574000,2.107000);
+positions[836] = Vec3(0.002000,2.128000,0.932000);
+positions[837] = Vec3(2.787000,2.375000,0.234000);
+positions[838] = Vec3(2.235000,1.852000,1.620000);
+positions[839] = Vec3(2.782000,1.642000,0.422000);
+positions[840] = Vec3(2.915000,1.760000,1.699000);
+positions[841] = Vec3(2.047000,2.178000,1.549000);
+positions[842] = Vec3(1.808000,1.878000,1.556000);
+positions[843] = Vec3(2.224000,2.043000,0.913000);
+positions[844] = Vec3(2.619000,2.611000,1.237000);
+positions[845] = Vec3(2.916000,2.726000,0.168000);
+positions[846] = Vec3(2.021000,2.833000,1.176000);
+positions[847] = Vec3(2.967000,2.308000,2.258000);
+positions[848] = Vec3(2.778000,2.270000,1.477000);
+positions[849] = Vec3(2.121000,1.834000,2.002000);
+positions[850] = Vec3(2.097000,2.752000,0.808000);
+positions[851] = Vec3(1.897000,0.566000,1.501000);
+positions[852] = Vec3(0.359000,2.802000,0.036000);
+positions[853] = Vec3(2.966000,2.454000,1.186000);
+positions[854] = Vec3(2.461000,2.964000,1.132000);
+positions[855] = Vec3(2.093000,1.821000,1.243000);
+positions[856] = Vec3(1.706000,2.659000,1.841000);
+positions[857] = Vec3(2.074000,1.709000,0.342000);
+positions[858] = Vec3(2.137000,2.894000,1.813000);
+positions[859] = Vec3(0.223000,2.293000,1.417000);
+positions[860] = Vec3(2.637000,0.007000,0.197000);
+positions[861] = Vec3(1.416000,0.050000,0.483000);
+positions[862] = Vec3(1.845000,2.250000,1.251000);
+positions[863] = Vec3(2.906000,0.034000,2.896000);
+positions[864] = Vec3(2.481000,0.204000,0.474000);
+positions[865] = Vec3(2.234000,2.051000,0.158000);
+positions[866] = Vec3(0.185000,2.453000,0.055000);
+positions[867] = Vec3(2.509000,0.048000,2.786000);
+positions[868] = Vec3(2.202000,2.206000,2.027000);
+positions[869] = Vec3(0.061000,2.367000,2.656000);
+positions[870] = Vec3(3.003000,2.755000,2.241000);
+positions[871] = Vec3(0.297000,2.131000,2.463000);
+positions[872] = Vec3(1.553000,0.429000,1.573000);
+positions[873] = Vec3(2.506000,1.832000,1.911000);
+positions[874] = Vec3(2.472000,1.814000,2.759000);
+positions[875] = Vec3(1.922000,1.563000,2.278000);
+positions[876] = Vec3(2.623000,2.666000,2.169000);
+positions[877] = Vec3(0.120000,1.834000,2.723000);
+positions[878] = Vec3(0.294000,0.103000,2.826000);
+positions[879] = Vec3(2.364000,2.821000,0.417000);
+positions[880] = Vec3(2.446000,1.734000,0.153000);
+positions[881] = Vec3(2.777000,2.037000,2.565000);
+positions[882] = Vec3(2.837000,2.477000,1.924000);
+positions[883] = Vec3(2.221000,1.961000,2.443000);
+positions[884] = Vec3(2.284000,2.895000,2.157000);
+positions[885] = Vec3(2.728000,2.880000,1.861000);
+positions[886] = Vec3(0.454000,2.080000,2.868000);
+positions[887] = Vec3(2.430000,2.790000,2.524000);
+positions[888] = Vec3(1.808000,2.213000,1.899000);
+positions[889] = Vec3(2.666000,0.053000,2.309000);
+positions[890] = Vec3(2.290000,2.408000,2.995000);
+positions[891] = Vec3(2.646000,2.592000,1.625000);
+positions[892] = Vec3(2.750000,2.508000,2.489000);
+positions[893] = Vec3(0.211000,1.753000,1.939000);

platforms/cuda/src/CudaIntegrationUtilities.cpp

@@ -123,15 +123,21 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
 
     // Record the set of constraints and how many constraints each atom is involved in.
 
-    int numConstraints = system.getNumConstraints();
-    vector<int> atom1(numConstraints);
-    vector<int> atom2(numConstraints);
-    vector<double> distance(numConstraints);
+    vector<int> atom1;
+    vector<int> atom2;
+    vector<double> distance;
     vector<int> constraintCount(context.getNumAtoms(), 0);
-    for (int i = 0; i < numConstraints; i++) {
-        system.getConstraintParameters(i, atom1[i], atom2[i], distance[i]);
-        constraintCount[atom1[i]]++;
-        constraintCount[atom2[i]]++;
+    for (int i = 0; i < system.getNumConstraints(); i++) {
+        int p1, p2;
+        double d;
+        system.getConstraintParameters(i, p1, p2, d);
+        if (system.getParticleMass(p1) != 0 || system.getParticleMass(p2) != 0) {
+            atom1.push_back(p1);
+            atom2.push_back(p2);
+            distance.push_back(d);
+            constraintCount[p1]++;
+            constraintCount[p2]++;
+        }
     }
 
     // Identify clusters of three atoms that can be treated with SETTLE.  First, for every

platforms/cuda/tests/TestCudaBrownianIntegrator.cpp

@@ -172,6 +172,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    BrownianIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -237,6 +268,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/cuda/tests/TestCudaCustomIntegrator.cpp

@@ -233,6 +233,43 @@ void testVelocityConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    CustomIntegrator integrator(0.002);
+    integrator.addPerDofVariable("oldx", 0);
+    integrator.addComputePerDof("v", "v+dt*f/m");
+    integrator.addComputePerDof("oldx", "x");
+    integrator.addComputePerDof("x", "x+dt*v");
+    integrator.addConstrainPositions();
+    integrator.addComputePerDof("v", "(x-oldx)/dt");
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities | State::Positions);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 /**
  * Test an integrator with an AndersenThermostat to see if updateContextState()
  * is being handled correctly.
@@ -724,6 +761,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testConstraints();
         testVelocityConstraints();
+        testConstrainedMasslessParticles();
         testWithThermostat();
         testMonteCarlo();
         testSum();

platforms/cuda/tests/TestCudaLangevinIntegrator.cpp

@@ -177,6 +177,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    LangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -241,6 +272,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/cuda/tests/TestCudaVariableLangevinIntegrator.cpp

@@ -172,6 +172,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -295,6 +326,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
         testArgonBox();
     }