/**
 * @file app_delphi.cpp
 * @brief Main function to generate the executable delphicpp
 *
 * @author Chuan Li, chuanli@clemson.edu
 *
 * delphicpp is the object-oriented C++ version of the program DelPhi, which takes as input a Brookhaven
 * database coordinate file format of a molecule or equivalent data for geometrical objects and/or charge
 * distributions and calculates the electrostatic potential in and around the system, using a finite
 * difference solution to the Poisson-Boltzmann equation. This can be done for any given concentration of
 * any two different salts. The system can be composed of different parts having different dielectrics.
 *
 * =================================================================================================
 *
 * [QUICK REFERENCE TO ACCESS SVN REPOSITORY ON SERVER consus.clemson.edu] \n
 * - TO SHOW LOG HISTORY FOR THIS PROJECT:
 *       svn log svn+ssh://username\@consus.clemson.edu/home/common/delphi_distr/delphicpp
 * - TO CHECKOUT THE LATEST REVERSION:
 *       svn co svn+ssh://username\@consus.clemson.edu/home/common/delphi_distr/delphicpp
 * - TO REMOVE/ADD FILES FROM A CHECKED-OUT PROJECT:
 *       svn rm  <file/directory name>
 *       svn add <file/directory name>
 * - TO REVIEW LOCAL CHANGES:
 *       svn status
 * - TO COMMIT CHANGES:
 *       export SVN_EDITOR=/usr/bin/vim
 *       svn commit
 * - TO UPDATE YOUR LOCAL FILES:
 *       svn update
 *
 *   (See "subversion user guide" created by S. Sarkar and C. Li for other svn options)
 *
 * =================================================================================================
 *
 **/

#include "../delphi/delphi_data.h"
#include "../delphi/delphi_datamarshal.h"
#include "../space/space.h"
#include "../solver/solver.h"
#include "../energy/energy.h"
#include "../site/site.h"

#include <time.h>
#include <math.h>
#include <iostream>
#include <cstddef>
#include <sstream>
#include <string>

#ifdef PARALLEL_OMP
#include <omp.h>
#endif

using namespace std;


int main(int argc, char *argv[])
{
    /*
     * bool values are inserted/extracted by their textual representation: either true or false, instead of integral values.
     * This flag can be unset with the noboolalpha manipulator.
     */


    size_t MAXWIDTH = 45;
    string info_string = " Info>";

#ifdef VERBOSE
    cout << boolalpha;
    cerr << boolalpha;
#endif

#ifdef DEVELOPER
    cout << fixed << setprecision(7); //cout.precision(15)
#else
    cout << fixed << setprecision(3); //cout.precision(7)
#endif

    try
    {
        CTimer * pTester = new CTimer; // record execution time

        pTester->start();

        shared_ptr<CTimer> pTimer(new CTimer); // record execution time

#ifndef PARALLEL_OMP
        cout << "\n             DelPhi is running on (single CPU) version" << endl;
#else
        int omp_num_threads, omp_thread_id;
        int num_threads_to_use = 0;
        const char* env_num_threads = getenv("OMP_NUM_THREADS");

        std::cout << "DEBUG: OMP_NUM_THREADS from getenv(): ";
        if (env_num_threads != nullptr) {
            std::cout << "'" << env_num_threads << "'" << std::endl;
            num_threads_to_use = atoi(env_num_threads);
        } else {
            std::cout << "(not set)" << std::endl;
        }
        std::cout << "DEBUG: Integer value parsed: " << num_threads_to_use << std::endl;


        if (num_threads_to_use > 0) {
            omp_set_num_threads(num_threads_to_use);
        } else {
            // Fallback: Use the number of available processors
            omp_set_num_threads(omp_get_num_procs());
        }
        omp_thread_id = omp_get_max_threads(); // this retreives the number of threads to use set earlier.
        if (0 == omp_thread_id)
        {
            omp_num_threads = omp_get_num_threads();
            cout << "\n       DelPhi is running on (OMP) version with " << omp_num_threads << " thread(s)" << endl;
        }
#endif

        //********************************************************************************//
        //                                                                                //
        //                  read the parameter, size, charge, PDB files                   //
        //                                                                                //
        //********************************************************************************//

        //cout << "DPLT IDELPHI Creating" << endl;
        //---------- a shared_ptr to an object of IDataContainer
        shared_ptr<IDataContainer> pDataContainer(new CDelphiData(argc, argv, pTimer));
        //cout << "DPLT IDELPHI Created" << endl;

#ifdef DEBUG_DELPHI_SPACE
        cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_SPACE] data container is written into file test_delphicpp_atbeginning.dat---\n\n"
            << "\033[0m";

        pDataContainer->showMap("test_delphicpp_atbeginning.dat");
#endif

        //pDataContainer->showMap("test_delphicpp_before_space.dat");

#ifdef DEVELOPER
        cout << "\n\n---------- delphicpp finishes IO on ";
        pTester->showTime();
        cout << "\n\n";
#endif

        //***************************** TO RUN OR NOT RUN MMPBSA *************************//
        bool& bDoMMPBSA(pDataContainer->getKey_Ref<bool>("doMMPBSA"));
        // If yes: then [space -> solver -> energy -> site] will be run for all frames
        // if no : then [space -> solver -> energy -> site] will be run for that one input structure (standard)


        // ************************* PARAMS FOR GAUSSIAN RUN ****************************//
        int& inhomo(pDataContainer->getKey_Ref<int>("inhomo"));
        const int& iGaussian(pDataContainer->getKey_constRef<int>("gaussian"));
        bool& logs(pDataContainer->getKey_Ref<bool>("logs"));
        bool& logg(pDataContainer->getKey_Ref<bool>("logg"));
        //********************************************************************************//
        //                                                                                //
        //    realize an object of CDelphiSpace class to construct molecular surfaces     //
        //                                                                                //
        //********************************************************************************//
        if (!bDoMMPBSA )
        {
            inhomo = 0;

            if (iGaussian == 1 && logs)
            {
                logg = true; //for gaussian
                inhomo = 1;
            }

            unique_ptr<CTimer> spaceTimer(new CTimer);
            spaceTimer->start();
            unique_ptr<IAbstractModule> pSpace(new CDelphiSpace(pDataContainer, pTimer));
            
            pSpace->run();
            pSpace.reset();

            spaceTimer->showElapseMessage("Total space-module runtime is");
            spaceTimer.reset();
#ifdef DEBUG_DELPHI_SPACE
            cout << "---[DEBUG_DELPHI_SPACE] data container is written into file test_delphicpp_aftersurf.dat---\n\n";

            pDataContainer->showMap("test_delphicpp_aftersurf.dat");
#endif


#ifdef DEVELOPER
            cout << "\n\n---------- delphicpp finishes SPACE class on ";
            pTester->showTime();
            cout << "\n\n";
            cout << endl;
#endif

#ifdef VERBOSE
            if ( !(iGaussian == 1 && inhomo == 0 && logs) )
                cout << left << setw(MAXWIDTH) << " Number of atom coordinates read" << " : "
                    << pDataContainer->getKey_constRef< delphi_integer > ("natom") << endl;
#endif

            if (pDataContainer->getKey_constRef<bool>("isolv"))
            {
                if ( !(iGaussian == 1 && inhomo == 0 && logs) )
                {
                    cout << left << setw(MAXWIDTH) << " Total number of assigned charges" << " : "
                        << pDataContainer->getKey_constRef< delphi_integer > ("nqass") << endl;
                    cout << left << setw(MAXWIDTH) << " Net assigned charge" << " : "
                        << right << setw(10) << pDataContainer->getKey_constRef< delphi_real > ("qnet") << endl;
                    cout << left << setw(MAXWIDTH) << " Assigned positive charge" << " : "
                        << right << setw(10) << pDataContainer->getKey_constRef< delphi_real > ("qplus") << endl;
                    cout << left << setw(MAXWIDTH) << " Centred at (gu)" << " : "
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqplus").nX << " " << right << setw(10)
                        << pDataContainer->getKey_constRef< SGrid<delphi_real> > ("cqplus").nY << " " << right << setw(10)
                        << pDataContainer->getKey_constRef< SGrid<delphi_real> > ("cqplus").nZ << endl;
                    cout << left << setw(MAXWIDTH) << " Assigned negative charge"
                        << " : " << pDataContainer->getKey_constRef< delphi_real > ("qmin") << endl;
                    cout << left << setw(MAXWIDTH) << " Centred at (gu)" << " : "
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqmin").nX << " " << right << setw(10)
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqmin").nY << " " << right << setw(10)
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqmin").nZ << endl;

#ifdef VERBOSE
                    cout << left << setw(MAXWIDTH) << " Number of dielectric boundary points" << " : "
                        << pDataContainer->getKey_constRef < delphi_integer > ("ibnum") << endl;
#endif

                    cout << endl;
                }

                if (pDataContainer->getKey_constRef<bool>("iexun") && 0 == pDataContainer->getKey_constRef < delphi_integer > ("ibnum"))
                    throw CNoBndyAndDielec(pTimer);

                //********************************************************************************//
                //                                                                                //
                //   realize an object of CDelphiSolver class to calculate potentials on grids    //
                //                                                                                //
                //********************************************************************************//

#ifdef DEBUG_DELPHI_SOLVER
                cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_SOLVER] data container is read from file test_fromsurf.dat---\n\n" << "\033[0m";

                pDataContainer->reset("test_fromsurf.dat");

                cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_SOLVER] data container is written into file test_delphicpp_beforeitr.dat---\n\n"
                    << "\033[0m";

                pDataContainer->showMap("test_delphicpp_beforeitr.dat");

#endif // DEBUG_DELPHI_SOLVER
            }

            unique_ptr<CTimer> solverTimer(new CTimer);
            solverTimer->start();
            unique_ptr<IAbstractModule> pSolver(new CDelphiSolver(pDataContainer, pTimer));

            pSolver->run();
            pSolver.reset();

            solverTimer->showElapseMessage("Total solver-module runtime is");
            solverTimer.reset();
#ifdef DEBUG_DELPHI_SOLVER
            cout << "---[DEBUG_DELPHI_SOLVER] data container is written into file test_delphicpp_afteritr.dat---\n\n";

            pDataContainer->showMap("test_delphicpp_afteritr.dat");
#endif // DEBUG_DELPHI_SOLVER

#ifdef DEVELOPER
            cout << "\n\n---------- delphicpp finishes SOLVER class on ";
            pTester->showTime();
            cout << "\n\n";
#endif

            //********************************************************************************//
            //                                                                                //
            //          realize an object of CDelphiEnergy class to calculate energies        //
            //                                                                                //
            //********************************************************************************//

#ifdef DEBUG_DELPHI_ENERGY
            cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_ENERGY] data container is read from file test_fromsurf.dat "
                << "and test_fromsolver.dat---\n\n" << "\033[0m";
            pDataContainer->reset("test_fromsurf.dat");
            pDataContainer->reset("test_fromsolver.dat");
            cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_SOLVER] data container is written into file test_delphicpp_beforeenergy.dat---\n\n"
                << "\033[0m";
            pDataContainer->showMap("test_delphicpp_beforeenergy.dat");
