//
//  io_topol.cpp
//  Created by Arghya Chakravorty on 2/09/19.
//
//  Contains functions to read topology of molecular systems
//  in various formats.

#include "io.h"

#define PROTON2AMBER_CHARGE_CONVERSIONFACTOR 18.2223

// AMBER2PROTON_CHARGE_CONVERSIONFACTOR = 1.0/PROTON2AMBER_CHARGE_CONVERSIONFACTOR
#define AMBER2PROTON_CHARGE_CONVERSIONFACTOR 0.05487781454591352

#define BAD_GRO -1
#define SUCCESS 0
#define BAD_PSF -1
#define BAD_PRMTOP -3
#define BAD_ATOM_LINE -2
#define BAD_DCD -1
#define MAXLENGTH 80
#define VECTORDIM 3
#define TPX_TAG_RELEASE "release"
#define NR_CBTDIHS 6  //(look for idef.h)
#define NR_RBDIHS 6   //(look for idef.h)
#define NR_FOURDIHS 4 //(look for idef.h)
#define NOTSET -12345 //(look for typedefs.h)
#define MAXNODES 256  //(look for tpxio.c)

// reading PSF file
delphi_integer CIO::readPSF(ifstream &ifPsfFileStream, vector<CSimAtom> &vctOfSimAtoms)
{
    // resize:: reset vctOfAtoms
    vctOfSimAtoms.resize(0);

    // local variables
    bool bIsValidPSF = false;
    bool bHasCheq = false;
    bool bHasExt = false;
    bool bHasXplor = false;
    bool bCanReadAtom = false;
    bool bPostProcess = true; // delphi io_pdb

    int natoms = -1;
    int numRemarks = -1;
    // I,LSEGID,LRESID,LRES,TYPE(I),IAC(I),CG(I),AMASS(I),IMOVE(I),ECH(I),EHA(I)

    //     standard format:
    //     (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,I4,1X,2G14.6,I8,2G14.6)
    //     (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8,2G14.6)  XPLOR
    //     expanded format EXT:
    //     (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,I4,1X,2G14.6,I8,2G14.6)
    //     (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,A4,1X,2G14.6,I8,2G14.6) XPLOR
    int atomId, atomResId, lastResId = -6000;
    char chrAtomResname[10];
    char chrAtomName[10];

    char chrAtInf[16]; // need only 15
    int iGoodAtInf;

    string strLine, strSubLine;
    string strAtomId, strResId, strAtomResname, strAtomName, strQ;
    string strFmt, strFmtName, strMisc;
    delphi_real fCharge;
    delphi_real fMass;

    int varCount = 0, i = 0;
    int strLineSize = 0;

    // read the psf file line by line
    while (!ifPsfFileStream.eof())
    {

        // get a line
        getline(ifPsfFileStream, strLine);

        // to ignore possible empty lines
        if (0 == removeSpace(strLine).compare(0, 1, ""))
            continue;

        /* PSF format info from: https://userguide.mdanalysis.org/1.0.0/formats/reference/psf.html
        CHEQ:
        I,LSEGID,LRESID,LRES,TYPE(I),IAC(I),CG(I),AMASS(I),IMOVE(I),ECH(I),EHA(I)

        standard format:
        (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,I4,1X,2G14.6,I8,2G14.6)
        (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8,2G14.6)  XPLOR
        expanded format EXT:
        (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,I4,1X,2G14.6,I8,2G14.6)
        (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,A4,1X,2G14.6,I8,2G14.6) XPLOR

        no CHEQ
        II,LSEGID,LRESID,LRES,TYPE(I),IAC(I),CG(I),AMASS(I),IMOVE(I)

        standard format:
        (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,I4,1X,2G14.6,I8)
        (I8,1X,A4,1X,A4,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8)  XPLOR
        expanded format EXT:
        (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,I4,1X,2G14.6,I8)
        (I10,1X,A8,1X,A8,1X,A8,1X,A8,1X,A4,1X,2G14.6,I8) XPLOR
        */
        string fmtCheq = "%8d %4s %4s %4s %4s %4d %14.6f%14.6f%8d%14.6f%14.6f";
        string fmtCheqXplor = "%8d %4s %4s %4s %4s %4s %14.6f%14.6f%8d%14.6f%14.6f";
        string fmtCheqExt = "%10d %8s %8s %8s %8s %4d %14.6f%14.6f%8d%14.6f%14.6f";
        string fmtCheqExtXplor = "%10d %8s %8s %8s %8s %4s %14.6f%14.6f%8d%14.6f%14.6f";

        string fmtNoCheq = "%8d %4s %4s %4s %4s %4d %14.6f%14.6f%8d";
        string fmtNoCheqXplor = "%8d %4s %4s %4s %4s %4s %14.6f%14.6f%8d";
        string fmtNoCheqExt = "%10d %8s %8s %8s %8s %4d %14.6f%14.6f%8d";
        string fmtNoCheqExtXplor = "%10d %8s %8s %8s %8s %4s %14.6f%14.6f%8d";
        // look for "PSF" in the very first line
        if ((0 == strLine.compare(0, 3, "PSF")) && (!bIsValidPSF))
        {
            bIsValidPSF = true;
            if (strLine.find("CHEQ") != string::npos)
                bHasCheq = true;
            if (strLine.find("EXT") != string::npos)
                bHasExt = true;
            if (strLine.find("XPLOR") != string::npos)
                bHasXplor = true;
            cout << " Found a good PSF file with CHEQ=";
            cout << boolalpha << bHasCheq << " EXT=" << bHasExt;
            cout << " XPLOR=" << boolalpha << bHasXplor << endl;
            if (bHasCheq)
            {
                if (bHasExt)
                {
                    if (bHasXplor)
                    {
                        strFmt = fmtCheqExtXplor;
                        strFmtName = "fmtCheqExtXplor";
                    }
                    else
                    {
                        strFmt = fmtCheqExt;
                        strFmtName = "fmtCheqExt";
                    }
                }
                else
                {
                    if (bHasXplor)
                    {
                        strFmt = fmtCheqXplor;
                        strFmtName = "fmtCheqXplor";
                    }
                    else
                    {
                        strFmt = fmtCheq;
                        strFmtName = "fmtCheq";
                    }
                }
            }
            else
            {
                if (bHasExt)
                {
                    if (bHasXplor)
                    {
                        strFmt = fmtNoCheqExtXplor;
                        strFmtName = "fmtNoCheqExt";
                    }
                    else
                    {
                        strFmt = fmtNoCheqExt;
                        strFmtName = "fmtNoCheqExtXplor";
                    }
                }
                else
                {
                    if (bHasXplor)
                    {
                        strFmt = fmtNoCheqXplor;
                        strFmtName = "fmtNoCheqXplor";
                    }
                    else
                    {
                        strFmt = fmtNoCheq;
                        strFmtName = "fmtNoCheq";
                    }
                }
            }
            cout << "ATOM fromat string '" << strFmt << "'" << endl;
            continue;
        }

        // if valid PSF
        if (bIsValidPSF)
        {
            // SKIP REMARKS STATEMENTS
            if (strstr(strLine.c_str(), "!NTITLE"))
            {
                sscanf(strLine.c_str(), "%d %*s", &numRemarks);
                // sscanf(strLine.c_str(),"%d %s", &numRemarks, varString);
                cout << " Expected " << numRemarks << " remarks" << endl;
                // skip those many lines
                while (varCount < numRemarks)
                {
                    getline(ifPsfFileStream, strLine);
                    varCount++;
                }

                varCount = 0; // reset
                continue;
            }

            // RECOGNIZE START OF ATOM BLOCK
            if ((!bCanReadAtom) && strstr(strLine.c_str(), "!NATOM"))
            {
                bCanReadAtom = true;
                sscanf(strLine.c_str(), "%d %*s", &natoms);
                cout << " Reading data for " << natoms << " atoms" << endl;
                vctOfSimAtoms.resize(natoms);
                cout << " A size of " << vctOfSimAtoms.size() << " has been assigned to vctOfSimAtoms " << endl;

                // reset number of residues to 0
                iResidueNum = 0;

                // add io_pdb post processing for ony-molecule inputs
                bPostProcess = true;
                continue;
            }

            // READ ATOM BLOCK
            if (bCanReadAtom)
            {
                // SNAP OUT IF START OF ANY OF THESE BLOCKS IS ENCOUNTERED
                if (strstr(strLine.c_str(), "!NBOND") || strstr(strLine.c_str(), "!NTHETA") || strstr(strLine.c_str(), "!NPHI") || strstr(strLine.c_str(), "!NIMPHI") || strstr(strLine.c_str(), "!NDON") || strstr(strLine.c_str(), "!NACC") || strstr(strLine.c_str(), "!NGRP") || strstr(strLine.c_str(), "!NCRTERM"))
                {
                    bCanReadAtom = false;
                    break;
                }

                /*
                 * strLineSize = sscanf(strLine.c_str(), "%d %*s %d %s %s %*s %*s %*f", &atomId, &atomResId, chrAtomResname, chrAtomName);
                 * string strAtomName(chrAtomName);
                 * string strAtomResname(chrAtomResname);
                 */
                if (strFmtName == "fmtNoCheq" || strFmtName == "fmtCheq" || strFmtName == "fmtNoCheqXplor" || strFmtName == "fmtCheqXplor")
                {
                    // "%8d %4s %4s %4s %4s %4d %14.6f%14.6f%8d"
                    strAtomId = strLine.substr(0, 8);
                    strAtomId = removeSpace(strAtomId);
                    atomId = atoi(strAtomId.c_str()); // atomId (I8)
                    strResId = strLine.substr(14, 4);
                    strResId = removeSpace(strResId);
                    atomResId = atoi(strResId.c_str()); // resid (A4)
                    strAtomResname = strLine.substr(19, 4);
                    strAtomResname = removeSpace(strAtomResname); // resname (A4)
                    strAtomName = strLine.substr(24, 4);
                    strAtomName = removeSpace(strAtomName); // atom name (A4);
                    strQ = strLine.substr(34, 14);
                    strQ = removeSpace(strQ);
                    fCharge = atof(strQ.c_str()); // charge (G14.6)
                }
                else if (strFmtName == "fmtNoCheqExt" || strFmtName == "fmtCheqExt" || strFmtName == "fmtNoCheqExtXplor" || strFmtName == "fmtCheqExtXplor")
                {
                    // "%10d %8s %8s %8s %8s %4d %14.6f%14.6f%8d"
                    strAtomId = strLine.substr(0, 10);
                    strAtomId = removeSpace(strAtomId);
                    atomId = atoi(strAtomId.c_str()); // atomId (I8)
                    strResId = strLine.substr(20, 8);
                    strResId = removeSpace(strResId);
                    atomResId = atoi(strResId.c_str()); // resid (A4)
                    strAtomResname = strLine.substr(29, 8);
                    strAtomResname = removeSpace(strAtomResname); // resname (A4)
                    strAtomName = strLine.substr(38, 8);
                    strAtomName = removeSpace(strAtomName); // atom name (A4);
                    strQ = strLine.substr(52, 14);
                    strQ = removeSpace(strQ);
                    fCharge = atof(strQ.c_str()); // charge (G14.6)
                }

                CSimAtom newAtom;
                newAtom.setIndex(atomId);
                newAtom.setAtomResid(atomResId);
                newAtom.setAtomName(strAtomName);
                newAtom.setAtomResname(strAtomResname);
                newAtom.setCharge(fCharge);

                iGoodAtInf = sprintf(chrAtInf, " %-4s %3s %1s%4d", strAtomName.c_str(), strAtomResname.c_str(), " ", atomResId);
                const string strAtInf(chrAtInf);
                newAtom.setAtInf(strAtInf);

                if (atomResId != lastResId)
                {
                    iResidueNum++;
                    lastResId = atomResId;
                }

                // all the atoms will belong to the same media/object
                vctiAtomMediaNum.push_back(iObjectMediaNum);

                vctOfSimAtoms[i] = newAtom;
                i++;
                continue;

            } // bCanReadAtom
        }
        else
        {
            return BAD_PSF;
        }

    } // while

    if (bPostProcess)
    {
        if (0 == vctfMediaEps.size())
            vctfMediaEps.push_back(0.0); // medeps(0)

        // cout << "You are not reading from an object file! \n";
        // cout << "Assuming having only molecules, and one medium \n";
        // CNotReadingObjectFile warning;

        vctfMediaEps.push_back(fDielec / fEPKT); // medeps(1)=repsin/epkt
        // dataobject(1,1) and dataobject(1,2)
        vctstrObject.push_back(string("is a molecule    0"));
        vctstrObject.push_back(string(" "));
        prgiObjectMediaNum.push_back(iObjectMediaNum);
    }

    iObjectNum = vctstrObject.size() / 2;
    iMediaNum = vctfMediaEps.size() - 1;

    // regardless of being a molecule or not, iatmmed has a field to say which is its medium
    if (0 != prgiObjectMediaNum.size())
    {
        for (unsigned i = 0; i < prgiObjectMediaNum.size(); i++)
            vctiAtomMediaNum.push_back(prgiObjectMediaNum[i]);

        prgiObjectMediaNum.clear();
    }

    cout << " Read parameters for a total of " << vctOfSimAtoms.size() << " atoms from PSF file" << endl;

#ifdef VERBOSE
    vector<CSimAtom>::iterator ait = vctOfSimAtoms.begin();
    while (ait != vctOfSimAtoms.end())
    {
        cout << " " << ait->getIndex() << " " << ait->getAtomName() << " " << ait->getAtomResname() << ait->getAtomResid() << endl;
        ait++;
    }
#endif

    return natoms;
} // END of readPSF

