source: pcp/src/init.c@ e6a122

Last change on this file since e6a122 was 774ae8, checked in by Frederik Heber <heber@…>, 18 years ago

WriteParameters() added
OutputIonCoordinates() uses WriteParameters() and additional parameter whether its MD or StructOpt
  • Property mode set to 100644
File size: 97.2 KB
RevLine 
[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 */
63static 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 */
77void 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,"/");
[64fa9e]108 if (dummy) Free(dummy, "CreateMainname: dummy");
[a0bcf1]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 */
126static 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] = UpDummy*NDIM_NDIM*NDIM;
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[i*NDIM+0],Lat->ReciBasis[i*NDIM+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 */
301static int TestG(struct Lattice *Lat, int g1, int g2, int g3, double *GSq, double G[NDIM]) {
302 G[0] = g1*Lat->ReciBasis[0]+g2*Lat->ReciBasis[3]+g3*Lat->ReciBasis[6];
303 G[1] = g1*Lat->ReciBasis[1]+g2*Lat->ReciBasis[4]+g3*Lat->ReciBasis[7];
304 G[2] = g1*Lat->ReciBasis[2]+g2*Lat->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 */
321static 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 */
341static 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 */
351static 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 */
362static 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 */
379static 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 */
410static 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(MaxPE*MaxG*(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
[64fa9e]522 Free(GlobalGArray, "HashGArray: GlobalGArray");
523 Free(GlobalGArray_sort, "HashGArray: GlobalGArray_sort");
[a0bcf1]524 for (i=0;i<NDIM;i++) {
[64fa9e]525 Free(G[i], "HashGArray: G[i]");
526 Free(G_send[i], "HashGArray: G_send[i]");
[a0bcf1]527 }
[64fa9e]528 Free(G, "HashGArray: G");
529 Free(G_send, "HashGArray: G_send");
[a0bcf1]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 */
538static 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(2*sizeof(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[2*Lev0->MaxDoubleG+1] < 0) Error(SomeError,"InitLevel0: Lev0->DoubleG[2*Lev->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 */
744static 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 */
800static 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->MaxG*Lev->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->MaxG*Lev->NUp1[0]*Lev->NUp1[1]*Lev->NUp1[2], "InitLatticeLevel: RTPosFactors1");
867 Lat->RT.PosFactor[RTPFRto0] = (fftw_complex *)
868 Malloc(sizeof(fftw_complex)*Lev->MaxG*Lev->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(2*sizeof(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[2*Lev->MaxDoubleG+1] < 0) Error(SomeError,"Lev->DoubleG[2*Lev->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 */
967static 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 */
1062static 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 */
1169void 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)
[961b75]1301 * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
[a0bcf1]1302 * \param critical necessity of this keyword being specified (optional, critical)
1303 * \return 1 - found, 0 - not found
1304 * \sa ReadMainParameterFile Calling routine
1305 */
[961b75]1306int ParseForParameter(int verbose, FILE *file, const char *const name, int sequential, int const xth, int const yth, enum value_type type, void *value, int repetition, enum necessity critical) {
[a0bcf1]1307 int i,j; // loop variables
1308 int me, length, maxlength = -1;
1309 long file_position = ftell(file); // mark current position
[961b75]1310 char *dummy2, *dummy1, *dummy, *free_dummy; // pointers in the line that is read in per step
1311
1312 if (repetition == 0)
1313 Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
1314
1315 // allocated memory and rank in multi-process word
[a0bcf1]1316 dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char),"ParseForParameter: MemAlloc for dummy1 failed");
1317 MPI_Comm_rank(MPI_COMM_WORLD, &me);
1318
1319 //fprintf(stderr,"(%i) Parsing for %s\n",me,name);
1320 int line = 0; // marks line where parameter was found
1321 int found = (type >= grid) ? 0 : (-yth + 1); // marks if yth parameter name was found
[961b75]1322 while((found != repetition)) {
[a0bcf1]1323 dummy1 = dummy = free_dummy;
1324 do {
1325 fgets(dummy1, 256, file); // Read the whole line
1326 if (feof(file)) {
1327 if ((critical) && (found == 0)) {
[64fa9e]1328 Free(free_dummy, "ParseForParameter: free_dummy");
[a0bcf1]1329 Error(InitReading, name);
1330 } else {
1331 //if (!sequential)
1332 fseek(file, file_position, SEEK_SET); // rewind to start position
[64fa9e]1333 Free(free_dummy, "ParseForParameter: free_dummy");
[a0bcf1]1334 return 0;
1335 }
1336 }
1337 line++;
1338 } while (((dummy1[0] == '#') || (dummy1[0] == '\n')) && dummy != NULL); // skip commentary and empty lines
1339
1340 // C++ getline removes newline at end, thus re-add
1341 if (strchr(dummy1,'\n') == NULL) {
1342 i = strlen(dummy1);
1343 dummy1[i] = '\n';
1344 dummy1[i+1] = '\0';
1345 }
1346
1347 if (dummy1 == NULL) {
1348 if (verbose) fprintf(stderr,"(%i) Error reading line %i\n",me,line);
1349 } else {
1350 //fprintf(stderr,"(%i) Reading next line %i: %s\n",me, line, dummy1);
1351 }
1352 // Seek for possible end of keyword on line if given ...
1353 if (name != NULL) {
1354 dummy = strchr(dummy1,'\t'); // set dummy on first tab or space which ever nearer
[961b75]1355 dummy2 = strchr(dummy1, ' '); // if not found seek for space
1356 if ((dummy != NULL) || (dummy2 != NULL)) {
1357 if ((dummy == NULL) || ((dummy2 != NULL) && (dummy2 < dummy)))
1358 dummy = dummy2;
[a0bcf1]1359 }
1360 if (dummy == NULL) {
1361 dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
1362 //fprintf(stderr,"(%i)Error: Cannot find tabs or spaces on line %i in search for %s\n", me, line, name);
[64fa9e]1363 //Free(free_dummy, "ParseForParameter: free_dummy");
[a0bcf1]1364 //Error(FileOpenParams, NULL);
1365 } else {
[961b75]1366 //fprintf(stderr,"(%i) found tab at %li\n",me,(long)((char *)dummy-(char *)dummy1));
[a0bcf1]1367 }
1368 } else dummy = dummy1;
1369 // ... and check if it is the keyword!
[961b75]1370 if ((name == NULL) || ((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0) && (dummy-dummy1 == strlen(name)))) {
[a0bcf1]1371 found++; // found the parameter!
1372 //fprintf(stderr,"(%i) found %s at line %i\n",me, name, line);
1373
[961b75]1374 if (found == repetition) {
[a0bcf1]1375 for (i=0;i<xth;i++) { // i = rows
1376 if (type >= grid) {
1377 // grid structure means that grid starts on the next line, not right after keyword
1378 dummy1 = dummy = free_dummy;
1379 do {
1380 fgets(dummy1, 256, file); // Read the whole line, skip commentary and empty ones
1381 if (feof(file)) {
1382 if ((critical) && (found == 0)) {
[64fa9e]1383 Free(free_dummy, "ParseForParameter: free_dummy");
[a0bcf1]1384 Error(InitReading, name);
1385 } else {
1386 //if (!sequential)
1387 fseek(file, file_position, SEEK_SET); // rewind to start position
[64fa9e]1388 Free(free_dummy, "ParseForParameter: free_dummy");
[a0bcf1]1389 return 0;
1390 }
1391 }
1392 line++;
1393 } while ((dummy1[0] == '#') || (dummy1[0] == '\n'));
1394 if (dummy1 == NULL){
1395 if (verbose) fprintf(stderr,"(%i) Error reading line %i\n", me, line);
1396 } else {
1397 //fprintf(stderr,"(%i) Reading next line %i: %s\n", me, line, dummy1);
1398 }
1399 } else { // simple int, strings or doubles start in the same line
1400 while ((*dummy == '\t') || (*dummy == ' ')) // skip interjacent tabs and spaces
1401 dummy++;
1402 }
1403 // C++ getline removes newline at end, thus re-add
1404 if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
1405 j = strlen(dummy1);
1406 dummy1[j] = '\n';
1407 dummy1[j+1] = '\0';
1408 }
1409
1410 int start = (type >= grid) ? 0 : yth-1 ;
1411 for (j=start;j<yth;j++) { // j = columns
1412 // check for lower triangular area and upper triangular area
1413 if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
1414 *((double *)value) = 0.0;
1415 //fprintf(stderr,"%f\t",*((double *)value));
1416 value += sizeof(double);
1417 } else {
1418 // otherwise we must skip all interjacent tabs and spaces and find next value
1419 dummy1 = dummy;
[961b75]1420 dummy = strchr(dummy1, '\t'); // seek for tab or space (space if appears sooner)
1421 dummy2 = strchr(dummy1, ' ');
1422 //fprintf(stderr,"dummy (%p) is %c\t dummy2 (%p) is %c\n", dummy, (dummy == NULL ? '-' : *dummy), dummy2, (dummy2 == NULL ? '-' : *dummy2));
1423 if ((dummy != NULL) || (dummy2 != NULL)) {
1424 if ((dummy == NULL) || ((dummy2 != NULL) && (dummy2 < dummy)))
1425 dummy = dummy2;
1426 //while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' '))) // skip some more tabs and spaces if necessary
1427 // dummy++;
[a0bcf1]1428 }
1429/* while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' '))) // skip some more tabs and spaces if necessary
1430 dummy++;*/
[961b75]1431 //fprintf(stderr,"dummy is %c\n", *dummy);
1432 dummy2 = strchr(dummy1, '#');
1433 if ((dummy == NULL) || ((dummy2 != NULL) && (dummy1 >= dummy2))) { // if still zero returned or in comment area ...
[a0bcf1]1434 dummy = strchr(dummy1, '\n'); // ... at line end then
[961b75]1435 if ((i < xth-1) && (type < 4)) { // check if xth value or not yet
1436 if (critical) {
1437 if (verbose) fprintf(stderr,"(%i)Error: EoL or comment at %i and still missing %i value(s) for parameter %s\n", me, line, yth-j, name);
1438 Free(free_dummy, "ParseForParameter: free_dummy");
1439 Error(InitReading, name);
1440 } else {
1441 //if (!sequential)
1442 fseek(file, file_position, SEEK_SET); // rewind to start position
1443 Free(free_dummy, "ParseForParameter: free_dummy");
1444 return 0;
1445 }
[a0bcf1]1446 //Error(FileOpenParams, NULL);
1447 }
1448 } else {
1449 //fprintf(stderr,"(%i) found tab at %i\n",me,(char *)dummy-(char *)free_dummy);
1450 }
1451 //fprintf(stderr,"(%i) value from %i to %i\n",me,(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
1452 switch(type) {
1453 case (row_int):
1454 *((int *)value) = atoi(dummy1);
1455 if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"(%i) %s = ",me, name);
1456 if (verbose) fprintf(stderr,"%i\t",*((int *)value));
1457 value += sizeof(int);
1458 break;
1459 case (row_double):
1460 case(grid):
1461 case(lower_trigrid):
1462 case(upper_trigrid):
1463 *((double *)value) = atof(dummy1);
1464 if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"(%i) %s = ",me, name);
1465 if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
1466 value += sizeof(double);
1467 break;
1468 case(double_type):
1469 *((double *)value) = atof(dummy1);
1470 if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s = %lg\n",me, name, *((double *) value));
1471 //value += sizeof(double);
1472 break;
1473 case(int_type):
1474 *((int *)value) = atoi(dummy1);
1475 if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s = %i\n",me, name, *((int *) value));
1476 //value += sizeof(int);
1477 break;
1478 default:
1479 case(string_type):
1480 if (value != NULL) {
1481 if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
1482 length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
1483 strncpy((char *)value, dummy1, length); // copy as much
1484 ((char *)value)[length] = '\0'; // and set end marker
1485 if ((verbose) && (i == xth-1)) fprintf(stderr,"(%i) %s is '%s' (%i chars)\n",me,name,((char *) value), length);
1486 //value += sizeof(char);
1487 } else {
1488 Error(SomeError, "ParseForParameter: given string array points to NULL!");
1489 }
1490 break;
1491 }
1492 }
1493 while (*dummy == '\t')
1494 dummy++;
1495 }
1496 }
1497 }
1498 }
1499 }
1500 if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
1501 if (!sequential) fseek(file, file_position, SEEK_SET); // rewind to start position
[961b75]1502 Free(free_dummy, "ParseForParameter: free_dummy");
[a0bcf1]1503 //fprintf(stderr, "(%i) End of Parsing\n\n",me);
1504
1505 return (found); // true if found, false if not
1506}
1507
1508/** Readin of main parameter file.
1509 * creates RunStruct and FileData, opens \a filename \n
1510 * Only for process 0: and then scans the whole thing.
1511 *
1512 * \param P Problem structure to be filled
1513 * \param filename const char array with parameter file name
1514 * \sa ParseForParameter
1515 */
1516void ReadParameters(struct Problem *const P, const char *const filename) {
1517 /* Oeffne Hauptparameterdatei */
1518 struct RunStruct *R = &P->R;
1519 struct FileData *F = &P->Files;
1520 FILE *file;
1521 int i, di, me;
[961b75]1522 //double BField[NDIM];
[a0bcf1]1523
1524 MPI_Comm_rank(MPI_COMM_WORLD, &me);
1525 if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)ReadMainParameterFile\n", me);
1526 file = fopen(filename, "r");
1527 if (file == NULL) {
1528 fprintf(stderr,"(%i)Error: Cannot open main parameter file: %s\n", me, filename);
1529 Error(FileOpenParams, NULL);
1530 } else if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)File is open: %s\n", me, filename);
1531
1532 /* Namen einlesen */
1533
1534 P->Files.filename = MallocString(MAXDUMMYSTRING,"ReadParameters: filename");
[961b75]1535 ParseForParameter(P->Call.out[ReadOut],file, "mainname", 0, 1, 1, string_type, P->Files.filename, 1, critical);
[a0bcf1]1536 //debug(P,"mainname");
1537 CreateMainname(P, filename);
1538 P->Files.default_path = MallocString(MAXDUMMYSTRING,"ReadParameters: default_path");
[961b75]1539 ParseForParameter(P->Call.out[ReadOut],file, "defaultpath", 0, 1, 1, string_type, P->Files.default_path, 1, critical);
[a0bcf1]1540 P->Files.pseudopot_path = MallocString(MAXDUMMYSTRING,"ReadParameters: pseudopot_path");
[961b75]1541 ParseForParameter(P->Call.out[ReadOut],file, "pseudopotpath", 0, 1, 1, string_type, P->Files.pseudopot_path, 1, critical);
1542 ParseForParameter(P->Call.out[ReadOut],file,"ProcPEGamma", 0, 1, 1, int_type, &(P->Par.proc[PEGamma]), 1, critical);
1543 ParseForParameter(P->Call.out[ReadOut],file,"ProcPEPsi", 0, 1, 1, int_type, &(P->Par.proc[PEPsi]), 1, critical);
[a0bcf1]1544 if (P->Call.proc[PEPsi]) { /* Kommandozeilenoption */
1545 P->Par.proc[PEPsi] = P->Call.proc[PEPsi];
1546 P->Par.proc[PEGamma] = P->Call.proc[PEGamma];
1547 }
1548 /*
1549 Run
1550 */
[961b75]1551 if (!ParseForParameter(P->Call.out[ReadOut],file,"Seed", 0, 1, 1, int_type, &(R->Seed), 1, optional))
[a0bcf1]1552 R->Seed = 1;
[961b75]1553 if (!ParseForParameter(P->Call.out[ReadOut],file,"WaveNr", 0, 1, 1, int_type, &(R->WaveNr), 1, optional))
[a0bcf1]1554 R->WaveNr = 0;
[961b75]1555 if (!ParseForParameter(P->Call.out[ReadOut],file,"Lambda", 0, 1, 1, double_type, &R->Lambda, 1, optional))
[a0bcf1]1556 R->Lambda = 1.;
1557
[961b75]1558 if(!ParseForParameter(P->Call.out[ReadOut],file,"DoOutOrbitals", 0, 1, 1, int_type, &F->DoOutOrbitals, 1, optional)) {
[a0bcf1]1559 F->DoOutOrbitals = 0;
1560 } else {
1561 if (F->DoOutOrbitals < 0) F->DoOutOrbitals = 0;
1562 if (F->DoOutOrbitals > 1) F->DoOutOrbitals = 1;
1563 }
[961b75]1564 ParseForParameter(P->Call.out[ReadOut],file,"DoOutVis", 0, 1, 1, int_type, &F->DoOutVis, 1, critical);
[a0bcf1]1565 if (F->DoOutVis < 0) F->DoOutVis = 0;
[961b75]1566 if (F->DoOutVis > 2) F->DoOutVis = 1;
1567 if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorPlane", 0, 1, 1, int_type, &R->VectorPlane, 1, optional))
[a0bcf1]1568 R->VectorPlane = -1;
1569 if (R->VectorPlane < -1 || R->VectorPlane > 2) {
1570 fprintf(stderr,"(%i) ! -1 <= R-VectorPlane <= 2: We simulate three-dimensional worlds! Disabling current vector plot: -1\n", P->Par.me);
1571 R->VectorPlane = -1;
1572 }
[961b75]1573 if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorCut", 0, 1, 1, double_type, &R->VectorCut, 1, optional))
[a0bcf1]1574 R->VectorCut = 0.;
[961b75]1575 ParseForParameter(P->Call.out[ReadOut],file,"DoOutMes", 0, 1, 1, int_type, &F->DoOutMes, 1, critical);
[a0bcf1]1576 if (F->DoOutMes < 0) F->DoOutMes = 0;
1577 if (F->DoOutMes > 1) F->DoOutMes = 1;
[961b75]1578 if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutCurr", 0, 1, 1, int_type, &F->DoOutCurr, 1, optional))
[a0bcf1]1579 F->DoOutCurr = 0;
[961b75]1580 if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutNICS", 0, 1, 1, int_type, &F->DoOutNICS, 1, optional))
1581 F->DoOutNICS = 0;
[a0bcf1]1582 if (F->DoOutCurr < 0) F->DoOutCurr = 0;
1583 if (F->DoOutCurr > 1) F->DoOutCurr = 1;
[961b75]1584 ParseForParameter(P->Call.out[ReadOut],file,"AddGramSch", 0, 1, 1, int_type, &R->UseAddGramSch, 1, critical);
[a0bcf1]1585 if (R->UseAddGramSch < 0) R->UseAddGramSch = 0;
1586 if (R->UseAddGramSch > 2) R->UseAddGramSch = 2;
1587 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)UseAddGramSch = %i\n",me,R->UseAddGramSch);
[961b75]1588 if(!ParseForParameter(P->Call.out[ReadOut],file,"CommonWannier", 0, 1, 1, int_type, &R->CommonWannier, 1, optional)) {
[a0bcf1]1589 R->CommonWannier = 0;
1590 } else {
1591 if (R->CommonWannier < 0) R->CommonWannier = 0;
1592 if (R->CommonWannier > 4) R->CommonWannier = 4;
1593 }
[961b75]1594 if(!ParseForParameter(P->Call.out[ReadOut],file,"SawtoothStart", 0, 1, 1, double_type, &P->Lat.SawtoothStart, 1, optional)) {
[a0bcf1]1595 P->Lat.SawtoothStart = 0.01;
1596 } else {
1597 if (P->Lat.SawtoothStart < 0.) P->Lat.SawtoothStart = 0.;
1598 if (P->Lat.SawtoothStart > 1.) P->Lat.SawtoothStart = 1.;
1599 }
[961b75]1600 if (!ParseForParameter(P->Call.out[ReadOut],file,"DoBrent", 0, 1, 1, int_type, &R->DoBrent, 1, optional)) {
[a0bcf1]1601 R->DoBrent = 0;
1602 } else {
1603 if (R->DoBrent < 0) R->DoBrent = 0;
1604 if (R->DoBrent > 1) R->DoBrent = 1;
1605 }
1606
1607
[961b75]1608 if (!ParseForParameter(P->Call.out[ReadOut],file,"MaxOuterStep", 0, 1, 1, int_type, &R->MaxOuterStep, 1, optional))
1609 R->MaxOuterStep = 0;
1610 if (!ParseForParameter(P->Call.out[ReadOut],file,"MaxStructOptStep", 0, 1, 1, int_type, &R->MaxStructOptStep, 1, optional)) {
1611 if (R->MaxOuterStep != 0) { // if StructOptStep not explicitely specified, then use OuterStep
1612 R->MaxStructOptStep = R->MaxOuterStep;
1613 } else {
1614 R->MaxStructOptStep = 100;
1615 }
1616 }
[a0bcf1]1617 // the following values are MD specific and should be dropped in further revisions
[961b75]1618 if (!ParseForParameter(P->Call.out[ReadOut],file,"Deltat", 0, 1, 1, double_type, &R->delta_t, 1, optional))
1619 R->delta_t = 1.;
1620 if (!ParseForParameter(P->Call.out[ReadOut],file,"OutVisStep", 0, 1, 1, int_type, &R->OutVisStep, 1, optional))
[a0bcf1]1621 R->OutVisStep = 10;
[961b75]1622 if (!ParseForParameter(P->Call.out[ReadOut],file,"OutSrcStep", 0, 1, 1, int_type, &R->OutSrcStep, 1, optional))
[a0bcf1]1623 R->OutSrcStep = 5;
[961b75]1624 if (!ParseForParameter(P->Call.out[ReadOut],file,"TargetTemp", 0, 1, 1, double_type, &R->TargetTemp, 1, optional))
[a0bcf1]1625 R->TargetTemp = 0;
1626
[961b75]1627 if (!ParseForParameter(P->Call.out[ReadOut],file,"EpsWannier", 0, 1, 1, double_type, &R->EpsWannier, 1, optional))
[a0bcf1]1628 R->EpsWannier = 1e-8;
1629
1630 // stop conditions
1631 R->t = 0;
1632 //if (R->MaxOuterStep <= 0) R->MaxOuterStep = 1;
1633 if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxOuterStep = %i\n",me,R->MaxOuterStep);
[961b75]1634 if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)Deltat = %g\n",me,R->delta_t);
1635 ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiStep", 0, 1, 1, int_type, &R->MaxPsiStep, 1, critical);
[a0bcf1]1636 if (R->MaxPsiStep <= 0) R->MaxPsiStep = 3;
1637 if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxPsiStep = %i\n",me,R->MaxPsiStep);
1638
[961b75]1639 ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStep", 0, 1, 1, int_type, &R->MaxMinStep, 1, critical);
1640 ParseForParameter(P->Call.out[ReadOut],file,"RelEpsTotalE", 0, 1, 1, double_type, &R->RelEpsTotalEnergy, 1, critical);
1641 ParseForParameter(P->Call.out[ReadOut],file,"RelEpsKineticE", 0, 1, 1, double_type, &R->RelEpsKineticEnergy, 1, critical);
1642 ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStopStep", 0, 1, 1, int_type, &R->MaxMinStopStep, 1, critical);
1643 ParseForParameter(P->Call.out[ReadOut],file,"MaxMinGapStopStep", 0, 1, 1, int_type, &R->MaxMinGapStopStep, 1, critical);
[a0bcf1]1644 if (R->MaxMinStep <= 0) R->MaxMinStep = R->MaxPsiStep;
1645 if (R->MaxMinStopStep < 1) R->MaxMinStopStep = 1;
1646 if (R->MaxMinGapStopStep < 1) R->MaxMinGapStopStep = 1;
1647 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);
1648
[961b75]1649 ParseForParameter(P->Call.out[ReadOut],file,"MaxInitMinStep", 0, 1, 1, int_type, &R->MaxInitMinStep, 1, critical);
1650 ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsTotalE", 0, 1, 1, double_type, &R->InitRelEpsTotalEnergy, 1, critical);
1651 ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsKineticE", 0, 1, 1, double_type, &R->InitRelEpsKineticEnergy, 1, critical);
1652 ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinStopStep", 0, 1, 1, int_type, &R->InitMaxMinStopStep, 1, critical);
1653 ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &R->InitMaxMinGapStopStep, 1, critical);
[a0bcf1]1654 if (R->MaxInitMinStep <= 0) R->MaxInitMinStep = R->MaxPsiStep;
1655 if (R->InitMaxMinStopStep < 1) R->InitMaxMinStopStep = 1;
1656 if (R->InitMaxMinGapStopStep < 1) R->InitMaxMinGapStopStep = 1;
1657 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);
1658
1659 // Unit cell and magnetic field
[961b75]1660 ParseForParameter(P->Call.out[ReadOut],file, "BoxLength", 0, 3, 3, lower_trigrid, &P->Lat.RealBasis, 1, critical); /* Lattice->RealBasis */
1661 ParseForParameter(P->Call.out[ReadOut],file, "DoPerturbation", 0, 1, 1, int_type, &R->DoPerturbation, 1, optional);
[a0bcf1]1662 if (!R->DoPerturbation) {
[961b75]1663 if (!ParseForParameter(P->Call.out[ReadOut],file, "BField", 0, 1, 3, grid, &R->BField, 1, optional)) { // old parameter for perturbation with field direction
[a0bcf1]1664 if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) No perturbation specified.\n", P->Par.me);
[961b75]1665 if (F->DoOutCurr || F->DoOutNICS) {
1666 if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) DoOutCurr/DoOutNICS = 1 though no DoPerturbation specified: setting DoOutCurr = 0\n", P->Par.me);
[a0bcf1]1667 F->DoOutCurr = 0;
[961b75]1668 F->DoOutNICS = 0;
[a0bcf1]1669 R->DoPerturbation = 0;
1670 }
1671 } else {
1672 R->DoPerturbation = 1;
1673 }
1674 }
1675 if (!P->Call.WriteSrcFiles && R->DoPerturbation) {
1676 if (P->Call.out[ReadOut]) fprintf(stderr,"(%i)Perturbation: Always writing source files (\"-w\").\n", P->Par.me);
1677 P->Call.WriteSrcFiles = 1;
1678 }
1679 RMatReci3(P->Lat.InvBasis, P->Lat.RealBasis);
1680
[961b75]1681 if (!ParseForParameter(P->Call.out[ReadOut],file,"DoFullCurrent", 0, 1, 1, int_type, &R->DoFullCurrent, 1, optional))
[a0bcf1]1682 R->DoFullCurrent = 0;
1683 if (R->DoFullCurrent < 0) R->DoFullCurrent = 0;
1684 if (R->DoFullCurrent > 2) R->DoFullCurrent = 2;
1685 if (R->DoPerturbation == 0) {
1686 if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) No Magnetic field: RunStruct#DoFullCurrent = 0\n", P->Par.me);
1687 R->DoFullCurrent = 0;
1688 }
1689
[961b75]1690 ParseForParameter(P->Call.out[ReadOut],file,"ECut", 0, 1, 1, double_type, &P->Lat.ECut, 1, critical);
[a0bcf1]1691 if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) ECut = %e -> %e Rydberg\n", me, P->Lat.ECut, 2.*P->Lat.ECut);
[961b75]1692 ParseForParameter(P->Call.out[ReadOut],file,"MaxLevel", 0, 1, 1, int_type, &P->Lat.MaxLevel, 1, critical);
1693 ParseForParameter(P->Call.out[ReadOut],file,"Level0Factor", 0, 1, 1, int_type, &P->Lat.Lev0Factor, 1, critical);
[a0bcf1]1694 if (P->Lat.Lev0Factor < 2) {
1695 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set Lev0Factor to 2!\n",me);
1696 P->Lat.Lev0Factor = 2;
1697 }
[961b75]1698 ParseForParameter(P->Call.out[ReadOut],file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
[a0bcf1]1699 if (di >= 0 && di < 2) {
1700 P->Lat.RT.Use = (enum UseRiemannTensor)di;
1701 } else {
1702 Error(SomeError, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
1703 }
1704 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!RiemannTensor = %i\n",me,P->Lat.RT.Use), critical;
1705 switch (P->Lat.RT.Use) {
1706 case UseNotRT:
1707 if (P->Lat.MaxLevel < 2) {
1708 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set MaxLevel to 2!\n",me);
1709 P->Lat.MaxLevel = 2;
1710 }
1711 P->Lat.LevRFactor = 2;
1712 P->Lat.RT.ActualUse = inactive;
1713 break;
1714 case UseRT:
1715 if (P->Lat.MaxLevel < 3) {
1716 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set MaxLevel to 3!\n",me);
1717 P->Lat.MaxLevel = 3;
1718 }
[961b75]1719 ParseForParameter(P->Call.out[ReadOut],file,"RiemannLevel", 0, 1, 1, int_type, &P->Lat.RT.RiemannLevel, 1, critical);
[a0bcf1]1720 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!RiemannLevel = %i\n",me,P->Lat.RT.RiemannLevel);
1721 if (P->Lat.RT.RiemannLevel < 2) {
1722 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set RiemannLevel to 2!\n",me);
1723 P->Lat.RT.RiemannLevel = 2;
1724 }
1725 if (P->Lat.RT.RiemannLevel > P->Lat.MaxLevel-1) {
1726 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set RiemannLevel to %i (MaxLevel-1)!\n",me,P->Lat.MaxLevel-1);
1727 P->Lat.RT.RiemannLevel = P->Lat.MaxLevel-1;
1728 }
[961b75]1729 ParseForParameter(P->Call.out[ReadOut],file,"LevRFactor", 0, 1, 1, int_type, &P->Lat.LevRFactor, 1, critical);
[a0bcf1]1730 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!LevRFactor = %i\n",me,P->Lat.LevRFactor);
1731 if (P->Lat.LevRFactor < 2) {
1732 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set LevRFactor to 2!\n",me);
1733 P->Lat.LevRFactor = 2;
1734 }
1735 P->Lat.Lev0Factor = 2;
1736 P->Lat.RT.ActualUse = standby;
1737 break;
1738 }
[961b75]1739 ParseForParameter(P->Call.out[ReadOut],file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
[a0bcf1]1740 if (di >= 0 && di < 2) {
1741 P->Lat.Psi.Use = (enum UseSpinType)di;
1742 } else {
1743 Error(SomeError, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
1744 }
1745 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiType = %i\n",me,P->Lat.Psi.Use);
1746 for (i=0; i<MaxPsiNoType; i++)
1747 P->Lat.Psi.GlobalNo[i] = 0;
1748 switch (P->Lat.Psi.Use) {
1749 case UseSpinDouble:
[961b75]1750 ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiDouble", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDouble], 1, critical);
1751 if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], 1, optional))
[a0bcf1]1752 P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
1753 if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
1754 R->DoUnOccupied = 1;
1755 else
1756 R->DoUnOccupied = 0;
1757 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!MaxPsiDouble = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoDouble],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
1758 //P->Lat.Psi.GlobalNo[PsiMaxNoDouble] += P->Lat.Psi.GlobalNo[PsiMaxAdd]; no more adding up on occupied ones
1759 P->Lat.Psi.GlobalNo[PsiMaxNo] = P->Lat.Psi.GlobalNo[PsiMaxNoDouble];
1760 break;
1761 case UseSpinUpDown:
1762 if (P->Par.proc[PEPsi] % 2) Error(SomeError,"UseSpinUpDown & P->Par.proc[PEGamma] % 2");
[961b75]1763 ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoUp", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoUp], 1, critical);
1764 ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoDown", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDown], 1, critical);
1765 if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], 1, optional))
[a0bcf1]1766 P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
1767 if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
1768 R->DoUnOccupied = 1;
1769 else
1770 R->DoUnOccupied = 0;
1771 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiMaxPsiUp = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoUp],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
1772 if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!PsiMaxPsiDown = %i + %i\n",me,P->Lat.Psi.GlobalNo[PsiMaxNoDown],P->Lat.Psi.GlobalNo[PsiMaxAdd]);
1773 //P->Lat.Psi.GlobalNo[PsiMaxNoUp] += P->Lat.Psi.GlobalNo[PsiMaxAdd];
1774 //P->Lat.Psi.GlobalNo[PsiMaxNoDown] += P->Lat.Psi.GlobalNo[PsiMaxAdd];
1775 if (abs(P->Lat.Psi.GlobalNo[PsiMaxNoUp]-P->Lat.Psi.GlobalNo[PsiMaxNoDown])>1) Error(SomeError,"|Up - Down| > 1");
1776 P->Lat.Psi.GlobalNo[PsiMaxNo] = P->Lat.Psi.GlobalNo[PsiMaxNoUp] + P->Lat.Psi.GlobalNo[PsiMaxNoDown];
1777 break;
1778 }
1779
1780 IonsInitRead(P, file);
1781 fclose(file);
1782}
1783
[774ae8]1784/** Writes parameters as a parsable config file to \a filename.
1785 * \param *P Problem at hand
1786 * \param filename name of config file to be written
1787 */
1788void WriteParameters(struct Problem *const P, const char *const filename) {
1789 struct Lattice *Lat = &P->Lat;
1790 struct RunStruct *R = &P->R;
1791 struct FileData *F = &P->Files;
1792 struct Ions *I = &P->Ion;
1793 FILE *output;
1794 double BorAng = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
1795 int i, it, nr = 0;
1796 if ((P->Par.me == 0) && (output = fopen(filename,"w"))) {
1797 //fprintf(stderr, "(%i) Writing config file %s\n",P->Par.me, filename);
1798
1799 fprintf(output,"# ParallelCarParinello - main configuration file - optimized geometry\n");
1800 fprintf(output,"\n");
1801 fprintf(output,"mainname\t%s\t# programm name (for runtime files)\n", F->mainname);
1802 fprintf(output,"defaultpath\t%s\t# where to put files during runtime\n", F->default_path);
1803 fprintf(output,"pseudopotpath\t%s\t# where to find pseudopotentials\n", F->pseudopot_path);
1804 fprintf(output,"\n");
1805 fprintf(output,"ProcPEGamma\t%i\t# for parallel computing: share constants\n", P->Par.proc[PEGamma]);
1806 fprintf(output,"ProcPEPsi\t%i\t# for parallel computing: share wave functions\n", P->Par.proc[PEPsi]);
1807 fprintf(output,"DoOutVis\t%i\t# Output data for OpenDX\n", F->DoOutVis);
1808 fprintf(output,"DoOutMes\t%i\t# Output data for measurements\n", F->DoOutMes);
1809 fprintf(output,"DoOutCurr\t%i\t# Ouput current density for OpenDx\n", F->DoOutCurr);
1810 fprintf(output,"DoPerturbation\t%i\t# Do perturbation calculate and determine susceptibility and shielding\n", R->DoPerturbation);
1811 fprintf(output,"DoFullCurrent\t%i\t# Do full perturbation\n", R->DoFullCurrent);
1812 fprintf(output,"DoOutOrbitals\t%i\t# Output all orbitals on end\n", F->DoOutOrbitals);
1813 fprintf(output,"CommonWannier\t%i\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center\n", R->CommonWannier);
1814 fprintf(output,"SawtoothStart\t%lg\t# Absolute value for smooth transition at cell border \n", Lat->SawtoothStart);
1815 fprintf(output,"VectorPlane\t%i\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot\n", R->VectorPlane);
1816 fprintf(output,"VectorCut\t%lg\t# Cut plane axis value\n", R->VectorCut);
1817 fprintf(output,"AddGramSch\t%i\t# Additional GramSchmidtOrtogonalization to be safe (1 - per MinStep, 2 - per PsiStep)\n", R->UseAddGramSch);
1818 fprintf(output,"Seed\t\t%i\t# initial value for random seed for Psi coefficients\n", R->Seed);
1819 fprintf(output,"DoBrent\t%i\t# Brent line search instead of CG with second taylor approximation\n", R->DoBrent);
1820 fprintf(output,"\n");
1821 fprintf(output,"MaxOuterStep\t%i\t# number of MolecularDynamics/Structure optimization steps\n", R->MaxOuterStep);
1822 fprintf(output,"Deltat\t%lg\t# time in atomic time per MD step\n", R->delta_t);
1823 fprintf(output,"OutVisStep\t%i\t# Output visual data every ...th step\n", R->OutVisStep);
1824 fprintf(output,"OutSrcStep\t%i\t# Output \"restart\" data every ..th step\n", R->OutSrcStep);
1825 fprintf(output,"TargetTemp\t%lg\t# Target temperature in Hartree\n", R->TargetTemp);
1826 fprintf(output,"Thermostat\t%s", I->ThermostatNames[I->Thermostat]);
1827 switch(I->Thermostat) {
1828 default:
1829 case 0: // None
1830 break;
1831 case 1: // Woodcock
1832 fprintf(output,"\t%i", R->ScaleTempStep);
1833 break;
1834 case 2: // Gaussian
1835 fprintf(output,"\t%i", R->ScaleTempStep);
1836 break;
1837 case 3: // Langevin
1838 fprintf(output,"\t%lg\t%lg", R->TempFrequency, R->alpha);
1839 break;
1840 case 4: // Berendsen
1841 fprintf(output,"\t%lg", R->TempFrequency);
1842 break;
1843 case 5: // NoseHoover
1844 fprintf(output,"\t%lg", R->HooverMass);
1845 break;
1846 }
1847 fprintf(output,"\t# Thermostat to be used with parameters\n");
1848 fprintf(output,"MaxPsiStep\t%i\t# number of Minimisation steps per state (0 - default)\n", R->MaxPsiStep);
1849 fprintf(output,"EpsWannier\t%lg\t# tolerance value for spread minimisation of orbitals\n", R->EpsWannier);
1850 fprintf(output,"\n");
1851 fprintf(output,"# Values specifying when to stop\n");
1852 fprintf(output,"MaxMinStep\t%i\t# Maximum number of steps\n", R->MaxMinStep);
1853 fprintf(output,"RelEpsTotalE\t%lg\t# relative change in total energy\n", R->RelEpsTotalEnergy);
1854 fprintf(output,"RelEpsKineticE\t%lg\t# relative change in kinetic energy\n", R->RelEpsKineticEnergy);
1855 fprintf(output,"MaxMinStopStep\t%i\t# check every ..th steps\n", R->MaxMinStopStep);
1856 fprintf(output,"MaxMinGapStopStep\t%i\t# check every ..th steps\n", R->MaxMinGapStopStep);
1857 fprintf(output,"\n");
1858 fprintf(output,"# Values specifying when to stop for INIT, otherwise same as above\n");
1859 fprintf(output,"MaxInitMinStep\t%i\t# Maximum number of steps\n", R->MaxInitMinStep);
1860 fprintf(output,"InitRelEpsTotalE\t%lg\t# relative change in total energy\n", R->InitRelEpsTotalEnergy);
1861 fprintf(output,"InitRelEpsKineticE\t%lg\t# relative change in kinetic energy\n", R->InitRelEpsKineticEnergy);
1862 fprintf(output,"InitMaxMinStopStep\t%i\t# check every ..th steps\n", R->InitMaxMinStopStep);
1863 fprintf(output,"InitMaxMinGapStopStep\t%i\t# check every ..th steps\n", R->InitMaxMinGapStopStep);
1864 fprintf(output,"\n");
1865 fprintf(output,"BoxLength\t\t\t# (Length of a unit cell)\n");
1866 fprintf(output,"%lg\t\n", Lat->RealBasis[0]*BorAng);
1867 fprintf(output,"%lg\t%lg\t\n", Lat->RealBasis[3]*BorAng, Lat->RealBasis[4]*BorAng);
1868 fprintf(output,"%lg\t%lg\t%lg\t\n", Lat->RealBasis[6]*BorAng, Lat->RealBasis[7]*BorAng, Lat->RealBasis[8]*BorAng);
1869 // FIXME
1870 fprintf(output,"\n");
1871 fprintf(output,"ECut\t\t%lg\t# energy cutoff for discretization in Hartrees\n", Lat->ECut/(Lat->LevelSizes[0]*Lat->LevelSizes[0]));
1872 fprintf(output,"MaxLevel\t%i\t# number of different levels in the code, >=2\n", P->Lat.MaxLevel);
1873 fprintf(output,"Level0Factor\t%i\t# factor by which node number increases from S to 0 level\n", P->Lat.Lev0Factor);
1874 fprintf(output,"RiemannTensor\t%i\t# (Use metric)\n", P->Lat.RT.Use);
1875 switch (P->Lat.RT.Use) {
1876 case 0: //UseNoRT
1877 break;
1878 case 1: // UseRT
1879 fprintf(output,"RiemannLevel\t%i\t# Number of Riemann Levels\n", P->Lat.RT.RiemannLevel);
1880 fprintf(output,"LevRFactor\t%i\t# factor by which node number increases from 0 to R level from\n", P->Lat.LevRFactor);
1881 break;
1882 }
1883 fprintf(output,"PsiType\t\t%i\t# 0 - doubly occupied, 1 - SpinUp,SpinDown\n", P->Lat.Psi.Use);
1884 // write out both types for easier changing afterwards
1885 // switch (PsiType) {
1886 // case 0:
1887 fprintf(output,"MaxPsiDouble\t%i\t# here: specifying both maximum number of SpinUp- and -Down-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoDouble]);
1888 // break;
1889 // case 1:
1890 fprintf(output,"PsiMaxNoUp\t%i\t# here: specifying maximum number of SpinUp-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoUp]);
1891 fprintf(output,"PsiMaxNoDown\t%i\t# here: specifying maximum number of SpinDown-states\n", P->Lat.Psi.GlobalNo[PsiMaxNoDown]);
1892 // break;
1893 // }
1894 fprintf(output,"AddPsis\t\t%i\t# Additional unoccupied Psis for bandgap determination\n", P->Lat.Psi.GlobalNo[PsiMaxAdd]);
1895 fprintf(output,"\n");
1896 fprintf(output,"RCut\t\t%lg\t# R-cut for the ewald summation\n", I->R_cut);
1897 fprintf(output,"StructOpt\t%i\t# Do structure optimization beforehand\n", I->StructOpt);
1898 fprintf(output,"IsAngstroem\t%i\t# 0 - Bohr, 1 - Angstroem\n", I->IsAngstroem);
1899 fprintf(output,"RelativeCoord\t0\t# whether ion coordinates are relative (1) or absolute (0)\n");
1900 fprintf(output,"MaxTypes\t%i\t# maximum number of different ion types\n", I->Max_Types);
1901 fprintf(output,"\n");
1902 // now elements ...
1903 fprintf(output,"# Ion type data (PP = PseudoPotential, Z = atomic number)\n");
1904 fprintf(output,"#Ion_TypeNr.\tAmount\tZ\tRGauss\tL_Max(PP)L_Loc(PP)IonMass\tName\tSymbol\n");
1905 for (i=0; i < I->Max_Types; i++)
1906 fprintf(output,"Ion_Type%i\t%i\t%i\t%lg\t%i\t%i\t%lg\t%s\t%s\n", i+1, I->I[i].Max_IonsOfType, I->I[i].Z, I->I[i].rgauss, I->I[i].l_max, I->I[i].l_loc, I->I[i].IonMass/Units2Electronmass, &I->I[i].Name[0], &I->I[i].Symbol[0]);
1907 // ... and each ion coordination
1908 fprintf(output,"#Ion_TypeNr._Nr.R[0]\tR[1]\tR[2]\tMoveType (0 MoveIon, 1 FixedIon) U[0]\tU[1]\tU[2] (velocities in %s)\n", (I->IsAngstroem) ? "Angstroem" : "atomic units");
1909 for (i=0; i < I->Max_Types; i++)
1910 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
1911 fprintf(output,"Ion_Type%i_%i\t%2.9f\t%2.9f\t%2.9f\t%i\t%2.9f\t%2.9f\t%2.9f\t# Number in molecule %i\n", i+1, it+1, I->I[i].R[0+NDIM*it]*BorAng, I->I[i].R[1+NDIM*it]*BorAng, I->I[i].R[2+NDIM*it]*BorAng, I->I[i].IMT[it], I->I[i].U[0+NDIM*it]*BorAng, I->I[i].U[1+NDIM*it]*BorAng, I->I[i].U[2+NDIM*it]*BorAng, ++nr);
1912 }
1913 fflush(output);
1914 fclose(output);
1915 }
1916}
1917
[a0bcf1]1918/** Main Initialization.
1919 * Massive calling of initializing super-functions init...():
1920 * -# InitOutputFiles()\n
1921 * Opening and Initializing of output measurement files.
1922 * -# InitLattice()\n
1923 * Initializing Lattice structure.
1924 * -# InitRunLevel()\n
1925 * Initializing RunLevel structure.
1926 * -# InitLatticeLevel0()\n
1927 * Initializing LatticeLevel structure.
1928 * -# InitOtherLevel()\n
1929 * Initialization of first to Lattice::MaxLevel LatticeLevel's.
1930 * -# InitRun()\n
1931 * Initialization of RunStruct structure.
1932 * -# InitPsis()\n
1933 * Initialization of wave functions as a whole.
1934 * -# InitDensity()\n
1935 * Initializing real and complex Density arrays.
1936 * -# either ReadSrcFiles() or InitPsisValues()\n
1937 * read old calculations from source file if desired or initialise with
1938 * random values InitPsisValue().
1939 * -# InitRiemannTensor()\n
1940 * Initializing RiemannTensor structure.
1941 * -# FirstInitGramSchData()\n
1942 * Sets up OnePsiElement as an MPI structure and subsequently initialises all of them.
1943 *
1944 * Now follow various (timed) tests:
1945 * -# GramSch()\n
1946 * (twice) to orthonormalize the random values
1947 * -# TestGramSch()\n
1948 * Test if orbitals now fulfill kronecker delta relation.
1949 * -# InitGradients()\n
1950 * Initializing Gradient structure, allocating its arrays.
1951 * -# InitDensityCalculation()\n
1952 * Calculate naive, real and complex densities and sum them up.
1953 * -# InitPseudoRead()\n
1954 * Read PseudoPot'entials from file.
1955 * -# InitEnergyArray()\n
1956 * Initialize Energy arrays.
1957 * -# InitPsiEnergyCalculation()\n
1958 * Calculate all psi energies and sum up.
1959 * -# CalculateDensityEnergy()\n
1960 * Calculate Density energy.
1961 * -# CalculateIonsEnergy()\n
1962 * Calculate ionic energy.
1963 * -# CalculateEwald()\n
1964 * Calculate ionic ewald summation.
1965 * -# EnergyAllReduce()\n
1966 * Gather energies from all processes and sum up.
1967 * -# EnergyOutput()\n
1968 * First output to file.
1969 * -# InitOutVisArray()\n
1970 * First output of visual for Opendx.
1971 *
1972 * \param *P Problem at hand
1973 */
1974void Init(struct Problem *const P) {
1975 struct RunStruct *R = &P->R;
1976 InitOutputFiles(P);
1977 InitLattice(P);
1978 InitRunLevel(P);
1979 InitLatticeLevel0(P);
1980 InitOtherLevel(P);
1981 InitPsiGroupNo(P);
1982 InitRun(P);
1983 InitPsis(P);
1984 InitDensity(P);
1985 if (P->Call.ReadSrcFiles) {
1986 if (ReadSrcIons(P) && P->Call.out[NormalOut])
1987 fprintf(stderr,"(%i) Ionic structure read successfully.\n", P->Par.me);
1988 else
1989 fprintf(stderr,"(%i) No readible Ionic structure found.\n", P->Par.me);
1990 } else {
1991 InitPsisValue(P, 0, P->Lat.Psi.TypeStartIndex[Perturbed_P0]);
1992 }
1993 InitRiemannTensor(P);
1994 if (P->Lat.RT.ActualUse == standby && R->LevRSNo == R->LevSNo) P->Lat.RT.ActualUse = active;
1995 CalculateRiemannTensorData(P);
1996 FirstInitGramSchData(P,&P->Lat.Psi);
1997 /* Test */
1998 SpeedMeasure(P, InitSimTime, StartTimeDo);
1999 SpeedMeasure(P, InitGramSchTime, StartTimeDo);
2000 GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
2001 ResetGramSch(P, &P->Lat.Psi);
2002 if (R->DoUnOccupied == 1) {
2003 ResetGramSchTagType(P, &P->Lat.Psi, UnOccupied, NotOrthogonal);
2004 GramSch(P, R->LevS, &P->Lat.Psi, Orthonormalize);
2005 }
2006 SpeedMeasure(P, InitGramSchTime, StopTimeDo);
2007 SpeedMeasure(P, InitSimTime, StopTimeDo);
2008 TestGramSch(P, R->LevS, &P->Lat.Psi, -1);
2009 ResetGramSchTagType(P, &P->Lat.Psi, UnOccupied, NotUsedToOrtho);
2010 InitGradients(P, &P->Grad);
2011 SpeedMeasure(P, InitSimTime, StartTimeDo);
2012 SpeedMeasure(P, InitDensityTime, StartTimeDo);
2013 InitDensityCalculation(P);
2014 SpeedMeasure(P, InitDensityTime, StopTimeDo);
2015 SpeedMeasure(P, InitSimTime, StopTimeDo);
2016 InitPseudoRead(P, P->Files.pseudopot_path);
2017 InitEnergyArray(P);
2018 SpeedMeasure(P, InitSimTime, StartTimeDo);
2019 InitPsiEnergyCalculation(P, Occupied); /* STANDARTLEVEL */
2020 CalculateDensityEnergy(P, 1); /* STANDARTLEVEL */
2021 CalculateIonsEnergy(P);
2022 CalculateEwald(P, 1);
2023 EnergyAllReduce(P);
2024/* SetCurrentMinState(P,UnOccupied);
2025 InitPsiEnergyCalculation(P,UnOccupied);
2026 CalculateGapEnergy(P);
2027 EnergyAllReduce(P);
2028 SetCurrentMinState(P,Occupied);*/
2029 SpeedMeasure(P, InitSimTime, StopTimeDo);
2030 R->LevS->Step++;
2031 EnergyOutput(P,0);
2032 InitOutVisArray(P);
[64fa9e]2033 Free(P->Lat.GArrayForSort, "Init: P->Lat.GArrayForSort");
[a0bcf1]2034 P->Speed.InitSteps++;
2035}
2036
Note: See TracBrowser for help on using the repository browser.