#endif // DEBUG_DELPHI_ENERGY

            unique_ptr<CTimer> energyTimer(new CTimer);
            energyTimer->start();

            unique_ptr<IAbstractModule> pEnergy(new CDelphiEnergy(pDataContainer, pTimer));
            pEnergy->run();
            pEnergy.reset();
            
            energyTimer->showElapseMessage("Total energy-module runtime is");
            energyTimer.reset();
#ifdef DEBUG_DELPHI_ENERGY
            cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_ENERGY] data container is written into file test_delphicpp_aftereng.dat---\n\n"
                << "\033[0m";
            pDataContainer->showMap("test_delphicpp_aftereng.dat");
#endif // DEBUG_DELPHI_ENERGY

#ifdef DEVELOPER
            cout << "\n\n---------- delphicpp finishes ENERGY class on ";
            pTester->showTime();
            cout << "\n\n";
#endif

            if (iGaussian == 1 && inhomo == 1 && logs) //second run for Gaussian
            {

                //pDataContainer.reset();

                //shared_ptr<IDataContainer> pDataContainer( new CDelphiData(argc,argv,pTimer) );
                inhomo = 0;

                unique_ptr < IAbstractModule > pSpace(new CDelphiSpace(pDataContainer, pTimer));
                pSpace->run();
                pSpace.reset();

#ifdef DEBUG_DELPHI_SPACE
                cout << "---[DEBUG_DELPHI_SPACE] data container is written into file test_delphicpp_aftersurf.dat---\n\n";
                pDataContainer->showMap("test_delphicpp_aftersurf.dat");
#endif

                unique_ptr < IAbstractModule > pSolver(new CDelphiSolver(pDataContainer, pTimer));
                pSolver->run();
                pSolver.reset();

#ifdef DEBUG_DELPHI_SOLVER
                cout << "---[DEBUG_DELPHI_SOLVER] data container is written into file test_delphicpp_afteritr.dat---\n\n";

                pDataContainer->showMap("test_delphicpp_afteritr.dat");
#endif
                unique_ptr < IAbstractModule > pEnergy(new CDelphiEnergy(pDataContainer, pTimer));
                pEnergy->run();
                pEnergy.reset();
            }

            //********************************************************************************//
            //                                                                                //
            //               realize an object of CSite class to write site info              //
            //                                                                                //
            //********************************************************************************//

