Skip to content

GitLab

Unverified Commit 3ef272df authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2623 from peastman/patches 

Changes to support CHARMM polarizable force field
parents 9e75b226 ac6f8942
Expand all Show whitespace changes
Inline Side-by-side
Showing with 168 additions and 185 deletions
wrappers/python/simtk/openmm/app/internal/compiled.pyx
View file @ 3ef272df
......@@ -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] = [templateAtoms.index(template.atoms[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, templateBondedTo, 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 templateBondedTo[templateAtoms[matches[bonded]]]
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, templateBondedTo, 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 @ 3ef272df
......@@ -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, False)
# 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.
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(residueNoEP.atoms(), matches):
templateAtomName = templatesNoEP[t].atoms[match].name
for atom, match in zip(atomsNoEP, matches):
templateAtomName = templateAtomsNoEP[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))
# Add the regular atoms.
......
wrappers/python/tests/TestPatches.py
View file @ 3ef272df
......@@ -32,6 +32,7 @@ class TestPatches(unittest.TestCase):
<AddExternalBond atomName="A"/>
<RemoveExternalBond atomName="C"/>
<ApplyToResidue name="RES"/>
<VirtualSite type="average2" siteName="A" atomName1="B" atomName2="B" weight1="0.8" weight2="0.2"/>
</Patch>
</Patches>
</ForceField>"""
......@@ -47,6 +48,7 @@ class TestPatches(unittest.TestCase):
self.assertEqual(1, len(patch.deletedBonds))
self.assertEqual(1, len(patch.addedExternalBonds))
self.assertEqual(1, len(patch.deletedExternalBonds))
self.assertEqual(1, len(patch.virtualSites))
self.assertEqual(1, len(ff._templatePatches))
self.assertEqual(1, len(ff._templatePatches['RES']))
self.assertEqual('A', patch.addedAtoms[0][0].name)
......@@ -67,6 +69,9 @@ class TestPatches(unittest.TestCase):
self.assertEqual(0, patch.addedExternalBonds[0].residue)
self.assertEqual('C', patch.deletedExternalBonds[0].name)
self.assertEqual(0, patch.deletedExternalBonds[0].residue)
self.assertEqual(0, patch.virtualSites[0][0].index)
self.assertEqual([1, 1], patch.virtualSites[0][0].atoms)
self.assertEqual([0.8, 0.2], patch.virtualSites[0][0].weights)
patch = list(ff._templatePatches['RES'])[0]
self.assertEqual('Test', patch[0])
self.assertEqual(0, patch[1])
......
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