# -*- 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,
)


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]

    def unpack_int64(self):
        return self._unpack_value(8, '>q')

    def unpack_uint64(self):
        return self._unpack_value(8, '>Q')

    def unpack_ushort(self):
        return self.unpack_uint()

    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')

    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()


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).
    """
    @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

    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

    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')

    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')

    def do_string(self):
        """
        Emulate gmx_fio_do_string
        """
        n = self.unpack_uint64()
        return self.unpack_fstring(n)


def ndo_int(data, n):
    """mimic of gmx_fio_ndo_real in gromacs"""
    return [data.unpack_int() for i in range(n)]


def ndo_real(data, n):
    """mimic of gmx_fio_ndo_real in gromacs"""
    return [data.unpack_real() for i in range(n)]


def do_rvec(data):
    return data.unpack_farray(setting.DIM, data.unpack_real)


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)]


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)]


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 "
        )


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}")


def read_tpxheader(data):
    """this function is now compatible with do_tpxheader in tpxio.cpp
    """
    # Last compatibility check with gromacs-2016
    ver_str = data.do_string()  # version string e.g. VERSION 4.0.5
    if not ver_str.startswith(b'VERSION'):
        raise ValueError('Input does not look like a TPR file.')
    precision = data.unpack_int()  # e.g. 4
    define_unpack_real(precision, data)
    fileVersion = data.unpack_int()  # version of tpx file
    fileVersion_err(fileVersion)

    # This is for backward compatibility with development version 77-79 where
    # the tag was, mistakenly, placed before the generation.
    # These versions are development versions between the 4.5 and 4.6 series.
    if 77 <= fileVersion <= 79:
        data.unpack_int()  # the value is 8, but haven't found the
        file_tag = data.do_string()

    fileGeneration = data.unpack_int() if fileVersion >= 26 else 0  # generation of tpx file, e.g. 17

    # Versions before 77 don't have the tag, set it to TPX_TAG_RELEASE file_tag
    # file_tag is used for comparing with tpx_tag. Only tpr files with a
    # tpx_tag from a lower or the same version of gromacs code can be parsed by
    # the tpxio.c

    file_tag = data.do_string() if fileVersion >= 81 else setting.TPX_TAG_RELEASE

    natoms = data.unpack_int()  # total number of atoms
    ngtc = data.unpack_int() if fileVersion >= 28 else 0  # number of groups for T-coupling

    if fileVersion < 62:
        # not sure what these two are for.
        data.unpack_int()  # idum
        data.unpack_real()  # rdum

    fep_state = data.unpack_int() if fileVersion >= 79 else 0

    # actually, it's lambda, not sure what is it. us lamb because lambda is a
    # keyword in python
    lamb = data.unpack_real()
    bIr = data.unpack_int()  # has input record or not
    bTop = data.unpack_int()  # has topology or not
    bX = data.unpack_int()  # has coordinates or not
    bV = data.unpack_int()  # has velocity or not
    bF = data.unpack_int()  # has force or not
    bBox = data.unpack_int()  # has box or not

    sizeOfTprBody = None
    if fileVersion >= setting.tpxv_AddSizeField and fileGeneration >= 27:
        sizeOfTprBody = data.unpack_int64()

    th = obj.TpxHeader(ver_str, precision, fileVersion, fileGeneration,
                       file_tag, natoms, ngtc, fep_state, lamb,
                       bIr, bTop, bX, bV, bF, bBox, sizeOfTprBody)
    return th


def extract_box_info(data, fver):
    box = ndo_rvec(data, setting.DIM)
    box_rel = ndo_rvec(data, setting.DIM)
    box_v = ndo_rvec(data, setting.DIM)

    return obj.Box(box, box_rel, box_v)


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


def do_symstr(data, symtab):
    # do_symstr: get a string based on index from the symtab
    ndx = data.unpack_int()
    return symtab[ndx]


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


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


def do_harm(data):
    data.unpack_real()  # rA
    data.unpack_real()  # krA
    data.unpack_real()  # rB
    data.unpack_real()  # krB


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


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


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)


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


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)


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)
    ]


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)


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


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