source: pcp/src/data.h@ cfaf42

#ifndef data_h
#define data_h

/** \file data.h
 * Defines most of the important data structures.
 * Such as structures containing command line options CallOptions, densities Density, energies Energy,
 * lattice information Lattice and levelwise LatticeLevel, for all wave functions Psis downto one Psi
 * OnePsiElement and OnePsiElementAddData, wave functions on the lattice level LevelPsi, the reciprocal
 * grid vector OneGData, filenames and paths FileData oder super structures containing the whole Problem
 * or data for the parallelisation ParallelSimulationData and the fft plans RPlans.
 *
  Project: ParallelCarParrinello
  Jan Hamaekers
  2000

  File: data.h
  $Id: data.h,v 1.102 2007-02-09 09:13:48 foo Exp $
*/
#include "mpi.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 "defs.h"
#include <signal.h>
/*extern volatile sig_atomic_t cpulim;*/
#define PROCDIM (2)                     //!< specifying the number in \ref ParallelDivideDims
#define MaxOutGroup 7           //!< specifying the number in \ref OutGroup
//! various verbose output groups, individually associable with levels (LeaderOut = process 0,...)
enum OutGroup { MinOut,         //!< minimum output
                NormalOut,      //!< normal output
                ValueOut,       //!< output of calculated value
                LeaderOut,      //!< output of process 0
                ReadOut,        //!< output of parsed values
                PsiOut,         //!< output of Psi-related numbers
                StepLeaderOut   //!< output of minimsation step and NIDensity
              };
//! enumerating output verbosity levels
enum OutType { OutNone,     //!< no -v specified (no output)
               OutNormal,   //!< -v (status information, minimum output, some calculate values)
               OutNormalP,  //!< -vv (status information, minimum output, more calculated values)
               OutMore,     //!< -vvv (steady output while minimising)
               OutMoreP,    //!< -vvvv (all processes)
               OutAll,      //!< not used
               OutAllP      //!< not used
             };

//! enumerating whether to parse source coefficients from file
enum ParseSrcDensities { DoNotParse,
                         DoReadAllSrcDensities,
                         DoReadOccupiedSrcDensities,
                         DoReadAndMinimise
};
             
//! enumerating usage status
enum UseType { NotInUse,  //!< is not in use
               InUse      //!< is in use
             };
//! enumerating the two groups in parallel divide&conquer strategy
enum ParallelDivideDims { PEGamma, //!< process in Psi group among coefficients are shared
                          PEPsi    //!< process groups among wave functions are shared
                        };
//! enumerating RiemannTensor usage
enum UseRiemannTensor { UseNotRT,   //!< don't use RiemannTensor calculus
                        UseRT       //!< use RiemannTensor calculus
                      };        
//! enumerating how to deal with spin
enum UseSpinType { UseSpinDouble,   //!< Spin is always double in each orbit (occupation number is 2)
                   UseSpinUpDown    //!< Treat each orbit differently for its SpinType#SpinUp and SpinType#SpinDown part
                 };     
//! SpinType is separated Up and Down or combined Double
enum SpinType { SpinDouble,   //!< Type is double, orbits always fully occupied
                SpinUp,       //!< Type of orbit spin is up
                SpinDown      //!< Type of orbit spin is down
              }; /*!< Double == 0 !!! */
#define MaxDensityTypes 22 //!< number of different densities
#define MaxInitDensityTypes 13 //!< cardinal number of density up to which all before need to be reseted in InitDensityCalculation()
//! enumerating the various density types
enum DensityTypes { ActualDensity,      //!< current density (in calculation)
                    TotalDensity,       //!< total density over all processes
                    TotalLocalDensity,  //!< total density stored in this process
                    TotalUpDensity,     //!< total density with spin up over all processes
                    TotalDownDensity,   //!< total density with spin down over all processes
                    CoreWaveDensity,    //!< density of the wave functions in the core range, see PseudoPot
                    HGcDensity,         //!< sum of gaussian and local PseudoPot'ential density
                    GapDensity,         //!< density of the all "unoccupied" additional wave functions
                    GapUpDensity,       //!< density of the all "unoccupied" additional wave functions
                    GapDownDensity,     //!< density of the all "unoccupied" additional wave functions
                    GapLocalDensity,    //!< density of the "unoccupied" additional wave functions stored in this process
                    TempDensity,        //!< temporal working array, used mostly in fft transforms
                    Temp2Density,        //!< another temporal working array
                    CurrentDensity0,     //!< current density[0][0] = d/dB_0 current[0]
                    CurrentDensity1,     //!< current density[1][0] = d/dB_0 current[1]
                    CurrentDensity2,     //!< current density[2][0] = d/dB_0 current[2]
                    CurrentDensity3,     //!< current density[0][1] = d/dB_1 current[0]
                    CurrentDensity4,     //!< current density[1][1] = d/dB_1 current[1]
                    CurrentDensity5,     //!< current density[2][1] = d/dB_1 current[2]
                    CurrentDensity6,     //!< current density[0][2] = d/dB_2 current[0]
                    CurrentDensity7,     //!< current density[1][2] = d/dB_2 current[1]
                    CurrentDensity8      //!< current density[2][2] = d/dB_2 current[2]
                  };    
