Skip to content

GitLab

Commit ffe0512f authored by peastman's avatar peastman
Browse files

addExtraParticles() works with patches

parent a752d258
Hide whitespace changes
Inline Side-by-side
Showing with 147 additions and 185 deletions
wrappers/python/simtk/openmm/app/forcefield.py
View file @ ffe0512f
......@@ -553,20 +553,30 @@ class ForceField(object):
class _SystemData(object):
"""Inner class used to encapsulate data about the system being created."""
def __init__(self):
def __init__(self, topology):
self.atomType = {}
self.atomParameters = {}
self.atomTemplateIndexes = {}
self.atoms = []
self.excludeAtomWith = []
self.atoms = list(topology.atoms())
self.excludeAtomWith = [[] for a in self.atoms]
self.virtualSites = {}
self.bonds = []
self.bonds = [ForceField._BondData(bond[0].index, bond[1].index) for bond in topology.bonds()]
self.angles = []
self.propers = []
self.impropers = []
self.atomBonds = []
self.atomBonds = [[] for a in self.atoms]
self.isAngleConstrained = []
self.constraints = {}
self.bondedToAtom = [set() for a in self.atoms]
# Record which atoms are bonded to each other atom
for i in range(len(self.bonds)):
bond = self.bonds[i]
self.bondedToAtom[bond.atom1].add(bond.atom2)
self.bondedToAtom[bond.atom2].add(bond.atom1)
self.atomBonds[bond.atom1].append(i)
self.atomBonds[bond.atom2].append(i)
def addConstraint(self, system, atom1, atom2, distance):
"""Add a constraint to the system, avoiding duplicate constraints."""
......@@ -883,8 +893,8 @@ class ForceField(object):
raise ValueError('%s: No parameters defined for atom type %s' % (self.forceName, t))
def _getResidueTemplateMatches(self, res, bondedToAtom, templateSignatures=None, ignoreExternalBonds=False):
"""Return the residue template matches, or None if none are found.
def _getResidueTemplateMatches(self, res, bondedToAtom, templateSignatures=None, ignoreExternalBonds=False, ignoreExtraParticles=False):
"""Return the templates that match a residue, or None if none are found.
Parameters
----------
......@@ -910,7 +920,7 @@ class ForceField(object):
if signature in templateSignatures:
allMatches = []
for t in templateSignatures[signature]:
match = compiled.matchResidueToTemplate(res, t, bondedToAtom, ignoreExternalBonds)
match = compiled.matchResidueToTemplate(res, t, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
if match is not None:
allMatches.append((t, match))
if len(allMatches) == 1:
......@@ -1112,83 +1122,20 @@ class ForceField(object):
"""
args['switchDistance'] = switchDistance
args['flexibleConstraints'] = flexibleConstraints
data = ForceField._SystemData()
data.atoms = list(topology.atoms())
for atom in data.atoms:
data.excludeAtomWith.append([])
data = ForceField._SystemData(topology)
rigidResidue = [False]*topology.getNumResidues()
# Make a list of all bonds
for bond in topology.bonds():
data.bonds.append(ForceField._BondData(bond[0].index, bond[1].index))
# Record which atoms are bonded to each other atom
bondedToAtom = []
for i in range(len(data.atoms)):
bondedToAtom.append(set())
data.atomBonds.append([])
for i in range(len(data.bonds)):
bond = data.bonds[i]
bondedToAtom[bond.atom1].add(bond.atom2)
bondedToAtom[bond.atom2].add(bond.atom1)
data.atomBonds[bond.atom1].append(i)
data.atomBonds[bond.atom2].append(i)
# Find the template matching each residue and assign atom types.
unmatchedResidues = []
for chain in topology.chains():
for res in chain.residues():
if res in residueTemplates:
tname = residueTemplates[res]
template = self._templates[tname]
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:
# Attempt to match one of the existing templates.
[template, matches] = self._getResidueTemplateMatches(res, bondedToAtom, ignoreExternalBonds=ignoreExternalBonds)
if matches is None:
unmatchedResidues.append(res)
else:
data.recordMatchedAtomParameters(res, template, matches)
if res.name == 'HOH':
# Determine whether this should be a rigid water.
if rigidWater is None and template is not None:
rigidResidue[res.index] = template.rigidWater
elif rigidWater:
rigidResidue[res.index] = True
# Try to apply patches to find matches for any unmatched residues.
if len(unmatchedResidues) > 0:
unmatchedResidues = _applyPatchesToMatchResidues(self, data, unmatchedResidues, bondedToAtom, ignoreExternalBonds)
# If we still haven't found a match for a residue, attempt to use residue template generators to create
# new templates (and potentially atom types/parameters).
for res in unmatchedResidues:
# A template might have been generated on an earlier iteration of this loop.
[template, matches] = self._getResidueTemplateMatches(res, bondedToAtom, ignoreExternalBonds=ignoreExternalBonds)
if matches is None:
# Try all generators.
for generator in self._templateGenerators:
if generator(self, res):
# This generator has registered a new residue template that should match.
[template, matches] = self._getResidueTemplateMatches(res, bondedToAtom, ignoreExternalBonds=ignoreExternalBonds)
if matches is None:
# Something went wrong because the generated template does not match the residue signature.
raise Exception('The residue handler %s indicated it had correctly parameterized residue %s, but the generated template did not match the residue signature.' % (generator.__class__.__name__, str(res)))
else:
# We successfully generated a residue template. Break out of the for loop.
break
if matches is None:
raise ValueError('No template found for residue %d (%s). %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
else:
data.recordMatchedAtomParameters(res, template, matches)
templateForResidue = self._matchAllResiduesToTemplates(data, topology, residueTemplates, ignoreExternalBonds)
for res in topology.residues():
if res.name == 'HOH':
# Determine whether this should be a rigid water.
if rigidWater is None:
rigidResidue[res.index] = templateForResidue[res.index].rigidWater
elif rigidWater:
rigidResidue[res.index] = True
# Create the System and add atoms
......@@ -1234,13 +1181,13 @@ class ForceField(object):
uniqueAngles = set()
for bond in data.bonds:
for atom in bondedToAtom[bond.atom1]:
for atom in data.bondedToAtom[bond.atom1]:
if atom != bond.atom2:
if atom < bond.atom2:
uniqueAngles.add((atom, bond.atom1, bond.atom2))
else:
uniqueAngles.add((bond.atom2, bond.atom1, atom))
for atom in bondedToAtom[bond.atom2]:
for atom in data.bondedToAtom[bond.atom2]:
if atom != bond.atom1:
if atom > bond.atom1:
uniqueAngles.add((bond.atom1, bond.atom2, atom))
......@@ -1252,13 +1199,13 @@ class ForceField(object):
uniquePropers = set()
for angle in data.angles:
for atom in bondedToAtom[angle[0]]:
for atom in data.bondedToAtom[angle[0]]:
if atom not in angle:
if atom < angle[2]:
uniquePropers.add((atom, angle[0], angle[1], angle[2]))
else:
uniquePropers.add((angle[2], angle[1], angle[0], atom))
for atom in bondedToAtom[angle[2]]:
for atom in data.bondedToAtom[angle[2]]:
if atom not in angle:
if atom > angle[0]:
uniquePropers.add((angle[0], angle[1], angle[2], atom))
......@@ -1268,8 +1215,8 @@ class ForceField(object):
# Make a list of all unique improper torsions
for atom in range(len(bondedToAtom)):
bondedTo = bondedToAtom[atom]
for atom in range(len(data.bondedToAtom)):
bondedTo = data.bondedToAtom[atom]
if len(bondedTo) > 2:
for subset in itertools.combinations(bondedTo, 3):
data.impropers.append((atom, subset[0], subset[1], subset[2]))
......@@ -1348,6 +1295,58 @@ class ForceField(object):
return sys
def _matchAllResiduesToTemplates(self, data, topology, residueTemplates, ignoreExternalBonds, ignoreExtraParticles=False):
"""Return a list of which template matches each residue in the topology, and assign atom types."""
templateForResidue = [None]*topology.getNumResidues()
unmatchedResidues = []
for chain in topology.chains():
for res in chain.residues():
if res in residueTemplates:
tname = residueTemplates[res]
template = self._templates[tname]
matches = compiled.matchResidueToTemplate(res, template, data.bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
if matches is None:
raise Exception('User-supplied template %s does not match the residue %d (%s)' % (tname, res.index+1, res.name))
else:
# Attempt to match one of the existing templates.
[template, matches] = self._getResidueTemplateMatches(res, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
if matches is None:
unmatchedResidues.append(res)
else:
data.recordMatchedAtomParameters(res, template, matches)
templateForResidue[res.index] = template
# Try to apply patches to find matches for any unmatched residues.
if len(unmatchedResidues) > 0:
unmatchedResidues = _applyPatchesToMatchResidues(self, data, unmatchedResidues, templateForResidue, data.bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
# If we still haven't found a match for a residue, attempt to use residue template generators to create
# new templates (and potentially atom types/parameters).
for res in unmatchedResidues:
# A template might have been generated on an earlier iteration of this loop.
[template, matches] = self._getResidueTemplateMatches(res, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
if matches is None:
# Try all generators.
for generator in self._templateGenerators:
if generator(self, res):
# This generator has registered a new residue template that should match.
[template, matches] = self._getResidueTemplateMatches(res, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
if matches is None:
# Something went wrong because the generated template does not match the residue signature.
raise Exception('The residue handler %s indicated it had correctly parameterized residue %s, but the generated template did not match the residue signature.' % (generator.__class__.__name__, str(res)))
else:
# We successfully generated a residue template. Break out of the for loop.
break
if matches is None:
raise ValueError('No template found for residue %d (%s). %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
else:
data.recordMatchedAtomParameters(res, template, matches)
templateForResidue[res.index] = template
return templateForResidue
def _findBondsForExclusions(data, sys):
"""Create a list of bonds to use when identifying exclusions."""
bondIndices = []
......@@ -1469,7 +1468,7 @@ def _createResidueSignature(elements):
return s
def _applyPatchesToMatchResidues(forcefield, data, residues, bondedToAtom, ignoreExternalBonds):
def _applyPatchesToMatchResidues(forcefield, data, residues, templateForResidue, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles):
"""Try to apply patches to find matches for residues."""
# Start by creating all templates than can be created by applying a combination of one-residue patches
# to a single template. The number of these is usually not too large, and they often cover a large fraction
......@@ -1495,11 +1494,12 @@ def _applyPatchesToMatchResidues(forcefield, data, residues, bondedToAtom, ignor
unmatchedResidues = []
for res in residues:
[template, matches] = forcefield._getResidueTemplateMatches(res, bondedToAtom, patchedTemplateSignatures, ignoreExternalBonds)
[template, matches] = forcefield._getResidueTemplateMatches(res, bondedToAtom, patchedTemplateSignatures, ignoreExternalBonds, ignoreExtraParticles)
if matches is None:
unmatchedResidues.append(res)
else:
data.recordMatchedAtomParameters(res, template, matches)
templateForResidue[res.index] = template
if len(unmatchedResidues) == 0:
return []
......@@ -1547,7 +1547,7 @@ def _applyPatchesToMatchResidues(forcefield, data, residues, bondedToAtom, ignor
for patchName in patches:
patch = forcefield._patches[patchName]
if patch.numResidues == clusterSize:
matchedClusters = _matchToMultiResiduePatchedTemplates(data, clusters, patch, patches[patchName], bondedToAtom, ignoreExternalBonds)
matchedClusters = _matchToMultiResiduePatchedTemplates(data, clusters, patch, patches[patchName], bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
for cluster in matchedClusters:
for residue in cluster:
unmatchedResidues.remove(residue)
......@@ -1597,21 +1597,21 @@ def _generatePatchedSingleResidueTemplates(template, patches, index, newTemplate
_generatePatchedSingleResidueTemplates(patchedTemplate, patches, index+1, newTemplates, newAlteredAtoms)
def _matchToMultiResiduePatchedTemplates(data, clusters, patch, residueTemplates, bondedToAtom, ignoreExternalBonds):
def _matchToMultiResiduePatchedTemplates(data, clusters, patch, residueTemplates, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles):
"""Apply a multi-residue patch to templates, then try to match them against clusters of residues."""
matchedClusters = []
selectedTemplates = [None]*patch.numResidues
_applyMultiResiduePatch(data, clusters, patch, residueTemplates, selectedTemplates, 0, matchedClusters, bondedToAtom, ignoreExternalBonds)
_applyMultiResiduePatch(data, clusters, patch, residueTemplates, selectedTemplates, 0, matchedClusters, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
return matchedClusters
def _applyMultiResiduePatch(data, clusters, patch, candidateTemplates, selectedTemplates, index, matchedClusters, bondedToAtom, ignoreExternalBonds):
def _applyMultiResiduePatch(data, clusters, patch, candidateTemplates, selectedTemplates, index, matchedClusters, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles):
"""This is called recursively to apply a multi-residue patch to all possible combinations of templates."""
if index < patch.numResidues:
for template in candidateTemplates[index]:
selectedTemplates[index] = template
_applyMultiResiduePatch(data, clusters, patch, candidateTemplates, selectedTemplates, index+1, matchedClusters, bondedToAtom, ignoreExternalBonds)
_applyMultiResiduePatch(data, clusters, patch, candidateTemplates, selectedTemplates, index+1, matchedClusters, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
else:
# We're at the deepest level of the recursion. We've selected a template for each residue, so apply the patch,
# then try to match it against clusters.
......@@ -1627,7 +1627,7 @@ def _applyMultiResiduePatch(data, clusters, patch, candidateTemplates, selectedT
for residues in itertools.permutations(cluster):
residueMatches = []
for residue, template in zip(residues, patchedTemplates):
matches = compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds)
matches = compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
if matches is None:
residueMatches = None
break
......@@ -2431,12 +2431,12 @@ class LennardJonesGenerator(object):
def registerNBFIX(self, parameters):
types = self.ff._findAtomTypes(parameters, 2)
if None not in types:
type1 = types[0][0]
type2 = types[1][0]
epsilon = _convertParameterToNumber(parameters['epsilon'])
sigma = _convertParameterToNumber(parameters['sigma'])
self.nbfixTypes[(type1, type2)] = [sigma, epsilon]
self.nbfixTypes[(type2, type1)] = [sigma, epsilon]
for type1 in types[0]:
for type2 in types[1]:
epsilon = _convertParameterToNumber(parameters['epsilon'])
sigma = _convertParameterToNumber(parameters['sigma'])
self.nbfixTypes[(type1, type2)] = [sigma, epsilon]
self.nbfixTypes[(type2, type1)] = [sigma, epsilon]
def registerLennardJones(self, parameters):
self.ljTypes.registerAtom(parameters)
......@@ -5570,11 +5570,6 @@ class AmoebaUreyBradleyGenerator(object):
generator.length.append(float(bond.attrib['d']))
generator.k.append(float(bond.attrib['k']))
else:
outputString = "AmoebaUreyBradleyGenerator: error getting types: %s %s %s" % (
bond.attrib['class1'], bond.attrib['class2'], bond.attrib['class3'])
raise ValueError(outputString)
#=============================================================================================
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
......
wrappers/python/simtk/openmm/app/internal/compiled.pyx
View file @ ffe0512f
......@@ -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) 2018 Stanford University and the Authors.
Portions copyright (c) 2018-2020 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -68,7 +68,7 @@ cdef class periodicDistance:
return sqrt(dx*dx + dy*dy + dz*dz)
def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds=False):
def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds=False, bint ignoreExtraParticles=False):
"""Determine whether a residue matches a template and return a list of corresponding atoms.
This is used heavily in ForceField.
......@@ -82,6 +82,8 @@ def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds
Enumerates which other atoms each atom is bonded to
ignoreExternalBonds : bool
If true, ignore external bonds when matching templates
ignoreExtraParticles : bool
If true, ignore extra particles (ones whose element is None) when matching templates
Returns
-------
......@@ -91,8 +93,18 @@ def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds
"""
cdef int numAtoms, i, j
atoms = list(res.atoms())
if ignoreExtraParticles:
atoms = [a for a in atoms if a.element is not None]
templateAtoms = [a for a in template.atoms if a.element is not None]
templateBondedTo = {}
for i, atom in enumerate(template.atoms):
if atom.element is not None:
templateBondedTo[atom] = [j for j in atom.bondedTo if template.atoms[j].element is not None]
else:
templateAtoms = template.atoms
templateBondedTo = dict((atom, atom.bondedTo) for atom in template.atoms)
numAtoms = len(atoms)
if numAtoms != len(template.atoms):
if numAtoms != len(templateAtoms):
return None
# Translate from global to local atom indices, and record the bonds for each atom.
......@@ -117,8 +129,8 @@ def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds
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)
for i, atom in enumerate(templateAtoms):
key = (atom.element, len(templateBondedTo[atom]), 0 if ignoreExternalBonds else atom.externalBonds)
if key not in templateTypeCount:
templateTypeCount[key] = 1
templateTypeCount[key] += 1
......@@ -130,11 +142,11 @@ def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds
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):
exactNameMatch = (atoms[i].element is None and any(atom.element is None and atom.name == atoms[i].name for atom in templateAtoms))
for j, atom in enumerate(templateAtoms):
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]):
if len(templateBondedTo[atom]) != len(bondedTo[i]):
continue
if not ignoreExternalBonds and atom.externalBonds != externalBonds[i]:
continue
......@@ -174,38 +186,38 @@ def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds
matches = numAtoms*[0]
hasMatch = numAtoms*[False]
if _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, 0):
if _findAtomMatches(templateAtoms, bondedTo, templateBondedTo, matches, hasMatch, candidates, 0):
return [matches[inverseSearchOrder[i]] for i in range(numAtoms)]
return None
def _getAtomMatchCandidates(template, bondedTo, matches, candidates, position):
def _getAtomMatchCandidates(templateAtoms, 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 templateAtoms[matches[bonded]].bondedTo
return candidates[position]
def _findAtomMatches(template, bondedTo, matches, hasMatch, candidates, int position):
def _findAtomMatches(templateAtoms, bondedTo, templateBondedTo, 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]
for i in _getAtomMatchCandidates(templateAtoms, bondedTo, matches, candidates, position):
atom = templateAtoms[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]))
allBondsMatch = all((bonded > position or matches[bonded] in templateBondedTo[atom] 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):
if _findAtomMatches(templateAtoms, bondedTo, templateBondedTo, matches, hasMatch, candidates, position+1):
return True
hasMatch[i] = False
return False
wrappers/python/simtk/openmm/app/modeller.py
View file @ ffe0512f
......@@ -1028,44 +1028,14 @@ class Modeller(object):
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.
templatesNoEP = {}
for resName, template in forcefield._templates.items():
if any(atom.element is None for atom in template.atoms):
index = 0
newIndex = {}
newTemplate = ForceField._TemplateData(resName)
for i, atom in enumerate(template.atoms):
if atom.element is not None:
newIndex[i] = index
index += 1
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]))
newTemplate.atoms[newIndex[b1]].bondedTo.append(newIndex[b2])
newTemplate.atoms[newIndex[b2]].bondedTo.append(newIndex[b1])
for b in template.externalBonds:
if b in newIndex:
newTemplate.externalBonds.append(newIndex[b])
templatesNoEP[template] = newTemplate
# Record which atoms are bonded to each other atom, with and without extra particles.
# Record which atoms are bonded to each other atom.
bondedToAtom = []
bondedToAtomNoEP = []
for atom in self.topology.atoms():
bondedToAtom.append(set())
bondedToAtomNoEP.append(set())
for atom1, atom2 in self.topology.bonds():
bondedToAtom[atom1.index].add(atom2.index)
bondedToAtom[atom2.index].add(atom1.index)
if atom1.element is not None and atom2.element is not None:
bondedToAtomNoEP[atom1.index].add(atom2.index)
bondedToAtomNoEP[atom2.index].add(atom1.index)
# If the force field has a DrudeForce, record the types of Drude particles and their parents since we'll
# need them for picking particle positions.
......@@ -1076,6 +1046,10 @@ class Modeller(object):
for type in force.typeMap:
drudeTypeMap[type] = force.typeMap[type][0]
# Identify the template to use for each residue.
templates = forcefield._matchAllResiduesToTemplates(ForceField._SystemData(self.topology), self.topology, {}, False, True)
# Create the new Topology.
newTopology = Topology()
......@@ -1087,16 +1061,8 @@ class Modeller(object):
newChain = newTopology.addChain(chain.id)
for residue in chain.residues():
newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
# Look for a matching template.
matchFound = False
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, ignoreExternalBonds) is not None:
matchFound = True
if matchFound:
template = templates[residue.index]
if len(template.atoms) == len(list(residue.atoms())):
# Just copy the residue over.
for atom in residue.atoms():
......@@ -1104,28 +1070,17 @@ class Modeller(object):
newAtoms[atom] = newAtom
newPositions.append(deepcopy(self.positions[atom.index]))
else:
# There's no matching template. Try to find one that matches based on everything except
# extra points.
template = None
residueNoEP = Residue(residue.name, residue.index, residue.chain, residue.id, residue.insertionCode)
residueNoEP._atoms = [atom for atom in residue.atoms() if atom.element is not None]
if signature in forcefield._templateSignatures:
for t in forcefield._templateSignatures[signature]:
if t in templatesNoEP:
matches = compiled.matchResidueToTemplate(residueNoEP, templatesNoEP[t], bondedToAtomNoEP, ignoreExternalBonds)
if matches is not None:
template = t;
# Record the corresponding atoms.
matchingAtoms = {}
for atom, match in zip(residueNoEP.atoms(), matches):
templateAtomName = templatesNoEP[t].atoms[match].name
for templateAtom in template.atoms:
if templateAtom.name == templateAtomName:
matchingAtoms[templateAtom] = atom
break
if template is None:
raise ValueError('Residue %d (%s) does not match any template defined by the ForceField.' % (residue.index+1, residue.name))
# Record the corresponding atoms.
matches = compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, True)
atomsNoEP = [a for a in residue.atoms() if a.element is not None]
templateAtomsNoEP = [a for a in template.atoms if a.element is not None]
matchingAtoms = {}
for atom, match in zip(atomsNoEP, matches):
templateAtomName = templateAtomsNoEP[match].name
for templateAtom in template.atoms:
if templateAtom.name == templateAtomName:
matchingAtoms[templateAtom] = atom
# Add the regular atoms.
......
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