/* 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
 */

#ifndef CT_RHF_FOCKBUILD_H
#define CT_RHF_FOCKBUILD_H

#include <string>
#include "CtBasisSet.h"
#include "CtAtomSet.h"
#include "CtMatrix.h"
// #include "CtDma.h"
#include "CtDftGrid.h" // for grid params.
#include "CtDftFunc.h"
#include "CtFockBuild.h"
#include "CtRhfOptions.h"
#include "CtTiming.h"
#include "CtIo.h"

namespace ct {

enum FFockBuilderFlags {
   FOCKBUILD_ClosedOnly //
};


struct FFockComponentBuilder : public FIntrusivePtrDest
{
   explicit FFockComponentBuilder(FLog &Log_, FTimerSet *pTimers_=0) : Energy(0.), m_Log(Log_), m_pTimers(pTimers_) {};

   virtual void Init(FWfDecl const &WfDecl_, FBasisSet *pOrbBasis_, FAtomSet const &Atoms_, FHfOptions const &Options_, FMemoryStack &Mem) = 0;
   // Accumulate fock matrix contribution to FockC/FockO.
   // Notes:
   //   - OrbC and OrbO are supposed to have their respective occupation numbers absorbed,
   //     such that the closed/open densities are obtained by DenC := OrbC x OrbC.T and DenO := OrbO x OrbO.T
   //   - Both orbitals and Fock matrices are expressed in terms of pBasis. pBasis may differ from pOrbBasis given in Init().
   virtual void AccFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FMemoryStack &Mem);
   virtual ~FFockComponentBuilder() = 0;

   virtual void AccGradient(FMatrixView Gradient, FMatrixView COccC, FMatrixView COccO, FMemoryStack &Mem);

   double Energy;
   virtual void PrintEnergyContribs();

   // returns true if this fock builder did anything here. Default implementation returns false (command ignored---no grids)
   virtual bool SwitchToRefineGrid(FDftGridParams const &NewGridParams);
protected:
   FLog
      &m_Log;
   FTimerSet
      *m_pTimers;
   FWfDecl
      m_WfDecl;
private:
   FFockComponentBuilder(FFockComponentBuilder const &); // not implemented
   void operator = (FFockComponentBuilder const &); // not implemented
};
typedef boost::intrusive_ptr<FFockComponentBuilder>
   FFockComponentBuilderPtr;


typedef std::vector<FFockComponentBuilderPtr>
   FFockComponentBuilderList;
void BuildFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FFockComponentBuilderList &BuilderList, FMemoryStack &Mem);



struct FDfJkOptions
{
   FDfJkOptions() : fExchFactor(1.0) {}
   double fExchFactor;
};


// build both DF Coulomb and Exchange
struct FFockComponentBuilderDfJk : public FFockComponentBuilder
{
   FFockComponentBuilderDfJk(FDfJkOptions const &Options, FLog &Log_, FTimerSet *pTimers_);
   void Init(FWfDecl const &WfDecl_, FBasisSet *pOrbBasis_, FAtomSet const &Atoms_, FHfOptions const &Options_, FMemoryStack &Mem);
   void AccFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FMemoryStack &Mem);
   void AccGradient(FMatrixView Gradient, FMatrixView COccC, FMatrixView COccO, FMemoryStack &Mem);
   void PrintEnergyContribs();
   ~FFockComponentBuilderDfJk();
public:
   FDfJkOptions
      Options;
   FAtomSet const
      *pAtoms;
   FBasisSetPtr
      pFitBasis;
   FBasisSet
      *pOrbBasis;
   FMatrixView
      Jcd;
   size_t
      nFit;
   double
      EnergyCoulomb, EnergyExch, EnergyXc;
};


// build DF Coulomb only
struct FFockComponentBuilderDfCoul : public FFockComponentBuilder
{
   FFockComponentBuilderDfCoul(FDfJkOptions const &Options, FLog &Log_);
   void Init(FWfDecl const &WfDecl_, FBasisSet *pOrbBasis_, FAtomSet const &Atoms_, FHfOptions const &Options_, FMemoryStack &Mem);
   void AccFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FMemoryStack &Mem);
   ~FFockComponentBuilderDfCoul();
public:
   FDfJkOptions
      Options;
   FAtomSet const
      *pAtoms;
   FBasisSetPtr
      pFitBasis;
   FBasisSet
      *pOrbBasis;
   FMatrixView
      Jcd;
   size_t
      nFit;
};


