Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

source: pcp/src/init.c@ 64ca279

Visit:

Last change on this file since 64ca279 was a0bcf1, checked in by Frederik Heber <heber@…>, 17 years ago
-initial commit -Minimum set of files needed from ESPACK SVN repository -Switch to three tantamount package parts instead of all relating to pcp (as at some time Ralf's might find inclusion as well)
Property mode set to `100644`
File size: 86.4 KB

Rev	Line
[a0bcf1]	1	/** \file init.c
	2	* Initialization and Readin of parameter files.
	3	* Initialization Init() is split up into various parts:
	4	* lattice InitLattice(), first lattice level InitLatticeLevel0(), remaining lattice level InitOtherLevel(),
	5	* wave functions InitPsis() and their values InitPsisValue() and distribution among the processes InitPsiGroupNo(),
	6	* Many parts are found in their respective files such as for densities density.c and energies or
	7	* GramSchmidt gramsch.c.\n
	8	* Reading of the parameter file is done in ReadParameters() by parsing the file with ParseForParameter(), CreateMainname()
	9	* separates the path from the main programme file.\n
	10	* There remain some special routines used during initialization: splitting grid vector index on level change SplitIndex(),
	11	* calculating the grid vector index GetIndex(), sort criteria on squared grid vector norm GetKeyOneG(), testing whether G's
	12	* norm is below the energy cutoff TestG, or discerning wethers it's stored only half as often due to symmetry TestGDouble().
	13	* Also due to density being calculated on a higher level there are certain factors that can be pulled out which can be
	14	* precalculated CreatePosFacForG().
	15	*
	16	Project: ParallelCarParrinello
	17	\author Jan Hamaekers, Frederik Heber
	18	\date 2000, 2006
	19
	20	File: init.c
	21	$Id: init.c,v 1.97 2007-04-10 14:34:08 foo Exp $
	22	*/
	23
	24	#include <stdlib.h>
	25	#include <stdio.h>
	26	#include <math.h>
	27	#include <string.h>
	28
	29	// use double precision fft when we have it
	30	#ifdef HAVE_CONFIG_H
	31	#include <config.h>
	32	#endif
	33
	34	#ifdef HAVE_DFFTW_H
	35	#include "dfftw.h"
	36	#else
	37	#include "fftw.h"
	38	#endif
	39
	40	#include "data.h"
	41	#include "helpers.h"
	42	#include "errors.h"
	43	#include "init.h"
	44	#include "mergesort2.h"
	45	#include "myfft.h"
	46	#include "gramsch.h"
	47	#include "density.h"
	48	#include "riemann.h"
	49	#include "output.h"
	50	#include "energy.h"
	51	#include "mymath.h"
	52	#include "pseudo.h"
	53	#include "ions.h"
	54
	55	/** Splits the index of a reciprocal grid vector by the number of nodes per dimension.
	56	* \param Index current index of the grid vector
	57	* \param nx number of nodes in x direction
	58	* \param nz number of nodes in z direction
	59	* \param *x return value (Index % nx)
	60	* \param *y return value (Index / nx / nz)
	61	* \param *z return value (Index / nx % nz)
	62	*/
	63	static void SplitIndex(int Index, const int nx, const int nz, int x, int y, int *z) {
	64	if ((nx == 0) \|\| (nz == 0)) fprintf(stderr,"SplitIndex: nx = %d, nz == %d",nx,nz);
	65	*x = Index % nx;
	66	Index /= nx;
	67	*z = Index % nz;
	68	*y = Index / nz;
	69	}
	70
	71	/** Prepare mainname (and mainpath).
	72	* Extracts the path specified from the config file on command line by hopping along the path delimiters,
	73	* this becomes Files::mainpath, including the filename itself it becomes Files::mainname.
	74	* \param *P Problem at hand
	75	* \param filename pointer to file name string
	76	*/
	77	void CreateMainname(struct Problem const P, const char const filename)
	78	{
	79	char *dummy;
	80	char *token;
	81	char *token2;
	82	const char delimiters[] = "/";
	83	int stop = 0;
	84	int setit = 0;
	85	dummy = (char) / Kopie von filename */
	86	Malloc(strlen(filename)+1, "PrecreateMainname - dummy");
	87	sprintf(dummy,"%s",filename);
	88	P->Files.mainname = (char) / Speicher fuer P->Files.mainname */
	89	Malloc(strlen(filename)+strlen(P->Files.filename)+2, "CreateMainname - P->Files.mainname");
	90	P->Files.mainname[0] = 0;
	91	/* Jetzt mit strtok mainname vorbereiten */
	92	if (dummy) {
	93	if (dummy[0] == '/') setit = 1; // remember there was delimiter at the very beginning
	94	token = strtok(dummy, delimiters);
	95	do {
	96	token2 = strtok(0, delimiters); // a subsequent call of strtok to the end of this path piece (between token and token2)
	97	if (token2) { // if it's not the last piece
	98	if (setit) strcat(P->Files.mainname,"/");
	99	setit = 1; // afterwards always replace the delimiters in the target string
	100	strcat(P->Files.mainname, token); // copy the path piece
	101	token = token2; // and go to next section
	102	} else {
	103	stop = 1;
	104	}
	105	} while (!stop);
	106	}
	107	if (setit) strcat(P->Files.mainname,"/");
	108	if (dummy) Free(dummy);
	109	P->Files.mainpath = (char) / Speicher fuer P->Files.mainname */
	110	Malloc(strlen(P->Files.mainname)+10, "CreateMainname - P->Files.mainpath");
	111	sprintf(P->Files.mainpath, "%s" ,P->Files.mainname);
	112	if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) CreateMainame: mainpath: %s\t", P->Par.me, P->Files.mainpath);
	113	//strcat(P->Files.mainname, P->Files.filename); /* mainname fertigstellen */
	114	strcpy(P->Files.mainname, P->Files.filename);
	115	if (P->Call.out[ReadOut]) fprintf(stderr,"mainname: %s\n", P->Files.mainname);
	116	}
	117
	118	/** Initializing Lattice structure.
	119	* Calculates unit cell's volume, (transposed) reciprocal base Lattice::ReciBasis
	120	* and square-rooted Lattice::ReciBasisSQ, preparation of LatticeLevel structure,
	121	* NFields depends on use of RiemannTensor, calculates the orthonormal reciprocal
	122	* base and thus how many mesh points for the ECut constraint are at least necessary.
	123	* \param *P Problem at hand
	124	* \return 0 - everything went alright, 1 - some problem (volume is zero)
	125	*/
	126	static int InitLattice(struct Problem *const P) {
	127	struct Lattice *const Lat = &P->Lat;
	128	double factor = 2.*PI;
	129	double h[NDIM_NDIM];
	130	double V[NDIM];
	131	double NormV;
	132	int AddNFactor = Lat->AddNFactor = P->Call.AddNFactor;
	133	int i,j,UpDummy;
	134	int MaxG[NDIM], N[NDIM], NFactor, Lev1N[NDIM];
	135	struct LatticeLevel *Lev0;
	136	int RL=-1;
	137	// Volume
	138	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitLattice\n", P->Par.me);
	139	Lat->Volume = RDET3(Lat->RealBasis);
	140	if (P->Call.out[ValueOut]) fprintf(stderr,"Volume = %g\n",Lat->Volume);
	141	if (Lat->Volume == 0.0) {
	142	return(1);
	143	}
	144	// square-rooted and normal reciprocal (transposed) Base
	145	for (i=0; i < NDIM_NDIM; i++) h[i] = Lat->RealBasis[i];
	146	RTranspose3(h);
	147	RMatReci3(Lat->ReciBasis, h);
	148	for (i=0; i < NDIM_NDIM; i++) Lat->ReciBasis[i] *= factor; // SM(Lat->ReciBasis, factor, NDIM_NDIM);
	149	for (i=0; i < NDIM; i++) {
	150	Lat->ReciBasisSQ[i] = RNORMSQ3(&(Lat->ReciBasis[i*NDIM]));
	151	Lat->RealBasisSQ[i] = RNORMSQ3(&(Lat->RealBasis[i*NDIM]));
	152	Lat->RealBasisQ[i] = sqrt(Lat->RealBasisSQ[i]);
	153	}
	154	// Prepares LatticeLevels
	155	Lat->Lev = (struct LatticeLevel *)
	156	Malloc(Lat->MaxLevel*sizeof(struct LatticeLevel),"InitLattice: Lat->Lev");
	157	for (i=0; i < Lat->MaxLevel; i++)
	158	Lat->Lev[i].Step = 0;
	159	Lev0 = Lat->Lev; // short for &(Lat->Lev[0])
	160	Lat->LevelSizes = (int *)
	161	Malloc(sizeof(int)*Lat->MaxLevel,"InitLattice: LevelSizes");
	162	Lat->LevelSizes[0] = P->Lat.Lev0Factor;
	163	for (i=1; i < Lat->MaxLevel-1; i++) Lat->LevelSizes[i] = 2;
	164	if (Lat->RT.Use == UseRT) {
	165	RL = Lat->RT.RiemannLevel;
	166	Lat->LevelSizes[RL-1] = P->Lat.LevRFactor;
	167	}
	168	Lat->LevelSizes[Lat->MaxLevel-1] = 1;
	169	// NFields are next ... depending on use of RiemannTensor
	170	Lat->MaxNoOfnFields = (int *)
	171	Malloc(sizeof(int)*Lat->MaxLevel,"InitLattice: MaxNoOfnFields");;
	172	Lat->NFields = (int **)
	173	Malloc(sizeof(int)Lat->MaxLevel,"InitLattice: NFields");
	174	switch (Lat->RT.Use) {
	175	case UseNotRT:
	176	for (i=0; i<Lat->MaxLevel; i++) {
	177	if (i==0) {
	178	Lat->MaxNoOfnFields[i] = 1;
	179	Lat->NFields[i] = (int *)
	180	Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
	181	Lat->NFields[i][FFTNF1] = 1;
	182	} else {
	183	Lat->MaxNoOfnFields[i] = 2;
	184	Lat->NFields[i] = (int *)
	185	Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
	186	Lat->NFields[i][FFTNF1] = 1;
	187	Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[i-1]Lat->LevelSizes[i-1]Lat->LevelSizes[i-1];
	188	}
	189	}
	190	break;
	191	case UseRT:
	192	for (i=0; i<Lat->MaxLevel; i++) {
	193	switch (i) {
	194	case 0:
	195	Lat->MaxNoOfnFields[i] = 2;
	196	Lat->NFields[i] = (int *)
	197	Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
	198	Lat->NFields[i][FFTNF1] = 1;
	199	Lat->NFields[i][FFTNF0Vec] = NDIM;
	200	break;
	201	case STANDARTLEVEL:
	202	Lat->MaxNoOfnFields[i] = 4;
	203	Lat->NFields[i] = (int *)
	204	Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
	205	Lat->NFields[i][FFTNF1] = 1;
	206	Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[0]Lat->LevelSizes[0]Lat->LevelSizes[0];
	207	Lat->NFields[i][FFTNFSVecUp] = Lat->NFields[i][FFTNFUp]*NDIM;
	208	Lat->NFields[i][FFTNFSVec] = NDIM;
	209	break;
	210	default:
	211	if (i == RL) {
	212	Lat->MaxNoOfnFields[i] = 5;
	213	Lat->NFields[i] = (int *)
	214	Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
	215	Lat->NFields[i][FFTNF1] = 1;
	216	Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[i-1]Lat->LevelSizes[i-1]Lat->LevelSizes[i-1];
	217	UpDummy = 1;
	218	for (j=RL-1; j >= STANDARTLEVEL; j--)
	219	UpDummy = Lat->LevelSizes[j]Lat->LevelSizes[j]*Lat->LevelSizes[j];
	220	/* UpDummy RL -> S */
	221	Lat->NFields[i][FFTNFRMatVecUpS] = UpDummyNDIM_NDIMNDIM;
	222
	223	UpDummy = 1;
	224	for (j=RL-1; j >= 0; j--)
	225	UpDummy = Lat->LevelSizes[j]Lat->LevelSizes[j]*Lat->LevelSizes[j];
	226	/* UpDummy RL -> 0 */
	227	Lat->NFields[i][FFTNFRMatUp0] = UpDummy*NDIM_NDIM;
	228	Lat->NFields[i][FFTNFRVecUp0] = UpDummy*NDIM;
	229	} else {
	230	Lat->MaxNoOfnFields[i] = 2;
	231	Lat->NFields[i] = (int *)
	232	Malloc(sizeof(int)*Lat->MaxNoOfnFields[i],"InitLattice: NFields");
	233	Lat->NFields[i][FFTNF1] = 1;
	234	Lat->NFields[i][FFTNFUp] = Lat->LevelSizes[i-1]Lat->LevelSizes[i-1]Lat->LevelSizes[i-1];
	235	}
	236	break;
	237	}
	238	}
	239	break;
	240	}
	241	if (P->Call.out[ValueOut]) {
	242	for (i=0; i<Lat->MaxLevel; i++) {
	243	fprintf(stderr,"NFields[%i] =",i);
	244	for (j=0; j<Lat->MaxNoOfnFields[i]; j++) {
	245	fprintf(stderr," %i",Lat->NFields[i][j]);
	246	}
	247	fprintf(stderr,"\n");
	248	}
	249	}
	250	// calculate mesh points per dimension from ECut constraints for 0th level
	251	for (i=0; i< NDIM; i++) { // by vector product creates orthogonal reciprocal base vectors, takes norm, ...
	252	VP3(V, &Lat->ReciBasis[NDIM((i+1)%NDIM)], &Lat->ReciBasis[NDIM((i+2)%NDIM)]);
	253	NormV = sqrt(RNORMSQ3(V));
	254	if (fabs(NormV) < MYEPSILON) fprintf(stderr,"InitLattice: NormV = %lg",NormV);
	255	for (j=0;j<NDIM;j++) V[j] /= NormV;
	256	Lat->ReciBasisO[i] = fabs(RSP3(V,&Lat->ReciBasis[NDIM*i]));
	257	Lat->ReciBasisOSQ[i] = Lat->ReciBasisO[i]*Lat->ReciBasisO[i];
	258	if (P->Call.out[ValueOut]) fprintf(stderr,"G[%i] = (%g ,%g ,%g) GSQ[%i] = %g GOSQ[%i] = %g\n",i,Lat->ReciBasis[iNDIM+0],Lat->ReciBasis[iNDIM+1],Lat->ReciBasis[i*NDIM+2],i,Lat->ReciBasisSQ[i],i,Lat->ReciBasisOSQ[i]);
	259	if (Lat->ReciBasisOSQ[i] < MYEPSILON) fprintf(stderr,"InitLattice: Lat->ReciBasisOSQ[i] = %lg",Lat->ReciBasisOSQ[i]);
	260	MaxG[i] = floor(sqrt(2.*Lat->ECut/Lat->ReciBasisOSQ[i]));
	261	N[i] = 2*MaxG[i]+2;
	262	NFactor = ( i==2 ? 2 : 1)P->Par.proc[PEGamma]; / TESTINIT */
	263	for (j=1; j < Lat->MaxLevel; j++) NFactor *= Lat->LevelSizes[j];
	264	if (AddNFactor == 0) fprintf(stderr,"InitLattice: AddNFactor = %d",AddNFactor);
	265	NFactor *= AddNFactor;
	266	Lev0->N[i] = NFactor*ceil((double)N[i]/(double)NFactor);
	267	Lev1N[i] = Lev0->N[i];
	268	Lev0->N[i] *= Lat->LevelSizes[0];
	269	}
	270	if (P->Call.out[LeaderOut] && P->Par.me == 0) fprintf(stderr,"MaxG[] = %i %i %i -> N[] = %i %i %i -> Lev1->N[] %i %i %i -> Lev0->N[] %i %i %i\n",MaxG[0],MaxG[1],MaxG[2],N[0],N[1],N[2],Lev1N[0],Lev1N[1],Lev1N[2],Lev0->N[0],Lev0->N[1],Lev0->N[2]);
	271	if (P->Files.MeOutMes)
	272	fprintf(P->Files.SpeedFile,"\nMaxG[] = %i %i %i -> N[] = %i %i %i -> Lev1->N[] %i %i %i -> Lev0->N[] %i %i %i\n",MaxG[0],MaxG[1],MaxG[2],N[0],N[1],N[2],Lev1N[0],Lev1N[1],Lev1N[2],Lev0->N[0],Lev0->N[1],Lev0->N[2]);
	273	Lat->E = &Lat->Energy[Occupied];
	274
	275	if (P->Call.out[ValueOut]) {
	276	fprintf(stderr,"(%i) Lat->Realbasis = \n",P->Par.me);
	277	for (i=0; i < NDIM_NDIM; i++) {
	278	fprintf(stderr,"%lg ",Lat->RealBasis[i]);
	279	if (i%NDIM == NDIM-1) fprintf(stderr,"\n");
	280	}
	281	fprintf(stderr,"(%i) Lat->RealbasisQ = ",P->Par.me);
	282	for (i=0; i < NDIM; i++)
	283	fprintf(stderr,"%lg ",Lat->RealBasisQ[i]);
	284	fprintf(stderr,"\n");
	285	}
	286
	287	return(0);
	288	}
	289
	290	/** Test if a reciprocal grid vector is still within the energy cutoff Lattice::ECut.
	291	* The vector is generated with the given coefficients, its norm calculated and compared
	292	* with twice the ECut.
	293	* \param *Lat Lattice containing the reciprocal basis vectors and ECut
	294	* \param g1 first coefficient
	295	* \param g2 second coefficient
	296	* \param g3 third coefficient
	297	* \param *GSq return value for norm of the vector
	298	* \param G return array of NDIM for the vector
	299	* \return 0 - vector ok, 1 - norm is bigger than twice the energy cutoff
	300	*/
	301	static int TestG(struct Lattice Lat, int g1, int g2, int g3, double GSq, double G[NDIM]) {
	302	G[0] = g1Lat->ReciBasis[0]+g2Lat->ReciBasis[3]+g3*Lat->ReciBasis[6];
	303	G[1] = g1Lat->ReciBasis[1]+g2Lat->ReciBasis[4]+g3*Lat->ReciBasis[7];
	304	G[2] = g1Lat->ReciBasis[2]+g2Lat->ReciBasis[5]+g3*Lat->ReciBasis[8];
	305	*GSq = RNORMSQ3(G);
	306	if (0.5(GSq) <= Lat->ECut) return (1);
	307	return (0);
	308	}
	309
	310	/** Specifies which DoubleGType a certain combination of coefficients is.
	311	* Basically, on the (gz=0)-plane the gy half is DoubleGUse, the lower DoubleGNotUse (including
	312	* the (gx>=0)-line)
	313	* \param gx x coefficient of reciprocal grid vector
	314	* \param gy y coefficient of reciprocal grid vector
	315	* \param gz z coefficient of reciprocal grid vector
	316	* \return DoubleG0 : if all coefficients are zero <br>
	317	* DoubleGUse : if gz is zero, and gx is strictly positive and gy positive or vice versa <br>
	318	* DoubleGNotUse: if gz is zero, and gx is (strictly) positive and gy (strictly) negative or vice versa <br>
	319	* DoubleGNot : if gz is not equal to zero
	320	*/
	321	static enum DoubleGType TestGDouble(int gx, int gy, int gz) {
	322	if (gx == 0 && gy == 0 && gz == 0) return DoubleG0;
	323	if (gz == 0) {
	324	if (gx > 0) {
	325	if (gy >= 0) return DoubleGUse;
	326	else return DoubleGNotUse;
	327	} else {
	328	if (gy > 0) return DoubleGUse;
	329	else return DoubleGNotUse;
	330	}
	331	}
	332	return DoubleGNot;
	333	}
	334
	335	/** Returns squared product of the reciprocal vector as a sort criteria.
	336	* \param *a OneGData structure array containing the grid vectors
	337	* \param i i-th reciprocal vector
	338	* \param *Args not used, specified for contingency
	339	* \return a[i].OneGData::GSq
	340	*/
	341	static double GetKeyOneG(void a, int i, void Args) {
	342	return(((struct OneGData *)a)[i].GSq);
	343	}
	344
	345	/** Returns component of reciprocal vector as a sort criteria.
	346	* \param *a OneGData structure array containing the grid vectors
	347	* \param i i-th reciprocal vector
	348	* \param *Args specifies index
	349	* \return a[i].OneGData::G[Args]
	350	*/
	351	static double GetKeyOneG2(void a, int i, void Args) {
	352	//fprintf(stderr,"G[%i] = %e\n", (int )Args,((struct OneGData )a)[i].G[(int *)Args] );
	353	return(((struct OneGData )a)[i].G[(int *)Args]);
	354	}
	355
	356	/** Copies each entry of OneGData from \a b to \a a.
	357	* \param *a destination reciprocal grid vector
	358	* \param i i-th grid vector of \a a
	359	* \param *b source reciprocal grid vector
	360	* \param j j-th vector of \a b
	361	*/
	362	static void CopyElementOneG(void a, int i, void b, int j) {
	363	((struct OneGData )a)[i].Index = ((struct OneGData )b)[j].Index;
	364	((struct OneGData )a)[i].GlobalIndex = ((struct OneGData )b)[j].GlobalIndex;
	365	((struct OneGData )a)[i].GSq = ((struct OneGData )b)[j].GSq;
	366	((struct OneGData )a)[i].G[0] = ((struct OneGData )b)[j].G[0];
	367	((struct OneGData )a)[i].G[1] = ((struct OneGData )b)[j].G[1];
	368	((struct OneGData )a)[i].G[2] = ((struct OneGData )b)[j].G[2];
	369	}
	370
	371	/** Calculates the index of a reciprocal grid vector by its components.
	372	* \param *plan
	373	* \param *LPlan LevelPlan structure
	374	* \param gx first component
	375	* \param gy second component
	376	* \param gz third component
	377	* \return index
	378	*/
	379	static int GetIndex(struct fft_plan_3d plan, struct LevelPlan LPlan, int gx, int gy, int gz) {
	380	int x = (gx < 0 ? gx+LPlan->N[0] : gx);
	381	int y = (gy < 0 ? gy+LPlan->N[1] : gy);
	382	int Y = y/LPlan->LevelFactor;
	383	int YY = y%LPlan->LevelFactor;
	384	int z = (gz < 0 ? gz+LPlan->N[2]/2+1 : gz);
	385	int Z = z/LPlan->LevelFactor;
	386	int ZZ = z%LPlan->LevelFactor;
	387	int pwIndex = Z+(plan->NLSR[2]/2+1)*Y;
	388	int i=0,pwy,pwz,iz;
	389	for (i=0; i < plan->LocalMaxpw; i++)
	390	if (plan->Localpw[i].Index == pwIndex) {
	391	pwz = i%(plan->NLSR[2]/2+1);
	392	if (pwz == plan->NLSR[2]/2) {
	393	iz = LPlan->N[2]/2;
	394	} else {
	395	iz = pwz*LPlan->LevelFactor+ZZ;
	396	}
	397	pwy = i/(plan->NLSR[2]/2+1);
	398	return (x+LPlan->N[0](iz+LPlan->nz(pwy*LPlan->LevelFactor+YY)));
	399	}
	400	return -1;
	401	}
	402
	403	/** Builds a hash table and sets a global index for GArray.
	404	* The global index is needed in order to initalise the wave functions randomly but solely dependent
	405	* on the seed value, not on the number of coefficient sharing processes!
	406	* \param *P Problem at hand
	407	* \param *Lat Pointer to Lattice structure
	408	* \param *Lev Pointer to LatticeLevel structure, containing GArray
	409	*/
	410	static void HashGArray(struct Problem P, struct Lattice Lat, struct LatticeLevel *Lev)
	411	{
	412	int index, i, g, dupes;
	413	int myPE, MaxPE;
	414	struct OneGData GlobalGArray, GlobalGArray_sort;
	415	double G, G_send;
	416	int base = 0;
	417	int MaxG = 0;
	418	int AllMaxG = 0;
	419	MPI_Comm_size(P->Par.comm_ST_Psi, &MaxPE); // Determines the size of the group associated with a communictor
	420	MPI_Comm_rank(P->Par.comm_ST_Psi, &myPE); // Determines the rank of the calling process in the communicator
	421
	422	// first gather GArray from all coefficient sharing processes into one array
	423	// get maximum number of G in one process
	424	for (i=0;i<MaxPE; i++) { // go through every process
	425	if (MaxG < Lev->AllMaxG[i])
	426	MaxG = Lev->AllMaxG[i];
	427	AllMaxG += Lev->AllMaxG[i];
	428	}
	429	GlobalGArray = (struct OneGData ) Malloc(AllMaxG(sizeof(struct OneGData)),"HashGArray: GlobalGArray");
	430	GlobalGArray_sort = (struct OneGData ) Malloc(AllMaxG(sizeof(struct OneGData)),"HashGArray: GlobalGArray");
	431	G = (double *) Malloc(NDIM(sizeof(double)),"HashGArray: G");
	432	G_send = (double *) Malloc(NDIM(sizeof(double)),"HashGArray: G_send");
	433	for (i=0;i<NDIM;i++) {
	434	G[i] = (double ) Malloc(MaxPEMaxG*(sizeof(double)),"HashGArray: G[i]");
	435	G_send[i] = (double ) Malloc(MaxG(sizeof(double)),"HashGArray: G_send[i]");
	436	}
	437	// fill send array
	438	for (index=0;index<NDIM;index++) {
	439	for (i=0;i<Lev->MaxG;i++)
	440	G_send[index][i] = Lev->GArray[i].G[index];
	441	for (;i<MaxG;i++)
	442	G_send[index][i] = 0.;
	443	}
	444
	445	// gather among all processes
	446	MPI_Allgather(G_send[0], MaxG, MPI_DOUBLE, G[0], MaxG, MPI_DOUBLE, P->Par.comm_ST_Psi);
	447	MPI_Allgather(G_send[1], MaxG, MPI_DOUBLE, G[1], MaxG, MPI_DOUBLE, P->Par.comm_ST_Psi);
	448	MPI_Allgather(G_send[2], MaxG, MPI_DOUBLE, G[2], MaxG, MPI_DOUBLE, P->Par.comm_ST_Psi);
	449	// fill gathered results into GlobalGArray
	450	//if (P->Call.out[LeaderOut])
	451	fprintf(stderr,"(%i) MaxPE %i / myPE %i, AllMaxG %i / MaxG %i / Lev->MaxG %i, TotalAllMaxG %i\n", P->Par.me, MaxPE, myPE, AllMaxG, MaxG, Lev->MaxG, Lev->TotalAllMaxG);
	452	base = 0;
	453	for (i=0;i<MaxPE; i++) { // go through every process
	454	for (g=0;g<Lev->AllMaxG[i];g++) { // through every G of above process
	455	for (index=0;index<NDIM;index++)
	456	GlobalGArray[base+g].G[index] = G[index][MaxG*i+g];
	457	GlobalGArray[base+g].GlobalIndex = i;
	458	GlobalGArray[base+g].Index = g;
	459	//if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray[%i+%i].G[index] = (%e,%e,%e), .index = %i, myPE = %i\n",P->Par.me, base, g, GlobalGArray[base+g].G[0], GlobalGArray[base+g].G[1], GlobalGArray[base+g].G[2], GlobalGArray[base+g].Index, GlobalGArray[base+g].GlobalIndex);
	460	}
	461	base += Lev->AllMaxG[i];
	462	}
	463
	464	// then sort array
	465	// sort by x coordinate
	466	index = 0;
	467	naturalmergesort(GlobalGArray,GlobalGArray_sort,0,AllMaxG-1,&GetKeyOneG2,&index,&CopyElementOneG);
	468	//if (P->Call.out[LeaderOut]) fprintf(stderr,"\n\n");
	469	base = 0;
	470	//for (g=0;g<AllMaxG;g++) // through every G of above process
	471	//if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray_sort[%i].G[index] = (%e,%e,%e), .index = %i, myPE= %i\n",P->Par.me, g, GlobalGArray_sort[g].G[0], GlobalGArray_sort[g].G[1], GlobalGArray_sort[g].G[2], GlobalGArray_sort[g].Index, GlobalGArray[base+g].GlobalIndex);
	472	// sort by y coordinate
	473	//if (P->Call.out[LeaderOut]) fprintf(stderr,"\n\n");
	474	index = 1;
	475	base = 0;
	476	for (g=0;g<AllMaxG;g++) { // through every G of above process
	477	if (GlobalGArray[base].G[index-1] != GlobalGArray[g].G[index-1]) {
	478	naturalmergesort(GlobalGArray_sort,GlobalGArray,base,g-1,&GetKeyOneG2,&index,&CopyElementOneG);
	479	base = g;
	480	}
	481	}
	482	if (base != AllMaxG-1) naturalmergesort(GlobalGArray_sort,GlobalGArray,base,AllMaxG-1,&GetKeyOneG2,&index,&CopyElementOneG); // last section
	483	base = 0;
	484	//for (g=0;g<AllMaxG;g++) // through every G of above process
	485	//if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray[%i].G[index] = (%e,%e,%e), .index = %i, myPE= %i\n",P->Par.me, g, GlobalGArray[g].G[0], GlobalGArray[g].G[1], GlobalGArray[g].G[2], GlobalGArray[g].Index, GlobalGArray[base+g].GlobalIndex);
	486	// sort by z coordinate
	487	index = 2;
	488	base = 0;
	489	for (g=0;g<AllMaxG;g++) { // through every G of above process
	490	if (GlobalGArray[base].G[index-1] != GlobalGArray[g].G[index-1]) {
	491	naturalmergesort(GlobalGArray,GlobalGArray_sort,base,g-1,&GetKeyOneG2,&index,&CopyElementOneG);
	492	base = g;
	493	}
	494	}
	495	if (base != AllMaxG-1) naturalmergesort(GlobalGArray,GlobalGArray_sort,base,AllMaxG-1,&GetKeyOneG2,&index,&CopyElementOneG); // last section
	496	base = 0;
	497	//for (g=0;g<AllMaxG;g++) // through every G of above process
	498	//if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) GlobalGArray_sort[%i].G = (%e,%e,%e), myPE = %i, index = %i\n", P->Par.me, g, GlobalGArray_sort[g].G[0], GlobalGArray_sort[g].G[1], GlobalGArray_sort[g].G[2], GlobalGArray_sort[g].GlobalIndex, GlobalGArray_sort[g].Index);
	499
	500
	501	// now index 'em through and build the hash table
	502	for (g=0;g<AllMaxG;g++) {// through every G of above process
	503	Lev->HashG[g].i = GlobalGArray_sort[g].Index;
	504	Lev->HashG[g].myPE = GlobalGArray_sort[g].GlobalIndex;
	505	//if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) HashG[%i]: i %i, myPE %i\n",P->Par.me, g, Lev->HashG[g].i, Lev->HashG[g].myPE);
	506	}
	507
	508	// find the dupes
	509	dupes = 0;
	510	for (g=0;g<AllMaxG-1;g++) // through every G of above process
	511	if ((GlobalGArray_sort[g].G[0] == GlobalGArray_sort[g+1].G[0]) &&
	512	(GlobalGArray_sort[g].G[1] == GlobalGArray_sort[g+1].G[1]) &&
	513	(GlobalGArray_sort[g].G[2] == GlobalGArray_sort[g+1].G[2])) {
	514	dupes++;
	515	if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) (i, myPE): (%i,%i) and (%i,%i) are duplicate!\n",P->Par.me, GlobalGArray_sort[g].Index, GlobalGArray_sort[g].GlobalIndex, GlobalGArray_sort[g+1].Index, GlobalGArray_sort[g+1].GlobalIndex);
	516	}
	517	if (P->Call.out[LeaderOut]) {
	518	if (dupes) fprintf(stderr,"(%i) %i duplicates in GArrays overall in Lev[%i]\n",P->Par.me, dupes, Lev->LevelNo);
	519	else fprintf(stderr,"(%i) There were no duplicate G vectors in Lev[%i] \n",P->Par.me, Lev->LevelNo);
	520	}
	521
	522	Free(GlobalGArray);
	523	Free(GlobalGArray_sort);
	524	for (i=0;i<NDIM;i++) {
	525	Free(G[i]);
	526	Free(G_send[i]);
	527	}
	528	Free(G);
	529	Free(G_send);
	530	}
	531
	532	/** Initialization of zeroth lattice level.
	533	* Creates each possible reciprocal grid vector, first counting them all (LatticeLevel::MaxG and LatticeLevel::MaxDoubleG),
	534	* then allocating memory in LatticeLevel::GArray, creating and storing them all. Afterwards, these are sorted.
	535	* LatticeLevel::MaxG is gathered from all processes and LatticeLevel::AllMaxG calculated
	536	\param P Problem at hand
	537	*/
	538	static void InitLatticeLevel0(struct Problem *const P) {
	539	struct Lattice *Lat = &P->Lat;
	540	struct LatticeLevel *Lev0 = Lat->Lev; // short for &(Lat->Lev[0])
	541	struct RPlans *RPs = &Lev0->Plan0;
	542	double GSq=0;
	543	double ECut = Lat->ECut;
	544	double G[NDIM];
	545	int x,ly,y,z,Index,d,Y,Z,iY,iZ,iy,iz,lz,X,i,yY,k,AllMaxG;
	546
	547	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitLatticeLevel0\n", P->Par.me);
	548	Lat->ECut = Lat->LevelSizes[0]Lat->LevelSizes[0];
	549	Lat->plan = fft_3d_create_plan(P, P->Par.comm_ST_Psi, Lat->ECut, Lat->MaxLevel, Lat->LevelSizes, Lat->ReciBasisOSQ, Lev0->N, Lat->MaxNoOfnFields, Lat->NFields);
	550	RPs->plan = &Lat->plan->Levplan[0];
	551	RPs->cdata = Lat->plan->local_data;
	552	RPs->rdata = (fftw_real *)Lat->plan->local_data;
	553	for (d=0; d < NDIM; d++) {
	554	Lev0->NUp[d] = 1;
	555	Lev0->NUp0[d] = 1;
	556	Lev0->NUp1[d] = 0;
	557	}
	558	Lev0->MaxN = Lev0->N[0]Lev0->N[1]Lev0->N[2];
	559	Lev0->PosFactorUp = NULL;
	560	Lev0->LevelNo = 0;
	561	Lev0->MaxG = 0;
	562	Lev0->MaxDoubleG = 0;
	563	Lev0->MaxNUp = 0;
	564	// goes through each possible reciprocal grid vector, counting them (for MaxG and MaxDoubleG)
	565	for (i = 0; i < Lat->plan->LocalMaxpw; i++) {
	566	Index = Lat->plan->Localpw[i].Index;
	567	Z = Index%(Lat->plan->NLSR[2]/2+1);
	568	Y = Index/(Lat->plan->NLSR[2]/2+1);
	569	for (ly = 0; ly < RPs->plan->LevelFactor; ly++) {
	570	yY = ly + RPs->plan->LevelFactor*Y;
	571	if (Z==Lat->plan->NLSR[2]/2) { // if it's the last possible one
	572	z = RPs->plan->N[2]/2;
	573	for (X=0; X < RPs->plan->N[0]; X++) {
	574	x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
	575	y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
	576	if (TestG(&P->Lat, x, y, z, &GSq, G)) Lev0->MaxG++;
	577	}
	578	} else {
	579	for (z=(Z)RPs->plan->LevelFactor; z < (Z+1)RPs->plan->LevelFactor; z++)
	580	for (X=0; X < RPs->plan->N[0]; X++) {
	581	x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
	582	y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
	583	if (TestG(&P->Lat, x, y, z, &GSq, G)) {
	584	switch (TestGDouble(x,y,z)) {
	585	case DoubleGNot:
	586	Lev0->MaxG++;
	587	break;
	588	case DoubleG0:
	589	Lev0->MaxG++;
	590	break;
	591	case DoubleGUse:
	592	Lev0->MaxG++;
	593	Lev0->MaxDoubleG++;
	594	break;
	595	case DoubleGNotUse:
	596	break;
	597	}
	598	}
	599	}
	600	}
	601	}
	602	}
	603	if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev0->MaxG %i MaxDoubleG %i\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], Lev0->MaxG, Lev0->MaxDoubleG);
	604	// then allocate memory for all these vectors
	605	Lev0->GArray = (struct OneGData *)
	606	Malloc(Lev0->MaxG*sizeof(struct OneGData),"InitLevel0: Lev0->GArray");
	607	Lat->GArrayForSort = (struct OneGData *)
	608	Malloc(Lev0->MaxG*sizeof(struct OneGData),"InitLevel0: Lat->GArrayForSort");
	609	if (Lev0->MaxDoubleG) {
	610	Lev0->DoubleG = (int *)
	611	Malloc(2sizeof(int)Lev0->MaxDoubleG,"InitLevel0: Lev0->DoubleG");
	612	} else {
	613	Lev0->DoubleG = NULL;
	614	}
	615	Lev0->MaxDoubleG = 0;
	616	Lev0->MaxG = 0;
	617	Lev0->ECut = 0;
	618	Lev0->G0 = -1;
	619	// and generate them again, this time noting them down in GArray
	620	for (i = 0; i < Lat->plan->LocalMaxpw; i++) {
	621	Index = Lat->plan->Localpw[i].Index;
	622	Z = Index%(Lat->plan->NLSR[2]/2+1);
	623	Y = Index/(Lat->plan->NLSR[2]/2+1);
	624	iZ = i%(Lat->plan->NLSR[2]/2+1);
	625	iY = i/(Lat->plan->NLSR[2]/2+1);
	626	for (ly = 0; ly < RPs->plan->LevelFactor; ly++) {
	627	yY = ly + RPs->plan->LevelFactor*Y;
	628	iy = ly + RPs->plan->LevelFactor*iY;
	629	if (Z==Lat->plan->NLSR[2]/2) {
	630	z = RPs->plan->N[2]/2;
	631	iz = RPs->plan->N[2]/2;
	632	for (X=0; X < RPs->plan->N[0]; X++) {
	633	x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
	634	y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
	635	if (TestG(&P->Lat, x, y, z, &GSq, G)) {
	636	Index = X+RPs->plan->N[0](iz+RPs->plan->nz(iy));
	637	Lev0->GArray[Lev0->MaxG].Index = Index;
	638	Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
	639	Lev0->GArray[Lev0->MaxG].GSq = GSq;
	640	Lev0->GArray[Lev0->MaxG].G[0] = G[0];
	641	Lev0->GArray[Lev0->MaxG].G[1] = G[1];
	642	Lev0->GArray[Lev0->MaxG].G[2] = G[2];
	643	if (GSq > Lev0->ECut) Lev0->ECut = GSq;
	644	Lev0->MaxG++;
	645	}
	646	}
	647	} else {
	648	for (lz=0; lz < RPs->plan->LevelFactor; lz++) {
	649	iz = iZ*RPs->plan->LevelFactor+lz;
	650	z = Z*RPs->plan->LevelFactor+lz;
	651	for (X=0; X < RPs->plan->N[0]; X++) {
	652	x = ( X >= RPs->plan->N[0]/2+1 ? X - RPs->plan->N[0] : X);
	653	y = ( yY >= RPs->plan->N[1]/2+1 ? yY - RPs->plan->N[1] : yY);
	654	if (TestG(&P->Lat, x, y, z, &GSq, G)) {
	655	switch (TestGDouble(x,y,z)) {
	656	case DoubleGNot:
	657	Index = X+RPs->plan->N[0](iz+RPs->plan->nz(iy));
	658	Lev0->GArray[Lev0->MaxG].Index = Index;
	659	Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
	660	Lev0->GArray[Lev0->MaxG].GSq = GSq;
	661	Lev0->GArray[Lev0->MaxG].G[0] = G[0];
	662	Lev0->GArray[Lev0->MaxG].G[1] = G[1];
	663	Lev0->GArray[Lev0->MaxG].G[2] = G[2];
	664	if (GSq > Lev0->ECut) Lev0->ECut = GSq;
	665	Lev0->MaxG++;
	666	break;
	667	case DoubleG0:
	668	Index = X+RPs->plan->N[0](iz+RPs->plan->nz(iy));
	669	Lev0->GArray[Lev0->MaxG].Index = Index;
	670	Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
	671	Lev0->GArray[Lev0->MaxG].GSq = GSq;
	672	Lev0->GArray[Lev0->MaxG].G[0] = G[0];
	673	Lev0->GArray[Lev0->MaxG].G[1] = G[1];
	674	Lev0->GArray[Lev0->MaxG].G[2] = G[2];
	675	if (GSq > Lev0->ECut) Lev0->ECut = GSq;
	676	Lev0->MaxG++;
	677	Lev0->G0 = Index;
	678	break;
	679	case DoubleGUse:
	680	Index = X+RPs->plan->N[0](iz+RPs->plan->nz(iy));
	681	Lev0->GArray[Lev0->MaxG].Index = Index;
	682	Lev0->GArray[Lev0->MaxG].GlobalIndex = 0;
	683	Lev0->GArray[Lev0->MaxG].GSq = GSq;
	684	Lev0->GArray[Lev0->MaxG].G[0] = G[0];
	685	Lev0->GArray[Lev0->MaxG].G[1] = G[1];
	686	Lev0->GArray[Lev0->MaxG].G[2] = G[2];
	687	if (GSq > Lev0->ECut) Lev0->ECut = GSq;
	688	Lev0->DoubleG[2*Lev0->MaxDoubleG] = Index;
	689	Lev0->DoubleG[2*Lev0->MaxDoubleG+1] = GetIndex(Lat->plan,RPs->plan,-x,-y,-z);
	690	if (Lev0->DoubleG[2Lev0->MaxDoubleG+1] < 0) Error(SomeError,"InitLevel0: Lev0->DoubleG[2Lev->MaxDoubleG+1]");
	691	Lev0->MaxG++;
	692	Lev0->MaxDoubleG++;
	693	break;
	694	case DoubleGNotUse:
	695	break;
	696	}
	697
	698	}
	699	}
	700	}
	701	}
	702	}
	703	}
	704	Lev0->ECut *= 0.5;
	705	naturalmergesort(Lev0->GArray,Lat->GArrayForSort,0,Lev0->MaxG-1,&GetKeyOneG,NULL,&CopyElementOneG);
	706	MPI_Allreduce(&Lev0->ECut , &ECut, 1, MPI_DOUBLE, MPI_MAX, P->Par.comm_ST_Psi );
	707	if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev0->ECut %g (%g) %g\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], Lev0->ECut, ECut, ECut - Lat->ECut);
	708	Lev0->ECut = Lat->ECut;
	709	if (ECut > Lat->ECut) Error(SomeError,"InitLevel0: ECut > Lat->ECut");
	710	Lev0->AllMaxG = (int *)
	711	Malloc(sizeof(int)*P->Par.Max_me_comm_ST_Psi,"InitLevel0: AllMaxG");
	712	MPI_Allgather ( &Lev0->MaxG, 1, MPI_INT,
	713	Lev0->AllMaxG, 1, MPI_INT, P->Par.comm_ST_Psi );
	714	AllMaxG = 0;
	715	for (k=0; k<P->Par.Max_me_comm_ST_Psi; k++) AllMaxG += Lev0->AllMaxG[k];
	716	if (P->Call.out[NormalOut]) {
	717	fprintf(stderr,"(%i) Lev0->AllMaxG %i -> RealAllMaxG %i\n", P->Par.mytid, AllMaxG, 2*AllMaxG-1);
	718	}
	719	Lev0->TotalAllMaxG = AllMaxG;
	720	Lev0->TotalRealAllMaxG = 2*AllMaxG-1;
	721	Lev0->HashG = (struct GDataHash *)
	722	Malloc(AllMaxG*sizeof(struct GDataHash),"InitLevel0: Lev0->HashG");
	723
	724	HashGArray(P,Lat,Lev0);
	725	}
	726
	727	/** Precalculate Position Factors LatticeLevel::PosFactorUp.
	728	* Due to the finer resolution of the densities certain transformations are necessary.
	729	* One position factor can be pulled out and pre-calculated
	730	* \f[
	731	* P_p(G) = \exp(i2\pi \sum^3_{j=1} \frac{g_j p_j}{2N_j} )
	732	* \f]
	733	* where \f$N_j\f$ is RPlans::plan::N, p is the position vector and G the reciprocal
	734	* grid vector, j goes up to NDIM (here 3)
	735	* \param *P Problem at hand
	736	* \param *Lat Lattice structure
	737	* \param *Lev LatticeLevel structure
	738	* \param *RPsUp RPlans structure
	739	* \param g number of reciprocal grid points
	740	* \param x number of mesh points in x direction (in the upper level)
	741	* \param y number of mesh points in y direction (in the upper level)
	742	* \param z number of mesh points in z direction (in the upper level)
	743	*/
	744	static void CreatePosFacForG(struct Problem P, struct Lattice Lat, const struct LatticeLevel
	745	Lev, const struct RPlans RPsUp, int g, int x, int y, int z)
	746	{
	747	double posfacreal[NDIM];
	748	struct RiemannTensor *RT = &Lat->RT;
	749	struct RunStruct *R = &P->R;
	750	int pos[NDIM];
	751	for (pos[0]=0; pos[0] < Lev->NUp[0]; pos[0]++) {
	752	for (pos[1]=0; pos[1] < Lev->NUp[1]; pos[1]++) {
	753	for (pos[2]=0; pos[2] < Lev->NUp[2]; pos[2]++) {
	754	posfacreal[0] = (x*pos[0])/(double)(RPsUp->plan->N[0]);
	755	posfacreal[1] = (y*pos[1])/(double)(RPsUp->plan->N[1]);
	756	posfacreal[2] = (z*pos[2])/(double)(RPsUp->plan->N[2]);
	757	Lev->PosFactorUp[pos[2]+Lev->NUp[2](pos[1]+Lev->NUp[1](pos[0]+Lev->NUp[0]g))].re = cos(2.PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
	758	Lev->PosFactorUp[pos[2]+Lev->NUp[2](pos[1]+Lev->NUp[1](pos[0]+Lev->NUp[0]g))].im = sin(2.PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
	759	}
	760	}
	761	}
	762	if (RT->Use == UseRT) {
	763	if (Lev->LevelNo == R->LevRNo) {
	764	for (pos[0]=0; pos[0] < RT->NUpLevRS[0]; pos[0]++) {
	765	for (pos[1]=0; pos[1] < RT->NUpLevRS[1]; pos[1]++) {
	766	for (pos[2]=0; pos[2] < RT->NUpLevRS[2]; pos[2]++) {
	767	posfacreal[0] = (x*pos[0])/(double)(Lat->Lev[R->LevRSNo].N[0]);
	768	posfacreal[1] = (y*pos[1])/(double)(Lat->Lev[R->LevRSNo].N[1]);
	769	posfacreal[2] = (z*pos[2])/(double)(Lat->Lev[R->LevRSNo].N[2]);
	770	RT->PosFactor[RTPFRtoS][pos[2]+RT->NUpLevRS[2](pos[1]+RT->NUpLevRS[1](pos[0]+RT->NUpLevRS[0]g))].re = cos(2.PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
	771	RT->PosFactor[RTPFRtoS][pos[2]+RT->NUpLevRS[2](pos[1]+RT->NUpLevRS[1](pos[0]+RT->NUpLevRS[0]g))].im = sin(2.PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
	772	}
	773	}
	774	}
	775	for (pos[0]=0; pos[0] < RT->NUpLevR0[0]; pos[0]++) {
	776	for (pos[1]=0; pos[1] < RT->NUpLevR0[1]; pos[1]++) {
	777	for (pos[2]=0; pos[2] < RT->NUpLevR0[2]; pos[2]++) {
	778	posfacreal[0] = (x*pos[0])/(double)(Lat->Lev[R->LevR0No].N[0]);
	779	posfacreal[1] = (y*pos[1])/(double)(Lat->Lev[R->LevR0No].N[1]);
	780	posfacreal[2] = (z*pos[2])/(double)(Lat->Lev[R->LevR0No].N[2]);
	781	RT->PosFactor[RTPFRto0][pos[2]+RT->NUpLevR0[2](pos[1]+RT->NUpLevR0[1](pos[0]+RT->NUpLevR0[0]g))].re = cos(2.PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
	782	RT->PosFactor[RTPFRto0][pos[2]+RT->NUpLevR0[2](pos[1]+RT->NUpLevR0[1](pos[0]+RT->NUpLevR0[0]g))].im = sin(2.PI*(posfacreal[0]+posfacreal[1]+posfacreal[2]));
	783	}
	784	}
	785	}
	786	}
	787	}
	788	}
	789
	790	/** Initialization of first to Lattice::MaxLevel LatticeLevel's.
	791	* Sets numbers and ratios of mesh points (maximum (LatticeLevel::MaxN, LatticeLevel::MaxNUp), per dimension
	792	* (LatticeLevel::N[], LatticeLevel::NUp[])) for the respective levels and its upper level. LatticeLevel::LevelNo and
	793	* RPlans are set from LevelPlan[], number of reciprocal grid vectors is recounted (LatticeLevel::MaxG,
	794	* LatticeLevel::MaxDoubleG), then they are created, stored (LatticeLevel::GArray) and their position factor
	795	* (LatticeLevel::PosFactorUp) recalculated. MaxG is gathered from all processes in LatticeLevel::AllMaxG and
	796	* LatticeLevel::TotalAllMaxG generated by summing these up.
	797	* \param *P Problem at hand
	798	* \sa InitLatticeLevel0()
	799	*/
	800	static void InitOtherLevel(struct Problem *const P) {
	801	struct Lattice *Lat = &P->Lat;
	802	struct LatticeLevel *Lev;
	803	struct LatticeLevel *LevUp;
	804	struct RPlans *RPs;
	805	struct RPlans *RPsUp;
	806	double GSq=0,ECut=0;
	807	int lx,x,ly,y,z,Index,Y,YY,Z,ZZ,lz,I,k,AllMaxG;
	808	int d,MaxGUp,i;
	809	for (i=1; i < P->Lat.MaxLevel; i++) { // for each level
	810	Lev = &Lat->Lev[i];
	811	LevUp = &Lat->Lev[i-1];
	812	RPs = &Lev->Plan0;
	813	RPsUp = &LevUp->Plan0;
	814	// set number of mesh points
	815	for (d=0; d < NDIM; d++) { // for each dimension
	816	Lev->NUp[d] = Lat->LevelSizes[i-1];
	817	Lev->NUp0[d] = LevUp->NUp0[d]*Lev->NUp[d];
	818	Lev->NUp1[d] = (i == STANDARTLEVEL ? 1 : LevUp->NUp1[d]*Lev->NUp[d]);
	819	Lev->N[d] = LevUp->N[d]/Lev->NUp[d];
	820	}
	821	Lev->MaxN = Lev->N[0]Lev->N[1]Lev->N[2];
	822	Lev->MaxNUp = Lev->NUp[0]Lev->NUp[1]Lev->NUp[2];
	823	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitLatticeLevel %i\tN[] = %i %i %i\n",P->Par.me, i,Lev->N[0],Lev->N[1],Lev->N[2]);
	824	// fft plan and its data
	825	RPs->plan = &Lat->plan->Levplan[i];
	826	RPs->cdata = Lat->plan->local_data;
	827	RPs->rdata = (fftw_real *)Lat->plan->local_data;
	828	Lev->LevelNo = i;
	829	// recount MaxG, MaxDoubleG by creating all g vectors
	830	Lev->MaxG=0;
	831	Lev->MaxDoubleG = 0;
	832	for (MaxGUp=0; MaxGUp < LevUp->MaxG; MaxGUp++) {
	833	SplitIndex(LevUp->GArray[MaxGUp].Index,RPsUp->plan->N[0],RPsUp->plan->nz,&lx,&ly,&lz);
	834	Y = ly/RPsUp->plan->LevelFactor;
	835	YY = ly%RPsUp->plan->LevelFactor;
	836	Z = lz/RPsUp->plan->LevelFactor;
	837	ZZ = lz%RPsUp->plan->LevelFactor;
	838	Index = Lat->plan->Localpw[Z+(Lat->plan->NLSR[2]/2+1)*(Y)].Index;
	839	ly = YY+RPsUp->plan->LevelFactor*(Index/(Lat->plan->NLSR[2]/2+1));
	840	z = ZZ+RPsUp->plan->LevelFactor*(Index%(Lat->plan->NLSR[2]/2+1));
	841	x = ( lx >= RPsUp->plan->N[0]/2+1 ? lx-RPsUp->plan->N[0] : lx);
	842	y = ( ly >= RPsUp->plan->N[1]/2+1 ? ly-RPsUp->plan->N[1] : ly);
	843	if ((x%Lev->NUp[0] == 0) && (y%Lev->NUp[1] == 0) && (z%Lev->NUp[2] == 0)) {
	844	switch (TestGDouble(x/Lev->NUp[0],y/Lev->NUp[1],z/Lev->NUp[2])) {
	845	case DoubleGNot:
	846	Lev->MaxG++;
	847	break;
	848	case DoubleG0:
	849	Lev->MaxG++;
	850	break;
	851	case DoubleGUse:
	852	Lev->MaxG++;
	853	Lev->MaxDoubleG++;
	854	break;
	855	case DoubleGNotUse:
	856	Lev->MaxG--;
	857	break;
	858	}
	859	}
	860	}
	861	if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev%i->MaxG %i MaxDoubleG %i\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], i, Lev->MaxG, Lev->MaxDoubleG);
	862	Lev->PosFactorUp = (fftw_complex ) Malloc(sizeof(fftw_complex)Lev->MaxGLev->NUp[0]Lev->NUp[1]*Lev->NUp[2], "InitLatticeLevel: PosFactors");
	863	if (Lat->RT.Use == UseRT) {
	864	if (Lev->LevelNo == Lat->RT.RiemannLevel) {
	865	Lat->RT.PosFactor[RTPFRtoS] = (fftw_complex *)
	866	Malloc(sizeof(fftw_complex)Lev->MaxGLev->NUp1[0]Lev->NUp1[1]Lev->NUp1[2], "InitLatticeLevel: RTPosFactors1");
	867	Lat->RT.PosFactor[RTPFRto0] = (fftw_complex *)
	868	Malloc(sizeof(fftw_complex)Lev->MaxGLev->NUp0[0]Lev->NUp0[1]Lev->NUp0[2], "InitLatticeLevel: RTPosFactors0");
	869	}
	870	}
	871	// now create the G vectors and store them
	872	Lev->GArray = (struct OneGData *)
	873	Malloc(Lev->MaxG*sizeof(struct OneGData),"InitLevel: Lev->GArray");
	874	if (Lev->MaxDoubleG) {
	875	Lev->DoubleG = (int *)
	876	Malloc(2sizeof(int)Lev->MaxDoubleG,"InitLevel: Lev0->DoubleG");
	877	} else {
	878	Lev->DoubleG = NULL;
	879	}
	880	Lev->MaxDoubleG = 0;
	881	Lev->MaxG=0;
	882	Lev->ECut = 0;
	883	Lev->G0 = -1;
	884	for (MaxGUp=0; MaxGUp < LevUp->MaxG; MaxGUp++) {
	885	SplitIndex(LevUp->GArray[MaxGUp].Index,RPsUp->plan->N[0],RPsUp->plan->nz,&lx,&ly,&lz);
	886	Index = lx/Lev->NUp[0]+RPs->plan->N[0](lz/Lev->NUp[2]+RPs->plan->nz(ly/Lev->NUp[1]));
	887	Y = ly/RPsUp->plan->LevelFactor;
	888	YY = ly%RPsUp->plan->LevelFactor;
	889	Z = lz/RPsUp->plan->LevelFactor;
	890	ZZ = lz%RPsUp->plan->LevelFactor;
	891	I = Lat->plan->Localpw[Z+(Lat->plan->NLSR[2]/2+1)*(Y)].Index;
	892	ly = YY+RPsUp->plan->LevelFactor*(I/(Lat->plan->NLSR[2]/2+1));
	893	z = ZZ+RPsUp->plan->LevelFactor*(I%(Lat->plan->NLSR[2]/2+1));
	894	x = ( lx >= RPsUp->plan->N[0]/2+1 ? lx-RPsUp->plan->N[0] : lx);
	895	y = ( ly >= RPsUp->plan->N[1]/2+1 ? ly-RPsUp->plan->N[1] : ly);
	896	if ((x%Lev->NUp[0] == 0) && (y%Lev->NUp[1] == 0) && (z%Lev->NUp[2] == 0)) {
	897	switch (TestGDouble(x/Lev->NUp[0],y/Lev->NUp[1],z/Lev->NUp[2])) {
	898	case DoubleGNot:
	899	TestG(&P->Lat, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2], &GSq, Lev->GArray[Lev->MaxG].G);
	900	Lev->GArray[Lev->MaxG].Index = Index;
	901	Lev->GArray[Lev->MaxG].GlobalIndex = 0;
	902	Lev->GArray[Lev->MaxG].GSq = GSq;
	903	if (GSq > Lev->ECut) Lev->ECut = GSq;
	904	CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
	905	Lev->MaxG++;
	906	break;
	907	case DoubleG0:
	908	TestG(&P->Lat, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2], &GSq, Lev->GArray[Lev->MaxG].G);
	909	Lev->GArray[Lev->MaxG].Index = Index;
	910	Lev->GArray[Lev->MaxG].GlobalIndex = 0;
	911	Lev->GArray[Lev->MaxG].GSq = GSq;
	912	if (GSq > Lev->ECut) Lev->ECut = GSq;
	913	CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
	914	Lev->MaxG++;
	915	Lev->G0 = Index;
	916	break;
	917	case DoubleGUse:
	918	TestG(&P->Lat, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2], &GSq, Lev->GArray[Lev->MaxG].G);
	919	Lev->GArray[Lev->MaxG].Index = Index;
	920	Lev->GArray[Lev->MaxG].GlobalIndex = 0;
	921	Lev->GArray[Lev->MaxG].GSq = GSq;
	922	if (GSq > Lev->ECut) Lev->ECut = GSq;
	923	CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
	924	Lev->DoubleG[2*Lev->MaxDoubleG] = Index;
	925	Lev->DoubleG[2*Lev->MaxDoubleG+1] = GetIndex(Lat->plan,RPs->plan,-(x/Lev->NUp[0]),-(y/Lev->NUp[1]),-(z/Lev->NUp[2]));
	926	if (Lev->DoubleG[2Lev->MaxDoubleG+1] < 0) Error(SomeError,"Lev->DoubleG[2Lev->MaxDoubleG+1]");
	927	CreatePosFacForG(P, Lat, Lev, RPsUp, Lev->MaxG, x/Lev->NUp[0], y/Lev->NUp[1], z/Lev->NUp[2]);
	928	Lev->MaxG++;
	929	Lev->MaxDoubleG++;
	930	break;
	931	case DoubleGNotUse:
	932	break;
	933	}
	934	}
	935	}
	936	Lev->ECut *= 0.5;
	937	MPI_Allreduce(&Lev->ECut , &ECut, 1, MPI_DOUBLE, MPI_MAX, P->Par.comm_ST_Psi );
	938	if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) (%i,%i) Lev%i->ECut %g (%g) %g\n", P->Par.mytid, P->Par.mypos[PEPsi], P->Par.mypos[PEGamma], i, Lev->ECut, ECut, ECut - Lat->ECut/(Lev->NUp0[0]*Lev->NUp0[0]));
	939	Lev->ECut = Lat->ECut/(Lev->NUp0[0]*Lev->NUp0[0]);
	940	if (ECut > Lev->ECut) Error(SomeError,"InitLevel: ECut > Lat->ECut");
	941	Lev->AllMaxG = (int *)
	942	Malloc(sizeof(int)*P->Par.Max_me_comm_ST_Psi,"InitLevel: AllMaxG");
	943	MPI_Allgather ( &Lev->MaxG, 1, MPI_INT,
	944	Lev->AllMaxG, 1, MPI_INT, P->Par.comm_ST_Psi );
	945	AllMaxG = 0;
	946	for (k=0; k<P->Par.Max_me_comm_ST_Psi; k++) AllMaxG += Lev->AllMaxG[k];
	947	if (P->Call.out[NormalOut]) {
	948	fprintf(stderr,"(%i) Lev%i->AllMaxG %i -> RealAllMaxG %i\n", P->Par.mytid, i, AllMaxG, 2*AllMaxG-1);
	949	}
	950	Lev->TotalAllMaxG = AllMaxG;
	951	Lev->TotalRealAllMaxG = 2*AllMaxG-1;
	952	Lev->HashG = (struct GDataHash *)
	953	Malloc(AllMaxG*sizeof(struct GDataHash),"InitLevel0: Lev0->HashG");
	954	HashGArray(P, Lat, Lev);
	955	}
	956	}
	957
	958	/** Sets the Spin group number for each process.
	959	* First of all, SpinUps and SpinDowns Psis are completely separated (they have their own communicator groups).
	960	* Psis::PsiGroup is then just the rank in this Psi group, where Psis::MaxPsiGroup is the number of participating
	961	* processes. Afterwards, depending on the SpinType Psis:PsiST is set and Psis::LocalNo calculated by modulo of
	962	* global number of Psis and the number of Psi groups, e.g. 18 psis, 5 processes, make 3 times 4 Psis per process
	963	* and 2 times 3 Psis per process. Finally, the maximum of all the local number of Psis is taken and stored in
	964	* Psis::AllMaxLocalNo.
	965	* \param *P Problem at hand
	966	*/
	967	static void InitPsiGroupNo(struct Problem *const P) {
	968	struct Lattice *Lat = &P->Lat;
	969	struct Psis *Psi = &Lat->Psi;
	970	int i, r, r1, NoPsis = 0, type;
	971	MPI_Status status;
	972	int AllMaxLocalNoUp, AllMaxLocalNoDown;
	973	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitPsisGroupNo\n", P->Par.me);
	974	Psi->MaxPsiGroup = P->Par.Max_my_color_comm_ST_Psi;
	975	Psi->PsiGroup = P->Par.my_color_comm_ST_Psi;
	976	// now discern between the different spin cases
	977	switch (Psi->Use) {
	978	case UseSpinDouble:
	979	NoPsis = Psi->GlobalNo[PsiMaxNoDouble];
	980	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i) PsiGlobalNo %i = %i + %i\n", P->Par.me, Psi->GlobalNo[PsiMaxNo], NoPsis, Psi->GlobalNo[PsiMaxAdd]);
	981	Psi->PsiST = SpinDouble; // set SpinType
	982	break;
	983	case UseSpinUpDown:
	984	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i) PsiGlobalNo %i = (%i+%i) + (%i+%i)\n", P->Par.me, Psi->GlobalNo[PsiMaxNo], Psi->GlobalNo[PsiMaxNoUp], Psi->GlobalNo[PsiMaxAdd], Psi->GlobalNo[PsiMaxNoDown], Psi->GlobalNo[PsiMaxAdd]);
	985	// depending on my spin type (up or down) ...
	986	switch (P->Par.my_color_comm_ST) {
	987	case ((int)SpinUp) - 1:
	988	Psi->PsiST = SpinUp;
	989	NoPsis = Psi->GlobalNo[PsiMaxNoUp];
	990	break;
	991	case ((int)SpinDown) -1:
	992	Psi->PsiST = SpinDown;
	993	NoPsis = Psi->GlobalNo[PsiMaxNoDown];
	994	break;
	995	}
	996	break;
	997	}
	998	Psi->LocalNo = (NoPsis + P->Lat.Psi.GlobalNo[PsiMaxAdd]) / Psi->MaxPsiGroup; // local number of Psis per process
	999	r = (NoPsis + P->Lat.Psi.GlobalNo[PsiMaxAdd]) % Psi->MaxPsiGroup; // calculate the rest
	1000	Psi->LocalNoAdd = (P->Lat.Psi.GlobalNo[PsiMaxAdd]) / Psi->MaxPsiGroup; // local number of Psis per process
	1001	r1 = (P->Lat.Psi.GlobalNo[PsiMaxAdd]) % Psi->MaxPsiGroup; // calculate the rest
	1002	if (P->R.DoUnOccupied && (r1 != 0)) Error(SomeError,"AddPsis must be multiple of PEPsi!\n");
	1003	if (Psi->PsiGroup < r) Psi->LocalNo++; // distribute the rest into the first processes
	1004	if (Psi->PsiGroup < r1) Psi->LocalNoAdd++; // distribute the rest into the first processes
	1005	Psi->TypeStartIndex[Occupied] = 0;
	1006	Psi->TypeStartIndex[UnOccupied] = Psi->LocalNo - Psi->LocalNoAdd;
	1007	if (P->R.DoPerturbation == 1) { // six times NoPsis for Perturbed_RxP and Perturbed_P
	1008	for (i=Perturbed_P0; i<= Perturbed_RxP2;i++) {
	1009	Psi->TypeStartIndex[i] = Psi->LocalNo;
	1010	Psi->LocalNo += (NoPsis) / Psi->MaxPsiGroup; // local number of Psis per process (just one array for all perturbed)
	1011	r = (NoPsis) % Psi->MaxPsiGroup; // calculate the rest
	1012	if (Psi->PsiGroup < r) Error(SomeError,"Psis must be multiple of PEPsi!\n");
	1013	//if (Psi->PsiGroup < r) Psi->LocalNo++; // distribute the rest into the first processes
	1014	}
	1015	} else {
	1016	for (i=Perturbed_P0; i<= Perturbed_RxP2;i++)
	1017	Psi->TypeStartIndex[i] = Psi->LocalNo;
	1018	}
	1019	Psi->TypeStartIndex[Extra] = Psi->LocalNo; // always set to end
	1020	Psi->TypeStartIndex[Extra+1] = Psi->LocalNo+1; // always set to end
	1021
	1022	if (P->Call.out[PsiOut])
	1023	for(type=Occupied;type<=Extra+1;type++)
	1024	fprintf(stderr,"(%i) TypeStartIndex[%i] = %i\n",P->Par.me, type,Psi->TypeStartIndex[type]);
	1025
	1026	if (P->Call.out[StepLeaderOut]) fprintf(stderr, "(%i) PsiLocalNo %i \n", P->Par.me, Psi->LocalNo);
	1027	// gather local number of Psis from all processes
	1028	Psi->RealAllLocalNo = (int *)
	1029	Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"FirstInitGramSchData: Psi->RealAllLocalNo");
	1030	MPI_Allgather ( &Psi->LocalNo, 1, MPI_INT,
	1031	Psi->RealAllLocalNo, 1, MPI_INT, P->Par.comm_ST_PsiT );
	1032	// determine maximum local number, depending on spin
	1033	switch (Psi->PsiST) {
	1034	case SpinDouble:
	1035	MPI_Allreduce(&Psi->LocalNo, &Psi->AllMaxLocalNo, 1,
	1036	MPI_INT, MPI_MAX, P->Par.comm_ST_PsiT );
	1037	break;
	1038	case SpinUp:
	1039	MPI_Allreduce(&Psi->LocalNo, &AllMaxLocalNoUp, 1,
	1040	MPI_INT, MPI_MAX, P->Par.comm_ST_PsiT );
	1041	MPI_Sendrecv( &AllMaxLocalNoUp, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag1,
	1042	&AllMaxLocalNoDown, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag2,
	1043	P->Par.comm_STInter, &status );
	1044	Psi->AllMaxLocalNo = MAX(AllMaxLocalNoUp, AllMaxLocalNoDown);
	1045	break;
	1046	case SpinDown:
	1047	MPI_Allreduce(&Psi->LocalNo, &AllMaxLocalNoDown, 1,
	1048	MPI_INT, MPI_MAX, P->Par.comm_ST_PsiT );
	1049	MPI_Sendrecv( &AllMaxLocalNoDown, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag2,
	1050	&AllMaxLocalNoUp, 1, MPI_INT, P->Par.me_comm_ST, AllMaxLocalNoTag1,
	1051	P->Par.comm_STInter, &status );
	1052	Psi->AllMaxLocalNo = MAX(AllMaxLocalNoUp, AllMaxLocalNoDown);
	1053	break;
	1054	}
	1055	}
	1056
	1057	/** Intialization of the wave functions.
	1058	* Allocates/sets OnePsiElementAddData elements to zero, allocates memory for coefficients and
	1059	* pointing to them from Lattice and LatticeLevel.
	1060	* \param *P Problem at hand
	1061	*/
	1062	static void InitPsis(struct Problem *const P) {
	1063	int i,j,l;
	1064	struct Lattice *Lat = &P->Lat;
	1065	struct RunStruct *R = &P->R;
	1066	struct LatticeLevel *Lev = R->InitLevS;
	1067	int CreateOld = P->R.UseOldPsi;
	1068	int NoPsis;
	1069	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitPsis\n", P->Par.me);
	1070	/// For Lattice::Psi memory for OnePsiElementAddData is allocated and entries zeroed.
	1071	Lat->Psi.AddData = (struct OnePsiElementAddData *)
	1072	Malloc((Lat->Psi.LocalNo+1)*sizeof(struct OnePsiElementAddData), "InitPsis: Psi->AddData");
	1073	for (i = 0; i <= Lat->Psi.LocalNo; i++) {
	1074	Lat->Psi.AddData[i].T = 0.0;
	1075	Lat->Psi.AddData[i].Lambda = 0.0;
	1076	Lat->Psi.AddData[i].Gamma = 0.0;
	1077	Lat->Psi.AddData[i].Step = 0;
	1078	Lat->Psi.AddData[i].WannierSpread = 0.;
	1079	for (j=0;j<NDIM;j++)
	1080	Lat->Psi.AddData[i].WannierCentre[j] = 0. ; //Lat->RealBasisQ[j]/2.;
	1081	}
	1082	// lambda contains H eigenvalues
	1083	switch (Lat->Psi.PsiST) {
	1084	default:
	1085	case SpinDouble:
	1086	Lat->Psi.NoOfPsis = Lat->Psi.GlobalNo[PsiMaxNoDouble];
	1087	Lat->Psi.NoOfTotalPsis = Lat->Psi.NoOfPsis + Lat->Psi.GlobalNo[PsiMaxAdd];
	1088	break;
	1089	case SpinUp:
	1090	Lat->Psi.NoOfPsis = Lat->Psi.GlobalNo[PsiMaxNoUp];
	1091	Lat->Psi.NoOfTotalPsis = Lat->Psi.NoOfPsis + Lat->Psi.GlobalNo[PsiMaxAdd];
	1092	break;
	1093	case SpinDown:
	1094	Lat->Psi.NoOfPsis = Lat->Psi.GlobalNo[PsiMaxNoDown];
	1095	Lat->Psi.NoOfTotalPsis = Lat->Psi.NoOfPsis + Lat->Psi.GlobalNo[PsiMaxAdd];
	1096	break;
	1097	};
	1098	if (P->Call.out[StepLeaderOut]) fprintf(stderr, "(%i) NoOfPsis %i \n", P->Par.me, Lat->Psi.NoOfPsis);
	1099	Lat->Psi.lambda = (double *) Malloc(sizeof(double )* Lat->Psi.NoOfTotalPsis, "InitPsis: Psi->Lambda");
	1100	Lat->Psi.Overlap = (double *) Malloc(sizeof(double )* Lat->Psi.NoOfPsis, "InitPsis: Psi->Overlap");
	1101	for (i=0;i<Lat->Psi.NoOfTotalPsis;i++)
	1102	Lat->Psi.lambda[i] = (double ) Malloc(sizeof(double) Lat->Psi.NoOfTotalPsis, "InitPsis: Psi->Lambda[i]");
	1103	for (i=0;i<Lat->Psi.NoOfPsis;i++)
	1104	Lat->Psi.Overlap[i] = (double ) Malloc(sizeof(double) Lat->Psi.NoOfPsis, "InitPsis: Psi->Overlap[i]");
	1105
	1106	/// Allocation for initial LevelPsi LatticeLevel::LPsi and LevelPsi::PsiDat
	1107	NoPsis = (Lat->Psi.TypeStartIndex[Perturbed_P1] - Lat->Psi.TypeStartIndex[Occupied]); //
	1108	if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) InitPsis(): Number of Psi arrays in memory %i\n", P->Par.me, NoPsis);
	1109	Lev->LPsi = (struct LevelPsi *)
	1110	Malloc(sizeof(struct LevelPsi),"InitPsis: LPsi");
	1111	Lev->LPsi->PsiDat = (fftw_complex *)
	1112	Malloc((NoPsis+3)sizeof(fftw_complex)Lev->MaxG, "InitPsis: LPsi->PsiDat"); // +3 because of TempPsi and TempPsi2 and extra
	1113	if (CreateOld)
	1114	Lev->LPsi->OldPsiDat = (fftw_complex *)
	1115	Malloc((NoPsis+1)sizeof(fftw_complex)Lev->MaxG, "InitPsis: LPsi->PsiDat");
	1116	else
	1117	Lev->LPsi->OldPsiDat = NULL;
	1118
	1119	/// Allocation for all other levels LatticeLevel::LPsi and LevelPsi::PsiDat, the latter initialised from initial level
	1120	/// here Lev[1].LPsi, because it's the biggest level for the Psis and they might get overwritten otherwise during the perturbed
	1121	/// load and transform-sequence, see ReadSrcPerturbedPsis()
	1122	for (i=1; i < Lat->MaxLevel; i++) { // for all levels
	1123	if (i > 1) Lat->Lev[i].LPsi = (struct LevelPsi *)
	1124	Malloc(sizeof(struct LevelPsi),"InitPsis: LPsi");
	1125	Lat->Lev[i].LPsi->PsiDat = Lev->LPsi->PsiDat;
	1126	Lat->Lev[i].LPsi->OldPsiDat = Lev->LPsi->OldPsiDat;
	1127	Lat->Lev[i].LPsi->LocalPsi = (fftw_complex **)
	1128	Malloc((Lat->Psi.LocalNo+1)sizeof(fftw_complex ),"InitPsi: LocalPsi");
	1129	/// Initialise each LevelPsi::LocalPsi by pointing it to its respective part of LevelPsi::PsiDat
	1130	for (j=0; j <= Lat->Psi.LocalNo; j++) { // for each local Psi
	1131	if (j < Lat->Psi.TypeStartIndex[Perturbed_P0]) { // if it's not a perturbed wave function
	1132	//if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Setting LocalPsi[%i] to %p\n",P->Par.me, j, &Lat->Lev[i].LPsi->PsiDat[j*Lat->Lev[1].MaxG]);
	1133	Lat->Lev[i].LPsi->LocalPsi[j] = &Lat->Lev[i].LPsi->PsiDat[j*Lat->Lev[1].MaxG];
	1134	} else if (j >= Lat->Psi.TypeStartIndex[Extra]) { // if it's the extra wave function
	1135	//if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Setting LocalPsi[%i] to %p\n",P->Par.me, j, &Lat->Lev[i].LPsi->PsiDat[(NoPsis)*Lat->Lev[1].MaxG]);
	1136	Lat->Lev[i].LPsi->LocalPsi[j] = &Lat->Lev[i].LPsi->PsiDat[(NoPsis)*Lat->Lev[1].MaxG];
	1137	} else { // ... a perturbed wave function
	1138	l = Lat->Psi.TypeStartIndex[Perturbed_P0]
	1139	+ ((j - Lat->Psi.TypeStartIndex[Perturbed_P0]) % (Lat->Psi.TypeStartIndex[Perturbed_P1] - Lat->Psi.TypeStartIndex[Perturbed_P0]));
	1140	//if (P->Call.out[PsiOut]) fprintf(stderr,"(%i) Setting LocalPsi[%i] to %p\n",P->Par.me, j, &Lat->Lev[i].LPsi->PsiDat[l*Lat->Lev[1].MaxG]);
	1141	Lat->Lev[i].LPsi->LocalPsi[j] = &Lat->Lev[i].LPsi->PsiDat[l*Lat->Lev[1].MaxG];
	1142	}
	1143	}
	1144	Lat->Lev[i].LPsi->TempPsi = &Lat->Lev[i].LPsi->PsiDat[(NoPsis+1)*Lat->Lev[1].MaxG];
	1145	Lat->Lev[i].LPsi->TempPsi2 = &Lat->Lev[i].LPsi->PsiDat[(NoPsis+2)*Lat->Lev[1].MaxG];
	1146
	1147	Lat->Lev[i].LPsi->OldLocalPsi = (fftw_complex **)
	1148	Malloc(((Lat->Psi.TypeStartIndex[Perturbed_P0] - Lat->Psi.TypeStartIndex[Occupied])+1)sizeof(fftw_complex ),"InitPsi: OldLocalPsi");
	1149	for (j=Lat->Psi.TypeStartIndex[Occupied]; j <= Lat->Psi.TypeStartIndex[UnOccupied]; j++) { // OldPsis are only needed during the Wannier procedure
	1150	if (CreateOld)
	1151	Lat->Lev[i].LPsi->OldLocalPsi[j] = &Lat->Lev[i].LPsi->OldPsiDat[j*Lat->Lev[1].MaxG];
	1152	else
	1153	Lat->Lev[i].LPsi->OldLocalPsi[j] = NULL;
	1154	}
	1155	}
	1156	Lat->Lev[0].LPsi = NULL;
	1157	}
	1158
	1159	/** Initialization of wave function coefficients.
	1160	* The random number generator's seed is set to 1, rand called numerously depending on the rank in the
	1161	* communicator, the number of LatticeLevel::AllMaxG in each process. Finally, for each local Psi and each
	1162	* reciprocal grid vector the real and imaginary part of the specific coefficient \f$c_{i,G}\f$ are
	1163	* initialized by a \f$1-\frac{1}{2} \frac{rand()}{\|G\|^2}\f$, where rand() yields a number between 0
	1164	* and the largest integer value.
	1165	* \param *P Problem at hand
	1166	* \param start The init starts at local this local wave function array ...
	1167	* \param end ... and ends at this local array, e.g. Lattice#Psis#LocalNo
	1168	*/
	1169	void InitPsisValue(struct Problem *const P, int start, int end) {
	1170	int i,j,k, index;
	1171	int myPE;
	1172	int NoBefore = 0;
	1173	struct Lattice *Lat = &P->Lat;
	1174	struct RunStruct *R = &P->R;
	1175	struct LatticeLevel *LevS = R->LevS;
	1176	MPI_Comm_rank(P->Par.comm_ST_Psi, &myPE); // Determines the rank of the calling process in the communicator
	1177	//char name[255];
	1178	//sprintf(name, ".Psis-%i.all", myPE);
	1179	//FILE *psi;
	1180	//OpenFile(P,&psi, name, "w",P->Call.out[ReadOut]);
	1181
	1182	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitPsisValue\n", P->Par.me);
	1183	for (i=0; i < P->Par.my_color_comm_ST_Psi; i++) {
	1184	NoBefore += Lat->Psi.RealAllLocalNo[i];
	1185	}
	1186	srand(R->Seed);
	1187	for (i=0; i<NoBefore; i++)
	1188	for (k=0; k<P->Par.Max_me_comm_ST_Psi; k++)
	1189	for (j=0; j<2*LevS->AllMaxG[k]; j++)
	1190	rand();
	1191	for (j=0;(j<end) && (j < Lat->Psi.LocalNo);j++) // for all local Psis
	1192	for (i=0; i < LevS->TotalAllMaxG; i++) { // for all coefficients
	1193	if ((LevS->HashG[i].myPE == myPE) && (j >= start)) {
	1194	index = LevS->HashG[i].i;
	1195	//if (!j) fprintf(stderr,"(%i) LevS->GArray[%i].GSq = %e\n",P->Par.me, index, LevS->GArray[index].GSq);
	1196	LevS->LPsi->LocalPsi[j][index].re = (1.-2.(rand()/(double)RAND_MAX))(LevS->GArray[index].GSq == 0.0 ? 1.0 : 1/LevS->GArray[index].GSq);
	1197	LevS->LPsi->LocalPsi[j][index].im = (1.-2.(rand()/(double)RAND_MAX))(LevS->GArray[index].GSq == 0.0 ? 0.0 : 1/LevS->GArray[index].GSq);
	1198	//if (!j) fprintf(stderr,"(%i) LevS->LPsi->LocalPsi[%i][%i/%i] = %e + i%e, GArray[%i].GSq = %e\n",myPE, j, index, i, LevS->LPsi->LocalPsi[j][index].re, LevS->LPsi->LocalPsi[j][index].im, index, LevS->GArray[index].GSq);
	1199	//fprintf(psi, "LevS->LPsi->LocalPsi[%i][%i] = %e + i%e\n", j, i, LevS->LPsi->LocalPsi[j][index].re, LevS->LPsi->LocalPsi[j][index].im);
	1200	} else {
	1201	rand();
	1202	rand();
	1203	}
	1204
	1205	/* for (j=0;j<Lat->Psi.LocalNo;j++) // for all local Psis
	1206	for (i=0; i < LevS->MaxG; i++) { // for all coefficients
	1207	if ((R->Seed == 0) && (j >= Lat->Psi.TypeStartIndex[Perturbed_P0])) {
	1208	if (i == 0) fprintf(stderr,"(%i) Initializing wave function %i with 1 over %i\n",P->Par.me, j, LevS->TotalAllMaxG);
	1209	LevS->LPsi->LocalPsi[j][i].re = 1./LevS->TotalAllMaxG;
	1210	LevS->LPsi->LocalPsi[j][i].im = 0.;
	1211	} else {
	1212	LevS->LPsi->LocalPsi[j][i].re = (1.-2.(rand()/(double)RAND_MAX))(LevS->GArray[i].GSq == 0.0 ? 1.0 : 1/LevS->GArray[i].GSq);
	1213	LevS->LPsi->LocalPsi[j][i].im = (1.-2.(rand()/(double)RAND_MAX))(LevS->GArray[i].GSq == 0.0 ? 0.0 : 1/LevS->GArray[i].GSq);
	1214	}*/
	1215
	1216	/*if ((j==0 && i==0) \|\| (j==1 && i==1) \|\| (j==2 && i==2) \|\| (j==3 && i==3))
	1217	LevS->LPsi->LocalPsi[j][i].re = 1.0;*/
	1218	/*
	1219	if ( j==0 && P->Par.me_comm_ST_PsiT == 0 &&
	1220	fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1221	fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1222	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1223	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1224	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1225	}
	1226	if ( j==1 && P->Par.me_comm_ST_PsiT == 0 &&
	1227	fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1228	fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1229	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1230	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1231	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1232	}
	1233	if ( j==2 && P->Par.me_comm_ST_PsiT == 0 &&
	1234	fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1235	fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1236	fabs(LevS->GArray[i].G[2] - 1.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1237	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1238	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1239	}
	1240	if ( j==3 && P->Par.me_comm_ST_PsiT == 0 &&
	1241	fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1242	fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1243	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1244	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1245	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1246	}
	1247	*/
	1248	/*
	1249	if ( j==0 && P->Par.me_comm_ST_PsiT == 0 &&
	1250	fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1251	fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1252	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1253	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1254	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1255	}
	1256	if ( j==1 && P->Par.me_comm_ST_PsiT == 0 &&
	1257	fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1258	fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1259	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1260	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1261	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1262	}
	1263	if ( j==0 && P->Par.me_comm_ST_PsiT == 1 &&
	1264	fabs(LevS->GArray[i].G[0] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1265	fabs(LevS->GArray[i].G[1] - 0.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1266	fabs(LevS->GArray[i].G[2] - 1.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1267	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1268	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1269	}
	1270	if ( j==1 && P->Par.me_comm_ST_PsiT == 1 &&
	1271	fabs(LevS->GArray[i].G[0] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1272	fabs(LevS->GArray[i].G[1] - 1.0*Lat->ReciBasis[0]) < MYEPSILON &&
	1273	fabs(LevS->GArray[i].G[2] - 0.0*Lat->ReciBasis[0]) < MYEPSILON ) {
	1274	LevS->LPsi->LocalPsi[j][i].re = 1.0;
	1275	fprintf(stdout, "(%i) %i g = %i\n", P->Par.me, j, i);
	1276	}
	1277	*/
	1278	}
	1279	//fclose(psi);
	1280	}
	1281
	1282	/* Init Lattice End */
	1283
	1284
	1285	/** Reads parameter from a parsed file.
	1286	* The file is either parsed for a certain keyword or if null is given for
	1287	* the value in row yth and column xth. If the keyword was necessity#critical,
	1288	* then an error is thrown and the programme aborted.
	1289	* \warning value is modified (both in contents and position)!
	1290	* \param verbose 1 - print found value to stderr, 0 - don't
	1291	* \param file file to be parsed
	1292	* \param name Name of value in file (at least 3 chars!)
	1293	* \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
	1294	* (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
	1295	* best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
	1296	* counted from this unresetted position!)
	1297	* \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
	1298	* \param yth In grid case specifying column number, otherwise the yth \a name matching line
	1299	* \param type Type of the Parameter to be read
	1300	* \param value address of the value to be read (must have been allocated)
	1301	* \param critical necessity of this keyword being specified (optional, critical)
	1302	* \return 1 - found, 0 - not found
	1303	* \sa ReadMainParameterFile Calling routine
	1304	*/
	1305	int ParseForParameter(int verbose, FILE file, const char const name, int sequential, int const xth, int const yth, enum value_type type, void *value, enum necessity critical) {
	1306	int i,j; // loop variables
	1307	int me, length, maxlength = -1;
	1308	long file_position = ftell(file); // mark current position
	1309	char dummy1, dummy, *free_dummy; // pointers in the line that is read in per step
	1310	dummy1 = free_dummy = (char ) Malloc(256 sizeof(char),"ParseForParameter: MemAlloc for dummy1 failed");
	1311	MPI_Comm_rank(MPI_COMM_WORLD, &me);
	1312
	1313	//fprintf(stderr,"(%i) Parsing for %s\n",me,name);
	1314
	1315	int line = 0; // marks line where parameter was found
	1316	int found = (type >= grid) ? 0 : (-yth + 1); // marks if yth parameter name was found
	1317	while((found != 1)) {
	1318	dummy1 = dummy = free_dummy;
	1319	do {
	1320	fgets(dummy1, 256, file); // Read the whole line
	1321	if (feof(file)) {
	1322	if ((critical) && (found == 0)) {
	1323	Free(free_dummy);
	1324	Error(InitReading, name);
	1325	} else {
	1326	//if (!sequential)
	1327	fseek(file, file_position, SEEK_SET); // rewind to start position
	1328	Free(free_dummy);
	1329	return 0;
	1330	}
	1331	}
	1332	line++;
	1333	} while (((dummy1[0] == '#') \|\| (dummy1[0] == '\n')) && dummy != NULL); // skip commentary and empty lines
	1334
	1335	// C++ getline removes newline at end, thus re-add
	1336	if (strchr(dummy1,'\n') == NULL) {
	1337	i = strlen(dummy1);
	1338	dummy1[i] = '\n';
	1339	dummy1[i+1] = '\0';
	1340	}
	1341
	1342	if (dummy1 == NULL) {
	1343	if (verbose) fprintf(stderr,"(%i) Error reading line %i\n",me,line);
	1344	} else {
	1345	//fprintf(stderr,"(%i) Reading next line %i: %s\n",me, line, dummy1);
	1346	}
	1347	// Seek for possible end of keyword on line if given ...
	1348	if (name != NULL) {
	1349	dummy = strchr(dummy1,'\t'); // set dummy on first tab or space which ever nearer
	1350	if (dummy == NULL) {
	1351	dummy = strchr(dummy1, ' '); // if not found seek for space
	1352	while ((dummy != NULL) && ((dummy == '\t') \|\| (dummy == ' '))) // skip some more tabs and spaces if necessary
	1353	dummy++;
	1354	}
	1355	if (dummy == NULL) {
	1356	dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
	1357	//fprintf(stderr,"(%i)Error: Cannot find tabs or spaces on line %i in search for %s\n", me, line, name);
	1358	//Free(free_dummy);
	1359	//Error(FileOpenParams, NULL);
	1360	} else {
	1361	//fprintf(stderr,"(%i) found tab at %i\n",me,(char )dummy-(char )dummy1);
	1362	}
	1363	} else dummy = dummy1;
	1364	// ... and check if it is the keyword!
	1365	if ((name == NULL) \|\| ((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0))) {
	1366	found++; // found the parameter!
	1367	//fprintf(stderr,"(%i) found %s at line %i\n",me, name, line);
	1368
	1369	if (found == 1) {
	1370	for (i=0;i<xth;i++) { // i = rows
	1371	if (type >= grid) {
	1372	// grid structure means that grid starts on the next line, not right after keyword
	1373	dummy1 = dummy = free_dummy;
	1374	do {
	1375	fgets(dummy1, 256, file); // Read the whole line, skip commentary and empty ones
	1376	if (feof(file)) {
	1377	if ((critical) && (found == 0)) {
	1378	Free(free_dummy);
	1379	Error(InitReading, name);
	1380	} else {
	1381	//if (!sequential)
	1382	fseek(file, file_position, SEEK_SET); // rewind to start position
	1383	Free(free_dummy);
	1384	return 0;
	1385	}
	1386	}
	1387	line++;
	1388	} while ((dummy1[0] == '#') \|\| (dummy1[0] == '\n'));
	1389	if (dummy1 == NULL){
	1390	if (verbose) fprintf(stderr,"(%i) Error reading line %i\n", me, line);
	1391	} else {
	1392	//fprintf(stderr,"(%i) Reading next line %i: %s\n", me, line, dummy1);
	1393	}
	1394	} else { // simple int, strings or doubles start in the same line
	1395	while ((dummy == '\t') \|\| (dummy == ' ')) // skip interjacent tabs and spaces
	1396	dummy++;
	1397	}
	1398	// C++ getline removes newline at end, thus re-add
	1399	if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
	1400	j = strlen(dummy1);
	1401	dummy1[j] = '\n';
	1402	dummy1[j+1] = '\0';
	1403	}
	1404
	1405	int start = (type >= grid) ? 0 : yth-1 ;
	1406	for (j=start;j<yth;j++) { // j = columns
	1407	// check for lower triangular area and upper triangular area
	1408	if ( ((i > j) && (type == upper_trigrid)) \|\| ((j > i) && (type == lower_trigrid))) {
	1409	((double )value) = 0.0;
	1410	//fprintf(stderr,"%f\t",((double )value));
	1411	value += sizeof(double);
	1412	} else {
	1413	// otherwise we must skip all interjacent tabs and spaces and find next value
	1414	dummy1 = dummy;
	1415	dummy = strchr(dummy1, '\t'); // seek for tab or space
	1416	if (dummy == NULL) {
	1417	dummy = strchr(dummy1, ' '); // if not found seek for space
	1418	while ((dummy != NULL) && ((dummy == '\t') \|\| (dummy == ' '))) // skip some more tabs and spaces if necessary
	1419	dummy++;
	1420	}
	1421	/* while ((dummy != NULL) && ((dummy == '\t') \|\| (dummy == ' '))) // skip some more tabs and spaces if necessary
	1422	dummy++;*/
	1423	if (dummy == NULL) { // if still zero returned ...
	1424	dummy = strchr(dummy1, '\n'); // ... at line end then
	1425	if ((j < yth-1) && (type < 4)) { // check if xth value or not yet
	1426	if (verbose) fprintf(stderr,"(%i)Error: EoL at %i and still missing %i value(s) for parameter %s\n", me, line, yth-j, name);
	1427	Free(free_dummy);
	1428	return 0;
	1429	//Error(FileOpenParams, NULL);
	1430	}
	1431	} else {
	1432	//fprintf(stderr,"(%i) found tab at %i\n",me,(char )dummy-(char )free_dummy);
	1433	}
	1434	//fprintf(stderr,"(%i) value from %i to %i\n",me,(char )dummy1-(char )free_dummy,(char )dummy-(char )free_dummy);
	1435	switch(type) {
	1436	case (row_int):
	1437	((int )value) = atoi(dummy1);
	1438	if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"(%i) %s = ",me, name);
	1439	if (verbose) fprintf(stderr,"%i\t",((int )value));
	1440	value += sizeof(int);
	1441	break;
	1442	case (row_double):
	1443	case(grid):
	1444	case(lower_trigrid):
	1445	case(upper_trigrid):
	1446	((double )value) = atof(dummy1);
	1447	if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"(%i) %s = ",me, name);
	1448	if (verbose) fprintf(stderr,"%lg\t",((double )value));
	1449	value += sizeof(double);
	1450	break;
	1451	case(double_type):
	1452	((double )value) = atof(dummy1);
	1453	if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s = %lg\n",me, name, ((double ) value));
	1454	//value += sizeof(double);
	1455	break;
	1456	case(int_type):
	1457	((int )value) = atoi(dummy1);
	1458	if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s = %i\n",me, name, ((int ) value));
	1459	//value += sizeof(int);
	1460	break;
	1461	default:
	1462	case(string_type):
	1463	if (value != NULL) {
	1464	if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
	1465	length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
	1466	strncpy((char *)value, dummy1, length); // copy as much
	1467	((char *)value)[length] = '\0'; // and set end marker
	1468	if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s is '%s' (%i chars)\n",me,name,((char *) value), length);
	1469	//value += sizeof(char);
	1470	} else {
	1471	Error(SomeError, "ParseForParameter: given string array points to NULL!");
	1472	}
	1473	break;
	1474	}
	1475	}
	1476	while (*dummy == '\t')
	1477	dummy++;
	1478	}
	1479	}
	1480	}
	1481	}
	1482	}
	1483	if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
	1484	Free(free_dummy);
	1485	if (!sequential) fseek(file, file_position, SEEK_SET); // rewind to start position
	1486	//fprintf(stderr, "(%i) End of Parsing\n\n",me);
	1487
	1488	return (found); // true if found, false if not
	1489	}
	1490
	1491	/** Readin of main parameter file.
	1492	* creates RunStruct and FileData, opens \a filename \n
	1493	* Only for process 0: and then scans the whole thing.
	1494	*
	1495	* \param P Problem structure to be filled
	1496	* \param filename const char array with parameter file name
	1497	* \sa ParseForParameter
	1498	*/
	1499	void ReadParameters(struct Problem const P, const char const filename) {
	1500	/* Oeffne Hauptparameterdatei */
	1501	struct RunStruct *R = &P->R;
	1502	struct FileData *F = &P->Files;
	1503	FILE *file;
	1504	int i, di, me;
	1505	double BField[NDIM_NDIM];
	1506
	1507	MPI_Comm_rank(MPI_COMM_WORLD, &me);
	1508	if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadMainParameterFile\n", me);
	1509	file = fopen(filename, "r");
	1510	if (file == NULL) {
	1511	fprintf(stderr,"(%i)Error: Cannot open main parameter file: %s\n", me, filename);
	1512	Error(FileOpenParams, NULL);
	1513	} else if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)File is open: %s\n", me, filename);
	1514
	1515	/* Namen einlesen */
	1516
	1517	P->Files.filename = MallocString(MAXDUMMYSTRING,"ReadParameters: filename");
	1518	ParseForParameter(P->Call.out[ReadOut],file, "mainname", 0, 1, 1, string_type, P->Files.filename, critical);
	1519	//debug(P,"mainname");
	1520	CreateMainname(P, filename);
	1521	P->Files.default_path = MallocString(MAXDUMMYSTRING,"ReadParameters: default_path");
	1522	ParseForParameter(P->Call.out[ReadOut],file, "defaultpath", 0, 1, 1, string_type, P->Files.default_path, critical);
	1523	P->Files.pseudopot_path = MallocString(MAXDUMMYSTRING,"ReadParameters: pseudopot_path");
	1524	ParseForParameter(P->Call.out[ReadOut],file, "pseudopotpath", 0, 1, 1, string_type, P->Files.pseudopot_path, critical);
	1525	ParseForParameter(P->Call.out[ReadOut],file,"ProcPEGamma", 0, 1, 1, int_type, &(P->Par.proc[PEGamma]), critical);
	1526	ParseForParameter(P->Call.out[ReadOut],file,"ProcPEPsi", 0, 1, 1, int_type, &(P->Par.proc[PEPsi]), critical);
	1527	if (P->Call.proc[PEPsi]) { /* Kommandozeilenoption */
	1528	P->Par.proc[PEPsi] = P->Call.proc[PEPsi];
	1529	P->Par.proc[PEGamma] = P->Call.proc[PEGamma];
	1530	}
	1531	/*
	1532	Run
	1533	*/
	1534	if (!ParseForParameter(P->Call.out[ReadOut],file,"Seed", 0, 1, 1, int_type, &(R->Seed), optional))
	1535	R->Seed = 1;
	1536	if (!ParseForParameter(P->Call.out[ReadOut],file,"WaveNr", 0, 1, 1, int_type, &(R->WaveNr), optional))
	1537	R->WaveNr = 0;
	1538	if (!ParseForParameter(P->Call.out[ReadOut],file,"Lambda", 0, 1, 1, double_type, &R->Lambda, optional))
	1539	R->Lambda = 1.;
	1540
	1541	if(!ParseForParameter(P->Call.out[ReadOut],file,"DoOutOrbitals", 0, 1, 1, int_type, &F->DoOutOrbitals, optional)) {
	1542	F->DoOutOrbitals = 0;
	1543	} else {
	1544	if (F->DoOutOrbitals < 0) F->DoOutOrbitals = 0;
	1545	if (F->DoOutOrbitals > 1) F->DoOutOrbitals = 1;
	1546	}
	1547	ParseForParameter(P->Call.out[ReadOut],file,"DoOutVis", 0, 1, 1, int_type, &F->DoOutVis, critical);
	1548	if (F->DoOutVis < 0) F->DoOutVis = 0;
	1549	if (F->DoOutVis > 1) F->DoOutVis = 1;
	1550	if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorPlane", 0, 1, 1, int_type, &R->VectorPlane, optional))
	1551	R->VectorPlane = -1;
	1552	if (R->VectorPlane < -1 \|\| R->VectorPlane > 2) {
	1553	fprintf(stderr,"(%i) ! -1 <= R-VectorPlane <= 2: We simulate three-dimensional worlds! Disabling current vector plot: -1\n", P->Par.me);
	1554	R->VectorPlane = -1;
	1555	}
	1556	if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorCut", 0, 1, 1, double_type, &R->VectorCut, optional))
	1557	R->VectorCut = 0.;
	1558	ParseForParameter(P->Call.out[ReadOut],file,"DoOutMes", 0, 1, 1, int_type, &F->DoOutMes, critical);
	1559	if (F->DoOutMes < 0) F->DoOutMes = 0;
	1560	if (F->DoOutMes > 1) F->DoOutMes = 1;
	1561	if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutCurr", 0, 1, 1, int_type, &F->DoOutCurr, optional))
	1562	F->DoOutCurr = 0;
	1563	if (F->DoOutCurr < 0) F->DoOutCurr = 0;
	1564	if (F->DoOutCurr > 1) F->DoOutCurr = 1;
	1565	ParseForParameter(P->Call.out[ReadOut],file,"AddGramSch", 0, 1, 1, int_type, &R->UseAddGramSch, critical);
	1566	if (R->UseAddGramSch < 0) R->UseAddGramSch = 0;
	1567	if (R->UseAddGramSch > 2) R->UseAddGramSch = 2;
	1568	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)UseAddGramSch = %i\n",me,R->UseAddGramSch);
	1569	if(!ParseForParameter(P->Call.out[ReadOut],file,"CommonWannier", 0, 1, 1, int_type, &R->CommonWannier, optional)) {
	1570	R->CommonWannier = 0;
	1571	} else {
	1572	if (R->CommonWannier < 0) R->CommonWannier = 0;
	1573	if (R->CommonWannier > 4) R->CommonWannier = 4;
	1574	}
	1575	if(!ParseForParameter(P->Call.out[ReadOut],file,"SawtoothStart", 0, 1, 1, double_type, &P->Lat.SawtoothStart, optional)) {
	1576	P->Lat.SawtoothStart = 0.01;
	1577	} else {
	1578	if (P->Lat.SawtoothStart < 0.) P->Lat.SawtoothStart = 0.;
	1579	if (P->Lat.SawtoothStart > 1.) P->Lat.SawtoothStart = 1.;
	1580	}
	1581	if (!ParseForParameter(P->Call.out[ReadOut],file,"DoBrent", 0, 1, 1, int_type, &R->DoBrent, optional)) {
	1582	R->DoBrent = 0;
	1583	} else {
	1584	if (R->DoBrent < 0) R->DoBrent = 0;
	1585	if (R->DoBrent > 1) R->DoBrent = 1;
	1586	}
	1587
	1588
	1589	ParseForParameter(P->Call.out[ReadOut],file,"MaxOuterStep", 0, 1, 1, int_type, &R->MaxOuterStep, critical);
	1590	// the following values are MD specific and should be dropped in further revisions
	1591	if (!ParseForParameter(P->Call.out[ReadOut],file,"Deltat", 0, 1, 1, double_type, &R->delta_t, optional))
	1592	R->delta_t = 0;
	1593	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutVisStep", 0, 1, 1, int_type, &R->OutVisStep, optional))
	1594	R->OutVisStep = 10;
	1595	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutSrcStep", 0, 1, 1, int_type, &R->OutSrcStep, optional))
	1596	R->OutSrcStep = 5;
	1597	if (!ParseForParameter(P->Call.out[ReadOut],file,"TargetTemp", 0, 1, 1, double_type, &R->TargetTemp, optional))
	1598	R->TargetTemp = 0;
	1599	if (!ParseForParameter(P->Call.out[ReadOut],file,"ScaleTempStep", 0, 1, 1, int_type, &R->ScaleTempStep, optional))
	1600	R->ScaleTempStep = 25;
	1601
	1602	if (!ParseForParameter(P->Call.out[ReadOut],file,"EpsWannier", 0, 1, 1, double_type, &R->EpsWannier, optional))
	1603	R->EpsWannier = 1e-8;
	1604
	1605	// stop conditions
	1606	R->t = 0;
	1607	//if (R->MaxOuterStep <= 0) R->MaxOuterStep = 1;
	1608	if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxOuterStep = %i\n",me,R->MaxOuterStep);
	1609	if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)t = %g\n",me,R->t);
	1610	ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiStep", 0, 1, 1, int_type, &R->MaxPsiStep, critical);
	1611	if (R->MaxPsiStep <= 0) R->MaxPsiStep = 3;
	1612	if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxPsiStep = %i\n",me,R->MaxPsiStep);
	1613
	1614	ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStep", 0, 1, 1, int_type, &R->MaxMinStep, critical);
	1615	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsTotalE", 0, 1, 1, double_type, &R->RelEpsTotalEnergy, critical);
	1616	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsKineticE", 0, 1, 1, double_type, &R->RelEpsKineticEnergy, critical);
	1617	ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStopStep", 0, 1, 1, int_type, &R->MaxMinStopStep, critical);
	1618	ParseForParameter(P->Call.out[ReadOut],file,"MaxMinGapStopStep", 0, 1, 1, int_type, &R->MaxMinGapStopStep, critical);
	1619	if (R->MaxMinStep <= 0) R->MaxMinStep = R->MaxPsiStep;
	1620	if (R->MaxMinStopStep < 1) R->MaxMinStopStep = 1;
	1621	if (R->MaxMinGapStopStep < 1) R->MaxMinGapStopStep = 1;
	1622	if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxMinStep = %i\tRelEpsTotalEnergy: %e\tRelEpsKineticEnergy %e\tMaxMinStopStep %i\n",me,R->MaxMinStep,R->RelEpsTotalEnergy,R->RelEpsKineticEnergy,R->MaxMinStopStep);
	1623
	1624	ParseForParameter(P->Call.out[ReadOut],file,"MaxInitMinStep", 0, 1, 1, int_type, &R->MaxInitMinStep, critical);
	1625	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsTotalE", 0, 1, 1, double_type, &R->InitRelEpsTotalEnergy, critical);
	1626	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsKineticE", 0, 1, 1, double_type, &R->InitRelEpsKineticEnergy, critical);
	1627	ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinStopStep", 0, 1, 1, int_type, &R->InitMaxMinStopStep, critical);
	1628	ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &R->InitMaxMinGapStopStep, critical);
	1629	if (R->MaxInitMinStep <= 0) R->MaxInitMinStep = R->MaxPsiStep;
	1630	if (R->InitMaxMinStopStep < 1) R->InitMaxMinStopStep = 1;
	1631	if (R->InitMaxMinGapStopStep < 1) R->InitMaxMinGapStopStep = 1;
	1632	if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)INIT:MaxMinStep = %i\tRelEpsTotalEnergy: %e\tRelEpsKineticEnergy %e\tMaxMinStopStep %i\tMaxMinGapStopStep %i\n",me,R->MaxInitMinStep,R->InitRelEpsTotalEnergy,R->InitRelEpsKineticEnergy,R->InitMaxMinStopStep,R->InitMaxMinStopStep);
	1633
	1634	// Unit cell and magnetic field
	1635	ParseForParameter(P->Call.out[ReadOut],file, "BoxLength", 0, 3, 3, lower_trigrid, &P->Lat.RealBasis, critical); /* Lattice->RealBasis */
	1636	ParseForParameter(P->Call.out[ReadOut],file, "DoPerturbation", 0, 1, 1, int_type, &R->DoPerturbation, optional);
	1637	if (!R->DoPerturbation) {
	1638	if (!ParseForParameter(P->Call.out[ReadOut],file, "BField", 0, 1, 3, grid, &BField, optional)) { // old parameter for perturbation with field direction
	1639	if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) No perturbation specified.\n", P->Par.me);
	1640	if (F->DoOutCurr) {
	1641	if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) DoOutCurr= 1 though no DoPerturbation specified: setting DoOutCurr = 0\n", P->Par.me);
	1642	F->DoOutCurr = 0;
	1643	R->DoPerturbation = 0;
	1644	}
	1645	} else {
	1646	R->DoPerturbation = 1;
	1647	}
	1648	}
	1649	if (!P->Call.WriteSrcFiles && R->DoPerturbation) {
	1650	if (P->Call.out[ReadOut]) fprintf(stderr,"(%i)Perturbation: Always writing source files (\"-w\").\n", P->Par.me);
	1651	P->Call.WriteSrcFiles = 1;
	1652	}
	1653	RMatReci3(P->Lat.InvBasis, P->Lat.RealBasis);
	1654
	1655	if (!ParseForParameter(P->Call.out[ReadOut],file,"DoFullCurrent", 0, 1, 1, int_type, &R->DoFullCurrent, optional))
	1656	R->DoFullCurrent = 0;
	1657	if (R->DoFullCurrent < 0) R->DoFullCurrent = 0;
	1658	if (R->DoFullCurrent > 2) R->DoFullCurrent = 2;
	1659	if (R->DoPerturbation == 0) {
	1660	if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) No Magnetic field: RunStruct#DoFullCurrent = 0\n", P->Par.me);
	1661	R->DoFullCurrent = 0;
	1662	}
	1663
	1664	ParseForParameter(P->Call.out[ReadOut],file,"ECut", 0, 1, 1, double_type, &P->Lat.ECut, critical);
	1665	if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) ECut = %e -> %e Rydberg\n", me, P->Lat.ECut, 2.*P->Lat.ECut);
	1666	ParseForParameter(P->Call.out[ReadOut],file,"MaxLevel", 0, 1, 1, int_type, &P->Lat.MaxLevel, critical);
	1667	ParseForParameter(P->Call.out[ReadOut],file,"Level0Factor", 0, 1, 1, int_type, &P->Lat.Lev0Factor, critical);
	1668	if (P->Lat.Lev0Factor < 2) {
	1669	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set Lev0Factor to 2!\n",me);
	1670	P->Lat.Lev0Factor = 2;
	1671	}
	1672	ParseForParameter(P->Call.out[ReadOut],file,"RiemannTensor", 0, 1, 1, int_type, &di, critical);
	1673	if (di >= 0 && di < 2) {
	1674	P->Lat.RT.Use = (enum UseRiemannTensor)di;
	1675	} else {
	1676	Error(SomeError, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
	1677	}
	1678	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!RiemannTensor = %i\n",me,P->Lat.RT.Use), critical;
	1679	switch (P->Lat.RT.Use) {
	1680	case UseNotRT:
	1681	if (P->Lat.MaxLevel < 2) {
	1682	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set MaxLevel to 2!\n",me);
	1683	P->Lat.MaxLevel = 2;
	1684	}
	1685	P->Lat.LevRFactor = 2;
	1686	P->Lat.RT.ActualUse = inactive;
	1687	break;
	1688	case UseRT:
	1689	if (P->Lat.MaxLevel < 3) {
	1690	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set MaxLevel to 3!\n",me);
	1691	P->Lat.MaxLevel = 3;
	1692	}
	1693	ParseForParameter(P->Call.out[ReadOut],file,"RiemannLevel", 0, 1, 1, int_type, &P->Lat.RT.RiemannLevel, critical);
	1694	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!RiemannLevel = %i\n",me,P->Lat.RT.RiemannLevel);
	1695	if (P->Lat.RT.RiemannLevel < 2) {
	1696	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set RiemannLevel to 2!\n",me);
	1697	P->Lat.RT.RiemannLevel = 2;
	1698	}
	1699	if (P->Lat.RT.RiemannLevel > P->Lat.MaxLevel-1) {
	1700	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set RiemannLevel to %i (MaxLevel-1)!\n",me,P->Lat.MaxLevel-1);
	1701	P->Lat.RT.RiemannLevel = P->Lat.MaxLevel-1;
	1702	}
	1703	ParseForParameter(P->Call.out[ReadOut],file,"LevRFactor", 0, 1, 1, int_type, &P->Lat.LevRFactor, critical);
	1704	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!LevRFactor = %i\n",me,P->Lat.LevRFactor);
	1705	if (P->Lat.LevRFactor < 2) {
	1706	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set LevRFactor to 2!\n",me);
	1707	P->Lat.LevRFactor = 2;
	1708	}
	1709	P->Lat.Lev0Factor = 2;
	1710	P->Lat.RT.ActualUse = standby;
	1711	break;
	1712	}
	1713	ParseForParameter(P->Call.out[ReadOut],file,"PsiType", 0, 1, 1, int_type, &di, critical);
	1714	if (di >= 0 && di < 2) {
	1715	P->Lat.Psi.Use = (enum UseSpinType)di;
	1716	} else {
	1717	Error(SomeError, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
	1718	}
	1719	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiType = %i\n",me,P->Lat.Psi.Use);
	1720	for (i=0; i<MaxPsiNoType; i++)
	1721	P->Lat.Psi.GlobalNo[i] = 0;
	1722	switch (P->Lat.Psi.Use) {
	1723	case UseSpinDouble:
	1724	ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiDouble", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDouble], critical);
	1725	if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], optional))
	1726	P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
	1727	if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
	1728	R->DoUnOccupied = 1;
	1729	else
	1730	R->DoUnOccupied = 0;
	1731	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!MaxPsiDouble = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoDouble],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
	1732	//P->Lat.Psi.GlobalNo[PsiMaxNoDouble] += P->Lat.Psi.GlobalNo[PsiMaxAdd]; no more adding up on occupied ones
	1733	P->Lat.Psi.GlobalNo[PsiMaxNo] = P->Lat.Psi.GlobalNo[PsiMaxNoDouble];
	1734	break;
	1735	case UseSpinUpDown:
	1736	if (P->Par.proc[PEPsi] % 2) Error(SomeError,"UseSpinUpDown & P->Par.proc[PEGamma] % 2");
	1737	ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoUp", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoUp], critical);
	1738	ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoDown", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDown], critical);
	1739	if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], optional))
	1740	P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
	1741	if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
	1742	R->DoUnOccupied = 1;
	1743	else
	1744	R->DoUnOccupied = 0;
	1745	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiMaxPsiUp = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoUp],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
	1746	if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiMaxPsiDown = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoDown],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
	1747	//P->Lat.Psi.GlobalNo[PsiMaxNoUp] += P->Lat.Psi.GlobalNo[PsiMaxAdd];
	1748	//P->Lat.Psi.GlobalNo[PsiMaxNoDown] += P->Lat.Psi.GlobalNo[PsiMaxAdd];
	1749	if (abs(P->Lat.Psi.GlobalNo[PsiMaxNoUp]-P->Lat.Psi.GlobalNo[PsiMaxNoDown])>1) Error(SomeError,"\|Up - Down\| > 1");
	1750	P->Lat.Psi.GlobalNo[PsiMaxNo] = P->Lat.Psi.GlobalNo[PsiMaxNoUp] + P->Lat.Psi.GlobalNo[PsiMaxNoDown];
	1751	break;
	1752	}
	1753
	1754	IonsInitRead(P, file);
	1755	fclose(file);
	1756	}
	1757
	1758	/** Main Initialization.
	1759	* Massive calling of initializing super-functions init...():
	1760	* -# InitOutputFiles()\n
	1761	* Opening and Initializing of output measurement files.
	1762	* -# InitLattice()\n
	1763	* Initializing Lattice structure.
	1764	* -# InitRunLevel()\n
	1765	* Initializing RunLevel structure.
	1766	* -# InitLatticeLevel0()\n
	1767	* Initializing LatticeLevel structure.
	1768	* -# InitOtherLevel()\n
	1769	* Initialization of first to Lattice::MaxLevel LatticeLevel's.
	1770	* -# InitRun()\n
	1771	* Initialization of RunStruct structure.
	1772	* -# InitPsis()\n
	1773	* Initialization of wave functions as a whole.
	1774	* -# InitDensity()\n
	1775	* Initializing real and complex Density arrays.
	1776	* -# either ReadSrcFiles() or InitPsisValues()\n
	1777	* read old calculations from source file if desired or initialise with
	1778	* random values InitPsisValue().
	1779	* -# InitRiemannTensor()\n
	1780	* Initializing RiemannTensor structure.
	1781	* -# FirstInitGramSchData()\n
	1782	* Sets up OnePsiElement as an MPI structure and subsequently initialises all of them.
	1783	*
	1784	* Now follow various (timed) tests:
	1785	* -# GramSch()\n
	1786	* (twice) to orthonormalize the random values
	1787	* -# TestGramSch()\n
	1788	* Test if orbitals now fulfill kronecker delta relation.
	1789	* -# InitGradients()\n
	1790	* Initializing Gradient structure, allocating its arrays.
	1791	* -# InitDensityCalculation()\n
	1792	* Calculate naive, real and complex densities and sum them up.
	1793	* -# InitPseudoRead()\n
	1794	* Read PseudoPot'entials from file.
	1795	* -# InitEnergyArray()\n
	1796	* Initialize Energy arrays.
	1797	* -# InitPsiEnergyCalculation()\n
	1798	* Calculate all psi energies and sum up.
	1799	* -# CalculateDensityEnergy()\n
	1800	* Calculate Density energy.
	1801	* -# CalculateIonsEnergy()\n
	1802	* Calculate ionic energy.
	1803	* -# CalculateEwald()\n
	1804	* Calculate ionic ewald summation.
	1805	* -# EnergyAllReduce()\n
	1806	* Gather energies from all processes and sum up.
	1807	* -# EnergyOutput()\n
	1808	* First output to file.
	1809	* -# InitOutVisArray()\n
	1810	* First output of visual for Opendx.
	1811	*
	1812	* \param *P Problem at hand
	1813	*/
	1814	void Init(struct Problem *const P) {
	1815	struct RunStruct *R = &P->R;
	1816	InitOutputFiles(P);
	1817	InitLattice(P);
	1818	InitRunLevel(P);
	1819	InitLatticeLevel0(P);
	1820	InitOtherLevel(P);
	1821	InitPsiGroupNo(P);
	1822	InitRun(P);
	1823	InitPsis(P);
	1824	InitDensity(P);
	1825	if (P->Call.ReadSrcFiles) {
	1826	if (ReadSrcIons(P) && P->Call.out[NormalOut])
	1827	fprintf(stderr,"(%i) Ionic structure read successfully.\n", P->Par.me);
	1828	else
	1829	fprintf(stderr,"(%i) No readible Ionic structure found.\n", P->Par.me);
	1830	} else {
	1831	InitPsisValue(P, 0, P->Lat.Psi.TypeStartIndex[Perturbed_P0]);
	1832	}
	1833	InitRiemannTensor(P);
	1834	if (P->Lat.RT.ActualUse == standby && R->LevRSNo == R->LevSNo) P->Lat.RT.ActualUse = active;
	1835	CalculateRiemannTensorData(P);
	1836	FirstInitGramSchData(P,&P->Lat.Psi);
	1837	/* Test */
	1838	SpeedMeasure(P, InitSimTime, StartTimeDo);
	1839	SpeedMeasure(P, InitGramSchTime, StartTimeDo);
	1840	GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
	1841	ResetGramSch(P, &P->Lat.Psi);
	1842	if (R->DoUnOccupied == 1) {
	1843	ResetGramSchTagType(P, &P->Lat.Psi, UnOccupied, NotOrthogonal);
	1844	GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
	1845	}
	1846	SpeedMeasure(P, InitGramSchTime, StopTimeDo);
	1847	SpeedMeasure(P, InitSimTime, StopTimeDo);
	1848	TestGramSch(P, R->LevS, &P->Lat.Psi, -1);
	1849	ResetGramSchTagType(P, &P->Lat.Psi, UnOccupied, NotUsedToOrtho);
	1850	InitGradients(P, &P->Grad);
	1851	SpeedMeasure(P, InitSimTime, StartTimeDo);
	1852	SpeedMeasure(P, InitDensityTime, StartTimeDo);
	1853	InitDensityCalculation(P);
	1854	SpeedMeasure(P, InitDensityTime, StopTimeDo);
	1855	SpeedMeasure(P, InitSimTime, StopTimeDo);
	1856	InitPseudoRead(P, P->Files.pseudopot_path);
	1857	InitEnergyArray(P);
	1858	SpeedMeasure(P, InitSimTime, StartTimeDo);
	1859	InitPsiEnergyCalculation(P, Occupied); /* STANDARTLEVEL */
	1860	CalculateDensityEnergy(P, 1); /* STANDARTLEVEL */
	1861	CalculateIonsEnergy(P);
	1862	CalculateEwald(P, 1);
	1863	EnergyAllReduce(P);
	1864	/* SetCurrentMinState(P,UnOccupied);
	1865	InitPsiEnergyCalculation(P,UnOccupied);
	1866	CalculateGapEnergy(P);
	1867	EnergyAllReduce(P);
	1868	SetCurrentMinState(P,Occupied);*/
	1869	SpeedMeasure(P, InitSimTime, StopTimeDo);
	1870	R->LevS->Step++;
	1871	EnergyOutput(P,0);
	1872	InitOutVisArray(P);
	1873	Free(P->Lat.GArrayForSort);
	1874	P->Speed.InitSteps++;
	1875	}
	1876

Note: See TracBrowser for help on using the repository browser.

Download in other formats: