source: pcp/src/output.c@ e24bbb

/** \file output.c
 * Output of forces and energies, Visuals (density and ions for OpenDX).
 * 
 * Herein all the functions concerning the output of data are gathered:
 * Initialization of the FileData structure InitOutVisArray() and the files themselves 
 * InitOutputFiles(),
 * Saving of a calculated state OutputVisSrcFiles() - 1. Psis OutputSrcPsiDensity(), 2. Ions
 * OutSrcIons() - or retrieving Psis ReadSrcPsiDensity() and Ions ReadSrcIons(),
 * Preparation (in case of RiemannTensor use) CalculateOutVisDensityPos(), OutVisPosRTransformPosNFRto0(),
 * Output of visual data (for OpenDx explorer) OutputVis() - 1. OpenDX files CreateDensityOutputGeneral(),
 * 2. electronic densitiy data OutVisDensity() (uses MPI sending density coefficients OutputOutVisDensity()
 * and receiving and writing CombineOutVisDensity()), 3. ion data OutVisIons() -
 * Closing of files CloseOutputFiles().
 * 
 * There are some more routines: OutputCurrentDensity() outputs the current density for each magnetic field
 * direction, OutputVisAllOrbital() outputs all orbitals of a certain minimsation type and
 * TestReadnWriteSrcDensity() checks whether a certain minimisation group can be written and read again
 * correctly.
 * 
 * There are some helpers that open files with certain filenames, making extensive usage of all
 * the suffixes defined in here: OpenFile(), OpenFileNo(), OpenFileNo2(), OpenFileNoNo() and
 * OpenFileNoPost().
 *
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers, Frederik Heber
 \date 2000, 2006

  File: output.c
  $Id: output.c,v 1.51.2.2 2007-04-21 12:55:50 foo Exp $
*/
#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<math.h>
//#include<sys/time.h>
#include <time.h>
#include<unistd.h>
#include"data.h"
#include"density.h"
#include"errors.h"
#include"gramsch.h"
#include"helpers.h"
#include "init.h"
#include"myfft.h"
#include"mymath.h"
#include"output.h"
#include"pdbformat.h"
#include"perturbed.h"


/* Konvention: Rueckgabe 0 einer Funktion, bedeutet keinen Fehler (entsprechend exitcode 0) */ 
/* Oeffnet Datei P->mainname+"..." mit what*/

/** Opens a file with FileData::mainname+suffix.
 * Both the path under which the main programme resides and the path to the config file are
 * tried subsequently.
 * \param *P Problem at hand
 * \param **file file pointer array
 * \param *suffix suffix after mainname
 * \param *what access parameter for the file: r, w, rw, r+, w+ ...
 * \param verbose 1 - print status, 0 - don't print anything
 * \return 0 - could not open file, 1 - file is open
 */
int OpenFile(struct Problem *P, FILE** file, const char* suffix, const char* what, int verbose)
{
  char* name; /* Zu erzeugender Dateiname */
  name = (char*)
    Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 1,"OpenFile");
  sprintf(name, "%s%s%s", P->Files.default_path, P->Files.mainname, suffix);
  *file = fopen(name, what);
  if (*file == NULL) {
    fprintf(stderr,"(%i) Normal access failed: name %s, what %s\n", P->Par.me, name, what);
    /* hier default file ausprobieren, falls nur gelesen werden soll! */
    if(*what == 'r') {
      name = (char*)
             Realloc(name,strlen(P->Files.default_path) + strlen(suffix)+1,"OpenFile");
      sprintf(name, "%s%s", P->Files.default_path,suffix);
      *file = fopen(name,what);
      if (*file != NULL) {
        if (verbose) fprintf(stderr,"(%i) Default file is open: %s\n",P->Par.me,name);
        Free(name, "OpenFile: name");
        return(1);
      } else {
        fprintf(stderr,"\n(%i)Error: Cannot open neither normal nor default file for reading: %s\n",P->Par.me,name);
        Free(name, "OpenFile: name");
        return(0);
      }
    } else {
      fprintf(stderr,"\n(%i)Error: Cannot open normal file for writing: %s\n",P->Par.me,name);
      Free(name, "OpenFile: name");
      return(0);
    }
  } else {
    if(verbose) fprintf(stderr,"(%i) File is open: %s\n",P->Par.me,name);
    Free(name, "OpenFile: name");
    return(1);
  }
}

/** Opens a file with FileData::mainname+suffix+"."+No (2 digits).
 * Both the path under which the main programme resides and the path to the config file are
 * tried subsequently.
 * \param *P Problem at hand
 * \param **file file pointer array
 * \param *suffix suffix after mainname
 * \param No the number with up to two digits
 * \param *what access parameter for the file: r, w, rw, r+, w+ ...
 * \param verbose 1 - print status, 0 - don't print anything
 * \return 0 - could not open file, 1 - file is open
 */
int OpenFileNo2(struct Problem *P, FILE** file, const char* suffix, int No, const char* what, int verbose)
{
  char* name; /* Zu erzeugender Dateiname */
  name = (char*)
    Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 4,"OpenFile");
  sprintf(name, "%s%s%s.%02i", P->Files.default_path, P->Files.mainname, suffix, No);
  *file = fopen(name, what);
  if (*file == NULL) {
    /* falls nur gelesen wird, auch default_path ausprobieren */
    if (*what == 'r') {
      name = (char*)
        Realloc(name,strlen(P->Files.default_path) + strlen(suffix) + 4,"OpenFileNo2");
      sprintf(name, "%s%s.%02i", P->Files.default_path, suffix, No);
      *file = fopen(name, what);
      if (*file != NULL) {
        if(verbose) fprintf(stderr,"(%i) Default file is open: %s\n", P->Par.me, name);
        Free(name, "OpenFileNo2: name");
        return(1);
      }
    }
    fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
    Free(name, "OpenFileNo2: name");
    return(0);
  } else {
    if(verbose) fprintf(stderr,"(%i) File is open: %s\n",P->Par.me, name);
    Free(name, "OpenFileNo2: name");
    return(1);
  }
}

/* Oeffnet Datei P->Files.mainname+"...".Nr(4stellig) mit what*/
/** Opens a file with FileData::mainname+suffix+"."+No (4 digits).
 * Both the path under which the main programme resides and the path to the config file are
 * tried subsequently.
 * \param *P Problem at hand
 * \param **file file pointer array
 * \param *suffix suffix after mainname
 * \param No the number with up to four digits
 * \param *what access parameter for the file: r, w, rw, r+, w+ ...
 * \param verbose 1 - print status, 0 - don't print anything
 * \return 0 - could not open file, 1 - file is open
 */
int OpenFileNo(struct Problem *P, FILE** file, const char* suffix, int No, const char* what, int verbose)
{
  char* name; /* Zu erzeugender Dateiname */
  name = (char*)
    Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 6,"OpenFileNo");
  sprintf(name, "%s%s%s.%04i", P->Files.default_path, P->Files.mainname, suffix, No);
  *file = fopen(name, what);
  if (*file == NULL) {
    /* falls nur gelesen wird, auch default_path ausprobieren */
    if (*what == 'r') {
      name = (char*)
        Realloc(name,strlen(P->Files.default_path) + strlen(suffix) + 6,"OpenFileNo");
      sprintf(name, "%s%s.%04i", P->Files.default_path, suffix, No);
      *file = fopen(name, what);
      if (*file != NULL) {
        if(verbose) fprintf(stderr,"(%i) Default file is open: %s\n", P->Par.me, name);
        Free(name, "OpenFileNo: name");
        return(1);
      }
    }
    fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
    Free(name, "OpenFileNo: name");
    return(0);
  } else {
    if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name);
    Free(name, "OpenFileNo: name");
    return(1);
  }
}

/* die nachfolgende Routine wird von seq und par benutzt!!! */
/* Oeffnet Datei P->Files.mainname+"...".No.postfix mit what*/

/** Opens a file with FileData::mainname+suffix+"."+No(4 digits)+postfix.
 * Only the path under which the main programme resides is tried.
 * \param *P Problem at hand
 * \param **file file pointer array
 * \param *suffix suffix after mainname
 * \param No the number with up to four digits
 * \param *postfix post-suffix at the end of filename
 * \param *what access parameter for the file: r, w, rw, r+, w+ ...
 * \param verbose 1 - print status, 0 - don't print anything
 * \return 0 - could not open file, 1 - file is open
 */
int OpenFileNoPost(struct Problem *P, FILE** file, const char* suffix, int No, const char* postfix, const char* what, int verbose)
{
  char* name; /* Zu erzeugender Dateiname */
  name = (char*)
    Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(postfix) + strlen(suffix) + 6,"OpenFileNoPost");
  sprintf(name, "%s%s%s.%04i%s", P->Files.default_path, P->Files.mainname, suffix, No, postfix);
  *file = fopen(name, what);
  if (*file == NULL) {
    fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
    Free(name, "OpenFileNoPost: name");
    return(0);
  } else {
    if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name);
    Free(name, "OpenFileNoPost: name");
    return(1);
  }
}

/** Opens a file with FileData::mainname+suffix+"."+No1(4 digits)+"."+No2(4 digits).
 * Only the path under which the main programme resides is tried.
 * \param *P Problem at hand
 * \param **file file pointer array
 * \param *suffix suffix after mainname
 * \param No1 first number with up to four digits
 * \param No2 second number with up to four digits
 * \param *what access parameter for the file: r, w, rw, r+, w+ ...
 * \param verbose 1 - print status, 0 - don't print anything
 * \return 0 - could not open file, 1 - file is open
 */
int OpenFileNoNo(struct Problem *P, FILE** file, const char* suffix, int No1, int No2, const char* what, int verbose)
{ /* zuerst suffix; No1: lfd, No2: procId */
  char *name;
  name = (char*)
    Malloc(strlen(P->Files.default_path) + strlen(P->Files.mainname) + strlen(suffix) + 5 + 6 + 1,"OpenFileNoNo");
  sprintf(name,"%s%s%s.%04i.%04i", P->Files.default_path, P->Files.mainname, suffix, No1, No2);
  *file = fopen(name, what);
  if(*file == NULL) {
    fprintf(stderr,"\n(%i)Error: Cannot open file: %s\n", P->Par.me, name);
    Free(name, "OpenFileNoNo: name");
    return(0);
  }
  if(verbose) fprintf(stderr,"(%i) File is open: %s\n", P->Par.me, name);
  Free(name, "OpenFileNoNo: name");
  return(1);
}

/** Transformation of wave function of Riemann level to zeroth before output as Vis.
 * \param *RT RiemannTensor structure, contains RiemannTensor::LatticeLevel
 * \param *source source wave function array
 * \param *dest destination wave function array
 * \param NF dimensional factor (NDIM, generally 3)
 */
static void OutVisPosRTransformPosNFRto0(const struct RiemannTensor *RT, fftw_real *source, fftw_real *dest, const int NF)
{
  struct LatticeLevel *Lev = RT->LevR;
  int es = Lev->NUp0[2]*NF;
  unsigned int cpyes = sizeof(fftw_real)*es;
  int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp0[0],Ny=Lev->NUp0[1];
  int lx,ly,z,Lx,Ly;
  for(lx=0; lx < nx; lx++)
    for(Lx=0; Lx < Nx; Lx++)
      for(ly=0; ly < ny; ly++)
        for(Ly=0; Ly < Ny; Ly++)
          for(z=0; z < nz; z++) {
            memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
                    &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
                    cpyes);
          }
}


/** prints Norm of each wave function to screen.
 * \param *out output destination (eg. stdout)
 * \param *P Problem at hand
 */
void OutputNorm (FILE *out, struct Problem *P) {
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct LatticeLevel *Lev = P->R.LevS;
  struct OnePsiElement *OnePsi, *LOnePsi;
  int i;
  
  for (i=0;i<Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT;i++) {
    OnePsi =&Psi->AllPsiStatus[i];
    if (OnePsi->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // local one
      LOnePsi = &Psi->LocalPsiStatus[OnePsi->MyLocalNo];
      fprintf(out,"(%i) Norm of Psi %i: %e\n", P->Par.me, OnePsi->MyGlobalNo, GramSchGetNorm2(P,Lev,Lev->LPsi->LocalPsi[OnePsi->MyLocalNo]));
    }
  }
}

