Merge pull request #2184 from peastman/cython

Optimizations to Python code

.travis.yml

@@ -132,9 +132,9 @@ before_install:
       ${AMDAPPSDK}/bin/x86_64/clinfo;
       sudo apt-get install -y libgl1-mesa-dev;
     fi
-  # Install swig for Python wrappers. However, testing CUDA and OpenCL, we
-  # skip the Python wrapper for speed. We're not using anaconda python,
-  # but this is a fast way to get an apparently functional precompiled
+  # Install SWIG and Cython for Python wrappers. However, testing CUDA and
+  # OpenCL, we skip the Python wrapper for speed. We're not using anaconda
+  # python, but this is a fast way to get an apparently functional precompiled
   # build of swig that's more modern than what's in apt.
   - if [[ "$OPENCL" == "false" && "$CUDA" == "false" && "$TRAVIS_OS_NAME" == "linux" ]]; then
       wget http://anaconda.org/omnia/swig/3.0.7/download/linux-64/swig-3.0.7-0.tar.bz2;
@@ -142,6 +142,7 @@ before_install:
       tar -xjvf swig-3.0.7-0.tar.bz2 -C $HOME/swig;
       export PATH=$HOME/swig/bin:$PATH;
       export SWIG_LIB=$HOME/swig/share/swig/3.0.7;
+      pip install cython;
     fi
   - if [[ "$OPENCL" == "false" && "$CUDA" == "false" && "$TRAVIS_OS_NAME" == "osx" ]]; then
       wget http://anaconda.org/omnia/swig/3.0.7/download/osx-64/swig-3.0.7-0.tar.bz2;
@@ -149,6 +150,7 @@ before_install:
       tar -xjvf swig-3.0.7-0.tar.bz2 -C $HOME/swig;
       export PATH=$HOME/swig/bin:$PATH;
       export SWIG_LIB=$HOME/swig/share/swig/3.0.7;
+      sudo pip install cython;
     fi
 
   - if [[ "$CUDA" == "true" ]]; then

appveyor.yml

@@ -11,6 +11,7 @@ install:
   - "set PATH=C:\\Program Files (x86)\\Git\\bin;%PATH%"
   - pip install pytest
   - pip install numpy
+  - pip install cython
 
 # Use cclash for compiler caching (experimental)
   - ps: wget https://github.com/inorton/cclash/releases/download/0.3.14/cclash-0.3.14.zip -OutFile cclash-0.3.14.zip

devtools/ci/jenkins/install.sh

@@ -26,7 +26,7 @@ create_conda_env() {
 }
 conda remove -yn ${CONDAENV} --all --quiet || true
 create_conda_env || create_conda_env # Crappy way to work around conda concurrency restrictions
-conda install -yn ${CONDAENV} numpy scipy pytest --quiet
+conda install -yn ${CONDAENV} numpy scipy pytest cython --quiet
 source activate ${CONDAENV} # enter our new environment
 
 # Build OpenMM

docs-source/usersguide/library.rst

@@ -455,6 +455,7 @@ CMake.
    * Python 2.7 or later (http://www.python.org)
    * SWIG (http://www.swig.org)
    * Doxygen (http://www.doxygen.org)
+   * Cython (https://cython.org)
 
 * For compiling the CPU platform, you need:
 
@@ -2555,69 +2556,6 @@ Python API
 **********
 
 
-Installing the Python API
-=========================
-
-There are currently two types of packages for installing the Python API.  One
-contains wrapper source code for Unix-type machines (including Linux and Mac
-operating systems).  You will need a C++ compiler to install it using this type
-of package.  The other type of installation package is a binary package which
-contains compiled wrapper code for Windows machines (no compilers are needed to
-install binary packages).
-
-Installing on Windows
----------------------
-
-OpenMM on Windows only works with Python 3.5, so make sure that version is
-installed before you try installing. For Python installation packages and
-instructions, go to http://python.org.  We suggest that you install Python using
-the default options.
-
-Double click on the Python API Installer icon, located in the top level
-directory for the OpenMM installation (by default, this is C:\Program
-Files\OpenMM).  This will install the OpenMM package into the Python
-installation area.  If you have more than one Python installation, you will be
-asked which Python to use—make sure to select Python 3.5.
-
-Installing on Linux and Mac
----------------------------
-
-Make sure you have Python 2.7 or later installed.  For Python installation
-packages and instructions, go to http://python.org.  If you do not have the
-correct Python version, install a valid version using the default options.  Most
-versions of Linux and Mac OS X have a suitable Python preinstalled.  You can
-check by typing “\ :code:`python` |--|\ :code:`version`\ ” in a terminal window.
-
-You must have a C++ compiler to install the OpenMM Python API.  If you are using
-a Mac, install Apple's Xcode development tools
-(http://developer.apple.com/TOOLS/Xcode) to get the needed compiler.  On other
-Unix-type systems, install gcc or clang.  We recommend clang, since it produces
-faster code than gcc.
-
-The install.sh script installs the Python API automatically as part of the
-installation process, so you probably already have it installed.  If for some
-reason you need to install it manually, you can do that with the
-:code:`setup.py` script included with OpenMM.  Before executing this script,
-you must set two environment variables: :code:`OPENMM_INCLUDE_PATH` must
-point to the directory containing OpenMM header files, and
-:code:`OPENMM_LIB_PATH` must point to the directory containing OpenMM library
-files.  Assuming OpenMM is installed in the default location
-(\ :code:`/usr/local/openmm`\ ), you would type the following commands.
-Note that if you are using the system Python (as opposed to a locally installed
-version), you may need to use the :code:`sudo` command when running
-:code:`python setup.py install`\ .
-::
-
-    export OPENMM_INCLUDE_PATH=/usr/local/openmm/include
-    export OPENMM_LIB_PATH=/usr/local/openmm/lib
-    python setup.py build
-    python setup.py install
-
-If you are compiling OpenMM from source, you can also install by building the
-“PythonInstall” target:
-
-:code:`make PythonInstall` OR :code:`sudo make PythonInstall`
-
 Mapping from the C++ API to the Python API
 ==========================================
 

wrappers/python/CMakeLists.txt

@@ -49,6 +49,7 @@ foreach(SUBDIR ${SUBDIRS})
     file(GLOB_RECURSE STAGING_INPUT_FILES1 RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*README.txt"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.py"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pyx"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.i"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ini"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.sh"

wrappers/python/setup.py

@@ -8,6 +8,7 @@ import sys
 import platform
 import numpy
 from distutils.core import setup
+from Cython.Build import cythonize
 
 MAJOR_VERSION_NUM='@OPENMM_MAJOR_VERSION@'
 MINOR_VERSION_NUM='@OPENMM_MINOR_VERSION@'
@@ -217,6 +218,7 @@ def buildKeywordDictionary(major_version_num=MAJOR_VERSION_NUM,
     if platform.system() != "Windows":
         extensionArgs["runtime_library_dirs"] = library_dirs
     setupKeywords["ext_modules"] = [Extension(**extensionArgs)]
+    setupKeywords["ext_modules"] += cythonize('simtk/openmm/app/internal/*.pyx', language='c++')
 
     outputString = ''
     firstTab     = 40

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

@@ -39,13 +39,13 @@ import xml.etree.ElementTree as etree
 import math
 from math import sqrt, cos
 from copy import deepcopy
-from heapq import heappush, heappop
 from collections import defaultdict
 import simtk.openmm as mm
 import simtk.unit as unit
 from . import element as elem
 from simtk.openmm.app import Topology
 from simtk.openmm.app.internal.singleton import Singleton
+from simtk.openmm.app.internal import compiled
 
 # Directories from which to load built in force fields.
 
@@ -239,8 +239,8 @@ class ForceField(object):
                 parentDir = os.path.dirname(file)
             else:
                 parentDir = ''
-            for include in tree.getroot().findall('Include'):
-                includeFile = include.attrib['file']
+            for included in tree.getroot().findall('Include'):
+                includeFile = included.attrib['file']
                 joined = os.path.join(parentDir, includeFile)
                 if os.path.isfile(joined):
                     includeFile = joined
@@ -875,7 +875,7 @@ class ForceField(object):
         if signature in templateSignatures:
             allMatches = []
             for t in templateSignatures[signature]:
-                match = _matchResidue(res, t, bondedToAtom, ignoreExternalBonds)
+                match = compiled.matchResidueToTemplate(res, t, bondedToAtom, ignoreExternalBonds)
                 if match is not None:
                     allMatches.append((t, match))
             if len(allMatches) == 1:
@@ -1004,7 +1004,7 @@ class ForceField(object):
             if signature in signatures:
                 # Signature is the same as an existing residue; check connectivity.
                 for check_residue in unique_unmatched_residues:
-                    matches = _matchResidue(check_residue, template, bondedToAtom, False)
+                    matches = compiled.matchResidueToTemplate(check_residue, template, bondedToAtom, False)
                     if matches is not None:
                         is_unique = False
             if is_unique:
@@ -1101,7 +1101,7 @@ class ForceField(object):
                 if res in residueTemplates:
                     tname = residueTemplates[res]
                     template = self._templates[tname]
-                    matches = _matchResidue(res, template, bondedToAtom, ignoreExternalBonds)
+                    matches = compiled.matchResidueToTemplate(res, template, bondedToAtom, ignoreExternalBonds)
                     if matches is None:
                         raise Exception('User-supplied template %s does not match the residue %d (%s)' % (tname, res.index+1, res.name))
                 else:
@@ -1356,144 +1356,6 @@ def _createResidueSignature(elements):
         s += element.symbol+str(count)
     return s
 
-def _matchResidue(res, template, bondedToAtom, ignoreExternalBonds=False):
-    """Determine whether a residue matches a template and return a list of corresponding atoms.
-
-    Parameters
-    ----------
-    res : Residue
-        The residue to check
-    template : _TemplateData
-        The template to compare it to
-    bondedToAtom : list
-        Enumerates which other atoms each atom is bonded to
-    ignoreExternalBonds : bool
-        If true, ignore external bonds when matching templates
-
-    Returns
-    -------
-    list
-        a list specifying which atom of the template each atom of the residue
-        corresponds to, or None if it does not match the template
-    """
-    atoms = list(res.atoms())
-    numAtoms = len(atoms)
-    if numAtoms != len(template.atoms):
-        return None
-
-    # Translate from global to local atom indices, and record the bonds for each atom.
-
-    renumberAtoms = {}
-    for i in range(numAtoms):
-        renumberAtoms[atoms[i].index] = i
-    bondedTo = []
-    externalBonds = []
-    for atom in atoms:
-        bonds = [renumberAtoms[x] for x in bondedToAtom[atom.index] if x in renumberAtoms]
-        bondedTo.append(bonds)
-        externalBonds.append(0 if ignoreExternalBonds else len([x for x in bondedToAtom[atom.index] if x not in renumberAtoms]))
-
-    # For each unique combination of element and number of bonds, make sure the residue and
-    # template have the same number of atoms.
-
-    residueTypeCount = {}
-    for i, atom in enumerate(atoms):
-        key = (atom.element, len(bondedTo[i]), externalBonds[i])
-        if key not in residueTypeCount:
-            residueTypeCount[key] = 1
-        residueTypeCount[key] += 1
-    templateTypeCount = {}
-    for i, atom in enumerate(template.atoms):
-        key = (atom.element, len(atom.bondedTo), 0 if ignoreExternalBonds else atom.externalBonds)
-        if key not in templateTypeCount:
-            templateTypeCount[key] = 1
-        templateTypeCount[key] += 1
-    if residueTypeCount != templateTypeCount:
-        return None
-
-    # Identify template atoms that could potentially be matches for each atom.
-
-    candidates = [[] for i in range(numAtoms)]
-    for i in range(numAtoms):
-        exactNameMatch = (atoms[i].element is None and any(atom.element is None and atom.name == atoms[i].name for atom in template.atoms))
-        for j, atom in enumerate(template.atoms):
-            if (atom.element is not None and atom.element != atoms[i].element) or (exactNameMatch and atom.name != atoms[i].name):
-                continue
-            if len(atom.bondedTo) != len(bondedTo[i]):
-                continue
-            if not ignoreExternalBonds and atom.externalBonds != externalBonds[i]:
-                continue
-            candidates[i].append(j)
-
-    # Find an optimal ordering for matching atoms.  This means 1) start with the one that has the fewest options,
-    # and 2) follow with ones that are bonded to an already matched atom.
-
-    searchOrder = []
-    atomsToOrder = set(range(numAtoms))
-    efficientAtomSet = set()
-    efficientAtomHeap = []
-    while len(atomsToOrder) > 0:
-        if len(efficientAtomSet) == 0:
-            fewestNeighbors = numAtoms+1
-            for i in atomsToOrder:
-                if len(candidates[i]) < fewestNeighbors:
-                    nextAtom = i
-                    fewestNeighbors = len(candidates[i])
-        else:
-            nextAtom = heappop(efficientAtomHeap)[1]
-            efficientAtomSet.remove(nextAtom)
-        searchOrder.append(nextAtom)
-        atomsToOrder.remove(nextAtom)
-        for i in bondedTo[nextAtom]:
-            if i in atomsToOrder:
-                if i not in efficientAtomSet:
-                    efficientAtomSet.add(i)
-                    heappush(efficientAtomHeap, (len(candidates[i]), i))
-    inverseSearchOrder = [0]*numAtoms
-    for i in range(numAtoms):
-        inverseSearchOrder[searchOrder[i]] = i
-    bondedTo = [[inverseSearchOrder[bondedTo[i][j]] for j in range(len(bondedTo[i]))] for i in searchOrder]
-    candidates = [candidates[i] for i in searchOrder]
-
-    # Recursively match atoms.
-
-    matches = numAtoms*[0]
-    hasMatch = numAtoms*[False]
-    if _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, 0):
-        return [matches[inverseSearchOrder[i]] for i in range(numAtoms)]
-    return None
-
-
-def _getAtomMatchCandidates(template, bondedTo, matches, candidates, position):
-    """Get a list of template atoms that are potential matches for the next atom."""
-    for bonded in bondedTo[position]:
-        if bonded < position:
-            # This atom is bonded to another one for which we already have a match, so only consider
-            # template atoms that *that* one is bonded to.
-            return template.atoms[matches[bonded]].bondedTo
-    return candidates[position]
-
-
-def _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, position):
-    """This is called recursively from inside _matchResidue() to identify matching atoms."""
-    if position == len(matches):
-        return True
-    for i in _getAtomMatchCandidates(template, bondedTo, matches, candidates, position):
-        atom = template.atoms[i]
-        if not hasMatch[i] and i in candidates[position]:
-            # See if the bonds for this identification are consistent
-
-            allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))
-            if allBondsMatch:
-                # This is a possible match, so try matching the rest of the residue.
-
-                matches[position] = i
-                hasMatch[i] = True
-                if _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, position+1):
-                    return True
-                hasMatch[i] = False
-    return False
-
 
 def _applyPatchesToMatchResidues(forcefield, data, residues, bondedToAtom, ignoreExternalBonds):
     """Try to apply patches to find matches for residues."""
@@ -1647,14 +1509,13 @@ def _applyMultiResiduePatch(data, clusters, patch, candidateTemplates, selectedT
         except:
             # This probably means the patch is inconsistent with another one that has already been applied,
             # so just ignore it.
-            raise
             return
         newlyMatchedClusters = []
         for cluster in clusters:
             for residues in itertools.permutations(cluster):
                 residueMatches = []
                 for residue, template in zip(residues, patchedTemplates):
-                    matches = _matchResidue(residue, template, bondedToAtom, ignoreExternalBonds)
+                    matches = compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds)
                     if matches is None:
                         residueMatches = None
                         break
@@ -5339,7 +5200,7 @@ class AmoebaGeneralizedKirkwoodGenerator(object):
         if (atomicNumber in bondiMap):
             radius = bondiMap[atomicNumber]
         else:
-            outputString = "Warning no Bondi radius for atom %s of %s %d using default value=%f" % (atom.name, atom.residue.name, atom.residue.index, radius)
+            outputString = "Warning no Bondi radius for atom %s of %s %d using default value" % (atom.name, atom.residue.name, atom.residue.index)
             raise ValueError( outputString )
 
         return radius
@@ -5439,6 +5300,7 @@ class AmoebaUreyBradleyGenerator(object):
 
     def __init__(self):
 
+        self.anglesForAtom2Type = defaultdict(list)
         self.types1 = []
         self.types2 = []
         self.types3 = []
@@ -5459,11 +5321,12 @@ class AmoebaUreyBradleyGenerator(object):
         for bond in element.findall('UreyBradley'):
             types = forceField._findAtomTypes(bond.attrib, 3)
             if None not in types:
-
+                index = len(generator.types1)
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.types3.append(types[2])
-
+                for t in types[1]:
+                    generator.anglesForAtom2Type[t].append(index)
                 generator.length.append(float(bond.attrib['d']))
                 generator.k.append(float(bond.attrib['k']))
 
@@ -5491,7 +5354,7 @@ class AmoebaUreyBradleyGenerator(object):
             type1 = data.atomType[data.atoms[angle[0]]]
             type2 = data.atomType[data.atoms[angle[1]]]
             type3 = data.atomType[data.atoms[angle[2]]]
-            for i in range(len(self.types1)):
+            for i in self.anglesForAtom2Type[type2]:
                 types1 = self.types1[i]
                 types2 = self.types2[i]
                 types3 = self.types3[i]

wrappers/python/simtk/openmm/app/internal/compiled.pyx

+"""
+compiled.pyx: Utility functions that are compiled with Cython for speed
+
+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) 2018 Stanford University and the Authors.
+Authors: Peter Eastman
+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.
+"""
+from __future__ import absolute_import
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+from heapq import heappush, heappop
+
+cdef extern from "math.h":
+    double round(double x)
+    double sqrt(double x)
+
+cdef class periodicDistance:
+    """This is a callable object that computes the distance between two points, taking
+    periodic boundary conditions into account.  This is used heavily in Modeller.addSolvent().
+    """
+    cdef double vectors[3][3]
+    cdef double invBoxSize[3]
+    
+    def __init__(self, boxVectors):
+        for i in range(3):
+            for j in range(3):
+                self.vectors[i][j] = boxVectors[i][j]
+            self.invBoxSize[i] = 1.0/boxVectors[i][i]
+
+    def __call__(self, pos1, pos2):
+        cdef double dx, dy, dz, scale1, scale2, scale3
+        dx = pos1[0]-pos2[0]
+        dy = pos1[1]-pos2[1]
+        dz = pos1[2]-pos2[2]
+        scale3 = round(dz*self.invBoxSize[2])
+        dx -= scale3*self.vectors[2][0]
+        dy -= scale3*self.vectors[2][1]
+        dz -= scale3*self.vectors[2][2]
+        scale2 = round(dy*self.invBoxSize[1])
+        dx -= scale2*self.vectors[1][0]
+        dy -= scale2*self.vectors[1][1]
+        scale1 = round(dx*self.invBoxSize[0])
+        dx -= scale1*self.vectors[0][0]
+        return sqrt(dx*dx + dy*dy + dz*dz)
+
+
+def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds=False):
+    """Determine whether a residue matches a template and return a list of corresponding atoms.
+    This is used heavily in ForceField.
+
+    Parameters
+    ----------
+    res : Residue
+        The residue to check
+    template : _TemplateData
+        The template to compare it to
+    bondedToAtom : list
+        Enumerates which other atoms each atom is bonded to
+    ignoreExternalBonds : bool
+        If true, ignore external bonds when matching templates
+
+    Returns
+    -------
+    list
+        a list specifying which atom of the template each atom of the residue
+        corresponds to, or None if it does not match the template
+    """
+    cdef int numAtoms, i, j
+    atoms = list(res.atoms())
+    numAtoms = len(atoms)
+    if numAtoms != len(template.atoms):
+        return None
+
+    # Translate from global to local atom indices, and record the bonds for each atom.
+
+    renumberAtoms = {}
+    for i in range(numAtoms):
+        renumberAtoms[atoms[i].index] = i
+    bondedTo = []
+    externalBonds = []
+    for atom in atoms:
+        bonds = [renumberAtoms[x] for x in bondedToAtom[atom.index] if x in renumberAtoms]
+        bondedTo.append(bonds)
+        externalBonds.append(0 if ignoreExternalBonds else len([x for x in bondedToAtom[atom.index] if x not in renumberAtoms]))
+
+    # For each unique combination of element and number of bonds, make sure the residue and
+    # template have the same number of atoms.
+
+    residueTypeCount = {}
+    for i, atom in enumerate(atoms):
+        key = (atom.element, len(bondedTo[i]), externalBonds[i])
+        if key not in residueTypeCount:
+            residueTypeCount[key] = 1
+        residueTypeCount[key] += 1
+    templateTypeCount = {}
+    for i, atom in enumerate(template.atoms):
+        key = (atom.element, len(atom.bondedTo), 0 if ignoreExternalBonds else atom.externalBonds)
+        if key not in templateTypeCount:
+            templateTypeCount[key] = 1
+        templateTypeCount[key] += 1
+    if residueTypeCount != templateTypeCount:
+        return None
+
+    # Identify template atoms that could potentially be matches for each atom.
+
+    candidates = [[] for i in range(numAtoms)]
+    cdef bint exactNameMatch
+    for i in range(numAtoms):
+        exactNameMatch = (atoms[i].element is None and any(atom.element is None and atom.name == atoms[i].name for atom in template.atoms))
+        for j, atom in enumerate(template.atoms):
+            if (atom.element is not None and atom.element != atoms[i].element) or (exactNameMatch and atom.name != atoms[i].name):
+                continue
+            if len(atom.bondedTo) != len(bondedTo[i]):
+                continue
+            if not ignoreExternalBonds and atom.externalBonds != externalBonds[i]:
+                continue
+            candidates[i].append(j)
+
+    # Find an optimal ordering for matching atoms.  This means 1) start with the one that has the fewest options,
+    # and 2) follow with ones that are bonded to an already matched atom.
+
+    searchOrder = []
+    atomsToOrder = set(range(numAtoms))
+    efficientAtomSet = set()
+    efficientAtomHeap = []
+    while len(atomsToOrder) > 0:
+        if len(efficientAtomSet) == 0:
+            fewestNeighbors = numAtoms+1
+            for i in atomsToOrder:
+                if len(candidates[i]) < fewestNeighbors:
+                    nextAtom = i
+                    fewestNeighbors = len(candidates[i])
+        else:
+            nextAtom = heappop(efficientAtomHeap)[1]
+            efficientAtomSet.remove(nextAtom)
+        searchOrder.append(nextAtom)
+        atomsToOrder.remove(nextAtom)
+        for i in bondedTo[nextAtom]:
+            if i in atomsToOrder:
+                if i not in efficientAtomSet:
+                    efficientAtomSet.add(i)
+                    heappush(efficientAtomHeap, (len(candidates[i]), i))
+    inverseSearchOrder = [0]*numAtoms
+    for i in range(numAtoms):
+        inverseSearchOrder[searchOrder[i]] = i
+    bondedTo = [[inverseSearchOrder[bondedTo[i][j]] for j in range(len(bondedTo[i]))] for i in searchOrder]
+    candidates = [candidates[i] for i in searchOrder]
+
+    # Recursively match atoms.
+
+    matches = numAtoms*[0]
+    hasMatch = numAtoms*[False]
+    if _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, 0):
+        return [matches[inverseSearchOrder[i]] for i in range(numAtoms)]
+    return None
+
+
+def _getAtomMatchCandidates(template, bondedTo, matches, candidates, position):
+    """Get a list of template atoms that are potential matches for the next atom."""
+    for bonded in bondedTo[position]:
+        if bonded < position:
+            # This atom is bonded to another one for which we already have a match, so only consider
+            # template atoms that *that* one is bonded to.
+            return template.atoms[matches[bonded]].bondedTo
+    return candidates[position]
+
+
+def _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, int position):
+    """This is called recursively from inside matchResidueToTemplate() to identify matching atoms."""
+    if position == len(matches):
+        return True
+    cdef int i
+    for i in _getAtomMatchCandidates(template, bondedTo, matches, candidates, position):
+        atom = template.atoms[i]
+        if not hasMatch[i] and i in candidates[position]:
+            # See if the bonds for this identification are consistent
+
+            allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))
+            if allBondsMatch:
+                # This is a possible match, so try matching the rest of the residue.
+
+                matches[position] = i
+                hasMatch[i] = True
+                if _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, position+1):
+                    return True
+                hasMatch[i] = False
+    return False