/*! enumerating density types which differ for DensityArray and DensityCArray, as follows:
  CoreWave   only DensityArray (real)
  ActualPsiDensity only DensityCArray (complex)
  HGcDensity       only DensityArray (real)
  HGDensity  only DensityCArray (complex)
 */

enum DoubleDensityTypes { HGDensity=HGcDensity,             //!< the local potential \f$V^H (G) + V^{ps,loc} (G)\f$
                          ActualPsiDensity=CoreWaveDensity  //!< the local psi coefficients
                        };
#define MaxPsiNoType 7  //!< maximum number of different wave function types
//! enumerating Psi numbering type
enum PsiNoType { PsiMaxNo,        //!< maximum number of wave functions
                 PsiMaxNoDouble,  //!< maximum number of wave functions in SpinType#SpinDouble case
                 PsiMaxNoUp,      //!< maximum number of wave functions with spin up
                 PsiMaxNoDown,    //!< maximum number of wave functions with spin down
                 PsiMaxAdd       //!< number of additional (unoccupied) wave functions helpful in improving the minimisation
               };
//! Enumerating different G vector types, due to the gamma point symmetry! (mirroring of Gs)
enum DoubleGType { DoubleGNot,    //!< not a reciprocal grid vector on (z=0)-plane (no symmetry can be utilized)
                   DoubleG0,      //!< is the (0,0,0) reciprocal grid vector
                   DoubleGUse,    //!< reciprocal grid vector which resides on the (z=0)-plane in the first and third sector and thus has symmetry \f$c_{i,G}^\ast = -c_{i,-G}\f$
                   DoubleGNotUse  //!< reciprocal grid vector which resides on the (z=0)-plane in the second and fourth sector 
                 };
//! Enumerating states of the current GramSchmidt-Orthonormalization of a wave function
enum PsiGramSchStatusType { NotOrthogonal,  //!< is not yet orthogonalized
                            IsOrthogonal,   //!< is orthohonal yet not normalized 
                            IsOrthonormal,  //!< is orthonormal
                            NotUsedToOrtho  //!< not touched during GramSch()
                          };
//! Enumerating what remains to do for this wave function
enum PsiGramSchToDoType { Orthonormalize,  //!< must still be orthogonalized and normalized
                          Orthogonalize    //!< must still be normalized
                        };
                        
//! Enumerating the ToDo status for this wave function
enum MinimisationStatusType { Minimalise,           //!< wave function ought to be minimalized
                              DoNotMinimalise,       //!< wave function is used to evaluate however not itself minimised
                              NotUsedToMinimalise   //!< ignore this wave function completely during minimisation
                            };


#define perturbations 6

//! Enumerating the Type of this wave function
enum PsiTypeTag { Occupied,     //!< normal occupied wave function
                  UnOccupied,    //!< additional unoccupied wave function foor gap energy
                  Perturbed_P0,    //!< perturbed wave function \f$|\varphi^{(p_0)}\rangle\f$ used to evaluate current density
                  Perturbed_P1,    //!< perturbed wave function \f$|\varphi^{(p_1)}\rangle\f$ used to evaluate current density
                  Perturbed_P2,    //!< perturbed wave function \f$|\varphi^{(p_2)}\rangle\f$ used to evaluate current density
                  Perturbed_RxP0,  //!< perturbed wave function \f$|\varphi^{(r\times p)_0}\rangle\f$ used to evaluate current density
                  Perturbed_RxP1,  //!< perturbed wave function \f$|\varphi^{(r\times p)_1}\rangle\f$ used to evaluate current density
                  Perturbed_RxP2,  //!< perturbed wave function \f$|\varphi^{(r\times p)_2}\rangle\f$ used to evaluate current density
                  Extra        //!< extra wave function (used for gradient calculation)
                };

//! Enumerating whether this is active or not (such as motion of ions, use of RiemannTensor, ...)
enum ModeType { inactive,   //!< generally deactivated
                active,     //!< activated
                standby     //!< temporarily deactivated
              };
