/** \file energy.c
 * Energy calculation, especially kinetic and electronic.
 * 
 * Contains routines for output of calculated values EnergyOutput(),
 * initialization of the Energy structure InitEnergyArray() or inital calculations
 * InitPsiEnergyCalculation(),
 * updating UpdateEnergyArray() and UpdatePsiEnergyCalculation() and of course
 * calculation of the various parts: Density CalculateDensityEnergy(), ionic
 * CalculateIonsEnergy(), electronic CalculatePsiEnergy(), kinetic with RiemannTensor
 * CalculateKineticEnergyUseRT() (calling EnergyGNSP()) and without
 * CalculateKineticEnergyNoRT() (using EnergyLapNSP()). CalculateGapEnergy() evaluates
 * would-be addition from unoccupied states to total energy.
 * Smaller functions such as EnergyAllReduce() gather intermediate results from all
 * processes and sum up to total value.  
 * \note Some others are stored in excor.c (exchange energies) and pseudo.c (density
 *       energies)
 * 
 * 
  Project: ParallelCarParrinello
  \author Jan Hamaekers
 \date 2000

  File: energy.c
  $Id: energy.c,v 1.68 2007/02/09 09:13:48 foo Exp $
*/

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include <string.h>

#include "data.h"
#include "errors.h"
#include "helpers.h"
#include "energy.h"
#include "riemann.h"
#include "excor.h"
#include "myfft.h"
#include "mymath.h"
//#include "output.h"
#include "run.h"

/** Calculates kinetic energy for a given wave function in reciprocal base.
 * \param *P Problem at hand
 * \param *Lev LatticeLevel structure
 * \param *LPsiDatA complex wave function coefficients \f$c_{i,G}\f$ (reciprocal base)
 * \param *T kinetic eigenvalue, herein result is additionally stored
 * \param Factor  Additional factor for return value (such as occupation number \f$f_i\f$ )
 * \return \f$E_{kin} = \frac{1}{2} f_i \sum_G |G|^2 |c_{i,G}|^2\f$
 * \sa used in CalculateKineticEnergyNoRT()
 */
static double EnergyGNSP(const struct Problem *P, const struct LatticeLevel *Lev, fftw_complex *LPsiDatA, double *T, const double Factor) {
  double LocalSP=0.0,PsiSP;
  int i,s = 0;
  struct OneGData *G = Lev->GArray;
  if (Lev->GArray[0].GSq == 0.0) {  // only for continuity (see other functions)
    LocalSP += G[0].GSq*LPsiDatA[0].re*LPsiDatA[0].re;
    s++;
  }
  /// Performs complex product for each reciprocal grid vector times twice this vector's norm.
  for (i=s; i < Lev->MaxG; i++) {
    LocalSP += G[i].GSq*2.*(LPsiDatA[i].re*LPsiDatA[i].re+LPsiDatA[i].im*LPsiDatA[i].im);  // factor 2 due to gamma point symmetry
  }
  /// Gathers partial results by summation from all other coefficients sharing processes.
  MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  /// Multiplies by \a Factor and returns result.
  *T = Factor*PsiSP;
  return(Factor*PsiSP);
}

/** Calculates kinetic energy for a given wave function in reciprocal base in case of Riemann tensor use.
 * \param *P Problem at hand
 * \param *RT RiemannTensor structure
 * \param *src replaces reciprocal grid vectors
 * \param *Lap replaces complex Psi wave function coefficients
 * \param *T kinetic eigenvalue, herein result is additionally stored
 * \param Factor Additional factor for return value
 * \sa used in CalculateKineticEnergyUseRT(), compare with EnergyGNSP()
 */
static double EnergyLapNSP(const struct Problem *P, const struct RiemannTensor *RT, fftw_complex *src, fftw_complex *Lap, double *T, const double Factor) {
  double LocalSP=0.0,PsiSP;
  int i;
  int s = 0;
  struct LatticeLevel *Lev = RT->LevS;
  if (Lev->GArray[0].GSq == 0.0) {
    LocalSP += src[0].re*Lap[0].re+src[0].im*Lap[0].im;
    s++;
  }
  for (i=s; i < Lev->MaxG; i++) {
    LocalSP += 2.*(src[i].re*Lap[i].re+src[i].im*Lap[i].im);
  }
  MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  *T = Factor*PsiSP;
  return(Factor*PsiSP);
}

/** Calculcates Kinetic energy without Riemann tensor.
 * \param *P Problem at hand
 * \param p local psi wave function number
 * \sa EnergyGNSP()
 */
static void CalculateKineticEnergyNoRT(struct Problem *P, const int p) 
{
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = Lat->E;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct LevelPsi *LPsi = LevS->LPsi;
  struct Psis *Psi = &Lat->Psi;
  const double PsiFactor = Psi->LocalPsiStatus[p].PsiFactor;
  fftw_complex *psi;
  psi = LPsi->LocalPsi[p];
  E->PsiEnergy[KineticEnergy][p] = EnergyGNSP(P, LevS, psi, &Psi->AddData[p].T, 0.5*PsiFactor);
}

