Bug fixes to reference implementation of DrudeLangevinIntegrator

plugins/drude/platforms/reference/src/ReferenceDrudeKernelFactory.cpp

@@ -52,9 +52,10 @@ extern "C" void initDrudeReferenceKernels() {
 }
 
 KernelImpl* ReferenceDrudeKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
+    ReferencePlatform::PlatformData& data = *static_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
     if (name == CalcDrudeForceKernel::Name())
         return new ReferenceCalcDrudeForceKernel(name, platform);
     if (name == IntegrateDrudeLangevinStepKernel::Name())
-        return new ReferenceIntegrateDrudeLangevinStepKernel(name, platform);
+        return new ReferenceIntegrateDrudeLangevinStepKernel(name, platform, data);
     throw OpenMMException((std::string("Tried to create kernel with illegal kernel name '")+name+"'").c_str());
 }

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp

@@ -35,6 +35,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "SimTKReference/ReferenceCCMAAlgorithm.h"
+#include "SimTKReference/ReferenceVirtualSites.h"
 #include <set>
 
 using namespace OpenMM;
@@ -256,22 +257,23 @@ void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system,
 }
 
 void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
+    vector<RealVec>& pos = extractPositions(context);
+    vector<RealVec>& vel = extractVelocities(context);
+    vector<RealVec>& force = extractForces(context);
+    
     // Update velocities of ordinary particles.
     
     const RealOpenMM vscale = exp(-integrator.getStepSize()*integrator.getFriction());
     const RealOpenMM fscale = (1-vscale)/integrator.getFriction();
     const RealOpenMM kT = BOLTZ*integrator.getTemperature();
     const RealOpenMM noisescale = sqrt(2*kT*integrator.getFriction())*sqrt(0.5*(1-vscale*vscale)/integrator.getFriction());
-    vector<RealVec>& pos = extractPositions(context);
-    vector<RealVec>& vel = extractVelocities(context);
-    vector<RealVec>& force = extractForces(context);
     for (int i = 0; i < (int) normalParticles.size(); i++) {
         int index = normalParticles[i];
         RealOpenMM invMass = particleInvMass[index];
         if (invMass != 0.0) {
             RealOpenMM sqrtInvMass = sqrt(invMass);
             for (int j = 0; j < 3; j++)
-                vel[i][j] = vscale*vel[i][j] + fscale*invMass*force[i][j] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+                vel[index][j] = vscale*vel[index][j] + fscale*invMass*force[index][j] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
         }
     }
     
@@ -284,18 +286,20 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
     for (int i = 0; i < (int) pairParticles.size(); i++) {
         int p1 = pairParticles[i].first;
         int p2 = pairParticles[i].second;
+        RealOpenMM mass1fract = pairInvTotalMass[i]/particleInvMass[p1];
+        RealOpenMM mass2fract = pairInvTotalMass[i]/particleInvMass[p2];
         RealOpenMM sqrtInvTotalMass = sqrt(pairInvTotalMass[i]);
         RealOpenMM sqrtInvReducedMass = sqrt(pairInvReducedMass[i]);
-        RealVec cmVel = (vel[p2]+vel[p1])*0.5;
+        RealVec cmVel = vel[p1]*mass1fract+vel[p2]*mass2fract;
         RealVec relVel = vel[p2]-vel[p1];
         RealVec cmForce = force[p1]+force[p2];
-        RealVec relForce = force[p2]*(pairInvTotalMass[i]/particleInvMass[p1]) - force[p1]*(pairInvTotalMass[i]/particleInvMass[p2]);
+        RealVec relForce = force[p2]*mass1fract - force[p1]*mass2fract;
         for (int j = 0; j < 3; j++) {
             cmVel[j] = vscale*cmVel[j] + fscale*pairInvTotalMass[i]*cmForce[j] + noisescale*sqrtInvTotalMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
             relVel[j] = vscaleDrude*relVel[j] + fscaleDrude*pairInvReducedMass[i]*relForce[j] + noisescaleDrude*sqrtInvReducedMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
         }
-        vel[p1] = cmVel-relVel*0.5;
-        vel[p2] = cmVel+relVel*0.5;
+        vel[p1] = cmVel-relVel*mass2fract;
+        vel[p2] = cmVel+relVel*mass1fract;
     }
 
     // Update the particle positions.
@@ -309,7 +313,8 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
     
     // Apply constraints.
     
-    constraints->apply(numParticles, pos, xPrime, particleInvMass);
+    if (constraints != NULL)
+        constraints->apply(numParticles, pos, xPrime, particleInvMass);
     
     // Record the constrained positions and velocities.
     
@@ -320,6 +325,9 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
             pos[i] = xPrime[i];
         }
     }
+    ReferenceVirtualSites::computePositions(context.getSystem(), pos);
+    data.time += integrator.getStepSize();
+    data.stepCount++;
 }
 
 double ReferenceIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
@@ -335,7 +343,7 @@ double ReferenceIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextIm
     double timeShift = 0.5*integrator.getStepSize();
     for (int i = 0; i < numParticles; ++i) {
         if (particleInvMass[i] > 0)
-            shiftedVel[i] = velData[i]+forceData[i]*(timeShift/particleInvMass[i]);
+            shiftedVel[i] = velData[i]+forceData[i]*(timeShift*particleInvMass[i]);
         else
             shiftedVel[i] = velData[i];
     }

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.h

