/** \file pseudo.c
 * PseudoPotentials to approximate core densities for plane waves.
 * Note: Pseudopotentials replace the electron-core-interaction potential \f$V^{El-K}\f$.
 * 
 * Herein are various functions dealing with the Pseudopotential approximation, which
 * calculate local and non-local form factors FormFacGauss(), FormFacLocPot(), FormFacNonLocPot(),
 * energies CalculateNonLocalEnergyNoRT(), CalculateIonNonLocalForce(), CalculateSomeEnergyAndHGNoRT()
 * and forces CalculateIonLocalForce(), CalculateIonNonLocalForce(),
 * also small functions for evaluating often used expressions: CalcExpiGR(), dexpiGRpkg(), Get_pkga_MinusG(),
 * CalculateCDfnl(), CalcSG(), CalcCoreCorrection() - some of these are updated in case of new waves or
 * finer mesh by calling UpdatePseudoToNewWaves(), ChangePseudoToLevUp(),
 * for allocation and parameter readin InitPseudoRead() and memory deallocation
 * RemovePseudoRead().
 * 
 * Missing so far are: CalculateAddNLPot()
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers
 \date 2000

  File: pseudo.c
  $Id: pseudo.c,v 1.51 2007-03-29 13:38:47 foo Exp $
*/

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


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

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

#include "data.h"
#include "errors.h"
#include "helpers.h"
#include "ions.h"
#include "init.h"
#include "mymath.h"
#include "myfft.h"
#include "pseudo.h"
#include "run.h"

/** Calculates structure factor \f$S_{I_s} (G)\f$.
 * \f[
 *    S_{I_s} (G) = \sum_{I_a} \exp(i G\cdot R_{I_s,I_a}) \qquad (4.14)
 * \f]
 * \param *P Problem at hand
 */
static void CalcSG(struct Problem *P) 
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  fftw_complex dS;
  double dSPGRi;
  int it,iot,g;
  struct OneGData *G = R->Lev0->GArray;
  for (it=0; it < I->Max_Types; it++)
    for (g=0; g < Lev0->MaxG; g++) {
      c_re(dS) = 0;
      c_im(dS) = 0;
      for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
      	dSPGRi = RSP3(G[g].G,&I->I[it].R[NDIM*iot]);
      	c_re(PP->ei[it][iot][g]) = cos(dSPGRi);
      	c_im(PP->ei[it][iot][g]) = -sin(dSPGRi);
      	c_re(dS) += c_re(PP->ei[it][iot][g]);
      	c_im(dS) += c_im(PP->ei[it][iot][g]);
      }
      c_re(I->I[it].SFactor[g]) = c_re(dS);
      c_im(I->I[it].SFactor[g]) = c_im(dS);
    }
}


/** Calculates \f$\exp(-i(G\cdot R)) = \cos(GR) - i\cdot\sin(GR)\f$.
 * Fills PseudoPot::expiGR array.
 * \param *P Problem at hand
 * \sa expiGR() for element retrieval
 */
static void CalcExpiGR(struct Problem *P) 
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  double dSPGRi;
  int it,iot,g;
  struct OneGData *G = LevS->GArray;
  for (it=0; it < I->Max_Types; it++)
    for (iot=0; iot < I->I[it].Max_IonsOfType; iot++)
      for (g=0; g < LevS->MaxG; g++) {
      	dSPGRi = RSP3(G[g].G,&I->I[it].R[NDIM*iot]);
      	c_re(PP->expiGR[it][iot][g]) = cos(dSPGRi);
      	c_im(PP->expiGR[it][iot][g]) = -sin(dSPGRi);
      }
}

/** Calculates \f$\exp(-i(G)R_{I_s,I_a})\f$.
 * \param *PP PseudoPot ential structure
 * \param it ion type \f$I_s\f$
 * \param iot ion number within certain type (ion of type) \f$I_a\f$
 * \param g reciprocal grid vector from GArray
 * \param *res complex return value
 * \sa CalcExpiGR() fills the array
 */
#ifdef HAVE_INLINE
inline static void expiGR(struct PseudoPot *PP, int it, int iot, int g, fftw_complex *res) {
#else
static void expiGR(struct PseudoPot *PP, int it, int iot, int g, fftw_complex *res) {
#endif
  c_re(res[0]) = c_re(PP->expiGR[it][iot][g]);
  c_im(res[0]) = c_im(PP->expiGR[it][iot][g]);
}   

/** Calculates \f$\exp(-i(G)R_{I_s,I_a})\phi_{I_s,l,m}^{ps,nl} (G)\f$.
 * \param *PP PseudoPot ential structure
 * \param it ion type \f$I_s\f$
 * \param iot ion number within certain type (ion of type), \f$I_a\f$
 * \param ilm m value for ion type
 * \param g reciprocal grid vector from GArray
 * \param *res complex return value
 */
#ifdef HAVE_INLINE
inline static void dexpiGRpkg(struct PseudoPot *PP, int it, int iot, int ilm, int g, fftw_complex *res) {
#else
static void dexpiGRpkg(struct PseudoPot *PP, int it, int iot, int ilm, int g, fftw_complex *res) {
#endif
  expiGR(PP, it, iot, g, res); 
  c_re(res[0]) = c_re(res[0])*PP->phi_ps_nl[it][ilm-1][g];
  c_im(res[0]) = -c_im(res[0])*PP->phi_ps_nl[it][ilm-1][g];
}

/** Calculates Gaussian form factor \f$\phi^{gauss}_{I_s}\f$.
 * \f$I_s\f$ ion types, V volume of cell, \f$r_I^{Gauss}\f$ is defined by the
 * core charge within via \f$n^{Gauss}\f$, G reciprocal grid vector 
 * \f[
 * 	\phi^{gauss}_{I_s} = - \frac{Z_{I_s}}{V} \exp(-\frac{1}{4} (r^{Gauss}_{I_s})^2 |G|^2)
 *  \qquad (4.23)
 * \f]
 * \param *P Problem at hand
 */
static void FormFacGauss(struct Problem *P) 
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  double dFac;
  int it, g;
  struct OneGData *G = Lev0->GArray;
  for (it=0; it < I->Max_Types; it++) {	// for each all types ...
    dFac = 0.25*I->I[it].rgauss*I->I[it].rgauss;
    for (g=0; g < Lev0->MaxG; g++) {	// for each reciprocal grid vector ...
      PP->FacGauss[it][g] = -PP->zval[it]/Lat->Volume*exp(-dFac*G[g].GSq);
    }
  }
}

/** Calculates local PseudoPot form factor \f$\phi_{I_s}^{ps,loc} (G)\f$.
 * \f$I_s\f$ ion types, \f$I_a\f$ ion number of type, V volume of cell, 
 * \f$r_I^{Gauss}\f$ is defined by the core charge within via \f$n^{Gauss}\f$,
 * G reciprocal grid vector
 * \f[
 * 	\phi_{I_s}^{ps,loc} (G) = \frac{4}{\pi} \int_0^\infty r^2 j_0(r|G|)
 * 		\Bigr ( \nu_{I_s}^{loc} (r) + \frac{Z_{I_s}}{r}erf\bigr ( \frac{r}
 * 		{r_{I_s}^{Gauss}} \bigl )\Bigl )
 *  \qquad (4.15)
 * \f]
 * \param *P Problem at hand
 * \note profiling says, sin()/cos() is biggest cpu hog
 */
static void FormFacLocPot(struct Problem *P) 
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  double vv,dint,darg;
  int I_s, ll, ir, g;                   // I_s = iontype, ir = ion radial pos, s = counter for g (0 treated specially), g index for Grid vectors GArray
  struct OneGData *G = Lev0->GArray;
  for (I_s=0; I_s < I->Max_Types; I_s++) {
    ll = I->I[I_s].l_loc-1;
    for (ir = 0; ir < PP->mmax[I_s]; ir++) {	 //fill integrand1 function array with its values: Z/r * erf(r/R_g)
      if (fabs(PP->r[I_s][ir]) < MYEPSILON) fprintf(stderr,"FormFacLocPot: PP->r[it][ir] = %lg\n",PP->r[I_s][ir]);
      if (fabs(I->I[I_s].rgauss) < MYEPSILON) fprintf(stderr,"FormFacLocPot: I->I[it].rgauss = %lg\n",I->I[I_s].rgauss);
      vv = -PP->zval[I_s]/PP->r[I_s][ir] * derf(PP->r[I_s][ir]/I->I[I_s].rgauss);
      PP->integrand1[ir] = PP->v_loc[I_s][ll][ir]-vv;
    }
    g=0;
    if (Lev0->GArray[0].GSq == 0.0) {	// if grid vectors contain zero vector, treat special
      for (ir = 0; ir < PP->mmax[I_s]; ir++) {	// fill the to be integrated function array with its values
		    PP->integrand[ir] = PP->integrand1[ir]*PP->r[I_s][ir]*PP->r[I_s][ir]*PP->r[I_s][ir];
      }
      dint = Simps(PP->mmax[I_s],PP->integrand,PP->clog[I_s]);	// do the integration
      PP->phi_ps_loc[I_s][0] = dint*4.*PI/Lat->Volume;
      g++;
    }
    for (; g < Lev0->MaxG; g++) { 	// pre-calculate for the rest of all the grid vectors
      for (ir = 0; ir < PP->mmax[I_s]; ir++) {   // r^3 * j_0 * dummy1
    		darg = PP->r[I_s][ir]*sqrt(G[g].GSq);
        if (fabs(darg) < MYEPSILON) fprintf(stderr,"FormFacLocPot: darg = %lg\n",darg);
        PP->integrand[ir] = PP->integrand1[ir]*sin(darg)/darg*PP->r[I_s][ir]*PP->r[I_s][ir]*PP->r[I_s][ir];
      }
      dint = Simps(PP->mmax[I_s],PP->integrand,PP->clog[I_s]);
      PP->phi_ps_loc[I_s][g] = dint*4.*PI/Lat->Volume;
    }
  }
}

