/* Copyright (c) 2015  Gerald Knizia
 * 
 * This file is part of the IboView program (see: http://www.iboview.org)
 * 
 * IboView is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 * 
 * IboView is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
 * 
 * Please see IboView documentation in README.txt for:
 * -- A list of included external software and their licenses. The included
 *    external software's copyright is not touched by this agreement.
 * -- Notes on re-distribution and contributions to/further development of
 *    the IboView software
 */

#include <math.h>
#include <stdexcept>
#include "Ir.h"
#include "CtDft.h"
#include "CxAlgebra.h"
#include "CxFortranInt.h"

namespace ct {


#define GRID_OBTAIN_GRADWT_THRD FORT_Extern(grid_obtain_gradwt_thrd,GRID_OBTAIN_GRADWT_THRD)
void GRID_OBTAIN_GRADWT_THRD(FINTARG icatom, FINTARG npt, double const (*pXyz)[3], double const *pGridWt, FINTARG jwt, double *pWtGrd, double *buf)
{
   assert_rt(0); // not implemented here.
}

// calculate buffer size for grid_obtain_gradwt_thrd; defined in dftgrid.f
#define GRID_OBTAIN_GRADWT_BUFSIZE FORT_Extern(grid_obtain_gradwt_bufsize,GRID_OBTAIN_GRADWT_BUFSIZE)
void GRID_OBTAIN_GRADWT_BUFSIZE(FORTINT *nbuf, FINTARG icatom, FINTARG npt, FINTARG jwt)
{
   assert_rt(0); // not implemented here.
}



// okay, now comes the ugly part:
// for what we do here, at the moment we do not define proper interfaces
// to grids or density functionals. There are simply a few *GLOBAL*
// objects to that degree, which are set externally. For the moment that
// is enough. Got other problems.

FXcFunctionalPtr
   g_pXcFunctional;


void DFTFUN_CXX(FINTARG fderiv, FINTARG open, FINTARG igrad, FINTARG nGridPt_, double const *wt,
         // density inputs.
         double const *rhoc, double const *rhoo, double const *sigmacc, double const *sigmaco, double const *sigmaoo, double const *tauc, double const *tauo, double const *upsilonc, double const *upsilono,
         // density functional outputs (summed over functionals)
         double *zk, double *vrhoc, double *vrhoo, double *vsigmacc, double *vsigmaco, double *vsigmaoo, double *vtauc, double *vtauo, double *vupsilonc, double *vupsilono,
         // energy outputs (split according to functions; one output per ndftu)
         double *dfu_energies,
         // temporary buffer (currently: len npt).
         double *tmp,
         // flags indicating whether tau/upsilon arrays were passed (to control our rather hacky way of scaling dftfun inputs..)
         FINTARG iHaveTau, FINTARG iHaveUpsilon)
{
   size_t
      nGridPt = nGridPt_;
   double
      *pZeros = 0;
   if ( open == 0 ) {
      // pure closed shell case -- code is written in such a way that these
      // quantities are never allocated in this case.
      pZeros = tmp;
      memset(pZeros, 0, sizeof(pZeros[0]) * nGridPt);
      rhoo = pZeros;
      sigmaco = pZeros;
      sigmaoo = pZeros;
      tauo = pZeros;
   }

   FXcFunctional
      *pXcFunctional = g_pXcFunctional.get();

   // copy input into linearized format expected by FXcFunctional.
   double
      *pIn,
      *pOut;
   if ( iHaveUpsilon )
      throw std::runtime_error("Upsilon-style Meta-GGAs not supported."  " DFTFUN_CXX_PROXY");
   uint
      nComp = 2;
   bool
      NeedSigma = pXcFunctional->NeedSigma(),
      NeedTau = pXcFunctional->NeedTau();
   if ( NeedSigma )
      nComp = 5;
   if ( NeedTau )
      nComp = 7;
   assert((bool)iHaveTau == NeedTau);
   assert(NeedSigma || !NeedTau); // tau implies sigma.

   size_t
      nInSt = nComp,
      nOutSt = nInSt + 1; // one for the energies, rest for derivatives wrt input
//    Mem.Alloc(pIn, nInSt*nGridPt);
//    Mem.Alloc(pOut, nOutSt*nGridPt);
   if ( fderiv == 0 )
      nOutSt = 1; // evaluate energy only, not potentials.
   pIn = tmp + nGridPt; // first one is for zeros in closed-shell case.
   pOut = pIn + nInSt * nGridPt;
   for ( size_t i = 0; i < nGridPt; ++ i ) {
//       pIn[i*nInSt + 0] = rhoc[i] + 1e-30;
      pIn[i*nInSt + 0] = rhoc[i] + 1e-15; // PBEC might blow up otherwise...
      pIn[i*nInSt + 1] = rhoo[i];
      if ( NeedSigma ) {
         pIn[i*nInSt + 2] = sigmacc[i];
         pIn[i*nInSt + 3] = sigmaco[i];
         pIn[i*nInSt + 4] = sigmaoo[i];
         if ( NeedTau ) {
            pIn[i*nInSt + 5] = tauc[i];
            pIn[i*nInSt + 6] = tauo[i];
         };
      }
   }

   memset(pOut, 0, sizeof(pOut[0])*nGridPt*nOutSt);
   pXcFunctional->Eval(pOut, nOutSt, pIn, nInSt, nGridPt, fderiv);

   // copy back outputs into target format.
   for ( size_t i = 0; i < nGridPt; ++ i )
      zk[i] = pOut[i*nOutSt];
   if ( fderiv != 0 ) {
      for ( size_t i = 0; i < nGridPt; ++ i ) {
         vrhoc[i] = pOut[i*nOutSt + 1];
         vrhoo[i] = pOut[i*nOutSt + 2];
         if ( NeedSigma ) {
            vsigmacc[i] = pOut[i*nOutSt + 3];
            vsigmaco[i] = pOut[i*nOutSt + 4];
            vsigmaoo[i] = pOut[i*nOutSt + 5];
            if ( NeedTau ) {
               vtauc[i] = pOut[i*nOutSt + 6];
               vtauo[i] = pOut[i*nOutSt + 7];
            }
         }
      }
   }
//    for ( size_t i = 0; i < nGridPt; ++ i ) {
//       double rho = pIn[i*nInSt];
//       double h = rho * 0.0001;
//       double dout[3] = {0};
//       double di[2] = {0};
//       di[0] = rho + h;
//       pXcFunctional->Eval(&dout[1], 1, &di[0], 2, 1, 0);
//       di[0] = rho - h;
//       pXcFunctional->Eval(&dout[0], 1, &di[0], 2, 1, 0);
//       vrhoc[i] = (dout[1] - dout[0])/(2*h);
//       vrhoo[i] = 0;
//       if ( isnan(vrhoc[i]) )
//          __asm__("int $0x03");
//    }

   dfu_energies[0] += Dot(wt, zk, nGridPt);
//    Mem.Free(pIn);
}




} // namespace ct
