/*
 * io_misc.cpp
 *
 *  Created on: Feb 4, 2014
 *      Author: chuan
 */

#include "io.h"
#ifdef ENABLE_TRAJECTORY_NC
#include "io_ncamber.h"
//#include "io_ncbase.h"
#endif

//-----------------------------------------------------------------------//
void CIO::remove_forbidden_chars(string& str)
{
    string::iterator it = str.begin();
    while ( it != str.end() )
    {
        int ascii = (int) *it;
        if ( (ascii < 33) || (ascii > 126) ) *it = NULL;
        it++;
    }
    return;
}

//-----------------------------------------------------------------------//
string CIO::decode_gro_string(ifstream& fileStream)
{
    int strSize = this->decode<int>(fileStream);
    
    // learnt aboiut this from unpack_fstring(n) in https://github.com/python/cpython/blob/master/Lib/xdrlib.py
    strSize = 4 * (int) ((strSize + 3)/4);

    char outChar[strSize];
    fileStream.read(reinterpret_cast<char *>(&outChar), strSize);
    
    string outStr(outChar);
    remove_forbidden_chars(outStr);
    
    return outStr;
} 

//-----------------------------------------------------------------------//
bool CIO::shouldRecordFrame(int& iFrameCurrent, const int& iFrameFirst, const int& iFrameLast, const int& iFrameStride)
{
    bool retFlag = ( (iFrameCurrent - iFrameFirst) % iFrameStride == 0 && iFrameCurrent >= iFrameFirst && iFrameCurrent <= iFrameLast );
    return retFlag;
}

//-----------------------------------------------------------------------//
string CIO::toUpperCase(const string& strLine)
{
   const string strAlphabet = "abcdefghijklmnopqrstuvwxyz";
   string strNewLine;
   locale loc;
   size_t found;
   
   for (size_t i = 0; i < strLine.size(); i++)
   {
      found = strAlphabet.find_first_of(strLine[i]);
      
      if (string::npos != found)
         strNewLine += toupper(strLine[i],loc);
      else
         strNewLine += strLine[i];   
   }
  
   return strNewLine;
}

//-----------------------------------------------------------------------//
string CIO::removeSpace(const string& strLine)
{
    string strNewLine;
    
    for (size_t i = 0; i < strLine.size(); i++)
       if (' ' != strLine[i]) strNewLine += strLine[i];
    
    return strNewLine;
}

//-- 2016_05_08 Lin Li: this function doesn't work well. This function is unneccessary -----------------------------------//
bool CIO::checkFileFormat(const string& strFile)
{
   bool bIsFormatted = true; // initially assume formatted file

   ifstream ifFileHandle;

   ifFileHandle.open(strFile.c_str());
    
   if (!ifFileHandle.is_open()) throw CUnknownIOFile(strFile);   
   
   char cTestChar[80];  
   
   ifFileHandle.read(cTestChar,80);

   string strASCI = "1234567890 .-+#,$asdfghjklzxcvbnmqwertyuiopASDFGHJKLZXCVBNMQWERTYUIOP)(}{][/";
   int iCount = 0;
   
   for (int i = 0; i < 80; i++)
      if (string::npos == strASCI.find(cTestChar[i]) ) iCount += 1;

   if (10 < iCount) bIsFormatted = false; // unformatted file
   
   ifFileHandle.close();   
   
   return bIsFormatted;
}

//-----------------------------------------------------------------------//
void CIO::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)
{
    // perform traj-topol compatibility checks
    if ( iTrajFormatIn == DCDTRAJ ) 
    {
        if ( !(iTopolFormatIn == PSFTOP) )
            throw CBadTrajTop(iTrajFormatIn, iTopolFormatIn);
    }
    else if (iTrajFormatIn == TRRTRAJ)
    {
        if ( iTopolFormatIn != TPRTOP )
            throw CBadTrajTop(iTrajFormatIn, iTopolFormatIn);
    }
	else if (iTrajFormatIn == MDCRDTRAJ) {
		if (iTopolFormatIn != PRMTOP)
			throw CBadTrajTop(iTrajFormatIn, iTopolFormatIn);
	} 
#ifdef ENABLE_TRAJECTORY_NC
   else if (iTrajFormatIn == NCTRAJ) {
		if (iTopolFormatIn != PRMTOP)
			throw CBadTrajTop(iTrajFormatIn, iTopolFormatIn);
	}
#endif


    // ensure that last frame is always greater than or equal to the first frame 
    if ( iFrameLast < iFrameFirst ) throw CBadFrameBounds(iFrameFirst, iFrameLast);

    int natomsTraj = -1;
    int natomsTopol = -1;
    int    iFound  = -1;
    delphi_real   fValue  = 0.0; // radius or charge
    delphi_real   fSigmaLJ = 0., fEpsilonLJ = 0., fVdwGamma = 0.;  // paramters for non-polar energy calculations
    string strAtInf,strAtom,strResidue,strResidueNum,strChain, strSub0, strSub1; 
    vector <delphi_real> vctfRadii;
    
	CIO::CParmtopReadHead prmtopReadResult;
    
	/*
     * -------------- TOPOLOGY -------------
     * READ ATOM INFORMATION FROM VARIOUS FORMATS OF TOPOLOGY FILES
     */
    // read PSF if you must
    if ( iTopolFormatIn == PSFTOP )
    {
        ifstream psfFileStream;
        psfFileStream.open(strTopolFile.c_str());

        natomsTopol = this->readPSF(psfFileStream, vctSimAtoms);
        psfFileStream.close();

    }
    
    // read TPR if you must
    else if ( iTopolFormatIn == TPRTOP )
    {
        ifstream tprFileStream;
        tprFileStream.open(strTopolFile.c_str());

        natomsTopol = this->readTPR(tprFileStream, vctSimAtoms);
        tprFileStream.close();

    }
    // read PRMTOP if you must
    else if (iTopolFormatIn == PRMTOP) {
	ifstream prmtopFileStream;
	prmtopFileStream.open(strTopolFile.c_str());

	prmtopReadResult = this->readPRMTOP(prmtopFileStream, vctSimAtoms);
	natomsTopol = (prmtopReadResult.readSuccess) ? vctSimAtoms.size() : 0;
	prmtopFileStream.close();
    }

    /*
     * read SIZ file
     * CHARGE VALUES will be read from the MD topology
     */
    this->readForceFile(strSizeFile); // read delphi size file
    // if (bSolvePB) this->readForceFile(strCrgFile);  // read delphi charge file

    if (!bExistRadiiInfo) // true for all types of traj-topol combo 
    {
        vctfRadii.resize(natomsTopol);
        for (delphi_integer iThisAtom = 0; iThisAtom < natomsTopol; iThisAtom++)
        {

            strAtInf = vctSimAtoms[iThisAtom].getAtInf(); strAtInf = toUpperCase(strAtInf);

            strAtom       = strAtInf.substr(0,5);  strAtom       = removeSpace(strAtom);
            strResidue    = strAtInf.substr(6,3);  strResidue    = removeSpace(strResidue);
            strChain      = strAtInf.substr(10,1); strChain      = removeSpace(strChain);
            strResidueNum = strAtInf.substr(11,4); strResidueNum = removeSpace(strResidueNum);
            
            // assign radius, searching for decreasingly specific specification ending with generic atom type
            // note all atoms must have an assignment         
            fValue = 0.0;

            iFound = FindRecord(strAtom,strResidue,strResidueNum,strChain,SIZEFILE,fValue);

            if (-1 == iFound)
            {
                CUnknownRadius warning(strAtInf);
            }

            //strSub0 = strAtom.substr(0,1); strSub1 = strAtom.substr(1,1);
            //if (1.0e-6 > fValue && 0 != strSub0.compare("H") &&  0 != strSub1.compare("H"))
            //    CZeroHeavyAtomRadius warning(strWarn);

            vctSimAtoms[iThisAtom].setRadius(fValue);
            vctfRadii[iThisAtom] = fValue;

        } //---------- end of loop over iThisAtom =0:iAtomNum-1 
    } //---------- end of if (!bExistRadiiInfo)
    

    /*
     * -------------- TRAJECTORY -------------
     */
    // read DCD if you must
    if ( iTrajFormatIn == DCDTRAJ )
    {
        ifstream dcdFileStream;
        dcdFileStream.open(strTrajFile.c_str(), ios::binary);

        natomsTraj = this->readDCD( dcdFileStream, vctFrames, vctfRadii, iFrameFirst, iFrameLast, iFrameStride );
        dcdFileStream.close();

    }

    // read TRR if you must
    else if ( iTrajFormatIn == TRRTRAJ )
    {
        ifstream trrFileStream;
        trrFileStream.open(strTrajFile.c_str(), ios::binary);

        natomsTraj = this->readTRR( trrFileStream, vctFrames, vctfRadii, iFrameFirst, iFrameLast, iFrameStride );
        trrFileStream.close();

    }

    // read MDCRD if you must
    else if (iTrajFormatIn == MDCRDTRAJ) 
    {
		ifstream mdcrdFileStream;
		mdcrdFileStream.open(strTrajFile.c_str());

		natomsTraj = this->readMDCRD(mdcrdFileStream, prmtopReadResult.prmtopHead, vctFrames, vctfRadii, iFrameFirst, iFrameLast, iFrameStride);
		mdcrdFileStream.close();
     }
    #ifdef ENABLE_TRAJECTORY_NC
    // read NETCDF if you must
    else if (iTrajFormatIn == NCTRAJ) 
    {
        CNCAmberTraj nctraj(strTrajFile.c_str());
        natomsTraj = nctraj.readAmberNC(prmtopReadResult.prmtopHead, vctFrames, vctfRadii, iFrameFirst, iFrameLast, iFrameStride);
        nctraj.NC_close();
    }
    #endif

    // CHECK IF IDENTICAL NUMBER OF ATOMS ARE READ FROM THE TOPOL AND THE TRAJECTIRY FILE //
    if (natomsTraj == natomsTopol)
    {
        iAtomNum = natomsTraj;
        #ifdef VERBOSE
        cout << " iomisc> atoms in traj = " << natomsTraj << " and atoms in top = " << natomsTopol << endl;
        #endif
    }
    else
        throw CAtomNumberMismatch(natomsTraj, natomsTopol); 


   if (bNonPolEng || bLJEng) 
   {
      this->readVdwFile(strVdwFile); // read delphi vdw param file
      for (delphi_integer iThisAtom = 0; iThisAtom < iAtomNum; iThisAtom++)
      {
         strAtInf = vctSimAtoms[iThisAtom].getAtInf(); strAtInf = toUpperCase(strAtInf);
         strAtom  = strAtInf.substr(0,5);  strAtom       = removeSpace(strAtom);

         // assign LJ and gamma parameters
         // note all atoms must have an assignment         
         iFound = FindVdwRecord(strAtom,fSigmaLJ, fEpsilonLJ, fVdwGamma);
      
         if (-1 == iFound)
         {
              CUnknownLJParam warning (strAtInf);   
         }

         vctSimAtoms[iThisAtom].setSigmaLJ(fSigmaLJ);
         vctSimAtoms[iThisAtom].setEpsilonLJ(fEpsilonLJ);
         vctSimAtoms[iThisAtom].setVdwGamma(fVdwGamma);
         
      } //---------- end of loop over iThisAtom =0:iAtomNum-1 
   } //---------- end of if (bNonPolEng)
   

    #ifdef VERBOSE
    for ( auto& satom : this->vctSimAtoms ) 
    {
        cout << left << setw(6) << satom.getAtomName() << " of " ;
        cout << left << setw(4) << satom.getAtomResname(); 
        cout << left << setw(6) << satom.getAtomResid();
        cout << " | ";
        cout << " Radius = " << right << setw(10) << satom.getRadius() << ", ";
        cout << " Charge = " << right << setw(10) << satom.getCharge() << ", ";
        if ( bNonPolEng || bLJEng )
	{
		cout << " EpsLJ = " << right << setw(10) << satom.getEpsilonLJ() << ", ";
		cout << " SigLJ = " << right << setw(10) << satom.getEpsilonLJ() << ", ";
		cout << " GamLJ = " << right << setw(10) << satom.getVdwGamma() << ", ";
	}
	cout << endl;
    }
    #endif

   return;

}// END OF FUNCTION


