Merge pull request #369 from peastman/rpmd

Created RPMDIntegrator::getTotalEnergy()

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

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -195,6 +195,11 @@ public:
      * and energies.  Group i will be included if (groups&(1<<i)) != 0.  The default value includes all groups.
      */
     State getState(int copy, int types, bool enforcePeriodicBox=false, int groups=0xFFFFFFFF);
+    /**
+     * Get the total energy of the ring polymer.  This includes the potential and kinetic energies of all copies,
+     * plus the potential energy of the harmonic springs that link copies together.
+     */
+    double getTotalEnergy();
     /**
      * Advance a simulation through time by taking a series of time steps.
      *

plugins/rpmd/openmmapi/src/RPMDIntegrator.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) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -34,6 +34,7 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/RpmdKernels.h"
+#include "SimTKOpenMMRealType.h"
 #include <cmath>
 #include <ctime>
 #include <string>
@@ -195,3 +196,27 @@ void RPMDIntegrator::step(int steps) {
         forcesAreValid = true;
     }
 }
+
+double RPMDIntegrator::getTotalEnergy() {
+    const System& system = owner->getSystem();
+    int numParticles = system.getNumParticles();
+    double energy = 0.0;
+    const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
+    const double wn = numCopies*BOLTZ*temperature/hbar;
+    State prevState = getState(numCopies-1, State::Positions);
+    for (int i = 0; i < numCopies; i++) {
+        // Add the energy of this copy.
+        
+        State state = getState(i, State::Positions | State::Energy);
+        energy += state.getKineticEnergy()+state.getPotentialEnergy();
+        
+        // Add the energy from the springs connecting it to the previous copy.
+        
+        for (int j = 0; j < numParticles; j++) {
+            Vec3 delta = state.getPositions()[j]-prevState.getPositions()[j];
+            energy += 0.5*wn*wn*system.getParticleMass(j)*delta.dot(delta);
+        }
+        prevState = state;
+    }
+    return energy;
+}

plugins/rpmd/platforms/cuda/tests/TestCudaRpmd.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) 2011-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -466,30 +466,9 @@ void testWithoutThermostat() {
     
     double initialEnergy;
     int numSteps = 100;
-    const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
-    const double wn = numCopies*BOLTZ*temperature/hbar;
-    const double springConstant = mass*wn*wn;
     for (int i = 0; i < numSteps; i++) {
         integ.step(1);
-        
-        // Sum the energies of all the copies.
-        
-        double energy = 0.0;
-        for (int j = 0; j < numCopies; j++) {
-            State state = integ.getState(j, State::Positions | State::Energy);
-            positions[j] = state.getPositions();
-            energy += state.getPotentialEnergy()+state.getKineticEnergy();
-        }
-        
-        // Add the energy from the springs connecting copies.
-        
-        for (int j = 0; j < numCopies; j++) {
-            int previous = (j == 0 ? numCopies-1 : j-1);
-            for (int k = 0; k < numParticles; k++) {
-                Vec3 delta = positions[j][k]-positions[previous][k];
-                energy += 0.5*springConstant*delta.dot(delta);
-            }
-        }
+        double energy = integ.getTotalEnergy();
         if (i == 0)
             initialEnergy = energy;
         else

plugins/rpmd/platforms/opencl/tests/TestOpenCLRpmd.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) 2011-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -467,30 +467,9 @@ void testWithoutThermostat() {
     
     double initialEnergy;
     int numSteps = 100;
-    const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
-    const double wn = numCopies*BOLTZ*temperature/hbar;
-    const double springConstant = mass*wn*wn;
     for (int i = 0; i < numSteps; i++) {
         integ.step(1);
-        
-        // Sum the energies of all the copies.
-        
-        double energy = 0.0;
-        for (int j = 0; j < numCopies; j++) {
-            State state = integ.getState(j, State::Positions | State::Energy);
-            positions[j] = state.getPositions();
-            energy += state.getPotentialEnergy()+state.getKineticEnergy();
-        }
-        
-        // Add the energy from the springs connecting copies.
-        
-        for (int j = 0; j < numCopies; j++) {
-            int previous = (j == 0 ? numCopies-1 : j-1);
-            for (int k = 0; k < numParticles; k++) {
-                Vec3 delta = positions[j][k]-positions[previous][k];
-                energy += 0.5*springConstant*delta.dot(delta);
-            }
-        }
+        double energy = integ.getTotalEnergy();
         if (i == 0)
             initialEnergy = energy;
         else

plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.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) 2011-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -350,30 +350,9 @@ void testWithoutThermostat() {
     
     double initialEnergy;
     int numSteps = 100;
-    const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
-    const double wn = numCopies*BOLTZ*temperature/hbar;
-    const double springConstant = mass*wn*wn;
     for (int i = 0; i < numSteps; i++) {
         integ.step(1);
-        
-        // Sum the energies of all the copies.
-        
-        double energy = 0.0;
-        for (int j = 0; j < numCopies; j++) {
-            State state = integ.getState(j, State::Positions | State::Energy);
-            positions[j] = state.getPositions();
-            energy += state.getPotentialEnergy()+state.getKineticEnergy();
-        }
-        
-        // Add the energy from the springs connecting copies.
-        
-        for (int j = 0; j < numCopies; j++) {
-            int previous = (j == 0 ? numCopies-1 : j-1);
-            for (int k = 0; k < numParticles; k++) {
-                Vec3 delta = positions[j][k]-positions[previous][k];
-                energy += 0.5*springConstant*delta.dot(delta);
-            }
-        }
+        double energy = integ.getTotalEnergy();
         if (i == 0)
             initialEnergy = energy;
         else

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -424,5 +424,6 @@ UNITS = {
 ("DrudeLangevinIntegrator", "getDrudeFriction") : ("1/unit.picosecond", ()),
 ("DrudeSCFIntegrator", "getMinimizationErrorTolerance") : ("unit.kilojoules_per_mole/unit.nanometer", ()),
 ("RPMDIntegrator", "getContractions") : (None, ()),
+("RPMDIntegrator", "getTotalEnergy") : ("unit.kilojoules_per_mole", ()),
 }