/**
 * @file io.h
 * @brief IO class for reading/writing various standard files
 *
 * @author Chuan Li, chuanli@clemson.edu
 */

#ifndef CIO_H_
#define CIO_H_

#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>   // std::string
#include <vector>   // std::vector
#include <locale>
#include <stdlib.h> // atof
#include <cmath>

#include "../interface/environment.h"
#include "../misc/misc_grid.h"
#include "io_datatype.h"
#include "io_exceptions.h"
#include "../delphi/delphi_constants.h"

using namespace std;

const int SIZEFILE   = 1;
const int CHARGEFILE = 2;
const int VDWFILE    = 3;
const int SIGMAGAUSSFILE = 4; // Added by Shailesh: Nov 7, 2019

const int STDPDB     = 10;
const int MODPDB     = 20;
const int PQRPDB     = 21;
const int GROPDB     = 22;
const int MOD4PDB    = 30;
const int PQR4PDB    = 31;

const int DCDTRAJ    = 101;
const int XTCTRAJ    = 102;
const int TRRTRAJ    = 103;
const int MDCRDTRAJ  = 104;
const int NCTRAJ     = 105;

const int PSFTOP     = 201;
const int TPRTOP     = 202;
const int PRMTOP     = 204;


/* CommPDB format is updated and used for communication with DelPhi and PrimePKA without read/write files
 * Updated: Dec 15, 2014 by Lin Wang
 */
const int CommPDB    = 40;


class CIO
{
   protected:
      /**
       * F95 var.: repsin
       */
      const delphi_real fDielec;

      /**
       * F95 var.: epkt
       */
      const delphi_real fEPKT;

      /*
       * ******************** miscellany functions defined in io_misc.cpp *******************
       */

      /**
       * function to convert a string to upper case
       * @param strLine input string
       * @return converted string in upper case
       */
      string toUpperCase(const string& strLine);

      /**
       * function to remove spaces from a string
       * @param strLine input string
       * @return converted string without spaces
       */
      string removeSpace(const string& strLine);

      /**
       * function to check if the input file is formatted or not
       * @param strFile name of the input file
       * @return true if the file formatted, false otherwise.
       */
      bool checkFileFormat(const string& strFile);

      /*
       * ******************* for reading atom force (radii/size and charge) file *******************
       * ********** these static attributes are declared here and defined in io_force.cpp **********
       */

      /**
       * F95 var.: nrmax
       * # of entries in radius file
       */
      static delphi_integer iRadiusNum;

      /**
       * F95 var.: radii(:)
       * a vector containing all radius read from .siz file
       */
      static vector<CForce> prgas;

      /**
       * F95 var.: nchrec
       * # of entries in charge file
       */
      static delphi_integer iCrgNum;

      /**
       * F95 var.: charge(:)
       * a vector containing all charges read from .crg file
       */
      static vector<CForce> prgac;

      /**
       * No F95 equivalent
       * # of entries in vdw file
       */
      static delphi_integer iVdwNum;

      /**
       * No F95 equivalent
       * var.: NON-POLAR PARAMS(:)
       * a vector containing all sigma-epsilon-gamma read from .vdw file
       */
      static vector<CForce> prgan;

      /**
       * F95 var.: None
       * # of entries in sigma-gaussian file
       */
      static delphi_integer iSggNum;

      /**
       * No F95 equivalent
       * var.: GAUSIAN-SIGMA PARAMS(:)
       * a vector containing all sigma-gaussian read from .sgg file
       */
      static vector<CForce> prgasg;

     /**
       * read a .siz or .crg file in non-PK format
       * @param ifFileStream input file stream
       * @param iFileType input file type
       */
      void readFileInNotPKFormat(ifstream& ifFileStream, const int& iFileType);

      /**
       * read a .siz or .crg file in PK format
       * @param ifFileStream input file stream
       * @param iFileType input file type
       */
      void readFileInPKFormat(ifstream& ifFileStream, const int& iFileType);

      /**
       * read a .vdw file 
       * @param ifFileStream input file stream
       */
      void readFileInVdwFormat(ifstream& ifFileStream);

#ifdef DEBUG_IO_FORCE
      /**
       * debug function to print vector storing values read from .siz or .crg file
       * @param iFileType input file type
       * @prgaf pointer pointing to the vector of values read from .siz or .crg file
       * @prgaf pointer now also will point to the vector of values read from .vdw file
       */
      void printForce(const int& iFileType, const vector<CForce>& prgaf);
      /** 
       * debug function to print vecot storing the values read from .siz/.crg/.sgg file
       * @param iFileType input file type
       */
      void printForce(const int& iFileType);
#endif