/* MPI Tags */
#define GramSchTag1 100   //!< Message consists of wave function coefficients needed during the Orthogonalization in GramSch()
#define GramSchTag2 101   //!< Message consists of orthogonal projections calculated during the Orthogonalization in GramSch()
#define GramSchTag3 102   //!< Message consists of wave function coefficients during Testing in TestGramSch()
#define InterTag1 110
#define DensityTag1 120   //!< Message on exchange of complex TotalUpDensity, see InitDensityCalculation()
#define DensityTag2 121   //!< Message on exchange of complex TotalDownDensity, see InitDensityCalculation()
#define DensityTag3 122   //!< Message on exchange of TotalUpDensity, see InitDensityCalculation()
#define DensityTag4 123   //!< Message on exchange of TotalDownDensity, see InitDensityCalculation()
#define DensityTag5 124   //!< Message on exchange of NIUpDensity, see InitDensityCalculation()
#define DensityTag6 125   //!< Message on exchange of NIDownDensity, see InitDensityCalculation()
#define DensityTag7 126   //!< Message on exchange of complex GapUpDensity, see InitDensityCalculation()
#define DensityTag8 127   //!< Message on exchange of complex GapDownDensity, see InitDensityCalculation()
#define DensityTag9 128   //!< Message on exchange of GapUpDensity, see InitDensityCalculation()
#define DensityTag0 129   //!< Message on exchange of GapDownDensity, see InitDensityCalculation()
#define OutputDensTag 140   //!< Message consists of FileData::work being sent to process 0 on output of visual data, see OutputOutVisDensity()
#define ReadSrcPsiTag 141 //!< Used during reading and send/recv of psi coefficients in ReadSrcPsiDensity()
#define OutputSrcPsiTag 142 //!< Message consists of Density::DensityArray[TempDensity] being sent to process 0 on output of source Psis (saving of current state), see OutputSrcPsiDensity()
#define ParseWannierTag 143 //!< Used during send/recv of wannier centres and spread in ParseWannierFile()
#define PlotRealDensityTag 145 //!< 4 doubles exchanged during PlotRealDensity() storing to file
#define AllMaxLocalNoTag1 150
#define AllMaxLocalNoTag2 151
#define HamiltonianTag 160  //!< Message consists of orbit coefficients sent/received during the setup of the hamiltonian matrix
#define HamiltonianTag2 161  //!< Message consists of orbit coefficients sent/received during the setup of the hamiltonian matrix
#define WannierTag1 170   //!< Message consists of orbital coefficients sent/received during iteration of wannier functions
#define WannierTag2 171   //!< Message consists of orbital coefficients sent/received during application of transformation matrix
#define WannierTag3 172   //!< Message consists of matrix elements A_ij sent/received during evaluating of operator matrices
#define WannierTag4 173   //!< Message consists of matrix elements B_ij sent/received during evaluating of operator matrices
#define OtherPsiTag1 180  //!< Message consists of OnePsiElement data from the Spinype#SpinUp group
#define OtherPsiTag2 181  //!< Message consists of OnePsiElement data from the Spinype#SpinDown group
#define StopTag1 190  //!< Message consists of stop flag from the Spinype#SpinUp group
#define StopTag2 191  //!< Message consists of stop flag from the Spinype#SpinDown group
#define EnergyTag1 200  //!< Message consists of Energy#AllUpPsiEnergy during EnergyAllReduce()
#define EnergyTag2 201  //!< Message consists of Energy#AllDownPsiEnergy during EnergyAllReduce()
#define EnergyTag3 203
#define EnergyTag4 204
#define EnergyTag5 205  //!< Message consists of partial Energy#TotalEnergy results, sent during EnergyAllReduce()
#define EnergyTag6 206  //!< Message consists of partial Energy#TotalEnergy results, sent during EnergyAllReduce()
#define OverlapTag 210  //!< Message consists of wave function coefficients sent during CalculatePerturbedOverlap()
#define PerturbedTag 211 //!< Message consists of wave function coefficients sent during Calculate1stPerturbedDerivative()
#define CurrentTag1 221 //!< Message consists of current density SpinType#SpinUp components sent during CalculateForce()
#define CurrentTag2 222 //!< Message consists of current density Spinype#SpinDown components sent during CalculateForce()
#define WannierCTag 230 //!< Wannier rotation exchange 
#define WannierSTag 231 //!< Wannier rotation exchange
#define WannierALTag 240  //!< Wannier index exchange 
#define WannierARTag 241  //!< Wannier index exchange


/** Options from command line.
 * This structure is generally filled with options from the command
 * line, scanned by \ref GetOptions()
 */
struct CallOptions {
  char* MainParameterFile;      //!< main parameter file
        char* ForcesFile;                                       //!< ForcesFile: is NULL if we don't parse forces from file and solve the ground state problem, otherwise don't solve and parse
  int debug;                    //!< (1) debug, (2) debug on errors
  int nicelevel;                //!< nice level for executable
  int Out;                                              //!< verbosity level
  int out[MaxOutGroup];           //!< print stderr-msg? for each group -v
  unsigned int alarm;           //!< set alarm after alarm seconds
  int proc[PROCDIM];                //!< process per gamma point, process per wave fct
  enum ParseSrcDensities ReadSrcFiles;  //!< whether to parse source coefficients from file or use minimisation
  int WriteSrcFiles;                //!< 0 - don't, 1 - write source file on exit with so far made calculations
  int AddNFactor;                             //!< ugly hack to fix read srcpsi
};

/** Structure containing file names and paths.
 * mainname, mainpath, path to config file, but also integers stating
 * whether outputfile (visual, measures) should be written and their
 * respective output file pointers
 */
