Checked in explicit solvent implementation for Cuda platform (cutoffs,...

Checked in explicit solvent implementation for Cuda platform (cutoffs, periodic boundary conditions, SETTLE)

openmmapi/include/GBSAOBCForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -40,6 +40,9 @@ namespace OpenMM {
 
 /**
  * This class implements an implicit solvation force using the GBSA-OBC model.
+ * <p>
+ * If the System also contains a NonbondedForce, this force will use the cutoffs
+ * and periodic boundary conditions specified in it.
  */
 
 class OPENMM_EXPORT GBSAOBCForce : public Force {

openmmapi/include/System.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -125,6 +125,14 @@ public:
     int getNumForces() const {
         return forces.size();
     }
+    /**
+     * Get a reference to one of the Forces in this System.
+     * 
+     * @param index  the index of the Force to get
+     */
+    const Force& getForce(int index) const {
+        return *forces[index];
+    }
     /**
      * Get a reference to one of the Forces in this System.
      *

openmmapi/include/internal/ForceImpl.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include <map>
+#include <string>
 #include <vector>
 
 namespace OpenMM {

platforms/cuda/cudpp/license.txt

+CUDA Data-Parallel Primitives Library (CUDPP) is the proprietary
+property of The Regents of the University of California ("The
+Regents") and NVIDIA Corporation ("NVIDIA").
+
+Copyright (c) 2007-2008 The Regents of the University of California, Davis
+campus and NVIDIA Corporation. All Rights Reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+  notice, this list of conditions and the following disclaimer. 
+
+- Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions and the following disclaimer in the
+  documentation and/or other materials provided with the distribution. 
+
+- Neither the name of The Regents, NVIDIA, nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+The end-user understands that the program was developed for research
+purposes and is advised not to rely exclusively on the program for any
+reason. No other rights or permissions are granted beyond what is stated 
+explictely herein.
+
+THE SOFTWARE PROVIDED IS ON AN "AS IS" BASIS, AND THE REGENTS, NVIDIA
+AND CONTRIBUTORS HAVE NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
+UPDATES, ENHANCEMENTS, OR MODIFICATIONS. THE REGENTS, NVIDIA AND
+CONTRIBUTORS SPECIFICALLY DISCLAIM ANY EXPRESS OR IMPLIED WARRANTIES,
+INCLUDING BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE REGENTS, NVIDIA OR CONTRIBUTORS BE LIABLE TO ANY
+PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, EXEMPLARY OR
+CONSEQUENTIAL DAMAGES, INCLUDING BUT NOT LIMITED TO PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES, LOSS OF USE, DATA OR PROFITS, OR
+BUSINESS INTERRUPTION, HOWEVER CAUSED AND UNDER ANY THEORY OF
+LIABILITY WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE AND ITS DOCUMENTATION, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
+
+If you do not agree to these terms, do not download or use the
+software.  This license may be modified only in a writing signed by
+authorized signatory of all parties.  For The Regents contact
+copyright@ucdavis.edu.
+
+Relating to funding received by the Regents-
+
+Acknowledgment: This material is based in part upon work supported by
+the Department of Energy under Award Numbers DE-FG02-04ER25609 and
+DE-FC02-06ER25777.
+
+Disclaimer: This report was prepared as an account of work sponsored
+by an agency of the United States Government. Neither the United
+States Government nor any agency thereof, nor any of their employees,
+makes any warranty, express or implied, or assumes any legal liability
+or responsibility for the accuracy, completeness, or usefulness of any
+information, apparatus, product, or process disclosed, or represents
+that its use would not infringe privately owned rights. Reference
+herein to any specific commercial product, process, or service by
+trade name, trademark, manufacturer, or otherwise does not necessarily
+constitute or imply its endorsement, recommendation, or favoring by
+the United States Government or any agency thereof. The views and
+opinions of authors expressed herein do not necessarily state or
+reflect those of the United States Government or any agency hereof.

platforms/cuda/cudpp/license_nv.txt

+This portion of the CUDA Data-Parallel Primitives Library (CUDPP) is the 
+proprietary property of NVIDIA Corporation ("NVIDIA").
+
+Copyright (c) 2007-2008 NVIDIA Corporation. All Rights Reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+- Redistributions of source code must retain the above copyright
+  notice, this list of conditions and the following disclaimer. 
+
+- Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions and the following disclaimer in the
+  documentation and/or other materials provided with the distribution. 
+
+- Neither the name of NVIDIA, nor the names of its contributors may be used to 
+  endorse or promote products derived from this software without specific prior 
+  written permission.
+
+The end-user understands that the program was developed for research
+purposes and is advised not to rely exclusively on the program for any
+reason. No other rights or permissions are granted beyond what is stated 
+explictely herein.
+
+THE SOFTWARE PROVIDED IS ON AN "AS IS" BASIS, AND NVIDIA
+AND CONTRIBUTORS HAVE NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
+UPDATES, ENHANCEMENTS, OR MODIFICATIONS. NVIDIA AND CONTRIBUTORS 
+SPECIFICALLY DISCLAIM ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING 
+BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, 
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT 
+SHALL NVIDIA OR CONTRIBUTORS BE LIABLE TO ANY PARTY FOR DIRECT, 
+INDIRECT, SPECIAL, INCIDENTAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES, 
+INCLUDING BUT NOT LIMITED TO PROCUREMENT OF SUBSTITUTE GOODS OR 
+SERVICES, LOSS OF USE, DATA OR PROFITS, OR BUSINESS INTERRUPTION, 
+HOWEVER CAUSED AND UNDER ANY THEORY OF LIABILITY WHETHER IN 
+CONTRACT, STRICT LIABILITY OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE AND ITS 
+DOCUMENTATION, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+If you do not agree to these terms, do not download or use the
+software.  This license may be modified only in a writing signed by
+authorized signatory of all parties.  
\ No newline at end of file

platforms/cuda/include/CudaPlatform.h

@@ -65,14 +65,16 @@ private:
 
 class CudaPlatform::PlatformData {
 public:
-    PlatformData(_gpuContext* gpu) : gpu(gpu), removeCM(false), useOBC(false), hasBonds(false), hasAngles(false),
-            hasPeriodicTorsions(false), hasRB(false), hasNonbonded(false), primaryKernel(NULL) {
+    PlatformData(_gpuContext* gpu) : gpu(gpu), removeCM(false), nonbondedMethod(0), hasBonds(false), hasAngles(false),
+            hasPeriodicTorsions(false), hasRB(false), hasNonbonded(false), primaryKernel(NULL), stepCount(0) {
     }
     _gpuContext* gpu;
     KernelImpl* primaryKernel;
-    bool removeCM, useOBC;
+    bool removeCM;
     bool hasBonds, hasAngles, hasPeriodicTorsions, hasRB, hasNonbonded;
+    int nonbondedMethod;
     int cmMotionFrequency;
+    int stepCount;
 };
 
 } // namespace OpenMM

platforms/cuda/sharedTarget/CMakeLists.txt

@@ -5,19 +5,18 @@ INCLUDE(${FINDCUDA_DIR}/FindCuda.cmake)
 INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
 LINK_DIRECTORIES(${CUDA_TARGET_LINK})
 FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
-    FILE(GLOB src_files ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src/*.cu  ${CMAKE_SOURCE_DIR}/platforms/cuda/src/*/*.cu)
+    FILE(GLOB src_files ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src/*.cu  ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src/*/*.cu)
     FOREACH(file ${src_files})
         FILE(RELATIVE_PATH file ${CMAKE_SOURCE_DIR}/platforms/cuda ${file})
         SET(SOURCE_FILES         ${SOURCE_FILES}         ${file})   #append
     ENDFOREACH(file)
-    CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/../${subdir}/include)
+    CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/include)
+    CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src)
 ENDFOREACH(subdir)
 
-CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/../src)
-
 CUDA_ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
 
-TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}_d optimized ${OPENMM_LIBRARY_NAME})
+TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}_d optimized ${OPENMM_LIBRARY_NAME} cudpp cutil)
 SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY")
 
 INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${SHARED_TARGET})

platforms/cuda/src/CudaKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -47,18 +47,24 @@ using namespace std;
 
 static void calcForces(OpenMMContextImpl& context, CudaPlatform::PlatformData& data) {
     _gpuContext* gpu = data.gpu;
-    if (data.useOBC) {
+    if (data.stepCount%100 == 0)
+        gpuReorderAtoms(gpu);
+    data.stepCount++;
+    kClearForces(gpu);
+    if (gpu->bIncludeGBSA) {
         gpu->bRecalculateBornRadii = true;
         kCalculateCDLJObcGbsaForces1(gpu);
         kReduceObcGbsaBornForces(gpu);
         kCalculateObcGbsaForces2(gpu);
     }
-    else {
-        kClearForces(gpu);
+    else if (data.hasNonbonded) {
         kCalculateCDLJForces(gpu);
     }
     kCalculateLocalForces(gpu);
+    if (gpu->bIncludeGBSA)
         kReduceBornSumAndForces(gpu);
+    else
+        kReduceForces(gpu);
 }
 
 static double calcEnergy(OpenMMContextImpl& context, System& system) {
@@ -259,7 +265,19 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             exclusionList[i] = vector<int>(exclusions[i].begin(), exclusions[i].end());
             exclusionList[i].push_back(i);
         }
-        gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList);
+        CudaNonbondedMethod method = NO_CUTOFF;
+        if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
+            gpuSetNonbondedCutoff(gpu, force.getCutoffDistance(), 78.3f);
+            method = CUTOFF;
+        }
+        if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
+            Vec3 boxVectors[3];
+            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+            gpuSetPeriodicBoxSize(gpu, boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
+            method = PERIODIC;
+        }
+        data.nonbondedMethod = method;
+        gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
     }
 
     // Initialize 1-4 nonbonded interactions.
@@ -312,7 +330,6 @@ void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCF
         scale[i] = (float) scalingFactor;
     }
     gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), particle, radius, scale);
-    data.useOBC = true;
 }
 
 void CudaCalcGBSAOBCForceKernel::executeForces(OpenMMContextImpl& context) {
@@ -361,7 +378,7 @@ static void initializeIntegration(const System& system, CudaPlatform::PlatformDa
         gpuSetRbDihedralParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>(),
                 vector<float>(), vector<float>(), vector<float>(), vector<float>());
     if (!data.hasNonbonded) {
-        gpuSetCoulombParameters(gpu, 138.935485f, vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<char>(), vector<vector<int> >());
+        gpuSetCoulombParameters(gpu, 138.935485f, vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<char>(), vector<vector<int> >(), NO_CUTOFF);
         gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<float>());
     }
     
@@ -401,6 +418,7 @@ void CudaIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, const Ve
     }
     kVerletUpdatePart1(gpu);
     kApplyFirstShake(gpu);
+    kApplyFirstSettle(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)
@@ -444,6 +462,7 @@ void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const
     }
     kUpdatePart1(gpu);
     kApplyFirstShake(gpu);