@@ -84,7 +84,8 @@ private:
  */
 class ReferenceIntegrateDrudeLangevinStepKernel : public IntegrateDrudeLangevinStepKernel {
 public:
-    ReferenceIntegrateDrudeLangevinStepKernel(std::string name, const Platform& platform) : IntegrateDrudeLangevinStepKernel(name, platform), constraints(NULL) {
+    ReferenceIntegrateDrudeLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) :
+        IntegrateDrudeLangevinStepKernel(name, platform), data(data), constraints(NULL) {
     }
     ~ReferenceIntegrateDrudeLangevinStepKernel();
     /**
@@ -110,6 +111,7 @@ public:
      */
     double computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator);
 private:
+    ReferencePlatform::PlatformData& data;
     const DrudeForce* drudeForce;
     std::vector<int> normalParticles;
     std::vector<std::pair<int, int> > pairParticles;

plugins/drude/platforms/reference/tests/TestReferenceDrudeLangevinIntegrator.cpp

@@ -38,6 +38,7 @@
 #include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
+#include "openmm/VirtualSite.h"
 #include "openmm/DrudeForce.h"
 #include "openmm/DrudeLangevinIntegrator.h"
 #include "SimTKUtilities/SimTKOpenMMUtilities.h"
@@ -67,7 +68,7 @@ void testSinglePair() {
     vector<Vec3> positions(2);
     positions[0] = Vec3(0, 0, 0);
     positions[1] = Vec3(0, 0, 0);
-    DrudeLangevinIntegrator integ(temperature, 10.0, temperatureDrude, 10.0, 0.003);
+    DrudeLangevinIntegrator integ(temperature, 20.0, temperatureDrude, 20.0, 0.003);
     Platform& platform = Platform::getPlatformByName("Reference");
     Context context(system, integ, platform);
     context.setPositions(positions);
@@ -93,9 +94,89 @@ void testSinglePair() {
     ASSERT_USUALLY_EQUAL_TOL(3*0.5*BOLTZ*temperatureDrude, keInternal/numSteps, 0.01);
 }
 
+void testWater() {
+    // Create a box of SWM4-NDP water molecules.  This involves constraints, virtual sites,
+    // and Drude particles.
+    
+    const int gridSize = 4;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const double spacing = 0.6;
+    const double boxSize = spacing*(gridSize+1);
+    const double temperature = 300.0;
+    const double temperatureDrude = 10.0;
+    System system;
+    NonbondedForce* nonbonded = new NonbondedForce();
+    DrudeForce* drude = new DrudeForce();
+    system.addForce(nonbonded);
+    system.addForce(drude);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setCutoffDistance(1.0);
+    for (int i = 0; i < numMolecules; i++) {
+        int startIndex = system.getNumParticles();
+        system.addParticle(15.6); // O
+        system.addParticle(0.4);  // D
+        system.addParticle(1.0);  // H1
+        system.addParticle(1.0);  // H2
+        system.addParticle(0.0);  // M
+        nonbonded->addParticle(1.71636, 0.318395, 0.21094*4.184);
+        nonbonded->addParticle(-1.71636, 1, 0);
+        nonbonded->addParticle(0.55733, 1, 0);
+        nonbonded->addParticle(0.55733, 1, 0);
+        nonbonded->addParticle(-1.11466, 1, 0);
+        for (int j = 0; j < 5; j++)
+            for (int k = 0; k < j; k++)
+                nonbonded->addException(startIndex+j, startIndex+k, 0, 1, 0);
+        system.addConstraint(startIndex, startIndex+2, 0.09572);
+        system.addConstraint(startIndex, startIndex+3, 0.09572);
+        system.addConstraint(startIndex+2, startIndex+3, 0.15139);
+        system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+    }
+    vector<Vec3> positions;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                Vec3 pos(i*spacing, j*spacing, k*spacing);
+                positions.push_back(pos);
+                positions.push_back(pos);
+                positions.push_back(pos+Vec3(0.09572, 0, 0));
+                positions.push_back(pos+Vec3(-0.023999, 0.092663, 0));
+                positions.push_back(pos);
+            }
+    
+    // Simulate it and check the temperature.
+    
+    DrudeLangevinIntegrator integ(temperature, 50.0, temperatureDrude, 50.0, 0.0005);
+    Platform& platform = Platform::getPlatformByName("Reference");
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    context.applyConstraints(1e-5);
+    
+    // Equilibrate.
+    
+    integ.step(500);
+    
+    // Compute the internal and center of mass temperatures.
+    
+    double ke = 0;
+    int numSteps = 2000;
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        ke += context.getState(State::Energy).getKineticEnergy();
+    }
+    ke /= numSteps;
+    int numStandardDof = 3*3*numMolecules-system.getNumConstraints();
+    int numDrudeDof = 3*numMolecules;
+    int numDof = numStandardDof+numDrudeDof;
+    double expectedTemp = (numStandardDof*temperature+numDrudeDof*temperatureDrude)/numDof;
+    ASSERT_USUALLY_EQUAL_TOL(expectedTemp, ke/(0.5*numDof*BOLTZ), 0.02);
+}
+
 int main() {
     try {
         testSinglePair();
+        testWater();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;