// Read PRMTOP file, return number of atoms on success otherwise -1
// Shailesh: April 12, 2019
CIO::CParmtopReadHead CIO::readPRMTOP(ifstream &ifPrmtopFileStream, vector<CSimAtom> &vctOfSimAtoms)
{
    // resize:: reset vctOfAtoms
    vctOfSimAtoms.resize(0);

    // local variables
    bool bIsValidPRMTOP = false;
    bool bCanReadAtom = false;
    bool bPostProcess = true; // delphi io_pdb

    int nAtoms = -1;
    int numRemarks = -1;

    int atomId, atomResId, lastResId = -6000;
    char chrAtomResname[10];
    char chrAtomName[10];
    char chrAtInf[16]; // need only 15
    int iGoodAtInf;

    string strLine, strSubLine;
    string strAtomId, strResId, strAtomResname, strAtomName, strQ;
    delphi_real fCharge;

    int varCount = 0, i = 0;
    int strLineSize = 0;

    CIO::CParmtopReadHead readResult;
    // CIO::SPrmtopHeadPointers prmtopPointers = &(readResult.prmtopHead);
    CIO::SPrmtopHeadPointers prmtopPointers = readResult.prmtopHead;
    readResult.readSuccess = false;

    vector<CIO::CPrmtopSectionDataFormat> sectonsToRead;
    sectonsToRead.push_back(CPrmtopSectionDataFormat("POINTERS", CIO::enumPrmtopDataType::INTEGER, 8, 10, 0, 0, 0));
    sectonsToRead.push_back(CPrmtopSectionDataFormat("ATOM_NAME", CIO::enumPrmtopDataType::CHAR_ARRAY, 4, 20, 0, 0, 0));
    sectonsToRead.push_back(CPrmtopSectionDataFormat("CHARGE", CIO::enumPrmtopDataType::DOUBLE, 16, 5, 0, 1, 3));
    sectonsToRead.push_back(CPrmtopSectionDataFormat("RESIDUE_LABEL", CIO::enumPrmtopDataType::CHAR_ARRAY, 4, 20, 0, 0, 0));
    sectonsToRead.push_back(CPrmtopSectionDataFormat("RESIDUE_POINTER", CIO::enumPrmtopDataType::INTEGER, 8, 10, 0, 0, 0));
    // sectonsToRead.push_back(sectionDataFormat("RADII", 'E', 16, 5, 0, 1, 3));

    vector<string> atomNames;
    vector<double> atomCharges;
    vector<string> residueNames;
    vector<int> residuePointer;

    CIO::CPrmtopSectionDataFormat currentSection("", CIO::enumPrmtopDataType::INTEGER, 1, 1, 0, 0, 0);
    bool is_a_flag_line = false;

    while (!ifPrmtopFileStream.eof())
    {
        // Do not read the line if `strLine` already contains %FLAG
        if (!is_a_flag_line)
        {
            // get a line
            getline(ifPrmtopFileStream, strLine);
#ifdef VERBOSE
            cout << " io_topol> strLine from prmtop : " << strLine << endl;
#endif

            // to ignore possible empty lines
            if (0 == removeSpace(strLine).compare(0, 1, ""))
                continue;
        }

        // check and proceed if line starts with %FLAG
        // it indicates start of a PRMTOP section
        if (strLine.rfind("%FLAG", 0) == 0)
        {
            string sectionName;
            sectionName = strLine.substr(6);
            sectionName = removeSpace(sectionName);
#ifdef VERBOSE
            cout << "io_topol> section name : |" << sectionName << "|" << endl;
#endif

            is_a_flag_line = false;
            string sectionCommentOptional;

            getline(ifPrmtopFileStream, sectionCommentOptional);
            while (strLine.rfind("%COMMENT", 0) == 0)
            {
                getline(ifPrmtopFileStream, sectionCommentOptional);
            }

            string sectionFormat;
            // The line does not match the %COMMENT flag therefore it must be %FORMAT line
            if (strLine.rfind("%FORMAT", 0) == 0)
            {
                sectionFormat = sectionCommentOptional;
            }

            // getline(ifPrmtopFileStream, sectionFormat);
            for (auto index = 0; index < sectonsToRead.size(); index++)
            {
                if (sectionName == sectonsToRead[index].sectionName)
                {
                    currentSection = sectonsToRead[index];
                    break;
                }
            }

        } // if( strLine.rfind("%FLAG", 0) == 0)

        if (currentSection.sectionName == "POINTERS")
        {
            int n_values_to_read = 31;
            string strTmpLine1, strTmpLine2, strTmpLine3, strTmpLine4;
            getline(ifPrmtopFileStream, strTmpLine1);
            getline(ifPrmtopFileStream, strTmpLine2);
            getline(ifPrmtopFileStream, strTmpLine3);
            getline(ifPrmtopFileStream, strTmpLine4);

#ifdef VERBOSE
            cout << "io_topol> strTmpLine1 : " << strTmpLine1 << endl;
            cout << "io_topol> strTmpLine2 : " << strTmpLine2 << endl;
            cout << "io_topol> strTmpLine3 : " << strTmpLine3 << endl;
            cout << "io_topol> strTmpLine4 : " << strTmpLine4 << endl;
#endif

            readResult.prmtopHead.natoms = atoi(removeSpace(strTmpLine1.substr(0, 8)).c_str());
            readResult.prmtopHead.ntypes = atoi(removeSpace(strTmpLine1.substr(8, 8)).c_str());
            readResult.prmtopHead.nbonh = atoi(removeSpace(strTmpLine1.substr(16, 8)).c_str());
            readResult.prmtopHead.mbona = atoi(removeSpace(strTmpLine1.substr(24, 8)).c_str());
            readResult.prmtopHead.ntheth = atoi(removeSpace(strTmpLine1.substr(32, 8)).c_str());
            readResult.prmtopHead.mtheta = atoi(removeSpace(strTmpLine1.substr(40, 8)).c_str());
            readResult.prmtopHead.nphih = atoi(removeSpace(strTmpLine1.substr(48, 8)).c_str());
            readResult.prmtopHead.mphia = atoi(removeSpace(strTmpLine1.substr(56, 8)).c_str());
            readResult.prmtopHead.nhparm = atoi(removeSpace(strTmpLine1.substr(64, 8)).c_str());
            readResult.prmtopHead.nparm = atoi(removeSpace(strTmpLine1.substr(72, 8)).c_str());

            readResult.prmtopHead.nnb = atoi(removeSpace(strTmpLine2.substr(0, 8)).c_str());
            readResult.prmtopHead.nres = atoi(removeSpace(strTmpLine2.substr(8, 8)).c_str());
            readResult.prmtopHead.nbona = atoi(removeSpace(strTmpLine2.substr(16, 8)).c_str());
            readResult.prmtopHead.ntheta = atoi(removeSpace(strTmpLine2.substr(24, 8)).c_str());
            readResult.prmtopHead.nphia = atoi(removeSpace(strTmpLine2.substr(32, 8)).c_str());
            readResult.prmtopHead.numbnd = atoi(removeSpace(strTmpLine2.substr(40, 8)).c_str());
            readResult.prmtopHead.numang = atoi(removeSpace(strTmpLine2.substr(48, 8)).c_str());
            readResult.prmtopHead.nptra = atoi(removeSpace(strTmpLine2.substr(56, 8)).c_str());
            readResult.prmtopHead.natyp = atoi(removeSpace(strTmpLine2.substr(64, 8)).c_str());
            readResult.prmtopHead.nphb = atoi(removeSpace(strTmpLine2.substr(72, 8)).c_str());

            readResult.prmtopHead.ifpert = atoi(removeSpace(strTmpLine3.substr(0, 8)).c_str());
            readResult.prmtopHead.nbper = atoi(removeSpace(strTmpLine3.substr(8, 8)).c_str());
            readResult.prmtopHead.ngper = atoi(removeSpace(strTmpLine3.substr(16, 8)).c_str());
            readResult.prmtopHead.ndper = atoi(removeSpace(strTmpLine3.substr(24, 8)).c_str());
            readResult.prmtopHead.mbper = atoi(removeSpace(strTmpLine3.substr(32, 8)).c_str());
            readResult.prmtopHead.mgper = atoi(removeSpace(strTmpLine3.substr(40, 8)).c_str());
            readResult.prmtopHead.mdper = atoi(removeSpace(strTmpLine3.substr(48, 8)).c_str());
            readResult.prmtopHead.ifbox = atoi(removeSpace(strTmpLine3.substr(56, 8)).c_str());
            readResult.prmtopHead.nmxrs = atoi(removeSpace(strTmpLine3.substr(64, 8)).c_str());
            readResult.prmtopHead.ifcap = atoi(removeSpace(strTmpLine3.substr(72, 8)).c_str());

            readResult.prmtopHead.numextra = atoi(removeSpace(strTmpLine4.substr(0, 8)).c_str());
            if (strTmpLine4.length() == 16)
            {
                readResult.prmtopHead.ncopy = atoi(removeSpace(strTmpLine4.substr(8, 8)).c_str());
            }

            currentSection.sectionName = "";
        }
        else if (currentSection.sectionName == "ATOM_NAME")
        {
            delphi_integer n_values_to_read = readResult.prmtopHead.natoms;
            string strTmpLine;
            try
            {
                while (n_values_to_read > 0)
                {
                    getline(ifPrmtopFileStream, strTmpLine);
                    delphi_integer values_in_line = 0;

                    if (n_values_to_read >= currentSection.valuesPerLine)
                        values_in_line = currentSection.valuesPerLine;
                    else
                        values_in_line = n_values_to_read;

                    for (auto temp = 0; temp < values_in_line; ++temp)
                    {
                        string strTemp = strTmpLine.substr(temp * currentSection.fieldWidth, currentSection.fieldWidth);
                        atomNames.push_back(removeSpace(strTemp));
                        n_values_to_read--;
                    }
                }

                currentSection.sectionName = "";
            }
            catch (exception &ex)
            {
                readResult.readSuccess = false;
                cout << "Found number of atoms (" << atomNames.size() << ") is less than expected (" << readResult.prmtopHead.natoms << ") in PRMTOP" << endl;
            } // try-catch
        }
        else if (currentSection.sectionName == "CHARGE")
        {
            delphi_integer n_values_to_read = readResult.prmtopHead.natoms;
            string strTmpLine;
            try
            {
                while (n_values_to_read > 0)
                {
                    getline(ifPrmtopFileStream, strTmpLine);
                    delphi_integer values_in_line = 0;
                    if (n_values_to_read >= currentSection.valuesPerLine)
                        values_in_line = currentSection.valuesPerLine;
                    else
                        values_in_line = n_values_to_read;

                    for (auto temp = 0; temp < values_in_line; ++temp)
                    {
                        string strValue = removeSpace(strTmpLine.substr(temp * currentSection.fieldWidth, currentSection.fieldWidth));
                        atomCharges.push_back(atof(strValue.c_str()) * AMBER2PROTON_CHARGE_CONVERSIONFACTOR);
                        n_values_to_read--;
                    }
                }
                currentSection.sectionName = "";
            }
            catch (exception &ex)
            {
                readResult.readSuccess = false;
                cout << "Found number of atom charges (" << atomCharges.size() << ") is less than expected (" << readResult.prmtopHead.natoms << ") in PRMTOP"
                     << endl;
            }
        } // CHARGE
        else if (currentSection.sectionName == "RESIDUE_LABEL")
        {
            delphi_integer n_values_to_read = readResult.prmtopHead.nres;
            string strTmpLine;
            try
            {
                while (n_values_to_read > 0)
                {
                    getline(ifPrmtopFileStream, strTmpLine);
                    delphi_integer values_in_line = 0;

                    if (n_values_to_read >= currentSection.valuesPerLine)
                        values_in_line = currentSection.valuesPerLine;
                    else
                        values_in_line = n_values_to_read;

                    for (auto temp = 0; temp < values_in_line; ++temp)
                    {
                        string strValue = removeSpace(strTmpLine.substr(temp * currentSection.fieldWidth, currentSection.fieldWidth));
                        residueNames.push_back(strValue);
                        n_values_to_read--;
                    }
                }
                currentSection.sectionName = "";
            }
            catch (exception &ex)
            {
                readResult.readSuccess = false;
                cout << "Found number of residue labels (" << residueNames.size() << ") is less than expected (" << readResult.prmtopHead.nres << ") in PRMTOP"
                     << endl;
            }
        } // RESIDUE_LABEL
        else if (currentSection.sectionName == "RESIDUE_POINTER")
        {
            delphi_integer n_values_to_read = readResult.prmtopHead.nres;
            string strTmpLine;
            while (n_values_to_read > 0)
            {
                getline(ifPrmtopFileStream, strTmpLine);
                delphi_integer values_in_line = 0;
                if (n_values_to_read >= currentSection.valuesPerLine)
                {
                    values_in_line = currentSection.valuesPerLine;
                }
                else
                {
                    values_in_line = n_values_to_read;
                }
                for (auto temp = 0; temp < values_in_line; ++temp)
                {
                    string strValue = removeSpace(strTmpLine.substr(temp * currentSection.fieldWidth, currentSection.fieldWidth));
                    residuePointer.push_back(atoi(strValue.c_str()));
                    n_values_to_read--;
                }
            }

            currentSection.sectionName = "";

        } // RESIDUE_POINTER

        else
        {
            // skip data for this section
            while (!ifPrmtopFileStream.eof())
            {
                getline(ifPrmtopFileStream, strLine);
                if (strLine.rfind("%FLAG", 0) == 0)
                {
                    is_a_flag_line = true;
                    break;
                }
            }
        } // else

    } // while (!ifPrmtopFileStream.eof())

    atomResId = 0;
    if (atomNames.size() == atomCharges.size())
    {
        for (auto tmp = 0; tmp < atomNames.size(); ++tmp)
        {
            atomId = tmp + 1; // Prmtop does not store atomid, it starts from 1 and changes with atomIndex
            // atomResId is not stored explicitly and it starts from 1 in prmtop
            // 1-based index of first-atom-of-residues in stored in `residuePointer`
            // there when `tmp+1` 1 based matches with first-atom-of-residue, resid is incremented by 1
            if (tmp + 1 == residuePointer[atomResId])
            {
                atomResId++;
            }
            strAtomResname = residueNames[atomResId - 1];
            strAtomName = atomNames[tmp];
            fCharge = atomCharges[tmp];

            CSimAtom newAtom;
            newAtom.setIndex(atomId);
            newAtom.setAtomResid(atomResId);
            newAtom.setAtomName(strAtomName);
            newAtom.setAtomResname(strAtomResname);
            newAtom.setCharge(fCharge);

            iGoodAtInf = sprintf(chrAtInf, " %-4s %3s %1s%4d", strAtomName.c_str(), strAtomResname.c_str(), " ", atomResId);
            const string strAtInf(chrAtInf);
            newAtom.setAtInf(strAtInf);

            if (atomResId != lastResId)
            {
                iResidueNum++;
                lastResId = atomResId;
            }

            // all the atoms will belong to the same media/object
            vctiAtomMediaNum.push_back(iObjectMediaNum);
            vctOfSimAtoms.push_back(newAtom);
        }
        readResult.readSuccess = true;
    }

    if (bPostProcess)
    {
        if (0 == vctfMediaEps.size())
            vctfMediaEps.push_back(0.0); // medeps(0)

        // cout << "You are not reading from an object file! \n";
        // cout << "Assuming having only molecules, and one medium \n";
        // CNotReadingObjectFile warning;

        vctfMediaEps.push_back(fDielec / fEPKT); // medeps(1)=repsin/epkt
        // dataobject(1,1) and dataobject(1,2)
        vctstrObject.push_back(string("is a molecule    0"));
        vctstrObject.push_back(string(" "));
        prgiObjectMediaNum.push_back(iObjectMediaNum);
    }

    iObjectNum = vctstrObject.size() / 2;
    iMediaNum = vctfMediaEps.size() - 1;

    // regardless of being a molecule or not, iatmmed has a field to say which is its medium
    if (0 != prgiObjectMediaNum.size())
    {
        for (unsigned i = 0; i < prgiObjectMediaNum.size(); i++)
            vctiAtomMediaNum.push_back(prgiObjectMediaNum[i]);

        prgiObjectMediaNum.clear();
    }

    nAtoms = vctOfSimAtoms.size();
    cout << " Read parameters for a total of " << vctOfSimAtoms.size() << " atoms from PRMTOP file" << endl;

#ifdef VERBOSE
    vector<CSimAtom>::iterator ait = vctOfSimAtoms.begin();
    while (ait != vctOfSimAtoms.end())
    {
        cout << " " << ait->getIndex() << " " << ait->getAtomName() << " " << ait->getAtomResname() << ait->getAtomResid() << endl;
        ait++;
    }
#endif

    return readResult;
} // END of readPRMTOP