struct FileData {
  char *mainname;        //!< full name of programme including path
  char *filename;        //!< name of programme
  char *mainpath;        //!< full path to programme
  char *default_path;    //!< path to default parameter file
  char *pseudopot_path;  //!< path ot pseudopotential files
  int *OutVisStep;       //!< Currently so and so often has visual data output occurred
  int MeOutVis;         //!< 1 - visual data are written to files, 0 - no output (for this process)
  int MeOutCurr;        //!< 1 - visual data of current density written to file, 0 - no output, no calculation (for this process)
  int MeOutMes;         //!< 1 - energy, forces and temperatures are written to files, 0 - no output (for this process)
  FILE *ForcesFile;                     //!< where the forces are written to
  FILE *EnergyFile;                     //!< where the energies are written to
  FILE *SpeedFile;                      //!< where the timing are written to
  FILE *HamiltonianFile;//!< where the explicit hamiltonian is written to (for retrieving Kohn-Sham-Eigenvalues)
  FILE *SpreadFile;     //!< where spread and wannier centers for each orbital are written to
  FILE *ReciSpreadFile; //!< where reciprocal spread for each orbital is written to
  FILE *MinimisationFile; //!< where TE, ATE, delta and various other elements from each minimisation steps are written to
  FILE *TemperatureFile;//!< where the temperatures are written to
  /* Arrays */
  int TotalSize;        //!< total number of coefficients (nodes) on all processes
  int LocalSizeR;       //!< local (total) number of real nodes, twice FileData::LocalSizeC
  int LocalSizeC;       //!< local (total) number of real nodes
  int MaxNUp;
  fftw_complex *PosC;   //!< complex coefficients array for the output of densities
  fftw_complex *PosTemp;
  fftw_complex *PosFactor;
  fftw_complex *work;   //!< working coefficients array for the output of densities
  fftw_real *PosR;      //!< real coefficients array for the output of densities
  enum ModeType *OutputPosType; //!< active - RiemannTensor is used, nonactive/standby - not
  int DoOutVis;         //!< 1 - visual data are written to files, 0 - no output (generally)
  int DoOutMes;         //!< 1 - energy and forces are written to files, 0 - no output (generally)
  int DoOutCurr;        //!< 1 - visual data of current density written to file, 0 - no output, no calculation
  int DoOutNICS;        //!< 1 - visual data of nuclear induced current shieldings written to file, 0 - no output, no calculation
  int DoOutOrbitals;    //!< 1 - output each orbital density, 0 - just total density
};

/** Structure containing info for MPI interface.
 * Such as process number, various communicators (interfacing groups in MPI), TID.
 * Always it's my (of this current process) process number within the group, within
 * the communicator, etc.
 */
struct ParallelSimulationData {
  int procs;                    //!< overall number of processes
  int proc[PROCDIM];            //!< number of processes per gamma point and per wave function
  int mypos[PROCDIM];           //!< I am the mypos-th process with the per gamma point or per wave function group
  int me;                       //!< my process id among all
  int mytid;                                                                            //!< My TID
  int partid;                   //!< Parent TID

  MPI_Comm world;               //!< Group: MPI_COMM_WORLD (including maybe callers)
  MPI_Comm comm;                //!< Group: all participating processe on this Problem at hand

  MPI_Comm comm_ST;             //!< Communicator for SpinDouble or SpinUpDown, thus either the whold world of one half of the world
  int me_comm_ST;               //!< my number within this communicator (going from 0 to Max_me_comm_ST-1)
  int my_color_comm_ST;         //!< either just 0 color(SpinDouble) or 0,1 (SpinUp - SpinDown)
  int Max_me_comm_ST;           //!< maximum number of processes in this SinUp/Down communicator, full number of total processes in SpinDouble, half the number in SpinUpDown
  int Max_my_color_comm_ST;     //!< maximum number regarding color (either 1 for SpinDouble or 2 for SpinUpDown)

  MPI_Comm comm_STInter;        //!< InterComm needed for SpinUp/Down case

  MPI_Comm comm_ST_Psi;         //!< SubComm for comm_ST Communicator Psi, these refer to the same wave function, sharing its coefficients on both grids
  int me_comm_ST_Psi;           //!< my number within this communicator
  int my_color_comm_ST_Psi;     //!< my number, either just one color(SpinDouble) or two (SpinUp - SpinDown)
  int Max_me_comm_ST_Psi;       //!< maximum number of processes in this Psi communicator
  int Max_my_color_comm_ST_Psi; //!< ??? maximum number regarding color (one or two)

  MPI_Comm comm_ST_PsiT;        //!< SubComm of comm_ST Communicator of processes that need to exchange wave funcions (GramSch())
  int me_comm_ST_PsiT;          //!< Transposed of comm_ST_Psi
  int my_color_comm_ST_PsiT;    //!< either just one color(SpinDouble) or two (SpinUp - SpinDown)
  int Max_me_comm_ST_PsiT;      //!< maximum number of processes in this GramSch communicator
  int Max_my_color_comm_ST_PsiT;//!< ??? maximum number regarding color (one or two)
};

struct RPlans {
  struct LevelPlan *plan;
  fftw_complex *cdata;
  fftw_real *rdata;
};

/** one reciprocal grid vector G.
 */
struct OneGData {
  int Index;       //!< index is needed to access elements in the upper Lev0, when densities are stored in reciprocal base, Density#DensityCArray
  int GlobalIndex;
  double GSq;              //!< squared norm of G vector
  double G[NDIM];        //!< integer cooefficients of reciprocal vector
};

struct GDataHash {
  int i;           //!< index running over each GArray from process zero to max.
  int myPE;        //!< to which process does it belong
};

/** Structure for the Psis per LatticeLevel.
 */
