source: pcp/src/grad.c@ 774ae8

/** \file grad.c
 * Gradient Calculation.
 * 
 * The most important structure is Gradient with its Gradient::GradientArray.
 * It contains the various stages of the calculated conjugate direction, beginning
 * from the Psis::LocalPsiStatus. Subsequently, CalculateCGGradient(), CalculatePreConGrad(),
 * CalculateConDir() and CalculateLineSearch() are called which use earlier stages to
 * calculate the latter ones: ActualGradient, PreConDirGradient, ConDirGradient and OldConDirGradient,
 * do the one-dimensional minimisation and finally update the coefficients of the wave function,
 * whilst these ought to remain orthogonal.
 * 
 * Initialization of the Gradient structure and arrays is done in InitGradients(), deallocation
 * in RemoveGradients().
 * 
 * Small functions calculate the scalar product between two gradient directions GradSP(), find
 * the minimal angle in the line search CalculateDeltaI() or evaluate complex expressions
 * CalculateConDirHConDir(), down to basic calculation of the gradient of the wave function
 * CalculateGradientNoRT(). CalculateHamiltonian() sets up the Hamiltonian for the Kohn-Sham
 * orbitals and calculates the eigenvalues.
 * 
 * Other functions update the coefficients after the position of the ion has changed - 
 * CalcAlphaForUpdateWaveAfterIonMove() and UpdateWaveAfterIonMove().
 * 
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers
 \date 2000

  File: grad.c
  $Id: grad.c,v 1.90 2007-03-29 13:38:04 foo Exp $
*/

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_heapsort.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_blas.h>

#include "mymath.h"
#include "data.h"
#include "errors.h"
#include "energy.h"
#include "helpers.h"
#include "myfft.h"
#include "density.h"
#include "gramsch.h"
#include "perturbed.h"
#include "pseudo.h"
#include "grad.h"
#include "run.h"
#include "wannier.h"

/** Initialization of Gradient Calculation.
 * Gradient::GradientArray[GraSchGradient] is set onto "extra" Psis::LocalPsi (needed
 * for orthogonal conjugate gradient directions!),
 * for the remaining GradientTypes memory is allocated and zerod.
 * \param *P Problem at hand
 * \param *Grad Gradient structure
 */
void InitGradients(struct Problem *P, struct Gradient *Grad)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct LatticeLevel *ILevS = R->InitLevS;
  struct Psis *Psi = &Lat->Psi;
  int i;
  Grad->GradientArray[GraSchGradient] = LevS->LPsi->LocalPsi[Psi->LocalNo];
  for (i=1; i<MaxGradientArrayTypes; i++) {
    Grad->GradientArray[i] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILevS->MaxG, "InitGradients");
    SetArrayToDouble0((double*)Grad->GradientArray[i], 2*ILevS->MaxG);
  }
}

/** Deallocates memory used in Gradient::GradientArray.
 * Gradient::GradientArray[GraSchGradient] is set to null.
 * \param *P Problem at hand
 * \param *Grad Gradient structure
 */
void RemoveGradients(struct Problem *P, struct Gradient *Grad)
{
  int i;
  P=P;
  Grad->GradientArray[GraSchGradient] = NULL;
  for (i=1; i<MaxGradientArrayTypes; i++)
    Free(Grad->GradientArray[i], "RemoveGradients: Grad->GradientArray[i]");
}

/** Calculates the scalar product of two wave functions with equal symmetry.
 * \param *P Problem at hand
 * \param *Lev LatticeLevel structure
 * \param *LPsiDatA first array of complex wave functions coefficients
 * \param *LPsiDatB second array of complex wave functions coefficients
 * \return complex scalar product of the two wave functions.
 * \warning MPI_Allreduce has to be done still, remember coefficients are shared among processes!
 */
double GradSP(struct Problem *P, struct LatticeLevel *Lev, const fftw_complex *LPsiDatA, const fftw_complex *LPsiDatB) {
  double LocalSP=0.0;
  int i = 0;
  if (Lev->GArray[0].GSq == 0.0) {
    LocalSP += LPsiDatA[0].re * LPsiDatB[0].re;
    i++;
  }
  for (; i < Lev->MaxG; i++) {
    LocalSP += 2.*(LPsiDatA[i].re * LPsiDatB[i].re + LPsiDatA[i].im * LPsiDatB[i].im); // factor 2 because of the symmetry arising from gamma point!
  }
  return(LocalSP);
}

/** Calculates the scalar product of two wave functions with unequal symmetry.
 * \param *P Problem at hand
 * \param *Lev LatticeLevel structure
 * \param *LPsiDatA first array of complex wave functions coefficients
 * \param *LPsiDatB second array of complex wave functions coefficients
 * \return complex scalar product of the two wave functions.
 * \warning MPI_Allreduce has to be done still, remember coefficients are shared among processes!
 */
double GradImSP(struct Problem *P, struct LatticeLevel *Lev, fftw_complex *LPsiDatA, fftw_complex *LPsiDatB) {
  double LocalSP=0.0;
  int i = 0;
  if (Lev->GArray[0].GSq == 0.0) {
    i++;
  }
  for (; i < Lev->MaxG; i++) {
    LocalSP += 2.*(LPsiDatA[i].re * LPsiDatB[i].im - LPsiDatA[i].im * LPsiDatB[i].re); // factor 2 because of the symmetry arising from gamma point!
  }
  return(LocalSP);
}

/** Calculation of gradient of wave function(without RiemannTensor).
 * We want to calculate the variational derivative \f$\frac{\delta}{\delta \psi^*_j} E[{\psi^*_j}]\f$,
 * which we need for the minimisation of the energy functional. In the plane wave basis this evaluates to
 * \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
 *    \qquad (section 4.3)
 * \f]
 * Therefore, DensityType#HGDensity is back-transformed, exchange correlation potential added up (CalculateXCPotentialNoRT()),
 * and vector multiplied with the wave function Psi (per node R), transformation into reciprocal grid space and we are back
 * in the plane wave basis set as desired. CalculateAddNLPot() adds the non-local potential. Finally, the kinetic term
 * is added up (simple in reciprocal space). The energy eigenvalue \a *Lambda is computed by summation over all reciprocal grid
 * vectors and in the meantime also the gradient direction (negative derivative) stored in \a *grad while the old values
 * are moved to \a *oldgrad. At last, \a *Lambda is summed up among all processes.
 * \f[
 *    \varphi_i^{(m)} (G) = \frac{1}{2} |G|^2 c_{i,G} + \underbrace{V^H (G) + V^{ps,loc} (G)}_{+V^{HG} (G)} + V^{ps,nl} (G) | \psi_i (G) \rangle
 *    \qquad (\textnormal{end of section 4.5.2})
 * \f]
 * \param *P Problem at hand, contains Lattice and RunStruct
 * \param *source source wave function coefficients
 * \param PsiFactor occupation number, see Psis::PsiFactor
 * \param *Lambda eigenvalue \f$\lambda_i = \langle \psi_i |h|\psi_i \rangle\f$, see OnePsiElementAddData::Lambda
 * \param ***fnl non-local form factors over ion type, ion of type and reciprocal grid vector, see PseudoPot::fnl
 * \param *grad return array for (negative) gradient coefficients \f$\varphi_i^{(m)} (G)\f$
 * \param *oldgrad return array for (negative) gradient coefficients \f$\varphi_i^{(m-1)} (G)\f$
 * \param *Hc_grad return array for complex coefficients of \f$H | \psi_i (G) \rangle\f$, see GradientArray#HcGradient
 * \sa CalculateConDirHConDir() - similar calculation for another explanation.
 */
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)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  struct Density *Dens0 = Lev0->Dens;
  struct  fft_plan_3d *plan = Lat->plan; 
  int Index, g;
  fftw_complex *work = Dens0->DensityCArray[TempDensity];
  fftw_real *HGcR = Dens0->DensityArray[HGcDensity];
  fftw_complex *HGcRC = (fftw_complex*)HGcR;
  fftw_complex *HGC = Dens0->DensityCArray[HGDensity];
  fftw_real *HGCR = (fftw_real *)HGC;
  fftw_complex *PsiC = Dens0->DensityCArray[ActualPsiDensity];
  fftw_real *PsiCR = (fftw_real *)PsiC;
  int nx,ny,nz,iS,i0,s;
  const int Nx = LevS->Plan0.plan->local_nx;
  const int Ny = LevS->Plan0.plan->N[1];
  const int Nz = LevS->Plan0.plan->N[2];
  const int NUpx = LevS->NUp[0];
  const int NUpy = LevS->NUp[1];
  const int NUpz = LevS->NUp[2];
  double lambda;
  const double HGcRCFactor = 1./LevS->MaxN;
  SpeedMeasure(P, LocTime, StartTimeDo);
  // back-transform HGDensity
  fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF1, HGC, work);
  // evaluate exchange potential with this density, add up onto HGCR
  //if (isnan(HGCR[0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGcR[%i]_%i = NaN!\n", 0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
  CalculateXCPotentialNoRT(P, HGCR); // add V^{xc} or \epsilon^{xc} upon V^H + V^{ps}
  for (nx=0;nx<Nx;nx++)  
    for (ny=0;ny<Ny;ny++) 
      for (nz=0;nz<Nz;nz++) {
        i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
        iS = nz+Nz*(ny+Ny*nx);
        HGcR[iS] = HGCR[i0]*PsiCR[i0]; /* Matrix Vector Mult */
        //if (isnan(HGCR[i0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGcR[%i]_%i = NaN!\n", i0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
      }
  fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGcRC, work);
  SpeedMeasure(P, LocTime, StopTimeDo);
  /* NonLocalPP */
  //SpeedMeasure(P, NonLocTime, StartTimeDo);
  CalculateAddNLPot(P, Hc_grad, fnl, PsiFactor);  // also resets Hc_grad beforehand
  //SetArrayToDouble0((double *)Hc_grad,2*LevS->MaxG);
  SpeedMeasure(P, NonLocTime, StopTimeDo);
  /* Lambda und Grad */
  for (g=0;g<LevS->MaxG;g++) {
    Index = LevS->GArray[g].Index; /* FIXME - factoren */
    Hc_grad[g].re += PsiFactor*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[g].GSq*source[g].re);
    Hc_grad[g].im += PsiFactor*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[g].GSq*source[g].im);
    //if (isnan(source[g].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): source_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); }  
  }
  if (grad != NULL) {
    lambda = 0.0;
    s = 0;
    if (LevS->GArray[0].GSq == 0.0) {
      lambda += Hc_grad[0].re*source[0].re + Hc_grad[0].im*source[0].im;
      oldgrad[0].re = grad[0].re; // store away old preconditioned steepest decent gradient
      oldgrad[0].im = grad[0].im;
      grad[0].re = -Hc_grad[0].re;
      grad[0].im = -Hc_grad[0].im;
      s++;
    }
    for (g=s;g<LevS->MaxG;g++) {
      lambda += 2.*(Hc_grad[g].re*source[g].re + Hc_grad[g].im*source[g].im);
      oldgrad[g].re = grad[g].re; // store away old preconditioned steepest decent gradient
      oldgrad[g].im = grad[g].im;
      grad[g].re = -Hc_grad[g].re;
      grad[g].im = -Hc_grad[g].im;
      //if (isnan(Hc_grad[g].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): Hc_grad_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); }
    }
    MPI_Allreduce ( &lambda, Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
    //fprintf(stderr,"(%i) Lambda = %e\n",P->Par.me, *Lambda);
  }
}