/** Output of current Psis state of source level RunStruct::LevS's LevelPsi::LocalPsi to file \ref suffixsrcpsidat.
 * In case of process 0, the doc file is written in a parsable way with minimisation type RunStruct#CurrentMin, level number
 * LatticeLevel#LevelNo, and number of grid nodes LatticeLevel#N.
 * The two (three) different Psis::SpinType's are discerned and in case of SpinUpDown separate data files opened.
 *
 * Then for every global wave function of the desired type each coefficient of the reciprocal grid is copied into
 * Density::DensityCArray[TempDensity], fouriertransformed, copied into Density::DensityArray[TempDensity].
 * 
 * As the file output is only handled by process 0, all other coefficient shares are sent within the Psi group to which
 * the current global wave function belongs to process 0 of the Psi group and from there on to process 0 of all via
 * MPI.
 * \param *P Problem at hand
 * \param type Current minimisation state
 * \note This serves as a backup file, when the process is terminated and one would like to restart it at the
 *       current calculation lateron, see ReadSrcPsiDensity(). Note also that it is not necessary to specify the
 *       same number of processes on later restart, any number may be used under the condition that the number of
 *       grid nodes match and that there 2 sharing wave functions in case of SpinUpDown.
 * \sa ReadSrcPsiDensity() - same for Reading the coefficients, TestReadnWriteSrcDensity() - checks both routines against
 *  each other
 */
void OutputSrcPsiDensity(struct Problem *P, enum PsiTypeTag type) 
{
  int i,j,k, Index, zahl, owner;
  struct Lattice *Lat = &P->Lat;
  //struct RunStruct *R = &P->R;
  struct fft_plan_3d *plan = Lat->plan;
  struct LatticeLevel *LevS = P->R.LevS;
  struct Psis *Psi = &Lat->Psi;
  fftw_complex *work = (fftw_complex *)LevS->Dens->DensityArray[TempDensity];
  double *destpsiR = (double *)LevS->Dens->DensityArray[TempDensity];
  fftw_real *srcpsiR = (fftw_real *)LevS->Dens->DensityCArray[TempDensity];
  fftw_complex *srcpsiC = (fftw_complex *)LevS->Dens->DensityCArray[TempDensity];
  FILE *SrcPsiData, *SrcPsiDoc;
  char suffixdat[255], suffixdoc[255];
  MPI_Status status;
  struct OnePsiElement *OnePsiA, *LOnePsiA;
  fftw_complex *LPsiDatA;
  int Num = 0, colorNo = 0;
  int Sent = 0, sent = 0;
 
  SpeedMeasure(P,ReadnWriteTime,StartTimeDo);
  sprintf(suffixdat, ".%.254s.L%i", P->R.MinimisationName[type], LevS->LevelNo);
  strncpy (suffixdoc, suffixdat, 255);
  // for the various spin cases, output the doc-file if it's process 0
  if (!(P->Par.me_comm_ST)) {
    switch (Lat->Psi.PsiST) {
    case SpinDouble:
      colorNo = 0;
      strncat (suffixdat, suffixsrcpsidat, 255-strlen(suffixdat)); 
      strncat (suffixdoc, suffixsrcpsidoc, 255-strlen(suffixdoc));
      Num =  Lat->Psi.GlobalNo[PsiMaxNoDouble];
      break;
    case SpinUp:
      colorNo = 0;
      strncat (suffixdat, suffixsrcpsiupdat, 255-strlen(suffixdat)); 
      strncat (suffixdoc, suffixsrcpsiupdoc, 255-strlen(suffixdoc)); 
      Num = Lat->Psi.GlobalNo[PsiMaxNoUp];
      break;
    case SpinDown:
      colorNo = 1;
      strncat (suffixdat, suffixsrcpsidowndat, 255-strlen(suffixdat)); 
      strncat (suffixdoc, suffixsrcpsidowndoc, 255-strlen(suffixdoc)); 
      Num =  Lat->Psi.GlobalNo[PsiMaxNoDown];
      break;
    }
    if (!OpenFileNo(P, &SrcPsiData, suffixdat, colorNo, "wb",P->Call.out[ReadOut])) // open SourcePsiData as write binary
      fprintf(stderr,"(%i) Error opening file with suffix %s for writing!\n",P->Par.me, suffixdat);
    if (!(P->Par.my_color_comm_ST_Psi)) { // if we are process 0
        if (!OpenFile(P, &SrcPsiDoc, suffixdoc, "w",P->Call.out[ReadOut])) // open the (text) doc file
        fprintf(stderr,"(%i) Error opening file with suffix %s for writing!\n",P->Par.me, suffixdoc);
        fprintf(SrcPsiDoc, "Mintype\t%i\n", (int)type);
      fprintf(SrcPsiDoc, "LevelNo\t%i\n", LevS->LevelNo);
      fprintf(SrcPsiDoc, "GridNodes\t%i\t%i\t%i\n", LevS->N[0], LevS->N[1], LevS->N[2]);
      fprintf(SrcPsiDoc, "PsiNo\t%i\t%i\n", Num, P->Lat.Psi.GlobalNo[PsiMaxAdd]);
      fprintf(SrcPsiDoc, "Epsilon\t%lg\t%lg\n", P->R.RelEpsTotalEnergy, P->R.RelEpsKineticEnergy);
        for (i = 0; i < P->Par.Max_my_color_comm_ST_Psi; i++) {
          fprintf(SrcPsiDoc, "\t%i", Lat->Psi.RealAllLocalNo[i]); // print number of local Psis
      }
      fprintf(SrcPsiDoc, "\n");
        fclose(SrcPsiDoc);
    }
  }
  
  // send/receive around and write share of coefficient array of each wave function 
  MPI_Allreduce(&sent, &Sent, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); // catch all at the starter line
  if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) me (%i/%i) \t Psi (%i/%i)\t PsiT (%i/%i)\n", P->Par.me, P->Par.me_comm_ST, P->Par.Max_me_comm_ST, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi, P->Par.me_comm_ST_PsiT, P->Par.Max_me_comm_ST_PsiT);
  k = -1;  // k is global PsiNo counter for the desired group
  for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) {  // go through all wave functions (plus the extra one for each process)
  OnePsiA = &Psi->AllPsiStatus[j];    // grab OnePsiA
    if (OnePsiA->PsiType == type) { // only take desired minimisation group
      k++;
      owner = 0; // notes down in process 0 of each psi group the owner of the next coefficient share
      //fprintf(stderr,"(%i) ST_Psi: OnePsiA %i\tP->Par.me %i\n", P->Par.me,OnePsiA->my_color_comm_ST_Psi,P->Par.my_color_comm_ST_Psi);
      if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // Belongs to my Psi group?
        LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
        LPsiDatA = LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo];
        SetArrayToDouble0((double *)srcpsiR, LevS->Dens->TotalSize*2);   // zero DensityCArray[TempDensity]
        for (i=0;i<LevS->MaxG;i++) {    // for each every unique G grid vector
          Index = LevS->GArray[i].Index;
          srcpsiC[Index].re = LPsiDatA[i].re;  // copy real value
          srcpsiC[Index].im = LPsiDatA[i].im;  // copy imaginary value
        }
        for (i=0; i<LevS->MaxDoubleG; i++) {  // also for every doubly appearing G vector (symmetry savings)
          srcpsiC[LevS->DoubleG[2*i+1]].re = srcpsiC[LevS->DoubleG[2*i]].re;
          srcpsiC[LevS->DoubleG[2*i+1]].im = -srcpsiC[LevS->DoubleG[2*i]].im;
        }
        // do an fft transform from complex to real on these srcPsiR    
        fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNF1, srcpsiC, work);
        
        for (i=0; i < LevS->Dens->LocalSizeR; i++)
          destpsiR[i] = (double)srcpsiR[i];
      } else 
         LOnePsiA = NULL;
  
      if (P->Par.me_comm_ST == 0) { // if we are process 0 of all, only we may access the file
        for (i=0; i<P->Par.Max_me_comm_ST_Psi;i++) { // for each share of the coefficient in the PsiGroup
          if (LOnePsiA == NULL) {   // if it's not local, receive coefficients from correct PsiGroup (process 0 within that group)
            if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, OnePsiA->my_color_comm_ST_Psi, OutputSrcPsiTag, P->Par.comm_ST_PsiT, &status ) != MPI_SUCCESS)
              Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
            MPI_Get_count(&status, MPI_DOUBLE, &zahl);
            if (zahl != LevS->Dens->LocalSizeR) // check number of elements
              fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, i, zahl, LevS->Dens->LocalSizeR);
            //else        
              //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, i);
          } else { // if it's local ...
            if (i != 0) { // but share of array not for us, receive from owner process within Psi group
              if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, i, OutputSrcPsiTag, P->Par.comm_ST_Psi, &status ) != MPI_SUCCESS)
                Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
              MPI_Get_count(&status, MPI_DOUBLE, &zahl);
              if (zahl != LevS->Dens->LocalSizeR) // check number of elements
                fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, i, zahl, LevS->Dens->LocalSizeR);
              //else        
                //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, i);
            } // otherwise it was our share already 
          }
          // store the final share on disc
          if ((zahl = fwrite(destpsiR, sizeof(double), (size_t)(LevS->Dens->LocalSizeR), SrcPsiData)) != (size_t)(LevS->Dens->LocalSizeR)) {
            fclose(SrcPsiData);
            //if (P->Par.me == 0)
              fprintf(stderr, "(%i)OutputSrcPsiDensity: only %i bytes written instead of expected %i\n", P->Par.me, zahl, LevS->Dens->LocalSizeR);
            Error(SomeError,"OutputSrcPsiDensity: fwrite Error");
          }
        }
      } else { // if we are not process 0  of all, we are but a deliverer
        if (LOnePsiA != NULL) { // send if it's local
          if (P->Par.me_comm_ST_Psi == 0) { // if we are process 0 in the group, send final share to our process 0
            for (owner = 0; owner < P->Par.Max_me_comm_ST_Psi; owner++) { // for all processes of our Psi group             
              if (owner != 0) { // still not our share of coefficients, receive from owner in our Psi group (increasing owner)
                if (MPI_Recv( destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, owner, OutputSrcPsiTag, P->Par.comm_ST_Psi, &status ) != MPI_SUCCESS)
                  Error(SomeError, "OutputSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
                MPI_Get_count(&status, MPI_DOUBLE, &zahl);
                if (zahl != LevS->Dens->LocalSizeR) // check number of elements
                  fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, owner, zahl, LevS->Dens->LocalSizeR);
                //else        
                  //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, owner);
              } else sent++; // only count sent if it was our share 
              
              if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, OutputSrcPsiTag, P->Par.comm_ST_PsiT) != MPI_SUCCESS)
                Error(SomeError, "OutputSrcPsiDensity: MPI_Send of loaded coefficients failed!");
              //else 
                //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Send to process %i in PsiT group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, 0, k);
            }
          } else {
            sent++;
            if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, OutputSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
              Error(SomeError, "OutputSrcPsiDensity: MPI_Send of loaded coefficients failed!");
            //else 
              //fprintf(stderr,"(%i)OutputSrcPsiDensity: MPI_Send to process %i in Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, 0, k);
          }
        }
        // otherwise we don't have anything to do with this
      }
    }
  }
  MPI_Allreduce(&sent, &Sent, 1, MPI_INT, MPI_SUM, P->Par.comm_ST); // catch all again at finish
  if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Out of %i shares %i had to be sent in total, %i from this process alone.\n", P->Par.me, P->Par.Max_me_comm_ST_Psi*Psi->NoOfPsis, Sent, sent);
  if (!(P->Par.me_comm_ST))
    fclose(SrcPsiData);
  SpeedMeasure(P,ReadnWriteTime,StopTimeDo);
}

