/** \file init.c
 * Initialization and Readin of parameter files.
 * Initialization Init() is split up into various parts:
 * lattice InitLattice(), first lattice level InitLatticeLevel0(), remaining lattice level InitOtherLevel(),
 * wave functions InitPsis() and their values InitPsisValue() and distribution among the processes InitPsiGroupNo(),
 * Many parts are found in their respective files such as for densities density.c and energies or
 * GramSchmidt gramsch.c.\n
 * Reading of the parameter file is done in ReadParameters() by parsing the file with ParseForParameter(), CreateMainname()
 * separates the path from the main programme file.\n
 * There remain some special routines used during initialization: splitting grid vector index on level change SplitIndex(),
 * calculating the grid vector index GetIndex(), sort criteria on squared grid vector norm GetKeyOneG(), testing whether G's
 * norm is below the energy cutoff TestG, or discerning wethers it's stored only half as often due to symmetry TestGDouble().
 * Also due to density being calculated on a higher level there are certain factors that can be pulled out which can be
 * precalculated CreatePosFacForG().
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers, Frederik Heber
 \date 2000, 2006

  File: init.c
  $Id: init.c,v 1.97 2007-04-10 14:34:08 foo Exp $
*/

#include <stdlib.h>
#include <stdio.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 "helpers.h"
#include "errors.h"
#include "init.h"
#include "mergesort2.h"
#include "myfft.h"
#include "gramsch.h"
#include "density.h"
#include "riemann.h"
#include "output.h"
#include "energy.h"
#include "mymath.h"
#include "pseudo.h"
#include "ions.h"

/** Splits the index of a reciprocal grid vector by the number of nodes per dimension.
 * \param Index current index of the grid vector
 * \param nx number of nodes in x direction
 * \param nz number of nodes in z direction
 * \param *x return value (Index % nx)
 * \param *y return value (Index / nx / nz)
 * \param *z return value (Index / nx % nz)
 */
static void SplitIndex(int Index, const int nx, const int nz, int *x, int *y, int *z) {
  if ((nx == 0) || (nz == 0)) fprintf(stderr,"SplitIndex: nx = %d, nz == %d",nx,nz);
  *x = Index % nx;
  Index /= nx;
  *z = Index % nz;
  *y = Index / nz;
}

/** Prepare mainname (and mainpath).
 * Extracts the path specified from the config file on command line by hopping along the path delimiters,
 * this becomes Files::mainpath, including the filename itself it becomes Files::mainname. 
 * \param *P Problem at hand
 * \param filename pointer to file name string
 */
void CreateMainname(struct Problem *const P, const char* const filename)
{
  char *dummy;
  char *token;
  char *token2;
  const char delimiters[] = "/";
  int stop = 0;
  int setit = 0;
  dummy = (char*) /* Kopie von filename */
    Malloc(strlen(filename)+1, "PrecreateMainname - dummy");
  sprintf(dummy,"%s",filename);
  P->Files.mainname = (char*) /* Speicher fuer P->Files.mainname */
    Malloc(strlen(filename)+strlen(P->Files.filename)+2, "CreateMainname - P->Files.mainname");
  P->Files.mainname[0] = 0;
  /* Jetzt mit strtok mainname vorbereiten */
  if (dummy) {
    if (dummy[0] == '/') setit = 1;   // remember there was delimiter at the very beginning
    token = strtok(dummy, delimiters);
    do {
      token2 = strtok(0, delimiters);   // a subsequent call of strtok to the end of this path piece (between token and token2)
      if (token2) {                     // if it's not the last piece
				if (setit) strcat(P->Files.mainname,"/");
				setit = 1;                      // afterwards always replace the delimiters in the target string
				strcat(P->Files.mainname, token); // copy the path piece
				token = token2;                 // and go to next section
      } else {
				stop = 1;
      }
    } while (!stop);
  }
  if (setit) strcat(P->Files.mainname,"/");
  if (dummy) Free(dummy, "CreateMainname: dummy");
  P->Files.mainpath = (char*) /* Speicher fuer P->Files.mainname */
    Malloc(strlen(P->Files.mainname)+10, "CreateMainname - P->Files.mainpath");
  sprintf(P->Files.mainpath, "%s" ,P->Files.mainname);
  if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) CreateMainame: mainpath: %s\t", P->Par.me, P->Files.mainpath);
  //strcat(P->Files.mainname, P->Files.filename); /* mainname fertigstellen */  
  strcpy(P->Files.mainname, P->Files.filename);
  if (P->Call.out[ReadOut]) fprintf(stderr,"mainname: %s\n", P->Files.mainname);
}

/** Initializing Lattice structure.
 * Calculates unit cell's volume, (transposed) reciprocal base Lattice::ReciBasis
 * and square-rooted Lattice::ReciBasisSQ, preparation of LatticeLevel structure,
 * NFields depends on use of RiemannTensor, calculates the orthonormal reciprocal
 * base and thus how many mesh points for the ECut constraint are at least necessary.
 * \param *P Problem at hand
 * \return 0 - everything went alright, 1 - some problem (volume is zero)
 */
static int InitLattice(struct Problem *const P) {
  struct Lattice *const Lat = &P->Lat;
  double factor = 2.*PI;
  double h[NDIM_NDIM];
  double V[NDIM];
  double NormV;
  int AddNFactor = Lat->AddNFactor = P->Call.AddNFactor;
  int i,j,UpDummy;
  int MaxG[NDIM], N[NDIM], NFactor, Lev1N[NDIM];
  struct LatticeLevel *Lev0;
  int RL=-1;
  // Volume
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitLattice\n", P->Par.me);
  Lat->Volume = RDET3(Lat->RealBasis);
  if (P->Call.out[ValueOut]) fprintf(stderr,"Volume = %g\n",Lat->Volume);
  if (Lat->Volume == 0.0) {
    return(1);
  }
  // square-rooted and normal reciprocal (transposed) Base
  for (i=0; i < NDIM_NDIM; i++) h[i] = Lat->RealBasis[i];
  RTranspose3(h);
  RMatReci3(Lat->ReciBasis, h);
  for (i=0; i < NDIM_NDIM; i++) Lat->ReciBasis[i] *= factor;		// SM(Lat->ReciBasis, factor, NDIM_NDIM);
  for (i=0; i < NDIM; i++) {
    Lat->ReciBasisSQ[i] = RNORMSQ3(&(Lat->ReciBasis[i*NDIM]));
    Lat->RealBasisSQ[i] = RNORMSQ3(&(Lat->RealBasis[i*NDIM]));
    Lat->RealBasisQ[i] = sqrt(Lat->RealBasisSQ[i]);
  }
	// Prepares LatticeLevels
  Lat->Lev = (struct LatticeLevel *)
    Malloc(Lat->MaxLevel*sizeof(struct LatticeLevel),"InitLattice: Lat->Lev");
  for (i=0; i < Lat->MaxLevel; i++)
    Lat->Lev[i].Step = 0;
  Lev0 = Lat->Lev;    // short for &(Lat->Lev[0])
  Lat->LevelSizes = (int *)
    Malloc(sizeof(int)*Lat->MaxLevel,"InitLattice: LevelSizes");
  Lat->LevelSizes[0] = P->Lat.Lev0Factor;
  for (i=1; i < Lat->MaxLevel-1; i++) Lat->LevelSizes[i] = 2; 
  if (Lat->RT.Use == UseRT) {
    RL = Lat->RT.RiemannLevel;
    Lat->LevelSizes[RL-1] = P->Lat.LevRFactor;
  }
  Lat->LevelSizes[Lat->MaxLevel-1] = 1;
  // NFields are next ... depending on use of RiemannTensor
  Lat->MaxNoOfnFields = (int *)
    Malloc(sizeof(int)*Lat->MaxLevel,"InitLattice: MaxNoOfnFields");;
  Lat->NFields = (int **)
    Malloc(sizeof(int*)*Lat->MaxLevel,"InitLattice: NFields");
  switch (Lat->RT.Use) {
  case UseNotRT:
    for (i=0; i<Lat->MaxLevel; i++) {
      if (i==0) {
				Lat->MaxNoOfnFields[i] = 1;
				Lat->NFields[i] = (int *)
				  Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
				Lat->NFields[i][FFTNF1] = 1;
      } else {
				Lat->MaxNoOfnFields[i] = 2;
				Lat->NFields[i] = (int *)
				  Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
				Lat->NFields[i][FFTNF1] = 1;
				Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[i-1]*Lat->LevelSizes[i-1]*Lat->LevelSizes[i-1];
      }
    }
    break;
  case UseRT:
    for (i=0; i<Lat->MaxLevel; i++) {
      switch (i) {
      case 0:
				Lat->MaxNoOfnFields[i] = 2;
				Lat->NFields[i] = (int *)
				  Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
				Lat->NFields[i][FFTNF1] = 1;
				Lat->NFields[i][FFTNF0Vec] = NDIM;
				break;
      case STANDARTLEVEL:
				Lat->MaxNoOfnFields[i] = 4;
				Lat->NFields[i] = (int *)
				  Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
				Lat->NFields[i][FFTNF1] = 1;
				Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[0]*Lat->LevelSizes[0]*Lat->LevelSizes[0];
				Lat->NFields[i][FFTNFSVecUp] = Lat->NFields[i][FFTNFUp]*NDIM;
				Lat->NFields[i][FFTNFSVec] = NDIM;
				break;
      default:
				if (i == RL) {
				  Lat->MaxNoOfnFields[i] = 5;
				  Lat->NFields[i] = (int *)
				    Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
				  Lat->NFields[i][FFTNF1] = 1;
				  Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[i-1]*Lat->LevelSizes[i-1]*Lat->LevelSizes[i-1];
				  UpDummy = 1;
				  for (j=RL-1; j >= STANDARTLEVEL; j--) 
				    UpDummy *= Lat->LevelSizes[j]*Lat->LevelSizes[j]*Lat->LevelSizes[j];  
				  /* UpDummy RL -> S */
				  Lat->NFields[i][FFTNFRMatVecUpS] = UpDummy*NDIM_NDIM*NDIM;
				  
				  UpDummy = 1;
				  for (j=RL-1; j >= 0; j--) 
				    UpDummy *= Lat->LevelSizes[j]*Lat->LevelSizes[j]*Lat->LevelSizes[j];
				  /* UpDummy RL -> 0 */
				  Lat->NFields[i][FFTNFRMatUp0] = UpDummy*NDIM_NDIM;
				  Lat->NFields[i][FFTNFRVecUp0] = UpDummy*NDIM;
				} else {
				  Lat->MaxNoOfnFields[i] = 2;
				  Lat->NFields[i] = (int *)
				    Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
				  Lat->NFields[i][FFTNF1] = 1;
				  Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[i-1]*Lat->LevelSizes[i-1]*Lat->LevelSizes[i-1];
				}
				break;
      }
    }
    break;
  }
  if (P->Call.out[ValueOut]) {
    for (i=0; i<Lat->MaxLevel; i++) {
      fprintf(stderr,"NFields[%i] =",i);
      for (j=0; j<Lat->MaxNoOfnFields[i]; j++) {
				fprintf(stderr," %i",Lat->NFields[i][j]);
      }
      fprintf(stderr,"\n");
    }
  }
  // calculate mesh points per dimension from ECut constraints for 0th level
  for (i=0; i< NDIM; i++) { // by vector product creates orthogonal reciprocal base vectors, takes norm, ...
    VP3(V, &Lat->ReciBasis[NDIM*((i+1)%NDIM)], &Lat->ReciBasis[NDIM*((i+2)%NDIM)]);
    NormV = sqrt(RNORMSQ3(V));
    if (fabs(NormV) < MYEPSILON) fprintf(stderr,"InitLattice: NormV = %lg",NormV);
    for (j=0;j<NDIM;j++) V[j] /= NormV;
    Lat->ReciBasisO[i] = fabs(RSP3(V,&Lat->ReciBasis[NDIM*i]));
    Lat->ReciBasisOSQ[i] = Lat->ReciBasisO[i]*Lat->ReciBasisO[i];
    if (P->Call.out[ValueOut]) fprintf(stderr,"G[%i] = (%g ,%g ,%g) GSQ[%i] = %g GOSQ[%i] = %g\n",i,Lat->ReciBasis[i*NDIM+0],Lat->ReciBasis[i*NDIM+1],Lat->ReciBasis[i*NDIM+2],i,Lat->ReciBasisSQ[i],i,Lat->ReciBasisOSQ[i]);
    if (Lat->ReciBasisOSQ[i] < MYEPSILON) fprintf(stderr,"InitLattice: Lat->ReciBasisOSQ[i] = %lg",Lat->ReciBasisOSQ[i]);
    MaxG[i] = floor(sqrt(2.*Lat->ECut/Lat->ReciBasisOSQ[i]));
    N[i] = 2*MaxG[i]+2;
    NFactor = ( i==2 ? 2 : 1)*P->Par.proc[PEGamma]; /* TESTINIT */
    for (j=1; j < Lat->MaxLevel; j++) NFactor *= Lat->LevelSizes[j];
    if (AddNFactor == 0) fprintf(stderr,"InitLattice: AddNFactor = %d",AddNFactor);
    NFactor *= AddNFactor;
    Lev0->N[i] = NFactor*ceil((double)N[i]/(double)NFactor);
    Lev1N[i] = Lev0->N[i];
    Lev0->N[i] *= Lat->LevelSizes[0];
  }
  if (P->Call.out[LeaderOut] && P->Par.me == 0) fprintf(stderr,"MaxG[] = %i %i %i -> N[] = %i %i %i -> Lev1->N[] %i %i %i -> Lev0->N[] %i %i %i\n",MaxG[0],MaxG[1],MaxG[2],N[0],N[1],N[2],Lev1N[0],Lev1N[1],Lev1N[2],Lev0->N[0],Lev0->N[1],Lev0->N[2]);
  if (P->Files.MeOutMes)
    fprintf(P->Files.SpeedFile,"\nMaxG[] = %i %i %i -> N[] = %i %i %i -> Lev1->N[] %i %i %i -> Lev0->N[] %i %i %i\n",MaxG[0],MaxG[1],MaxG[2],N[0],N[1],N[2],Lev1N[0],Lev1N[1],Lev1N[2],Lev0->N[0],Lev0->N[1],Lev0->N[2]);
  Lat->E = &Lat->Energy[Occupied];

  if (P->Call.out[ValueOut]) {
    fprintf(stderr,"(%i) Lat->Realbasis =  \n",P->Par.me);
    for (i=0; i < NDIM_NDIM; i++) {
      fprintf(stderr,"%lg ",Lat->RealBasis[i]);
      if (i%NDIM == NDIM-1) fprintf(stderr,"\n"); 
    }
    fprintf(stderr,"(%i) Lat->RealbasisQ =  ",P->Par.me);
    for (i=0; i < NDIM; i++)
      fprintf(stderr,"%lg ",Lat->RealBasisQ[i]);
    fprintf(stderr,"\n");
  }
  
  return(0);
}