      /**
       * function to get index of a specified atom in the record
       * @param strAtom atom name
       * @param strResidue residue name
       * @param strResidueNum residue number
       * @param strChain chain name
       * @param iRecordNum total number of records
       * @param prgaf pointer to the vector of records
       * @return -1 if no record is found, index otherwise.
       */
      delphi_integer FindRecordIndex(const string& strAtom,const string& strResidue,const string& strResidueNum,
                                     const string& strChain,const delphi_integer& iRecordNum,vector<CForce>* prgaf);
      
      /**
       * function to get index of a specified atom in the record
       * @param strAtom atom name
       * @param iRecordNum total number of records
       * @param prgaf pointer to the vector of records
       * @return -1 if no record is found, index otherwise.
       */
      delphi_integer FindVdwRecordIndex(const string& strAtom,const delphi_integer& iRecordNum,vector<CForce>* prgaf);

      /**
       * function to find a particular record matching the specified atom
       * @param strAtom atom name
       * @param strResidue residue name
       * @param strResidueNum residue number
       * @param strChain chain name
       * @param iFileType file type
       * @param fValue value (charge or radius) associated with this atom
       * @return -1 if unfounded, index otherwise
       */
      delphi_integer FindRecord(const string& strAtom,const string& strResidue,const string& strResidueNum,
                                const string& strChain,const int& iFileType, delphi_real& fValue);

      
      /**
       * function to find a particular record matching the specified atom
       * @param strAtom atom name
       * @param strResidue residue name
       * @param strResidueNum residue number
       * @param strChain chain name
       * @param iFileType file type
       * @param fValue value (sigma-gaussian) associated with this atom
       * @return -1 if unfounded, index otherwise
       */
      delphi_integer FindSggRecord(const string& strAtom,const string& strResidue,const string& strResidueNum,
                                const string& strChain,const int& iFileType, delphi_real& fValue);


      /**
       * function to find a particular record matching the specified atom
       * @param strAtom atom name
       * @param fVdwRmin value associated with this atom
       * @param fVdwEpsilon value associated with this atom
       * @param fVdwGamma value associated with this atom
       * @return -1 if unfounded, index otherwise
       */
      delphi_integer FindVdwRecord(const string& strAtom,  delphi_real& fVdwRmin, delphi_real& fVdwEpsilon, delphi_real& fVdwGamma);

      //void checkCharge(const bool& bSurfCrgInSite);

      /*
       * ********** functions defined in io_pdb.cpp for reading PDB file in various formats ********
       */

      /**
       * F95 var.: imedianumb
       * number of objects read from .pdb file
       */
      delphi_integer iObjectMediaNum;

      /**
       * function to read standard .pdb file
       * @param ifPdbFileStream input .pdb file stream
       */
      void readStdPdbFile(ifstream& ifPdbFileStream);

      /**
       * function to read modified .pdb file (type = 1)
       * @param ifPdbFileStream input .pdb file stream
       */
      void readModFile1(ifstream& ifPdbFileStream);

      /**
       * function to read modified .pdb file (type = 4)
       * @param ifPdbFileStream input .pdb file stream
       */
      void readModFile4(ifstream& ifPdbFileStream);

      /**
       * function to read .pqr file
       * @param ifPdbFileStream input .pqr file stream
       */
      void readPqrFile(ifstream& ifPdbFileStream);

      /**
       * function to read mod4 pdb file
       * @param ifPdbFileStream input .pdb file stream
       */
      void readMod4File(ifstream& ifPdbFileStream);

      /**
       * function to read pqr4 file
       * @param ifPdbFileStream input .pdb file stream
       */
      void readPqr4File(ifstream& ifPdbFileStream);

      /**
       * function to read gromacs' .gro file
       * @param ifPdbFileStream input .gro file stream
       */
      void readGroFile(ifstream& ifPdbFileStream);
      
      /**
       * function to read unformatted pdb file
       * @param strPdbFile input .pdb file name
       * @param ifPdbFileStream input .pdb file stream
       * @param bPostProcess flag for post processing
       */
      void readUnformattedPdb(const string& strPdbFile,ifstream& ifPdbFileStream,bool& bPostProcess);

      void commVecPDB(const vector<string>& strCommPDB);

#ifdef DEBUG_IO_PDB
      /**
       * debug function to print values read from .pdb file
       */
      void printPDB();
#endif

      /**
       * F95 var.: iatrad
       * whether there is radius info
       */
      bool    bExistRadiiInfo;

      /**
       * F95 var.: tmpiatmmed(nobjectmax)
       * vector containing internal media-number per object
       */
      vector<delphi_integer> prgiObjectMediaNum;

   public:

      /**
       * F95 var.: nmedia
       * # of media
       */
      delphi_integer    iMediaNum;

