Skip to content

Commit

Permalink
Bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@peastman
peastman committed Jul 13, 2024
1 parent ddb9fca commit 7a45e9d
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 30 deletions.
41 changes: 26 additions & 15 deletions pdbfixer/pdbfixer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -963,12 +963,23 @@ def findMissingAtoms(self):
missingAtoms = {}
missingTerminals = {}

# Record which other atoms each atom is bonded to.
# Determine which atoms have an external bond to another residue.

bondedTo = defaultdict(set)
hasExternal = defaultdict(bool)
for atom1, atom2 in self.topology.bonds():
bondedTo[atom1.name].add(atom2)
bondedTo[atom2.name].add(atom1)
if atom1.residue != atom2.residue:
hasExternal[(atom1.residue, atom1.name)] = True
hasExternal[(atom2.residue, atom2.name)] = True
for chain in self.topology.chains():
chainResidues = list(chain.residues())
for residue in chain.residues():
atomNames = [atom.name for atom in residue.atoms()]
if all(name in atomNames for name in ['C', 'O', 'CA']):
# We'll be adding peptide bonds.
if residue != chainResidues[0]:
hasExternal[(residue, 'N')] = True
if residue != chainResidues[-1]:
hasExternal[(residue, 'C')] = True

# Loop over residues.

Expand All @@ -981,23 +992,23 @@ def findMissingAtoms(self):
# If an atom is marked as terminal only, and if it is bonded to any atom that has an external bond
# to another residue, we need to omit that atom and any other terminal-only atom bonded to it.

bondedTo = defaultdict(set)
for atom1, atom2 in template.topology.bonds():
bondedTo[atom1].add(atom2)
bondedTo[atom2].add(atom1)
skip = set()
hasExternal = set()
for atom in residue.atoms():
if any(a.residue != residue for a in bondedTo[atom.name]):
hasExternal.add(atom.name)
for atom, terminal in zip(template.topology.atoms(), template.terminal):
if terminal:
for atom2 in bondedTo[atom.name]:
if atom2.name in hasExternal:
skip.add(atom.name)
for atom2 in bondedTo[atom]:
if hasExternal[(residue, atom2.name)]:
skip.add(atom)
for atom, terminal in zip(template.topology.atoms(), template.terminal):
if terminal:
for atom2 in bondedTo[atom.name]:
if atom2.name in skip:
skip.add(atom.name)
for atom2 in bondedTo[atom]:
if atom2 in skip:
skip.add(atom)
atomNames = set(atom.name for atom in residue.atoms())
templateAtoms = [atom for atom in template.topology.atoms() if atom.name not in skip]
templateAtoms = [atom for atom in template.topology.atoms() if atom not in skip]
if nucleic and residue == chainResidues[0] and (chain.index, 0) not in self.missingResidues:
templateAtoms = [atom for atom in templateAtoms if atom.name not in ('P', 'OP1', 'OP2')]

Expand Down
30 changes: 15 additions & 15 deletions pdbfixer/tests/test_build_and_simulate.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@

from __future__ import print_function
import pdbfixer
import openmm
from pdbfixer.pdbfixer import PDBFixer
from openmm import app

import os
import os.path
Expand Down Expand Up @@ -97,40 +97,38 @@ def test_build_and_simulate():
]

# DEBUG: A few small test cases.
pdbcodes_to_build = ['110D', '116D', '117D', '118D', '134D', '135D', '136D', '138D', '143D', '148D', '151D', '152D', '159D', '177D', '17RA', '183D', '184D', '186D', '187D', '188D', '189D', '1A11', '1A13', '1A1P', '1A3P', '1A51', '1A60', '1A83', '1A9L', '1AAF', '1AB1', '1ABZ', '1AC7', '1ACW', '1AD7', '1ADX', '1AFP', '1AFT', '1AFX', '1AG7', '1AGG', '1AGL', '1AGT', '1AHL', '1AIE', '1AJ1', '1AJF', '1AJJ', '1AJU', '1AKG', '1AKX', '1AL1', '1ALE', '1ALF', '1ALG', '1AM0', '1AMB', '1AMC', '1AML', '1ANP', '1ANR', '1ANS', '1AOO']
pdbcodes_to_build = ['110D', '116D', '117D', '118D', '134D', '135D', '136D', '138D', '143D', '148D', '151D', '152D', '177D', '17RA', '183D', '184D', '186D', '187D', '188D', '189D', '1A11', '1A13', '1A1P', '1A3P', '1A51', '1A60', '1A83', '1A9L', '1AAF', '1AB1', '1ABZ', '1AC7', '1ACW', '1AD7', '1ADX', '1AFP', '1AFT', '1AFX', '1AG7', '1AGG', '1AGL', '1AGT', '1AHL', '1AIE', '1AJ1', '1AJF', '1AJJ', '1AKG', '1AL1', '1ALE', '1ALF', '1ALG', '1AM0', '1AMB', '1AMC', '1AML', '1ANP', '1ANR', '1ANS', '1AOO', '1BH7', '1BX8', '1CEK']

# Don't simulate any.
pdbcodes_to_simulate = []

# Keep track of list of failures.
failures = list()

forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
for pdbcode in pdbcodes_to_build:
print("------------------------------------------------")
print(pdbcode)

output_pdb_filename = 'output.pdb'

# PDB setup parameters.
# TODO: Try several combinations?
from openmm import unit
pH = 7.0
ionic = 50.0 * unit.millimolar
box = 10.0 * unit.angstrom
positiveIon = 'Na+'
negativeIon = 'Cl-'

outfile = tempfile.NamedTemporaryFile(mode='w', delete=False)
output_pdb_filename = outfile.name

timeout_seconds = 30
watchdog = Watchdog(timeout_seconds)
build_successful = False
try:
from pdbfixer.pdbfixer import PDBFixer
from openmm import app
try:
stage = "Creating PDBFixer..."
fixer = PDBFixer(pdbid=pdbcode)
stage = "Deleting hydrogens..."
if pdbcode in ['135D', '136D', '177D', '1A83', '1AGG', '1AJ1']:
# These input files include extra hydrogens that aren't supported by the force field.
# To avoid problems, delete all pre-existing hydrogens.
modeller = app.Modeller(fixer.topology, fixer.positions)
modeller.delete([a for a in fixer.topology.atoms() if a.element == app.element.hydrogen])
fixer.topology = modeller.topology
fixer.positions = modeller.positions
stage = "Finding missing residues..."
fixer.findMissingResidues()
stage = "Finding nonstandard residues..."
Expand All @@ -147,6 +145,8 @@ def test_build_and_simulate():
fixer.addMissingHydrogens(pH)
stage = "Writing PDB file..."
app.PDBFile.writeFile(fixer.topology, fixer.positions, outfile)
stage = "Create System..."
forcefield.createSystem(fixer.topology)
stage = "Done."
outfile.close()
build_successful = True
Expand Down

0 comments on commit 7a45e9d

Please sign in to comment.