/** Calculcates Kinetic energy with Riemann tensor.
 * CalculateRiemannLaplaceS() is called and Lattice->RT.RTLaplaceS used instead
 * of Lattice::Psi.
 * \param *P Problem at hand
 * \param p local psi wave function number
 * \sa EnergyLapNSP()
 */
static void CalculateKineticEnergyUseRT(struct Problem *P, const int p)
{
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = Lat->E;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct LevelPsi *LPsi = LevS->LPsi;
  struct Psis *Psi = &Lat->Psi;
  const double PsiFactor = Psi->LocalPsiStatus[p].PsiFactor;
  fftw_complex *psi;
  psi = LPsi->LocalPsi[p];
  CalculateRiemannLaplaceS(Lat, &Lat->RT, psi, Lat->RT.RTLaplaceS);
  E->PsiEnergy[KineticEnergy][p] = EnergyLapNSP(P, &Lat->RT, psi, Lat->RT.RTLaplaceS, &Psi->AddData[p].T, 0.5*PsiFactor);
  //fprintf(stderr,"CalculateKineticEnergyUseRT: PsiEnergy[KineticEnergy] = %lg\n",E->PsiEnergy[KineticEnergy][p]);
}

/** Sums up calculated energies from all processes, calculating Energy::TotalEnergy.
 * Via MPI_Allreduce is Energy::AllLocalPsiEnergy summed up and the sum
 * redistributed back to all processes for the various spin cases. In
 * case of SpinType#SpinUp or SpinType#SpinDown the "other" energy (e.g. down
 * for an up process) is via MPI_SendRecv communicated, and lateron
 * Energy#AllTotalPsiEnergy as the sum of two states generated in each process.<br>
 * Finally, the Energy::TotalEnergy[0] is the sum of Energy::AllTotalPsiEnergy[],
 * Energy::AllTotalDensityEnergy[], - Energy::AllTotalIonsEnergy[GaussSelfEnergy]
 * and Energy::AllTotalIonsEnergy[EwaldEnergy].
 * Again, as in UpdateEnergyArray(), this is split up completely for PsiTypeTag#UnOccupied
 * and PsiTypeTag#Occupied.
 * \param *P Problem at hand
 */