      /**
       * F95 var.: nobject
       * # of objects
       */
      delphi_integer    iObjectNum;

      /**
       * F95 var.: natom
       * # of atoms
       */
      delphi_integer    iAtomNum;

      /**
       * F95 var.: resnummax
       * maximum residue number
       */
      delphi_integer    iResidueNum;

      /**
       * F95 var.: ionlymol
       * true if there are only molecules in the system (no objects)
       */
      bool       bOnlyMolecule;

      /**
       * F95 var.: delphipdb(Natom)
       * array of structure to store info read from pdb file
       */
      vector<CAtomPdb> vctapAtomPdb;

      /**
       * F95 var.: medeps(0:nmediamax)
       * vector containing correspondence media<->epsilon/epkt
       */
      vector<delphi_real>    vctfMediaEps;

      /**
       * F95 var.: dataobject(nobjectmax,2)
       * vector containing string with object data, and pre-elab data
       * changed it to vctstrObject(2*nobjectmax)
       */
      vector<string>  vctstrObject;

      /**
       * F95 var.: iatmmed(Natom+Nobjectmax)
       * vector containing internal media-number per atom and object
       */
      vector<delphi_integer> vctiAtomMediaNum;

      /**
       * ************ MMPBSA/ TRAJECTORY ANALYSIS ****************
       */

      // vector containing the frames read from a trajectory
      vector<CFrame> vctFrames;

      // vector containing objects of Class CSimAtom which contains the information about the atoms
      // of the simulation
      vector<CSimAtom> vctSimAtoms;

      // should you swap bytes if the endianness of the source and host machine are different
      bool bSwapBytes;

      // ENUM for GROMACS energy function types
      enum enFuncType 
      {
        F_BONDS ,  F_G96BONDS ,  F_MORSE ,
        F_CUBICBONDS ,  F_CONNBONDS , F_HARMONIC ,
        F_FENEBONDS ,  F_TABBONDS ,  F_TABBONDSNC ,
        F_RESTRBONDS , F_ANGLES ,  F_G96ANGLES ,
        F_RESTRANGLES ,  F_LINEAR_ANGLES ,  F_CROSS_BOND_BONDS ,
        F_CROSS_BOND_ANGLES ,  F_UREY_BRADLEY ,  F_QUARTIC_ANGLES ,
        F_TABANGLES ,  F_PDIHS ,    F_RBDIHS ,
        F_RESTRDIHS ,  F_CBTDIHS ,  F_FOURDIHS ,
        F_IDIHS , F_PIDIHS ,  F_TABDIHS ,
        F_CMAP ,  F_GB12 ,  F_GB13 ,
        F_GB14 ,  F_GBPOL ,  F_NPSOLVATION ,
        F_LJ14 ,  F_COUL14 , F_LJC14_Q ,
        F_LJC_PAIRS_NB ,  F_LJ ,     F_BHAM ,
        F_LJ_LR , F_BHAM_LR ,  F_DISPCORR ,
        F_COUL_SR ,  F_COUL_LR ,  F_RF_EXCL ,
        F_COUL_RECIP ,  F_LJ_RECIP ,  F_DPD ,
        F_POLARIZATION ,  F_WATER_POL ,     F_THOLE_POL ,
        F_ANHARM_POL ,  F_POSRES ,  F_FBPOSRES ,
        F_DISRES ,     F_DISRESVIOL ,  F_ORIRES ,
        F_ORIRESDEV ,  F_ANGRES ,  F_ANGRESZ ,
        F_DIHRES ,  F_DIHRESVIOL ,  F_CONSTR ,
        F_CONSTRNC ,  F_SETTLE , F_VSITE2 ,
        F_VSITE3 ,  F_VSITE3FD ,  F_VSITE3FAD ,
        F_VSITE3OUT , F_VSITE4FD ,  F_VSITE4FDN ,
        F_VSITEN ,  F_COM_PULL ,  F_EQM ,
        F_EPOT ,  F_EKIN ,  F_ETOT ,
        F_ECONSERVED ,  F_TEMP , F_VTEMP_NOLONGERUSED ,
        F_PDISPCORR ,  F_PRES ,  F_DVDL_CONSTR ,
        F_DVDL , F_DKDL ,  F_DVDL_COUL ,
        F_DVDL_VDW ,  F_DVDL_BONDED ,  F_DVDL_RESTRAINT ,
        F_DVDL_TEMPERATURE ,  F_NRE
      };

