/*
 * mpi_setCrg.cpp
 *
 *  Created on: Aug 26, 2019
 *      Author: chuanli
 */

#include "solver.h"

//-----------------------------------------------------------------------//
void CDelphiSolver::mpi_setCrg()
{
    //++++++++++ INPUT:
    //const vector<SGridValue<delphi_real> >& chrgv2 = prggvCrg2Grid;

    //const vector<SGrid<delphi_integer> >& iepsmp = prgigEpsMap;
    SGrid<delphi_integer> iepsmp_val;

    vector<char>& idebmap = prgbDielecMap;
    const int& idirectalg = iDirectEpsMap;
    //const vector<delphi_real>& medeps = prgfMediaEps;
    //vector<delphi_integer>::const_iterator nqgrdtonqass = prgiCrg2GridMap.cbegin();
    //const vector<delphi_real>& atmeps = prgfAtomEps;

    //++++++++++ LOCAL:
    int kb1, kb2, kb3, i, j, k, m, ix, iy, iz;
    delphi_integer iw;
    delphi_real cg1, cg2, cg3, re, deb;
    SGrid < delphi_integer > igCoord;
    vector < delphi_real > phimap(iGrid * iGrid * iGrid, 0.0); // phimap(1:igrid,1:igrid,1:igrid) = 0.0

    delphi_integer inearest[6], iext, ibgp, cont, itmp, ico;
    delphi_real temp;
    vector<delphi_real> gchrgd, gchrgtmp;
    SDoubleGridValue fdbGridVal;
    vector < delphi_integer > gchrg2;
    vector < delphi_real > vecttemp(6, 0.0); // vecttemp
    vector < vector < delphi_real >> gaussianChargeDielecTemp;
    vector < delphi_real > gaussianChargeDensityTemp;

    //++++++++++ OUTPUT:
    delphi_integer& icount1a = iCrgedGridEven;
    delphi_integer& icount1b = iCrgedGridSum;
    delphi_integer& ibc = iCrgBndyGridNum;
    vector < delphi_integer > &iqpos = prgiCrgPose;
    vector<delphi_real>& qval = prgfCrgValA, &gval = prgfCrgValG, &gchrg = prgfGridCrg;
    vector < SGrid<delphi_integer> > &gchrgp = prgigGridCrgPose;
    vector < SDoubleGridValue > &cgbp = prgdgvCrgBndyGrid;

    iext = 0;
    ibgp = 0;

    size_t mpi_size = pdc->sizeofGlobal1D<SGridValue<delphi_real>>("chrgv2");
    SGridValue<delphi_real> mpi_chrgv2;

    delphi_integer mpi_nqgrdtonqass_indx, mpi_nqgrdtonqass;
    delphi_real mpi_atmeps;

    for (size_t itr = 0; itr < mpi_size; itr++)
    {
        mpi_chrgv2 = pdc->readGlobalVector1D< SGridValue<delphi_real> >("chrgv2", itr);

        kb1 = mpi_chrgv2.nGrid.nX;
        kb2 = mpi_chrgv2.nGrid.nY;
        kb3 = mpi_chrgv2.nGrid.nZ;

        for ( ix = 0; ix <= 1; ix++ )
        {
            i = kb1 + ix;
            cg1 = kb1 - mpi_chrgv2.nGrid.nX + 1 - ix;
            for ( iy = 0; iy <= 1; iy++ )
            {
                j = kb2 + iy;
                cg2 = kb2 - mpi_chrgv2.nGrid.nY + 1 - iy;
                for ( iz = 0; iz <= 1; iz++ )
                {
                    k = kb3 + iz;
                    cg3 = kb3 - mpi_chrgv2.nGrid.nZ + 1 - iz;
                    re = abs(cg1 * cg2 * cg3) * (mpi_chrgv2.nValue);
                    phimap[(k - 1) * iGrid * iGrid + (j - 1) * iGrid + (i - 1)] += re;
                }
            }
        }
    }

    for ( k = 2; k < iGrid; k++ )
    {
        for ( j = 2; j < iGrid; j++ )
        {
            for ( i = 2; i < iGrid; i++ )
            {
                m = (k - 1) * iGrid * iGrid + (j - 1) * iGrid + (i - 1);
                if ( abs(phimap[m]) > fZero ) // phimap(i,j,k).ne.0.
                {
                    igCoord.nX = i;
                    igCoord.nY = j;
                    igCoord.nZ = k;
                    gchrgp.push_back(igCoord);
                    gchrg.push_back(phimap[m]);
                    phimap[m] = 0.0;

                    //---------- determine how many charged grid points are even (odd in Fortran)
                    if ( 0 != optSum < delphi_integer > (igCoord) % 2 )
                        icount1a += 1;
                }
            }
        }
    }

    mpi_nqgrdtonqass_indx = 0;

    for (size_t itr = 0; itr < mpi_size; itr++)
    {
        mpi_chrgv2 = pdc->readGlobalVector1D< SGridValue<delphi_real> >("chrgv2", itr);

        kb1 = mpi_chrgv2.nGrid.nX;
        kb2 = mpi_chrgv2.nGrid.nY;
        kb3 = mpi_chrgv2.nGrid.nZ;

        for ( ix = 0; ix <= 1; ix++ )
        {
            i = kb1 + ix;
            cg1 = kb1 -  mpi_chrgv2.nGrid.nX + 1 - ix;
            for ( iy = 0; iy <= 1; iy++ )
            {
                j = kb2 + iy;
                cg2 = kb2 -  mpi_chrgv2.nGrid.nY + 1 - iy;
                for ( iz = 0; iz <= 1; iz++ )
                {
                    k = kb3 + iz;
                    cg3 = kb3 -  mpi_chrgv2.nGrid.nZ + 1 - iz;
                    re = abs(cg1 * cg2 * cg3) * ( mpi_chrgv2.nValue);
                    deb = 0.0;
                    if ( idebmap[(k - 1) * iGrid * iGrid + (j - 1) * iGrid + (i - 1)] ) deb = 1.0;

                    mpi_nqgrdtonqass = pdc->readGlobalVector1D< delphi_integer >("nqgrdtonqass", mpi_nqgrdtonqass_indx);
                    mpi_atmeps       = pdc->readGlobalVector1D< delphi_real >("atmeps", mpi_nqgrdtonqass - 1);
                    phimap[(k - 1) * iGrid * iGrid + (j - 1) * iGrid + (i - 1)] += re / (6.0 * mpi_atmeps + fDebFct * deb); // atmeps[*nqgrdtonqass-1]
                }
            }
        }

        mpi_nqgrdtonqass_indx++; // ic1=nqgrdtonqass(ig)
    }

    icount1b = gchrg.size();

    for ( vector<SGrid<delphi_integer> >::iterator it = gchrgp.begin(); it != gchrgp.end(); ++it )
    {
        ix = it->nX;
        iy = it->nY;
        iz = it->nZ;
        iw = (iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1);
        gchrgtmp.push_back(phimap[iw]);
    }

    //vector<delphi_real>().swap(phimap); // remove the vector phimap. No need in below

    //---------------------------- Begin of wrtcrg.f -----------------------------//
    if ( bGridCrgOut )
    {
        string strCrgFile = "crg.dat"; // charge file name
        ofstream ofCrgFileStream;
        ofCrgFileStream.open(strCrgFile.c_str());

        ofCrgFileStream << fixed << setprecision(3);

        ofCrgFileStream << "DELPHI OUTPUT FILE: GRID CHARGE" << endl;
        ofCrgFileStream << "FORMAT NUMBER= 1" << endl;
        ofCrgFileStream << "NUMBER OF CHARGES=" << setw(8) << right << icount1b << endl;

        delphi_real mpi_medeps;

        for ( int i = 0; i <= iMediaNum; i++ )
        {
            mpi_medeps = pdc->readGlobalVector1D< delphi_real >("medeps", i);
            ofCrgFileStream << "DIELECTRIC IN MEDIUM NUMBER " << setw(3) << right << i << " : " << setw(8)
                            << right << mpi_medeps * fEPKT << endl;

        }
        ofCrgFileStream << "GRID SCALE=" << fScale << endl;

        for ( int i = 0; i < icount1b; i++ )
        {
            ofCrgFileStream << scientific << setprecision(17) << setw(25) << right << gchrg[i];
            ofCrgFileStream << fixed << setw(8) << right << gchrgp[i].nX << setw(8)
                            << right << gchrgp[i].nY << setw(8) << right << gchrgp[i].nZ << endl;
        }

        ofCrgFileStream.close();
    }
    //---------------------------- End of wrtcrg.f -----------------------------//

    //---------- set up odd/even pointer array, to be used in making qval and iqpos
    i = 1;
    j = icount1a + 1;
    for ( vector<SGrid<delphi_integer> >::iterator it = gchrgp.begin(); it != gchrgp.end(); ++it )
    {
        if ( 0 != optSum < delphi_integer > (*it) % 2 )
        {
            gchrg2.push_back(i);
            i++;
        }
        else
        {
            gchrg2.push_back(j);
            j++;
        }
    }

    //---------- determine denominator at all charged grid points
    if ( 0 != ibc ) ibc = 0;
    ico = 0;
    i = 0;

    for ( vector<SGrid<delphi_integer> >::iterator it = gchrgp.begin(); it != gchrgp.end(); ++it )
    {
        ix = it->nX;
        iy = it->nY;
        iz = it->nZ;

        if ( iGaussian == 0 )
        {
            iext = 0;
            ibgp = 0;

            iw          = (iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1);
            iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);
            inearest[0] = iepsmp_val.nX / iEpsDim;
            inearest[1] = iepsmp_val.nY / iEpsDim;
            inearest[2] = iepsmp_val.nZ / iEpsDim;

            iw          = (iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1 - 1);
            iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);
            inearest[3] = iepsmp_val.nX / iEpsDim;

            iw          = (iz - 1) * iGrid * iGrid + (iy - 1 - 1) * iGrid + (ix - 1);
            iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);
            inearest[4] = iepsmp_val.nY / iEpsDim;

            iw          = (iz - 1 - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1);
            iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);
            inearest[5] = iepsmp_val.nZ / iEpsDim;

            if ( 0 == inearest[0] )
                iext = 1;

            if ( inearest[0] != inearest[5] )
                ibgp = 1;

            for ( cont = 1; cont < 6; cont++ )
            {
                if ( 0 == inearest[cont] )
                    iext = 1;
                if ( inearest[cont - 1] != inearest[cont] )
                    ibgp = 1;
            }
        } //iGaussian == 0 |--> 2-eps model

        deb = 0.0;

        if ( idebmap[(iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1)] )
            deb = 1.0;

        vector < delphi_real > dbrow_original;
        delphi_real gridDensity;
        if (iGaussian==1)
        {
            //gridDensity = gaussianDensityMap[(iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1)];
            iw = (iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1);
            gridDensity = pdc->readGlobalVector1D< delphi_real >("gdensity", iw);
        }

        if ( 0 == idirectalg )
        {
            throw CDirectEpsilonMap(idirectalg);
        }
        else
        {
            if ( iGaussian == 0 )
            {
                temp = 0.0;
                delphi_real mpi_medeps;

                iw          = (iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1);
                iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);

                itmp        = iepsmp_val.nX / iEpsDim;
                mpi_medeps  = pdc->readGlobalVector1D< delphi_real >("medeps", itmp);
                temp       += mpi_medeps;
                itmp        = iepsmp_val.nY / iEpsDim;
                temp       += mpi_medeps;
                itmp        = iepsmp_val.nZ / iEpsDim;
                temp       += mpi_medeps;

                iw          = (iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1 - 1);
                iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);
                itmp        = iepsmp_val.nX / iEpsDim;
                mpi_medeps  = pdc->readGlobalVector1D< delphi_real >("medeps", itmp);
                temp       += mpi_medeps;

                iw          = (iz - 1) * iGrid * iGrid + (iy - 1 - 1) * iGrid + (ix - 1);
                iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);
                itmp        = iepsmp_val.nY / iEpsDim;
                mpi_medeps  = pdc->readGlobalVector1D< delphi_real >("medeps", itmp);
                temp       += mpi_medeps;

                iw          = (iz - 1 - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix - 1);
                iepsmp_val  = pdc->readGlobalVector1D< SGrid<delphi_integer> >("iepsmp", iw);
                itmp        = iepsmp_val.nZ / iEpsDim;
                mpi_medeps  = pdc->readGlobalVector1D< delphi_real >("medeps", itmp);
                temp       += mpi_medeps;

                temp += fDebFct * deb;

                gchrgd.push_back(temp);
            }
            else if ( iGaussian == 1 )
            {
                if ( uniformdiel )
                {
                    temp = fEpsIn * 6;
                    gchrgd.push_back(temp + fDebFct * deb);
                }
                else
                {
                    SGrid<delphi_real> mpi_gepsmp2;

                    temp        = 0;

                    iw          = (ix - 1) * iGrid * iGrid + (iy - 1) * iGrid + (iz - 1);
                    mpi_gepsmp2 = pdc->readGlobalVector1D< SGrid<delphi_real> >("gepsmp2", iw);

                    temp       += mpi_gepsmp2.nX;
                    vecttemp[0] = mpi_gepsmp2.nX;

                    temp       += mpi_gepsmp2.nY;
                    vecttemp[1] = mpi_gepsmp2.nY;

                    temp       += mpi_gepsmp2.nZ;
                    vecttemp[2] = mpi_gepsmp2.nZ;

                    iw          = (ix - 1 - 1) * iGrid * iGrid + (iy - 1) * iGrid + (iz - 1);
                    mpi_gepsmp2 = pdc->readGlobalVector1D< SGrid<delphi_real> >("gepsmp2", iw);

                    temp       += mpi_gepsmp2.nX;
                    vecttemp[3] = mpi_gepsmp2.nX;

                    iw          = (ix - 1) * iGrid * iGrid + (iy - 1 - 1) * iGrid + (iz - 1);
                    mpi_gepsmp2 = pdc->readGlobalVector1D< SGrid<delphi_real> >("gepsmp2", iw);

                    temp       += mpi_gepsmp2.nY;
                    vecttemp[4] = mpi_gepsmp2.nY;

                    iw          = (ix - 1) * iGrid * iGrid + (iy - 1) * iGrid + (iz - 1 - 1);
                    mpi_gepsmp2 = pdc->readGlobalVector1D< SGrid<delphi_real> >("gepsmp2", iw);

                    temp       += mpi_gepsmp2.nZ;
                    vecttemp[5] = mpi_gepsmp2.nZ;

                    temp = temp / fEPKT; //temp=temp/epkt
                    gchrgd.push_back(temp + fDebFct * deb); //gchrgd(i)=temp+ debfct*deb

                    vector < delphi_real > dbrow_original;
                    dbrow_original.push_back(vecttemp[3]);
                    dbrow_original.push_back(vecttemp[0]);
                    dbrow_original.push_back(vecttemp[4]);
                    dbrow_original.push_back(vecttemp[1]);
                    dbrow_original.push_back(vecttemp[5]);
                    dbrow_original.push_back(vecttemp[2]);

                    gaussianChargeDielecTemp.push_back(dbrow_original);
                    gaussianChargeDensityTemp.push_back(gridDensity);
                }
            }
        }

        if ( 1 == ibgp || 1 == iext )
        {
            ibc += 1;

            if ( 0 == ibgp ) ico += 1;

            fdbGridVal.fgCoord.nX = (delphi_real) ix;
            fdbGridVal.fgCoord.nY = (delphi_real) iy;
            fdbGridVal.fgCoord.nZ = (delphi_real) iz;
            fdbGridVal.fVal1 = gchrgtmp[i] * f4Pi * fScale;
            fdbGridVal.fVal2 = gchrg2[i];

            cgbp.push_back(fdbGridVal);
        }

        i++;
    } //----- end of for (vector< SGrid<delphi_integer> >::iterator it = gchrgp.begin(); it != gchrgp.end(); ++it)

    #ifdef VERBOSE
    cout << " Info> Number of charged boundary grid points = " << ibc << endl;
    if (0 != ibc) CCrgedPtsInSolution waring(ico);
    #endif

    //---------- make qval, fpoh term so potentials will be in kt/e
    qval.assign(icount1b, 0.0);
    gval.assign(icount1b, 0.0);
    vecttemp.assign(6, 0.0);
    gaussianChargeDielec.assign(icount1b, vecttemp);
    gaussianChargeDensity.assign(icount1b, 0.0);

    for ( i = 0; i < icount1b; i++ )
    {
        j = gchrg2[i];
        qval[j - 1] = gchrg[i] * (f4Pi * fScale / gchrgd[i]);
        gval[j - 1] = gchrg[i];
        if ( iGaussian != 0 )
        {
            gaussianChargeDielec[j - 1] = gaussianChargeDielecTemp[i];
            gaussianChargeDensity[j - 1] = gaussianChargeDensityTemp[i];
        }
    }

    vector<delphi_real>().swap(gchrgd); // remove gchrgd, need in below

    iqpos.assign(icount1b, 0);
    for ( i = 0; i < icount1b; i++ )
    {
        j = gchrg2[i];
        ix = gchrgp[i].nX;
        iy = gchrgp[i].nY;
        iz = gchrgp[i].nZ;
        iw = (iz - 1) * iGrid * iGrid + (iy - 1) * iGrid + (ix + 1);
        iqpos[j - 1] = iw / 2;
    }
}