/** Tests whether writing and successant reading of coefficient array is working correctly.
 * The local wave function array is written to a disc (\sa OutputSrcPsiDensity()), the by \a wavenr
 * specified coefficient array copied to OldPsiDat and afterwards read again via ReadSrcPsiDensity().
 * Comparison per process of each local coefficient shows incorrect read or writes.
 * \param *P Problem at hand
 * \param type minimisation type array to test for read and write
 * \return 1 - successful, 0 - test failed 
 */
int TestReadnWriteSrcDensity(struct Problem *P, enum PsiTypeTag type)
{
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  int i,k;
  fftw_complex *destpsiC, *srcpsiC;
  
  //fprintf(stderr,"(%i)TestReadnWriteSrcDensity\n",P->Par.me); 
  // write whole array of type to disc
  OutputSrcPsiDensity(P,type);
  debug(P,"array written");
  
  // copy specified array to OldPsiDat
  for (k=0; k < Lat->Psi.LocalNo; k++)  // for every local wave function of type, copy coefficients
    if (Psi->LocalPsiStatus[k].PsiType == type) {  // ... yet only for given type
      srcpsiC = LevS->LPsi->LocalPsi[k];
      destpsiC = LevS->LPsi->OldLocalPsi[k - Psi->TypeStartIndex[type]];
      for (i=0;i<LevS->MaxG;i++) {    // for each every unique G grid vector
        destpsiC[i].re = srcpsiC[i].re;  // copy real value
        destpsiC[i].im = srcpsiC[i].im;  // copy imaginary value
      }
    }
  debug(P,"array copied");
    
  // read whole array again
  if (!ReadSrcPsiDensity(P,type,0,R->LevSNo))
    return 0;
  debug(P,"array read");
  
  // compare with copied array
  for (k=0; k < Lat->Psi.LocalNo; k++)  // for every local wave function of type, compare coefficients
    if (Psi->LocalPsiStatus[k].PsiType == type) {  // ... yet only for given type
      srcpsiC = LevS->LPsi->LocalPsi[k];
      destpsiC = LevS->LPsi->OldLocalPsi[k - Psi->TypeStartIndex[type]];
      for (i=0;i<LevS->MaxG;i++) // for each every unique G grid vector
        if ((fabs(destpsiC[i].re - srcpsiC[i].re) >= MYEPSILON) ||(fabs(destpsiC[i].im - srcpsiC[i].im) >= MYEPSILON)) {
          fprintf(stderr,"(%i)TestReadnWriteSrcDensity: First difference at index %i - %lg+i%lg against loaded %lg+i%lg\n",P->Par.me, i, srcpsiC[i].re, srcpsiC[i].im,destpsiC[i].re,destpsiC[i].im); 
          return 0; 
        }
    }
  debug(P,"array compared");
  fprintf(stderr,"(%i)TestReadnWriteSrcDensity: OK!\n",P->Par.me); 
  return 1;
}


/** Read Psis state from an earlier run.
 * The doc file is opened, mesh sizes LatticeLevel::N[], global number of Psis read and checked against the known
 * values in the inital level RunStruct::LevS. 
 * Note, only process 0 handles the files, all read coefficients are transfered to their respective owners via MPI
 * afterwards. Here, process 0 of a certain Psi group is used as a transferer for the coefficients of the other processes
 * in this Psi group. He receives them all from process 0 and sends them onward accordingly. The complete set of
 * coefficients on the real grid for one wave function in the Psi group are transformed into complex wave function
 * coefficients by the usual fft procedure (see ChangeToLevUp()).
 * 
 * \param *P Problem at hand
 * \param type minimisation type to read
 * \param test whether to just test for presence of files (1) or really read them (0)
 * \param LevSNo level number to be read
 * \note This is the counterpart to OutputSrcPsiDensity().
 * \return 1 - success, 0 - failure
 * \note It is not necessary to specify the same number of processes on later restart, any number may be used under
 *       the condition that the number of grid nodes match and that there at least 2 processes sharing wave functions
 *       in case of SpinUpDown.
 * \sa OutputSrcPsiDensity() - same for writing the coefficients, TestReadnWriteSrcDensity() - checks both routines against
 *  each other
 */
int ReadSrcPsiDensity(struct Problem *P, enum PsiTypeTag type, int test, int LevSNo)
{
  int i, j, k, Index, owner;
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct fft_plan_3d *plan = Lat->plan;
  struct LatticeLevel *LevS = R->LevS;  // keep open for LevelNo read from file
  struct Psis *Psi = &Lat->Psi;
  //struct Energy *E = Lat->E;
  fftw_complex *work;
  double *destpsiR;
  fftw_real *srcpsiR;
  fftw_complex *srcpsiC;
  FILE *SrcPsiData, *SrcPsiDoc;
  int N[NDIM], GlobalNo[2];
  int LevelNo, readnr=0;
  int zahl, signal = test ? 1 : 2; // 0 - ok, 1 - test failed, 2 - throw Error
  char suffixdat[255], suffixdoc[255];
  int read_type, Num = 0, colorNo = 0;
  char spin[20];
  double Eps[2];
  MPI_Status status;
  struct OnePsiElement *OnePsiA, *LOnePsiA;
  int Recv=0, recv=0;
  
  SpeedMeasure(P,ReadnWriteTime,StartTimeDo);
  sprintf(suffixdat, ".%.254s.L%i", P->R.MinimisationName[type], LevSNo);
  strncpy (suffixdoc, suffixdat, 255);
  // Depending on Psis::SpinType the source psi doc file is opened and header written
  switch (Lat->Psi.PsiST) {
  case SpinDouble:    
    colorNo = 0;
    strncat (suffixdat, suffixsrcpsidat, 255-strlen(suffixdat)); 
    strncat (suffixdoc, suffixsrcpsidoc, 255-strlen(suffixdoc));
    strncpy (spin, "GlobalNoSpinDouble", 20); 
    Num = Lat->Psi.GlobalNo[PsiMaxNoDouble];
    break;
  case SpinUp:
    colorNo = 0;
    strncat (suffixdat, suffixsrcpsiupdat, 255-strlen(suffixdat)); 
    strncat (suffixdoc, suffixsrcpsiupdoc, 255-strlen(suffixdoc)); 
    strncpy (spin, "GlobalNoSpinUp", 20); 
    Num = Lat->Psi.GlobalNo[PsiMaxNoUp];
    break;
  case SpinDown:
    colorNo = 1;
    strncat (suffixdat, suffixsrcpsidowndat, 255-strlen(suffixdat)); 
    strncat (suffixdoc, suffixsrcpsidowndoc, 255-strlen(suffixdoc)); 
    strncpy (spin, "GlobalNoSpinDown", 20); 
    Num = Lat->Psi.GlobalNo[PsiMaxNoDown];
    break;
  }
  // open doc file ...
  if (!(P->Par.me_comm_ST)) {
    if (!OpenFile(P, &SrcPsiDoc, suffixdoc, "r", test ? 0 : P->Call.out[ReadOut])) { // open doc file
      debug(P,"ReadSrcPsiDensity: doc file pointer NULL\n");
      if (test) {
        if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
          Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
        return 0;
      }
      if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
        Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
      Error(SomeError,"ReadSrcPsiDensity: cannot open doc file!");
    }
    // ... and parse critical ...
    readnr += ParseForParameter(0,SrcPsiDoc,"Mintype",0,1,1,int_type,(int *)&read_type, 1, test ? optional : critical);
    readnr += ParseForParameter(0,SrcPsiDoc,"LevelNo",0,1,1,int_type,&LevelNo,1,  test ? optional : critical);
    readnr += 3*ParseForParameter(0,SrcPsiDoc,"GridNodes",0,3,1,row_int,&N[0], 1, test ? optional : critical);
    readnr += 2*ParseForParameter(0,SrcPsiDoc,"PsiNo",0,2,1,row_int,&GlobalNo[0], 1, test ? optional : critical);
    // and optional items ...
    if (ParseForParameter(0,SrcPsiDoc,"Epsilon",0,2,1,row_double,&Eps[0],1,optional))
      if ((P->Call.ReadSrcFiles == 1) && ((Eps[1] < R->RelEpsKineticEnergy) || (Eps[0] < R->RelEpsTotalEnergy))) {
        //fprintf(stderr,"(%i) Eps %lg %lg\tRelEps %lg %lg\n", P->Par.me, Eps[0], Eps[1], R->RelEpsTotalEnergy, R->RelEpsKineticEnergy);
        fprintf(stderr,"(%i) NOTE: Doing minimization after source file parsing due to smaller specified epsilon stop conditions.\n",P->Par.me); 
        P->Call.ReadSrcFiles = 2; // do minimisation even after loading
      }
    if (readnr != 7) { // check number of items read
      debug(P, "ReadSrcPsiDensity: too few doc items in file\n");
      fclose(SrcPsiDoc);
      if (test) {
        if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
          Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
        return 0;
      }
      fprintf(stderr,"ReadSrcPsiDensity: Only %i items read!\n",readnr);
      if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
        Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
      Error(SomeError, "ReadSrcPsiDensity: read error");
    }
    // check if levels match
    if (LevSNo != LevelNo) { 
      debug(P,"ReadSrcPsiDensity: mismatching levels\n");
      fclose(SrcPsiDoc);
      if (test) {
        if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
          Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
        return 0;
      }
      if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
        Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
      Error(SomeError,"ReadSrcPsiDensity: Mismatching levels!");
    } else {
      LevS = &P->Lat.Lev[LevelNo];
    }
    
    // check if systems in memory and file match
    if ((read_type != R->CurrentMin) || (N[0] !=  LevS->N[0]) || (N[1] !=  LevS->N[1]) || (N[2] !=  LevS->N[2]) || (GlobalNo[0] != Num) || (GlobalNo[1] != Lat->Psi.GlobalNo[PsiMaxAdd])) { // check read system
      debug(P,"ReadSrcPsiDensity: srcpsi file does not fit to system\n");
      fclose(SrcPsiDoc);
      if (test) {
        fprintf(stderr,"(%i) Min\t N(x,y,z)\tPsiNo+AddNo\n file: %s\t %i %i %i\t %i + %i\nsystem: %s\t %d %d %d\t %d + %d\n",P->Par.me, R->MinimisationName[read_type], N[0], N[1], N[2], GlobalNo[0], GlobalNo[1], R->MinimisationName[R->CurrentMin], LevS->N[0] , LevS->N[1], LevS->N[2], Lat->Psi.GlobalNo[PsiMaxNoDouble], Lat->Psi.GlobalNo[PsiMaxAdd]);
        if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
          Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
        return 0;
      }
      fprintf(stderr,"ReadSrcPsiDensity: Type %i != CurrentMin %i || N[0] %i != %i || N[1] %i != %i || N[2] %i != %i || %s %i + %i != %i + %i\n", read_type, R->CurrentMin, N[0], LevS->N[0], N[1], LevS->N[1], N[2], LevS->N[2], spin, GlobalNo[0], GlobalNo[1], Num, P->Lat.Psi.GlobalNo[PsiMaxAdd]);
      if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
        Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
      Error(SomeError,"ReadSrcPsiDensity: srcpsi file does not fit to system");
    }
    signal = 0; // everything went alright, signal ok
    if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
      Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
  } else {  // others wait for signal from root process
    if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
      Error(SomeError,"ReadSrcPsiDensity: Bcast of signal failed\n");
    if (signal == 1)
      return 0;
    else if (signal == 2)
      Error(SomeError, "ReadSrcPsiDensity: Something went utterly wrong, see root process");
    else if (P->Call.out[PsiOut])
      fprintf(stderr,"(%i) ReadSrcPsiDensity: Everything went alright so far\n", P->Par.me);
  }
  if (MPI_Bcast(&LevelNo,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
    Error(SomeError,"ReadSrcPsiDensity: Bcast of LevelNo failed\n");
  LevS = &P->Lat.Lev[LevelNo];
  //if (!test) fprintf(stderr,"(%i) LevelSNo %i\n", P->Par.me, LevS->LevelNo);

  if (!test) {
    // set some pointers for work to follow
    work = (fftw_complex *)LevS->Dens->DensityArray[TempDensity];
    destpsiR = (double *)LevS->Dens->DensityArray[TempDensity];
    srcpsiR = (fftw_real *)LevS->Dens->DensityCArray[TempDensity];
    srcpsiC = (fftw_complex *)LevS->Dens->DensityCArray[TempDensity];
  
    // read share of coefficient array for each wave function and send/receive around
    owner = 0;
    MPI_Allreduce (&recv, &Recv, 1, MPI_INT, MPI_SUM, P->Par.comm_ST);
    //fprintf(stderr,"(%i) me (%i/%i) \t Psi (%i/%i)\t PsiT (%i/%i)\n", P->Par.me, P->Par.me_comm_ST, P->Par.Max_me_comm_ST, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi, P->Par.me_comm_ST_PsiT, P->Par.Max_me_comm_ST_PsiT);
    k = -1;  // k is global PsiNo counter for the desired group
    for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) {  // go through all wave functions (plus the extra one for each process)
      OnePsiA = &Psi->AllPsiStatus[j];    // grab OnePsiA
      if (OnePsiA->PsiType == type) { // only take desired minimisation group
        k++;
        //fprintf(stderr,"(%i) ST_Psi: OnePsiA %i\tP->Par.my_color_comm_ST_Psi %i\n", P->Par.me,OnePsiA->my_color_comm_ST_Psi,P->Par.my_color_comm_ST_Psi);
        if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // Belongs to my Psi group?
           LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
        else 
           LOnePsiA = NULL;
    
        if (P->Par.me_comm_ST == 0) { // if we are process 0 of all, we may access file
          if (!OpenFileNo(P, &SrcPsiData, suffixdat, colorNo, "r", test ? 0 : P->Call.out[ReadOut])) {
            Error(SomeError,"ReadSrcPsiDensity: cannot open data file");
          }        
          for (i=P->Par.Max_me_comm_ST_Psi-1; i>=0;i--) { // load coefficients
            fseek( SrcPsiData, (long)((k*N[0]*N[1]*N[2]+i*((long)LevS->Dens->LocalSizeR))*sizeof(double)), SEEK_SET); // seek to beginning of this process' coefficients
            // readin 
            if ((zahl = fread(destpsiR, sizeof(double), (size_t)(LevS->Dens->LocalSizeR), SrcPsiData)) != (size_t)LevS->Dens->LocalSizeR) {
              fclose(SrcPsiData);
              fprintf(stderr, "(%i)ReadSrcPsiDensity: only %i bytes read instead of expected %i\n", P->Par.me, zahl, LevS->Dens->LocalSizeR);
              Error(SomeError,"ReadSrcPsiDensity: fread Error");
            }
            if (LOnePsiA == NULL) {   // if it's not local, send away coefficients to correct PsiGroup (process 0 within that group)
              if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, OnePsiA->my_color_comm_ST_Psi, ReadSrcPsiTag, P->Par.comm_ST_PsiT) != MPI_SUCCESS)
                Error(SomeError, "ReadSrcPsiDensity: MPI_Send of loaded coefficients failed!");
              //else 
                //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i of loaded coefficients GlobalNo %i, owner %i succeeded!\n", P->Par.me, OnePsiA->my_color_comm_ST_Psi, k, i); 
            } else { // if it's local ...
              if (i != 0) { // but share of array not for us, send to owner process within Psi group
                if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, i, ReadSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
                  Error(SomeError, "ReadSrcPsiDensity: MPI_Send within Psi group of loaded coefficients failed!");
                //else 
                  //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i within Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, i, k);
              } // otherwise it was our share already 
            }
          }
        } else {
          if (LOnePsiA != NULL) { // receive
            if (P->Par.me_comm_ST_Psi == 0) { // if we are process 0 in the group, receive share from process 0 of all
              for (owner = P->Par.Max_me_comm_ST_Psi -1; owner >=0; owner--) { // for all processes of our Psi group
                if (MPI_Recv(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, ReadSrcPsiTag, P->Par.comm_ST_PsiT, &status) != MPI_SUCCESS)
                  Error(SomeError, "ReadSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
                MPI_Get_count(&status, MPI_DOUBLE, &zahl);
                if (zahl != LevS->Dens->LocalSizeR) // check number of elements
                  fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv from process 0 of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, owner, zahl, LevS->Dens->LocalSizeR);
                //else        
                  //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv from process 0 of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, owner);
                
                if (owner != 0) { // not our share of coefficients, send to owner in our Psi group
                  if (MPI_Send(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, owner, ReadSrcPsiTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
                    Error(SomeError, "ReadSrcPsiDensity: MPI_Send within Psi group of loaded coefficients failed!");
                  //else 
                    //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Send to process %i within Psi group of loaded coefficients GlobalNo %i succeeded!\n", P->Par.me, owner, k);
                } else recv++;
              }
              // otherwise it's our share!
            } else {  // our share within Psi Group not belonging to process 0 of all
              recv++;
              if (MPI_Recv(destpsiR, LevS->Dens->LocalSizeR, MPI_DOUBLE, 0, ReadSrcPsiTag, P->Par.comm_ST_Psi, &status) != MPI_SUCCESS)
                Error(SomeError, "ReadSrcPsiDensity: MPI_Recv of loaded coefficients failed!");
              MPI_Get_count(&status, MPI_DOUBLE, &zahl);
              if (zahl != LevS->Dens->LocalSizeR) // check number of elements
                fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i failed: Too few coefficients - %i instead of %i!\n", P->Par.me, k, P->Par.me_comm_ST_Psi, zahl, LevS->Dens->LocalSizeR);
              //else        
                //fprintf(stderr,"(%i)ReadSrcPsiDensity: MPI_Recv of loaded coefficients of GlobalNo %i, owner %i succeeded!\n", P->Par.me, k, P->Par.me_comm_ST_Psi);
            } 
          }
          // otherwise we don't have anything to do with this
        }
        
        if (LOnePsiA != NULL) {
          SetArrayToDouble0((double *)srcpsiR, LevS->Dens->TotalSize*2);
          for (i=0; i < LevS->Dens->LocalSizeR; i++)  // copy dest to src
            srcpsiR[i] = (fftw_real)destpsiR[i];
    
          fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, srcpsiC, work); // fft transform
          //if (P->Call.out[PsiOut]) 
          //fprintf(stderr,"(%i) LevSNo %i\t LocalPsi %p\n", P->Par.me, LevS->LevelNo, LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo]);
          for (i=0;i<LevS->MaxG;i++) {  // and copy into wave functions coefficients
            Index = LevS->GArray[i].Index;
            LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo][i].re = srcpsiC[Index].re/LevS->MaxN;
            LevS->LPsi->LocalPsi[LOnePsiA->MyLocalNo][i].im = srcpsiC[Index].im/LevS->MaxN;
          }
        }
        if ((P->Par.me_comm_ST == 0) && (SrcPsiData != NULL)) fclose(SrcPsiData);
      }
    }
    MPI_Allreduce (&recv, &Recv, 1, MPI_INT, MPI_SUM, P->Par.comm_ST);
    if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Out of %i shares %i had to be received in total, %i from this process alone.\n", P->Par.me, P->Par.Max_me_comm_ST_Psi*Psi->NoOfPsis, Recv, recv);
    SpeedMeasure(P,ReadnWriteTime,StopTimeDo);
  }
  return 1; // everything went well till the end
}