      /**
       * constructor I using default Dielectric constant and epkt values
       */
      CIO():fDielec(4.0),fEPKT(167100.9162872952/297.3342119)
      {
#ifdef DEBUG_OBJECT
         cout << endl;
         cout << "****************************************************************\n";
         cout << "*                     CIO is constructed                       *\n";
         cout << "****************************************************************\n";
#endif

         //iRadiusNum      = 0;
         //iCrgNum         = 0;
         iMediaNum       = 1;
         iObjectNum      = 1;
         iAtomNum        = 0;
         iObjectMediaNum = 1;
         bOnlyMolecule   = true;
         bExistRadiiInfo = false;
         iResidueNum     = 0;
         bSwapBytes      = false;
      };

      /**
       * constructor II using user-specified Dielectric constant and epkt values
       */
      CIO(delphi_real fDielecIn,delphi_real fEPKTIn):fDielec(fDielecIn),fEPKT(fEPKTIn)
      {
#ifdef DEBUG_OBJECT
         cout << endl;
         cout << "****************************************************************\n";
         cout << "*                     CIO is constructed                       *\n";
         cout << "****************************************************************\n";
#endif

         //iRadiusNum      = 0;
         //iCrgNum         = 0;
         iMediaNum       = 1;
         iObjectNum      = 1;
         iAtomNum        = 0;
         iObjectMediaNum = 1;
         bOnlyMolecule   = true;
         bExistRadiiInfo = false;
         iResidueNum     = 0;
         bSwapBytes      = false;
      };

      /**
       * destructor
       */
      ~CIO()
      {
#ifdef DEBUG_OBJECT
         cout << endl;
         cout << "****************************************************************\n";
         cout << "*                      CIO is destroyed                        *\n";
         cout << "****************************************************************\n";
#endif
      };


      /**
       * function for reading atom force (radii/size and charge) file
       * @param strFile file name
       */
      void readForceFile(const string& strFile);
      
      /**
       * function for reading atom vdw (sigma-epsilon-gamma) file
       * @param strFile file name
       */
      void readVdwFile(const string& strFile);

      /**
       * function for reading PDB file in various formats
       * @param strPdbFile pdb file name
       * @param iPdbFormat pdb file format
       * @param bPdbUnformat true if unformatted pdb file
       */
      void readPdbFile(const string& strPdbFile,const int& iPdbFormat,const bool& bPdbUnformat, const vector<string>& strCommPDB);

      /**
       * function for writing unformatted pdb file
       * @param strPdbFile pdb file name
       */
      void writeUnformatPdb(const string& strPdbFile);

      /**
       * function for writing modified pdb file
       * @param strPdbFile pdb file name
       * @param iModPdbFormatOut file format
       */
      void writeModifiedPdb(const string& strPdbFile,const int& iModPdbFormatOut);

      /*
       * miscellany
       */

      /**
       * function to set DelPhi-style atom list vctapAtomPdb(iAtomNum), i.e, (delphipdb(natom))
       * @param bSolvePB true if solving PB equation
       * @param bMultiSigmaGauss true if multi-sigma-gaussian used
       * @param bSurfCrgInSite true if writing surface chagre info into site file
       * @param strSizeFile name of .siz file
       * @param strCrgFile name of .crg file
       * @param strSggFile name of .sgg (sigma-gaussian: a .crg format) file
       * @param fSigma the single-sigma-gaussian value supplied by user or (default=1.0 set in delphi_datamarshel_setDefaults.cpp)
       * @param strPdbFile name of .pdb file
       * @param iPdbFormat .pdb file format
       * @param bPdbUnformat true if unformatted pdb file
       * @param bLJEng true if computing Lennard Jonees energy 
       * @param bNonPolEng true if computing non-polar energies
       * @param strVdwFile name of .vdw file
       */
      void setDelphiAtom(const bool& bSolvePB,       const bool& bUseGaussian,  const bool& bMultiSigmaGauss,     
                         const bool& bSurfCrgInSite, const string& strSizeFile, const string& strCrgFile,         
                         const string& strSggFile,   const delphi_real& fSigma, const string& strPdbFile,         
                         const int& iPdbFormat,      const bool& bPdbUnformat,  const vector<string>& strCommPDB, 
                         const bool& bLJEng,         const bool& bNonPolEng,    const string& strVdwFile);

      // void setDelphiAtom(const bool& bSolvePB,const bool& bSurfCrgInSite,const string& strSizeFile,
      //                    const string& strCrgFile,const string& strPdbFile,const int& iPdbFormat,
      //                    const bool& bPdbUnformat, const vector<string>& strCommPDB,
      //                    const bool& bLJEng, const bool& bNonPolEng, const string& strVdwFile);

      /**
       * function for reading FRC file
       * @param strFrcFile name of frc file
       * @param fgOffCenter offset from box center
       * @param fScale scale used for calculations
       * @return position of box center
       */
      SGrid<delphi_real> readFrcFile(const string& strFrcFile,const SGrid<delphi_real>& fgOffCenter,
                                     const delphi_real& fScale);

