Python Topology.atoms方法代码示例

# 需要导入模块: from simtk.openmm.app import Topology [as 别名]
# 或者: from simtk.openmm.app.Topology import atoms [as 别名]
def generateTopologyFromOEMol(molecule):
    """
    Generate an OpenMM Topology object from an OEMol molecule.

    Parameters
    ----------
    molecule : openeye.oechem.OEMol
        The molecule from which a Topology object is to be generated.

    Returns
    -------
    topology : simtk.openmm.app.Topology
        The Topology object generated from `molecule`.

    """
    # Create a Topology object with one Chain and one Residue.
    from simtk.openmm.app import Topology
    topology = Topology()
    chain = topology.addChain()
    resname = molecule.GetTitle()
    residue = topology.addResidue(resname, chain)

    # Create atoms in the residue.
    for atom in molecule.GetAtoms():
        name = atom.GetName()
        element = Element.getByAtomicNumber(atom.GetAtomicNum())
        atom = topology.addAtom(name, element, residue)

    # Create bonds.
    atoms = { atom.name : atom for atom in topology.atoms() }
    for bond in molecule.GetBonds():
        topology.addBond(atoms[bond.GetBgn().GetName()], atoms[bond.GetEnd().GetName()])

    return topology

# 需要导入模块: from simtk.openmm.app import Topology [as 别名]
# 或者: from simtk.openmm.app.Topology import atoms [as 别名]
    def __init__(self, file):
        """Load a prmtop file."""
        top = Topology()
        ## The Topology read from the prmtop file
        self.topology = top

        # Load the prmtop file

        prmtop = amber_file_parser.PrmtopLoader(file)
        self._prmtop = prmtop

        # Add atoms to the topology

        PDBFile._loadNameReplacementTables()
        lastResidue = None
        c = top.addChain()
        for index in range(prmtop.getNumAtoms()):
            resNumber = prmtop.getResidueNumber(index)
            if resNumber != lastResidue:
                lastResidue = resNumber
                resName = prmtop.getResidueLabel(iAtom=index).strip()
                if resName in PDBFile._residueNameReplacements:
                    resName = PDBFile._residueNameReplacements[resName]
                r = top.addResidue(resName, c)
                if resName in PDBFile._atomNameReplacements:
                    atomReplacements = PDBFile._atomNameReplacements[resName]
                else:
                    atomReplacements = {}
            atomName = prmtop.getAtomName(index).strip()
            if atomName in atomReplacements:
                atomName = atomReplacements[atomName]

            # Try to guess the element.

            upper = atomName.upper()
            if upper.startswith('CL'):
                element = elem.chlorine
            elif upper.startswith('NA'):
                element = elem.sodium
            elif upper.startswith('MG'):
                element = elem.magnesium
            else:
                try:
                    element = elem.get_by_symbol(atomName[0])
                except KeyError:
                    element = None
            top.addAtom(atomName, element, r)

        # Add bonds to the topology

        atoms = list(top.atoms())
        for bond in prmtop.getBondsWithH():
            top.addBond(atoms[bond[0]], atoms[bond[1]])
        for bond in prmtop.getBondsNoH():
            top.addBond(atoms[bond[0]], atoms[bond[1]])

        # Set the periodic box size.

        if prmtop.getIfBox():
            top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer)

