Bug fixes and simplifications to integration

platforms/cuda/tests/TestCudaLangevinIntegrator.cpp

@@ -126,7 +126,7 @@ void testTemperature() {
     }
     ke /= 1000;
     double expected = 0.5*numParticles*3*BOLTZ*temp;
-    ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
+    ASSERT_EQUAL_TOL(expected, ke, 3/std::sqrt(1000.0));
 }
 
 void testConstraints() {

platforms/opencl/src/OpenCLIntegrationUtilities.cpp

@@ -67,12 +67,11 @@ struct OpenCLIntegrationUtilities::ShakeCluster {
 };
 
 OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, const System& system) : context(context),
-        oldPos(NULL), posDelta(NULL), settleAtoms(NULL), settleParams(NULL), shakeAtoms(NULL), shakeParams(NULL),
+        posDelta(NULL), settleAtoms(NULL), settleParams(NULL), shakeAtoms(NULL), shakeParams(NULL),
         random(NULL), randomSeed(NULL), randomPos(NULL) {
     // Create workspace arrays.
 
     posDelta = new OpenCLArray<mm_float4>(context, context.getPaddedNumAtoms(), "posDelta");
-    oldPos = new OpenCLArray<mm_float4>(context, context.getPaddedNumAtoms(), "oldPos");
 
     // Create kernels for enforcing constraints.
 
@@ -249,8 +248,6 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
 OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
     if (posDelta != NULL)
         delete posDelta;
-    if (oldPos != NULL)
-        delete oldPos;
     if (settleAtoms != NULL)
         delete settleAtoms;
     if (settleParams != NULL)
@@ -265,13 +262,13 @@ OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
         delete randomSeed;
 }
 