/** Creates the density \ref suffixdensdoc and \ref suffixdensdx files for OpenDx.
 * Opens \ref suffixdensdoc, fills (pos&data file name, byte order, max mesh points, matrix alignment, steps)
 * and closes it.
 * Opens \ref suffixdensdx, then for every FileData::OutVisStep a describing structure for DX is written and
 * the file closed again.
 * \param *P Problem at hand
 * \param me my process number in communicator Psi (0 - do nothing, else - do)
 */
static void CreateDensityOutputGeneral(struct Problem *P, const int me)
{
  FILE *DensityDoc, *DensityDx;
  char *posname, *datname; 
  struct LatticeLevel *Lev = &P->Lat.Lev[STANDARTLEVEL];
  unsigned int i, MaxPoints, N[NDIM];
  double *RB = P->Lat.RealBasis;
  if (me) return;
  N[0] = Lev->N[0]*Lev->NUp[0];
  N[1] = Lev->N[1]*Lev->NUp[1];
  N[2] = Lev->N[2]*Lev->NUp[2];
  MaxPoints = (N[0]+1)*(N[1]+1)*(N[2]+1);
  posname = (char*)
    Malloc(strlen(P->Files.mainname) + strlen(suffixdenspos) + 1,"OpenFile");
  sprintf(posname, "%s%s", P->Files.mainname, suffixdenspos);
  datname = (char*)
    Malloc(strlen(P->Files.mainname) + strlen(suffixdensdat) + 1,"OpenFile");
  sprintf(datname, "%s%s", P->Files.mainname, suffixdensdat);
  // write doc file
  OpenFile(P, &DensityDoc, suffixdensdoc, "w",P->Call.out[ReadOut]);
  fprintf(DensityDoc,"DensityPositions file = %s.####\n", posname);
  fprintf(DensityDoc,"DensityData file = %s.####\n", datname);
  fprintf(DensityDoc,"format = ieee float (Bytes %lu) %s\n",(unsigned long) sizeof(float),msb);
  fprintf(DensityDoc,"points = %i\n", MaxPoints);
  fprintf(DensityDoc,"majority = row\n");
  fprintf(DensityDoc,"TimeSeries = %i\n",P->Files.OutVisStep+1);
  fclose(DensityDoc);
  // write DX file
  OpenFile(P, &DensityDx, suffixdensdx, "w",P->Call.out[ReadOut]);
  for (i=0; i < (unsigned int)P->Files.OutVisStep+1; i++) { // for every OutVis step
    if (i==0) { 
      fprintf(DensityDx,"object \"gridcon\" class gridconnections counts %i %i %i\n\n",(N[0]+1),(N[1]+1),(N[2]+1));
      if (P->Files.OutputPosType[i] != active)
        fprintf(DensityDx, "object \"posdens\" class gridpositions counts %i %i %i\norigin 0.0  0.0  0.0\ndelta  %f  %f  %f\ndelta  %f  %f  %f\ndelta  %f  %f  %f\n\n",
                (N[0]+1),(N[1]+1),(N[2]+1),
                (float)(RB[0]/N[0]),(float)(RB[1]/N[0]),(float)RB[2]/N[0],
                (float)(RB[3]/N[1]),(float)(RB[4]/N[1]),(float)RB[5]/N[1],
                (float)(RB[6]/N[2]),(float)(RB[7]/N[2]),(float)RB[8]/N[2]);
      }
      if (P->Files.OutputPosType[i] == active) {
        fprintf(DensityDx,
              "object \"pos.%04u\" class array type float rank 1 shape 3 items %i %s binary\n",i,MaxPoints,msb); 
        fprintf(DensityDx,"data file %s.%04u,0\n",posname,i);
      }
      fprintf(DensityDx,"attribute \"dep\" string \"positions\"\n");
      fprintf(DensityDx,"# %lu - %lu Bytes\n\n",MaxPoints*i*(unsigned long)sizeof(float)*NDIM,MaxPoints*(i+1)*(unsigned long)sizeof(float)*NDIM-1);
  
      fprintf(DensityDx,"object \"dat.%04u\" class array type float rank 0 items %i %s binary\n",i,MaxPoints,msb);
      fprintf(DensityDx,"data file %s.%04u,0\n",datname,i);
      fprintf(DensityDx,"attribute \"dep\" string \"positions\"\n");
      fprintf(DensityDx,"# %lu - %lu Bytes\n\n",MaxPoints*i*(unsigned long)sizeof(float),MaxPoints*(i+1)*(unsigned long)sizeof(float)-1);
  
      fprintf(DensityDx,"object \"obj.%04u\" class field\n",i);
      if (P->Files.OutputPosType[i] == active)
        fprintf(DensityDx,"component \"positions\" \"pos.%04i\"\n",i);
      if (P->Files.OutputPosType[i] != active)
        fprintf(DensityDx,"component \"positions\" \"posdens\"\n");
      fprintf(DensityDx,"component \"connections\" \"gridcon\"\n");
      fprintf(DensityDx,"component \"data\" \"dat.%04i\"\n",i);
    }
  fprintf(DensityDx,"\nobject \"series\" class series\n");
  for (i=0; i < (unsigned int)P->Files.OutVisStep+1; i++) 
    fprintf(DensityDx,"member %i \"obj.%04u\" position %f\n",i,i,(float)i);
  fprintf(DensityDx,"end\n");

  fclose(DensityDx);
  Free(posname, "CreateDensityOutputGeneral: posname");
  Free(datname, "CreateDensityOutputGeneral: datname");
}