/** Calculates the Hamiltonian matrix in Kohn-Sham basis set and outputs to file.
 * In a similiar manner to TestGramSch() all necessary coefficients are retrieved from the
 * other processes and \f$\langle Psi_i |  H | Psi_j \rangle\f$ = \f$\lambda_{i,j}\f$ calculated.
 * \param *P Problem at hand
 * \sa CalculateGradientNoRT() - same calculation (of \f$\lambda_i\f$), yet only for diagonal elements
 */
void CalculateHamiltonian(struct Problem *P) {
  struct Lattice *Lat = &P->Lat;
  struct FileData *F = &P->Files;
  struct RunStruct *R = &P->R;
  struct Psis *Psi = &Lat->Psi;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  //struct LevelPsi *LPsi = LevS->LPsi;
  struct Density *Dens0 = Lev0->Dens;
  int Index;
  struct  fft_plan_3d *plan = Lat->plan; 
  fftw_complex *work = Dens0->DensityCArray[TempDensity];
  fftw_real *HGcR = Dens0->DensityArray[HGcDensity];
  fftw_complex *HGcRC = (fftw_complex*)HGcR;
  fftw_complex *HGC = Dens0->DensityCArray[HGDensity];
  fftw_real *HGCR = (fftw_real *)HGC;
  fftw_complex *PsiC = Dens0->DensityCArray[ActualPsiDensity];
  fftw_real *PsiCR = (fftw_real *)PsiC;
  fftw_complex ***fnl;
  fftw_complex *Hc_grad = P->Grad.GradientArray[HcGradient];
  int nx,ny,nz,iS,i0,s;
  const int Nx = LevS->Plan0.plan->local_nx;
  const int Ny = LevS->Plan0.plan->N[1];
  const int Nz = LevS->Plan0.plan->N[2];
  const int NUpx = LevS->NUp[0];
  const int NUpy = LevS->NUp[1];
  const int NUpz = LevS->NUp[2];
  double lambda, Lambda;
  const double HGcRCFactor = 1./LevS->MaxN;
  double pot_xc, Pot_xc;
  int PsiFactorA = 1;
  int PsiFactorB = 1;
  int NoPsis = 0;
  int NoOfPsis = 0;

  if (P->Call.out[LeaderOut])
    fprintf(stderr, "(%i)Setting up Hamiltonian for Level %d...",P->Par.me, LevS->LevelNo);
  if (P->Call.out[StepLeaderOut])
    fprintf(stderr, "\n(%i) H = \n",P->Par.me);

  ApplyTotalHamiltonian(P,LevS->LPsi->LocalPsi[0],Hc_grad, P->PP.fnl[0], 1., 1); // update HGDensity->
  
  // go through all pairs of orbitals
  int i,k,l,m,j=0,RecvSource;
  MPI_Status status;
  struct OnePsiElement *OnePsiB, *OnePsiA, *LOnePsiA, *LOnePsiB;
  int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  fftw_complex *LPsiDatA=NULL, *LPsiDatB;
  
  l = -1;
  switch (Psi->PsiST) {
    case SpinDouble:
      NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoDouble];
      NoOfPsis += Psi->GlobalNo[PsiMaxAdd];      
      break;
    case SpinUp:
      NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoUp];
      NoOfPsis += Psi->GlobalNo[PsiMaxAdd];      
      break;
    case SpinDown:
      NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoDown];
      NoOfPsis += Psi->GlobalNo[PsiMaxAdd];      
      break;
  };
  gsl_matrix *H = gsl_matrix_calloc(NoOfPsis,NoOfPsis);
  for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) {  // go through all wave functions
    //fprintf(stderr,"(%i) GlobalNo: %d\tLocalNo: %d\n", P->Par.me,Psi->AllPsiStatus[i].MyGlobalNo,Psi->AllPsiStatus[i].MyLocalNo);
    OnePsiA = &Psi->AllPsiStatus[i];    // grab OnePsiA
    if (OnePsiA->PsiType <= UnOccupied) {   // drop the extra and perturbed ones
      l++;
      if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
         LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
      else 
         LOnePsiA = NULL;
      if (LOnePsiA == NULL) {   // if it's not local ... receive it from respective process into TempPsi
        RecvSource = OnePsiA->my_color_comm_ST_Psi;
        MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, HamiltonianTag, P->Par.comm_ST_PsiT, &status );
        LPsiDatA=LevS->LPsi->TempPsi;
      } else {                  // .. otherwise send it to all other processes (Max_me... - 1)
        for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++)
          if (k != OnePsiA->my_color_comm_ST_Psi) 
            MPI_Send( LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, k, HamiltonianTag, P->Par.comm_ST_PsiT);
        LPsiDatA=LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo];
      } // LPsiDatA is now set to the coefficients of OnePsi either stored or MPI_Received

      if (LOnePsiA != NULL) { // if fnl (!) are locally available!
        // Trick is: fft wave function by moving it to higher level and there take only each ...th element into account
        PsiFactorA = 1; //Psi->LocalPsiStatus[OnePsiA->MyLocalNo].PsiFactor;
        CalculateOneDensityR(Lat, LevS, Dens0, LPsiDatA, Dens0->DensityArray[ActualDensity], R->FactorDensityR*PsiFactorA, 1);
        // note: factor is not used when storing result in DensityCArray[ActualPsiDensity] in CalculateOneDensityR()!
        for (nx=0;nx<Nx;nx++)  
          for (ny=0;ny<Ny;ny++) 
            for (nz=0;nz<Nz;nz++) {
              i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
              iS = nz+Nz*(ny+Ny*nx);
              HGcR[iS] = HGCR[i0]*PsiCR[i0]; /* Matrix Vector Mult */
            }
        // V^{HG} (R) + V^{XC} (R) | \Psi_b (R) \rangle
        fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGcRC, work);
        // V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle          
        /* NonLocalPP */
        fnl = P->PP.fnl[OnePsiA->MyLocalNo];    // fnl only locally available!
        CalculateAddNLPot(P, Hc_grad, fnl, PsiFactorA); // sets Hc_grad to zero beforehand
        // V^{ps,nl} (G)| \Psi_b (G) + [V^{HG} (G) + V^{XC} (G)] | \Psi_b (G) \rangle          
        /* Lambda und Grad */
        s = 0;
        if (LevS->GArray[0].GSq == 0.0) {
          Index = LevS->GArray[0].Index; /* FIXME - factoren */
          Hc_grad[0].re += PsiFactorA*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[0].GSq*LPsiDatA[0].re);
          Hc_grad[0].im += PsiFactorA*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[0].GSq*LPsiDatA[0].im);
          s++;
        }
        for (;s<LevS->MaxG;s++) {
          Index = LevS->GArray[s].Index; /* FIXME - factoren */
          Hc_grad[s].re += PsiFactorA*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[s].GSq*LPsiDatA[s].re);
          Hc_grad[s].im += PsiFactorA*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[s].GSq*LPsiDatA[s].im);
        }
        // 1/2 |G|^2 + V^{ps,nl} (G) + V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle          
      } else {
        SetArrayToDouble0((double *)Hc_grad,2*LevS->MaxG);     
      }

      m = -1;
      // former border of following loop: Psi->GlobalNo[Psi->PsiST+1]+Psi->GlobalNo[Psi->PsiST+4]+P->Par.Max_me_comm_ST_PsiT
      for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) {  // go through all wave functions
        OnePsiB = &Psi->AllPsiStatus[j];    // grab OnePsiA
        if (OnePsiB->PsiType <= UnOccupied) {   // drop the extra ones
          m++;
          if (OnePsiB->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
             LOnePsiB = &Psi->LocalPsiStatus[OnePsiB->MyLocalNo];
          else 
             LOnePsiB = NULL;
          if (LOnePsiB == NULL) {   // if it's not local ... receive it from respective process into TempPsi
            RecvSource = OnePsiB->my_color_comm_ST_Psi;
            MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, HamiltonianTag, P->Par.comm_ST_PsiT, &status );
            LPsiDatB=LevS->LPsi->TempPsi;
          } else {                  // .. otherwise send it to all other processes (Max_me... - 1)
            for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++)
              if (k != OnePsiB->my_color_comm_ST_Psi) 
                MPI_Send( LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, k, HamiltonianTag, P->Par.comm_ST_PsiT);
            LPsiDatB=LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo];
          } // LPsiDatB is now set to the coefficients of OnePsi either stored or MPI_Received

          PsiFactorB = 1; //Psi->LocalPsiStatus[OnePsiB->MyLocalNo].PsiFactor;
          lambda = 0; //PsiFactorB * GradSP(P, LevS, LPsiDatB, Hc_grad); 
          Lambda = 0; 
          s=0;
          if (LevS->GArray[0].GSq == 0.0) {
            lambda += GSL_REAL(gsl_complex_mul(gsl_complex_conjugate(convertComplex(LPsiDatB[s])), convertComplex(Hc_grad[s])));
            s++;
          }
          for (; s < LevS->MaxG; s++) {
            lambda += GSL_REAL(gsl_complex_mul(gsl_complex_conjugate(convertComplex(LPsiDatB[s])), convertComplex(Hc_grad[s]))); 
            lambda += GSL_REAL(gsl_complex_mul(convertComplex(LPsiDatB[s]), gsl_complex_conjugate(convertComplex(Hc_grad[s]))));  // c(-G) = c^\ast (G) due to gamma point symmetry !
          } // due to the symmetry the resulting lambda is both real and symmetric in l,m (i,j) !
          lambda *= PsiFactorB;
          // \langle \Psi_a | 1/2 |G|^2 + V^{ps,nl} (G) + V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle
          MPI_Allreduce ( &lambda, &Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST);  // gather from all(!) (coeffs- and Psis-sharing) groups

          //if (P->Par.me == 0) {
          // and Output
            if (P->Call.out[StepLeaderOut])
              fprintf(stderr,"%2.5lg\t",Lambda);
            if (P->Par.me == 0) fprintf(F->HamiltonianFile, "%lg\t", Lambda);
          //}
          if (m < NoOfPsis && l < NoOfPsis) {
            //fprintf(stderr,"Index (%i/%i,%i/%i): %e \n",l,NoOfPsis,m,NoOfPsis, Lambda);
            gsl_matrix_set(H,l,m,Lambda);
          } else
            fprintf(stderr,"Index (%i,%i) out of range %i\n",l,m,NoOfPsis);
        }
        if ((P->Par.me == 0) && (j == NoOfPsis-1))
          fprintf(F->HamiltonianFile, "\n");        
        if ((P->Call.out[StepLeaderOut]) && (j == NoOfPsis-1))
            fprintf(stderr,"\n");       
      }
    }
  }
  if (P->Par.me == 0) {
    fprintf(F->HamiltonianFile, "\n");
    fflush(F->HamiltonianFile);
  }
  if ((P->Call.out[LeaderOut]) && (!P->Call.out[StepLeaderOut])) {
    fprintf(stderr,"finished.\n");
  }
  
  // storing H matrix for latter use in perturbed calculation
  for (i=0;i<Psi->NoOfPsis;i++)
    for (j=0;j<Psi->NoOfPsis;j++)
      Psi->lambda[i][j] = gsl_matrix_get(H,i,j);

  // retrieve EW from H
  gsl_vector *eval = gsl_vector_alloc(NoOfPsis);
  gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc(NoOfPsis);
  gsl_eigen_symm(H, eval, w);
  gsl_eigen_symm_free(w);
  gsl_matrix_free(H);
  
  // print eigenvalues
  if(P->Call.out[LeaderOut]) {
    fprintf(stderr,"(%i) Unsorted Eigenvalues:", P->Par.me);
    for (i=0;i<NoOfPsis;i++)
      fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
    fprintf(stderr,"\n");
  }
  gsl_sort_vector(eval);  // sort eigenvalues
  // print eigenvalues
  if(P->Call.out[LeaderOut]) {
    fprintf(stderr,"(%i) All Eigenvalues:", P->Par.me);
    for (i=0;i<NoOfPsis;i++)
      fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
    fprintf(stderr,"\n");
    fprintf(stderr,"(%i) KS Eigenvalues:", P->Par.me);
    for (i=0;i<NoOfPsis;i++)
      fprintf(stderr,"\t%lg",Psi->AddData[i].Lambda);
    fprintf(stderr,"\n");
  }
  
  // print difference between highest occupied and lowest unoccupied
  if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0) {
    Lat->Energy[UnOccupied].bandgap = gsl_vector_get(eval, Psi->NoOfPsis) - gsl_vector_get(eval, Psi->NoOfPsis-1);
    if (P->Call.out[NormalOut])
      fprintf(stderr,"(%i) Band gap: \\epsilon_{%d+%d}(%d+1) - \\epsilon_{%d+%d}(%d) = %lg Ht\n", P->Par.me, NoPsis,P->Lat.Psi.GlobalNo[PsiMaxAdd],NoPsis, NoPsis, P->Lat.Psi.GlobalNo[PsiMaxAdd], NoPsis, Lat->Energy[UnOccupied].bandgap);
  }
  // now, calculate \int v_xc (r) n(r) d^3r
  SetArrayToDouble0((double *)HGCR,2*Dens0->TotalSize);
  CalculateXCPotentialNoRT(P, HGCR); /* FIXME - kann man halbieren */
  pot_xc = 0;
  Pot_xc = 0;
  for (i=0;i<Dens0->LocalSizeR;i++)
    pot_xc += HGCR[i]*Dens0->DensityArray[TotalDensity][i];
  MPI_Allreduce ( &pot_xc, &Pot_xc, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  if (P->Call.out[StepLeaderOut]) {
    fprintf(stderr,"(%i) Exchange-Correlation potential energy: \\int v_xc[n] (r) n(r) d^3r = %lg\n", P->Par.me, HGcRCFactor*Pot_xc*Lat->Volume);
  }
  gsl_vector_free(eval);
}