      /** OLD OLD OLD OLD OLD OLD OLD OLD OLD
       * function for writing EPS file
       * @param iAtomNumIn number of atoms
       * @param iObjectNumIn number of objects
       * @param iGrid number of grid points on each direction
       * @param fScale scale used for calculations
       * @param fgBoxCenter position of box center
       * @param vctigEpsMap vector of epsilon map
       * @param vctbDielecMap vector of dielectric constants map
       * @param strEpsFile output eps file name
       */
      // void writeEpsMap(const delphi_integer& iAtomNumIn,const delphi_integer& iObjectNumIn,
      //                  const delphi_integer& iGrid,const delphi_real& fScale,const SGrid<delphi_real>& fgBoxCenter,
      //                  const vector< SGrid<delphi_integer> >& vctigEpsMap,const vector<bool>& vctbDielecMap,
      //                  const string& strEpsFile);

    /** function for writing EPS file for 2-dielectric model
      * only CUBE format
      * @param iGrid number of grid points on each direction
      * @param repsout solvent dielctric value
      * @param repsin  solute dielectric value
      * @param fScale scale used for calculations
      * @param fgBoxCenter the coordinates of the box center
      * @param vctbDielecMap vector of dielectric constants map
      * @param strEpsFile output eps file name
      */

      void writeHomoEpsMap(const delphi_integer& iGrid,
                            const delphi_real&    repsout,
                            const delphi_real&    repsin,
                            const delphi_real& fScale,
                            const SGrid<delphi_real>& fgBoxCenter,
                            const vector<char>& vctbDielecMap,
                            const string& strEpsFile);

      /** function for writing EPS file for gaussian model
        * only CUBE format
        * @param iGrid number of grid points on each direction
        * @param repsout solvent dielctric value
        * @param fScale scale used for calculations
        * @param vct_cepsmap vector of dielectric constants map
        * @param strEpsFile output eps file name
        */

        void writeGaussEpsMap(const delphi_integer& iGrid,
                              const delphi_real&    repsout,
                              const delphi_real& fScale,
                              const SGrid<delphi_real>& fgBoxCenter,
                              const vector< SGrid<delphi_real> >& vctigEpsMap,
                              const string& strEpsFile);

      /* function for writing any grid-based map in CUBE format
	   * @param iGrid number of grid points on each direction
	   * @param fScale scale used for calculations
	   * @param fgBoxCenter position of box center
	   * @param vctMap vector containing the grid data in T datatype. The data will be printed as float.
	   * @param strOutFile output file name
	   */
	   template <class T> void writeMapInCube(const delphi_integer& iGrid,
								const delphi_real& fScale,
								const SGrid<delphi_real>& fgBoxCenter,
								const vector<T>& vctMap,
								const string& strOutFile);

      /*
       * for reading/writing PHI file
       */

      /*
       *for writing ENERGY file
       */

      /*
       * for writing HSURF2 file
       */

      /*
       *for writing SURFEN file
       */

      /*
       * for writing SCRG file
       */
        
      /*
       *
       * ************ MMPBSA/ TRAJECTORY ANALYSIS ****************
       *
       *
       * function to collect information about coordinates of atoms from a MD simulation trajectory
       * and their properties from the corresponding topology file.
       * @param strTrajFile name of trajectory file (.dcd, .xtc, .trr)
       * @param iTrajFormatIn trajectory file type
       * @param strTopolFile name of topology file (.psf, .tpr)
       * @param iTopolFormatIn topology file type
       * @param strSizeFile name of .siz file
       * @param strCrgFile name of .crg file
       * @param strVdwFile name of .vdw file
       * @param bSolvePB true if solving PB equation
       * @param bLJEng true if computing Lennard Jonees energy 
       * @param bNonPolEng true if computing non-polar energies
       * @param iFrameFirst index of the first frame to be read (0-indexed)
       * @param iFrameLast  index of the last frame to be read (0-indexed)
       * @param iFrameStride number of frames to skip between two consecutively read frames
       */
      void setDelphiTraj(const string& strTrajFile, int& iTrajFormatIn,
                         const string& strTopolFile, int& iTopolFormatIn,
                         const string& strSizeFile,
                         const string& strCrgFile,
                         const string& strVdwFile,
                         const bool& bSolvePB,
                         const bool& bLJEng,
			 const bool& bNonPolEng, 
			 const int& iFrameFirst,
                         const int& iFrameLast,
                         const int& iFrameStride);


