Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Improved method of selecting positions for added residues #198

Closed
wants to merge 1 commit into from
Closed

[WIP] Improved method of selecting positions for added residues #198

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 55 additions & 3 deletions pdbfixer/pdbfixer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.

Portions copyright (c) 2013-2015 Stanford University and the Authors.
Portions copyright (c) 2013-2020 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:

Expand Down Expand Up @@ -374,6 +374,7 @@ def _addAtomsToTopology(self, heavyAtomsOnly, omitUnknownMolecules):
addedAtomMap = {}
addedOXT = []
residueCenters = [self._computeResidueCenter(res).value_in_unit(unit.nanometers) for res in self.topology.residues()]*unit.nanometers
addedResidues = []
for chain in self.topology.chains():
if omitUnknownMolecules and not any(residue.name in self.templates for residue in chain.residues()):
continue
Expand All @@ -397,7 +398,7 @@ def _addAtomsToTopology(self, heavyAtomsOnly, omitUnknownMolecules):
startPosition = endPosition+outward
loopDirection = None
firstIndex = int(residue.id)-len(insertHere)
self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, loopDirection, residue, newAtoms, newPositions, firstIndex)
addedResidues += self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, loopDirection, residue, newAtoms, newPositions, firstIndex)

# Create the new residue and add existing heavy atoms.

Expand Down Expand Up @@ -452,7 +453,7 @@ def _addAtomsToTopology(self, heavyAtomsOnly, omitUnknownMolecules):
outward *= len(insertHere)*0.5*unit.nanometer/norm
endPosition = startPosition+outward
firstIndex = int(residue.id)+1
self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, None, residue, newAtoms, newPositions, firstIndex)
addedResidues += self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, None, residue, newAtoms, newPositions, firstIndex)
newResidue = list(newChain.residues())[-1]
if newResidue.name in proteinResidues:
terminalsToAdd = ['OXT']
Expand Down Expand Up @@ -484,6 +485,8 @@ def _addAtomsToTopology(self, heavyAtomsOnly, omitUnknownMolecules):

# Return the results.

if len(addedResidues) > 0:
newPositions = self._optimizedAddedResiduePositions(newTopology, newPositions, addedResidues)
return (newTopology, newPositions, newAtoms, existingAtomMap)

def _computeResidueCenter(self, residue):
Expand All @@ -507,6 +510,7 @@ def _addMissingResiduesToChain(self, chain, residueNames, startPosition, endPosi

# Add the residues.

addedResidues = []
for i, residueName in enumerate(residueNames):
template = self.templates[residueName]

Expand All @@ -523,6 +527,7 @@ def _addMissingResiduesToChain(self, chain, residueNames, startPosition, endPosi
# Create the new residue.

newResidue = chain.topology.addResidue(residueName, chain, "%d" % ((firstIndex+i)%10000))
addedResidues.append(newResidue)
fraction = (i+1.0)/(numResidues+1.0)
translate = startPosition + (endPosition-startPosition)*fraction + loopHeight*math.sin(fraction*math.pi)*loopDirection
templateAtoms = list(template.topology.atoms())
Expand All @@ -533,6 +538,53 @@ def _addMissingResiduesToChain(self, chain, residueNames, startPosition, endPosi
newAtoms.append(newAtom)
templatePosition = template.positions[atom.index].value_in_unit(unit.nanometer)
newPositions.append(mm.Vec3(*np.dot(rotate, templatePosition))*unit.nanometer+translate)
return addedResidues

def _optimizedAddedResiduePositions(self, topology, positions, addedResidues):
"""Optimize the positions of added residues by minimizing a coarse grained model."""
# Find the position and size of each residue.

residueCenters = []
residueSizes = []
positions = positions.value_in_unit(unit.nanometer)
addedResidues = set(addedResidues)
for res in topology.residues():
numAtoms = len(list(res.atoms()))
center = unit.sum([positions[atom.index] for atom in res.atoms()])/numAtoms
radius = np.sqrt(unit.sum([unit.norm(positions[atom.index]-center)**2 for atom in res.atoms()])/numAtoms) + 0.1
residueCenters.append(center)
residueSizes.append(radius)

# Construct a coarse grained model for optimizing the positions.

system = mm.System()
nb = mm.CustomNonbondedForce('1/((r/(radius1+radius2))^4+1)')
nb.addPerParticleParameter('radius')
system.addForce(nb)
bonded = mm.HarmonicBondForce()
system.addForce(bonded)
prevRes = None
for res, center, radius in zip(topology.residues(), residueCenters, residueSizes):
system.addParticle(1.0 if res in addedResidues else 0.0)
nb.addParticle([radius])
if prevRes is not None and res.chain == prevRes.chain:
dist = unit.norm(center-residueCenters[prevRes.index])
bonded.addBond(res.index, prevRes.index, dist, 100.0)
prevRes = res

# Run the optimization and update the positions of added residues.

integrator = mm.VerletIntegrator(1*unit.femtosecond)
context = mm.Context(system, integrator)
context.setPositions(residueCenters)
mm.LocalEnergyMinimizer.minimize(context, tolerance=0.1)
newCenters = context.getState(getPositions=True).getPositions().value_in_unit(unit.nanometer)
newPositions = positions[:]
for res in addedResidues:
delta = newCenters[res.index] - residueCenters[res.index]
for atom in res.atoms():
newPositions[atom.index] += delta
return newPositions*unit.nanometer

def removeChains(self, chainIndices=None, chainIds=None):
"""Remove a set of chains from the structure.
Expand Down