/*
 * delphi_datamarshal_updateTrajParameters.cpp
 *
 *  Created on: Feb 01, 2019
 *      Author: Arghya 
 */

#include "delphi_datamarshal.h"


//-----------------------------------------------------------------------//
void CDelphiDataMarshal::updateTrajParameters()
{
    // ++++++ update parameters after reading the parameter file ++++++//
    //*****************************************************************//
    //                                                                 //
    //        perform updates after reading the parameter file         //
    //                                                                 //
    //*************** For trajectory analyses/MMPBSA ******************//
#ifdef ENABLE_TRAJECTORY
    if (-1 == vctfProbeRadius[1])
        vctfProbeRadius[1] = vctfProbeRadius[0];

    /*
     * parameters set at the end of subroutine rdprm
     */
    if (0 > fExDielec || 0 > fInDielec)
    {
        CMinusDielec warning(fExDielec, fInDielec);

        fExDielec = abs(fExDielec);
        fInDielec = abs(fInDielec);
    }

    delphi_real fZ1Plus  = vctiValence1[0]; // valence 1
    delphi_real fZ1Minus = vctiValence1[1];
    delphi_real fZ2Plus  = vctiValence2[0]; // valence 2
    delphi_real fZ2Minus = vctiValence2[1];

    /*
     * concentration of positive ion
     */
    delphi_real fZ1PlusConcentrate = vctfSalt[0]*fZ1Minus;
    delphi_real fZ2PlusConcentrate = vctfSalt[1]*fZ2Minus;

    fIonStrength = (fZ1PlusConcentrate*fZ1Plus*(fZ1Plus+fZ1Minus) + fZ2PlusConcentrate*fZ2Plus*(fZ2Plus+fZ2Minus))/2.0;

    /*
     * coefficients in Taylor series of the charge concentration apart from n! (order >=1)
     * (chuan 2012Apr24) Correct coefficients in Taylor series. NOT in compact form just for clean math formula
     */
    fTaylorCoeff1 = -2.0*fIonStrength;

    fTaylorCoeff2 =  ( fZ1PlusConcentrate*fZ1Plus *pow(fZ1Plus, 2) -
                       fZ1PlusConcentrate*fZ1Minus*pow(fZ1Minus,2) +
                       fZ2PlusConcentrate*fZ2Plus *pow(fZ2Plus, 2) -
                       fZ2PlusConcentrate*fZ2Minus*pow(fZ2Minus,2) )/2.0;

    fTaylorCoeff3 = -( fZ1PlusConcentrate*fZ1Plus *pow(fZ1Plus, 3) +
                       fZ1PlusConcentrate*fZ1Minus*pow(fZ1Minus,3) +
                       fZ2PlusConcentrate*fZ2Plus *pow(fZ2Plus, 3) +
                       fZ2PlusConcentrate*fZ2Minus*pow(fZ2Minus,3) )/6.0;

    fTaylorCoeff4 =  ( fZ1PlusConcentrate*fZ1Plus *pow(fZ1Plus, 4) -
                       fZ1PlusConcentrate*fZ1Minus*pow(fZ1Minus,4) +
                       fZ2PlusConcentrate*fZ2Plus *pow(fZ2Plus, 4) -
                       fZ2PlusConcentrate*fZ2Minus*pow(fZ2Minus,4) )/24.0;

    fTaylorCoeff5 = -( fZ1PlusConcentrate*fZ1Plus *pow(fZ1Plus, 5) +
                       fZ1PlusConcentrate*fZ1Minus*pow(fZ1Minus,5) +
                       fZ2PlusConcentrate*fZ2Plus *pow(fZ2Plus, 5) +
                       fZ2PlusConcentrate*fZ2Minus*pow(fZ2Minus,5) )/120.0;

    /*
     * convert ionic strength to debye length
     */
     if (fZero < fIonStrength)
     {
         delphi_real fDebyeFactor = 0.01990076478*sqrt(fTemper*fExDielec);
         fDebyeLength = fDebyeFactor/sqrt(fIonStrength);

         if (0 < iNonIterateNum) bNonlinearEng = true;
     }
     else
     {
         bIonsEng = false;
         fDebyeLength = 1.0e6;
     }

	/*
	 * epkt assignment as a function of temperature
	 */
	fEPKT = dEPK/fTemper;

	/*
	 * set epsin and epsout (= epkt adjusted dielectrics such that all distances are in angstroms, charges in e)
	 */
	fEpsIn  = fInDielec/fEPKT;
	fEpsOut = fExDielec/fEPKT;


	string strASCI = "1234567890 .-+#,$asdfghjklzxcvbnmqwertyuiopASDFGHJKLZXCVBNMQWERTYUIOP)(}{][/";

	ifstream inTrajHandle, inTopolHandle;

	inTrajHandle.open(strTrajFile.c_str());
	if ( !inTrajHandle.is_open() ) throw CUnknownFile(strTrajFile);
	inTrajHandle.close();

	inTopolHandle.open(strTopolFile.c_str());
	if ( !inTopolHandle.is_open() ) throw CUnknownFile(strTopolFile);
	inTopolHandle.close();


	//*****************************************************************//
	//                                                                 //
	//        read size, charge, topology and trajectory files         //
	//                                                                 //
	//*****************************************************************//
	unique_ptr<CIO> pIO(new CIO(fInDielec,fEPKT)); // smart unique_ptr

	pIO->setDelphiTraj(strTrajFile, iTrajFormatIn, strTopolFile, iTopolFormatIn, strSizeFile, strCrgFile, strVdwFile, bSolvePB, bLJEng, bNonPolEng, iFrameFirst, iFrameLast, iFrameStride);
	
	iMediaNum        = pIO->iMediaNum;         // nmedia
	iObjectNum       = pIO->iObjectNum;        // nobject
	iAtomNum         = pIO->iAtomNum;          // natom
	iResidueNum      = pIO->iResidueNum;       // resnummax
	bOnlyMolecule    = pIO->bOnlyMolecule;     // ionlymol (true)

	vctapAtomPdb.clear();                      // keep it clean here
	vctfMediaEps     = pIO->vctfMediaEps;      // medeps(0:nmediamax)
	vctstrObject     = pIO->vctstrObject;      // dataobject(nobjectmax,2)
	vctiAtomMediaNum = pIO->vctiAtomMediaNum;  // iatmmed(Natom+Nobjectmax)
	
	iNumFrames       = pIO->vctFrames.size();  // nframes
	vctSimFrames     = pIO->vctFrames;
	vctSimAtoms      = pIO->vctSimAtoms;

	fEpsIn = pIO->vctfMediaEps[1];
	
	if (1 < iMediaNum) iDirectEpsMap = 1;
	
	cout << " number of objects = " << iObjectNum << endl;
	// ---------------------- extremes, gsize, scales, perfils for each frame --------------- //
	
    iMoleculeNum = 1;    // numbmol
    fMaxRadius   = 0.01; // rdmx
	

	if (-1 > iAtomNum) throw CNoAtomsInMolecule(0);

	for ( CSimAtom& tmpAtom: vctSimAtoms )
		fMaxRadius = max(fMaxRadius, tmpAtom.getRadius()); 

	for ( CFrame& tmpFrame: vctSimFrames )
	{
		// assign box centers
		if (bIsAcent)
			tmpFrame.gfFrameBoxCenter = gfAcent;
		else
			tmpFrame.gfFrameBoxCenter = tmpFrame.gfFrameMidPoint - gfOffCenter/fScale;

		// compute side with max dimension
		SGrid<delphi_real> gVec1 = 2.0*optABS<delphi_real>(tmpFrame.gfFrameMaxCoordinate - tmpFrame.gfFrameBoxCenter);
		delphi_real fMaxVal1 = optMax<delphi_real>(gVec1);

		SGrid<delphi_real> gVec2 = 2.0*optABS<delphi_real>(tmpFrame.gfFrameMinCoordinate - tmpFrame.gfFrameBoxCenter);
		delphi_real fMaxVal2 = optMax<delphi_real>(gVec2);

		fMaxDimension = (fMaxVal1 > fMaxVal2) ? fMaxVal1 : fMaxVal2;

		// iFrameGrid, fFrameScale and fFramePercentageFill
		if (0 == iGrid)
		{
			tmpFrame.fFrameScale = fScale;
			tmpFrame.fFramePercentageFill = fPercentageFill;

			if (fZero > abs(tmpFrame.fFrameScale - 10000.0) )          tmpFrame.fFrameScale = 2.0;
			if (fZero > abs(tmpFrame.fFramePercentageFill - 10000.0) ) tmpFrame.fFramePercentageFill = 80.0;

			tmpFrame.iFrameGrid = tmpFrame.fFrameScale * 100.0/tmpFrame.fFramePercentageFill * fMaxDimension;
		}
		else if (fZero > abs(fScale-10000.0) )
		{
			tmpFrame.iFrameGrid = iGrid;
			tmpFrame.fFramePercentageFill = fPercentageFill;
			if (fZero > abs(tmpFrame.fFramePercentageFill - 10000.0) )
			{
				tmpFrame.fFrameScale          = 2.0;
				tmpFrame.fFramePercentageFill = 100.0 * fMaxDimension * tmpFrame.fFrameScale/( tmpFrame.iFrameGrid - 1 );
			}
			else
			{
				tmpFrame.fFrameScale = (tmpFrame.iFrameGrid - 1) * tmpFrame.fFramePercentageFill/(100.0 * fMaxDimension);
			}
		}
		else
		{
			tmpFrame.iFrameGrid = iGrid;
			tmpFrame.fFrameScale = fScale;
			tmpFrame.fFramePercentageFill = 100.0*fMaxDimension * tmpFrame.fFrameScale/(tmpFrame.iFrameGrid - 1);
		}

		if (0 == tmpFrame.iFrameGrid%2) tmpFrame.iFrameGrid += 1;

		// check if within bounds
		SGrid<delphi_real> gLeftBndy  = tmpFrame.gfFrameBoxCenter - (1.0 / tmpFrame.fFrameScale)*(tmpFrame.iFrameGrid + 1)*0.5;
		SGrid<delphi_real> gRightBndy = tmpFrame.gfFrameBoxCenter + (1.0 / tmpFrame.fFrameScale)*(tmpFrame.iFrameGrid + 1)*0.5;
		if (optORLT(tmpFrame.gfFrameMinCoordinate, gLeftBndy) || optORGT(tmpFrame.gfFrameMaxCoordinate, gRightBndy)) 
			CSystemOutsideBox warning;


		/*
		* convert atom coordinates from angstroms to grid units
		*/
		for ( SGrid<delphi_real>& pos : tmpFrame.coordinates )
		{
			tmpFrame.aCoordinates.push_back( pos );
			tmpFrame.gCoordinates.push_back( (pos - tmpFrame.gfFrameBoxCenter) * tmpFrame.fFrameScale + (delphi_real)((tmpFrame.iFrameGrid + 1)/2) );
		}

		// assign debye length
		tmpFrame.fDebyeNum = (100.0/tmpFrame.fFramePercentageFill-1.0) * fMaxDimension/fDebyeLength;
	} //---------------- for tmpFrame
	
	// uniform dielectric is set to false (though relevant for objects only, i have tried to be consistent)
	vctfMediaEps[0] = fEpsOut;
	bUniformDielec = true;
	for ( delphi_integer i = 0; i < iMediaNum; i++ )
	{
		if (vctfMediaEps[i] != vctfMediaEps[i+1]) bUniformDielec = false;
	}

	/*
	 * new updates in c++ for more rigid check of input parameter values
	 */
	if (false == bAutoConverge && 0 == iLinIterateNum && 0 == iNonIterateNum)
		throw CBadAutoConvergence(bAutoConverge);

	showTrajParameters();

#ifdef VERBOSE
	int iFrame = 0;
	for ( CFrame& tmpFrame: vctSimFrames ) 
	{
		cout << "------------------- Frame " << iFrame << " ---------------------" << endl;
		// tmpFrame.dumpFrame();
		iFrame++;
	}
#endif

	pIO.reset();
#endif
}