void EnergyAllReduce(struct Problem *P)
{
  struct Psis *Psi = &P->Lat.Psi;
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = P->Lat.E;
  struct OnePsiElement *OnePsiB;
  MPI_Status status;
  int i,j,m;
  double lambda, overlap, lambda2, overlap2; //, part1, part2, tmp;
  MPI_Allreduce(E->AllLocalDensityEnergy, E->AllTotalDensityEnergy, MAXALLDENSITYENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi );
  switch (Psi->PsiST) {
  case SpinDouble:
    MPI_Allreduce(E->AllLocalPsiEnergy, E->AllTotalPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); 
    break;
  case SpinUp:
    MPI_Allreduce(E->AllLocalPsiEnergy, E->AllUpPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); 
    MPI_Sendrecv(E->AllUpPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1,
		 E->AllDownPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2,
		 P->Par.comm_STInter, &status );
    for (i=0; i< MAXALLPSIENERGY; i++)
      E->AllTotalPsiEnergy[i] = E->AllUpPsiEnergy[i] + E->AllDownPsiEnergy[i];
    break;
  case SpinDown:
    MPI_Allreduce(E->AllLocalPsiEnergy, E->AllDownPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Sendrecv(E->AllDownPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2,
		 E->AllUpPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1,
		 P->Par.comm_STInter, &status );
    for (i=0; i< MAXALLPSIENERGY; i++)
      E->AllTotalPsiEnergy[i] = E->AllUpPsiEnergy[i] + E->AllDownPsiEnergy[i];
    break;
  }
  switch (R->CurrentMin) {
    case Occupied:
      E->TotalEnergy[0] = E->AllTotalPsiEnergy[KineticEnergy];
      E->TotalEnergy[0] += E->AllTotalPsiEnergy[NonLocalEnergy]
      ;
      E->TotalEnergy[0] += E->AllTotalDensityEnergy[CorrelationEnergy]+E->AllTotalDensityEnergy[ExchangeEnergy];
      E->TotalEnergy[0] += E->AllTotalDensityEnergy[PseudoEnergy];
      E->TotalEnergy[0] += E->AllTotalDensityEnergy[HartreePotentialEnergy];
      E->TotalEnergy[0] -= E->AllTotalIonsEnergy[GaussSelfEnergy];
      E->TotalEnergy[0] += E->AllTotalIonsEnergy[EwaldEnergy];
//      fprintf(stderr,"(%i) Tot %1.5f ---\t Diff %1.5f\t ---\t Kin %1.5f\tNLErg %1.5f ---\t Corr %1.5f\tXC %1.5f\t Pseudo %1.5f\tHP %1.5f\n", P->Par.me,
//        E->TotalEnergy[0],(E->TotalEnergy[1]-E->TotalEnergy[0])/fabs(E->TotalEnergy[0]),
//        E->AllTotalPsiEnergy[KineticEnergy],E->AllTotalPsiEnergy[NonLocalEnergy],
//        E->AllTotalDensityEnergy[CorrelationEnergy],E->AllTotalDensityEnergy[ExchangeEnergy],
//        E->AllTotalDensityEnergy[PseudoEnergy],E->AllTotalDensityEnergy[HartreePotentialEnergy]);
      break;
    case UnOccupied:
      E->TotalEnergy[0]  = E->AllTotalPsiEnergy[KineticEnergy]+E->AllTotalPsiEnergy[NonLocalEnergy];
      E->TotalEnergy[0] += E->AllTotalDensityEnergy[CorrelationEnergy]+E->AllTotalDensityEnergy[ExchangeEnergy];
      E->TotalEnergy[0] += E->AllTotalDensityEnergy[PseudoEnergy]+E->AllTotalDensityEnergy[HartreePotentialEnergy];
      break;
    case Perturbed_P0:
    case Perturbed_P1:
    case Perturbed_P2:
    case Perturbed_RxP0:
    case Perturbed_RxP1:
    case Perturbed_RxP2:
      E->TotalEnergy[0]  = E->AllTotalPsiEnergy[Perturbed1_0Energy]+E->AllTotalPsiEnergy[Perturbed0_1Energy];
      E->parts[0]  = E->AllTotalPsiEnergy[Perturbed1_0Energy]+E->AllTotalPsiEnergy[Perturbed0_1Energy];
      //fprintf(stderr,"(%i) summed single: AllTotalPsiEnergy[total] = %lg + %lg = %lg\n",P->Par.me,E->AllTotalPsiEnergy[Perturbed1_0Energy],E->AllTotalPsiEnergy[Perturbed0_1Energy], E->TotalEnergy[0]);
      //E->TotalEnergy[0] = 0.;
      m = -1;
      E->parts[1] = E->parts[2] = 0.;
      for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) {  // go through all wave functions: \sum_k
        OnePsiB = &Psi->AllPsiStatus[i];    // grab OnePsiB
        if (OnePsiB->PsiType == R->CurrentMin) {   // drop all but occupied types
          m++;  // increase m if it is occupied wave function
          lambda = overlap = 0.;
          if (OnePsiB->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // local?
            //fprintf(stderr,"(%i) Lambda[%i] = %lg\n",P->Par.me, m,Lat->Psi.AddData[i].Lambda);
            lambda = Lat->Psi.AddData[OnePsiB->MyLocalNo].Lambda * Lat->Psi.LocalPsiStatus[OnePsiB->MyLocalNo].PsiFactor;
            for (j=0;j<Psi->NoOfPsis;j++) { // over all Psis: \sum_l
              //lambda -= Lat->Psi.lambda[j][i - Psi->TypeStartIndex[R->CurrentMin]]*GradSP(P,R->LevS,R->LevS->LPsi->LocalPsi[j + Psi->TypeStartIndex[R->CurrentMin]],R->LevS->LPsi->LocalPsi[i]);
              overlap -= Lat->Psi.lambda[j][m]*Lat->Psi.Overlap[j][m];
            }
            // send results to other processes
          }
          //fprintf(stderr,"(%i) before MPI: \tlambda[%i] = %lg\t overlap[%i] = %lg\n",P->Par.me,m,lambda,m,overlap);
          MPI_Allreduce( &lambda, &lambda2, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); // lambda is just local
          MPI_Allreduce( &overlap, &overlap2, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); // we just summed up over local m's!
          //fprintf(stderr,"(%i) after MPI: \t lambda[%i] = %lg\t overlap[%i] = %lg\n",P->Par.me, OnePsiB->MyGlobalNo,lambda2, m, overlap2);
          E->TotalEnergy[0] += lambda2;
          E->TotalEnergy[0] += overlap2;
          E->parts[1] += lambda2;
          E->parts[2] += overlap2;
        }
      }
      
      switch(Psi->PsiST) {
        default:
            lambda = overlap = 0.;
            break;
        case SpinUp:
            MPI_Sendrecv(&E->parts[1], 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1,
             &lambda, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2, P->Par.comm_STInter, &status );
            E->TotalEnergy[0] += lambda;
            MPI_Sendrecv(&E->parts[2],1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
             &overlap, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4, P->Par.comm_STInter, &status );
            E->TotalEnergy[0] += overlap;
            break;
        case SpinDown:
            MPI_Sendrecv(&E->parts[1], 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2,
             &lambda, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1, P->Par.comm_STInter, &status );
            E->TotalEnergy[0] += lambda;
            MPI_Sendrecv(&E->parts[2],1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
             &overlap, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3, P->Par.comm_STInter, &status );
            E->TotalEnergy[0] += overlap;
            break;
      }
      //fprintf(stderr,"(%i) summed single: lambda[total] = %lg\t overlap[total] = %lg\n",P->Par.me,E->parts[1]+lambda, E->parts[2]+overlap);
      break;
    default:
      break;
  }
}

/** Initializes Energy Array.
 * Arrays are malloc'd and set to zero, then InitExchangeCorrelationEnergy()
 * called.
 * \param *P Problem at hand
 */
