RPMD works correctly with a barostat

plugins/rpmd/platforms/cuda/src/CudaRpmdKernels.cpp

@@ -252,7 +252,18 @@ void CudaIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
                 &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
         cu.executeKernel(copyToContextKernel, copyToContextArgs, cu.getNumAtoms());
         context.computeVirtualSites();
+        Vec3 initialBox[3];
+        context.getPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
         context.updateContextState();
+        Vec3 finalBox[3];
+        context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
+        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2]) {
+            // A barostat was applied during updateContextState().  Adjust the particle positions in all the
+            // other copies to match this one.
+            
+            void* args[] = {&positions->getDevicePointer(), &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
+            cu.executeKernel(translateKernel, args, cu.getNumAtoms());
+        }
         context.calcForcesAndEnergy(true, false, groupsNotContracted);
         void* copyFromContextArgs[] = {&cu.getForce().getDevicePointer(), &forces->getDevicePointer(), &cu.getVelm().getDevicePointer(),
                 &velocities->getDevicePointer(), &cu.getPosq().getDevicePointer(), &positions->getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};

plugins/rpmd/platforms/cuda/src/kernels/rpmd.cu

@@ -217,7 +217,7 @@ extern "C" __global__ void copyDataFromContext(long long* srcForce, long long* d
 }
 
 /**
- * Update atom positions so all copies are offset by the same number of periodic box widths.
+ * Atom positions in one copy have been modified.  Apply the same offsets to all the other copies.
  */
 extern "C" __global__ void applyCellTranslations(mixed4* posq, real4* movedPos, int* order, int movedCopy) {
     for (int particle = blockIdx.x*blockDim.x+threadIdx.x; particle < NUM_ATOMS; particle += blockDim.x*gridDim.x) {

plugins/rpmd/platforms/cuda/tests/TestCudaRpmd.cpp

@@ -38,6 +38,7 @@
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
+#include "openmm/MonteCarloBarostat.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
@@ -496,6 +497,80 @@ void testWithoutThermostat() {
     }
 }
 
+void testWithBarostat() {
+    const int gridSize = 3;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int numParticles = numMolecules*2;
+    const int numCopies = 10;
+    const double spacing = 2.0;
+    const double cutoff = 3.0;
+    const double boxSize = spacing*(gridSize+1);
+    const double temperature = 300.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setForceGroup(1);
+    nonbonded->setReciprocalSpaceForceGroup(2);
+    system.addForce(nonbonded);
+    system.addForce(new MonteCarloBarostat(0.5, temperature));
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(-0.2, 0.2, 0.2);
+        nonbonded->addParticle(0.2, 0.2, 0.2);
+        nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
+        bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
+    }
+    RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("CUDA");
+    Context context(system, integ, platform);
+    for (int copy = 0; copy < numCopies; copy++) {
+        for (int i = 0; i < gridSize; i++)
+            for (int j = 0; j < gridSize; j++)
+                for (int k = 0; k < gridSize; k++) {
+                    Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
+                    int index = k+gridSize*(j+gridSize*i);
+                    positions[2*index] = pos;
+                    positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
+                }
+        integ.setPositions(copy, positions);
+    }
+
+    // Check the temperature.
+    
+    const int numSteps = 500;
+    integ.step(100);
+    vector<double> ke(numCopies, 0.0);
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        vector<State> state(numCopies);
+        for (int j = 0; j < numCopies; j++)
+            state[j] = integ.getState(j, State::Velocities, true);
+        for (int j = 0; j < numParticles; j++) {
+            for (int k = 0; k < numCopies; k++) {
+                Vec3 v = state[k].getVelocities()[j];
+                ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
+            }
+        }
+    }
+    double meanKE = 0.0;
+    for (int i = 0; i < numCopies; i++)
+        meanKE += ke[i];
+    meanKE /= numSteps*numCopies;
+    double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
+    ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDCudaKernelFactories();
@@ -507,6 +582,7 @@ int main(int argc, char* argv[]) {
         testVirtualSites();
         testContractions();
         testWithoutThermostat();
+        testWithBarostat();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.cpp

@@ -268,7 +268,18 @@ void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
         copyToContextKernel.setArg<cl_int>(5, i);
         cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
         context.computeVirtualSites();
+        Vec3 initialBox[3];
+        context.getPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
         context.updateContextState();
+        Vec3 finalBox[3];
+        context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
+        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2]) {
+            // A barostat was applied during updateContextState().  Adjust the particle positions in all the
+            // other copies to match this one.
+            
+            translateKernel.setArg<cl_int>(3, i);
+            cl.executeKernel(translateKernel, cl.getNumAtoms());
+        }
         context.calcForcesAndEnergy(true, false, groupsNotContracted);
         copyFromContextKernel.setArg<cl_int>(7, i);
         cl.executeKernel(copyFromContextKernel, cl.getNumAtoms());

plugins/rpmd/platforms/opencl/src/kernels/rpmd.cl

