Skip to content

GitLab

Commit 2a6a5edc authored by peastman's avatar peastman
Browse files

Optimizations to Python code

parent 76780eef
Hide whitespace changes
Inline Side-by-side
Showing with 267 additions and 183 deletions
wrappers/python/CMakeLists.txt
View file @ 2a6a5edc
......@@ -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
View file @ 2a6a5edc
......@@ -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
View file @ 2a6a5edc
......@@ -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 0 → 100644
View file @ 2a6a5edc
"""
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
View file @ 2a6a5edc
......@@ -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
View file @ 2a6a5edc
......@@ -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)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment