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

/** \file density.c
 * Density Calculation.
 * 
 * The basic routines are CalculateOneDensityR() and CalculateOneDensityC(), which do the
 * level up from wave function coefficients in the current level to real density in the upper
 * and from density in the current (upper) to complex wave function coefficients.
 * Initialization of necessary structure such as Density is done in InitDensity(), where in
 * InitDensityCalculation() the density arrays are allocated, zeroed, and the first density
 * calculated and stored. The Density is updated after every minimisation step via 
 * UpdateDensityCalculation.
 * Small routines such as DensityRTransformPos() - performing coefficient permutation before
 * the Levelup fft - CalculateNativeIntDens() - calculating \a naively the density and
 * ControlNativeDensity - checking current electronic density against number of electrons -
 * help in the calculation of these and keeping them valid.
 * One last important function is ChangePsiAndDensToLevUp(), which performs similarly to
 * CalculateOneDensityR(), implementing afterwards also CalculateOneDensityC. It basically
 * calculates the density in the upper level from the coefficients and there retrieves the
 * coefficients by back transformation, thus doing the LevelUp for the Psis and the Density.
 * 
 * 
  Project: ParallelCarParrinello
  \author Jan Hamaekers
  \date 2000

  File: density.c
  $Id: density.c,v 1.70 2007-03-30 15:29:18 foo Exp $
*/

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

#include"data.h"
#include"errors.h"
#include"helpers.h"
#include"density.h"
#include"myfft.h"
#include"mymath.h"
#include"density.h"
#include"perturbed.h"


/** Marks one array as in use and not to be overwritten.
 * \param *Dens pointer to Density structure
 * \param DensityType DensityType
 * \param re_im 0 - Density#DensityArray, 1 - Density#DensityCArray
 */ 
void LockDensityArray(struct Density *Dens, enum UseType DensityType, enum complex re_im) {
  enum UseType *marker;
  if (DensityType >=0 && DensityType < MaxDensityTypes) {
    if (re_im == real)
      marker = &Dens->LockArray[DensityType];
    else
      marker = &Dens->LockCArray[DensityType];
    if (*marker == InUse) {
      fprintf(stderr,"Regarding %i array of type %i (0=real, 1=imag): %i\n",DensityType,re_im, *marker);
      Error(SomeError,"Array is already locked!");
    } else {
      *marker = InUse;
      //fprintf(stderr,"Locking %i (0=real, 1=imag) array of type %i ... \n",re_im,DensityType);
    }
  } else {
    fprintf(stderr,"LockDensityArray(): ILLEGAL Type %i ... \n",DensityType);
  }
}

/** Marks one array as not in use.
 * \param *Dens pointer to Density structure
 * \param DensityType DensityType
 * \param re_im 0 - Density#DensityArray, 1 - Density#DensityCArray
 */ 
void UnLockDensityArray(struct Density *Dens, enum UseType DensityType, enum complex re_im) {
  enum UseType *marker;
  if (DensityType >=0 && DensityType < MaxDensityTypes) {
    if (re_im == real)
      marker = &Dens->LockArray[DensityType];
    else
      marker = &Dens->LockCArray[DensityType];
    if (*marker == NotInUse) {
      fprintf(stderr,"Regarding %i array of type %i (0=real, 1=imag): %i\n",DensityType,re_im, *marker);
      Error(SomeError,"Array was not locked!");
    } else {
      *marker = NotInUse;
      //fprintf(stderr,"UnLocking %i (0=real, 1=imag) array of type %i ... \n",re_im,DensityType);
    }
  } else {
    fprintf(stderr,"LockDensityArray(): ILLEGAL Type %i ... \n",DensityType);
  }
}

/** Calculation of naive initial density.
 * \a *densR is summed up over all points on the grid (divided by RiemannTensor::DensityR if used),
 * the sum times \a factor is summed up among all processes and divived over the number of points on
 * the mesh LatticeLevel::MaxN returned.
 * \param *P Problem at hand, contains Lattice::RiemannTensor
 * \param *Lev LatticeLevel structure
 * \param *densR \f$n_i\f$, real density array for every point i on the mesh
 * \param factor additional factor (such as volume element)
 * \return \f$\frac{\sum^N_i n_i}{N}\f$
 */
double CalculateNativeIntDens(const struct Problem *P, const struct LatticeLevel *Lev, fftw_real *densR, const double factor) 
{
  int i;
  double Sum = 0.0, res;
  fftw_real *RTFactor = NULL;
  switch (P->Lat.RT.ActualUse) {
  case inactive:
  case standby:
    for (i=0; i < Lev->Dens->LocalSizeR; i++) 
      Sum  += densR[i];
    break;
  case active:
    RTFactor = P->Lat.RT.DensityR[RTADetPreRT];
    for (i=0; i < Lev->Dens->LocalSizeR; i++) 
      Sum  += densR[i]/fabs(RTFactor[i]);
    break;
  }
  Sum *= factor;
  MPI_Allreduce(&Sum , &res, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi );
  return(res/Lev->MaxN);
}

/** Permutation of wave function coefficients before LevelUp FFT.
 * \f[
 *    N^{up}_z (z + N_z (L_y + N^{up}_y (y + N_y ( L_x + N^{up}_x \cdot x))))
 *    \quad \rightarrow \quad
 *    N^{up}_z ( L_y + N^{up}_y ( L_x + N^{up}_x (z + N_z (y + N_y \cdot x))))
 * \f]
 * In the first part each mesh point (x,y,z) is a 2x2x2 cube, whereas in the last
 * case there is a 2x2x2 alignment of the cubes of all mesh points (x,y,z). (Here
 * assumed that LevelUp doubles the mesh's finesse)
 * \param *Lev LatticeLevel, contains LevelPlan::local_nx and LatticeLevel::N
 * \param *source source real coefficients array
 * \param *dest destination real coefficients array
 * \sa CalculateOneDensityR() - used there.
 */
void DensityRTransformPos(const struct LatticeLevel *Lev, fftw_real *source, fftw_real *dest)
{
  int es = Lev->NUp[2];
  unsigned int cpyes = sizeof(fftw_real)*es;
  int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp[0],Ny=Lev->NUp[1];
  int lx,ly,z,Lx,Ly;
  for(lx=0; lx < nx; lx++)
    for(Lx=0; Lx < Nx; Lx++)
      for(ly=0; ly < ny; ly++)
        for(Ly=0; Ly < Ny; Ly++)
          for(z=0; z < nz; z++) {
            memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
                    &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
                    cpyes);
          }
}

/*int GetOtherIndex(const struct fft_plan_3d *plan, const int Level, int Index)
{
  int N = Index;
  int x,yz,Y,Z,YZ,y,z;
  int Nx = plan->Levplan[Level].nx;
  int Ny = plan->Levplan[Level].local_ny ;
  int Nz = plan->Levplan[Level].nz;
  int Nyz = plan->Levplan[Level].ny*Nz;
  int NZ = plan->NLSR[2]/2+1;
  int YY;
  x = N%Nx;
  N /= Nx;
  z = N%Nz;
  y = N/Nz;
  Y = (y+plan->myPE*Ny)/plan->Levplan[Level].LevelFactor;
  YY = (y+plan->myPE*Ny)%plan->Levplan[Level].LevelFactor;
  Z = z/plan->Levplan[Level].LevelFactor;
  YZ = Z+NZ*Y;
  YZ = plan->pw[YZ].Index;
  Z = (YZ%NZ)*plan->Levplan[Level].LevelFactor+(z)%plan->Levplan[Level].LevelFactor;
  yz = Z+Nz*(((YZ/NZ)*plan->Levplan[Level].LevelFactor+YY));
  return(yz+Nyz*x);
}
*/