-void OpenCLIntegrationUtilities::applyConstraints(double tol, OpenCLArray<mm_float4>& oldPositions, OpenCLArray<mm_float4>& positionDeltas, OpenCLArray<mm_float4>& newDeltas) {
+void OpenCLIntegrationUtilities::applyConstraints(double tol) {
     if (settleAtoms != NULL) {
         settleKernel.setArg<cl_int>(0, settleAtoms->getSize());
         settleKernel.setArg<cl_float>(1, tol);
-        settleKernel.setArg<cl::Buffer>(2, oldPositions.getDeviceBuffer());
-        settleKernel.setArg<cl::Buffer>(3, positionDeltas.getDeviceBuffer());
-        settleKernel.setArg<cl::Buffer>(4, newDeltas.getDeviceBuffer());
+        settleKernel.setArg<cl::Buffer>(2, context.getPosq().getDeviceBuffer());
+        settleKernel.setArg<cl::Buffer>(3, posDelta->getDeviceBuffer());
+        settleKernel.setArg<cl::Buffer>(4, posDelta->getDeviceBuffer());
         settleKernel.setArg<cl::Buffer>(5, context.getVelm().getDeviceBuffer());
         settleKernel.setArg<cl::Buffer>(6, settleAtoms->getDeviceBuffer());
         settleKernel.setArg<cl::Buffer>(7, settleParams->getDeviceBuffer());
@@ -280,9 +277,9 @@ void OpenCLIntegrationUtilities::applyConstraints(double tol, OpenCLArray<mm_flo
     if (shakeAtoms != NULL) {
         shakeKernel.setArg<cl_int>(0, shakeAtoms->getSize());
         shakeKernel.setArg<cl_float>(1, tol);
-        shakeKernel.setArg<cl::Buffer>(2, oldPositions.getDeviceBuffer());
-        shakeKernel.setArg<cl::Buffer>(3, positionDeltas.getDeviceBuffer());
-        shakeKernel.setArg<cl::Buffer>(4, newDeltas.getDeviceBuffer());
+        shakeKernel.setArg<cl::Buffer>(2, context.getPosq().getDeviceBuffer());
+        shakeKernel.setArg<cl::Buffer>(3, posDelta->getDeviceBuffer());
+        shakeKernel.setArg<cl::Buffer>(4, posDelta->getDeviceBuffer());
         shakeKernel.setArg<cl::Buffer>(5, shakeAtoms->getDeviceBuffer());
         shakeKernel.setArg<cl::Buffer>(6, shakeParams->getDeviceBuffer());
         context.executeKernel(shakeKernel, shakeAtoms->getSize());
@@ -299,7 +296,7 @@ void OpenCLIntegrationUtilities::initRandomNumberGenerator(unsigned int randomNu
     // Create the random number arrays.
 
     lastSeed = randomNumberSeed;
-    random = new OpenCLArray<mm_float4>(context, 32*context.getNumAtoms(), "random");
+    random = new OpenCLArray<mm_float4>(context, 32*context.getPaddedNumAtoms(), "random");
     randomSeed = new OpenCLArray<mm_int4>(context, context.getNumThreadBlocks()*OpenCLContext::ThreadBlockSize, "randomSeed");
     randomPos = random->getSize();
 
@@ -328,7 +325,5 @@ int OpenCLIntegrationUtilities::prepareRandomNumbers(int numValues) {
     randomKernel.setArg<cl::Buffer>(2, randomSeed->getDeviceBuffer());
     context.executeKernel(randomKernel, random->getSize());
     randomPos = numValues;
-    vector<mm_float4> r(random->getSize());
-    random->download(r);
     return 0;
 }

platforms/opencl/src/OpenCLIntegrationUtilities.h

@@ -47,12 +47,6 @@ public:
     OpenCLArray<mm_float4>& getPosDelta() {
         return *posDelta;
     }
-    /**
-     * Get the array which contains positions from before the current time step.
-     */
-    OpenCLArray<mm_float4>& getOldPos() {
-        return *oldPos;
-    }
     /**
      * Get the array which contains random values.
      */
@@ -63,11 +57,8 @@ public:
      * Apply constraints to the atom positions.
      *
      * @param tol             the constraint tolerance
-     * @param oldPositions    the atom positions from before the current time step
-     * @param positionDeltas  the offsets being added to atom positions
-     * @param newDeltas       the array to store constrained deltas into
      */
-    void applyConstraints(double tol, OpenCLArray<mm_float4>& oldPositions, OpenCLArray<mm_float4>& positionDeltas, OpenCLArray<mm_float4>& newDeltas);
+    void applyConstraints(double tol);
     /**
      * Initialize the random number generator.
      */
@@ -85,7 +76,6 @@ private:
     cl::Kernel shakeKernel;
     cl::Kernel randomKernel;
     OpenCLArray<mm_float4>* posDelta;
-    OpenCLArray<mm_float4>* oldPos;
     OpenCLArray<mm_int4>* settleAtoms;
     OpenCLArray<mm_float2>* settleParams;
     OpenCLArray<mm_int4>* shakeAtoms;

platforms/opencl/src/OpenCLKernels.cpp

@@ -820,12 +820,11 @@ void OpenCLIntegrateVerletStepKernel::execute(ContextImpl& context, const Verlet
     kernel1.setArg<cl::Buffer>(3, cl.getVelm().getDeviceBuffer());
     kernel1.setArg<cl::Buffer>(4, cl.getForce().getDeviceBuffer());
     kernel1.setArg<cl::Buffer>(5, integration.getPosDelta().getDeviceBuffer());
-    kernel1.setArg<cl::Buffer>(6, integration.getOldPos().getDeviceBuffer());
     cl.executeKernel(kernel1, numAtoms);
 
     // Apply constraints.
 
-    integration.applyConstraints(integrator.getConstraintTolerance(), integration.getOldPos(), integration.getPosDelta(), integration.getPosDelta());
+    integration.applyConstraints(integrator.getConstraintTolerance());
 
     // Call the second integration kernel.
 
@@ -869,6 +868,7 @@ void OpenCLIntegrateLangevinStepKernel::initialize(const System& system, const L
 void OpenCLIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
     OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilties();
     int numAtoms = cl.getNumAtoms();
+    int numThreads = cl.getNumThreadBlocks()*cl.ThreadBlockSize;
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
     double stepSize = integrator.getStepSize();
@@ -940,21 +940,20 @@ void OpenCLIntegrateLangevinStepKernel::execute(ContextImpl& context, const Lang
     // Call the first integration kernel.
 
     kernel1.setArg<cl_int>(0, numAtoms);
-    kernel1.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
-    kernel1.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
-    kernel1.setArg<cl::Buffer>(3, cl.getForce().getDeviceBuffer());
-    kernel1.setArg<cl::Buffer>(4, integration.getPosDelta().getDeviceBuffer());
-    kernel1.setArg<cl::Buffer>(5, params->getDeviceBuffer());
-    kernel1.setArg(6, params->getSize()*sizeof(cl_float), NULL);
-    kernel1.setArg<cl::Buffer>(7, xVector->getDeviceBuffer());
-    kernel1.setArg<cl::Buffer>(8, vVector->getDeviceBuffer());
-    kernel1.setArg<cl::Buffer>(9,integration.getRandom().getDeviceBuffer());
-    kernel1.setArg<cl_uint>(10, integration.prepareRandomNumbers(2*cl.getPaddedNumAtoms()));
+    kernel1.setArg<cl::Buffer>(1, cl.getVelm().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(2, cl.getForce().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(3, integration.getPosDelta().getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(4, params->getDeviceBuffer());
+    kernel1.setArg(5, params->getSize()*sizeof(cl_float), NULL);
+    kernel1.setArg<cl::Buffer>(6, xVector->getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(7, vVector->getDeviceBuffer());
+    kernel1.setArg<cl::Buffer>(8,integration.getRandom().getDeviceBuffer());
+    kernel1.setArg<cl_uint>(9, integration.prepareRandomNumbers(2*numThreads));
     cl.executeKernel(kernel1, numAtoms);
 
     // Apply constraints.
 
-    integration.applyConstraints(integrator.getConstraintTolerance(), integration.getOldPos(), integration.getPosDelta(), integration.getPosDelta());
+    integration.applyConstraints(integrator.getConstraintTolerance());
 
     // Call the second integration kernel.
 
@@ -966,12 +965,12 @@ void OpenCLIntegrateLangevinStepKernel::execute(ContextImpl& context, const Lang
     kernel2.setArg<cl::Buffer>(5, xVector->getDeviceBuffer());
     kernel2.setArg<cl::Buffer>(6, vVector->getDeviceBuffer());
     kernel2.setArg<cl::Buffer>(7,integration.getRandom().getDeviceBuffer());
-    kernel2.setArg<cl_uint>(8, integration.prepareRandomNumbers(2*cl.getPaddedNumAtoms()));
+    kernel2.setArg<cl_uint>(8, integration.prepareRandomNumbers(2*numThreads));
     cl.executeKernel(kernel2, numAtoms);
 
     // Reapply constraints.
 
-    integration.applyConstraints(integrator.getConstraintTolerance(), integration.getOldPos(), integration.getPosDelta(), cl.getPosq());
+    integration.applyConstraints(integrator.getConstraintTolerance());
 
     // Call the third integration kernel.
 

platforms/opencl/src/kernels/langevin.cl

@@ -4,9 +4,9 @@ enum {EM, EM_V, DOverTauC, TauOneMinusEM_V, TauDOverEMMinusOne, V, X, Yv, Yx, Fi
  * Perform the first step of Langevin integration.
  */
 
-__kernel void integrateLangevinPart1(int numAtoms, __global float4* posq, __global float4* velm, __global float4* force,
-        __global float4* posDelta, __global float* paramBuffer, __local float* params, __global float4* xVector,
-        __global float4* vVector, __global float4* random, unsigned int randomIndex) {
+__kernel void integrateLangevinPart1(int numAtoms, __global float4* velm, __global float4* force, __global float4* posDelta,
+        __global float* paramBuffer, __local float* params, __global float4* xVector, __global float4* vVector,
+        __global float4* random, unsigned int randomIndex) {
 
     // Load the parameters into local memory for faster access.
 
@@ -47,7 +47,7 @@ __kernel void integrateLangevinPart2(int numAtoms, __global float4* velm, __glob
     while (index < numAtoms) {
         float4 delta = posDelta[index];
         float4 velocity = velm[index];
-        float sqrtInvMass = 0.0f;//sqrt(velocity.w);
+        float sqrtInvMass = sqrt(velocity.w);
         velocity.xyz = delta.xyz*params[OneOverFix1];
         float4 xmh = (float4) (vVector[index].xyz*params[TauDOverEMMinusOne] + sqrtInvMass*params[Yx]*random[randomIndex].xyz, 0.0f);
         randomIndex += get_global_size(0);

platforms/opencl/src/kernels/verlet.cl

@@ -2,12 +2,11 @@
  * Perform the first step of verlet integration.
  */
 
-__kernel void integrateVerletPart1(int numAtoms, float dt, __global float4* posq, __global float4* velm, __global float4* force, __global float4* posDelta, __global float4* oldPos) {
+__kernel void integrateVerletPart1(int numAtoms, float dt, __global float4* posq, __global float4* velm, __global float4* force, __global float4* posDelta) {
     int index = get_global_id(0);
     while (index < numAtoms) {
         float4 pos = posq[index];
         float4 velocity = velm[index];
-        oldPos[index] = pos;
         velocity.xyz += force[index].xyz*dt*velocity.w;
         pos.xyz = velocity.xyz*dt;
         posDelta[index] = pos;

platforms/opencl/tests/TestOpenCLLangevinIntegrator.cpp

@@ -126,7 +126,7 @@ void testTemperature() {
     }
     ke /= 1000;
     double expected = 0.5*numParticles*3*BOLTZ*temp;
-    ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
+    ASSERT_EQUAL_TOL(expected, ke, 3/std::sqrt(1000.0));
 }
 
 void testConstraints() {

platforms/opencl/tests/TestOpenCLSettle.cpp

@@ -39,7 +39,7 @@
 #include "OpenCLPlatform.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
-#include "openmm/VerletIntegrator.h"
+#include "openmm/LangevinIntegrator.h"
 #include "../src/sfmt/SFMT.h"
 #include <iostream>
 #include <vector>
@@ -54,7 +54,7 @@ void testConstraints() {
     const double temp = 100.0;
     OpenCLPlatform platform;
     System system;
-    VerletIntegrator integrator(0.001);
+    LangevinIntegrator integrator(temp, 2.0, 0.001);
     integrator.setConstraintTolerance(1e-5);
     NonbondedForce* forceField = new NonbondedForce();
     for (int i = 0; i < numMolecules; ++i) {

platforms/reference/tests/TestReferenceLangevinIntegrator.cpp

@@ -126,7 +126,7 @@ void testTemperature() {
     }
     ke /= 1000;
     double expected = 0.5*numParticles*3*BOLTZ*temp;
-    ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
+    ASSERT_EQUAL_TOL(expected, ke, 3/std::sqrt(1000.0));
 }
 
 void testConstraints() {