# -*- 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 Lesser GNU Public Licence, v2.1 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
#
# TPR parser and tpr support module
# Copyright (c) 2011 Zhuyi Xue
# Released under the Lesser GNU Public Licence, v2.1+
"""
Utilities for the TPRParser
===========================
Function calling order::
(TPRParser.py call do_mtop)
do_mtop -> do_symtab
-> do_ffparams -> do_iparams
-> do_moltype -> do_atoms -> do_atom
-> do_resinfo
-> do_ilists
-> do_block
-> do_blocka
-> do_molblock
Then compose the stuffs in the format :class:`MDAnalysis.Universe` reads in.
The module also contains the :func:`do_inputrec` to read the TPR header with.
"""
import numpy as np
from mda_xdrlib import xdrlib
import struct
from . import obj
from . import setting
from ..base import squash_by
from ...core.topology import Topology
from ...core.topologyattrs import (
Atomids,
Atomnames,
Atomtypes,
Masses,
Charges,
Elements,
Resids,
Resnames,
Moltypes,
Molnums,
Segids,
ChainIDs,
Bonds,
Angles,
Dihedrals,
Impropers,
)
[docs]
class TPXUnpacker(xdrlib.Unpacker):
"""
Extend the standard XDR unpacker for the specificity of TPX files.
"""
def __init__(self, data):
super().__init__(data)
self._buf = self.get_buffer()
# The parent class uses a dunder attribute to store the
# cursor position. This property makes it easier to manipulate
# this attribute that is otherwise "protected".
@property
def _pos(self):
return self.get_position()
@_pos.setter
def _pos(self, value):
self.set_position(value)
def _unpack_value(self, item_size, struct_template):
start_position = self._pos
end_position = self._pos = start_position + item_size
content = self._buf[start_position:end_position]
if len(content) != item_size:
raise EOFError
return struct.unpack(struct_template, content)[0]
[docs]
def unpack_int64(self):
return self._unpack_value(8, '>q')
[docs]
def unpack_uint64(self):
return self._unpack_value(8, '>Q')
[docs]
def unpack_ushort(self):
return self.unpack_uint()
[docs]
def unpack_uchar(self):
# TPX files prior to gromacs 2020 (tpx version < 119) use unsigned ints
# (4 bytes) instead of unsigned chars.
return self._unpack_value(4, '>I')
[docs]
def do_string(self):
"""
Emulate gmx_fio_do_string
gmx_fio_do_string reads a string from a XDR file. On the contrary to the
python unpack_string, gmx_fio_do_string reads the size as an unsigned
integer before reading the actual string.
See <gromacs-2016-src>/src/gromacs/fileio/gmx_system_xdr.c:454
"""
self.unpack_int()
return self.unpack_string()
[docs]
class TPXUnpacker2020(TPXUnpacker):
"""
Unpacker for TPX files from and later than gromacs 2020.
A new implementation of the serializer (InMemorySerializer), introduced in
gromacs 2020, changes le meaning of some types in the file body (the header
keep using the previous implementation of the serializer).
"""
[docs]
@classmethod
def from_unpacker(cls, unpacker):
new_unpacker = cls(unpacker.get_buffer())
new_unpacker._pos = unpacker.get_position()
if hasattr(unpacker, 'unpack_real'):
if unpacker.unpack_real == unpacker.unpack_float:
new_unpacker.unpack_real = new_unpacker.unpack_float
elif unpacker.unpack_real == unpacker.unpack_double:
new_unpacker.unpack_real = new_unpacker.unpack_double
else:
raise ValueError("Unrecognized precision")
return new_unpacker
[docs]
def unpack_fstring(self, n):
if n < 0:
raise ValueError('Size of fstring cannot be negative.')
start_position = self._pos
end_position = self._pos = start_position + n
if end_position > len(self._buf):
raise EOFError
content = self._buf[start_position:end_position]
return content
[docs]
def unpack_ushort(self):
# The InMemorySerializer implements ushort according to the XDR standard
# on the contrary to the IO serializer.
return self._unpack_value(2, '>H')
[docs]
def unpack_uchar(self):
# The InMemorySerializer implements uchar according to the XDR standard
# on the contrary to the IO serializer.
return self._unpack_value(1, '>B')
[docs]
def do_string(self):
"""
Emulate gmx_fio_do_string
"""
n = self.unpack_uint64()
return self.unpack_fstring(n)
[docs]
def ndo_int(data, n):
"""mimic of gmx_fio_ndo_real in gromacs"""
return [data.unpack_int() for i in range(n)]
[docs]
def ndo_real(data, n):
"""mimic of gmx_fio_ndo_real in gromacs"""
return [data.unpack_real() for i in range(n)]
[docs]
def do_rvec(data):
return data.unpack_farray(setting.DIM, data.unpack_real)
[docs]
def ndo_rvec(data, n):
"""mimic of gmx_fio_ndo_rvec in gromacs"""
return [data.unpack_farray(setting.DIM, data.unpack_real) for i in range(n)]
[docs]
def ndo_ivec(data, n):
"""mimic of gmx_fio_ndo_rvec in gromacs"""
return [data.unpack_farray(setting.DIM, data.unpack_int) for i in range(n)]
[docs]
def fileVersion_err(fver):
if fver not in setting.SUPPORTED_VERSIONS:
raise NotImplementedError(
f"Your tpx version is {fver}, which this parser does not support, yet "
)
[docs]
def define_unpack_real(prec, data):
"""Define an unpack_real method of data based on the float precision used"""
if prec == 4:
data.unpack_real = data.unpack_float
elif prec == 8:
data.unpack_real = data.unpack_double
else:
raise ValueError(f"unsupported precision: {prec}")
[docs]
def do_mtop(data, fver, tpr_resid_from_one=False, precision=4):
# mtop: the topology of the whole system
symtab = do_symtab(data)
do_symstr(data, symtab) # system_name
ff_params = do_ffparams(data, fver) # params
nmoltype = data.unpack_int()
moltypes = [] # non-gromacs
for i in range(nmoltype):
moltype = do_moltype(data, symtab, fver)
moltypes.append(moltype)
nmolblock = data.unpack_int()
mtop = obj.Mtop(nmoltype, moltypes, nmolblock)
bonds = []
angles = []
dihedrals = []
impropers = []
atomids = []
segids = []
chainIDs = []
resids = []
resnames = []
atomnames = []
atomtypes = []
moltypes = []
molnums = []
charges = []
masses = []
elements = []
atom_start_ndx = 0
res_start_ndx = 0
molnum = 0
for i in range(mtop.nmolblock):
# molb_type is just an index for moltypes/molecule_types
mb = do_molblock(data)
# segment is made to correspond to the molblock as in gromacs, the
# naming is kind of arbitrary
molblock = mtop.moltypes[mb.molb_type].name.decode('utf-8')
segid = f"seg_{i}_{molblock}"
chainID = molblock[14:] if molblock[:14] == "Protein_chain_" else molblock
for j in range(mb.molb_nmol):
mt = mtop.moltypes[mb.molb_type] # mt: molecule type
for atomkind in mt.atomkinds:
atomids.append(atomkind.id + atom_start_ndx)
segids.append(segid)
chainIDs.append(chainID)
resids.append(atomkind.resid + res_start_ndx)
resnames.append(atomkind.resname.decode())
atomnames.append(atomkind.name.decode())
atomtypes.append(atomkind.type.decode())
moltypes.append(molblock)
molnums.append(molnum)
charges.append(atomkind.charge)
masses.append(atomkind.mass)
elements.append(atomkind.element_symbol)
molnum += 1
# remap_ method returns [blah, blah, ..] or []
bonds.extend(mt.remap_bonds(atom_start_ndx))
angles.extend(mt.remap_angles(atom_start_ndx))
dihedrals.extend(mt.remap_dihe(atom_start_ndx))
impropers.extend(mt.remap_impr(atom_start_ndx))
atom_start_ndx += mt.number_of_atoms()
res_start_ndx += mt.number_of_residues()
atomids = Atomids(np.array(atomids, dtype=np.int32))
atomnames = Atomnames(np.array(atomnames, dtype=object))
atomtypes = Atomtypes(np.array(atomtypes, dtype=object))
charges = Charges(np.array(charges, dtype=np.float32))
masses = Masses(np.array(masses, dtype=np.float32))
moltypes = np.array(moltypes, dtype=object)
molnums = np.array(molnums, dtype=np.int32)
segids = np.array(segids, dtype=object)
chainIDs = np.array(chainIDs, dtype=object)
resids = np.array(resids, dtype=np.int32)
if tpr_resid_from_one:
resids += 1
resnames = np.array(resnames, dtype=object)
(residx, new_resids,
(new_resnames,
new_moltypes,
new_molnums,
perres_segids
)
) = squash_by(resids,
resnames,
moltypes,
molnums,
segids)
residueids = Resids(new_resids)
residuenames = Resnames(new_resnames)
residue_moltypes = Moltypes(new_moltypes)
residue_molnums = Molnums(new_molnums)
segidx, perseg_segids = squash_by(perres_segids)[:2]
segids = Segids(perseg_segids)
chainIDs = ChainIDs(chainIDs)
top = Topology(
len(atomids),
len(new_resids),
len(perseg_segids),
attrs=[
atomids,
atomnames,
atomtypes,
charges,
masses,
residueids,
residuenames,
residue_moltypes,
residue_molnums,
segids,
chainIDs,
],
atom_resindex=residx,
residue_segindex=segidx,
)
top.add_TopologyAttr(Bonds([bond for bond in bonds if bond]))
top.add_TopologyAttr(Angles([angle for angle in angles if angle]))
top.add_TopologyAttr(Dihedrals([dihedral for dihedral in dihedrals
if dihedral]))
top.add_TopologyAttr(Impropers([improper for improper in impropers
if improper]))
if any(elements):
elements = Elements(np.array(elements, dtype=object))
top.add_TopologyAttr(elements)
# NOTE: the tpr striding code below serves the
# purpose of placing us in a suitable "seek" position
# in the binary file such that we are well placed
# for reading coords and velocities if they are present
# the order of operations is based on an analysis of
# the C++ code in do_mtop() function in the GROMACS
# source at:
# src/gromacs/fileio/tpxio.cpp
# TODO: expand tpx version support for striding to
# the coordinates
atnr = ff_params.atnr
if fver >= 58:
# TODO: the following value is important, and not sure
# how to access programmatically yet...
# from GMX source code:
# api/legacy/include/gromacs/topology/topology_enums.h
# worst case scenario we hard code it based on
# tpx/GMX version?
# for details, see gh-4873
SimulationAtomGroupType_size = 10
n_atoms = data.unpack_int()
if fver < 116:
for i in range(3 * atnr * int(precision / 4)):
data.unpack_int()
if fver < 129:
# NOTE: speculative, Tyler did this by
# inspecting binary data and relative file offsets
# for older tpr files relative to more recent;
# the stride skip appears related to `atnr` in GMX source, the
# number of non-bonded atom types
for i in range(atnr + 1):
data.unpack_int()
if 58 < fver <= 83:
for i in range(atnr * int(precision / 4) * 2):
data.unpack_int()
if fver > 83:
interm = data.unpack_uchar()
if 83 < fver < 116:
for i in range(2 * atnr):
data.unpack_int()
if fver > 58:
ngrid = data.unpack_int()
grid_spacing = data.unpack_int()
n_elements = grid_spacing ** 2
for i in range(ngrid):
for j in range(n_elements):
ndo_real(data, 4)
for i in range(SimulationAtomGroupType_size):
group_size = data.unpack_int()
ndo_int(data, group_size)
n_group_names = data.unpack_int()
for i in range(n_group_names):
data.unpack_int()
for i in range(SimulationAtomGroupType_size):
n_grp_numbers = data.unpack_int()
if n_grp_numbers != 0:
for i in range(n_grp_numbers):
data.unpack_uchar()
if fver > 119:
im_excl_grp_size = data.unpack_int()
ndo_int(data, im_excl_grp_size)
# TODO: why is this needed?
# NOTE: this `mystery_n` value is at least sometimes used
# for tpx >= 137, but is often 0 otherwise
mystery_n = data.unpack_int()
if fver >= 137:
# NOTE: Tyler observed that the first float32
# positional coordinate was offset by `mystery_n`
# 32-bit units in a sample tpx 137 tpr file
# by empirical testing
for i in range(mystery_n):
data.unpack_int()
return top
[docs]
def do_symstr(data, symtab):
# do_symstr: get a string based on index from the symtab
ndx = data.unpack_int()
return symtab[ndx]
[docs]
def do_symtab(data):
symtab_nr = data.unpack_int() # number of symbols
symtab = []
for i in range(symtab_nr):
j = data.do_string()
symtab.append(j)
return symtab
[docs]
def do_ffparams(data, fver):
atnr = data.unpack_int()
ntypes = data.unpack_int()
functype = ndo_int(data, ntypes)
reppow = data.unpack_double() if fver >= 66 else 12.0
fudgeQQ = data.unpack_real()
# mimicing the c code,
# remapping the functype due to inconsistency in different versions
for i in range(len(functype)):
for k in setting.ftupd:
# j[0]: tpx_version, j[1] funtype
if fver < k[0] and functype[i] >= k[1]:
functype[i] += 1
do_iparams(data, functype, fver)
params = obj.Params(atnr, ntypes, functype, reppow, fudgeQQ)
return params
[docs]
def do_harm(data):
data.unpack_real() # rA
data.unpack_real() # krA
data.unpack_real() # rB
data.unpack_real() # krB
[docs]
def do_iparams(data, functypes, fver):
# Not all elif cases in this function has been used and tested
for k, i in enumerate(functypes):
if i in [
setting.F_ANGLES, setting.F_G96ANGLES,
setting.F_BONDS, setting.F_G96BONDS,
setting.F_HARMONIC, setting.F_IDIHS
]:
do_harm(data)
elif i in [setting.F_RESTRANGLES]:
data.unpack_real() # harmonic.rA
data.unpack_real() # harmonic.krA
if fver >= 134:
data.unpack_real() # harmonic.rB
data.unpack_real() # harmonic.krB
elif i in [setting.F_LINEAR_ANGLES]:
data.unpack_real() # linangle.klinA
data.unpack_real() # linangle.aA
data.unpack_real() # linangle.klinB
data.unpack_real() # linangle.aB);
elif i in [setting.F_FENEBONDS]:
data.unpack_real() # fene.bm
data.unpack_real() # fene.kb
elif i in [setting.F_RESTRBONDS]:
data.unpack_real() # restraint.lowA
data.unpack_real() # restraint.up1A
data.unpack_real() # restraint.up2A
data.unpack_real() # restraint.kA
data.unpack_real() # restraint.lowB
data.unpack_real() # restraint.up1B
data.unpack_real() # restraint.up2B
data.unpack_real() # restraint.kB
elif i in [
setting.F_TABBONDS, setting.F_TABBONDSNC,
setting.F_TABANGLES, setting.F_TABDIHS,
]:
data.unpack_real() # tab.kA
data.unpack_int() # tab.table
data.unpack_real() # tab.kB
elif i in [setting.F_CROSS_BOND_BONDS]:
data.unpack_real() # cross_bb.r1e
data.unpack_real() # cross_bb.r2e
data.unpack_real() # cross_bb.krr
elif i in [setting.F_CROSS_BOND_ANGLES]:
data.unpack_real() # cross_ba.r1e
data.unpack_real() # cross_ba.r2e
data.unpack_real() # cross_ba.r3e
data.unpack_real() # cross_ba.krt
elif i in [setting.F_UREY_BRADLEY]:
data.unpack_real() # u_b.theta
data.unpack_real() # u_b.ktheta
data.unpack_real() # u_b.r13
data.unpack_real() # u_b.kUB
if fver >= 79:
data.unpack_real() # u_b.thetaB
data.unpack_real() # u_b.kthetaB
data.unpack_real() # u_b.r13B
data.unpack_real() # u_b.kUBB
elif i in [setting.F_QUARTIC_ANGLES]:
data.unpack_real() # qangle.theta
ndo_real(data, 5) # qangle.c
elif i in [setting.F_BHAM]:
data.unpack_real() # bham.a
data.unpack_real() # bham.b
data.unpack_real() # bham.c
elif i in [setting.F_MORSE]:
data.unpack_real() # morse.b0
data.unpack_real() # morse.cb
data.unpack_real() # morse.beta
if fver >= 79:
data.unpack_real() # morse.b0B
data.unpack_real() # morse.cbB
data.unpack_real() # morse.betaB
elif i in [setting.F_CUBICBONDS]:
data.unpack_real() # cubic.b0g
data.unpack_real() # cubic.kb
data.unpack_real() # cubic.kcub
elif i in [setting.F_CONNBONDS]:
pass
elif i in [setting.F_POLARIZATION]:
data.unpack_real() # polarize.alpha
elif i in [setting.F_ANHARM_POL]:
data.unpack_real() # anharm_polarize.alpha
data.unpack_real() # anharm_polarize.drcut
data.unpack_real() # anharm_polarize.khyp
elif i in [setting.F_WATER_POL]:
data.unpack_real() # wpol.al_x
data.unpack_real() # wpol.al_y
data.unpack_real() # wpol.al_z
data.unpack_real() # wpol.rOH
data.unpack_real() # wpol.rHH
data.unpack_real() # wpol.rOD
elif i in [setting.F_THOLE_POL]:
data.unpack_real() # thole.a
data.unpack_real() # thole.alpha1
data.unpack_real() # thole.alpha2
if fver < setting.tpxv_RemoveTholeRfac:
data.unpack_real() # thole.rfac
elif i in [setting.F_LJ]:
data.unpack_real() # lj_c6
data.unpack_real() # lj_c9
elif i in [setting.F_LJ14]:
data.unpack_real() # lj14_c6A
data.unpack_real() # lj14_c12A
data.unpack_real() # lj14_c6B
data.unpack_real() # lj14_c12B
elif i in [setting.F_LJC14_Q]:
data.unpack_real() # ljc14.fqq
data.unpack_real() # ljc14.qi
data.unpack_real() # ljc14.qj
data.unpack_real() # ljc14.c6
data.unpack_real() # ljc14.c12
elif i in [setting.F_LJC_PAIRS_NB]:
data.unpack_real() # ljcnb.qi
data.unpack_real() # ljcnb.qj
data.unpack_real() # ljcnb.c6
data.unpack_real() # ljcnb.c12
elif i in [
setting.F_PIDIHS, setting.F_ANGRES,
setting.F_ANGRESZ, setting.F_PDIHS,
]:
data.unpack_real() # pdihs_phiA
data.unpack_real() # pdihs_cpA
data.unpack_real() # pdihs_phiB
data.unpack_real() # pdihs_cpB
data.unpack_int() # pdihs_mult
elif i in [setting.F_RESTRDIHS]:
data.unpack_real() # pdihs.phiA
data.unpack_real() # pdihs.cpA
if fver >= 134:
data.unpack_real() # pdihs.phiB
data.unpack_real() # pdihs.cpB
elif i in [setting.F_DISRES]:
data.unpack_int() # disres.label
data.unpack_int() # disres.type
data.unpack_real() # disres.low
data.unpack_real() # disres.up1
data.unpack_real() # disres.up2
data.unpack_real() # disres.kfac
elif i in [setting.F_ORIRES]:
data.unpack_int() # orires.ex
data.unpack_int() # orires.label
data.unpack_int() # orires.power
data.unpack_real() # orires.c
data.unpack_real() # orires.obs
data.unpack_real() # orires.kfac
elif i in [setting.F_DIHRES]:
if fver < 72:
data.unpack_int() # idum
data.unpack_int() # idum
data.unpack_real() # dihres.phiA
data.unpack_real() # dihres.dphiA
data.unpack_real() # dihres.kfacA
if fver >= 72:
data.unpack_real() # dihres.phiB
data.unpack_real() # dihres.dphiB
data.unpack_real() # dihres.kfacB
elif i in [setting.F_POSRES]:
do_rvec(data) # posres.pos0A
do_rvec(data) # posres.fcA
do_rvec(data) # posres.pos0B
do_rvec(data) # posres.fcB
elif i in [setting.F_FBPOSRES]:
data.unpack_int() # fbposres.geom
do_rvec(data) # fbposres.pos0
data.unpack_real() # fbposres.r
data.unpack_real() # fbposres.k
elif i in [setting.F_CBTDIHS]:
ndo_real(data, setting.NR_CBTDIHS) # cbtdihs.cbtcA
if fver >= 134:
ndo_real(data, setting.NR_CBTDIHS) # cbtdihs.cbtcB
elif i in [setting.F_RBDIHS]:
ndo_real(data, setting.NR_RBDIHS) # iparams_rbdihs_rbcA
ndo_real(data, setting.NR_RBDIHS) # iparams_rbdihs_rbcB
elif i in [setting.F_FOURDIHS]:
# Fourier dihedrals
ndo_real(data, setting.NR_RBDIHS) # rbdihs.rbcA
ndo_real(data, setting.NR_RBDIHS) # rbdihs.rbcB
elif i in [setting.F_CONSTR, setting.F_CONSTRNC]:
data.unpack_real() # dA
data.unpack_real() # dB
elif i in [setting.F_SETTLE]:
data.unpack_real() # settle.doh
data.unpack_real() # settle.dhh
elif i in [setting.F_VSITE1]:
pass
elif i in [setting.F_VSITE2, setting.F_VSITE2FD]:
data.unpack_real() # vsite.a
elif i in [setting.F_VSITE3, setting.F_VSITE3FD, setting.F_VSITE3FAD]:
data.unpack_real() # vsite.a
data.unpack_real() # vsite.b
elif i in [setting.F_VSITE3OUT, setting.F_VSITE4FD, setting.F_VSITE4FDN]:
data.unpack_real() # vsite.a
data.unpack_real() # vsite.b
data.unpack_real() # vsite.c
elif i in [setting.F_VSITEN]:
data.unpack_int() # vsiten.n
data.unpack_real() # vsiten.a
elif i in [setting.F_GB12, setting.F_GB13, setting.F_GB14]:
# /* We got rid of some parameters in version 68 */
if fver < 68:
data.unpack_real() # rdum
data.unpack_real() # rdum
data.unpack_real() # rdum
data.unpack_real() # rdum
data.unpack_real() # gb.sar
data.unpack_real() # gb.st
data.unpack_real() # gb.pi
data.unpack_real() # gb.gbr
data.unpack_real() # gb.bmlt
elif i in [setting.F_CMAP]:
data.unpack_int() # cmap.cmapA
data.unpack_int() # cmap.cmapB
elif i in [setting.F_ENNPOT]:
# TODO: handle the new neural network potentials
# supported from GROMACS 2025.0
pass
else:
raise NotImplementedError(f"unknown functype: {i}")
return
[docs]
def do_moltype(data, symtab, fver):
molname = do_symstr(data, symtab)
# key info: about atoms
atoms_obj = do_atoms(data, symtab, fver)
#### start: MDAnalysis specific
atomkinds = []
for k, a in enumerate(atoms_obj.atoms):
atomkinds.append(obj.AtomKind(
k,
atoms_obj.atomnames[k],
atoms_obj.type[k],
a.resind,
atoms_obj.resnames[a.resind],
a.m,
a.q,
a.atomnumber,
))
#### end: MDAnalysis specific
# key info: about bonds, angles, dih, improp dih.
ilists = do_ilists(data, fver)
#### start: MDAnalysis specific
# these may not be available for certain molecules, e.g. tip4p
bonds, angles, dihs, impr = [], [], [], []
for ilist in ilists:
if ilist.nr > 0:
ik_obj = obj.InteractionKind(*ilist.ik)
ias = ilist.iatoms
# the following if..elif..else statement needs to be updated as new
# type of interactions become interested
if ik_obj.name in ['BONDS', 'G96BONDS', 'MORSE', 'CUBICBONDS',
'CONNBONDS', 'HARMONIC', 'FENEBONDS',
'RESTRAINTPOT', 'CONSTR', 'CONSTRNC',
'TABBONDS', 'TABBONDSNC']:
bonds += list(ik_obj.process(ias))
elif ik_obj.name in ['ANGLES', 'G96ANGLES', 'CROSS_BOND_BOND',
'CROSS_BOND_ANGLE', 'UREY_BRADLEY', 'QANGLES',
'RESTRANGLES', 'TABANGLES']:
angles += list(ik_obj.process(ias))
elif ik_obj.name in ['PDIHS', 'RBDIHS', 'RESTRDIHS', 'CBTDIHS',
'FOURDIHS', 'TABDIHS']:
dihs += list(ik_obj.process(ias))
elif ik_obj.name in ['IDIHS', 'PIDIHS']:
impr += list(ik_obj.process(ias))
elif ik_obj.name == 'SETTLE':
# SETTLE interactions are optimized triangular constraints for
# water molecules. They should be counted as a pair of bonds
# between the oxygen and the hydrogens. In older versions of
# the TPR format only specifies the index of the oxygen and
# assumes that the next two atoms are the hydrogens.
if len(ias) == 2:
# Old format. Only the first atom is specified.
base_atom = ias[1]
bonds += [
[base_atom, base_atom + 1],
[base_atom, base_atom + 2],
]
else:
all_settle = ik_obj.process(ias)
for settle in all_settle:
base_atom = settle[0]
bonds += [
[settle[0], settle[1]],
[settle[0], settle[2]],
]
else:
# other interaction types are not interested at the moment
pass
bonds = bonds if bonds else None
angles = angles if angles else None
dihs = dihs if dihs else None
impr = impr if impr else None
moltype = obj.MoleculeKind(molname, atomkinds, bonds, angles, dihs, impr)
#### end: MDAnalysis specific
# info in do_block and do_blocka is not interesting, but has to be parsed
# here so that all moltypes can be parsed properly
do_block(data)
do_blocka(data)
return moltype
[docs]
def do_atoms(data, symtab, fver):
nr = data.unpack_int() # number of atoms in a particular molecule
nres = data.unpack_int() # number of residues in a particular molecule
atoms = []
for i in range(nr):
A = do_atom(data, fver)
atoms.append(A)
# do_strstr
atomnames = [symtab[i] for i in ndo_int(data, nr)]
type = [symtab[i] for i in ndo_int(data, nr)] # e.g. opls_111
typeB = [symtab[i] for i in ndo_int(data, nr)]
resnames = do_resinfo(data, symtab, fver, nres)
return obj.Atoms(atoms, nr, nres, type, typeB, atomnames, resnames)
[docs]
def do_resinfo(data, symtab, fver, nres):
if fver < 63:
resnames = [symtab[i] for i in ndo_int(data, nres)]
else:
resnames = []
for i in range(nres):
resnames.append(symtab[data.unpack_int()])
data.unpack_int()
data.unpack_uchar()
return resnames
[docs]
def do_atom(data, fver):
m = data.unpack_real() # mass
q = data.unpack_real() # charge
mB = data.unpack_real()
qB = data.unpack_real()
tp = data.unpack_ushort() # type is a keyword in python
typeB = data.unpack_ushort()
ptype = data.unpack_int() # regular atom, virtual site or others
resind = data.unpack_int() # index of residue
atomnumber = data.unpack_int() # index of atom type
return obj.Atom(m, q, mB, qB, tp, typeB, ptype, resind, atomnumber)
[docs]
def do_ilists(data, fver):
nr = [] # number of ilist
iatoms = [] # atoms involved in a particular interaction type
pos = []
for j in range(setting.F_NRE): # total number of energies (i.e. interaction types)
bClear = False
for k in setting.ftupd:
if fver < k[0] and j == k[1]:
bClear = True
if bClear:
nr.append(0)
iatoms.append(None)
else:
# do_ilist
n = data.unpack_int()
nr.append(n)
l_ = []
for i in range(n):
l_.append(data.unpack_int())
iatoms.append(l_)
return [
obj.Ilist(n, it, i)
for n, it, i in zip(nr, setting.interaction_types, iatoms)
]
[docs]
def do_molblock(data):
molb_type = data.unpack_int()
molb_nmol = data.unpack_int() # number of moles in the molblock
molb_natoms_mol = data.unpack_int() # number of atoms in a molecule
molb_nposres_xA = data.unpack_int()
if molb_nposres_xA > 0:
ndo_rvec(data, molb_nposres_xA)
molb_nposres_xB = data.unpack_int() # The number of posres coords for top B
if molb_nposres_xB > 0:
ndo_rvec(data, molb_nposres_xB)
return obj.Molblock(molb_type, molb_nmol, molb_natoms_mol,
molb_nposres_xA, molb_nposres_xB)
[docs]
def do_block(data):
block_nr = data.unpack_int() # for cgs: charge groups
# starting or ending atom indices, based on which cgs are grouped
ndo_int(data, block_nr + 1)
return do_block
[docs]
def do_blocka(data):
block_nr = data.unpack_int() # No. of atoms with excls
block_nra = data.unpack_int() # total times fo appearance of atoms for excls
ndo_int(data, block_nr + 1)
ndo_int(data, block_nra)
return block_nr, block_nra