wrappers/python/simtk/openmm/app/internal/pdbstructure.py

@@ -6,7 +6,7 @@ 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) 2012-2015 Stanford University and the Authors.
+Portions copyright (c) 2012-2018 Stanford University and the Authors.
 Authors: Christopher M. Bruns
 Contributors: Peter Eastman
 
@@ -158,28 +158,37 @@ class PdbStructure(object):
         for pdb_line in input_stream:
             if not isinstance(pdb_line, str):
                 pdb_line = pdb_line.decode('utf-8')
+            command = pdb_line[:6]
             # Look for atoms
-            if (pdb_line.find("ATOM  ") == 0) or (pdb_line.find("HETATM") == 0):
+            if command == "ATOM  " or command == "HETATM":
                 self._add_atom(Atom(pdb_line, self, self.extraParticleIdentifier))
+            elif command == "CONECT":
+                atoms = [int(pdb_line[6:11])]
+                for pos in (11,16,21,26):
+                    try:
+                        atoms.append(int(pdb_line[pos:pos+5]))
+                    except:
+                        pass
+                self._current_model.connects.append(atoms)
             # Notice MODEL punctuation, for the next level of detail
             # in the structure->model->chain->residue->atom->position hierarchy
-            elif (pdb_line.find("MODEL") == 0):
+            elif pdb_line[:5] == "MODEL":
                 model_number = int(pdb_line[10:14])
                 self._add_model(Model(model_number))
                 self._reset_atom_numbers()
                 self._reset_residue_numbers()