/** Test if a reciprocal grid vector is still within the energy cutoff Lattice::ECut.
 * The vector is generated with the given coefficients, its norm calculated and compared
 * with twice the ECut.
 * \param *Lat Lattice containing the reciprocal basis vectors and ECut
 * \param g1 first coefficient
 * \param g2 second coefficient
 * \param g3 third coefficient
 * \param *GSq return value for norm of the vector
 * \param G return array of NDIM for the vector 
 * \return 0 - vector ok, 1 -  norm is bigger than twice the energy cutoff
 */
static int TestG(struct Lattice *Lat, int g1, int g2, int g3, double *GSq, double G[NDIM]) {
  G[0] = g1*Lat->ReciBasis[0]+g2*Lat->ReciBasis[3]+g3*Lat->ReciBasis[6];
  G[1] = g1*Lat->ReciBasis[1]+g2*Lat->ReciBasis[4]+g3*Lat->ReciBasis[7];
  G[2] = g1*Lat->ReciBasis[2]+g2*Lat->ReciBasis[5]+g3*Lat->ReciBasis[8];
  *GSq = RNORMSQ3(G);
  if (0.5*(*GSq) <= Lat->ECut) return (1);
  return (0);
}

/** Specifies which DoubleGType a certain combination of coefficients is.
 * Basically, on the (gz=0)-plane the gy half is DoubleGUse, the lower DoubleGNotUse (including
 * the (gx>=0)-line)
 * \param gx x coefficient of reciprocal grid vector
 * \param gy y coefficient of reciprocal grid vector
 * \param gz z coefficient of reciprocal grid vector
 * \return DoubleG0     : if all coefficients are zero <br>
 *         DoubleGUse   : if gz is zero, and gx is strictly positive and gy positive or vice versa <br>
 *         DoubleGNotUse: if gz is zero, and gx is (strictly) positive and gy (strictly) negative or vice versa <br>
 *         DoubleGNot   : if gz is not equal to zero
 */
static enum DoubleGType TestGDouble(int gx, int gy, int gz) { 
  if (gx == 0 && gy == 0 && gz == 0) return DoubleG0;
  if (gz == 0) {
    if (gx > 0) {
      if (gy >= 0) return DoubleGUse;
      else return DoubleGNotUse;
    } else {
      if (gy > 0) return DoubleGUse;
      else return DoubleGNotUse;
    }
  }
  return DoubleGNot;
}

/** Returns squared product of the reciprocal vector as a sort criteria.
 * \param *a OneGData structure array containing the grid vectors
 * \param i i-th reciprocal vector
 * \param *Args not used, specified for contingency
 * \return a[i].OneGData::GSq
 */
static double GetKeyOneG(void *a, int i, void *Args) {
  return(((struct OneGData *)a)[i].GSq);
}

/** Returns component of reciprocal vector as a sort criteria.
 * \param *a OneGData structure array containing the grid vectors
 * \param i i-th reciprocal vector
 * \param *Args specifies index
 * \return a[i].OneGData::G[Args]
 */
static double GetKeyOneG2(void *a, int i, void *Args) {
  //fprintf(stderr,"G[%i] = %e\n", *(int *)Args,((struct OneGData *)a)[i].G[*(int *)Args] );
  return(((struct OneGData *)a)[i].G[*(int *)Args]);
}

/** Copies each entry of OneGData from \a b to \a a.
 * \param *a destination reciprocal grid vector
 * \param i i-th grid vector of \a a
 * \param *b source reciprocal grid vector
 * \param j j-th vector of \a b
 */
static void CopyElementOneG(void *a, int i, void *b, int j) {
  ((struct OneGData *)a)[i].Index = ((struct OneGData *)b)[j].Index;
  ((struct OneGData *)a)[i].GlobalIndex = ((struct OneGData *)b)[j].GlobalIndex;
  ((struct OneGData *)a)[i].GSq = ((struct OneGData *)b)[j].GSq;
  ((struct OneGData *)a)[i].G[0] = ((struct OneGData *)b)[j].G[0];
  ((struct OneGData *)a)[i].G[1] = ((struct OneGData *)b)[j].G[1];
  ((struct OneGData *)a)[i].G[2] = ((struct OneGData *)b)[j].G[2];
}

/** Calculates the index of a reciprocal grid vector by its components.
 * \param *plan
 * \param *LPlan LevelPlan structure
 * \param gx first component
 * \param gy second component
 * \param gz third component
 * \return index
 */
static int GetIndex(struct fft_plan_3d *plan, struct LevelPlan *LPlan, int gx, int gy, int gz) {
  int x = (gx < 0 ? gx+LPlan->N[0] : gx);
  int y = (gy < 0 ? gy+LPlan->N[1] : gy);
  int Y = y/LPlan->LevelFactor;
  int YY = y%LPlan->LevelFactor;
  int z = (gz < 0 ? gz+LPlan->N[2]/2+1 : gz);
  int Z = z/LPlan->LevelFactor;
  int ZZ = z%LPlan->LevelFactor;
  int pwIndex = Z+(plan->NLSR[2]/2+1)*Y;
  int i=0,pwy,pwz,iz;
  for (i=0; i < plan->LocalMaxpw; i++) 
    if (plan->Localpw[i].Index == pwIndex) {
      pwz = i%(plan->NLSR[2]/2+1);
      if (pwz == plan->NLSR[2]/2) {
      	iz = LPlan->N[2]/2;
      } else {
      	iz = pwz*LPlan->LevelFactor+ZZ;
      }
      pwy = i/(plan->NLSR[2]/2+1);
      return (x+LPlan->N[0]*(iz+LPlan->nz*(pwy*LPlan->LevelFactor+YY)));
    }
  return -1;
}

/** Builds a hash table and sets a global index for GArray.
 * The global index is needed in order to initalise the wave functions randomly but solely dependent
 * on the seed value, not on the number of coefficient sharing processes! The creation of the G vectors
 * is checked whether there were duplicates or not by sorting the vectors.
 * \param *P Problem at hand
 * \param *Lat Pointer to Lattice structure
 * \param *Lev Pointer to LatticeLevel structure, containing GArray
 */
static void HashGArray(struct Problem *P, struct Lattice *Lat, struct LatticeLevel *Lev) 
{
  int index, i, g, dupes;
  int myPE, MaxPE;
  struct OneGData *GlobalGArray, *GlobalGArray_sort;
  double **G, **G_send;
  int base = 0;
  int MaxG = 0;
  int AllMaxG = 0;
  MPI_Comm_size(P->Par.comm_ST_Psi, &MaxPE); // Determines the size of the group associated with a communictor
  MPI_Comm_rank(P->Par.comm_ST_Psi, &myPE);  // Determines the rank of the calling process in the communicator  
  
  // first gather GArray from all coefficient sharing processes into one array
  // get maximum number of G in one process
  for (i=0;i<MaxPE; i++) {  // go through every process
    if (MaxG < Lev->AllMaxG[i]) 
      MaxG = Lev->AllMaxG[i];
    AllMaxG += Lev->AllMaxG[i];
  } 
  GlobalGArray = (struct OneGData *) Malloc(AllMaxG*(sizeof(struct OneGData)),"HashGArray: GlobalGArray");
  GlobalGArray_sort = (struct OneGData *) Malloc(AllMaxG*(sizeof(struct OneGData)),"HashGArray: GlobalGArray");
  G = (double **) Malloc(NDIM*(sizeof(double)),"HashGArray: G");
  G_send = (double **) Malloc(NDIM*(sizeof(double)),"HashGArray: G_send");
  for (i=0;i<NDIM;i++) {
    G[i] = (double *) Malloc(MaxPE*MaxG*(sizeof(double)),"HashGArray: G[i]");
    G_send[i] = (double *) Malloc(MaxG*(sizeof(double)),"HashGArray: G_send[i]");
  }
  // fill send array
  for (index=0;index<NDIM;index++) {
    for (i=0;i<Lev->MaxG;i++)
      G_send[index][i] = Lev->GArray[i].G[index];      
    for (;i<MaxG;i++)
      G_send[index][i] = 0.;
  }

  // gather among all processes
  MPI_Allgather(G_send[0], MaxG, MPI_DOUBLE, G[0], MaxG, MPI_DOUBLE, P->Par.comm_ST_Psi);
  MPI_Allgather(G_send[1], MaxG, MPI_DOUBLE, G[1], MaxG, MPI_DOUBLE, P->Par.comm_ST_Psi);
  MPI_Allgather(G_send[2], MaxG, MPI_DOUBLE, G[2], MaxG, MPI_DOUBLE, P->Par.comm_ST_Psi);
  // fill gathered results into GlobalGArray
  //if (P->Call.out[LeaderOut]) 
  fprintf(stderr,"(%i) MaxPE %i / myPE %i, AllMaxG %i / MaxG %i / Lev->MaxG %i, TotalAllMaxG %i, ECut %lg\n", P->Par.me, MaxPE, myPE, AllMaxG, MaxG, Lev->MaxG, Lev->TotalAllMaxG, Lev->ECut/4.);   // Ecut is in Rydberg
  base = 0;
  for (i=0;i<MaxPE; i++) {  // go through every process 
    for (g=0;g<Lev->AllMaxG[i];g++) { // through every G of above process
      for (index=0;index<NDIM;index++)
        GlobalGArray[base+g].G[index] = G[index][MaxG*i+g];
      GlobalGArray[base+g].GlobalIndex = i;
      GlobalGArray[base+g].Index = g;
      //if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray[%i+%i].G[index] = (%e,%e,%e), .index = %i, myPE = %i\n",P->Par.me, base, g, GlobalGArray[base+g].G[0], GlobalGArray[base+g].G[1], GlobalGArray[base+g].G[2], GlobalGArray[base+g].Index, GlobalGArray[base+g].GlobalIndex);
    }
    base += Lev->AllMaxG[i];
  }

  // then sort array
  // sort by x coordinate
  index = 0;
  naturalmergesort(GlobalGArray,GlobalGArray_sort,0,AllMaxG-1,&GetKeyOneG2,&index,&CopyElementOneG);
  //if (P->Call.out[LeaderOut]) fprintf(stderr,"\n\n");
  base = 0;
  //for (g=0;g<AllMaxG;g++) // through every G of above process
    //if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray_sort[%i].G[index] = (%e,%e,%e), .index = %i, myPE= %i\n",P->Par.me, g, GlobalGArray_sort[g].G[0], GlobalGArray_sort[g].G[1], GlobalGArray_sort[g].G[2], GlobalGArray_sort[g].Index, GlobalGArray[base+g].GlobalIndex);
  // sort by y coordinate
  //if (P->Call.out[LeaderOut]) fprintf(stderr,"\n\n");
  index = 1;
  base = 0;
  for (g=0;g<AllMaxG;g++) { // through every G of above process
    if (GlobalGArray[base].G[index-1] != GlobalGArray[g].G[index-1]) {
      naturalmergesort(GlobalGArray_sort,GlobalGArray,base,g-1,&GetKeyOneG2,&index,&CopyElementOneG);
      base = g;
    } 
  }
  if (base != AllMaxG-1) naturalmergesort(GlobalGArray_sort,GlobalGArray,base,AllMaxG-1,&GetKeyOneG2,&index,&CopyElementOneG); // last section
  base = 0;
  //for (g=0;g<AllMaxG;g++) // through every G of above process
    //if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray[%i].G[index] = (%e,%e,%e), .index = %i, myPE= %i\n",P->Par.me, g, GlobalGArray[g].G[0], GlobalGArray[g].G[1], GlobalGArray[g].G[2], GlobalGArray[g].Index, GlobalGArray[base+g].GlobalIndex);
  // sort by z coordinate
  index = 2;
  base = 0;
  for (g=0;g<AllMaxG;g++) { // through every G of above process
    if (GlobalGArray[base].G[index-1] != GlobalGArray[g].G[index-1]) {
      naturalmergesort(GlobalGArray,GlobalGArray_sort,base,g-1,&GetKeyOneG2,&index,&CopyElementOneG);
      base = g;
    } 
  }
  if (base != AllMaxG-1) naturalmergesort(GlobalGArray,GlobalGArray_sort,base,AllMaxG-1,&GetKeyOneG2,&index,&CopyElementOneG); // last section
  base = 0;
  //for (g=0;g<AllMaxG;g++) // through every G of above process
    //if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray_sort[%i].G = (%e,%e,%e), myPE = %i, index = %i\n", P->Par.me, g, GlobalGArray_sort[g].G[0], GlobalGArray_sort[g].G[1], GlobalGArray_sort[g].G[2], GlobalGArray_sort[g].GlobalIndex, GlobalGArray_sort[g].Index);

  
  // now index 'em through and build the hash table
  for (g=0;g<AllMaxG;g++) {// through every G of above process
    Lev->HashG[g].i = GlobalGArray_sort[g].Index;
    Lev->HashG[g].myPE = GlobalGArray_sort[g].GlobalIndex;
    //if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) HashG[%i]: i %i, myPE %i\n",P->Par.me, g, Lev->HashG[g].i, Lev->HashG[g].myPE);
  }
  
  // find the dupes
  dupes = 0;
  for (g=0;g<AllMaxG-1;g++) // through every G of above process
    if ((GlobalGArray_sort[g].G[0] == GlobalGArray_sort[g+1].G[0]) &&
        (GlobalGArray_sort[g].G[1] == GlobalGArray_sort[g+1].G[1]) &&
        (GlobalGArray_sort[g].G[2] == GlobalGArray_sort[g+1].G[2])) {
      dupes++;
      if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) (i, myPE): (%i,%i) and (%i,%i) are duplicate!\n",P->Par.me, GlobalGArray_sort[g].Index, GlobalGArray_sort[g].GlobalIndex, GlobalGArray_sort[g+1].Index, GlobalGArray_sort[g+1].GlobalIndex);    
    }
  if (P->Call.out[LeaderOut]) {
    if (dupes) fprintf(stderr,"(%i) %i duplicates in GArrays overall in Lev[%i]\n",P->Par.me, dupes, Lev->LevelNo);
    else fprintf(stderr,"(%i) There were no duplicate G vectors in Lev[%i] \n",P->Par.me, Lev->LevelNo);
  }
    
  Free(GlobalGArray, "HashGArray: GlobalGArray");
  Free(GlobalGArray_sort, "HashGArray: GlobalGArray_sort");
  for (i=0;i<NDIM;i++) {
    Free(G[i], "HashGArray: G[i]");
    Free(G_send[i], "HashGArray: G_send[i]");
  }
  Free(G, "HashGArray: G");
  Free(G_send, "HashGArray: G_send");
}

/** Initialization of zeroth lattice level.
 * Creates each possible reciprocal grid vector, first counting them all (LatticeLevel::MaxG and LatticeLevel::MaxDoubleG),
 * then allocating memory in LatticeLevel::GArray, creating and storing them all. Afterwards, these are sorted.
 * LatticeLevel::MaxG is gathered from all processes and LatticeLevel::AllMaxG calculated
 *\param *P Problem at hand
 */
