Merge pull request #1522 from peastman/drude

Fixed energy error for Drude particle at same location as its parent atom

platforms/cpu/src/CpuNonbondedForce.cpp

@@ -386,13 +386,15 @@ void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex
                             float inverseR = 1/r;
                             float chargeProdOverR = scaledChargeI*posq[4*j+3]*inverseR;
                             float dEdR = chargeProdOverR*inverseR*inverseR;
-                            dEdR = dEdR * (erfAlphaR-(float)TWO_OVER_SQRT_PI*alphaR*(float)exp(-alphaR*alphaR));
+                            dEdR = dEdR * (erfAlphaR-TWO_OVER_SQRT_PI*alphaR*(float)exp(-alphaR*alphaR));
                             fvec4 result = deltaR*dEdR;
                             (fvec4(forces+4*i)-result).store(forces+4*i);
                             (fvec4(forces+4*j)+result).store(forces+4*j);
                             if (includeEnergy)
                                 threadEnergy[threadIndex] -= chargeProdOverR*erfAlphaR;
                         }
+                        else if (includeEnergy)
+                           threadEnergy[threadIndex] -= alphaEwald*TWO_OVER_SQRT_PI*scaledChargeI*posq[4*j+3];
                     }
                 }
             }

platforms/cuda/src/kernels/coulombLennardJones.cu

@@ -25,8 +25,10 @@
             tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
             tempEnergy += -prefactor*erfAlphaR;
         }
-        else
+        else {
             includeInteraction = false;
+            tempEnergy -= TWO_OVER_SQRT_PI*EWALD_ALPHA*138.935456f*posq1.w*posq2.w;
+        }
     }
     else {
 #if HAS_LENNARD_JONES

platforms/opencl/src/OpenCLContext.cpp

@@ -179,10 +179,8 @@ OpenCLContext::OpenCLContext(const System& system, int platformIndex, int device
         compilationDefines["WORK_GROUP_SIZE"] = intToString(ThreadBlockSize);
         if (platformVendor.size() >= 5 && platformVendor.substr(0, 5) == "Intel")
             defaultOptimizationOptions = "";
-        else if (platformVendor == "Apple")
-            defaultOptimizationOptions = "-cl-mad-enable -cl-no-signed-zeros";
         else
-            defaultOptimizationOptions = "-cl-fast-relaxed-math";
+            defaultOptimizationOptions = "-cl-mad-enable -cl-no-signed-zeros";
         supports64BitGlobalAtomics = (device.getInfo<CL_DEVICE_EXTENSIONS>().find("cl_khr_int64_base_atomics") != string::npos);
         supportsDoublePrecision = (device.getInfo<CL_DEVICE_EXTENSIONS>().find("cl_khr_fp64") != string::npos);
         if ((useDoublePrecision || useMixedPrecision) && !supportsDoublePrecision)

platforms/opencl/src/kernels/coulombLennardJones.cl

@@ -30,8 +30,10 @@
             tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
             tempEnergy += -prefactor*erfAlphaR;
         }
-        else
+        else {
             includeInteraction = false;
+            tempEnergy -= TWO_OVER_SQRT_PI*EWALD_ALPHA*138.935456f*posq1.w*posq2.w;
+        }
     }
     else {
 #if HAS_LENNARD_JONES

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -420,6 +420,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
     // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
 
     RealOpenMM totalExclusionEnergy = 0.0f;
+    const double TWO_OVER_SQRT_PI = 2/sqrt(PI_M);
     for (int i = 0; i < numberOfAtoms; i++)
         for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
             if (*iter > i) {
@@ -446,12 +447,15 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
                    // accumulate energies
 
                    realSpaceEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR));
+               }
+               else {
+                   realSpaceEwaldEnergy = (RealOpenMM) (alphaEwald*TWO_OVER_SQRT_PI*ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]);
+               }
 
-                   totalExclusionEnergy += realSpaceEwaldEnergy;
-                   if (energyByAtom) {
-                       energyByAtom[ii] -= realSpaceEwaldEnergy;
-                       energyByAtom[jj] -= realSpaceEwaldEnergy;
-                   }
+               totalExclusionEnergy += realSpaceEwaldEnergy;
+               if (energyByAtom) {
+                   energyByAtom[ii] -= realSpaceEwaldEnergy;
+                   energyByAtom[jj] -= realSpaceEwaldEnergy;
                }
             }
         }

plugins/drude/platforms/cuda/tests/TestCudaDrudeLangevinIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -179,6 +179,66 @@ void testWater() {
     ASSERT_USUALLY_EQUAL_TOL(expectedTemp, ke/(0.5*numDof*BOLTZ), 0.02);
 }
 