-            elif (pdb_line.find("ENDMDL") == 0):
+            elif command == "ENDMDL":
                 self._current_model._finalize()
                 if not self.load_all_models:
                     break
-            elif (pdb_line.find("END") == 0):
+            elif pdb_line[:3] == "END":
                 self._current_model._finalize()
                 if not self.load_all_models:
                     break
-            elif (pdb_line.find("TER") == 0 and pdb_line.split()[0] == "TER"):
+            elif pdb_line[:3] == "TER" and pdb_line.split()[0] == "TER":
                 self._current_model._current_chain._add_ter_record()
                 self._reset_residue_numbers()
-            elif (pdb_line.find("CRYST1") == 0):
+            elif command == "CRYST1":
                 a_length = float(pdb_line[6:15])*0.1
                 b_length = float(pdb_line[15:24])*0.1
                 c_length = float(pdb_line[24:33])*0.1
@@ -187,20 +196,12 @@ class PdbStructure(object):
                 beta = float(pdb_line[40:47])*math.pi/180.0
                 gamma = float(pdb_line[47:54])*math.pi/180.0
                 self._periodic_box_vectors = computePeriodicBoxVectors(a_length, b_length, c_length, alpha, beta, gamma)
-            elif (pdb_line.find("CONECT") == 0):
-                atoms = [int(pdb_line[6:11])]
-                for pos in (11,16,21,26):
-                    try:
-                        atoms.append(int(pdb_line[pos:pos+5]))
-                    except:
-                        pass
-                self._current_model.connects.append(atoms)
-            elif (pdb_line.find("SEQRES") == 0):
+            elif command == "SEQRES":
                 chain_id = pdb_line[11]
                 if len(self.sequences) == 0 or chain_id != self.sequences[-1].chain_id:
                     self.sequences.append(Sequence(chain_id))
                 self.sequences[-1].residues.extend(pdb_line[19:].split())