/** Calculates the OutVis density of the RiemannTensor level.
 * The usual pattern arises when a density is fftransformed:
 * -# over all grid vectors up to MaxG
 * -# over all doubly grid vectors up to MaxDoubleG
 * -# call to fft_3d_complex_to_real()
 * 
 * In this case here followed by call to OutVisPosRTransformPosNFRto0() and finally FileData::work
 * is copied to FileData::PosR.
 * \param *Lat Lattice structure, containing Lattice::plan and LatticeLevel
 * \param *F FileData structure, containing FileData::PosC, FileData::PosR, FileData::work, FileData::Totalsize, FileData::LocalSizeR
 */
static void CalculateOutVisDensityPos(struct Lattice *Lat, struct FileData *F/*, const double FactorC_R, const double FactorR_C*/) 
{
  struct  fft_plan_3d *plan = Lat->plan;
  struct RiemannTensor *RT = &Lat->RT;
  struct LatticeLevel *LevR = RT->LevR;
  fftw_complex *destC = F->PosC;
  fftw_real *destR = F->PosR;
  fftw_complex *work = F->work;
  fftw_real *workR = (fftw_real*)work;
  fftw_complex *PosFactor = F->PosFactor;
  fftw_complex *posfac, *destpos, *destRCS, *destRCD;
  fftw_complex *coeff;
  fftw_complex source;
  int i, Index, pos, n;
  int NF = NDIM, MaxNUp = F->MaxNUp, TotalSize = F->TotalSize, LocalSizeR = F->LocalSizeR;
  SetArrayToDouble0((double *)destC, TotalSize*2);
  for (i=0;i < LevR->MaxG;i++) {
    Index = LevR->GArray[i].Index;
    posfac = &PosFactor[MaxNUp*i];
    destpos = &destC[MaxNUp*Index*NF];
    coeff = &RT->Coeff[i];
    for (pos=0; pos < MaxNUp; pos++) {
      for (n=0; n < NF; n++) {
        source.re = coeff[n].re;
        source.im = coeff[n].im;
        destpos[n+NF*pos].re = source.re*posfac[pos].re-source.im*posfac[pos].im;
              destpos[n+NF*pos].im = source.re*posfac[pos].im+source.im*posfac[pos].re;
      }
    }
  }
  for (i=0; i < LevR->MaxDoubleG; i++) {
    destRCS = &destC[LevR->DoubleG[2*i]*MaxNUp*NF];
    destRCD = &destC[LevR->DoubleG[2*i+1]*MaxNUp*NF];
    for (pos=0; pos < MaxNUp; pos++) {
      for (n=0; n < NF; n++) {
        destRCD[n+NF*pos].re =  destRCS[n+NF*pos].re;
        destRCD[n+NF*pos].im = -destRCS[n+NF*pos].im;
      }
    }
  }
  fft_3d_complex_to_real(plan, LevR->LevelNo, FFTNFRVecUp0, destC, work);
  OutVisPosRTransformPosNFRto0(RT, destR, workR, NF);
  memcpy(destR,workR,sizeof(fftw_real)*LocalSizeR);
}

/** Prepare Density::DensityArray for output.
 * Into FileData::work subsequently each node (all z, all y, all x) is written as \f$\log(1+x)\f$,
 * where x is Density::DensityArray[TotalDensity]. In the end result is send to process 0 (yet not
 * received here, see CombineOutVisArray()). In  case of RiemannTensor use, some more complex calculations
 * are made: FileData::PosR is used, the coefficient offset'ed to the current node and the log taken there.
 * \param *P Problem at hand
 * \param myPE this ranks process in the Psi communcator ParallelSimulationData::me_comm_ST_Psi
 * \param *srcdens Pointer to DensityArray which is to be displayed
 */
static void OutputOutVisDensity(struct Problem *P, const int myPE, fftw_real *srcdens) 
{
  int N[NDIM], n[NDIM], pos[NDIM];
  int destpos = 0;
  double fac[NDIM], posd[NDIM];
  float posf[NDIM+1];
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *Lev0 = &Lat->Lev[0];
  fftw_real *srcpos = P->Files.PosR;
  //fftw_real *srcdens = Lev0->Dens->DensityArray[ActualDensity]; //[TotalDensity]; trick to display single density
  float *dest = (float *)P->Files.work;
  int Nx = Lev0->Plan0.plan->N[0];
  int i;
  double min, max;

  N[0] = Lev0->Plan0.plan->local_nx;
  N[1] = Lev0->Plan0.plan->N[1];
  N[2] = Lev0->Plan0.plan->N[2];
  
  max = min = srcdens[0];
  for (i=1;i<P->R.Lev0->Dens->LocalSizeR;i++) {
    if (srcdens[i] < min) min = srcdens[i];
    if (srcdens[i] > max) max = srcdens[i];
  }
  if (P->Call.out[PsiOut]) fprintf(stderr,"(%i)OutputOutVisDensity: min %e\tmax %e\n",P->Par.me, min, max);
  
  // go through all nodes
  for (n[0]=0; n[0] < N[0]; n[0]++) {
    pos[0] = (n[0] == N[0] ? 0 : n[0]);
    for (n[1]=0; n[1] <= N[1]; n[1]++) {
      pos[1] = (n[1] == N[1] ? 0 : n[1]);
      for (n[2]=0; n[2] <= N[2]; n[2]++) {
        pos[2] = (n[2] == N[2] ? 0 : n[2]);
        // depending on RiemannTensor use, fill FileData::work
        switch (Lat->RT.ActualUse) {
        case inactive:
        case standby:
          if ((srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]) > 0.)
            dest[destpos] = log(1.0+(srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]));
          else 
            dest[destpos] = 0.;
          destpos++;
          break;
        case active:
          posf[0] = srcpos[0+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))];
          posf[1] = srcpos[1+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))];
          posf[2] = srcpos[2+NDIM*(pos[2]+N[2]*(pos[1]+N[1]*pos[0]))];
          fac[0] = ((n[0]+N[0]*myPE)/(double)Nx);
          fac[1] = (n[1]/(double)N[1]);
          fac[2] = (n[2]/(double)N[2]); 
          RMat33Vec3(posd, Lat->RealBasis, fac);
          posf[0] += posd[0];
          posf[1] += posd[1];
          posf[2] += posd[2];
          if ((srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]) > 0.)
            posf[3] = log(1.0+(srcdens[pos[2]+N[2]*(pos[1]+N[1]*pos[0])]));
          else 
            posf[3] = 0.;
          dest[destpos+0] = posf[0];
          dest[destpos+1] = posf[1];
          dest[destpos+2] = posf[2];
          dest[destpos+3] = posf[3];
          destpos += 4;
          break;
        }
      }
    }
  }
  if (myPE) MPI_Send(dest, destpos, MPI_FLOAT, 0, OutputDensTag, P->Par.comm_ST_Psi);
}

/** Combines prepared electronic Psis density and output to file.
 * If we are process 0, open file suffixdensdat (only when RiemannTensor is used) and suffixdenspos, receive
 * FileData::work logarithmic coefficients sent by the other processes in OutputOutVisDensity(), go through all
 * nodes and save the coefficient to file - again depending on RiemannTensor use - followed by FileData::PosTemp
 * (for y and z nodes), close file(s).
 * \param *P Problem at hand 
 * \param me this ranks process in the Psi communcator ParallelSimulationData::me_comm_ST_Psi
 * \param Maxme number of processes in this Psi communcator ParallelSimulationData::Max_me_comm_ST_Psi
 */