struct LevelPsi {
  fftw_complex **LocalPsi;  //!< are the current complex coefficients of Psis, stored locally in this process (OnePsiElement::LocalNo), initialised by PsiDat, used for orthonormlization of non-Psis
  fftw_complex **OldLocalPsi;//!< are the old complex coefficients of Psis, stored locally in this process (OnePsiElement::LocalNo), initialised by OldPsiDat. Needed to UpdateWavesAfterIonMove() and ComputeMLWF()
  fftw_complex *PsiDat;     //!< contains the initial complex coefficients of all wave functions (even temporary ones)
  fftw_complex *OldPsiDat;  //!< contains the old complex coefficients of all wave functions (even temporary ones)
  fftw_complex *TempPsi;    //!< temporal array used in GramSch(), size of LatticeLevel:MaxG, always the Psi to be orthogonalized
  fftw_complex *TempPsi2;   //!< pointer to a Psi used in GramSch(), always the Psi already orthogonalized
};

/** Structure containing one wave function.
 */
struct OnePsiElement {
  int me_comm_ST_Psi;       //!< this wave function belongs to this process in the communicator Psi
  int my_color_comm_ST_Psi; //!< this wave function belongs to this process in the communicator PsiT
  int MyLocalNo;            //!< continuing number among the Psis that are local to this process
  int MyGlobalNo;           //!< continuing global number among all Psis
  int/*enum PsiGramSchStatusType*/ PsiGramSchStatus;  //!< Status of how far the Gram-Schmidt-Orthonormalization has gone yet for this Psi
  enum MinimisationStatusType MinimisationStatus;     //!< Status of how far the Gram-Schmidt-Orthonormalization has gone yet for this Psi
  enum PsiTypeTag PsiType;                            //!< what type of wave function: normal (occupied), extra (gradient) or unoccupied Psi
  double PsiFactor;         //!< occupation number, hard-coded (2.0 in SpinDouble case or 1.0 in SpinUpDown case)
  double PsiReciNorm2;      //!< reciprocal norm of Psi being the Sum over all reciprocal grid vectors.
  int DoBrent;              //!< signals switching to brent iteration in line search
};

/** Additional Data to a wave function OnePsiElement.
 * Contains energy, kinetic eigenvalues, performed minimalization steps.
 */
struct OnePsiElementAddData {
  double Lambda;    //!< energy eigenvalue \f$\lambda_i = \langle \psi_i^{(m)}|H|\psi_i^{(m)} \rangle\f$
  double T;         //!< kinetic eigenvalue \f$\langle \psi_i | \ \frac{1}{2} \nabla^2 | psi_i \rangle\f$
  double Gamma;     //!< scalar product between the gradient of the wave function and the preconditioned one at this \ref Step
  int Step;         //!< holds count of currently made minimalization steps
  double WannierCentre[NDIM]; //!< NDIM coordinates of wannier centre of this orbital
  double WannierSpread;       //!< current spread of this orbital
};


/** Structure about all wave functions a.k.a orbits.
 * Containing variables such as SpinType, LocalNo, PsiGroup, OnePsiElement array and
 * various densities.
 */
struct Psis {
  enum UseSpinType Use;                  //!< Use Spin Up/Down or not (Double)
  enum SpinType PsiST;                   //!< Spin Type: Up, Down or Double
  int GlobalNo[MaxPsiNoType];//!< global number of Psis for the four plus one cases: Max, MaxDouble, MaxUp, MaxDown, MaxAdd
  int LocalNo;                                       //!< number of occupied and unoccupied Psis that are local (ly accessible) in this process
  int LocalNoAdd;            //!< number of unoccupied Psis that are local (ly accessible) in this process (spinup/-down or double)
  int NoOfPsis;              //!< contains SpinType-dependent sum of GlobalNo, giving number of wave functions in the current minimisation group (except for PsiTypeTag#UnOccupied of course)
  int NoOfTotalPsis;         //!< contains SpinType-dependent sum of GlobalNo with added PsiTypeTag#UnOccupied states
  int TypeStartIndex[Extra+2];//!< index array where the respective PsiTypeTag type starts in LocalPsiStatus
  int AllMaxLocalNo;        //!< maximum local number of Psis on one process (some processes have one more due to modulo != 0 on sharing)
  int MaxPsiGroup;          //!< number of processes among which the Psis are divided
  int PsiGroup;             //!< rank of this process in the communicator ParallelSimulationData::comm_ST_Psi
  int MaxPsiOfType;         //!< overall number of Psis
  int *AllLocalNo;          //!< array over all processes in the GramSch() communicator ParallelSimulationData::comm_ST_PsiT containing their number of local Psis, AllLocalNo[i] = RealAllLocalNo[i] + 1
  int *RealAllLocalNo;      //!< array over all processes in the GramSch() communicator ParallelSimulationData::comm_ST_PsiT containing their number of local Psis, i-th process: Psi->LocalNo = RealAllLocalNo[i]
  int MyStartNo;            //!< at which Psis do the locally accessible of this process start (going to MyStartNo+LocalNo)
  struct OnePsiElement *AllPsiStatus;         //!< array over all PsiTypeTag's of all wave functions, yet without coefficients (see Psis::LocalPsiStatus)
  struct OnePsiElement *AllPsiStatusForSort;  //!< temporary array of same size as *AllPsiStatus, only used in naturalmergesort()
  struct OnePsiElement *LocalPsiStatus;       //!< array over all PsiTypeTag's of the local wave functions which are accessible to/stored in this process
  int *TempSendA;           //!< In GramSch(): Holds count to which process in GramSch group a local wave function has been sent
  double NIDensity[Extra];                                //!< Density over all
  double NIDensityUp[Extra];                      //!< Density of SpinUp
  double NIDensityDown[Extra];                  //!< Density of SpinDown
  int *AllActualLocalPsiNo;
  int *AllOldActualLocalPsiNo;
  struct OnePsiElementAddData *AddData;   //!< some additional local data, such as energy eigenvalue
  double **lambda;            //!< contains \f$\lambda_{kl} = \langle \varphi_k^{(0)} | H^{(0)} | \varphi_l^{(0)} \rangle\f$, sa CalculateHamiltonian()
  double **Overlap;            //!< contains \f$S_{kl} = \langle \varphi_k^{(1)} | \varphi_l^{(1)} \rangle\f$, sa CalculatePerturbedOverlap()
};

