Index: pcp/src/energy.c
===================================================================
--- pcp/src/energy.c	(revision a0bcf15381c41adb281e0a5d9293f02515ba70c0)
+++ pcp/src/energy.c	(revision a0bcf15381c41adb281e0a5d9293f02515ba70c0)
@@ -0,0 +1,752 @@
+/** \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);
+      }
+    }
+  }
+}