static void CombineOutVisDensity(struct Problem *P, const int me, const int Maxme)
{
  int i,n[NDIM], N[NDIM];
  float posf[NDIM+1];
  float *source = (float *)P->Files.work;
  double posd[NDIM], fac[NDIM];
  float *Temp = (float *)P->Files.PosTemp;
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *Lev0 = &Lat->Lev[0];
  float step = P->Files.OutVisStep;
  int No=0, destpos;
  FILE *DensityData, *DensityPos;
  int Nx = Lev0->Plan0.plan->N[0]+1;
  MPI_Status status;
  if (me) return; // if we are process 0!
  N[0] = Lev0->Plan0.plan->local_nx;
  N[1] = Lev0->Plan0.plan->N[1]+1;
  N[2] = Lev0->Plan0.plan->N[2]+1;
  // Open respective file depending on RiemannTensor use
  switch (Lat->RT.ActualUse) {
  case active:
    OpenFileNo(P, &DensityPos, suffixdenspos, (int)step, "wb",P->Call.out[ReadOut]);
  case inactive:
  case standby:
    OpenFileNo(P, &DensityData, suffixdensdat, (int)step, "wb",P->Call.out[ReadOut]);
    break;
  }
  // for all processes in the communicator
  for (i=0; i< Maxme; i++) {
    if (i) {  // if process != 0, receive from this process
      /*      MPI_Probe( i, OutputDensTag, P->Par.comm_ST_Psi, &status );*/
      switch (Lat->RT.ActualUse) {
      case inactive:
      case standby:
        MPI_Recv( source, N[0]*N[1]*N[2], MPI_FLOAT, i, OutputDensTag, P->Par.comm_ST_Psi, &status );
        break;
      case active:
        MPI_Recv( source, N[0]*N[1]*N[2]*4, MPI_FLOAT, i, OutputDensTag, P->Par.comm_ST_Psi, &status );
        break;
      }
    }
    destpos = 0;
    // go through all nodes and save the coefficient to file DensityData, depending on RiemannTensor
    for (n[0]=0; n[0] < N[0]; n[0]++) {
      for (n[1]=0; n[1] < N[1]; n[1]++) {
        for (n[2]=0; n[2] < N[2]; n[2]++) {
          switch (Lat->RT.ActualUse) {
          case inactive:
          case standby:
            posf[3] = source[destpos];
            destpos++;
            (void)fwrite(&posf[3], sizeof(float), (size_t)(1), DensityData);
            No++;
            if (i==0 && n[0] == 0) 
              Temp[(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[3];
            break;
          case active:
            posf[0] = source[destpos+0];
            posf[1] = source[destpos+1];
            posf[2] = source[destpos+2];
            posf[3] = source[destpos+3];
            destpos += 4;
            (void)fwrite(posf, sizeof(float), (size_t)(NDIM), DensityPos);
            (void)fwrite(&posf[3], sizeof(float), (size_t)(1), DensityData);
            No++;
            if (i==0 && n[0] == 0) {
              fac[0] = ((n[0]+N[0]*i)/(double)(Nx-1));
              fac[1] = (n[1]/(double)(N[1]-1));
              fac[2] = (n[2]/(double)(N[2]-1)); 
              RMat33Vec3(posd, Lat->RealBasis, fac);
              posf[0] -= posd[0];
              posf[1] -= posd[1];
              posf[2] -= posd[2];
              fac[0] = ((Nx-1)/(double)(Nx-1));
              fac[1] = (n[1]/(double)(N[1]-1));
              fac[2] = (n[2]/(double)(N[2]-1)); 
              RMat33Vec3(posd, Lat->RealBasis, fac);
              posf[0] += posd[0];
              posf[1] += posd[1];
              posf[2] += posd[2];
              Temp[0+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[0]; 
              Temp[1+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[1]; 
              Temp[2+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[2]; 
              Temp[3+(NDIM+1)*(n[2]+N[2]*(n[1]+N[1]*n[0]))] = posf[3]; 
            }
            break;
          }
        }
      }
    }
  }
  n[0] = N[0];
  for (n[1]=0; n[1] < N[1]; n[1]++) {
    for (n[2]=0; n[2] < N[2]; n[2]++) {
      switch (Lat->RT.ActualUse) {
      case inactive:
      case standby:
        (void)fwrite(&Temp[n[2]+N[2]*(n[1])], sizeof(float), (size_t)(1), DensityData);
        No++;
        break;
      case active:
        (void)fwrite(&Temp[(NDIM+1)*(n[2]+N[2]*(n[1]))], sizeof(float), (size_t)(NDIM), DensityPos);
        (void)fwrite(&Temp[(NDIM+1)*(n[2]+N[2]*(n[1]))+3], sizeof(float), (size_t)(1), DensityData);
        No++;
        break;
      }
    }
  }
  if (No != Nx*N[1]*N[2]) Error(SomeError,"CombineOutVisDensity: No != points");
  switch (Lat->RT.ActualUse) {
  case active:
    fclose(DensityPos);
  case inactive:
  case standby:
    fclose(DensityData);
    break;
  }
}

/** Main output electronic Psis density for OpenDX.
 * If FileData::MeOutVis is set, calls OutputOutVisDensity() followed by CombineOutVisDensity().
 * Beforehand CalculateOutVisDensityPos() is called if RiemannTensor is used.
 * \param *P Problem at hand 
 * \param *src_dens Pointer to DensityArray which is to be displayed
 */
static void OutVisDensity(struct Problem *P, fftw_real *src_dens) 
{
  if (!P->Files.MeOutVis) return;
  if (P->Lat.RT.ActualUse == active) CalculateOutVisDensityPos(&P->Lat, &P->Files/*, P->Lat.FactorDensityR, P->Lat.FactorDensityC*/); 
  OutputOutVisDensity(P, P->Par.me_comm_ST_Psi, src_dens);
  /* Achtung hier: P->Files.work (RT->TempC, Dens->DensityCArray[TempDensity]) fuer myPE == 0 nicht veraendern !!! */
  CombineOutVisDensity(P, P->Par.me_comm_ST_Psi, P->Par.Max_me_comm_ST_Psi);
} 

/** Opening and Initializing of output measurement files.
 * If this is process 0, opens and writes top line of FileData::ForcesFile, FileData::EnergyFile.
 * and sets FileData::MeOutVis and FileData::MeOutMes (if output desired) to 1, otherwise 0.
 * \param *P Problem at hand
 */
void InitOutputFiles(struct Problem *P) 
{
  struct FileData *F = &P->Files;
  F->ForcesFile = NULL;
  F->EnergyFile = NULL;
  F->HamiltonianFile = NULL;
  F->MinimisationFile = NULL;
  F->SpreadFile = NULL;
  F->ReciSpreadFile = NULL;
  F->TemperatureFile = NULL;
  // process 0 ?
  F->MeOutVis = ((P->Par.my_color_comm_ST == 0 && P->Par.my_color_comm_ST_Psi == 0 && F->DoOutVis) ? 1 : 0);
  F->MeOutCurr = ((P->Par.my_color_comm_ST == 0 && P->Par.my_color_comm_ST_Psi == 0 && F->DoOutCurr) ? 1 : 0);
  F->MeOutMes = ((P->Par.me == 0 && F->DoOutMes) ? 1 : 0);
  OpenFile(P, &F->HamiltonianFile, suffixhamiltonianall, "w",P->Call.out[ReadOut]);

  if (!F->MeOutMes) return;
  OpenFile(P, &F->ForcesFile, suffixforcesall, "w",P->Call.out[ReadOut]);
  if (F->ForcesFile == NULL) fprintf(stderr,"Error opening ForcesFile\n");
  // write header of forces file
  fprintf(F->ForcesFile, "%s\t%s\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t%s\t%s\t%s\t\t%s\t\t%s\n",
          "Type", "No",
          "Pos0", "Pos1", "Pos2",
          "Total0", "Total1", "Total2",
          "Local0", "Local1", "Local2",
          "NLocal0", "NLocal1", "NLocal2", 
    "Magnetic0", "Magnetic1", "Magnetic2", 
          "Ewald0", "Ewald1", "Ewald2");
  OpenFile(P, &F->EnergyFile, suffixenergyall, "w",P->Call.out[ReadOut]);
  if (F->EnergyFile == NULL) fprintf(stderr,"Error opening EnergyFile\n");
  // write header of energy file
  if (P->R.DoUnOccupied) {
    fprintf(F->EnergyFile, "%s\t\t%s\t\t%s\t%s\t\t%s\t%s\t\t%s\t%s\t%s\t\t%s\t\t%s\t%s\t\t%s\t\t%s\t\t%s\n",
          "Time",
          "Total",
      "Total+Gap",
          "Kinetic", "NonLocal",
      "GapPsi",
          "Correlation", "Exchange",
          "Pseudo", "Hartree",
      "GapDensity",
          "-Gauss",
          "Ewald",
          "IonKin",
          "ETotal");
  } else {
    fprintf(F->EnergyFile, "%s\t\t%s\t\t%s\t\t%s\t%s\t%s\t%s\t\t%s\t\t%s\t\t%s\t\t%s\t\t%s\n",
      "Time",
      "Total",
      "Kinetic", "NonLocal",
      "Correlation", "Exchange",
      "Pseudo", "Hartree",
      "-Gauss",
      "Ewald",
      "IonKin",
      "ETotal");
  }
  OpenFile(P, &F->MinimisationFile, suffixminall, "w",P->Call.out[ReadOut]);
  if (F->MinimisationFile == NULL) fprintf(stderr,"Error opening MinimsationFile\n");
  // write header of minimsation file
  fprintf(F->MinimisationFile, "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n","Step", "Psi", "PsiStep", "TE", "ATE", "delta", "dEdt0", "ddEdtt0");

  OpenFile(P, &F->TemperatureFile, suffixtempall, "w",P->Call.out[ReadOut]);
  if (F->TemperatureFile == NULL) fprintf(stderr,"Error opening TemperatureFile\n");
}

/** Closes all output measurement files.
 * Closes FileData::ForcesFile and FileData::EnergyFile
 * \param *P Problem at hand
 */
void CloseOutputFiles(struct Problem *P) 
{
  struct FileData *F = &P->Files;
  if (!F->MeOutMes) return;
  if (F->ForcesFile != NULL) fclose(F->ForcesFile); // only they've been opened (thus not pointing to NULL)
  if (F->EnergyFile != NULL) fclose(F->EnergyFile);
  if (F->HamiltonianFile != NULL) fclose(F->HamiltonianFile);
  if (F->MinimisationFile != NULL) fclose(F->MinimisationFile);
  if (F->SpreadFile != NULL) fclose(F->SpreadFile);
}

/** Initialization of Problem::FileData structure for visual output.
 * If this is process 0 (and OutVis desired), allocate memory for FileData::PosTemp, FileData::work,
 * set the entries of FileData all to their corresponding values from RiemannTensor,
 * FileData::OutVisStep to zero.
 * \param *P Problem at hand
 */
void InitOutVisArray(struct Problem *P) 
{
  struct FileData *F = &P->Files;
  struct Lattice *Lat = &P->Lat;
  struct RiemannTensor *RT = &Lat->RT;
  struct LatticeLevel *Lev0 = &Lat->Lev[0];
  F->OutputPosType = NULL;
  if (!F->MeOutVis) return;
  F->OutVisStep = 0;
  F->PosTemp = (fftw_complex *)
    Malloc(sizeof(float)*(Lev0->Plan0.plan->N[1]+1)*(Lev0->Plan0.plan->N[2]+1)*
           ((Lat->RT.Use != UseRT ? 0 : NDIM)+1), "InitOutVisArray: TempC");
  F->work = (fftw_complex *)Lev0->Dens->DensityCArray[TempDensity];
  if (Lat->RT.Use != UseRT) return;
  F->TotalSize = RT->TotalSize[RTAIRT]/NDIM;
  F->LocalSizeR = RT->LocalSizeR[RTAiGcg];
  F->LocalSizeC = RT->LocalSizeC[RTAiGcg];
  F->MaxNUp = RT->MaxNUp[RTPFRto0];
  F->PosC = RT->DensityC[RTAiGcg];
  F->PosR = (fftw_real *)F->PosC;
  F->work = RT->TempC;
  F->PosFactor = RT->PosFactor[RTPFRto0];
}

static const char suffixionfor[] = ".ions.force";     //!< Suffix for ion forces file
static const char suffixionZ[] = ".ions.datZ";        //!< Suffix for ion datZ file
static const char suffixionpos[] = ".ions.pos";       //!< Suffix for ion positions file
static const char suffixiondx[] = ".ions.dx";         //!< Suffix for ions dx file
static const char suffixiondoc[] = ".ions.doc";       //!< Suffix for ions doc file
static const char suffixsrciondoc[] = ".srcion.doc";  //!< Suffix for state ions doc file
static const char suffixsrciondat[] = ".srcion.data"; //!< Suffix for state ions data file

/** Output current ions state.
 * If this is process0, open file suffixsrciondoc for writing, output Ions::Max_Types and
 * Ions::Max_IonsOfType of each type - each in a new line - closes it.
 * Then opens suffixsrciondat for binary writing, outputs Lattice:RealBasis vectors and
 * position IonType::R and speed IonType::U, closes it.
 * \param *P Problem at hand
 * \note This is the ionic counterpart to the elecontric OutputSrcPsiDensity(), storing a so far made
 *       calculation to file.
 */
static void OutSrcIons(struct Problem *P) 
{
  struct Ions *I = &P->Ion;
  double *U, *pos;
  double data[2*NDIM];
  int is,ia,i;
  FILE *SrcIonDoc, *SrcIonData;

  if (!(P->Par.me == 0)) return;
  
  // output of ion types and numbers per type
  OpenFile(P, &SrcIonDoc, suffixsrciondoc, "w",P->Call.out[ReadOut]);
  fprintf(SrcIonDoc,"%i\n", I->Max_Types);
  for (is=0; is < I->Max_Types; is++)
    fprintf(SrcIonDoc,"%i\n", I->I[is].Max_IonsOfType);
  fclose(SrcIonDoc);

  OpenFile(P, &SrcIonData, suffixsrciondat, "wb",P->Call.out[ReadOut]);
  (void)fwrite(P->Lat.RealBasis, sizeof(double), (size_t)(NDIM_NDIM), SrcIonData);
  for (is=0; is < I->Max_Types; is++) {
    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
      U = &I->I[is].U[NDIM*ia];
      pos = &I->I[is].R[NDIM*ia];
      for (i=0;i<NDIM;i++) {
        data[i] = pos[i];
        data[i+NDIM] = U[i];
      }
      (void)fwrite(&data, sizeof(double), (size_t)(2*NDIM), SrcIonData);
    }
  }
  fclose(SrcIonData);
}

/** Read ions state from a file.
 * Reads the suffixsrciondoc file and checks it against the current state in Ions regarding
 * IonType::MaxTypes and IonType::Max_IonsOfType, closes it.
 * Afterwards, opens suffixsrciondat for binary reading, retrieves the basis checking it against
 * Problem::Lattice::RealBasis. If ok, reads position IonType::R and speed IonType::U, closes it.
 * \param *P Problem at hand
 * \return 1 - successful, 0 - failure
 * \note This is the ionic counterpart to the electronic ReadSrcPsiDensity(), see also OutSrcIons().
 */
int ReadSrcIons(struct Problem *P) 
{
  struct Ions *I = &P->Ion;
  double *U, *pos;
  double data[2*NDIM];
  int is,ia,i;
  int Max_Types;
  int *Max_IonsOfType = NULL;
  double RealBasis[NDIM_NDIM];
  FILE *SrcIonDoc, *SrcIonData;
  // read the doc file and check
  if (OpenFile(P, &SrcIonDoc, suffixsrciondoc, "r",P->Call.out[ReadOut])) {
    if (fscanf(SrcIonDoc,"%i", &Max_Types) != 1)
      //Error(SomeError, "ReadSrcIons: read error");
      return 0;
    if (Max_Types != I->Max_Types)
      //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, MaxTypes");
      return 0;
    Max_IonsOfType = (int *) Malloc(Max_Types*sizeof(int), "ReadSrcIons: Max_IonsOfType");
    for (is=0; is < Max_Types; is++) {
      if (fscanf(SrcIonDoc,"%i", &Max_IonsOfType[is]) != 1)
        //Error(SomeError, "ReadSrcIons: read error");
      return 0;
      if (Max_IonsOfType[is] != I->I[is].Max_IonsOfType)
        //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, Max_IonsOfType");
      return 0;
    }
    fclose(SrcIonDoc);
    // read basis, then positions and speeds of ions
    if (OpenFile(P, &SrcIonData, suffixsrciondat, "rb",P->Call.out[ReadOut])) {
      if (fread(RealBasis, sizeof(double), (size_t)(NDIM_NDIM), SrcIonData) != NDIM_NDIM)
        //Error(SomeError, "ReadSrcIons: read error");
        return 0;
      for (i=0; i < NDIM_NDIM; i++)
        if (RealBasis[i] != P->Lat.RealBasis[i])
          //Error(SomeError, "ReadSrcIons: srcion file does not fit to system, RealBasis");
          return 0;
      for (is=0; is < I->Max_Types; is++) {
        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
          if (fread(&data, sizeof(double), (size_t)(2*NDIM), SrcIonData) != 2*NDIM)
                //Error(SomeError, "ReadSrcIons: read error");
            return 0;
          U = &I->I[is].U[NDIM*ia];
          pos = &I->I[is].R[NDIM*ia];
          for (i=0;i<NDIM;i++) {
                pos[i] = data[i];
                U[i] = data[i+NDIM];
          }
        }
      }
      fclose(SrcIonData);
    }
    if (Max_IonsOfType != NULL) Free(Max_IonsOfType, "Max_IonsOfType");
  }
  return 1;
}

/** Output of ion doc, dx, forces and positions file for OpenDX.
 * If this is process 0,
 * open, fill and close IonDoc file suffixiondoc,
 * open, fill for each FileData::OutVisStep and close IonDX file suffixiondx
 * for every
 * open suffixionfor, suffixionpos (and suffixionZ in case of only FileData::OutVisStep), fill
 * them with the ion forces IonType::FIon and positions IonType::R of each type and each ion per type,
 * close them all.
 * \param *P Problem at hand
 */
static void OutVisIons(struct Problem *P) 
{
  struct Ions *I = &P->Ion;
  struct FileData *F = &P->Files;
  int i,is,ia;
  double *fion, *pos;
  float data[6];      // holds temporarily twice NDIM values as write buffer
  int Z;
  char *datnamef, *datnameZ, *posname;
  FILE *IonsDataF, *IonsDataZ, *IonsPos, *IonsDoc, *IonsDx;
  if (!P->Files.MeOutVis && P->Par.me == 0) return;
  // generate file names
  datnamef = (char*)
    malloc(strlen(P->Files.mainname)+strlen(suffixionfor) + 1);
  sprintf(datnamef, "%s%s", P->Files.mainname, suffixionfor);
  datnameZ = (char*)
    malloc(strlen(P->Files.mainname)+strlen(suffixionZ) + 1);
  sprintf(datnameZ, "%s%s", P->Files.mainname, suffixionZ);
  posname = (char*)
    malloc(strlen(P->Files.mainname)+strlen(suffixionpos) + 1);
  sprintf(posname, "%s%s", P->Files.mainname, suffixionpos);
  // open, fill and close doc file
  if (OpenFile(P, &IonsDoc, suffixiondoc, "w",P->Call.out[ReadOut])) {
    fprintf(IonsDoc,"IonsPos file = %s.####\n", posname);
    fprintf(IonsDoc,"IonsForce file = %s.####\n", datnamef);
    fprintf(IonsDoc,"format = ieee float (Bytes %lu) %s = Force(3)\n",(unsigned long) sizeof(float),msb);
    fprintf(IonsDoc,"IonsZ file = %s.####\n", datnameZ);
    fprintf(IonsDoc,"format = int (Bytes %lu) %s = Z(1)\n",(unsigned long) sizeof(int),msb);
    fprintf(IonsDoc,"points = %i\n", I->Max_TotalIons);
    fprintf(IonsDoc,"majority = row\n");
    fprintf(IonsDoc,"TimeSeries = %i\n",F->OutVisStep+1);
    fclose(IonsDoc);
  }
  // open dx file and fill it with each output step, close it
  if (OpenFile(P, &IonsDx, suffixiondx, "w",P->Call.out[ReadOut])) {
    for (i=0; i < F->OutVisStep+1; i++) {
      if (i==0) {
        /*      fprintf(IonsDx,"object \"ioncon\" class array type int rank 1 shape 2 items 0 data follows\nattribute \"element type\" string \"lines\"\nattribute \"ref\" string \"positions\"\n\n"); */
        fprintf(IonsDx,"object \"iondatZ\" class array type int rank 0 items %i %s binary\n",I->Max_TotalIons,msb);
        fprintf(IonsDx,"data file %s,0\n",datnameZ);
        fprintf(IonsDx,"attribute \"dep\" string \"positions\"\n\n");
      }
      
      fprintf(IonsDx,"object \"ionpos.%04i\" class array type float rank 1 shape 3 items %i %s binary\n",i,I->Max_TotalIons,msb);
      fprintf(IonsDx,"data file %s.%04i,0\n\n",posname,i);
      
      fprintf(IonsDx,"object \"iondatF.%04i\" class array type float rank 1 shape 3 items %i %s binary\n",i,I->Max_TotalIons,msb);
      fprintf(IonsDx,"data file %s.%04i,0\n",datnamef,i);
      fprintf(IonsDx,"attribute \"dep\" string \"positions\"\n\n");
      
      fprintf(IonsDx,"object \"ionobjF.%04i\" class field\n",i);
      fprintf(IonsDx,"component \"positions\" \"ionpos.%04i\"\n",i);
      fprintf(IonsDx,"component \"data\" \"iondatF.%04i\"\n",i);
      /*fprintf(IonsDx,"component \"connections\" \"ioncon\"\n\n");*/
      
      fprintf(IonsDx,"object \"ionobjZ.%04i\" class field\n",i);
      fprintf(IonsDx,"component \"positions\" \"ionpos.%04i\"\n",i);
      fprintf(IonsDx,"component \"data\" \"iondatZ\"\n");
      /*    fprintf(IonsDx,"component \"connections\" \"ioncon\"\n\n");*/
    }
    fprintf(IonsDx,"\nobject \"ionseriesF\" class series\n");
    for (i=0; i < F->OutVisStep+1; i++) 
      fprintf(IonsDx,"member %i \"ionobjF.%04i\" position %f\n",i,i,(float)i);
    fprintf(IonsDx,"\nobject \"ionseriesZ\" class series\n");
    for (i=0; i < F->OutVisStep+1; i++) 
      fprintf(IonsDx,"member %i \"ionobjZ.%04i\" position %f\n",i,i,(float)i);
    fprintf(IonsDx,"end\n");
    fclose(IonsDx);
  }
  free(datnamef);
  free(datnameZ);
  free(posname);
  // open IonForces, IonZ and IonPosition file, write forces respectively positions for each ion of each type, close them
  if (OpenFileNo(P, &IonsDataF, suffixionfor, F->OutVisStep, "wb",P->Call.out[ReadOut])) {
    if (F->OutVisStep == 0) 
      OpenFile(P, &IonsDataZ, suffixionZ, "wb",P->Call.out[ReadOut]);
    if (OpenFileNo(P, &IonsPos, suffixionpos, F->OutVisStep, "wb",P->Call.out[ReadOut])) {
      for (is=0; is < I->Max_Types; is++) {
        for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
          fion = &I->I[is].FIon[NDIM*ia];
          pos = &I->I[is].R[NDIM*ia];
          for (i=0;i<3;i++) {
                data[i+3] = fion[i];
                data[i] = pos[i];
          }
          Z = I->I[is].Z;
          if (fwrite(&data[0],sizeof(float), (size_t)(3),IonsPos) != 3)
            Error(FileOpenParams, "Error writing ion positions!");
          if (F->OutVisStep == 0) (void)fwrite(&Z,sizeof(int), (size_t)(1),IonsDataZ);
          if (fwrite(&data[3],sizeof(float), (size_t)(3),IonsDataF) != 3)
            Error(FileOpenParams, "Error writing ionic forces!");
        }
      }
      fclose(IonsPos);
    }
    if (F->OutVisStep == 0) 
      fclose(IonsDataZ);
    fclose(IonsDataF);
  }
}

/** Output electronic density and ion state files with so far made calculations.
 * If CallOptions::WriteSrcFiles is set, OutputSrcPsiDensity() and OutSrcIons() are called.
 * \param *P Problem at hand
 * \param type which PsiTypeTag should be put to file
 */
void OutputVisSrcFiles(struct Problem *P, enum PsiTypeTag type)
{
  if (P->Call.WriteSrcFiles) {
    if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Writing %s srcpsi to disk\n", P->Par.me, P->R.MinimisationName[type]);
    OutputSrcPsiDensity(P, type);
    if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Writing srcion to disk\n", P->Par.me);
    OutSrcIons(P);
  } 
  // if (!P->Files.MeOutVis) return;
  // P->Files.OutputPosType = 
  //   Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+1),"OutputVis");
  // P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
  // CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
  // OutVisDensity(P);
  // OutVisIons(P);
  // if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
  // /* P->Files.OutVisStep++; Genau ebend nicht hochzaehlen - wird immer ueberschrieben */
}

/** Main output total electronic density and ion data for OpenDX.
 * Calls subsequently preparing CreateDensityOutputGeneral(), then output of electronic
 * densities OutVisDensity() and ion data OutVisIons(), increasing finally FileData::OutVisStep.
 * \param *P Problem at hand
  * \param *srcdens Pointer to DensityArray which is to be displayed
 * \note Output is made only RunStruct::OutVisStep steps and if FileData::MeOutVis is set.
 */
void OutputVis(struct Problem *P, fftw_real *srcdens)
{
  if (!P->Files.MeOutVis) return;
  P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+1),"OutputVis");
  P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
  
  CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
  OutVisDensity(P, srcdens);
  OutVisIons(P);
  if(P->Call.out[MinOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
  P->Files.OutVisStep++;
}

/** Output of each current density component for OpenDX.
 * \param *P Problem at hand
 * \note Output is made only RunStruct::OutVisStep steps and if FileData::MeOutVis is set.
 */
void OutputCurrentDensity(struct Problem *P)
{
  int index, i, r;
  fftw_real *density = P->R.Lev0->Dens->DensityArray[ActualDensity];
  fftw_real *CurrentDensity[NDIM*NDIM];
  if (!P->Files.MeOutCurr) return;

  P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+(1)*NDIM),"OutputVis");
  for (i=0;i<(1)*NDIM;i++)
    P->Files.OutputPosType[P->Files.OutVisStep+i] = P->Lat.RT.ActualUse;
  if(P->Call.out[PsiOut]) fprintf(stderr,"(%i) OutVisStep %i, OutputPosType %p\n",P->Par.me, P->Files.OutVisStep, P->Files.OutputPosType);
  
  // due to preprocessor values we can't put the following stuff into a loop
  CurrentDensity[0] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity0];
  CurrentDensity[1] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity1];
  CurrentDensity[2] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity2];
  CurrentDensity[3] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity3];
  CurrentDensity[4] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity4];
  CurrentDensity[5] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity5];
  CurrentDensity[6] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity6];
  CurrentDensity[7] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity7];
  CurrentDensity[8] = (fftw_real *) P->R.Lev0->Dens->DensityArray[CurrentDensity8];
  
  // output current density, not vector component
  for (index=0;index<NDIM;index++) {
    // evaluate euclidian norm for given B component
    SetArrayToDouble0((double *)density,P->R.Lev0->Dens->TotalSize*2); // reset
    for (r=0;r<P->R.Lev0->Dens->LocalSizeR;r++) {
      for (i=0;i<NDIM;i++)
        density[r] += CurrentDensity[i + index*NDIM][r]*CurrentDensity[i + index*NDIM][r];
      density[r] = sqrt(density[r]);      
    }
    // output
    CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
    OutVisDensity(P, density);
    OutVisIons(P);
    if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
    P->Files.OutVisStep++;
  }
}