# 需要导入模块: from simtk.openmm.app import Topology [as 别名]
# 或者: from simtk.openmm.app.Topology import atoms [as 别名]
class OpenMMAmberParm(AmberParm):
    """
    OpenMM-compatible subclass of AmberParm. This object should still work with
    the ParmEd API while also being compatible with OpenMM's environment
    """
   
    # Define default force groups for all of the bonded terms. This allows them
    # to be turned on and off selectively. This is a way to implement per-term
    # energy decomposition to compare individual components

    BOND_FORCE_GROUP = 0
    ANGLE_FORCE_GROUP = 1
    DIHEDRAL_FORCE_GROUP = 2
    NONBONDED_FORCE_GROUP = 3
    GB_FORCE_GROUP = 3

    def openmm_LJ(self):
        """
        Same as fill_LJ, except it uses 0.5 for the LJ radius for H-atoms with
        no vdW parameters (per OpenMM's standard)

        Returns:
            list, list : The 1st list is the list of all Rmin/2 terms. The
                         2nd is the list of all epsilon (or well depth) terms.
        """
        LJ_radius = []  # empty LJ_radii so it can be re-filled
        LJ_depth = []   # empty LJ_depths so it can be re-filled
        one_sixth = 1 / 6    # we need to raise some numbers to the 1/6th power

        ntypes = self.pointers['NTYPES']
        acoef = self.parm_data['LENNARD_JONES_ACOEF']
        bcoef = self.parm_data['LENNARD_JONES_BCOEF']

        for i in range(ntypes):
            lj_index = self.parm_data["NONBONDED_PARM_INDEX"][ntypes*i+i] - 1
            if acoef[lj_index] < 1.0e-10:
                LJ_radius.append(0.5)
                LJ_depth.append(0)
            else:
                factor = (2 * acoef[lj_index] / bcoef[lj_index])
                LJ_radius.append(pow(factor, one_sixth) * 0.5)
                LJ_depth.append(bcoef[lj_index] / 2 / factor)
      
        # Now check that we haven't modified any off-diagonals, since that will
        # not work with OpenMM
        for i in range(ntypes):
            for j in range(ntypes):
                idx = self.parm_data['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
                rij = LJ_radius[i] + LJ_radius[j]
                wdij = sqrt(LJ_depth[i] * LJ_depth[j])
                a = acoef[idx]
                b = bcoef[idx]
                if a == 0 or b == 0:
                    if a != 0 or b != 0 or (wdij != 0 and rij != 0):
                        raise OpenMMError('Off-diagonal LJ modifications '
                                          'detected. These are incompatible '
                                          'with the OpenMM API')
                elif (abs((a - (wdij * rij**12)) / a) > 1e-6 or
                      abs((b - (2 * wdij * rij**6)) / b) > 1e-6):
                    raise OpenMMError(
                            'Off-diagonal LJ modifications detected. These are '
                            'incompatible with the OpenMM API. Acoef=%s; '
                            'computed=%s. Bcoef=%s; computed=%s' %
                            (acoef, wdij*rij**12, bcoef, 2*wdij*rij**6)
                    )

        return LJ_radius, LJ_depth

    def openmm_14_LJ(self):
        """
        Returns the radii and depths for the LJ interactions between 1-4 pairs.
        For Amber topology files this is the same as the normal LJ parameters,
        but is done here so that OpenMMChamberParm can inherit and override this
        behavior without having to duplicate all of the system creation code.
        """
        return self.openmm_LJ()

    @property
    def topology(self):
        """
        The OpenMM Topology object. Cached when possible, but any changes to the
        topology object lists results in the topology being deleted and rebuilt
        """
        # If anything changed, rebuild the topology
        if not self._topology_changed():
            try:
                return self._topology
            except AttributeError:
                pass
        else:
            self.remake_parm()
      
        self._topology = Topology()

        # Add all of the atoms to the topology file in the same chain
        chain = self._topology.addChain()
        last_residue = None
        for i, atm in enumerate(self.atom_list):
            resnum = atm.residue.idx
            if last_residue != resnum:
#.........这里部分代码省略.........

# 需要导入模块: from simtk.openmm.app import Topology [as 别名]
# 或者: from simtk.openmm.app.Topology import atoms [as 别名]
    def addHydrogens(self, forcefield, pH=7.0, variants=None, platform=None):
        """Add missing hydrogens to the model.

        Some residues can exist in multiple forms depending on the pH and properties of the local environment.  These
        variants differ in the presence or absence of particular hydrogens.  In particular, the following variants
        are supported:

        Aspartic acid:
            ASH: Neutral form with a hydrogen on one of the delta oxygens
            ASP: Negatively charged form without a hydrogen on either delta oxygen

        Cysteine:
            CYS: Neutral form with a hydrogen on the sulfur
            CYX: No hydrogen on the sulfur (either negatively charged, or part of a disulfide bond)

        Glutamic acid:
            GLH: Neutral form with a hydrogen on one of the epsilon oxygens
            GLU: Negatively charged form without a hydrogen on either epsilon oxygen

        Histidine:
            HID: Neutral form with a hydrogen on the ND1 atom
            HIE: Neutral form with a hydrogen on the NE2 atom
            HIP: Positively charged form with hydrogens on both ND1 and NE2

        Lysine:
            LYN: Neutral form with two hydrogens on the zeta nitrogen
            LYS: Positively charged form with three hydrogens on the zeta nitrogen

        The variant to use for each residue is determined by the following rules:

        1. The most common variant at the specified pH is selected.
        2. Any Cysteine that participates in a disulfide bond uses the CYX variant regardless of pH.
        3. For a neutral Histidine residue, the HID or HIE variant is selected based on which one forms a better hydrogen bond.

        You can override these rules by explicitly specifying a variant for any residue.  Also keep in mind that this
        function will only add hydrogens.  It will never remove ones that are already present in the model, regardless
        of the specified pH.

        Definitions for standard amino acids and nucleotides are built in.  You can call loadHydrogenDefinitions() to load
        additional definitions for other residue types.

        Parameters:
         - forcefield (ForceField) the ForceField to use for determining the positions of hydrogens
         - pH (float=7.0) the pH based on which to select variants
         - variants (list=None) an optional list of variants to use.  If this is specified, its length must equal the number
           of residues in the model.  variants[i] is the name of the variant to use for residue i (indexed starting at 0).
           If an element is None, the standard rules will be followed to select a variant for that residue.
         - platform (Platform=None) the Platform to use when computing the hydrogen atom positions.  If this is None,
           the default Platform will be used.
        Returns: a list of what variant was actually selected for each residue, in the same format as the variants parameter
        """
        # Check the list of variants.

        residues = list(self.topology.residues())
        if variants is not None:
            if len(variants) != len(residues):
                raise ValueError("The length of the variants list must equal the number of residues")
        else:
            variants = [None]*len(residues)
        actualVariants = [None]*len(residues)

        # Load the residue specifications.

        if not Modeller._hasLoadedStandardHydrogens:
            Modeller.loadHydrogenDefinitions(os.path.join(os.path.dirname(__file__), 'data', 'hydrogens.xml'))

        # Make a list of atoms bonded to each atom.

        bonded = {}
        for atom in self.topology.atoms():
            bonded[atom] = []
        for atom1, atom2 in self.topology.bonds():
            bonded[atom1].append(atom2)
            bonded[atom2].append(atom1)

        # Define a function that decides whether a set of atoms form a hydrogen bond, using fairly tolerant criteria.

        def isHbond(d, h, a):
            if norm(d-a) > 0.35*nanometer:
                return False
            deltaDH = h-d
            deltaHA = a-h
            deltaDH /= norm(deltaDH)
            deltaHA /= norm(deltaHA)
            return acos(dot(deltaDH, deltaHA)) < 50*degree

        # Loop over residues.

        newTopology = Topology()
        newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
        newAtoms = {}
        newPositions = []*nanometer
        newIndices = []
        acceptors = [atom for atom in self.topology.atoms() if atom.element in (elem.oxygen, elem.nitrogen)]
        for chain in self.topology.chains():
            newChain = newTopology.addChain()
            for residue in chain.residues():
                newResidue = newTopology.addResidue(residue.name, newChain)
                isNTerminal = (residue == chain._residues[0])
                isCTerminal = (residue == chain._residues[-1])
#.........这里部分代码省略.........