/** Structure containing the various (radially discretized) densities.
 * Densities are always calculated on one level higher than the respective wave functions as the fftransformation then
 * becomes exact (fourier basis).
 */
struct Density {
  int TotalSize;    // total number of nodes
  int LocalSizeC;   // number of nodes in radial mesh of the complex density (locally accessible)
  int LocalSizeR;   // number of nodes in radial mesh of the density (locally accessible)
  enum UseType DensityTypeUse[MaxDensityTypes]; //!< Density is used or not (e.g. Up/Down not used in SpinType::SpinDouble case)
  fftw_real *DensityArray[MaxDensityTypes];     //!< density array (wave functions summed and squared, R-dependent) (radially discretized), \sa CalculateOneDensityR()
  fftw_complex *DensityCArray[MaxDensityTypes]; //!< density array (wave functions summed and squared, G-dependent), \sa CalculateOneDensityC()
  enum UseType LockArray[MaxDensityTypes];               //!< 1 - real density array is currently in use, 0 - array may be locked for sole
  enum UseType LockCArray[MaxDensityTypes];              //!< 1 - complex density array is currently in use, 0 - array may be locked for sole
};

enum complex {
              real,   //!< number is real
              imag    //!< number is complex
             };
/** Structure containing one lattice level.
 * Containing Level number LevelNo, maximum grid mesh points normal MaxN and reciprocal MaxG, factors
 * between this one and upper levels, array of the grid vectors GArray, wave function on this level
 * LPsi and densities Dens.
 */
struct LatticeLevel {
  int LevelNo;                //!< current number of this level
  int MaxN;                   //!< number of grid points
  int N[NDIM];                //!< number of grid points for this level in each dimension NDIM
  int NUp[NDIM];              //!< ratio of number of mesh points between this and the upper level for each dimension NDIM
  int NUp0[NDIM];             //!< ratio of number of mesh points this and the 0th level for each dimension NDIM
  int NUp1[NDIM];             //!< ratio of number of mesh points and the first level for each dimension NDIM
  int MaxNUp;                 //!< ratio of grid points between this and the upper level overall (all dimensions together)
  int MaxG;                                                                             //!< number of reciprocal grid vectors
  int MaxDoubleG;             //!< number of reciprocal grid vectors which are double by gamma point symmetry (only one is stored)
  int *AllMaxG;               //!< number of reciprocal grid vectors for each process
  int TotalAllMaxG;           //!< number of reciprocal grid (complex) vectors for all process
  int TotalRealAllMaxG;       //!< number of reciprocal grid (real) vectors for all process
  double ECut;                                                          //!< Cutoff energy (discretizes possible reciprocal vectors), here: maximum norm of reciprocal grid vector
  struct RPlans Plan0;
  struct OneGData *GArray;              //!< Array of reciprocal grid vectors
  struct GDataHash *HashG;    //!< Hash-array with global index for each G
  int *DoubleG;               //!< Array of indices of the doubly appearing grid vectors, with twice entries per vector: normal index and the inverse index (-x,-y,-z)
  int G0;                     //!< OneGData::Index of the reciprocal grid vector with zero components
  fftw_complex *PosFactorUp;  //!< precalculated position factor \f$P_p(G) \f$when going from wave coefficients to density, huge array of \ref MaxNUp times \ref MaxG, see CreatePosFacForG()
  struct LevelPsi *LPsi;
  struct Density *Dens;       //!< Density on this level
  int Step;                   //!< holds count of Molecular Dynamics steps
};


//! List of NField types, for this level, for the upper level
//! Enumerating entries in nFields for current and upper level
enum FFTNFields { FFTNF1,   //!< Number of nFields in current level
                  FFTNFUp   //!< Number of nFields in upper level
                };
enum FFTNFieldsS {                FFTNFSVecUp=2, 
                                  FFTNFSVec
                 };
enum FFTNFieldsR {                FFTNFRMatUp0=2,
                                  FFTNFRMatVecUpS, 
                                  FFTNFRVecUp0
                 };
enum FFTNFields0 {        FFTNF0Vec=1
                 };

#define MAXRTPOSFAC 2   //!< maximum number of different RiemannTensor position factors
//! Enumerating to which level the position factor relates
enum RTPosFacType { RTPFRtoS,   //!< position factor is relative to STANDARDLEVEL
                    RTPFRto0    //!< position factor is relative to topmost level
                   };