/** Calculates the real density of one wave function orbit.
 * In order to evaluate densities from the wave function coefficients there is a reoccurring pattern:
 * First of all, the (finer) cubic grid of coefficients must be built. Therefore we go through all reciprocal
 * grid vectors G up to LatticeLevel::MaxG (N), through all possible nodes up to LatticeLevel::MaxNUp (8) in the
 * upper level, multiply each complex coefficient \a source[G] with the pulled out LatticeLevel::PosFactorUp[pos]
 * and store in the temporal array Density::DensityArray[DensityTypes#TempDensity].
 * The symmetry condition \f$c_{i,G}^\ast = - c_{i,G}\f$ on the (z=0)-plane up to LatticeLevel::MaxDoubleG
 * is applied.
 * Then follows the complex to real fftransformation fft_3d_complex_to_real() and finally a permutation by
 * calling DensityRTransformPos().\n
 * The resulting density is stored in Density::DensityCArray[DensityTypes::ActualPsiDensity] if desired. In any
 * case the return array is \a *destR, each entry multiplied by \a Factor and in RiemannTensor case by an
 * additional \a RTFactor.
 * \param *Lat Lattice structure
 * \param *Lev LatticeLevel structure
 * \param *Dens Density structure, contains Density::DensityCArray and Density::DensityArray
 * \param *source complex wave function coefficients
 * \param *destR where the density is returned
 * \param Factor returned density is multiplied by this factor (e.g. occupation number)
 * \param StorePsi Whether the Psi should be stored in Density::DensityCArray[DoubleDensityTypes#ActualPsiDensity] or not
 */
void CalculateOneDensityR(const struct Lattice *Lat, const struct LatticeLevel *Lev, const struct Density *Dens, const fftw_complex *source, fftw_real *destR, const double Factor, int StorePsi) 
{
  int Index,i,pos;
  struct  fft_plan_3d *plan = Lat->plan; 
  fftw_complex *work = (fftw_complex *)Dens->DensityCArray[TempDensity];
  fftw_complex *tempdestRC = (fftw_complex *)Dens->DensityArray[TempDensity];
  fftw_complex *posfac, *destpos, *destRCS, *destRCD;
  fftw_real *RTFactor = NULL;
  fftw_real *Psic = (fftw_real *)Dens->DensityCArray[ActualPsiDensity];
  if (Lat->RT.ActualUse == active)
    RTFactor = Lat->RT.DensityR[RTADetPreRT];
  SetArrayToDouble0((double *)tempdestRC, Dens->TotalSize*2);   // set array to zero
  // the density in the upper level is created by filling the 2x2x2 block with the coefficient times the position factor
  for (i=0;i<Lev->MaxG;i++) {
    Index = Lev->GArray[i].Index;
    posfac = &Lev->PosFactorUp[Lev->MaxNUp*i];      // jump to the specific 2x2x2 block for the position factor (taken from upper level!)
    destpos = &tempdestRC[Lev->MaxNUp*Index];       // jump to the speficif 2x2x2 block for the destination values (resorted with Index)
    for (pos=0; pos < Lev->MaxNUp; pos++) {     // complex:  destpos = source * postfac
      destpos[pos].re = source[i].re*posfac[pos].re-source[i].im*posfac[pos].im;
      destpos[pos].im = source[i].re*posfac[pos].im+source[i].im*posfac[pos].re;
    }
  }
  // symmetry condition is applied to the (z=0)-plane
  for (i=0; i<Lev->MaxDoubleG; i++) {
    destRCS = &tempdestRC[Lev->DoubleG[2*i]*Lev->MaxNUp];
    destRCD = &tempdestRC[Lev->DoubleG[2*i+1]*Lev->MaxNUp];
    for (pos=0; pos < Lev->MaxNUp; pos++) {
      //if (isnan(destRCD[pos].re)) fprintf(stderr,"CalculateOneDensityR: destRCD.re NaN!\n");
      //if (isnan(destRCD[pos].im)) fprintf(stderr,"CalculateOneDensityR: destRCD.im NaN!\n");
      destRCD[pos].re =  destRCS[pos].re;
      destRCD[pos].im = -destRCS[pos].im;
    }
  } 
  fft_3d_complex_to_real(plan, Lev->LevelNo, FFTNFUp, tempdestRC, work);
  DensityRTransformPos(Lev,(fftw_real*)tempdestRC,destR);
  
  if (StorePsi)
    for (i=0; i < Dens->LocalSizeR; i++)  {
      Psic[i] = destR[i]; /* FIXME RT */
      //if (isnan(Psic[i])) fprintf(stderr,"CalculateOneDensityR: Psic NaN!\n");
    }

  switch (Lat->RT.ActualUse) {
  case inactive:
  case standby:
    for (i=0; i < Dens->LocalSizeR; i++) 
      destR[i]  *= destR[i]*Factor;
    break;
  case active:
    for (i=0; i < Dens->LocalSizeR; i++) 
      destR[i]  *= destR[i]*Factor*fabs(RTFactor[i]);
    break;
  }  
}

/** Calculates the real density of one wave function orbit with negative symmetry.
 * In order to evaluate densities from the wave function coefficients there is a reoccurring pattern:
 * First of all, the (finer) cubic grid of coefficients must be built. Therefore we go through all reciprocal
 * grid vectors G up to LatticeLevel::MaxG (N), through all possible nodes up to LatticeLevel::MaxNUp (8) in the
 * upper level, multiply each complex coefficient \a source[G] with the pulled out LatticeLevel::PosFactorUp[pos]
 * and store in the temporal array Density::DensityArray[DensityTypes#TempDensity].
 * The symmetry condition \f$c_{i,G}^\ast = - c_{i,G}\f$ on the (z=0)-plane up to LatticeLevel::MaxDoubleG
 * is applied.
 * Then follows the complex to real fftransformation fft_3d_complex_to_real() and finally a permutation by
 * calling DensityRTransformPos().\n
 * The resulting density is stored in Density::DensityCArray[DensityTypes::ActualPsiDensity] if desired. In any
 * case the return array is \a *destR, each entry multiplied by \a Factor and in RiemannTensor case by an
 * additional \a RTFactor.
 * \param *Lat Lattice structure
 * \param *Lev LatticeLevel structure
 * \param *Dens Density structure, contains Density::DensityCArray and Density::DensityArray
 * \param *source complex wave function coefficients
 * \param *destR where the density is returned
 * \param Factor returned density is multiplied by this factor (e.g. occupation number)
 * \param StorePsi Whether the Psi should be stored in Density::DensityCArray[DoubleDensityTypes#ActualPsiDensity] or not
 */
void CalculateOneDensityImR(const struct Lattice *Lat, const struct LatticeLevel *Lev, const struct Density *Dens, fftw_complex *source, fftw_real *destR, const double Factor, int StorePsi) 
{
  int Index,i,pos;
  struct  fft_plan_3d *plan = Lat->plan; 
  fftw_complex *work = (fftw_complex *)Dens->DensityCArray[TempDensity];
  fftw_complex *tempdestRC = (fftw_complex *)Dens->DensityArray[TempDensity];
  fftw_complex *posfac, *destpos, *destRCS, *destRCD;
  fftw_real *RTFactor = NULL;
  fftw_real *Psic = (fftw_real *)Dens->DensityCArray[ActualPsiDensity];
  if (Lat->RT.ActualUse == active)
    RTFactor = Lat->RT.DensityR[RTADetPreRT];
  SetArrayToDouble0((double *)tempdestRC, Dens->TotalSize*2);   // set array to zero
  // the density in the upper level is created by filling the 2x2x2 block with the coefficient times the position factor
  for (i=0;i<Lev->MaxG;i++) {
    Index = Lev->GArray[i].Index;
    posfac = &Lev->PosFactorUp[Lev->MaxNUp*i];      // jump to the specific 2x2x2 block for the position factor (taken from upper level!)
    destpos = &tempdestRC[Lev->MaxNUp*Index];       // jump to the speficif 2x2x2 block for the destination values (resorted with Index)
    for (pos=0; pos < Lev->MaxNUp; pos++) {     // complex:  destpos = source * postfac
      destpos[pos].re = source[i].re*posfac[pos].re-source[i].im*posfac[pos].im;
      destpos[pos].im = source[i].re*posfac[pos].im+source[i].im*posfac[pos].re;
    }
  }
  // symmetry condition is applied to the (z=0)-plane
  for (i=0; i<Lev->MaxDoubleG; i++) {
    destRCS = &tempdestRC[Lev->DoubleG[2*i]*Lev->MaxNUp];
    destRCD = &tempdestRC[Lev->DoubleG[2*i+1]*Lev->MaxNUp];
    for (pos=0; pos < Lev->MaxNUp; pos++) {
      //if (isnan(destRCD[pos].re)) fprintf(stderr,"CalculateOneDensityR: destRCD.re NaN!\n");
      //if (isnan(destRCD[pos].im)) fprintf(stderr,"CalculateOneDensityR: destRCD.im NaN!\n");
      destRCD[pos].re = -destRCS[pos].re;
      destRCD[pos].im =  destRCS[pos].im;
    }
  } 
  fft_3d_complex_to_real(plan, Lev->LevelNo, FFTNFUp, tempdestRC, work);
  DensityRTransformPos(Lev,(fftw_real*)tempdestRC,destR);
  
  if (StorePsi)
    for (i=0; i < Dens->LocalSizeR; i++)  {
      Psic[i] = destR[i]; /* FIXME RT */
      //if (isnan(Psic[i])) fprintf(stderr,"CalculateOneDensityR: Psic NaN!\n");
    }

  switch (Lat->RT.ActualUse) {
  case inactive:
  case standby:
    for (i=0; i < Dens->LocalSizeR; i++) 
      destR[i]  *= destR[i]*Factor;
    break;
  case active:
    for (i=0; i < Dens->LocalSizeR; i++) 
      destR[i]  *= destR[i]*Factor*fabs(RTFactor[i]);
    break;
  }  
}

