Merge pull request #1 from openmm/master

Update from upstream

.travis.yml

@@ -75,7 +75,7 @@ matrix:
 
     - sudo: false
       dist: trusty
-      python: 2.7_with_system_site_packages
+      python: 3.6
       env: ==STATIC_LIB==
            OPENCL=false
            CUDA=false
@@ -85,8 +85,8 @@ matrix:
 
     - sudo: false
       dist: trusty
-      python: 2.7_with_system_site_packages
-      env: ==PYTHON_2==
+      python: 3.6
+      env: ==PYTHON_3_6==
            OPENCL=false
            CUDA=false
            CC=$CCACHE/clang
@@ -95,9 +95,9 @@ matrix:
            CMAKE_FLAGS="-DOPENMM_GENERATE_API_DOCS=ON"
 
     - sudo: false
-      dist: trusty
-      python: 3.4
-      env: ==PYTHON_3==
+      dist: xenial
+      python: 3.8
+      env: ==PYTHON_3_8==
            OPENCL=false
            CUDA=false
            CC=$CCACHE/gcc

openmmapi/include/openmm/VirtualSite.h

@@ -112,7 +112,7 @@ public:
      * @param particle3    the index of the third particle
      * @param weight1      the weight factor (between 0 and 1) for the first particle
      * @param weight2      the weight factor (between 0 and 1) for the second particle
-     * @param weight2      the weight factor (between 0 and 1) for the third particle
+     * @param weight3      the weight factor (between 0 and 1) for the third particle
      */
     ThreeParticleAverageSite(int particle1, int particle2, int particle3, double weight1, double weight2, double weight3);
     /**

platforms/reference/include/ReferenceTabulatedFunction.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) 2014-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2014-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,6 +35,7 @@
 #include "openmm/TabulatedFunction.h"
 #include "openmm/internal/windowsExport.h"
 #include "lepton/CustomFunction.h"
+#include <memory>
 #include <vector>
 
 namespace OpenMM {
@@ -146,6 +147,22 @@ private:
     std::vector<double> values;
 };
 
+/**
+ * This is a lightweight wrapper around an immutable CustomFunction.  It makes
+ * cloning very inexpensive since nothing needs to be copied except a single
+ * pointer.
+ */
+class OPENMM_EXPORT SharedFunctionWrapper : public Lepton::CustomFunction {
+public:
+    SharedFunctionWrapper(std::shared_ptr<const CustomFunction> pointer);
+    int getNumArguments() const;
+    double evaluate(const double* arguments) const;
+    double evaluateDerivative(const double* arguments, const int* derivOrder) const;
+    CustomFunction* clone() const;
+private:
+    std::shared_ptr<const CustomFunction> pointer;
+};
+
 } // namespace OpenMM
 
 #endif /*OPENMM_REFERENCETABULATEDFUNCTION_H_*/