-            elif (pdb_line.find("MODRES") == 0):
+            elif command == "MODRES":
                 self.modified_residues.append(ModifiedResidue(pdb_line[16], int(pdb_line[18:22]), pdb_line[12:15].strip(), pdb_line[24:27].strip()))
         self._finalize()
 
@@ -785,11 +786,11 @@ class Atom(object):
         except:
             occupancy = 1.0
         try:
-            temperature_factor = float(pdb_line[60:66]) * unit.angstroms * unit.angstroms
+            temperature_factor = unit.Quantity(float(pdb_line[60:66]), unit.angstroms**2)
         except:
-            temperature_factor = 0.0 * unit.angstroms * unit.angstroms
+            temperature_factor = unit.Quantity(0.0, unit.angstroms**2)
         self.locations = {}
-        loc = Atom.Location(alternate_location_indicator, Vec3(x,y,z) * unit.angstroms, occupancy, temperature_factor, self.residue_name_with_spaces)
+        loc = Atom.Location(alternate_location_indicator, unit.Quantity(Vec3(x,y,z), unit.angstroms), occupancy, temperature_factor, self.residue_name_with_spaces)
         self.locations[alternate_location_indicator] = loc
         self.default_location_id = alternate_location_indicator
         # segment id, element_symbol, and formal_charge are not always present

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