      /**
       * function that returns a boolean flag to indicate if a frame from the trajectory should be read 
       * when provided with the boundary indices and stride.
       * @param iFrameCurrent index of the current frame
       * @param iFrameFirst index of the first frame to be read (0-indexed)
       * @param iFrameLast  index of the last frame to be read (0-indexed)
       * @param iFrameStride number of frames to skip between two consecutively read frames
       * @return true/false
       */
      bool shouldRecordFrame(int& iFrameCurrent, const int& iFrameFirst, const int& iFrameLast, const int& iFrameStride );

      /* 
       * function to unpack the next few bytes corresponding to the size of T
       * from a binary input stream.
       * @param fileStream input file stream
       * @return decoded variable of type T
       */
      template <class T> T decode(ifstream& fileStream);

      /*
       * function to unpack the next few bytes into a string in a way 
       * gromacs does it. Every string starts with an int that indicates the
       * length of the string that is to be unoacked.
       * @param fileStream input file stream
       * @return decoded string
       */
      string decode_gro_string(ifstream& fileStream);

      /*
       * funciton that removes the forbidden characters in a string (32 < ASCII < 127)
       * that might be present when unpacking a string from a binary stream.
       * @param str
       */
      void remove_forbidden_chars(string& str);

      /*
       * function to swap bytes if endianness of the host and source machnies are different.
       * @param input : variable whose bytes are to be swapped
       * @return Of class T whose bytes are reversed wrt to that of the input.
       */
      template <class T> T swapbytes(T input);

      /*
       * function that will skip a number of bytes corresponding to a certain datatype T
       * @param fileStream input file stream
       * @param numInstances number of instabces of datatype T placed consecutively that is to be skipped.
       */
      template <class T> void skip_n(ifstream& fileStream, int numInstances);
      
      /*
       * **********     functions defined in io_dcd.cpp for     ********  
       * ********** reading simulation frames from NAMD/CHARMM  ********
       * **********          DCD trajectory file                ********
       */

      /**
       * function to read blocks sizes NAMD/CHARMM DCD format trajectory file
       * Block sizes are at the beginning of each block in the DCD files and 
       * the same block size is at the end. They work like parenthesis. 
       * @param fileStream : .dcd file stream
       * @param blockSize : to read these many bytes
       */
      void readBlockSize(ifstream& fileStream, int& blockSize);

      /**
       * function to skip some number of bytes in the input AMD/CHARMM DCD format trajectory file
       * @param fileStream : .dcd file stream
       * @param nbytes : to skip these many bytes
       */
      void skipBytes(ifstream& fileStream, int& nbytes);

      /**
       * function to read a NAMD/CHARMM DCD format trajectory file
       * and gather coordinates of atoms in the simulation and make objects of Class CFrame.
       * @param ifDcdFileStream : .dcd file stream
       * @param vctOfFrames : vector of CFrame that contains coordinates of atoms in each snapshot.
       * @param vctAtomRadii : vector of radius values per atoms
       * @param iFrameFirst index of the first frame to be read (0-indexed)
       * @param iFrameLast  index of the last frame to be read (0-indexed)
       * @param iFrameStride number of frames to skip between two consecutively read frames
       * return -1 if reading wasnt successful which is possible due to bad DCD file.
       * return `# of atoms` if reading was successful.
       */
      delphi_integer readDCD(ifstream& ifDcdFileStream, 
                             vector<CFrame>& vctOfFrames, 
                             vector<delphi_real>& vctAtomRadii,
                             const int& iFrameFirst, 
                             const int& iFrameLast, 
                             const int& iFrameStride);
      /*
       * ********** functions defined in io_psf.cpp for rading atom properties from a NAMD/CHARMM PSF topology file ********
       */

      /**
       * function to read atom props from a NAMD/CHARMM PSF file.
       * It reads the atom and residue info of each atom. The corresponding trajectory (.dcd)
       * file is expected to have the coordinates of the atoms in the order of their appearence 
       * in this topology file. 
       * @param ifPsfFileStream : .psf file stream
       * @param vctOfSimAtoms  : vector of Class CSimAtom that contains for every atom in the simulation the attributes contained in the PSF file.
       * return -1 if reading wasnt successful which is possible due to incorrect format of PSF file.
       * return `# of atoms` if reading was successful.
       */
      delphi_integer readPSF(ifstream& ifPsfFileStream, vector<CSimAtom>& vctOfSimAtoms);

      
	/* ***************** BEGIN AMBER PRMTOP READING CODES ******************
	 * **********     functions defined in io_ambertraj.cpp for     ********
	 * **********        reading simulation frames from AMBER       ********
	 * **********          MDCRD trajectory file                    ********
	 * ********** Author :: Shailesh Kumar Panday
	 * ********** Date   :: Apr 14, 2019
	 **+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++**/

