Cleaning up the RPMD code

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

@@ -168,7 +168,7 @@ protected:
 private:
     double temperature, friction;
     int numCopies, randomNumberSeed;
-    bool forcesAreValid;
+    bool forcesAreValid, hasSetPosition, hasSetVelocity;
     ContextImpl* context;
     Context* owner;
     Kernel kernel;

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -42,7 +42,7 @@ using std::string;
 using std::vector;
 
 RPMDIntegrator::RPMDIntegrator(int numCopies, double temperature, double frictionCoeff, double stepSize) :
-        owner(NULL), numCopies(numCopies), forcesAreValid(false) {
+        owner(NULL), numCopies(numCopies), forcesAreValid(false), hasSetPosition(false), hasSetVelocity(false) {
     setTemperature(temperature);
     setFriction(frictionCoeff);
     setStepSize(stepSize);
@@ -71,10 +71,12 @@ vector<string> RPMDIntegrator::getKernelNames() {
 
 void RPMDIntegrator::setPositions(int copy, const vector<Vec3>& positions) {
     dynamic_cast<IntegrateRPMDStepKernel&>(kernel.getImpl()).setPositions(copy, positions);
+    hasSetPosition = true;
 }
 
 void RPMDIntegrator::setVelocities(int copy, const vector<Vec3>& velocities) {
     dynamic_cast<IntegrateRPMDStepKernel&>(kernel.getImpl()).setVelocities(copy, velocities);
+    hasSetVelocity = true;
 }
 
 State RPMDIntegrator::getState(int copy, int types) {
@@ -83,6 +85,20 @@ State RPMDIntegrator::getState(int copy, int types) {
 }
 
 void RPMDIntegrator::step(int steps) {
+    if (!hasSetPosition) {
+        // Initialize the positions from the context.
+        
+        State s = context->getOwner().getState(State::Positions);
+        for (int i = 0; i < numCopies; i++)
+            setPositions(i, s.getPositions());
+    }
+    if (!hasSetVelocity) {
+        // Initialize the velocities from the context.
+        
+        State s = context->getOwner().getState(State::Velocities);
+        for (int i = 0; i < numCopies; i++)
+            setVelocities(i, s.getVelocities());
+    }
     for (int i = 0; i < steps; ++i) {
         dynamic_cast<IntegrateRPMDStepKernel&>(kernel.getImpl()).execute(*context, *this, forcesAreValid);
         forcesAreValid = true;

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

@@ -79,22 +79,6 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     vector<RealVec>& pos = extractPositions(context);
     vector<RealVec>& vel = extractVelocities(context);
     vector<RealVec>& f = extractForces(context);
-    if (!hasSetPosition) {
-        // Initialize the positions from the context.
-        
-        for (int i = 0; i < numCopies; i++)
-            for (int j = 0; j < numParticles; j++)
-                positions[i][j] = pos[j];
-        hasSetPosition = true;
-    }
-    if (!hasSetVelocity) {
-        // Initialize the velocities from the context.
-        
-        for (int i = 0; i < numCopies; i++)
-            for (int j = 0; j < numParticles; j++)
-                velocities[i][j] = vel[j];
-        hasSetVelocity = true;
-    }
     
     // Loop over copies and compute the force on each one.
     
@@ -109,24 +93,24 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
 
     // Apply the PILE-L thermostat.
     
-    vector<t_complex> p(numCopies);
+    vector<t_complex> v(numCopies);
     vector<t_complex> q(numCopies);
     const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
-    const RealOpenMM scale = 1.0/sqrt(numCopies);
+    const RealOpenMM scale = 1.0/sqrt((RealOpenMM) numCopies);
     const RealOpenMM nkT = numCopies*BOLTZ*integrator.getTemperature();
     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++) {
-        const RealOpenMM c3_0 = c2_0*sqrt(system.getParticleMass(particle)*nkT);
+        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++)
-                p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
-            fftpack_exec_1d(fft, FFTPACK_FORWARD, &p[0], &p[0]);
+                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.
 
-            p[0].re = p[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+            v[0].re = v[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
 
             // Use critical damping white noise for the remaining modes.
             
@@ -135,16 +119,16 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
                 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(system.getParticleMass(particle)*nkT);
+                const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle));
                 RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
                 RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
-                p[k] = p[k]*c1 + t_complex(rand1, rand2);
+                v[k] = v[k]*c1 + t_complex(rand1, rand2);
                 if (k < numCopies-k)
-                    p[numCopies-k] = p[numCopies-k]*c1 + t_complex(rand1, -rand2);                
+                    v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
             }
-            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &p[0], &p[0]);
+            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &v[0], &v[0]);
             for (int k = 0; k < numCopies; k++)
-                velocities[k][particle][component] = scale*p[k].re/system.getParticleMass(particle);
+                velocities[k][particle][component] = scale*v[k].re;
         }
     }
 