// reading TPR file
int CIO::readTPR(ifstream &ifTprHandler, vector<CSimAtom> &vctOfSimAtoms)
{

    // collect size (total number of bytes)
    ifTprHandler.seekg(0, ios_base::end);
    long int totalSize = ifTprHandler.tellg();
    ifTprHandler.seekg(0, ios_base::beg);

    int numAtoms = -1;

    // local variables
    string strAtomResname, strAtomName;
    delphi_real fCharge;
    int atomId, atomResId, lastResId = -6000;
    bool bPostProcess = true; // delphi io_pdb

    int precision, fileVersion, iBox, num_gtc;
    int num_moltypes;

    vector<delphi_real> vctAtomCharges;
    vector<delphi_integer> vctAtomResIds;
    vector<string> vctAtomNames;
    vector<string> vctResNames;

    while (ifTprHandler.tellg() != totalSize)
    {

        // read header
        readTPRHeader(ifTprHandler, numAtoms, precision, fileVersion, iBox, num_gtc);
        cout << " Number of atoms =  " << numAtoms << endl;
        cout << " Number of t-couple groups = " << num_gtc << endl;
        cout << " ---------------- HEADER FINSIHED ----------------- " << endl;

        // box information (skip bytes)
        readBox(ifTprHandler, fileVersion, precision, iBox);
        cout << " -------------- BOX FINISHED ------------------ " << endl;

        // t_couple_groups (Skip bytes)
        read_tCouplGroups(ifTprHandler, fileVersion, precision, num_gtc);
        cout << " -------------- T-Couple FINISHED ------------------ " << endl;

        // fetch symbol table (IMPORTATN)
        vector<string> symTable = readSymTable(ifTprHandler);
        int num_symbols = (int)symTable.size();

        // get symstr (skip)
        skip_n<int>(ifTprHandler, 1);

        /**************************
         *     FF PARAMS          *
         **************************/
        skip_n<int>(ifTprHandler, 1); // atnr

        if (fileVersion < 57)
            skip_n<int>(ifTprHandler, 1); // idum

        // ff params (skip carefully)
        readFFParams(ifTprHandler, fileVersion, precision);

        // moltypes
        num_moltypes = decode<int>(ifTprHandler);
        cout << " Total number of moltypes in the system : " << num_moltypes << endl;

        // create arrays that will store atom and residue count for each of these molecules
        int atomCounts[num_moltypes];
        int resCounts[num_moltypes];

        // for each moltype, do a certain set of operations
        for (int iMol = 0; iMol < num_moltypes; iMol++)
        {

            // do the equivalent of `do_moltype` function

            // start
            if (fileVersion >= 57)
            {
                // get symstr (skip)
                skip_n<int>(ifTprHandler, 1);
            }

            int mol_natoms = decode<int>(ifTprHandler);
            atomCounts[iMol] = mol_natoms;
            int mol_nres = decode<int>(ifTprHandler);
            resCounts[iMol] = mol_nres;

            cout << " Found a molecule with " << mol_natoms << " atoms and " << mol_nres << " residues." << endl;

            if (fileVersion < 57)
                // group number (skip)
                skip_n<int>(ifTprHandler, 1);

            // do the equivalent of `do_atoms`
            vctAtomCharges.resize(mol_natoms);
            vctAtomResIds.resize(mol_natoms);
            vctAtomNames.resize(mol_natoms);
            vctResNames.resize(mol_nres);

            // do the equivalent of `do_atoms`
            if (mol_natoms > 0)
            {

                delphi_real fCharge = 0.0;
                uint uiType;

                for (int iAtom = 0; iAtom < mol_natoms; iAtom++)
                {
                    // mass (skip)
                    skip_nreal(ifTprHandler, 1, precision);

                    // charge (IMPORTANT)
                    if (precision == sizeof(float))
                        fCharge = decode<float>(ifTprHandler);
                    else if (precision == sizeof(double))
                        fCharge = (float)decode<double>(ifTprHandler);

                    vctAtomCharges[iAtom] = fCharge;

                    // mass B (skip)
                    skip_nreal(ifTprHandler, 1, precision);

                    // charge B (skip)
                    skip_nreal(ifTprHandler, 1, precision);

                    // type of atom
                    uiType = decode<uint>(ifTprHandler);

                    // type of atom B
                    skip_n<uint>(ifTprHandler, 1);

                    // ptype : 'Atom' or 'vitual site', etc.
                    skip_n<int>(ifTprHandler, 1);

                    // resId
                    vctAtomResIds[iAtom] = decode<int>(ifTprHandler);

                    if (fileVersion >= 52)
                    {
                        // skip atom number
                        skip_n<int>(ifTprHandler, 1);
                    }

                    // skipping parts where group numbers are assigned to atoms
                    // not necessary here.

                } // for iAtom 0 -> mol_natoms

                // fetch atom names
                int idx;
                string atomName;
                for (int iAtom = 0; iAtom < mol_natoms; iAtom++)
                {
                    idx = decode<int>(ifTprHandler);
                    atomName = symTable[idx];
                    vctAtomNames[iAtom] = atomName;

                } // for iAtom 0 -> mol_natoms

                // atomtype assignment (skip)
                if (fileVersion > 20)
                {
                    skip_n<int>(ifTprHandler, mol_natoms);

                    // type B assignment
                    skip_n<int>(ifTprHandler, mol_natoms);
                }

            } // if mol_natoms > 0

            // do the equivalent of `do_resinfo`
            if (mol_nres > 0)
            {
                int idx;
                string resName;
                for (int iRes = 0; iRes < mol_nres; iRes++)
                {
                    if (fileVersion >= 63)
                    {
                        idx = decode<int>(ifTprHandler);
                        resName = symTable[idx];
                        vctResNames[iRes] = resName;

                        // skip resid because we already have it availabale separately
                        // from a diferet function
                        skip_n<int>(ifTprHandler, 2);
                    }

                } // for iRes : 0 -> mol_nres
            }     // if mol_nres > 0

            // do the equivalent of `do_ilists`
            for (int iEnerType = 0; iEnerType < F_NRE; iEnerType++)
            {
                // do the equivalent of `do_ilist`
                if (fileVersion < 44)
                    skip_n<int>(ifTprHandler, 1); // idum

                // length of atom list that contains atoms which participate
                // in a certain kind of interaction.
                // Will be important when doing MM part
                int lenAtomList = decode<int>(ifTprHandler);

                // skip details for now
                skip_n<int>(ifTprHandler, lenAtomList);
            }

            // do the equivalent of `do_block`
            if (fileVersion < 44)
                skip_n<int>(ifTprHandler, MAXNODES);

            int lenBlock = decode<int>(ifTprHandler);

            int dum_nra = 0;
            if (fileVersion < 51)
                dum_nra = decode<int>(ifTprHandler);

            // skip `lenBlock + 1` real quantities
            skip_nreal(ifTprHandler, lenBlock + 1, precision);

            if (fileVersion < 51 && dum_nra > 0)
                skip_n<int>(ifTprHandler, dum_nra);

            // do the equivalent of `do_blocka`
            int nr = decode<int>(ifTprHandler);
            int nra = decode<int>(ifTprHandler);

            skip_n<int>(ifTprHandler, nr + 1);
            skip_n<int>(ifTprHandler, nra);

        } // iMol: 0 -> num_molTypes

        // we should have all the important information by this point.
        // SO EXIT!!!!
        break;

    } // while ifTprHandler.tellg() != totalSize
    cout << " iotopol> Finished reading the TPR file via a while loop" << endl;

    // now create instances of CSimAtoms using the information
    // obtained from reading the binary TPR file.
    vctOfSimAtoms.resize(numAtoms);
    this->iResidueNum = 0;
    char chrAtInf[16]; // need only 15
    int iGoodAtInf;

    for (int idx = 0; idx < numAtoms; idx++)
    {

        atomId = idx + 1;
        atomResId = vctAtomResIds[idx] + 1;
        strAtomResname = vctResNames[vctAtomResIds[idx]];
        strAtomResname = removeSpace(strAtomResname);
        strAtomName = vctAtomNames[idx];
        strAtomName = removeSpace(strAtomName);
        fCharge = vctAtomCharges[idx];

        CSimAtom newAtom;
        newAtom.setIndex(atomId);
        newAtom.setAtomResid(atomResId);
        newAtom.setAtomName(strAtomName);
        newAtom.setAtomResname(strAtomResname);
        newAtom.setCharge(fCharge);

        iGoodAtInf = sprintf(chrAtInf, " %-4s %3s %1s%4d", strAtomName.c_str(), strAtomResname.c_str(), " ", atomResId);
        const string strAtInf(chrAtInf);
        newAtom.setAtInf(strAtInf);

        if (atomResId != lastResId)
        {
            this->iResidueNum++;
            lastResId = atomResId;
        }

        // all the atoms will belong to the same media/object
        vctiAtomMediaNum.push_back(iObjectMediaNum);

        vctOfSimAtoms[idx] = newAtom;

#ifdef VERBOSE
        cout << "ATOM  ";
        cout << setw(5) << right << idx + 1;
        cout << " " << setw(4) << vctAtomNames[idx];
        cout << " " << setw(3) << vctResNames[vctAtomResIds[idx]];
        cout << " " << setw(4) << 1 + vctAtomResIds[idx];
        cout << " " << setw(10) << vctAtomCharges[idx];
        cout << endl;
#endif
    } // idx : 0 -> numAtoms

    if (bPostProcess)
    {
        if (0 == vctfMediaEps.size())
            vctfMediaEps.push_back(0.0); // medeps(0)

        // cout << "You are not reading from an object file! \n";
        // cout << "Assuming having only molecules, and one medium \n";
        // CNotReadingObjectFile warning;

        vctfMediaEps.push_back(fDielec / fEPKT); // medeps(1)=repsin/epkt
        // dataobject(1,1) and dataobject(1,2)
        vctstrObject.push_back(string("is a molecule    0"));
        vctstrObject.push_back(string(" "));
        prgiObjectMediaNum.push_back(iObjectMediaNum);
    }

    iObjectNum = vctstrObject.size() / 2;
    iMediaNum = vctfMediaEps.size() - 1;

    // regardless of being a molecule or not, iatmmed has a field to say which is its medium
    if (0 != prgiObjectMediaNum.size())
    {
        for (unsigned i = 0; i < prgiObjectMediaNum.size(); i++)
            vctiAtomMediaNum.push_back(prgiObjectMediaNum[i]);

        prgiObjectMediaNum.clear();
    }

    cout << " Read parameters for a total of " << vctOfSimAtoms.size() << " atoms from GROMACS TPR file" << endl;

#ifdef VERBOSE
    vector<CSimAtom>::iterator ait = vctOfSimAtoms.begin();
    while (ait != vctOfSimAtoms.end())
    {
        cout << " " << ait->getIndex() << " " << ait->getAtomName() << " " << ait->getAtomResname() << ait->getAtomResid() << endl;
        ait++;
    }
#endif
    return numAtoms;
}

