Source code for parmed.structure

"""
This module contains the core base class for all of the chemical structures with
various topological and force field features.
"""
from __future__ import absolute_import, division

import logging
import math
import os
from collections import defaultdict
from copy import copy

import numpy as np

from . import unit as u
from . import residue
from .constants import DEG_TO_RAD, SMALL
from .exceptions import ParameterError
from .geometry import (STANDARD_BOND_LENGTHS_SQUARED, box_lengths_and_angles_to_vectors,
                       box_vectors_to_lengths_and_angles, distance2)
from .topologyobjects import (AcceptorDonor, Angle, Atom, AtomList, Bond, ChiralFrame, Cmap,
                              Dihedral, DihedralType, DihedralTypeList, ExtraPoint, Group, Improper,
                              MultipoleFrame, NonbondedException, NoUreyBradley, OutOfPlaneBend,
                              OutOfPlaneExtraPointFrame, PiTorsion, ResidueList, StretchBend,
                              ThreeParticleExtraPointFrame, TorsionTorsion, TrackedList,
                              TrigonalAngle, TwoParticleExtraPointFrame, UnassignedAtomType,
                              UreyBradley, Link)
from .utils import PYPY, find_atom_pairs, tag_molecules
from .utils.decorators import needs_openmm
from .utils.six import integer_types, iteritems, string_types
from .utils.six.moves import range, zip
from .vec3 import Vec3

# Try to import the OpenMM modules
try:
    from simtk.openmm import app
    from simtk import openmm as mm
    from simtk.openmm.app.internal.unitcell import reducePeriodicBoxVectors
except ImportError:  # pragma: no cover
    app = mm = None  # pragma: no cover

LOGGER = logging.getLogger(__name__)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Private attributes and methods

def _strip_box_units(args):
    new_args = []
    for arg in args:
        # Handle 3 types of arguments here: units, regular numbers, and
        # iterables. Iterables have units removed recursively. But protect
        # against strings, since you get infinite recursion that way
        if u.is_quantity(arg):
            if arg.unit.is_compatible(u.degree):
                new_args.append(arg.value_in_unit(u.degree))
            else:
                new_args.append(arg.value_in_unit(u.angstroms))
        elif isinstance(arg, string_types):
            raise TypeError('Unit cell cannot have strings')
        else:
            try:
                iter(arg)
                arg = _strip_box_units(arg)
            except TypeError:
                pass
            new_args.append(arg)
    return new_args

def _bondi(atom):
    if atom.atomic_number == 6: return 1.7
    if atom.atomic_number == 1: return 1.2
    if atom.atomic_number == 7: return 1.55
    if atom.atomic_number == 14: return 2.1
    if atom.atomic_number == 15: return 1.85
    if atom.atomic_number == 16: return 1.8
    return 1.5

def _mbondi(atom):
    if atom.atomic_number == 1:
        bondeds = atom.bond_partners
        if bondeds[0].atomic_number in (6, 7):
            return 1.3
        if bondeds[0].atomic_number in (8, 16):
            return 0.8
        return 1.2
    return _bondi(atom)

def _mbondi2(atom):
    if atom.atomic_number == 1:
        if atom.bond_partners[0].atomic_number == 7:
            return 1.3
        return 1.2
    return _bondi(atom)

def _mbondi3(atom):
    if atom.residue.name in ('GLU', 'ASP', 'GL4', 'AS4'):
        if atom.name.startswith('OE') or atom.name.startswith('OD'):
            return 1.4
    elif atom.residue.name == 'ARG':
        if atom.name.startswith('HH') or atom.name.startswith('HE'):
            return 1.17
    if atom.name == 'OXT':
        return 1.4
    return _mbondi2(atom)

@needs_openmm
def _gb_rad_screen(atom, model):
    """
    Gets the default GB parameters for a given atom according to a specific
    Generalized Born model

    Parameters
    ----------
    atom : :class:`Atom`
        The atom to get the default GB parameters for
    model : ``app.HCT, app.OBC1, or app.OBC2``
        The GB model to get the default parameters for (app.GBn and app.GBn2 are
        already handled in Structure._get_gb_parameters)

    Returns
    -------
    radius, screen [,alpha, beta, gamma] : ``float, float [,float, float, float]``
        The intrinsic radius of the atom and the screening factor of the atom.
        If the model is GBn2, alpha, beta, and gamma parameters are also
        returned
    """
    if model in (app.OBC1, app.OBC2):
        rad = _mbondi2(atom)
    else:
        rad = _mbondi(atom)
    if atom.atomic_number == 1: return rad, 0.85
    if atom.atomic_number == 6: return rad, 0.72
    if atom.atomic_number == 7: return rad, 0.79
    if atom.atomic_number == 8: return rad, 0.85
    if atom.atomic_number == 9: return rad, 0.88
    if atom.atomic_number == 15: return rad, 0.86
    if atom.atomic_number == 16: return rad, 0.96
    return rad, 0.8

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