/** Output each orbital in a certain step order for OpenDX.
 * \param *P Problem at hand
 * \param offset from which step do we start
 * \param increment by which increment do we advance step-wise
 * \param type Only PsiTypeTag orbitals are displayed 
 */
void OutputVisAllOrbital(struct Problem *P, int offset, int increment, enum PsiTypeTag type) {
  struct Lattice *Lat = &P->Lat;
  struct Psis *Psi = &Lat->Psi;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevS = R->LevS;
  struct LatticeLevel *Lev0 = R->Lev0;
  struct Density *Dens0 = Lev0->Dens;
  struct OnePsiElement *OnePsiA, *LOnePsiA;
  MPI_Status status;
  int ElementSize = (sizeof(fftw_complex) / sizeof(double));
  fftw_complex *LPsiDatA;  
  int i, p, RecvSource;
  int Num = Psi->NoOfPsis;

  if (!P->Files.MeOutVis) return;

  P->Files.OutputPosType = (enum ModeType *) Realloc(P->Files.OutputPosType,sizeof(enum ModeType)*(P->Files.OutVisStep+(Num)),"OutputVis");
  
  P->Files.OutVisStep += offset;
  P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;  
  for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) {  // go through all wave functions (plus the extra one for each process)
    OnePsiA = &Psi->AllPsiStatus[i];    // grab the desired OnePsiA
    if (OnePsiA->PsiType == type) {   // drop if extra one
      if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
         LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
      else 
         LOnePsiA = NULL;
      if (LOnePsiA == NULL) {   // if it's not local ... receive it from respective process into TempPsi
        RecvSource = OnePsiA->my_color_comm_ST_Psi;
        MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, WannierTag1, P->Par.comm_ST_PsiT, &status );
        LPsiDatA=LevS->LPsi->TempPsi;
      } else {                  // .. otherwise send it to all other processes (Max_me... - 1)
        for (p=0;p<P->Par.Max_me_comm_ST_PsiT;p++)
          if (p != OnePsiA->my_color_comm_ST_Psi) 
            MPI_Send( LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, p, WannierTag1, P->Par.comm_ST_PsiT);
        LPsiDatA=LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo];
      } // LPsiDatA is now set to the coefficients of OnePsi either stored or MPI_Received

      if (P->Files.MeOutVis) {
        P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;
        // recalculate density for the specific wave function ...
        CalculateOneDensityR(Lat, LevS, Dens0, LPsiDatA, Dens0->DensityArray[ActualDensity], R->FactorDensityR, 0);
        // ... and output (wherein ActualDensity is used instead of TotalDensity)
        //OutputVis(P);
        CreateDensityOutputGeneral(P, P->Par.me_comm_ST_Psi);
        OutVisDensity(P, Dens0->DensityArray[ActualDensity]);
        OutVisIons(P);
        if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Written OutVisStep %i to disk\n", P->Par.me, P->Files.OutVisStep);
        P->Files.OutVisStep+=increment;
        P->Files.OutputPosType[P->Files.OutVisStep] = P->Lat.RT.ActualUse;  
      }
    }
  }
}