+void testForceEnergyConsistency() {
+    // Create a box of polarizable particles.
+    
+    const int gridSize = 3;
+    const int numAtoms = 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;
+    vector<Vec3> positions;
+    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::PME);
+    nonbonded->setCutoffDistance(1.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(0.9);
+    nonbonded->setEwaldErrorTolerance(5e-5);
+    for (int i = 0; i < numAtoms; i++) {
+        int startIndex = system.getNumParticles();
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(1.0, 0.3, 1.0);
+        nonbonded->addParticle(-1.0, 0.3, 1.0);
+        nonbonded->addException(startIndex, startIndex+1, 0, 1, 0);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.0, 0.001, 1, 1);
+    }
+    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);
+            }
+    
+    // Simulate it and check that force and energy remain consistent.
+    
+    DrudeLangevinIntegrator integ(temperature, 50.0, temperatureDrude, 50.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("CUDA");
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State prevState;
+    for (int i = 0; i < 100; i++) {
+        State state = context.getState(State::Energy | State::Forces | State::Positions);
+        if (i > 0) {
+            double expectedEnergyChange = 0;
+            for (int j = 0; j < system.getNumParticles(); j++) {
+                Vec3 delta = state.getPositions()[j]-prevState.getPositions()[j];
+                expectedEnergyChange -= 0.5*(state.getForces()[j]+prevState.getForces()[j]).dot(delta);
+            }
+            ASSERT_EQUAL_TOL(expectedEnergyChange, state.getPotentialEnergy()-prevState.getPotentialEnergy(), 0.05);
+        }
+        prevState = state;
+        integ.step(1);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerDrudeCudaKernelFactories();
@@ -186,6 +246,7 @@ int main(int argc, char* argv[]) {
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("Precision", string(argv[1]));
         testSinglePair();
         testWater();
+        testForceEnergyConsistency();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/drude/platforms/opencl/tests/TestOpenCLDrudeLangevinIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -179,6 +179,66 @@ void testWater() {
     ASSERT_USUALLY_EQUAL_TOL(expectedTemp, ke/(0.5*numDof*BOLTZ), 0.02);
 }
 
+void testForceEnergyConsistency() {
+    // Create a box of polarizable particles.
+    
+    const int gridSize = 3;
+    const int numAtoms = 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;
+    vector<Vec3> positions;
+    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::PME);
+    nonbonded->setCutoffDistance(1.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(0.9);
+    nonbonded->setEwaldErrorTolerance(5e-5);
+    for (int i = 0; i < numAtoms; i++) {
+        int startIndex = system.getNumParticles();
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(1.0, 0.3, 1.0);
+        nonbonded->addParticle(-1.0, 0.3, 1.0);
+        nonbonded->addException(startIndex, startIndex+1, 0, 1, 0);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.0, 0.001, 1, 1);
+    }
+    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);
+            }
+    
+    // Simulate it and check that force and energy remain consistent.
+    
+    DrudeLangevinIntegrator integ(temperature, 50.0, temperatureDrude, 50.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("OpenCL");
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State prevState;
+    for (int i = 0; i < 100; i++) {
+        State state = context.getState(State::Energy | State::Forces | State::Positions);
+        if (i > 0) {
+            double expectedEnergyChange = 0;
+            for (int j = 0; j < system.getNumParticles(); j++) {
+                Vec3 delta = state.getPositions()[j]-prevState.getPositions()[j];
+                expectedEnergyChange -= 0.5*(state.getForces()[j]+prevState.getForces()[j]).dot(delta);
+            }
+            ASSERT_EQUAL_TOL(expectedEnergyChange, state.getPotentialEnergy()-prevState.getPotentialEnergy(), 0.05);
+        }
+        prevState = state;
+        integ.step(1);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerDrudeOpenCLKernelFactories();
@@ -186,6 +246,7 @@ int main(int argc, char* argv[]) {
             Platform::getPlatformByName("OpenCL").setPropertyDefaultValue("Precision", string(argv[1]));
         testSinglePair();
         testWater();
+        testForceEnergyConsistency();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

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

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -179,11 +179,72 @@ void testWater() {
     ASSERT_USUALLY_EQUAL_TOL(expectedTemp, ke/(0.5*numDof*BOLTZ), 0.03);
 }
 
+void testForceEnergyConsistency() {
+    // Create a box of polarizable particles.
+    
+    const int gridSize = 3;
+    const int numAtoms = 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;
+    vector<Vec3> positions;
+    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::PME);
+    nonbonded->setCutoffDistance(1.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(0.9);
+    nonbonded->setEwaldErrorTolerance(5e-5);
+    for (int i = 0; i < numAtoms; i++) {
+        int startIndex = system.getNumParticles();
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        nonbonded->addParticle(1.0, 0.3, 1.0);
+        nonbonded->addParticle(-1.0, 0.3, 1.0);
+        nonbonded->addException(startIndex, startIndex+1, 0, 1, 0);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.0, 0.001, 1, 1);
+    }
+    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);
+            }
+    
+    // Simulate it and check that force and energy remain consistent.
+    
+    DrudeLangevinIntegrator integ(temperature, 50.0, temperatureDrude, 50.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("Reference");
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State prevState;
+    for (int i = 0; i < 100; i++) {
+        State state = context.getState(State::Energy | State::Forces | State::Positions);
+        if (i > 0) {
+            double expectedEnergyChange = 0;
+            for (int j = 0; j < system.getNumParticles(); j++) {
+                Vec3 delta = state.getPositions()[j]-prevState.getPositions()[j];
+                expectedEnergyChange -= 0.5*(state.getForces()[j]+prevState.getForces()[j]).dot(delta);
+            }
+            ASSERT_EQUAL_TOL(expectedEnergyChange, state.getPotentialEnergy()-prevState.getPotentialEnergy(), 0.05);
+        }
+        prevState = state;
+        integ.step(1);
+    }
+}
+
 int main() {
     try {
         registerDrudeReferenceKernelFactories();
         testSinglePair();
         testWater();
+        testForceEnergyConsistency();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;