static void InitLatticeLevel0(struct Problem *const P) {
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *Lev0 = Lat->Lev;     // short for &(Lat->Lev[0])
  struct RPlans *RPs = &Lev0->Plan0;
  double GSq=0;
  double ECut = Lat->ECut;
  double G[NDIM];
  int x,ly,y,z,Index,d,Y,Z,iY,iZ,iy,iz,lz,X,i,yY,k,AllMaxG;

  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitLatticeLevel0\n", P->Par.me);
  Lat->ECut *= Lat->LevelSizes[0]*Lat->LevelSizes[0];
  Lat->plan = fft_3d_create_plan(P, P->Par.comm_ST_Psi, Lat->ECut, Lat->MaxLevel, Lat->LevelSizes, Lat->ReciBasisOSQ, Lev0->N, Lat->MaxNoOfnFields, Lat->NFields); 
  RPs->plan = &Lat->plan->Levplan[0];
  RPs->cdata = Lat->plan->local_data;
  RPs->rdata = (fftw_real *)Lat->plan->local_data;
  for (d=0; d < NDIM; d++) {
    Lev0->NUp[d] = 1;
    Lev0->NUp0[d] = 1;
    Lev0->NUp1[d] = 0;
  }
  Lev0->MaxN = Lev0->N[0]*Lev0->N[1]*Lev0->N[2];
  Lev0->PosFactorUp = NULL;
  Lev0->LevelNo = 0;
  Lev0->MaxG = 0;
  Lev0->MaxDoubleG = 0;
  Lev0->MaxNUp = 0;
  // goes through each possible reciprocal grid vector, counting them (for MaxG and MaxDoubleG)
  for (i = 0; i < Lat->plan->LocalMaxpw; i++) {
    Index = Lat->plan->Localpw[i].Index;
    Z = Index%(Lat->plan->NLSR[2]/2+1);
    Y = Index/(Lat->plan->NLSR[2]/2+1);
    for (ly = 0; ly < RPs->plan->LevelFactor; ly++) {
      yY = ly + RPs->plan->LevelFactor*Y;
      if (Z==Lat->plan->NLSR[2]/2) {    // if it's the last possible one
				z = RPs->plan->N[2]/2;
				for (X=0; X < RPs->plan->N[0]; X++) {
				  x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
				  y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
				  if (TestG(&P->Lat, x, y, z, &GSq, G)) Lev0->MaxG++;
				}
      } else {
				for (z=(Z)*RPs->plan->LevelFactor; z < (Z+1)*RPs->plan->LevelFactor; z++)
				  for (X=0; X < RPs->plan->N[0]; X++) {
				    x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
				    y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
				    if (TestG(&P->Lat, x, y, z, &GSq, G)) {
				      switch (TestGDouble(x,y,z)) {
				      case DoubleGNot:
								Lev0->MaxG++;
								break;
				      case DoubleG0:
								Lev0->MaxG++;
								break;
				      case DoubleGUse:
								Lev0->MaxG++;
								Lev0->MaxDoubleG++;
								break;
				      case DoubleGNotUse:
								break;
				      }
				    }
				  }
      }
    }
  }
  if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev0->MaxG %i MaxDoubleG %i\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], Lev0->MaxG, Lev0->MaxDoubleG);
  // then allocate memory for all these vectors
  Lev0->GArray = (struct OneGData *)
    Malloc(Lev0->MaxG*sizeof(struct OneGData),"InitLevel0: Lev0->GArray");
  Lat->GArrayForSort = (struct OneGData *)
    Malloc(Lev0->MaxG*sizeof(struct OneGData),"InitLevel0: Lat->GArrayForSort");
  if (Lev0->MaxDoubleG) {
    Lev0->DoubleG = (int *)
      Malloc(2*sizeof(int)*Lev0->MaxDoubleG,"InitLevel0: Lev0->DoubleG");
  } else {
    Lev0->DoubleG = NULL;
  }
  Lev0->MaxDoubleG = 0;
  Lev0->MaxG = 0;
  Lev0->ECut = 0;
  Lev0->G0 = -1;
  // and generate them again, this time noting them down in GArray
  for (i = 0; i < Lat->plan->LocalMaxpw; i++) {
    Index = Lat->plan->Localpw[i].Index;
    Z = Index%(Lat->plan->NLSR[2]/2+1);
    Y = Index/(Lat->plan->NLSR[2]/2+1);
    iZ = i%(Lat->plan->NLSR[2]/2+1);
    iY = i/(Lat->plan->NLSR[2]/2+1);
    for (ly = 0; ly < RPs->plan->LevelFactor; ly++) {
      yY = ly + RPs->plan->LevelFactor*Y;
      iy = ly + RPs->plan->LevelFactor*iY; 
      if (Z==Lat->plan->NLSR[2]/2) {
				z = RPs->plan->N[2]/2;
				iz = RPs->plan->N[2]/2;
				for (X=0; X < RPs->plan->N[0]; X++) {
				  x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
				  y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
				  if (TestG(&P->Lat, x, y, z, &GSq, G)) {
				    Index = X+RPs->plan->N[0]*(iz+RPs->plan->nz*(iy));
				    Lev0->GArray[Lev0->MaxG].Index = Index;
            Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
				    Lev0->GArray[Lev0->MaxG].GSq = GSq;
				    Lev0->GArray[Lev0->MaxG].G[0] = G[0];
				    Lev0->GArray[Lev0->MaxG].G[1] = G[1];
				    Lev0->GArray[Lev0->MaxG].G[2] = G[2];
				    if (GSq > Lev0->ECut) Lev0->ECut = GSq;
				    Lev0->MaxG++;
				  }
				}
      } else {
				for (lz=0; lz < RPs->plan->LevelFactor; lz++) {
				  iz = iZ*RPs->plan->LevelFactor+lz;
				  z = Z*RPs->plan->LevelFactor+lz;
				  for (X=0; X < RPs->plan->N[0]; X++) {
				    x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
				    y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
				    if (TestG(&P->Lat, x, y, z, &GSq, G)) {
				      switch (TestGDouble(x,y,z)) {
				      case DoubleGNot:
								Index = X+RPs->plan->N[0]*(iz+RPs->plan->nz*(iy));
								Lev0->GArray[Lev0->MaxG].Index = Index;
                Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
								Lev0->GArray[Lev0->MaxG].GSq = GSq;
								Lev0->GArray[Lev0->MaxG].G[0] = G[0];
								Lev0->GArray[Lev0->MaxG].G[1] = G[1];
								Lev0->GArray[Lev0->MaxG].G[2] = G[2];
								if (GSq > Lev0->ECut) Lev0->ECut = GSq;
								Lev0->MaxG++;
								break;
				      case DoubleG0:
								Index = X+RPs->plan->N[0]*(iz+RPs->plan->nz*(iy));
								Lev0->GArray[Lev0->MaxG].Index = Index;
                Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
								Lev0->GArray[Lev0->MaxG].GSq = GSq;
								Lev0->GArray[Lev0->MaxG].G[0] = G[0];
								Lev0->GArray[Lev0->MaxG].G[1] = G[1];
								Lev0->GArray[Lev0->MaxG].G[2] = G[2];
								if (GSq > Lev0->ECut) Lev0->ECut = GSq;
								Lev0->MaxG++;
								Lev0->G0 = Index;
								break;
				      case DoubleGUse:
								Index = X+RPs->plan->N[0]*(iz+RPs->plan->nz*(iy));
								Lev0->GArray[Lev0->MaxG].Index = Index;
                Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
								Lev0->GArray[Lev0->MaxG].GSq = GSq;
								Lev0->GArray[Lev0->MaxG].G[0] = G[0];
								Lev0->GArray[Lev0->MaxG].G[1] = G[1];
								Lev0->GArray[Lev0->MaxG].G[2] = G[2];
								if (GSq > Lev0->ECut) Lev0->ECut = GSq;
								Lev0->DoubleG[2*Lev0->MaxDoubleG] = Index;
								Lev0->DoubleG[2*Lev0->MaxDoubleG+1] = GetIndex(Lat->plan,RPs->plan,-x,-y,-z);
								if (Lev0->DoubleG[2*Lev0->MaxDoubleG+1] < 0) Error(SomeError,"InitLevel0: Lev0->DoubleG[2*Lev->MaxDoubleG+1]");
								Lev0->MaxG++;
								Lev0->MaxDoubleG++;
								break;
				      case DoubleGNotUse:
								break;
				      }
				      
				    }
				  }
				}
      }
    }
  }
  Lev0->ECut *= 0.5;
  naturalmergesort(Lev0->GArray,Lat->GArrayForSort,0,Lev0->MaxG-1,&GetKeyOneG,NULL,&CopyElementOneG);
  MPI_Allreduce(&Lev0->ECut , &ECut, 1, MPI_DOUBLE, MPI_MAX, P->Par.comm_ST_Psi );
  if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev0->ECut %g (%g) %g\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], Lev0->ECut, ECut, ECut - Lat->ECut);
  Lev0->ECut = Lat->ECut;
  if (ECut > Lat->ECut) Error(SomeError,"InitLevel0: ECut > Lat->ECut");
  Lev0->AllMaxG = (int *)
    Malloc(sizeof(int)*P->Par.Max_me_comm_ST_Psi,"InitLevel0: AllMaxG");
  MPI_Allgather ( &Lev0->MaxG, 1, MPI_INT,
		  Lev0->AllMaxG, 1, MPI_INT, P->Par.comm_ST_Psi );
  AllMaxG = 0;
  for (k=0; k<P->Par.Max_me_comm_ST_Psi; k++) AllMaxG += Lev0->AllMaxG[k]; 
  if (P->Call.out[NormalOut]) {
    fprintf(stderr,"(%i) Lev0->AllMaxG %i -> RealAllMaxG %i\n", P->Par.mytid, AllMaxG, 2*AllMaxG-1);
  }
  Lev0->TotalAllMaxG = AllMaxG;
  Lev0->TotalRealAllMaxG = 2*AllMaxG-1;
  Lev0->HashG = (struct GDataHash *)
    Malloc(AllMaxG*sizeof(struct GDataHash),"InitLevel0: Lev0->HashG");
    
  HashGArray(P,Lat,Lev0);
}

/** Precalculate Position Factors LatticeLevel::PosFactorUp.
 * Due to the finer resolution of the densities certain transformations are necessary.
 * One position factor can be pulled out and pre-calculated
 * \f[
 *    P_p(G) = \exp(i2\pi \sum^3_{j=1} \frac{g_j p_j}{2N_j} )
 * \f]
 * where \f$N_j\f$ is RPlans::plan::N, p is the position vector and G the reciprocal
 * grid vector, j goes up to NDIM (here 3)
 * \param *P Problem at hand
 * \param *Lat Lattice structure
 * \param *Lev LatticeLevel structure
 * \param *RPsUp RPlans structure
 * \param g number of reciprocal grid points
 * \param x number of mesh points in x direction (in the upper level)
 * \param y number of mesh points in y direction (in the upper level)
 * \param z number of mesh points in z direction (in the upper level)
 */