// reading TPR header
void CIO::readTPRHeader(ifstream &fileStream, int &numAtoms, int &precision, int &fileVersion, int &iBox, int &n_gtc)
{

    int varint;
    string varstring;

    // skip the first integer
    skip_n<int>(fileStream, 1);

    // check for need to reverse bytes
    varint = decode<int>(fileStream);
    cout << " Strsize before swapping : " << varint << endl;
    if (varint > MAXLENGTH)
    {
        varint = swapbytes<int>(varint);
        cout << " Strsize after swapping : " << varint << endl;
        if (varint <= MAXLENGTH)
        {
            // endianness difference detected
            bSwapBytes = true;
        }
    } // check for need to reverse bytes done

    // place the seeker at the 4 byte position
    fileStream.seekg(sizeof(int), ios_base::beg);

    // version
    string version = decode_gro_string(fileStream);
    cout << " Version : " << version << " of length " << strlen(version.c_str()) << endl;
    cout << " Present seeker position at : " << fileStream.tellg() << endl;

    // precision
    precision = decode<int>(fileStream);
    cout << " Precision : " << precision << endl;

    // file version
    fileVersion = decode<int>(fileStream);
    cout << " File Version : " << fileVersion << endl;

    // file tag
    string fileTag;
    if (fileVersion >= 77 && fileVersion <= 79)
    {
        fileTag = decode_gro_string(fileStream);
        cout << " File Tag : " << fileTag << " of length " << strlen(fileTag.c_str()) << endl;
    }

    // generation
    int fileGen = -1;
    if (fileVersion >= 26)
    {
        fileGen = decode<int>(fileStream);
        cout << " File generation : " << fileGen << endl;
    }
    else
        fileGen = 0;

    if (fileVersion >= 81)
    {
        varint = decode<int>(fileStream);
        fileTag = decode_gro_string(fileStream);
        cout << " File tag for version > 81 : " << fileTag << " of length " << strlen(fileTag.c_str()) << endl;
    }
    else if (fileVersion < 77)
    {
        cout << " will do nothing " << endl;
    }

    // number of atoms
    numAtoms = decode<int>(fileStream);

    // # of t-coupl groups
    n_gtc = decode<int>(fileStream);

    // idum and rdum (skip)
    if (fileVersion < 62)
        skip_n<int>(fileStream, 2);

    // fep state (skip)
    if (fileVersion >= 79)
        skip_n<int>(fileStream, 1);

    // lambda (skip)
    skip_nreal(fileStream, 1, precision);

    // bIr (input record?), bTop (topology?), bX (coords?), bV (velocities?), bF (forces?) (bool)
    // skip 5 ints
    skip_n<int>(fileStream, 5);

    // has Box
    iBox = decode<int>(fileStream);

    return;
}

