#ifndef grad_h
#define grad_h
/** \file grad.h
 * Header file for \ref grad.c
 * 
 * Contains declarations of the functions implemented in \ref grad.c and
 * enumeration of the various GradientTypes which is used in the Gradient
 * structure herein discerning the various stages of finding the conjugate
 * gradient direction.
 * 
 * 
  Project: ParallelCarParrinello
  Jan Hamaekers
  2000

  File: grad.h
  $Id: grad.h,v 1.21 2007/02/09 09:22:05 foo Exp $
*/


// use double precision fft when we have it
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef HAVE_DFFTW_H
#include "dfftw.h"
#else
#include "fftw.h"
#endif


#define MaxGradientArrayTypes 9 //!< number of different stages in the conjugate direction
/// Enumerating the different Gradient types during the construction of the conjugate gradient
enum GradientTypes { GraSchGradient,    //!< set onto extra local Psi in order to orthonormalise CG by GramSch() (working array)
                     ActualGradient,    //!< contains orthonormalized direction of steepest descent of step m \f$\zeta^{(m)}_i = - (H-\lambda_i)\psi^{(m)}_i\f$
                     OldActualGradient, //!< contains orthonormalized direction of steepest descent of step m-1 \f$\zeta^{(m-1)}_i = - (H-\lambda_i)\psi^{(m-1)}_i\f$
                     HcGradient,        //!< contains gradient of wave function \f$\langle \chi_G| -\frac{1}{2}\nabla^2 + \underbrace{V^H + V^{ps,loc} + V^{XC}}_{V^{loc}} + V^{ps,nl} | \psi_i \rangle\f$
                     H1cGradient,        //!< contains gradient of wave function \f$\langle {\cal H}^{(1)} | \psi_i^{(0)} \rangle\f$
                     PreConGradient,    //!< contains preconditioned, orthonormalized direction of steepest descent \f$M \tilde{\zeta}_i\f$
                     ConDirGradient,    //!< conjugate gradient direction of current step m, \f$\tilde{\tilde{\zeta}}_i^{(m)}\f$
                     OldConDirGradient, //!< old conjugate gradient direction of last step m-1, \f$\tilde{\tilde{\zeta}}_i^{(m-1)}\f$
                     TempGradient       //!< temporary working gradient array
                   }; 

/** Structure containing calculated conjugate gradient at its various stages.
 */
struct Gradient {
  fftw_complex *GradientArray[MaxGradientArrayTypes]; //!< Array containing the conjugate direction: initial, steepest descent, preconditioned, conjugate, ...
};

void InitGradients(struct Problem *P, struct Gradient *Grad);
void RemoveGradients(struct Problem *P, struct Gradient *Grad);
void CalculateGradientNoRT(struct Problem *P, fftw_complex *source, const double PsiFactor, double *Lambda, fftw_complex ***fnl, fftw_complex *grad, fftw_complex *oldgrad, fftw_complex *Hc_grad);
void CalculateNewWave(struct Problem *P, double *x);
void CalculateConDirHConDir(struct Problem *P, fftw_complex *ConDir, const double PsiFactor, double *ConDirHConDir, double *HartreeddEddt0, double *XCddEddt0);
double CalculateDeltaI(double E0, double dEdt0, double ddEddt0, double *E, double *dE, double *ddE, double *dcos, double *dsin);
void CalculateHamiltonian(struct Problem *P);
void CalculateUnOccupied(struct Problem *P);
double GradSP(struct Problem *P, struct LatticeLevel *Lev, const fftw_complex *LPsiDatA, const fftw_complex *LPsiDatB);
double GradImSP(struct Problem *P, struct LatticeLevel *Lev, fftw_complex *LPsiDatA, fftw_complex *LPsiDatB);
double UpdateWaveAfterIonMove(struct Problem *P);
#endif