//-----------------------------------------------------------------------//
void CIO::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)
{
   int    iFound  = -1;
   delphi_real   fValue  = 0.0; // radius or charge
   delphi_real   fSigmaLJ = 0., fEpsilonLJ = 0., fVdwGamma = 0.;  // paramters for non-polar energy calculations
   string strAtInf,strSub0,strSub1,strAtom,strResidue,strResidueNum,strChain; 

#ifdef VERBOSE
   cout << "\nassigning charges and radii... \n\n";
#endif
   /*
    * only standard PDB requires to read charge and size from separate files. Otherwise, charges and sizes are already
    * included in the input PDB file.
    */
   if (STDPDB == iPdbFormat || GROPDB == iPdbFormat)
   {
      this->readForceFile(strSizeFile); // read delphi size file
      if (bSolvePB) this->readForceFile(strCrgFile);  // read delphi charge file
      if (bMultiSigmaGauss) this->readForceFile(strSggFile); // read sigma-gausian file
   }
   else
   {
      CReadSizeCrgFromPDB warning;
   }
    
   this->readPdbFile(strPdbFile,iPdbFormat,bPdbUnformat, strCommPDB); // read pdb file

   if (!bExistRadiiInfo) // pdb file, such as std pdb or gromacs' gro, does not have radii info
   {
      for (delphi_integer iThisAtom = 0; iThisAtom < iAtomNum; iThisAtom++)
      {
         strAtInf = vctapAtomPdb[iThisAtom].getAtInf(); strAtInf = toUpperCase(strAtInf);

         /*
          * for certain pdb files such as 1AB1.pdb, position 5 (maybe and 9) are for special purpose
          */
         //strAtom       = strAtInf.substr(0,6);  strAtom       = removeSpace(strAtom);
         //strResidue    = strAtInf.substr(6,4);  strResidue    = removeSpace(strResidue);
         //strChain      = strAtInf.substr(10,1); strChain      = removeSpace(strChain);
         //strResidueNum = strAtInf.substr(11,4); strResidueNum = removeSpace(strResidueNum);

         strAtom       = strAtInf.substr(0,5);  strAtom       = removeSpace(strAtom);
         strResidue    = strAtInf.substr(6,3);  strResidue    = removeSpace(strResidue);
         strChain      = strAtInf.substr(10,1); strChain      = removeSpace(strChain);
         strResidueNum = strAtInf.substr(11,4); strResidueNum = removeSpace(strResidueNum);

         // assign radius, searching for decreasingly specific specification ending with generic atom type
         // note all atoms must have an assignment         
         fValue = 0.0;

         iFound = FindRecord(strAtom,strResidue,strResidueNum,strChain,SIZEFILE,fValue);
      
         if (-1 == iFound)
         {
            //throw  CUnknownRadius(strAtInf); // need stop here
            CUnknownRadius warning(strAtInf);
         }

         strSub0 = strAtom.substr(0,1); strSub1 = strAtom.substr(1,1);
         if (1.0e-6 > fValue && 0 != strSub0.compare("H") &&  0 != strSub1.compare("H"))
            CZeroHeavyAtomRadius warning(strAtInf);
      
         vctapAtomPdb[iThisAtom].setRadius(fValue);
         
         fValue = 0.0;
         
         if (bSolvePB)
         {
            iFound = FindRecord(strAtom,strResidue,strResidueNum,strChain,CHARGEFILE,fValue);
            if (-1 == iFound) CUnknownCharge warning(strAtInf); // need stop here
            vctapAtomPdb[iThisAtom].setCharge(fValue);
         }
      } //---------- end of loop over iThisAtom =0:iAtomNum-1 
   } //---------- end of if (!bExistRadiiInfo)
   if(bUseGaussian)
   {
      for (delphi_integer iThisAtom = 0; iThisAtom < iAtomNum; iThisAtom++)
      {
         
         strAtInf = vctapAtomPdb[iThisAtom].getAtInf(); strAtInf = toUpperCase(strAtInf);

         /*
          * for certain pdb files such as 1AB1.pdb, position 5 (maybe and 9) are for special purpose
          */

         strAtom       = strAtInf.substr(0,5);  strAtom       = removeSpace(strAtom);
         strResidue    = strAtInf.substr(6,3);  strResidue    = removeSpace(strResidue);
         strChain      = strAtInf.substr(10,1); strChain      = removeSpace(strChain);
         strResidueNum = strAtInf.substr(11,4); strResidueNum = removeSpace(strResidueNum);
         
         fValue = 0.0;
         if (bMultiSigmaGauss)
         {
            iFound = FindSggRecord(strAtom,strResidue,strResidueNum,strChain,SIGMAGAUSSFILE,fValue);
            if (-1 == iFound) 
	    {
	        vctapAtomPdb[iThisAtom].setSigmaGaussian(fSigma);
                CUnknownSigmaGaussian warning(strAtInf, fSigma); // need stop here
	    }
            else 
            {
                vctapAtomPdb[iThisAtom].setSigmaGaussian(fValue);
            }
         } else {
           vctapAtomPdb[iThisAtom].setSigmaGaussian(fSigma);
         }
      }
   }

   // check charge 
   // inc_setrcmod_chkcrg but combined the two cases of isitsf = .true. and .false.
   string strLastResidue    = "    "; // a4
   string strLastResidueNum = "    "; // a4
   delphi_real   fCrgSum,fError; // total net charge

   fCrgSum = 0.0;
   
   for (delphi_integer iThisAtom = 0; iThisAtom < iAtomNum; iThisAtom++)
   {
      fValue = vctapAtomPdb[iThisAtom].getCharge();
         
      strAtInf      = vctapAtomPdb[iThisAtom].getAtInf();
      strResidue    = strAtInf.substr(6,4); 
      strResidueNum = strAtInf.substr(11,4);
         
      if ( 0 != strResidue.compare(strLastResidue) || 0 != strResidueNum.compare(strLastResidueNum) )
      {
         if (1.0e-4 < abs(fCrgSum))
         {
            fError = abs(fCrgSum) - 1.0;
               
            if (1.0e-4 < abs(fError)) 
               CNonZeroNetCrg waring(strLastResidue,strLastResidueNum,fCrgSum);
         }
         
         strLastResidue    = strResidue;
         strLastResidueNum = strResidueNum;
         fCrgSum           = fValue; 
      }
      else
         fCrgSum += fValue;
          

      // specific section to find maximum residue number. Goes here when isitsf = .true.
      if (bSurfCrgInSite)
      {
         delphi_integer iNumber =  atoi( strResidueNum.c_str() );
         if (iNumber > iResidueNum) iResidueNum = iNumber;
      }
   }
  
   if (1.0e-4 < abs(fCrgSum))
   {
      fError = abs(fCrgSum) - 1.0;
      if (1.0e-4 < abs(fError))
         CNonZeroNetCrg waring(strLastResidue,strLastResidueNum,fCrgSum);
   }  

   if (bNonPolEng || bLJEng) 
   {
      this->readVdwFile(strVdwFile); // read delphi vdw param file
      for (delphi_integer iThisAtom = 0; iThisAtom < iAtomNum; iThisAtom++)
      {
         strAtInf = vctapAtomPdb[iThisAtom].getAtInf(); strAtInf = toUpperCase(strAtInf);
         strAtom  = strAtInf.substr(0,5);  strAtom       = removeSpace(strAtom);

         // assign LJ and gamma parameters
         // note all atoms must have an assignment         
         iFound = FindVdwRecord(strAtom,fSigmaLJ, fEpsilonLJ, fVdwGamma);
      
         if (-1 == iFound)
         {
              CUnknownLJParam warning (strAtInf);   
         }

         vctapAtomPdb[iThisAtom].setSigmaLJ(fSigmaLJ);
         vctapAtomPdb[iThisAtom].setEpsilonLJ(fEpsilonLJ);
         vctapAtomPdb[iThisAtom].setVdwGamma(fVdwGamma);
         
      } //---------- end of loop over iThisAtom =0:iAtomNum-1 
   } //---------- end of if (bNonPolEng)
#ifdef DEBUG_IO_PDB
   printPDB();
#endif

}                        
