void InitEnergyArray(struct Problem *P) 
{
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = Lat->E;
  struct Psis *Psi = &Lat->Psi;
  int i, type, oldtype = P->R.CurrentMin;
  for (type=Occupied;type<Extra;type++) {
    SetCurrentMinState(P, type);
    for (i=0; i<MAXALLPSIENERGY; i++) {
      Lat->Energy[type].PsiEnergy[i] = (double *)
        Malloc(sizeof(double)*Psi->LocalNo,"InitEnergyCalculations: PsiKineticEnergy");
      SetArrayToDouble0(Lat->Energy[type].PsiEnergy[i], Psi->LocalNo);
    }
    SetArrayToDouble0(E->AllLocalDensityEnergy, MAXALLDENSITYENERGY);
    SetArrayToDouble0(E->AllTotalDensityEnergy, MAXALLDENSITYENERGY);
    SetArrayToDouble0(E->AllTotalIonsEnergy, MAXALLIONSENERGY);
    SetArrayToDouble0(E->TotalEnergy, MAXOLD);
  }
  InitExchangeCorrelationEnergy(P, &P->ExCo);
  SetCurrentMinState(P, oldtype);
}

/** Shifts stored total energy values up to MAXOLD.
 * Energy::PsiEnergy, Energy::AllLocalDensityEnergy, Energy::AllTotalDensityEnergy
 * are all set to zero. Energy::TotalEnergy is shifted by one to make
 * place for next value. This is done either for occupied or for the gap energies,
 * stated by \a IncType.
 * \param *P Problem at hand
 */
void UpdateEnergyArray(struct Problem *P) 
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct Energy *E = Lat->E;
  int i;
  for (i=0; i<MAXALLPSIENERGY; i++)
    E->PsiEnergy[i][R->OldActualLocalPsiNo] = 0.0;
  SetArrayToDouble0(E->AllLocalDensityEnergy, MAXALLDENSITYENERGY);
  SetArrayToDouble0(E->AllTotalDensityEnergy, MAXALLDENSITYENERGY);
  for (i=MAXOLD-1; i > 0; i--)
    Lat->Energy[R->CurrentMin].TotalEnergy[i] = Lat->Energy[R->CurrentMin].TotalEnergy[i-1];
}

/** Calculates ion energy.
 * CalculateGaussSelfEnergyNoRT() is called.
 * \param *P Problem at hand
 * \warning CalculateIonsUseRT not implemented
 */
void CalculateIonsEnergy(struct Problem *P)
{
  struct Lattice *Lat = &P->Lat;
  switch (Lat->RT.ActualUse) {
  case inactive:
  case standby:
    CalculateGaussSelfEnergyNoRT(P, &P->PP, &P->Ion);
    break;
  case active:
    Error(SomeError, "CalculateIonsUseRT: Not implemented");
    break;
  }
}

/** Calculates density energy.
 * Depending on RT::ActualUse CalculateXCEnergyNoRT(), CalculateSomeEnergyAndHGNoRT() 
 * and CalculateGapEnergy() are called within SpeedMeasure() brackets.
 * Afterwards, CalculateGradientNoRT() for RunStruct::ActualLocalPsiNo.
 * \param *P Problem at hand
 * \param UseInitTime whether SpeedMeasure uses InitLocTime (Init) or LocTime (Update)
 * \warning CalculateSomeEnergyAndHGUseRT not implemented
 * \warning call EnergyAllReduce() afterwards!
 */
void CalculateDensityEnergy(struct Problem *P, const int UseInitTime)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  switch (Lat->RT.ActualUse) {
  case inactive:
  case standby:
    SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
    CalculateXCEnergyNoRT(P);
    CalculateSomeEnergyAndHGNoRT(P);
    SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
    CalculateGradientNoRT(P, 
			  LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], 
			  Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor,
			  &Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda, 
			  P->PP.fnl[R->ActualLocalPsiNo], 
			  P->Grad.GradientArray[ActualGradient],  
        P->Grad.GradientArray[OldActualGradient],  
			  P->Grad.GradientArray[HcGradient]);
      break;
  case active:
    CalculateXCEnergyUseRT(P);
    /* FIXME - Hier fehlt noch was */
    Error(SomeError, "CalculateSomeEnergyAndHGUseRT: Not implemented");
    break;
  }
  /* Achtung danach noch EnergyAllReduce aufrufen */
}

/** Calculation of an additional total energy from unoccupied states.
 * By including additional (unoccupied) orbitals into the minimisation process,
 * which also stand orthogonal on all others, however are minimised in the field
 * of occupied ones (without altering it in a self-consistent manner), a first
 * approximation to the band gap energy can be obtained. The density part of this
 * virtual energy is stored in the PsiTypeTag#UnOccupied part of Lattice#Energy.
 * Note that the band gap is approximated via the eigenvalues of the additional
 * orbitals. This energy here is purely fictious.
 * \param *P Problem at hand
 * \note No RiemannTensor use so far implemented. And mostly copy&paste from
 *        CalculateXCEnergyNoRT() and CalculateSomeEnergyAndHGNoRT().
 * \bug Understand CoreCorrection part in \f$E_{XC}\f$
 */