#ifdef DEBUG_DELPHI_SITE
            cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_SITE] data container is read from file test_fromsurf.dat, "
                << "test_fromsolver.dat and test_fromenergy.dat---\n\n" << "\033[0m";
            pDataContainer->reset("test_fromsurf.dat");
            pDataContainer->reset("test_fromsolver.dat");
            pDataContainer->reset("test_fromenergy.dat");
            cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_SOLVER] data container is written into file test_delphicpp_beforesite.dat---\n\n"
                << "\033[0m";
            pDataContainer->showMap("test_delphicpp_beforesite.dat");
#endif // DEBUG_DELPHI_SITE

            unique_ptr<CSite> pSite(new CSite(pDataContainer, pTimer));

            if (pDataContainer->getKey_constRef<bool>("isite"))
            {
                int iisitsf = 0;
                if (pDataContainer->getKey_Ref<bool>("isitsf")) iisitsf = 1;
                pSite->writeSite(iisitsf);
            }

            if (pDataContainer->getKey_constRef<bool>("phiwrt")) pSite->writePhi();
            if (pDataContainer->getKey_constRef < delphi_real > ("radipz") > 0.) pSite->wirtePAnalysis();

            pSite.reset();

#ifdef DEBUG_DELPHI_SITE
            cerr << "\n\033[1;35m" << "---[DEBUG_DELPHI_SITE] data container is written into file test_delphicpp_atend.dat---\n\n"
                << "\033[0m";
            pDataContainer->showMap("test_delphicpp_atend.dat");
#endif // DEBUG_DELPHI_SITE

#ifdef DEVELOPER
            cout << "\n\n---------- delphicpp finishes SITE class on ";
            pTester->showTime();
            cout << "\n\n";