/** Calculation of direction of steepest descent.
 * The Gradient has been evaluated in CalculateGradient(). Now the SD-direction
 * \f$\zeta_i = - (H-\lambda_i)\psi_i\f$ is stored in 
 * Gradient::GradientArray[GradientTypes#GraSchGradient].
 * Afterwards, the "extra" local Psi (which is the above gradient) is orthogonalized
 * under exclusion of the to be minimialized Psi to ascertain Orthonormality of this
 * gradient with respect to all other states/orbits. 
 * \param *P Problem at hand, contains Lattice and RunStruct
 * \param *source source wave function coefficients, \f$\psi_i\f$
 * \param *grad gradient coefficients, see CalculateGradient(), on return contains conjugate gradient coefficients, \f$\varphi_i^{(m)}\f$
 * \param *Lambda eigenvalue \f$\lambda_i\f$ = \f$\langle \psi_i^{(m)}|H|\psi_i^{(m)} \rangle\f$
 * \warning coefficients in \a *grad are overwritten
 */
static void CalculateCGGradient(struct Problem *P, fftw_complex *source, fftw_complex *grad, double Lambda)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct Psis *Psi = &Lat->Psi; 
  struct LatticeLevel *LevS = R->LevS;
  fftw_complex *GradOrtho = P->Grad.GradientArray[GraSchGradient];
  int g;
  int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  for (g=0;g<LevS->MaxG;g++) {
    GradOrtho[g].re = grad[g].re+Lambda*source[g].re;
    GradOrtho[g].im = grad[g].im+Lambda*source[g].im;
  }
  // include the ExtraPsi (which is the GraSchGradient!)
  SetGramSchExtraPsi(P, Psi, NotOrthogonal);
  // exclude the minimised Psi
  SetGramSchActualPsi(P, Psi, NotUsedToOrtho);
  SpeedMeasure(P, GramSchTime, StartTimeDo);
  // makes conjugate gradient orthogonal to all other orbits
  //fprintf(stderr,"CalculateCGGradient: GramSch() for extra orbital\n");
  GramSch(P, LevS, Psi, Orthogonalize);
  SpeedMeasure(P, GramSchTime, StopTimeDo);
  SetGramSchActualPsi(P, Psi, IsOrthonormal);
  memcpy(grad, GradOrtho, ElementSize*LevS->MaxG*sizeof(double));
}

