Added option to omit PILE thermostat from RPMD

plugins/rpmd/openmmapi/include/openmm/RPMDIntegrator.h

@@ -47,6 +47,9 @@ namespace OpenMM {
  * springs to form a ring.  This allows certain quantum mechanical effects to be efficiently
  * simulated.
  * 
+ * By default this Integrator applies a PILE thermostat to the system to simulate constant
+ * temperature dynamics.  You can disable the thermostat by calling setApplyThermostat(false).
+ * 
  * Because this Integrator simulates many copies of the System at once, it must be used
  * differently from other Integrators.  Instead of setting positions and velocities by
  * calling methods of the Context, you should use the corresponding methods of the Integrator
@@ -127,6 +130,18 @@ public:
     void setFriction(double coeff) {
         friction = coeff;
     }
+    /**
+     * Get whether a thermostat is applied to the system.
+     */
+    bool getApplyThermostat() const {
+        return applyThermostat;
+    }
+    /**
+     * Set whether a thermostat is applied to the system.
+     */
+    void setApplyThermostat(bool apply) {
+        applyThermostat = apply;
+    }
     /**
      * Get the random number seed.  See setRandomNumberSeed() for details.
      */
@@ -213,6 +228,7 @@ protected:
 private:
     double temperature, friction;
     int numCopies, randomNumberSeed;
+    bool applyThermostat;
     std::map<int, int> contractions;
     bool forcesAreValid, hasSetPosition, hasSetVelocity, isFirstStep;
     Kernel kernel;

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -42,7 +42,7 @@ using namespace OpenMM;
 using namespace std;
 
 RPMDIntegrator::RPMDIntegrator(int numCopies, double temperature, double frictionCoeff, double stepSize, const map<int, int>& contractions) :
-        numCopies(numCopies), contractions(contractions), forcesAreValid(false), hasSetPosition(false), hasSetVelocity(false), isFirstStep(true) {
+        numCopies(numCopies), applyThermostat(true), contractions(contractions), forcesAreValid(false), hasSetPosition(false), hasSetVelocity(false), isFirstStep(true) {
     setTemperature(temperature);
     setFriction(frictionCoeff);
     setStepSize(stepSize);
@@ -51,7 +51,7 @@ RPMDIntegrator::RPMDIntegrator(int numCopies, double temperature, double frictio
 }
 
 RPMDIntegrator::RPMDIntegrator(int numCopies, double temperature, double frictionCoeff, double stepSize) :
-        numCopies(numCopies), forcesAreValid(false), hasSetPosition(false), hasSetVelocity(false), isFirstStep(true) {
+        numCopies(numCopies), applyThermostat(true), forcesAreValid(false), hasSetPosition(false), hasSetVelocity(false), isFirstStep(true) {
     setTemperature(temperature);
     setFriction(frictionCoeff);
     setStepSize(stepSize);

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

@@ -205,6 +205,7 @@ void CudaIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegr
     void* frictionPtr = (useDoublePrecision ? (void*) &friction : (void*) &frictionFloat);
     int randomIndex = integration.prepareRandomNumbers(numParticles*numCopies);
     void* pileArgs[] = {&velocities->getDevicePointer(), &integration.getRandom().getDevicePointer(), &randomIndex, dtPtr, kTPtr, frictionPtr};
+    if (integrator.getApplyThermostat())
         cu.executeKernel(pileKernel, pileArgs, numParticles*numCopies, workgroupSize);
 
     // Update positions and velocities.
@@ -223,8 +224,10 @@ void CudaIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegr
 
     // Apply the PILE-L thermostat again.
 
+    if (integrator.getApplyThermostat()) {
         randomIndex = integration.prepareRandomNumbers(numParticles*numCopies);
         cu.executeKernel(pileKernel, pileArgs, numParticles*numCopies, workgroupSize);
+    }
 
     // Update the time and step count.
 

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

@@ -431,6 +431,71 @@ void testContractions() {
     ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
 }
 
+void testWithoutThermostat() {
+    const int numParticles = 20;
+    const int numCopies = 10;
+    const double temperature = 300.0;
+    const double mass = 2.0;
+    
+    // Create a chain of particles.
+    
+    System system;
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(mass);
+        if (i > 0)
+            bonds->addBond(i-1, i, 1.0, 1000.0);
+    }
+    RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
+    integ.setApplyThermostat(false);
+    Platform& platform = Platform::getPlatformByName("CUDA");
+    Context context(system, integ, platform);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<vector<Vec3> > positions(numCopies);
+    for (int i = 0; i < numCopies; i++) {
+        positions[i].resize(numParticles);
+        for (int j = 0; j < numParticles; j++)
+            positions[i][j] = Vec3(0.95*j, 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
+        integ.setPositions(i, positions[i]);
+    }
+    
+    // Simulate it and see if the energy remains constant.
+    
+    double initialEnergy;
+    int numSteps = 100;
+    const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
+    const double wn = numCopies*BOLTZ*temperature/hbar;
+    const double springConstant = mass*wn*wn;
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        
+        // Sum the energies of all the copies.
+        
+        double energy = 0.0;
+        for (int j = 0; j < numCopies; j++) {
+            State state = integ.getState(j, State::Positions | State::Energy);
+            positions[j] = state.getPositions();
+            energy += state.getPotentialEnergy()+state.getKineticEnergy();
+        }
+        
+        // Add the energy from the springs connecting copies.
+        
+        for (int j = 0; j < numCopies; j++) {
+            int previous = (j == 0 ? numCopies-1 : j-1);
+            for (int k = 0; k < numParticles; k++) {
+                Vec3 delta = positions[j][k]-positions[previous][k];
+                energy += 0.5*springConstant*delta.dot(delta);
+            }
+        }
+        if (i == 0)
+            initialEnergy = energy;
+        else
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 1e-4);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDCudaKernelFactories();
@@ -441,6 +506,7 @@ int main(int argc, char* argv[]) {
         testCMMotionRemoval();
         testVirtualSites();
         testContractions();
+        testWithoutThermostat();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

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

@@ -223,6 +223,7 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte
         stepKernel.setArg<cl_float>(4, (cl_float) (integrator.getTemperature()*BOLTZ));
         velocitiesKernel.setArg<cl_float>(2, (cl_float) dt);
     }
+    if (integrator.getApplyThermostat())
         cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize);
 
     // Update positions and velocities.
@@ -238,8 +239,10 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte
 
     // Apply the PILE-L thermostat again.
 
+    if (integrator.getApplyThermostat()) {
         pileKernel.setArg<cl_uint>(2, integration.prepareRandomNumbers(numParticles*numCopies));
         cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize);
+    }
 
     // Update the time and step count.
 

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

@@ -432,6 +432,71 @@ void testContractions() {
     ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
 }
 
+void testWithoutThermostat() {
+    const int numParticles = 20;
+    const int numCopies = 10;
+    const double temperature = 300.0;
+    const double mass = 2.0;
+    
+    // Create a chain of particles.
+    
+    System system;
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(mass);
+        if (i > 0)
+            bonds->addBond(i-1, i, 1.0, 1000.0);
+    }
+    RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
+    integ.setApplyThermostat(false);
+    Platform& platform = Platform::getPlatformByName("OpenCL");
+    Context context(system, integ, platform);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<vector<Vec3> > positions(numCopies);
+    for (int i = 0; i < numCopies; i++) {
+        positions[i].resize(numParticles);
+        for (int j = 0; j < numParticles; j++)
+            positions[i][j] = Vec3(0.95*j, 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
+        integ.setPositions(i, positions[i]);
+    }
+    
+    // Simulate it and see if the energy remains constant.
+    
+    double initialEnergy;
+    int numSteps = 100;
+    const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
+    const double wn = numCopies*BOLTZ*temperature/hbar;
+    const double springConstant = mass*wn*wn;
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        
+        // Sum the energies of all the copies.
+        
+        double energy = 0.0;
+        for (int j = 0; j < numCopies; j++) {
+            State state = integ.getState(j, State::Positions | State::Energy);
+            positions[j] = state.getPositions();
+            energy += state.getPotentialEnergy()+state.getKineticEnergy();
+        }
+        
+        // Add the energy from the springs connecting copies.
+        
+        for (int j = 0; j < numCopies; j++) {
+            int previous = (j == 0 ? numCopies-1 : j-1);
+            for (int k = 0; k < numParticles; k++) {
+                Vec3 delta = positions[j][k]-positions[previous][k];
+                energy += 0.5*springConstant*delta.dot(delta);
+            }
+        }
+        if (i == 0)
+            initialEnergy = energy;
+        else
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 1e-4);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDOpenCLKernelFactories();
@@ -442,6 +507,7 @@ int main(int argc, char* argv[]) {
         testCMMotionRemoval();
         testVirtualSites();
         testContractions();
+        testWithoutThermostat();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

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

@@ -135,6 +135,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     const RealOpenMM twown = 2.0*nkT/hbar;
     const RealOpenMM c1_0 = exp(-halfdt*integrator.getFriction());
     const RealOpenMM c2_0 = sqrt(1.0-c1_0*c1_0);
+    if (integrator.getApplyThermostat()) {
         for (int particle = 0; particle < numParticles; particle++) {
             if (system.getParticleMass(particle) == 0.0)
                 continue;
@@ -167,6 +168,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
                     velocities[k][particle][component] = scale*v[k].re;
             }
         }
+    }
 
     // Update velocities.
     
@@ -220,6 +222,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
 
     // Apply the PILE-L thermostat again.
     
+    if (integrator.getApplyThermostat()) {
         for (int particle = 0; particle < numParticles; particle++) {
             if (system.getParticleMass(particle) == 0.0)
                 continue;
@@ -252,6 +255,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
                     velocities[k][particle][component] = scale*v[k].re;
             }
         }
+    }
     
     // Update the time.
     

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

@@ -313,12 +313,78 @@ void testContractions() {
     ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
 }
 
+void testWithoutThermostat() {
+    const int numParticles = 20;
+    const int numCopies = 10;
+    const double temperature = 300.0;
+    const double mass = 2.0;
+    
+    // Create a chain of particles.
+    
+    System system;
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(mass);
+        if (i > 0)
+            bonds->addBond(i-1, i, 1.0, 1000.0);
+    }
+    RPMDIntegrator integ(numCopies, temperature, 1.0, 0.001);
+    integ.setApplyThermostat(false);
+    Platform& platform = Platform::getPlatformByName("Reference");
+    Context context(system, integ, platform);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<vector<Vec3> > positions(numCopies);
+    for (int i = 0; i < numCopies; i++) {
+        positions[i].resize(numParticles);
+        for (int j = 0; j < numParticles; j++)
+            positions[i][j] = Vec3(0.95*j, 0.01*genrand_real2(sfmt), 0.01*genrand_real2(sfmt));
+        integ.setPositions(i, positions[i]);
+    }
+    
+    // Simulate it and see if the energy remains constant.
+    
+    double initialEnergy;
+    int numSteps = 100;
+    const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
+    const double wn = numCopies*BOLTZ*temperature/hbar;
+    const double springConstant = mass*wn*wn;
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        
+        // Sum the energies of all the copies.
+        
+        double energy = 0.0;
+        for (int j = 0; j < numCopies; j++) {
+            State state = integ.getState(j, State::Positions | State::Energy);
+            positions[j] = state.getPositions();
+            energy += state.getPotentialEnergy()+state.getKineticEnergy();
+        }
+        
+        // Add the energy from the springs connecting copies.
+        
+        for (int j = 0; j < numCopies; j++) {
+            int previous = (j == 0 ? numCopies-1 : j-1);
+            for (int k = 0; k < numParticles; k++) {
+                Vec3 delta = positions[j][k]-positions[previous][k];
+                energy += 0.5*springConstant*delta.dot(delta);
+            }
+        }
+        if (i == 0)
+            initialEnergy = energy;
+        else
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 1e-4);
+    }
+}
+
 int main() {
     try {
         testFreeParticles();
         testCMMotionRemoval();
         testVirtualSites();
         testContractions();
+        testWithoutThermostat();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;