Merge pull request #1234 from rmcgibbo/enforcePeriodicBox

getState enforcePeriodicBox with triclinic unit cells

openmmapi/src/Context.cpp

@@ -114,21 +114,15 @@ State Context::getState(int types, bool enforcePeriodicBox, int groups) const {
                 center *= 1.0/molecules[i].size();
 
                 // Find the displacement to move it into the first periodic box.
-
-                int xcell = (int) floor(center[0]/periodicBoxSize[0][0]);
-                int ycell = (int) floor(center[1]/periodicBoxSize[1][1]);
-                int zcell = (int) floor(center[2]/periodicBoxSize[2][2]);
-                double dx = xcell*periodicBoxSize[0][0];
-                double dy = ycell*periodicBoxSize[1][1];
-                double dz = zcell*periodicBoxSize[2][2];
+                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]);
 
                 // Translate all the particles in the molecule.
-                
                 for (int j = 0; j < (int) molecules[i].size(); j++) {
                     Vec3& pos = positions[molecules[i][j]];
-                    pos[0] -= dx;
-                    pos[1] -= dy;
-                    pos[2] -= dz;
+                    pos -= diff;
                 }
             }
         }

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -128,21 +128,15 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
             center *= 1.0/molecules[i].size();
 
             // Find the displacement to move it into the first periodic box.
-
-            int xcell = (int) floor(center[0]/periodicBoxSize[0][0]);
-            int ycell = (int) floor(center[1]/periodicBoxSize[1][1]);
-            int zcell = (int) floor(center[2]/periodicBoxSize[2][2]);
-            double dx = xcell*periodicBoxSize[0][0];
-            double dy = ycell*periodicBoxSize[1][1];
-            double dz = zcell*periodicBoxSize[2][2];
+            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]);
 
             // Translate all the particles in the molecule.
-
             for (int j = 0; j < (int) molecules[i].size(); j++) {
                 Vec3& pos = positions[molecules[i][j]];
-                pos[0] -= dx;
-                pos[1] -= dy;
-                pos[2] -= dz;
+                pos -= diff;
             }
         }
 
@@ -170,7 +164,7 @@ double RPMDIntegrator::computeKineticEnergy() {
 
 void RPMDIntegrator::step(int steps) {
     if (context == NULL)
-        throw OpenMMException("This Integrator is not bound to a context!");    
+        throw OpenMMException("This Integrator is not bound to a context!");
     if (!hasSetPosition) {
         // Initialize the positions from the context.
 

tests/TestEnforcePeriodicBox.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
+ * Authors: Robert McGibbon                                                   *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/Platform.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+
+using namespace OpenMM;
+using namespace std;
+
+void testTruncatedOctahedron() {
+    const int numMolecules = 5;
+    const int numParticles = numMolecules*2;
+    const float cutoff = 2.0;
+    Vec3 a(6.7929, 0, 0);
+    Vec3 b(-2.264163559406279, 6.404455775962287, 0);
+    Vec3 c(-2.264163559406279, -3.2019384603140684, 5.54658849047036);
+        
+    System system;
+    system.setDefaultPeriodicBoxVectors(a, b, c);
+    NonbondedForce* force = new NonbondedForce();
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    
+    force->setCutoffDistance(cutoff);
+    force->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        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+1] = positions[2*i] + Vec3(1.0, 0.0, 0.0);
+        system.addConstraint(2*i, 2*i+1, 1.0);
+    }
+    system.addForce(force);
+    
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    context.setPositions(positions);
+    State initialState = context.getState(State::Positions | State::Energy, true);
+    double initialEnergy = initialState.getPotentialEnergy();
+
+    context.setState(initialState);
+    State finalState = context.getState(State::Positions | State::Energy, true);
+    double finalEnergy = finalState.getPotentialEnergy();
+
+    ASSERT_EQUAL_TOL(initialEnergy, finalEnergy, 1e-4);
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        testTruncatedOctahedron();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}