static void CreatePosFacForG(struct Problem *P, struct Lattice *Lat, const struct LatticeLevel
 *Lev, const struct RPlans *RPsUp, int g, int x, int y, int z)
{
  double posfacreal[NDIM];
  struct RiemannTensor *RT = &Lat->RT;
  struct RunStruct *R = &P->R; 
  int pos[NDIM];
  for (pos[0]=0; pos[0] < Lev->NUp[0]; pos[0]++) {
    for (pos[1]=0; pos[1] < Lev->NUp[1]; pos[1]++) {
      for (pos[2]=0; pos[2] < Lev->NUp[2]; pos[2]++) {
      	posfacreal[0] = (x*pos[0])/(double)(RPsUp->plan->N[0]);
      	posfacreal[1] = (y*pos[1])/(double)(RPsUp->plan->N[1]);
      	posfacreal[2] = (z*pos[2])/(double)(RPsUp->plan->N[2]);
      	Lev->PosFactorUp[pos[2]+Lev->NUp[2]*(pos[1]+Lev->NUp[1]*(pos[0]+Lev->NUp[0]*g))].re = cos(2.*PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
      	Lev->PosFactorUp[pos[2]+Lev->NUp[2]*(pos[1]+Lev->NUp[1]*(pos[0]+Lev->NUp[0]*g))].im = sin(2.*PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
      }
    }
  }
  if (RT->Use == UseRT) {
    if (Lev->LevelNo == R->LevRNo) {
      for (pos[0]=0; pos[0] < RT->NUpLevRS[0]; pos[0]++) {
      	for (pos[1]=0; pos[1] < RT->NUpLevRS[1]; pos[1]++) {
      	  for (pos[2]=0; pos[2] < RT->NUpLevRS[2]; pos[2]++) {
      	    posfacreal[0] = (x*pos[0])/(double)(Lat->Lev[R->LevRSNo].N[0]);
      	    posfacreal[1] = (y*pos[1])/(double)(Lat->Lev[R->LevRSNo].N[1]);
      	    posfacreal[2] = (z*pos[2])/(double)(Lat->Lev[R->LevRSNo].N[2]);
      	    RT->PosFactor[RTPFRtoS][pos[2]+RT->NUpLevRS[2]*(pos[1]+RT->NUpLevRS[1]*(pos[0]+RT->NUpLevRS[0]*g))].re = cos(2.*PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
      	    RT->PosFactor[RTPFRtoS][pos[2]+RT->NUpLevRS[2]*(pos[1]+RT->NUpLevRS[1]*(pos[0]+RT->NUpLevRS[0]*g))].im = sin(2.*PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
      	  }
      	}
      }
      for (pos[0]=0; pos[0] < RT->NUpLevR0[0]; pos[0]++) {
      	for (pos[1]=0; pos[1] < RT->NUpLevR0[1]; pos[1]++) {
      	  for (pos[2]=0; pos[2] < RT->NUpLevR0[2]; pos[2]++) {
      	    posfacreal[0] = (x*pos[0])/(double)(Lat->Lev[R->LevR0No].N[0]);
      	    posfacreal[1] = (y*pos[1])/(double)(Lat->Lev[R->LevR0No].N[1]);
      	    posfacreal[2] = (z*pos[2])/(double)(Lat->Lev[R->LevR0No].N[2]);
      	    RT->PosFactor[RTPFRto0][pos[2]+RT->NUpLevR0[2]*(pos[1]+RT->NUpLevR0[1]*(pos[0]+RT->NUpLevR0[0]*g))].re = cos(2.*PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
      	    RT->PosFactor[RTPFRto0][pos[2]+RT->NUpLevR0[2]*(pos[1]+RT->NUpLevR0[1]*(pos[0]+RT->NUpLevR0[0]*g))].im = sin(2.*PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
      	  }
      	}
      }
    }
  }
}

/** Initialization of first to Lattice::MaxLevel LatticeLevel's.
 * Sets numbers and ratios  of mesh points (maximum (LatticeLevel::MaxN, LatticeLevel::MaxNUp), per dimension
 * (LatticeLevel::N[], LatticeLevel::NUp[])) for the respective levels and its upper level. LatticeLevel::LevelNo and
 * RPlans are set from LevelPlan[], number of reciprocal grid vectors is recounted (LatticeLevel::MaxG,
 * LatticeLevel::MaxDoubleG), then they are created, stored (LatticeLevel::GArray) and their position factor
 * (LatticeLevel::PosFactorUp) recalculated. MaxG is gathered from all processes in LatticeLevel::AllMaxG and
 * LatticeLevel::TotalAllMaxG generated by summing these up.
 * \param *P Problem at hand
 * \sa InitLatticeLevel0()
 */
static void InitOtherLevel(struct Problem *const P) {
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *Lev;
  struct LatticeLevel *LevUp;
  struct RPlans *RPs;
  struct RPlans *RPsUp;
  double GSq=0,ECut=0;
  int lx,x,ly,y,z,Index,Y,YY,Z,ZZ,lz,I,k,AllMaxG;
  int d,MaxGUp,i;
  for (i=1; i < P->Lat.MaxLevel; i++) {   // for each level
    Lev = &Lat->Lev[i];
    LevUp = &Lat->Lev[i-1];
    RPs = &Lev->Plan0;
    RPsUp = &LevUp->Plan0;
    // set number of mesh points
    for (d=0; d < NDIM; d++) {  // for each dimension
      Lev->NUp[d] = Lat->LevelSizes[i-1];
      Lev->NUp0[d] = LevUp->NUp0[d]*Lev->NUp[d];
      Lev->NUp1[d] = (i == STANDARTLEVEL ? 1 : LevUp->NUp1[d]*Lev->NUp[d]);
      Lev->N[d] = LevUp->N[d]/Lev->NUp[d];
    }
    Lev->MaxN = Lev->N[0]*Lev->N[1]*Lev->N[2];
    Lev->MaxNUp = Lev->NUp[0]*Lev->NUp[1]*Lev->NUp[2];
    if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitLatticeLevel %i\tN[] = %i %i %i\n",P->Par.me, i,Lev->N[0],Lev->N[1],Lev->N[2]);
    // fft plan and its data
    RPs->plan = &Lat->plan->Levplan[i];
    RPs->cdata = Lat->plan->local_data;
    RPs->rdata = (fftw_real *)Lat->plan->local_data;
    Lev->LevelNo = i;
    // recount MaxG, MaxDoubleG by creating all g vectors
    Lev->MaxG=0;
    Lev->MaxDoubleG = 0;
    for (MaxGUp=0; MaxGUp < LevUp->MaxG; MaxGUp++) {
      SplitIndex(LevUp->GArray[MaxGUp].Index,RPsUp->plan->N[0],RPsUp->plan->nz,&lx,&ly,&lz);
      Y = ly/RPsUp->plan->LevelFactor;
      YY = ly%RPsUp->plan->LevelFactor;
      Z = lz/RPsUp->plan->LevelFactor;
      ZZ = lz%RPsUp->plan->LevelFactor;
      Index = Lat->plan->Localpw[Z+(Lat->plan->NLSR[2]/2+1)*(Y)].Index;
      ly = YY+RPsUp->plan->LevelFactor*(Index/(Lat->plan->NLSR[2]/2+1));
      z = ZZ+RPsUp->plan->LevelFactor*(Index%(Lat->plan->NLSR[2]/2+1));
      x = ( lx >= RPsUp->plan->N[0]/2+1 ? lx-RPsUp->plan->N[0] : lx);
      y = ( ly >= RPsUp->plan->N[1]/2+1 ? ly-RPsUp->plan->N[1] : ly);
      if ((x%Lev->NUp[0] == 0) && (y%Lev->NUp[1] == 0) && (z%Lev->NUp[2] == 0)) {
				switch (TestGDouble(x/Lev->NUp[0],y/Lev->NUp[1],z/Lev->NUp[2])) {
				case DoubleGNot:
				  Lev->MaxG++;
				  break;
				case DoubleG0:
				  Lev->MaxG++;
				  break;
				case DoubleGUse:
				  Lev->MaxG++; 
				  Lev->MaxDoubleG++;
				  break;
				case DoubleGNotUse:
				  Lev->MaxG--;
				  break;
				}
      }
    }
    if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev%i->MaxG %i MaxDoubleG %i\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], i, Lev->MaxG, Lev->MaxDoubleG);
    Lev->PosFactorUp = (fftw_complex *) Malloc(sizeof(fftw_complex)*Lev->MaxG*Lev->NUp[0]*Lev->NUp[1]*Lev->NUp[2], "InitLatticeLevel: PosFactors");
    if (Lat->RT.Use == UseRT) {
      if (Lev->LevelNo == Lat->RT.RiemannLevel) {	
				Lat->RT.PosFactor[RTPFRtoS] = (fftw_complex *)
				  Malloc(sizeof(fftw_complex)*Lev->MaxG*Lev->NUp1[0]*Lev->NUp1[1]*Lev->NUp1[2], "InitLatticeLevel: RTPosFactors1");
				Lat->RT.PosFactor[RTPFRto0] = (fftw_complex *)
				  Malloc(sizeof(fftw_complex)*Lev->MaxG*Lev->NUp0[0]*Lev->NUp0[1]*Lev->NUp0[2], "InitLatticeLevel: RTPosFactors0");
      }
    }
    // now create the G vectors and store them
    Lev->GArray = (struct OneGData *)
      Malloc(Lev->MaxG*sizeof(struct OneGData),"InitLevel: Lev->GArray");
    if (Lev->MaxDoubleG) {
      Lev->DoubleG = (int *)
				Malloc(2*sizeof(int)*Lev->MaxDoubleG,"InitLevel: Lev0->DoubleG");
    } else {
      Lev->DoubleG = NULL;
    }
    Lev->MaxDoubleG = 0;
    Lev->MaxG=0;
    Lev->ECut = 0;
    Lev->G0 = -1;
    for (MaxGUp=0; MaxGUp < LevUp->MaxG; MaxGUp++) {
      SplitIndex(LevUp->GArray[MaxGUp].Index,RPsUp->plan->N[0],RPsUp->plan->nz,&lx,&ly,&lz);
      Index = lx/Lev->NUp[0]+RPs->plan->N[0]*(lz/Lev->NUp[2]+RPs->plan->nz*(ly/Lev->NUp[1]));
      Y = ly/RPsUp->plan->LevelFactor;
      YY = ly%RPsUp->plan->LevelFactor;
      Z = lz/RPsUp->plan->LevelFactor;
      ZZ = lz%RPsUp->plan->LevelFactor;
      I = Lat->plan->Localpw[Z+(Lat->plan->NLSR[2]/2+1)*(Y)].Index;
      ly = YY+RPsUp->plan->LevelFactor*(I/(Lat->plan->NLSR[2]/2+1));
      z = ZZ+RPsUp->plan->LevelFactor*(I%(Lat->plan->NLSR[2]/2+1));
      x = ( lx >= RPsUp->plan->N[0]/2+1 ? lx-RPsUp->plan->N[0] : lx);
      y = ( ly >= RPsUp->plan->N[1]/2+1 ? ly-RPsUp->plan->N[1] : ly);
      if ((x%Lev->NUp[0] == 0) && (y%Lev->NUp[1] == 0) && (z%Lev->NUp[2] == 0)) {
				switch (TestGDouble(x/Lev->NUp[0],y/Lev->NUp[1],z/Lev->NUp[2])) {
				case DoubleGNot:
				  TestG(&P->Lat, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2], &GSq, Lev->GArray[Lev->MaxG].G);
				  Lev->GArray[Lev->MaxG].Index = Index;
          Lev->GArray[Lev->MaxG].GlobalIndex = 0;
				  Lev->GArray[Lev->MaxG].GSq = GSq;
				  if (GSq > Lev->ECut) Lev->ECut = GSq;
				  CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
				  Lev->MaxG++;
				  break;
				case DoubleG0:
				  TestG(&P->Lat, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2], &GSq, Lev->GArray[Lev->MaxG].G);
				  Lev->GArray[Lev->MaxG].Index = Index;
          Lev->GArray[Lev->MaxG].GlobalIndex = 0;
				  Lev->GArray[Lev->MaxG].GSq = GSq;
				  if (GSq > Lev->ECut) Lev->ECut = GSq;
				  CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
				  Lev->MaxG++;
				  Lev->G0 = Index;
				  break;
				case DoubleGUse:
				  TestG(&P->Lat, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2], &GSq, Lev->GArray[Lev->MaxG].G);
				  Lev->GArray[Lev->MaxG].Index = Index;
          Lev->GArray[Lev->MaxG].GlobalIndex = 0;
				  Lev->GArray[Lev->MaxG].GSq = GSq;
				  if (GSq > Lev->ECut) Lev->ECut = GSq;
				  CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
				  Lev->DoubleG[2*Lev->MaxDoubleG] = Index;
				  Lev->DoubleG[2*Lev->MaxDoubleG+1] = GetIndex(Lat->plan,RPs->plan,-(x/Lev->NUp[0]),-(y/Lev->NUp[1]),-(z/Lev->NUp[2]));
				  if (Lev->DoubleG[2*Lev->MaxDoubleG+1] < 0) Error(SomeError,"Lev->DoubleG[2*Lev->MaxDoubleG+1]");
				  CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
				  Lev->MaxG++;
				  Lev->MaxDoubleG++;
				  break;
				case DoubleGNotUse:
				  break;
				}
      }
    }
    Lev->ECut *= 0.5;
    MPI_Allreduce(&Lev->ECut , &ECut, 1, MPI_DOUBLE, MPI_MAX, P->Par.comm_ST_Psi );
    if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev%i->ECut %g (%g) %g\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], i, Lev->ECut, ECut, ECut - Lat->ECut/(Lev->NUp0[0]*Lev->NUp0[0]));
    Lev->ECut = Lat->ECut/(Lev->NUp0[0]*Lev->NUp0[0]);
    if (ECut > Lev->ECut) Error(SomeError,"InitLevel: ECut > Lat->ECut");
    Lev->AllMaxG = (int *)
      Malloc(sizeof(int)*P->Par.Max_me_comm_ST_Psi,"InitLevel: AllMaxG");
    MPI_Allgather ( &Lev->MaxG, 1, MPI_INT,
		    Lev->AllMaxG, 1, MPI_INT, P->Par.comm_ST_Psi );
    AllMaxG = 0;
    for (k=0; k<P->Par.Max_me_comm_ST_Psi; k++) AllMaxG += Lev->AllMaxG[k]; 
    if (P->Call.out[NormalOut]) {
      fprintf(stderr,"(%i) Lev%i->AllMaxG %i -> RealAllMaxG %i\n", P->Par.mytid, i, AllMaxG, 2*AllMaxG-1);
    }
    Lev->TotalAllMaxG = AllMaxG;
    Lev->TotalRealAllMaxG = 2*AllMaxG-1;
    Lev->HashG = (struct GDataHash *)
      Malloc(AllMaxG*sizeof(struct GDataHash),"InitLevel0: Lev0->HashG");
    HashGArray(P, Lat, Lev);
  }
}

/** Sets the Spin group number for each process.
 * First of all, SpinUps and SpinDowns Psis are completely separated (they have their own communicator groups).
 * Psis::PsiGroup is then just the rank in this Psi group, where Psis::MaxPsiGroup is the number of participating
 * processes. Afterwards, depending on the SpinType Psis:PsiST is set and Psis::LocalNo calculated by modulo of
 * global number of Psis and the number of Psi groups, e.g. 18 psis, 5 processes, make 3 times 4 Psis per process
 * and 2 times 3 Psis per process. Finally, the maximum of all the local number of Psis is taken and stored in
 * Psis::AllMaxLocalNo.
 * \param *P Problem at hand
 */
static void InitPsiGroupNo(struct Problem *const P) {
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  int i, r, r1, NoPsis = 0, type;
  MPI_Status status;
  int AllMaxLocalNoUp, AllMaxLocalNoDown;
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitPsisGroupNo\n", P->Par.me);
  Psi->MaxPsiGroup = P->Par.Max_my_color_comm_ST_Psi;
  Psi->PsiGroup = P->Par.my_color_comm_ST_Psi;
  // now discern between the different spin cases
  switch (Psi->Use) {
  case UseSpinDouble:
    NoPsis = Psi->GlobalNo[PsiMaxNoDouble];
    if (P->Call.out[NormalOut]) fprintf(stderr, "(%i) PsiGlobalNo %i = %i + %i\n", P->Par.me, Psi->GlobalNo[PsiMaxNo], NoPsis, Psi->GlobalNo[PsiMaxAdd]); 
    Psi->PsiST = SpinDouble;    // set SpinType
    break;
  case UseSpinUpDown:
    if (P->Call.out[NormalOut]) fprintf(stderr, "(%i) PsiGlobalNo %i = (%i+%i) + (%i+%i)\n", P->Par.me, Psi->GlobalNo[PsiMaxNo], Psi->GlobalNo[PsiMaxNoUp], Psi->GlobalNo[PsiMaxAdd], Psi->GlobalNo[PsiMaxNoDown], Psi->GlobalNo[PsiMaxAdd]);
    // depending on my spin type (up or down) ...
    switch (P->Par.my_color_comm_ST) {
    case ((int)SpinUp) - 1:
      Psi->PsiST = SpinUp;
      NoPsis = Psi->GlobalNo[PsiMaxNoUp];
      break;
    case ((int)SpinDown) -1:
      Psi->PsiST = SpinDown;
      NoPsis = Psi->GlobalNo[PsiMaxNoDown];
      break;
    }
    break;
  }
  Psi->LocalNo = (NoPsis + P->Lat.Psi.GlobalNo[PsiMaxAdd]) / Psi->MaxPsiGroup;  // local number of Psis per process
  r = (NoPsis + P->Lat.Psi.GlobalNo[PsiMaxAdd]) % Psi->MaxPsiGroup;     // calculate the rest
  Psi->LocalNoAdd = (P->Lat.Psi.GlobalNo[PsiMaxAdd]) / Psi->MaxPsiGroup;  // local number of Psis per process
  r1 = (P->Lat.Psi.GlobalNo[PsiMaxAdd]) % Psi->MaxPsiGroup;     // calculate the rest
  if (P->R.DoUnOccupied && (r1 != 0)) Error(SomeError,"AddPsis must be multiple of PEPsi!\n");
  if (Psi->PsiGroup < r) Psi->LocalNo++;      // distribute the rest into the first processes
  if (Psi->PsiGroup < r1) Psi->LocalNoAdd++;      // distribute the rest into the first processes
  Psi->TypeStartIndex[Occupied] = 0;
  Psi->TypeStartIndex[UnOccupied] = Psi->LocalNo - Psi->LocalNoAdd;
  if (P->R.DoPerturbation == 1) { // six times NoPsis for Perturbed_RxP and Perturbed_P
    for (i=Perturbed_P0; i<= Perturbed_RxP2;i++) {
      Psi->TypeStartIndex[i] = Psi->LocalNo;
      Psi->LocalNo += (NoPsis) / Psi->MaxPsiGroup;  // local number of Psis per process (just one array for all perturbed)
      r = (NoPsis) % Psi->MaxPsiGroup;     // calculate the rest
      if (Psi->PsiGroup < r) Error(SomeError,"Psis must be multiple of PEPsi!\n");
      //if (Psi->PsiGroup < r) Psi->LocalNo++;      // distribute the rest into the first processes
    }
  } else {
    for (i=Perturbed_P0; i<= Perturbed_RxP2;i++)
      Psi->TypeStartIndex[i] = Psi->LocalNo;
  }
  Psi->TypeStartIndex[Extra] = Psi->LocalNo; // always set to end
  Psi->TypeStartIndex[Extra+1] = Psi->LocalNo+1; // always set to end
  
  if (P->Call.out[PsiOut])
    for(type=Occupied;type<=Extra+1;type++)
      fprintf(stderr,"(%i) TypeStartIndex[%i] = %i\n",P->Par.me, type,Psi->TypeStartIndex[type]);

  if (P->Call.out[StepLeaderOut]) fprintf(stderr, "(%i) PsiLocalNo %i \n", P->Par.me, Psi->LocalNo);
  // gather local number of Psis from all processes
  Psi->RealAllLocalNo = (int *)
    Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"FirstInitGramSchData: Psi->RealAllLocalNo");
  MPI_Allgather ( &Psi->LocalNo, 1, MPI_INT,
		  Psi->RealAllLocalNo, 1, MPI_INT, P->Par.comm_ST_PsiT );
  // determine maximum local number, depending on spin
  switch (Psi->PsiST) {
  case SpinDouble:
    MPI_Allreduce(&Psi->LocalNo, &Psi->AllMaxLocalNo, 1, 
		  MPI_INT, MPI_MAX, P->Par.comm_ST_PsiT );
    break;
  case SpinUp:
    MPI_Allreduce(&Psi->LocalNo, &AllMaxLocalNoUp, 1, 
		  MPI_INT, MPI_MAX, P->Par.comm_ST_PsiT );
    MPI_Sendrecv( &AllMaxLocalNoUp, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag1, 
		  &AllMaxLocalNoDown, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag2, 
		  P->Par.comm_STInter, &status );
    Psi->AllMaxLocalNo = MAX(AllMaxLocalNoUp, AllMaxLocalNoDown);
    break;
  case SpinDown:
    MPI_Allreduce(&Psi->LocalNo, &AllMaxLocalNoDown, 1, 
		  MPI_INT, MPI_MAX, P->Par.comm_ST_PsiT );
    MPI_Sendrecv( &AllMaxLocalNoDown, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag2, 
		  &AllMaxLocalNoUp, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag1, 
		  P->Par.comm_STInter, &status );
    Psi->AllMaxLocalNo = MAX(AllMaxLocalNoUp, AllMaxLocalNoDown);
    break;
  }
}