platforms/reference/src/ReferenceTabulatedFunction.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) 2014-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2014-2019 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -56,19 +56,22 @@ using namespace std;
 using Lepton::CustomFunction;
 
 extern "C" OPENMM_EXPORT CustomFunction* createReferenceTabulatedFunction(const TabulatedFunction& function) {
+    CustomFunction* fn;
     if (dynamic_cast<const Continuous1DFunction*>(&function) != NULL)
-        return new ReferenceContinuous1DFunction(dynamic_cast<const Continuous1DFunction&>(function));
-    if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)
-        return new ReferenceContinuous2DFunction(dynamic_cast<const Continuous2DFunction&>(function));
-    if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL)
-        return new ReferenceContinuous3DFunction(dynamic_cast<const Continuous3DFunction&>(function));
-    if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL)
-        return new ReferenceDiscrete1DFunction(dynamic_cast<const Discrete1DFunction&>(function));
-    if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL)
-        return new ReferenceDiscrete2DFunction(dynamic_cast<const Discrete2DFunction&>(function));
-    if (dynamic_cast<const Discrete3DFunction*>(&function) != NULL)
-        return new ReferenceDiscrete3DFunction(dynamic_cast<const Discrete3DFunction&>(function));
-    throw OpenMMException("createReferenceTabulatedFunction: Unknown function type");
+        fn = new ReferenceContinuous1DFunction(dynamic_cast<const Continuous1DFunction&>(function));
+    else if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)
+        fn = new ReferenceContinuous2DFunction(dynamic_cast<const Continuous2DFunction&>(function));
+    else if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL)
+        fn = new ReferenceContinuous3DFunction(dynamic_cast<const Continuous3DFunction&>(function));
+    else if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL)
+        fn = new ReferenceDiscrete1DFunction(dynamic_cast<const Discrete1DFunction&>(function));
+    else if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL)
+        fn = new ReferenceDiscrete2DFunction(dynamic_cast<const Discrete2DFunction&>(function));
+    else if (dynamic_cast<const Discrete3DFunction*>(&function) != NULL)
+        fn = new ReferenceDiscrete3DFunction(dynamic_cast<const Discrete3DFunction&>(function));
+    else
+        throw OpenMMException("createReferenceTabulatedFunction: Unknown function type");
+    return new SharedFunctionWrapper(shared_ptr<const CustomFunction>(fn));
 }
 
 ReferenceContinuous1DFunction::ReferenceContinuous1DFunction(const Continuous1DFunction& function) : function(function) {
@@ -298,3 +301,22 @@ double ReferenceDiscrete3DFunction::evaluateDerivative(const double* arguments,
 CustomFunction* ReferenceDiscrete3DFunction::clone() const {
     return new ReferenceDiscrete3DFunction(function);
 }
+
+SharedFunctionWrapper::SharedFunctionWrapper(shared_ptr<const CustomFunction> pointer) : pointer(pointer) {
+}
+
+int SharedFunctionWrapper::getNumArguments() const {
+    return pointer->getNumArguments();
+}
+
+double SharedFunctionWrapper::evaluate(const double* arguments) const {
+    return pointer->evaluate(arguments);
+}
+
+double SharedFunctionWrapper::evaluateDerivative(const double* arguments, const int* derivOrder) const {
+    return pointer->evaluateDerivative(arguments, derivOrder);
+}
+
+CustomFunction* SharedFunctionWrapper::clone() const {
+    return new SharedFunctionWrapper(pointer);
+}

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

@@ -563,20 +563,27 @@ class GromacsTopFile(object):
                     if atomName in atomReplacements:
                         atomName = atomReplacements[atomName]
 
-                    # Try to guess the element.
-
-                    upper = atomName.upper()
-                    if upper.startswith('CL'):
-                        element = elem.chlorine
-                    elif upper.startswith('NA'):
-                        element = elem.sodium
-                    elif upper.startswith('MG'):
-                        element = elem.magnesium
+                    # Try to determine the element.
+
+                    atomicNumber = self._atomTypes[fields[1]][2]
+                    if atomicNumber is None:
+                        # Try to guess the element from the name.
+                        upper = atomName.upper()
+                        if upper.startswith('CL'):
+                            element = elem.chlorine
+                        elif upper.startswith('NA'):
+                            element = elem.sodium
+                        elif upper.startswith('MG'):
+                            element = elem.magnesium
+                        else:
+                            try:
+                                element = elem.get_by_symbol(atomName[0])
+                            except KeyError:
+                                element = None
+                    elif atomicNumber == '0':
+                        element = None
                     else:
-                        try:
-                            element = elem.get_by_symbol(atomName[0])
-                        except KeyError:
-                            element = None
+                        element = elem.Element.getByAtomicNumber(int(atomicNumber))
                     atoms.append(top.addAtom(atomName, element, r))
 
                 # Add bonds to the topology

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

@@ -1004,7 +1004,7 @@ class Modeller(object):
         del context
         return actualVariants
 
-    def addExtraParticles(self, forcefield):
+    def addExtraParticles(self, forcefield, ignoreExternalBonds=False):
         """Add missing extra particles to the model that are required by a force
         field.
 
@@ -1023,6 +1023,10 @@ class Modeller(object):
         ----------
         forcefield : ForceField
             the ForceField defining what extra particles should be present
+        ignoreExternalBonds : boolean=False
+            If true, ignore external bonds when matching residues to templates.
+            This is useful when the Topology represents one piece of a larger
+            molecule, so chains are not terminated properly.
         """
         # Create copies of all residue templates that have had all extra points removed.
 
@@ -1090,7 +1094,7 @@ class Modeller(object):
                 signature = _createResidueSignature([atom.element for atom in residue.atoms()])
                 if signature in forcefield._templateSignatures:
                     for t in forcefield._templateSignatures[signature]:
-                        if compiled.matchResidueToTemplate(residue, t, bondedToAtom, False) is not None:
+                        if compiled.matchResidueToTemplate(residue, t, bondedToAtom, ignoreExternalBonds) is not None:
                             matchFound = True
                 if matchFound:
                     # Just copy the residue over.
@@ -1109,7 +1113,7 @@ class Modeller(object):
                     if signature in forcefield._templateSignatures:
                         for t in forcefield._templateSignatures[signature]:
                             if t in templatesNoEP:
-                                matches = compiled.matchResidueToTemplate(residueNoEP, templatesNoEP[t], bondedToAtomNoEP, False)
+                                matches = compiled.matchResidueToTemplate(residueNoEP, templatesNoEP[t], bondedToAtomNoEP, ignoreExternalBonds)
                                 if matches is not None:
                                     template = t;
                                     # Record the corresponding atoms.

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

@@ -48,6 +48,11 @@ try:
     have_gzip = True
 except: have_gzip = False
 
+try:
+    import numpy
+    have_numpy = True
+except: have_numpy = False
+
 class SimulatedTempering(object):
     """SimulatedTempering implements the simulated tempering algorithm for accelerated sampling.
     
@@ -214,7 +219,10 @@ class SimulatedTempering(object):
                     # Rescale the velocities.
                     
                     scale = math.sqrt(self.temperatures[j]/self.temperatures[self.currentTemperature])
-                    velocities = [v*scale for v in state.getVelocities().value_in_unit(unit.nanometers/unit.picoseconds)]
+                    if have_numpy:
+                        velocities = scale*state.getVelocities(asNumpy=True).value_in_unit(unit.nanometers/unit.picoseconds)
+                    else:
+                        velocities = [v*scale for v in state.getVelocities().value_in_unit(unit.nanometers/unit.picoseconds)]
                     self.simulation.context.setVelocities(velocities)
 
                     # Select this temperature.

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

@@ -196,10 +196,11 @@ class Simulation(object):
             while stepsToGo > 10:
                 self.integrator.step(10) # Only take 10 steps at a time, to give Python more chances to respond to a control-c.
                 stepsToGo -= 10
+                self.currentStep += 10
                 if endTime is not None and datetime.now() >= endTime:
                     return
             self.integrator.step(stepsToGo)
-            self.currentStep += nextSteps
+            self.currentStep += stepsToGo
             if anyReport:
                 # One or more reporters are ready to generate reports.  Organize them into three
                 # groups: ones that want wrapped positions, ones that want unwrapped positions,

wrappers/python/tests/TestGromacsTopFile.py

@@ -172,6 +172,13 @@ class TestGromacsTopFile(unittest.TestCase):
 
         top = GromacsTopFile('systems/bnz.top')
         gro = GromacsGroFile('systems/bnz.gro')
+        for atom in top.topology.atoms():
+            if atom.name.startswith('C'):
+                self.assertEqual(elem.carbon, atom.element)
+            elif atom.name.startswith('H'):
+                self.assertEqual(elem.hydrogen, atom.element)
+            else:
+                self.assertIsNone(atom.element)
         system = top.createSystem()
 
         self.assertEqual(26, system.getNumParticles())

wrappers/python/tests/TestSimulation.py

@@ -146,6 +146,12 @@ class TestSimulation(unittest.TestCase):
 
         self.assertTrue(endTime >= startTime+timedelta(seconds=5))
         self.assertTrue(endTime < startTime+timedelta(seconds=10))
+        
+        # Check that the time and step count are consistent.
+        
+        time = simulation.context.getState().getTime().value_in_unit(picoseconds)
+        expectedTime = simulation.currentStep*integrator.getStepSize().value_in_unit(picoseconds)
+        self.assertAlmostEqual(expectedTime, time)
 
         # Load the checkpoint and state and make sure they are both correct.