source: pcp/src/init.c@ 90c027

Last change on this file since 90c027 was f5586e, checked in by Frederik Heber <heber@…>, 18 years ago

RealBasisQ[] removed

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