/** Intialization of the wave functions.
 * Allocates/sets OnePsiElementAddData elements to zero, allocates memory for coefficients and
 * pointing to them from Lattice and LatticeLevel.
 * \param *P Problem at hand
 */
static void InitPsis(struct Problem *const P) {
  int i,j,l;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *Lev = R->InitLevS;
  int CreateOld = P->R.UseOldPsi;
  int NoPsis;
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitPsis\n", P->Par.me);
  /// For Lattice::Psi memory for OnePsiElementAddData is allocated and entries zeroed.
  Lat->Psi.AddData = (struct OnePsiElementAddData *) 
    Malloc((Lat->Psi.LocalNo+1)*sizeof(struct OnePsiElementAddData), "InitPsis: Psi->AddData");
  for (i = 0; i <= Lat->Psi.LocalNo; i++) {
    Lat->Psi.AddData[i].T = 0.0;
    Lat->Psi.AddData[i].Lambda = 0.0;
    Lat->Psi.AddData[i].Gamma = 0.0;
    Lat->Psi.AddData[i].Step = 0;
    Lat->Psi.AddData[i].WannierSpread = 0.;
    for (j=0;j<NDIM;j++)
      Lat->Psi.AddData[i].WannierCentre[j] = 0. ;
  }
  // lambda contains H eigenvalues
  switch (Lat->Psi.PsiST) {
    default:
    case SpinDouble:
      Lat->Psi.NoOfPsis = Lat->Psi.GlobalNo[PsiMaxNoDouble];
      Lat->Psi.NoOfTotalPsis = Lat->Psi.NoOfPsis + Lat->Psi.GlobalNo[PsiMaxAdd];      
      break;
    case SpinUp:
      Lat->Psi.NoOfPsis = Lat->Psi.GlobalNo[PsiMaxNoUp];
      Lat->Psi.NoOfTotalPsis = Lat->Psi.NoOfPsis + Lat->Psi.GlobalNo[PsiMaxAdd];      
      break;
    case SpinDown:
      Lat->Psi.NoOfPsis = Lat->Psi.GlobalNo[PsiMaxNoDown];
      Lat->Psi.NoOfTotalPsis = Lat->Psi.NoOfPsis + Lat->Psi.GlobalNo[PsiMaxAdd];      
      break;
  };
  if (P->Call.out[StepLeaderOut]) fprintf(stderr, "(%i) NoOfPsis %i \n", P->Par.me, Lat->Psi.NoOfPsis);
  Lat->Psi.lambda = (double **) Malloc(sizeof(double *)* Lat->Psi.NoOfTotalPsis, "InitPsis: Psi->Lambda");
  Lat->Psi.Overlap = (double **) Malloc(sizeof(double *)* Lat->Psi.NoOfPsis, "InitPsis: Psi->Overlap");
  for (i=0;i<Lat->Psi.NoOfTotalPsis;i++)
    Lat->Psi.lambda[i] = (double *) Malloc(sizeof(double)* Lat->Psi.NoOfTotalPsis, "InitPsis: Psi->Lambda[i]");
  for (i=0;i<Lat->Psi.NoOfPsis;i++)
    Lat->Psi.Overlap[i] = (double *) Malloc(sizeof(double)* Lat->Psi.NoOfPsis, "InitPsis: Psi->Overlap[i]");

  /// Allocation for initial LevelPsi LatticeLevel::LPsi and LevelPsi::PsiDat
  NoPsis = (Lat->Psi.TypeStartIndex[Perturbed_P1] - Lat->Psi.TypeStartIndex[Occupied]); // 
  if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) InitPsis(): Number of Psi arrays in memory %i\n", P->Par.me, NoPsis);
  Lev->LPsi = (struct LevelPsi *)
    Malloc(sizeof(struct LevelPsi),"InitPsis: LPsi");
  Lev->LPsi->PsiDat = (fftw_complex *)
    Malloc((NoPsis+3)*sizeof(fftw_complex)*Lev->MaxG, "InitPsis: LPsi->PsiDat");  // +3 because of TempPsi and TempPsi2 and extra
  if (CreateOld) 
    Lev->LPsi->OldPsiDat = (fftw_complex *)
      Malloc((NoPsis+1)*sizeof(fftw_complex)*Lev->MaxG, "InitPsis: LPsi->PsiDat");
  else
    Lev->LPsi->OldPsiDat = NULL;
  
  /// Allocation for all other levels LatticeLevel::LPsi and LevelPsi::PsiDat, the latter initialised from initial level
  /// here Lev[1].LPsi, because it's the biggest level for the Psis and they might get overwritten otherwise during the perturbed
  /// load and transform-sequence, see ReadSrcPerturbedPsis()
  for (i=1; i < Lat->MaxLevel; i++) { // for all levels
    if (i > 1) Lat->Lev[i].LPsi = (struct LevelPsi *)
		                                Malloc(sizeof(struct LevelPsi),"InitPsis: LPsi");
    Lat->Lev[i].LPsi->PsiDat = Lev->LPsi->PsiDat;
    Lat->Lev[i].LPsi->OldPsiDat = Lev->LPsi->OldPsiDat;
    Lat->Lev[i].LPsi->LocalPsi = (fftw_complex **)
      Malloc((Lat->Psi.LocalNo+1)*sizeof(fftw_complex *),"InitPsi: LocalPsi");
    /// Initialise each LevelPsi::LocalPsi by pointing it to its respective part of LevelPsi::PsiDat
    for (j=0; j <= Lat->Psi.LocalNo; j++) { // for each local Psi
      if (j < Lat->Psi.TypeStartIndex[Perturbed_P0]) { // if it's not a perturbed wave function
        //if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Setting LocalPsi[%i] to %p\n",P->Par.me, j, &Lat->Lev[i].LPsi->PsiDat[j*Lat->Lev[1].MaxG]);
        Lat->Lev[i].LPsi->LocalPsi[j] = &Lat->Lev[i].LPsi->PsiDat[j*Lat->Lev[1].MaxG];
      } else if (j >= Lat->Psi.TypeStartIndex[Extra]) { // if it's the extra wave function
        //if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Setting LocalPsi[%i] to %p\n",P->Par.me, j, &Lat->Lev[i].LPsi->PsiDat[(NoPsis)*Lat->Lev[1].MaxG]);
        Lat->Lev[i].LPsi->LocalPsi[j] = &Lat->Lev[i].LPsi->PsiDat[(NoPsis)*Lat->Lev[1].MaxG];
      } else { // ... a perturbed wave function
        l = Lat->Psi.TypeStartIndex[Perturbed_P0] 
          + ((j - Lat->Psi.TypeStartIndex[Perturbed_P0]) % (Lat->Psi.TypeStartIndex[Perturbed_P1] - Lat->Psi.TypeStartIndex[Perturbed_P0]));
        //if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Setting LocalPsi[%i] to %p\n",P->Par.me, j, &Lat->Lev[i].LPsi->PsiDat[l*Lat->Lev[1].MaxG]);
        Lat->Lev[i].LPsi->LocalPsi[j] = &Lat->Lev[i].LPsi->PsiDat[l*Lat->Lev[1].MaxG];
      }
    }
    Lat->Lev[i].LPsi->TempPsi = &Lat->Lev[i].LPsi->PsiDat[(NoPsis+1)*Lat->Lev[1].MaxG];
    Lat->Lev[i].LPsi->TempPsi2 = &Lat->Lev[i].LPsi->PsiDat[(NoPsis+2)*Lat->Lev[1].MaxG];

    Lat->Lev[i].LPsi->OldLocalPsi = (fftw_complex **)
      Malloc(((Lat->Psi.TypeStartIndex[Perturbed_P0] - Lat->Psi.TypeStartIndex[Occupied])+1)*sizeof(fftw_complex *),"InitPsi: OldLocalPsi");
    for (j=Lat->Psi.TypeStartIndex[Occupied]; j < Lat->Psi.TypeStartIndex[Perturbed_P0]; j++) {  // OldPsis are only needed during the Wannier procedure
      if (CreateOld)
      	Lat->Lev[i].LPsi->OldLocalPsi[j] = &Lat->Lev[i].LPsi->OldPsiDat[j*Lat->Lev[1].MaxG];
      else 
      	Lat->Lev[i].LPsi->OldLocalPsi[j] = NULL;
    }
  }
  Lat->Lev[0].LPsi = NULL;
}

/** Initialization of wave function coefficients.
 * The random number generator's seed is set to 1, rand called numerously depending on the rank in the
 * communicator, the number of LatticeLevel::AllMaxG in each process. Finally, for each local Psi and each
 * reciprocal grid vector the real and imaginary part of the specific coefficient \f$c_{i,G}\f$ are
 * initialized by a \f$1-\frac{1}{2} \frac{rand()}{|G|^2}\f$, where rand() yields a number between 0
 * and the largest integer value.
 * \param *P Problem at hand
 * \param start The init starts at local this local wave function array ...
 * \param end ... and ends at this local array, e.g. Lattice#Psis#LocalNo
 */
void InitPsisValue(struct Problem *const P, int start, int end) {
  int i,j,k, index;
  int myPE;
  int NoBefore = 0;
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  MPI_Comm_rank(P->Par.comm_ST_Psi, &myPE);  // Determines the rank of the calling process in the communicator
  //char name[255];
  //sprintf(name, ".Psis-%i.all", myPE);
  //FILE *psi;
  //OpenFile(P,&psi, name, "w",P->Call.out[ReadOut]);  

  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitPsisValue\n", P->Par.me);
  for (i=0; i < P->Par.my_color_comm_ST_Psi; i++) {
    NoBefore += Lat->Psi.RealAllLocalNo[i];
  }
  srand(R->Seed);
  for (i=0; i<NoBefore; i++)
    for (k=0; k<P->Par.Max_me_comm_ST_Psi; k++)
      for (j=0; j<2*LevS->AllMaxG[k]; j++)
	     rand();
  for (j=0;(j<end) && (j < Lat->Psi.LocalNo);j++)  // for all local Psis
    for (i=0; i < LevS->TotalAllMaxG; i++) {  // for all coefficients
      if ((LevS->HashG[i].myPE == myPE) && (j >= start)) {
        index = LevS->HashG[i].i;
        //if (!j) fprintf(stderr,"(%i) LevS->GArray[%i].GSq = %e\n",P->Par.me, index, LevS->GArray[index].GSq);
        LevS->LPsi->LocalPsi[j][index].re = (1.-2.*(rand()/(double)RAND_MAX))*(LevS->GArray[index].GSq == 0.0 ? 1.0 : 1/LevS->GArray[index].GSq);
        LevS->LPsi->LocalPsi[j][index].im = (1.-2.*(rand()/(double)RAND_MAX))*(LevS->GArray[index].GSq == 0.0 ? 0.0 : 1/LevS->GArray[index].GSq);
        //if (!j) fprintf(stderr,"(%i) LevS->LPsi->LocalPsi[%i][%i/%i] = %e + i%e,  GArray[%i].GSq = %e\n",myPE, j, index, i, LevS->LPsi->LocalPsi[j][index].re, LevS->LPsi->LocalPsi[j][index].im, index, LevS->GArray[index].GSq);
        //fprintf(psi, "LevS->LPsi->LocalPsi[%i][%i] = %e + i%e\n", j, i, LevS->LPsi->LocalPsi[j][index].re, LevS->LPsi->LocalPsi[j][index].im);
      } else {
        rand();
        rand();
      }

/*  for (j=0;j<Lat->Psi.LocalNo;j++)  // for all local Psis
    for (i=0; i < LevS->MaxG; i++) {  // for all coefficients
      if ((R->Seed == 0) && (j >= Lat->Psi.TypeStartIndex[Perturbed_P0])) {
        if (i == 0) fprintf(stderr,"(%i) Initializing wave function %i with 1 over %i\n",P->Par.me, j, LevS->TotalAllMaxG);
        LevS->LPsi->LocalPsi[j][i].re = 1./LevS->TotalAllMaxG;
        LevS->LPsi->LocalPsi[j][i].im = 0.;
      } else {      
        LevS->LPsi->LocalPsi[j][i].re = (1.-2.*(rand()/(double)RAND_MAX))*(LevS->GArray[i].GSq == 0.0 ? 1.0 : 1/LevS->GArray[i].GSq);
        LevS->LPsi->LocalPsi[j][i].im = (1.-2.*(rand()/(double)RAND_MAX))*(LevS->GArray[i].GSq == 0.0 ? 0.0 : 1/LevS->GArray[i].GSq);
      }*/

      /*if ((j==0 && i==0) || (j==1 && i==1) || (j==2 && i==2) || (j==3 && i==3))
	LevS->LPsi->LocalPsi[j][i].re = 1.0;*/
      /*
      if ( j==0 && P->Par.me_comm_ST_PsiT == 0 &&
	   fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i); 
      }
      if ( j==1 && P->Par.me_comm_ST_PsiT == 0 &&
	   fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i); 
      }
      if ( j==2 && P->Par.me_comm_ST_PsiT == 0 &&
	   fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[2] - 1.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i); 
      }
      if ( j==3 && P->Par.me_comm_ST_PsiT == 0 &&
	   fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	   fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0; 
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i); 
      }
      */
      /*
	if ( j==0 && P->Par.me_comm_ST_PsiT == 0 &&
	fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i); 
	}
	if ( j==1 && P->Par.me_comm_ST_PsiT == 0 &&
	fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i); 
	}
	if ( j==0 && P->Par.me_comm_ST_PsiT == 1 &&
	fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[2] - 1.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i); 
	}
	if ( j==1 && P->Par.me_comm_ST_PsiT == 1 &&
	fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON && 
	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	LevS->LPsi->LocalPsi[j][i].re = 1.0; 
	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	}
      */
    }
    //fclose(psi);
}

