Skip to content

GitLab

Commit dfdf8328 authored by Peter Eastman's avatar Peter Eastman
Browse files

Implementing addHydrogens()

parent 4ae53b55
Hide whitespace changes
Inline Side-by-side
Showing with 671 additions and 1 deletion
wrappers/python/simtk/openmm/app/data/hydrogens.xml 0 → 100644
View file @ dfdf8328
<Residues>
<Residue name="A">
<H name="H1'" parent="C1'"/>
<H name="H2" parent="C2"/>
<H name="H2'" parent="C2'"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H61" parent="N6"/>
<H name="H62" parent="N6"/>
<H name="H8" parent="C8"/>
<H name="HO2'" parent="O2'"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="ACE">
<H name="H" parent="C"/>
<H name="H1" parent="CH3"/>
<H name="H2" parent="CH3"/>
<H name="H3" parent="CH3"/>
</Residue>
<Residue name="AIB">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HB11" parent="CB1"/>
<H name="HB12" parent="CB1"/>
<H name="HB13" parent="CB1"/>
<H name="HB21" parent="CB2"/>
<H name="HB22" parent="CB2"/>
<H name="HB23" parent="CB2"/>
</Residue>
<Residue name="ALA">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB1" parent="CB"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
</Residue>
<Residue name="ARG">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD2" parent="CD"/>
<H name="HD3" parent="CD"/>
<H name="HE" parent="NE"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
<H name="HH11" parent="NH1"/>
<H name="HH12" parent="NH1"/>
<H name="HH21" parent="NH2"/>
<H name="HH22" parent="NH2"/>
</Residue>
<Residue name="ASN">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD21" parent="ND2"/>
<H name="HD22" parent="ND2"/>
</Residue>
<Residue name="ASP">
<Variant name="ASP"/>
<Variant name="ASH"/>
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD2" parent="OD2" maxph="4.4" variant="ASH"/>
</Residue>
<Residue name="C">
<H name="H1'" parent="C1'"/>
<H name="H2'" parent="C2'"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H41" parent="N4"/>
<H name="H42" parent="N4"/>
<H name="H5" parent="C5"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H6" parent="C6"/>
<H name="HO2'" parent="O2'"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="CYS">
<Variant name="CYS"/>
<Variant name="CYX"/>
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HG" parent="SG" maxph="8.5" variant="CYS"/>
</Residue>
<Residue name="DA">
<H name="H1'" parent="C1'"/>
<H name="H2" parent="C2"/>
<H name="H2'" parent="C2'"/>
<H name="H2''" parent="C2'"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H61" parent="N6"/>
<H name="H62" parent="N6"/>
<H name="H8" parent="C8"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="DC">
<H name="H1'" parent="C1'"/>
<H name="H2'" parent="C2'"/>
<H name="H2''" parent="C2'"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H41" parent="N4"/>
<H name="H42" parent="N4"/>
<H name="H5" parent="C5"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H6" parent="C6"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="DG">
<H name="H1" parent="N1"/>
<H name="H1'" parent="C1'"/>
<H name="H2'" parent="C2'"/>
<H name="H2''" parent="C2'"/>
<H name="H21" parent="N2"/>
<H name="H22" parent="N2"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H8" parent="C8"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="DT">
<H name="H1'" parent="C1'"/>
<H name="H2'" parent="C2'"/>
<H name="H2''" parent="C2'"/>
<H name="H3" parent="N3"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H6" parent="C6"/>
<H name="H71" parent="C7"/>
<H name="H72" parent="C7"/>
<H name="H73" parent="C7"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="FOR">
<H name="H1" parent="C"/>
<H name="H2" parent="C"/>
</Residue>
<Residue name="G">
<H name="H1" parent="N1"/>
<H name="H1'" parent="C1'"/>
<H name="H2'" parent="C2'"/>
<H name="H21" parent="N2"/>
<H name="H22" parent="N2"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H8" parent="C8"/>
<H name="HO2'" parent="O2'"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="GLN">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HE21" parent="NE2"/>
<H name="HE22" parent="NE2"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
</Residue>
<Residue name="GLU">
<Variant name="GLU"/>
<Variant name="GLH"/>
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HE2" parent="OE2" maxph="4.4" variant="GLH"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
</Residue>
<Residue name="GLY">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA2" parent="CA"/>
<H name="HA3" parent="CA"/>
</Residue>
<Residue name="HIS">
<Variant name="HID"/>
<Variant name="HIE"/>
<Variant name="HIP"/>
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD1" parent="ND1" variant="HID,HIP"/>
<H name="HD2" parent="CD2"/>
<H name="HE1" parent="CE1"/>
<H name="HE2" parent="NE2" variant="HIE,HIP"/>
</Residue>
<Residue name="HOH">
<H name="H1" parent="O"/>
<H name="H2" parent="O"/>
</Residue>
<Residue name="ILE">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB" parent="CB"/>
<H name="HD11" parent="CD1"/>
<H name="HD12" parent="CD1"/>
<H name="HD13" parent="CD1"/>
<H name="HG12" parent="CG1"/>
<H name="HG13" parent="CG1"/>
<H name="HG21" parent="CG2"/>
<H name="HG22" parent="CG2"/>
<H name="HG23" parent="CG2"/>
</Residue>
<Residue name="LEU">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD11" parent="CD1"/>
<H name="HD12" parent="CD1"/>
<H name="HD13" parent="CD1"/>
<H name="HD21" parent="CD2"/>
<H name="HD22" parent="CD2"/>
<H name="HD23" parent="CD2"/>
<H name="HG" parent="CG"/>
</Residue>
<Residue name="LYS">
<Variant name="LYS"/>
<Variant name="LYN"/>
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD2" parent="CD"/>
<H name="HD3" parent="CD"/>
<H name="HE2" parent="CE"/>
<H name="HE3" parent="CE"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
<H name="HZ1" parent="NZ"/>
<H name="HZ2" parent="NZ"/>
<H name="HZ3" parent="NZ" maxph="10.0" variant="LYS"/>
</Residue>
<Residue name="MET">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HE1" parent="CE"/>
<H name="HE2" parent="CE"/>
<H name="HE3" parent="CE"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
</Residue>
<Residue name="NH2">
<H name="HN1" parent="N"/>
<H name="HN2" parent="N"/>
</Residue>
<Residue name="NME">
<H name="H" parent="N"/>
<H name="H1" parent="C"/>
<H name="H2" parent="C"/>
<H name="H3" parent="C"/>
<H name="HN2" parent="N"/>
</Residue>
<Residue name="ORN">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD2" parent="CD"/>
<H name="HD3" parent="CD"/>
<H name="HE1" parent="NE"/>
<H name="HE2" parent="NE"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
</Residue>
<Residue name="PCA">
<H name="H" parent="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
</Residue>
<Residue name="PHE">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD1" parent="CD1"/>
<H name="HD2" parent="CD2"/>
<H name="HE1" parent="CE1"/>
<H name="HE2" parent="CE2"/>
<H name="HZ" parent="CZ"/>
</Residue>
<Residue name="PRO">
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD2" parent="CD"/>
<H name="HD3" parent="CD"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
</Residue>
<Residue name="SER">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HG" parent="OG"/>
</Residue>
<Residue name="THR">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB" parent="CB"/>
<H name="HG1" parent="OG1"/>
<H name="HG21" parent="CG2"/>
<H name="HG22" parent="CG2"/>
<H name="HG23" parent="CG2"/>
</Residue>
<Residue name="TRP">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD1" parent="CD1"/>
<H name="HE1" parent="NE1"/>
<H name="HE3" parent="CE3"/>
<H name="HH2" parent="CH2"/>
<H name="HZ2" parent="CZ2"/>
<H name="HZ3" parent="CZ3"/>
</Residue>
<Residue name="TYR">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB2" parent="CB"/>
<H name="HB3" parent="CB"/>
<H name="HD1" parent="CD1"/>
<H name="HD2" parent="CD2"/>
<H name="HE1" parent="CE1"/>
<H name="HE2" parent="CE2"/>
<H name="HH" parent="OH"/>
</Residue>
<Residue name="U">
<H name="H1'" parent="C1'"/>
<H name="H2'" parent="C2'"/>
<H name="H3" parent="N3"/>
<H name="H3'" parent="C3'"/>
<H name="H4'" parent="C4'"/>
<H name="H5" parent="C5"/>
<H name="H5'" parent="C5'"/>
<H name="H5''" parent="C5'"/>
<H name="H6" parent="C6"/>
<H name="HO2'" parent="O2'"/>
<H name="HO3'" parent="O3'"/>
<H name="HOP2" parent="OP2"/>
<H name="HOP3" parent="OP3"/>
</Residue>
<Residue name="UNK">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB1" parent="CB"/>
<H name="HB2" parent="CB"/>
<H name="HG1" parent="CG"/>
<H name="HG2" parent="CG"/>
<H name="HG3" parent="CG"/>
</Residue>
<Residue name="VAL">
<H name="H" parent="N"/>
<H name="H2" parent="N" terminal="N"/>
<H name="H3" parent="N" terminal="N"/>
<H name="HA" parent="CA"/>
<H name="HB" parent="CB"/>
<H name="HG11" parent="CG1"/>
<H name="HG12" parent="CG1"/>
<H name="HG13" parent="CG1"/>
<H name="HG21" parent="CG2"/>
<H name="HG22" parent="CG2"/>
<H name="HG23" parent="CG2"/>
</Residue>
</Residues>
wrappers/python/simtk/openmm/app/modeller.py
View file @ dfdf8328
...@@ -5,12 +5,14 @@ __author__ = "Peter Eastman" ...@@ -5,12 +5,14 @@ __author__ = "Peter Eastman"
__version__ = "1.0" __version__ = "1.0"
from simtk.openmm.app import Topology, PDBFile from simtk.openmm.app import Topology, PDBFile
from simtk.openmm.app.forcefield import HAngles
from simtk.openmm.vec3 import Vec3 from simtk.openmm.vec3 import Vec3
from simtk.openmm import NonbondedForce from simtk.openmm import System, Context, NonbondedForce, VerletIntegrator
from simtk.unit import nanometer, molar, elementary_charge, amu, gram, liter, sqrt, is_quantity from simtk.unit import nanometer, molar, elementary_charge, amu, gram, liter, sqrt, is_quantity
import element as elem import element as elem
import os import os
import random import random
import xml.etree.ElementTree as etree
from copy import deepcopy from copy import deepcopy
from math import ceil, floor from math import ceil, floor
...@@ -23,6 +25,8 @@ class Modeller(object): ...@@ -23,6 +25,8 @@ class Modeller(object):
and getPositions() to get the results. and getPositions() to get the results.
""" """
_residueHydrogens = None
def __init__(self, topology, positions): def __init__(self, topology, positions):
"""Create a new Modeller object """Create a new Modeller object
...@@ -353,3 +357,219 @@ class Modeller(object): ...@@ -353,3 +357,219 @@ class Modeller(object):
newTopology.setUnitCellDimensions(deepcopy(box)*nanometer) newTopology.setUnitCellDimensions(deepcopy(box)*nanometer)
self.topology = newTopology self.topology = newTopology
self.positions = newPositions self.positions = newPositions
class _ResidueData:
"""Inner class used to encapsulate data about the hydrogens for a residue."""
def __init__(self, name):
self.name = name
self.variants = []
self.hydrogens = []
class _Hydrogen:
"""Inner class used to encapsulate data about a hydrogen atom."""
def __init__(self, name, parent, maxph, variants, terminal):
self.name = name
self.parent = parent
self.maxph = maxph
self.variants = variants
self.terminal = terminal
def addHydrogens(self, forcefield, pH=7.0, variants=None):
"""Add missing hydrogens to the model.
Some residues can exist in multiple forms depending on the pH and properties of the local environment. These
variants differ in the presence or absence of particular hydrogens. In particular, the following variants
are supported:
Aspartic acid:
ASH: Neutral form with a hydrogen on one of the delta oxygens
ASP: Negatively charged form without a hydrogen on either delta oxygen
Cysteine:
CYS: Neutral form with a hydrogen on the sulfur
CYX: No hydrogen on the sulfur (either negatively charged, or part of a disulfide bond)
Glutamic acid:
GLH: Neutral form with a hydrogen on one of the epsilon oxygens
GLU: Negatively charged form without a hydrogen on either epsilon oxygen
Histidine:
HID: Neutral form with a hydrogen on the ND1 atom
HIE: Neutral form with a hydrogen on the NE2 atom
HIP: Positively charged form with hydrogens on both ND1 and NE2
Lysine:
LYN: Neutral form with two hydrogens on the zeta nitrogen
LYS: Positively charged form with three hydrogens on the zeta nitrogen
The variant to use for each residue is determined by the following rules:
1. The most common variant at the specified pH is selected.
2. Any Cysteine that participates in a disulfide bond uses the CYX variant regardless of pH.
3. For a neutral Histidine residue, the HID or HIE variant is selected based on which one forms a better hydrogen bond.
You can override these rules by explicitly specifying a variant for any residue. Also keep in mind that this
function will only add hydrogens. It will never remove ones that are already present in the model, regardless
of the specified pH.
Parameters:
- forcefield (ForceField) the ForceField to use for determining the positions of hydrogens
- pH (float=7.0) the pH based on which to select variants
- variants (list=None) an optional list of variants to use. If this is specified, its length must equal the number
of residues in the model. variants[i] is the name of the variant to use for residue i (indexed starting at 0).
If an element is None, the standard rules will be followed to select a variant for that residue.
"""
# Check the list of variants.
residues = list(self.topology.residues())
if variants is not None:
if len(variants) != len(residues):
raise ValueError("The length of the variants list must equal the number of residues")
else:
variants = [None]*len(residues)
# Load the residue specifications.
if Modeller._residueHydrogens is None:
Modeller._residueHydrogens = {}
tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', 'hydrogens.xml'))
infinity = float('Inf')
for residue in tree.getroot().findall('Residue'):
resName = residue.attrib['name']
data = Modeller._ResidueData(resName)
Modeller._residueHydrogens[resName] = data
for variant in residue.findall('Variant'):
data.variants.append(variant.attrib['name'])
for hydrogen in residue.findall('H'):
maxph = infinity
if 'maxph' in hydrogen.attrib:
maxph = float(hydrogen.attrib['maxph'])
atomVariants = None
if 'variant' in hydrogen.attrib:
atomVariants = hydrogen.attrib['variant'].split(',')
terminal = None
if 'terminal' in hydrogen.attrib:
terminal = hydrogen.attrib['terminal']
data.hydrogens.append(Modeller._Hydrogen(hydrogen.attrib['name'], hydrogen.attrib['parent'], maxph, atomVariants, terminal))
# Make a list of atoms bonded to each atom.
bonded = {}
for atom1, atom2 in self.topology.bonds():
if atom1 not in bonded:
bonded[atom1] = []
if atom2 not in bonded:
bonded[atom2] = []
bonded[atom1].append(atom2)
bonded[atom2].append(atom1)
# Loop over residues.
newTopology = Topology()
newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
newAtoms = {}
newPositions = []*nanometer
newIndices = []
acceptors = [atom for atom in self.topology.atoms() if atom.element in (elem.oxygen, elem.nitrogen)]
for chain in self.topology.chains():
newChain = newTopology.addChain()
for residue in chain.residues():
newResidue = newTopology.addResidue(residue.name, newChain)
isNTerminal = (residue == chain._residues[0])
isCTerminal = (residue == chain._residues[-1])
if residue.name in Modeller._residueHydrogens:
# Add hydrogens. First select which variant to use.
spec = Modeller._residueHydrogens[residue.name]
variant = variants[residue.index]
if variant is None:
if residue.name == 'CYS':
# If this is part of a disulfide, use CYX.
sulfur = [atom for atom in residue.atoms() if atom.element == elem.sulfur]
if len(sulfur) == 1 and any((atom.residue != residue for atom in bonded[sulfur[0]])):
variant = 'CYX'
if residue.name == 'HIS' and pH > 6.5:
variant = 'HID' # Fix this!!!!!!!!!!!!!!!!!!!
elif residue.name == 'HIS':
variant = 'HIP'
if variant is not None and variant not in spec.variants:
raise ValueError('Illegal variant for %s residue: %s' % (residue.name, variant))
# Make a list of hydrogens that should be present in the residue.
parents = [atom for atom in residue.atoms() if atom.element != elem.hydrogen]
parentNames = [atom.name for atom in parents]
hydrogens = [h for h in spec.hydrogens if (variant is None and pH <= h.maxph) or (h.variants is None and pH <= h.maxph) or (h.variants is not None and variant in h.variants)]
hydrogens = [h for h in hydrogens if h.terminal is None or (isNTerminal and h.terminal == 'N') or (isCTerminal and h.terminal == 'C')]
hydrogens = [h for h in hydrogens if h.parent in parentNames]
# Loop over atoms in the residue, adding them to the new topology along with required hydrogens.
for parent in residue.atoms():
# Add the atom.
newAtom = newTopology.addAtom(parent.name, parent.element, newResidue)
newAtoms[parent] = newAtom
newPositions.append(deepcopy(self.positions[parent.index]))
if parent in parents:
# Match expected hydrogens with existing ones and find which ones need to be added.
existing = [atom for atom in bonded[parent] if atom.element == elem.hydrogen]
expected = [h for h in hydrogens if h.parent == parent.name]
if len(existing) < len(expected):
# Try to match up existing hydrogens to expected ones.
matches = []
for e in existing:
match = [h for h in expected if h.name == e.name]
if len(match) > 0:
matches.append(match[0])
expected.remove(match[0])
else:
matches.append(None)
# If any hydrogens couldn't be matched by name, just match them arbitrarily.
for i in range(len(matches)):
if matches[i] is None:
matches[i] = expected[-1]
expected.remove(expected[-1])
# Add the missing hydrogens.
for h in expected:
newH = newTopology.addAtom(h.name, elem.hydrogen, newResidue)
newIndices.append(newH.index)
if len(expected) == 1:
delta = Vec3(0, 0, 0)*nanometer
for other in bonded[parent]:
delta += self.positions[parent.index]-self.positions[other.index]
else:
delta = Vec3(random.random(), random.random(), random.random())*nanometer
delta *= 0.1*nanometer/sqrt(delta[0]*delta[0]+delta[1]*delta[1]+delta[2]*delta[2])
newPositions.append(self.positions[parent.index]+delta)
newTopology.addBond(newAtom, newH)
else:
# Just copy over the residue.
for atom in residue.atoms():
newAtom = newTopology.addAtom(atom.name, atom.element, newResidue)
newAtoms[atom] = newAtom
newPositions.append(deepcopy(self.positions[atom.index]))
for bond in self.topology.bonds():
if bond[0] in newAtoms and bond[1] in newAtoms:
newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
# The hydrogens were added at random positions. Now use the ForceField to fix them up.
system = forcefield.createSystem(newTopology, constraints=HAngles)
system2 = System()
for i in range(system.getNumParticles()):
system2.addParticle(system.getParticleMass(i))
atoms = list(newTopology.atoms())
for i in range(system.getNumConstraints()):
p1, p2, distance = system.getConstraintParameters(i)
if atoms[p1].element == elem.hydrogen or atoms[p2].element == elem.hydrogen:
system2.addConstraint(p1, p2, distance)
context = Context(system2, VerletIntegrator(0.0))
context.setPositions(newPositions)
context.applyConstraints(0.0001)
self.topology = newTopology
self.positions = context.getState(getPositions=True).getPositions()
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