void CalculateGapEnergy(struct Problem *P) {
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct Energy *E = Lat->E;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  struct Density *Dens = Lev0->Dens;
  struct ExCor *EC = &P->ExCo;
  struct Ions *I = &P->Ion;
  double SumEc = 0;
  double rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx=0.0;
  double Factor = R->XCEnergyFactor/Lev0->MaxN;
  fftw_complex vp,rhoe,rhoe_gap;
  fftw_complex rp,rhog;
  double SumEHP,Fac,SumEH,SumEG,SumEPS;
  int g,s=0,it,Index;
  int i;
  fftw_complex *HG = Dens->DensityCArray[HGDensity];
  SetArrayToDouble0((double *)HG,2*Dens->TotalSize);

  SpeedMeasure(P, GapTime, StartTimeDo);
  // === Calculate exchange and correlation energy ===
  for (i = 0; i < Dens->LocalSizeR; i++) {  // for each node in radial mesh
    // put (corecorrected) densities in p, pUp, pDown  
    p = Dens->DensityArray[TotalDensity][i];
    switch (Psi->PsiST) {
    case SpinDouble:
      pUp = 0.5*p; 
      pDown = 0.5*p;
      break;
    case SpinUp:  
    case SpinDown:
      pUp = Dens->DensityArray[TotalUpDensity][i];
      pDown = Dens->DensityArray[TotalDownDensity][i];
      break;
    default:
      ;
    }
    if ((p < 0) || (pUp < 0) || (pDown < 0)) {
      /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
      p = 0.0; 
      pUp = 0.0;
      pDown = 0.0;
    }
    rs = Calcrs(EC,p);
    rsUp = Calcrs(EC,pUp);
    rsDown = Calcrs(EC,pDown);
    zeta = CalcZeta(EC,pUp,pDown);
    // Calculation with the densities and summation
    SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown);
    SumEc += CalcSECr(EC, rs, zeta, p);
  }
  E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
  E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx;

  // === Calculate electrostatic and local pseudopotential energy ===
  SumEHP = 0.0;
  SumEH = 0.0;
  SumEG = 0.0;
  SumEPS = 0.0;
  // calculates energy of local pseudopotential
  if (Lev0->GArray[0].GSq == 0.0) {
    Index = Lev0->GArray[0].Index;
    c_re(vp) = 0.0;
    c_im(vp) = 0.0;
    for (it = 0; it < I->Max_Types; it++) {
      c_re(vp) += (c_re(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]); 
      c_im(vp) += (c_im(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
    }
    SumEPS += (c_re(Dens->DensityCArray[GapDensity][Index])*c_re(vp) +
         c_im(Dens->DensityCArray[GapDensity][Index])*c_im(vp))*R->HGcFactor; 
    s++;
  }
  for (g=s; g < Lev0->MaxG; g++) {
    Index = Lev0->GArray[g].Index;
    c_re(vp) = 0.0;
    c_im(vp) = 0.0;
    c_re(rp) = 0.0;
    c_im(rp) = 0.0;
    Fac = 2.*PI/(Lev0->GArray[g].GSq);  // 4*.. -> 2*...: 1/2 here ! (due to lacking dependence of first n^{tot} (r) on gap psis)
    for (it = 0; it < I->Max_Types; it++) {
      c_re(vp) += (c_re(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]); 
      c_im(vp) += (c_im(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
      c_re(rp) += (c_re(I->I[it].SFactor[g])*PP->FacGauss[it][g]); 
      c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
    }
    c_re(rhog) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp);
    c_im(rhog) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp);
    c_re(rhoe) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
    c_im(rhoe) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
    c_re(rhoe_gap) = c_re(Dens->DensityCArray[GapDensity][Index])*R->HGcFactor;
    c_im(rhoe_gap) = c_im(Dens->DensityCArray[GapDensity][Index])*R->HGcFactor;
    //c_re(PP->VCoulombc[g]) = Fac*c_re(rhog);
    //c_im(PP->VCoulombc[g]) = -Fac*c_im(rhog);
    c_re(HG[Index]) = c_re(vp) + Fac * (c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp));
    c_im(HG[Index]) = c_im(vp) + Fac * (c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp));
    /*if (P->first) */
    //SumEG += Fac*(c_re(rp)*c_re(rp)+c_im(rp)*c_im(rp));             // Coulomb energy
    // changed: took out gaussian part! rhog -> rhoe
    SumEHP += 2*Fac*(c_re(rhoe)*c_re(rhoe_gap)+c_im(rhoe)*c_im(rhoe_gap));    // E_ES first part
    //SumEH += Fac*(c_re(rhoe)*c_re(rhoe)+c_im(rhoe)*c_im(rhoe));
    SumEPS += 2.*(c_re(Dens->DensityCArray[GapDensity][Index])*c_re(vp)+
      c_im(Dens->DensityCArray[GapDensity][Index])*c_im(vp))*R->HGcFactor;
  }
  //
  for (i=0; i<Lev0->MaxDoubleG; i++) {
    HG[Lev0->DoubleG[2*i+1]].re =  HG[Lev0->DoubleG[2*i]].re;
    HG[Lev0->DoubleG[2*i+1]].im = -HG[Lev0->DoubleG[2*i]].im;
  }
  /*if (P->first) */
  //E->AllLocalDensityEnergy[GaussEnergy] = SumEG*Lat->Volume;
  E->AllLocalDensityEnergy[PseudoEnergy] = SumEPS*Lat->Volume;   // local pseudopotential
  E->AllLocalDensityEnergy[HartreePotentialEnergy] = SumEHP*Lat->Volume;   // Gauss n^Gauss and Hartree n potential
  //E->AllLocalDensityEnergy[HartreeEnergy] = SumEH*Lat->Volume;
  SpeedMeasure(P, GapTime, StopTimeDo); 
  
  // === Calculate Gradient ===
  CalculateGradientNoRT(P, 
      LevS->LPsi->LocalPsi[R->ActualLocalPsiNo], 
      Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor,
      &Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda, 
      P->PP.fnl[R->ActualLocalPsiNo], 
      P->Grad.GradientArray[ActualGradient],  
      P->Grad.GradientArray[OldActualGradient],  
      P->Grad.GradientArray[HcGradient]);
}