// read(skip) FF PARAMS
void CIO::readFFParams(ifstream &fileStream, int &fileVersion, int &precision)
{
    // function types
    int num_fTypes = decode<int>(fileStream);
    int funcTypes[num_fTypes];
    int iFunc = 0;
    while (iFunc < num_fTypes)
    {
        funcTypes[iFunc] = decode<int>(fileStream);
        iFunc++;
    }

    // reppow (skip double)
    if (fileVersion >= 66)
    {
        skip_n<double>(fileStream, 1);
    }

    // fudge QQ (skip)
    skip_nreal(fileStream, 1, precision);

    // Gromacs checks wether all the fucntions defined are compatible with the code.
    // this is for backward compatibility and the numbers are being altered from old
    // to new numbers for funcTypes

    // read parameters for all the function types read from the TPR file
    // (basically skip them because we dont need them now)
    for (int i = 0; i < num_fTypes; i++)
    {
        iFunc = funcTypes[i];
        if ((iFunc == F_ANGLES) || (iFunc == F_G96ANGLES) || (iFunc == F_BONDS) || (iFunc == F_G96BONDS) || (iFunc == F_HARMONIC) || (iFunc == F_IDIHS))
            skip_nreal(fileStream, 4, precision);

        else if (iFunc == F_RESTRANGLES)
            skip_nreal(fileStream, 2, precision);

        else if (iFunc == F_LINEAR_ANGLES)
            skip_nreal(fileStream, 4, precision);

        else if (iFunc == F_FENEBONDS)
            skip_nreal(fileStream, 2, precision);

        else if (iFunc == F_RESTRBONDS)
            skip_nreal(fileStream, 8, precision);

        else if ((iFunc == F_TABBONDS) || (iFunc == F_TABBONDSNC) || (iFunc == F_TABANGLES) || (iFunc == F_TABDIHS))
        {
            skip_nreal(fileStream, 1, precision);
            skip_n<int>(fileStream, 1);
            skip_nreal(fileStream, 1, precision);
        }
        else if (iFunc == F_CROSS_BOND_BONDS)
            skip_nreal(fileStream, 3, precision);

        else if (iFunc == F_CROSS_BOND_ANGLES)
            skip_nreal(fileStream, 3, precision);

        else if (iFunc == F_UREY_BRADLEY)
        {
            skip_nreal(fileStream, 4, precision);
            if (fileVersion >= 79)
                skip_nreal(fileStream, 4, precision);
        }
        else if (iFunc == F_QUARTIC_ANGLES)
            skip_nreal(fileStream, 1 + 5, precision);

        else if (iFunc == F_BHAM)
            skip_nreal(fileStream, 3, precision);

        else if (iFunc == F_MORSE)
        {
            skip_nreal(fileStream, 3, precision);
            if (fileVersion >= 79)
                skip_nreal(fileStream, 3, precision);
            else
                fake_a_death();
        }

        else if (iFunc == F_CUBICBONDS)
            skip_nreal(fileStream, 3, precision);

        else if (iFunc == F_CONNBONDS)
            fake_a_death();

        else if (iFunc == F_POLARIZATION)
            skip_nreal(fileStream, 1, precision);

        else if (iFunc == F_ANHARM_POL)
            skip_nreal(fileStream, 3, precision);

        else if (iFunc == F_WATER_POL)
        {
            if (fileVersion < 31)
            {
                fake_a_death();
            }
            else
                skip_nreal(fileStream, 6, precision);
        }

        else if (iFunc == F_THOLE_POL)
            skip_nreal(fileStream, 4, precision);

        else if (iFunc == F_LJ)
            skip_nreal(fileStream, 2, precision);

        else if (iFunc == F_LJ14)
            skip_nreal(fileStream, 4, precision);

        else if (iFunc == F_LJC14_Q)
            skip_nreal(fileStream, 5, precision);

        else if (iFunc == F_LJC_PAIRS_NB)
            skip_nreal(fileStream, 4, precision);

        else if ((iFunc == F_PDIHS) || (iFunc == F_PIDIHS) || (iFunc == F_ANGRES) || (iFunc == F_ANGRESZ))
        {
            skip_nreal(fileStream, 2, precision);
            if (((iFunc == F_ANGRES) || (iFunc == F_ANGRESZ)) && fileVersion < 42)
            {
                skip_nreal(fileStream, 2, precision);
            }
            else
            {
                skip_nreal(fileStream, 2, precision);
                skip_n<int>(fileStream, 1);
            }
        }

        else if (iFunc == F_RESTRDIHS)
            skip_nreal(fileStream, 2, precision);

        else if (iFunc == F_DISRES)
        {
            skip_n<int>(fileStream, 2);
            skip_nreal(fileStream, 2, precision);
        }

        else if (iFunc == F_ORIRES)
        {
            skip_n<int>(fileStream, 3);
            skip_nreal(fileStream, 3, precision);
        }

        else if (iFunc == F_DIHRES)
        {
            if (fileVersion < 82)
                skip_n<int>(fileStream, 2);

            skip_nreal(fileStream, 3, precision);
            if (fileVersion >= 82)
            {
                skip_nreal(fileStream, 3, precision);
            }
            else
            {
                fake_a_death();
            }
        }

        else if (iFunc == F_POSRES)
        {
            // POSITIONS
            skip_nreal(fileStream, 3, precision);

            // forces
            skip_nreal(fileStream, 3, precision);

            if (fileVersion < 27)
            {
                fake_a_death();
            }
            else
            {
                skip_nreal(fileStream, 3, precision); //(B)
                skip_nreal(fileStream, 3, precision); //(B)
            }
        }

        else if (iFunc == F_FBPOSRES)
        {
            skip_n<int>(fileStream, 1);
            skip_nreal(fileStream, 3, precision); // vector
            skip_nreal(fileStream, 2, precision);
        }

        else if (iFunc == F_CBTDIHS)
            skip_nreal(fileStream, NR_CBTDIHS, precision);

        else if (iFunc == F_RBDIHS)
        {
            skip_nreal(fileStream, NR_RBDIHS, precision);
            if (fileVersion > 25)
                skip_nreal(fileStream, NR_RBDIHS, precision);
        }

        else if (iFunc == F_FOURDIHS)
            skip_nreal(fileStream, 2 * NR_RBDIHS, precision);

        else if ((iFunc == F_CONSTR) || (iFunc == F_CONSTRNC))
            skip_nreal(fileStream, 2, precision);

        else if (iFunc == F_SETTLE)
            skip_nreal(fileStream, 2, precision);

        else if (iFunc == F_VSITE2)
            skip_nreal(fileStream, 1, precision);

        else if ((iFunc == F_VSITE3) || (iFunc == F_VSITE3FD) || (iFunc == F_VSITE3FAD))
            skip_nreal(fileStream, 2, precision);

        else if ((iFunc == F_VSITE3OUT) || (iFunc == F_VSITE4FD) || (iFunc == F_VSITE4FDN))
            skip_nreal(fileStream, 3, precision);

        else if (iFunc == F_VSITEN)
        {
            skip_n<int>(fileStream, 1);
            skip_nreal(fileStream, 1, precision);
        }

        else if ((iFunc == F_GB12) || (iFunc == F_GB13) || (iFunc == F_GB14))
        {
            if (fileVersion < 68)
                skip_nreal(fileStream, 4, precision);

            skip_nreal(fileStream, 5, precision);
        }

        else if (iFunc == F_CMAP)
            skip_n<int>(fileStream, 2);

        else
            cout << " unknown function type with id " << iFunc << " in the set of func with ids in range of 0 to " << num_fTypes << endl;

    } // iFunc: 0 -> num_fTypes

    return;

} // END of readFFParams