@@ -35,8 +35,9 @@ __author__ = "Peter Eastman"
 __version__ = "1.0"
 
 from simtk.openmm.app import Topology, PDBFile, ForceField
-from simtk.openmm.app.forcefield import HAngles, AllBonds, CutoffNonPeriodic, CutoffPeriodic, _createResidueSignature, _matchResidue, DrudeGenerator
+from simtk.openmm.app.forcefield import HAngles, AllBonds, CutoffNonPeriodic, CutoffPeriodic, _createResidueSignature, DrudeGenerator
 from simtk.openmm.app.topology import Residue
+from simtk.openmm.app.internal import compiled
 from simtk.openmm.vec3 import Vec3
 from simtk.openmm import System, Context, NonbondedForce, CustomNonbondedForce, HarmonicBondForce, HarmonicAngleForce, VerletIntegrator, LangevinIntegrator, LocalEnergyMinimizer
 from simtk.unit import nanometer, molar, elementary_charge, amu, gram, liter, degree, sqrt, acos, is_quantity, dot, norm, kilojoules_per_mole
@@ -416,14 +417,9 @@ class Modeller(object):
             positions = self.positions.value_in_unit(nanometer)[:]
         cells = _CellList(positions, maxCutoff, vectors, True)
 
-        # Define a function to compute the distance between two points, taking periodic boundary conditions into account.
+        # Create a function to compute the distance between two points, taking periodic boundary conditions into account.
 