@@ -160,26 +144,26 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
         for (int component = 0; component < 3; component++) {
             for (int k = 0; k < numCopies; k++) {
                 q[k] = t_complex(scale*positions[k][particle][component], 0.0);
-                p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
+                v[k] = t_complex(scale*velocities[k][particle][component], 0.0);
             }
             fftpack_exec_1d(fft, FFTPACK_FORWARD, &q[0], &q[0]);
-            fftpack_exec_1d(fft, FFTPACK_FORWARD, &p[0], &p[0]);
-            q[0] += p[0]*(dt/system.getParticleMass(particle));
+            fftpack_exec_1d(fft, FFTPACK_FORWARD, &v[0], &v[0]);
+            q[0] += v[0]*dt;
             for (int k = 1; k < numCopies; k++) {
                 const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
                 const RealOpenMM wt = wk*dt;
                 const RealOpenMM coswt = cos(wt);
                 const RealOpenMM sinwt = sin(wt);
                 const RealOpenMM wm = wk*system.getParticleMass(particle);
-                const t_complex pprime = p[k]*coswt - q[k]*(wm*sinwt); // Advance momentum from t to t+dt
-                q[k] = p[k]*(sinwt/wm) + q[k]*coswt; // Advance position from t to t+dt
-                p[k] = pprime;
+                const t_complex vprime = v[k]*coswt - q[k]*(wk*sinwt); // Advance velocity from t to t+dt
+                q[k] = v[k]*(sinwt/wk) + q[k]*coswt; // Advance position from t to t+dt
+                v[k] = vprime;
             }
             fftpack_exec_1d(fft, FFTPACK_BACKWARD, &q[0], &q[0]);
-            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &p[0], &p[0]);
+            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &v[0], &v[0]);
             for (int k = 0; k < numCopies; k++) {
                 positions[k][particle][component] = scale*q[k].re;
-                velocities[k][particle][component] = scale*p[k].re/system.getParticleMass(particle);
+                velocities[k][particle][component] = scale*v[k].re;
             }
         }
     }
@@ -206,15 +190,15 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     // Apply the PILE-L thermostat again.
     
     for (int particle = 0; particle < numParticles; particle++) {
-        const RealOpenMM c3_0 = c2_0*sqrt(system.getParticleMass(particle)*nkT);
+        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++)
-                p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
-            fftpack_exec_1d(fft, FFTPACK_FORWARD, &p[0], &p[0]);
+                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.
 
-            p[0].re = p[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+            v[0].re = v[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
 
             // Use critical damping white noise for the remaining modes.
             
@@ -223,16 +207,16 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
                 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(system.getParticleMass(particle)*nkT);
+                const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle));
                 RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
                 RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
-                p[k] = p[k]*c1 + t_complex(rand1, rand2);
+                v[k] = v[k]*c1 + t_complex(rand1, rand2);
                 if (k < numCopies-k)
-                    p[numCopies-k] = p[numCopies-k]*c1 + t_complex(rand1, -rand2);                
+                    v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
             }
-            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &p[0], &p[0]);
+            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &v[0], &v[0]);
             for (int k = 0; k < numCopies; k++)
-                velocities[k][particle][component] = scale*p[k].re/system.getParticleMass(particle);
+                velocities[k][particle][component] = scale*v[k].re;
         }
     }
     
@@ -247,22 +231,6 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
 //    vector<RealVec>& pos = extractPositions(context);
 //    vector<RealVec>& vel = extractVelocities(context);
 //    vector<RealVec>& f = extractForces(context);
-//    if (!hasSetPosition) {
-//        // Initialize the positions from the context.
-//        
-//        for (int i = 0; i < numCopies; i++)
-//            for (int j = 0; j < numParticles; j++)
-//                positions[i][j] = pos[j];
-//        hasSetPosition = true;
-//    }
-//    if (!hasSetVelocity) {
-//        // Initialize the velocities from the context.
-//        
-//        for (int i = 0; i < numCopies; i++)
-//            for (int j = 0; j < numParticles; j++)
-//                velocities[i][j] = vel[j];
-//        hasSetVelocity = true;
-//    }
 //    
 //    // Loop over copies and compute the force on each one.
 //    
@@ -356,14 +324,12 @@ void ReferenceIntegrateRPMDStepKernel::setPositions(int copy, const vector<Vec3>
     int numParticles = positions[copy].size();
     for (int i = 0; i < numParticles; i++)
         positions[copy][i] = pos[i];
-    hasSetPosition = true;
 }
 
 void ReferenceIntegrateRPMDStepKernel::setVelocities(int copy, const vector<Vec3>& vel) {
     int numParticles = velocities[copy].size();
     for (int i = 0; i < numParticles; i++)
         velocities[copy][i] = vel[i];
-    hasSetVelocity = true;
 }
 
 void ReferenceIntegrateRPMDStepKernel::copyToContext(int copy, ContextImpl& context) const {

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

@@ -46,7 +46,7 @@ namespace OpenMM {
 class ReferenceIntegrateRPMDStepKernel : public IntegrateRPMDStepKernel {
 public:
     ReferenceIntegrateRPMDStepKernel(std::string name, const Platform& platform) :
-            IntegrateRPMDStepKernel(name, platform), fft(NULL), hasSetPosition(false), hasSetVelocity(false) {
+            IntegrateRPMDStepKernel(name, platform), fft(NULL) {
     }
     ~ReferenceIntegrateRPMDStepKernel();
     /**
@@ -82,7 +82,6 @@ private:
     std::vector<std::vector<RealVec> > velocities;
     std::vector<std::vector<RealVec> > forces;
     fftpack* fft;
-    bool hasSetPosition, hasSetVelocity;
 };
 
 } // namespace OpenMM