#define MAXRTARRAYS 8   //!< maximum number of different RiemannTensor arrays
//! Enumerating different RiemannTensor arrays
enum RTArrayType { RTADetPreRT, 
                   RTAPreA, 
                   RTAA, 
                   RTAIRT, 
                   RTARTA, 
                   RTAARTA, 
                   RTAiGcg, 
                   RTAcg 
                 };

struct RiemannTensor {
  int RiemannLevel;
  int NUpLevRS[NDIM];
  int NUpLevR0[NDIM];
  enum UseRiemannTensor Use;
  enum ModeType ActualUse;
  fftw_complex *Coeff;
  struct LatticeLevel *LevR; /* RiemannLevel */
  struct LatticeLevel *LevS; /* StandartLevel */
  struct LatticeLevel *Lev0; /*  0 Level */
  int MaxNUp[MAXRTPOSFAC];
  int TotalSize[MAXRTARRAYS];
  int LocalSizeC[MAXRTARRAYS];
  int LocalSizeR[MAXRTARRAYS];
  int NFields[MAXRTARRAYS];
  fftw_complex *PosFactor[MAXRTPOSFAC];
  fftw_complex *DensityC[MAXRTARRAYS];
  fftw_real *DensityR[MAXRTARRAYS];
  fftw_complex *RTLaplaceS;
  fftw_complex *RTLaplace0;
  fftw_complex *TempC;
  size_t TempTotalSize;
};

#define MAXALLPSIENERGY 4  //!< number of different wave function energies

//! Enumerating energy types of wave function
enum AllPsiEnergyTypes { KineticEnergy,   //!< kinetic energy
                         NonLocalEnergy,  //!< non-local pseudo potential energy
                         Perturbed1_0Energy,  //!< \f$\langle \varphi_l^{(1)} | H^{(1)} | \varphi_l^{(0)} \rangle \f$
                         Perturbed0_1Energy   //!< \f$\langle \varphi_l^{(0)} | H^{(1)} | \varphi_l^{(1)} \rangle \f$
                       }; 
#define MAXALLDENSITYENERGY 6  //!< number of different density energies
//! Enumerating density energy types
enum AllDensityEnergyTypes { CorrelationEnergy,       //!< correlation energy \f$E_C\f$
                             ExchangeEnergy,          //!< exchange energy \f$E_X\f$
                                               GaussEnergy,             //!< gaussian energy \f$E_{gauss}\f$
                             PseudoEnergy,            //!< pseudopotential energy \f$E_{ps}\f$
                                               HartreePotentialEnergy,  //!< Hartree potential energy including gaussian density \f$E_H\f$
                             HartreeEnergy           //!< Hartree energy \f$E_H\f$
                            };
#define MAXALLIONSENERGY 2  //!< number of different ion energies
//! Enumerating different ion energy types
enum AllIonsEnergyTypes { GaussSelfEnergy,  //!< Gaussian self energy due to charge
                          EwaldEnergy       //!< core to core energy summed over all super cells by ewald summation
                        }; 
#define MAXOLD 2                        //!< how many older calculated values are archived in various arrays

/** Energy structure.
 * Contains Total and local (i.e. within one process) energies of the
 * three classes: wave function, density and ion, also the first and
 * second time derivative, dating back to MAXOLD
 */
struct Energy {
  double AllTotalPsiEnergy[MAXALLPSIENERGY];    //!< Total energy in SpinType#SpinDouble case
  double AllLocalPsiEnergy[MAXALLPSIENERGY];    //!< Calculated energies of the local process, summed up via MPI_Allreduce to Energy::AllUpPsiEnergy, Energy::AllDownPsiEnergy or Energy::AllTotalPsiEnergy respectively
  double AllUpPsiEnergy[MAXALLPSIENERGY];       //!< Total energy in SpinType#SpinUp case
  double AllDownPsiEnergy[MAXALLPSIENERGY];     //!< Total energy in SpinType#SpinDown case
  double *PsiEnergy[MAXALLPSIENERGY];
  double AllLocalDensityEnergy[MAXALLDENSITYENERGY];  //!< local total energy of electron density within one process
  double AllTotalDensityEnergy[MAXALLDENSITYENERGY];  //!< Total energy resulting from electron density
  double AllTotalIonsEnergy[MAXALLIONSENERGY];  //!< Total energy of the Ions, Ewald and gaussian
  double TotalEnergy[MAXOLD];                   //!< Total energy as the sum of AllTotalPsiEnergy, AllTotalDensityEnergy and AllTotalIonsEnergy
  double TotalEnergyOuter[MAXOLD];
  double TotalEnergyFixed[MAXOLD];
  double delta[MAXOLD];
  double dEdt0[MAXOLD];
  double ddEddt0[MAXOLD];
  double ATE[MAXOLD];
  double parts[3];
  double bandgap;
  //double homolumo;
};

/** Lattice structure.
 * containing real, inverted, reciprocal basis (squared and/or orthonormal), Volume,
 * ECut, LatticeLevel structure, Psis structure , the fft_plan_3d structure and
 * Energy structure
 */