[docs]class Structure(object): """ A chemical structure composed of atoms, bonds, angles, torsions, and other topological features Attributes ---------- atoms : :class:`AtomList` List of all atoms in the structure residues : :class:`ResidueList` List of all residues in the structure bonds : :class:`TrackedList` (:class:`Bond`) List of all bonds in the structure angles : :class:`TrackedList` (:class:`Angle`) List of all angles in the structure dihedrals : :class:`TrackedList` (:class:`Dihedral`) List of all dihedrals in the structure rb_torsions : :class:`TrackedList` (:class:`Dihedral`) List of all Ryckaert-Bellemans torsions in the structure urey_bradleys : :class:`TrackedList` (:class:`UreyBradley`) List of all Urey-Bradley angle bends in the structure impropers : :class:`TrackedList` (:class:`Improper`) List of all CHARMM-style improper torsions in the structure cmaps : :class:`TrackedList` (:class:`Cmap`) List of all CMAP objects in the structure trigonal_angles : :class:`TrackedList` (:class:`TrigonalAngle`) List of all AMOEBA-style trigonal angles in the structure out_of_plane_bends : :class:`TrackedList` (:class:`OutOfPlaneBends`) List of all AMOEBA-style out-of-plane bending angles pi_torsions : :class:`TrackedList` (:class:`PiTorsion`) List of all AMOEBA-style pi-torsion angles stretch_bends : :class:`TrackedList` (:class:`StretchBend`) List of all AMOEBA-style stretch-bend compound bond/angle terms torsion_torsions : :class:`TrackedList` (:class:`TorsionTorsion`) List of all AMOEBA-style coupled torsion-torsion terms chiral_frames : :class:`TrackedList` (:class:`ChiralFrame`) List of all AMOEBA-style chiral frames defined in the structure multipole_frames : :class:`TrackedList` (:class:`MultipoleFrame`) List of all AMOEBA-style multipole frames defined in the structure adjusts : :class:`TrackedList` (:class:`NonbondedException`) List of all nonbonded pair-exception rules acceptors : :class:`TrackedList` (:class:`AcceptorDonor`) List of all H-bond acceptors, if that information is present donors : :class:`TrackedList` (:class:`AcceptorDonor`) List of all H-bond donors, if that information is present groups : :class:`TrackedList` (:class:`Group`) List of all CHARMM-style GROUP objects (whatever those are used for) box : ``list of 6 floats`` Box dimensions (a, b, c, alpha, beta, gamma) for the unit cell. If no box is defined, `box` is set to `None` space_group : ``str`` The space group of the structure (default is "P 1") nrexcl : ``int`` The number of bonds away that an atom can be in order to be excluded from the direct nonbonded summation title : ``str`` Cosmetic only, it is an arbitrary title assigned to the system. Default value is an empty string positions : u.Quantity(list(Vec3), u.angstroms) Unit-bearing atomic coordinates. If not all atoms have coordinates, this property is None coordinates : np.ndarray of shape (nframes, natom, 3) If no coordinates are set, this is set to None. The first frame will match the coordinates present on the atoms. symmetry : :class:`Symmetry` if no symmetry is set, this is set to None. links : :class:`TrackedList` (:class:`Link`) The list of Link definitions for this Structure Notes ----- This class also has a handful of type lists for each of the attributes above (excluding `atoms`, `residues`, `chiral_frames`, and `multipole_frames`). They are all TrackedList instances that are designed to hold the relevant parameter type. The list is: bond_types, angle_types, dihedral_types, urey_bradley_types, improper_types, cmap_types, trigonal_angle_types, out_of_plane_bend_types, pi_torsion_types, stretch_bend_types, torsion_torsion_types, adjust_types dihedral_types _may_ be a list of :class:`DihedralType` instances, since torsion profiles are often represented by a Fourier series with multiple terms """ # Force groups assigned to each type of force BOND_FORCE_GROUP = 0 ANGLE_FORCE_GROUP = 1 DIHEDRAL_FORCE_GROUP = 2 UREY_BRADLEY_FORCE_GROUP = 3 IMPROPER_FORCE_GROUP = 4 CMAP_FORCE_GROUP = 5 TRIGONAL_ANGLE_FORCE_GROUP = 6 OUT_OF_PLANE_BEND_FORCE_GROUP = 7 PI_TORSION_FORCE_GROUP = 8 STRETCH_BEND_FORCE_GROUP = 9 TORSION_TORSION_FORCE_GROUP = 10 NONBONDED_FORCE_GROUP = 11 RB_TORSION_FORCE_GROUP = 12 #===================================================
[docs] def __init__(self): # Topological object lists self.atoms = AtomList() self.residues = ResidueList() self.bonds = TrackedList() self.angles = TrackedList() self.dihedrals = TrackedList() self.rb_torsions = TrackedList() self.urey_bradleys = TrackedList() self.impropers = TrackedList() self.cmaps = TrackedList() self.trigonal_angles = TrackedList() self.out_of_plane_bends = TrackedList() self.pi_torsions = TrackedList() self.stretch_bends = TrackedList() self.torsion_torsions = TrackedList() self.chiral_frames = TrackedList() self.multipole_frames = TrackedList() self.adjusts = TrackedList() # Extraneous information stored in CHARMM PSF files... not used as far # as I can tell for much of anything self.acceptors = TrackedList() self.donors = TrackedList() self.groups = TrackedList() # Parameter type lists self.bond_types = TrackedList() self.angle_types = TrackedList() self.dihedral_types = TrackedList() self.urey_bradley_types = TrackedList() self.improper_types = TrackedList() self.rb_torsion_types = TrackedList() self.cmap_types = TrackedList() self.trigonal_angle_types = TrackedList() self.out_of_plane_bend_types = TrackedList() self.pi_torsion_types = TrackedList() self.stretch_bend_types = TrackedList() self.torsion_torsion_types = TrackedList() self.adjust_types = TrackedList() self.links = TrackedList() self._box = None self._coordinates = None self.space_group = "P 1" self.unknown_functional = False self.nrexcl = 3 self.title = '' self._combining_rule = 'lorentz' self.symmetry = None
#=================================================== def __repr__(self): natom = len(self.atoms) nres = len(self.residues) nextra = sum([isinstance(a, ExtraPoint) for a in self.atoms]) retstr = ['<%s %d atoms' % (type(self).__name__, natom)] if nextra > 0: retstr.append(' [%d EPs]' % nextra) retstr.append('; %d residues' % nres) nbond = len(self.bonds) retstr.append('; %d bonds' % nbond) if self.box is not None: retstr.append('; PBC') if (abs(self.box[3]-90) > SMALL or abs(self.box[4]-90) > SMALL or abs(self.box[5]-90) > SMALL): retstr.append(' (triclinic)') else: retstr.append(' (orthogonal)') # Just assume that if the first bond has a defined type, so does # everything else... we don't want __repr__ to be super expensive if len(self.bonds) > 0 and self.bonds[0].type is not None: retstr.append('; parametrized>') else: retstr.append('; NOT parametrized>') return ''.join(retstr) #===================================================
[docs] def add_atom(self, atom, resname, resnum, chain='', inscode='', segid=''): """ Adds a new atom to the Structure, adding a new residue to `residues` if it has a different name or number as the last residue added and adding it to the `atoms` list. Parameters ---------- atom : :class:`Atom` The atom to add to this residue list resname : ``str`` The name of the residue this atom belongs to resnum : ``int`` The number of the residue this atom belongs to chain : ``str`` The chain ID character for this residue inscode : ``str`` The insertion code ID character for this residue (it is stripped) segid : ``str`` The segment identifier for this residue (it is stripped) Notes ----- If the residue name and number differ from the last residue in this list, a new residue is added and the atom is added to that residue """ self.residues.add_atom(atom, resname, resnum, chain, inscode, segid) self.atoms.append(atom)
#===================================================
[docs] def add_atom_to_residue(self, atom, residue): """ Adds a new atom to the Structure at the end if the given residue Parameters ---------- atom : :class:`Atom` The atom to add to the system residue : :class:`Residue` The residue to which to add this atom. It MUST be part of this Structure instance already or a ValueError is raised Notes ----- This atom is added at the end of the residue and is inserted into the `atoms` list in such a way that all residues are composed of atoms contiguous in the atoms list. For large systems, this may be a relatively expensive operation """ # Make sure residue belongs to this list if residue.list is not self.residues: raise ValueError('Residue is not part of the structure') last_atom = residue.atoms[-1] residue.add_atom(atom) # Special-case if this is going to be the last atom if not self.atoms or last_atom is self.atoms[-1]: self.atoms.append(atom) else: self.atoms.insert(last_atom.idx + 1, atom)
#=================================================== def __copy__(self): """ A deep copy of the Structure """ return self.copy(type(self)) #===================================================
[docs] def copy(self, cls, split_dihedrals=False): """ Makes a copy of the current structure as an instance of a specified subclass Parameters ---------- cls : Structure subclass The returned object is a copy of this structure as a `cls` instance split_dihedrals : ``bool`` If True, then the Dihedral entries will be split up so that each one is paired with a single DihedralType (rather than a DihedralTypeList) Returns ------- *cls* instance The instance of the Structure subclass `cls` with a copy of the current Structure's topology information """ c = cls() for atom in self.atoms: res = atom.residue a = copy(atom) c.add_atom(a, res.name, res.number, res.chain, res.insertion_code, res.segid) # Now copy all of the types for bt in self.bond_types: c.bond_types.append(copy(bt)) c.bond_types.claim() for at in self.angle_types: c.angle_types.append(copy(at)) c.angle_types.claim() if split_dihedrals: ndt = 0 mapdt = {} for idt, dt in enumerate(self.dihedral_types): if hasattr(dt, '__iter__'): for t in dt: c.dihedral_types.append(copy(t)) mapdt.setdefault(idt, []).append(ndt) ndt += 1 else: mapdt.setdefault(idt, []).append(ndt) ndt += 1 c.dihedral_types.append(copy(dt)) else: for dt in self.dihedral_types: c.dihedral_types.append(copy(dt)) c.dihedral_types.claim() for ut in self.urey_bradley_types: c.urey_bradley_types.append(copy(ut)) c.urey_bradley_types.claim() for it in self.improper_types: c.improper_types.append(copy(it)) c.improper_types.claim() for rt in self.rb_torsion_types: c.rb_torsion_types.append(copy(rt)) c.rb_torsion_types.claim() for ct in self.cmap_types: c.cmap_types.append(copy(ct)) c.cmap_types.claim() for ta in self.trigonal_angle_types: c.trigonal_angle_types.append(copy(ta)) c.trigonal_angle_types.claim() for ot in self.out_of_plane_bend_types: c.out_of_plane_bend_types.append(copy(ot)) c.out_of_plane_bend_types.claim() for pt in self.pi_torsion_types: c.pi_torsion_types.append(copy(pt)) c.pi_torsion_types.claim() for st in self.stretch_bend_types: c.stretch_bend_types.append(copy(st)) c.stretch_bend_types.claim() for tt in self.torsion_torsion_types: c.torsion_torsion_types.append(copy(tt)) c.torsion_torsion_types.claim() for at in self.adjust_types: c.adjust_types.append(copy(at)) c.adjust_types.claim() # Now create the topological objects atoms = c.atoms for b in self.bonds: c.bonds.append( Bond(atoms[b.atom1.idx], atoms[b.atom2.idx], b.type and c.bond_types[b.type.idx]) ) c.bonds[-1].funct = b.funct for a in self.angles: c.angles.append( Angle(atoms[a.atom1.idx], atoms[a.atom2.idx], atoms[a.atom3.idx], a.type and c.angle_types[a.type.idx]) ) c.angles[-1].funct = a.funct if split_dihedrals: for d in self.dihedrals: if hasattr(d.type, '__iter__'): for i in range(len(d.type)): ie = d.ignore_end or i < len(d.type) - 1 ti = mapdt[d.type.idx][i] c.dihedrals.append( Dihedral(atoms[d.atom1.idx], atoms[d.atom2.idx], atoms[d.atom3.idx], atoms[d.atom4.idx], improper=d.improper, ignore_end=ie, type=c.dihedral_types[ti]) ) c.dihedrals[-1]._funct = d._funct else: ti = mapdt[d.type.idx][0] c.dihedrals.append( Dihedral(atoms[d.atom1.idx], atoms[d.atom2.idx], atoms[d.atom3.idx], atoms[d.atom4.idx], improper=d.improper, ignore_end=d.ignore_end, type=d.type and c.dihedral_types[ti]) ) c.dihedrals[-1]._funct = d._funct else: for d in self.dihedrals: if d.type is None: typ = None else: typ = c.dihedral_types[d.type.idx] c.dihedrals.append( Dihedral(atoms[d.atom1.idx], atoms[d.atom2.idx], atoms[d.atom3.idx], atoms[d.atom4.idx], improper=d.improper, ignore_end=d.ignore_end, type=typ) ) c.dihedrals[-1]._funct = d._funct for ub in self.urey_bradleys: if ub.type is NoUreyBradley: typ = NoUreyBradley else: typ = ub.type and c.urey_bradley_types[ub.type.idx] c.urey_bradleys.append( UreyBradley(atoms[ub.atom1.idx], atoms[ub.atom2.idx], typ) ) for i in self.impropers: c.impropers.append( Improper(atoms[i.atom1.idx], atoms[i.atom2.idx], atoms[i.atom3.idx], atoms[i.atom4.idx], i.type and c.improper_types[i.type.idx]) ) for r in self.rb_torsions: c.rb_torsions.append( Dihedral(atoms[r.atom1.idx], atoms[r.atom2.idx], atoms[r.atom3.idx], atoms[r.atom4.idx], type=r.type and c.rb_torsion_types[r.type.idx]) ) c.rb_torsions[-1]._funct = r._funct for cm in self.cmaps: c.cmaps.append( Cmap(atoms[cm.atom1.idx], atoms[cm.atom2.idx], atoms[cm.atom3.idx], atoms[cm.atom4.idx], atoms[cm.atom5.idx], cm.type and c.cmap_types[cm.type.idx]) ) c.cmaps[-1].funct = cm.funct for t in self.trigonal_angles: c.trigonal_angles.append( TrigonalAngle(atoms[t.atom1.idx], atoms[t.atom2.idx], atoms[t.atom3.idx], atoms[t.atom4.idx], t.type and c.trigonal_angle_types[t.type.idx]) ) for o in self.out_of_plane_bends: c.out_of_plane_bends.append( OutOfPlaneBend(atoms[o.atom1.idx], atoms[o.atom2.idx], atoms[o.atom3.idx], atoms[o.atom4.idx], o.type and c.out_of_plane_bend_types[o.type.idx]) ) for p in self.pi_torsions: c.pi_torsions.append( PiTorsion(atoms[p.atom1.idx], atoms[p.atom2.idx], atoms[p.atom3.idx], atoms[p.atom4.idx], atoms[p.atom5.idx], atoms[p.atom6.idx], p.type and c.pi_torsion_types[p.type.idx]) ) for s in self.stretch_bends: c.stretch_bends.append( StretchBend(atoms[s.atom1.idx], atoms[s.atom2.idx], atoms[s.atom3.idx], s.type and c.stretch_bend_types[s.type.idx]) ) for t in self.torsion_torsions: c.torsion_torsions.append( TorsionTorsion(atoms[t.atom1.idx], atoms[t.atom2.idx], atoms[t.atom3.idx], atoms[t.atom4.idx], atoms[t.atom5.idx], t.type and c.torsion_torsion_types[t.type.idx]) ) for ch in self.chiral_frames: c.chiral_frames.append( ChiralFrame(atoms[ch.atom1.idx], atoms[ch.atom2.idx], ch.chirality) ) for m in self.multipole_frames: c.multipole_frames.append( MultipoleFrame(atoms[m.atom.idx], m.frame_pt_num, m.vectail, m.vechead, m.nvec) ) for a in self.adjusts: c.adjusts.append( NonbondedException(atoms[a.atom1.idx], atoms[a.atom2.idx], a.type and c.adjust_types[a.type.idx]) ) for a in self.acceptors: c.acceptors.append( AcceptorDonor(atoms[a.atom1.idx], atoms[a.atom2.idx]) ) for d in self.donors: c.donors.append( AcceptorDonor(atoms[d.atom1.idx], atoms[d.atom2.idx]) ) for g in self.groups: c.groups.append(Group(atoms[g.atom.idx], g.type, g.move)) for l in self.links: c.links.append( Link(atoms[l.atom1.idx], atoms[l.atom2.idx], l.length, l.symmetry_op1, l.symmetry_op2) ) c._box = copy(self._box) c._coordinates = copy(self._coordinates) c.combining_rule = self.combining_rule # Transfer TER cards for r1, r2 in zip(c.residues, self.residues): r1.ter = r2.ter return c
#===================================================
[docs] def to_dataframe(self): """ Generates a DataFrame from the current Structure's atomic properties Returns ------- df : DataFrame DataFrame with all atomic properties See Also -------- :func:`parmed.utils.pandautils.create_dataframe` """ from parmed.utils.pandautils import create_dataframe return create_dataframe(self)
#===================================================
[docs] def load_dataframe(self, df): """ Loads atomic properties from an input DataFrame Parameters ---------- df : pandas.DataFrame A pandas DataFrame with atomic properties that will be used to set the properties on the current list of atoms See Also -------- :func:`parmed.utils.pandautils.load_dataframe` """ from parmed.utils.pandautils import load_dataframe return load_dataframe(self, df)
#===================================================
[docs] def is_changed(self): """ Determines if any of the topology has changed for this structure """ return (self.atoms.changed or self.residues.changed or self.bonds.changed or self.trigonal_angles.changed or self.dihedrals.changed or self.urey_bradleys.changed or self.impropers.changed or self.cmaps.changed or self.angles.changed or self.out_of_plane_bends.changed or self.pi_torsions.changed or self.stretch_bends.changed or self.torsion_torsions.changed or self.chiral_frames.changed or self.multipole_frames.changed or self.adjusts.changed or self.acceptors.changed or self.donors.changed or self.groups.changed or self.bond_types.changed or self.angle_types.changed or self.dihedral_types.changed or self.urey_bradley_types.changed or self.cmap_types.changed or self.improper_types.changed or self.adjust_types.changed or self.trigonal_angle_types.changed or self.out_of_plane_bends.changed or self.stretch_bend_types.changed or self.torsion_torsion_types.changed or self.pi_torsion_types.changed or self.rb_torsions.changed or self.rb_torsion_types.changed)
#===================================================
[docs] def unchange(self): """ Toggles all lists so that they do not indicate any changes """ self.atoms.changed = False self.residues.changed = False self.bonds.changed = False self.angles.changed = False self.dihedrals.changed = False self.urey_bradleys.changed = False self.impropers.changed = False self.cmaps.changed = False self.trigonal_angles.changed = False self.out_of_plane_bends.changed = False self.pi_torsions.changed = False self.stretch_bends.changed = False self.torsion_torsions.changed = False self.chiral_frames.changed = False self.multipole_frames.changed = False self.adjusts.changed = False self.acceptors.changed = False self.donors.changed = False self.groups.changed = False # Parameter type lists self.bond_types.changed = False self.angle_types.changed = False self.dihedral_types.changed = False self.urey_bradley_types.changed = False self.improper_types.changed = False self.cmap_types.changed = False self.trigonal_angle_types.changed = False self.out_of_plane_bend_types.changed = False self.pi_torsion_types.changed = False self.stretch_bend_types.changed = False self.torsion_torsion_types.changed = False self.adjust_types.changed = False
#===================================================
[docs] def prune_empty_terms(self): """ Looks through all of the topological lists and gets rid of terms in which at least one of the atoms is None or has an `idx` attribute set to -1 (indicating that it has been removed from the `atoms` atom list) """ self._prune_empty_bonds() self._prune_empty_angles() self._prune_empty_dihedrals() self._prune_empty_rb_torsions() self._prune_empty_ureys() self._prune_empty_impropers() self._prune_empty_cmaps() self._prune_empty_trigonal_angles() self._prune_empty_out_of_plane_bends() self._prune_empty_pi_torsions() self._prune_empty_stretch_bends() self._prune_empty_torsion_torsions() self._prune_empty_chiral_frames() self._prune_empty_multipole_frames() self._prune_empty_adjusts()
#===================================================
[docs] def update_dihedral_exclusions(self): """ Nonbonded exclusions and exceptions have the following priority: bond -> angle -> dihedral Since bonds and angles are completely excluded, any ring systems in which two atoms are attached by a bond or angle as well as a dihedral should be completely excluded as the bond and angle exclusion rules take precedence. If a Bond or Angle was _added_ to the structure between a pair of atoms previously connected only by a dihedral term, it's possible that those two atoms have both an exclusion *and* an exception defined. The result of this scenario is that sander and pmemd will happily compute an energy, _including_ the 1-4 nonbonded terms between atoms now connected by a bond or an Angle. OpenMM, on the other hand, will complain about an exception specified multiple times. This method scans through all of the dihedrals in which `ignore_end` is `False` and turns it to `True` if the two end atoms are in the bond or angle partners arrays """ set14 = set() deferred_dihedrals = [] # to work around pmemd tossing pn=0 dihedrals for dihedral in self.dihedrals: if dihedral.ignore_end : continue if (dihedral.atom1 in dihedral.atom4.bond_partners or dihedral.atom1 in dihedral.atom4.angle_partners): dihedral.ignore_end = True elif (dihedral.atom1.idx, dihedral.atom4.idx) in set14: # Avoid double counting of 1-4 in a six-membered ring dihedral.ignore_end = True elif isinstance(dihedral.type, DihedralType) and dihedral.type.per == 0: # This needs to be done to work around a "feature" in pmemd where periodicity=0 # torsions are thrown away. Since AmberParm never has any DihedralTypeList, we only # need to defer periodicity=0 torsions for DihedralType torsions. deferred_dihedrals.append(dihedral) else: set14.add((dihedral.atom1.idx, dihedral.atom4.idx)) set14.add((dihedral.atom4.idx, dihedral.atom1.idx)) # Only keep ignore_end = False if we *must*; i.e., if the exclusion it would have # added was not added by anybody else for dihedral in deferred_dihedrals: if (dihedral.atom1.idx, dihedral.atom4.idx) in set14: dihedral.ignore_end = True
#===================================================
[docs] def strip(self, selection): """ Deletes a subset of the atoms corresponding to an atom-based selection. Parameters ---------- selection : :class:`AmberMask`, ``str``, or ``iterable`` This is the selection of atoms that will be deleted from this structure. If it is a string, it will be interpreted as an AmberMask. If it is an AmberMask, it will be converted to a selection of atoms. If it is an iterable, it must be the same length as the `atoms` list. """ from parmed.amber import AmberMask if isinstance(selection, AmberMask): if selection.parm is not self: raise TypeError('passed mask does not belong to Structure') sel = selection.Selection() elif isinstance(selection, string_types): sel = AmberMask(self, selection).Selection() else: try: sel = list(selection) except TypeError: raise TypeError('Selection not a supported type [%s]' % type(selection).__name__) if len(sel) != len(self.atoms): raise ValueError('Selection iterable wrong length') atomlist = sorted([i for i, s in enumerate(sel) if s]) for i in reversed(atomlist): del self.atoms[i] self.prune_empty_terms() self.residues.prune() self.unchange() # Slice out coordinates if present if self._coordinates is not None: if PYPY: # numpypy does not currently support advanced indexing it seems self._coordinates = np.array( # pragma: no cover [[[x, y, z] for i, (x, y, z) in enumerate(crd) if sel[i]==0] for crd in self._coordinates] ) else: # Pure numpy is faster in CPython, so do that when we can self._coordinates = self._coordinates[:, np.array(sel)==0]
#===================================================
[docs] def assign_bonds(self, *reslibs): """ Assigns bonds to all atoms based on the provided residue template libraries. Atoms whose names are *not* in the templates, as well as those residues for whom no template is found, is assigned to bonds based on distances. Parameters ---------- reslibs : dict{str: ResidueTemplate} Any number of residue template libraries. By default, assign_bonds knows about the standard amino acid, RNA, and DNA residues. """ # Import here to avoid circular references from parmed.modeller import StandardBiomolecularResidues # Build a composite dict of all residue templates all_residues = copy(StandardBiomolecularResidues) for lib in reslibs: all_residues.update(lib) # Walk through every residue and assign bonds from the templates unassigned_residues = set() unassigned_atoms = set() cysteine_sg = set() for res in self.residues: templ = _res_in_templlib(res, all_residues) if templ is None: # Don't have a template for this residue. Keep a note and go on unassigned_atoms.update(res.atoms) unassigned_residues.add(res) continue resatoms = {a.name: a for a in res.atoms} for a in res.atoms: if a.name not in templ.map: unassigned_atoms.add(a) continue # We matched our template atom. Transfer the element # information, as it is more reliable a.atomic_number = templ.map[a.name].atomic_number for bp in templ.map[a.name].bond_partners: if (bp.name in resatoms and resatoms[bp.name] not in a.bond_partners): if a not in resatoms[bp.name].bond_partners: self.bonds.append(Bond(a, resatoms[bp.name])) # Now go through each residue and assign heads and tails. This walks # through the residues and bonds residue i's tail with residue j's head. # So there is nothing to do for the last residue. Omit it from the loop # so self.residues[i+1] never raises an IndexError for i, res in enumerate(self.residues[:-1]): templ = _res_in_templlib(res, all_residues) # TER cards and changing chains prevents bonding to the next residue if res.ter or (res.chain and (res.chain != self.residues[i+1].chain)): continue ntempl = _res_in_templlib(self.residues[i+1], all_residues) if templ is None and ntempl is None: # Any cross-link here picked up later continue if templ is None: if ntempl.head is None: continue # Next residue doesn't bond to the previous one # See if any atom in templ is close enough to bond to the head # atom of the next residue's template for head in self.residues[i+1].atoms: if head.name == ntempl.head.name: break else: LOGGER.warning('Could not find the head atom of the next template! Bond ' 'pattern may be wrong, which could lead to extra TER cards in a ' 'PDB file') continue # head atom not found! for a in res.atoms: maxdist = STANDARD_BOND_LENGTHS_SQUARED[(a.atomic_number, head.atomic_number)] if distance2(a, head) < maxdist: if a not in head.bond_partners: self.bonds.append(Bond(a, head)) break continue if templ.tail is None: continue # This residue does not bond with the next one if ntempl is None: # See if any atom in the next residue is bonding distance away # from my tail for tail in res.atoms: if tail.name == templ.tail.name: break else: continue # tail not found for a in self.residues[i+1].atoms: maxdist = STANDARD_BOND_LENGTHS_SQUARED[(a.atomic_number, tail.atomic_number)] if distance2(a, tail) < maxdist: if a not in tail.bond_partners: self.bonds.append(Bond(a, tail)) break continue if ntempl.head is None: continue # Next residue does not bond with this one # We have templates for both atoms, and both have a head and a tail for tail in res.atoms: if tail.name == templ.tail.name: break else: continue # head could not be found for head in self.residues[i+1].atoms: if head.name == ntempl.head.name: break else: continue # tail could not be found if head not in tail.bond_partners: self.bonds.append(Bond(head, tail)) # Now time to find bonds based on distance. First thing to do is to find # all CYS residues and pull out those that have an SG atom with only 1 # bond partner, since those are *very* commonly bonded by # post-translational modifications. for res in self.residues: if len(res.name) != 3 or not residue.AminoAcidResidue.has(res.name): continue if residue.AminoAcidResidue.get(res.name).abbr != 'CYS': continue # Cysteine! Find the SG and add it to unassigned atoms if it only # has 1 atom bonded to it for a in res.atoms: if a.name == 'SG' and len(a.bond_partners) < 2: unassigned_atoms.add(a) cysteine_sg.add(a) break # Can't do distance-based bond determination without coordinates if self.coordinates is None: return # OK, now go through all atoms that have not been assigned. If an atom # has not been assigned, but its residue HAS been assigned, then that # means the name is wrong. A likely culprit is bad nomenclature for # atoms. So in this case (and this case only), look for bond partners in # an 'assigned' residue (but only the same residue we are part of). One # thing this misses is when the head or tail atom is misnamed. mindist = math.sqrt(max(STANDARD_BOND_LENGTHS_SQUARED.values())) pairs = find_atom_pairs(self, mindist, unassigned_atoms) for atom in unassigned_atoms: for partner in pairs[atom.idx]: maxdist = STANDARD_BOND_LENGTHS_SQUARED[(atom.atomic_number, partner.atomic_number)] if (distance2(atom, partner) < maxdist and atom not in partner.bond_partners): self.bonds.append(Bond(atom, partner)) # Now look through all atoms in this template if it's a template # that's already been assigned. If it's already in an unassigned # residue, then all that residue's atoms were in unassigned_atoms, # so there's no need to look through the residue again. Also, we # added cysteine SG atoms to the unassigned atoms list to see if it # is involved in post-translational bonds. But there *was* a # template match for cysteine and SG was paired correctly. So don't # try to assign it to bonds with other atoms in the cysteine here. if atom.residue in unassigned_residues or atom in cysteine_sg: continue for partner in atom.residue.atoms: if partner is atom: continue maxdist = STANDARD_BOND_LENGTHS_SQUARED[(atom.atomic_number, partner.atomic_number)] if (distance2(atom, partner) < maxdist and atom not in partner.bond_partners): self.bonds.append(Bond(atom, partner))
# All reasonable bonds have now been added #===================================================
[docs] def visualize(self, *args, **kwargs): """Use nglview for visualization. This only works with Jupyter notebook and require to install `nglview` Examples -------- >>> import parmed as pmd >>> parm = pmd.download_PDB('1tsu') >>> parm.visualize() Parameters ---------- args and kwargs : positional and keyword arguments given to nglview, optional """ if self.coordinates is None: raise ValueError('coordinates must not be None') from nglview import show_parmed return show_parmed(self, *args, **kwargs)
#=================================================== def __getitem__(self, selection): """ Allows extracting a single atom from the structure or a slice of atoms as a new Structure instance. The following syntaxes are allowed: - struct[str] : str is interpreted as an Amber selection mask - struct[[sel,[sel,]],sel] : sel can be a list of indices (or strings for chain selection) an integer, or a slice Parameters ---------- selection : str, or slice|iter|str, slice|iter|int, slice|iter|int Which atoms to select for the new structure Returns ------- struct : :class:`Structure` If more than one atom is selected, the resulting return value is a new structure containing all selected atoms and all valence terms (and types) in which all involved atoms are present or atom : :class:`Atom` If the selection is a single integer, a tuple of two integers, or a tuple of a string and two integers. This is equivalent to selecting a particular atom, a particular atom from a particular residue, or a particular atom from a particular residue from a particular chain, respectively. All other selections return a Structure instance, even if that selection happens to only select a single atom. Notes ----- When selecting more than 1 atom, this is a costly operation. It is currently implemented by making a full copy of the object and then stripping the unused atoms. Raises ------ ``TypeError`` : if selection (or any part of the selection) is not an allowed type ``ValueError`` : if the selection is a boolean-like list and its length is not the same as the number of atoms in the system """ if isinstance(selection, integer_types): return self.atoms[selection] selection = self._get_selection_array(selection) # _get_selection_array might have returned either an Atom or None if # there is no selection. Handle those and return, or continue if we got # a list of atoms if selection is None: return type(self)() elif isinstance(selection, Atom): return selection sumsel = sum(selection) if sumsel == 0: # No atoms selected. Return empty type return type(self)() # The cumulative sum of selection will give our index + 1 of each # selected atom into the new structure scan = [selection[0]] for i in range(1, len(selection)): scan.append(scan[i-1] + selection[i]) # Zero-out the unselected atoms scan = [x * y for x, y in zip(scan, selection)] # Copy all parameters struct = type(self)() for i, atom in enumerate(self.atoms): if not selection[i]: continue res = atom.residue if res.number == 0: num = res.idx else: num = res.number struct.add_atom(copy(atom), res.name, num, res.chain, res.insertion_code, res.segid) def copy_valence_terms(oval, otyp, sval, styp, attrlist): """ Copies the valence terms from one list to another; oval=Other VALence; otyp=Other TYPe; sval=Self VALence; styp=Self TYPe; attrlist=ATTRibute LIST (atom1, atom2, ...) """ otypcp = [copy(typ) for typ in styp] used_types = [False for typ in otypcp] for val in sval: ats = [getattr(val, attr) for attr in attrlist] # Make sure all of our atoms in this valence term is "selected" indices = [scan[at.idx] for at in ats if isinstance(at, Atom)] if not all(indices): continue # Add the type if applicable kws = dict() if hasattr(val, 'type') and val.type is NoUreyBradley: kws['type'] = NoUreyBradley # special-case singleton elif otypcp and val.type is not None: kws['type'] = otypcp[val.type.idx] used_types[val.type.idx] = True for i, at in enumerate(ats): if isinstance(at, Atom): ats[i] = struct.atoms[scan[at.idx]-1] oval.append(type(val)(*ats, **kws)) if hasattr(val, 'funct'): oval[-1].funct = val.funct # Now tack on the "new" types copied from `other` for used, typ in zip(used_types, otypcp): if used: otyp.append(typ) if hasattr(otyp, 'claim'): otyp.claim() copy_valence_terms(struct.bonds, struct.bond_types, self.bonds, self.bond_types, ['atom1', 'atom2']) copy_valence_terms(struct.angles, struct.angle_types, self.angles, self.angle_types, ['atom1', 'atom2', 'atom3']) copy_valence_terms(struct.dihedrals, struct.dihedral_types, self.dihedrals, self.dihedral_types, ['atom1', 'atom2', 'atom3', 'atom4', 'improper', 'ignore_end']) copy_valence_terms(struct.rb_torsions, struct.rb_torsion_types, self.rb_torsions, self.rb_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'improper', 'ignore_end']) copy_valence_terms(struct.urey_bradleys, struct.urey_bradley_types, self.urey_bradleys, self.urey_bradley_types, ['atom1', 'atom2']) copy_valence_terms(struct.impropers, struct.improper_types, self.impropers, self.improper_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(struct.cmaps, struct.cmap_types, self.cmaps, self.cmap_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) copy_valence_terms(struct.trigonal_angles, struct.trigonal_angle_types, self.trigonal_angles, self.trigonal_angle_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(struct.out_of_plane_bends, struct.out_of_plane_bend_types, self.out_of_plane_bends, self.out_of_plane_bend_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(struct.pi_torsions, struct.pi_torsion_types, self.pi_torsions, self.pi_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5', 'atom6']) copy_valence_terms(struct.stretch_bends, struct.stretch_bend_types, self.stretch_bends, self.stretch_bend_types, ['atom1', 'atom2', 'atom3']) copy_valence_terms(struct.torsion_torsions, struct.torsion_torsion_types, self.torsion_torsions, self.torsion_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) copy_valence_terms(struct.chiral_frames, [], self.chiral_frames, [], ['atom1', 'atom2', 'chirality']) copy_valence_terms(struct.multipole_frames, [], self.multipole_frames, [], ['atom', 'frame_pt_num', 'vectail', 'vechead', 'nvec']) copy_valence_terms(struct.adjusts, struct.adjust_types, self.adjusts, self.adjust_types, ['atom1', 'atom2']) copy_valence_terms(struct.donors, [], self.donors, [], ['atom1', 'atom2']) copy_valence_terms(struct.acceptors, [], self.acceptors, [], ['atom1', 'atom2']) copy_valence_terms(struct.groups, [], self.groups, [], ['atom', 'type', 'move']) struct._box = self._box struct.symmetry = self.symmetry struct.space_group = self.space_group return struct def _get_selection_array(self, selection): """ Private method to convert a selection into an array -- common use for viewing and regular selection Parameters ---------- selection : selector (slice, tuple, ... etc.) The selection given to the [] operator """ from parmed.amber import AmberMask # Now we select a subset of atoms. Convert "selection" into a natom list # with 0s and 1s, depending on what the input selection is if isinstance(selection, string_types): mask = AmberMask(self, selection) selection = mask.Selection() elif isinstance(selection, slice): sel = [0 for a in self.atoms] for idx in list(range(len(self.atoms)))[selection]: sel[idx] = 1 selection = sel elif isinstance(selection, tuple) and len(selection) in (2, 3): # This is a selection of chains, and/or residues, and/or atoms if len(selection) == 2: # Select just residues and atoms. Special-case a single-atom # selection for speed -- orders of magnitude improvement in # efficiency ressel, atomsel = selection if isinstance(ressel, integer_types) and isinstance(atomsel, integer_types): return self.residues[ressel][atomsel] has_chain = False elif len(selection) == 3: chainsel, ressel, atomsel = selection chainmap = defaultdict(TrackedList) for r in self.residues: chainmap[r.chain].append(r) if (isinstance(chainsel, string_types) and isinstance(ressel, integer_types) and isinstance(atomsel, integer_types)): # Special-case single-atom selection for efficiency chainmap = dict(chainmap) # no longer defaultdict try: return chainmap[chainsel][ressel][atomsel] except KeyError: raise IndexError('No chain %s in Structure' % chainsel) if isinstance(chainsel, string_types): if chainsel in chainmap: chainset = set([chainsel]) else: # The selected chain is not in the system. Cannot # possibly select *any* atoms. Bail now for speed return None elif isinstance(chainsel, slice): # Build an ordered set of chains chains = [self.residues[0].chain] for res in self.residues: if res.chain != chains[-1]: chains.append(res.chain) chainset = set(chains[chainsel]) else: chainset = set(chainsel) for chain in chainset: if chain in chainmap: break else: # No requested chain is present. Bail now for speed return None has_chain = True # Residue selection can either be by name or index if isinstance(ressel, slice): resset = set(list(range(len(self.residues)))[ressel]) elif isinstance(ressel, string_types) or isinstance(ressel, integer_types): resset = set([ressel]) else: resset = set(ressel) if isinstance(atomsel, slice): atomset = set(list(range(len(self.atoms)))[atomsel]) elif isinstance(atomsel, string_types) or isinstance(atomsel, integer_types): atomset = set([atomsel]) else: atomset = set(atomsel) if has_chain: try: # To allow us to index the atoms residues from their list of # chains, temporarily have the chainmap lists claim the # residues. This must be reversed or the Structure will be # broken for chain_name, chain in iteritems(chainmap): chain.claim() selection = [ (a.residue.chain in chainset) and (a.residue.name in resset or a.residue.idx in resset) and (a.name in atomset or a.idx-a.residue[0].idx in atomset) for a in self.atoms ] finally: # Make sure that the residues list *always* reclaims its # contents self.residues.claim() else: selection = [ (a.name in atomset or a.idx-a.residue[0].idx in atomset) and (a.residue.name in resset or a.residue.idx in resset) for a in self.atoms ] else: # Assume it is an iterable. If it is the same length as the atoms, # it is a boolean mask array. Otherwise, it is a list of atom # indices to select sel = [0 for atom in self.atoms] selection = list(selection) if len(selection) == len(self.atoms): for i, val in enumerate(selection): if val: sel[i] = 1 elif len(selection) > len(self.atoms): raise ValueError('Selection iterable is too long') else: try: for val in selection: sel[val] = 1 except IndexError: raise ValueError('Selected atom out of range') selection = sel # Make sure we have an integer array of 0s and 1s return [int(bool(x)) for x in selection] #=================================================== @property def view(self): """ Returns an indexable object that can be indexed like a standard Structure, but returns a *view* rather than a copy See Also -------- Structure.__getitem__ """ return _StructureViewerCreator(self) #===================================================
[docs] def split(self): """ Split the current Structure into separate Structure instances for each unique molecule. A molecule is defined as all atoms connected by a graph of covalent bonds. Returns ------- [structs, counts] : list of (:class:`Structure`, list) tuples List of all molecules in the order that they appear first in the parent structure accompanied by the list of the molecule numbers in which that molecule appears in the Structure """ tag_molecules(self) mollist = [atom.marked for atom in self.atoms] nmol = max(mollist) structs = [] counts = [] res_molecules = dict() # Divvy up the atoms into their separate molecules molatoms = [[] for i in range(nmol)] for atom in self.atoms: molatoms[atom.marked-1].append(atom) for i in range(nmol): sel = molatoms[i] involved_residues = set(atom.residue.idx for atom in sel) # Shortcut -- keep names of single-residue molecules and the number # of atoms present in those residues in a dict and use that to see # if two molecules are unique -- this gives drastic speedup for # systems with many duplicated single-residue molecules # (i.e., just about every solvent out there) if len(involved_residues) == 1: res = sel[0].residue names = tuple(a.name for a in res) charges = tuple('%.6f' % a.charge for a in res) rmins = tuple('%.6f' % a.rmin for a in res) epsilons = tuple('%.6f' % a.epsilon for a in res) oneres_key = (res.name, len(res), names, charges, rmins, epsilons) if oneres_key in res_molecules: counts[res_molecules[oneres_key]].add(i) continue is_duplicate = False for j, struct in enumerate(structs): if len(struct.atoms) == len(sel): for a1, a2 in zip(struct.atoms, sel): assert None not in (a1.residue, a2.residue), 'Residues must all be set' if a1.residue.name != a2.residue.name: break if a1.name != a2.name: break if '%.6f' % a1.charge != '%.6f' % a2.charge: break if '%.6f' % a1.rmin != '%.6f' % a2.rmin: break if '%.6f' % a1.epsilon != '%.6f' % a2.epsilon: break else: counts[j].add(i) is_duplicate = True break if not is_duplicate: if len(involved_residues) == 1: # Cache single-residue molecules res_molecules[oneres_key] = len(structs) mol = self[[atom.marked == i+1 for atom in self.atoms]] structs.append(mol) counts.append(set([i])) return list(zip(structs, counts))
#===================================================
[docs] def save(self, fname, format=None, overwrite=False, **kwargs): """ Saves the current Structure in the requested file format. Supported formats can be specified explicitly or determined by file-name extension. The following formats are supported, with the recognized suffix and ``format`` keyword shown in parentheses: - PDB (.pdb, pdb) - PDBx/mmCIF (.cif, cif) - PQR (.pqr, pqr) - Amber topology file (.prmtop/.parm7, amber) - CHARMM PSF file (.psf, psf) - CHARMM coordinate file (.crd, charmmcrd) - Gromacs topology file (.top, gromacs) - Gromacs GRO file (.gro, gro) - Mol2 file (.mol2, mol2) - Mol3 file (.mol3, mol3) - Amber ASCII restart (.rst7/.inpcrd/.restrt, rst7) - Amber NetCDF restart (.ncrst, ncrst) Parameters ---------- fname : str or file-like object Name of the file or file-like object to save. If ``format`` is ``None`` (see below), the file type will be determined based on the filename extension. If ``fname`` is file-like object, ``format`` must be provided. If the type cannot be determined, a ValueError is raised. format : str, optional The case-insensitive keyword specifying what type of file ``fname`` should be saved as. If ``None`` (default), the file type will be determined from filename extension of ``fname`` overwrite : bool, optional If True, allow the target file to be overwritten. Otherwise, an IOError is raised if the file exists. Default is False kwargs : keyword-arguments Remaining arguments are passed on to the file writing routines that are called by this function Raises ------ ValueError if either filename extension or ``format`` are not recognized TypeError if the structure cannot be converted to the desired format for whatever reason IOError if the file cannot be written either because it exists and ``overwrite`` is ``False``, the filesystem is read-only, or write permissions are not granted for the user """ from parmed import amber, charmm, formats, gromacs extmap = { '.pdb' : 'PDB', '.pqr' : 'PQR', '.cif' : 'CIF', '.pdbx' : 'CIF', '.parm7' : 'AMBER', '.prmtop' : 'AMBER', '.psf' : 'PSF', '.top' : 'GROMACS', '.gro' : 'GRO', '.mol2' : 'MOL2', '.mol3' : 'MOL3', '.crd' : 'CHARMMCRD', '.rst7' : 'RST7', '.inpcrd' : 'RST7', '.restrt' : 'RST7', '.ncrst' : 'NCRST', } # Basically everybody uses atom type names instead of type indexes. So # convert to atom type names and switch back if need be if not hasattr(fname, 'write'): if os.path.exists(fname) and not overwrite: raise IOError('%s exists; not overwriting' % fname) else: if format is None: raise RuntimeError('Must provide supported format if using file-like object') all_ints = True for atom in self.atoms: if (isinstance(atom.type, integer_types) and atom.atom_type is not UnassignedAtomType): atom.type = str(atom.atom_type) else: all_ints = False try: if format is not None: format = format.upper() else: base, ext = os.path.splitext(fname) if ext in ('.bz2', '.gz'): ext = os.path.splitext(base)[1] try: format = extmap[ext] except KeyError: raise ValueError('Could not determine file type of %s' % fname) # Dispatch if format == 'PDB': self.write_pdb(fname, **kwargs) elif format == 'CIF': self.write_cif(fname, **kwargs) elif format == 'PQR': formats.PQRFile.write(self, fname, **kwargs) elif format == 'PSF': s = charmm.CharmmPsfFile.from_structure(self) s.write_psf(fname, **kwargs) elif format == 'GRO': gromacs.GromacsGroFile.write(self, fname, **kwargs) elif format == 'MOL2': formats.Mol2File.write(self, fname, **kwargs) elif format == 'MOL3': formats.Mol2File.write(self, fname, mol3=True, **kwargs) elif format == 'GROMACS': s = gromacs.GromacsTopologyFile.from_structure(self) s.write(fname, **kwargs) elif format == 'CHARMMCRD': charmm.CharmmCrdFile.write(self, fname, **kwargs) elif format == 'AMBER': if (self.trigonal_angles or self.out_of_plane_bends or self.torsion_torsions or self.pi_torsions or self.stretch_bends or self.chiral_frames or self.multipole_frames): s = amber.AmoebaParm.from_structure(self) s.write_parm(fname, **kwargs) elif self.urey_bradleys or self.impropers or self.cmaps: s = amber.ChamberParm.from_structure(self) s.write_parm(fname, **kwargs) else: try: s = amber.AmberParm.from_structure(self) except TypeError as e: if 'Cannot translate exceptions' in str(e): s = amber.ChamberParm.from_structure(self) else: raise # pragma: no cover s.write_parm(fname, **kwargs) elif format in ('RST7', 'NCRST'): rst7 = amber.Rst7(natom=len(self.atoms), **kwargs) rst7.coordinates = self.coordinates rst7.vels = self.velocities rst7.box = self.box rst7.write(fname, netcdf=(format == 'NCRST')) else: raise ValueError('No file type matching %s' % format) finally: if all_ints: for atom in self.atoms: atom.type = int(atom.atom_type)
#===================================================
[docs] def join_dihedrals(self): """ Joins multi-term torsions into a single term and makes all of the parameters DihedralTypeList instances. If any dihedrals are *already* DihedralTypeList instances, or any are not parametrized, or there are no dihedral_types, this method returns without doing anything """ if not self.dihedral_types: return # nothing to do if any(isinstance(t, DihedralTypeList) for t in self.dihedral_types): return # already done if any(d.type is None for d in self.dihedrals): return # Not fully parametrized dihedrals_to_delete = list() dihedrals_processed = dict() new_dihedral_types = TrackedList() for i, d in enumerate(self.dihedrals): if d.atom1 < d.atom4: key = (d.atom1, d.atom2, d.atom3, d.atom4) else: key = (d.atom4, d.atom3, d.atom2, d.atom1) if key in dihedrals_processed: dihedrals_processed[key].append(d.type) dihedrals_to_delete.append(i) else: dihedrals_processed[key] = dtl = DihedralTypeList() dtl.append(d.type) new_dihedral_types.append(dtl) d.type = dtl # Now drop the new dihedral types into place self.dihedral_types = new_dihedral_types # Remove the "duplicate" dihedrals for i in reversed(dihedrals_to_delete): self.dihedrals[i].delete() del self.dihedrals[i]
#=================================================== @property def combining_rule(self): return self._combining_rule @combining_rule.setter def combining_rule(self, thing): if thing not in ('lorentz', 'geometric'): raise ValueError("combining_rule must be 'lorentz' or 'geometric'") self._combining_rule = thing #=================================================== @property @needs_openmm def topology(self): """ The OpenMM Topology object. Cached when possible, but any changes to the Structure instance results in the topology being deleted and rebuilt Notes ----- This function calls ``prune_empty_terms`` if any topology lists have changed. """ if not self.is_changed(): try: return self._topology except AttributeError: pass else: self.prune_empty_terms() self.unchange() self._topology = top = app.Topology() chain = top.addChain() try: last_chain = self.residues[0].chain last_residue = None last_omm_residue = None except IndexError: # Empty topology return self._topology # Add the atoms for i, atom in enumerate(self.atoms): # See if we need to add a new residue if atom.residue is not last_residue: # See if we need a new chain if last_chain != atom.residue.chain: last_chain = atom.residue.chain chain = top.addChain() last_residue = atom.residue last_omm_residue = top.addResidue(atom.residue.name, chain) try: elem = app.element.Element.getByAtomicNumber(atom.atomic_number) except KeyError: elem = None top.addAtom(atom.name, elem, last_omm_residue) # Add the bonds atoms = list(top.atoms()) for bond in self.bonds: top.addBond(atoms[bond.atom1.idx], atoms[bond.atom2.idx]) # Set the unit cell dimensions if self.box is not None: top.setPeriodicBoxVectors( reducePeriodicBoxVectors( box_lengths_and_angles_to_vectors(*self.box) ) ) return top #=================================================== @property def positions(self): """ A list of 3-element Quantity tuples of dimension length representing the atomic positions for every atom in the system. If set with unitless numbers, those numbers are assumed to be in angstroms. If any atoms do not have coordinates, this is simply ``None``. """ try: return [Vec3(a.xx,a.xy,a.xz) for a in self.atoms] * u.angstroms except AttributeError: return None @positions.setter def positions(self, value): """ A list of 3-element Quantity tuples of dimension length representing the atomic positions for every atom in the system. If set with unitless numbers, those numbers are assumed to be in angstroms. Raises ------ ValueError if the positions are not either a (natom, 3)-shape iterable of iterables or a (natom*3)-length iterable. """ if u.is_quantity(value): value = value.value_in_unit(u.angstroms) # See if the array is flattened if len(value) == len(self.atoms): # It had better all be 3-length iterables for i, atom in enumerate(self.atoms): atom.xx, atom.xy, atom.xz = value[i] elif len(value) == 3 * len(self.atoms): for i, atom in enumerate(self.atoms): i3 = i * 3 atom.xx, atom.xy, atom.xz = value[i3:i3+3] else: raise ValueError('Wrong shape for position array') @property def coordinates(self): try: coords = [[a.xx, a.xy, a.xz] for a in self.atoms] except AttributeError: return None return np.array(coords) @coordinates.setter def coordinates(self, value): """ Setting coordinates will also set xx, xy, and xz on the atoms """ if value is None: # Wipe out coordinates self._coordinates = None for atom in self.atoms: try: del atom.xx, atom.xy, atom.xz except AttributeError: pass else: if u.is_quantity(value): value = value.value_in_unit(u.angstroms) value = list(value) coords = np.array(value, dtype=np.float64, copy=False, subok=True) coords = coords.reshape((-1, len(self.atoms), 3)) if len(coords) > 0: for a, xyz in zip(self.atoms, coords[0]): a.xx, a.xy, a.xz = xyz self._coordinates = coords else: # We set coordinates to an empty iterable. Wipe them out here for a in self.atoms: try: del a.xx, a.xy, a.xz except AttributeError: pass self._coordinates = None
[docs] def get_coordinates(self, frame='all'): """ In some cases, multiple conformations may be stored in the Structure. This function retrieves a particular frame's coordinates Parameters ---------- frame : int or 'all', optional The frame number whose coordinates should be retrieved. Default is 'all' Returns ------- coords : np.ndarray, shape([#,] natom, 3) or None If frame is 'all', all coordinates are returned with shape (#, natom, 3). Otherwise the requested frame is returned with shape (natom, 3). If no coordinates exist and 'all' is requested, None is returned Raises ------ IndexError if there are fewer than ``frame`` coordinates """ try: coords = [[a.xx, a.xy, a.xz] for a in self.atoms] except AttributeError: coords = None else: coords = np.array(coords) # Wipe out our existing coordinates if we think anything might have changed if self._coordinates is not None: if coords is None: self._coordinates = None elif coords.shape != self._coordinates.shape[1:]: self._coordinates = None elif np.abs(coords - self._coordinates[0]).max() > SMALL: self._coordinates = None if frame == 'all': if self._coordinates is not None: return self._coordinates elif coords is not None: return coords.reshape((1, len(self.atoms), 3)) return None elif self._coordinates is None: if frame == 0 and coords is not None: return coords # Requested NOT the first frame raise IndexError('No coordinate frames present') return self._coordinates[frame]
@property def box(self): if self._box is None: return None return self._box[0] @box.setter def box(self, value): if value is None: self._box = None else: if isinstance(value, np.ndarray): box = value else: box = _strip_box_units(list(value)) box = np.array(box, dtype=np.float64, copy=False, subok=True) if box.shape != (6,): if len(box.shape) != 2 or box.shape[-1] != 6: raise ValueError('Box information must be 6 floats') self._box = box.reshape((-1, 6))
[docs] def get_box(self, frame='all'): """ In some cases, multiple conformations may be stored in the Structure. This function retrieves a particular frame's unit cell (box) dimensions Parameters ---------- frame : int or 'all', optional The frame number whose unit cell should be retrieved. Default is 'all' Returns ------- box : np.ndarray, shape([#,] 6) or None If frame is 'all', all unit cells are returned with shape (#, 6). Otherwise the requested frame is returned with shape (6,). If no unit cell exist and 'all' is requested, None is returned Raises ------ IndexError if there are fewer than ``frame`` unit cell dimensions """ if self._box is None: if frame != 'all': raise IndexError('No unit cell frames present') return None elif frame == 'all': return self._box else: return self._box[frame]
#=================================================== @property def velocities(self): """ A (natom, 3)-shape numpy array with atomic velocities for every atom in the system (in units of angstrom/picosecond), or None if there are no velocities """ try: return np.array([[a.vx, a.vy, a.vz] for a in self.atoms]) except AttributeError: return None @velocities.setter def velocities(self, value): """ A list of 3-element Quantity tuples of dimension length representing the atomic velocities for every atom in the system """ if u.is_quantity(value): value = value.value_in_unit(u.angstroms/u.picoseconds) if value is None: for atom in self.atoms: try: del atom.vx, atom.vy, atom.vz except AttributeError: pass else: value = np.array(value, copy=False).reshape((-1,len(self.atoms),3)) for atom, xyz in zip(self.atoms, value[0]): atom.vx, atom.vy, atom.vz = xyz #===================================================
[docs] def has_NBFIX(self): """ Returns whether or not any pairs of atom types have their LJ interactions modified by an NBFIX definition Returns ------- has_nbfix : bool If True, at least two atom types have NBFIX mod definitions """ typemap = dict() for a in self.atoms: if a.atom_type is UnassignedAtomType: continue typemap[str(a.atom_type)] = a.atom_type # Now we have a map of all atom types that we have defined in our # system. Look through all of the atom types and see if any of their # NBFIX definitions are also keys in typemap for key, type in iteritems(typemap): for key in type.nbfix: if key in typemap: return True return False
#=================================================== @property def box_vectors(self): """ 3, 3-element tuple of unit cell vectors that are Quantity objects of dimension length """ if self._box is None: return None return box_lengths_and_angles_to_vectors(*self.box) @box_vectors.setter def box_vectors(self, value): """ 3, 3-element tuple of unit cell vectors that are Quantity objects of dimension length """ (a, b, c), (A, B, G) = box_vectors_to_lengths_and_angles(*value) a = a.value_in_unit(u.angstroms) b = b.value_in_unit(u.angstroms) c = c.value_in_unit(u.angstroms) A = A.value_in_unit(u.degrees) B = B.value_in_unit(u.degrees) G = G.value_in_unit(u.degrees) self._box = np.array([[a, b, c, A, B, G]], dtype=np.float64) #===================================================
[docs] @needs_openmm def createSystem(self, nonbondedMethod=None, nonbondedCutoff=8.0*u.angstroms, switchDistance=0.0*u.angstroms, constraints=None, rigidWater=True, implicitSolvent=None, implicitSolventKappa=None, implicitSolventSaltConc=0.0*u.moles/u.liters, temperature=298.15*u.kelvin, soluteDielectric=1.0, solventDielectric=78.5, useSASA=False, removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005, flexibleConstraints=True, verbose=False, splitDihedrals=False): """ Construct an OpenMM System representing the topology described by the prmtop file. Parameters ---------- nonbondedMethod : cutoff method This is the cutoff method. It can be either the NoCutoff, CutoffNonPeriodic, CutoffPeriodic, PME, or Ewald objects from the simtk.openmm.app namespace nonbondedCutoff : float or distance Quantity The nonbonded cutoff must be either a floating point number (interpreted as nanometers) or a Quantity with attached units. This is ignored if nonbondedMethod is NoCutoff. switchDistance : float or distance Quantity The distance at which the switching function is turned on for van der Waals interactions. This is ignored when no cutoff is used, and no switch is used if switchDistance is 0, negative, or greater than the cutoff constraints : None, app.HBonds, app.HAngles, or app.AllBonds Which type of constraints to add to the system (e.g., SHAKE). None means no bonds are constrained. HBonds means bonds with hydrogen are constrained rigidWater : bool=True If True, water is kept rigid regardless of the value of constraints. A value of False is ignored if constraints is not None. implicitSolvent : None, app.HCT, app.OBC1, app.OBC2, app.GBn, app.GBn2 The Generalized Born implicit solvent model to use. implicitSolventKappa : float or 1/distance Quantity = None This is the Debye kappa property related to modeling saltwater conditions in GB. It should have units of 1/distance (1/nanometers is assumed if no units present). A value of None means that kappa will be calculated from implicitSolventSaltConc (below) implicitSolventSaltConc : float or amount/volume Quantity=0 moles/liter If implicitSolventKappa is None, the kappa will be computed from the salt concentration. It should have units compatible with mol/L temperature : float or temperature Quantity = 298.15 kelvin This is only used to compute kappa from implicitSolventSaltConc soluteDielectric : float=1.0 The dielectric constant of the protein interior used in GB solventDielectric : float=78.5 The dielectric constant of the water used in GB useSASA : bool=False If True, use the ACE non-polar solvation model. Otherwise, use no SASA-based nonpolar solvation model. removeCMMotion : bool=True If True, the center-of-mass motion will be removed periodically during the simulation. If False, it will not. hydrogenMass : float or mass quantity = None If not None, hydrogen masses will be changed to this mass and the difference subtracted from the attached heavy atom (hydrogen mass repartitioning) ewaldErrorTolerance : float=0.0005 When using PME or Ewald, the Ewald parameters will be calculated from this value flexibleConstraints : bool=True If False, the energies and forces from the constrained degrees of freedom will NOT be computed. If True, they will (but those degrees of freedom will *still* be constrained). verbose : bool=False If True, the progress of this subroutine will be printed to stdout splitDihedrals : bool=False If True, the dihedrals will be split into two forces -- proper and impropers. This is primarily useful for debugging torsion parameter assignments. Notes ----- This function calls prune_empty_terms if any Topology lists have changed """ if self.unknown_functional: raise ParameterError('Cannot createSystem: unknown functional') # Establish defaults if nonbondedMethod is None: nonbondedMethod = app.NoCutoff system = mm.System() # Make sure periodic simulations have a box if nonbondedMethod in (app.CutoffPeriodic, app.PME, app.Ewald): if self.box is None: raise ValueError('No periodic boundary conditions detected') # Do hydrogen mass repartitioning if necessary masses = [atom.mass for atom in self.atoms] if hydrogenMass is not None: if u.is_quantity(hydrogenMass): hydrogenMass = hydrogenMass.value_in_unit(u.dalton) if hydrogenMass <= 0: raise ValueError('Hydrogen mass must be positive') for atom in self.atoms: if atom.element == 1: heavy_atom = None for a2 in atom.bond_partners: if a2.element != 1: heavy_atom = a2 break if heavy_atom is not None: masses[atom.idx] = hydrogenMass masses[heavy_atom.idx] -= hydrogenMass - atom.mass for mass in masses: system.addParticle(mass) self.omm_add_constraints(system, constraints, rigidWater) # Prune empty terms if we have changed if self.is_changed(): self.prune_empty_terms() self.unchange() # Add the various types of forces LOGGER.info('Adding bonds...') self._add_force_to_system(system, self.omm_bond_force(constraints, rigidWater, flexibleConstraints) ) LOGGER.info('Adding angles...') self._add_force_to_system(system, self.omm_angle_force(constraints, flexibleConstraints)) LOGGER.info('Adding dihedrals...') self._add_force_to_system(system, self.omm_dihedral_force(splitDihedrals)) LOGGER.info('Adding Ryckaert-Bellemans torsions...') self._add_force_to_system(system, self.omm_rb_torsion_force()) LOGGER.info('Adding Urey-Bradleys...') self._add_force_to_system(system, self.omm_urey_bradley_force()) LOGGER.info('Adding improper torsions...') self._add_force_to_system(system, self.omm_improper_force()) LOGGER.info('Adding CMAP torsions...') self._add_force_to_system(system, self.omm_cmap_force()) LOGGER.info('Adding trigonal angle terms...') self._add_force_to_system(system, self.omm_trigonal_angle_force()) LOGGER.info('Adding out-of-plane bends...') self._add_force_to_system(system, self.omm_out_of_plane_bend_force()) LOGGER.info('Adding pi-torsions...') self._add_force_to_system(system, self.omm_pi_torsion_force()) LOGGER.info('Adding stretch-bends...') self._add_force_to_system(system, self.omm_stretch_bend_force()) LOGGER.info('Adding torsion-torsions...') self._add_force_to_system(system, self.omm_torsion_torsion_force()) LOGGER.info('Adding Nonbonded force...') if implicitSolvent is not None: rf_dielc = 1.0 else: rf_dielc = 78.5 self._add_force_to_system(system, self.omm_nonbonded_force(nonbondedMethod, nonbondedCutoff, switchDistance, ewaldErrorTolerance, rf_dielc) ) if implicitSolvent is not None: LOGGER.info('Adding GB force...') self._add_force_to_system(system, self.omm_gbsa_force(implicitSolvent, nonbondedMethod, nonbondedCutoff, soluteDielectric, solventDielectric, implicitSolventKappa, implicitSolventSaltConc, temperature, useSASA) ) if removeCMMotion: system.addForce(mm.CMMotionRemover()) if self.box is not None: system.setDefaultPeriodicBoxVectors(*reducePeriodicBoxVectors(self.box_vectors)) self.omm_set_virtual_sites(system) return system
#===================================================
[docs] @needs_openmm def omm_add_constraints(self, system, constraints, rigidWater): """ Adds constraints to a given system Parameters ---------- system : mm.System The OpenMM system for which constraints should be added constraints : None, app.HBonds, app.AllBonds, or app.HAngles Which kind of constraints should be used rigidWater : bool If True, water bonds are constrained regardless of whether constrains is None """ if constraints is None and not rigidWater: return if constraints not in (None, app.HBonds, app.AllBonds, app.HAngles): raise ValueError("Unrecognized constraints option (%s)" % constraints) length_conv = u.angstrom.conversion_factor_to(u.nanometer) constraint_bond_set = set() constraint_angle_set = set() is_water = _settler(self) if constraints is app.AllBonds or constraints is app.HAngles: for bond in self.bonds: # Skip all extra points... don't constrain those if isinstance(bond.atom1, ExtraPoint): continue if isinstance(bond.atom2, ExtraPoint): continue constraint_bond_set.add(frozenset((bond.atom1.idx, bond.atom2.idx))) elif constraints is app.HBonds: for bond in self.bonds: if isinstance(bond.atom1, ExtraPoint): continue if isinstance(bond.atom2, ExtraPoint): continue if bond.atom1.element == 1 or bond.atom2.element == 1: constraint_bond_set.add(frozenset((bond.atom1.idx, bond.atom2.idx))) if rigidWater: for bond in self.bonds: if isinstance(bond.atom1, ExtraPoint): continue if isinstance(bond.atom2, ExtraPoint): continue if is_water[bond.atom1.residue.idx]: constraint_bond_set.add(frozenset((bond.atom1.idx, bond.atom2.idx))) # Add bond constraints for bond in self.bonds: if frozenset((bond.atom1.idx, bond.atom2.idx)) in constraint_bond_set: system.addConstraint(bond.atom1.idx, bond.atom2.idx, bond.type.req * length_conv) if constraints is app.HAngles: for angle in self.angles: numH = 0 if angle.atom1.element == 1: numH += 1 if angle.atom3.element == 1: numH += 1 if numH == 2 or (numH == 1 and angle.atom2.element == 8): constraint_angle_set.add((angle.atom1.idx, angle.atom2.idx, angle.atom3.idx)) if rigidWater: for angle in self.angles: if is_water[angle.atom1.residue.idx]: constraint_angle_set.add((angle.atom1.idx, angle.atom2.idx, angle.atom3.idx)) # Add angle constraints for angle in self.angles: if (angle.atom1.idx, angle.atom2.idx, angle.atom3.idx) in constraint_angle_set: if frozenset((angle.atom1.idx, angle.atom3.idx)) in constraint_bond_set: continue # Constrain this angle l1 = l2 = None for bond in angle.atom2.bonds: if bond in angle and angle.atom1 in bond: l1 = bond.type.req * length_conv elif bond in angle and angle.atom3 in bond: l2 = bond.type.req * length_conv # Law of cosines to find the constraint distance if l1 is None or l2 is None: continue # no bonds found... cost = math.cos(angle.type.theteq * DEG_TO_RAD) length = math.sqrt(l1 * l1 + l2 * l2 - 2 * l1 * l2 * cost) system.addConstraint(angle.atom1.idx, angle.atom3.idx, length)
#===================================================
[docs] @needs_openmm def omm_set_virtual_sites(self, system): """ Sets the virtual sites in a given OpenMM `System` object from the extra points defined in this system Parameters ---------- system : mm.System The system for which the virtual sites will be set. All particles must have already been added to this System before calling this method """ if system.getNumParticles() != len(self.atoms): raise ValueError('OpenMM System does not correspond to Structure') for atom in self.atoms: if not isinstance(atom, ExtraPoint): continue # This is a virtual site... get its frame type typ = atom.frame_type weights = typ.get_weights() refatoms = typ.get_atoms() if isinstance(typ, TwoParticleExtraPointFrame): a1, a2 = refatoms w1, w2 = weights system.setVirtualSite(atom.idx, mm.TwoParticleAverageSite(a1.idx, a2.idx, w1, w2)) elif isinstance(typ, ThreeParticleExtraPointFrame): a1, a2, a3 = refatoms w1, w2, w3 = weights system.setVirtualSite( atom.idx, mm.ThreeParticleAverageSite(a1.idx, a2.idx, a3.idx, w1, w2, w3) ) elif isinstance(typ, OutOfPlaneExtraPointFrame): a1, a2, a3 = refatoms w1, w2, w3 = weights system.setVirtualSite( atom.idx, mm.OutOfPlaneSite(a1.idx, a2.idx, a3.idx, w1, w2, w3) )
#===================================================
[docs] @needs_openmm def omm_bond_force(self, constraints=None, rigidWater=True, flexibleConstraints=True): """ Creates an OpenMM Bond Force object (or AmoebaBondForce if the bonds are for an Amoeba-parametrized system) Parameters ---------- constraints : None, app.HBonds, app.AllBonds, or app.HAngles The types of constraints that are on the system. If flexibleConstraints is False, then the constrained bonds will not be added to the resulting Force rigidWater : bool=True Should water-H bonds be constrained regardless of `constraints`? flexibleConstraints : bool=True If True, all bonds are added to the force regardless of `constraints` Returns ------- force HarmonicBondForce (or AmoebaBondForce if this is an Amoeba system), or None if there are no bonds to add """ if not flexibleConstraints and constraints in (app.HAngles, app.AllBonds) or not self.bonds: return None # No bonds to add length_conv = u.angstroms.conversion_factor_to(u.nanometers) _ambfrc = u.kilocalorie_per_mole/u.angstrom**2 _ommfrc = u.kilojoule_per_mole/u.nanometer**2 frc_conv = _ambfrc.conversion_factor_to(_ommfrc) # See if we need to add Amoeba bonds or regular bonds if hasattr(self.bond_types, 'degree') and hasattr(self.bond_types, 'coeffs'): force = mm.AmoebaBondForce() force.setAmoebaGlobalBondCubic(self.bond_types.coeffs[3]/length_conv) force.setAmoebaGlobalBondQuartic(self.bond_types.coeffs[4]/length_conv**2) else: force = mm.HarmonicBondForce() force.setForceGroup(self.BOND_FORCE_GROUP) # Add the bonds if rigidWater: is_water = _settler(self) else: is_water = [False for r in self.residues] for bond in self.bonds: if 1 in (bond.atom1.element, bond.atom2.element) and ( not flexibleConstraints and constraints is app.HBonds): continue if not flexibleConstraints and is_water[bond.atom1.residue.idx]: # is_water is False even for waters if rigidWater is False, so # we can rely on is_water to be correct for our needs # regardless continue if bond.type is None: raise ParameterError('Cannot find necessary parameters') force.addBond(bond.atom1.idx, bond.atom2.idx, bond.type.req*length_conv, 2*bond.type.k*frc_conv) # Done adding the force if force.getNumBonds() == 0: return None return force
#===================================================
[docs] @needs_openmm def omm_angle_force(self, constraints=None, flexibleConstraints=True): """ Creates an OpenMM HarmonicAngleForce object (or AmoebaAngleForce if the angles are for an Amoeba-parametrized system) Parameters ---------- constraints : None, app.HBonds, app.AllBonds, or app.HAngles The types of constraints that are on the system. If flexibleConstraints is False, then the constrained bonds will not be added to the resulting Force flexibleConstraints : bool=True If True, all bonds are added to the force regardless of `constraints` Returns ------- force HarmonicAngleForce (or AmoebaAngleForce if this is an Amoeba system), or None if there are no angles to add """ if not self.angles: return None frc_conv = u.kilocalories.conversion_factor_to(u.kilojoules) if (hasattr(self.angle_types, 'degree') and hasattr(self.angle_types, 'coeffs')): c = self.angle_types.coeffs force = mm.AmoebaAngleForce() force.setAmoebaGlobalAngleCubic(c[3]) force.setAmoebaGlobalAngleQuartic(c[4]) force.setAmoebaGlobalAnglePentic(c[5]) force.setAmoebaGlobalAngleSextic(c[6]) else: force = mm.HarmonicAngleForce() force.setForceGroup(self.ANGLE_FORCE_GROUP) for angle in self.angles: num_h = (angle.atom1.element == 1) + (angle.atom3.element == 1) if constraints is app.HAngles and (num_h == 2 or (num_h == 1 and angle.atom2.element == 8) and not flexibleConstraints): continue # pragma: no cover if angle.type is None: raise ParameterError('Cannot find angle parameters') force.addAngle(angle.atom1.idx, angle.atom2.idx, angle.atom3.idx, angle.type.theteq*DEG_TO_RAD, 2*angle.type.k*frc_conv) if force.getNumAngles() == 0: return None return force
#===================================================
[docs] @needs_openmm def omm_dihedral_force(self, split=False): """ Creates the OpenMM PeriodicTorsionForce modeling dihedrals Parameters ---------- split : bool, optional, default=False If True, separate PeriodicTorsionForce instances with the propers in the first and impropers in the second return item. If no impropers or propers are present, the instances with zero terms are not returned. Returns ------- PeriodicTorsionForce[, PeriodicTorsionForce] Or None if no torsions are present in this system """ if not self.dihedrals: return None frc_conv = u.kilocalories.conversion_factor_to(u.kilojoules) proper = mm.PeriodicTorsionForce() improper = mm.PeriodicTorsionForce() proper.setForceGroup(self.DIHEDRAL_FORCE_GROUP) improper.setForceGroup(self.IMPROPER_FORCE_GROUP) for tor in self.dihedrals: if tor.type is None: raise ParameterError('Cannot find torsion parameters') if split and tor.improper: force = improper else: force = proper if isinstance(tor.type, DihedralTypeList): for typ in tor.type: force.addTorsion(tor.atom1.idx, tor.atom2.idx, tor.atom3.idx, tor.atom4.idx, int(typ.per), typ.phase*DEG_TO_RAD, typ.phi_k*frc_conv) else: force.addTorsion(tor.atom1.idx, tor.atom2.idx, tor.atom3.idx, tor.atom4.idx, int(tor.type.per), tor.type.phase*DEG_TO_RAD, tor.type.phi_k*frc_conv) if proper.getNumTorsions() == 0: return improper if improper.getNumTorsions() == 0: return proper return proper, improper
#===================================================
[docs] @needs_openmm def omm_rb_torsion_force(self): """ Creates the OpenMM RBTorsionForce for Ryckaert-Bellemans torsions Returns ------- RBTorsionForce Or None if no torsions are present in this system """ if not self.rb_torsions: return None conv = u.kilocalories.conversion_factor_to(u.kilojoules) force = mm.RBTorsionForce() force.setForceGroup(self.RB_TORSION_FORCE_GROUP) for tor in self.rb_torsions: if tor.type is None: raise ParameterError('Cannot find R-B torsion parameters') force.addTorsion(tor.atom1.idx, tor.atom2.idx, tor.atom3.idx, tor.atom4.idx, tor.type.c0*conv, tor.type.c1*conv, tor.type.c2*conv, tor.type.c3*conv, tor.type.c4*conv, tor.type.c5*conv) return force
#===================================================
[docs] @needs_openmm def omm_urey_bradley_force(self): """ Creates the OpenMM Urey-Bradley force Returns ------- HarmonicBondForce Or None, if no urey-bradleys are present """ if not self.urey_bradleys: return None length_conv = u.angstroms.conversion_factor_to(u.nanometers) _ambfrc = u.kilocalorie_per_mole/u.angstrom**2 _ommfrc = u.kilojoule_per_mole/u.nanometer**2 frc_conv = _ambfrc.conversion_factor_to(_ommfrc) force = mm.HarmonicBondForce() force.setForceGroup(self.UREY_BRADLEY_FORCE_GROUP) for urey in self.urey_bradleys: if urey.type is None: raise ParameterError('Cannot find urey-bradley parameters') force.addBond(urey.atom1.idx, urey.atom2.idx, urey.type.req*length_conv, 2*urey.type.k*frc_conv) return force
#===================================================
[docs] @needs_openmm def omm_improper_force(self): """ Creates the OpenMM improper torsion force (quadratic bias) Returns ------- CustomTorsionForce With the formula k*(phi-phi0)^2, or None if there are no impropers """ if not self.impropers: return None frc_conv = u.kilocalories.conversion_factor_to(u.kilojoules) energy_function = 'k*dtheta_torus^2;' energy_function += 'dtheta_torus = dtheta - floor(dtheta/(2*pi)+0.5)*(2*pi);' energy_function += 'dtheta = theta - theta0;' energy_function += 'pi = %f;' % math.pi force = mm.CustomTorsionForce(energy_function) force.addPerTorsionParameter('k') force.addPerTorsionParameter('theta0') force.setForceGroup(self.IMPROPER_FORCE_GROUP) for imp in self.impropers: if imp.type is None: raise ParameterError('Cannot find improper torsion parameters') force.addTorsion(imp.atom1.idx, imp.atom2.idx, imp.atom3.idx, imp.atom4.idx, (imp.type.psi_k*frc_conv, imp.type.psi_eq*DEG_TO_RAD)) return force
#===================================================
[docs] @needs_openmm def omm_cmap_force(self): """ Creates the OpenMM CMAP torsion force Returns ------- CMAPTorsionForce Or None, if no CMAP terms are present """ if not self.cmaps: return None frc_conv = u.kilocalories.conversion_factor_to(u.kilojoules) force = mm.CMAPTorsionForce() force.setForceGroup(self.CMAP_FORCE_GROUP) # First get the list of cmap maps we're going to use. Just store the IDs # so we have simple integer comparisons to do later cmap_type_list = [] cmap_map = dict() for cmap in self.cmaps: if cmap.type is None: raise ParameterError('Cannot find CMAP torsion parameters') if not id(cmap.type) in cmap_type_list: ct = cmap.type cmap_type_list.append(id(ct)) # OpenMM takes the correction maps in the range 0 to 360 # degrees, but we store it in -180 -- 180 (and in the transpose # of what OpenMM expects). grid = ct.grid.switch_range().T m = force.addMap(ct.resolution, [x*frc_conv for x in grid]) cmap_map[id(ct)] = m # Now add all of the cmaps for cmap in self.cmaps: force.addTorsion(cmap_map[id(cmap.type)], cmap.atom1.idx, cmap.atom2.idx, cmap.atom3.idx, cmap.atom4.idx, cmap.atom2.idx, cmap.atom3.idx, cmap.atom4.idx, cmap.atom5.idx) return force
#===================================================
[docs] @needs_openmm def omm_nonbonded_force(self, nonbondedMethod=None, nonbondedCutoff=8*u.angstroms, switchDistance=0*u.angstroms, ewaldErrorTolerance=0.0005, reactionFieldDielectric=78.5): """ Creates the OpenMM NonbondedForce instance Parameters ---------- nonbondedMethod : cutoff method This is the cutoff method. It can be either the NoCutoff, CutoffNonPeriodic, CutoffPeriodic, PME, or Ewald objects from the simtk.openmm.app namespace nonbondedCutoff : float or distance Quantity The nonbonded cutoff must be either a floating point number (interpreted as nanometers) or a Quantity with attached units. This is ignored if nonbondedMethod is NoCutoff. switchDistance : float or distance Quantity The distance at which the switching function is turned on for van der Waals interactions. This is ignored when no cutoff is used, and no switch is used if switchDistance is 0, negative, or greater than the cutoff ewaldErrorTolerance : float=0.0005 When using PME or Ewald, the Ewald parameters will be calculated from this value reactionFieldDielectric : float=78.5 If the nonbondedMethod is CutoffPeriodic or CutoffNonPeriodic, the region beyond the cutoff is treated using a reaction field method with this dielectric constant. It should be set to 1 if another implicit solvent model is being used (e.g., GB) Returns ------- NonbondedForce This just implements the very basic NonbondedForce with the typical charge-charge and 12-6 Lennard-Jones interactions with the Lorentz-Berthelot combining rules. Notes ----- Subclasses of Structure for which this nonbonded treatment is inadequate should override this method to implement what is needed. If nrexcl is set to 3 and no exception parameters are stored in the adjusts list, the 1-4 interactions are determined from the list of dihedrals """ if not self.atoms: return None length_conv = u.angstrom.conversion_factor_to(u.nanometer) ene_conv = u.kilocalories.conversion_factor_to(u.kilojoules) force = mm.NonbondedForce() force.setForceGroup(self.NONBONDED_FORCE_GROUP) if u.is_quantity(nonbondedCutoff): nonbondedCutoff = nonbondedCutoff.value_in_unit(u.nanometers) if nonbondedMethod is None or nonbondedMethod is app.NoCutoff: force.setNonbondedMethod(mm.NonbondedForce.NoCutoff) elif nonbondedMethod is app.CutoffNonPeriodic: force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic) force.setCutoffDistance(nonbondedCutoff) elif nonbondedMethod is app.CutoffPeriodic: force.setNonbondedMethod(mm.NonbondedForce.CutoffPeriodic) force.setCutoffDistance(nonbondedCutoff) elif nonbondedMethod is app.PME: force.setNonbondedMethod(mm.NonbondedForce.PME) force.setCutoffDistance(nonbondedCutoff) force.setEwaldErrorTolerance(ewaldErrorTolerance) elif nonbondedMethod is app.Ewald: force.setNonbondedMethod(mm.NonbondedForce.Ewald) force.setCutoffDistance(nonbondedCutoff) force.setEwaldErrorTolerance(ewaldErrorTolerance) else: raise ValueError('Unrecognized nonbondedMethod (%s)' % nonbondedMethod) force.setReactionFieldDielectric(reactionFieldDielectric) # Now add the particles sigma_scale = length_conv * 2 * 2**(-1/6) for atom in self.atoms: force.addParticle(atom.charge, atom.sigma*length_conv, abs(atom.epsilon*ene_conv)) # Add exclusions from the bond graph out to nrexcl-1 bonds away (atoms # nrexcl bonds away will be exceptions defined later) def exclude_to(origin, atom, level, end): if level >= end: return for partner in atom.bond_partners: if partner is origin: continue if atom is not origin: force.addException(origin.idx, atom.idx, 0.0, 0.5, 0.0, True) # Exclude EP children, too for child in origin.children: force.addException(partner.idx, child.idx, 0.0, 0.5, 0.0, True) for child in partner.children: force.addException(child.idx, origin.idx, 0.0, 0.5, 0.0, True) for child in origin.children: for child2 in partner.children: force.addException(child.idx, child2.idx, 0.0, 0.5, 0.0, True) exclude_to(origin, partner, level+1, end) for atom in self.atoms: exclude_to(atom, atom, 0, self.nrexcl) # Add the exceptions from the dihedral list IFF no explicit exceptions # (or *adjusts*) have been specified. If dihedral.ignore_end is False, a # 1-4 with the appropriate scaling factor is used as the exception. sigma_scale = 2**(-1/6) * length_conv if self.combining_rule == 'lorentz': comb_sig = lambda sig1, sig2: 0.5 * (sig1 + sig2) elif self.combining_rule == 'geometric': comb_sig = lambda sig1, sig2: math.sqrt(sig1 * sig2) if not self.adjusts: for dih in self.dihedrals + self.rb_torsions: if dih.ignore_end: continue if isinstance(dih.type, DihedralTypeList): scee = scnb = 0 i = 0 while (scee == 0 or scnb == 0) and i < len(dih.type): scee = dih.type[i].scee scnb = dih.type[i].scnb i += 1 if scee == 0 or scnb == 0: raise ValueError('Detected scaling constants of 0 for ' 'dihedral containing 1-4 info!') else: scee = dih.type.scee scnb = dih.type.scnb try: rij, wdij, rij14, wdij14 = dih.atom1.atom_type.nbfix[str(dih.atom4.atom_type)] except (KeyError, AttributeError): epsprod = abs(dih.atom1.epsilon_14 * dih.atom4.epsilon_14) epsprod = math.sqrt(epsprod) * ene_conv / scnb sigprod = comb_sig(dih.atom1.sigma_14, dih.atom4.sigma_14) sigprod *= length_conv else: epsprod = wdij14 * ene_conv / scnb sigprod = rij14 * length_conv * sigma_scale chgprod = dih.atom1.charge * dih.atom4.charge / scee force.addException(dih.atom1.idx, dih.atom4.idx, chgprod, sigprod, epsprod, True) for child in dih.atom1.children: epsprod = abs(child.epsilon_14 * dih.atom4.epsilon_14) epsprod = math.sqrt(epsprod) * ene_conv / scnb sigprod = comb_sig(child.sigma_14, dih.atom4.sigma_14) sigprod *= length_conv chgprod = (child.charge * dih.atom4.charge) / scee force.addException(child.idx, dih.atom4.idx, chgprod, sigprod, epsprod, True) for child in dih.atom4.children: epsprod = abs(child.epsilon_14 * dih.atom1.epsilon_14) epsprod = math.sqrt(epsprod) * ene_conv / scnb sigprod = comb_sig(child.sigma_14, dih.atom1.sigma_14) sigprod *= length_conv chgprod = child.charge * dih.atom1.charge / scee force.addException(child.idx, dih.atom1.idx, chgprod, sigprod, epsprod, True) for c1 in dih.atom1.children: for c2 in dih.atom2.children: epsprod = abs(c1.epsilon_14 * c2.epsilon_14) epsprod = math.sqrt(epsprod) * ene_conv / scnb sigprod = comb_sig(c1.sigma_14, c2.sigma_14) sigprod *= length_conv chgprod = c1.charge * c2.charge / scee force.addException(c1.idx, c2.idx, chgprod, sigprod, epsprod, True) # Allow our specific exceptions (in adjusts) to override anything that # came before for pair in self.adjusts: chgprod = pair.atom1.charge * pair.atom2.charge * pair.type.chgscale force.addException(pair.atom1.idx, pair.atom2.idx, chgprod, pair.type.sigma*length_conv, pair.type.epsilon*ene_conv, True) # Any exclusion partners we already added will zero-out existing # exclusions/exceptions for atom in self.atoms: for a2 in atom.exclusion_partners: force.addException(atom.idx, a2.idx, 0.0, 0.5, 0.0, True) if switchDistance and nonbondedMethod is not app.NoCutoff: if u.is_quantity(switchDistance): switchDistance = switchDistance.value_in_unit(u.nanometers) if 0 < switchDistance < nonbondedCutoff: force.setUseSwitchingFunction(True) force.setSwitchingDistance(switchDistance) # Now see if we have any NBFIXes that we need to implement via a # CustomNonbondedForce if self.combining_rule == 'geometric': if self.has_NBFIX(): return (force, self._omm_nbfixed_force(force, nonbondedMethod), self._omm_geometric_force(force, nonbondedMethod)) else: return force, self._omm_geometric_force(force, nonbondedMethod) if self.has_NBFIX(): return force, self._omm_nbfixed_force(force, nonbondedMethod) return force
#=================================================== def _omm_nbfixed_force(self, nonbfrc, nonbondedMethod): """ Private method for creating a CustomNonbondedForce with a lookup table. This should not be called by users -- you have been warned. Parameters ---------- nonbfrc : NonbondedForce NonbondedForce for the "standard" nonbonded interactions. This will be modified (specifically, L-J ixns will be zeroed) nonbondedMethod : Nonbonded Method (e.g., NoCutoff, PME, etc.) The nonbonded method to apply here. Ewald and PME will be interpreted as CutoffPeriodic for the CustomNonbondedForce. Returns ------- force : CustomNonbondedForce The L-J force with NBFIX implemented as a lookup table """ length_conv = u.angstroms.conversion_factor_to(u.nanometers) ene_conv = u.kilocalories.conversion_factor_to(u.kilojoules) # We need a CustomNonbondedForce to implement the NBFIX functionality. # First derive the type lookup tables lj_idx_list = [0 for atom in self.atoms] lj_radii, lj_depths = [], [] num_lj_types = 0 lj_type_list = [] for i, atom in enumerate(self.atoms): atype = atom.atom_type if lj_idx_list[i]: continue # already assigned num_lj_types += 1 lj_idx_list[i] = num_lj_types ljtype = (atype.rmin, abs(atype.epsilon)) lj_type_list.append(atype) lj_radii.append(atype.rmin) lj_depths.append(abs(atype.epsilon)) for j in range(i+1, len(self.atoms)): atype2 = self.atoms[j].atom_type if lj_idx_list[j] > 0: continue # already assigned if atype2 is atype: lj_idx_list[j] = num_lj_types elif not atype.nbfix: # Only non-NBFIXed atom types can be compressed ljtype2 = (atype2.rmin, abs(atype2.epsilon)) if ljtype == ljtype2: lj_idx_list[j] = num_lj_types # Now everything is assigned. Create the A- and B-coefficient arrays acoef = [0 for i in range(num_lj_types*num_lj_types)] bcoef = acoef[:] for i in range(num_lj_types): for j in range(num_lj_types): namej = lj_type_list[j].name try: rij, wdij, rij14, wdij14 = lj_type_list[i].nbfix[namej] except KeyError: rij = (lj_radii[i] + lj_radii[j]) * length_conv wdij = math.sqrt(lj_depths[i] * lj_depths[j]) * ene_conv else: rij *= length_conv wdij *= ene_conv rij6 = rij**6 acoef[i+num_lj_types*j] = math.sqrt(wdij) * rij6 bcoef[i+num_lj_types*j] = 2 * wdij * rij6 force = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r2*r2*r2; r2=r^2; ' 'a=acoef(type1, type2); b=bcoef(type1, type2)') force.addTabulatedFunction('acoef', mm.Discrete2DFunction(num_lj_types, num_lj_types, acoef)) force.addTabulatedFunction('bcoef', mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef)) force.addPerParticleParameter('type') force.setForceGroup(self.NONBONDED_FORCE_GROUP) if (nonbondedMethod is app.PME or nonbondedMethod is app.Ewald or nonbondedMethod is app.CutoffPeriodic): force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic) elif nonbondedMethod is app.NoCutoff: force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff) elif nonbondedMethod is app.CutoffNonPeriodic: force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic) else: raise AssertionError('Unrecognized nonbonded method [%s]' % nonbondedMethod) # Add the particles for i in lj_idx_list: force.addParticle((i-1,)) # Now wipe out the L-J parameters in the nonbonded force for i in range(nonbfrc.getNumParticles()): chg, sig, eps = nonbfrc.getParticleParameters(i) nonbfrc.setParticleParameters(i, chg, 0.5, 0.0) # Now transfer the exclusions for ii in range(nonbfrc.getNumExceptions()): i, j, qq, ss, ee = nonbfrc.getExceptionParameters(ii) force.addExclusion(i, j) # Now transfer the other properties (cutoff, switching function, etc.) force.setUseLongRangeCorrection(True) if nonbondedMethod is app.NoCutoff: force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff) elif nonbondedMethod is app.CutoffNonPeriodic: force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic) elif nonbondedMethod in (app.PME, app.Ewald, app.CutoffPeriodic): force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic) else: raise AssertionError('Unsupported nonbonded method %s' % nonbondedMethod) force.setCutoffDistance(nonbfrc.getCutoffDistance()) if nonbfrc.getUseSwitchingFunction(): force.setUseSwitchingFunction(True) force.setSwitchingDistance(nonbfrc.getSwitchingDistance()) return force #=================================================== def _omm_geometric_force(self, nonbfrc, nonbondedMethod): """ Private method for creating a CustomNonbondedForce with a lookup table. This should not be called by users -- you have been warned. Parameters ---------- nonbfrc : NonbondedForce NonbondedForce for the "standard" nonbonded interactions. This will be modified (specifically, L-J ixns will be zeroed) nonbondedMethod : Nonbonded Method (e.g., NoCutoff, PME, etc.) The nonbonded method to apply here. Ewald and PME will be interpreted as CutoffPeriodic for the CustomNonbondedForce. Returns ------- force : CustomNonbondedForce The L-J force with NBFIX implemented as a lookup table """ length_conv = u.angstroms.conversion_factor_to(u.nanometers) ene_conv = u.kilocalories.conversion_factor_to(u.kilojoules) # We need a CustomNonbondedForce to implement the geometric combining # rules force = mm.CustomNonbondedForce('epsilon1*epsilon2*(sigr6^2-sigr6); sigr6=sigr2*sigr2*sigr2; ' 'sigr2=(sigc/r)^2; sigc=sigma1*sigma2') force.addPerParticleParameter('epsilon') force.addPerParticleParameter('sigma') force.setForceGroup(self.NONBONDED_FORCE_GROUP) if (nonbondedMethod is app.PME or nonbondedMethod is app.Ewald or nonbondedMethod is app.CutoffPeriodic): force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic) elif nonbondedMethod is app.NoCutoff: force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff) elif nonbondedMethod is app.CutoffNonPeriodic: force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic) else: raise AssertionError('Unrecognized nonbonded method [%s]' % nonbondedMethod) # Add the particles for atom in self.atoms: eps = math.sqrt(atom.epsilon*ene_conv) * 2 sig = math.sqrt(atom.sigma*length_conv) force.addParticle((eps, sig)) # Now wipe out the L-J parameters in the nonbonded force for i in range(nonbfrc.getNumParticles()): chg, sig, eps = nonbfrc.getParticleParameters(i) nonbfrc.setParticleParameters(i, chg, 0.5, 0.0) # Now transfer the exclusions for ii in range(nonbfrc.getNumExceptions()): i, j, qq, ss, ee = nonbfrc.getExceptionParameters(ii) force.addExclusion(i, j) # Now transfer the other properties (cutoff, switching function, etc.) force.setUseLongRangeCorrection(True) if nonbondedMethod is app.NoCutoff: force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff) elif nonbondedMethod is app.CutoffNonPeriodic: force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic) elif nonbondedMethod in (app.PME, app.Ewald, app.CutoffPeriodic): force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic) else: raise AssertionError('Unsupported nonbonded method %s' % nonbondedMethod) force.setCutoffDistance(nonbfrc.getCutoffDistance()) if nonbfrc.getUseSwitchingFunction(): force.setUseSwitchingFunction(True) force.setSwitchingDistance(nonbfrc.getSwitchingDistance()) return force #===================================================
[docs] @needs_openmm def omm_gbsa_force(self, implicitSolvent, nonbondedMethod=None, nonbondedCutoff=30.0*u.angstroms, soluteDielectric=1.0, solventDielectric=78.5, implicitSolventKappa=None, implicitSolventSaltConc=0.0*u.moles/u.liter, temperature=298.15*u.kelvin, useSASA=True): """ Creates a Generalized Born force for running implicit solvent calculations Parameters ---------- implicitSolvent : app.HCT, app.OBC1, app.OBC2, app.GBn, app.GBn2 The Generalized Born implicit solvent model to use. nonbondedMethod : cutoff method This is the cutoff method. It can be either the NoCutoff, CutoffNonPeriodic, CutoffPeriodic, PME, or Ewald objects from the simtk.openmm.app namespace. Default is NoCutoff nonbondedCutoff : float or distance Quantity The nonbonded cutoff must be either a floating opint number (interpreted as nanometers) or a Quantity with attached units. This is ignored if nonbondedMethod is NoCutoff implicitSolventKappa : float or 1/distance Quantity = None This is the Debye kappa property related to modeling saltwater conditions in GB. It should have units of 1/distance (1/nanometers is assumed if no units present). A value of None means that kappa will be calculated from implicitSolventSaltConc (below) implicitSolventSaltConc : float or amount/volume Quantity=0 moles/liter If implicitSolventKappa is None, the kappa will be computed from the salt concentration. It should have units compatible with mol/L temperature : float or temperature Quantity = 298.15 kelvin This is only used to compute kappa from implicitSolventSaltConc soluteDielectric : float=1.0 The dielectric constant of the protein interior used in GB solventDielectric : float=78.5 The dielectric constant of the water used in GB """ from simtk.openmm.app.internal.customgbforces import (GBSAHCTForce, GBSAOBC1Force, GBSAOBC2Force, GBSAGBnForce, GBSAGBn2Force) try: from simtk.openmm.app.internal.customgbforces import convertParameters except ImportError: # Unnecessary in newer versions of OpenMM convertParameters = lambda params, choice: params if implicitSolvent is None: return None if useSASA: sasa = 'ACE' else: sasa = None if nonbondedMethod is None: nonbondedMethod = app.NoCutoff if implicitSolvent not in (app.HCT, app.OBC1, app.OBC2, app.GBn, app.GBn2): raise ValueError('Unrecognized implicit solvent model') gb_parms = convertParameters(self._get_gb_parameters(implicitSolvent), str(implicitSolvent)) if implicitSolventKappa is None: if u.is_quantity(implicitSolventSaltConc): sc = implicitSolventSaltConc.value_in_unit(u.moles/u.liter) implicitSolventSaltConc = sc if u.is_quantity(temperature): temperature = temperature.value_in_unit(u.kelvin) # The constant is 1 / sqrt(eps_0 * kB / (2*NA*q^2*1000)) where NA is # Avogadro's number, eps_0 is the permittivity of free space, q is # the charge (this # matches Amber's conversion factor) implicitSolventKappa = 50.33355 * math.sqrt(implicitSolventSaltConc / solventDielectric / temperature) # Multiply by 0.73 to account for ion exclusions, and multiply by 10 # to convert to 1/nm from 1/angstroms implicitSolventKappa *= 7.3 elif u.is_quantity(implicitSolventKappa): implicitSolventKappa = implicitSolventKappa.value_in_unit(u.nanometer**-1) if nonbondedMethod is app.NoCutoff: cutoff = None elif u.is_quantity(nonbondedCutoff): cutoff = nonbondedCutoff.value_in_unit(u.nanometers) else: cutoff = nonbondedCutoff if implicitSolvent is app.HCT: force = GBSAHCTForce(solventDielectric, soluteDielectric, sasa, cutoff, kappa=implicitSolventKappa) elif implicitSolvent is app.OBC1: force = GBSAOBC1Force(solventDielectric, soluteDielectric, sasa, cutoff, kappa=implicitSolventKappa) elif implicitSolvent is app.OBC2: force = GBSAOBC2Force(solventDielectric, soluteDielectric, sasa, cutoff, kappa=implicitSolventKappa) elif implicitSolvent is app.GBn: force = GBSAGBnForce(solventDielectric, soluteDielectric, sasa, cutoff, kappa=implicitSolventKappa) elif implicitSolvent is app.GBn2: force = GBSAGBn2Force(solventDielectric, soluteDielectric, sasa, cutoff, kappa=implicitSolventKappa) else: raise AssertionError('Unexpected implicit solvent model... should not be here') for atom, parms in zip(self.atoms, gb_parms): force.addParticle([atom.charge] + list(parms)) try: force.finalize() except AttributeError: # only new versions of omm require calling finalize pass # Set cutoff method if nonbondedMethod is app.NoCutoff: force.setNonbondedMethod(mm.CustomGBForce.NoCutoff) elif nonbondedMethod is app.CutoffNonPeriodic: force.setNonbondedMethod(mm.CustomGBForce.CutoffNonPeriodic) force.setCutoffDistance(cutoff) else: # cutoff periodic (PME, CutoffPeriodic, Ewald) force.setNonbondedMethod(mm.CustomGBForce.CutoffPeriodic) force.setCutoffDistance(cutoff) force.setForceGroup(self.NONBONDED_FORCE_GROUP) return force
#=================================================== # Amoeba-specific forces below
[docs] @needs_openmm def omm_trigonal_angle_force(self): """ Creates the Amoeba trigonal-angle force Returns ------- AmoebaInPlaneAngleForce The trigonal in-plane Angle force """ if not self.trigonal_angles: return None frc_conv = u.kilocalories.conversion_factor_to(u.kilojoules) if (not hasattr(self.trigonal_angle_types, 'degree') or not hasattr(self.trigonal_angle_types, 'coeffs')): raise ParameterError('Do not have the trigonal angle force table parameters') force = mm.AmoebaInPlaneAngleForce() c = self.trigonal_angle_types.coeffs force.setAmoebaGlobalInPlaneAngleCubic(c[3]) force.setAmoebaGlobalInPlaneAngleQuartic(c[4]) force.setAmoebaGlobalInPlaneAnglePentic(c[5]) force.setAmoebaGlobalInPlaneAngleSextic(c[6]) force.setForceGroup(self.TRIGONAL_ANGLE_FORCE_GROUP) for ang in self.trigonal_angles: if ang.type is None: raise ParameterError('Missing trigonal angle parameters') force.addAngle(ang.atom1.idx, ang.atom2.idx, ang.atom3.idx, ang.atom4.idx, ang.type.theteq, ang.type.k*frc_conv) return force
#===================================================
[docs] @needs_openmm def omm_out_of_plane_bend_force(self): """ Creates the Amoeba out-of-plane bend force Returns ------- AmoebaOutOfPlaneBendForce The out-of-plane bend Angle force """ if not self.out_of_plane_bends: return None frc_conv = u.kilocalories.conversion_factor_to(u.kilojoules) if (not hasattr(self.out_of_plane_bend_types, 'degree') or not hasattr(self.out_of_plane_bend_types, 'coeffs')): raise ParameterError('Do not have the trigonal angle force table parameters') force = mm.AmoebaOutOfPlaneBendForce() c = self.out_of_plane_bend_types.coeffs force.setAmoebaGlobalOutOfPlaneBendCubic(c[3]) force.setAmoebaGlobalOutOfPlaneBendQuartic(c[4]) force.setAmoebaGlobalOutOfPlaneBendPentic(c[5]) force.setAmoebaGlobalOutOfPlaneBendSextic(c[6]) force.setForceGroup(self.OUT_OF_PLANE_BEND_FORCE_GROUP) for ang in self.out_of_plane_bends: if ang.type is None: raise ParameterError('Missing out-of-plane bend parameters') force.addOutOfPlaneBend(ang.atom1.idx, ang.atom2.idx, ang.atom3.idx, ang.atom4.idx, 2*ang.type.k*frc_conv) return force
#===================================================
[docs] @needs_openmm def omm_pi_torsion_force(self): """ Creates the Amoeba pi-torsion force Returns ------- AmoebaPiTorsionForce The pi-torsion force """ if not self.pi_torsions: return None frc_conv = u.kilocalories.conversion_factor_to(u.kilojoules) force = mm.AmoebaPiTorsionForce() force.setForceGroup(self.PI_TORSION_FORCE_GROUP) for ang in self.pi_torsions: if ang.type is None: raise ParameterError('Missing pi-torsion parameters') force.addPiTorsion(ang.atom1.idx, ang.atom2.idx, ang.atom3.idx, ang.atom4.idx, ang.atom5.idx, ang.atom6.idx, ang.type.phi_k*frc_conv) return force
#===================================================
[docs] @needs_openmm def omm_stretch_bend_force(self): """ Create the OpenMM Amoeba stretch-bend force for this system Returns ------- AmoebaStretchBendForce The stretch-bend force containing all terms in this system """ if not self.stretch_bends: return None # Conversion factor taken from pyopenmm/processTinkerForceField.py frc_conv = math.pi / 180 * 41.84 # 4.184 * 10 length_conv = u.angstroms.conversion_factor_to(u.nanometers) force = mm.AmoebaStretchBendForce() force.setForceGroup(self.STRETCH_BEND_FORCE_GROUP) for strbnd in self.stretch_bends: if strbnd.type is None: raise ParameterError("Missing stretch-bend parameters") force.addStretchBend(strbnd.atom1.idx, strbnd.atom2.idx, strbnd.atom3.idx, strbnd.type.req1*length_conv, strbnd.type.req2*length_conv, strbnd.type.theteq*math.pi/180, strbnd.type.k1*frc_conv, strbnd.type.k2*frc_conv) return force
#===================================================
[docs] @needs_openmm def omm_torsion_torsion_force(self): """ Create the OpenMM Amoeba coupled-torsion (CMAP) force Returns ------- AmoebaTorsionTorsionForce The torsion-torsion (CMAP) force with all coupled-torsion parameters for this system """ if not self.torsion_torsions: return None # Not implemented yet... raise NotImplementedError("Torsion-torsions not yet implemented")
#=================================================== def _prune_empty_bonds(self): """ Gets rid of any empty bonds """ for i in reversed(range(len(self.bonds))): bond = self.bonds[i] if bond.atom1 is None and bond.atom2 is None: del self.bonds[i] elif bond.atom1.idx == -1 or bond.atom2.idx == -1: bond.delete() del self.bonds[i] #=================================================== def _prune_empty_angles(self): """ Gets rid of any empty angles """ for i in reversed(range(len(self.angles))): angle = self.angles[i] if (angle.atom1 is None and angle.atom2 is None and angle.atom3 is None): del self.angles[i] elif -1 in (angle.atom1.idx, angle.atom2.idx, angle.atom3.idx): angle.delete() del self.angles[i] #=================================================== def _prune_empty_dihedrals(self, dlist='dihedrals'): """ Gets rid of any empty dihedrals """ for i in reversed(range(len(getattr(self, dlist)))): dihed = getattr(self, dlist)[i] if (dihed.atom1 is None and dihed.atom2 is None and dihed.atom3 is None and dihed.atom4 is None): del getattr(self, dlist)[i] elif -1 in (dihed.atom1.idx, dihed.atom2.idx, dihed.atom3.idx, dihed.atom4.idx): dihed.delete() del getattr(self, dlist)[i] #=================================================== def _prune_empty_rb_torsions(self): """ Gets rid of any empty R-B torsions """ self._prune_empty_dihedrals('rb_torsions') #=================================================== def _prune_empty_ureys(self): """ Gets rid of any empty Urey-Bradley terms """ for i in reversed(range(len(self.urey_bradleys))): ub = self.urey_bradleys[i] if ub.atom1 is None and ub.atom2 is None: del self.urey_bradleys[i] elif ub.atom1.idx == -1 or ub.atom2.idx == -1: ub.delete() del self.urey_bradleys[i] #=================================================== def _prune_empty_impropers(self): """ Gets rid of any empty improper torsions """ for i in reversed(range(len(self.impropers))): imp = self.impropers[i] if imp.atom1 is None and imp.atom2 is None and imp.atom3 is None and imp.atom4 is None: del self.impropers[i] elif -1 in (imp.atom1.idx, imp.atom2.idx, imp.atom3.idx, imp.atom4.idx): imp.delete() del self.impropers[i] #=================================================== def _prune_empty_cmaps(self): """ Gets rid of any empty CMAP terms """ for i in reversed(range(len(self.cmaps))): cmap = self.cmaps[i] if (cmap.atom1 is None and cmap.atom2 is None and cmap.atom3 is None and cmap.atom4 is None and cmap.atom5 is None): del self.cmaps[i] elif -1 in (cmap.atom1.idx, cmap.atom2.idx, cmap.atom3.idx, cmap.atom4.idx, cmap.atom5.idx): cmap.delete() del self.cmaps[i] #=================================================== def _prune_empty_trigonal_angles(self): """ Gets rid of any empty trigonal angles """ self._prune_empty_dihedrals('trigonal_angles') #=================================================== def _prune_empty_out_of_plane_bends(self): """ Gets rid of any empty out-of-plane bends """ self._prune_empty_dihedrals('out_of_plane_bends') #=================================================== def _prune_empty_pi_torsions(self): """ Gets rid of any empty pi-torsions """ for i in reversed(range(len(self.pi_torsions))): pit = self.pi_torsions[i] if (pit.atom1 is None and pit.atom2 is None and pit.atom3 is None and pit.atom4 is None and pit.atom5 is None and pit.atom6 is None): del self.pi_torsions[i] elif -1 in (pit.atom1.idx, pit.atom2.idx, pit.atom3.idx, pit.atom4.idx, pit.atom5.idx, pit.atom6.idx): # Not stored anywhere, no need to call delete() del self.pi_torsions[i] #=================================================== def _prune_empty_stretch_bends(self): """ Gets rid of any empty stretch-bend terms """ for i in reversed(range(len(self.stretch_bends))): sb = self.stretch_bends[i] if sb.atom1 is None and sb.atom2 is None and sb.atom3 is None: del self.stretch_bends[i] elif sb.atom1.idx == -1 or sb.atom2.idx == -1 or sb.atom3.idx == -1: # Not stored anywhere, no need to call delete() del self.stretch_bends[i] #=================================================== def _prune_empty_torsion_torsions(self): """ Gets rid of any empty torsion-torsion terms """ for i in reversed(range(len(self.torsion_torsions))): tt = self.torsion_torsions[i] if (tt.atom1 is None and tt.atom2 is None and tt.atom3 is None and tt.atom4 is None and tt.atom5 is None): del self.torsion_torsions[i] elif -1 in (tt.atom1.idx, tt.atom2.idx, tt.atom3.idx, tt.atom4.idx, tt.atom5.idx): tt.delete() del self.torsion_torsions[i] #=================================================== def _prune_empty_chiral_frames(self): """ Gets rid of any empty chiral frame terms """ for i in reversed(range(len(self.chiral_frames))): cf = self.chiral_frames[i] if cf.atom1 is None or cf.atom2 is None: del self.chiral_frames[i] elif cf.atom1.idx == -1 or cf.atom2.idx == -1: del self.chiral_frames[i] #=================================================== def _prune_empty_multipole_frames(self): """ Gets rid of any empty multipole frame terms """ for i in reversed(range(len(self.multipole_frames))): mf = self.multipole_frames[i] if mf.atom is None or mf.atom.idx == -1: del self.multipole_frames[i] #=================================================== def _prune_empty_adjusts(self): """ Gets rid of any empty nonbonded exception adjustments """ for i in reversed(range(len(self.adjusts))): adj = self.adjusts[i] if adj.atom1 is None or adj.atom2 is None: del self.adjusts[i] elif adj.atom1.idx == -1 or adj.atom2.idx == -1: del self.adjusts[i] #=================================================== @needs_openmm def _get_gb_parameters(self, implicitSolvent): """ Gets the GB parameters for the requested GB model used by OpenMM Parameters ---------- implicitSolvent : app.HCT, app.OBC1, app.OBC2, app.GBn, or app.GBn2 The object specifying a particular GB model in OpenMM Returns ------- parameters : list of float List of parameters for the requested GB model """ if implicitSolvent is app.GBn: screen = [0.5 for atom in self.atoms] radii = [atom.solvent_radius for atom in self.atoms] for i, atom in enumerate(self.atoms): if atom.element == 6: screen[i] = 0.48435382330 elif atom.element == 1: screen[i] = 1.09085413633 elif atom.element == 7: screen[i] = 0.700147318409 elif atom.element == 8: screen[i] = 1.06557401132 elif atom.element == 16: screen[i] = 0.602256336067 if radii[i] == 0: radii[i] = _bondi(atom) elif implicitSolvent is app.GBn2: # Add non-optimized values as defaults alpha = [1.0 for i in self.atoms] beta = [0.8 for i in self.atoms] gamma = [4.85 for i in self.atoms] screen = [0.5 for i in self.atoms] radii = [atom.solvent_radius for atom in self.atoms] for i, atom in enumerate(self.atoms): if atom.element == 6: screen[i] = 1.058554 alpha[i] = 0.733756 beta[i] = 0.506378 gamma[i] = 0.205844 elif atom.element == 1: screen[i] = 1.425952 alpha[i] = 0.788440 beta[i] = 0.798699 gamma[i] = 0.437334 elif atom.element == 7: screen[i] = 0.733599 alpha[i] = 0.503364 beta[i] = 0.316828 gamma[i] = 0.192915 elif atom.element == 8: screen[i] = 1.061039 alpha[i] = 0.867814 beta[i] = 0.876635 gamma[i] = 0.387882 elif atom.element == 16: screen[i] = -0.703469 alpha[i] = 0.867814 beta[i] = 0.876635 gamma[i] = 0.387882 if not radii[i]: radii[i] = _mbondi3(atom) else: radii = [atom.solvent_radius for atom in self.atoms] screen = [atom.screen for atom in self.atoms] for i, atom in enumerate(self.atoms): if not radii[i] or not screen[i]: # Replace with defaults radii[i], screen[i] = _gb_rad_screen(atom, implicitSolvent) length_conv = u.angstrom.conversion_factor_to(u.nanometer) radii = [x * length_conv for x in radii] if implicitSolvent is app.GBn2: return list(zip(radii, screen, alpha, beta, gamma)) return list(zip(radii, screen)) #=================================================== def __str__(self): if hasattr(self, 'name') and self.name: return self.name return repr(self) #=================================================== # Support concatenating and replicating the Structure def __add__(self, other): cp = copy(self) cp += other return cp def __iadd__(self, other): if not isinstance(other, Structure): return NotImplemented # Cache my coordinates, since changing the structure will destroy the # coordinate list mycrd = self.get_coordinates('all') # The basic approach taken here is to extend the atom list, then scan # through all of the valence terms in `other`, adding them to the # corresponding arrays of `self`, using an offset to look into the atom # array. The types are done a little differently. The types are all # copied into a "temporary" array so they have the same index as the # original Structure `other`, and as the valence terms are created the # reference to that type is added. Afterwards, those types are tacked on # to the end of the types list on `self` aoffset = len(self.atoms) try: roffset = self.residues[-1].number + 1 except IndexError: # No residues... must be an empty structure roffset = 0 for atom in other.atoms: res = atom.residue self.add_atom(copy(atom), res.name, res.idx+roffset, res.chain, res.insertion_code, res.segid) def copy_valence_terms(oval, otyp, sval, styp, attrlist): """ Copies the valence terms from one list to another; oval=Other VALence; otyp=Other TYPe; sval=Self VALence; styp=Self TYPe; attrlist=ATTRibute LIST (atom1, atom2, ...) """ otypcp = [copy(typ) for typ in otyp] for val in oval: ats = [getattr(val, attr) for attr in attrlist] # Replace any atoms with the ones from self.atoms for i, at in enumerate(ats): if isinstance(at, Atom): ats[i] = self.atoms[at.idx+aoffset] kws = dict() if hasattr(val, 'type') and val.type is NoUreyBradley: kws['type'] = NoUreyBradley # special-case singleton elif otypcp and val.type is not None: kws['type'] = otypcp[val.type.idx] sval.append(type(val)(*ats, **kws)) if hasattr(val, 'funct'): sval[-1].funct = val.funct # Now tack on the "new" types copied from `other` styp.extend(otypcp) if hasattr(styp, 'claim'): styp.claim() copy_valence_terms(other.bonds, other.bond_types, self.bonds, self.bond_types, ['atom1', 'atom2']) copy_valence_terms(other.angles, other.angle_types, self.angles, self.angle_types, ['atom1', 'atom2', 'atom3']) copy_valence_terms(other.dihedrals, other.dihedral_types, self.dihedrals, self.dihedral_types, ['atom1', 'atom2', 'atom3', 'atom4', 'improper', 'ignore_end']) copy_valence_terms(other.rb_torsions, other.rb_torsion_types, self.rb_torsions, self.rb_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'improper', 'ignore_end']) copy_valence_terms(other.urey_bradleys, other.urey_bradley_types, self.urey_bradleys, self.urey_bradley_types, ['atom1', 'atom2']) copy_valence_terms(other.impropers, other.improper_types, self.impropers, self.improper_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(other.cmaps, other.cmap_types, self.cmaps, self.cmap_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) copy_valence_terms(other.trigonal_angles, other.trigonal_angle_types, self.trigonal_angles, self.trigonal_angle_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(other.out_of_plane_bends, other.out_of_plane_bend_types, self.out_of_plane_bends, self.out_of_plane_bend_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(other.pi_torsions, other.pi_torsion_types, self.pi_torsions, self.pi_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5', 'atom6']) copy_valence_terms(other.stretch_bends, other.stretch_bend_types, self.stretch_bends, self.stretch_bend_types, ['atom1', 'atom2', 'atom3']) copy_valence_terms(other.torsion_torsions, other.torsion_torsion_types, self.torsion_torsions, self.torsion_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) copy_valence_terms(other.chiral_frames, [], self.chiral_frames, [], ['atom1', 'atom2', 'chirality']) copy_valence_terms(other.multipole_frames, [], self.multipole_frames, [], ['atom', 'frame_pt_num', 'vectail', 'vechead', 'nvec']) copy_valence_terms(other.adjusts, other.adjust_types, self.adjusts, self.adjust_types, ['atom1', 'atom2']) copy_valence_terms(other.donors, [], self.donors, [], ['atom1', 'atom2']) copy_valence_terms(other.acceptors, [], self.acceptors, [], ['atom1', 'atom2']) copy_valence_terms(other.groups, [], self.groups, [], ['atom', 'type', 'move']) if mycrd is None or other._coordinates is None: self._coordinates = None else: ocrd = other.get_coordinates('all') nframes = min(ocrd.shape[0], mycrd.shape[0]) self._coordinates = np.concatenate((mycrd[:nframes,:,:], ocrd[:nframes,:,:]), axis=1) return self def __mul__(self, ncopies): """ Replicates the current Structure `ncopies` times """ cp = copy(self) return cp.__imul__(ncopies, self) __rmul__ = __mul__ def __imul__(self, ncopies, other=None): if not isinstance(ncopies, integer_types): return NotImplemented # The basic approach here is similar to what we used in __iadd__, except # we don't have to extend the type arrays at all -- we just point to the # same one that the first copy pointed to. def copy_valence_terms(oval, aoffset, sval, styp, attrlist): """ Copies the valence terms from one list to another; oval=Other VALence; otyp=Other TYPe; sval=Self VALence; styp=Self TYPe; attrlist=ATTRibute LIST (atom1, atom2, ...) """ for val in oval: ats = [getattr(val, attr) for attr in attrlist] # Replace any atoms with the ones from self.atoms for i, at in enumerate(ats): if isinstance(at, Atom): ats[i] = self.atoms[at.idx+aoffset] kws = dict() if hasattr(val, 'type') and val.type is NoUreyBradley: kws['type'] = NoUreyBradley # special-case singleton elif styp and val.type is not None: kws['type'] = styp[val.type.idx] sval.append(type(val)(*ats, **kws)) if hasattr(val, 'funct'): sval[-1].funct = val.funct if other is None: other = copy(self) for i in range(ncopies-1): aoffset = len(self.atoms) roffset = self.residues[-1].number + 1 for atom in other.atoms: res = atom.residue self.add_atom(copy(atom), res.name, res.idx+roffset, res.chain, res.insertion_code, res.segid) copy_valence_terms(other.bonds, aoffset, self.bonds, self.bond_types, ['atom1', 'atom2']) copy_valence_terms(other.angles, aoffset, self.angles, self.angle_types, ['atom1', 'atom2', 'atom3']) copy_valence_terms(other.dihedrals, aoffset, self.dihedrals, self.dihedral_types, ['atom1', 'atom2', 'atom3', 'atom4', 'improper', 'ignore_end']) copy_valence_terms(other.rb_torsions, aoffset, self.rb_torsions, self.rb_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'improper', 'ignore_end']) copy_valence_terms(other.urey_bradleys, aoffset, self.urey_bradleys, self.urey_bradley_types, ['atom1', 'atom2']) copy_valence_terms(other.impropers, aoffset, self.impropers, self.improper_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(other.cmaps, aoffset, self.cmaps, self.cmap_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) copy_valence_terms(other.trigonal_angles, aoffset, self.trigonal_angles, self.trigonal_angle_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(other.out_of_plane_bends, aoffset, self.out_of_plane_bends, self.out_of_plane_bend_types, ['atom1', 'atom2', 'atom3', 'atom4']) copy_valence_terms(other.pi_torsions, aoffset, self.pi_torsions, self.pi_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5', 'atom6']) copy_valence_terms(other.stretch_bends, aoffset, self.stretch_bends, self.stretch_bend_types, ['atom1', 'atom2', 'atom3']) copy_valence_terms(other.torsion_torsions, aoffset, self.torsion_torsions, self.torsion_torsion_types, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) copy_valence_terms(other.chiral_frames, aoffset, self.chiral_frames, [], ['atom1', 'atom2', 'chirality']) copy_valence_terms(other.multipole_frames, aoffset, self.multipole_frames, [], ['atom', 'frame_pt_num', 'vectail', 'vechead', 'nvec']) copy_valence_terms(other.adjusts, aoffset, self.adjusts, self.adjust_types, ['atom1', 'atom2']) copy_valence_terms(other.donors, aoffset, self.donors, [], ['atom1', 'atom2']) copy_valence_terms(other.acceptors, aoffset, self.acceptors, [], ['atom1', 'atom2']) copy_valence_terms(other.groups, aoffset, self.groups, [], ['atom', 'type', 'move']) if self._coordinates is not None: self._coordinates = np.hstack([self._coordinates for i in range(ncopies)]) return self #=================================================== # For the bool-ness of Structure. An empty structure evaluates to # boolean-false, but this means that the structure must have no atoms, # residues, topological features, or parameter types at all. def __bool__(self): return bool(self.atoms or self.residues or self.bonds or self.angles or self.dihedrals or self.impropers or self.rb_torsions or self.cmaps or self.torsion_torsions or self.stretch_bends or self.out_of_plane_bends or self.trigonal_angles or self.torsion_torsions or self.pi_torsions or self.urey_bradleys or self.chiral_frames or self.multipole_frames or self.adjusts or self.acceptors or self.donors or self.groups or self.bond_types or self.angle_types or self.dihedral_types or self.urey_bradley_types or self.improper_types or self.rb_torsion_types or self.cmap_types or self.trigonal_angle_types or self.out_of_plane_bend_types or self.pi_torsion_types or self.torsion_torsion_types or self.adjust_types) __nonzero__ = __bool__ # for Python 2 #=================================================== @staticmethod def _add_force_to_system(system, force): """ Adds an OpenMM force to a system IFF the force is not None """ if force is None: return if isinstance(force, tuple) or isinstance(force, list): # It's possible we got multiple forces to add for f in force: system.addForce(f) return system.addForce(force) #=================================================== def __getstate__(self): """ Serializes a structure """ retdict = dict(residues=self.residues, bond_types=self.bond_types, angle_types=self.angle_types, dihedral_types=self.dihedral_types, urey_bradley_types=self.urey_bradley_types, improper_types=self.improper_types, rb_torsion_types=self.rb_torsion_types, cmap_types=self.cmap_types, trigonal_angle_types=self.trigonal_angle_types, out_of_plane_bend_types=self.out_of_plane_bend_types, pi_torsion_types=self.pi_torsion_types, stretch_bend_types=self.stretch_bend_types, torsion_torsion_types=self.torsion_torsion_types, adjust_types=self.adjust_types, groups=self.groups, _coordinates=self._coordinates, _box=self._box, nrexcl=self.nrexcl, _combining_rule=self._combining_rule, unknown_functional=self.unknown_functional, space_group=self.space_group, ) def idx(thing): return thing.idx if thing is not None else None retdict['bonds'] = [(b.atom1.idx, b.atom2.idx, idx(b.type)) for b in self.bonds] retdict['angles'] = [(a.atom1.idx, a.atom2.idx, a.atom3.idx, idx(a.type)) for a in self.angles] retdict['dihedrals'] = [(d.atom1.idx, d.atom2.idx, d.atom3.idx, d.atom4.idx, d.improper, d.ignore_end, idx(d.type)) for d in self.dihedrals] retdict['impropers'] = [(d.atom1.idx, d.atom2.idx, d.atom3.idx, d.atom4.idx, idx(d.type)) for d in self.impropers] retdict['rb_torsions'] = [(d.atom1.idx, d.atom2.idx, d.atom3.idx, d.atom4.idx, idx(d.type)) for d in self.rb_torsions] retdict['urey_bradleys'] = [(u.atom1.idx, u.atom2.idx, idx(u.type)) for u in self.urey_bradleys] retdict['cmaps'] = [(c.atom1.idx, c.atom2.idx, c.atom3.idx, c.atom4.idx, c.atom5.idx, idx(c.type)) for c in self.cmaps] retdict['trigonal_angles'] = [(t.atom1.idx, t.atom2.idx, t.atom3.idx, t.atom4.idx, idx(t.type)) for t in self.trigonal_angles] retdict['out_of_plane_bends'] = [(o.atom1.idx, o.atom2.idx, o.atom3.idx, o.atom4.idx, idx(o.type)) for o in self.out_of_plane_bends] retdict['pi_torsions'] = [(p.atom1.idx, p.atom2.idx, p.atom3.idx, p.atom4.idx, p.atom5.idx, p.atom6.idx, idx(p.type)) for p in self.pi_torsions] retdict['stretch_bends'] = [(s.atom1.idx, s.atom2.idx, s.atom3.idx, idx(s.type)) for s in self.stretch_bends] retdict['torsion_torsions'] = [(t.atom1.idx, t.atom2.idx, t.atom3.idx, t.atom4.idx, t.atom5.idx, idx(t.type)) for t in self.torsion_torsions] retdict['chiral_frames'] = [(c.atom1.idx, c.atom2.idx, c.chirality) for c in self.chiral_frames] retdict['multipole_frames'] = [(f.atom.idx, f.frame_pt_num, f.vectail, f.vechead, f.nvec) for f in self.multipole_frames] retdict['adjusts'] = [(e.atom1.idx, e.atom2.idx, idx(e.type)) for e in self.adjusts] retdict['acceptors'] = [(a.atom1.idx, a.atom2.idx) for a in self.acceptors] retdict['donors'] = [(d.atom1.idx, d.atom2.idx) for d in self.donors] retdict['exclusions'] = [tuple(e.idx for e in a._exclusion_partners) for a in self.atoms] # Now the metadata stuff, if applicable for key in ('experimental', 'journal', 'authors', 'keywords', 'doi', 'pmid', 'journal_authors', 'volume', 'title', 'year', 'resolution', 'related_entries'): try: retdict[key] = getattr(self, key) except AttributeError: continue return retdict def __setstate__(self, d): # Assign the type lists we need to just copy over, and make sure the # attributes are claimed. for attr in ('residues', 'bond_types', 'angle_types', 'dihedral_types', 'urey_bradley_types', 'improper_types', 'rb_torsion_types', 'cmap_types', 'trigonal_angle_types', 'adjust_types', 'out_of_plane_bend_types', 'pi_torsion_types', 'groups', 'stretch_bend_types', 'torsion_torsion_types'): setattr(self, attr, d[attr]) getattr(self, attr).claim() # Assign the possible metadata for key in ('experimental', 'journal', 'authors', 'keywords', 'doi', 'pmid', 'journal_authors', 'volume', 'title', 'year', 'resolution', 'related_entries', '_coordinates', '_box', 'nrexcl', '_combining_rule', 'unknown_functional', 'space_group'): if key in d: setattr(self, key, d[key]) self.atoms = AtomList() for r in self.residues: self.atoms.extend(r.atoms) def assign_type(typelist, idx): if idx is None: return None return typelist[idx] # Set the topology arrays self.bonds = TrackedList( Bond(self.atoms[it[0]], self.atoms[it[1]], type=assign_type(self.bond_types, it[2])) for it in d['bonds'] ) self.angles = TrackedList( Angle(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], type=assign_type(self.angle_types, it[3])) for it in d['angles'] ) self.dihedrals = TrackedList( Dihedral(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], improper=it[4], ignore_end=it[5], type=assign_type(self.dihedral_types, it[6])) for it in d['dihedrals'] ) self.impropers = TrackedList( Improper(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], type=assign_type(self.improper_types, it[4])) for it in d['impropers'] ) self.urey_bradleys = TrackedList( UreyBradley(self.atoms[it[0]], self.atoms[it[1]], type=assign_type(self.urey_bradley_types, it[2])) for it in d['urey_bradleys'] ) self.rb_torsions = TrackedList( Dihedral(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], type=assign_type(self.rb_torsion_types, it[4])) for it in d['rb_torsions'] ) self.cmaps = TrackedList( Cmap(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], self.atoms[it[4]], type=assign_type(self.cmap_types, it[5])) for it in d['cmaps'] ) self.trigonal_angles = TrackedList( TrigonalAngle(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], type=assign_type(self.trigonal_angle_types, it[4])) for it in d['trigonal_angles'] ) self.out_of_plane_bends = TrackedList( OutOfPlaneBend(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], type=assign_type(self.out_of_plane_bend_types, it[4])) for it in d['out_of_plane_bends'] ) self.pi_torsions = TrackedList( PiTorsion(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], self.atoms[it[4]], self.atoms[it[5]], type=assign_type(self.pi_torsion_types, it[6])) for it in d['pi_torsions'] ) self.stretch_bends = TrackedList( StretchBend(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], type=assign_type(self.stretch_bend_types, it[3])) for it in d['stretch_bends'] ) self.torsion_torsions = TrackedList( TorsionTorsion(self.atoms[it[0]], self.atoms[it[1]], self.atoms[it[2]], self.atoms[it[3]], self.atoms[it[4]], type=assign_type(self.torsion_torsion_types, it[5])) for it in d['torsion_torsions'] ) self.chiral_frames = TrackedList( ChiralFrame(self.atoms[it[0]], self.atoms[it[1]], it[2]) for it in d['chiral_frames'] ) self.multipole_frames = TrackedList( MultipoleFrame(self.atoms[it[0]], *it[1:]) for it in d['multipole_frames'] ) self.adjusts = TrackedList( NonbondedException(self.atoms[it[0]], self.atoms[it[1]], assign_type(self.adjust_types, it[2])) for it in d['adjusts'] ) self.acceptors = TrackedList( AcceptorDonor(self.atoms[it[0]], self.atoms[it[1]]) for it in d['acceptors'] ) self.donors = TrackedList( AcceptorDonor(self.atoms[it[0]], self.atoms[it[1]]) for it in d['donors'] ) # Transfer the exclusions for atom, excl in zip(self.atoms, d['exclusions']): for idx in excl: atom.exclude(self.atoms[idx])
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class _StructureViewerCreator(object): """ Class responsible for creating a StructureView when indexed """ def __init__(self, struct): self.struct = struct def __getitem__(self, selection): struct = self.struct if isinstance(selection, integer_types): return struct.atoms[selection] view = StructureView() selection = struct._get_selection_array(selection) # _get_selection_array might have returned either an Atom or None if # there is no selection. Handle those and return, or continue if we got # a list of atoms if selection is None: return view elif isinstance(selection, Atom): return selection sumsel = sum(selection) if sumsel == 0: # No atoms selected. Return None return view # The cumulative sum of selection will give our index + 1 of each # selected atom into the new structure scan = [selection[0]] for i in range(1, len(selection)): scan.append(scan[i-1] + selection[i]) # Zero-out the unselected atoms scan = [x * y for x, y in zip(scan, selection)] # Copy all parameters sel_res = set() for i, atom in enumerate(struct.atoms): if not selection[i]: continue view.atoms.append(atom) if atom.residue in sel_res: continue view.residues.append(atom.residue) sel_res.add(atom.residue) def add_valence_terms(oval, sval, attrlist): """ Adds the valence terms from one list to another; oval=Other VALence; otyp=Other TYPe; sval=Self VALence; styp=Self TYPe; attrlist=ATTRibute LIST (atom1, atom2, ...) """ for val in sval: ats = [getattr(val, attr) for attr in attrlist] # Make sure all of our atoms in this valence term is "selected" indices = [scan[at.idx] for at in ats if isinstance(at, Atom)] if not all(indices): continue # Add the type if applicable oval.append(val) add_valence_terms(view.bonds, struct.bonds, ['atom1', 'atom2']) add_valence_terms(view.angles, struct.angles, ['atom1', 'atom2', 'atom3']) add_valence_terms(view.dihedrals, struct.dihedrals, ['atom1', 'atom2', 'atom3', 'atom4']) add_valence_terms(view.rb_torsions, struct.rb_torsions, ['atom1', 'atom2', 'atom3', 'atom4']) add_valence_terms(view.urey_bradleys, struct.urey_bradleys, ['atom1', 'atom2']) add_valence_terms(view.impropers, struct.impropers, ['atom1', 'atom2', 'atom3', 'atom4']) add_valence_terms(view.cmaps, struct.cmaps, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) add_valence_terms(view.trigonal_angles, struct.trigonal_angles, ['atom1', 'atom2', 'atom3', 'atom4']) add_valence_terms(view.out_of_plane_bends, struct.out_of_plane_bends, ['atom1', 'atom2', 'atom3', 'atom4']) add_valence_terms(view.pi_torsions, struct.pi_torsions, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5', 'atom6']) add_valence_terms(view.stretch_bends, struct.stretch_bends, ['atom1', 'atom2', 'atom3']) add_valence_terms(view.torsion_torsions, struct.torsion_torsions, ['atom1', 'atom2', 'atom3', 'atom4', 'atom5']) add_valence_terms(view.chiral_frames, struct.chiral_frames, ['atom1', 'atom2']) add_valence_terms(view.multipole_frames, struct.multipole_frames, ['atom']) add_valence_terms(view.adjusts, struct.adjusts, ['atom1', 'atom2']) add_valence_terms(view.donors, struct.donors, ['atom1', 'atom2']) add_valence_terms(view.acceptors, struct.acceptors, ['atom1', 'atom2']) return view #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[docs]class StructureView(object): """ A view of a Structure. In many cases, this can serve as a duck-typed Structure and it has many of the same attributes. However, none of its lists *own* their objects, and the lists of atoms, residues, and parameters/topologies are regular lists, rather than TrackedList instances. Therefore, the indexes correspond to the indexes from the *original* Structure from which this view was taken. Furthermore, there are no "type" lists, as they would be exactly equivalent to the type lists of the parent Structure instance. Attributes ---------- atoms : :class:`AtomList` List of all atoms in the structure residues : :class:`ResidueList` List of all residues in the structure bonds : :class:`TrackedList` (:class:`Bond`) List of all bonds in the structure angles : :class:`TrackedList` (:class:`Angle`) List of all angles in the structure dihedrals : :class:`TrackedList` (:class:`Dihedral`) List of all dihedrals in the structure rb_torsions : :class:`TrackedList` (:class:`RBTorsion`) List of all Ryckaert-Bellemans torsions in the structure urey_bradleys : :class:`TrackedList` (:class:`UreyBradley`) List of all Urey-Bradley angle bends in the structure impropers : :class:`TrackedList` (:class:`Improper`) List of all CHARMM-style improper torsions in the structure cmaps : :class:`TrackedList` (:class:`Cmap`) List of all CMAP objects in the structure trigonal_angles : :class:`TrackedList` (:class:`TrigonalAngle`) List of all AMOEBA-style trigonal angles in the structure out_of_plane_bends : :class:`TrackedList` (:class:`OutOfPlaneBends`) List of all AMOEBA-style out-of-plane bending angles pi_torsions : :class:`TrackedList` (:class:`PiTorsion`) List of all AMOEBA-style pi-torsion angles stretch_bends : :class:`TrackedList` (:class:`StretchBend`) List of all AMOEBA-style stretch-bend compound bond/angle terms torsion_torsions : :class:`TrackedList` (:class:`TorsionTorsion`) List of all AMOEBA-style coupled torsion-torsion terms chiral_frames : :class:`TrackedList` (:class:`ChiralFrame`) List of all AMOEBA-style chiral frames defined in the structure multipole_frames : :class:`TrackedList` (:class:`MultipoleFrame`) List of all AMOEBA-style multipole frames defined in the structure adjusts : :class:`TrackedList` (:class:`NonbondedException`) List of all nonbonded pair-exception rules acceptors : :class:`TrackedList` (:class:`AcceptorDonor`) List of all H-bond acceptors, if that information is present donors : :class:`TrackedList` (:class:`AcceptorDonor`) List of all H-bond donors, if that information is present positions : u.Quantity(list(Vec3), u.angstroms) Unit-bearing atomic coordinates. If not all atoms have coordinates, this property is None coordinates : np.ndarray of shape (nframes, natom, 3) If no coordinates are set, this is set to None. The first frame will match the coordinates present on the atoms. """ def __init__(self): self.atoms = [] self.residues = [] self.bonds = [] self.angles = [] self.dihedrals = [] self.rb_torsions = [] self.urey_bradleys = [] self.impropers = [] self.cmaps = [] self.trigonal_angles = [] self.out_of_plane_bends = [] self.pi_torsions = [] self.stretch_bends = [] self.torsion_torsions = [] self.chiral_frames = [] self.multipole_frames = [] self.adjusts = [] self.acceptors = [] self.donors = []
[docs] def to_dataframe(self): """ Generates a DataFrame from the current Structure's atomic properties Returns ------- df : DataFrame DataFrame with all atomic properties See Also -------- :func:`parmed.utils.pandautils.create_dataframe` """ from parmed.utils.pandautils import create_dataframe return create_dataframe(self)
[docs] def load_dataframe(self, df): """ Loads atomic properties from an input DataFrame Parameters ---------- df : pandas.DataFrame A pandas DataFrame with atomic properties that will be used to set the properties on the current list of atoms See Also -------- :func:`parmed.utils.pandautils.load_dataframe` """ from parmed.utils.pandautils import load_dataframe return load_dataframe(self, df)
@property def coordinates(self): try: return np.array([[a.xx, a.xy, a.xz] for a in self.atoms]) except AttributeError: return None @property def positions(self): try: return [Vec3(a.xx,a.xy,a.xz) for a in self.atoms] * u.angstroms except AttributeError: return None def __bool__(self): return bool(self.atoms or self.residues or self.bonds or self.angles or self.dihedrals or self.impropers or self.rb_torsions or self.cmaps or self.torsion_torsions or self.stretch_bends or self.out_of_plane_bends or self.trigonal_angles or self.torsion_torsions or self.pi_torsions or self.urey_bradleys or self.chiral_frames or self.multipole_frames or self.adjusts or self.acceptors or self.donors) # Ignore types here because it's just a view def __repr__(self): natom = len(self.atoms) nres = len(self.residues) nextra = sum([isinstance(a, ExtraPoint) for a in self.atoms]) retstr = ['<%s %d atoms' % (type(self).__name__, natom)] if nextra > 0: retstr.append(' [%d EPs]' % nextra) retstr.append('; %d residues' % nres) nbond = len(self.bonds) retstr.append('; %d bonds>' % nbond) return ''.join(retstr) __nonzero__ = __bool__ # For Python 2 def __iter__(self): return iter(self.atoms)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def _settler(parm): """ Identifies residues that can have SETTLE applied to it """ is_water = [] for r in parm.residues: na = sum(1 for a in r if not isinstance(a, ExtraPoint)) if na != 3: is_water.append(False) continue # Make sure we are not bonded to any other residues for a in r: for a2 in a.bond_partners: if a2.residue is not r: is_water.append(False) break else: # this atom is OK continue # We only hit here if we hit the break above. So break here break else: # Don't check elements, since they may not *always* be accurate. # It's more likely that a 3-atom residue not bonded to any other # residue is water. is_water.append(True) assert len(is_water) == len(parm.residues), 'Incorrect length' return is_water def _res_in_templlib(res, lib): """ Returns the residue template inside lib that matches res """ if res.name in lib: return lib[res.name] if len(res.name) == 3 and residue.AminoAcidResidue.has(res.name): return lib[residue.AminoAcidResidue.get(res.name).abbr] if residue.DNAResidue.has(res.name): return lib[residue.DNAResidue.get(res.name).abbr] if residue.RNAResidue.has(res.name) and residue.RNAResidue.get(res.name).abbr != 'T': return lib[residue.RNAResidue.get(res.name).abbr] # Not present return None