More bug fixes to handling of random numbers by CustomIntegrator

platforms/cuda/src/CudaKernels.cpp

@@ -4842,7 +4842,8 @@ void CudaIntegrateCustomStepKernel::prepareForComputation(ContextImpl& context,
         for (int step = 1; step < numSteps; step++) {
             if (needsForces[step] || needsEnergy[step])
                 continue;
-            if (stepType[step-1] == CustomIntegrator::ComputeGlobal && stepType[step] == CustomIntegrator::ComputeGlobal)
+            if (stepType[step-1] == CustomIntegrator::ComputeGlobal && stepType[step] == CustomIntegrator::ComputeGlobal &&
+                    !usesVariable(expression[step], "uniform") && !usesVariable(expression[step], "gaussian"))
                 merged[step] = true;
             if (stepType[step-1] == CustomIntegrator::ComputePerDof && stepType[step] == CustomIntegrator::ComputePerDof)
                 merged[step] = true;
@@ -4867,9 +4868,9 @@ void CudaIntegrateCustomStepKernel::prepareForComputation(ContextImpl& context,
                     numUniform += numAtoms*usesVariable(expression[j], "uniform");
                     compute << "{\n";
                     if (numGaussian > 0)
-                        compute << "float4 gaussian = gaussianValues[gaussianIndex];\n";
+                        compute << "float4 gaussian = gaussianValues[gaussianIndex+index];\n";
                     if (numUniform > 0)
-                        compute << "float4 uniform = uniformValues[uniformIndex];\n";
+                        compute << "float4 uniform = uniformValues[uniformIndex+index];\n";
                     for (int i = 0; i < 3; i++)
                         compute << createPerDofComputation(stepType[j] == CustomIntegrator::ComputePerDof ? variable[j] : "", expression[j], i, integrator, forceName[j], energyName[j]);
                     if (variable[j] == "x") {
@@ -4889,9 +4890,9 @@ void CudaIntegrateCustomStepKernel::prepareForComputation(ContextImpl& context,
                         }
                     }
                     if (numGaussian > 0)
-                        compute << "gaussianIndex += blockDim.x*gridDim.x;\n";
+                        compute << "gaussianIndex += NUM_ATOMS;\n";
                     if (numUniform > 0)
-                        compute << "uniformIndex += blockDim.x*gridDim.x;\n";
+                        compute << "uniformIndex += NUM_ATOMS;\n";
                     compute << "}\n";
                 }
                 map<string, string> replacements;

platforms/cuda/src/kernels/customIntegratorPerDof.cu

@@ -34,12 +34,10 @@ inline __device__ mixed4 convertFromDouble4(double4 a) {
 extern "C" __global__ void computePerDof(real4* __restrict__ posq, real4* __restrict__ posqCorrection, mixed4* __restrict__ posDelta,
         mixed4* __restrict__ velm, const long long* __restrict__ force, const mixed2* __restrict__ dt, const mixed* __restrict__ globals,
         const mixed* __restrict__ params, mixed* __restrict__ sum, const float4* __restrict__ gaussianValues,
-        unsigned int gaussianIndex, const float4* __restrict__ uniformValues, const real* __restrict__ energy
+        unsigned int gaussianBaseIndex, const float4* __restrict__ uniformValues, const real* __restrict__ energy
         PARAMETER_ARGUMENTS) {
     mixed stepSize = dt[0].y;
     int index = blockIdx.x*blockDim.x+threadIdx.x;
-    gaussianIndex += index;
-    int uniformIndex = index;
     const double forceScale = 1.0/0xFFFFFFFF;
     while (index < NUM_ATOMS) {
 #ifdef LOAD_POS_AS_DELTA
@@ -51,6 +49,8 @@ extern "C" __global__ void computePerDof(real4* __restrict__ posq, real4* __rest
         double4 f = make_double4(forceScale*force[index], forceScale*force[index+PADDED_NUM_ATOMS], forceScale*force[index+PADDED_NUM_ATOMS*2], 0.0);
         double mass = 1.0/velocity.w;
         if (velocity.w != 0.0) {
+            int gaussianIndex = gaussianBaseIndex;
+            int uniformIndex = 0;
             COMPUTE_STEP
         }
         index += blockDim.x*gridDim.x;

platforms/cuda/tests/TestCudaCustomIntegrator.cpp

@@ -651,6 +651,72 @@ void testRespa() {
     }
 }
 
+/**
+ * Make sure random numbers are computed correctly when steps get merged.
+ */
+void testMergedRandoms() {
+    const int numParticles = 10;
+    const int numSteps = 10;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomIntegrator integrator(0.1);
+    integrator.addPerDofVariable("dofUniform1", 0);
+    integrator.addPerDofVariable("dofUniform2", 0);
+    integrator.addPerDofVariable("dofGaussian1", 0);
+    integrator.addPerDofVariable("dofGaussian2", 0);
+    integrator.addGlobalVariable("globalUniform1", 0);
+    integrator.addGlobalVariable("globalUniform2", 0);
+    integrator.addGlobalVariable("globalGaussian1", 0);
+    integrator.addGlobalVariable("globalGaussian2", 0);
+    integrator.addComputePerDof("dofUniform1", "uniform");
+    integrator.addComputePerDof("dofUniform2", "uniform");
+    integrator.addComputePerDof("dofGaussian1", "gaussian");
+    integrator.addComputePerDof("dofGaussian2", "gaussian");
+    integrator.addComputeGlobal("globalUniform1", "uniform");
+    integrator.addComputeGlobal("globalUniform2", "uniform");
+    integrator.addComputeGlobal("globalGaussian1", "gaussian");
+    integrator.addComputeGlobal("globalGaussian2", "gaussian");
+    Context context(system, integrator, platform);
+    
+    // See if the random numbers are computed correctly.
+    
+    vector<Vec3> values1, values2;
+    for (int i = 0; i < numSteps; i++) {
+        integrator.step(1);
+        integrator.getPerDofVariable(0, values1);
+        integrator.getPerDofVariable(1, values2);
+        for (int i = 0; i < numParticles; i++)
+            for (int j = 0; j < 3; j++) {
+                double v1 = values1[i][j];
+                double v2 = values2[i][j];
+                ASSERT(v1 >= 0 && v1 < 1);
+                ASSERT(v2 >= 0 && v2 < 1);
+                ASSERT(v1 != v2);
+            }
+        integrator.getPerDofVariable(2, values1);
+        integrator.getPerDofVariable(3, values2);
+        for (int i = 0; i < numParticles; i++)
+            for (int j = 0; j < 3; j++) {
+                double v1 = values1[i][j];
+                double v2 = values2[i][j];
+                ASSERT(v1 >= -10 && v1 < 10);
+                ASSERT(v2 >= -10 && v2 < 10);
+                ASSERT(v1 != v2);
+            }
+        double v1 = integrator.getGlobalVariable(0);
+        double v2 = integrator.getGlobalVariable(1);
+        ASSERT(v1 >= 0 && v1 < 1);
+        ASSERT(v2 >= 0 && v2 < 1);
+        ASSERT(v1 != v2);
+        v1 = integrator.getGlobalVariable(2);
+        v2 = integrator.getGlobalVariable(3);
+        ASSERT(v1 >= -10 && v1 < 10);
+        ASSERT(v2 >= -10 && v2 < 10);
+        ASSERT(v1 != v2);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -666,6 +732,7 @@ int main(int argc, char* argv[]) {
         testPerDofVariables();
         testForceGroups();
         testRespa();
+        testMergedRandoms();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/src/OpenCLKernels.cpp

@@ -5070,7 +5070,8 @@ void OpenCLIntegrateCustomStepKernel::prepareForComputation(ContextImpl& context
         for (int step = 1; step < numSteps; step++) {
             if (needsForces[step] || needsEnergy[step])
                 continue;
-            if (stepType[step-1] == CustomIntegrator::ComputeGlobal && stepType[step] == CustomIntegrator::ComputeGlobal)
+            if (stepType[step-1] == CustomIntegrator::ComputeGlobal && stepType[step] == CustomIntegrator::ComputeGlobal &&
+                    !usesVariable(expression[step], "uniform") && !usesVariable(expression[step], "gaussian"))
                 merged[step] = true;
             if (stepType[step-1] == CustomIntegrator::ComputePerDof && stepType[step] == CustomIntegrator::ComputePerDof)
                 merged[step] = true;
@@ -5095,9 +5096,9 @@ void OpenCLIntegrateCustomStepKernel::prepareForComputation(ContextImpl& context
                     numUniform += numAtoms*usesVariable(expression[j], "uniform");
                     compute << "{\n";
                     if (numGaussian > 0)
-                        compute << "float4 gaussian = gaussianValues[gaussianIndex];\n";
+                        compute << "float4 gaussian = gaussianValues[gaussianIndex+index];\n";
                     if (numUniform > 0)
-                        compute << "float4 uniform = uniformValues[uniformIndex];\n";
+                        compute << "float4 uniform = uniformValues[uniformIndex+index];\n";
                     for (int i = 0; i < 3; i++)
                         compute << createPerDofComputation(stepType[j] == CustomIntegrator::ComputePerDof ? variable[j] : "", expression[j], i, integrator, forceName[j], energyName[j]);
                     if (variable[j] == "x") {
@@ -5121,9 +5122,9 @@ void OpenCLIntegrateCustomStepKernel::prepareForComputation(ContextImpl& context
                         }
                     }
                     if (numGaussian > 0)
-                        compute << "gaussianIndex += get_global_size(0);\n";
+                        compute << "gaussianIndex += NUM_ATOMS;\n";
                     if (numUniform > 0)
-                        compute << "uniformIndex += get_global_size(0);\n";
+                        compute << "uniformIndex += NUM_ATOMS;\n";
                     compute << "}\n";
                 }
                 map<string, string> replacements;

platforms/opencl/src/kernels/customIntegratorPerDof.cl

@@ -26,12 +26,10 @@ void storePos(__global real4* restrict posq, __global real4* restrict posqCorrec
 __kernel void computePerDof(__global real4* restrict posq, __global real4* restrict posqCorrection, __global mixed4* restrict posDelta,
         __global mixed4* restrict velm, __global const real4* restrict force, __global const mixed2* restrict dt, __global const mixed* restrict globals,
         __global const mixed* restrict params, __global mixed* restrict sum, __global const float4* restrict gaussianValues,
-        unsigned int gaussianIndex, __global const float4* restrict uniformValues, __global const real* restrict energy
+        unsigned int gaussianBaseIndex, __global const float4* restrict uniformValues, __global const real* restrict energy
         PARAMETER_ARGUMENTS) {
     mixed stepSize = dt[0].y;
     int index = get_global_id(0);
-    gaussianIndex += index;
-    int uniformIndex = index;
     while (index < NUM_ATOMS) {
 #ifdef LOAD_POS_AS_DELTA
         mixed4 position = loadPos(posq, posqCorrection, index)+posDelta[index];
@@ -42,6 +40,8 @@ __kernel void computePerDof(__global real4* restrict posq, __global real4* restr
         real4 f = force[index];
         mixed mass = 1/velocity.w;
         if (velocity.w != 0.0) {
+            int gaussianIndex = gaussianBaseIndex;
+            int uniformIndex = 0;
             COMPUTE_STEP
         }
         index += get_global_size(0);

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -651,6 +651,72 @@ void testRespa() {
     }
 }
 
+/**
+ * Make sure random numbers are computed correctly when steps get merged.
+ */
+void testMergedRandoms() {
+    const int numParticles = 10;
+    const int numSteps = 10;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomIntegrator integrator(0.1);
+    integrator.addPerDofVariable("dofUniform1", 0);
+    integrator.addPerDofVariable("dofUniform2", 0);
+    integrator.addPerDofVariable("dofGaussian1", 0);
+    integrator.addPerDofVariable("dofGaussian2", 0);
+    integrator.addGlobalVariable("globalUniform1", 0);
+    integrator.addGlobalVariable("globalUniform2", 0);
+    integrator.addGlobalVariable("globalGaussian1", 0);
+    integrator.addGlobalVariable("globalGaussian2", 0);
+    integrator.addComputePerDof("dofUniform1", "uniform");
+    integrator.addComputePerDof("dofUniform2", "uniform");
+    integrator.addComputePerDof("dofGaussian1", "gaussian");
+    integrator.addComputePerDof("dofGaussian2", "gaussian");
+    integrator.addComputeGlobal("globalUniform1", "uniform");
+    integrator.addComputeGlobal("globalUniform2", "uniform");
+    integrator.addComputeGlobal("globalGaussian1", "gaussian");
+    integrator.addComputeGlobal("globalGaussian2", "gaussian");
+    Context context(system, integrator, platform);
+    
+    // See if the random numbers are computed correctly.
+    
+    vector<Vec3> values1, values2;
+    for (int i = 0; i < numSteps; i++) {
+        integrator.step(1);
+        integrator.getPerDofVariable(0, values1);
+        integrator.getPerDofVariable(1, values2);
+        for (int i = 0; i < numParticles; i++)
+            for (int j = 0; j < 3; j++) {
+                double v1 = values1[i][j];
+                double v2 = values2[i][j];
+                ASSERT(v1 >= 0 && v1 < 1);
+                ASSERT(v2 >= 0 && v2 < 1);
+                ASSERT(v1 != v2);
+            }
+        integrator.getPerDofVariable(2, values1);
+        integrator.getPerDofVariable(3, values2);
+        for (int i = 0; i < numParticles; i++)
+            for (int j = 0; j < 3; j++) {
+                double v1 = values1[i][j];
+                double v2 = values2[i][j];
+                ASSERT(v1 >= -10 && v1 < 10);
+                ASSERT(v2 >= -10 && v2 < 10);
+                ASSERT(v1 != v2);
+            }
+        double v1 = integrator.getGlobalVariable(0);
+        double v2 = integrator.getGlobalVariable(1);
+        ASSERT(v1 >= 0 && v1 < 1);
+        ASSERT(v2 >= 0 && v2 < 1);
+        ASSERT(v1 != v2);
+        v1 = integrator.getGlobalVariable(2);
+        v2 = integrator.getGlobalVariable(3);
+        ASSERT(v1 >= -10 && v1 < 10);
+        ASSERT(v2 >= -10 && v2 < 10);
+        ASSERT(v1 != v2);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -666,6 +732,7 @@ int main(int argc, char* argv[]) {
         testPerDofVariables();
         testForceGroups();
         testRespa();
+        testMergedRandoms();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;