+    kApplyFirstSettle(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)
@@ -451,6 +470,7 @@ void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const
     }
     kUpdatePart2(gpu);
     kApplySecondShake(gpu);
+    kApplySecondSettle(gpu);
 }
 
 CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
@@ -479,6 +499,7 @@ void CudaIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const
     }
     kBrownianUpdatePart1(gpu);
     kApplyFirstShake(gpu);
+    kApplyFirstSettle(gpu);
     if (data.removeCM) {
         int step = (int) (context.getTime()/stepSize);
         if (step%data.cmMotionFrequency == 0)

platforms/cuda/staticTarget/CMakeLists.txt

@@ -5,16 +5,15 @@ INCLUDE(${FINDCUDA_DIR}/FindCuda.cmake)
 INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
 LINK_DIRECTORIES(${CUDA_TARGET_LINK})
 FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
-    FILE(GLOB src_files ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src/*.cu  ${CMAKE_SOURCE_DIR}/platforms/cuda/src/*/*.cu)
+    FILE(GLOB src_files ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src/*.cu  ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src/*/*.cu)
     FOREACH(file ${src_files})
         FILE(RELATIVE_PATH file ${CMAKE_SOURCE_DIR}/platforms/cuda ${file})
         SET(SOURCE_FILES         ${SOURCE_FILES}         ${file})   #append
     ENDFOREACH(file)
-    CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/../${subdir}/include)
+    CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/include)
+    CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/${subdir}/src)
 ENDFOREACH(subdir)
 
-CUDA_INCLUDE_DIRECTORIES(BEFORE ${CMAKE_SOURCE_DIR}/platforms/cuda/../src)
-
 CUDA_ADD_LIBRARY(${STATIC_TARGET} STATIC ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
 
 # required for getting OPENMM_EXPORT to be set correctly in 'class OPENMM_EXPORT CudaStreamFactory', ...
@@ -23,6 +22,6 @@ CUDA_ADD_LIBRARY(${STATIC_TARGET} STATIC ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES}
 SET( CUDA_STATIC_COMPILE_FLAG "-DOPENMM_BUILDING_STATIC_LIBRARY -DOPENMM_USE_STATIC_LIBRARIES")
 SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES COMPILE_FLAGS ${CUDA_STATIC_COMPILE_FLAG})
 
-TARGET_LINK_LIBRARIES(${STATIC_TARGET} debug ${OPENMM_LIBRARY_NAME}_static_d optimized ${OPENMM_LIBRARY_NAME}_static)
+TARGET_LINK_LIBRARIES(${STATIC_TARGET} debug ${OPENMM_LIBRARY_NAME}_static_d optimized ${OPENMM_LIBRARY_NAME}_static cudpp cutil)
 
 INSTALL_TARGETS(/lib RUNTIME_DIRECTORY /lib ${STATIC_TARGET})

platforms/cuda/tests/TestCudaGBSAOBCForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -36,6 +36,7 @@
 #include "../../../tests/AssertionUtilities.h"
 #include "OpenMMContext.h"
 #include "CudaPlatform.h"
+#include "ReferencePlatform.h"
 #include "GBSAOBCForce.h"
 #include "System.h"
 #include "LangevinIntegrator.h"