/** Calculates the energy of a wave function psi.
 * Calculates kinetc and non local energies, taking care of possible
 * Riemann tensor usage.
 * \param *P Problem at hand
 * \param p number of wave function Psi
 * \param UseInitTime Whether speed measuring uses InitLocTime (Init) or LocTime (Update)
 * \sa EnergyGNSP(), EnergyLapNSP()
 * \warning call EnergyAllReduce() afterwards!
 */
void CalculatePsiEnergy(struct Problem *P, const int p, const int UseInitTime)
{
  struct Lattice *Lat = &(P->Lat);
  switch (Lat->RT.ActualUse) {
  case inactive:
  case standby:
    SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
    CalculateKineticEnergyNoRT(P, p);
    SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
    SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StartTimeDo);
    CalculateNonLocalEnergyNoRT(P, p);
    SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StopTimeDo);
    break;
  case active:
    SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
    CalculateKineticEnergyUseRT(P, p);
    SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
    SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StartTimeDo);
    CalculateNonLocalEnergyUseRT(P, p);
    SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StopTimeDo);
    break;
  }
  /* Achtung danach noch EnergyAllReduce aufrufen */
}

/** Initialization of wave function energy calculation.
 * For every wave function CalculatePsiEnergy() is called,
 * Energy::AllLocalPsiEnergy is set to sum of these.
 * \param *P Problem at hand
 * \param type minimisation type PsiTypeTag to calculate psi energy for
 * \warning call EnergyAllReduce() afterwards!
 */
void InitPsiEnergyCalculation(struct Problem *P, enum PsiTypeTag type)
{
  struct Lattice *Lat = &(P->Lat);
  int p,i, oldtype = P->R.CurrentMin;
  SetCurrentMinState(P,type);
  for (p=0; p < Lat->Psi.LocalNo; p++) { 
    if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
      CalculatePsiEnergy(P, p, 1);
  }
  for (i=0; i< MAXALLPSIENERGY; i++) {
    Lat->E->AllLocalPsiEnergy[i] = 0.0;
    for (p=0; p < Lat->Psi.LocalNo; p++)
      if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
        Lat->E->AllLocalPsiEnergy[i] += Lat->E->PsiEnergy[i][p];
  }   
  SetCurrentMinState(P,oldtype);
  /* Achtung danach noch EnergyAllReduce aufrufen */
}

/** Updating wave function energy.
 * CalculatePsiEnergy() is called for RunStruct::OldActualLocalPsiNo,
 * Energy::AllLocalPsiEnergy is set to sum of all Psi energies.
 * \param *P Problem at hand
 * \warning call EnergyAllReduce() afterwards!
 */
void UpdatePsiEnergyCalculation(struct Problem *P)
{
  struct Lattice *Lat = &(P->Lat);
  struct RunStruct *R = &P->R;
  int p = R->OldActualLocalPsiNo, i;
  if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
    CalculatePsiEnergy(P, p, 0);
  for (i=0; i< MAXALLPSIENERGY; i++) {
    Lat->E->AllLocalPsiEnergy[i] = 0.0;
    for (p=0; p < P->Lat.Psi.LocalNo; p++)
      if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
        Lat->E->AllLocalPsiEnergy[i] += Lat->E->PsiEnergy[i][p];
  }   
  /* Achtung danach noch EnergyAllReduce aufrufen */
}

/** Output of energies.
 * Prints either Time, PsiEnergy(Up/Down), DensityEnergy, IonsEnergy, KineticEnergy,
 * TotalEnergy and Temperature to screen or<br>
 * Time, TotalEnergy, NonLocalEnergy, CorrelationEnergy, ExchangeEnergy,
 * PseudoEnergy, HartreePotentialEnergy, GaussSelfEnergy, EwaldEnergy,
 * KineticEnergy and TotalEnergy+KineticEnergy to the energy file
 * \param *P Problem at hand
 * \param doout 1 - print stuff to stderr, 0 - put compact form into energy file
 */