/** Determines whether its a symmetric or antisymmetric non-local pseudopotential form factor pkga[][value][].
 * \param value pgga l value
 * \return 1 - symmetric (0,4,5,6,7,8) , -1 - antisymmetric (1,2,3)
 */
inline static int Get_pkga_MinusG(const int value)
{ 
  switch (value) {
  case 0: return(1);
  case 1:
  case 2:
  case 3: return(-1);
  case 4:
  case 5:
  case 6:
  case 7:
  case 8: return(1);
  default:
    Error(SomeError,"Get_pkga_MinusG: Unknown Value");
  }
  return 0;
}

/** Calculates non-local PseudoPot form factor \f$\phi_{I_s,l,m}^{ps,nl} (G)\f$ and non-local weight \f$w^{nl}_{I_s,l}\f$.
 * \f$I_s\f$ ion types, \f$I_a\f$ ion number of type, V volume of cell, 
 * \f$j_l\f$ Bessel function of l-th order,
 * G reciprocal grid vector
 * \f[
 * 	\phi_{I_s,l,m}^{ps,nl} (G) = \sqrt{\frac{4\pi}{2l+1}} \int_0^\infty r^2 j_l(|G|r)
 * 		\Delta \nu^{nl}_{I_s} (r) R_{I_s,l}(r) y_{lm} (\vartheta_G, \varphi_G) dr
 *  \qquad (4.17)
 * \f]
 * \f[
 *  w^{nl}_{I_s,l} = \frac{4\pi}{V} (2l+1) \Bigr ( \int^\infty_0 r^2 R_{I_s,l} (r)
 *      \Delta \nu^{nl}_{I_s,l} (r) R_{I_s,l} (r) dr \Bigl )^{-1} 
 *  \qquad (4.18)
 * \f]
 * \param *P Problem at hand
 * \note see e.g. MathWorld for the spherical harmonics in cartesian coordinates
 */