@@ -55,12 +56,12 @@ void testSingleParticle() {
     System system(1, 0);
     system.setParticleMass(0, 2.0);
     LangevinIntegrator integrator(0, 0.1, 0.01);
-    GBSAOBCForce* forceField = new GBSAOBCForce(1);
-    NonbondedForce* standard = new NonbondedForce(1, 0);
-    forceField->setParticleParameters(0, 0.5, 0.15, 1);
-    standard->setParticleParameters(0, 0.5, 1, 0);
-    system.addForce(forceField);
-    system.addForce(standard);
+    GBSAOBCForce* gbsa = new GBSAOBCForce(1);
+    NonbondedForce* nonbonded = new NonbondedForce(1, 0);
+    gbsa->setParticleParameters(0, 0.5, 0.15, 1);
+    nonbonded->setParticleParameters(0, 0.5, 1, 0);
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(1);
     positions[0] = Vec3(0, 0, 0);
@@ -68,45 +69,120 @@ void testSingleParticle() {
     State state = context.getState(State::Energy);
     double bornRadius = 0.15-0.009; // dielectric offset
     double eps0 = EPSILON0;
-    double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
+    double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
     double extendedRadius = bornRadius+0.14; // probe radius
     double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from?  Just copied it from CpuImplicitSolvent.cpp
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
 }
 
-void testForce() {
-    CudaPlatform platform;
-    const int numParticles = 10;
+void testCutoffAndPeriodic() {
+    CudaPlatform cuda;
+    System system(2, 0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* gbsa = new GBSAOBCForce(2);
+    NonbondedForce* nonbonded = new NonbondedForce(2, 0);
+    gbsa->setParticleParameters(0, -1, 0.15, 1);
+    nonbonded->setParticleParameters(0, -1, 1, 0);
+    gbsa->setParticleParameters(1, 1, 0.15, 1);
+    nonbonded->setParticleParameters(1, 1, 1, 0);
+    const double cutoffDistance = 3.0;
+    const double boxSize = 10.0;
+    nonbonded->setCutoffDistance(cutoffDistance);
+    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+
+    // Calculate the forces for both cutoff and periodic with two different atom positions.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    OpenMMContext context(system, integrator, cuda);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state2 = context.getState(State::Forces);
+    positions[1][0]+= boxSize;
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state3 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state4 = context.getState(State::Forces);
+
+    // All forces should be identical, exception state3 which should be zero.
+
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state2.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state2.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state4.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state4.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[0], Vec3(0, 0, 0), 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
+}
+
+void testForce(int numParticles, NonbondedForce::NonbondedMethod method) {
+    CudaPlatform cuda;
+    ReferencePlatform reference;
     System system(numParticles, 0);
     LangevinIntegrator integrator(0, 0.1, 0.01);
-    GBSAOBCForce* forceField = new GBSAOBCForce(numParticles);
-    NonbondedForce* standard = new NonbondedForce(numParticles, 0);
+    GBSAOBCForce* gbsa = new GBSAOBCForce(numParticles);
+    NonbondedForce* nonbonded = new NonbondedForce(numParticles, 0);
     for (int i = 0; i < numParticles; ++i) {
         double charge = i%2 == 0 ? -1 : 1;
-        forceField->setParticleParameters(i, charge, 0.15, 1);
-        standard->setParticleParameters(i, charge, 1, 0);
+        gbsa->setParticleParameters(i, charge, 0.15, 1);
+        nonbonded->setParticleParameters(i, charge, 1, 0);
     }
-    system.addForce(forceField);
-    system.addForce(standard);
-    OpenMMContext context(system, integrator, platform);
+    nonbonded->setNonbondedMethod(method);
+    nonbonded->setCutoffDistance(3.0);
+    int grid = (int) floor(0.5+pow(numParticles, 1.0/3.0));
+    if (method == NonbondedForce::CutoffPeriodic) {
+        double boxSize = (grid+1)*2.0;
+        nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    }
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
+    OpenMMContext context(system, integrator, cuda);
+    OpenMMContext refContext(system, integrator, reference);
     
-    // Set random positions for all the particles.
+    // Set random (but uniformly distributed) positions for all the particles.
     
     vector<Vec3> positions(numParticles);
     init_gen_rand(0);
+    for (int i = 0; i < grid; i++)
+        for (int j = 0; j < grid; j++)
+            for (int k = 0; k < grid; k++)
+                positions[i*grid*grid+j*grid+k] = Vec3(i*2.0, j*2.0, k*2.0);
     for (int i = 0; i < numParticles; ++i)
-        positions[i] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
+        positions[i] = positions[i] + Vec3(0.5*genrand_real2(), 0.5*genrand_real2(), 0.5*genrand_real2());
     context.setPositions(positions);
+    refContext.setPositions(positions);
     State state = context.getState(State::Forces | State::Energy);
+    State refState = refContext.getState(State::Forces | State::Energy);
 
-    // Take a small step in the direction of the energy gradient.
+    // Make sure the Cuda and Reference platforms agree.
 
     double norm = 0.0;
+    double diff = 0.0;
     for (int i = 0; i < numParticles; ++i) {
         Vec3 f = state.getForces()[i];
         norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+        Vec3 delta = f-refState.getForces()[i];
+        diff += delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
     }
     norm = std::sqrt(norm);
+    diff = std::sqrt(diff);
+    ASSERT_EQUAL_TOL(0.0, diff, 0.001*norm); 
+
+    // Take a small step in the direction of the energy gradient.  (This doesn't work with cutoffs, since the energy
+    // changes discontinuously at the cutoff distance.)
+
+    if (method == NonbondedForce::NoCutoff)
+    {
         const double delta = 1e-2;
         double step = delta/norm;
         for (int i = 0; i < numParticles; ++i) {
@@ -119,13 +195,19 @@ void testForce() {
         // See whether the potential energy changed by the expected amount.
 
         State state2 = context.getState(State::Energy);
-    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 0.02)
+        ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 1e-3*abs(state.getPotentialEnergy()));
+    }
 }
 
 int main() {
     try {
         testSingleParticle();
-        testForce();
+        testCutoffAndPeriodic();
+        for (int i = 2; i < 11; i++) {
+            testForce(i*i*i, NonbondedForce::NoCutoff);
+            testForce(i*i*i, NonbondedForce::CutoffNonPeriodic);
+            testForce(i*i*i, NonbondedForce::CutoffPeriodic);
+        }
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/cuda/tests/TestCudaNonbondedForce.cpp

@@ -36,11 +36,16 @@
 #include "../../../tests/AssertionUtilities.h"
 #include "OpenMMContext.h"
 #include "CudaPlatform.h"
+#include "ReferencePlatform.h"
 #include "HarmonicBondForce.h"
 #include "NonbondedForce.h"
 #include "System.h"
 #include "LangevinIntegrator.h"
+#include "VerletIntegrator.h"
+#include "internal/OpenMMContextImpl.h"
+#include "kernels/gputypes.h"
 #include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
 #include <iostream>
 #include <vector>
 
@@ -211,6 +216,8 @@ void testCutoff14() {
     system.addForce(bonds);
     NonbondedForce* nonbonded = new NonbondedForce(5, 2);
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, 0.0);
+    nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
     const double cutoff = 3.5;
     nonbonded->setCutoffDistance(cutoff);
     system.addForce(nonbonded);
@@ -316,14 +323,219 @@ void testPeriodic() {
     ASSERT_EQUAL_TOL(2*138.935485*(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 = 10.0;
+    const double tol = 1e-3;
+    CudaPlatform cuda;
+    ReferencePlatform reference;
+    System system(numParticles, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce(numParticles, 0);
+    HarmonicBondForce* bonds = new HarmonicBondForce(numMolecules);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            nonbonded->setParticleParameters(2*i, 1.0, 0.2, 0.1);
+            nonbonded->setParticleParameters(2*i+1, 1.0, 0.1, 0.1);
+        }
+        else {
+            nonbonded->setParticleParameters(2*i, 1.0, 0.2, 0.2);
+            nonbonded->setParticleParameters(2*i+1, 1.0, 0.1, 0.2);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(), boxSize*genrand_real2(), boxSize*genrand_real2());
+        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(), genrand_real2(), genrand_real2());
+        velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
+        bonds->setBondParameters(i, 2*i, 2*i+1, 1.0, 0.1);
+    }
+
+    // Try with cutoffs but not periodic boundary conditions, and make sure the Cuda and Reference
+    // platforms agree.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.addForce(nonbonded);
+    system.addForce(bonds);
+    OpenMMContext cudaContext(system, integrator, cuda);
+    OpenMMContext referenceContext(system, integrator, reference);
+    cudaContext.setPositions(positions);
+    cudaContext.setVelocities(velocities);
+    referenceContext.setPositions(positions);
+    referenceContext.setVelocities(velocities);
+    State cudaState = cudaContext.getState(State::Positions | State::Velocities | State::Forces);
+    State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(cudaState.getPositions()[i], referenceState.getPositions()[i], tol);
+        ASSERT_EQUAL_VEC(cudaState.getVelocities()[i], referenceState.getVelocities()[i], tol);
+        ASSERT_EQUAL_VEC(cudaState.getForces()[i], referenceState.getForces()[i], tol);
+    }
+
+    // Now do the same thing with periodic boundary conditions.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    cudaContext.reinitialize();
+    referenceContext.reinitialize();
+    cudaContext.setPositions(positions);
+    cudaContext.setVelocities(velocities);
+    referenceContext.setPositions(positions);
+    referenceContext.setVelocities(velocities);
+    cudaState = cudaContext.getState(State::Positions | State::Velocities | State::Forces);
+    referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(cudaState.getPositions()[i], referenceState.getPositions()[i], tol);
+        ASSERT_EQUAL_VEC(cudaState.getVelocities()[i], referenceState.getVelocities()[i], tol);
+        ASSERT_EQUAL_VEC(cudaState.getForces()[i], referenceState.getForces()[i], tol);
+    }
+}
+
+void testBlockInteractions(bool periodic) {
+    const int blockSize = 32;
+    const int numBlocks = 100;
+    const int numParticles = blockSize*numBlocks;
+    const double cutoff = 1.0;
+    const double boxSize = (periodic ? 5.1 : 1.1);
+    CudaPlatform cuda;
+    System system(numParticles, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce(numParticles, 0);
+    vector<Vec3> positions(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numParticles; i++) {
+        nonbonded->setParticleParameters(i, 1.0, 0.2, 0.2);
+        positions[i] = Vec3(boxSize*(3*genrand_real2()-1), boxSize*(3*genrand_real2()-1), boxSize*(3*genrand_real2()-1));
+    }
+    nonbonded->setNonbondedMethod(periodic ? NonbondedForce::CutoffPeriodic : NonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+    OpenMMContext context(system, integrator, cuda);
+    context.setPositions(positions);
+    State state = context.getState(State::Positions | State::Velocities | State::Forces);
+    OpenMMContextImpl* contextImpl = *reinterpret_cast<OpenMMContextImpl**>(&context);
+    CudaPlatform::PlatformData& data = *static_cast<CudaPlatform::PlatformData*>(contextImpl->getPlatformData());
+    
+    // Verify that the bounds of each block were calculated correctly.
+
+    data.gpu->psPosq4->Download();
+    data.gpu->psGridBoundingBox->Download();
+    data.gpu->psGridCenter->Download();
+    for (int i = 0; i < numBlocks; i++) {
+        float4 gridSize = data.gpu->psGridBoundingBox->_pSysData[i];
+        float4 center = data.gpu->psGridCenter->_pSysData[i];
+        if (periodic) {
+            ASSERT(gridSize.x < 0.5*boxSize);
+            ASSERT(gridSize.y < 0.5*boxSize);
+            ASSERT(gridSize.z < 0.5*boxSize);
+        }
+        float minx = 0.0, maxx = 0.0, miny = 0.0, maxy = 0.0, minz = 0.0, maxz = 0.0, radius = 0.0;
+        for (int j = 0; j < blockSize; j++) {
+            float4 pos = data.gpu->psPosq4->_pSysData[i*blockSize+j];
+            float dx = pos.x-center.x;
+            float dy = pos.y-center.y;
+            float dz = pos.z-center.z;
+            if (periodic) {
+                dx -= floor(0.5+dx/boxSize)*boxSize;
+                dy -= floor(0.5+dy/boxSize)*boxSize;
+                dz -= floor(0.5+dz/boxSize)*boxSize;
+            }
+            ASSERT(abs(dx) < gridSize.x+TOL);
+            ASSERT(abs(dy) < gridSize.y+TOL);
+            ASSERT(abs(dz) < gridSize.z+TOL);
+            minx = min(minx, dx);
+            maxx = max(maxx, dx);
+            miny = min(miny, dy);
+            maxy = max(maxy, dy);
+            minz = min(minz, dz);
+            maxz = max(maxz, dz);
+        }
+        ASSERT_EQUAL_TOL(-minx, gridSize.x, TOL);
+        ASSERT_EQUAL_TOL(maxx, gridSize.x, TOL);
+        ASSERT_EQUAL_TOL(-miny, gridSize.y, TOL);
+        ASSERT_EQUAL_TOL(maxy, gridSize.y, TOL);
+        ASSERT_EQUAL_TOL(-minz, gridSize.z, TOL);
+        ASSERT_EQUAL_TOL(maxz, gridSize.z, TOL);
+    }
+
+    // Verify that interactions were identified correctly.
+
+    data.gpu->psInteractionCount->Download();
+    int numWithInteractions = data.gpu->psInteractionCount->_pSysData[0];
+    vector<bool> hasInteractions(data.gpu->sim.workUnits, false);
+    data.gpu->psInteractingWorkUnit->Download();
+    data.gpu->psInteractionFlag->Download();
+    const unsigned int atoms = data.gpu->sim.paddedNumberOfAtoms;
+    const unsigned int grid = data.gpu->grid;
+    const unsigned int dim = (atoms+(grid-1))/grid;
+    for (int i = 0; i < numWithInteractions; i++) {
+        unsigned int workUnit = data.gpu->psInteractingWorkUnit->_pSysData[i];
+        unsigned int x = (workUnit >> 17);
+        unsigned int y = ((workUnit >> 2) & 0x7fff);
+        int tile = (x > y ? x+y*dim-y*(y+1)/2 : y+x*dim-x*(x+1)/2);
+        hasInteractions[tile] = true;
+
+        // Make sure this tile really should have been flagged based on bounding volumes.
+
+        float4 gridSize1 = data.gpu->psGridBoundingBox->_pSysData[x];
+        float4 gridSize2 = data.gpu->psGridBoundingBox->_pSysData[y];
+        float4 center1 = data.gpu->psGridCenter->_pSysData[x];
+        float4 center2 = data.gpu->psGridCenter->_pSysData[y];
+        float dx = center1.x-center2.x;
+        float dy = center1.y-center2.y;
+        float dz = center1.z-center2.z;
+        if (periodic) {
+            dx -= floor(0.5+dx/boxSize)*boxSize;
+            dy -= floor(0.5+dy/boxSize)*boxSize;
+            dz -= floor(0.5+dz/boxSize)*boxSize;
+        }
+        dx = max(0.0f, abs(dx)-gridSize1.x-gridSize2.x);
+        dy = max(0.0f, abs(dy)-gridSize1.y-gridSize2.y);
+        dz = max(0.0f, abs(dz)-gridSize1.z-gridSize2.z);
+        ASSERT(sqrt(dx*dx+dy*dy+dz*dz) < cutoff+TOL);
+    }
+
+    // Check the tiles that did not have interactions to make sure all atoms are beyond the cutoff.
+
+    data.gpu->psWorkUnit->Download();
+    for (int i = 0; i < hasInteractions.size(); i++)
+        if (!hasInteractions[i]) {
+            unsigned int workUnit = data.gpu->psWorkUnit->_pSysData[i];
+            unsigned int x = (workUnit >> 17);
+            unsigned int y = ((workUnit >> 2) & 0x7fff);
+            for (int atom1 = 0; atom1 < blockSize; ++atom1) {
+                float4 pos1 = data.gpu->psPosq4->_pSysData[x*blockSize+atom1];
+                for (int atom2 = 0; atom2 < blockSize; ++atom2) {
+                    float4 pos2 = data.gpu->psPosq4->_pSysData[y*blockSize+atom2];
+                    float dx = pos1.x-pos2.x;
+                    float dy = pos1.y-pos2.y;
+                    float dz = pos1.z-pos2.z;
+                    if (periodic) {
+                        dx -= floor(0.5+dx/boxSize)*boxSize;
+                        dy -= floor(0.5+dy/boxSize)*boxSize;
+                        dz -= floor(0.5+dz/boxSize)*boxSize;
+                    }
+                    ASSERT(dx*dx+dy*dy+dz*dz > cutoff*cutoff);
+                }
+            }
+        }
+}
+
 int main() {
     try {
         testCoulomb();
         testLJ();
         testExclusionsAnd14();
-//        testCutoff();
-//        testCutoff14();
-//        testPeriodic();
+        testCutoff();
+        testCutoff14();
+        testPeriodic();
+        testLargeSystem();
+        testBlockInteractions(false);
+        testBlockInteractions(true);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/cuda/tests/TestCudaSettle.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-2009 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 Cuda implementation of the SETTLE algorithm.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "OpenMMContext.h"
+#include "CudaPlatform.h"
+#include "NonbondedForce.h"
+#include "System.h"
+#include "LangevinIntegrator.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testConstraints() {
+    const int numMolecules = 10;
+    const int numParticles = numMolecules*3;
+    const int numConstraints = numMolecules*3;
+    const double temp = 100.0;
+    CudaPlatform platform;
+    System system(numParticles, numConstraints);
+    LangevinIntegrator integrator(temp, 2.0, 0.001);
+    integrator.setConstraintTolerance(1e-5);
+    NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
+    for (int i = 0; i < numMolecules; ++i) {
+        system.setParticleMass(i*3, 16.0);
+        system.setParticleMass(i*3+1, 1.0);
+        system.setParticleMass(i*3+2, 1.0);
+        forceField->setParticleParameters(i*3, -0.82, 0.317, 0.65);
+        forceField->setParticleParameters(i*3+1, 0.41, 1.0, 0.0);
+        forceField->setParticleParameters(i*3+2, 0.41, 1.0, 0.0);
+        system.setConstraintParameters(i*3, i*3, i*3+1, 0.1);
+        system.setConstraintParameters(i*3+1, i*3, i*3+2, 0.1);
+        system.setConstraintParameters(i*3+2, i*3+1, i*3+2, 0.163);
+    }
+    system.addForce(forceField);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numMolecules; ++i) {
+        positions[i*3] = Vec3((i%4)*0.4, (i/4)*0.4, 0);
+        positions[i*3+1] = positions[i*3]+Vec3(0.1, 0, 0);
+        positions[i*3+2] = positions[i*3]+Vec3(-0.03333, 0.09428, 0);
+        velocities[i*3] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
+        velocities[i*3+1] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
+        velocities[i*3+2] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
+    }
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+
+    // Simulate it and see whether the constraints remain satisfied.
+
+    for (int i = 0; i < 1000; ++i) {
+        integrator.step(1);
+        State state = context.getState(State::Positions | State::Forces);
+        for (int j = 0; j < numConstraints; ++j) {
+            int particle1, particle2;
+            double distance;
+            system.getConstraintParameters(j, particle1, particle2, distance);
+            Vec3 p1 = state.getPositions()[particle1];
+            Vec3 p2 = state.getPositions()[particle2];
+            double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
+            ASSERT_EQUAL_TOL(distance, dist, 1e-5);
+        }
+    }
+}
+
+int main() {
+    try {
+        testConstraints();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cuda/tests/TestCudaStreams.cpp

@@ -35,7 +35,9 @@
 
 #include "../../../tests/AssertionUtilities.h"
 #include "CudaPlatform.h"
+#include "OpenMMContext.h"
 #include "Stream.h"
+#include "VerletIntegrator.h"
 #include <iostream>
 
 using namespace OpenMM;
@@ -47,7 +49,10 @@ template <class T, int WIDTH>
 void testStream(Stream::DataType type, T scale) {
     const int size = 100;
     CudaPlatform platform;
-    OpenMMContextImpl* impl = 0; // We need to pass an OpenMMContextImpl to createStream, but it will be ignored.
+    System system(size, 0);
+    VerletIntegrator integrator(0.01);
+    OpenMMContext context(system, integrator, platform);
+    OpenMMContextImpl* impl = *reinterpret_cast<OpenMMContextImpl**>(&context);
     Stream stream = platform.createStream("", size, type, *impl);
     const int length = size*WIDTH;
     T array[size*WIDTH+1];

platforms/reference/src/gbsa/CpuObc.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -29,6 +29,7 @@
 #include "../SimTKUtilities/SimTKOpenMMLog.h"
 #include "../SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "CpuObc.h"
+#include "../SimTKReference/ReferenceForce.h"
 #include <math.h>
 
 /**---------------------------------------------------------------------------------------
@@ -265,12 +266,14 @@ int CpuObc::computeBornRadii( RealOpenMM** atomCoordinates, RealOpenMM* bornRadi
 
          if( atomJ != atomI ){
 
-            RealOpenMM deltaX          = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
-            RealOpenMM deltaY          = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
-            RealOpenMM deltaZ          = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
- 
-            RealOpenMM r2              = deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ;
-            RealOpenMM r               = SQRT( r2 );
+            RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+            if (_obcParameters->getPeriodic())
+                ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
+            else
+                ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
+            RealOpenMM r               = deltaR[ReferenceForce::RIndex];
+            if (_obcParameters->getUseCutoff() && r > _obcParameters->getCutoffDistance())
+                continue;
             RealOpenMM offsetRadiusJ   = atomicRadii[atomJ] - dielectricOffset; 
             RealOpenMM scaledRadiusJ   = offsetRadiusJ*scaledRadiusFactor[atomJ];
             RealOpenMM rScaledRadiusJ  = r + scaledRadiusJ;
@@ -417,15 +420,17 @@ int CpuObc::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoo
       RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
       for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
 
-         // 3 FLOP
-
-         RealOpenMM deltaX             = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
-         RealOpenMM deltaY             = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
-         RealOpenMM deltaZ             = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
- 
-         // 5 FLOP
-
-         RealOpenMM r2                 = deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ;
+         RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+         if (_obcParameters->getPeriodic())
+             ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
+         else
+             ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
+         if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
+             continue;
+         RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
+         RealOpenMM deltaX             = deltaR[ReferenceForce::XIndex];
+         RealOpenMM deltaY             = deltaR[ReferenceForce::YIndex];
+         RealOpenMM deltaZ             = deltaR[ReferenceForce::ZIndex];
 
          // 3 FLOP
 
@@ -957,15 +962,17 @@ FILE* logFile = NULL;
       RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
       for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
 
-         // 3 FLOP
-
-         RealOpenMM deltaX             = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
-         RealOpenMM deltaY             = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
-         RealOpenMM deltaZ             = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
- 
-         // 5 FLOP
-
-         RealOpenMM r2                 = deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ;
+         RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+         if (_obcParameters->getPeriodic())
+             ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
+         else
+             ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
+         if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
+             continue;
+         RealOpenMM r2                 = deltaR[ReferenceForce::R2Index];
+         RealOpenMM deltaX             = deltaR[ReferenceForce::XIndex];
+         RealOpenMM deltaY             = deltaR[ReferenceForce::YIndex];
+         RealOpenMM deltaZ             = deltaR[ReferenceForce::ZIndex];
 
          // 3 FLOP
 

platforms/reference/src/gbsa/ObcParameters.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -101,7 +101,7 @@ const std::string ObcParameters::ParameterFileName = std::string( "params.agb" )
 
    --------------------------------------------------------------------------------------- */
 
-ObcParameters::ObcParameters( int numberOfAtoms, ObcParameters::ObcType obcType ) : ImplicitSolventParameters( numberOfAtoms ){
+ObcParameters::ObcParameters( int numberOfAtoms, ObcParameters::ObcType obcType ) : ImplicitSolventParameters( numberOfAtoms ), cutoff(false), periodic(false) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -656,3 +656,84 @@ int ObcParameters::isNotReady( void ) const {
    return errors;
 }
 
+/**---------------------------------------------------------------------------------------
+
+     Set the force to use a cutoff.
+
+     @param distance            the cutoff distance
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+int ObcParameters::setUseCutoff( RealOpenMM distance ) {
+
+    cutoff = true;
+    cutoffDistance = distance;
+    return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get whether to use a cutoff.
+
+     --------------------------------------------------------------------------------------- */
+
+bool ObcParameters::getUseCutoff() {
+    return cutoff;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get the cutoff distance.
+
+     --------------------------------------------------------------------------------------- */
+
+RealOpenMM ObcParameters::getCutoffDistance() {
+    return cutoffDistance;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     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
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+int ObcParameters::setPeriodic( RealOpenMM* boxSize ) {
+
+    assert(cutoff);
+    assert(boxSize[0] >= 2.0*cutoffDistance);
+    assert(boxSize[1] >= 2.0*cutoffDistance);
+    assert(boxSize[2] >= 2.0*cutoffDistance);
+    periodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+    return SimTKOpenMMCommon::DefaultReturn;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get whether to use periodic boundary conditions.
+
+     --------------------------------------------------------------------------------------- */
+
+bool ObcParameters::getPeriodic() {
+    return periodic;
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Get the periodic box dimension
+
+     --------------------------------------------------------------------------------------- */
+
+const RealOpenMM* ObcParameters::getPeriodicBox() {
+    return periodicBoxSize;
+}

platforms/reference/src/gbsa/ObcParameters.h

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -55,6 +55,14 @@ class ObcParameters : public ImplicitSolventParameters {
       int _ownScaledRadiusFactors;
       RealOpenMM* _scaledRadiusFactors;
 
+      // cutoff and periodic boundary conditions
+      
+      bool cutoff;
+      bool periodic;
+      RealOpenMM periodicBoxSize[3];
+      RealOpenMM cutoffDistance;
+
+
       /**---------------------------------------------------------------------------------------
       
          Set solvent dielectric (Simbios) 
@@ -286,6 +294,64 @@ class ObcParameters : public ImplicitSolventParameters {
        
       int isNotReady( void ) const;
 
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use a cutoff.
+
+         @param distance            the cutoff distance
+
+         @return SimTKOpenMMCommon::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int setUseCutoff( RealOpenMM distance );
+
+      /**---------------------------------------------------------------------------------------
+
+         Get whether to use a cutoff.
+
+         --------------------------------------------------------------------------------------- */
+
+      bool getUseCutoff();
+
+      /**---------------------------------------------------------------------------------------
+
+         Get the cutoff distance.
+
+         --------------------------------------------------------------------------------------- */
+
+      RealOpenMM getCutoffDistance();
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use periodic boundary conditions.  This requires that a cutoff has
+         already been set, and the smallest side of the periodic box is at least twice the cutoff
+         distance.
+
+         @param boxSize             the X, Y, and Z widths of the periodic box
+
+         @return SimTKOpenMMCommon::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int setPeriodic( RealOpenMM* boxSize );
+
+      /**---------------------------------------------------------------------------------------
+
+         Get whether to use periodic boundary conditions.
+
+         --------------------------------------------------------------------------------------- */
+
+      bool getPeriodic();
+
+      /**---------------------------------------------------------------------------------------
+
+         Get the periodic box dimension
+
+         --------------------------------------------------------------------------------------- */
+
+      const RealOpenMM* getPeriodicBox();
+
 };
    
 /**---------------------------------------------------------------------------------------

platforms/reference/tests/TestReferenceGBSAOBCForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -39,6 +39,7 @@
 #include "GBSAOBCForce.h"
 #include "System.h"
 #include "LangevinIntegrator.h"
+#include "NonbondedForce.h"
 #include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
 #include "../src/sfmt/SFMT.h"
 #include <iostream>
@@ -70,6 +71,56 @@ void testSingleParticle() {
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
 }
 
+void testCutoffAndPeriodic() {
+    ReferencePlatform platform;
+    System system(2, 0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* gbsa = new GBSAOBCForce(2);
+    NonbondedForce* nonbonded = new NonbondedForce(2, 0);
+    gbsa->setParticleParameters(0, -1, 0.15, 1);
+    nonbonded->setParticleParameters(0, -1, 1, 0);
+    gbsa->setParticleParameters(1, 1, 0.15, 1);
+    nonbonded->setParticleParameters(1, 1, 1, 0);
+    const double cutoffDistance = 3.0;
+    const double boxSize = 10.0;
+    nonbonded->setCutoffDistance(cutoffDistance);
+    nonbonded->setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+
+    // Calculate the forces for both cutoff and periodic with two different atom positions.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    OpenMMContext context(system, integrator, platform);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state2 = context.getState(State::Forces);
+    positions[1][0]+= boxSize;
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state3 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state4 = context.getState(State::Forces);
+
+    // All forces should be identical, exception state3 which should be zero.
+
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state2.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state2.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state4.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state4.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[0], Vec3(0, 0, 0), 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
+}
+
 void testForce() {
     ReferencePlatform platform;
     const int numParticles = 10;
@@ -116,6 +167,7 @@ void testForce() {
 int main() {
     try {
         testSingleParticle();
+        testCutoffAndPeriodic();
         testForce();
     }
     catch(const exception& e) {