struct Lattice {
  double RealBasis[NDIM_NDIM];  //!< Coefficients of the basis vectors
  double RealBasisSQ[NDIM];     //!< squared Norm of each basis vector
  double RealBasisCenter[NDIM]; //!< center of the unit cell ((0.5,0.5,0.5) transformed by RealBasis)
  double InvBasis[NDIM_NDIM];   //!< Matrix-wise inverted basis vectors
  double ReciBasis[NDIM_NDIM];  //!< Coefficients of the transposed(!), inverse basis "matrix" (i.e. reciprocal basis)
  double ReciBasisSQ[NDIM];     //!< Norm of each reciprocal basis vectors
  double ReciBasisO[NDIM];      //!< Measure of how orthonormal each basis vectors is to the others
  double ReciBasisOSQ[NDIM];    //!< Square of the orthonormality measure ReciBasisO[]
  double Volume;                //!< volume as the determinant of the matrix of the basis vectors
  double ECut;                  //!< Energy cutoff, limiting number of reciprocal wave vectors
  double SawtoothStart;
  int Lev0Factor;               //!< LevelSizes of 0th level
  int LevRFactor;               //!< LevelSizes of (RL-1)th level (upper of RiemannLevel)
  int MaxLevel;                 //!< Maximum number of levels
  int AddNFactor;
  int *LevelSizes;              //!< Factor for all levels, stating the size in comparsion to the one below (most often factor 2)
  int *MaxNoOfnFields;          //!< maximum number of entries in NFields array per level
  int **NFields;
  struct LatticeLevel *Lev;     //!< Array of different LatticeLevels, up to MaxLevel
  struct OneGData *GArrayForSort;
  struct Psis Psi;              //!< Wave function structure residing within the lattice
  struct RiemannTensor RT;
  struct fft_plan_3d *plan;
  struct Energy Energy[Extra];  //!< Energy structure for this particular lattice setting
  struct Energy *E;             //!< pointer to current Energy structure within minimisation scheme
};

#define MAXTIMETYPES 19                 //!< specifying number of TimeTypes
//! Enumerating timing groups, each having its own timer \warning LevSMaxG must be the last one */
enum TimeTypes { SimTime,           //!< simulation time
                 InitSimTime,       //!< simulation time during which initialisation of calculations is performed
                 InitTime,          //!< initialization time at programme start
                             InitGramSchTime,   //!< initialization time for Gram-Schmidt-Orthogonalization
                 InitLocTime,       //!< initialization time for local pseudopotential calculations
                 InitNonLocTime,    //!< initialization time for non-local pseudopotential calculations
                 InitDensityTime,   //!< initialization time for density calculations                 
                             GramSchTime,       //!< time spent orthonormalizing the orbits
                 LocTime,           //!< time spent evaluating local pseudopotential
                 NonLocTime,        //!< time spent evaluating non-local pseudopotential
                 DensityTime,       //!< time spent evaluating density
                             LocFTime,          //!< time spent evaluating local form factors
                 NonLocFTime,       //!< time spent evaluating non-local form factors
                 EwaldTime,         //!< time spent evaluating ewald summation over ion core-to-core forces
                 GapTime,           //!< time spent evaluating the gap energy
                 CurrDensTime,      //!< time spent evaluating the current density
                 WannierTime,       //!< time spent localizing the orbitals
                 ReadnWriteTime,    //!< time spent during reading and writing of wave function coefficients
                             LevSMaxG};         //!< number of reciprocal grid vectors to calculate standard deviation, must always be the last one!

//! enumerating if time is started or stopped
enum TimeDoTypes { StartTimeDo,   //!< start timer
                   StopTimeDo     //!< stop timer
                 };


/** Array structure for timekeeping.
 * It consists of arrays for each timing group(!) for
 * steps, averages, stored times, deviations, min and max
 */
struct SpeedStruct {
  double SpeedStep[MAXTIMETYPES];
  double time1[MAXTIMETYPES];                           //!< time before the time before current time
  double time2[MAXTIMETYPES];                           //!< time before current time
  double time[MAXTIMETYPES];                            //!< current time
  double average[MAXTIMETYPES];                 //!< storing average
  double stddev[MAXTIMETYPES];                  //!< standard deviation
  double min[MAXTIMETYPES];                                     //!< minimum time
  double max[MAXTIMETYPES];                                     //!< maximum time
  int InitSteps;
  int LevUpSteps;
  int Steps;                        //!< holds count of made Molecular dynamic steps
};

struct Problem;
#include "ions.h"
#include "pseudo.h"
#include "excor.h"
#include "grad.h"
#include "run.h"

/** Physical problem superstructure.
 * This structure contains all the vital information specifying the
 * problem at hand, such as the command line and parameter file options,
 * also very important RunStruct for mpi groups and SpeedStruct for time measuring,
 * or simply all other structures.
 * \sa GetOptions() and ReadParameters()
 */
struct Problem {
  struct CallOptions Call;
  struct FileData Files;
  struct ParallelSimulationData Par;
  struct Lattice Lat;
  struct Ions Ion;
  struct PseudoPot PP;
  struct ExCor ExCo;
  struct Gradient Grad;
  struct RunStruct R;
  struct SpeedStruct Speed;
};


#endif