Fixed bug in Langevin integration when used with constraints

platforms/cuda/src/CudaKernels.cpp

@@ -876,6 +876,7 @@ void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langev
     kApplySecondShake(gpu);
     kApplySecondSettle(gpu);
     kApplySecondCCMA(gpu);
+    kSetVelocitiesFromPositions(gpu);
     data.time += stepSize;
     data.stepCount++;
 }
@@ -995,6 +996,7 @@ void CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, cons
     kApplySecondShake(gpu);
     kApplySecondSettle(gpu);
     kApplySecondCCMA(gpu);
+    kSetVelocitiesFromPositions(gpu);
     gpu->psStepSize->Download();
     data.time += (*gpu->psStepSize)[0].y;
     if ((*gpu->psStepSize)[0].y == maxStepSize)

platforms/cuda/src/kernels/cudaKernels.h

@@ -62,6 +62,7 @@ extern void kApplySecondSettle(gpuContext gpu);
 extern void kLangevinUpdatePart1(gpuContext gpu);
 extern void kLangevinUpdatePart2(gpuContext gpu);
 extern void kSelectLangevinStepSize(gpuContext gpu, float maxTimeStep);
+extern void kSetVelocitiesFromPositions(gpuContext gpu);
 extern void kVerletUpdatePart1(gpuContext gpu);
 extern void kVerletUpdatePart2(gpuContext gpu);
 extern void kSelectVerletStepSize(gpuContext gpu, float maxTimeStep);

platforms/cuda/src/kernels/kLangevinUpdate.cu

@@ -169,3 +169,36 @@ void kSelectLangevinStepSize(gpuContext gpu, float maxTimeStep)
     kSelectLangevinStepSize_kernel<<<1, gpu->sim.update_threads_per_block, sizeof(float)*gpu->sim.update_threads_per_block>>>(maxTimeStep);
     LAUNCHERROR("kSelectLangevinStepSize");
 }
+
+__global__
+#if (__CUDA_ARCH__ >= 200)
+__launch_bounds__(GF1XX_UPDATE_THREADS_PER_BLOCK, 1)
+#elif (__CUDA_ARCH__ >= 130)
+__launch_bounds__(GT2XX_UPDATE_THREADS_PER_BLOCK, 1)
+#else
+__launch_bounds__(G8X_UPDATE_THREADS_PER_BLOCK, 1)
+#endif
+void kSetVelocitiesFromPositions_kernel()
+{
+    float2 stepSize = cSim.pStepSize[0];
+    float oneOverDt = 2.0f/(stepSize.x+stepSize.y);
+    unsigned int pos = threadIdx.x;
+    while (pos < cSim.atoms)
+    {
+        float4 posq = cSim.pPosq[pos];
+        float4 oldPosq = cSim.pOldPosq[pos];
+        float4 velm = cSim.pVelm4[pos];
+        velm.x = oneOverDt*(posq.x-oldPosq.x);
+        velm.y = oneOverDt*(posq.y-oldPosq.y);
+        velm.z = oneOverDt*(posq.z-oldPosq.z);
+        cSim.pVelm4[pos] = velm;
+        pos += blockDim.x * gridDim.x;
+    }
+}
+
+void kSetVelocitiesFromPositions(gpuContext gpu)
+{
+//    printf("kSetVelocitiesFromPositions\n");
+    kSetVelocitiesFromPositions_kernel<<<1, gpu->sim.update_threads_per_block, sizeof(float)*gpu->sim.update_threads_per_block>>>();
+    LAUNCHERROR("kSetVelocitiesFromPositions");
+}

platforms/cuda/src/kernels/kLangevinUpdate.h

@@ -126,7 +126,14 @@ void kLangevinUpdatePart2_kernel()
 {
     unsigned int pos            = threadIdx.x + blockIdx.x * blockDim.x;
     unsigned int rpos           = cSim.pRandomPosition[blockIdx.x];
-    float dt = cSim.pStepSize[0].y;
+    __shared__ float dt;
+    if (threadIdx.x == 0)
+    {
+        dt = cSim.pStepSize[0].y;
+        if (pos == 0)
+            cSim.pStepSize[0].x = dt;
+    }
+    __syncthreads();
 #ifdef REMOVE_CM
     extern __shared__ float3 sCM[];
     float3 CM                   = {0.0f, 0.0f, 0.0f};

platforms/opencl/src/OpenCLKernels.cpp

@@ -3023,6 +3023,8 @@ void OpenCLIntegrateLangevinStepKernel::execute(ContextImpl& context, const Lang
         kernel1.setArg<cl::Buffer>(5, integration.getRandom().getDeviceBuffer());
         kernel2.setArg<cl::Buffer>(0, cl.getPosq().getDeviceBuffer());
         kernel2.setArg<cl::Buffer>(1, integration.getPosDelta().getDeviceBuffer());
+        kernel2.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
+        kernel2.setArg<cl::Buffer>(3, integration.getStepSize().getDeviceBuffer());
     }
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
@@ -3224,6 +3226,8 @@ void OpenCLIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
         kernel1.setArg<cl::Buffer>(5, integration.getRandom().getDeviceBuffer());
         kernel2.setArg<cl::Buffer>(0, cl.getPosq().getDeviceBuffer());
         kernel2.setArg<cl::Buffer>(1, integration.getPosDelta().getDeviceBuffer());
+        kernel2.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
+        kernel2.setArg<cl::Buffer>(3, integration.getStepSize().getDeviceBuffer());
         selectSizeKernel.setArg<cl::Buffer>(4, integration.getStepSize().getDeviceBuffer());
         selectSizeKernel.setArg<cl::Buffer>(5, cl.getVelm().getDeviceBuffer());
         selectSizeKernel.setArg<cl::Buffer>(6, cl.getForce().getDeviceBuffer());

platforms/opencl/src/kernels/langevin.cl

@@ -27,13 +27,17 @@ __kernel void integrateLangevinPart1(__global float4* velm, __global float4* for
  * Perform the second step of Langevin integration.
  */
 
-__kernel void integrateLangevinPart2(__global float4* posq, __global float4* posDelta) {
+__kernel void integrateLangevinPart2(__global float4* posq, __global float4* posDelta, __global float4* velm, __global float2* dt) {
+    float invStepSize = 1.0f/dt[0].y;
     int index = get_global_id(0);
     while (index < NUM_ATOMS) {
         float4 pos = posq[index];
         float4 delta = posDelta[index];
+        float4 vel = velm[index];
         pos.xyz += delta.xyz;
+        vel.xyz = invStepSize*delta.xyz;
         posq[index] = pos;
+        velm[index] = vel;
         index += get_global_size(0);
     }
 }

platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp

@@ -293,11 +293,12 @@ int ReferenceStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoo
 
    // copy xPrime -> atomCoordinates
 
-   for( int ii = 0; ii < numberOfAtoms; ii++ ){
-      atomCoordinates[ii][0] = xPrime[ii][0];
-      atomCoordinates[ii][1] = xPrime[ii][1];
-      atomCoordinates[ii][2] = xPrime[ii][2];
-   }
+   RealOpenMM invStepSize = 1.0/getDeltaT();
+   for (int i = 0; i < numberOfAtoms; ++i)
+       for (int j = 0; j < 3; ++j) {
+           velocities[i][j] = invStepSize*(xPrime[i][j]-atomCoordinates[i][j]);
+           atomCoordinates[i][j] = xPrime[i][j];
+       }
 
    incrementTimeStep();