# -*- 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
#
# TPR parser and tpr support module
# Copyright (c) 2011 Zhuyi Xue
# Released under the GNU Public Licence, v2
"""
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.
"""
from __future__ import absolute_import
from six.moves import range
import warnings
import numpy as np
import xdrlib
import struct
from . import obj
from . import setting as S
from ..base import squash_by
from ...core.topology import Topology
from ...core.topologyattrs import (
Atomids,
Atomnames,
Atomtypes,
Masses,
Charges,
Resids,
Resnames,
Moltypes,
Molnums,
Segids,
Bonds,
Angles,
Dihedrals,
Impropers
)
[docs]class TPXUnpacker(xdrlib.Unpacker):
"""
Extend the standard XDR unpacker for the specificity of TPX files.
"""
def __init__(self, data):
# Using super here would be ideal, but it does not work on python2 as
# TPXUnpacker is a `classobj` and not a `type`.
# super().__init__(data)
xdrlib.Unpacker.__init__(self, 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".
# The use of this property works well in python3, but fails in python2
# where direct access to the mangled attribute seems required.
@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):
# Ideally, we should use the _pos attribute, but it present some
# unexpected behaviour on python2 where the position of the cursor is
# not kept in sync between the method defined in TPXUnpacker and the
# methods defined in the base class.
# start_position = self._pos
# end_position = self._pos = start_position + item_size
start_position = self._Unpacker__pos # pylint: disable=access-member-before-definition
end_position = self._pos = self._Unpacker__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._buf)
new_unpacker._pos = new_unpacker._Unpacker__pos = unpacker._Unpacker__pos
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._Unpacker__pos # pylint: disable=access-member-before-definition
end_position = self._pos = self._Unpacker__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(S.DIM, data.unpack_real)
[docs]def ndo_rvec(data, n):
"""mimic of gmx_fio_ndo_rvec in gromacs"""
return [data.unpack_farray(S.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(S.DIM, data.unpack_int) for i in range(n)]
[docs]def fileVersion_err(fver):
if fver not in S.SUPPORTED_VERSIONS:
raise NotImplementedError(
"Your tpx version is {0}, which this parser does not support, yet ".format(
fver))
[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("unsupported precision: {0}".format(prec))
[docs]def do_mtop(data, fver, tpr_resid_from_one=False):
# mtop: the topology of the whole system
symtab = do_symtab(data)
do_symstr(data, symtab) # system_name
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 = []
resids = []
resnames = []
atomnames = []
atomtypes = []
moltypes = []
molnums = []
charges = []
masses = []
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 = "seg_{0}_{1}".format(i, 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)
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)
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()
# not useful here
# data.unpack_int() # mtop_natoms
# do_atomtypes(data)
# mtop_ffparams_cmap_grid_ngrid = 0
# mtop_ffparams_cmap_grid_grid_spacing = 0.1
# mtop_ffparams_cmap_grid_cmapdata = 'NULL'
# do_groups(data, symtab)
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)
resids = np.array(resids, dtype=np.int32)
if tpr_resid_from_one:
resids += 1
else:
warnings.warn("TPR files index residues from 0. "
"From MDAnalysis version 2.0, resids will start "
"at 1 instead. If you wish to keep indexing "
"resids from 0, please set "
"`tpr_resid_from_one=False` as a keyword argument "
"when you create a new Topology or Universe.",
category=DeprecationWarning)
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)
top = Topology(len(atomids), len(new_resids), len(perseg_segids),
attrs=[atomids, atomnames, atomtypes,
charges, masses,
residueids, residuenames,
residue_moltypes, residue_molnums,
segids],
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]))
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()
if fver < 57:
data.unpack_int() # idum
ntypes = data.unpack_int()
functype = ndo_int(data, ntypes)
reppow = data.unpack_double() if fver >= 66 else 12.0
if fver >= 57:
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 S.ftupd:
# j[0]: tpx_version, j[1] funtype
if fver < k[0] and functype[i] >= k[1]:
functype[i] += 1
# parameters for different functions, None returned for now since not sure
# what is iparams
iparams = do_iparams(data, functype, fver)
params = obj.Params(atnr, ntypes, functype, reppow, fudgeQQ, iparams)
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 [
S.F_ANGLES, S.F_G96ANGLES,
S.F_BONDS, S.F_G96BONDS,
S.F_HARMONIC, S.F_IDIHS
]:
do_harm(data)
elif i in [S.F_RESTRANGLES]:
data.unpack_real() # harmonic.rA
data.unpack_real() # harmonic.krA
elif i in [S.F_LINEAR_ANGLES]:
data.unpack_real() # linangle.klinA
data.unpack_real() # linangle.aA
data.unpack() # linangle.klinB
data.unpack() # linangle.aB);
elif i in [S.F_FENEBONDS]:
data.unpack_real() # fene.bm
data.unpack_real() # fene.kb
elif i in [S.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 [S.F_TABBONDS, S.F_TABBONDSNC, S.F_TABANGLES, S.F_TABDIHS]:
data.unpack_real() # tab.kA
data.unpack_int() # tab.table
data.unpack_real() # tab.kB
elif i in [S.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 [S.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 [S.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 [S.F_QUARTIC_ANGLES]:
data.unpack_real() # qangle.theta
ndo_real(data, 5) # qangle.c
elif i in [S.F_BHAM]:
data.unpack_real() # bham.a
data.unpack_real() # bham.b
data.unpack_real() # bham.c
elif i in [S.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 [S.F_CUBICBONDS]:
data.unpack_real() # cubic.b0
data.unpack_real() # cubic.kb
data.unpack_real() # cubic.kcub
elif i in [S.F_CONNBONDS]:
pass
elif i in [S.F_POLARIZATION]:
data.unpack_real() # polarize.alpha
elif i in [S.F_ANHARM_POL]:
data.unpack_real() # anharm_polarize.alpha
data.unpack_real() # anharm_polarize.drcut
data.unpack_real() # anharm_polarize.khyp
elif i in [S.F_WATER_POL]:
if fver < 31:
fver_err(fver)
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 [S.F_THOLE_POL]:
data.unpack_real() # thole.a
data.unpack_real() # thole.alpha1
data.unpack_real() # thole.alpha2
data.unpack_real() # thole.rfac
elif i in [S.F_LJ]:
data.unpack_real() # lj_c6
data.unpack_real() # lj_c9
elif i in [S.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 [S.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 [S.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 [S.F_PIDIHS, S.F_ANGRES, S.F_ANGRESZ, S.F_PDIHS]:
data.unpack_real() # pdihs_phiA
data.unpack_real() # pdihs_cpA
if (i == S.F_ANGRES or i == S.F_ANGRESZ) and fver < 42:
data.unpack_real() # harmonic.rB
data.unpack_real() # harmonic.krB
else:
data.unpack_real() # pdihs_phiB
data.unpack_real() # pdihs_cpB
data.unpack_int() # pdihs_mult
elif i in [S.F_RESTRDIHS]:
data.unpack_real() # pdihs.phiA
data.unpack_real() # pdihs.cpA
elif i in [S.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 [S.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 [S.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 [S.F_POSRES]:
do_rvec(data) # posres.pos0A
do_rvec(data) # posres.fcA
if fver < 27:
fver_err(fver)
else:
do_rvec(data) # posres.pos0B
do_rvec(data) # posres.fcB
elif i in [S.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 [S.F_CBTDIHS]:
ndo_real(data, S.NR_CBTDIHS) # cbtdihs.cbtcA
elif i in [S.F_RBDIHS]:
ndo_real(data, S.NR_RBDIHS) # iparams_rbdihs_rbcA
if fver >= 25:
ndo_real(data, S.NR_RBDIHS) # iparams_rbdihs_rbcB
elif i in [S.F_FOURDIHS]:
# Fourier dihedrals
ndo_real(data, S.NR_RBDIHS) # rbdihs.rbcA
ndo_real(data, S.NR_RBDIHS) # rbdihs.rbcB
elif i in [S.F_CONSTR, S.F_CONSTRNC]:
data.unpack_real() # dA
data.unpack_real() # dB
elif i in [S.F_SETTLE]:
data.unpack_real() # settle.doh
data.unpack_real() # settle.dhh
elif i in [S.F_VSITE2, S.F_VSITE2FD]:
data.unpack_real() # vsite.a
elif i in [S.F_VSITE3, S.F_VSITE3FD, S.F_VSITE3FAD]:
data.unpack_real() # vsite.a
data.unpack_real() # vsite.b
elif i in [S.F_VSITE3OUT, S.F_VSITE4FD, S.F_VSITE4FDN]:
data.unpack_real() # vsite.a
data.unpack_real() # vsite.b
data.unpack_real() # vsite.c
elif i in [S.F_VSITEN]:
data.unpack_int() # vsiten.n
data.unpack_real() # vsiten.a
elif i in [S.F_GB12, S.F_GB13, S.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 [S.F_CMAP]:
data.unpack_int() # cmap.cmapA
data.unpack_int() # cmap.cmapB
else:
raise NotImplementedError("unknown functype: {0}".format(i))
return
[docs]def do_moltype(data, symtab, fver):
if fver >= 57:
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))
#### 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 oxigen 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 interested, 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
if fver < 57:
fver_err(fver)
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)]
if fver <= 20:
fver_err(fver)
else:
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)
if fver < 57:
fver_err(fver)
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
if fver >= 52:
atomnumber = data.unpack_int() # index of atom type
if fver < 23 or fver < 39 or fver < 57:
fver_err(fver)
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(S.F_NRE): # total number of energies (i.e. interaction types)
bClear = False
for k in S.ftupd:
if fver < k[0] and j == k[1]:
bClear = True
if bClear:
nr.append(0)
iatoms.append(None)
else:
if fver < 44:
fver_err(fver)
# 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, S.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
##############UTILS FOR INFORMATION NOT INTERESTED AT THE MOMENT###############
[docs]def do_grps(data): # pragma: no cover
grps_nr = []
myngrps = ngrps = S.egcNR # remind of version inconsistency
for j in range(ngrps):
if j < myngrps:
v = data.unpack_int()
grps_nr.append(v)
ndo_int(data, v)
return grps_nr
[docs]def do_groups(data, symtab): # pragma: no cover
do_grps(data)
ngrpname = data.unpack_int()
# do_strstr, list of indices of group name: e.g. System, Protein,
# Protein-H, etc. use symtab[i] to check
ndo_int(data, ngrpname)
ngrpnr = []
grpnr = []
for i in range(S.egcNR):
x = data.unpack_int()
ngrpnr.append(x)
if x == 0:
grpnr.append(None)
else:
l_ = []
for i in range(x):
l_.append(data.unpack_uint())
grpnr.append(l_)
# print ngrpnr
# print [len(i) for i in grpnr if i]
return
[docs]def do_atomtypes(data): # pragma: no cover
at_nr = data.unpack_int() # at: atomtype
at_radius = ndo_real(data, at_nr)
at_vol = ndo_real(data, at_nr)
at_surftens = ndo_real(data, at_nr)
at_atomnumber = ndo_int(data, at_nr)
return at_radius, at_vol, at_surftens, at_atomnumber