Created Cuda implementation of AndersenThermostat. Tweaked test cases for BrownianIntegrator.

platforms/cuda/src/CudaKernelFactory.cpp

@@ -55,8 +55,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaIntegrateLangevinStepKernel(name, platform, data);
     if (name == IntegrateBrownianStepKernel::Name())
         return new CudaIntegrateBrownianStepKernel(name, platform, data);
-//    if (name == ApplyAndersenThermostatKernel::Name())
-//        return new CudaApplyAndersenThermostatKernel(name, platform);
+    if (name == ApplyAndersenThermostatKernel::Name())
+        return new CudaApplyAndersenThermostatKernel(name, platform, data);
     if (name == CalcKineticEnergyKernel::Name())
         return new CudaCalcKineticEnergyKernel(name, platform);
     if (name == RemoveCMMotionKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -471,15 +471,31 @@ void CudaIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const
     }
     kBrownianUpdatePart2(gpu);
 }
-//
-//CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
-//}
-//
-//void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
-//}
-//
-//void CudaApplyAndersenThermostatKernel::execute(OpenMMContextImpl& context) {
-//}
+
+CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
+}
+
+void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
+    prevStepSize = -1.0;
+}
+
+void CudaApplyAndersenThermostatKernel::execute(OpenMMContextImpl& context) {
+    _gpuContext* gpu = data.gpu;
+    double temperature = context.getParameter(AndersenThermostat::Temperature);
+    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency);
+    double stepSize = context.getIntegrator().getStepSize();
+    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
+        // Initialize the GPU parameters.
+        
+        gpuSetAndersenThermostatParameters(gpu, temperature, frequency*stepSize);
+        gpuSetConstants(gpu);
+        kGenerateRandoms(gpu);
+        prevTemp = temperature;
+        prevFrequency = frequency;
+        prevStepSize = stepSize;
+    }
+    kCalculateAndersenThermostat(gpu);
+}
 
 void CudaCalcKineticEnergyKernel::initialize(const System& system) {
     int numParticles = system.getNumParticles();

platforms/cuda/src/CudaKernels.h

@@ -331,32 +331,32 @@ private:
     CudaPlatform::PlatformData& data;
     double prevTemp, prevFriction, prevStepSize;
 };
-//
-///**
-// * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
-// */
-//class CudaApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
-//public:
-//    CudaApplyAndersenThermostatKernel(std::string name, const Platform& platform) : ApplyAndersenThermostatKernel(name, platform), thermostat(0) {
-//    }
-//    ~CudaApplyAndersenThermostatKernel();
-//    /**
-//     * Initialize the kernel.
-//     * 
-//     * @param system     the System this kernel will be applied to
-//     * @param thermostat the AndersenThermostat this kernel will be used for
-//     */
-//    void initialize(const System& system, const AndersenThermostat& thermostat);
-//    /**
-//     * Execute the kernel.
-//     * 
-//     * @param context    the context in which to execute this kernel
-//     */
-//    void execute(OpenMMContextImpl& context);
-//private:
-//    CudaAndersenThermostat* thermostat;
-//    RealOpenMM* masses;
-//};
+
+/**
+ * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
+ */
+class CudaApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
+public:
+    CudaApplyAndersenThermostatKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data) : ApplyAndersenThermostatKernel(name, platform), data(data) {
+    }
+    ~CudaApplyAndersenThermostatKernel();
+    /**
+     * Initialize the kernel.
+     * 
+     * @param system     the System this kernel will be applied to
+     * @param thermostat the AndersenThermostat this kernel will be used for
+     */
+    void initialize(const System& system, const AndersenThermostat& thermostat);
+    /**
+     * Execute the kernel.
+     * 
+     * @param context    the context in which to execute this kernel
+     */
+    void execute(OpenMMContextImpl& context);
+private:
+    CudaPlatform::PlatformData& data;
+    double prevTemp, prevFrequency, prevStepSize;
+};
 
 /**
  * This kernel is invoked to calculate the kinetic energy of the system.

platforms/cuda/src/CudaPlatform.cpp

@@ -49,7 +49,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
-//    registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
+    registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
     registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
 }

platforms/cuda/tests/TestCudaBrownianIntegrator.cpp

@@ -90,7 +90,7 @@ void testSingleBond() {
 void testTemperature() {
     const int numParticles = 8;
     const int numBonds = numParticles-1;
-    const double temp = 100.0;
+    const double temp = 10.0;
     CudaPlatform platform;
     System system(numParticles, 0);
     BrownianIntegrator integrator(temp, 2.0, 0.01);
@@ -98,7 +98,7 @@ void testTemperature() {
     for (int i = 0; i < numParticles; ++i)
         system.setParticleMass(i, 2.0);
     for (int i = 0; i < numBonds; ++i)
-        forceField->setBondParameters(i, i, i+1, 1.0, i+1);
+        forceField->setBondParameters(i, i, i+1, 1.0, 5.0);
     system.addForce(forceField);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(numParticles);
@@ -108,12 +108,12 @@ void testTemperature() {
     
     // Let it equilibrate.
     
-    integrator.step(1000);
+    integrator.step(10000);
     
     // Now run it for a while and see if the temperature is correct.
     
     double pe = 0.0;
-    const int steps = 10000;
+    const int steps = 50000;
     for (int i = 0; i < steps; ++i) {
         State state = context.getState(State::Energy);
         pe += state.getPotentialEnergy();
@@ -121,8 +121,7 @@ void testTemperature() {
     }
     pe /= steps;
     double expected = 0.5*numBonds*BOLTZ*temp;
-    expected *= sqrt(2.0);
-    ASSERT_EQUAL_TOL(expected, pe, 4*expected/std::sqrt(steps));
+    ASSERT_EQUAL_TOL(expected, pe, 20*expected/std::sqrt(steps));
 }
 
 void testConstraints() {

platforms/reference/tests/TestReferenceBrownianIntegrator.cpp

@@ -87,7 +87,7 @@ void testSingleBond() {
 void testTemperature() {
     const int numParticles = 8;
     const int numBonds = numParticles-1;
-    const double temp = 100.0;
+    const double temp = 10.0;
     ReferencePlatform platform;
     System system(numParticles, 0);
     BrownianIntegrator integrator(temp, 2.0, 0.01);
@@ -96,7 +96,7 @@ void testTemperature() {
         system.setParticleMass(i, 2.0);
     }
     for (int i = 0; i < numBonds; ++i)
-        forceField->setBondParameters(i, i, i+1, 1.0, i+1);
+        forceField->setBondParameters(i, i, i+1, 1.0, 5.0);
     system.addForce(forceField);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(numParticles);
@@ -106,7 +106,7 @@ void testTemperature() {
     
     // Let it equilibrate.
     
-    integrator.step(1000);
+    integrator.step(10000);
     
     // Now run it for a while and see if the temperature is correct.
     
@@ -119,8 +119,7 @@ void testTemperature() {
     }
     pe /= steps;
     double expected = 0.5*numBonds*BOLTZ*temp;
-    expected *= sqrt(2.0);
-    ASSERT_EQUAL_TOL(expected, pe, 4*expected/std::sqrt(steps));
+    ASSERT_EQUAL_TOL(expected, pe, 20*expected/std::sqrt(steps));
 }
 
 void testConstraints() {