/** Preconditioning for the gradient calculation.
 * A suitable matrix for preconditioning is [TPA89]
 * \f[
 *    M_{G,G'} = \delta_{G,G'} \frac{27+18x+12x^2+8x^3}{27+18x+12x^2+8x^3+16x^4} \qquad (5.3)
 * \f]
 *  where \f$x= \frac{\langle \xi_G|-\frac{1}{2}\nabla^2|\xi_G \rangle}
 *    {\langle \psi_i^{(m)}|-\frac{1}{2}\nabla^2|\psi_i^{(m)} \rangle}\f$.
 * (This simple form is due to the diagonal dominance of the kinetic term in Hamiltonian, and
 * x is much simpler to evaluate in reciprocal space!)
 * With this one the higher plane wave states are being transformed such that they
 * nearly equal the error vector and converge at an equal rate.
 * So that the new steepest descent direction is \f$M \tilde{\zeta}_i\f$.
 * Afterwards, the new direction is orthonormalized via GramSch() and the "extra"
 * local Psi in a likewise manner to CalculateCGGradient().
 * 
 * \param *P Problem at hand, contains LatticeLevel and RunStruct
 * \param *grad array of steepest descent coefficients, \f$\varphi_i^{(m)}\f$
 * \param *PCgrad return array for preconditioned steepest descent coefficients, \f$\tilde{\varphi}_i^{(m)}\f$
 * \param T kinetic eigenvalue \f$\langle \psi_i | - \frac{1}{2} \nabla^2 | \psi_i \rangle\f$
 */
static void CalculatePreConGrad(struct Problem *P, fftw_complex *grad, fftw_complex *PCgrad, double T)
{
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  int g;
  double x, K, dK;
  struct Psis *Psi = &P->Lat.Psi; 
  fftw_complex *PCOrtho = P->Grad.GradientArray[GraSchGradient];
  int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  if (fabs(T) < MYEPSILON) T = 1;
  for (g=0;g<LevS->MaxG;g++) {
    x = LevS->GArray[g].GSq;
    x /= T;
    K = (((8*x+12)*x +18)*x+27);
    dK = ((((16*x+8)*x +12)*x+18)*x+27);
    K /= dK;
    c_re(PCOrtho[g]) = K*c_re(grad[g]);
    c_im(PCOrtho[g]) = K*c_im(grad[g]);  
  }
  SetGramSchExtraPsi(P, Psi, NotOrthogonal);
  SpeedMeasure(P, GramSchTime, StartTimeDo);
  // preconditioned direction is orthogonalized
  //fprintf(stderr,"CalculatePreConGrad: GramSch() for extra orbital\n");
  GramSch(P, LevS, Psi, Orthogonalize);
  SpeedMeasure(P, GramSchTime, StopTimeDo);
  memcpy(PCgrad, PCOrtho, ElementSize*LevS->MaxG*sizeof(double));
}

/** Calculation of conjugate direction.
 * The new conjugate direction is according to the algorithm of Fletcher-Reeves:
 * \f[
 *  \varphi^{(m)}_i = 
 *    \tilde{\eta}^{(m)}_i + \frac{\langle \tilde{\eta}^{(m)}_i | \tilde{\varphi}_i^{(m)} \rangle}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} \varphi_i^{(m-1)}  
 * \f]
 * Instead we use the algorithm of Polak-Ribiere:
 * \f[
 *  \varphi^{(m)}_i = 
 *    \tilde{\eta}^{(m)}_i + \frac{\langle \tilde{\eta}^{(m)}_i | \bigl (\tilde{\varphi}_i^{(m)} \rangle - \tilde{\varphi}_i^{(m-1)} \rangle \bigr)}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} \varphi_i^{(m-1)}  
 * \f]
 * with an additional automatic restart:
 * \f[
 *  \varphi^{(m)}_i = \tilde{\eta}^{(m)}_i \quad \textnormal{if}\ \frac{\langle \tilde{\eta}^{(m)}_i | \bigl (\tilde{\varphi}_i^{(m)} \rangle - \tilde{\varphi}_i^{(m-1)} \rangle \bigr)}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} <= 0
 * \f]
 * That's why ConDirGradient and OldConDirGradient are stored, the scalar product - summed up among all processes - 
 * is saved in OnePsiAddData::Gamma. The conjugate direction is again orthonormalized by GramSch() and stored
 * in \a *ConDir.
 * \param *P Problem at hand, containing LatticeLevel and RunStruct
 * \param step minimisation step m
 * \param *GammaDiv scalar product \f$\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle\f$, on calling expects (m-1), on return contains (m)
 * \param *grad array for coefficients of steepest descent direction \f$\tilde{\zeta}_i^{(m)}\f$
 * \param *oldgrad array for coefficients of steepest descent direction of \a step-1 \f$\tilde{\zeta}_i^{(m-1)}\f$
 * \param *PCgrad array for coefficients of preconditioned steepest descent direction \f$\tilde{\eta}_i^{(m)}\f$
 * \param *ConDir return array for coefficients of conjugate array of \a step, \f$\tilde{\varphi}_i^{(m)}\f$ = \f$\frac{\tilde{\varphi}_i^{(m)}} {||\tilde{\varphi}_i^{(m)}||}\f$
 * \param *ConDir_old array for coefficients of conjugate array of \a step-1, \f$\tilde{\varphi}_i^{(m-1)}\f$
 */
