/* 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 IR_H
#define IR_H

#ifndef M_PI
   #define M_PI 3.14159265358979323846
#endif

#ifdef _DEBUG
   #include <ostream> // for debug printing of shell information.
   #include <stdio.h> // for printf
#endif

#include <stddef.h> // for size_t

#ifdef MOLPRO
   #define INCLUDE_OPTIONALS
#endif

#include "CxDefs.h"
#include "CxMemoryStack.h"
#include "IrBoysFn.h"

namespace ir {
   // low level shell data structure the integral drivers work with
   // (for simplifying talking to existing fortran and C++ programs!)
   struct FRawShell
   {
      unsigned
         // angular momentum, number of primitive exponents, number of contractions
         l, nExp, nCo;
      double const
         // pExp[iExp], iExp = 0..nExp-1: Exponent of primitive #i.
         *pExp,
         // pCo[iExp + nExp*iCo]: coefficient of primitive #iExp in contraction #iCo
         //
         // Note: Contraction coefficients are stored with respect to raw,
         // unnormalized Gaussians. Use the RawGaussNorm() function
         // to convert between contraction coefficients in terms of normalized
         // Gaussians (as used in libraries, for example), and raw coefficients
         // (as used in most programs)
         *pCo,
         // pCenter[i]: components 0..2: x,y,z position in space.
         *vCen;
      double const
         // Can be 0. If provided: Screening information.
         // [0]: largest range of any contracted function in the shell,
         // [1+iExp]: range of primitive #iExp
         // [1+nExp+iCo]: range of contraction #iCo
         *pRange;

      FRawShell() {}
      FRawShell(unsigned l_, double const *pExp_, unsigned nExp_, double const *pCo_, unsigned nCo_,
                double const *pvCen_, double const *pRange_=0)
         : l(l_), nExp(nExp_), nCo(nCo_), pExp(pExp_), pCo(pCo_),
           vCen(pvCen_), pRange(pRange_)
      {}

      // return number of spherical components per contraction
      inline unsigned nSh() const { return (2*l+1); }
      // return number of functions represented by the shell
      inline unsigned nFn() const { return (2*l+1) * nCo; }
      // return contraction coefficient of primitive iExp in contraction iCo.
      inline double fCo(unsigned iExp, unsigned iCo) const { assert(iCo < nCo && iExp < nExp); return pCo[iExp + nExp*iCo]; }

      inline double MaxCoRange() const { assert(pRange); return pRange[0]; }
      inline double CoRange(unsigned iCo) const { assert(pRange && iCo < nCo); return pRange[1+nExp+iCo]; }
      inline double ExpRange(unsigned iExp) const { assert(pRange && iExp < nExp); return pRange[1+iExp]; }
   };


   // integral kernels represent scalar functions G(r_12), to be used in matrix
   // elements like \int A(r_1) G(r_1 - r_2) C(r_2) d^3r_1 d^3r_2.
   struct FIntegralKernel
   {
      virtual double MaxRange() const = 0;
      // evaluate Gm(rho,T) for m = 0...MaxM  (MaxM included), where T = rho |P-Q|^2.
      // All output data is multiplied by Factor.
      virtual void EvalGm(double *pOut, double rho, double T, unsigned MaxM, double Factor) const = 0;
      virtual ~FIntegralKernel() = 0;
   };

   // Kernel for G(r_12) = delta(r_12)
   struct FOverlapKernel : public FIntegralKernel
   {
      double MaxRange() const; // override
      void EvalGm(double *pOut, double rho, double T, unsigned MaxM, double Factor) const; // override
      ~FOverlapKernel();
   };

   // Kernel for G(r_12) = 1/r_12
   struct FCoulombKernel : public FIntegralKernel
   {
      double MaxRange() const; // override
      void EvalGm(double *pOut, double rho, double T, unsigned MaxM, double Factor) const; // override
      ~FCoulombKernel();
   };


   extern FOverlapKernel g_IrOverlapKernel;
   extern FCoulombKernel g_IrCoulombKernel;

   // evaluate 2-electron, 3-center integrals
   //   (ab|kernel|c)
   void EvalInt2e3c(double *pOut, size_t *Strides,
      FRawShell const *pA, FRawShell const *pB, FRawShell const *pCs, size_t nC,
      double Prefactor, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);
   // evaluate c-contraction of 2-electron, 3-center integral (ab|kernel|c)
   //   Out[a,b,V] := \sum_c (ab|kernel|c) pIn[c,V]
   void EvalInt2e3c_ContractC(double *pOut, size_t *Strides, double const *pIn, size_t StrideIn, size_t nVecIn,
      FRawShell const *pA, FRawShell const *pB, FRawShell const *pCs, size_t nC,
      double Prefactor, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);

   // evaluate 2-electron, 2-center integrals
   //   (a|kernel|c)
   void EvalInt2e2c(double *pOut, size_t StrideA, size_t StrideC, FRawShell const *pA, FRawShell const *pC,
      double Prefactor, bool Add, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);
   // evaluate 2-electron, 2-center integrals
   //   (a|kernel Laplace^n|c)
   void EvalInt2e2c_LaplaceC(double *pOut, size_t StrideA, size_t StrideC, FRawShell const *pA, FRawShell const *pC,
      double Prefactor, bool Add, unsigned LaplaceOrder, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);
   // evaluate 1st derivative (with respect to atomic coordinates) of 2-electron, 2-center integrals:
   //   d/d[Eu] (a|kernel Laplace^n|c)
   // (E=A,C, u=x,y,z)
   void EvalInt2e2c1d_LaplaceC(double *pOutAxyz, double *pOutCxyz, size_t StrideA, size_t StrideC, size_t StrideDeriv, FRawShell const *pA, FRawShell const *pC,
      double Prefactor, bool Add, unsigned LaplaceOrder, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);

   // Compute 1st derivative integrals
   //
   //       d/d[Eu]  (a b|kernel|c)
   //
   // Where Eu (E=A,B, u=x,y,z) denotes the positions of the atomic centers on
   // which the basis functions pA and pB sit.
   //
   // Notes:
   // - The output of the derivative components is organized as
   //
   //       pOut[... + iComp*Strides[3]]
   //
   //   with iComp = 0...5 in order: d/d[Ax], d/d[Ay], d/d[Az], d/d[Bx], d/d[By], d/d[Bz].
   //
   // - This function *DOES NOT* explicitly compute the derivatives with respect
   //   to C. These can be obtained from the A and B derivatives via
   //   translational invariance:
   //
   //       (d/d[Ax] + d/d[Bx] + d/d[Cx]) I = 0.
   //
   // - TODO: consider extending this function to indirect addressing of output
   //   components by passing eiter iCenterA/iCenterB as input or a double
   //   **pOut denoting the output addresses directly. Currently I always
   //   seem to be storing stuff to buffers and munching on it afterwards.
   void EvalInt2e3c1d(double *pOut, size_t *Strides,
      FRawShell const *pA, FRawShell const *pB, FRawShell const *pCs, size_t nC,
      double Prefactor, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);
   void EvalInt2e3c2d(double *pOut, size_t *Strides,
      FRawShell const *pA, FRawShell const *pB, FRawShell const *pCs, size_t nC,
      double Prefactor, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);
   // evaluate laplace op commutator integral:  ([Nabla a]b|c) - (a [Nabla b]|c)
   void EvalInt2e3c_kcomm(double *pOut, size_t *Strides,
      FRawShell const *pA, FRawShell const *pB, FRawShell const *pCs, size_t nC,
      double Prefactor, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);
   void EvalInt2e3c1d_ContractC(double *pOut, size_t *Strides, double const *pIn, size_t StrideIn, size_t nVecIn,
      FRawShell const *pA, FRawShell const *pB, FRawShell const *pCs, size_t nC,
      double Prefactor, FIntegralKernel const *pKernel, ct::FMemoryStack &Mem);

   // TODO: make a function with similar interface which will work if functions are not all on the same center.
   // pBfPos then goes away. Screening interface could also be adjusted to IR-style.
   void EvalShellGroupOnGrid(double *pOut, size_t nCompStride, size_t nBfStride,
      size_t *pMap, size_t &nMap, size_t iFnBase,
      double const *pGridPos, size_t GridStride, size_t nGridPts,
      FRawShell const *pFirstBf, FRawShell const *pLastBf, double const *pBfPos,
      unsigned DerivOrder, double ThrOrb, double LogThrOrb, ct::FMemoryStack &Mem);

   // For a given 3x3 unitary rotation matrix R[3][3], construct the matrix T[lc,l'c'] which makes
   //    S[l,c](R r) = \sum_{c',l'} T[lc,l'c'] S[l',c'](r)
   // for all 3-vectors 'r' (i.e., compute the transformation between solid harmonic components
   // induced by rotations in real space).
   //
   // Notes:
   //  - The matrix elements T[lc,l'c'] for l != l' vanish.
   //  - T[lc,l'c'] is stored at pOut[iSlcX(l,c) + nSlmX(MaxL) * iSlcX(l',c')].
   //    nStride is set to nSlmX(MaxL); pOut is allocated on Mem by this routine.
   //  - Order of trafo not checked. If it doesn't work, use transpose(T) instead of T.
   void EvalSlcXRotationTrafo(double *&pOut, size_t &nStride, unsigned MaxL, double const *R, ct::FMemoryStack &Mem);

#ifdef _DEBUG
void IrPrintMatrixGen(std::ostream &xout, double *pData, size_t nRows, size_t iRowSt, size_t nCols, size_t iColSt, std::string const &Caption);
void operator << (std::ostream &xout, FRawShell const &rs);
#endif // _DEBUG

} // namespace ir

#endif // IR_H