/** Calculates the complex density for one wave function.
 * Copies given real coefficients from \a *srcR into \a *destC and performs the real to
 * complex LevelUp-Transformation onto it, where Density::DensityCArray[TempDensity] is used
 * as a temporal storage. The \a factor is later used on every density point on the grid.
 * \param *Lat Lattice structure
 * \param *Lev LatticeLevel structure
 * \param *Dens Density structure
 * \param *srcR one real wave function from which density is produced
 * \param *destC complex density array to be calculated
 * \param Factor
 */
void CalculateOneDensityC(const struct Lattice *Lat, const struct LatticeLevel *Lev, const struct Density *Dens, const fftw_real *srcR, fftw_complex *destC, const double Factor) 
{
  int i;
  struct LatticeLevel *LevUp = &Lat->Lev[Lev->LevelNo-1];
  struct  fft_plan_3d *plan = Lat->plan; 
  fftw_complex *work = (fftw_complex *)Dens->DensityCArray[TempDensity];
  fftw_real *destCR = (fftw_real *)destC;
  double factor = Factor/LevUp->MaxN;
  // copy coefficients into destC
  memcpy(destCR,srcR,sizeof(fftw_real)*Dens->LocalSizeR);
  // fftransform
  fft_3d_real_to_complex(plan, LevUp->LevelNo, FFTNF1, destC, work);
  // and multiply by given factor
  for (i=0; i<Dens->LocalSizeC; i++) {    
    destC[i].re *= factor;
    destC[i].im *= factor;
  }
}

/** Initialization of Density calculation.
 * Set Density::DensityArray's to zero.\n 
 * Now for the standard level LevS for each local Psi calculate the real density times volume times
 * occupation number (Psis::PsiFactor), build native initial density by adding up occupation numbers,
 * and sum up DensityTypes::TotalLocalDensity to total DensityTypes::ActualDensity.
 * The complex density is calculated from this (real) total local array, and depending on the SpinType
 * results are gathered from each process by summation: real and complex DensityTypes::TotalLocalDensity
 * and native initial density Psis::NIDensity.\n
 * For the topmost level CalculateNativeIntDens() is called and the result compared to Psis::NIDensity.
 * \param *P Problem at hand
 */