/* Init Lattice End */


/** Reads parameter from a parsed file.
 * The file is either parsed for a certain keyword or if null is given for
 * the value in row yth and column xth. If the keyword was necessity#critical,
 * then an error is thrown and the programme aborted.
 * \warning value is modified (both in contents and position)!
 * \param verbose 1 - print found value to stderr, 0 - don't
 * \param file file to be parsed
 * \param	name Name of value in file (at least 3 chars!)
 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning 
 * 				(if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
 * 				 best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
 * 				 counted from this unresetted position!)
 * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
 * \param yth In grid case specifying column number, otherwise the yth \a name matching line
 * \param type Type of the Parameter to be read
 * \param value address of the value to be read (must have been allocated)
 * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
 * \param critical necessity of this keyword being specified (optional, critical)
 * \return 1 - found, 0 - not found
 * \sa ReadMainParameterFile Calling routine
 */
int ParseForParameter(int verbose, FILE *file, const char *const name, int sequential, int const xth, int const yth, enum value_type type, void *value, int repetition, enum necessity critical) {
	int i,j;	// loop variables
	int me, length, maxlength = -1;
	long file_position = ftell(file);	// mark current position
	char *dummy2, *dummy1, *dummy, *free_dummy;		// pointers in the line that is read in per step

  if (repetition == 0)
    Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");

  // allocated memory and rank in multi-process word
	dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char),"ParseForParameter: MemAlloc for dummy1 failed");
	MPI_Comm_rank(MPI_COMM_WORLD, &me);

	//fprintf(stderr,"(%i) Parsing for %s\n",me,name);	
	int line = 0;	// marks line where parameter was found
	int found = (type >= grid) ? 0 : (-yth + 1);	// marks if yth parameter name was found
	while((found != repetition)) {
		dummy1 = dummy = free_dummy;
		do {
			fgets(dummy1, 256, file);	// Read the whole line
      if (feof(file)) {
        if ((critical) && (found == 0)) {
          Free(free_dummy, "ParseForParameter: free_dummy");
          Error(InitReading, name);
        } else {
          //if (!sequential)
          fseek(file, file_position, SEEK_SET);  // rewind to start position
          Free(free_dummy, "ParseForParameter: free_dummy");          
          return 0;
        }
      }
			line++;
		} while (((dummy1[0] == '#') || (dummy1[0] == '\n')) && dummy != NULL); // skip commentary and empty lines
    
    // C++ getline removes newline at end, thus re-add
    if (strchr(dummy1,'\n') == NULL) {
      i = strlen(dummy1);
      dummy1[i] = '\n';
      dummy1[i+1] = '\0';
    }
    
		if (dummy1 == NULL) {
			if (verbose) fprintf(stderr,"(%i) Error reading line %i\n",me,line);
		} else {
			//fprintf(stderr,"(%i) Reading next line %i: %s\n",me, line, dummy1);
		}
		// Seek for possible end of keyword on line if given ...
		if (name != NULL) {
			dummy = strchr(dummy1,'\t');	// set dummy on first tab or space which ever nearer
      dummy2 = strchr(dummy1, ' ');  // if not found seek for space
      if ((dummy != NULL) || (dummy2 != NULL)) {
        if ((dummy == NULL) || ((dummy2 != NULL) && (dummy2 < dummy))) 
          dummy = dummy2;
      }
			if (dummy == NULL) {
				dummy = strchr(dummy1, '\n');	// set on line end then (whole line = keyword)
				//fprintf(stderr,"(%i)Error: Cannot find tabs or spaces on line %i in search for %s\n", me, line, name);
        //Free(free_dummy, "ParseForParameter: free_dummy");
		    //Error(FileOpenParams, NULL);			
			} else {
				//fprintf(stderr,"(%i) found tab at %li\n",me,(long)((char *)dummy-(char *)dummy1));
			}
		} else dummy = dummy1;
		// ... and check if it is the keyword!
		if ((name == NULL) || ((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0) && (dummy-dummy1 == strlen(name)))) {
			found++; // found the parameter!
			//fprintf(stderr,"(%i) found %s at line %i\n",me, name, line);		
			
			if (found == repetition) {
				for (i=0;i<xth;i++) {	// i = rows
					if (type >= grid) {
            // grid structure means that grid starts on the next line, not right after keyword
						dummy1 = dummy = free_dummy;
						do {
							fgets(dummy1, 256, file);	// Read the whole line, skip commentary and empty ones
              if (feof(file)) {
                if ((critical) && (found == 0)) {
                  Free(free_dummy, "ParseForParameter: free_dummy");
                  Error(InitReading, name);
                } else {
                  //if (!sequential) 
                  fseek(file, file_position, SEEK_SET);  // rewind to start position
                  Free(free_dummy, "ParseForParameter: free_dummy");
                  return 0;
                }
              }
							line++;
						} while ((dummy1[0] == '#') || (dummy1[0] == '\n'));
						if (dummy1 == NULL){
							if (verbose) fprintf(stderr,"(%i) Error reading line %i\n", me, line);
						} else {
							//fprintf(stderr,"(%i) Reading next line %i: %s\n", me, line, dummy1);
						}
					} else { // simple int, strings or doubles start in the same line
						while ((*dummy == '\t') || (*dummy == ' '))		// skip interjacent tabs and spaces
							dummy++;
					}
          // C++ getline removes newline at end, thus re-add
          if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
            j = strlen(dummy1);
            dummy1[j] = '\n';
            dummy1[j+1] = '\0';
          }
	
					int start = (type >= grid) ? 0 : yth-1 ;
					for (j=start;j<yth;j++) { // j = columns
						// check for lower triangular area and upper triangular area
						if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
							*((double *)value) = 0.0;
							//fprintf(stderr,"%f\t",*((double *)value));
							value += sizeof(double);
						} else {
							// otherwise we must skip all interjacent tabs and spaces and find next value
							dummy1 = dummy;
							dummy = strchr(dummy1, '\t');	// seek for tab or space (space if appears sooner)
              dummy2 = strchr(dummy1, ' ');
              //fprintf(stderr,"dummy (%p) is %c\t dummy2 (%p) is %c\n", dummy, (dummy == NULL ? '-' : *dummy), dummy2, (dummy2 == NULL ? '-' : *dummy2));
              if ((dummy != NULL) || (dummy2 != NULL)) {
                if ((dummy == NULL) || ((dummy2 != NULL) && (dummy2 < dummy))) 
                  dummy = dummy2;
                //while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
                //  dummy++;
              }
/*							while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))		// skip some more tabs and spaces if necessary
								dummy++;*/
              //fprintf(stderr,"dummy is %c\n", *dummy);
              dummy2 = strchr(dummy1, '#');
              if ((dummy == NULL) || ((dummy2 != NULL) && (dummy1 >= dummy2))) { // if still zero returned or in comment area ...
								dummy = strchr(dummy1, '\n');	 // ... at line end then
								if ((i < xth-1) && (type < 4)) {	// check if xth value or not yet
                  if (critical) {
                    if (verbose) fprintf(stderr,"(%i)Error: EoL or comment at %i and still missing %i value(s) for parameter %s\n", me, line, yth-j, name);
                    Free(free_dummy, "ParseForParameter: free_dummy");
                    Error(InitReading, name);
                  } else {
                    //if (!sequential) 
                    fseek(file, file_position, SEEK_SET);  // rewind to start position
                    Free(free_dummy, "ParseForParameter: free_dummy");
                    return 0;
                  }
							    //Error(FileOpenParams, NULL);			
								}
							} else {
								//fprintf(stderr,"(%i) found tab at %i\n",me,(char *)dummy-(char *)free_dummy);
							}
							//fprintf(stderr,"(%i) value from %i to %i\n",me,(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
							switch(type) {
                case (row_int):
                  *((int *)value) = atoi(dummy1);
                  if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"(%i) %s = ",me, name);
                  if (verbose) fprintf(stderr,"%i\t",*((int *)value));
                    value += sizeof(int);
                  break;
                case (row_double): 
								case(grid):
								case(lower_trigrid):
								case(upper_trigrid):
									*((double *)value) = atof(dummy1);
									if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"(%i) %s = ",me, name);
									if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
									 value += sizeof(double);
									break;
								case(double_type):
									*((double *)value) = atof(dummy1);
								  if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s = %lg\n",me, name, *((double *) value));
									//value += sizeof(double);
									break;
								case(int_type):
									*((int *)value) = atoi(dummy1);
								  if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s = %i\n",me, name, *((int *) value));
									//value += sizeof(int);
									break;
								default:
								case(string_type):
                  if (value != NULL) {
                    if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
                    length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
  									strncpy((char *)value, dummy1, length);  // copy as much
                    ((char *)value)[length] = '\0';  // and set end marker
  									if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s is '%s' (%i chars)\n",me,name,((char *) value), length);
  									//value += sizeof(char);
                  } else {
                    Error(SomeError, "ParseForParameter: given string array points to NULL!");
                  }
								break;
							}
						}
						while (*dummy == '\t')
							dummy++;
					}
				}
			}
		}
	}	
  if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
	if (!sequential) fseek(file, file_position, SEEK_SET);	// rewind to start position
  Free(free_dummy, "ParseForParameter: free_dummy");
	//fprintf(stderr, "(%i) End of Parsing\n\n",me);
	
	return (found);	// true if found, false if not
}

/** Readin of main parameter file.
 * creates RunStruct and FileData, opens \a filename \n
 * Only for process 0: and then scans the whole thing.
 * 
 * \param P Problem structure to be filled
 * \param filename const char array with parameter file name
 * \sa ParseForParameter
 */
