BrownianIntegrator scales the friction coefficient by the mass

openmmapi/include/openmm/BrownianIntegrator.h

@@ -48,7 +48,7 @@ public:
      * Create a BrownianIntegrator.
      * 
      * @param temperature    the temperature of the heat bath (in Kelvin)
-     * @param frictionCoeff  the friction coefficient which couples the system to the heat bath
+     * @param frictionCoeff  the friction coefficient which couples the system to the heat bath, measured in 1/ps
      * @param stepSize       the step size with which to integrator the system (in picoseconds)
      */
     BrownianIntegrator(double temperature, double frictionCoeff, double stepSize);

platforms/cuda/src/kernels/kBrownianUpdate.cu

@@ -62,11 +62,14 @@ __global__ void kBrownianUpdatePart1_kernel()
         float4 random4a         = cSim.pRandom4a[rpos + pos];
         float4 apos             = cSim.pPosq[pos];
         float4 force            = cSim.pForce4[pos];
+        float invMass           = cSim.pVelm4[pos].w;
+        float forceScale        = cSim.tauDeltaT*invMass;
+        float noiseScale        = cSim.noiseAmplitude*sqrt(invMass);
 
         cSim.pOldPosq[pos]      = apos;
-        apos.x                  = force.x*cSim.tauDeltaT + cSim.noiseAmplitude*random4a.x;
-        apos.y                  = force.y*cSim.tauDeltaT + cSim.noiseAmplitude*random4a.y;
-        apos.z                  = force.z*cSim.tauDeltaT + cSim.noiseAmplitude*random4a.z;
+        apos.x                  = force.x*forceScale + noiseScale*random4a.x;
+        apos.y                  = force.y*forceScale + noiseScale*random4a.y;
+        apos.z                  = force.z*forceScale + noiseScale*random4a.z;
         cSim.pPosqP[pos]        = apos;
         pos                    += blockDim.x * gridDim.x;
     }

platforms/cuda/tests/TestCudaBrownianIntegrator.cpp

@@ -71,7 +71,7 @@ void testSingleBond() {
     
     // This is simply an overdamped harmonic oscillator, so compare it to the analytical solution.
     
-    double rate = 2*1.0/0.1;
+    double rate = 2*1.0/(0.1*2.0);
     for (int i = 0; i < 1000; ++i) {
         State state = context.getState(State::Positions | State::Velocities);
         double time = state.getTime();

platforms/opencl/src/OpenCLKernels.cpp

@@ -2680,7 +2680,8 @@ void OpenCLIntegrateBrownianStepKernel::execute(ContextImpl& context, const Brow
         hasInitializedKernels = true;
         kernel1.setArg<cl::Buffer>(2, cl.getForce().getDeviceBuffer());
         kernel1.setArg<cl::Buffer>(3, integration.getPosDelta().getDeviceBuffer());
-        kernel1.setArg<cl::Buffer>(4,integration.getRandom().getDeviceBuffer());
+        kernel1.setArg<cl::Buffer>(4, cl.getVelm().getDeviceBuffer());
+        kernel1.setArg<cl::Buffer>(5, integration.getRandom().getDeviceBuffer());
         kernel2.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
         kernel2.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
         kernel2.setArg<cl::Buffer>(3, integration.getPosDelta().getDeviceBuffer());
@@ -2700,7 +2701,7 @@ void OpenCLIntegrateBrownianStepKernel::execute(ContextImpl& context, const Brow
 
     // Call the first integration kernel.
 
-    kernel1.setArg<cl_uint>(5, integration.prepareRandomNumbers(cl.getPaddedNumAtoms()));
+    kernel1.setArg<cl_uint>(6, integration.prepareRandomNumbers(cl.getPaddedNumAtoms()));
     cl.executeKernel(kernel1, numAtoms);
 
     // Apply constraints.

platforms/opencl/src/kernels/brownian.cl

@@ -3,10 +3,11 @@
  */
 
 __kernel void integrateBrownianPart1(float tauDeltaT, float noiseAmplitude, __global float4* force,
-        __global float4* posDelta, __global float4* random, unsigned int randomIndex) {
+        __global float4* posDelta, __global float4* velm, __global float4* random, unsigned int randomIndex) {
     randomIndex += get_global_id(0);
     for (int index = get_global_id(0); index < NUM_ATOMS; index += get_global_size(0)) {
-        posDelta[index] = (float4) (tauDeltaT*force[index].xyz + noiseAmplitude*random[randomIndex].xyz, 0.0f);
+        float invMass = velm[index].w;
+        posDelta[index] = (float4) (tauDeltaT*invMass*force[index].xyz + noiseAmplitude*sqrt(invMass)*random[randomIndex].xyz, 0.0f);
         randomIndex += get_global_size(0);
     }
 }

platforms/opencl/tests/TestOpenCLBrownianIntegrator.cpp

@@ -71,7 +71,7 @@ void testSingleBond() {
 
     // This is simply an overdamped harmonic oscillator, so compare it to the analytical solution.
 
-    double rate = 2*1.0/0.1;
+    double rate = 2*1.0/(0.1*2.0);
     for (int i = 0; i < 1000; ++i) {
         State state = context.getState(State::Positions | State::Velocities);
         double time = state.getTime();

platforms/reference/src/SimTKReference/ReferenceBrownianDynamics.cpp

@@ -202,7 +202,7 @@ int ReferenceBrownianDynamics::update( int numberOfAtoms, RealOpenMM** atomCoord
    const RealOpenMM forceScale = getDeltaT()/getFriction();
    for (int i = 0; i < numberOfAtoms; ++i) {
        for (int j = 0; j < 3; ++j) {
-           xPrime[i][j] = atomCoordinates[i][j] + forceScale*forces[i][j] + noiseAmplitude*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+           xPrime[i][j] = atomCoordinates[i][j] + forceScale*inverseMasses[i]*forces[i][j] + noiseAmplitude*SQRT(inverseMasses[i])*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
        }
    }
    ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();

platforms/reference/tests/TestReferenceBrownianIntegrator.cpp

@@ -68,7 +68,7 @@ void testSingleBond() {
     
     // This is simply an overdamped harmonic oscillator, so compare it to the analytical solution.
     
-    double rate = 2*1.0/0.1;
+    double rate = 2*1.0/(0.1*2.0);
     for (int i = 0; i < 1000; ++i) {
         State state = context.getState(State::Positions | State::Velocities);
         double time = state.getTime();