// a dummy useless function that just "passes"
void CIO::fake_a_death()
{
    return;
}

// skip numInstances of real objects based on the precision of the gromacs file
void CIO::skip_nreal(ifstream &fileStream, int numInstances, int real_precision)
{
    if (real_precision == sizeof(float))
        skip_n<float>(fileStream, numInstances);
    else if (real_precision == sizeof(double))
        skip_n<double>(fileStream, numInstances);
    else
        return;
}

// read/skip box information in TPR files
void CIO::readBox(ifstream &fileStream, int &fileVersion, int &precision, int &iBox)
{
    // skip various box information
    if ((bool)iBox)
    {
        // skip 3 x 3 * (precision) bytes
        skip_nreal(fileStream, VECTORDIM * VECTORDIM, precision);

        // box_rel (skip)
        if (fileVersion >= 51)
        {
            skip_nreal(fileStream, VECTORDIM * VECTORDIM, precision);
        }

        // box_v
        skip_nreal(fileStream, VECTORDIM * VECTORDIM, precision);
        // mdum
        if (fileVersion < 56)
        {
            skip_nreal(fileStream, VECTORDIM * VECTORDIM, precision);
        }

    } // iBox
} // END of readBox

// read/skip information about t_couple groups
void CIO::read_tCouplGroups(ifstream &fileStream, int &fileVersion, int &precision, int &n_gtc)
{

    if (fileVersion < 69)
    {
        skip_nreal(fileStream, n_gtc, precision);
    }

    // redundant call (?)
    skip_nreal(fileStream, n_gtc, precision);

    return;

} // END of read_tCouplGroups

// fetch all the symbols used in a GROMACS TPR file
vector<string> CIO::readSymTable(ifstream &fileStream)
{
    // number of symbols
    int num_symbols = decode<int>(fileStream);
    cout << " Expecting " << num_symbols << " symbols in the topology." << endl;

    // create a vector of symbols (strings)
    vector<string> vctSymbols(num_symbols, " ");

    int iSymbol = 0;
    while (iSymbol < num_symbols)
    {
        skip_n<int>(fileStream, 1);
        vctSymbols[iSymbol] = decode_gro_string(fileStream);
        iSymbol++;
    }

    /* print the symbols
     iSymbol = 0;
     while (iSymbol < vctSymbols.size())
     {
     cout << iSymbol << " : " << vctSymbols[iSymbol] << right << setw(10) << " of length " << strlen(vctSymbols[iSymbol].c_str()) << endl;
     iSymbol++;
     }
     */
    return vctSymbols;

} // END of readSymTable