@@ -204,7 +204,7 @@ __kernel void copyDataFromContext(__global real4* srcForce, __global real4* dstF
 }
 
 /**
- * Update atom positions so all copies are offset by the same number of periodic box widths.
+ * Atom positions in one copy have been modified.  Apply the same offsets to all the other copies.
  */
 __kernel void applyCellTranslations(__global mixed4* posq, __global real4* movedPos, __global int* order, int movedCopy) {
     for (int particle = get_global_id(0); particle < NUM_ATOMS; particle += get_global_size(0)) {

plugins/rpmd/platforms/opencl/tests/TestOpenCLRpmd.cpp

@@ -38,6 +38,7 @@
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
+#include "openmm/MonteCarloBarostat.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
@@ -497,6 +498,80 @@ void testWithoutThermostat() {
     }
 }
 
+void testWithBarostat() {
+    const int gridSize = 3;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int numParticles = numMolecules*2;
+    const int numCopies = 10;
+    const double spacing = 2.0;
+    const double cutoff = 3.0;
+    const double boxSize = spacing*(gridSize+1);
+    const double temperature = 300.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setForceGroup(1);
+    nonbonded->setReciprocalSpaceForceGroup(2);
+    system.addForce(nonbonded);
+    system.addForce(new MonteCarloBarostat(0.5, temperature));
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(-0.2, 0.2, 0.2);
+        nonbonded->addParticle(0.2, 0.2, 0.2);
+        nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
+        bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
+    }
+    RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("OpenCL");
+    Context context(system, integ, platform);
+    for (int copy = 0; copy < numCopies; copy++) {
+        for (int i = 0; i < gridSize; i++)
+            for (int j = 0; j < gridSize; j++)
+                for (int k = 0; k < gridSize; k++) {
+                    Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
+                    int index = k+gridSize*(j+gridSize*i);
+                    positions[2*index] = pos;
+                    positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
+                }
+        integ.setPositions(copy, positions);
+    }
+
+    // Check the temperature.
+    
+    const int numSteps = 500;
+    integ.step(100);
+    vector<double> ke(numCopies, 0.0);
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        vector<State> state(numCopies);
+        for (int j = 0; j < numCopies; j++)
+            state[j] = integ.getState(j, State::Velocities, true);
+        for (int j = 0; j < numParticles; j++) {
+            for (int k = 0; k < numCopies; k++) {
+                Vec3 v = state[k].getVelocities()[j];
+                ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
+            }
+        }
+    }
+    double meanKE = 0.0;
+    for (int i = 0; i < numCopies; i++)
+        meanKE += ke[i];
+    meanKE /= numSteps*numCopies;
+    double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
+    ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDOpenCLKernelFactories();
@@ -508,6 +583,7 @@ int main(int argc, char* argv[]) {
         testVirtualSites();
         testContractions();
         testWithoutThermostat();
+        testWithBarostat();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp

@@ -275,7 +275,22 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
         pos = positions[i];
         vel = velocities[i];
         context.computeVirtualSites();
+        Vec3 initialBox[3];
+        context.getPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
         context.updateContextState();
+        Vec3 finalBox[3];
+        context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
+        if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2]) {
+            // A barostat was applied during updateContextState().  Adjust the particle positions in all the
+            // other copies to match this one.
+            
+            for (int j = 0; j < numParticles; j++) {
+                Vec3 delta = pos[j]-positions[i][j];
+                for (int k = 0; k < totalCopies; k++)
+                    if (k != i)
+                        positions[k][j] += delta;
+            }
+        }
         positions[i] = pos;
         velocities[i] = vel;
         context.calcForcesAndEnergy(true, false, groupsNotContracted);

plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.cpp

@@ -37,6 +37,7 @@
 #include "openmm/CMMotionRemover.h"
 #include "openmm/Context.h"
 #include "openmm/HarmonicBondForce.h"
+#include "openmm/MonteCarloBarostat.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
@@ -378,6 +379,80 @@ void testWithoutThermostat() {
     }
 }
 
+void testWithBarostat() {
+    const int gridSize = 3;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int numParticles = numMolecules*2;
+    const int numCopies = 5;
+    const double spacing = 2.0;
+    const double cutoff = 3.0;
+    const double boxSize = spacing*(gridSize+1);
+    const double temperature = 300.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    nonbonded->setForceGroup(1);
+    nonbonded->setReciprocalSpaceForceGroup(2);
+    system.addForce(nonbonded);
+    system.addForce(new MonteCarloBarostat(0.5, temperature));
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(-0.2, 0.2, 0.2);
+        nonbonded->addParticle(0.2, 0.2, 0.2);
+        nonbonded->addException(2*i, 2*i+1, 0, 1, 0);
+        bonds->addBond(2*i, 2*i+1, 1.0, 10000.0);
+    }
+    RPMDIntegrator integ(numCopies, temperature, 50.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("Reference");
+    Context context(system, integ, platform);
+    for (int copy = 0; copy < numCopies; copy++) {
+        for (int i = 0; i < gridSize; i++)
+            for (int j = 0; j < gridSize; j++)
+                for (int k = 0; k < gridSize; k++) {
+                    Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
+                    int index = k+gridSize*(j+gridSize*i);
+                    positions[2*index] = pos;
+                    positions[2*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
+                }
+        integ.setPositions(copy, positions);
+    }
+
+    // Check the temperature.
+    
+    const int numSteps = 500;
+    integ.step(100);
+    vector<double> ke(numCopies, 0.0);
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        vector<State> state(numCopies);
+        for (int j = 0; j < numCopies; j++)
+            state[j] = integ.getState(j, State::Velocities, true);
+        for (int j = 0; j < numParticles; j++) {
+            for (int k = 0; k < numCopies; k++) {
+                Vec3 v = state[k].getVelocities()[j];
+                ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
+            }
+        }
+    }
+    double meanKE = 0.0;
+    for (int i = 0; i < numCopies; i++)
+        meanKE += ke[i];
+    meanKE /= numSteps*numCopies;
+    double expectedKE = 0.5*numCopies*numParticles*3*BOLTZ*temperature;
+    ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
+}
+
 int main() {
     try {
         testFreeParticles();
@@ -385,6 +460,7 @@ int main() {
         testVirtualSites();
         testContractions();
         testWithoutThermostat();
+        testWithBarostat();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;