Merge pull request #162 from peastman/master

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,7 +205,8 @@ 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};
-    cu.executeKernel(pileKernel, pileArgs, numParticles*numCopies, workgroupSize);
+    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.
 
-    randomIndex = integration.prepareRandomNumbers(numParticles*numCopies);
-    cu.executeKernel(pileKernel, pileArgs, numParticles*numCopies, workgroupSize);
+    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,7 +223,8 @@ 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);
     }
-    cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize);
+    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.
 
-    pileKernel.setArg<cl_uint>(2, integration.prepareRandomNumbers(numParticles*numCopies));
-    cl.executeKernel(pileKernel, numParticles*numCopies, workgroupSize);
+    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,36 +135,38 @@ 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);
-    for (int particle = 0; particle < numParticles; particle++) {
-        if (system.getParticleMass(particle) == 0.0)
-            continue;
-        const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
-        for (int component = 0; component < 3; component++) {
-            for (int k = 0; k < numCopies; k++)
-                v[k] = t_complex(scale*velocities[k][particle][component], 0.0);
-            fftpack_exec_1d(fft, FFTPACK_FORWARD, &v[0], &v[0]);
-            
-            // Apply a local Langevin thermostat to the centroid mode.
+    if (integrator.getApplyThermostat()) {
+        for (int particle = 0; particle < numParticles; particle++) {
+            if (system.getParticleMass(particle) == 0.0)
+                continue;
+            const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
+            for (int component = 0; component < 3; component++) {
+                for (int k = 0; k < numCopies; k++)
+                    v[k] = t_complex(scale*velocities[k][particle][component], 0.0);
+                fftpack_exec_1d(fft, FFTPACK_FORWARD, &v[0], &v[0]);
 
-            v[0].re = v[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                // Apply a local Langevin thermostat to the centroid mode.
 
-            // Use critical damping white noise for the remaining modes.
-            
-            for (int k = 1; k <= numCopies/2; k++) {
-                const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
-                const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
-                const RealOpenMM c1 = exp(-2.0*wk*halfdt);
-                const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
-                const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle));
-                RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-                RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
-                v[k] = v[k]*c1 + t_complex(rand1, rand2);
-                if (k < numCopies-k)
-                    v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
+                v[0].re = v[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+
+                // Use critical damping white noise for the remaining modes.
+
+                for (int k = 1; k <= numCopies/2; k++) {
+                    const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
+                    const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
+                    const RealOpenMM c1 = exp(-2.0*wk*halfdt);
+                    const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
+                    const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle));
+                    RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                    RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
+                    v[k] = v[k]*c1 + t_complex(rand1, rand2);
+                    if (k < numCopies-k)
+                        v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
+                }
+                fftpack_exec_1d(fft, FFTPACK_BACKWARD, &v[0], &v[0]);
+                for (int k = 0; k < numCopies; k++)
+                    velocities[k][particle][component] = scale*v[k].re;
             }
-            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &v[0], &v[0]);
-            for (int k = 0; k < numCopies; k++)
-                velocities[k][particle][component] = scale*v[k].re;
         }
     }
 
@@ -220,36 +222,38 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
 
     // Apply the PILE-L thermostat again.
     
-    for (int particle = 0; particle < numParticles; particle++) {
-        if (system.getParticleMass(particle) == 0.0)
-            continue;
-        const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
-        for (int component = 0; component < 3; component++) {
-            for (int k = 0; k < numCopies; k++)
-                v[k] = t_complex(scale*velocities[k][particle][component], 0.0);
-            fftpack_exec_1d(fft, FFTPACK_FORWARD, &v[0], &v[0]);
-            
-            // Apply a local Langevin thermostat to the centroid mode.
+    if (integrator.getApplyThermostat()) {
+        for (int particle = 0; particle < numParticles; particle++) {
+            if (system.getParticleMass(particle) == 0.0)
+                continue;
+            const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
+            for (int component = 0; component < 3; component++) {
+                for (int k = 0; k < numCopies; k++)
+                    v[k] = t_complex(scale*velocities[k][particle][component], 0.0);
+                fftpack_exec_1d(fft, FFTPACK_FORWARD, &v[0], &v[0]);
 
-            v[0].re = v[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                // Apply a local Langevin thermostat to the centroid mode.
 
-            // Use critical damping white noise for the remaining modes.
-            
-            for (int k = 1; k <= numCopies/2; k++) {
-                const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
-                const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
-                const RealOpenMM c1 = exp(-2.0*wk*halfdt);
-                const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
-                const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle));
-                RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-                RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
-                v[k] = v[k]*c1 + t_complex(rand1, rand2);
-                if (k < numCopies-k)
-                    v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
+                v[0].re = v[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+
+                // Use critical damping white noise for the remaining modes.
+
+                for (int k = 1; k <= numCopies/2; k++) {
+                    const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
+                    const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
+                    const RealOpenMM c1 = exp(-2.0*wk*halfdt);
+                    const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
+                    const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle));
+                    RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                    RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
+                    v[k] = v[k]*c1 + t_complex(rand1, rand2);
+                    if (k < numCopies-k)
+                        v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
+                }
+                fftpack_exec_1d(fft, FFTPACK_BACKWARD, &v[0], &v[0]);
+                for (int k = 0; k < numCopies; k++)
+                    velocities[k][particle][component] = scale*v[k].re;
             }
-            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &v[0], &v[0]);
-            for (int k = 0; k < numCopies; k++)
-                velocities[k][particle][component] = scale*v[k].re;
         }
     }
     

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;