void CalculatePerturbedEnergy(struct Problem *P, int l, int DoGradient) 
{
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct Energy *E = Lat->E;
  struct PseudoPot *PP = &P->PP;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  int state = R->CurrentMin;
  int l_normal = Psi->TypeStartIndex[Occupied] + (l - Psi->TypeStartIndex[state]);  // offset l to \varphi_l^{(0)}
  int ActNum = l - Psi->TypeStartIndex[state] + Psi->TypeStartIndex[1] * Psi->LocalPsiStatus[l].my_color_comm_ST_Psi;
  int g, i, m, j;
  double lambda, Lambda;
  double RElambda10, RELambda10;
  double RElambda01, RELambda01;
  double IMlambda10, IMLambda10;
  double IMlambda01, IMLambda01;
  fftw_complex *source = LevS->LPsi->LocalPsi[l];
  fftw_complex *grad = P->Grad.GradientArray[ActualGradient];
  fftw_complex *Hc_grad = P->Grad.GradientArray[HcGradient];
  fftw_complex *H1c_grad = P->Grad.GradientArray[H1cGradient];
  fftw_complex *TempPsi_0 = H1c_grad;
  fftw_complex *TempPsi_1 = P->Grad.GradientArray[TempGradient];
  fftw_complex *varphi_1, *varphi_0;
  struct OnePsiElement *OnePsiB, *LOnePsiB;
  fftw_complex *LPsiDatB=NULL;
  int ElementSize = (sizeof(fftw_complex) / sizeof(double)), RecvSource;
  MPI_Status status;
  
  // ============ Calculate H^(0) psi^(1) =============================
  SetArrayToDouble0((double *)Hc_grad,2*R->InitLevS->MaxG);
  ApplyTotalHamiltonian(P,source,Hc_grad, PP->fnl[l], 1, 1);
    
  // ============ ENERGY FUNCTIONAL Evaluation  PART 1 ================  
  varphi_0 = LevS->LPsi->LocalPsi[l_normal];
  varphi_1 = LevS->LPsi->LocalPsi[l];
  switch (state) {
    case Perturbed_P0:
      CalculatePerturbationOperator_P(P,varphi_0,TempPsi_0,0,0); //  \nabla_0 | \varphi_l^{(0)} \rangle
      CalculatePerturbationOperator_P(P,varphi_1,TempPsi_1,0,0); //  \nabla_0 | \varphi_l^{(1)} \rangle
      break;
    case Perturbed_P1:
      CalculatePerturbationOperator_P(P,varphi_0,TempPsi_0,1,0); //  \nabla_1 | \varphi_l^{(0)} \rangle
      CalculatePerturbationOperator_P(P,varphi_1,TempPsi_1,1,0); //  \nabla_1 | \varphi_l^{(1)} \rangle
      break;
    case Perturbed_P2:
      CalculatePerturbationOperator_P(P,varphi_0,TempPsi_0,2,0); //  \nabla_1 | \varphi_l^{(0)} \rangle
      CalculatePerturbationOperator_P(P,varphi_1,TempPsi_1,2,0); //  \nabla_1 | \varphi_l^{(1)} \rangle
      break;
    case Perturbed_RxP0:
      CalculatePerturbationOperator_RxP(P,varphi_0,TempPsi_0,l,0,0); //  r \times \nabla | \varphi_l^{(0)} \rangle
      CalculatePerturbationOperator_RxP(P,varphi_1,TempPsi_1,l,0,0); //  r \times \nabla | \varphi_l^{(1)} \rangle
      //CalculatePerturbationOperator_R(P,varphi_0,TempPsi_0,0,l); //  r \times \nabla | \varphi_l^{(0)} \rangle
      //CalculatePerturbationOperator_R(P,varphi_1,TempPsi_1,0,l); //  r \times \nabla | \varphi_l^{(1)} \rangle
      break;
    case Perturbed_RxP1:
      CalculatePerturbationOperator_RxP(P,varphi_0,TempPsi_0,l,1,0); //  r \times \nabla | \varphi_l^{(0)} \rangle
      CalculatePerturbationOperator_RxP(P,varphi_1,TempPsi_1,l,1,0); //  r \times \nabla | \varphi_l^{(1)} \rangle
      //CalculatePerturbationOperator_R(P,varphi_0,TempPsi_0,1,l); //  r \times \nabla | \varphi_l^{(0)} \rangle
      //CalculatePerturbationOperator_R(P,varphi_1,TempPsi_1,1,l); //  r \times \nabla | \varphi_l^{(1)} \rangle
      break;
    case Perturbed_RxP2:
      CalculatePerturbationOperator_RxP(P,varphi_0,TempPsi_0,l,2,0); //  r \times \nabla | \varphi_l^{(0)} \rangle
      CalculatePerturbationOperator_RxP(P,varphi_1,TempPsi_1,l,2,0); //  r \times \nabla | \varphi_l^{(1)} \rangle
      //CalculatePerturbationOperator_R(P,varphi_0,TempPsi_0,2,l); //  r \times \nabla | \varphi_l^{(0)} \rangle
      //CalculatePerturbationOperator_R(P,varphi_1,TempPsi_1,2,l); //  r \times \nabla | \varphi_l^{(1)} \rangle
      break;
    default:
      fprintf(stderr,"(%i) CalculatePerturbedEnergy called whilst not within perturbation run: CurrentMin = %i !\n",P->Par.me, R->CurrentMin);
      break;
  }

  // ============ GRADIENT and EIGENVALUE Evaluation  Part 1==============
  lambda = 0.0;
  if ((DoGradient) && (grad != NULL)) {
    g = 0;
    if (LevS->GArray[0].GSq == 0.0) {
      lambda += Hc_grad[0].re*source[0].re;
      grad[0].re = -(Hc_grad[0].re + TempPsi_0[0].re);
      grad[0].im = -(Hc_grad[0].im + TempPsi_0[0].im);
      g++;
    }
    for (;g<LevS->MaxG;g++) {
      lambda += 2.*(Hc_grad[g].re*source[g].re + Hc_grad[g].im*source[g].im);
      grad[g].re = -(Hc_grad[g].re + TempPsi_0[g].re);
      grad[g].im = -(Hc_grad[g].im + TempPsi_0[g].im);
    }

    m = -1;
    for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) {  // go through all wave functions
      OnePsiB = &Psi->AllPsiStatus[j];    // grab OnePsiB
      if (OnePsiB->PsiType == state) {   // drop all but the ones of current min state
        m++;  // increase m if it is type-specific wave function
        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, PerturbedTag, P->Par.comm_ST_PsiT, &status );
          LPsiDatB=LevS->LPsi->TempPsi;
        } else {                  // .. otherwise send it to all other processes (Max_me... - 1)
          for (i=0;i<P->Par.Max_me_comm_ST_PsiT;i++)
            if (i != OnePsiB->my_color_comm_ST_Psi)
              MPI_Send( LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, i, PerturbedTag, 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
  
        g = 0;
        if (LevS->GArray[0].GSq == 0.0) { // perform the summation
          grad[0].re += Lat->Psi.lambda[ActNum][m]*LPsiDatB[0].re;
          grad[0].im += Lat->Psi.lambda[ActNum][m]*LPsiDatB[0].im;
          g++;
        }
        for (;g<LevS->MaxG;g++) {
          grad[g].re += Lat->Psi.lambda[ActNum][m]*LPsiDatB[g].re;
          grad[g].im += Lat->Psi.lambda[ActNum][m]*LPsiDatB[g].im;
        }
      }
    }
  } else {
    g = 0;
    if (LevS->GArray[0].GSq == 0.0) {
      lambda += Hc_grad[0].re*source[0].re;
      g++;
    }
    for (;g<LevS->MaxG;g++)
      lambda += 2.*(Hc_grad[g].re*source[g].re + Hc_grad[g].im*source[g].im);
  }
  MPI_Allreduce ( &lambda, &Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  //fprintf(stderr,"(%i) Lambda[%i] = %lg\n",P->Par.me, l, Lambda);
  Lat->Psi.AddData[l].Lambda = Lambda; 

  // ============ ENERGY FUNCTIONAL Evaluation  PART 2 ================
  // varphi_1 jas negative symmetry, returning TemPsi_0 from CalculatePerturbedOperator also, thus real part of scalar product
  // "-" due to purely imaginary wave function is on left hand side, thus becomes complex conjugated: i -> -i
  // (-i goes into pert. op., "-" remains when on right hand side)  
  IMlambda10 =  GradImSP(P,LevS,varphi_1,TempPsi_0) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
  IMlambda01 = -GradImSP(P,LevS,varphi_0,TempPsi_1) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
  RElambda10 =  GradSP(P,LevS,varphi_1,TempPsi_0) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
  RElambda01 = -GradSP(P,LevS,varphi_0,TempPsi_1) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
  
  MPI_Allreduce ( &RElambda10, &RELambda10, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  MPI_Allreduce ( &RElambda01, &RELambda01, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  MPI_Allreduce ( &IMlambda10, &IMLambda10, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  MPI_Allreduce ( &IMlambda01, &IMLambda01, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
  
  E->PsiEnergy[Perturbed1_0Energy][l] = RELambda10;
  E->PsiEnergy[Perturbed0_1Energy][l] = RELambda01;
  if (P->Par.me == 0) {
    fprintf(stderr,"RE.Lambda10[%i] = %lg\t RE.Lambda01[%i] = %lg\n", l, RELambda10, l, RELambda01);
    fprintf(stderr,"IM.Lambda10[%i] = %lg\t IM.Lambda01[%i] = %lg\n", l, IMLambda10, l, IMLambda01);
  }
}

void CalculatePerturbationOperator_RxP(struct Problem *P, fftw_complex *source, fftw_complex *dest, int l, int index, int symmetry)
{
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  struct Lattice *Lat = &P->Lat;
  struct fft_plan_3d *plan = Lat->plan;
  struct Density *Dens0 = Lev0->Dens;
  fftw_complex *tempdestRC =  Dens0->DensityCArray[CurrentDensity0];
  fftw_real *tempdestR = (fftw_real *) tempdestRC;
  fftw_complex *tempdestRC2 =  Dens0->DensityCArray[CurrentDensity1];
  fftw_real *tempdestR2 = (fftw_real *) tempdestRC2;
  fftw_complex *work =  Dens0->DensityCArray[TempDensity];
  fftw_complex *PsiC = (fftw_complex *) Dens0->DensityCArray[ActualPsiDensity];
  fftw_real *PsiCR = (fftw_real *) PsiC;
  fftw_real *RealPhiR = (fftw_real *) Dens0->DensityArray[TempDensity];
  fftw_real *RealPhiR2 = (fftw_real *) Dens0->DensityArray[Temp2Density];
  fftw_complex *posfac, *destsnd, *destrcv, *destsnd2, *destrcv2;
  double x[NDIM], fac[NDIM], *WCentre;
  int n[NDIM], N0, n0, g, Index, pos, iS, i0;

  // init pointers and values
  int myPE = P->Par.me_comm_ST_Psi;
  int N[NDIM], NUp[NDIM];
  N[0] = LevS->Plan0.plan->N[0];
  N[1] = LevS->Plan0.plan->N[1];
  N[2] = LevS->Plan0.plan->N[2];
  NUp[0] = LevS->NUp[0]; 
  NUp[1] = LevS->NUp[1]; 
  NUp[2] = LevS->NUp[2]; 
  double FFTFactor = 1./LevS->MaxN; 
  N0 = LevS->Plan0.plan->local_nx;
  WCentre = Lat->Psi.AddData[l].WannierCentre;
  
  // blow up source coefficients
  SetArrayToDouble0((double *)tempdestRC ,Dens0->TotalSize*2);
  SetArrayToDouble0((double *)tempdestRC2,Dens0->TotalSize*2);
  SetArrayToDouble0((double *)RealPhiR ,Dens0->TotalSize*2);
  SetArrayToDouble0((double *)RealPhiR2,Dens0->TotalSize*2);
  SetArrayToDouble0((double *)PsiC,Dens0->TotalSize*2);
  for (g=0; g<LevS->MaxG; g++) {
    Index = LevS->GArray[g].Index;
    posfac = &LevS->PosFactorUp[LevS->MaxNUp*g];
    destrcv = &tempdestRC[LevS->MaxNUp*Index]; 
    destrcv2 = &tempdestRC2[LevS->MaxNUp*Index]; 
    for (pos=0; pos<LevS->MaxNUp; pos++) {
      destrcv [pos].re = LevS->GArray[g].G[cross(index,1)]*(( source[g].im)*posfac[pos].re-(-source[g].re)*posfac[pos].im);
      destrcv [pos].im = LevS->GArray[g].G[cross(index,1)]*(( source[g].im)*posfac[pos].im+(-source[g].re)*posfac[pos].re);
      destrcv2[pos].re = LevS->GArray[g].G[cross(index,3)]*(( source[g].im)*posfac[pos].re-(-source[g].re)*posfac[pos].im);
      destrcv2[pos].im = LevS->GArray[g].G[cross(index,3)]*(( source[g].im)*posfac[pos].im+(-source[g].re)*posfac[pos].re);
    }
  }
  for (g=0; g<LevS->MaxDoubleG; g++) {
    destsnd  = &tempdestRC [LevS->DoubleG[2*g]*LevS->MaxNUp];
    destrcv  = &tempdestRC [LevS->DoubleG[2*g+1]*LevS->MaxNUp];
    destsnd2 = &tempdestRC2[LevS->DoubleG[2*g]*LevS->MaxNUp];
    destrcv2 = &tempdestRC2[LevS->DoubleG[2*g+1]*LevS->MaxNUp];
    for (pos=0; pos<LevS->MaxNUp; pos++) {
      destrcv [pos].re =  destsnd [pos].re;
      destrcv [pos].im = -destsnd [pos].im;
      destrcv2[pos].re =  destsnd2[pos].re;
      destrcv2[pos].im = -destsnd2[pos].im;
    }    
  }  
  
  // fourier transform blown up wave function
  fft_3d_complex_to_real(plan,LevS->LevelNo, FFTNFUp, tempdestRC , work);
  DensityRTransformPos(LevS,tempdestR ,RealPhiR );
  fft_3d_complex_to_real(plan,LevS->LevelNo, FFTNFUp, tempdestRC2, work);
  DensityRTransformPos(LevS,tempdestR2,RealPhiR2);

  //fft_Psi(P,source,RealPhiR ,cross(index,1),6);
  //fft_Psi(P,source,RealPhiR2,cross(index,3),6);
  
  // for every point on the real grid multiply with component of position vector
  for (n0=0; n0<N0; n0++)
    for (n[1]=0; n[1]<N[1]; n[1]++)
      for (n[2]=0; n[2]<N[2]; n[2]++) {
        n[0] = n0 + N0 * myPE;
        fac[0] = (double)(n[0])/(double)((N[0]));
        fac[1] = (double)(n[1])/(double)((N[1]));
        fac[2] = (double)(n[2])/(double)((N[2]));
        RMat33Vec3(x,Lat->RealBasis,fac);
        iS = n[2] + N[2]*(n[1] + N[1]*n0);  // mind splitting of x axis due to multiple processes
        i0 = n[2]*NUp[2]+N[2]*NUp[2]*(n[1]*NUp[1]+N[1]*NUp[1]*n0*NUp[0]);
        //PsiCR[iS] = ((double)n[cross(index,0)]/(double)N[cross(index,0)]*sqrt(Lat->RealBasisSQ[cross(index,0)]) - WCentre[cross(index,0)])*RealPhiR[i0] - ((double)n[cross(index,2)]/(double)N[cross(index,2)]*sqrt(Lat->RealBasisQ[cross(index,2)]) - WCentre[cross(index,2)])*RealPhiR2[i0];  
        PsiCR[iS] = 
            MinImageConv(Lat,x[cross(index,0)],WCentre[cross(index,0)],cross(index,0)) * RealPhiR [i0] 
          - MinImageConv(Lat,x[cross(index,2)],WCentre[cross(index,2)],cross(index,2)) * RealPhiR2[i0]; 
        //tmp += truedist(Lat,x[index_r],WCentre[index_r],index_r) * RealPhiR[i0];
        //tmp += sawtooth(P,truedist(Lat,x[index_r],WCentre[index_r],index_r), index_r)*RealPhiR[i0];
        //(Fehler mit falschem Ort ist vor dieser Stelle!): ueber result =  RealPhiR[i0] * (x[index_r]) * RealPhiR[i0]; gecheckt
      }
  
  // inverse fourier transform 
  fft_3d_real_to_complex(plan,LevS->LevelNo, FFTNF1, PsiC, work);
  
  // copy to destination array
  for (g=0; g<LevS->MaxG; g++) {
    Index = LevS->GArray[g].Index;
    dest[g].re = ( PsiC[Index].re)*FFTFactor; 
    dest[g].im = ( PsiC[Index].im)*FFTFactor;
    if (LevS->GArray[g].GSq == 0)
      dest[g].im = 0; // imaginary of G=0 is zero
  }
}