# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
"""
ParmEd topology parser
======================
Converts a `ParmEd <https://parmed.github.io/ParmEd/html>`_ :class:`parmed.structure.Structure` into a :class:`MDAnalysis.core.Topology`.
Example
-------
If you want to use an MDAnalysis-written ParmEd structure for simulation
in ParmEd, you need to first read your files with ParmEd to include the
necessary topology parameters. ::
>>> import parmed as pmd
>>> import MDAnalysis as mda
>>> from MDAnalysis.tests.datafiles import PRM7_ala2, RST7_ala2
>>> prm = pmd.load_file(PRM7_ala2, RST7_ala2)
>>> prm
<AmberParm 3026 atoms; 1003 residues; 3025 bonds; PBC (orthogonal); parametrized>
We can then convert this to an MDAnalysis structure, select only the
protein atoms, and then convert it back to ParmEd. ::
>>> u = mda.Universe(prm)
>>> u
<Universe with 3026 atoms>
>>> prot = u.select_atoms('protein')
>>> prm_prot = prot.convert_to('PARMED')
>>> prm_prot
<Structure 23 atoms; 2 residues; 22 bonds; PBC (orthogonal); parametrized>
From here you can create an OpenMM simulation system and minimize the energy. ::
>>> import simtk.openmm as mm
>>> import simtk.openmm.app as app
>>> from parmed import unit as u
>>> system = prm_prot.createSystem(nonbondedMethod=app.NoCutoff,
... constraints=app.HBonds,
... implicitSolvent=app.GBn2)
>>> integrator = mm.LangevinIntegrator(
... 300*u.kelvin, # Temperature of heat bath
... 1.0/u.picoseconds, # Friction coefficient
... 2.0*u.femtoseconds, # Time step
... )
>>> sim = app.Simulation(prm_prot.topology, system, integrator)
>>> sim.context.setPositions(prm_prot.positions)
>>> sim.minimizeEnergy(maxIterations=500)
Now you can continue on and run a simulation, if you wish.
Classes
-------
.. autoclass:: ParmEdParser
:members:
:inherited-members:
"""
from __future__ import absolute_import, division
from six.moves import range
from six import raise_from
import logging
import functools
from math import ceil
import numpy as np
from ..lib.util import openany
from .guessers import guess_atom_element
from .base import TopologyReaderBase, change_squash
from .tables import Z2SYMB
from ..core.topologyattrs import (
Atomids,
Atomnames,
AltLocs,
ChainIDs,
Atomtypes,
Occupancies,
Tempfactors,
Elements,
Masses,
Charges,
Resids,
Resnums,
Resnames,
Segids,
GBScreens,
SolventRadii,
NonbondedIndices,
RMins,
Epsilons,
RMin14s,
Epsilon14s,
Bonds,
UreyBradleys,
Angles,
Dihedrals,
Impropers,
CMaps
)
from ..core.topology import Topology
logger = logging.getLogger("MDAnalysis.topology.ParmEdParser")
def squash_identical(values):
if len(values) == 1:
return values[0]
else:
return tuple(values)
[docs]class ParmEdParser(TopologyReaderBase):
"""
For ParmEd structures
"""
format = 'PARMED'
@staticmethod
def _format_hint(thing):
"""Can this Parser read object *thing*?
.. versionadded:: 1.0.0
"""
try:
import parmed as pmd
except ImportError: # if no parmed, probably not parmed
return False
else:
return isinstance(thing, pmd.Structure)
[docs] def parse(self, **kwargs):
"""Parse PARMED into Topology
Returns
-------
MDAnalysis *Topology* object
"""
structure = self.filename
#### === ATOMS === ####
names = []
masses = []
charges = []
types = []
atomic_numbers = []
serials = []
resnames = []
resids = []
chainids = []
segids = []
altLocs = []
bfactors = []
occupancies = []
screens = []
solvent_radii = []
nonbonded_indices = []
rmins = []
epsilons = []
rmin14s = []
epsilon14s = []
for atom in structure.atoms:
names.append(atom.name)
masses.append(atom.mass)
charges.append(atom.charge)
types.append(atom.type)
atomic_numbers.append(atom.atomic_number)
serials.append(atom.number)
resnames.append(atom.residue.name)
resids.append(atom.residue.number)
chainids.append(atom.residue.chain)
segids.append(atom.residue.segid)
altLocs.append(atom.altloc)
bfactors.append(atom.bfactor)
occupancies.append(atom.occupancy)
screens.append(atom.screen)
solvent_radii.append(atom.solvent_radius)
nonbonded_indices.append(atom.nb_idx)
rmins.append(atom.rmin)
epsilons.append(atom.epsilon)
rmin14s.append(atom.rmin_14)
epsilon14s.append(atom.epsilon_14)
attrs = []
n_atoms = len(names)
elements = []
for z, name in zip(atomic_numbers, names):
try:
elements.append(Z2SYMB[z])
except KeyError:
elements.append(guess_atom_element(name))
# Make Atom TopologyAttrs
for vals, Attr, dtype in (
(names, Atomnames, object),
(masses, Masses, np.float32),
(charges, Charges, np.float32),
(types, Atomtypes, object),
(elements, Elements, object),
(serials, Atomids, np.int32),
(chainids, ChainIDs, object),
(altLocs, AltLocs, object),
(bfactors, Tempfactors, np.float32),
(occupancies, Occupancies, np.float32),
(screens, GBScreens, np.float32),
(solvent_radii, SolventRadii, np.float32),
(nonbonded_indices, NonbondedIndices, np.int32),
(rmins, RMins, np.float32),
(epsilons, Epsilons, np.float32),
(rmin14s, RMin14s, np.float32),
(epsilon14s, Epsilon14s, np.float32),
):
attrs.append(Attr(np.array(vals, dtype=dtype)))
resids = np.array(resids, dtype=np.int32)
resnames = np.array(resnames, dtype=object)
chainids = np.array(chainids, dtype=object)
segids = np.array(segids, dtype=object)
residx, (resids, resnames, chainids, segids) = change_squash(
(resids, resnames, chainids, segids), (resids, resnames, chainids, segids))
n_residues = len(resids)
attrs.append(Resids(resids))
attrs.append(Resnums(resids.copy()))
attrs.append(Resnames(resnames))
segidx, (segids,) = change_squash((segids,), (segids,))
n_segments = len(segids)
attrs.append(Segids(segids))
#### === OTHERS === ####
bond_values = {}
bond_types = []
bond_orders = []
ub_values = {}
ub_types = []
angle_values = {}
angle_types = []
dihedral_values = {}
dihedral_types = []
improper_values = {}
improper_types = []
cmap_values = {}
cmap_types = []
for bond in structure.bonds:
idx = (bond.atom1.idx, bond.atom2.idx)
if idx not in bond_values:
bond_values[idx] = ([bond], [bond.order])
else:
bond_values[idx][0].append(bond)
bond_values[idx][1].append(bond.order)
try:
bond_values, values = zip(*list(bond_values.items()))
except ValueError:
bond_values, bond_types, bond_orders = [], [], []
else:
bond_types, bond_orders = zip(*values)
bond_types = list(map(squash_identical, bond_types))
bond_orders = list(map(squash_identical, bond_orders))
attrs.append(Bonds(bond_values, types=bond_types, guessed=False,
order=bond_orders))
for pmdlist, na, values, types in (
(structure.urey_bradleys, 2, ub_values, ub_types),
(structure.angles, 3, angle_values, angle_types),
(structure.dihedrals, 4, dihedral_values, dihedral_types),
(structure.impropers, 4, improper_values, improper_types),
(structure.cmaps, 5, cmap_values, cmap_types),
):
for p in pmdlist:
atoms = ['atom{}'.format(i) for i in range(1, na+1)]
idx = tuple(getattr(p, a).idx for a in atoms)
if idx not in values:
values[idx] = [p]
else:
values[idx].append(p)
for dct, Attr in (
(ub_values, UreyBradleys),
(angle_values, Angles),
(dihedral_values, Dihedrals),
(improper_values, Impropers),
(cmap_values, CMaps),
):
try:
vals, types = zip(*list(dct.items()))
except ValueError:
vals, types = [], []
types = list(map(squash_identical, types))
attrs.append(Attr(vals, types=types, guessed=False, order=None))
top = Topology(n_atoms, n_residues, n_segments,
attrs=attrs,
atom_resindex=residx,
residue_segindex=segidx)
return top