OpenAtom  Version1.5a
rhoSlab.C
Go to the documentation of this file.
1 /*****************************************************************************
2  * $Source$
3  * $Author$
4  * $Date$
5  * $Revision$
6  *****************************************************************************/
7 
8 ///////////////////////////////////////////////////////////////////////////////=
9 ///////////////////////////////////////////////////////////////////////////////c
10 ///////////////////////////////////////////////////////////////////////////////=
11 /** \file rhoSlab.C
12  * Add functions to allow application programmers to initialize these and
13  * the corresponding functions in Charm++ to call these functions with
14  * appropriate parameters
15  */
16 ///////////////////////////////////////////////////////////////////////////////=
17 
18 #include "utility/util.h"
19 #include "debug_flags.h"
20 #include <cmath>
21 #include "main/cpaimd.h"
24 
25 ///////////////////////////////////////////////////////////////////////////////=
26 
27 extern CkVec <CProxy_FFTcache> UfftCacheProxy;
28 extern Config config;
29 extern int nstates;
30 extern int sizeX;
31 
32 ///////////////////////////////////////////////////////////////////////////////=
33 ///////////////////////////////////////////////////////////////////////////////=
34 ///////////////////////////////////////////////////////////////////////////////c
35 ///////////////////////////////////////////////////////////////////////////////=
36 RhoRealSlab::~RhoRealSlab(){
37 
38  fftw_free(densityC);
39  fftw_free(VksC);
40  fftw_free(rhoIRXC);
41  fftw_free(rhoIRYC);
42  fftw_free(rhoIRZC);
43  fftw_free(VksHartC);
44 }
45 ///////////////////////////////////////////////////////////////////////////////=
46 
47 
48 ///////////////////////////////////////////////////////////////////////////////=
49 ///////////////////////////////////////////////////////////////////////////////c
50 ///////////////////////////////////////////////////////////////////////////////=
51 /* This gets called at the end of the RealSpaceDensity Constructor */
52 ///////////////////////////////////////////////////////////////////////////////=
53 void initRhoRealSlab(RhoRealSlab *rho_rs, int xdim, int ydim, int zdim,
54  int xdimA, int ydimA, int myIndexX, int myIndexY,
55  int rhoRsubplanes)
56 ///////////////////////////////////////////////////////////////////////////////=
57  {//begin routine
58 ///////////////////////////////////////////////////////////////////////////////=
59 
60  rho_rs->rhoRsubplanes = rhoRsubplanes;
61  rho_rs->sizeX = xdim;
62  rho_rs->sizeY = ydim;
63  rho_rs->sizeZ = zdim;
64 
65  rho_rs->size = (rho_rs->sizeX+2) * (rho_rs->sizeY);
66  rho_rs->trueSize = (rho_rs->sizeX) * (rho_rs->sizeY);
67 
68  int sizenow = rho_rs->size;
69  int csizenow = sizenow/2;
70 
71  complex *dummy;
72  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
73  rho_rs->VksC = dummy;
74  rho_rs->Vks = reinterpret_cast<double*> (dummy);
75 
76  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
77  rho_rs->densityC = dummy;
78  rho_rs->density = reinterpret_cast<double*> (dummy);
79 
80  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
81  rho_rs->rhoIRXC = dummy;
82  rho_rs->rhoIRX = reinterpret_cast<double*> (dummy);
83 
84  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
85  rho_rs->rhoIRYC = dummy;
86  rho_rs->rhoIRY = reinterpret_cast<double*> (dummy);
87 
88  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
89  rho_rs->rhoIRZC = dummy;
90  rho_rs->rhoIRZ = reinterpret_cast<double*> (dummy);
91 
92  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
93  rho_rs->VksHartC = dummy;
94  rho_rs->VksHart = reinterpret_cast<double*> (dummy);
95 
96  // if you have an extra transpose, you need a little more memory
97  // to receive messages asynchronously from other elements
98 
99  rho_rs->csizeInt = 0;
100  rho_rs->rsizeInt = 0;
101  if(rhoRsubplanes>1){
102 
103  csizenow = xdimA*ydimA;
104  rho_rs->csizeInt = csizenow;
105  rho_rs->rsizeInt = 2*csizenow;
106 
107  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
108  rho_rs->rhoIRXCint = dummy;
109  rho_rs->rhoIRXint = reinterpret_cast<double*> (dummy);
110 
111  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
112  rho_rs->rhoIRYCint = dummy;
113  rho_rs->rhoIRYint = reinterpret_cast<double*> (dummy);
114 
115  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
116  rho_rs->rhoIRZCint = dummy;
117  rho_rs->rhoIRZint = reinterpret_cast<double*> (dummy);
118 
119  dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
120  rho_rs->VksHartCint = dummy;
121  rho_rs->VksHartint = reinterpret_cast<double*> (dummy);
122 
123  }//endif
124 
125 ///////////////////////////////////////////////////////////////////////////////=
126  }//end routine
127 ///////////////////////////////////////////////////////////////////////////////=
128 
129 
130 ///////////////////////////////////////////////////////////////////////////////=
131 ///////////////////////////////////////////////////////////////////////////////c
132 ///////////////////////////////////////////////////////////////////////////////=
133 RhoGSlab::~RhoGSlab(){
134 
135  if(Rho !=NULL)fftw_free(Rho);
136  if(divRhoX !=NULL)fftw_free(divRhoX);
137  if(divRhoY !=NULL)fftw_free(divRhoY);
138  if(divRhoZ !=NULL)fftw_free(divRhoZ);
139  if(packedRho!=NULL)fftw_free(packedRho);
140  if(packedVks!=NULL)fftw_free(packedVks);
141 
142  if(runs!=NULL) delete [] runs;
143  if(k_x !=NULL) fftw_free(k_x);
144  if(k_y !=NULL) fftw_free(k_y);
145  if(k_z !=NULL) fftw_free(k_z);
146  if(perdCorr!=NULL)fftw_free(perdCorr);
147 
148 }
149 ///////////////////////////////////////////////////////////////////////////////=
150 
151 ///////////////////////////////////////////////////////////////////////////////=
152 ///////////////////////////////////////////////////////////////////////////////c
153 ///////////////////////////////////////////////////////////////////////////////=
154 void RhoGSlab::pup(PUP::er &p) {
155  // local bools so we only bother with non null objects
156  p|iperd;
157  p|sizeX;
158  p|sizeY;
159  p|sizeZ;
160  p|runsToBeSent;
161  p|numRuns;
162  p|numLines;
163  p|numFull;
164  p|numPoints;
165  p|ehart_ret;
166  p|eext_ret;
167  p|ewd_ret;
168  p|size;
169  p|nPacked;
170 
171  bool RhoMake =false;
172  bool divRhoXMake =false;
173  bool divRhoYMake =false;
174  bool divRhoZMake =false;
175  bool packedRhoMake=false;
176  bool packedVksMake=false;
177  bool VksMake =false;
178  if(!p.isUnpacking()){// create flags for each array in pup
179  RhoMake = (Rho !=NULL) ? true :false;
180  divRhoXMake = (divRhoX !=NULL) ? true :false;
181  divRhoYMake = (divRhoY !=NULL) ? true :false;
182  divRhoZMake = (divRhoZ !=NULL) ? true :false;
183  packedRhoMake = (packedRho!=NULL) ? true :false;
184  packedVksMake = (packedVks!=NULL) ? true :false;
185  VksMake = (Vks !=NULL) ? true :false;
186  }//endif
187 
188  p|RhoMake;
189  p|divRhoXMake;
190  p|divRhoYMake;
191  p|divRhoZMake;
192  p|packedRhoMake;
193  p|packedVksMake;
194  p|VksMake;
195 
196  if(p.isUnpacking()){
197  if(RhoMake)
198  Rho = (complex *)fftw_malloc(numFull*sizeof(complex));
199  else
200  Rho = NULL;
201  if(divRhoXMake)
202  divRhoX = (complex *)fftw_malloc(numFull*sizeof(complex));
203  else
204  divRhoX = NULL;
205  if(divRhoYMake)
206  divRhoY = (complex *)fftw_malloc(numFull*sizeof(complex));
207  else
208  divRhoY = NULL;
209  if(divRhoZMake)
210  divRhoZ = (complex *)fftw_malloc(numFull*sizeof(complex));
211  else
212  divRhoZ = NULL;
213  if(packedRhoMake)
214  packedRho = (complex *)fftw_malloc(nPacked*sizeof(complex));
215  else
216  packedRho = NULL;
217  if(packedVksMake)
218  packedVks = (complex *)fftw_malloc(nPacked*sizeof(complex));
219  else
220  packedVks = NULL;
221  if(VksMake)
222  Vks = (complex *)fftw_malloc(numFull*sizeof(complex));
223  else
224  Vks = NULL;
225  runs = new RunDescriptor[numRuns];
226  k_x = (int *)fftw_malloc(numPoints*sizeof(int));
227  k_y = (int *)fftw_malloc(numPoints*sizeof(int));
228  k_z = (int *)fftw_malloc(numPoints*sizeof(int));
229  perdCorr = NULL;
230  if(iperd!=3){perdCorr = (double *)fftw_malloc(numPoints*sizeof(double));}
231  }//endif : unpacking malloc
232 
233  if(RhoMake) PUParray(p,Rho,numFull);
234  if(divRhoXMake) PUParray(p,divRhoX,numFull);
235  if(divRhoYMake) PUParray(p,divRhoY,numFull);
236  if(divRhoZMake) PUParray(p,divRhoZ,numFull);
237  if(packedRhoMake) PUParray(p,packedRho,nPacked);
238  if(packedVksMake) PUParray(p,packedVks,nPacked);
239  if(VksMake) PUParray(p,Vks,numFull);
240 
241  PUParray(p,runs,numRuns);
242  PUParray(p,k_x,numPoints);
243  PUParray(p,k_y,numPoints);
244  PUParray(p,k_z,numPoints);
245  if(iperd!=3){PUParray(p,perdCorr,numPoints);}
246 
247 //------------------------------------------------------------------------------
248  }//end intense pupping experience
249 ///////////////////////////////////////////////////////////////////////////////=
250 
251 
252 ///////////////////////////////////////////////////////////////////////////////=
253 /// packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz
254 ///////////////////////////////////////////////////////////////////////////////=
255 ///////////////////////////////////////////////////////////////////////////////c
256 ///////////////////////////////////////////////////////////////////////////////=
257 
258 void RhoGSlab::divRhoGdot(double *hmati, double tpi,complex *tmpRho){
259 
260 ///////////////////////////////////////////////////////////////////////////////=
261 
262  int nfftz = sizeZ;
263 
264 ///////////////////////////////////////////////////////////////////////////////=
265 ///
266 
267  bzero(divRhoX,sizeof(complex)*numFull);
268  bzero(divRhoY,sizeof(complex)*numFull);
269  bzero(divRhoZ,sizeof(complex)*numFull);
270  double gx,gy,gz;
271 
272  int koff = 0;
273  for (int r = 0,l=0; r < numRuns; r+=2,l++) {
274 
275  int joff = l*nfftz + runs[r].z;
276  for (int i=0,j=joff,k=koff; i<runs[r].length; i++,j++,k++) {
277  gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
278  gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
279  gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
280  complex tmp = (tmpRho[k].multiplyByi())*(-1.0);
281  divRhoX[j] = tmp*gx;
282  divRhoY[j] = tmp*gy;
283  divRhoZ[j] = tmp*gz;
284  }//endfor
285  koff += runs[r].length;
286 
287  int r1=r+1;
288  joff = l*nfftz + runs[r1].z;
289  for (int i=0,j=joff,k=koff; i<runs[r1].length; i++,j++,k++) {
290  gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
291  gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
292  gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
293  complex tmp = (tmpRho[k].multiplyByi())*(-1.0);
294  divRhoX[j] = tmp*gx;
295  divRhoY[j] = tmp*gy;
296  divRhoZ[j] = tmp*gz;
297  }//endfor
298  koff += runs[r1].length;
299 
300  }//endfor
301 
302  CkAssert(numPoints == koff);
303 
304 //------------------------------------------------------------------------------
305  }//end routine
306 ///////////////////////////////////////////////////////////////////////////////=
307 
308 
309 ///////////////////////////////////////////////////////////////////////////////=
310 /// packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz
311 ///////////////////////////////////////////////////////////////////////////////=
312 ///////////////////////////////////////////////////////////////////////////////c
313 ///////////////////////////////////////////////////////////////////////////////=
314 
315 void RhoGSlab::createWhiteByrd(double *hmati, double tpi){
316 
317 ///////////////////////////////////////////////////////////////////////////////=
318 
319  int nfftz = sizeZ;
320  double gx,gy,gz;
321  complex tmp;
322 
323 ///////////////////////////////////////////////////////////////////////////////=
324 
325  complex *whitebyrd = divRhoX; // zeroing done carefully inside loop
326  // so that we can save memory by reusing divRhoX
327  int koff = 0;
328  for (int r = 0,l=0; r < numRuns; r+=2,l++) {
329 
330  int joff1 = l*nfftz + runs[r].z;
331  for (int i=0,j=joff1,k=koff; i<runs[r].length; i++,j++,k++) {
332  gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
333  gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
334  gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
335  tmp = divRhoX[j]*gx + divRhoY[j]*gy + divRhoZ[j]*gz;
336  whitebyrd[j] = tmp.multiplyByi()*(-1.0);
337  }//endfor
338  koff += runs[r].length;
339 
340  int r1=r+1;
341  int joff2 = l*nfftz + runs[r1].z;
342  for (int i=0,j=joff2,k=koff; i<runs[r1].length; i++,j++,k++) {
343  gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
344  gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
345  gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
346  tmp = divRhoX[j]*gx + divRhoY[j]*gy + divRhoZ[j]*gz;
347  whitebyrd[j] = tmp.multiplyByi()*(-1.0);
348  }//endfor
349  koff += runs[r1].length;
350 
351  int joff3 = joff2+runs[r1].length;
352  for(int j=joff3;j<joff1;j++){whitebyrd[j]=0.0;}
353 
354  }//endfor
355 
356  CkAssert(numPoints == koff);
357 //------------------------------------------------------------------------------
358  }//end routine
359 ///////////////////////////////////////////////////////////////////////////////=
360 
361 
362 ///////////////////////////////////////////////////////////////////////////////=
363 ///////////////////////////////////////////////////////////////////////////////c
364 ///////////////////////////////////////////////////////////////////////////////=
365 /*
366  * In G-space our representation uses wrapping of coordinates.
367  * When reading from the files, if x < 0, x += sizeX;
368  */
369 ///////////////////////////////////////////////////////////////////////////////=
370 
372 
373 ///////////////////////////////////////////////////////////////////////=
374 /// Construct the k-vectors
375 
376  k_x = (int *)fftw_malloc(numPoints*sizeof(int));
377  k_y = (int *)fftw_malloc(numPoints*sizeof(int));
378  k_z = (int *)fftw_malloc(numPoints*sizeof(int));
379 
380  int r, i, dataCovered = 0;
381  int x, y, z;
382  for (r = 0; r < numRuns; r++) { // 2*number of lines z
383  x = runs[r].x;
384  if (x > sizeX/2) x -= sizeX;
385  y = runs[r].y;
386  if (y > sizeY/2) y -= sizeY;
387  z = runs[r].z;
388  if (z > sizeZ/2) z -= sizeZ;
389 
390  for (i = 0; i < runs[r].length; i++) { //pts in lines of z
391  k_x[dataCovered] = x;
392  k_y[dataCovered] = y;
393  k_z[dataCovered] = (z+i);
394  dataCovered++;
395  }//endfor
396  }//endfor
397 
398  CkAssert(dataCovered == numPoints);
399 
400 ///////////////////////////////////////////////////////////////////////////////=
401 /// Set the return values
402 
403  *n = numPoints;
404 
405 ///////////////////////////////////////////////////////////////////////////////=
406  }//end routine
407 ///////////////////////////////////////////////////////////////////////////////=
408 
409 ///////////////////////////////////////////////////////////////////////////////=
410 ///////////////////////////////////////////////////////////////////////////////c
411 ///////////////////////////////////////////////////////////////////////////////=
412 void RhoRealSlab::pup(PUP::er &p) {
413 
414  p|sizeZ; // fftX size
415  p|sizeY; // fftYZ size
416  p|sizeZ; // fftZ size
417  p|size; // plane size for fftw : bigger
418  p|trueSize;// the real plane size
419  p|exc_ret; // energies
420  p|muxc_ret;
421  p|exc_gga_ret;
422  p|rhoRsubplanes;
423  p|csizeInt;
424  p|rsizeInt;
425 
426  if(p.isUnpacking()){
427  int csize = size/2;
428  VksC = (complex*) fftw_malloc(csize*sizeof(complex));
429  Vks = reinterpret_cast<double*> (VksC);
430  densityC = (complex*) fftw_malloc(csize*sizeof(complex));
431  density = reinterpret_cast<double*> (densityC);
432  rhoIRXC = (complex*) fftw_malloc(csize*sizeof(complex));
433  rhoIRX = reinterpret_cast<double*> (rhoIRXC);
434  rhoIRYC = (complex*) fftw_malloc(csize*sizeof(complex));
435  rhoIRY = reinterpret_cast<double*> (rhoIRYC);
436  rhoIRZC = (complex*) fftw_malloc(csize*sizeof(complex));
437  rhoIRZ = reinterpret_cast<double*> (rhoIRZC);
438  VksHartC = (complex*) fftw_malloc(csize*sizeof(complex));
439  VksHart = reinterpret_cast<double*> (VksHartC);
440  if(rhoRsubplanes>1){
441  rhoIRXCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
442  rhoIRXint = reinterpret_cast<double*> (rhoIRXCint);
443  rhoIRYCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
444  rhoIRYint = reinterpret_cast<double*> (rhoIRYCint);
445  rhoIRZCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
446  rhoIRZint = reinterpret_cast<double*> (rhoIRZCint);
447  VksHartCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
448  VksHartint = reinterpret_cast<double*> (VksHartCint);
449  }
450  }
451  PUParray(p,Vks, size);
452  PUParray(p,density,size);
453  PUParray(p,rhoIRX, size);
454  PUParray(p,rhoIRY, size);
455  PUParray(p,rhoIRZ, size);
456  PUParray(p,VksHart,size);
457 
458  if(rhoRsubplanes==1){
459  PUParray(p,rhoIRXint, rsizeInt);
460  PUParray(p,rhoIRYint, rsizeInt);
461  PUParray(p,rhoIRZint, rsizeInt);
462  PUParray(p,VksHartint,rsizeInt);
463  }//endif
464 
465 }
466 ///////////////////////////////////////////////////////////////////////////////=
467 
468 
469 //////////////////////////////////////////////////////////////////////////////
470 //////////////////////////////////////////////////////////////////////////////
471 //////////////////////////////////////////////////////////////////////////////
472 void RhoRealSlab::uPackScale(double *uPackData, double *PackData,double scale){
473 
474  for(int i=0;i<size;i++){uPackData[i] = PackData[i]*scale;}
475 
476 }//end routine
477 //////////////////////////////////////////////////////////////////////////////
478 
479 
480 //////////////////////////////////////////////////////////////////////////////
481 //////////////////////////////////////////////////////////////////////////////
482 //////////////////////////////////////////////////////////////////////////////
483 void RhoRealSlab::scale(double *uPackData,double scale){
484 
485  for(int i=0;i<size;i++){uPackData[i] *= scale;}
486 
487 }//end routine
488 //////////////////////////////////////////////////////////////////////////////
489 
490 //////////////////////////////////////////////////////////////////////////////
491 //////////////////////////////////////////////////////////////////////////////
492 //////////////////////////////////////////////////////////////////////////////
493 void RhoRealSlab::uPackScaleShrink(double *uPackData,double *PackData,double scale){
494 
495  for(int i=0,i2=0;i<sizeY;i++,i2+=2){
496  for(int j=i*sizeX;j<(i+1)*sizeX;j++){
497  uPackData[j] =PackData[(j+i2)]*scale;
498  }//endfor
499  }//endfor
500 
501 }//end routine
502 //////////////////////////////////////////////////////////////////////////////
503 
504 
505 //////////////////////////////////////////////////////////////////////////////
506 //////////////////////////////////////////////////////////////////////////////
507 //////////////////////////////////////////////////////////////////////////////
508 void RhoRealSlab::uPackShrink(double *uPackData,double *PackData){
509 
510  for(int i=0,i2=0;i<sizeY;i++,i2+=2){
511  for(int j=i*sizeX;j<(i+1)*sizeX;j++){
512  uPackData[j] =PackData[(j+i2)];
513  }//endfor
514  }//endfor
515 
516 }//end routine
517 //////////////////////////////////////////////////////////////////////////////
518 
519 //////////////////////////////////////////////////////////////////////////////
520 //////////////////////////////////////////////////////////////////////////////
521 //////////////////////////////////////////////////////////////////////////////
522 void RhoRealSlab::uPackScaleGrow(double *uPackData,double *PackData,double scale){
523  for(int i=0,i2=0;i<sizeY;i++,i2+=2){
524  for(int j=i*sizeX;j<(i+1)*sizeX;j++){
525  uPackData[(j+i2)] =PackData[j]*scale;
526  }//endfor
527  int k = (i+1)*sizeX+i2;
528  uPackData[k] =0.0;
529  uPackData[(k+1)]=0.0;
530  }//endfor
531 
532 }//end routine
533 //////////////////////////////////////////////////////////////////////////////
void setKVectors(int *n)
Definition: rhoSlab.C:371
Add type declarations for simulationConstants class (readonly vars) and once class for each type of o...
void createWhiteByrd(double *, double)
= packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz
Definition: rhoSlab.C:315
== Index logic for lines of constant x,y in gspace.
Definition: RunDescriptor.h:22
Some basic data structures and the array map classes are defined here.
void divRhoGdot(double *, double, complex *)
= packed g-space of size numPoints is expanded to numFull =numRuns/2*nfftz
Definition: rhoSlab.C:258
Useful debugging flags.