static void CalculateConDir(struct Problem *P, int step, double *GammaDiv,
                            fftw_complex *grad,  
          fftw_complex *oldgrad,  
                            fftw_complex *PCgrad,
                            fftw_complex *ConDir, fftw_complex *ConDir_old) {
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct Psis *Psi = &P->Lat.Psi; 
  fftw_complex *Ortho = P->Grad.GradientArray[GraSchGradient];
  int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  int g;
  double S[3], Gamma, GammaDivOld = *GammaDiv;
  double LocalSP[2], PsiSP[2];
  LocalSP[0] = GradSP(P, LevS, PCgrad, grad);
  LocalSP[1] = GradSP(P, LevS, PCgrad, oldgrad);

  MPI_Allreduce ( LocalSP, PsiSP, 2, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  *GammaDiv = S[0] = PsiSP[0];
  S[1] = GammaDivOld;
  S[2] = PsiSP[1]; 
  if (step) { // only in later steps is the scalar product used, but always condir stored in oldcondir and Ortho (working gradient)
    if (fabs(S[1]) < MYEPSILON) fprintf(stderr,"CalculateConDir: S[1] = %lg\n",S[1]);
    Gamma = (S[0]-S[2])/S[1];  // the final CG beta coefficient (Polak-Ribiere)
    if (fabs(S[1]) < MYEPSILON) {
      if (fabs((S[0]-S[2])) < MYEPSILON)  
        Gamma = 1.0;
      else
        Gamma = 0.0;
    }
    if (Gamma < 0) { // Polak-Ribiere with automatic restart: gamma = max {0, gamma}
      Gamma = 0.;
      if (P->Call.out[StepLeaderOut]) fprintf(stderr,"(%i) Polak-Ribiere: Restarting iteration by setting gamma = 0.\n", P->Par.me);
    }
    for (g=0; g < LevS->MaxG; g++) {
      c_re(ConDir[g]) = c_re(PCgrad[g]) + Gamma*c_re(ConDir_old[g]);
      c_im(ConDir[g]) = c_im(PCgrad[g]) + Gamma*c_im(ConDir_old[g]);
      c_re(ConDir_old[g]) = c_re(ConDir[g]);
      c_im(ConDir_old[g]) = c_im(ConDir[g]);
      c_re(Ortho[g]) = c_re(ConDir[g]);
      c_im(Ortho[g]) = c_im(ConDir[g]);
    }
  } else {
    Gamma = 0.0;
    for (g=0; g < LevS->MaxG; g++) {
      c_re(ConDir[g]) = c_re(PCgrad[g]);
      c_im(ConDir[g]) = c_im(PCgrad[g]);
      c_re(ConDir_old[g]) = c_re(ConDir[g]);
      c_im(ConDir_old[g]) = c_im(ConDir[g]);
      c_re(Ortho[g]) = c_re(ConDir[g]);
      c_im(Ortho[g]) = c_im(ConDir[g]);
    }
  }
  // orthonormalize
  SetGramSchExtraPsi(P, Psi, NotOrthogonal);
  SpeedMeasure(P, GramSchTime, StartTimeDo);
  //fprintf(stderr,"CalculateConDir: GramSch() for extra orbital\n");
  GramSch(P, LevS, Psi, Orthonormalize);
  SpeedMeasure(P, GramSchTime, StopTimeDo);
  memcpy(ConDir, Ortho, ElementSize*LevS->MaxG*sizeof(double));
}

/** Calculation of \f$\langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\widetilde{\widetilde{\varphi}}_i^{(m)} \rangle\f$.
 * Evaluates most of the parts of the second derivative of the energy functional \f$\frac{\delta^2 E}{\delta \Theta^2}|_{\Theta=0}\f$,
 * most importantly the energy expection value of the conjugate direction.
 * \f[
 *    \langle \widetilde{\widetilde{\varphi}}_i^{(m)} (G) | -\frac{1}{2} \nabla^2 + \underbrace{V^H (G) + V^{ps,loc} (G)}_{=V^{HG}(G)} + V^{ps,nl} | \widetilde{\widetilde{\varphi}} (G) \rangle
 * \f]
 * The conjugate direction coefficients \f$\varphi_i^{(m)} (G)\f$ are fftransformed to \f$\varphi_i^{(m)} (R)\f$,
 * times the DoubleDensityTypes#HGDensity \f$V^{HG} (R) \f$we have the first part,
 * then \f$\rho (R) = 2\frac{f_i}{V} \Re (\varphi_i^{(m)} (R) \cdot \psi_i (R))\f$is evaluated,
 * HGDensity is fftransformed to \f$V^{HG} | \varphi_i^{(m)} (G) \rangle\f$,
 * CalculateAddNLPot() adds non-local potential \f$V^{ps,nl} (G)\f$ to it,
 * kinetic term is added in the reciprocal base as well \f$-\frac{1}{2} |G|^2\f$,
 * and at last scalar product between it and the conjugate direction, and done.
 *
 * Meanwhile the fftransformed ConDir density \f$\rho (R)\f$ is used for the ddEddt0s \f$A_H\f$, \f$A_{XC}\f$.
 * In case of RunStruct::DoCalcCGGauss in the following expression additionally \f$\hat{n}^{Gauss}\f$
 * is evaluted, otherwise left out:
 * \f[
 *    A^H = \frac{4\pi}{V} \sum_G \frac{|\frac{\rho(G)}{V^2 N} + \hat{n}^{Gauss}|}{|G|^2} 
 *    \qquad (\textnormal{section 5.1, line search})
 * \f]
 * \param *P Problem at hand
 * \param *ConDir conjugate direction \f$\tilde{\tilde{\varphi}}_i^{(m)}\f$
 * \param PsiFactor occupation number, see Psis::PsiFactor
 * \param *ConDirHConDir first return value: energy expectation value for conjugate direction
 * \param *HartreeddEddt0 second return value: second derivative of Hartree energy \f$A_H\f$
 * \param *XCddEddt0 third return value: second derivative of exchange correlation energy \f$A_{XC}\f$
 * \warning The MPI_Allreduce for the scalar product in the end has not been done!
 * \sa CalculateLineSearch() - used there
 */
void CalculateConDirHConDir(struct Problem *P, fftw_complex *ConDir, const double PsiFactor,
                                   double *ConDirHConDir, double *HartreeddEddt0, double *XCddEddt0)
{
  struct Lattice *Lat = &P->Lat;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  struct Density *Dens0 = Lev0->Dens;
  int Index, g, i, pos, s=0;
  struct  fft_plan_3d *plan = Lat->plan; 
  fftw_complex *work = Dens0->DensityCArray[TempDensity];
  fftw_real *PsiCD = (fftw_real *)work;
  fftw_complex *PsiCDC = work;
  fftw_real *PsiR = (fftw_real *)Dens0->DensityCArray[ActualPsiDensity];
  fftw_complex *tempdestRC = (fftw_complex *)Dens0->DensityArray[TempDensity];
  fftw_real *HGCDR = Dens0->DensityArray[HGcDensity];
  fftw_complex *HGCDRC = (fftw_complex*)HGCDR;
  fftw_complex *HGC = Dens0->DensityCArray[HGDensity];
  fftw_real *HGCR = (fftw_real *)HGC;
  fftw_complex *posfac, *destpos, *destRCS, *destRCD;
  const double HartreeFactorR = 2.*R->FactorDensityR*PsiFactor;
  double HartreeFactorC = 4.*PI*Lat->Volume;
  double HartreeFactorCGauss = 1./Lev0->MaxN; //R->FactorDensityC*R->HGcFactor/Lev0->MaxN; // same expression as first and second factor cancel each other
  int nx,ny,nz,iS,i0;
  const int Nx = LevS->Plan0.plan->local_nx;
  const int Ny = LevS->Plan0.plan->N[1];
  const int Nz = LevS->Plan0.plan->N[2];
  const int NUpx = LevS->NUp[0];
  const int NUpy = LevS->NUp[1];
  const int NUpz = LevS->NUp[2];
  const double HGcRCFactor = 1./LevS->MaxN;
  fftw_complex *HConDir = P->Grad.GradientArray[TempGradient];
  int DoCalcCGGauss = R->DoCalcCGGauss;
  fftw_complex rp, rhog;
  int it;
  struct Ions *I = &P->Ion;

  /* ConDir H ConDir */
  // (G)->(R): retrieve density from wave functions into PsiCD in the usual manner (yet within the same level!) (see CalculateOneDensityR())
  SetArrayToDouble0((double *)tempdestRC, Dens0->TotalSize*2);
  for (i=0;i<LevS->MaxG;i++) {
    Index = LevS->GArray[i].Index;
    posfac = &LevS->PosFactorUp[LevS->MaxNUp*i];
    destpos = &tempdestRC[LevS->MaxNUp*Index];
    for (pos=0; pos < LevS->MaxNUp; pos++) {
      destpos[pos].re = ConDir[i].re*posfac[pos].re-ConDir[i].im*posfac[pos].im;
      destpos[pos].im = ConDir[i].re*posfac[pos].im+ConDir[i].im*posfac[pos].re;
    }
  }
  for (i=0; i<LevS->MaxDoubleG; i++) {
    destRCS = &tempdestRC[LevS->DoubleG[2*i]*LevS->MaxNUp];
    destRCD = &tempdestRC[LevS->DoubleG[2*i+1]*LevS->MaxNUp];
    for (pos=0; pos < LevS->MaxNUp; pos++) {
      destRCD[pos].re =  destRCS[pos].re;
      destRCD[pos].im = -destRCS[pos].im;
    }
  }
  fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNFUp, tempdestRC, work);
  DensityRTransformPos(LevS,(fftw_real*)tempdestRC, PsiCD);
  // HGCR contains real(!) HGDensity from CalculateGradientNoRT() !
  for (nx=0;nx<Nx;nx++)  
    for (ny=0;ny<Ny;ny++) 
      for (nz=0;nz<Nz;nz++) {
        i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
        iS = nz+Nz*(ny+Ny*nx);
        HGCDR[iS] = HGCR[i0]*PsiCD[i0]; /* Matrix Vector Mult */
      }
   // now we have V^HG + V^XC | \varphi \rangle
  /*
    Hartree & Exchange Correlaton ddt0
   */
  for (i=0; i<Dens0->LocalSizeR; i++) 
    PsiCD[i] *= HartreeFactorR*PsiR[i];   // is rho(r) = PsiFactor times \psi_i is in PsiR times 1/V times \varphi_i^{(m)}
  *XCddEddt0 = CalculateXCddEddt0NoRT(P, PsiCD);    // calcs A^{XC}
  fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF1, PsiCDC, tempdestRC);
  s = 0;
  *HartreeddEddt0 = 0.0;
  if (Lev0->GArray[0].GSq == 0.0)
    s++;
  if (DoCalcCGGauss) {  // this part is wrong, n^Gauss is not part of HartreeddEddt0 (only appears in 2*(<phi_i| H |phi_i> - ...)) 
    for (g=s; g < Lev0->MaxG; g++) {
      Index = Lev0->GArray[g].Index;
      rp.re = 0.0; rp.im = 0.0;
      for (it = 0; it < I->Max_Types; it++) {
        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(PsiCDC[Index])*HartreeFactorCGauss+c_re(rp);
      c_im(rhog) = c_im(PsiCDC[Index])*HartreeFactorCGauss+c_im(rp);
      if (Lev0->GArray[g].GSq < MYEPSILON) fprintf(stderr,"CalculateConDirHConDir: Lev0->GArray[g].GSq = %lg",Lev0->GArray[g].GSq);
      *HartreeddEddt0 += (c_re(rhog)*c_re(rhog)+c_im(rhog)*c_im(rhog))/(Lev0->GArray[g].GSq);
    }
  } else {
    HartreeFactorC *= HartreeFactorCGauss*HartreeFactorCGauss;
    for (g=s; g < Lev0->MaxG; g++) {
      Index = Lev0->GArray[g].Index;
      if (Lev0->GArray[g].GSq < MYEPSILON) fprintf(stderr,"CalculateConDirHConDir: Lev0->GArray[g].GSq = %lg",Lev0->GArray[g].GSq);
      *HartreeddEddt0 += (c_re(PsiCDC[Index])*c_re(PsiCDC[Index])+c_im(PsiCDC[Index])*c_im(PsiCDC[Index]))/(Lev0->GArray[g].GSq);
    }
  }
  *HartreeddEddt0 *= HartreeFactorC;
  /* fertig mit ddEddt0 */
  
  fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGCDRC, work);
  CalculateCDfnl(P, ConDir, PP->CDfnl); // calculate needed non-local form factors
  CalculateAddNLPot(P, HConDir, PP->CDfnl, PsiFactor); // add non-local potential, sets to zero beforehand
  //SetArrayToDouble0((double *)HConDir,2*LevS->MaxG);
  // now we have V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle
  for (g=0;g<LevS->MaxG;g++) {  // factor by one over volume and add kinetic part
    Index = LevS->GArray[g].Index; /* FIXME - factoren */
    HConDir[g].re += PsiFactor*(HGCDRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[g].GSq*ConDir[g].re);
    HConDir[g].im += PsiFactor*(HGCDRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[g].GSq*ConDir[g].im);
  }
   // now, T^{kin} + V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle
  *ConDirHConDir = GradSP(P, LevS, ConDir, HConDir);
   // finally, \langle \varphi | V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle which is returned
}