	enum enumPrmtopDataType 
	{
		INTEGER = 1, DOUBLE = 3, CHAR_ARRAY = 5
	};

	struct SPrmtopHeadPointers {
		delphi_integer natoms;
		delphi_integer ntypes;
		delphi_integer nbonh;
		delphi_integer mbona;
		delphi_integer ntheth;
		delphi_integer mtheta;
		delphi_integer nphih;
		delphi_integer mphia;
		delphi_integer nhparm;
		delphi_integer nparm;

		delphi_integer nnb;
		delphi_integer nres;
		delphi_integer nbona;
		delphi_integer ntheta;
		delphi_integer nphia;
		delphi_integer numbnd;
		delphi_integer numang;
		delphi_integer nptra;
		delphi_integer natyp;
		delphi_integer nphb;

		delphi_integer ifpert;
		delphi_integer nbper;
		delphi_integer ngper;
		delphi_integer ndper;
		delphi_integer mbper;
		delphi_integer mgper;
		delphi_integer mdper;
		delphi_integer ifbox;
		delphi_integer nmxrs;
		delphi_integer ifcap;

		delphi_integer numextra;
		delphi_integer ncopy;
	};

	class CPrmtopSectionDataFormat {
	public:
		string sectionName;
		enumPrmtopDataType dataType;
		unsigned int fieldWidth;
		unsigned int valuesPerLine;
		unsigned int leadingZeros;
		unsigned int fieldsBeforePoint;
		unsigned int fieldsForExponent;

		CPrmtopSectionDataFormat(string name, enumPrmtopDataType dataType,
				unsigned int width,
				unsigned int valuesPerLine, unsigned int leadingZeros,
				unsigned int fieldsBeforePoint,
				unsigned int filedsForExponent) {
			this->sectionName = name;
			this->dataType = dataType;
			this->fieldWidth = width;
			this->valuesPerLine = valuesPerLine;
			this->leadingZeros = leadingZeros;
			this->fieldsBeforePoint = fieldsBeforePoint;
			this->fieldsForExponent = fieldsForExponent;
		}
	};

	// delphi_integer readSectionOfPRMTOP(ifstream& ifPsfFileStream,
	//		CPrmtopSectionDataFormat& currentSection, vector<T> sectionData);

	class CParmtopReadHead {
	public:
		SPrmtopHeadPointers prmtopHead;
		bool readSuccess;
	};
	/**
	 * function to read atom props from a AMBER PRMTOP file.
	 * It reads the atom and residue info of each atom. The corresponding trajectory (.mdcrd)
	 * file is expected to have the coordinates of the atoms in the order of their appearence
	 * in this topology file.
	 * @param ifPrmtopFileStream : .prmtop file stream
	 * @param vctOfSimAtoms  : vector of Class CSimAtom that contains for every atom in the simulation the attributes contained in the PSF file.
	 * return -1 if reading wasnt successful which is possible due to incorrect format of PRMTOP file.
	 * return `# of atoms` if reading was successful.
	 */
	CParmtopReadHead readPRMTOP(ifstream& ifPsfFileStream,
			vector<CSimAtom>& vctOfSimAtoms);

	delphi_integer readMDCRD(ifstream& ifMdcrdFileStream,
			SPrmtopHeadPointers prmtopHeadPointer,
			vector<CFrame>& vctOfFrames, vector<delphi_real>& vctAtomRadii,
			const int& iFrameFirst, const int& iFrameLast,
			const int& iFrameStride);

	/* *******************************  END  *******************************
	 * **********      AMBER trajectory functions declaration       ********
	 *---------------------------------------------------------------------*/
      /* **********     functions defined in io_trr.cpp for     ********  
       * ********** reading simulation frames from GROMACS      ********
       * **********          TRR trajectory file                ********
       */

      /**
       * function to read a  gromacs TRR format trajectory file
       * and gather coordinates of atoms in the simulation and make objects of Class CFrame.
       * @param ifTrrFileStream : .trr file stream
       * @param vctOfFrames : vector of CFrame that contains coordinates of atoms in each snapshot.
       * @param vctAtomRadii : vector of radius values per atoms
       * @param iFrameFirst index of the first frame to be read (0-indexed)
       * @param iFrameLast  index of the last frame to be read (0-indexed)
       * @param iFrameStride number of frames to skip between two consecutively read frames
       * return -1 if reading wasnt successful which is possible due to bad DCD file.
       * return `# of atoms` if reading was successful.
       */
      delphi_integer readTRR(ifstream& ifTrrFileStream, 
                             vector<CFrame>& vctOfFrames, 
                             vector<delphi_real>& vctAtomRadii, 
                             const int& iFrameFirst, 
                             const int& iFrameLast, 
                             const int& iFrameStride);
    