// build df-coulomb and auxiliary expanded xc. Cache all integrals in memory.
// struct FFockComponentBuilderDfCoulXcCachedImpl;
struct FFockComponentBuilderDfCoulXcCached : public FFockComponentBuilder
{
   FFockComponentBuilderDfCoulXcCached(FDfJkOptions const &Options, FDftGridParams const &GridParams_, std::string const &XcFunctionalName, FLog &Log_, FTimerSet *pTimers_);
   void Init(FWfDecl const &WfDecl_, FBasisSet *pOrbBasis_, FAtomSet const &Atoms_, FHfOptions const &Options_, FMemoryStack &Mem);
   void AccFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FMemoryStack &Mem);
   void AccGradient(FMatrixView Gradient, FMatrixView COccC, FMatrixView COccO, FMemoryStack &Mem);
   ~FFockComponentBuilderDfCoulXcCached();
   void PrintEnergyContribs();
   bool SwitchToRefineGrid(FDftGridParams const &NewGridParams);
protected:
   std::string
      XcFunctionalName;
   FDfJkOptions
      JkOptions;
   FAtomSet const
      *pAtoms;
   FBasisSetPtr
      pFitBasis;
   FBasisSet
      *pOrbBasis;
   FRawBasis
      *pOrbBasisRaw,
      *pFitBasisRaw;
   FMatrixView
      Jcd;
   size_t
      nAo,
      nFit;
   FDftGridParams
      GridParams;
   FDftGridPtr
      pDftGrid;
   FXcFunctionalPtr
      pXcFn;
   TArray<double>
      // format: nFit x [(nAo * (nAo+1))/2].
      // This version is for *small* cases and stores *all* integrals.
      Int3ixStorage;
   FMatrixView
      Int3ix;
//    FFockComponentBuilderDfCoulXcCachedImpl
//       *p;
   TArray<double>
      m_LastDensity;
   bool
      AuxExpandXc;

   double
      EnergyCoulomb, EnergyXc, fElecTotal, fElecTotalAnalytic;
   void Make1ixDensity(FMatrixView Jgamma, FMatrixView AuxDen, FMatrixView COccC, FMatrixView COccO, FMemoryStack &Mem);
private:
   FFockComponentBuilderDfCoulXcCached(FFockComponentBuilderDfCoulXcCached const&); // not implemented
   void operator = (FFockComponentBuilderDfCoulXcCached const &); // not implemented
};




// build IAO local exchange
struct FFockComponentBuilderDfLx : public FFockComponentBuilder
{
   FFockComponentBuilderDfLx(FDfJkOptions const &Options_, FMatrixView S1_, FMatrixView S1cd_, FLog &Log_);
   void Init(FWfDecl const &WfDecl_, FBasisSet *pOrbBasis_, FAtomSet const &Atoms_, FHfOptions const &Options_, FMemoryStack &Mem);
   void AccFock(FMatrixView &FockC, FMatrixView &FockO, FBasisSet *pBasis, FMatrixView const &COccC, FMatrixView const &COccO, uint Flags, FMemoryStack &Mem);
   ~FFockComponentBuilderDfLx();
public:
   FDfJkOptions
      Options;
   FAtomSet const
      *pAtoms;
   FBasisSetPtr
      pFitBasis,
      pMinBasis,
      pIaoBasis;
   FBasisSet
      *pOrbBasis;
   FMatrixView
      Jcd,  // cholesky decomposition of full JKFIT basis (may not be required, depending on the mode of use)
      S1,   // overlap matrix of orbital basis
      S1cd; // choleksy decomposition of orbital basis.
};


double *FormIntMNF(ir::FRawShell const &ShF,
   ir::FIntegralKernel *pIntKernel, FRawBasis const *pOrbBasisA, FRawBasis const *pOrbBasisB, FRawBasis const *pFitBasis, FMemoryStack &Mem, FMatrixView ScrDen, double fThr);
double *FormIntMNF(ir::FRawShell const &ShF,
   ir::FIntegralKernel *pIntKernel, FRawBasis const *pOrbBasis, FRawBasis const *pFitBasis, FMemoryStack &Mem, FMatrixView ScrDen, double fThr);


} // namespace ct

#endif // CT_RHF_FOCKBUILD_H