static void FormFacNonLocPot(struct Problem *P) 
{
  int I_s,ir,j,ll,il,ml,g,i;
  double delta_v,arg,ag,fac,cosx,cosy,cosz,dint,sqrt3=sqrt(3.0);
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct OneGData *G = LevS->GArray;
  for (I_s=0; I_s < I->Max_Types; I_s++) { // through all ion types
    ml = -1;        // m value
    for (il=0; il < I->I[I_s].l_max; il++) { // through all possible l values
      ll = I->I[I_s].l_loc-1;                // index for the local potential
      if (il != ll) {
        // calculate non-local weights
				for (ir = 0; ir < PP->mmax[I_s]; ir++) {
				  delta_v = PP->v_loc[I_s][il][ir]-PP->v_loc[I_s][ll][ir];   // v_l - v_l_loc
				  PP->integrand2[I_s][ir] = delta_v*PP->R[I_s][il][ir]*PP->r[I_s][ir]*PP->r[I_s][ir];    // for form factor
				  PP->integrand1[ir] = delta_v*PP->R[I_s][il][ir]*PP->R[I_s][il][ir]*PP->r[I_s][ir];    // for weight
				} 
				dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
        if (fabs(dint) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: dint[%d] = %lg\n",I_s,dint);
				if (il+1 < 4) {
				  for(j =1; j <= 2*(il+1)-1; j++) {
				    PP->wNonLoc[I_s][ml+j]=(2.0*(il+1)-1.0)*4.*PI/Lat->Volume/dint;  // store weight
				  }
				}
        // calculate non-local form factors
				for (g=0; g < LevS->MaxG; g++) { 
				  arg = sqrt(G[g].GSq);
				  switch (il) {     // switch on the order of the needed bessel function, note: only implemented till second order
				  case 0: 
				    if (arg == 0.0) {   // arg = 0 ...
				      for (ir=0; ir < PP->mmax[I_s]; ir++) {
								PP->integrand1[ir] = PP->integrand2[I_s][ir];
				      }
				    } else {    // ... or evaluate J_0
				      for (ir=0; ir < PP->mmax[I_s]; ir++) {
								fac=arg*PP->r[I_s][ir];
                if (fabs(fac) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: fac[%d][%d] = %lg\n",I_s,ir,fac);
								fac=sin(fac)/fac; // CHECKED: j_0 (x)
								PP->integrand1[ir] = fac*PP->integrand2[I_s][ir];
				      }
				    }
				    dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
				    PP->phi_ps_nl[I_s][ml+1][g]=dint;      // m = 0
				    break;
				  case 1: 
				    if (arg == 0.0) {   // arg = 0 ...
				      for (j=1; j <= 3; j++) {
								PP->phi_ps_nl[I_s][ml+j][g]=0.0;
				      }
				    } else {   // ... or evaluate J_1
				      for (ir=0; ir < PP->mmax[I_s]; ir++) {
								fac=arg*PP->r[I_s][ir];
                if (fabs(fac) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: fac[%d][%d] = %lg\n",I_s,ir,fac);
								fac=(sin(fac)/fac-cos(fac))/fac;  // CHECKED: j_1(x)
								PP->integrand1[ir] = fac*PP->integrand2[I_s][ir];
				      }
				      dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
				      PP->phi_ps_nl[I_s][ml+1][g]=dint*G[g].G[0]/arg; // m = -1
				      PP->phi_ps_nl[I_s][ml+2][g]=dint*G[g].G[1]/arg; // m = +1
				      PP->phi_ps_nl[I_s][ml+3][g]=dint*G[g].G[2]/arg; // m = 0
				    } 
				    break;
				  case 2:
				    if (arg == 0.0) {   // arg = 0 ...
				      for (j=1; j <= 5; j++) {
								PP->phi_ps_nl[I_s][ml+j][g]=0.0;
				      }
				    } else {   // ... or evaluate J_2
				      for (ir=0; ir < PP->mmax[I_s]; ir++) {
								ag=arg*PP->r[I_s][ir];
                if (fabs(ag) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: ag[%d][%d] = %lg\n",I_s,ir,ag);
								fac=(3.0/(ag*ag)-1.0)*sin(ag)-3.0/ag*cos(ag);
								fac=fac/ag;   // CHECKED: j_2(x)
								PP->integrand1[ir] = fac*PP->integrand2[I_s][ir];
				      }
				      dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
              if (fabs(arg) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: arg = %lg\n",arg);
				      cosx = G[g].G[0]/arg;
				      cosy = G[g].G[1]/arg;
				      cosz = G[g].G[2]/arg;
				      PP->phi_ps_nl[I_s][ml+1][g]=dint*0.5*(3.0*cosz*cosz-1.0);    // m = 0
				      PP->phi_ps_nl[I_s][ml+2][g]=dint*sqrt3*cosz*cosx; // m = +1
				      PP->phi_ps_nl[I_s][ml+3][g]=dint*sqrt3*cosz*cosy; // m = -1
				      PP->phi_ps_nl[I_s][ml+4][g]=dint*sqrt3*cosx*cosy; // m = -2
				      PP->phi_ps_nl[I_s][ml+5][g]=dint*sqrt3/2.0*(cosx*cosx-cosy*cosy); // m = +2
				    }
				    break;
          default:
            Error(SomeError, "FormFacNonLocPot: Order of sperhical Bessel function not implemented!");
            break;
				  }
          // phi_ps_nl_a removed, not needed, just a mem hog
//				  for (j=1; j<=2*(il+1)-1;j++) {
//				    PP->phi_ps_nl_a[I_s][ml+j][g] = PP->phi_ps_nl[I_s][ml+j][g];
//				  }
				}
				ml += (2*(il+1)-1);
      }
    }
    for (i=0; i < I->I[I_s].l_max*I->I[I_s].l_max-2*I->I[I_s].l_loc+1; i++) 
      if(P->Call.out[ReadOut]) 
				fprintf(stderr,"wnl:(%i,%i) = %g\n",I_s,i,PP->wNonLoc[I_s][i]);
  }
}

/** Calculates partial core form factors \f$\phi^{pc}_{I_s}(G)\f$ and \f$\hat{n}^{pc} (G)\f$.
 * \f[
 *    \phi^{pc}_{I_s}(G) = \frac{4\pi}{V} \int_0^\infty r^2 j_0(|G|r) n^{pc}_{I_s} (r) dr \qquad (4.21)
 * \f]
 * \f[
 *    \hat{n}^{pc} = \sum_{I_s} S(G) \phi^{pc}_{I_s}(G)
 * \f]
 * \param *P Problem at hand
 * \param Create 1 - calculate form factors, 0 - don't
 */
static void CalcCoreCorrection(struct Problem *P, const int Create)
{
  int i, s = 0;
  int it,g,Index,ir;
  double arg;
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct OneGData *G = Lev0->GArray;
  struct Density *Dens = Lev0->Dens;
  struct fft_plan_3d *plan = Lat->plan; 
  fftw_complex *work = (fftw_complex *)Dens->DensityCArray[TempDensity];
  fftw_complex *CoreWave = (fftw_complex *)Dens->DensityArray[CoreWaveDensity];
  // if desired the form factors are calculated
  if (Create) {
    for (it=0; it < I->Max_Types; it++) {
      if (I->I[it].corecorr == CoreCorrected) {
      	s=0;
      	if (Lev0->GArray[0].GSq == 0.0) {
      	  for (ir=0; ir < PP->mmax[it]; ir++) {
      	    PP->integrand[ir] = 4.*PI*PP->corewave[it][ir]*PP->r[it][ir]*PP->r[it][ir]*PP->r[it][ir]/Lat->Volume;
      	  }
      	  PP->formfCore[it][0] = Simps(PP->mmax[it], PP->integrand, PP->clog[it]);
      	  s++;
      	}
      	for (g=s; g <  Lev0->MaxG; g++) {
      	  for (ir=0; ir < PP->mmax[it]; ir++) {
      	    arg=PP->r[it][ir]*sqrt(G[g].GSq);
            if (fabs(arg) < MYEPSILON) fprintf(stderr,"CalcCoreCorrection: arg = %lg\n",arg);
      	    PP->integrand1[ir] = PP->integrand[ir]*sin(arg)/arg;
      	  }
      	  PP->formfCore[it][g] = Simps(PP->mmax[it], PP->integrand1, PP->clog[it]); 
      	} 
      }
    }
  }
  // here the density is evaluated, note again that it has to be enfolded due to the finer resolution
  SetArrayToDouble0((double *)CoreWave,2*Dens->TotalSize);
  for (g=0; g < Lev0->MaxG; g++) {
    Index = Lev0->GArray[g].Index;
    for (it = 0; it < I->Max_Types; it++) {
      if (I->I[it].corecorr == CoreCorrected) {
      	c_re(CoreWave[Index]) += PP->formfCore[it][g]*c_re(I->I[it].SFactor[g]);
      	c_im(CoreWave[Index]) += PP->formfCore[it][g]*c_im(I->I[it].SFactor[g]);
      }
    }
  }
  for (i=0; i<Lev0->MaxDoubleG; i++) {
    CoreWave[Lev0->DoubleG[2*i+1]].re =  CoreWave[Lev0->DoubleG[2*i]].re;
    CoreWave[Lev0->DoubleG[2*i+1]].im = -CoreWave[Lev0->DoubleG[2*i]].im;
  }
  fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF1, CoreWave, work);
}

/** Reads data for Pseudopotentials.
 * PseudoPot structure is allocated and intialized, depending on the given Z values in the 
 * config file, corresponding pseudopotential files are read (generated from FHMD98 package),
 * basically calls ChangePseudoToLevUp() (only it's implemented by hand).
 * \param *P Problem at hand
 * \param *pseudopot_path Path to pseudopot files
 * \sa ParseForParameter() for parsing, RemovePseudoRead() for deallocation
 */
void InitPseudoRead(struct Problem *P, char *pseudopot_path) 
{
  char cpiInputFileName[255];
  FILE *cpiInputFile;
  int i,it,ib,il,m,j,ia;
  int count = 0;
  double dummycorewave = 0., dummyr = 0.;
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *ILev0 = R->InitLev0;
  struct LatticeLevel *ILevS = R->InitLevS;
  PP->Mmax = 0;
  PP->core = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->rc = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->rcl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
  PP->al = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
  PP->bl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
  PP->mmax = (int *) Malloc(sizeof(int)*I->Max_Types, "InitPseudoRead: ");
  PP->clog = (double *) Malloc(sizeof(double)*I->Max_Types, "InitPseudoRead: ");
  PP->R = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
  PP->r = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->integrand2 = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->v_loc = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
  PP->zval = (double *)Malloc(sizeof(double)*I->Max_Types, "InitPseudoRead: ");
  PP->nang = (int *) Malloc(sizeof(int)*I->Max_Types, "InitPseudoRead: ");
  PP->FacGauss = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->phi_ps_loc = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->wNonLoc = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  //PP->phi_ps_nl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
  PP->phi_ps_nl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
  PP->lm_end = (int *) Malloc(sizeof(int)*I->Max_Types, "InitPseudoRead: ");
  PP->ei = (fftw_complex ***) Malloc(sizeof(fftw_complex**)*I->Max_Types, "InitPseudoRead: ");
  PP->expiGR = (fftw_complex ***) Malloc(sizeof(fftw_complex**)*I->Max_Types, "InitPseudoRead: ");
  PP->VCoulombc = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
  PP->corewave = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->formfCore = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
  PP->lm_endmax = 0;
  PP->corecorr = NotCoreCorrected;
  for (it = 0; it < I->Max_Types; it++) {
    PP->lm_end[it] = I->I[it].l_max*I->I[it].l_max+1-2*I->I[it].l_loc;
    if (PP->lm_endmax < PP->lm_end[it]) PP->lm_endmax = PP->lm_end[it];  
    //PP->phi_ps_nl[it] =  (double **) Malloc(sizeof(double*)*PP->lm_end[it], "InitPseudoRead: ");
    PP->phi_ps_nl[it] = (double **) Malloc(sizeof(double*)*PP->lm_end[it], "InitPseudoRead: ");
    PP->ei[it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
    PP->expiGR[it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
    for (ia=0;ia<I->I[it].Max_IonsOfType;ia++) {
      PP->ei[it][ia] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
      PP->expiGR[it][ia] =  (fftw_complex *) Malloc(sizeof(fftw_complex)*ILevS->MaxG, "InitPseudoRead: ");
    }
    for (il=0;il<PP->lm_end[it];il++) {
      //PP->phi_ps_nl[it][il] = (double *) Malloc(sizeof(double)*ILevS->MaxG, "InitPseudoRead: ");
      PP->phi_ps_nl[it][il] = (double *) Malloc(sizeof(double)*ILevS->MaxG, "InitPseudoRead: ");
    }
    PP->wNonLoc[it] = (double *) Malloc(sizeof(double)*PP->lm_end[it], "InitPseudoRead: ");
    sprintf(cpiInputFileName,"%s/pseudo.%02i",pseudopot_path,I->I[it].Z); 
    cpiInputFile = fopen(cpiInputFileName,"r");
    if (cpiInputFile == NULL)
      Error(SomeError, "Cant't find pseudopot path or file");
    else
    	if (P->Call.out[ReadOut]) fprintf(stderr,"%s was found\n",cpiInputFileName);
    int zeile = 1;
    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->zval[it], 1, critical);
    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, int_type, &PP->nang[it], 1, critical);
    I->TotalZval += PP->zval[it]*(I->I[it].Max_IonsOfType);
    PP->core[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
    PP->rc[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->core[it][0], 1, critical);
    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->rc[it][0], 1, critical);
    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->core[it][1], 1, critical);
    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->rc[it][1], 1, critical);
    PP->rcl[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
    PP->al[it] =  (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
    PP->bl[it] =  (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
    PP->FacGauss[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
    PP->phi_ps_loc[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
    I->I[it].SFactor = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
    for (il=0; il < 3; il++) {
      PP->rcl[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
      PP->al[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
      PP->bl[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
      for (ib = 0; ib < 3; ib++) {
		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->rcl[it][il][ib], 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->al[it][il][ib], 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 1, 3, zeile, double_type, &PP->bl[it][il][ib], 1, critical);
				if (PP->rcl[it][il][ib] < MYEPSILON)
				  PP->rcl[it][il][ib] = 0.0;
				else
				  PP->rcl[it][il][ib] = 1./sqrt(PP->rcl[it][il][ib]);
      }
    }
    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
    if (PP->mmax[it] > PP->Mmax) PP->Mmax = PP->mmax[it];
    PP->clog[it] = log(PP->clog[it]);
    PP->r[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
    PP->integrand2[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
    PP->R[it] = (double **) Malloc(sizeof(double*)*PP->nang[it], "InitPseudoRead: ");
    PP->v_loc[it] = (double **) Malloc(sizeof(double*)*PP->nang[it], "InitPseudoRead: ");
    for (il=0; il < PP->nang[it]; il++) {
      PP->R[it][il] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
      PP->v_loc[it][il] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
      for (j=0;j< PP->mmax[it];j++) {
		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &m, 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->r[it][j], 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->R[it][il][j], 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->v_loc[it][il][j], 1, critical);
      }
      if (il < PP->nang[it]-1) {
		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
				PP->clog[it] = log(PP->clog[it]);
      }
    }
    I->I[it].corecorr = NotCoreCorrected;
    count = 0;
    count += ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &dummyr, 1, optional);
    count += ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &dummycorewave, 1, optional);
    //fprintf(stderr,"(%i) %lg %lg\n",P->Par.me,dummyr,dummycorewave);
    if (count == 2) {
      PP->r[it][0] = dummyr;
      PP->corewave[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
      PP->formfCore[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
      I->I[it].corecorr = CoreCorrected;
      PP->corecorr = CoreCorrected;
      PP->corewave[it][0] = dummycorewave/(4.*PI);
	    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
	    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
      for (j=1;j < PP->mmax[it]; j++) {
		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->r[it][j], 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->corewave[it][j], 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
		    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
				PP->corewave[it][j] /= (4.*PI);
      }
    } else {
      PP->corewave[it] = NULL;
      PP->formfCore[it] = NULL;
    }
    fclose(cpiInputFile);
  }
  PP->fnl = (fftw_complex ****) Malloc(sizeof(fftw_complex***)*(Lat->Psi.LocalNo+1), "InitPseudoRead: ");
  for (i=0; i< Lat->Psi.LocalNo+1; i++) {
    PP->fnl[i] = (fftw_complex ***) Malloc(sizeof(fftw_complex**)*I->Max_Types, "InitPseudoRead: ");
    for (it=0; it < I->Max_Types; it++) {
      PP->fnl[i][it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*PP->lm_end[it], "InitPseudoRead: ");
      for (il =0; il < PP->lm_end[it]; il++) 
				PP->fnl[i][it][il] = (fftw_complex *) Malloc(sizeof(fftw_complex)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
    }
  }
  if (fabs(I->TotalZval - P->Lat.Psi.NIDensity[Occupied]) >= MYEPSILON) {
    if (P->Par.me_comm_ST == 0) // instead of P->Par.me, thus differing charge density output also for SpinUp/Down from both processes!
      fprintf(stderr, "TotalZval %i != NIDensity[0] %g eps (%g >= %g)\n",I->TotalZval, P->Lat.Psi.NIDensity[Occupied], fabs(I->TotalZval - P->Lat.Psi.NIDensity[Occupied]),MYEPSILON);
    Error(SomeError,"Readparameters: charged system is not implemented yet");
  }
  PP->CDfnl = PP->fnl[Lat->Psi.LocalNo];
  PP->dfnl = (fftw_complex *) Malloc(sizeof(fftw_complex)*I->Max_Max_IonsOfType, "InitPseudoRead: ");
  PP->rr = (double *) Malloc(sizeof(double)*PP->lm_endmax, "InitPseudoRead: ");
  PP->t = (fftw_complex *) Malloc(sizeof(fftw_complex)*PP->lm_endmax, "InitPseudoRead: ");
  PP->integrand = (double *) Malloc(sizeof(double)*PP->Mmax, "InitPseudoRead: ");
  PP->integrand1 = (double *) Malloc(sizeof(double)*PP->Mmax, "InitPseudoRead: ");
  for (it=0;it < I->Max_Types; it++) {
    if (I->I[it].corecorr == CoreCorrected) {
      for (j=0; j < PP->mmax[it]; j++) {
				PP->integrand[j] = 4.*PI*PP->corewave[it][j]*PP->r[it][j]*PP->r[it][j]*PP->r[it][j];
      }
      if(P->Call.out[ReadOut]) fprintf(stderr,"Ion[%i] is core corrected: core charge = %g\n", it, Simps(PP->mmax[it], PP->integrand, PP->clog[it]));
    }
  }
  PP->fac1sqrt2PI = 1./sqrt(2.*PI);
  PP->fac1sqrtPI = 1./sqrt(PI);
  SpeedMeasure(P, InitSimTime, StartTimeDo);
  SpeedMeasure(P, InitLocTime, StartTimeDo);
  CalcSG(P);
  CalcExpiGR(P);
  FormFacGauss(P);
  FormFacLocPot(P);
  SpeedMeasure(P, InitLocTime, StopTimeDo);
  SpeedMeasure(P, InitNonLocTime, StartTimeDo);
  FormFacNonLocPot(P); 
  SpeedMeasure(P, InitNonLocTime, StopTimeDo);
  SpeedMeasure(P, InitLocTime, StartTimeDo);
  CalcCoreCorrection(P, 1);
  SpeedMeasure(P, InitLocTime, StopTimeDo);
  SpeedMeasure(P, InitSimTime, StopTimeDo);
}

/** Updates Pseudopotentials due to change of mesh width.
 * Calls CalcSG(), CalcExpiGR(), recalculates the form factors by calling
 * FormFacGauss(), FormFacLocPot(), FormFacLonLocPot() and finally
 * CalcCoreCorrection(). All speed-measured in InitLocTime respectively
 * InitNonLocTime.
 * \param *P Problem at hand
 */
inline void ChangePseudoToLevUp(struct Problem *P)
{
  SpeedMeasure(P, InitLocTime, StartTimeDo);
  CalcSG(P);
  CalcExpiGR(P);
  FormFacGauss(P);
  FormFacLocPot(P);
  SpeedMeasure(P, InitLocTime, StopTimeDo);
  SpeedMeasure(P, InitNonLocTime, StartTimeDo);
  FormFacNonLocPot(P); 
  SpeedMeasure(P, InitNonLocTime, StopTimeDo);
  SpeedMeasure(P, InitLocTime, StartTimeDo);
  CalcCoreCorrection(P, 1);
  SpeedMeasure(P, InitLocTime, StopTimeDo);
}

/** Updates Pseudopotentials due to the newly minimized wave functions.
 * Calls CalcSG(), CalcExpiGR() and CalcCoreCorrection().
 * \param *P Problem at hand
 */
inline void UpdatePseudoToNewWaves(struct Problem *P)
{
  SpeedMeasure(P, InitLocTime, StartTimeDo);
  CalcSG(P);
  CalcExpiGR(P);
  CalcCoreCorrection(P, 0);
  SpeedMeasure(P, InitLocTime, StopTimeDo); 
}

/** Frees memory allocated from Readin.
 * All memory is free'd that was allocated in InitPseudoRead()
 * \param *P Problem at hand
 * \sa InitPseudoRead()
 */
void RemovePseudoRead(struct Problem *P)
{
  int i,it,il,ia;
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  Free(PP->integrand1, "RemovePseudoRead: PP->integrand1");
  Free(PP->integrand, "RemovePseudoRead: PP->integrand");
  Free(PP->dfnl, "RemovePseudoRead: PP->dfnl"); 
  Free(PP->rr, "RemovePseudoRead: PP->rr");
  Free(PP->t, "RemovePseudoRead: PP->t");
  for (i=0; i< Lat->Psi.LocalNo+1; i++) {
    for (it=0; it < I->Max_Types; it++) {
      for (il =0; il < PP->lm_end[it]; il++) 
      	Free(PP->fnl[i][it][il], "RemovePseudoRead: PP->fnl[i][it][il]");
      Free(PP->fnl[i][it], "RemovePseudoRead: PP->fnl[i][it]");
    }
    Free(PP->fnl[i], "RemovePseudoRead: PP->fnl[i]"); 
  }
  for (it=0; it < I->Max_Types; it++) {
    if (I->I[it].corecorr == CoreCorrected) {
      Free(PP->corewave[it], "RemovePseudoRead: PP->corewave[it]");
      Free(PP->formfCore[it], "RemovePseudoRead: PP->formfCore[it]");
    }
  }
  for (it=0; it < I->Max_Types; it++) {
    for (il=0; il < PP->nang[it]; il++) {
      Free(PP->R[it][il], "RemovePseudoRead: PP->R[it][il]");
      Free(PP->v_loc[it][il], "RemovePseudoRead: PP->v_loc[it][il]");
    }
    for (il=0; il < 3; il++) {
      Free(PP->rcl[it][il], "RemovePseudoRead: PP->rcl[it][il]");
      Free(PP->al[it][il], "RemovePseudoRead: PP->al[it][il]");
      Free(PP->bl[it][il], "RemovePseudoRead: PP->bl[it][il]");
    }
    for (il=0;il<PP->lm_end[it];il++) {
      //Free(PP->phi_ps_nl[it][il], "RemovePseudoRead: PP->phi_ps_nl[it][il]");
      Free(PP->phi_ps_nl[it][il], "RemovePseudoRead: PP->phi_ps_nl[it][il]");
    }
    for (ia=0;ia<I->I[it].Max_IonsOfType;ia++) {
      Free(PP->ei[it][ia], "RemovePseudoRead: PP->ei[it][ia]");
      Free(PP->expiGR[it][ia], "RemovePseudoRead: PP->expiGR[it][ia]");
    }
    Free(PP->ei[it], "RemovePseudoRead: PP->ei[it]");
    Free(PP->expiGR[it], "RemovePseudoRead: PP->expiGR[it]");
    //Free(PP->phi_ps_nl[it], "RemovePseudoRead: PP->phi_ps_nl[it]");
    Free(PP->phi_ps_nl[it], "RemovePseudoRead: PP->phi_ps_nl[it]");
    Free(PP->R[it], "RemovePseudoRead: PP->R[it]");
    Free(PP->v_loc[it], "RemovePseudoRead: PP->v_loc[it]");
    Free(PP->r[it], "RemovePseudoRead: PP->r[it]");
    Free(PP->integrand2[it], "RemovePseudoRead: PP->integrand2[it]");
    Free(PP->rcl[it], "RemovePseudoRead: PP->rcl[it]");
    Free(PP->al[it], "RemovePseudoRead: PP->al[it]");
    Free(PP->bl[it], "RemovePseudoRead: PP->bl[it]");
    Free(PP->FacGauss[it], "RemovePseudoRead: PP->FacGauss[it]");
    Free(PP->phi_ps_loc[it], "RemovePseudoRead: PP->phi_ps_loc[it]");
    Free(PP->core[it], "RemovePseudoRead: PP->core[it]");
    Free(PP->rc[it], "RemovePseudoRead: PP->rc[it]");
    Free(PP->wNonLoc[it], "RemovePseudoRead: PP->wNonLoc[it]");
  }
  //Free(PP->phi_ps_nl, "RemovePseudoRead: PP->phi_ps_nl");
  Free(PP->phi_ps_nl, "RemovePseudoRead: PP->phi_ps_nl");
  Free(PP->R, "RemovePseudoRead: PP->R");
  Free(PP->v_loc, "RemovePseudoRead: PP->v_loc");
  Free(PP->r, "RemovePseudoRead: PP->r");
  Free(PP->integrand2, "RemovePseudoRead: PP->integrand2");
  Free(PP->rcl, "RemovePseudoRead: PP->rcl");
  Free(PP->al, "RemovePseudoRead: PP->al");
  Free(PP->bl, "RemovePseudoRead: PP->bl");
  Free(PP->FacGauss, "RemovePseudoRead: PP->FacGauss");
  Free(PP->phi_ps_loc, "RemovePseudoRead: PP->phi_ps_loc");
  Free(PP->core, "RemovePseudoRead: PP->core");
  Free(PP->rc, "RemovePseudoRead: PP->rc");
  Free(PP->wNonLoc, "RemovePseudoRead: PP->wNonLoc");
  Free(PP->ei, "RemovePseudoRead: PP->ei");
  Free(PP->expiGR, "RemovePseudoRead: PP->expiGR");
  Free(PP->corewave, "RemovePseudoRead: PP->corewave");
  Free(PP->formfCore, "RemovePseudoRead: PP->formfCore");
  Free(PP->fnl, "RemovePseudoRead: PP->fnl"); 
  Free(PP->VCoulombc, "RemovePseudoRead: PP->VCoulombc");
  Free(PP->lm_end, "RemovePseudoRead: PP->lm_end");
  Free(PP->zval, "RemovePseudoRead: PP->zval");
  Free(PP->nang, "RemovePseudoRead: PP->nang");
  Free(PP->mmax, "RemovePseudoRead: PP->mmax");
  Free(PP->clog, "RemovePseudoRead: PP->clog");
}

/* Calculate Energy and Forces */

/** Calculates Gauss (self) energy term of ions \f$E_{self}\f$.
 * \f[
 *    E_{self} = \sum_{I_s,I_a} \frac{1}{\sqrt{2\pi}} \frac{Z_{I_s}^2}{r_{I_s}^{Gauss}} \qquad (4.10 last part)
 * \f]
 * \param *P Problem at hand
 * \param *PP PseudoPot ential structure
 * \param *I Ions structure
 */
void CalculateGaussSelfEnergyNoRT(struct Problem *P, struct PseudoPot *PP, struct Ions *I) 
{
  double SumE = 0.0;
  int it;
  for (it=0; it < I->Max_Types; it++) {
    if (I->I[it].rgauss < MYEPSILON) fprintf(stderr,"CalculateGaussSelfEnergyNoRT: I->I[it].rgauss = %lg\n",I->I[it].rgauss);
    SumE += PP->zval[it]*PP->zval[it]/I->I[it].rgauss*I->I[it].Max_IonsOfType;
  }
  P->Lat.E->AllTotalIonsEnergy[GaussSelfEnergy] = SumE*PP->fac1sqrt2PI; 
}

/** Calculates non local pseudopotential energy \f$E_{ps,nl,i}\f$.
 * First, calculates non-local form factors 
 * \f[
 *    f^{nl}_{i,I_s,I_a,l,m} = \sum_G \exp(-i(G R_{I_s,I_a}) 
 *      \phi_{I_s,l,m}^{ps,nl} (G) c_{i,G}
 * \f]
 * and afterwards evaluates with these
 * \f[
 *    E_{ps,nl,i} = f_{i} \sum_{I_s,I_a} \sum_{l,m} w^{nl}_{I_s,l} | f^{nl}_{i,I_s,I_a,l,m} |^2  \qquad (4.16)
 * \f]
 * \param *P Problem at hand
 * \param i Energy of i-th Psi
 */
void CalculateNonLocalEnergyNoRT(struct Problem *P, const int i)
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS; 
  int it,iot,ilm,ig,s=0;
  double sf,enl;//,enlm;
  fftw_complex PPdexpiGRpkg;
  fftw_complex *dfnl = PP->dfnl;
  fftw_complex *LocalPsi = LevS->LPsi->LocalPsi[i];   // i-th wave function coefficients
  const double PsiFactor = Lat->Psi.LocalPsiStatus[i].PsiFactor;
  int MinusGFactor;
  enl=0.0;
  for (it=0; it < I->Max_Types; it++) { // through all ion types
    for (ilm=1; ilm <=PP->lm_end[it]; ilm++) {  // over all possible m values
      MinusGFactor = Get_pkga_MinusG(ilm-1);
      SetArrayToDouble0((double *)dfnl, 2*I->I[it].Max_IonsOfType);
      for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {   // for each ion per type
    		s=0;
    		if (LevS->GArray[0].GSq == 0.0) {
    		  dexpiGRpkg(PP, it, iot, ilm, 0, &PPdexpiGRpkg);
    		  c_re(dfnl[iot]) += (c_re(LocalPsi[0])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[0])*c_im(PPdexpiGRpkg));
    		  c_im(dfnl[iot]) += (c_re(LocalPsi[0])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[0])*c_re(PPdexpiGRpkg));
    		  s++;
    		}
    		for (ig=s; ig < LevS->MaxG; ig++) {
    		  dexpiGRpkg(PP, it, iot, ilm, ig, &PPdexpiGRpkg);
    		  c_re(dfnl[iot]) += (c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg));
    		  c_im(dfnl[iot]) += (c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg));
    		  /* Minus G - due to inherent symmetry at gamma point! */
    		  c_re(dfnl[iot]) += MinusGFactor*(c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg));
    		  c_im(dfnl[iot]) -= MinusGFactor*(c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg));
    		}
      }
      // Allreduce as the coefficients are spread over many processes
      MPI_Allreduce( dfnl, PP->fnl[i][it][ilm-1], 2*I->I[it].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
    }
  }
  for (it=0; it < I->Max_Types; it++) {
    for (ilm=1; ilm <=PP->lm_end[it]; ilm++) {
      for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
    		//enlm = 0.0;
    		sf = c_re(PP->fnl[i][it][ilm-1][iot])*c_re(PP->fnl[i][it][ilm-1][iot])+c_im(PP->fnl[i][it][ilm-1][iot])*c_im(PP->fnl[i][it][ilm-1][iot]);
    		//enlm = (PsiFactor*sf);
    		enl += PsiFactor*sf*PP->wNonLoc[it][ilm-1]; 
      }
    }
  } 
  Lat->Energy[R->CurrentMin].PsiEnergy[NonLocalEnergy][i] = enl;
}


/** Calculates non local pseudopotential energy \f$E_{ps,nl,i}\f$ using Riemann tensor.
 * \param *P Problem at hand
 * \param i value
 * \note not implemented
 */
void CalculateNonLocalEnergyUseRT(struct Problem *P, const int i)
{
  Error(SomeError, "CalculateNonLocalEnergyUseRT: Not implemented");
}

/* AddVICGP aus scp */ /* HG wird geloescht und neu gesetzt */

/** Calculates Gauss, Pseudopotential, Hartreepotential and Hartree energies, also
 * PseudoPot::VCoulombc \f$V^H (G)\f$ and Density::DensityCArray[DoubleDensityTypes#HGDensity].
 * \f[
 *    E_{Gauss} = V \sum_{G\neq 0} \frac{4\pi}{|G|^2} |\overbrace{\sum_{I_s} S_{I_s} (G) \phi^{Gauss}_{I_s} (G)}^{\hat{n}^{Gauss}(G)}|^2 \qquad (\textnormal{just before 4.22})
 * \f]
 * \f[
 *    \widetilde{E}_{ps,loc} = 2V \sum_G \sum_{I_s} S_{I_s} (G) \phi^{ps,loc}_{I_s} (G) \hat{n}^\ast (G) \quad (4.13)
 * \f]
 * \f[
 *    E_{Hpot} = V \sum_{G\neq 0} \frac{4\pi}{|G|^2} |\hat{n}(G)+\hat{n}^{Gauss}(G)|^2  \qquad (\textnormal{first part 4.10})
 * \f]
 * \f[
 *    E_{Hartree} = V \sum_{G\neq 0} \frac{4\pi}{|G|^2} |\hat{n}(G)|^2
 * \f]
 * \f[
 *    V^H (G) = \frac{4\pi}{|G|^2} (\hat{n}(G)+\hat{n}^{Gauss}(G))  \qquad (\textnormal{section 4.3.2})
 * \f]
 * \param *P Problem at hand
 * \note There are some factor 2 discrepancies to formulas due to gamma point symmetry!
 */
void CalculateSomeEnergyAndHGNoRT(struct Problem *P) 
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct Lattice *Lat = &P->Lat;
  struct Energy *E = Lat->E;
  fftw_complex vp,rp,rhog,rhoe;
  double SumEHP,Fac,SumEH,SumEG,SumEPS;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct Density *Dens = Lev0->Dens;
  int g,s=0,it,Index,i;
  fftw_complex *HG = Dens->DensityCArray[HGDensity];
  SetArrayToDouble0((double *)HG,2*Dens->TotalSize);
  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]);
    }
    c_re(HG[Index]) = c_re(vp);
    SumEPS += (c_re(Dens->DensityCArray[TotalDensity][Index])*c_re(vp) +
         c_im(Dens->DensityCArray[TotalDensity][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 = 4.*PI/(Lev0->GArray[g].GSq);
    for (it = 0; it < I->Max_Types; it++) {
      c_re(vp) += (c_re(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);  // V^ps,loc
      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]);    // n^gauss
      c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
    } // rp = n^{Gauss)(G)
    c_re(rhog) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp);  // n + n^gauss = n^~
    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;
    // rhog = n(G) + n^{Gauss}(G), rhoe = n(G)
    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(rhog);
    c_im(HG[Index]) = c_im(vp)+Fac*c_im(rhog);
    /*if (P->first) */
    SumEG += Fac*(c_re(rp)*c_re(rp)+c_im(rp)*c_im(rp));             // Gauss energy
    SumEHP += Fac*(c_re(rhog)*c_re(rhog)+c_im(rhog)*c_im(rhog));    // E_ES first part
    SumEH += Fac*(c_re(rhoe)*c_re(rhoe)+c_im(rhoe)*c_im(rhoe));     // Hartree energy
    SumEPS += 2.*(c_re(Dens->DensityCArray[TotalDensity][Index])*c_re(vp)+
		  c_im(Dens->DensityCArray[TotalDensity][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;
  E->AllLocalDensityEnergy[HartreePotentialEnergy] = SumEHP*Lat->Volume;
  E->AllLocalDensityEnergy[HartreeEnergy] = SumEH*Lat->Volume;
  /* ReduceAllEnergy danach noch aufrufen */
}

/** Calculates \f$f^{nl}_{i,I_s,I_a,l,m}\f$ for conjugate direction Psis.
 * \param *P Problem at hand
 * \param *ConDir array of complex coefficients of the conjugate direction
 * \param ***CDfnl return array [ion type, lm value, ion of type]
 * \sa CalculateConDirHConDir() - used there
 */
void CalculateCDfnl(struct Problem *P, fftw_complex *ConDir, fftw_complex ***CDfnl)
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct RunStruct *R = &P->R;
  const struct LatticeLevel *LevS = R->LevS;
  int it,iot,ilm,ig,MinusGFactor,s;
  fftw_complex PPdexpiGRpkg;
  fftw_complex *dfnl = PP->dfnl;
  for (it=0; it < I->Max_Types; it++) {
    for (ilm=1; ilm <=PP->lm_end[it]; ilm++) {
      MinusGFactor = Get_pkga_MinusG(ilm-1);
      SetArrayToDouble0((double *)dfnl, 2*I->I[it].Max_IonsOfType);
      for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
      	s=0;
      	if (LevS->GArray[0].GSq == 0.0) {
      	  dexpiGRpkg(PP, it, iot, ilm, 0, &PPdexpiGRpkg);
      	  c_re(dfnl[iot]) += (c_re(ConDir[0])*c_re(PPdexpiGRpkg)-c_im(ConDir[0])*c_im(PPdexpiGRpkg));
      	  c_im(dfnl[iot]) += (c_re(ConDir[0])*c_im(PPdexpiGRpkg)+c_im(ConDir[0])*c_re(PPdexpiGRpkg));
      	  s++;
      	}
      	for (ig=s; ig < LevS->MaxG; ig++) {
      	  dexpiGRpkg(PP, it, iot, ilm, ig, &PPdexpiGRpkg);
      	  c_re(dfnl[iot]) += (c_re(ConDir[ig])*c_re(PPdexpiGRpkg)-c_im(ConDir[ig])*c_im(PPdexpiGRpkg));
      	  c_im(dfnl[iot]) += (c_re(ConDir[ig])*c_im(PPdexpiGRpkg)+c_im(ConDir[ig])*c_re(PPdexpiGRpkg));
      	  /* Minus G */
      	  c_re(dfnl[iot]) += MinusGFactor*(c_re(ConDir[ig])*c_re(PPdexpiGRpkg)-c_im(ConDir[ig])*c_im(PPdexpiGRpkg));
      	  c_im(dfnl[iot]) -= MinusGFactor*(c_re(ConDir[ig])*c_im(PPdexpiGRpkg)+c_im(ConDir[ig])*c_re(PPdexpiGRpkg));
      	}
      }
      MPI_Allreduce( dfnl, CDfnl[it][ilm-1], 2*I->I[it].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
    }
  }
}

/** Return non-local potential \f$V^{ps,nl}(G)|\Psi_{i} \rangle\f$.
 * \a *HNL is set to zero and the non-local potential added.
 * \param *P Problem at hand
 * \param *HNL return array of real coefficients
 * \param ***fnl non-local form factors for specific wave function, see PseudoPot::fnl
 * \param PsiFactor occupation number of respective wave function
 * \sa CalculcateGradientNoRT() - used there in gradient calculation
 *     CalculateConDirHConDir() - used there with conjugate directions (different non-local form factors!)
 */
void CalculateAddNLPot(struct Problem *P, fftw_complex *HNL, fftw_complex ***fnl, double PsiFactor) 
{
  struct PseudoPot *PP = &P->PP;
  struct Ions *I = &P->Ion;
  struct RunStruct *R = &P->R;
  const struct LatticeLevel *LevS = R->LevS;
  int it,iot,ig,ilm;
  fftw_complex dexpiGR, dpkg, tt1, tt2, tt3, tt4, tt;
  double rr1=0, rr2=0, rr3=0, rr4=0;
  double *rr = PP->rr;
  fftw_complex *t = PP->t;
  SetArrayToDouble0((double *)HNL,2*LevS->MaxG);
  for (it=0; it < I->Max_Types; it++) {   //over all types
    if (PP->lm_end[it] == 4) {
      rr1 = PP->wNonLoc[it][0];
      rr2 = PP->wNonLoc[it][1];
      rr3 = PP->wNonLoc[it][2];
      rr4 = PP->wNonLoc[it][3];
    } else {
      for (ilm=0; ilm < PP->lm_end[it]; ilm++) {
      	rr[ilm] = PP->wNonLoc[it][ilm];
      }
    }
    for (iot =0; iot < I->I[it].Max_IonsOfType; iot++) {  // over all ions of this type
      if (PP->lm_end[it] == 4) {  // over all lm-values
      	c_re(tt1) = -c_re(fnl[it][0][iot])*rr1;
      	c_im(tt1) = -c_im(fnl[it][0][iot])*rr1;
      	c_re(tt2) = -c_re(fnl[it][1][iot])*rr2;
      	c_im(tt2) = -c_im(fnl[it][1][iot])*rr2;
      	c_re(tt3) = -c_re(fnl[it][2][iot])*rr3;
      	c_im(tt3) = -c_im(fnl[it][2][iot])*rr3;
      	c_re(tt4) = -c_re(fnl[it][3][iot])*rr4;	
      	c_im(tt4) = -c_im(fnl[it][3][iot])*rr4;	
      	for (ig=0; ig < LevS->MaxG; ig++) {
      	  expiGR(PP, it, iot, ig, &dexpiGR);
      	  c_re(dpkg) = PP->phi_ps_nl[it][0][ig]*c_re(tt1)+PP->phi_ps_nl[it][1][ig]*c_re(tt2)+PP->phi_ps_nl[it][2][ig]*c_re(tt3)+PP->phi_ps_nl[it][3][ig]*c_re(tt4); 
      	  c_im(dpkg) = PP->phi_ps_nl[it][0][ig]*c_im(tt1)+PP->phi_ps_nl[it][1][ig]*c_im(tt2)+PP->phi_ps_nl[it][2][ig]*c_im(tt3)+PP->phi_ps_nl[it][3][ig]*c_im(tt4); 
      	  c_re(HNL[ig]) -= (c_re(dexpiGR)*c_re(dpkg)-c_im(dexpiGR)*c_im(dpkg))*PsiFactor;
      	  c_im(HNL[ig]) -= (c_re(dexpiGR)*c_im(dpkg)+c_im(dexpiGR)*c_re(dpkg))*PsiFactor;
      	}
      } else {
    		for (ilm=0; ilm < PP->lm_end[it]; ilm++) { // over all lm-values
    		  c_re(t[ilm]) = -c_re(fnl[it][ilm][iot])*rr[ilm];
    		  c_im(t[ilm]) = -c_im(fnl[it][ilm][iot])*rr[ilm];
    		}
    		if (PP->lm_end[it] > 0) {
    		  for (ig=0; ig < LevS->MaxG; ig++) {
    		    c_re(tt) = c_re(t[0])*PP->phi_ps_nl[it][0][ig];
    		    c_im(tt) = c_im(t[0])*PP->phi_ps_nl[it][0][ig];	
    		    for (ilm=1; ilm < PP->lm_end[it]; ilm++) {
    		      c_re(tt) += c_re(t[ilm])*PP->phi_ps_nl[it][ilm][ig];
    		      c_im(tt) += c_im(t[ilm])*PP->phi_ps_nl[it][ilm][ig];
    		    }
    		    expiGR(PP,it,iot,ig,&dexpiGR);
    		    c_re(HNL[ig]) -= (c_re(dexpiGR)*c_re(tt)-c_im(dexpiGR)*c_im(tt))*PsiFactor;
    		    c_im(HNL[ig]) -= (c_re(dexpiGR)*c_im(tt)+c_im(dexpiGR)*c_re(tt))*PsiFactor;
    		  }
    		}
      }
    }
  }
  if (LevS->GArray[0].GSq == 0.0)
    HNL[0].im = 0.0;
}

/** Calculates the local force acting on each ion.
 * \param *P Problem at hand
 */
void CalculateIonLocalForce(struct Problem *P) 
{
  int is,ia,g2,Index,s,i;
  double *G;
  struct Lattice *L = &P->Lat;
  struct Ions *I = &P->Ion;
  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  fft_plan_3d *plan = L->plan; 
  fftw_complex *work = Dens->DensityCArray[TempDensity];
  fftw_complex *HGC = Dens->DensityCArray[HGDensity];
  fftw_real *HGCR = (fftw_real *)HGC;
  double ForceFac = L->Volume;
  double force, *dsum;
  fftw_complex cv;
  
  // precalc V_XC(r) and fft
  //debug(P,"precalc V_XC\n");
  SetArrayToDouble0((double *)HGCR, Dens->TotalSize*2);
  CalculateXCPotentialNoRT(P, HGCR);
  fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGC, work);
  //debug(P,"precalc V_XC finished\n");
  for (is=0; is < I->Max_Types; is++) {
    SetArrayToDouble0(I->FTemp, NDIM*I->I[is].Max_IonsOfType);
    for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) {
      dsum = &I->FTemp[ia*NDIM];
      s = 0;
      if (Lev0->GArray[0].GSq == 0.0)
    		s++;
      for (g2=s; g2 < Lev0->MaxG; g2++) {
    		Index = Lev0->GArray[g2].Index;
    		G = Lev0->GArray[g2].G;
    		c_re(cv) = c_re(PP->VCoulombc[g2])*PP->FacGauss[is][g2]+PP->phi_ps_loc[is][g2]*c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
    		c_im(cv) = c_im(PP->VCoulombc[g2])*PP->FacGauss[is][g2]-PP->phi_ps_loc[is][g2]*c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
        if (I->I[is].corecorr == CoreCorrected) {
          // this summand was missing before, is not thoroughly tested
          //debug(P,"V_XC summand in CalculateIonLocalForce()\n");
          c_re(cv) +=  (PP->formfCore[is][g2]*c_re(HGC[Index]))*R->HGcFactor;
          c_im(cv) += -(PP->formfCore[is][g2]*c_im(HGC[Index]))*R->HGcFactor;
          //fprintf(stderr, "(%i) FormfCore %lg\tHGC %lg+i%lg\n",P->Par.me, PP->formfCore[is][g2],c_re(HGC[Index]),c_im(HGC[Index]));
        }
    		force = c_re(PP->ei[is][ia][g2])*c_im(cv)+c_im(PP->ei[is][ia][g2])*c_re(cv);
    		dsum[0] += 2.*(-G[0]*force); //2.* ... why? Nowhere from derivative!
    		dsum[1] += 2.*(-G[1]*force); //2.*
    		dsum[2] += 2.*(-G[2]*force); //2.*
      }
    }
    MPI_Allreduce( I->FTemp, I->I[is].FIonL, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
    for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) 
      for (i=0; i<NDIM; i++)
				I->I[is].FIonL[i+NDIM*ia] *= ForceFac;
  }
}

/** Calculates the non-local force acting on ion.
 * \param *P Problem at hand
 */
void CalculateIonNonLocalForce(struct Problem *P) 
{
  double *G;
  double fnl[NDIM], AllLocalfnl[NDIM], AllUpfnl[NDIM], AllDownfnl[NDIM], AllTotalfnl[NDIM];
  double Fac;
  struct Lattice *L = &P->Lat;
  struct Psis *Psi = &L->Psi;
  struct Ions *I = &P->Ion;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  int MinusGFactor;
  fftw_complex PPdexpiGRpkg, sf_a, sf_s, *LocalPsi;
  int is,ia,ilm,ig,i,s,d;
  MPI_Status status;
  for (is=0; is < I->Max_Types; is++) {
    SetArrayToDouble0(I->I[is].FIonNL, NDIM*I->I[is].Max_IonsOfType);
    for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) {
      for (d=0; d<NDIM; d++)
		    fnl[d] = 0.0;
      for (i=0; i < L->Psi.LocalNo; i++) if (Psi->LocalPsiStatus[i].PsiType == P->R.CurrentMin) {
    		LocalPsi = LevS->LPsi->LocalPsi[i]; 
    		for (ilm=1; ilm <=PP->lm_end[is]; ilm++) {
    		  MinusGFactor = Get_pkga_MinusG(ilm-1);
    		  Fac = Psi->LocalPsiStatus[i].PsiFactor*PP->wNonLoc[is][ilm-1];
    		  c_re(sf_s) = c_re(PP->fnl[i][is][ilm-1][ia]);
    		  c_im(sf_s) = c_im(PP->fnl[i][is][ilm-1][ia]);
    		  s=0;
    		  if (LevS->GArray[0].GSq == 0.0) {
    		    s++;
    		  }
    		  for (ig=s; ig < LevS->MaxG; ig++) {
    		    dexpiGRpkg(PP, is, ia, ilm, ig, &PPdexpiGRpkg);
    		    G = LevS->GArray[ig].G;    
    		    for (d=0; d <NDIM; d++) {
    		      c_re(sf_a) = (c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]);
    		      c_im(sf_a) = (c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]);
    		      /* Minus G */
    		      c_re(sf_a) -= MinusGFactor*(c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]);
    		      c_im(sf_a) += MinusGFactor*(c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]);
    		      fnl[d] += Fac*(c_re(sf_a)*c_im(sf_s)-c_im(sf_a)*c_re(sf_s));
    		    }
    		  }
    		}
      } 
      MPI_Allreduce( fnl, AllLocalfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
      switch (Psi->PsiST) {
	      case SpinDouble:
    			MPI_Allreduce(AllLocalfnl, AllTotalfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); 
		      break;
	      case SpinUp:
    			MPI_Allreduce(AllLocalfnl, AllUpfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); 
    			MPI_Sendrecv(AllUpfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
			     AllDownfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
			     P->Par.comm_STInter, &status );
      		for (d=0; d< NDIM; d++)
      		  AllTotalfnl[d] = AllUpfnl[d] + AllDownfnl[d];
      		break;
	      case SpinDown:
    			MPI_Allreduce(AllLocalfnl, AllDownfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    			MPI_Sendrecv(AllDownfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
    			     AllUpfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
    			     P->Par.comm_STInter, &status );
      		for (d=0; d < NDIM; d++)
      		  AllTotalfnl[d] = AllUpfnl[d] + AllDownfnl[d];
      		break;
    		}
      for (d=0; d<NDIM;d++)
    		I->I[is].FIonNL[d+NDIM*ia] -= AllTotalfnl[d]*2.0;
    }
  }  
}