-        def periodicDistance(pos1, pos2):
-            delta = pos1-pos2
-            delta -= vectors[2]*floor(delta[2]*invBox[2]+0.5)
-            delta -= vectors[1]*floor(delta[1]*invBox[1]+0.5)
-            delta -= vectors[0]*floor(delta[0]*invBox[0]+0.5)
-            return norm(delta)
+        periodicDistance = compiled.periodicDistance(vectors)
 
         # Find the list of water molecules to add.
 
@@ -990,7 +986,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 _matchResidue(residue, t, bondedToAtom, False) is not None:
+                        if compiled.matchResidueToTemplate(residue, t, bondedToAtom, False) is not None:
                             matchFound = True
                 if matchFound:
                     # Just copy the residue over.
@@ -1009,7 +1005,7 @@ class Modeller(object):
                     if signature in forcefield._templateSignatures:
                         for t in forcefield._templateSignatures[signature]:
                             if t in templatesNoEP:
-                                matches = _matchResidue(residueNoEP, templatesNoEP[t], bondedToAtomNoEP, False)
+                                matches = compiled.matchResidueToTemplate(residueNoEP, templatesNoEP[t], bondedToAtomNoEP, False)
                                 if matches is not None:
                                     template = t;
                                     # Record the corresponding atoms.