void EnergyOutput(struct Problem *P, int doout)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  struct Ions *I = &P->Ion;
  struct Psis *Psi = &Lat->Psi;
  struct FileData *F = &P->Files;
  FILE *eout = NULL;
  FILE *convergence = NULL;
  int i, it;
  if (P->Call.out[LeaderOut] && (P->Par.me == 0) && doout) {
    fprintf(stderr, "Time T:  :\t%e\n",P->R.t);
    fprintf(stderr, "PsiEnergy    :");
    for (i=0; i<MAXALLPSIENERGY; i++) {
      fprintf(stderr," %e",Lat->Energy[Occupied].AllTotalPsiEnergy[i]);
    }
    fprintf(stderr,"\n");
    if (Psi->PsiST != SpinDouble) {
      fprintf(stderr, "PsiEnergyUp  :");
      for (i=0; i<MAXALLPSIENERGY; i++) {
				fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllUpPsiEnergy[i]);
      }
      fprintf(stderr,"\n");
      fprintf(stderr, "PsiEnergyDown:");
      for (i=0; i<MAXALLPSIENERGY; i++) {
				fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllDownPsiEnergy[i]);
      }
      fprintf(stderr,"\n");
    }
    fprintf(stderr, "DensityEnergy:");
    for (i=0; i<MAXALLDENSITYENERGY; i++) {
      fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllTotalDensityEnergy[i]);
    }
    fprintf(stderr,"\n");
    fprintf(stderr, "IonsEnergy   :");
    for (i=0; i<MAXALLIONSENERGY; i++) {
      fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllTotalIonsEnergy[i]);
    }
    fprintf(stderr,"\n");
    fprintf(stderr, "TotalEnergy  :\t%e\n",Lat->Energy[Occupied].TotalEnergy[0]);
    if (R->DoPerturbation)
      fprintf(stderr, "PerturbedEnergy :\t%e\t%e\t%e\t%e\t%e\t%e\n",Lat->Energy[Perturbed_P0].TotalEnergy[0], Lat->Energy[Perturbed_P1].TotalEnergy[0], Lat->Energy[Perturbed_P2].TotalEnergy[0], Lat->Energy[Perturbed_RxP0].TotalEnergy[0], Lat->Energy[Perturbed_RxP1].TotalEnergy[0], Lat->Energy[Perturbed_RxP2].TotalEnergy[0]);
    if (R->DoUnOccupied) {
      fprintf(stderr, "GapEnergy  :\t%e\t%e\t%e\t%e\n",Lat->Energy[UnOccupied].AllTotalPsiEnergy[KineticEnergy],Lat->Energy[UnOccupied].AllTotalPsiEnergy[NonLocalEnergy], Lat->Energy[UnOccupied].AllTotalDensityEnergy[PseudoEnergy],Lat->Energy[UnOccupied].AllTotalDensityEnergy[HartreePotentialEnergy]);
      fprintf(stderr, "TotalGapEnergy  :\t%e\n",Lat->Energy[UnOccupied].TotalEnergy[0]);
      fprintf(stderr, "BandGap :\t%e\n", Lat->Energy[UnOccupied].bandgap);
    }
    fprintf(stderr, "IonKinEnergy  :\t%e\n",P->Ion.EKin);
    fprintf(stderr, "Temperature  :\t%e\n",P->Ion.ActualTemp);
  }
  if (P->Files.MeOutMes == 0 || doout == 0) return;
  eout = F->EnergyFile;
  if (eout != NULL) {
    if (R->DoUnOccupied) {
      fprintf(eout, "%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",
    	  P->R.t,
    	  Lat->Energy[Occupied].TotalEnergy[0], Lat->Energy[UnOccupied].TotalEnergy[0],
    	  Lat->Energy[Occupied].AllTotalPsiEnergy[KineticEnergy], Lat->Energy[Occupied].AllTotalPsiEnergy[NonLocalEnergy],
        Lat->Energy[UnOccupied].AllTotalPsiEnergy[KineticEnergy]+Lat->Energy[UnOccupied].AllTotalPsiEnergy[NonLocalEnergy],
    	  Lat->Energy[Occupied].AllTotalDensityEnergy[CorrelationEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[ExchangeEnergy],
    	  Lat->Energy[Occupied].AllTotalDensityEnergy[PseudoEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[HartreePotentialEnergy],
        Lat->Energy[UnOccupied].AllTotalDensityEnergy[PseudoEnergy]+Lat->Energy[UnOccupied].AllTotalDensityEnergy[HartreePotentialEnergy],
    	  -Lat->Energy[Occupied].AllTotalIonsEnergy[GaussSelfEnergy],
    	  Lat->Energy[Occupied].AllTotalIonsEnergy[EwaldEnergy],
    	  P->Ion.EKin,
    	  Lat->Energy[Occupied].TotalEnergy[0]+P->Ion.EKin);
    } else {
      fprintf(eout, "%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",
        P->R.t,
        Lat->Energy[Occupied].TotalEnergy[0],
        Lat->Energy[Occupied].AllTotalPsiEnergy[KineticEnergy], Lat->Energy[Occupied].AllTotalPsiEnergy[NonLocalEnergy],
        Lat->Energy[Occupied].AllTotalDensityEnergy[CorrelationEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[ExchangeEnergy],
        Lat->Energy[Occupied].AllTotalDensityEnergy[PseudoEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[HartreePotentialEnergy],
        -Lat->Energy[Occupied].AllTotalIonsEnergy[GaussSelfEnergy],
        Lat->Energy[Occupied].AllTotalIonsEnergy[EwaldEnergy],
        P->Ion.EKin,
        Lat->Energy[Occupied].TotalEnergy[0]+P->Ion.EKin);
    }
    fflush(eout);
  }
  
  // Output for testing convergence
  if (!P->Par.me) {
    if (R->DoPerturbation)  {
      if ((convergence = fopen("pcp.convergence.csv","r")) != NULL) { // only prepend with header if file doesn't exist yet
        fclose(convergence);
        convergence = fopen("pcp.convergence.csv","a");
      } else {
        convergence = fopen("pcp.convergence.csv","w");
        if (convergence != NULL) {
          fprintf(convergence,"Ecut\ta.u.\tMaxG\tN0\tN1\tN2\tChi00\tChi11\tChi22\tP0\tP1\tP2\tRxP0\tRxP1\tRxP2\tTE\t");
          for (it=0; it < I->Max_Types; it++) {  // over all types
            fprintf(convergence,"Sigma00_%s\tSigma11_%s\tSigma22_%s\t",I->I[it].Symbol,I->I[it].Symbol,I->I[it].Symbol);
            fprintf(convergence,"Sigma00_%s_rezi\tSigma11_%s_rezi\tSigma22_%s_rezi\t",I->I[it].Symbol,I->I[it].Symbol,I->I[it].Symbol);
          }
          fprintf(convergence,"Volume\tN^3\tSawStart\t");
          fprintf(convergence,"Perturbedrxp-0\tPerturbedrxp-1\tPerturbedrxp-2\t");
          fprintf(convergence,"Perturbedp-0\tPerturbedp-1\tPerturbedp-2\n");
        }
      }
      if (convergence != NULL) {
        fprintf(convergence,"%e\t%e\t%i\t%i\t%i\t%i\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t",     
          Lat->ECut,
          Lat->RealBasisQ[0], LevS->TotalAllMaxG, Lev0->Plan0.plan->N[0], Lev0->Plan0.plan->N[1], Lev0->Plan0.plan->N[2],      
          I->I[0].chi[0+0*NDIM], I->I[0].chi[1+1*NDIM], I->I[0].chi[2+2*NDIM], 
          Lat->Energy[Perturbed_P0].TotalEnergy[0], Lat->Energy[Perturbed_P1].TotalEnergy[0], Lat->Energy[Perturbed_P2].TotalEnergy[0], 
          Lat->Energy[Perturbed_RxP0].TotalEnergy[0], Lat->Energy[Perturbed_RxP1].TotalEnergy[0],Lat->Energy[Perturbed_RxP2].TotalEnergy[0],
          Lat->Energy[Occupied].TotalEnergy[0]);
        for (it=0; it < I->Max_Types; it++) {  // over all types
          //fprintf(convergence,"%e\t%e\t%e\t%e\t%e\t%e\t", 
          fprintf(convergence,"%e\t%e\t%e\t", 
          //I->I[it].sigma[0][0+0*NDIM], I->I[it].sigma[0][1+1*NDIM], I->I[it].sigma[0][2+2*NDIM], 
          I->I[it].sigma_rezi[0][0+0*NDIM], I->I[it].sigma_rezi[0][1+1*NDIM], I->I[it].sigma_rezi[0][2+2*NDIM]);
        }
        fprintf(convergence,"%e\t%i\t%e\t",Lat->Volume, Lev0->Plan0.plan->N[0]*Lev0->Plan0.plan->N[1]*Lev0->Plan0.plan->N[2], 
        P->Lat.SawtoothStart);
        fprintf(convergence,"%e\t%e\t%e\t",Lat->E[Perturbed_RxP0].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_RxP1].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_RxP2].AllTotalPsiEnergy[Perturbed1_0Energy]);
        fprintf(convergence,"%e\t%e\t%e\n",Lat->E[Perturbed_P0].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_P1].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_P2].AllTotalPsiEnergy[Perturbed1_0Energy]);
        fclose(convergence);
      }
    } else {
      if ((convergence = fopen("pcp.convergence.csv","r")) != NULL) { // only prepend with header if file doesn't exist yet
        fclose(convergence);
        convergence = fopen("pcp.convergence.csv","a");
      } else {
        convergence = fopen("pcp.convergence.csv","w");
        if (convergence != NULL)
          fprintf(convergence,"Ecut\ta.u.\tMaxG\tN0\tN1\tN2\tTE\tVolume\tN^3\tSawStart\n");
      }
      if (convergence != NULL) {
        fprintf(convergence,"%e\t%e\t%i\t%i\t%i\t%i\t%e\t%e\t%i\t%e\n",     
          Lat->ECut,
          Lat->RealBasisQ[0], LevS->TotalAllMaxG, Lev0->Plan0.plan->N[0], Lev0->Plan0.plan->N[1], Lev0->Plan0.plan->N[2],      
          Lat->Energy[Occupied].TotalEnergy[0],
          Lat->Volume, Lev0->Plan0.plan->N[0]*Lev0->Plan0.plan->N[1]*Lev0->Plan0.plan->N[2], 
          P->Lat.SawtoothStart);
        fclose(convergence);
      }
    }
  }
}