#endif

            pDataContainer.reset();

            pTimer->exit();
            pTimer.reset();

            delete pTester;
        } 
        else if ( bDoMMPBSA )
        {
            /** 
             * ---------------------------------------------
             * setup space variables to provide framewise 
             * ---------------------------------------------
             *  igrid
             *  perfil
             *  scale
             *  Min coordinate
             *  Max Coordinate
             *  BoxCenter
             *  Coordinate range
             *  geom center
             *  delphipdb atoms
             */

            delphi_integer&         iGrid1D(pDataContainer->getKey_Ref<int>("igrid"));
            delphi_real&            fPerfil(pDataContainer->getKey_Ref<delphi_real>("perfil"));
            delphi_real&            fScale(pDataContainer->getKey_Ref<delphi_real>("scale"));
            SGrid<delphi_real>&     gfMinCoordinate( pDataContainer->getKey_Ref< SGrid<delphi_real> >("cmin"));
            SGrid<delphi_real>&     gfMaxCoordinate( pDataContainer->getKey_Ref< SGrid<delphi_real> >("cmax"));
            SGrid<delphi_real>&     gfBoxCenter( pDataContainer->getKey_Ref< SGrid<delphi_real> >("oldmid"));
            SGrid<delphi_real>&     gfCoordinateRange( pDataContainer->getKey_Ref< SGrid<delphi_real> >("cran"));
            SGrid<delphi_real>&     gfGeometricCenter( pDataContainer->getKey_Ref< SGrid<delphi_real> >("pmid"));
            vector<CAtomPdb>&       vctapAtomPdb( pDataContainer->getKey_Ref< vector<CAtomPdb> >("delphipdb") );

            // angstromg and grid coord
            vector< SGrid<delphi_real> >& vctgfAtomCoordA(pDataContainer->getKey_Ref< vector<SGrid<delphi_real>> >("xn1"));
            vector< SGrid<delphi_real> >& vctgfAtomCoordG(pDataContainer->getKey_Ref< vector<SGrid<delphi_real>> >("xn2"));

            // extemes of the frame
            vector< SExtrema<delphi_real> >& vctefExtrema(pDataContainer->getKey_Ref< vector< SExtrema<delphi_real> > >("limobject"));

            /**
             * ---------------------------------------------
             * setup solver variables to provide framewise
             * ---------------------------------------------
             * iDielecBndySum
             * iCrgedGridSum
             * vctfGridCrg
             * vctgiGridCrgPose
             * iCrgBdyGrid
             * vctdgvCrgBndyGrid
             * vctfPhiMap
             * vctfPhiMap_Pre
             */
            delphi_integer&                 iDielecBndySum(pDataContainer->getKey_Ref<delphi_integer> ("icount2b"));
            delphi_integer&                 iCrgedGridSum(pDataContainer->getKey_Ref<delphi_integer> ("icount1b"));
            vector<delphi_real>&            vctfGridCrg(pDataContainer->getKey_Ref<vector<delphi_real>> ("gchrg"));
            vector<SGrid<delphi_integer>>&  vctgiGridCrgPose(pDataContainer->getKey_Ref<vector< SGrid<delphi_integer> >>("gchrgp"));
            delphi_integer&                 iCrgBdyGrid(pDataContainer->getKey_Ref<delphi_integer>("ibc"));
            vector<SDoubleGridValue>&       vctdgvCrgBndyGrid(pDataContainer->getKey_Ref<vector<SDoubleGridValue>>("cgbp"));
            vector<delphi_real>&            vctfPhiMap(pDataContainer->getKey_Ref<vector<delphi_real>>("phimap"));
            vector<delphi_real>&            vctfPhiMap_Pre(pDataContainer->getKey_Ref<vector<delphi_real>>("phimap_pre"));

            /*
             * ---------------------------------------------
             * simulation atoms and trajectory
             * ---------------------------------------------
             */
            vector<CFrame>&    vctSimFrames(pDataContainer->getKey_Ref< vector<CFrame> >("snapshots"));    
            vector<CSimAtom>&  vctSimAtoms(pDataContainer->getKey_Ref< vector<CSimAtom> >("trajatoms"));    

            /*
             * -----------------------------------------------------------
             * setup energy variables to determine 
             *  - which energies to print
             *  - what those energies are
             * -----------------------------------------------------------
             */
            bool&   bGridEng(pDataContainer->getKey_Ref<bool>("logg")); 
            bool&   bSolvEng(pDataContainer->getKey_Ref<bool>("logs"));
            bool&   bIonsEng(pDataContainer->getKey_Ref<bool>("logions")); 
            bool&   bCoulombEng(pDataContainer->getKey_Ref<bool>("logc"));
            bool&   bLJEng(pDataContainer->getKey_Ref<bool>("loglj"));
            bool&   bNonPolEng(pDataContainer->getKey_Ref<bool>("lognp"));

            delphi_real& fEngGrid(pDataContainer->getKey_Ref<delphi_real>("ergg"));
            delphi_real& fEngCorrect(pDataContainer->getKey_Ref<delphi_real>("ergs"));
            delphi_real& fEngIons(pDataContainer->getKey_Ref<delphi_real>("ergions"));
            delphi_real& fEngCoul(pDataContainer->getKey_Ref<delphi_real>("ergc"));
            delphi_real& fEngLJ(pDataContainer->getKey_Ref<delphi_real>("erglj"));
            delphi_real& fEngNonPolCav_SA(pDataContainer->getKey_Ref<delphi_real>("ergncs"));
            delphi_real& fEngNonPolCav_Vol(pDataContainer->getKey_Ref<delphi_real>("ergncv"));
            delphi_real& fEngNonPolCav(pDataContainer->getKey_Ref<delphi_real>("ergnc"));

            // prepare the header of the energy output file (default: "energy.dat")
            const string& strEnergyFile(pDataContainer->getKey_constRef < string > ("nrgnam"));
            ofstream ofEnergyFile(strEnergyFile.c_str());
            ofEnergyFile << setw(15) << left << "FRAME";
            if (bGridEng)        ofEnergyFile << setw(15) << left << "GRID";
            if (bSolvEng)        ofEnergyFile << setw(15) << left << "SOLV";
            if (bIonsEng)        ofEnergyFile << setw(15) << left << "IONIC";
            if (bCoulombEng)     ofEnergyFile << setw(15) << left << "COULOMBIC";
            if (bLJEng)          ofEnergyFile << setw(15) << left << "LJ";
            if (bNonPolEng)      ofEnergyFile << setw(15) << left << "CAVITY_SA" << setw(15) << left << "CAVITY_VOL" << setw(15) << left << "NONPOLAR";
            ofEnergyFile << endl;
            ofEnergyFile.close();

            /*
             * setup vctapAtomPdb with all the information from topology except for coordinates
             * which will be updtaesd using trajectories.
             */
            vctapAtomPdb.clear();
            vctapAtomPdb.resize(pDataContainer->getKey_Ref<delphi_integer>("natom"));

            for (int idx = 0; idx < vctapAtomPdb.size(); idx++)
            {
                CSimAtom idxAtom = vctSimAtoms[idx];

                vctapAtomPdb[idx].setAtInf(idxAtom.getAtInf());
                vctapAtomPdb[idx].setRadius(idxAtom.getRadius());
                vctapAtomPdb[idx].setCharge(idxAtom.getCharge());
                vctapAtomPdb[idx].setAtResname(idxAtom.getAtomResname());

                // if LJ or non-polar energy calculations are turned on
                vctapAtomPdb[idx].setSigmaLJ(idxAtom.getSigmaLJ());
                vctapAtomPdb[idx].setEpsilonLJ(idxAtom.getEpsilonLJ());
                vctapAtomPdb[idx].setVdwGamma(idxAtom.getVdwGamma());

            }

            // gaussian based PB
            if (iGaussian == 1 && logs)
            {
                logg = true; //for gaussian
                inhomo = 1;
            }

            /*
             * iterate over the coordinates of all the atoms provided in each frame
             * and perform space, solver, energy, site on all each of these frames.
             */
            int iFrame = 0;
            vector<CFrame>::iterator itFrame = vctSimFrames.begin();
            while ( itFrame != vctSimFrames.end() )
            {
                cout << fixed << setprecision(7); //cout.precision(15)

                cout << " ********************************************************** " << endl;
                cout << "                       Frame" << setw(5) << right << iFrame << endl;
                cout << " ********************************************************** " << endl;

#ifdef VERBOSE
                // itFrame->dumpFrame();
#endif

                iGrid1D = itFrame->iFrameGrid;
                fPerfil = itFrame->fFramePercentageFill;
                fScale  = itFrame->fFrameScale;

                gfMinCoordinate = itFrame->gfFrameMinCoordinate;
                gfMaxCoordinate = itFrame->gfFrameMaxCoordinate;
                gfCoordinateRange = itFrame->gfFrameCoordinateRange;

                gfGeometricCenter = itFrame->gfFrameMidPoint;
                gfBoxCenter = itFrame->gfFrameBoxCenter;

                vctgfAtomCoordA = itFrame->aCoordinates;
                vctgfAtomCoordG = itFrame->gCoordinates;

                /*
                 * update coordinate of each element of vctapAtomPdb 
                 * based on the information present in the current frame
                 */
                for (int idx = 0; idx < itFrame->numAtoms; idx++)
                    vctapAtomPdb[idx].setPose(itFrame->coordinates[idx].nX, itFrame->coordinates[idx].nY, itFrame->coordinates[idx].nZ);

#ifdef VERBOSE
                string strPDB = "frame" + to_string(iFrame) + ".pqr";
                ofstream ofPDB(strPDB.c_str());

                for (int idx = 0; idx < itFrame->numAtoms; idx++)
                {
                    ofPDB << setw(6)  << left << "ATOM";
                    ofPDB << setw(5)  << right  << idx+1;
                    ofPDB << setw(15) << left << vctSimAtoms[idx].getAtInf();
                    ofPDB << "    ";

                    ofPDB << fixed << setprecision(3);
                    ofPDB << setw(8) << right << vctapAtomPdb[idx].getPose().nX;
                    ofPDB << setw(8) << right << vctapAtomPdb[idx].getPose().nY;
                    ofPDB << setw(8) << right << vctapAtomPdb[idx].getPose().nZ;

                    ofPDB << fixed << setprecision(4);
                    ofPDB << setw(8) << right << vctapAtomPdb[idx].getCharge();
                    ofPDB << setw(7) << right << vctapAtomPdb[idx].getRadius();
                    ofPDB << endl;
                }
                ofPDB.close();
#endif

                inhomo = 0;
                if (iGaussian == 1 && logs)
                {
                    logg = true; //for gaussian
                    inhomo = 1;
                }

                cout << left << setw(MAXWIDTH) << " Grid size " << " : " << iGrid1D << endl;
                cout << left << setw(MAXWIDTH) << " Percentage filling " << " : "  << fPerfil << endl;
                cout << left << setw(MAXWIDTH) << " Scale " << " : " << fScale << endl;
                cout << left << setw(MAXWIDTH) << " xmin,xmax (A)"   << " : " << gfMinCoordinate.nX << " " << right << setw(10) << gfMaxCoordinate.nX << endl;
                cout << left << setw(MAXWIDTH) << " ymin,ymax (A)"   << " : " << gfMinCoordinate.nY << " " << right << setw(10) << gfMaxCoordinate.nY << endl;
                cout << left << setw(MAXWIDTH) << " zmin,zmax (A)"   << " : " << gfMinCoordinate.nZ << " " << right << setw(10) << gfMaxCoordinate.nZ << endl;
                cout << left << setw(MAXWIDTH) << " x,y,z range (A)" << " : " << gfCoordinateRange.nX << " "
                    << right << setw(10) << gfCoordinateRange.nY << " " << right << setw(10) << gfCoordinateRange.nZ << endl;
                cout << left << setw(MAXWIDTH) << " System geometric center (A)" << " : " << gfGeometricCenter.nX << " "
                    << right << setw(10) << gfGeometricCenter.nY << " "
                    << right << setw(10) << gfGeometricCenter.nZ << endl;
                cout << left << setw(MAXWIDTH) << " Grid box is centered (A)" << " : " << gfBoxCenter.nX << " "
                    << right << setw(10) << gfBoxCenter.nY << " "
                    << right << setw(10) << gfBoxCenter.nZ << endl;


                unique_ptr<IAbstractModule> pSpace(new CDelphiSpace(pDataContainer, pTimer));
                pSpace->run();
                pSpace.reset();

#ifdef DEVELOPER
                cout << "\n\n---------- delphicpp finishes SPACE class on ";
                pTester->showTime();
                cout << "\n\n";
                cout << endl;
#endif

                if ( !(iGaussian == 1 && inhomo == 0 && logs) )
                {
                    cout << left << setw(MAXWIDTH) << " Net assigned charge" << " : "
                        << right << setw(10) << pDataContainer->getKey_constRef< delphi_real > ("qnet") << endl;
                    cout << left << setw(MAXWIDTH) << " Centred at (gu)" << " : "
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqplus").nX << " " << right << setw(10)
                        << pDataContainer->getKey_constRef< SGrid<delphi_real> > ("cqplus").nY << " " << right << setw(10)
                        << pDataContainer->getKey_constRef< SGrid<delphi_real> > ("cqplus").nZ << endl;
                    cout << left << setw(MAXWIDTH) << " Centred at (gu)" << " : "
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqmin").nX << " " << right << setw(10)
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqmin").nY << " " << right << setw(10)
                        << pDataContainer->getKey_constRef < SGrid<delphi_real> > ("cqmin").nZ << endl;

                    cout << left << setw(MAXWIDTH) << " Number of dielectric boundary points" << " : "
                        << pDataContainer->getKey_constRef < delphi_integer > ("ibnum") << endl;

                    cout << endl;
                }

                if (pDataContainer->getKey_constRef<bool>("iexun") && 
                        0 == pDataContainer->getKey_constRef < delphi_integer > ("ibnum"))
                    throw CNoBndyAndDielec(pTimer);

                //********************************************************************************//
                //                                                                                //
                //   realize an object of CDelphiSolver class to calculate potentials on grids   //
                //                                                                                //
                //********************************************************************************//

                // reset solver variables for each frame
                iDielecBndySum = 0;
                iCrgedGridSum  = 0;
                vector<delphi_real>().swap(vctfGridCrg);
                vector<SGrid<delphi_integer>>().swap(vctgiGridCrgPose);
                iCrgBdyGrid    = 0;
                vector<SDoubleGridValue>().swap(vctdgvCrgBndyGrid);
                vector<delphi_real>().swap(vctfPhiMap);
                vector<delphi_real>().swap(vctfPhiMap_Pre);

                unique_ptr<IAbstractModule> pSolver(new CDelphiSolver(pDataContainer, pTimer));
                pSolver->run();
                pSolver.reset();

#ifdef DEVELOPER
                cout << "\n\n---------- delphicpp finishes SOLVER class on ";
                pTester->showTime();
                cout << "\n\n";
#endif

                //********************************************************************************//
                //                                                                                //
                //          realize an object of CDelphiEnergy class to calculate energies        //
                //                                                                                //
                //********************************************************************************//
                unique_ptr<IAbstractModule> pEnergy(new CDelphiEnergy(pDataContainer, pTimer));
                pEnergy->run();
                pEnergy.reset();


#ifdef DEVELOPER
                cout << "\n\n---------- delphicpp finishes ENERGY class on ";
                pTester->showTime();
                cout << "\n\n";
#endif

                if (iGaussian == 1 && inhomo == 1 && logs) //second run for Gaussian
                {

                    //pDataContainer.reset();

                    //shared_ptr<IDataContainer> pDataContainer( new CDelphiData(argc,argv,pTimer) );
                    inhomo = 0;

                    unique_ptr < IAbstractModule > pSpace(new CDelphiSpace(pDataContainer, pTimer));
                    pSpace->run();
                    pSpace.reset();

#ifdef DEBUG_DELPHI_SPACE
                    cout << "---[DEBUG_DELPHI_SPACE] data container is written into file test_delphicpp_aftersurf_run2_frame" + to_string(iFrame) + ".dat---\n\n";
                    pDataContainer->showMap("test_delphicpp_aftersurf_run2_frame" + to_string(iFrame) + ".dat");
#endif

                    unique_ptr < IAbstractModule > pSolver(new CDelphiSolver(pDataContainer, pTimer));
                    pSolver->run();
                    pSolver.reset();

#ifdef DEBUG_DELPHI_SOLVER
                    cout << "---[DEBUG_DELPHI_SOLVER] data container is written into file test_delphicpp_afteritr_run2_frame" + to_string(iFrame) + ".dat---\n\n";

                    pDataContainer->showMap("test_delphicpp_afteritr_run2_frame" + to_string(iFrame) + ".dat");
#endif
                    unique_ptr < IAbstractModule > pEnergy(new CDelphiEnergy(pDataContainer, pTimer));
                    pEnergy->run();
                    pEnergy.reset();
                }

                // After energy calculations are complete, 
                // write the requested energy terms 
                // into the file
                ofEnergyFile.open(strEnergyFile.c_str(), ofstream::app);
                ofEnergyFile << setw(15) << left << iFrame;

                ofEnergyFile << fixed << setprecision(3);
                if (bGridEng)        ofEnergyFile << setw(15) << left << fEngGrid; 
                if (bSolvEng)        ofEnergyFile << setw(15) << left << fEngCorrect;
                if (bIonsEng)        ofEnergyFile << setw(15) << left << fEngIons;
                if (bCoulombEng)     ofEnergyFile << setw(15) << left << fEngCoul;
                if (bLJEng)          ofEnergyFile << setw(15) << left << fEngLJ;
                if (bNonPolEng)      ofEnergyFile << setw(15) << left << fEngNonPolCav_SA << setw(15) << left << fEngNonPolCav_Vol << setw(15) << left << fEngNonPolCav;
                ofEnergyFile << endl;
                ofEnergyFile.close();

                //********************************************************************************//
                //                                                                                //
                //               realize an object of CSite class to write site info              //
                //                                                                                //
                //********************************************************************************//

                unique_ptr<CSite> pSite(new CSite(pDataContainer, pTimer));

                if (pDataContainer->getKey_constRef<bool>("isite"))
                {
                    int iisitsf = 0;
                    if (pDataContainer->getKey_Ref<bool>("isitsf")) iisitsf = 1;
                    pSite->writeSite(iisitsf);
                }

                if (pDataContainer->getKey_constRef<bool>("phiwrt")) pSite->writePhi();
                if (pDataContainer->getKey_constRef < delphi_real > ("radipz") > 0.) pSite->wirtePAnalysis();

                pSite.reset();


#ifdef DEVELOPER
                cout << "\n\n---------- delphicpp finishes SITE class on ";
                pTester->showTime();
                cout << "\n\n";
#endif

                // point to the next frame
                itFrame++;

                // incr the frame index
                iFrame++;

            }// for each frame iterator (itFrame)

            pDataContainer.reset();

            pTimer->exit();
            pTimer.reset();

            delete pTester;

        } // ---------- end of try block

    } // ---------- end of try block
    catch (CException&)
    {
        cerr << "\n\n ......... PROGRAM ABORTS WITH AN EXCEPTION AND "
            << CWarning::iWarningNum << " WARNING(S) ........\n\n";

        return 1;
    }

    cout << "\n\n .........  PROGRAM EXITS SUCCESSFULLY : WITH TOTAL "
        << CWarning::iWarningNum << " WARNING(S) ........\n\n";

    int numWarnings = 0;
    if (iMaxWarn >= CWarning::iWarningNum) 
    {
        numWarnings = CWarning::iWarningNum;
    } 
    else 
    {
        numWarnings = iMaxWarn;
    }

    cout << " .....................  SHOWING " << numWarnings << " OUT OF " 
        <<  CWarning::iWarningNum << " WARNING(S) ..................." << endl;


    // ARGO: modified the method of printing warnings based on the number requested by the user
    // Specifically, the buffer is made to have "|" char which act as delimiters.
    // The buffer is split by that delimiter.

    string warningStr = cwarn.rdbuf()->str()  + "|";
    int warnCount = -1;
    std::size_t current, previous = 0;

    current = warningStr.find_first_of("|");
    while (current != std::string::npos ) 
    {
        if (warnCount == iMaxWarn) break; 
        cout << warningStr.substr(previous, current - previous) << endl;
        previous = current + 1;
        current = warningStr.find_first_of("|", previous);
        warnCount++;
    }

    cerr << "\033[0m";
    cout.unsetf(ios_base::floatfield); // return to cout default notation

    return 0;

}