@@ -1383,12 +1379,22 @@ class Modeller(object):
         context = Context(system, integrator)
         context.setPositions(mergedPositions)
         LocalEnergyMinimizer.minimize(context, 10.0, 30)
+        try:
+            import numpy as np
+            hasNumpy = True
+            proteinPosArray = np.array(proteinPos)
+            scaledProteinPosArray = np.array(scaledProteinPos)
+        except:
+            hasNumpy = False
         for i in range(50):
             weight1 = i/49.0
             weight2 = 1.0-weight1
-            mergedPositions = context.getState(getPositions=True).getPositions().value_in_unit(nanometer)
-            for j in range(len(proteinPos)):
-                mergedPositions[j+numMembraneParticles] = (weight1*proteinPos[j] + weight2*scaledProteinPos[j])
+            mergedPositions = context.getState(getPositions=True).getPositions(asNumpy=hasNumpy).value_in_unit(nanometer)
+            if hasNumpy:
+                mergedPositions[numMembraneParticles:] = weight1*proteinPosArray + weight2*scaledProteinPosArray
+            else:
+                for j in range(len(proteinPos)):
+                    mergedPositions[j+numMembraneParticles] = (weight1*proteinPos[j] + weight2*scaledProteinPos[j])
             context.setPositions(mergedPositions)
             integrator.step(20)