1 /* -*- c++ -*- */ 2 #ifndef GRID_H 3 #define GRID_H 4 5 #include "Array.h" 6 #include "Vector3D.h" 7 #include "mathutilities.h" 8 #include "FFTComplex.h" 9 #include "Report.h" 10 #include "ScalarStructure.h" 11 12 namespace ProtoMol { 13 //_________________________________________________________________ Grid 14 /** 15 * A simple Grid class using T as interpolation scheme, assuming periodic 16 * boundary conditions. 17 */ 18 19 20 template<class TInterpolation> 21 class Grid { 22 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 23 // Typedef & const 24 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 25 private: 26 struct Int3D {int x; int y; int z;}; 27 /// 3d interpolation 28 struct Interpolation3D { 29 TInterpolation x; 30 TInterpolation y; 31 TInterpolation z; Interpolation3DInterpolation3D32 Interpolation3D(){}; Interpolation3DInterpolation3D33 Interpolation3D(unsigned int order):x(order),y(order),z(order){} 34 }; 35 36 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 37 // Constructors, destructors, assignment 38 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 39 public: 40 Grid(); 41 ~Grid(); 42 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 43 // New methods of class Grid 44 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 45 public: 46 47 48 void anterpolateCharge(Real q, const Vector3D& coord, unsigned int index); fftBack()49 void fftBack(){myFFT.backward();} 50 Real scalarSum(ScalarStructure* energies); 51 Real scalarSum(ScalarStructure* energies, unsigned int block, unsigned int nBlocks); fftForward()52 void fftForward(){myFFT.forward();} 53 void interpolateForce(Real q, unsigned int index, Vector3D& force); 54 55 void initialize(Real width, Real length, Real height, Real alpha, 56 unsigned int nx, unsigned int ny, unsigned int nz, 57 unsigned int interOrder, 58 unsigned int atomCount); 59 60 void clear(); getQ(Real * & begin,Real * & end)61 void getQ(Real*& begin, Real*& end) {begin=&(myQ.begin()->re);end=&(myQ.end()->re);} 62 void print(); 63 private: 64 void dftmod(unsigned int order, unsigned int n, Real* interpolation, Real* interpolationMod); 65 66 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 67 // My data members 68 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 69 private: 70 Array<zomplex,3> myQ; 71 Array<zomplex,3> myQTmp; 72 unsigned int myNX; 73 unsigned int myNY; 74 unsigned int myNZ; 75 int myNXOffset; 76 int myNYOffset; 77 int myNZOffset; 78 Real myWidth; 79 Real myLength; 80 Real myHeight; 81 Real myWidthr; 82 Real myLengthr; 83 Real myHeightr; 84 Real myV; 85 Real myHX; 86 Real myHY; 87 Real myHZ; 88 Real myHXr; 89 Real myHYr; 90 Real myHZr; 91 Real myAlpha; 92 std::vector<Int3D> myScaledParticleIntPositions; 93 std::vector<Interpolation3D> myInterpolations; 94 Real* myInerpolationModX; 95 Real* myInerpolationModY; 96 Real* myInerpolationModZ; 97 Real* myExpX; 98 Real* myExpY; 99 Real* myExpZ; 100 unsigned int myInterOrder; 101 FFTComplex myFFT; 102 unsigned int myAtomCount; 103 Real myFac; 104 }; 105 //______________________________________________________________________ INLINES 106 template<class TInterpolation> anterpolateCharge(Real q,const Vector3D & coord,unsigned int index)107 inline void Grid<TInterpolation>::anterpolateCharge(Real q, const Vector3D& coord, unsigned int index) { 108 Real x = coord.x*myHXr; 109 Real y = coord.y*myHYr; 110 Real z = coord.z*myHZr; 111 while(x < 0.0) x += myNX; 112 while(x >= myNX) x -= myNX; 113 while(y < 0.0) y += myNY; 114 while(y >= myNY) y -= myNY; 115 while(z < 0.0) z += myNZ; 116 while(z >= myNZ) z -= myNZ; 117 int intX = (int)x; 118 int intY = (int)y; 119 int intZ = (int)z; 120 int i0 = intX+myNXOffset; 121 int j0 = intY+myNYOffset; 122 int k0 = intZ+myNZOffset; 123 myScaledParticleIntPositions[index].x = i0; 124 myScaledParticleIntPositions[index].y = j0; 125 myScaledParticleIntPositions[index].z = k0; 126 myInterpolations[index].x.set(x-intX); 127 myInterpolations[index].y.set(y-intY); 128 myInterpolations[index].z.set(z-intZ); 129 Real* thetaX = myInterpolations[index].x.theta; 130 Real* thetaY = myInterpolations[index].y.theta; 131 Real* thetaZ = myInterpolations[index].z.theta; 132 for(unsigned int i=0;i<myInterOrder;i++){ 133 Real a = q*thetaX[i]; 134 int i1 = (i+i0) % myNX; 135 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 136 Array<zomplex,3>::RefArray<2> rQX = myQ[i1]; 137 #else 138 RefArray<zomplex,2> rQX = myQ[i1]; 139 #endif 140 for(unsigned int j=0;j<myInterOrder;j++){ 141 Real ab = a*thetaY[j]; 142 int j1 = (j+j0) % myNY; 143 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 144 Array<zomplex,3>::RefArray<1> rQXY = rQX[j1]; 145 #else 146 RefArray<zomplex,1> rQXY = rQX[j1]; 147 #endif 148 for(unsigned int k=0;k<myInterOrder;k++){ 149 int k1 = (k+k0) % myNZ; 150 rQXY[k1].re += ab*thetaZ[k]; 151 } 152 } 153 } 154 } 155 156 157 158 template<class TInterpolation> interpolateForce(Real q,unsigned int index,Vector3D & force)159 inline void Grid<TInterpolation>::interpolateForce(Real q, unsigned int index, Vector3D& force){ 160 int i0 = myScaledParticleIntPositions[index].x; 161 int j0 = myScaledParticleIntPositions[index].y; 162 int k0 = myScaledParticleIntPositions[index].z; 163 Real fx = 0.0; 164 Real fy = 0.0; 165 Real fz = 0.0; 166 Real* thetaX = myInterpolations[index].x.theta; 167 Real* thetaY = myInterpolations[index].y.theta; 168 Real* thetaZ = myInterpolations[index].z.theta; 169 Real* dThetaX = myInterpolations[index].x.dTheta; 170 Real* dThetaY = myInterpolations[index].y.dTheta; 171 Real* dThetaZ = myInterpolations[index].z.dTheta; 172 for(unsigned int i=0;i<myInterOrder;i++){ 173 int i1 = (i+i0)%myNX; 174 for(unsigned int j=0;j<myInterOrder;j++){ 175 int j1 = (j+j0)%myNY; 176 Real xij = dThetaX[i]*thetaY[j]; 177 Real yij = thetaX[i]*dThetaY[j]; 178 Real zij = thetaX[i]*thetaY[j]; 179 for(unsigned int k=0;k<myInterOrder;k++){ 180 int k1 = (k+k0)%myNZ; 181 Real term = myQ[i1][j1][k1].re; 182 fx -= term*xij*thetaZ[k]; 183 fy -= term*yij*thetaZ[k]; 184 fz -= term*zij*dThetaZ[k]; 185 } 186 } 187 } 188 force += Vector3D(fx*myHXr*q,fy*myHYr*q,fz*myHZr*q); 189 } 190 191 template<class TInterpolation> Grid()192 Grid<TInterpolation>::Grid(): 193 myNX(0), 194 myNY(0), 195 myNZ(0), 196 myNXOffset(0), 197 myNYOffset(0), 198 myNZOffset(0), 199 myWidth(0.0), 200 myLength(0.0), 201 myHeight(0.0), 202 myWidthr(0.0), 203 myLengthr(0.0), 204 myHeightr(0.0), 205 myV(0.0), 206 myHX(0.0), 207 myHY(0.0), 208 myHZ(0.0), 209 myHXr(0.0), 210 myHYr(0.0), 211 myHZr(0.0), 212 myAlpha(0.0), 213 myInerpolationModX(NULL), 214 myInerpolationModY(NULL), 215 myInerpolationModZ(NULL), 216 myExpX(NULL), 217 myExpY(NULL), 218 myExpZ(NULL), 219 myInterOrder(0), 220 myAtomCount(0), 221 myFac(0.0){ 222 } 223 224 template<class TInterpolation> ~Grid()225 Grid<TInterpolation>::~Grid(){ 226 if(myInerpolationModX != NULL) delete [] myInerpolationModX; 227 if(myInerpolationModY != NULL) delete [] myInerpolationModY; 228 if(myInerpolationModZ != NULL) delete [] myInerpolationModZ; 229 if(myExpX != NULL) delete [] myExpX; 230 if(myExpY != NULL) delete [] myExpY; 231 if(myExpZ != NULL) delete [] myExpZ; 232 } 233 234 235 236 template<class TInterpolation> scalarSum(ScalarStructure * energies)237 Real Grid<TInterpolation>::scalarSum(ScalarStructure* energies){ 238 Real energy = 0.0; 239 Real virialxx = 0.0; 240 Real virialxy = 0.0; 241 Real virialxz = 0.0; 242 Real virialyy = 0.0; 243 Real virialyz = 0.0; 244 Real virialzz = 0.0; 245 bool doVirial = energies->virial(); 246 bool doMolVirial = energies->molecularVirial(); 247 Real piVr = 1.0/(M_PI*myV); 248 int count = 0; 249 for (unsigned int i = 0; i < myNX; i++){ 250 int i0 = i <= myNX/2 ? i : i-myNX; 251 Real mi = i0*myWidthr; 252 Real ex = myExpX[i]*piVr; 253 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 254 Array<zomplex,3>::RefArray<2> rQX = myQ[i]; 255 #else 256 RefArray<zomplex,2> rQX = myQ[i]; 257 #endif 258 for (unsigned int j = 0 ; j < myNY; j++){ 259 int j0 = j <= myNY/2 ? j : j-myNY; 260 Real interpolationModXY = myInerpolationModX[i]*myInerpolationModY[j]; 261 Real mj = j0*myLengthr; 262 Real mij = mi*mi + mj*mj; 263 Real exy = ex*myExpY[j]; 264 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 265 Array<zomplex,3>::RefArray<1> rQXY = rQX[j]; 266 #else 267 RefArray<zomplex,1> rQXY = rQX[j]; 268 #endif 269 for (unsigned int k = (i!=0 || j!=0 ? 0:1); k <= myNZ/2; k++){ 270 count++; 271 int k0 = k <= myNZ/2 ? k : k-myNZ; 272 //Vector3D mHat(i0*myWidthr,j0*myLengthr,k0*myHeightr); 273 //Real mHatSquared = mHat.normSquared(); 274 //Real theta = myInerpolationModX[i]*myInerpolationModY[j]*myInerpolationModZ[k]*exp(-fac*mHatSquared)/(mHatSquared*M_PI*myV); 275 Real mk = k0*myHeightr; 276 Real mHatSquared = mij+mk*mk; 277 Real theta = interpolationModXY*myInerpolationModZ[k]*exy*myExpZ[k]/mHatSquared; 278 279 // Energy 280 Real q = power<2>(rQXY[k].re)+power<2>(rQXY[k].im); 281 //Real q = power<2>(myQ[i][j][k].re)+power<2>(myQ[i][j][k].im); 282 Real e = q*theta; 283 Real v = 2.0*(1.0/mHatSquared + myFac); 284 285 // Symmetric 286 if(k > 0 && ((k != myNZ/2) || (myNZ & 1))){ 287 e *= 2.0; 288 zomplex& w = myQ[(myNX-i)%myNX][(myNY-j)%myNY][(myNZ-k)%myNZ]; 289 w.re *= theta; 290 w.im *= theta; 291 } 292 293 energy += e; 294 295 // Virial 296 if(doVirial){ 297 virialxx += e *(1.0 - v * mi*mi); 298 virialxy -= e *(v * mi*mj); 299 virialxz -= e *(v * mi*mk); 300 virialyy += e *(1.0 - v * mj*mj); 301 virialyz -= e *(v * mj*mk); 302 virialzz += e *(1.0 - v * mk*mk); 303 } 304 305 rQXY[k].re *= theta; 306 rQXY[k].im *= theta; 307 //myQ[i][j][k].re *= theta; 308 //myQ[i][j][k].im *= theta; 309 } 310 } 311 } 312 // Just clear (0,0,0) since we did this not the nested loop. 313 myQ[0][0][0].re = 0.0; 314 myQ[0][0][0].im = 0.0; 315 316 if(doVirial){ 317 (*energies)[ScalarStructure::VIRIALXX] += virialxx*0.5; 318 (*energies)[ScalarStructure::VIRIALXY] += virialxy*0.5; 319 (*energies)[ScalarStructure::VIRIALXZ] += virialxz*0.5; 320 (*energies)[ScalarStructure::VIRIALYX] += virialxy*0.5; 321 (*energies)[ScalarStructure::VIRIALYY] += virialyy*0.5; 322 (*energies)[ScalarStructure::VIRIALYZ] += virialyz*0.5; 323 (*energies)[ScalarStructure::VIRIALZX] += virialxz*0.5; 324 (*energies)[ScalarStructure::VIRIALZY] += virialyz*0.5; 325 (*energies)[ScalarStructure::VIRIALZZ] += virialzz*0.5; 326 } 327 328 // Molecular Virial 329 if(doMolVirial) { 330 (*energies)[ScalarStructure::MOLVIRIALXX] += virialxx*0.5; 331 (*energies)[ScalarStructure::MOLVIRIALXY] += virialxy*0.5; 332 (*energies)[ScalarStructure::MOLVIRIALXZ] += virialxz*0.5; 333 (*energies)[ScalarStructure::MOLVIRIALYX] += virialxy*0.5; 334 (*energies)[ScalarStructure::MOLVIRIALYY] += virialyy*0.5; 335 (*energies)[ScalarStructure::MOLVIRIALYZ] += virialyz*0.5; 336 (*energies)[ScalarStructure::MOLVIRIALZX] += virialxz*0.5; 337 (*energies)[ScalarStructure::MOLVIRIALZY] += virialyz*0.5; 338 (*energies)[ScalarStructure::MOLVIRIALZZ] += virialzz*0.5; 339 } 340 return energy*0.5; 341 } 342 343 template<class TInterpolation> scalarSum(ScalarStructure * energies,unsigned int block,unsigned int nBlocks)344 Real Grid<TInterpolation>::scalarSum(ScalarStructure* energies, unsigned int block, unsigned int nBlocks){ 345 Real energy = 0.0; 346 Real virialxx = 0.0; 347 Real virialxy = 0.0; 348 Real virialxz = 0.0; 349 Real virialyy = 0.0; 350 Real virialyz = 0.0; 351 Real virialzz = 0.0; 352 bool doVirial = energies->virial(); 353 bool doMolVirial = energies->molecularVirial(); 354 Real piVr = 1.0/(M_PI*myV); 355 356 myQTmp.resize(ArraySizes(myNX)(myNY)(myNZ)); 357 int m = myQ.size(); 358 zomplex *q = myQ.begin(); 359 zomplex *t = myQTmp.begin(); 360 for(int i=0;i<m;i++){ 361 t[i].re = q[i].re; 362 t[i].im = q[i].im; 363 q[i].re = 0.0; 364 q[i].im = 0.0; 365 } 366 367 int nx = myNX; 368 int ny = myNY; 369 int nz = myNZ/2+1; 370 371 int nyz = ny*nz; 372 int n = nx*nyz; 373 int sn = (n*block)/nBlocks + (block==0?1:0); // Add 1 to skip i,j,k == 0 374 int en = (n*(block+1))/nBlocks - 1; 375 376 int count = 0; 377 int size = en-sn+1; 378 if(size == 0) 379 return 0.0; 380 381 int k = sn % nz; 382 int j = (sn / nz) % ny; 383 int i = (sn / nyz); 384 int ez = (en % nz)+1; 385 int ey = ((en / nz) % ny)+1; 386 int ex = (en / nyz)+1; 387 if(j < ey-1) 388 ez = nz; 389 if(i < ex-1){ 390 ey = ny; 391 ez = nz; 392 } 393 394 for (; i < ex; i++,j=0){ 395 int i0 = i <= static_cast<int>(myNX/2) ? i : i-myNX; 396 Real mi = i0*myWidthr; 397 Real ex = myExpX[i]*piVr; 398 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 399 Array<zomplex,3>::RefArray<2> rQX = myQ[i]; 400 Array<zomplex,3>::RefArray<2> rTQX = myQTmp[i]; 401 #else 402 RefArray<zomplex,2> rQX = myQ[i]; 403 RefArray<zomplex,2> rTQX = myQTmp[i]; 404 #endif 405 for (; j < ey; j++,k=0){ 406 int j0 = j <= static_cast<int>(myNY/2) ? j : j-myNY; 407 Real interpolationModXY = myInerpolationModX[i]*myInerpolationModY[j]; 408 Real mj = j0*myLengthr; 409 Real mij = mi*mi + mj*mj; 410 Real exy = ex*myExpY[j]; 411 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 412 Array<zomplex,3>::RefArray<1> rQXY = rQX[j]; 413 Array<zomplex,3>::RefArray<1> rTQXY = rTQX[j]; 414 #else 415 RefArray<zomplex,1> rQXY = rQX[j]; 416 RefArray<zomplex,1> rTQXY = rTQX[j]; 417 #endif 418 //for (unsigned int k = (i!=0 || j!=0 ? 0:1); k < myNZ; k++){ 419 for (; k < ez; k++){ 420 int k0 = k <= static_cast<int>(myNZ/2) ? k : k-myNZ; 421 //Vector3D mHat(i0*myWidthr,j0*myLengthr,k0*myHeightr); 422 //Real mHatSquared = mHat.normSquared(); 423 //Real theta = myInerpolationModX[i]*myInerpolationModY[j]*myInerpolationModZ[k]*exp(-fac*mHatSquared)/(mHatSquared*M_PI*myV); 424 Real mk = k0*myHeightr; 425 Real mHatSquared = mij+mk*mk; 426 Real theta = interpolationModXY*myInerpolationModZ[k]*exy*myExpZ[k]/mHatSquared; 427 428 // Energy 429 Real q = power<2>(rTQXY[k].re)+power<2>(rTQXY[k].im); 430 //Real q = power<2>(myQ[i][j][k].re)+power<2>(myQ[i][j][k].im); 431 Real e = q*theta; 432 Real v = 2.0*(1.0/mHatSquared + myFac); 433 434 // Symmetric 435 if(k > 0 && ((k != static_cast<int>(myNZ/2)) || (myNZ & 1))){ 436 e *= 2.0; 437 myQ[(myNX-i)%myNX][(myNY-j)%myNY][(myNZ-k)%myNZ].re = myQTmp[(myNX-i)%myNX][(myNY-j)%myNY][(myNZ-k)%myNZ].re * theta; 438 myQ[(myNX-i)%myNX][(myNY-j)%myNY][(myNZ-k)%myNZ].im = myQTmp[(myNX-i)%myNX][(myNY-j)%myNY][(myNZ-k)%myNZ].im * theta; 439 } 440 441 energy += e; 442 443 // Virial 444 if(doVirial){ 445 virialxx += e *(1.0 - v * mi*mi); 446 virialxy -= e *(v * mi*mj); 447 virialxz -= e *(v * mi*mk); 448 virialyy += e *(1.0 - v * mj*mj); 449 virialyz -= e *(v * mj*mk); 450 virialzz += e *(1.0 - v * mk*mk); 451 } 452 // 453 rQXY[k].re = rTQXY[k].re * theta; 454 rQXY[k].im = rTQXY[k].im * theta; 455 //myQ[i][j][k].re *= theta; 456 //myQ[i][j][k].im *= theta; 457 458 count++; 459 if(count >= size){ 460 i = myNX; 461 j = myNY; 462 k = myNZ; 463 } 464 } 465 } 466 } 467 // Just clear (0,0,0) since we did this not the nested loop. 468 myQ[0][0][0].re = 0.0; 469 myQ[0][0][0].im = 0.0; 470 471 if(doVirial){ 472 (*energies)[ScalarStructure::VIRIALXX] += virialxx*0.5; 473 (*energies)[ScalarStructure::VIRIALXY] += virialxy*0.5; 474 (*energies)[ScalarStructure::VIRIALXZ] += virialxz*0.5; 475 (*energies)[ScalarStructure::VIRIALYX] += virialxy*0.5; 476 (*energies)[ScalarStructure::VIRIALYY] += virialyy*0.5; 477 (*energies)[ScalarStructure::VIRIALYZ] += virialyz*0.5; 478 (*energies)[ScalarStructure::VIRIALZX] += virialxz*0.5; 479 (*energies)[ScalarStructure::VIRIALZY] += virialyz*0.5; 480 (*energies)[ScalarStructure::VIRIALZZ] += virialzz*0.5; 481 } 482 483 // Molecular Virial 484 if(doMolVirial) { 485 (*energies)[ScalarStructure::MOLVIRIALXX] += virialxx*0.5; 486 (*energies)[ScalarStructure::MOLVIRIALXY] += virialxy*0.5; 487 (*energies)[ScalarStructure::MOLVIRIALXZ] += virialxz*0.5; 488 (*energies)[ScalarStructure::MOLVIRIALYX] += virialxy*0.5; 489 (*energies)[ScalarStructure::MOLVIRIALYY] += virialyy*0.5; 490 (*energies)[ScalarStructure::MOLVIRIALYZ] += virialyz*0.5; 491 (*energies)[ScalarStructure::MOLVIRIALZX] += virialxz*0.5; 492 (*energies)[ScalarStructure::MOLVIRIALZY] += virialyz*0.5; 493 (*energies)[ScalarStructure::MOLVIRIALZZ] += virialzz*0.5; 494 } 495 return energy*0.5; 496 } 497 498 template<class TInterpolation> initialize(Real width,Real length,Real height,Real alpha,unsigned int nx,unsigned int ny,unsigned int nz,unsigned int interOrder,unsigned int atomCount)499 void Grid<TInterpolation>::initialize(Real width, Real length, Real height, Real alpha, 500 unsigned int nx, unsigned int ny, unsigned int nz, 501 unsigned int interOrder, 502 unsigned int atomCount){ 503 if(!myQ.resize(ArraySizes(nx)(ny)(nz))) 504 report << error <<"[Grid<>::initialize] Could not allocate memory for Q[" 505 <<nx<<"]["<<ny<<"]["<<nz<<"]."<<endr; 506 myNX = nx; 507 myNY = ny; 508 myNZ = nz; 509 myNXOffset = -((int)interOrder-1)/2 + (int)nx; 510 myNYOffset = -((int)interOrder-1)/2 + (int)ny; 511 myNZOffset = -((int)interOrder-1)/2 + (int)nz; 512 myWidth = width; 513 myLength = length; 514 myHeight = height; 515 myWidthr = 1.0/width; 516 myLengthr = 1.0/length; 517 myHeightr = 1.0/height; 518 myV = width*length*height; 519 myHX = width/nx; 520 myHY = length/ny; 521 myHZ = height/nz; 522 myHXr = nx/width; 523 myHYr = ny/length; 524 myHZr = nz/height; 525 myAlpha = alpha; 526 myInterOrder = interOrder; 527 myAtomCount = atomCount; 528 myFac = M_PI*M_PI/(myAlpha*myAlpha); 529 530 // Precompute exp(-pi^2m^2/alpha^2) 531 if(myExpX != NULL) 532 delete [] myExpX; 533 myExpX = new Real[nx]; 534 if(myExpY != NULL) 535 delete [] myExpY; 536 myExpY = new Real[ny]; 537 if(myExpZ != NULL) 538 delete [] myExpZ; 539 myExpZ = new Real[nz]; 540 541 for (unsigned int i = 0; i < myNX; i++){ 542 int i0 = i <= myNX/2 ? i : i-myNX; 543 myExpX[i] = exp(-myFac*power<2>(i0*myWidthr)); 544 } 545 for (unsigned int j = 0 ; j < myNY; j++){ 546 int j0 = j <= myNY/2 ? j : j-myNY; 547 myExpY[j] = exp(-myFac*power<2>(j0*myLengthr)); 548 } 549 550 //for (unsigned int k = 0; k < myNZ; k++){ 551 for (unsigned int k = 0; k<= myNZ/2 ; k++){ 552 int k0 = k <= myNZ/2 ? k : k-myNZ; 553 myExpZ[k] = exp(-myFac*power<2>(k0*myHeightr)); 554 } 555 556 557 // Precompute the mod TInterpolation, B(m1,m2,m3) 558 if(myInerpolationModX != NULL) 559 delete [] myInerpolationModX; 560 myInerpolationModX = new Real[nx]; 561 if(myInerpolationModY != NULL) 562 delete [] myInerpolationModY; 563 myInerpolationModY = new Real[ny]; 564 if(myInerpolationModZ != NULL) 565 delete [] myInerpolationModZ; 566 myInerpolationModZ = new Real[nz]; 567 568 TInterpolation interpolation = TInterpolation(myInterOrder,0.0); 569 dftmod(myInterOrder,nx,interpolation.theta,myInerpolationModX); 570 dftmod(myInterOrder,ny,interpolation.theta,myInerpolationModY); 571 dftmod(myInterOrder,nz,interpolation.theta,myInerpolationModZ); 572 //for(unsigned int i=0;i<nx;i++) 573 // report << plain <<"myInerpolationModX["<<i<<"]:"<<1.0/myInerpolationModX[i]<<endr; 574 //for(unsigned int i=0;i<ny;i++) 575 // report << plain <<"myInerpolationModY["<<i<<"]:"<<1.0/myInerpolationModY[i]<<endr; 576 //for(unsigned int i=0;i<nz;i++) 577 // report << plain <<"myInerpolationModZ["<<i<<"]:"<<1.0/myInerpolationModZ[i]<<endr; 578 579 580 581 // Resize the vector data members 582 myScaledParticleIntPositions.resize(myAtomCount); 583 myInterpolations.resize(myAtomCount,Interpolation3D(myInterOrder)); 584 585 myFFT.initialize(myNX,myNY,myNZ,&myQ[0][0][0]); 586 } 587 588 template<class TInterpolation> dftmod(unsigned int order,unsigned int n,Real * interpolation,Real * interpolationMod)589 void Grid<TInterpolation>::dftmod(unsigned int order, unsigned int n, Real* interpolation, Real* interpolationMod){ 590 for(unsigned int i=0;i<n;i++){ 591 Real sumCos = 0.0; 592 Real sumSin = 0.0; 593 for(unsigned int j=0;j<order;j++){ 594 Real x = M_PI*2.0*i*j/(Real)n; 595 sumCos += interpolation[j]*cos(x); 596 sumSin += interpolation[j]*sin(x); 597 } 598 interpolationMod[i] = 1.0/(sumCos*sumCos + sumSin*sumSin); 599 } 600 } 601 602 template<class TInterpolation> clear()603 void Grid<TInterpolation>::clear(){ 604 int n = myQ.size(); 605 zomplex *q = myQ.begin(); 606 for(int i=0;i<n;i++){ 607 q[i].re = 0.0; 608 q[i].im = 0.0; 609 } 610 } 611 612 template<class TInterpolation> print()613 void Grid<TInterpolation>::print(){ 614 Real q = 0.0; 615 report << plain; 616 for (unsigned int i = 0; i < myNX; i++){ 617 for (unsigned int j = 0 ; j < myNY; j++){ 618 report << "Q["<<i<<"]["<<j<<"][0-"<<myNZ-1<<"] : "; 619 for (unsigned int k = 0; k < myNZ; k++){ 620 report << "("<<myQ[i][j][k].re <<","<<myQ[i][j][k].im<<")"; 621 q += myQ[i][j][k].re; 622 } 623 report <<std::endl; 624 } 625 } 626 report <<"Sum Q :"<<q<<endr; 627 } 628 } 629 #endif 630