Merge branch 'master' of https://github.com/peastman/openmm

openmmapi/include/openmm/Context.h

@@ -247,6 +247,15 @@ public:
      * @param stream    an input stream the checkpoint data should be read from
      */
     void loadCheckpoint(std::istream& stream);
+    /**
+     * Get a description of how the particles in the system are grouped into molecules.  Two particles are in the
+     * same molecule if they are connected by constraints or bonds, where every Force object can define bonds
+     * in whatever way are appropriate to that force.
+     *
+     * Each element lists the indices of all particles in a single molecule.  Every particle is guaranteed to
+     * belong to exactly one molecule.
+     */
+    const std::vector<std::vector<int> >& getMolecules() const;
 private:
     friend class Force;
     friend class Platform;

openmmapi/src/Context.cpp

@@ -278,3 +278,7 @@ void Context::loadCheckpoint(istream& stream) {
 ContextImpl& Context::getImpl() {
     return *impl;
 }
+
+const vector<vector<int> >& Context::getMolecules() const {
+    return impl->getMolecules();
+}

platforms/cpu/tests/TestCpuCustomNonbondedForce.cpp

@@ -843,6 +843,10 @@ void testInteractionGroupLongRangeCorrection() {
 
 int main() {
     try {
+        if (!CpuPlatform::isProcessorSupported()) {
+            cout << "CPU is not supported.  Exiting." << endl;
+            return 0;
+        }
         testSimpleExpression();
         testParameters();
         testExclusions();

platforms/reference/include/ReferenceTabulatedFunction.h

@@ -42,7 +42,7 @@ namespace OpenMM {
 /**
  * Given a TabulatedFunction, wrap it in an appropriate subclass of Lepton::CustomFunction.
  */
-extern "C" Lepton::CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function);
+extern "C" OPENMM_EXPORT Lepton::CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function);
 
 /**
  * This class adapts a Continuous1DFunction into a Lepton::CustomFunction.

platforms/reference/src/ReferenceTabulatedFunction.cpp

@@ -48,7 +48,7 @@ using namespace OpenMM;
 using namespace std;
 using Lepton::CustomFunction;
 
-extern "C" CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function) {
+extern "C" OPENMM_EXPORT CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function) {
     if (dynamic_cast<const Continuous1DFunction*>(&function) != NULL)
         return new ReferenceContinuous1DFunction(dynamic_cast<const Continuous1DFunction&>(function));
     if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)

wrappers/python/simtk/openmm/app/__init__.py

@@ -25,6 +25,7 @@ from modeller import Modeller
 from statedatareporter import StateDataReporter
 from element import Element
 from desmonddmsfile import DesmondDMSFile
+from checkpointreporter import CheckpointReporter
 
 # Enumerated values
 

wrappers/python/simtk/openmm/app/checkpointreporter.py

+"""
+checkpointreporter.py: Saves checkpoint files for a simulation
+
+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) 2014 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.
+"""
+__author__ = "Robert McGibbon"
+__version__ = "1.0"
+
+import simtk.openmm as mm
+__all__ = ['CheckpointReporter']
+
+
+class CheckpointReporter(object):
+    """CheckpointReporter saves periodic checkpoints of a simulation.
+    The checkpoints will overwrite one another -- only the last checkpoint
+    will be saved in the file.
+
+    To use it, create a CheckpointReporter, then add it to the Simulation's
+    list of reporters. To load a checkpoint file and continue a simulation,
+    use the following recipe:
+
+    >>> with open('checkput.chk', 'rb') as f:
+    >>>     simulation.context.loadCheckpoint(f.read())
+
+    Notes:
+    A checkpoint contains not only publicly visible data such as the particle
+    positions and velocities, but also internal data such as the states of
+    random number generators.  Ideally, loading a checkpoint should restore the
+    Context to an identical state to when it was written, such that continuing
+    the simulation will produce an identical trajectory.  This is not strictly
+    guaranteed to be true, however, and should not be relied on.  For most
+    purposes, however, the internal state should be close enough to be
+    reasonably considered equivalent.
+
+    A checkpoint contains data that is highly specific to the Context from
+    which it was created. It depends on the details of the System, the
+    Platform being used, and the hardware and software of the computer it was
+    created on.  If you try to load it on a computer with different hardware,
+    or for a System that is different in any way, loading is likely to fail.
+    Checkpoints created with different versions of OpenMM are also often
+    incompatible.  If a checkpoint cannot be loaded, that is signaled by
+    throwing an exception.
+
+    """
+    def __init__(self, file, reportInterval):
+        """Create a CheckpointReporter.
+
+        Parameters:
+         - file (string or open file object) The file to write to. Any current
+           contents will be overwritten.
+         - reportInterval (int) The interval (in time steps) at which to write checkpoints
+        """
+
+        self._reportInterval = reportInterval
+        if isinstance(file, basestring):
+            self._own_handle = true
+            self._out = open(file, 'w+b', 0)
+        else:
+            self._out = file
+            self._own_handle = False
+
+    def describeNextReport(self, simulation):
+        """Get information about the next report this object will generate.
+
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+
+        Returns: A five element tuple.  The first element is the number of steps until the
+        next report.  The remaining elements specify whether that report will require
+        positions, velocities, forces, and energies respectively.
+        """
+        steps = self._reportInterval - simulation.currentStep%self._reportInterval
+        return (steps, False, False, False, False)
+
+    def report(self, simulation, state):
+        """Generate a report.
+
+        Parameters:
+         - simulation (Simulation) The Simulation to generate a report for
+         - state (State) The current state of the simulation
+        """
+        self._out.seek(0)
+        chk = simulation.context.createCheckpoint()
+        self._out.write(chk)
+        self._out.truncate()
+        self._out.flush()
+
+    def __del__(self):
+        if self._own_handle:
+            self._out.close()
+

wrappers/python/simtk/openmm/app/modeller.py

@@ -39,6 +39,7 @@ from simtk.openmm.app.topology import Residue
 from simtk.openmm.vec3 import Vec3
 from simtk.openmm import System, Context, NonbondedForce, CustomNonbondedForce, HarmonicBondForce, HarmonicAngleForce, VerletIntegrator, LocalEnergyMinimizer
 from simtk.unit import nanometer, molar, elementary_charge, amu, gram, liter, degree, sqrt, acos, is_quantity, dot, norm
+import simtk.unit as unit
 import element as elem
 import os
 import random
@@ -811,7 +812,7 @@ class Modeller(object):
         else:
             context = Context(system, VerletIntegrator(0.0), platform)
         context.setPositions(newPositions)
-        LocalEnergyMinimizer.minimize(context)
+        LocalEnergyMinimizer.minimize(context, 1.0, 50)
         self.topology = newTopology
         self.positions = context.getState(getPositions=True).getPositions()
         del context
@@ -843,7 +844,9 @@ class Modeller(object):
                     if atom.element is not None:
                         newIndex[i] = index
                         index += 1
-                        newTemplate.atoms.append(ForceField._TemplateAtomData(atom.name, atom.type, atom.element))
+                        newAtom = ForceField._TemplateAtomData(atom.name, atom.type, atom.element)
+                        newAtom.externalBonds = atom.externalBonds
+                        newTemplate.atoms.append(newAtom)
                 for b1, b2 in template.bonds:
                     if b1 in newIndex and b2 in newIndex:
                         newTemplate.bonds.append((newIndex[b1], newIndex[b2]))
@@ -968,7 +971,7 @@ class Modeller(object):
                                 # and hope that energy minimization will fix it.
 
                                 knownPositions = [x for x in templateAtomPositions if x is not None]
-                                position = sum(knownPositions)/len(knownPositions)
+                                position = unit.sum(knownPositions)/len(knownPositions)
                             newPositions.append(position*nanometer)
         for bond in self.topology.bonds():
             if bond[0] in newAtoms and bond[1] in newAtoms:

wrappers/python/simtk/openmm/app/statedatareporter.py

@@ -201,9 +201,12 @@ class StateDataReporter(object):
         if self._density:
             values.append((self._totalMass/volume).value_in_unit(unit.gram/unit.item/unit.milliliter))
         if self._speed:
-            elapsedDays = (clockTime-self._initialClockTime)/86400
+            elapsedDays = (clockTime-self._initialClockTime)/86400.0
             elapsedNs = (state.getTime()-self._initialSimulationTime).value_in_unit(unit.nanosecond)
-            values.append('%.3g' % (elapsedNs/elapsedDays))
+            if elapsedDays > 0.0:
+                values.append('%.3g' % (elapsedNs/elapsedDays))
+            else:
+                values.append('--')
         if self._remainingTime:
             elapsedSeconds = clockTime-self._initialClockTime
             elapsedSteps = simulation.currentStep-self._initialSteps

wrappers/python/src/swig_doxygen/OpenMM.i

@@ -13,10 +13,10 @@ See https://simtk.org/home/pyopenmm for details"
 
 %module (docstring=DOCSTRING) openmm
 
-%include "typemaps.i"
 %include "factory.i"
 %include "std_string.i"
 %include "std_iostream.i"
+%include "typemaps.i"
 
 %include "std_map.i"
 %include "std_pair.i"

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -297,6 +297,7 @@ UNITS = {
 ("AmoebaWcaDispersionForce",              "getShctd")                                      :  ( None, ()),
 
 ("Context", "getParameter") : (None, ()),
+("Context", "getMolecules") : (None, ()),
 ("CMAPTorsionForce", "getMapParameters") : (None, ()),
 ("CMAPTorsionForce", "getTorsionParameters") : (None, ()),
 ("CMMotionRemover", "getFrequency") : (None, ()),

wrappers/python/src/swig_doxygen/swig_lib/python/typemaps.i

@@ -156,3 +156,36 @@
   $3 = &tempC;
 }
 
+%typemap(out) std::string OpenMM::Context::createCheckpoint{
+    // createCheckpoint returns a bytes object
+    $result = PyBytes_FromStringAndSize($1.c_str(), $1.length());
+}
+
+
+%typemap(in) std::string {
+    // if we have a C++ method that takes in a std::string, we're most happy to accept
+    // a python bytes object. But if the user passes in a unicode object we'll try
+    // to recover by encoding it to UTF-8 bytes
+    PyObject* temp = NULL;
+    char* c_str = NULL;
+    Py_ssize_t len = 0;
+
+    if (PyUnicode_Check($input)) {
+        temp = PyUnicode_AsUTF8String($input);
+        if (temp == NULL) {
+            SWIG_exception_fail(SWIG_TypeError, "'utf-8' codec can't decode byte");
+        }
+        PyBytes_AsStringAndSize(temp, &c_str, &len);
+        Py_XDECREF(temp);
+    } else if (PyBytes_Check($input)) {
+        PyBytes_AsStringAndSize($input, &c_str, &len);
+    } else {
+         SWIG_exception_fail(SWIG_TypeError, "argument must be str or bytes");
+    }
+
+    if (c_str == NULL) {
+        SWIG_exception_fail(SWIG_TypeError, "argument must be str or bytes");
+    }
+
+    $1 = std::string(c_str, len);
+}
\ No newline at end of file

wrappers/python/tests/TestBytes.py

+import unittest
+import simtk.openmm as mm
+
+
+class TestBytes(unittest.TestCase):
+    def test_createCheckpoint(self):
+        system = mm.System()
+        system.addParticle(1.0)
+        refPositions = [(0,0,0)]
+
+        context = mm.Context(system, mm.VerletIntegrator(0))
+        context.setPositions(refPositions)
+        chk = context.createCheckpoint()
+        # check that the return value of createCheckpoint is of type bytes (non-unicode)
+        assert isinstance(chk, bytes)
+
+        # set the positions to something random then reload the checkpoint, and
+        # make sure that the positions get restored correctly
+        context.setPositions([(12345, 12345, 123451)])
+        context.loadCheckpoint(chk)
+        newPositions = context.getState(getPositions=True).getPositions()._value
+
+        assert newPositions == refPositions
+
+        # try encoding the checkpoint in utf-8. OpenMM should be able to handle this too
+        context.setPositions([(12345, 12345, 123451)])
+        context.loadCheckpoint(chk.decode('utf-8'))
+        newPositions = context.getState(getPositions=True).getPositions()._value
+
+        assert newPositions == refPositions
+
+
+if __name__ == '__main__':
+    unittest.main()
+
+    

wrappers/python/tests/TestCheckpointReporter.py

+import unittest
+import tempfile
+import numpy as np
+from simtk.openmm import app
+import simtk.openmm as mm
+from simtk import unit
+
+
+class TestCheckpointReporter(unittest.TestCase):
+    def setUp(self):
+        with open('systems/alanine-dipeptide-implicit.pdb') as f:
+            pdb = app.PDBFile(f)
+        forcefield = app.ForceField('amber99sbildn.xml')
+        system = forcefield.createSystem(pdb.topology,
+            nonbondedMethod=app.CutoffNonPeriodic, nonbondedCutoff=1.0*unit.nanometers,
+            constraints=app.HBonds)
+        self.simulation = app.Simulation(pdb.topology, system, mm.VerletIntegrator(0.002*unit.picoseconds))
+        self.simulation.context.setPositions(pdb.positions)
+
+    def test_1(self):
+        file = tempfile.NamedTemporaryFile()
+        self.simulation.reporters.append(app.CheckpointReporter(file, 1))
+        self.simulation.step(1)
+
+        # get the current positions
+        positions = self.simulation.context.getState(getPositions=True).getPositions(asNumpy=True)._value
+        # now set the positions into junk...
+        self.simulation.context.setPositions(np.random.random(positions.shape))
+        # then reload the right positions from the checkpoint
+        with open(file.name, 'rb') as f:
+            self.simulation.context.loadCheckpoint(f.read())
+        file.close()
+
+        newPositions = self.simulation.context.getState(getPositions=True).getPositions(asNumpy=True)._value
+        np.testing.assert_array_equal(positions, newPositions)
+
+if __name__ == '__main__':
+    unittest.main()