/** Find the minimal \f$\Theta_{min}\f$ for the line search.
 * We find it as
 * \f[ 
 *  \widetilde{\Theta}_{min} = \frac{1}{2} \tan^{-1} \Bigl ( -\frac{ \frac{\delta E}{\delta \Theta}_{|\Theta =0} }
 *    { \frac{1}{2} \frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0} } \Bigr )
 * \f]
 * Thus, with given \f$E(0)\f$, \f$\frac{\delta E}{\delta \Theta}_{|\Theta =0}\f$ and \f$\frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0}\f$
 * \f$\Theta_{min}\f$is calculated and the unknowns A, B and C (see CalculateLineSearch()) obtained. So that our
 * approximating formula \f$\widetilde{E}(\Theta) = A + B cos(2\Theta) + CSin(2\Theta)\f$ is fully known and we
 * can evaluate the energy functional and its derivative at \f$\Theta_{min} \f$.
 * \param E0 energy eigenvalue at \f$\Theta = 0\f$: \f$A= E(0) - B\f$
 * \param dEdt0 first derivative of energy eigenvalue at \f$\Theta = 0\f$: \f$C = \frac{1}{2} \frac{\delta E}{\delta \Theta}|_{\Theta=0}\f$
 * \param ddEddt0 first derivative of energy eigenvalue at \f$\Theta = 0\f$: \f$B = -\frac{1}{4} \frac{\delta^2 E}{\delta \Theta^2}|_{\Theta=0}\f$
 * \param *E returned energy functional at \f$\Theta_{min}\f$
 * \param *dE returned first derivative energy functional at \f$\Theta_{min}\f$
 * \param *ddE returned second derivative energy functional at \f$\Theta_{min}\f$
 * \param *dcos cosine of the found minimal angle
 * \param *dsin sine of the found minimal angle
 * \return \f$\Theta_{min}\f$
 * \sa CalculateLineSearch() - used there
 */
double CalculateDeltaI(double E0, double dEdt0, double ddEddt0, double *E, double *dE, double *ddE, double *dcos, double *dsin) 
{
  double delta, dcos2, dsin2, A, B, Eavg;
  if (fabs(ddEddt0) < MYEPSILON) fprintf(stderr,"CalculateDeltaI: ddEddt0 = %lg",ddEddt0);
  delta = atan(-2.*dEdt0/ddEddt0)/2.0; //.5*
  *dcos = cos(delta);
  *dsin = sin(delta);
  dcos2 = cos(delta*2.0); //2.*
  dsin2 = sin(delta*2.0); //2.*
  A = -ddEddt0/4.;
  B = dEdt0/2.;
  //fprintf(stderr,"CalculateDeltaI: A = %lg, B = %lg\n", A, B);
  Eavg = E0-A;
  *E = Eavg + A*dcos2 + B*dsin2;
  *dE = -2.*A*dsin2+2.*B*dcos2;
  *ddE = -4.*A*dcos2-4.*B*dsin2;
  if (*ddE < 0) { // if it's a maximum (indicated by positive second derivative at found theta)
    if (delta < 0) {  // go the minimum - the two are always separated by PI/2 by the ansatz
      delta += PI/2.0; ///2.
    } else {
      delta -= PI/2.0; ///2.
    }
    *dcos = cos(delta);
    *dsin = sin(delta);
    dcos2 = cos(delta*2.0); //2.*
    dsin2 = sin(delta*2.0); //2.*
    *E = Eavg + A*dcos2 + B*dsin2;
    *dE = -2.*A*dsin2+2.*B*dcos2;
    *ddE = -4.*A*dcos2-4.*B*dsin2;
  }
  return(delta);
}


/** Line search: One-dimensional minimisation along conjugate direction.
 * We want to minimize the energy functional along the calculated conjugate direction.
 * States \f$\psi_i^{(m)}(\Theta) = \widetilde{\psi}_i^{(m)} \cos(\Theta) + \widetilde{\widetilde{\varphi}}_i^{(m)}
 * \sin(\Theta)\f$ are normalized and orthogonal for every \f$\Theta \in {\cal R}\f$. Thus we
 * want to find the minimal \f$\Theta_{min}\f$. Intenting to evaluate the energy functional as
 * seldom as possible, we approximate it by:
 * \f$\widetilde{E}(\Theta) = A + B cos(2\Theta) + CSin(2\Theta)\f$
 * where the functions shall match the energy functional and its first and second derivative at
 * \f$\Theta = 0\f$, fixing the unknowns A, B and C.\n
 * \f$\frac{\delta E}{\delta \Theta}_{|\Theta =0} = 
 *   2 \Re \bigl ( \langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\psi_i^{(m)} \rangle \bigr )\f$
 * and 
 * \f$\frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0} = 
 *  2\bigl( \langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\widetilde{\widetilde{\varphi}}_i^{(m)} \rangle 
 *  - \langle \psi_i^{(m)}|H|\psi_i^{(m)} \rangle
 *  + \underbrace{\int_\Omega \nu_H[\rho](r)\rho(r)dr}_{A_H} 
 *  + \underbrace{\int_\Omega \rho(r)^2 \frac{\delta V^{xc}}{\delta n}(r) dr}_{A_{xc}}\f$. \n
 * \f$A_H\f$ and \f$A_{xc}\f$ can be calculated from the Hartree respectively the exchange
 * correlation energy.\n
 * Basically, this boils down to using CalculateConDirHConDir() (we already have OnePsiAddData::Lambda),
 * in evaluation of the energy functional and its derivatives as shown in the above formulas 
 * and calling with these CalculateDeltaI() which finds the above \f$\Theta_{min} \f$. \n
 * The coefficients of the current wave function are recalculated, thus we have the new wave function,
 * the total energy and its derivatives printed to stderr.\n
 * 
 * \param *P Problem at hand
 * \param *source current minimised wave function
 * \param *ConDir conjugate direction,  \f$\tilde{\tilde{\varphi}}_i^{(m)}\f$
 * \param *Hc_grad complex coefficients of \f$H | \psi_i (G) \rangle\f$, see GradientArray#HcGradient
 * \param x if NULL do approximative line search by second order taylor expansion, otherwise use given
 *        parameter between -M_PI...M_PI. Is simply put through to CalculateLineSearch()
 * \param PsiFactor occupation number
 */