void InitDensityCalculation(struct Problem *P)
{
  int p,d,i, type;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LevelPsi *LPsi = LevS->LPsi;
  struct Density *Dens = Lev0->Dens;
  struct Psis *Psi = &Lat->Psi;
  double factorDR = R->FactorDensityR;
  double factorDC = R->FactorDensityC;
  double NIDensity=0., NIDensityUp=0., NIDensityDown=0.;
  fftw_real *Temp1R, *Temp2R, *Temp3R;
  fftw_complex *Temp1C, *Temp2C, *Temp3C;
  MPI_Status status;
  
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitDensityCalculation\n", P->Par.me);
  // set arrays to zero if allocated
  for (d=0; d < MaxInitDensityTypes; d++) {
    if (Dens->DensityCArray[d] != NULL)
      SetArrayToDouble0((double *)Dens->DensityCArray[d], Dens->TotalSize*2);
    if (Dens->DensityArray[d] != NULL)
      SetArrayToDouble0((double *)Dens->DensityArray[d], Dens->TotalSize*2);
  }
  
  // for each local Psi calculate the real density times volume times occupation number (PsiFactor), build naive initial density by adding up occupation number and sum up local density to total one
  for (p=Psi->LocalNo-1; p >= 0; p--) {
    //if (isnan(Dens->DensityArray[ActualDensity][0])) { fprintf(stderr,"(%i) WARNING in InitDensityCalculation(): ActualDensity_%i[%i] = NaN!\n", P->Par.me, p, 0); Error(SomeError, "NaN-Fehler!"); }
    //if (P->Call.out[StepLeaderOut]) fprintf(stderr,"NIDensity: PsiType[%d] = %d\n",p, Psi->LocalPsiStatus[p].PsiType);
    if (R->CurrentMin == Psi->LocalPsiStatus[p].PsiType) NIDensity += Psi->LocalPsiStatus[p].PsiFactor;
    if (R->CurrentMin == Psi->LocalPsiStatus[p].PsiType || Psi->LocalPsiStatus[p].PsiType == Occupied) {
      CalculateOneDensityR(Lat, LevS, Dens, LPsi->LocalPsi[p], Dens->DensityArray[ActualDensity], factorDR*Psi->LocalPsiStatus[p].PsiFactor, (p==Psi->TypeStartIndex[R->CurrentMin]));
      if (Psi->LocalPsiStatus[p].PsiType == Occupied)
        Temp1R = Dens->DensityArray[TotalLocalDensity];
      else
        Temp1R = Dens->DensityArray[GapLocalDensity];
      Temp2R = Dens->DensityArray[ActualDensity];
      for (i=0; i < Dens->LocalSizeR; i++)
        Temp1R[i] += Temp2R[i];
    }
  }
  
  CalculateOneDensityR(Lat, LevS, Dens, LPsi->LocalPsi[Psi->TypeStartIndex[R->CurrentMin]], Dens->DensityArray[ActualDensity], factorDR*Psi->LocalPsiStatus[Psi->TypeStartIndex[R->CurrentMin]].PsiFactor, (Psi->TypeStartIndex[R->CurrentMin]));
  /* CalculateOneDensityC(Lat, LevS, Dens, Dens->DensityArray[ActualDensity], Dens->DensityCArray[ActualDensity], factorDC); */
  // calculate complex density from the (real) total local array
  CalculateOneDensityC(Lat, LevS, Dens, Dens->DensityArray[TotalLocalDensity], Dens->DensityCArray[TotalLocalDensity], factorDC);
  CalculateOneDensityC(Lat, LevS, Dens, Dens->DensityArray[GapLocalDensity], Dens->DensityCArray[GapLocalDensity], factorDC);
  // depending on SpinType gather results: real and complex DensityTypes::TotalLocalDensity and native initial density
  switch (Psi->PsiST) {
  case SpinDouble:
    //fprintf(stderr,"(%i) InitDensity SpinDouble\n",P->Par.me);
    if (R->CurrentMin == Occupied) { // occupied states
      MPI_Allreduce(Dens->DensityArray[TotalLocalDensity] ,Dens->DensityArray[TotalDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[TotalLocalDensity] ,Dens->DensityCArray[TotalDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    } else {// unoccupied
      MPI_Allreduce(Dens->DensityArray[GapLocalDensity] ,Dens->DensityArray[GapDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[GapLocalDensity] ,Dens->DensityCArray[GapDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    }
    
    // naive density
    MPI_Allreduce(&NIDensity ,&Psi->NIDensity[R->CurrentMin] , 1  , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    break;
  case SpinUp:
    //fprintf(stderr,"(%i) InitDensity SpinUp\n",P->Par.me);
    if (R->CurrentMin == Occupied) { // occupied states
      MPI_Allreduce(Dens->DensityArray[TotalLocalDensity] ,Dens->DensityArray[TotalUpDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[TotalLocalDensity] ,Dens->DensityCArray[TotalUpDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
        MPI_Sendrecv( Dens->DensityCArray[TotalUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag1, 
          Dens->DensityCArray[TotalDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag2, 
          P->Par.comm_STInter, &status );
        MPI_Sendrecv( Dens->DensityArray[TotalUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag3, 
          Dens->DensityArray[TotalDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag4, 
          P->Par.comm_STInter, &status );
      Temp1R = Dens->DensityArray[TotalDensity];
      Temp2R = Dens->DensityArray[TotalUpDensity];
      Temp3R = Dens->DensityArray[TotalDownDensity];
      Temp1C = Dens->DensityCArray[TotalDensity];
      Temp2C = Dens->DensityCArray[TotalUpDensity];
      Temp3C = Dens->DensityCArray[TotalDownDensity];
    } else { // unoccupied or perturbed states
      MPI_Allreduce(Dens->DensityArray[GapLocalDensity] ,Dens->DensityArray[GapUpDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[GapLocalDensity] ,Dens->DensityCArray[GapUpDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
      MPI_Sendrecv( Dens->DensityCArray[GapUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag7, 
        Dens->DensityCArray[GapDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag8, 
        P->Par.comm_STInter, &status );
      MPI_Sendrecv( Dens->DensityArray[GapUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag9, 
        Dens->DensityArray[GapDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag0, 
        P->Par.comm_STInter, &status );
      Temp1R = Dens->DensityArray[GapDensity];
      Temp2R = Dens->DensityArray[GapUpDensity];
      Temp3R = Dens->DensityArray[GapDownDensity];
      Temp1C = Dens->DensityCArray[GapDensity];
      Temp2C = Dens->DensityCArray[GapUpDensity];
      Temp3C = Dens->DensityCArray[GapDownDensity];
    }

    for (i=0; i < Dens->LocalSizeR; i++)
      Temp1R[i] = Temp2R[i]+Temp3R[i];
    for (i=0; i < Dens->LocalSizeC; i++) {
      Temp1C[i].re = Temp2C[i].re+Temp3C[i].re;
      Temp1C[i].im = Temp2C[i].im+Temp3C[i].im;
    }

    // we are up, so have naive spinup density
    MPI_Allreduce(&NIDensity ,&Psi->NIDensityUp[R->CurrentMin] ,1 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Sendrecv( Psi->NIDensityUp, Extra, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag5, 
                  Psi->NIDensityDown, Extra, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag6, 
                  P->Par.comm_STInter, &status );
    for (type=Occupied;type<Extra;type++)
      Psi->NIDensity[type] = Psi->NIDensityUp[type] + Psi->NIDensityDown[type];
    break;
  case SpinDown:
    //fprintf(stderr,"(%i) InitDensity SpinDown\n",P->Par.me);
    if (R->CurrentMin == Occupied) { // occupied states
      MPI_Allreduce(Dens->DensityArray[TotalLocalDensity] ,Dens->DensityArray[TotalDownDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[TotalLocalDensity] ,Dens->DensityCArray[TotalDownDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
      MPI_Sendrecv( Dens->DensityCArray[TotalDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag2, 
        Dens->DensityCArray[TotalUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag1, 
        P->Par.comm_STInter, &status );
      MPI_Sendrecv( Dens->DensityArray[TotalDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag4, 
        Dens->DensityArray[TotalUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag3, 
        P->Par.comm_STInter, &status );

      Temp1R = Dens->DensityArray[TotalDensity];
      Temp2R = Dens->DensityArray[TotalUpDensity];
      Temp3R = Dens->DensityArray[TotalDownDensity];
      Temp1C = Dens->DensityCArray[TotalDensity];
      Temp2C = Dens->DensityCArray[TotalUpDensity];
      Temp3C = Dens->DensityCArray[TotalDownDensity];
    } else { // unoccupied or perturbed states
      MPI_Allreduce(Dens->DensityArray[GapLocalDensity] ,Dens->DensityArray[GapDownDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[GapLocalDensity] ,Dens->DensityCArray[GapDownDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
      MPI_Sendrecv( Dens->DensityCArray[GapDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag8, 
        Dens->DensityCArray[GapUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag7, 
        P->Par.comm_STInter, &status );
      MPI_Sendrecv( Dens->DensityArray[GapDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag0, 
        Dens->DensityArray[GapUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag9, 
        P->Par.comm_STInter, &status );

      Temp1R = Dens->DensityArray[GapDensity];
      Temp2R = Dens->DensityArray[GapUpDensity];
      Temp3R = Dens->DensityArray[GapDownDensity];
      Temp1C = Dens->DensityCArray[GapDensity];
      Temp2C = Dens->DensityCArray[GapUpDensity];
      Temp3C = Dens->DensityCArray[GapDownDensity];
    }

    for (i=0; i < Dens->LocalSizeR; i++)
      Temp1R[i] = Temp2R[i]+Temp3R[i];
    for (i=0; i < Dens->LocalSizeC; i++) {
      Temp1C[i].re = Temp2C[i].re+Temp3C[i].re;
      Temp1C[i].im = Temp2C[i].im+Temp3C[i].im;
    }

    // we are up, so have naive spinup density
    MPI_Allreduce(&NIDensity ,&Psi->NIDensityDown[R->CurrentMin] , 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);

    MPI_Sendrecv( Psi->NIDensityDown, Extra, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag6, 
      Psi->NIDensityUp, Extra, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag5, 
      P->Par.comm_STInter, &status );
    for (type=Occupied;type<Extra;type++)
      Psi->NIDensity[type] = Psi->NIDensityUp[type] + Psi->NIDensityDown[type];
    break;
  }
  switch (R->CurrentMin) {
    case Occupied:
      NIDensity = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[TotalDensity], 1./factorDR);
      if (Psi->PsiST != SpinDouble) {
        NIDensityUp = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[TotalUpDensity], 1./factorDR);
        NIDensityDown = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[TotalDownDensity], 1./factorDR);
      }
      break;
    default:
    case UnOccupied:
      NIDensity = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[GapDensity], 1./factorDR);
      if (Psi->PsiST != SpinDouble) {
        NIDensityUp = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[GapUpDensity], 1./factorDR);
        NIDensityDown = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[GapDownDensity], 1./factorDR);
      }
      break;
  }
  if (P->Call.out[LeaderOut] && (P->Par.me == 0)) {
    if (fabs(NIDensity - Psi->NIDensity[R->CurrentMin]) >= MYEPSILON) {
      fprintf(stderr, "NativeDensity(L(%i)[%i]) = %g != %g eps (%g >= %g)\n",R->LevSNo,R->CurrentMin,NIDensity,Psi->NIDensity[R->CurrentMin],fabs(NIDensity - Psi->NIDensity[R->CurrentMin]),MYEPSILON);
    } else {
      fprintf(stderr, "NativeDensity(L(%i)[%i]) = %g Ok !\n",R->LevSNo,R->CurrentMin,NIDensity);
    }
    if (Psi->PsiST != SpinDouble) {
      if (fabs(NIDensityUp - Psi->NIDensityUp[R->CurrentMin]) >= MYEPSILON) {
                                fprintf(stderr, "NativeDensityUp(L(%i)[%i]) = %g != %g eps (%g >= %g)\n",R->LevSNo,R->CurrentMin,NIDensityUp,Psi->NIDensityUp[R->CurrentMin],fabs(NIDensityUp - Psi->NIDensityUp[R->CurrentMin]),MYEPSILON);
      } else {
                                fprintf(stderr, "NativeDensityUp(L(%i)[%i]) = %g Ok !\n",R->LevSNo,R->CurrentMin,NIDensityUp);
      }
      if (fabs(NIDensityDown - Psi->NIDensityDown[R->CurrentMin]) >= MYEPSILON) {
                                fprintf(stderr, "NativeDensityDown(L(%i)[%i]) = %g != %g eps (%g >= %g)\n",R->LevSNo,R->CurrentMin,NIDensityDown,Psi->NIDensityDown[R->CurrentMin],fabs(NIDensityDown - Psi->NIDensityDown[R->CurrentMin]),MYEPSILON);
      } else {
                                fprintf(stderr, "NativeDensityDown(L(%i)[%i]) = %g Ok !\n",R->LevSNo,R->CurrentMin,NIDensityDown);
      }
    }
  }
}

/** Recalculate real and complex density after the minimisation of a wave function.
 * In Density::DensityArray DensityTypes#ActualDensity is subtracted from DensityTypes#TotalLocalDensity, 
 * recalculated by CalculateOneDensityR() for the last minimised Psi and the new value re-added.\n
 * Afterwards, CalculateOneDensityR() for current minimised Psi is called and stored
 * in Density::Types#ActualDensity.\n
 * And CalculateOneDensityC() called to evaluate
 * Density::DensityCArray[DensityTypes#TotalLocalDensity].\n
 * Finally, Density::DensityArray and Density::DensityCArray are summed up among the processes
 * into DensityTypes#TotalDensity.
 *
 * \param *P Problem at hand
 * \sa InitDensityCalculation() - very similar procedure.
 * \note ActualDensity is expected to be the old one before minimisation took place, as it is subtracted,
 * the density recalculated (for the newly minimsed wave function) and added again.
 */
void UpdateDensityCalculation(struct Problem *P)
{
  int i;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LevelPsi *LPsi = LevS->LPsi;
  struct Density *Dens = Lev0->Dens;
  struct Psis *Psi = &Lat->Psi;
  double factorDR = R->FactorDensityR;
  double factorDC = R->FactorDensityC;
  const int p = R->ActualLocalPsiNo;
  const int p_old = R->OldActualLocalPsiNo;
  fftw_real *Temp1R, *Temp2R, *Temp3R;
  fftw_complex *Temp1C, *Temp2C, *Temp3C;
  MPI_Status status;
  
  // subtract actual density of last minimised Psi from total
  if (R->CurrentMin == Psi->LocalPsiStatus[p_old].PsiType || Psi->LocalPsiStatus[p_old].PsiType == Occupied) {
    if (Psi->LocalPsiStatus[p_old].PsiType == Occupied)
      Temp1R = Dens->DensityArray[TotalLocalDensity];
    else
      Temp1R = Dens->DensityArray[GapLocalDensity];
    Temp2R = Dens->DensityArray[ActualDensity];
    for (i=0; i < Dens->LocalSizeR; i++)
      Temp1R[i] -= Temp2R[i];
  }
    
  //SetArrayToDouble0((double *)Dens->DensityArray[ActualDensity], Dens->TotalSize*2);   // zero actual density (of the current wave function)
  // calculate density of last minimised Psi (but only store if same as current one)
  CalculateOneDensityR(Lat, LevS, Dens, LPsi->LocalPsi[p_old], Dens->DensityArray[ActualDensity], factorDR*Psi->LocalPsiStatus[p_old].PsiFactor, (p == p_old));
  //if (isnan(Dens->DensityArray[ActualDensity][0])) { fprintf(stderr,"(%i) WARNING in UpdateDensityCalculation(): ActualDensity_%i[%i] = NaN!\n", P->Par.me, p_old, 0); Error(SomeError, "NaN-Fehler!"); }
  
  // calculate actual density for current minimised Psi and add up again
  if (R->CurrentMin == Psi->LocalPsiStatus[p_old].PsiType || Psi->LocalPsiStatus[p_old].PsiType == Occupied) {
    if (Psi->LocalPsiStatus[p_old].PsiType == Occupied)
      Temp1R = Dens->DensityArray[TotalLocalDensity];
    else
      Temp1R = Dens->DensityArray[GapLocalDensity];
    Temp2R = Dens->DensityArray[ActualDensity];
    for (i=0; i < Dens->LocalSizeR; i++)
      Temp1R[i] += Temp2R[i];
  }

  if (p != p_old) { // recalc actual density if step to next wave function was made
    SetArrayToDouble0((double *)Dens->DensityArray[ActualDensity], Dens->TotalSize*2);  // set to zero
    CalculateOneDensityR(Lat, LevS, Dens, LPsi->LocalPsi[p], Dens->DensityArray[ActualDensity], factorDR*Psi->LocalPsiStatus[p].PsiFactor, 1);
  }
  
  // calculate complex density from the (real) total local array
  CalculateOneDensityC(Lat, LevS, Dens, Dens->DensityArray[TotalLocalDensity], Dens->DensityCArray[TotalLocalDensity], factorDC);
  CalculateOneDensityC(Lat, LevS, Dens, Dens->DensityArray[GapLocalDensity], Dens->DensityCArray[GapLocalDensity], factorDC);
  // depending on SpinType gather results: real and complex DensityTypes::TotalLocalDensity and native initial density
  switch (Psi->PsiST) {
  case SpinDouble:
     //fprintf(stderr,"(%i) InitDensity SpinDouble\n",P->Par.me);
    if (Psi->LocalPsiStatus[p_old].PsiType == Occupied) { // occupied states
      MPI_Allreduce(Dens->DensityArray[TotalLocalDensity] ,Dens->DensityArray[TotalDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[TotalLocalDensity] ,Dens->DensityCArray[TotalDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    } else { // unoccupied or perturbed
      MPI_Allreduce(Dens->DensityArray[GapLocalDensity] ,Dens->DensityArray[GapDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[GapLocalDensity] ,Dens->DensityCArray[GapDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    }
    break;
  case SpinUp:
    //fprintf(stderr,"(%i) InitDensity SpinUp\n",P->Par.me);
    if (Psi->LocalPsiStatus[p_old].PsiType == Occupied) { // occupied states
      MPI_Allreduce(Dens->DensityArray[TotalLocalDensity] ,Dens->DensityArray[TotalUpDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[TotalLocalDensity] ,Dens->DensityCArray[TotalUpDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
      MPI_Sendrecv( Dens->DensityCArray[TotalUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag1, 
        Dens->DensityCArray[TotalDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag2, 
        P->Par.comm_STInter, &status );
      MPI_Sendrecv( Dens->DensityArray[TotalUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag3, 
        Dens->DensityArray[TotalDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag4, 
        P->Par.comm_STInter, &status );
      Temp1R = Dens->DensityArray[TotalDensity];
      Temp2R = Dens->DensityArray[TotalUpDensity];
      Temp3R = Dens->DensityArray[TotalDownDensity];
      Temp1C = Dens->DensityCArray[TotalDensity];
      Temp2C = Dens->DensityCArray[TotalUpDensity];
      Temp3C = Dens->DensityCArray[TotalDownDensity];
    } else { // unoccupied or perturbed
      MPI_Allreduce(Dens->DensityArray[GapLocalDensity] ,Dens->DensityArray[GapUpDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[GapLocalDensity] ,Dens->DensityCArray[GapUpDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
      MPI_Sendrecv( Dens->DensityCArray[GapUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag7, 
        Dens->DensityCArray[GapDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag8, 
        P->Par.comm_STInter, &status );
      MPI_Sendrecv( Dens->DensityArray[GapUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag9, 
        Dens->DensityArray[GapDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag0, 
        P->Par.comm_STInter, &status );
      Temp1R = Dens->DensityArray[GapDensity];
      Temp2R = Dens->DensityArray[GapUpDensity];
      Temp3R = Dens->DensityArray[GapDownDensity];
      Temp1C = Dens->DensityCArray[GapDensity];
      Temp2C = Dens->DensityCArray[GapUpDensity];
      Temp3C = Dens->DensityCArray[GapDownDensity];
    }
  
    for (i=0; i < Dens->LocalSizeR; i++)
      Temp1R[i] = Temp2R[i]+Temp3R[i];
    for (i=0; i < Dens->LocalSizeC; i++) {
      Temp1C[i].re = Temp2C[i].re+Temp3C[i].re;
      Temp1C[i].im = Temp2C[i].im+Temp3C[i].im;
    }
    break;
  case SpinDown:
    //fprintf(stderr,"(%i) InitDensity SpinDown\n",P->Par.me);
    if (Psi->LocalPsiStatus[p_old].PsiType == Occupied) { // occupied states
      MPI_Allreduce(Dens->DensityArray[TotalLocalDensity] ,Dens->DensityArray[TotalDownDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[TotalLocalDensity] ,Dens->DensityCArray[TotalDownDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
      MPI_Sendrecv( Dens->DensityCArray[TotalDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag2, 
        Dens->DensityCArray[TotalUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag1, 
        P->Par.comm_STInter, &status );
      MPI_Sendrecv( Dens->DensityArray[TotalDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag4, 
        Dens->DensityArray[TotalUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag3, 
        P->Par.comm_STInter, &status );
      Temp1R = Dens->DensityArray[TotalDensity];
      Temp2R = Dens->DensityArray[TotalUpDensity];
      Temp3R = Dens->DensityArray[TotalDownDensity];
      Temp1C = Dens->DensityCArray[TotalDensity];
      Temp2C = Dens->DensityCArray[TotalUpDensity];
      Temp3C = Dens->DensityCArray[TotalDownDensity];
    } else { // unoccupied or perturbed
      MPI_Allreduce(Dens->DensityArray[GapLocalDensity] ,Dens->DensityArray[GapDownDensity] ,Dens->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      MPI_Allreduce(Dens->DensityCArray[GapLocalDensity] ,Dens->DensityCArray[GapDownDensity] ,Dens->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
      // exchange up and down densities
      MPI_Sendrecv( Dens->DensityCArray[GapDownDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag8, 
        Dens->DensityCArray[GapUpDensity], Dens->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag7, 
        P->Par.comm_STInter, &status );
      MPI_Sendrecv( Dens->DensityArray[GapDownDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag0, 
        Dens->DensityArray[GapUpDensity], Dens->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag9, 
        P->Par.comm_STInter, &status );
      Temp1R = Dens->DensityArray[GapDensity];
      Temp2R = Dens->DensityArray[GapUpDensity];
      Temp3R = Dens->DensityArray[GapDownDensity];
      Temp1C = Dens->DensityCArray[GapDensity];
      Temp2C = Dens->DensityCArray[GapUpDensity];
      Temp3C = Dens->DensityCArray[GapDownDensity];
    }

    for (i=0; i < Dens->LocalSizeR; i++)
      Temp1R[i] = Temp2R[i]+Temp3R[i];
    for (i=0; i < Dens->LocalSizeC; i++) {
      Temp1C[i].re = Temp2C[i].re+Temp3C[i].re;
      Temp1C[i].im = Temp2C[i].im+Temp3C[i].im;
    }
    break;
  }
  //if (isnan(c_re(Dens->DensityCArray[TotalDensity][0]))) { fprintf(stderr,"(%i) WARNING in UpdateDensityCalculation(): TotalDensity[%i] = NaN!\n", P->Par.me, 0); Error(SomeError, "NaN-Fehler!"); }
  //if (isnan(c_re(Dens->DensityCArray[GapDensity][0]))) { fprintf(stderr,"(%i) WARNING in UpdateDensityCalculation(): GapDensity[%i] = NaN!\n", P->Par.me, 0); Error(SomeError, "NaN-Fehler!"); }
}

/** Prepares wave functions and densities for the LevelUp.
 * Density::DensityArray's of all DensityTypes of the new level are set to zero. The density is
 * calculated in the usual manner - see CalculateOneDensityR() - from the wave functions for the
 * upper level. Then we have the real density in the upper level and by back-transformation
 * via fft_3d_real_to_complex() we gain the wave function coefficient on the finer mesh.\n
 * Finally, the new ActualDensity is calculated, summed up for all local coefficients and
 * gathered afterwards from all processes - for the various spin cases - into the
 * DensityTypes#TotalDensity.
 * \param *P Problem at hand
 */
void ChangePsiAndDensToLevUp(struct Problem *P)
{
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct LatticeLevel *LevS = R->LevS;
  struct Density *Dens = Lev0->Dens;
  struct LatticeLevel *LevUp = &Lat->Lev[Lev0->LevelNo-1];
  struct Density *DensUp = LevUp->Dens;
  struct Psis *Psi = &Lat->Psi;
  int Index,i,pos;
  struct  fft_plan_3d *plan = Lat->plan; 
  fftw_complex *work = (fftw_complex *)Dens->DensityCArray[TempDensity];
  fftw_complex *tempdestRC = (fftw_complex *)Dens->DensityArray[TempDensity];
  fftw_complex *posfac, *destpos, *destRCS, *destRCD;
  fftw_complex *source, *source0;
  double factorC = 1./Lev0->MaxN;
  double factorDC = R->FactorDensityC;
  double factorDR = R->FactorDensityR;
  int p,g,d;
  MPI_Status status;
  fftw_real *Temp1R, *Temp2R, *Temp3R, *Temp4R, *Temp5R, *Temp6R;
  fftw_complex *Temp1C, *Temp2C, *Temp3C, *Temp4C, *Temp5C, *Temp6C;
  
  // set density arrays of the new level for all types to zero
  for (d=0; d < (R->DoPerturbation ? MaxDensityTypes : MaxInitDensityTypes); d++) {
    if (DensUp->DensityCArray[d] != NULL)
      SetArrayToDouble0((double *)DensUp->DensityCArray[d], Dens->TotalSize*2);
    if (DensUp->DensityArray[d] != NULL)
      SetArrayToDouble0((double *)DensUp->DensityArray[d], Dens->TotalSize*2);
  }
  // for all local Psis do the usual transformation (completing coefficients for all grid vectors, fft, permutation)
  for (p=Psi->LocalNo-1; p >= 0; p--) {
    source = LevS->LPsi->LocalPsi[p];
    source0 = Lev0->LPsi->LocalPsi[p];
    SetArrayToDouble0((double *)tempdestRC, Dens->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];
      //if (isnan(source[i].re)) { fprintf(stderr,"(%i) WARNING in ChangePsiAndDensToLevUp(): source_%i[%i] = NaN!\n", P->Par.me, p, i); Error(SomeError, "NaN-Fehler!"); }
      for (pos=0; pos < LevS->MaxNUp; pos++) {
        destpos[pos].re = source[i].re*posfac[pos].re-source[i].im*posfac[pos].im;
        destpos[pos].im = source[i].re*posfac[pos].im+source[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,(fftw_real *)Dens->DensityCArray[ActualPsiDensity]);
    // now we have density in the upper level, fft back to complex and store it as wave function coefficients
    fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF1, Dens->DensityCArray[ActualPsiDensity], work);
    for (g=0; g < Lev0->MaxG; g++) {
      Index = Lev0->GArray[g].Index;
      source0[g].re = Dens->DensityCArray[ActualPsiDensity][Index].re*factorC; 
      source0[g].im = Dens->DensityCArray[ActualPsiDensity][Index].im*factorC; 
      //if (isnan(source0[g].re)) { fprintf(stderr,"(%i) WARNING in ChangePsiAndDensToLevUp(): source0_%i[%i] = NaN!\n", P->Par.me, p, g); Error(SomeError, "NaN-Fehler!"); }
    }
    if (Lev0->GArray[0].GSq == 0.0)
      source0[g].im = 0.0;
    // calculate real actual density in the new level and sum all up into TotalLocalDensity
    CalculateOneDensityR(Lat, Lev0, DensUp, source0, DensUp->DensityArray[ActualDensity], factorDR*Psi->LocalPsiStatus[p].PsiFactor, 0);
    if (R->CurrentMin == Psi->LocalPsiStatus[p].PsiType || Psi->LocalPsiStatus[p].PsiType == Occupied) {
      if (Psi->LocalPsiStatus[p].PsiType == Occupied)
        Temp1R = DensUp->DensityArray[TotalLocalDensity];
      else
        Temp1R = DensUp->DensityArray[GapLocalDensity];
      Temp2R = DensUp->DensityArray[ActualDensity];
      for (i=0; i < DensUp->LocalSizeR; i++)
        Temp1R[i] += Temp2R[i];
    }
  }
  /*if (R->CurrentMin > UnOccupied) // taken out as it makes no sense as Dens->PertMixed is updated not DensUp->... !
    CalculatePertMixedDensity(P); // updates DensityArray[PertMixedDensity] and changes ActualDensity
  */
  source0 = Lev0->LPsi->LocalPsi[R->ActualLocalPsiNo];
  CalculateOneDensityR(Lat, Lev0, DensUp, source0, DensUp->DensityArray[ActualDensity], factorDR*Psi->LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor, 1); 
  
  CalculateOneDensityC(Lat, Lev0, DensUp, DensUp->DensityArray[TotalLocalDensity], DensUp->DensityCArray[TotalLocalDensity], factorDC);
  CalculateOneDensityC(Lat, Lev0, DensUp, DensUp->DensityArray[GapLocalDensity], DensUp->DensityCArray[GapLocalDensity], factorDC);
  // gather results from all processes for the various spin cases
  switch (Psi->PsiST) {
  case SpinDouble:
     // occupied states
    MPI_Allreduce(DensUp->DensityArray[TotalLocalDensity] ,DensUp->DensityArray[TotalDensity] ,DensUp->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Allreduce(DensUp->DensityCArray[TotalLocalDensity] ,DensUp->DensityCArray[TotalDensity] ,DensUp->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
     // unoccupied
    MPI_Allreduce(DensUp->DensityArray[GapLocalDensity] ,DensUp->DensityArray[GapDensity] ,DensUp->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Allreduce(DensUp->DensityCArray[GapLocalDensity] ,DensUp->DensityCArray[GapDensity] ,DensUp->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    break;
  case SpinUp:
    //fprintf(stderr,"(%i) InitDensity SpinUp\n",P->Par.me);
    // occupied states
    MPI_Allreduce(DensUp->DensityArray[TotalLocalDensity] ,DensUp->DensityArray[TotalUpDensity] ,DensUp->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Allreduce(DensUp->DensityCArray[TotalLocalDensity] ,DensUp->DensityCArray[TotalUpDensity] ,DensUp->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    // unoccupied or perturbed states
    MPI_Allreduce(DensUp->DensityArray[GapLocalDensity] ,DensUp->DensityArray[GapUpDensity] ,DensUp->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Allreduce(DensUp->DensityCArray[GapLocalDensity] ,DensUp->DensityCArray[GapUpDensity] ,DensUp->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);

    // exchange up and down densities
    MPI_Sendrecv( DensUp->DensityCArray[TotalUpDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag1, 
      DensUp->DensityCArray[TotalDownDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag2, 
      P->Par.comm_STInter, &status );
    MPI_Sendrecv( DensUp->DensityArray[TotalUpDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag3, 
      DensUp->DensityArray[TotalDownDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag4, 
      P->Par.comm_STInter, &status );
    MPI_Sendrecv( DensUp->DensityCArray[GapUpDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag7, 
      DensUp->DensityCArray[GapDownDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag8, 
      P->Par.comm_STInter, &status );
    MPI_Sendrecv( DensUp->DensityArray[GapUpDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag9, 
      DensUp->DensityArray[GapDownDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag0, 
      P->Par.comm_STInter, &status );

    Temp1R = DensUp->DensityArray[TotalDensity];
    Temp2R = DensUp->DensityArray[TotalUpDensity];
    Temp3R = DensUp->DensityArray[TotalDownDensity];
    Temp4R = DensUp->DensityArray[GapDensity];
    Temp5R = DensUp->DensityArray[GapUpDensity];
    Temp6R = DensUp->DensityArray[GapDownDensity];
    for (i=0; i < DensUp->LocalSizeR; i++) {
      Temp1R[i] = Temp2R[i]+Temp3R[i];
      Temp4R[i] = Temp5R[i]+Temp6R[i];
    }
    Temp1C = DensUp->DensityCArray[TotalDensity];
    Temp2C = DensUp->DensityCArray[TotalUpDensity];
    Temp3C = DensUp->DensityCArray[TotalDownDensity];
    Temp4C = DensUp->DensityCArray[GapDensity];
    Temp5C = DensUp->DensityCArray[GapUpDensity];
    Temp6C = DensUp->DensityCArray[GapDownDensity];
    for (i=0; i < DensUp->LocalSizeC; i++) {
      Temp1C[i].re = Temp2C[i].re+Temp3C[i].re;
      Temp1C[i].im = Temp2C[i].im+Temp3C[i].im;
      Temp4C[i].re = Temp5C[i].re+Temp6C[i].re;
      Temp4C[i].im = Temp5C[i].im+Temp6C[i].im;
    }
    break;
  case SpinDown:
    //fprintf(stderr,"(%i) InitDensity SpinDown\n",P->Par.me);
    // occupied states
    MPI_Allreduce(DensUp->DensityArray[TotalLocalDensity] ,DensUp->DensityArray[TotalDownDensity] ,DensUp->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Allreduce(DensUp->DensityCArray[TotalLocalDensity] ,DensUp->DensityCArray[TotalDownDensity] ,DensUp->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    // unoccupied or perturbed states
    MPI_Allreduce(DensUp->DensityArray[GapLocalDensity] ,DensUp->DensityArray[GapDownDensity] ,DensUp->LocalSizeR , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
    MPI_Allreduce(DensUp->DensityCArray[GapLocalDensity] ,DensUp->DensityCArray[GapDownDensity] ,DensUp->LocalSizeC*2 , MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);

    // exchange up and down densities
    MPI_Sendrecv( DensUp->DensityCArray[TotalDownDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag2, 
      DensUp->DensityCArray[TotalUpDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag1, 
      P->Par.comm_STInter, &status );
    MPI_Sendrecv( DensUp->DensityArray[TotalDownDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag4, 
      DensUp->DensityArray[TotalUpDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag3, 
      P->Par.comm_STInter, &status );
    MPI_Sendrecv( DensUp->DensityCArray[GapDownDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag8, 
      DensUp->DensityCArray[GapUpDensity], DensUp->LocalSizeC*2, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag7, 
      P->Par.comm_STInter, &status );
    MPI_Sendrecv( DensUp->DensityArray[GapDownDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag0, 
      DensUp->DensityArray[GapUpDensity], DensUp->LocalSizeR, MPI_DOUBLE, P->Par.me_comm_ST, DensityTag9, 
      P->Par.comm_STInter, &status );

    Temp1R = DensUp->DensityArray[TotalDensity];
    Temp2R = DensUp->DensityArray[TotalUpDensity];
    Temp3R = DensUp->DensityArray[TotalDownDensity];
    Temp4R = DensUp->DensityArray[GapDensity];
    Temp5R = DensUp->DensityArray[GapUpDensity];
    Temp6R = DensUp->DensityArray[GapDownDensity];
    for (i=0; i < DensUp->LocalSizeR; i++) {
      Temp1R[i] = Temp2R[i]+Temp3R[i];
      Temp4R[i] = Temp5R[i]+Temp6R[i];
    }
    Temp1C = DensUp->DensityCArray[TotalDensity];
    Temp2C = DensUp->DensityCArray[TotalUpDensity];
    Temp3C = DensUp->DensityCArray[TotalDownDensity];
    Temp4C = DensUp->DensityCArray[GapDensity];
    Temp5C = DensUp->DensityCArray[GapUpDensity];
    Temp6C = DensUp->DensityCArray[GapDownDensity];
    for (i=0; i < DensUp->LocalSizeC; i++) {
      Temp1C[i].re = Temp2C[i].re+Temp3C[i].re;
      Temp1C[i].im = Temp2C[i].im+Temp3C[i].im;
      Temp4C[i].re = Temp5C[i].re+Temp6C[i].re;
      Temp4C[i].im = Temp5C[i].im+Temp6C[i].im;
    }
    break;
  }
}

/** Checks for deviation from Psis::NIDensity.
 * Calculates NIDensity (sum over DensityTypes#TotalDensity and DensityTypes#PertMixedDensity over R)
 * via CalculateNativeIntDens() and checks if it is within MYEPSILON to Psis::NIDensity. Note that sum must
 * be equal to that over DensityTypes#TotalDensity alone, as summed up DensityTypes#PertMixedDensity must
 * be zero (charge conservation).
 * \param *P Problem at hand
 */
void ControlNativeDensity(struct Problem *P)
{
  double NIDensity=0., NIDensityUp=0., NIDensityDown=0.;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct Density *Dens = Lev0->Dens;
  struct Psis *Psi = &Lat->Psi;
  double factorDR = R->FactorDensityR;

  switch (R->CurrentMin) {
    case Occupied:
      NIDensity = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[TotalDensity], 1./factorDR);
      if (Psi->PsiST != SpinDouble) {
        NIDensityUp = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[TotalUpDensity], 1./factorDR);
        NIDensityDown = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[TotalDownDensity], 1./factorDR);
      }
      break;
    default:
    case UnOccupied:
      NIDensity = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[GapDensity], 1./factorDR);
      if (Psi->PsiST != SpinDouble) {
        NIDensityUp = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[GapUpDensity], 1./factorDR);
        NIDensityDown = CalculateNativeIntDens(P, Lev0, Dens->DensityArray[GapDownDensity], 1./factorDR);
      }
      break;
  }
  if (P->Call.out[LeaderOut] && (P->Par.me == 0)) {
    if (fabs(NIDensity - Psi->NIDensity[R->CurrentMin]) >= MYEPSILON) {
      fprintf(stderr, "NativeDensity(L(%i)[%i]) = %g != %g eps (%g >= %g)\n",R->LevSNo,R->CurrentMin,NIDensity,Psi->NIDensity[R->CurrentMin],fabs(NIDensity - Psi->NIDensity[R->CurrentMin]),MYEPSILON);
    } else {
      if (P->Call.out[StepLeaderOut]) fprintf(stderr, "NativeDensity(L(%i) = %g Ok !\n",R->LevSNo,NIDensity);
    }
    if (Psi->PsiST != SpinDouble) {
        if (fabs(NIDensityUp - Psi->NIDensityUp[R->CurrentMin]) >= MYEPSILON) {
          fprintf(stderr, "NativeDensityUp(L(%i)[%i]) = %g != %g eps (%g >= %g)\n",R->LevSNo,R->CurrentMin,NIDensityUp,Psi->NIDensityUp[R->CurrentMin],fabs(NIDensityUp - Psi->NIDensityUp[R->CurrentMin]),MYEPSILON);
        } else {
          if (P->Call.out[StepLeaderOut]) fprintf(stderr, "NativeDensityUp(L(%i) = %g Ok !\n",R->LevSNo,NIDensityUp);
        }
        if (fabs(NIDensityDown - Psi->NIDensityDown[R->CurrentMin]) >= MYEPSILON) {
          fprintf(stderr, "NativeDensityDown(L(%i)[%i]) = %g != %g eps (%g >= %g)\n",R->LevSNo,R->CurrentMin,NIDensityDown,Psi->NIDensityDown[R->CurrentMin],fabs(NIDensityDown - Psi->NIDensityDown[R->CurrentMin]),MYEPSILON);
        } else {
          if (P->Call.out[StepLeaderOut]) fprintf(stderr, "NativeDensityDown(L(%i) = %g Ok !\n",R->LevSNo,NIDensityDown);
        }
    }
  }
}

/** Initialization of Density structure.
 * Sets initial values for Density structure such as Density::TotalSize,
 * Density::LocalSizeC and Density::LocalSizeR, defines on SpinUpDown/SpinDouble which
 * densities shall be calculated in Density::DensityTypeUse, allocates in UseType::InUse case
 * memory for the density arrays, sets them to zero, otherwise to NULL.
 * \param *P Problem at hand
 */
void InitDensity(struct Problem *const P) {
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *Lev, *LevDown;
  struct Density *Dens;
  int i,j,d;
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitDensity\n", P->Par.me);
  for (i = 0; i < Lat->MaxLevel; i++) {
    Lev = &Lat->Lev[i];
    if (i<Lat->MaxLevel-1) {
      LevDown = &Lat->Lev[i+1];
      Lev->Dens = (struct Density *)
             Malloc(sizeof(struct Density),"InitDensity: Dens");
      Dens = Lev->Dens;
      Dens->TotalSize = LevDown->Plan0.plan->LocalSizeC*LevDown->MaxNUp;
      Dens->LocalSizeC = Lev->Plan0.plan->LocalSizeC;
      Dens->LocalSizeR = Lev->Plan0.plan->LocalSizeR;
    } else {
      Lev->Dens = (struct Density *)
             Malloc(sizeof(struct Density),"InitDensity: Dens");
      Dens = Lev->Dens;
      Dens->TotalSize = Lev->Plan0.plan->LocalSizeC;
      Dens->LocalSizeC = Lev->Plan0.plan->LocalSizeC;
      Dens->LocalSizeR = Lev->Plan0.plan->LocalSizeR;
      /*
        Lev->Dens = NULL;
        Dens = Lev->Dens;
      */
    }
    if (i == 0) {
      for (d=0; d < (P->R.DoPerturbation == 1 ? MaxDensityTypes : MaxInitDensityTypes); d++) {
        switch (Lat->Psi.Use) {
        case UseSpinDouble:
          switch (d) {
          case ActualDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalLocalDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalUpDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalDownDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CoreWaveDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case HGcDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case GapDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case GapUpDensity:
            Dens->DensityTypeUse[d] = InUse;  // are needed for Wannier as well!
            break;
          case GapDownDensity:
            Dens->DensityTypeUse[d] = InUse;  // are needed for Wannier as well!
            break;
          case GapLocalDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;         
          case CurrentDensity0:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity1:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity2:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity3:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity4:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity5:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity6:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity7:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity8:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TempDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case Temp2Density:
            Dens->DensityTypeUse[d] = InUse;
            break;
          }
          break;
        case UseSpinUpDown:
          switch (d) {
          case ActualDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalLocalDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalUpDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TotalDownDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CoreWaveDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case HGcDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case GapDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case GapUpDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case GapDownDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case GapLocalDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;         
          case CurrentDensity0:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity1:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity2:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity3:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity4:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity5:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity6:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity7:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case CurrentDensity8:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case TempDensity:
            Dens->DensityTypeUse[d] = InUse;
            break;
          case Temp2Density:
            Dens->DensityTypeUse[d] = InUse;
            break;
          }
          break;
        }
        if (Dens->DensityTypeUse[d] == InUse) {
          Dens->DensityArray[d] = (fftw_real *)
            Malloc(sizeof(fftw_complex)*Dens->TotalSize,"InitDensity: DensityArray");
          for (j=0; j < Dens->TotalSize*2; j++) 
            Dens->DensityArray[d][j] = 0.0;
          Dens->DensityCArray[d] = (fftw_complex *)
            Malloc(sizeof(fftw_complex)*Dens->TotalSize,"InitDensity: DensityCArray");
          for (j=0; j < Dens->TotalSize; j++) {
            Dens->DensityCArray[d][j].re = 0.0;
            Dens->DensityCArray[d][j].im = 0.0;
          }
        } else {
          Dens->DensityArray[d] = NULL;
          Dens->DensityCArray[d] = NULL;
        }
        Dens->LockArray[d] = NotInUse;
        Dens->LockCArray[d] = NotInUse;
      }
    } else {
      if (Dens != NULL)
        for (d=0; d < (P->R.DoPerturbation == 1 ? MaxDensityTypes : MaxInitDensityTypes); d++) {
          Dens->DensityTypeUse[d] = Lat->Lev[i-1].Dens->DensityTypeUse[d];
          Dens->DensityArray[d] = Lat->Lev[i-1].Dens->DensityArray[d];
          Dens->DensityCArray[d] = Lat->Lev[i-1].Dens->DensityCArray[d];
        Dens->LockArray[d] = NotInUse;
        Dens->LockCArray[d] = NotInUse;
        }
    }
  }
}