void ReadParameters(struct Problem *const  P, const char *const filename) {
  /* Oeffne Hauptparameterdatei */
  struct RunStruct *R = &P->R;
  struct FileData *F = &P->Files;
  FILE *file;
  int i, di, me;
  //double BField[NDIM];
  
  MPI_Comm_rank(MPI_COMM_WORLD, &me);
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadMainParameterFile\n", me);
  file = fopen(filename, "r");
  if (file == NULL) {
    fprintf(stderr,"(%i)Error: Cannot open main parameter file: %s\n", me, filename);
    Error(FileOpenParams, NULL);
  } else if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)File is open: %s\n", me, filename);

  /* Namen einlesen */

  P->Files.filename = MallocString(MAXDUMMYSTRING,"ReadParameters: filename");
	ParseForParameter(P->Call.out[ReadOut],file, "mainname", 0, 1, 1, string_type, P->Files.filename, 1, critical);
  //debug(P,"mainname");
  CreateMainname(P, filename);
  P->Files.default_path = MallocString(MAXDUMMYSTRING,"ReadParameters: default_path");
  ParseForParameter(P->Call.out[ReadOut],file, "defaultpath", 0, 1, 1, string_type, P->Files.default_path, 1, critical);
  P->Files.pseudopot_path = MallocString(MAXDUMMYSTRING,"ReadParameters: pseudopot_path");
  ParseForParameter(P->Call.out[ReadOut],file, "pseudopotpath", 0, 1, 1, string_type, P->Files.pseudopot_path, 1, critical);
  ParseForParameter(P->Call.out[ReadOut],file,"ProcPEGamma", 0, 1, 1, int_type, &(P->Par.proc[PEGamma]), 1, critical);
  ParseForParameter(P->Call.out[ReadOut],file,"ProcPEPsi", 0, 1, 1, int_type, &(P->Par.proc[PEPsi]), 1, critical);
  if (P->Call.proc[PEPsi]) { /* Kommandozeilenoption */
    P->Par.proc[PEPsi]  = P->Call.proc[PEPsi];
    P->Par.proc[PEGamma]  = P->Call.proc[PEGamma];
  }
  /*
    Run
  */
  if (!ParseForParameter(P->Call.out[ReadOut],file,"Seed", 0, 1, 1, int_type, &(R->Seed), 1, optional))
    R->Seed = 1;
  if (!ParseForParameter(P->Call.out[ReadOut],file,"WaveNr", 0, 1, 1, int_type, &(R->WaveNr), 1, optional))
    R->WaveNr = 0;
  if (!ParseForParameter(P->Call.out[ReadOut],file,"Lambda", 0, 1, 1, double_type, &R->Lambda, 1, optional))
    R->Lambda = 1.;

  if(!ParseForParameter(P->Call.out[ReadOut],file,"DoOutOrbitals", 0, 1, 1, int_type, &F->DoOutOrbitals, 1, optional)) {
    F->DoOutOrbitals = 0;
  } else {
    if (F->DoOutOrbitals < 0) F->DoOutOrbitals = 0;
    if (F->DoOutOrbitals > 1) F->DoOutOrbitals = 1;
  }
  ParseForParameter(P->Call.out[ReadOut],file,"DoOutVis", 0, 1, 1, int_type, &F->DoOutVis, 1, critical);
  if (F->DoOutVis < 0) F->DoOutVis = 0;
  if (F->DoOutVis > 2) F->DoOutVis = 1;
  if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorPlane", 0, 1, 1, int_type, &R->VectorPlane, 1, optional))
    R->VectorPlane = -1;
  if (R->VectorPlane < -1 || R->VectorPlane > 2) {
    fprintf(stderr,"(%i) ! -1 <= R-VectorPlane <= 2: We simulate three-dimensional worlds! Disabling current vector plot: -1\n", P->Par.me);
    R->VectorPlane = -1;
  }
  if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorCut", 0, 1, 1, double_type, &R->VectorCut, 1, optional))
    R->VectorCut = 0.;
  ParseForParameter(P->Call.out[ReadOut],file,"DoOutMes", 0, 1, 1, int_type, &F->DoOutMes, 1, critical);
  if (F->DoOutMes < 0) F->DoOutMes = 0;
  if (F->DoOutMes > 1) F->DoOutMes = 1;
  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutCurr", 0, 1, 1, int_type, &F->DoOutCurr, 1, optional))
    F->DoOutCurr = 0;
  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutNICS", 0, 1, 1, int_type, &F->DoOutNICS, 1, optional))
    F->DoOutNICS = 0;
  if (F->DoOutCurr < 0) F->DoOutCurr = 0;
  if (F->DoOutCurr > 1) F->DoOutCurr = 1; 
  ParseForParameter(P->Call.out[ReadOut],file,"AddGramSch", 0, 1, 1, int_type, &R->UseAddGramSch, 1, critical);
  if (R->UseAddGramSch < 0) R->UseAddGramSch = 0;
  if (R->UseAddGramSch > 2) R->UseAddGramSch = 2;
  if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)UseAddGramSch = %i\n",me,R->UseAddGramSch);
  if(!ParseForParameter(P->Call.out[ReadOut],file,"CommonWannier", 0, 1, 1, int_type, &R->CommonWannier, 1, optional)) {
    R->CommonWannier = 0;
  } else {
    if (R->CommonWannier < 0) R->CommonWannier = 0;
    if (R->CommonWannier > 4) R->CommonWannier = 4;
  }
  if(!ParseForParameter(P->Call.out[ReadOut],file,"SawtoothStart", 0, 1, 1, double_type, &P->Lat.SawtoothStart, 1, optional)) {
    P->Lat.SawtoothStart = 0.01;
  } else {
    if (P->Lat.SawtoothStart < 0.) P->Lat.SawtoothStart = 0.;
    if (P->Lat.SawtoothStart > 1.) P->Lat.SawtoothStart = 1.;
  }
  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoBrent", 0, 1, 1, int_type, &R->DoBrent, 1, optional)) {
    R->DoBrent = 0;
  } else {
    if (R->DoBrent < 0) R->DoBrent = 0;
    if (R->DoBrent > 1) R->DoBrent = 1;
  }
   
  
  if (!ParseForParameter(P->Call.out[ReadOut],file,"MaxOuterStep", 0, 1, 1, int_type, &R->MaxOuterStep, 1, optional))
    R->MaxOuterStep = 0;
  if (!ParseForParameter(P->Call.out[ReadOut],file,"MaxStructOptStep", 0, 1, 1, int_type, &R->MaxStructOptStep, 1, optional)) {
    if (R->MaxOuterStep != 0) { // if StructOptStep not explicitely specified, then use OuterStep
      R->MaxStructOptStep = R->MaxOuterStep;
    } else {
      R->MaxStructOptStep = 100;
    }
  }
  // the following values are MD specific and should be dropped in further revisions
	if (!ParseForParameter(P->Call.out[ReadOut],file,"Deltat", 0, 1, 1, double_type, &R->delta_t, 1, optional))
    R->delta_t = 1.;
	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutVisStep", 0, 1, 1, int_type, &R->OutVisStep, 1, optional))
    R->OutVisStep = 10;
	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutSrcStep", 0, 1, 1, int_type, &R->OutSrcStep, 1, optional))
    R->OutSrcStep = 5;
	if (!ParseForParameter(P->Call.out[ReadOut],file,"TargetTemp", 0, 1, 1, double_type, &R->TargetTemp, 1, optional))
    R->TargetTemp = 0;
    
  if (!ParseForParameter(P->Call.out[ReadOut],file,"EpsWannier", 0, 1, 1, double_type, &R->EpsWannier, 1, optional))
    R->EpsWannier = 1e-8;
  
  // stop conditions
  R->t = 0;
  //if (R->MaxOuterStep <= 0) R->MaxOuterStep = 1;
  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxOuterStep = %i\n",me,R->MaxOuterStep);
  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)Deltat = %g\n",me,R->delta_t);
  ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiStep", 0, 1, 1, int_type, &R->MaxPsiStep, 1, critical);
  if (R->MaxPsiStep <= 0) R->MaxPsiStep = 3;
  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxPsiStep = %i\n",me,R->MaxPsiStep);
  
  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStep", 0, 1, 1, int_type, &R->MaxMinStep, 1, critical);
	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsTotalE", 0, 1, 1, double_type, &R->RelEpsTotalEnergy, 1, critical);
	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsKineticE", 0, 1, 1, double_type, &R->RelEpsKineticEnergy, 1, critical);
  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStopStep", 0, 1, 1, int_type, &R->MaxMinStopStep, 1, critical);
  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinGapStopStep", 0, 1, 1, int_type, &R->MaxMinGapStopStep, 1, critical);
  if (R->MaxMinStep <= 0) R->MaxMinStep = R->MaxPsiStep;
  if (R->MaxMinStopStep < 1) R->MaxMinStopStep = 1;
  if (R->MaxMinGapStopStep < 1) R->MaxMinGapStopStep = 1;
  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxMinStep = %i\tRelEpsTotalEnergy: %e\tRelEpsKineticEnergy %e\tMaxMinStopStep %i\n",me,R->MaxMinStep,R->RelEpsTotalEnergy,R->RelEpsKineticEnergy,R->MaxMinStopStep);
  
  ParseForParameter(P->Call.out[ReadOut],file,"MaxInitMinStep", 0, 1, 1, int_type, &R->MaxInitMinStep, 1, critical);
	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsTotalE", 0, 1, 1, double_type, &R->InitRelEpsTotalEnergy, 1, critical);
	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsKineticE", 0, 1, 1, double_type, &R->InitRelEpsKineticEnergy, 1, critical);
  ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinStopStep", 0, 1, 1, int_type, &R->InitMaxMinStopStep, 1, critical);
  ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &R->InitMaxMinGapStopStep, 1, critical);
  if (R->MaxInitMinStep <= 0) R->MaxInitMinStep = R->MaxPsiStep;
  if (R->InitMaxMinStopStep < 1) R->InitMaxMinStopStep = 1;
  if (R->InitMaxMinGapStopStep < 1) R->InitMaxMinGapStopStep = 1;
  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)INIT:MaxMinStep = %i\tRelEpsTotalEnergy: %e\tRelEpsKineticEnergy %e\tMaxMinStopStep %i\tMaxMinGapStopStep %i\n",me,R->MaxInitMinStep,R->InitRelEpsTotalEnergy,R->InitRelEpsKineticEnergy,R->InitMaxMinStopStep,R->InitMaxMinStopStep);
  
  // Unit cell and magnetic field
  ParseForParameter(P->Call.out[ReadOut],file, "BoxLength", 0, 3, 3, lower_trigrid, &P->Lat.RealBasis, 1, critical); /* Lattice->RealBasis */
  ParseForParameter(P->Call.out[ReadOut],file, "DoPerturbation", 0, 1, 1, int_type, &R->DoPerturbation, 1, optional);
  if (!R->DoPerturbation) {
    if (!ParseForParameter(P->Call.out[ReadOut],file, "BField", 0, 1, 3, grid, &R->BField, 1, optional)) { // old parameter for perturbation with field direction
      if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) No perturbation specified.\n", P->Par.me);
      if (F->DoOutCurr || F->DoOutNICS) {
        if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) DoOutCurr/DoOutNICS = 1 though no DoPerturbation specified: setting DoOutCurr = 0\n", P->Par.me);
        F->DoOutCurr = 0;
        F->DoOutNICS = 0;
        R->DoPerturbation = 0;
      }
    } else {
      R->DoPerturbation = 1;
    }
  }
  if (!P->Call.WriteSrcFiles && R->DoPerturbation) {
    if (P->Call.out[ReadOut]) fprintf(stderr,"(%i)Perturbation: Always writing source files (\"-w\").\n", P->Par.me);
    P->Call.WriteSrcFiles = 1;
  }
  RMatReci3(P->Lat.InvBasis,  P->Lat.RealBasis);

  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoFullCurrent", 0, 1, 1, int_type, &R->DoFullCurrent, 1, optional))
    R->DoFullCurrent = 0;
  if (R->DoFullCurrent < 0) R->DoFullCurrent = 0;
  if (R->DoFullCurrent > 2) R->DoFullCurrent = 2; 
  if (R->DoPerturbation == 0) {
    if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) No Magnetic field: RunStruct#DoFullCurrent = 0\n", P->Par.me);
    R->DoFullCurrent = 0;
  }

	ParseForParameter(P->Call.out[ReadOut],file,"ECut", 0, 1, 1, double_type, &P->Lat.ECut, 1, critical);
  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) ECut = %e -> %e Rydberg\n", me, P->Lat.ECut, 2.*P->Lat.ECut);
  ParseForParameter(P->Call.out[ReadOut],file,"MaxLevel", 0, 1, 1, int_type, &P->Lat.MaxLevel, 1, critical);
  ParseForParameter(P->Call.out[ReadOut],file,"Level0Factor", 0, 1, 1, int_type, &P->Lat.Lev0Factor, 1, critical);
  if (P->Lat.Lev0Factor < 2) {
    if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set Lev0Factor to 2!\n",me);
    P->Lat.Lev0Factor = 2;
  }
  ParseForParameter(P->Call.out[ReadOut],file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
  if (di >= 0 && di < 2) {
    P->Lat.RT.Use = (enum UseRiemannTensor)di;
  } else {
    Error(SomeError, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT"); 
  }
  if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!RiemannTensor = %i\n",me,P->Lat.RT.Use), critical;
  switch (P->Lat.RT.Use) {
	  case UseNotRT:
	    if (P->Lat.MaxLevel < 2) {
	      if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set MaxLevel to 2!\n",me);
	      P->Lat.MaxLevel = 2;
	    }
	    P->Lat.LevRFactor = 2;
	    P->Lat.RT.ActualUse = inactive;
	    break;
	  case UseRT:
	    if (P->Lat.MaxLevel < 3) {
	      if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set MaxLevel to 3!\n",me);
	      P->Lat.MaxLevel = 3;
	    }
	    ParseForParameter(P->Call.out[ReadOut],file,"RiemannLevel", 0, 1, 1, int_type, &P->Lat.RT.RiemannLevel, 1, critical);
		  if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!RiemannLevel = %i\n",me,P->Lat.RT.RiemannLevel);
	    if (P->Lat.RT.RiemannLevel < 2) {
	      if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set RiemannLevel to 2!\n",me);
	      P->Lat.RT.RiemannLevel = 2;
	    } 
	    if (P->Lat.RT.RiemannLevel > P->Lat.MaxLevel-1) {
	      if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set RiemannLevel to %i (MaxLevel-1)!\n",me,P->Lat.MaxLevel-1);
	      P->Lat.RT.RiemannLevel = P->Lat.MaxLevel-1;
	    }
	    ParseForParameter(P->Call.out[ReadOut],file,"LevRFactor", 0, 1, 1, int_type, &P->Lat.LevRFactor, 1, critical);
	    if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!LevRFactor = %i\n",me,P->Lat.LevRFactor);
	    if (P->Lat.LevRFactor < 2) {
	      if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set LevRFactor to 2!\n",me);
	      P->Lat.LevRFactor = 2;
	    } 
	    P->Lat.Lev0Factor = 2;
	    P->Lat.RT.ActualUse = standby;
	    break;
  }
  ParseForParameter(P->Call.out[ReadOut],file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
  if (di >= 0 && di < 2) {
    P->Lat.Psi.Use = (enum UseSpinType)di;
  } else {
    Error(SomeError, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown"); 
  }
  if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiType = %i\n",me,P->Lat.Psi.Use);
  for (i=0; i<MaxPsiNoType; i++)
    P->Lat.Psi.GlobalNo[i] = 0;  
  switch (P->Lat.Psi.Use) {
  case UseSpinDouble:
  	ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiDouble", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDouble], 1, critical);
    if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], 1, optional))
      P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
    if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
      R->DoUnOccupied = 1;
    else 
      R->DoUnOccupied = 0;
    if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!MaxPsiDouble = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoDouble],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
    //P->Lat.Psi.GlobalNo[PsiMaxNoDouble] += P->Lat.Psi.GlobalNo[PsiMaxAdd];  no more adding up on occupied ones
    P->Lat.Psi.GlobalNo[PsiMaxNo] = P->Lat.Psi.GlobalNo[PsiMaxNoDouble]; 
    break;
  case UseSpinUpDown:
    if (P->Par.proc[PEPsi] % 2) Error(SomeError,"UseSpinUpDown & P->Par.proc[PEGamma] % 2");
    ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoUp", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoUp], 1, critical);
    ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoDown", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDown], 1, critical);
    if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], 1, optional))
      P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
    if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
      R->DoUnOccupied = 1;
    else 
      R->DoUnOccupied = 0;
    if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiMaxPsiUp = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoUp],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
    if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiMaxPsiDown = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoDown],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
    //P->Lat.Psi.GlobalNo[PsiMaxNoUp] += P->Lat.Psi.GlobalNo[PsiMaxAdd];
    //P->Lat.Psi.GlobalNo[PsiMaxNoDown] += P->Lat.Psi.GlobalNo[PsiMaxAdd];
    if (abs(P->Lat.Psi.GlobalNo[PsiMaxNoUp]-P->Lat.Psi.GlobalNo[PsiMaxNoDown])>1) Error(SomeError,"|Up - Down| > 1");
    P->Lat.Psi.GlobalNo[PsiMaxNo] = P->Lat.Psi.GlobalNo[PsiMaxNoUp] + P->Lat.Psi.GlobalNo[PsiMaxNoDown];
    break;
  }
  
  IonsInitRead(P, file);
  fclose(file);
}

/** Writes parameters as a parsable config file to \a filename.
 * \param *P Problem at hand
 * \param filename name of config file to be written
 */
