Merge pull request #1275 from peastman/enforceperiodic

Fixed bug in getState() when enforcing periodic boundary conditions

openmmapi/src/Context.cpp

@@ -115,9 +115,9 @@ State Context::getState(int types, bool enforcePeriodicBox, int groups) const {
 
                 // Find the displacement to move it into the first periodic box.
                 Vec3 diff;
-                diff -= periodicBoxSize[0]*static_cast<int>(center[0]/periodicBoxSize[0][0]);
-                diff -= periodicBoxSize[1]*static_cast<int>(center[1]/periodicBoxSize[1][1]);
-                diff -= periodicBoxSize[2]*static_cast<int>(center[2]/periodicBoxSize[2][2]);
+                diff += periodicBoxSize[2]*floor(center[2]/periodicBoxSize[2][2]);
+                diff += periodicBoxSize[1]*floor((center[1]-diff[1])/periodicBoxSize[1][1]);
+                diff += periodicBoxSize[0]*floor((center[0]-diff[0])/periodicBoxSize[0][0]);
 
                 // Translate all the particles in the molecule.
                 for (int j = 0; j < (int) molecules[i].size(); j++) {

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -129,9 +129,9 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
 
             // Find the displacement to move it into the first periodic box.
             Vec3 diff;
-            diff -= periodicBoxSize[0]*static_cast<int>(center[0]/periodicBoxSize[0][0]);
-            diff -= periodicBoxSize[1]*static_cast<int>(center[1]/periodicBoxSize[1][1]);
-            diff -= periodicBoxSize[2]*static_cast<int>(center[2]/periodicBoxSize[2][2]);
+            diff += periodicBoxSize[2]*floor(center[2]/periodicBoxSize[2][2]);
+            diff += periodicBoxSize[1]*floor((center[1]-diff[1])/periodicBoxSize[1][1]);
+            diff += periodicBoxSize[0]*floor((center[0]-diff[0])/periodicBoxSize[0][0]);
 
             // Translate all the particles in the molecule.
             for (int j = 0; j < (int) molecules[i].size(); j++) {

tests/TestEnforcePeriodicBox.cpp

@@ -41,7 +41,7 @@ using namespace OpenMM;
 using namespace std;
 
 void testTruncatedOctahedron() {
-    const int numMolecules = 5;
+    const int numMolecules = 50;
     const int numParticles = numMolecules*2;
     const float cutoff = 2.0;
     Vec3 a(6.7929, 0, 0);
@@ -63,7 +63,7 @@ void testTruncatedOctahedron() {
         system.addParticle(1.0);
         force->addParticle(-1, 0.2, 0.2);
         force->addParticle(1, 0.2, 0.2);
-        positions[2*i] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        positions[2*i] = a*(5*genrand_real2(sfmt)-2) + b*(5*genrand_real2(sfmt)-2) + c*(5*genrand_real2(sfmt)-2);
         positions[2*i+1] = positions[2*i] + Vec3(1.0, 0.0, 0.0);
         system.addConstraint(2*i, 2*i+1, 1.0);
     }
@@ -73,6 +73,15 @@ void testTruncatedOctahedron() {
     Context context(system, integrator, Platform::getPlatformByName("Reference"));
     context.setPositions(positions);
     State initialState = context.getState(State::Positions | State::Energy, true);
+    for (int i = 0; i < numMolecules; i++) {
+        Vec3 center = (initialState.getPositions()[2*i]+initialState.getPositions()[2*i+1])*0.5;
+        ASSERT(center[0] >= 0.0);
+        ASSERT(center[1] >= 0.0);
+        ASSERT(center[2] >= 0.0);
+        ASSERT(center[0] <= a[0]);
+        ASSERT(center[1] <= b[1]);
+        ASSERT(center[2] <= c[2]);
+    }
     double initialEnergy = initialState.getPotentialEnergy();
 
     context.setState(initialState);