00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "util.h"
00019 #include "../../include/debug_flags.h"
00020 #include <math.h>
00021 #include "cpaimd.h"
00022 #include "groups.h"
00023 #include "fftCacheSlab.h"
00024 #include "CP_State_Plane.h"
00025
00026
00027
00028 extern CProxy_FFTcache fftCacheProxy;
00029 extern Config config;
00030 extern int nstates;
00031 extern int sizeX;
00032 extern CProxy_AtomsGrp atomsGrpProxy;
00033 extern CProxy_CPcharmParaInfoGrp scProxy;
00034
00035
00036
00037
00038
00039 RhoRealSlab::~RhoRealSlab(){
00040
00041 fftw_free(densityC);
00042 fftw_free(VksC);
00043 fftw_free(rhoIRXC);
00044 fftw_free(rhoIRYC);
00045 fftw_free(rhoIRZC);
00046 fftw_free(VksHartC);
00047 }
00048
00049
00050
00051
00052
00053
00054
00055
00056 void initRhoRealSlab(RhoRealSlab *rho_rs, int xdim, int ydim, int zdim,
00057 int xdimA, int ydimA, int myIndexX, int myIndexY,
00058 int rhoRsubplanes)
00059
00060 {
00061
00062
00063 rho_rs->rhoRsubplanes = rhoRsubplanes;
00064 rho_rs->sizeX = xdim;
00065 rho_rs->sizeY = ydim;
00066 rho_rs->sizeZ = zdim;
00067
00068 rho_rs->size = (rho_rs->sizeX+2) * (rho_rs->sizeY);
00069 rho_rs->trueSize = (rho_rs->sizeX) * (rho_rs->sizeY);
00070
00071 int sizenow = rho_rs->size;
00072 int csizenow = sizenow/2;
00073
00074 complex *dummy;
00075 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00076 rho_rs->VksC = dummy;
00077 rho_rs->Vks = reinterpret_cast<double*> (dummy);
00078
00079 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00080 rho_rs->densityC = dummy;
00081 rho_rs->density = reinterpret_cast<double*> (dummy);
00082
00083 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00084 rho_rs->rhoIRXC = dummy;
00085 rho_rs->rhoIRX = reinterpret_cast<double*> (dummy);
00086
00087 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00088 rho_rs->rhoIRYC = dummy;
00089 rho_rs->rhoIRY = reinterpret_cast<double*> (dummy);
00090
00091 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00092 rho_rs->rhoIRZC = dummy;
00093 rho_rs->rhoIRZ = reinterpret_cast<double*> (dummy);
00094
00095 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00096 rho_rs->VksHartC = dummy;
00097 rho_rs->VksHart = reinterpret_cast<double*> (dummy);
00098
00099
00100
00101
00102 rho_rs->csizeInt = 0;
00103 rho_rs->rsizeInt = 0;
00104 if(rhoRsubplanes>1){
00105
00106 csizenow = xdimA*ydimA;
00107 rho_rs->csizeInt = csizenow;
00108 rho_rs->rsizeInt = 2*csizenow;
00109
00110 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00111 rho_rs->rhoIRXCint = dummy;
00112 rho_rs->rhoIRXint = reinterpret_cast<double*> (dummy);
00113
00114 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00115 rho_rs->rhoIRYCint = dummy;
00116 rho_rs->rhoIRYint = reinterpret_cast<double*> (dummy);
00117
00118 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00119 rho_rs->rhoIRZCint = dummy;
00120 rho_rs->rhoIRZint = reinterpret_cast<double*> (dummy);
00121
00122 dummy = (complex*) fftw_malloc(csizenow*sizeof(complex));
00123 rho_rs->VksHartCint = dummy;
00124 rho_rs->VksHartint = reinterpret_cast<double*> (dummy);
00125
00126 }
00127
00128
00129 }
00130
00131
00132
00133
00134
00135
00136 RhoGSlab::~RhoGSlab(){
00137
00138 if(Rho !=NULL)fftw_free(Rho);
00139 if(divRhoX !=NULL)fftw_free(divRhoX);
00140 if(divRhoY !=NULL)fftw_free(divRhoY);
00141 if(divRhoZ !=NULL)fftw_free(divRhoZ);
00142 if(packedRho!=NULL)fftw_free(packedRho);
00143 if(packedVks!=NULL)fftw_free(packedVks);
00144
00145 if(runs!=NULL) delete [] runs;
00146 if(k_x !=NULL) fftw_free(k_x);
00147 if(k_y !=NULL) fftw_free(k_y);
00148 if(k_z !=NULL) fftw_free(k_z);
00149
00150 }
00151
00152
00153
00154
00155
00156 void RhoGSlab::pup(PUP::er &p) {
00157
00158 p|sizeX;
00159 p|sizeY;
00160 p|sizeZ;
00161 p|runsToBeSent;
00162 p|numRuns;
00163 p|numLines;
00164 p|numFull;
00165 p|numPoints;
00166 p|ehart_ret;
00167 p|eext_ret;
00168 p|ewd_ret;
00169 p|size;
00170 p|nPacked;
00171
00172 bool RhoMake =false;
00173 bool divRhoXMake =false;
00174 bool divRhoYMake =false;
00175 bool divRhoZMake =false;
00176 bool packedRhoMake=false;
00177 bool packedVksMake=false;
00178 bool VksMake =false;
00179 if(!p.isUnpacking()){
00180 RhoMake = (Rho !=NULL) ? true :false;
00181 divRhoXMake = (divRhoX !=NULL) ? true :false;
00182 divRhoYMake = (divRhoY !=NULL) ? true :false;
00183 divRhoZMake = (divRhoZ !=NULL) ? true :false;
00184 packedRhoMake = (packedRho!=NULL) ? true :false;
00185 packedVksMake = (packedVks!=NULL) ? true :false;
00186 VksMake = (Vks !=NULL) ? true :false;
00187 }
00188
00189 p|RhoMake;
00190 p|divRhoXMake;
00191 p|divRhoYMake;
00192 p|divRhoZMake;
00193 p|packedRhoMake;
00194 p|packedVksMake;
00195 p|VksMake;
00196
00197 if(p.isUnpacking()){
00198 if(RhoMake)
00199 Rho = (complex *)fftw_malloc(numFull*sizeof(complex));
00200 else
00201 Rho = NULL;
00202 if(divRhoXMake)
00203 divRhoX = (complex *)fftw_malloc(numFull*sizeof(complex));
00204 else
00205 divRhoX = NULL;
00206 if(divRhoYMake)
00207 divRhoY = (complex *)fftw_malloc(numFull*sizeof(complex));
00208 else
00209 divRhoY = NULL;
00210 if(divRhoZMake)
00211 divRhoZ = (complex *)fftw_malloc(numFull*sizeof(complex));
00212 else
00213 divRhoZ = NULL;
00214 if(packedRhoMake)
00215 packedRho = (complex *)fftw_malloc(nPacked*sizeof(complex));
00216 else
00217 packedRho = NULL;
00218 if(packedVksMake)
00219 packedVks = (complex *)fftw_malloc(nPacked*sizeof(complex));
00220 else
00221 packedVks = NULL;
00222 if(VksMake)
00223 Vks = (complex *)fftw_malloc(numFull*sizeof(complex));
00224 else
00225 Vks = NULL;
00226 runs = new RunDescriptor[numRuns];
00227 k_x = (int *)fftw_malloc(numPoints*sizeof(int));
00228 k_y = (int *)fftw_malloc(numPoints*sizeof(int));
00229 k_z = (int *)fftw_malloc(numPoints*sizeof(int));
00230
00231 }
00232
00233 if(RhoMake) PUParray(p,Rho,numFull);
00234 if(divRhoXMake) PUParray(p,divRhoX,numFull);
00235 if(divRhoYMake) PUParray(p,divRhoY,numFull);
00236 if(divRhoZMake) PUParray(p,divRhoZ,numFull);
00237 if(packedRhoMake) PUParray(p,packedRho,nPacked);
00238 if(packedVksMake) PUParray(p,packedVks,nPacked);
00239 if(VksMake) PUParray(p,Vks,numFull);
00240
00241 PUParray(p,runs,numRuns);
00242 PUParray(p,k_x,numPoints);
00243 PUParray(p,k_y,numPoints);
00244 PUParray(p,k_z,numPoints);
00245
00246
00247 }
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 void RhoGSlab::divRhoGdot(double *hmati, double tpi,complex *tmpRho){
00258
00259
00260
00261 int nfftz = sizeZ;
00262
00263
00264
00265
00266 bzero(divRhoX,sizeof(complex)*numFull);
00267 bzero(divRhoY,sizeof(complex)*numFull);
00268 bzero(divRhoZ,sizeof(complex)*numFull);
00269 double gx,gy,gz;
00270
00271 int koff = 0;
00272 for (int r = 0,l=0; r < numRuns; r+=2,l++) {
00273
00274 int joff = l*nfftz + runs[r].z;
00275 for (int i=0,j=joff,k=koff; i<runs[r].length; i++,j++,k++) {
00276 gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
00277 gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
00278 gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
00279 complex tmp = (tmpRho[k].multiplyByi())*(-1.0);
00280 divRhoX[j] = tmp*gx;
00281 divRhoY[j] = tmp*gy;
00282 divRhoZ[j] = tmp*gz;
00283 }
00284 koff += runs[r].length;
00285
00286 int r1=r+1;
00287 joff = l*nfftz + runs[r1].z;
00288 for (int i=0,j=joff,k=koff; i<runs[r1].length; i++,j++,k++) {
00289 gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
00290 gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
00291 gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
00292 complex tmp = (tmpRho[k].multiplyByi())*(-1.0);
00293 divRhoX[j] = tmp*gx;
00294 divRhoY[j] = tmp*gy;
00295 divRhoZ[j] = tmp*gz;
00296 }
00297 koff += runs[r1].length;
00298
00299 #ifdef CMK_VERSION_BLUEGENE
00300 if(r % 40==0){CmiNetworkProgress();}
00301 #endif
00302
00303 }
00304
00305 CkAssert(numPoints == koff);
00306
00307 #ifdef CMK_VERSION_BLUEGENE
00308 CmiNetworkProgress();
00309 #endif
00310
00311
00312 }
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 void RhoGSlab::createWhiteByrd(double *hmati, double tpi){
00323
00324
00325
00326 int nfftz = sizeZ;
00327 double gx,gy,gz;
00328 complex tmp;
00329
00330
00331
00332 complex *whitebyrd = divRhoX;
00333
00334 int koff = 0;
00335 for (int r = 0,l=0; r < numRuns; r+=2,l++) {
00336
00337 int joff1 = l*nfftz + runs[r].z;
00338 for (int i=0,j=joff1,k=koff; i<runs[r].length; i++,j++,k++) {
00339 gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
00340 gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
00341 gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
00342 tmp = divRhoX[j]*gx + divRhoY[j]*gy + divRhoZ[j]*gz;
00343 whitebyrd[j] = tmp.multiplyByi()*(-1.0);
00344 }
00345 koff += runs[r].length;
00346
00347 int r1=r+1;
00348 int joff2 = l*nfftz + runs[r1].z;
00349 for (int i=0,j=joff2,k=koff; i<runs[r1].length; i++,j++,k++) {
00350 gx = tpi*(k_x[k]*hmati[1] + k_y[k]*hmati[2] + k_z[k]*hmati[3]);
00351 gy = tpi*(k_x[k]*hmati[4] + k_y[k]*hmati[5] + k_z[k]*hmati[6]);
00352 gz = tpi*(k_x[k]*hmati[7] + k_y[k]*hmati[8] + k_z[k]*hmati[9]);
00353 tmp = divRhoX[j]*gx + divRhoY[j]*gy + divRhoZ[j]*gz;
00354 whitebyrd[j] = tmp.multiplyByi()*(-1.0);
00355 }
00356 koff += runs[r1].length;
00357
00358 int joff3 = joff2+runs[r1].length;
00359 for(int j=joff3;j<joff1;j++){whitebyrd[j]=0.0;}
00360
00361 #ifdef CMK_VERSION_BLUEGENE
00362 if(r % 40==0){CmiNetworkProgress();}
00363 #endif
00364
00365 }
00366
00367 CkAssert(numPoints == koff);
00368
00369 #ifdef CMK_VERSION_BLUEGENE
00370 CmiNetworkProgress();
00371 #endif
00372
00373
00374
00375 }
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 void RhoGSlab::setKVectors(int *n){
00389
00390
00391
00392
00393 k_x = (int *)fftw_malloc(numPoints*sizeof(int));
00394 k_y = (int *)fftw_malloc(numPoints*sizeof(int));
00395 k_z = (int *)fftw_malloc(numPoints*sizeof(int));
00396
00397 int r, i, dataCovered = 0;
00398 int x, y, z;
00399 for (r = 0; r < numRuns; r++) {
00400 x = runs[r].x;
00401 if (x > sizeX/2) x -= sizeX;
00402 y = runs[r].y;
00403 if (y > sizeY/2) y -= sizeY;
00404 z = runs[r].z;
00405 if (z > sizeZ/2) z -= sizeZ;
00406
00407 for (i = 0; i < runs[r].length; i++) {
00408 k_x[dataCovered] = x;
00409 k_y[dataCovered] = y;
00410 k_z[dataCovered] = (z+i);
00411 dataCovered++;
00412 }
00413 }
00414
00415 CkAssert(dataCovered == numPoints);
00416
00417
00418
00419
00420 *n = numPoints;
00421
00422
00423 }
00424
00425
00426
00427
00428
00429 void RhoRealSlab::pup(PUP::er &p) {
00430
00431 p|sizeZ;
00432 p|sizeY;
00433 p|sizeZ;
00434 p|size;
00435 p|trueSize;
00436 p|exc_ret;
00437 p|muxc_ret;
00438 p|exc_gga_ret;
00439 p|rhoRsubplanes;
00440 p|csizeInt;
00441 p|rsizeInt;
00442
00443 if(p.isUnpacking()){
00444 int csize = size/2;
00445 VksC = (complex*) fftw_malloc(csize*sizeof(complex));
00446 Vks = reinterpret_cast<double*> (VksC);
00447 densityC = (complex*) fftw_malloc(csize*sizeof(complex));
00448 density = reinterpret_cast<double*> (densityC);
00449 rhoIRXC = (complex*) fftw_malloc(csize*sizeof(complex));
00450 rhoIRX = reinterpret_cast<double*> (rhoIRXC);
00451 rhoIRYC = (complex*) fftw_malloc(csize*sizeof(complex));
00452 rhoIRY = reinterpret_cast<double*> (rhoIRYC);
00453 rhoIRZC = (complex*) fftw_malloc(csize*sizeof(complex));
00454 rhoIRZ = reinterpret_cast<double*> (rhoIRZC);
00455 VksHartC = (complex*) fftw_malloc(csize*sizeof(complex));
00456 VksHart = reinterpret_cast<double*> (VksHartC);
00457 if(rhoRsubplanes>1){
00458 rhoIRXCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
00459 rhoIRXint = reinterpret_cast<double*> (rhoIRXCint);
00460 rhoIRYCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
00461 rhoIRYint = reinterpret_cast<double*> (rhoIRYCint);
00462 rhoIRZCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
00463 rhoIRZint = reinterpret_cast<double*> (rhoIRZCint);
00464 VksHartCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
00465 VksHartint = reinterpret_cast<double*> (VksHartCint);
00466 }
00467 }
00468 PUParray(p,Vks, size);
00469 PUParray(p,density,size);
00470 PUParray(p,rhoIRX, size);
00471 PUParray(p,rhoIRY, size);
00472 PUParray(p,rhoIRZ, size);
00473 PUParray(p,VksHart,size);
00474
00475 if(rhoRsubplanes==1){
00476 PUParray(p,rhoIRXint, rsizeInt);
00477 PUParray(p,rhoIRYint, rsizeInt);
00478 PUParray(p,rhoIRZint, rsizeInt);
00479 PUParray(p,VksHartint,rsizeInt);
00480 }
00481
00482 }
00483
00484
00485
00486
00487
00488
00489 void RhoRealSlab::uPackScale(double *uPackData, double *PackData,double scale){
00490
00491 for(int i=0;i<size;i++){uPackData[i] = PackData[i]*scale;}
00492
00493 }
00494
00495
00496
00497
00498
00499
00500 void RhoRealSlab::scale(double *uPackData,double scale){
00501
00502 for(int i=0;i<size;i++){uPackData[i] *= scale;}
00503
00504 }
00505
00506
00507
00508
00509
00510 void RhoRealSlab::uPackScaleShrink(double *uPackData,double *PackData,double scale){
00511
00512 for(int i=0,i2=0;i<sizeY;i++,i2+=2){
00513 for(int j=i*sizeX;j<(i+1)*sizeX;j++){
00514 uPackData[j] =PackData[(j+i2)]*scale;
00515 }
00516 }
00517
00518 }
00519
00520
00521
00522
00523
00524
00525 void RhoRealSlab::uPackShrink(double *uPackData,double *PackData){
00526
00527 for(int i=0,i2=0;i<sizeY;i++,i2+=2){
00528 for(int j=i*sizeX;j<(i+1)*sizeX;j++){
00529 uPackData[j] =PackData[(j+i2)];
00530 }
00531 }
00532
00533 }
00534
00535
00536
00537
00538
00539 void RhoRealSlab::uPackScaleGrow(double *uPackData,double *PackData,double scale){
00540 for(int i=0,i2=0;i<sizeY;i++,i2+=2){
00541 for(int j=i*sizeX;j<(i+1)*sizeX;j++){
00542 uPackData[(j+i2)] =PackData[j]*scale;
00543 }
00544 int k = (i+1)*sizeX+i2;
00545 uPackData[k] =0.0;
00546 uPackData[(k+1)]=0.0;
00547 }
00548
00549 }
00550
00551