void WriteParameters(struct Problem *const  P, const char *const filename) {
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct FileData *F = &P->Files;
  struct Ions *I = &P->Ion;
  FILE *output;
  double BorAng = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
  int i, it, nr = 0;
  if ((P->Par.me == 0) && (output = fopen(filename,"w"))) {
    //fprintf(stderr, "(%i) Writing config file %s\n",P->Par.me, filename);
    
    fprintf(output,"# ParallelCarParinello - main configuration file - optimized geometry\n");
    fprintf(output,"\n");
    fprintf(output,"mainname\t%s\t# programm name (for runtime files)\n", F->mainname);
    fprintf(output,"defaultpath\t%s\t# where to put files during runtime\n", F->default_path);
    fprintf(output,"pseudopotpath\t%s\t# where to find pseudopotentials\n", F->pseudopot_path);
    fprintf(output,"\n");
    fprintf(output,"ProcPEGamma\t%i\t# for parallel computing: share constants\n", P->Par.proc[PEGamma]);
    fprintf(output,"ProcPEPsi\t%i\t# for parallel computing: share wave functions\n", P->Par.proc[PEPsi]);
    fprintf(output,"DoOutVis\t%i\t# Output data for OpenDX\n", F->DoOutVis);
    fprintf(output,"DoOutMes\t%i\t# Output data for measurements\n", F->DoOutMes);
    fprintf(output,"DoOutCurr\t%i\t# Ouput current density for OpenDx\n", F->DoOutCurr);
    fprintf(output,"DoPerturbation\t%i\t# Do perturbation calculate and determine susceptibility and shielding\n", R->DoPerturbation);
    fprintf(output,"DoFullCurrent\t%i\t# Do full perturbation\n", R->DoFullCurrent);
    fprintf(output,"DoOutOrbitals\t%i\t# Output all orbitals on end\n", F->DoOutOrbitals);
    fprintf(output,"CommonWannier\t%i\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center\n", R->CommonWannier);
    fprintf(output,"SawtoothStart\t%lg\t# Absolute value for smooth transition at cell border \n", Lat->SawtoothStart);
    fprintf(output,"VectorPlane\t%i\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot\n", R->VectorPlane);
    fprintf(output,"VectorCut\t%lg\t# Cut plane axis value\n", R->VectorCut);
    fprintf(output,"AddGramSch\t%i\t# Additional GramSchmidtOrtogonalization to be safe (1 - per MinStep, 2 - per PsiStep)\n", R->UseAddGramSch);
    fprintf(output,"Seed\t\t%i\t# initial value for random seed for Psi coefficients\n", R->Seed);
    fprintf(output,"DoBrent\t%i\t# Brent line search instead of CG with second taylor approximation\n", R->DoBrent);
    fprintf(output,"\n");
    fprintf(output,"MaxOuterStep\t%i\t# number of MolecularDynamics/Structure optimization steps\n", R->MaxOuterStep);
    fprintf(output,"Deltat\t%lg\t# time in atomic time per MD step\n", R->delta_t);
    fprintf(output,"OutVisStep\t%i\t#  Output visual data every ...th step\n", R->OutVisStep);
    fprintf(output,"OutSrcStep\t%i\t# Output \"restart\" data every ..th step\n", R->OutSrcStep);
    fprintf(output,"TargetTemp\t%lg\t# Target temperature in Hartree\n", R->TargetTemp);
    fprintf(output,"Thermostat\t%s", I->ThermostatNames[I->Thermostat]);
    switch(I->Thermostat) {
      default:
      case 0: // None
        break;
      case 1: // Woodcock
        fprintf(output,"\t%i", R->ScaleTempStep);
        break;
      case 2: // Gaussian
        fprintf(output,"\t%i", R->ScaleTempStep);
        break;
      case 3: // Langevin
        fprintf(output,"\t%lg\t%lg", R->TempFrequency, R->alpha);
        break;
      case 4: // Berendsen
        fprintf(output,"\t%lg", R->TempFrequency);
        break;
      case 5: // NoseHoover
        fprintf(output,"\t%lg", R->HooverMass);
        break;
    }
    fprintf(output,"\t# Thermostat to be used with parameters\n");
    fprintf(output,"MaxPsiStep\t%i\t# number of Minimisation steps per state (0 - default)\n", R->MaxPsiStep);
    fprintf(output,"EpsWannier\t%lg\t# tolerance value for spread minimisation of orbitals\n", R->EpsWannier);
    fprintf(output,"\n");
    fprintf(output,"# Values specifying when to stop\n");
    fprintf(output,"MaxMinStep\t%i\t# Maximum number of steps\n", R->MaxMinStep);
    fprintf(output,"RelEpsTotalE\t%lg\t# relative change in total energy\n", R->RelEpsTotalEnergy);
    fprintf(output,"RelEpsKineticE\t%lg\t# relative change in kinetic energy\n", R->RelEpsKineticEnergy);
    fprintf(output,"MaxMinStopStep\t%i\t# check every ..th steps\n", R->MaxMinStopStep);
    fprintf(output,"MaxMinGapStopStep\t%i\t# check every ..th steps\n", R->MaxMinGapStopStep);
    fprintf(output,"\n");
    fprintf(output,"# Values specifying when to stop for INIT, otherwise same as above\n");
    fprintf(output,"MaxInitMinStep\t%i\t# Maximum number of steps\n", R->MaxInitMinStep);
    fprintf(output,"InitRelEpsTotalE\t%lg\t# relative change in total energy\n", R->InitRelEpsTotalEnergy);
    fprintf(output,"InitRelEpsKineticE\t%lg\t# relative change in kinetic energy\n", R->InitRelEpsKineticEnergy);
    fprintf(output,"InitMaxMinStopStep\t%i\t# check every ..th steps\n", R->InitMaxMinStopStep);
    fprintf(output,"InitMaxMinGapStopStep\t%i\t# check every ..th steps\n", R->InitMaxMinGapStopStep);
    fprintf(output,"\n");
    fprintf(output,"BoxLength\t\t\t# (Length of a unit cell)\n");
    fprintf(output,"%lg\t\n", Lat->RealBasis[0]*BorAng);
    fprintf(output,"%lg\t%lg\t\n", Lat->RealBasis[3]*BorAng, Lat->RealBasis[4]*BorAng);
    fprintf(output,"%lg\t%lg\t%lg\t\n", Lat->RealBasis[6]*BorAng, Lat->RealBasis[7]*BorAng, Lat->RealBasis[8]*BorAng);
    // FIXME
    fprintf(output,"\n");
    fprintf(output,"ECut\t\t%lg\t# energy cutoff for discretization in Hartrees\n", Lat->ECut/(Lat->LevelSizes[0]*Lat->LevelSizes[0]));
    fprintf(output,"MaxLevel\t%i\t# number of different levels in the code, >=2\n", P->Lat.MaxLevel);
    fprintf(output,"Level0Factor\t%i\t# factor by which node number increases from S to 0 level\n", P->Lat.Lev0Factor);
    fprintf(output,"RiemannTensor\t%i\t# (Use metric)\n", P->Lat.RT.Use);
    switch (P->Lat.RT.Use) {
      case 0: //UseNoRT
        break;
      case 1: // UseRT
        fprintf(output,"RiemannLevel\t%i\t# Number of Riemann Levels\n", P->Lat.RT.RiemannLevel);
        fprintf(output,"LevRFactor\t%i\t# factor by which node number increases from 0 to R level from\n", P->Lat.LevRFactor);
        break;
    }
    fprintf(output,"PsiType\t\t%i\t# 0 - doubly occupied, 1 - SpinUp,SpinDown\n", P->Lat.Psi.Use);
    // write out both types for easier changing afterwards
  //  switch (PsiType) {
  //    case 0:
        fprintf(output,"MaxPsiDouble\t%i\t# here: specifying both maximum number of SpinUp- and -Down-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoDouble]);
  //      break;
  //    case 1:
        fprintf(output,"PsiMaxNoUp\t%i\t# here: specifying maximum number of SpinUp-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoUp]);
        fprintf(output,"PsiMaxNoDown\t%i\t# here: specifying maximum number of SpinDown-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoDown]);
  //      break;
  //  }
    fprintf(output,"AddPsis\t\t%i\t# Additional unoccupied Psis for bandgap determination\n", P->Lat.Psi.GlobalNo[PsiMaxAdd]);
    fprintf(output,"\n");
    fprintf(output,"RCut\t\t%lg\t# R-cut for the ewald summation\n", I->R_cut);
    fprintf(output,"StructOpt\t%i\t# Do structure optimization beforehand\n", I->StructOpt);
    fprintf(output,"IsAngstroem\t%i\t# 0 - Bohr, 1 - Angstroem\n", I->IsAngstroem);
    fprintf(output,"RelativeCoord\t0\t# whether ion coordinates are relative (1) or absolute (0)\n");
    fprintf(output,"MaxTypes\t%i\t# maximum number of different ion types\n", I->Max_Types);
    fprintf(output,"\n");
    // now elements ...
    fprintf(output,"# Ion type data (PP = PseudoPotential, Z = atomic number)\n");
    fprintf(output,"#Ion_TypeNr.\tAmount\tZ\tRGauss\tL_Max(PP)L_Loc(PP)IonMass\tName\tSymbol\n");
    for (i=0; i < I->Max_Types; i++)
      fprintf(output,"Ion_Type%i\t%i\t%i\t%lg\t%i\t%i\t%lg\t%s\t%s\n", i+1, I->I[i].Max_IonsOfType, I->I[i].Z, I->I[i].rgauss, I->I[i].l_max, I->I[i].l_loc, I->I[i].IonMass/Units2Electronmass, &I->I[i].Name[0], &I->I[i].Symbol[0]);
    // ... and each ion coordination
    fprintf(output,"#Ion_TypeNr._Nr.R[0]\tR[1]\tR[2]\tMoveType (0 MoveIon, 1 FixedIon) U[0]\tU[1]\tU[2] (velocities in %s)\n", (I->IsAngstroem) ? "Angstroem" : "atomic units");
    for (i=0; i < I->Max_Types; i++)
      for (it=0;it<I->I[i].Max_IonsOfType;it++) {
        fprintf(output,"Ion_Type%i_%i\t%2.9f\t%2.9f\t%2.9f\t%i\t%2.9f\t%2.9f\t%2.9f\t# Number in molecule %i\n", i+1, it+1, I->I[i].R[0+NDIM*it]*BorAng, I->I[i].R[1+NDIM*it]*BorAng, I->I[i].R[2+NDIM*it]*BorAng, I->I[i].IMT[it], I->I[i].U[0+NDIM*it]*BorAng, I->I[i].U[1+NDIM*it]*BorAng, I->I[i].U[2+NDIM*it]*BorAng, ++nr);
      }
    fflush(output);
    fclose(output);
  }
}

/** Main Initialization.
 * Massive calling of initializing super-functions init...():
 * -# InitOutputFiles()\n
 *    Opening and Initializing of output measurement files.
 * -# InitLattice()\n
 *    Initializing Lattice structure.
 * -# InitRunLevel()\n
 *    Initializing RunLevel structure.
 * -# InitLatticeLevel0()\n
 *    Initializing LatticeLevel structure.
 * -# InitOtherLevel()\n
 *    Initialization of first to Lattice::MaxLevel LatticeLevel's.
 * -# InitRun()\n
 *    Initialization of RunStruct structure.
 * -# InitPsis()\n
 *    Initialization of wave functions as a whole.
 * -# InitDensity()\n
 *    Initializing real and complex Density arrays.
 * -# either ReadSrcFiles() or InitPsisValues()\n
 *    read old calculations from source file if desired or initialise with
 *    random values InitPsisValue().
 * -# InitRiemannTensor()\n
 *    Initializing RiemannTensor structure.
 * -# FirstInitGramSchData()\n
 *    Sets up OnePsiElement as an MPI structure and subsequently initialises all of them.
 * 
 * Now follow various (timed) tests:
 * -# GramSch()\n
 *    (twice) to orthonormalize the random values
 * -# TestGramSch()\n
 *    Test if orbitals now fulfill kronecker delta relation.
 * -# InitGradients()\n
 *    Initializing Gradient structure, allocating its arrays.
 * -# InitDensityCalculation()\n
 *    Calculate naive, real and complex densities and sum them up.
 * -# InitPseudoRead()\n
 *    Read PseudoPot'entials from file.
 * -# InitEnergyArray()\n
 *    Initialize Energy arrays.
 * -# InitPsiEnergyCalculation()\n
 *    Calculate all psi energies and sum up.
 * -# CalculateDensityEnergy()\n
 *    Calculate Density energy.
 * -# CalculateIonsEnergy()\n
 *    Calculate ionic energy.
 * -# CalculateEwald()\n
 *    Calculate ionic ewald summation.
 * -# EnergyAllReduce()\n
 *    Gather energies from all processes and sum up.
 * -# EnergyOutput()\n
 *    First output to file.
 * -# InitOutVisArray()\n
 *    First output of visual for Opendx.
 * 
 * \param *P Problem at hand
 */
void Init(struct Problem *const P) {
  struct RunStruct *R = &P->R;
  InitOutputFiles(P);
  InitLattice(P);
  InitRunLevel(P);
  InitLatticeLevel0(P);
  InitOtherLevel(P);
  InitPsiGroupNo(P);
  InitRun(P);
  InitPsis(P);
  InitDensity(P);
  if (P->Call.ReadSrcFiles) {
    if (ReadSrcIons(P)) {
      if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) Ionic structure read successfully.\n", P->Par.me);
    } else {
      if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) No readible Ionic structure found.\n", P->Par.me);
    }
  } else {
    InitPsisValue(P, 0, P->Lat.Psi.TypeStartIndex[Perturbed_P0]);
  }
  InitRiemannTensor(P);
  if (P->Lat.RT.ActualUse == standby && R->LevRSNo == R->LevSNo) P->Lat.RT.ActualUse = active;
  CalculateRiemannTensorData(P);
  FirstInitGramSchData(P,&P->Lat.Psi);
  /* Test */
  SpeedMeasure(P, InitSimTime, StartTimeDo);
  SpeedMeasure(P, InitGramSchTime, StartTimeDo);
  GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);  
  ResetGramSch(P, &P->Lat.Psi);
  if (R->DoUnOccupied == 1) {
    ResetGramSchTagType(P, &P->Lat.Psi, UnOccupied, NotOrthogonal);
    GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
  }
  SpeedMeasure(P, InitGramSchTime, StopTimeDo);
  SpeedMeasure(P, InitSimTime, StopTimeDo);  
  TestGramSch(P, R->LevS, &P->Lat.Psi, -1);
  ResetGramSchTagType(P, &P->Lat.Psi, UnOccupied, NotUsedToOrtho);
  InitGradients(P, &P->Grad);
  SpeedMeasure(P, InitSimTime, StartTimeDo);
  SpeedMeasure(P, InitDensityTime, StartTimeDo);
  InitDensityCalculation(P);
  SpeedMeasure(P, InitDensityTime, StopTimeDo);
  SpeedMeasure(P, InitSimTime, StopTimeDo);
  InitPseudoRead(P, P->Files.pseudopot_path);
  InitEnergyArray(P);
  SpeedMeasure(P, InitSimTime, StartTimeDo);
  InitPsiEnergyCalculation(P, Occupied); /* STANDARTLEVEL */
  CalculateDensityEnergy(P, 1); /* STANDARTLEVEL */
  CalculateIonsEnergy(P);
  CalculateEwald(P, 1);
  EnergyAllReduce(P);
/*  SetCurrentMinState(P,UnOccupied);
  InitPsiEnergyCalculation(P,UnOccupied);
  CalculateGapEnergy(P);
  EnergyAllReduce(P);
  SetCurrentMinState(P,Occupied);*/
  SpeedMeasure(P, InitSimTime, StopTimeDo);
  R->LevS->Step++;
  EnergyOutput(P,0);
  InitOutVisArray(P);
  Free(P->Lat.GArrayForSort, "Init: P->Lat.GArrayForSort");
  P->Speed.InitSteps++;
}