static void CalculateLineSearch(struct Problem *P, 
                                fftw_complex *source,
                                fftw_complex *ConDir,
                                fftw_complex *Hc_grad, const double PsiFactor, const double *x)
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct FileData *F = &P->Files;
  //struct LatticeLevel *Lev0 = R->Lev0;
  //struct Density *Dens = Lev0->Dens;
  //struct Psis *Psi = &P->Lat.Psi;
  struct Energy *En = Lat->E; 
  //int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  double dEdt0, ddEddt0, ConDirHConDir, Eavg, A, B;//, sourceHsource;
  double E0, E, dE, ddE, delta, dcos, dsin, dcos2, dsin2;
  double EI, dEI, ddEI, deltaI, dcosI, dsinI;
  double HartreeddEddt0, XCddEddt0;
  double d[4],D[4], Diff;
  //double factorDR = R->FactorDensityR;
  int g, i;
  //int ActNum;

  // ========= dE / dt | 0 ============
  //fprintf(stderr,"(%i) |ConDir| = %lg, |Hc_grad| = %lg\n", P->Par.me, sqrt(GradSP(P, LevS, ConDir, ConDir)), sqrt(GradSP(P, LevS, Hc_grad, Hc_grad)));
  d[0] = 2.*GradSP(P, LevS, ConDir, Hc_grad);

  // ========== ddE / ddt | 0 =========
  CalculateConDirHConDir(P, ConDir, PsiFactor, &d[1], &d[2], &d[3]);
  
  MPI_Allreduce ( &d, &D, 4, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  dEdt0 = D[0];
  for (i=MAXOLD-1; i > 0; i--)
    En->dEdt0[i] = En->dEdt0[i-1];
  En->dEdt0[0] = dEdt0;
  ConDirHConDir = D[1];
  HartreeddEddt0 = D[2];
  XCddEddt0 = D[3];
  ddEddt0 = 0.0;
  switch (R->CurrentMin) {
    case Occupied:
      ddEddt0 = HartreeddEddt0+XCddEddt0;    // these terms drop out in unoccupied variation due to lacking dependence of n^{tot}
    case UnOccupied:
      ddEddt0 += 2.*(ConDirHConDir - Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda);
      break;
    default:
    ddEddt0 = 0.0;
      break;
  }
  //if (R->CurrentMin <= UnOccupied)
    //fprintf(stderr,"ConDirHConDir = %lg\tLambda = %lg\t Hartree = %lg\t XC = %lg\n",ConDirHConDir,Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda,HartreeddEddt0,XCddEddt0);

  for (i=MAXOLD-1; i > 0; i--)
    En->ddEddt0[i] = En->ddEddt0[i-1];
  En->ddEddt0[0] = ddEddt0;
  E0 = En->TotalEnergy[0];

  if (x == NULL) {
    // delta
    //if (isnan(E0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): E0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
    //if (isnan(dEdt0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): dEdt0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
    //if (isnan(ddEddt0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): ddEddt0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
    deltaI = CalculateDeltaI(E0, dEdt0, ddEddt0, 
                           &EI, &dEI, &ddEI, &dcosI, &dsinI);
  } else {
    deltaI = *x;
    dcosI = cos(deltaI);
    dsinI = sin(deltaI);
    dcos2 = cos(deltaI*2.0); //2.*
    dsin2 = sin(deltaI*2.0); //2.*
    A = -ddEddt0/4.;
    B = dEdt0/2.;
    //fprintf(stderr,"CalculateDeltaI: A = %lg, B = %lg\n", A, B);
    Eavg = E0-A;
    EI = Eavg + A*dcos2 + B*dsin2;
    dEI = -2.*A*dsin2+2.*B*dcos2;
    ddEI = -4.*A*dcos2-4.*B*dsin2;
  }
  delta = deltaI; E = EI; dE = dEI; ddE = ddEI; dcos = dcosI; dsin = dsinI; 
  if (R->ScanPotential >= R->ScanAtStep) {
    delta = (R->ScanPotential-(R->ScanAtStep-1))*R->ScanDelta;
    dcos = cos(delta);
    dsin = sin(delta);
    //fprintf(stderr,"(%i) Scaning Potential form with delta %lg\n", P->Par.me, delta);
  }
  // shift energy delta values
  for (i=MAXOLD-1; i > 0; i--) {
    En->delta[i] = En->delta[i-1];
    En->ATE[i] = En->ATE[i-1]; 
  }
  // store new one
  En->delta[0] = delta;
  En->ATE[0] = E;

  if (R->MinStep != 0)
    Diff = fabs(En->TotalEnergy[1] - E0)/(En->TotalEnergy[1] - E0)*fabs((E0 - En->ATE[1])/E0);
  else
    Diff = 0.;
  
  // rotate local Psi and ConDir according to angle Theta
  for (g = 0; g < LevS->MaxG; g++) {  // Here all coefficients are updated for the new found wave function
    //if (isnan(ConDir[g].re)) { fprintf(stderr,"WARNGING: CalculateLineSearch(): ConDir_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); }
    c_re(source[g]) = c_re(source[g])*dcos + c_re(ConDir[g])*dsin;
    c_im(source[g]) = c_im(source[g])*dcos + c_im(ConDir[g])*dsin;
  }
  // print to stderr
  if (P->Call.out[StepLeaderOut] && x == NULL) {
    if (1)
    switch (R->CurrentMin) {
      case UnOccupied:
        fprintf(stderr, "(%i,%i,%i)S(%i,%i,%i):\tTGE: %e\tATGE: %e\t Diff: %e\t --- d: %e\tdEdt0: %e\tddEddt0: %e\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, Lat->Energy[UnOccupied].TotalEnergy[0], E, Diff, delta, dEdt0, ddEddt0);
        break;
      default:
        fprintf(stderr, "(%i,%i,%i)S(%i,%i,%i):\tTE: %e\tATE: %e\t Diff: %e\t --- d: %e\tdEdt0: %e\tddEddt0: %e\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, Lat->Energy[Occupied].TotalEnergy[0], E, Diff, delta, dEdt0, ddEddt0);
        break;
    }
  } 
  if ((P->Par.me == 0) && (R->ScanFlag || (R->ScanPotential < R->ScanAtStep))) {
    fprintf(F->MinimisationFile, "%i\t%i\t%i\t%e\t%e\t%e\t%e\t%e\n",R->MinStep, R->ActualLocalPsiNo, R->PsiStep, E0, E, delta, dEdt0, ddEddt0);
    fflush(F->MinimisationFile);
  }
}

/** Find the new minimised wave function.
 * Note, the gradient wave the function is expected in GradientTypes#ActualGradient. Then,
 * subsequently calls (performing thereby a conjugate gradient minimisation algorithm
 * of Fletchers&Reeves):
 * - CalculateCGGradient()\n
 *    find direction of steepest descent: GradientTypes#ActualGradient
 * - CalculatePreConGrad()\n
 *    precondition steepest descent direction: GradientTypes#PreConGradient
 * - CalculateConDir()\n
 *    calculate conjugate direction: GradientTypes#ConDirGradient and GradientTypes#OldConDirGradient
 * - CalculateLineSearch()\n
 *    perform the minimisation in this one-dimensional subspace: LevelPsi::LocalPsi[RunStruct::ActualLocalPsiNo]
 * \note CalculateGradient() is not re-done.
 * \param *P Problem at hand
 * \param x if NULL do approximative line search by second order taylor expansion, otherwise use given paramter between -M_PI...M_PI
 */
void CalculateNewWave(struct Problem *P, double *x) 
{
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  const int i = R->ActualLocalPsiNo;
  const int ConDirStep = R->PsiStep;
  if (((!R->ScanPotential) || (R->ScanPotential < R->ScanAtStep)) && (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent != 1)) { // only evaluate CG once
    //fprintf(stderr,"(%i) Evaluating ConDirGradient from SD direction...\n", P->Par.me);
    CalculateCGGradient(P, LevS->LPsi->LocalPsi[i], 
                      P->Grad.GradientArray[ActualGradient], 
                      Lat->Psi.AddData[i].Lambda);
    //fprintf(stderr,"CalculateCGGradient: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[ActualGradient], P->Grad.GradientArray[ActualGradient]));
    //if (isnan(P->Grad.GradientArray[ActualGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): ActualGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
    CalculatePreConGrad(P, P->Grad.GradientArray[ActualGradient], 
                      P->Grad.GradientArray[PreConGradient], 
                      Lat->Psi.AddData[i].T);
    //if (isnan(P->Grad.GradientArray[PreConGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): PreConGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
    //fprintf(stderr,"CalculatePreConGrad: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[PreConGradient], P->Grad.GradientArray[PreConGradient]));
    CalculateConDir(P, ConDirStep, &Lat->Psi.AddData[i].Gamma,
                  P->Grad.GradientArray[ActualGradient],
        P->Grad.GradientArray[OldActualGradient],
                  P->Grad.GradientArray[PreConGradient],
                  P->Grad.GradientArray[ConDirGradient],
        P->Grad.GradientArray[OldConDirGradient]);
    //if (isnan(P->Grad.GradientArray[ConDirGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): ConDirGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
    //fprintf(stderr,"CalculateConDir: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[ConDirGradient], P->Grad.GradientArray[ConDirGradient]));
  } //else 
    //TestForOrth(P,R->LevS,P->Grad.GradientArray[GraSchGradient]);
  //fprintf(stderr,"(%i) DoBrent[%i] = %i\n", P->Par.me, R->ActualLocalPsiNo, Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent);
  if (R->ScanPotential) R->ScanPotential++;  // increment by one
  CalculateLineSearch(P, 
                      LevS->LPsi->LocalPsi[i],
                      P->Grad.GradientArray[ConDirGradient],
                      P->Grad.GradientArray[HcGradient],
                      Lat->Psi.LocalPsiStatus[i].PsiFactor, x);
  //if (isnan(LevS->LPsi->LocalPsi[i][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): {\\wildetilde \\varphi}_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
  SetGramSchExtraPsi(P, &Lat->Psi, NotUsedToOrtho);
  //TestForOrth(P,R->LevS,P->Grad.GradientArray[GraSchGradient]);
}


/** Calculates the parameter alpha for UpdateWaveAfterIonMove()
 * Takes ratios of differences of old and older positions with old and current,
 * respectively with old and older. Alpha is then the ratio between the two:
 * \f[
 *    \alpha = \frac{v_x[old] \cdot v_x[now] + v_y[old] \cdot v_y[now] + v_z[old] \cdot v_z[now]}{v_x[old]^2 + v_y[old]^2 + v_z[old]^2}
 * \f]
 * Something like a rate of change in the ion velocity.
 * \param *P Problem at hand
 * \return alpha
 */
static double CalcAlphaForUpdateWaveAfterIonMove(struct Problem *P) 
{
  struct Ions *I = &P->Ion;
  struct RunStruct *Run = &P->R;
  double alpha=0;
  double Sum1=0;
  double Sum2=0;
  int is,ia;
  double *R,*R_old,*R_old_old;
  Run->AlphaStep++;
  if (Run->AlphaStep <= 1) return(0.0);
  for (is=0; is < I->Max_Types; is++) {
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
      if (I->I[is].IMT[ia] == MoveIon) {
                                R = &I->I[is].R[NDIM*ia];
                                R_old = &I->I[is].R_old[NDIM*ia];
                                R_old_old = &I->I[is].R_old_old[NDIM*ia];  
                                Sum1 += (R_old[0]-R_old_old[0])*(R[0]-R_old[0]);
                                Sum1 += (R_old[1]-R_old_old[1])*(R[1]-R_old[1]);
                                Sum1 += (R_old[2]-R_old_old[2])*(R[2]-R_old[2]);
                                Sum2 += (R_old[0]-R_old_old[0])*(R_old[0]-R_old_old[0]);
                                Sum2 += (R_old[1]-R_old_old[1])*(R_old[1]-R_old_old[1]);
                                Sum2 += (R_old[2]-R_old_old[2])*(R_old[2]-R_old_old[2]);
      }
    }
  }
  if (fabs(Sum2) > MYEPSILON) {
    alpha = Sum1/Sum2;
  } else {
    fprintf(stderr,"(%i) Sum2 <= MYEPSILON\n", P->Par.me);
    alpha = 0.0;
  }
  //if (isnan(alpha)) { fprintf(stderr,"(%i) CalcAlphaForUpdateWaveAfterIonMove: alpha is NaN!\n", P->Par.me); Error(SomeError, "NaN-Fehler!");}
  return (alpha);
}

/** Updates Wave functions after UpdateIonsR.
 * LocalPsi coefficients are shifted by the difference between LevelPsi::LocalPsi and LevelPsi::OldLocalPsi
 * multiplied by alpha, new value stored in old one.
 * Afterwards all local wave functions are set to NotOrthogonal and then orthonormalized via
 * GramSch().
 * 
 * \param *P Problem at hand
 * \sa CalcAlphaForUpdateWaveAfterIonMove() - alpha is calculated here
 */
double UpdateWaveAfterIonMove(struct Problem *P) 
#if 0
static double CalculateDeltaII(double E0, double dEdt0, double E1, double deltastep, double *E, double *dE, double *ddE, double *dcos, double *dsin) 
{
  double delta, dcos2, dsin2, A, B, Eavg;
  if (fabs(1-cos(2*deltastep) < MYEPSILON) fprintf(stderr,"CalculateDeltaII: (1-cos(2*deltastep) = %lg",(1-cos(2*deltastep));
  A = (E0-E1+0.5*dEdt0*sin(2*deltastep))/(1-cos(2*deltastep));
  B = 0.5*dEdt0;
  Eavg = E0-A;
  if (fabs(A) < MYEPSILON) fprintf(stderr,"CalculateDeltaII: A = %lg",A);
  delta = 0.5*atan(B/A);
  *dcos = cos(delta);
  *dsin = sin(delta);
  dcos2 = cos(2.*delta);
  dsin2 = sin(2.*delta);
  *E = Eavg + A*dcos2 + B*dsin2;
  *dE = -2.*A*dsin2+2.*B*dcos2;
  *ddE = -4.*A*dcos2-4.*B*dsin2;
  if (*ddE < 0) {
    if (delta < 0) {
      delta += PI/2.;
    } else {
      delta -= PI/2.;
    }
    *dcos = cos(delta);
    *dsin = sin(delta);
    dcos2 = cos(2.*delta);
    dsin2 = sin(2.*delta);
    *E = Eavg + A*dcos2 + B*dsin2;
    *dE = -2.*A*dsin2+2.*B*dcos2;
    *ddE = -4.*A*dcos2-4.*B*dsin2;
  }
  return(delta);
}
#endif
{
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  double alpha=CalcAlphaForUpdateWaveAfterIonMove(P);
  int i,g,p,ResetNo=0;
  fftw_complex *source, *source_old;
  struct OnePsiElement *OnePsi = NULL;
  //fprintf(stderr, "%lg\n", alpha);
  if (P->R.AlphaStep > 1 && alpha != 0.0) {
    for (i=0; i < Psi->LocalNo; i++) {
      source = LevS->LPsi->LocalPsi[i];
      source_old = LevS->LPsi->OldLocalPsi[i];
      for (g=0; g < LevS->MaxG; g++) {
                                c_re(source[g]) = c_re(source[g])+alpha*(c_re(source[g])-c_re(source_old[g]));
                                c_im(source[g]) = c_im(source[g])+alpha*(c_im(source[g])-c_im(source_old[g]));
                                c_re(source_old[g])=c_re(source[g]);
                                c_im(source_old[g])=c_im(source[g]);
      }
    }
  } else {
    for (i=0; i < Psi->LocalNo; i++) {
      source = LevS->LPsi->LocalPsi[i];
      source_old = LevS->LPsi->OldLocalPsi[i];
      for (g=0; g < LevS->MaxG; g++) {
                                c_re(source_old[g])=c_re(source[g]);
                                c_im(source_old[g])=c_im(source[g]);
      }
    }
  }
  //fprintf(stderr, "%lg\n", alpha);
  if (alpha != 0.0) {
    for (p=0; p< Psi->LocalNo; p++) {
      Psi->LocalPsiStatus[p].PsiGramSchStatus = (Psi->LocalPsiStatus[p].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho;
      //if (R->CurrentMin > UnOccupied) fprintf(stderr,"Setting L-Status of %i to %i\n",p, Psi->LocalPsiStatus[p].PsiGramSchStatus);
    }
    ResetNo=0;
    for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
      for (p=ResetNo; p < ResetNo+Psi->AllLocalNo[i]-1; p++) {
                                Psi->AllPsiStatus[p].PsiGramSchStatus = (Psi->LocalPsiStatus[p].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho;
        //if (R->CurrentMin > UnOccupied) fprintf(stderr,"Setting A-Status of %i to %i\n",p, Psi->AllPsiStatus[p].PsiGramSchStatus);
      }
      ResetNo += Psi->AllLocalNo[i];
      OnePsi = &Psi->AllPsiStatus[ResetNo-1]; // extra wave function is set to NotUsed
      //if (R->CurrentMin > UnOccupied) fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo-1, NotUsedToOrtho);
      OnePsi->PsiGramSchStatus =  (int)NotUsedToOrtho;
    }
    SpeedMeasure(P, InitGramSchTime, StartTimeDo);
    //fprintf(stderr,"UpdateWaveAfterIonMove: ResetGramSch() for %i orbital\n",p);
    GramSch(P, LevS, Psi, Orthonormalize);
    SpeedMeasure(P, InitGramSchTime, StartTimeDo);
  }
  return(alpha);
}