      //reading TPR function
      int readTPR(ifstream& ifTprHandler, vector<CSimAtom>& vctOfSimAtoms);
      
      // skip float/double
      void skip_nreal(ifstream& fileStream, int numInstances, int real_precision);

      // a dummy useless function that just "passes"
      void fake_a_death();

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

      // read/skip box information
      void readBox(ifstream& fileStream, int& iBox, int& fileVersion, int& precision);

      // read/skip t_coupl_group info
      void read_tCouplGroups(ifstream& fileStream, int& fileVersion, int& precision, int& n_gtc);

      // read/fetch symbol table
      vector<string> readSymTable(ifstream& fileStream);

      // read(skip) ff params
      void readFFParams(ifstream& fileStream, int& fileVersion, int& precision);

};

/****************************************************
 *  DEFINING TEMPLATE FUNCTIONS OF TEH IO CLASS
 ****************************************************
 */

template <class T> T CIO::decode(ifstream& fileStream)
{
    T out;
    fileStream.read( reinterpret_cast<char *> (&out), sizeof(T) );
    if (this->bSwapBytes) out = this->swapbytes<T>(out);
    return out;

} // END OF FUNCTION


template <class T> T CIO::swapbytes(T input)
{
    T swapped;
    char * inputRepr = reinterpret_cast<char *>(&input);
    char * outputRepr = reinterpret_cast<char *>(&swapped);

    for (int k = 0; k < sizeof(T); k++ ) outputRepr[k] = inputRepr[sizeof(T) - 1 - k];
    return swapped;

} //END OF FUNCTION

template <class T> void CIO::skip_n(ifstream& fileStream, int numInstances)
{
    int ii = 0;
    while ( ii < numInstances )
    {
        fileStream.seekg(sizeof(T), ios_base::cur);
        ii++;
    }
    return;
} //END OF FUNCTION

/*
 * Function that writes any grid-based map in cube format.
 * the map must be supplied as a vector which stores data 
 * as instances of class or tyepname T.
 */
template <class T> void CIO::writeMapInCube(const delphi_integer& iGrid,
                                            const delphi_real& fScale,
                                            const SGrid<delphi_real>& fgBoxCenter,
                                            const vector<T>& vctMap,
                                            const string& strOutFile)
{

  int ix, iy, iz, iw, l;

  delphi_real coeff=0.5291772108;
  delphi_real stepsize=1.0/fScale;
  SGrid<delphi_real> origin = (fgBoxCenter-stepsize*(iGrid-1)/2.0)/coeff;

  ofstream outfile(strOutFile.c_str());

  outfile     << setw(10) << fixed << setprecision(6) << fScale << setw(6) << iGrid
              << setw(10) << setprecision(6) << fgBoxCenter.nX
              << setw(10) << setprecision(6) << fgBoxCenter.nY
              << setw(10) << setprecision(6) << fgBoxCenter.nZ
              << endl;
  outfile << "Gaussian cube format map generated by Delphi" << endl;

  outfile << fixed   << setprecision(6);
  outfile << setw(5) << right << 1     << setw(14) << right << origin.nX       << setw(14) << right          << origin.nY << setw(14) << right << origin.nZ << endl;
  outfile << setw(5) << right << iGrid << setw(14) << right << stepsize/coeff  << setw(14) << right          << 0.0       << setw(14) << right << 0.0       << endl;
  outfile << setw(5) << right << iGrid << setw(14) << right << 0.0 << setw(14) << right    << stepsize/coeff << setw(14)  << right    << 0.0   << endl;
  outfile << setw(5) << right << iGrid << setw(14) << right << 0.0 << setw(14) << right    << 0.0            << setw(14)  << right    << stepsize/coeff     << endl;
  outfile << setw(5) << right << 1     << setw(14) << right << 0.0 << setw(14) << right    << 0.0            << setw(14)  << right    << 0.0   << setw(14)  << right << 0.0 << endl;

  outfile.unsetf(ios_base::floatfield); // return to outfile default notation

  outfile.precision(5);


  for (ix = 1; ix <= iGrid; ix++)
  {
    for (iy = 1; iy <= iGrid; iy++)
    {
      l = 0;
      for (iz = 1; iz <= iGrid; iz++)
      {
        iw = (iz-1)*iGrid*iGrid + (iy-1)*iGrid + (ix-1);

        outfile << setw(13) << scientific << (float)vctMap[iw] << " " ;
        l++;
        if(l == 6)
        {
          outfile << endl;
          l=0;
        }

      }
      outfile << endl;
    }
  }

  outfile.close();

}
#endif // CIO_H_