/** Read source files containing stored calculations.
 * Calls ReadSrcPsiDensity() and ReadSrcIons().
 * \param *P Problem at hand
 */
void ReadSrcFiles(struct Problem *P)
{
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadSrcPsiDensity\n", P->Par.me);
  ReadSrcPsiDensity(P, Occupied, 0, P->R.LevSNo);
  ReadSrcPsiDensity(P, UnOccupied, 0, P->R.LevSNo);
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadSrcIons\n", P->Par.me);
  ReadSrcIons(P);
}

/** Plots a cut plane of the real density of one wave function.
 * \param *P Problem at hand
 * \param index index of axis (vector orthogonal to plane)
 * \param node node specifying where to cut at the given axis
 * \param wavenr global number of wave function
 * \param *density density array to plot
 * \sa PlotVectorPlane() - very similar 
 */
void PlotSrcPlane(struct Problem *P, int index, double node, int wavenr, fftw_real *density) 
{
  struct RunStruct *R = &P->R;
  struct Lattice *Lat = &P->Lat;
  struct LatticeLevel *Lev0 = R->Lev0;
  const int myPE = P->Par.me_comm_ST_Psi;
  char filename[255], spin[12];
  char *suchpointer;
  FILE *PlotFile = NULL;
  time_t seconds; 
  

  // open file
  time(&seconds); // get current time

  switch (Lat->Psi.PsiST) {
  case SpinDouble:    
    sprintf(&filename[0], ".psi%i_cut%i.csv", wavenr, index);
    strncat(spin,"SpinDouble",12);
    break;
  case SpinUp:
    sprintf(&filename[0], ".psi%i_cut%i_up.csv", wavenr, index);
    strncat(spin,"SpinUp",12);
    break;
  case SpinDown:
    sprintf(&filename[0], ".psi%i_cut%i_down.csv", wavenr, index);
    strncat(spin,"SpinDown",12);
    break;
  }
  
  if (!myPE) { // only process 0 writes to file
    OpenFile(P, &PlotFile, filename, "w", P->Call.out[ReadOut]);
    strcpy(filename, ctime(&seconds));
    suchpointer = strchr(filename, '\n');
    if (suchpointer != NULL)
      *suchpointer = '\0';
    if (PlotFile != NULL) { 
      fprintf(PlotFile,"# Psi %i, type %s (real density) plot of plane perpendicular to direction e_%i at node %lg, seed %i, config %s, run on %s, #cpus %i", wavenr, spin, index, node, R->Seed, P->Files.default_path, filename, P->Par.Max_me_comm_ST_Psi);
      fprintf(PlotFile,"\n");
    } else { Error(SomeError, "PlotSrcPlane: Opening Plot File"); }
  }
    
  // plot density
  PlotRealDensity(P, Lev0, PlotFile, index, node, density, density);

  if (PlotFile != NULL) { 
    // close file
    fclose(PlotFile);
  }
}

/** plots a cut plane of a given 3d real density.
 * \param *P Problem at hand, contains pointer to Lattice structure
 * \param *Lev LatticeLevel of the real density
 * \param PlotFile file pointer (already open and valid)
 * \param index index of lattice axis
 * \param n_orth position on lattice axis where to cut
 * \param *density1 first real density array
 * \param *density2 second real density array (point to \a *density1 if not needed)
 */
void PlotRealDensity(struct Problem *P, struct LatticeLevel *Lev, FILE *PlotFile, int index, double n_orth, fftw_real *density1, fftw_real *density2)
{
  struct Lattice *Lat = &P->Lat;
  int n[NDIM], n0;
  int N[NDIM];
  N[0] = Lev->Plan0.plan->N[0];
  N[1] = Lev->Plan0.plan->N[1];
  N[2] = Lev->Plan0.plan->N[2];
  const int N0 = Lev->Plan0.plan->local_nx;
  const int myPE = P->Par.me_comm_ST_Psi;
  double fac[NDIM], x[NDIM];
  int i0, i = 0;
  int PE, zahl;
  double *buffer;
  MPI_Status status;
  int sizes[P->Par.Max_me_comm_ST_Psi], c0, c1;
  double nodes[NDIM], node[NDIM];
  
  for(i=0;i<NDIM;i++) {
    nodes[i] = (i == index) ? n_orth : 0.;
    node[i] = 0.;
  }
  RMat33Vec3(node, Lat->ReciBasis, nodes);  // transform cartesian coordinates into cell coordinates [0,1]^3
  for(i=0;i<NDIM;i++) // now N^3 within node range of discrete grid
    node[i] = (int)(node[i]*N[i]/(2.*PI));
  fprintf(stderr,"(%i) n_orth %lg, index %i converted to plane offset vector (%lg, %lg, %lg).\n", P->Par.me, n_orth, index, node[0], node[1], node[2]);
  
  switch (index) {
    case 0:
        zahl = 4*N[1]*N[2];
      break;
    case 1:
        zahl = 4*N0*N[2];
      break;
    case 2:
        zahl = 4*N0*N[1];
      break;
  } 
  fprintf(stderr,"(%i) buffer size %i\n", P->Par.me, zahl);
  buffer = Malloc(sizeof(double)*zahl,"PlotRealDensity: buffer");
  
  c0 = cross(index,0);
  c1 = cross(index,1);
  // then for every point on the grid in real space ...
  i=0;
  for (n0=0;n0<N0;n0++)  // only local points on x axis
    for (n[1]=0;n[1]<N[1];n[1]++)
      for (n[2]=0;n[2]<N[2];n[2]++) {
        n[0]=n0 + N0*myPE; // global relative coordinate: due to partitoning of x-axis in PEPGamma>1 case
        if (n[index] == (int)node[index]) { // only on the correct plane orthogonal to desired axis and at desired node ...
          fac[0] = (double)n[0]/(double)N[0];
          fac[1] = (double)n[1]/(double)N[1];
          fac[2] = (double)n[2]/(double)N[2];
          RMat33Vec3(x, Lat->RealBasis, fac); // relative coordinate times basis matrix gives absolute ones
          i0 = n[2]+N[2]*(n[1]+N[1]*n0);  // index to local density array

          buffer[i++] = x[c0];  // fill buffer
          buffer[i++] = x[c1];
          buffer[i++] = density1[i0];
          buffer[i++] = density2[i0];
          if (i > zahl) Error(SomeError, "PlotRealDensity: buffer too small!");
        }
      }
  // exchange sizes of each buffer
  MPI_Allgather(&i, 1, MPI_INT, sizes, 1, MPI_INT, P->Par.comm_ST_Psi);   
  if (myPE == 0) {
    for (PE=0; PE < P->Par.Max_me_comm_ST_Psi; PE++) {
      if (PE != 0) {  
        // receive them
        if (MPI_Recv(buffer, sizes[PE], MPI_DOUBLE, PE, PlotRealDensityTag, P->Par.comm_ST_Psi, &status) != MPI_SUCCESS)
          Error(SomeError, "PlotRealDensity: MPI_Recv failure!");
        MPI_Get_count(&status, MPI_DOUBLE, &zahl);
        if (zahl != sizes[PE]) 
          Error(SomeError, "PlotRealDensity: received unexpected amount of elements!");
      }
      //write them: local one (still in buffer) and received ones
      for (i0 = 0; i0 < sizes[PE];) {
        fprintf(PlotFile,"%e", buffer[(i0)++]);
        if ((i0 % 4) == 0) {
          fprintf(PlotFile,"\n");
        } else {
          fprintf(PlotFile,"\t");
        }
      }
    }
  } else { // send them
    if (MPI_Send(buffer, i, MPI_DOUBLE, 0, PlotRealDensityTag, P->Par.comm_ST_Psi) != MPI_SUCCESS)
      Error(SomeError, "PlotRealDensity: MPI_Send failure!"); 
  }
  Free(buffer, "PlotRealDensity: buffer");
}