MD++: My Documents/Codes/MD++.svn/trunk/src/ewald.cpp Source File

00001 /*
00002   ewald.cpp
00003   by Hark Lee and Wei Cai, hblee@stanford.edu, caiwei@stanford.edu
00004   Last Modified : Thu Oct 18 18:20:05 2007
00005 
00006   FUNCTION  :  Ewald Summation (classical Ewald and Particle Mesh Ewald)
00007 
00008   To Do:
00009       1.  CE Virial stress still incorrect (perhaps from Self energy contribution)
00010       2.  CG relax dies early when both Y and Oxygen vacancy exist
00011       3.  PME still needs to be included
00012 
00013   Note:
00014       This file is included by md.cpp.  This file cannot be compiled by itself.
00015 */
00016 
00017 /********************************************************************/
00018 /* Ewald Summation for Coulomb Interaction */
00019 /********************************************************************/
00020 
00021 //###################################################################
00022 //  CE Functions
00023 //###################################################################
00024 
00025 //-------------------------------------------------------------------
00026 // Warning
00027 // Note that if you want to switch conj_fixbox (from 1 to 0) or (from 
00028 // 0 to 1) in the middle of simulation, CE should be reinitialized.
00029 //-------------------------------------------------------------------
00030 
00031 //-------------------------------------------------------------------
00032 // Classical Ewald
00033 //-------------------------------------------------------------------
00034 void MDFrame::CE()
00035 {
00036     int i;
00037     CE_Real();
00038     CE_Rec();
00039 
00040     _EPOT_Ewald = _EPOT_Ewald_Real + _EPOT_Ewald_Rec + _EPOT_Ewald_Self;
00041     _VIRIAL_Ewald = _VIRIAL_Ewald_Real + _VIRIAL_Ewald_Rec;
00042 
00043     for(i=0;i<_NP;i++)
00044       {
00045         _F_Ewald[i] = _F_Ewald_Real[i] + _F_Ewald_Rec[i];
00046         _EPOT_IND_Ewald[i] = _EPOT_IND_Ewald_Real[i] + _EPOT_IND_Ewald_Rec[i];
00047       }
00048 }
00049 
00050 //-------------------------------------------------------------------
00051 // Compute Real Part of Ewald Summation
00052 //-------------------------------------------------------------------
00053 void MDFrame::CE_Real()
00054 {
00055     int i,j,ipt,jpt,isp,jsp; // integers for atom counting, atom indices, atom species
00056     double UC,r,r2,QQ,EXP;   // U = potential | r = radius | r2 = r^2 | ri6 = r^-6
00057     Vector3 sij, rij, fijC; // sij : scaled coordinates | rij : coordinates | fij : forces
00058 
00059     DUMP(HIG"Coulomb Potential (real)"NOR);
00060 
00061     // start time of real part
00062     CE_sREAL=clock();
00063     
00064     // refresh neighbor list
00065     refreshneighborlist();
00066     
00067     // initialize force, potential and virials
00068     _EPOT_Ewald_Real = 0.0;
00069     _VIRIAL_Ewald_Real.clear();
00070 
00071     for(i=0;i<_NP;i++)
00072     {
00073         _F_Ewald_Real[i].clear();
00074         _EPOT_IND_Ewald_Real[i]=0.0;
00075     }
00076     
00077     //-------------------------------------------------------------------
00078     //       Calculating Short Range Potential and Ewald Real Part
00079     //-------------------------------------------------------------------
00080     for(ipt=0;ipt<_NP;ipt++) // for all particles
00081     {
00082         /* if fixed == -1, simply ignore this atom */
00083         if(fixed[ipt]==-1) continue;
00084         for(j=0;j<nn[ipt];j++) // for all neighbors
00085         {
00086             jpt=nindex[ipt][j]; // get the atom index of j-th neighbor
00087             /* if fixed == -1, simply ignore this atom */
00088             if(fixed[jpt]==-1) continue;
00089             
00090             if(ipt>=jpt) continue; // only if i<j -> avoid redundancy
00091             sij=_SR[jpt]-_SR[ipt]; // scaled interatomics distance
00092             sij.subint(); // apply periodic boundary condition (PBC)
00093             rij=_H*sij; // get the interatomic distance vector (with PBC)
00094             r2=rij.norm2(); // get the square of interatomic distance
00095             r=sqrt(r2); // interatomic distance
00096             
00097             if(r<=Ewald_Rcut) // if the interatomic distance is shorter than the cut-off radius
00098             {
00099                 // get the species indices
00100                 isp=species[ipt]; // species of ipt-atom
00101                 jsp=species[jpt]; // species of jpt-atom
00102                 QQ=P_CLMB * _ATOMCHARGE[isp]*_ATOMCHARGE[jsp]; // get the charge term
00103                 
00104                 // calculate potential
00105                 UC=(erfc(Ewald_Alpha*r))/r*QQ;
00106                 
00107                 // calculate force
00108                 EXP=(2.0/M_SQRT_PI)*Ewald_Alpha*exp(-SQR(Ewald_Alpha*r))*QQ;
00109                 fijC=rij*((UC+EXP)/r2);
00110                 
00111                 // add pair terms to the accumulative sum
00112                 _F_Ewald_Real[ipt]-=fijC; _F_Ewald_Real[jpt]+=fijC;
00113 
00114                 _EPOT_Ewald_Real+=UC;
00115 
00116                 _EPOT_IND_Ewald_Real[ipt]+=UC*0.5;
00117                 _EPOT_IND_Ewald_Real[jpt]+=UC*0.5;
00118                 _VIRIAL_Ewald_Real.addnvv(1.,fijC,rij);
00119             }
00120         }
00121     }
00122     
00123     // end time of real part
00124     CE_eREAL=clock();
00125 }
00126 
00127 //-------------------------------------------------------------------
00128 // Compute Reciprocal Part of Ewald Summation
00129 //-------------------------------------------------------------------
00130 void MDFrame::CE_Rec()
00131 {    
00132     int i, ik, kal, kbl;
00133     int ka, kb, kc, np_real;
00134     double k2, kx, ky, kz;
00135     double pfac, tmp_H_L1, tmp_H_L2, tmp_H_L3;
00136     double UE, EXP;
00137     EComplex sincosk;
00138 
00139      // start time of reciprocal part
00140     CE_sREC=clock();
00141     
00142     // initialize force, potential and virials
00143    _EPOT_Ewald_Rec = 0.0;
00144     _VIRIAL_Ewald_Rec.clear();
00145 
00146     for(i=0;i<_NP;i++)
00147     {
00148         _F_Ewald_Rec[i].clear();
00149         _EPOT_IND_Ewald_Rec[i]=0.0;
00150     }
00151     
00152     // recalculate the system volume (in case the box is changing)
00153     if(conj_fixbox==0)
00154       {
00155         tmp_H_L1=sqrt(SQR(_H[0][0])+SQR(_H[1][0])+SQR(_H[2][0]));
00156         tmp_H_L2=sqrt(SQR(_H[0][1])+SQR(_H[1][1])+SQR(_H[2][1]));
00157         tmp_H_L3=sqrt(SQR(_H[0][2])+SQR(_H[1][2])+SQR(_H[2][2]));
00158         Ewald_cell_V=_H.det();        
00159       }
00160     else
00161       {
00162         // just to avoid warning
00163         tmp_H_L1=Ewald_H_L1;
00164         tmp_H_L2=Ewald_H_L2;
00165         tmp_H_L3=Ewald_H_L3;
00166       }
00167 
00168     // precalculate exp(iKx*rx), exp(iKy*ry), exp(iKz*rz)
00169     CE_fillSinCosTables();
00170     
00171     kal=-1; // an impossible value
00172     kbl=0;  // Just to avoid a warning
00173     for(ik=0;ik<CE_nKV;ik++) // for all K-vectors
00174       {
00175         // kx, ky, kz calculated in enumKV()
00176         ka=CE_KV[ik].a;
00177         kb=CE_KV[ik].b;
00178         kc=CE_KV[ik].c;
00179         
00180         if(conj_fixbox==0)
00181           {
00182             kx=CE_KV[ik].kx*Ewald_H_L1/tmp_H_L1;
00183             ky=CE_KV[ik].ky*Ewald_H_L2/tmp_H_L2;
00184             kz=CE_KV[ik].kz*Ewald_H_L3/tmp_H_L3;
00185             k2=kx*kx+ky*ky+kz*kz;
00186           }
00187         else
00188           {
00189             kx=CE_KV[ik].kx;
00190             ky=CE_KV[ik].ky;
00191             kz=CE_KV[ik].kz;
00192             k2=CE_KV[ik].k2;
00193           }
00194         
00195         //quick sckk [ Qj*exp(iKr) ] calculation
00196         if(kal!=ka || kbl!=kb)  // Re-calculate sc_axby
00197           {
00198             kal=ka;
00199             kbl=kb;
00200             if(kb>=0) // if ky>=0
00201               for(i=0;i<_NP;i++) // for all atoms
00202                 EComplex::Mul(CE_scx[ka][i], CE_scy[kb][i], CE_sc_axby[i]);
00203             else      // if ky<0
00204               for(i=0;i<_NP;i++) // for all atoms
00205                 EComplex::MulConjugate(CE_scx[ka][i], CE_scy[-kb][i], CE_sc_axby[i]);
00206           }
00207         // initialize structure factor 
00208         sincosk.Re=0.0;
00209         sincosk.Im=0.0;
00210         
00211         if(kc>=0) // if kz>=0
00212           for(i=0;i<_NP;i++) 
00213             {
00214               /* if fixed == -1, simply ignore this atom */
00215               if(fixed[i]==-1) continue;
00216               
00217               EComplex::Mul(CE_sc_axby[i], CE_scz[kc][i], CE_sckk[i]);
00218               CE_sckk[i]*=P_SQRT_CLMB*_ATOMCHARGE[species[i]];
00219               sincosk+=CE_sckk[i];
00220             }
00221         else      // if kz<0
00222           for(i=0;i<_NP;i++)
00223             {
00224               /* if fixed == -1, simply ignore this atom */
00225               if(fixed[i]==-1) continue;
00226               
00227               EComplex::MulConjugate(CE_sc_axby[i], CE_scz[-kc][i], CE_sckk[i]);
00228               CE_sckk[i]*=P_SQRT_CLMB*_ATOMCHARGE[species[i]];
00229               sincosk+=CE_sckk[i];
00230             }
00231         
00232         // calculate the reciprocal energy pre-factor
00233         EXP=((8.0*M_PI)/Ewald_cell_V)*exp(-k2/(4.0*SQR(Ewald_Alpha)))/k2;
00234         
00235         // calculate the reciprocal force
00236         for(i=0;i<_NP;i++)
00237           {
00238             /* if fixed == -1, simply ignore this atom */
00239             if(fixed[i]==-1) continue;
00240             
00241             double ad=EXP*(CE_sckk[i].Im*sincosk.Re-CE_sckk[i].Re*sincosk.Im);
00242             _F_Ewald_Rec[i].x+=ad*((double)kx);
00243             _F_Ewald_Rec[i].y+=ad*((double)ky);
00244             _F_Ewald_Rec[i].z+=ad*((double)kz);
00245           }
00246         pfac=1.0/k2+1/(4.0*SQR(Ewald_Alpha));
00247         UE=EXP/2.0*sincosk.Norm2();
00248         
00249         // Potential Energy Calculation        
00250         _EPOT_Ewald_Rec+=UE;
00251         
00252         /* Temporary: distribute rec energy equally to all atoms */
00253         np_real = 0;
00254         for(i=0;i<_NP;i++)
00255           if(fixed[i]!=-1) np_real++;
00256         
00257         for(i=0;i<_NP;i++)
00258           {        
00259             /* if fixed == -1, simply ignore this atom */
00260             if(fixed[i]==-1) continue;
00261             
00262             _EPOT_IND_Ewald_Rec[i]=UE/np_real;
00263           }
00264         
00265         // Virial Stress Calculation
00266         _VIRIAL_Ewald_Rec[0][0]+=(UE*(1.0-2.0*(kx*kx)*pfac)); //xx
00267         _VIRIAL_Ewald_Rec[1][1]+=(UE*(1.0-2.0*(ky*ky)*pfac)); //yy
00268         _VIRIAL_Ewald_Rec[2][2]+=(UE*(1.0-2.0*(kz*kz)*pfac)); //zz
00269         _VIRIAL_Ewald_Rec[0][1]+=(UE*(-2.0*(kx*ky)*pfac));  //xy
00270         _VIRIAL_Ewald_Rec[1][0]+=(UE*(-2.0*(kx*ky)*pfac));  //yx
00271         _VIRIAL_Ewald_Rec[0][2]+=(UE*(-2.0*(kx*kz)*pfac));  //xz
00272         _VIRIAL_Ewald_Rec[2][0]+=(UE*(-2.0*(kx*kz)*pfac));  //zx
00273         _VIRIAL_Ewald_Rec[1][2]+=(UE*(-2.0*(ky*kz)*pfac));  //yz
00274         _VIRIAL_Ewald_Rec[2][1]+=(UE*(-2.0*(ky*kz)*pfac));  //zy
00275       }
00276     //INFO("[CE] Reciprocal Energy Calculated by CE = "<<_EPOT_Ewald_Rec);
00277     
00278     // start time of reciprocal part
00279     CE_eREC=clock();
00280 }
00281 
00282 //-------------------------------------------------------------------
00283 // Initialize Ewald  (which will call CE_init() or PME_init
00284 //-------------------------------------------------------------------
00285 void MDFrame::Ewald_init()
00286 {
00287     Matrix33 hinv;
00288     double tmp1, tmp2, tmp3, tmp4, tmp_RLIST;
00289     int i;
00290     
00291     // get the system volume
00292     Ewald_cell_V=_H.det();
00293 
00294     // Calculate inverse matrix and length of lattice vectors
00295     hinv = _H.inv();
00296     Ewald_H_L1=sqrt(SQR(_H[0][0])+SQR(_H[1][0])+SQR(_H[2][0]));
00297     Ewald_H_L2=sqrt(SQR(_H[0][1])+SQR(_H[1][1])+SQR(_H[2][1]));
00298     Ewald_H_L3=sqrt(SQR(_H[0][2])+SQR(_H[1][2])+SQR(_H[2][2]));
00299 
00300     // Allocate real space force
00301     Realloc(_F_Ewald,     Vector3,_NP);
00302     Realloc(_F_Ewald_Real,Vector3,_NP);
00303     Realloc(_F_Ewald_Rec ,Vector3,_NP);
00304 
00305     Realloc(_EPOT_IND_Ewald,     double,_NP);
00306     Realloc(_EPOT_IND_Ewald_Real,double,_NP);
00307     Realloc(_EPOT_IND_Ewald_Rec ,double,_NP);
00308     
00309     bindvar("F_Ewald",_F_Ewald,DOUBLE);
00310     bindvar("F_Ewald_Real",_F_Ewald_Real,DOUBLE);
00311     bindvar("F_Ewald_Rec", _F_Ewald_Rec, DOUBLE);
00312     
00313     bindvar("EPOT_IND_Ewald",_EPOT_IND_Ewald,DOUBLE);
00314     bindvar("EPOT_IND_Ewald_Real",_EPOT_IND_Ewald_Real,DOUBLE);
00315     bindvar("EPOT_IND_Ewald_Rec", _EPOT_IND_Ewald_Rec, DOUBLE);
00316     
00317     //-------------------------------------------------------------------    
00318     // Detemine the Ewald accuracy
00319     //-------------------------------------------------------------------
00320     // Ewaldprecision=sqrt(-ln(error*Ewaldprecision^2))
00321     //-------------------------------------------------------------------
00322     // sqrt(-ln(10^-4)) = 3.03485425877029
00323     // sqrt(-ln(10^-5)) = 3.39307021220756
00324     // sqrt(-ln(10^-6)) = 3.71692218884984
00325     // sqrt(-ln(10^-7)) = 4.01473481701573
00326     // sqrt(-ln(10^-8)) = 4.29193205257869
00327     // sqrt(-ln(10^-9)) = 4.55228138815544
00328     // sqrt(-ln(10^-10))= 4.79852591218808
00329     //-------------------------------------------------------------------
00330 
00331     if(Ewald_CE_or_PME==0)
00332     { /* CE: classical Ewald */
00333         if( Ewald_option_Alpha==0 )
00334         {
00335             // determine Alpha (for CE only)
00336             if( (Ewald_time_ratio==0) || (Ewald_precision==0))
00337             {
00338                 ERROR("Ewald_time_ratio or Ewald_precision is Not Specified -> Check Your Script File");
00339             }
00340             // get the optimal alpha
00341             Ewald_Alpha = M_SQRT_PI*cbrt(sqrt(Ewald_time_ratio*_NP)/Ewald_cell_V);
00342             Ewald_Rcut = Ewald_precision / Ewald_Alpha;
00343             
00344             if (_SKIN <= 0) _SKIN = 1.0;
00345             tmp_RLIST = Ewald_Rcut +_SKIN;
00346             tmp1=Ewald_H_L1/3.0;
00347             tmp2=Ewald_H_L2/3.0;
00348             tmp3=Ewald_H_L3/3.0;
00349             tmp4=tmp_RLIST;
00350             
00351             if(tmp4>tmp1){tmp4=tmp1;}
00352             if(tmp4>tmp2){tmp4=tmp2;}
00353             if(tmp4>tmp3){tmp4=tmp3;}
00354             
00355             if(tmp_RLIST>tmp4)
00356             {
00357                 WARNING("PP_RC Computed from Optimum Alpha is Too Big");
00358                 WARNING("Automatically Setting _RLIST Equal to 1/3 of Minimum Cell Length");
00359                 Ewald_Rcut = (tmp4*0.95) - _SKIN;
00360             }
00361         }
00362     
00363         if(Ewald_option_Alpha==1)
00364         {
00365             if( (Ewald_Rcut==0) || (Ewald_precision==0))
00366             {
00367                 ERROR("Ewald_Rcut or Ewald_precision is Not Specified -> Check Your Script File");
00368             }
00369             // get the controlled alpha
00370             Ewald_Alpha = Ewald_precision / Ewald_Rcut;
00371         }
00372     }
00373     else
00374     { /* PME: particle mesh Ewald */        
00375         // check if Alpha has non-zero value
00376         if((Ewald_option_Alpha==2)&&(Ewald_Alpha==0))
00377         {
00378             WARNING("Ewald Parameter Alpha is Not Specified -> Check Your Script File");
00379         }
00380     }
00381 
00382     
00383     // Update Verlet list cut-off
00384     _RLIST = Ewald_Rcut + _SKIN;
00385     
00386     // calculate sum of square charge
00387     Ewald_qSQRsum=0.0;
00388     for(i=0;i<_NP;i++)
00389     {
00390         /* if fixed == -1, simply ignore this atom */
00391         if(fixed[i]==-1) continue;
00392         Ewald_qSQRsum+=SQR(P_SQRT_CLMB*_ATOMCHARGE[species[i]]);
00393     }
00394 
00395     // compute the self-interaction energy
00396     _EPOT_Ewald_Self = - Ewald_Alpha / M_SQRT_PI * Ewald_qSQRsum;
00397 
00398     if(Ewald_CE_or_PME == 0)
00399     {
00400         CE_init();
00401     }
00402     else
00403     {
00404         PME_init();
00405     }
00406 }
00407 
00408 //-------------------------------------------------------------------
00409 // Initialize Classical Ewald 
00410 //-------------------------------------------------------------------
00411 void MDFrame::CE_init()
00412 {
00413     // Allocating sc_axby and sckk
00414     Realloc(CE_sckk,EComplex,(_NP*2));
00415     CE_sc_axby=CE_sckk+_NP;
00416 
00417     // determine the cutoff radius
00418     CE_Kc = 2*Ewald_Alpha*Ewald_precision; // reciprocal-space
00419       
00420     // Calculate the K-vectors
00421     CE_enumKV();
00422 
00423     // Print Out Ewald Information
00424     INFO("------------------------------------------------------------------------------");
00425     INFO("------------------------------------------------------------------------------");
00426     INFO("              Classical Ewald Parameters [Long Range Interaction] ");
00427     INFO("------------------------------------------------------------------------------");
00428     INFO("    Ewald Fourier Part Cutoff Radius (Kc)         = "<<CE_Kc);
00429     INFO("    Ewald Gaussian Factor (Alpha)                 = "<<Ewald_Alpha);
00430     INFO("    Total Number of K-Vectors                     = "<<CE_nKV);
00431     INFO("    Max Number of K-Vectors                       = "<<CE_MKV);
00432     INFO("    Max K-Points in A-direction                   = "<<CE_Mka);
00433     INFO("    Max K-Points in B-direction                   = "<<CE_Mkb);
00434     INFO("    Max K-Points in C-direction                   = "<<CE_Mkc);
00435     INFO("------------------------------------------------------------------------------");
00436     INFO("------------------------------------------------------------------------------");
00437 }
00438 
00439 //-------------------------------------------------------------------
00440 //  Clear Classical Ewald
00441 //-------------------------------------------------------------------
00442 void MDFrame::CE_clear()
00443 {
00444     // free variables used in the Classical Ewald summation    
00445     if(CE_scx) { free(CE_scx[0]); free(CE_scx); }
00446     if(CE_sckk){ free(CE_sckk); free(CE_KV); }
00447 }
00448 
00449 //-------------------------------------------------------------------
00450 // Filling Sin & Cos Table
00451 //-------------------------------------------------------------------
00452 // calculate exp(-iKx*rx), exp(-iKy*ry), exp(-iKz*rz) 
00453 // for all atoms and all K-vectors, should be updated every time step
00454 //-------------------------------------------------------------------
00455 void MDFrame::CE_fillSinCosTables()
00456 {
00457     int i,ki;
00458     EComplex *sc1;
00459     for(i=0;i<_NP;i++)
00460     {
00461         // when K-vector is zero vector (kx=0, ky=0, kz=0)
00462         CE_scx[0][i].Re=1.0; CE_scx[0][i].Im=0.0;
00463         CE_scy[0][i].Re=1.0; CE_scy[0][i].Im=0.0;
00464         CE_scz[0][i].Re=1.0; CE_scz[0][i].Im=0.0;
00465     }
00466     if(CE_Mka>0) 
00467     {
00468         // calculate the exp(-ikx*rx) when kx=1 (for all atoms)
00469         sc1=CE_scx[1];
00470         for(i=0;i<_NP;i++)
00471             sc1[i].ExpI(2.0*M_PI*_SR[i].x);
00472         // calculate all of the term for kx>=2 (for all atoms)
00473         for(ki=2;ki<=CE_Mka;ki++) for(i=0;i<_NP;i++)
00474             EComplex::Mul(CE_scx[ki-1][i], sc1[i], CE_scx[ki][i]);
00475     }
00476     if(CE_Mkb>0) 
00477     {
00478         // calculate the exp(-iky*ry) when ky=1 (for all atoms)
00479         sc1=CE_scy[1];
00480         for(i=0;i<_NP;i++)
00481             sc1[i].ExpI(2.0*M_PI*_SR[i].y);
00482         // calculate all of the term for ky>=2 (for all atoms)
00483         for(ki=2;ki<=CE_Mkb;ki++) for(i=0;i<_NP;i++)
00484             EComplex::Mul(CE_scy[ki-1][i], sc1[i], CE_scy[ki][i]);
00485     }
00486     if(CE_Mkc>0) 
00487     {
00488         // calculate the exp(-ikz*rz) when kz=1 (for all atoms)
00489         sc1=CE_scz[1];
00490         for(i=0;i<_NP;i++)
00491             sc1[i].ExpI(2.0*M_PI*_SR[i].z);
00492         // calculate all of the terms when kz=>2 (for all atoms)
00493         for(ki=2;ki<=CE_Mkc;ki++) for(i=0;i<_NP;i++)
00494             EComplex::Mul(CE_scz[ki-1][i], sc1[i], CE_scz[ki][i]);
00495     }
00496 }
00497 
00498 //-------------------------------------------------------------------
00499 // Enumerate K-Vectors
00500 //-------------------------------------------------------------------
00501 void MDFrame::CE_enumKV()
00502 {
00503     //------------------------------------------------------------------- 
00504     //  In this function, following variables are prepared
00505     //  Mka, Mkb, Mkc; MKV, nKV, KV
00506     //-------------------------------------------------------------------
00507     const double INHOM=1.1; //The expected error in k-point estimation
00508     
00509     EComplex *p;
00510     int i;
00511     double kx, ky, kz, k2;
00512     int nka, nkb, nkc, ika, ikb, ikc;
00513         
00514     // Kc = (2*pi/Li)*ni and ni = Kc*Li/(2*pi)
00515     // So we can multiply Kc and 1/Gi to get the integer ni
00516     nka=(int)(CE_Kc*Ewald_H_L1/(2.0*M_PI));
00517     nkb=(int)(CE_Kc*Ewald_H_L2/(2.0*M_PI));
00518     nkc=(int)(CE_Kc*Ewald_H_L3/(2.0*M_PI));
00519     
00520     // Get the max-number of K-vectors (half sphere)
00521     ika=(int)(INHOM/3.0*CE_Kc*CE_Kc*CE_Kc*Ewald_cell_V/(2.0*M_PI*2.0*M_PI));
00522 
00523     // Allocate K-vectors
00524     if(CE_MKV < ika) {CE_MKV=ika; Realloc(CE_KV,K_rec,CE_MKV);}
00525     
00526     // find K-vectors
00527     CE_nKV=0;
00528     CE_Mka=CE_Mkb=CE_Mkc=0;
00529     for(ika=0;ika<nka+1;ika++)
00530         for(ikb= (ika==0? 0:-nkb);ikb<=nkb;ikb++)
00531             for(ikc=((ika==0 && ikb==0)?1:-nkc);ikc<=nkc;ikc++)
00532             {
00533                 //Notice that this (kx,ky,kz) is row vector, mul in left (??)
00534                 kx=(2.0*M_PI)*(ika/Ewald_H_L1);
00535                 ky=(2.0*M_PI)*(ikb/Ewald_H_L2);
00536                 kz=(2.0*M_PI)*(ikc/Ewald_H_L3);
00537                 k2=kx*kx+ky*ky+kz*kz;
00538                 if(k2 > CE_Kc*CE_Kc) continue;
00539                 CE_KV[CE_nKV].kx=kx;
00540                 CE_KV[CE_nKV].ky=ky;
00541                 CE_KV[CE_nKV].kz=kz;
00542                 CE_KV[CE_nKV].k2=k2;
00543                 CE_KV[CE_nKV].a=ika;
00544                 CE_KV[CE_nKV].b=ikb;
00545                 CE_KV[CE_nKV].c=ikc;
00546                 CE_nKV++;
00547                 /*
00548                 if(CE_nKV==10000)
00549                   {
00550                     INFO("_H_L1 = "<<Ewald_H_L1);
00551                     INFO("_H_L2 = "<<Ewald_H_L2);
00552                     INFO("_H_L3 = "<<Ewald_H_L3);
00553                     INFO("kx = "<<kx);
00554                     INFO("ky = "<<ky);
00555                     INFO("kz = "<<kz);
00556                     INFO("k2 = "<<k2);
00557                   }
00558                 */
00559                 
00560                 if(CE_nKV >= CE_MKV) INFO("Error : enumKV() -> Out of KV space");
00561                 if(CE_Mka < ika) CE_Mka=ika;
00562                 if(CE_Mkb < abs(ikb)) CE_Mkb=abs(ikb); // Mkb can be negative
00563                 if(CE_Mkc < abs(ikc)) CE_Mkc=abs(ikc); // Mkc can be negative
00564             }
00565     if(CE_scx==NULL) p=0;
00566     else p=CE_scx[0]; // assign p to scx[0]
00567     Realloc(p,EComplex,(_NP*(CE_Mka+CE_Mkb+CE_Mkc+3)));
00568     INFO("CE - CE_enumKV() : Number of Sincos Table Entries = "<<(_NP*(CE_Mka+CE_Mkb+CE_Mkc+3)));
00569     if(p==NULL) INFO("Error : CE_enumKV() -> Out of memory allocating sin/cos tables");
00570     Realloc(CE_scx,EComplex*,(CE_Mka+CE_Mkb+CE_Mkc+3)); // need to study more this part
00571     if(CE_scx==NULL) INFO("Error : CE_enumKV() -> Out of memory allocating sin/cos tables");
00572     CE_scy=CE_scx+(CE_Mka+1);
00573     CE_scz=CE_scy+(CE_Mkb+1);
00574     
00575     CE_scx[0]=p; 
00576     for(i=1;i<=CE_Mka;i++) CE_scx[i]=CE_scx[i-1]+_NP;
00577     CE_scy[0]=CE_scx[0]+(CE_Mka+1)*_NP;
00578     for(i=1;i<=CE_Mkb;i++) CE_scy[i]=CE_scy[i-1]+_NP;
00579     CE_scz[0]=CE_scy[0]+(CE_Mkb+1)*_NP;
00580     for(i=1;i<=CE_Mkc;i++) CE_scz[i]=CE_scz[i-1]+_NP;
00581 }
00582 
00583 
00584 //###################################################################
00585 //  PME Functions
00586 //###################################################################
00587 
00588 //-------------------------------------------------------------------
00589 // Warning
00590 // Note that if you want to switch conj_fixbox (from 1 to 0) or (from 
00591 // 0 to 1) in the middle of simulation, PME should be reinitialized.
00592 //-------------------------------------------------------------------
00593 
00594 //-------------------------------------------------------------------
00595 // Particle Mesh Ewald
00596 //-------------------------------------------------------------------
00597 void MDFrame::PME()
00598 {
00599     int i;
00600     CE_Real();
00601     PME_Rec();
00602 
00603     _EPOT_Ewald = _EPOT_Ewald_Real + _EPOT_Ewald_Rec + _EPOT_Ewald_Self;
00604     _VIRIAL_Ewald = _VIRIAL_Ewald_Real + _VIRIAL_Ewald_Rec;
00605 
00606     for(i=0;i<_NP;i++)
00607     {
00608         _F_Ewald[i] = _F_Ewald_Real[i] + _F_Ewald_Rec[i];
00609         _EPOT_IND_Ewald[i] = _EPOT_IND_Ewald_Real[i] + _EPOT_IND_Ewald_Rec[i];
00610     }
00611 }
00612 
00613 //-------------------------------------------------------------------
00614 //  Initialize Particle Mesh Ewald
00615 //-------------------------------------------------------------------
00616 void MDFrame::PME_init()
00617 {    
00618     int m1, m2, m3, tmp, tmp_index; // k-space indices
00619     ComplexType PME_BC_b1, PME_BC_b2, PME_BC_b3; // for the B matrix
00620     double tmp1, tmp2; // temporary dummy variable
00621 
00622     // k-space mesh setting
00623     PME_K23=PME_K2*PME_K3;
00624     PME_K123=PME_K1*PME_K23; // equal to the number of data points
00625 
00626     PME_K1d=((double)PME_K1);
00627     PME_K2d=((double)PME_K2);
00628     PME_K3d=((double)PME_K3);
00629     PME_K123d=((double)PME_K123);
00630     
00631     PME_K3S=(PME_K3/2)+1;
00632     PME_K23S=PME_K2*PME_K3S;
00633     PME_K123S=PME_K1*PME_K2*PME_K3S;
00634         
00635     // allocate array BC,Q,IQ,CONV, etc.
00636     Realloc(PME_B,double,PME_K123);
00637     Realloc(PME_C,double,PME_K123);
00638     Realloc(PME_BC,double,PME_K123);
00639     Realloc(PME_Q,double,(2*PME_K123S));
00640     Realloc(PME_IQ,ComplexType,PME_K123S);
00641     Realloc(PME_CONV,double,(2*PME_K123S));
00642     
00643     // scaled coefficients _UR
00644     Realloc(_UR,Vector3,_NP);
00645     
00646     // Mn(k+1) and factorial arrays
00647     Realloc(PME_bsp_Mk,double,(PME_bsp_n-1)); // caution : type double total (n-1) entries
00648     Realloc(PME_bsp_fac,double,(PME_bsp_n+1)); // caution : type integer total (n+1) entries
00649     
00650     // start and end indices of k-space mesh point
00651     Realloc(PME_m1s,int,_NP);
00652     Realloc(PME_m2s,int,_NP);
00653     Realloc(PME_m3s,int,_NP);
00654     
00655     // effective k-space indices
00656     tmp=PME_bsp_n*_NP;
00657     Realloc(PME_m1,int,PME_bsp_n);
00658     Realloc(PME_m2,int,PME_bsp_n);
00659     Realloc(PME_m3,int,PME_bsp_n);
00660     Realloc(PME_m1mod,int,tmp);
00661     Realloc(PME_m2mod,int,tmp);
00662     Realloc(PME_m3mod,int,tmp);
00663     
00664     // (u-k),Mn(u-k),Mn-1(u-k),Mn-1(u-k-1) arrays
00665     Realloc(PME_x1,double,tmp);
00666     Realloc(PME_x2,double,tmp);
00667     Realloc(PME_x3,double,tmp);
00668     Realloc(PME_MU1,double,tmp);
00669     Realloc(PME_MU2,double,tmp);
00670     Realloc(PME_MU3,double,tmp);
00671     Realloc(PME_d1MU1,double,PME_bsp_n);
00672     Realloc(PME_d1MU2,double,PME_bsp_n);
00673     Realloc(PME_d1MU3,double,PME_bsp_n);
00674     Realloc(PME_d2MU1,double,PME_bsp_n);
00675     Realloc(PME_d2MU2,double,PME_bsp_n);
00676     Realloc(PME_d2MU3,double,PME_bsp_n);
00677 
00678     // prepare data which will be used in the calculation
00679     PME_cal_bsp_fac();
00680     PME_cal_bsp_Mk();
00681     
00682     // calculate matrix B and C
00683     for (m1=0;m1<PME_K1;m1++) // from 0 to (K1-1)
00684     {
00685         for (m2=0;m2<PME_K2;m2++) // from 0 to (K2-1)
00686         {
00687             for (m3=0;m3<PME_K3;m3++) // from 0 to (K3-1)
00688             {
00689                 tmp_index=m3+m2*PME_K3+m1*PME_K23;
00690 
00691                 // calculate each B and C component
00692                 PME_cal_bsp_Bm(m1,PME_K1); PME_BC_b1[0]=PME_R_bsp_Bm[0];  PME_BC_b1[1]=PME_R_bsp_Bm[1];
00693                 PME_cal_bsp_Bm(m2,PME_K2); PME_BC_b2[0]=PME_R_bsp_Bm[0];  PME_BC_b2[1]=PME_R_bsp_Bm[1];
00694                 PME_cal_bsp_Bm(m3,PME_K3); PME_BC_b3[0]=PME_R_bsp_Bm[0];  PME_BC_b3[1]=PME_R_bsp_Bm[1];
00695                 tmp1=PME_cal_bsp_Cm(m1,m2,m3);
00696                 tmp2=(SQR(PME_BC_b1[0])+SQR(PME_BC_b1[1]))*(SQR(PME_BC_b2[0])+SQR(PME_BC_b2[1]))*(SQR(PME_BC_b3[0])+SQR(PME_BC_b3[1]));
00697                 
00698                 // compute BC matrix entries
00699                 PME_B[tmp_index]=tmp2;
00700                 PME_C[tmp_index]=tmp1;
00701                 PME_BC[tmp_index]=tmp1*tmp2;
00702             }
00703         }
00704     }
00705 
00706     // prepare fft plans
00707 #if _USEFFTW
00708     PME_fft_plan1=fftw_plan_dft_r2c_3d(PME_K1,PME_K2,PME_K3,PME_Q,PME_IQ,FFTW_ESTIMATE);
00709     PME_fft_plan2=fftw_plan_dft_c2r_3d(PME_K1,PME_K2,PME_K3,PME_IQ,PME_CONV,FFTW_ESTIMATE);
00710 #else
00711     PME_fft_plan1 = 0;
00712     PME_fft_plan2 = 0;
00713 #endif
00714     
00715 
00716     // Print Out PME Information
00717     INFO("------------------------------------------------------------------------------");
00718     INFO("------------------------------------------------------------------------------");
00719     INFO("            Particle Mesh Ewald Parameters [Long Range Interaction] ");
00720     INFO("------------------------------------------------------------------------------");
00721     INFO("    Real Space Cutoff Radius (Rc)                 = "<<Ewald_Rcut);
00722     INFO("    Ewald Gaussian Factor (Alpha)                 = "<<Ewald_Alpha);
00723     INFO("    Mesh Points in X-direction (PME_K1)           = "<<PME_K1);
00724     INFO("    Mesh Points in Y-direction (PME_K2)           = "<<PME_K2);
00725     INFO("    Mesh Points in Z-direction (PME_K3)           = "<<PME_K3);
00726     INFO("    Order of B-spline Interpolation (PME_bsp_n)   = "<<PME_bsp_n);
00727     INFO("------------------------------------------------------------------------------");
00728     INFO("------------------------------------------------------------------------------");
00729 }
00730 
00731 //-------------------------------------------------------------------
00732 //  Calculate Factorials
00733 //-------------------------------------------------------------------
00734 void MDFrame::PME_cal_bsp_fac()
00735 {
00736     int i, j, mult;
00737     
00738     // calculate factorials (0! to PME_bsp_n!)
00739     PME_bsp_fac[0]=1.0; // 0!=1
00740     for (i=1;i<=PME_bsp_n;i++) // from 1 to n
00741     {
00742         mult=1;
00743         for (j=1;j<=i;j++)
00744         {
00745             mult=mult*j;
00746         }
00747         PME_bsp_fac[i]=((double)mult);
00748         //INFO("[PME] Calculating Factorial => "<<i<<"! = "<<mult);
00749     }
00750 }
00751 
00752 //-------------------------------------------------------------------
00753 //  Calculate Mn(n,k+1)
00754 //-------------------------------------------------------------------
00755 void MDFrame::PME_cal_bsp_Mk()
00756 {
00757     int i;
00758     
00759     // calculate Mn(k+1)
00760     for (i=0;i<(PME_bsp_n-1);i++) // from 0 to n-2
00761     {
00762         PME_bsp_Mk[i]=PME_cal_bsp_Mu(PME_bsp_n,((double)(i+1))); // Mn(n,i+1) => from 1 to n-1
00763         //INFO("[PME] bsp_Mk["<<i<<"] = "<<PME_bsp_Mk[i]);
00764     }
00765 }
00766 
00767 //-------------------------------------------------------------------
00768 //  Calculate Mn(n,u)
00769 //-------------------------------------------------------------------
00770 double MDFrame::PME_cal_bsp_Mu(int n, double u)
00771 {
00772     int i,k;
00773     double bsp_sign, bsp_Mu, bsp_Muk, tmp;
00774 
00775     // calculate Mn(u-k)
00776     bsp_sign=1.0; // start from +1
00777     bsp_Mu=0.0; // initialize sum quantity
00778     for (k=0;k<=n;k++) // from 0 to n
00779     {
00780         tmp=u-((double)k); // convert k to double and calculate (u-k)
00781         if(tmp>0.0) // only when (u-k)>0
00782         {
00783             // calculate (u-k)^(n-1)
00784             bsp_Muk=1.0;
00785             for (i=1;i<n;i++) // from 1 to n-1
00786             {
00787                 bsp_Muk*=tmp;
00788             }
00789             // calculate n!/[k!(n-k)!]*(u-k)^(n-1) and sum it up
00790             bsp_Mu+=(bsp_sign*PME_bsp_fac[n]/(PME_bsp_fac[k]*PME_bsp_fac[n-k])*bsp_Muk);
00791             bsp_sign*=-1.0;
00792         }
00793     }
00794     bsp_Mu/=PME_bsp_fac[n-1]; // divide by (n-1)!
00795     return bsp_Mu;
00796 }
00797 
00798 //-------------------------------------------------------------------
00799 //  Calculate bi(mi)
00800 //-------------------------------------------------------------------
00801 void MDFrame::PME_cal_bsp_Bm(int mi, int Ki)
00802 {
00803     int k;
00804     double tmp;
00805     ComplexType bsp_bm, bsp_bms, bsp_bmn;
00806     
00807     // initialize the summation dummy complex number
00808     bsp_bms[0]=0.0;
00809     bsp_bms[1]=0.0;
00810     
00811     // calculate the B matrix coefficient Bi(mi)
00812     // compute denominator c+di
00813     for (k=0;k<(PME_bsp_n-1);k++) // from 0 to (n-2)
00814     {
00815         tmp=2.0*M_PI*(((double)(mi*k))/((double)Ki));
00816         bsp_bms[0]+=PME_bsp_Mk[k]*cos(tmp);
00817         bsp_bms[1]+=PME_bsp_Mk[k]*sin(tmp);
00818     }
00819     // compute numerator a+bi
00820     tmp=2.0*M_PI*(((double)((PME_bsp_n-1)*mi))/((double)Ki));
00821     bsp_bmn[0]=cos(tmp);
00822     bsp_bmn[1]=sin(tmp);
00823     
00824     // compute c^2+d^2
00825     tmp=SQR(bsp_bms[0])+SQR(bsp_bms[1]);
00826     
00827     // compute (a+bi)/(c+di) = [(ac+bd)+(bc-ad)i]/(c^2+d^2)
00828     bsp_bm[0]=(bsp_bmn[0]*bsp_bms[0]+bsp_bmn[1]*bsp_bms[1])/tmp; //(ac+bd)
00829     bsp_bm[1]=(bsp_bmn[1]*bsp_bms[0]-bsp_bmn[0]*bsp_bms[1])/tmp; //(bc-ad)
00830     PME_R_bsp_Bm[0]=bsp_bm[0];
00831     PME_R_bsp_Bm[1]=bsp_bm[1];
00832 }
00833 
00834 //-------------------------------------------------------------------
00835 //  Calculate C(m1,m2,m3)
00836 //-------------------------------------------------------------------
00837 double MDFrame::PME_cal_bsp_Cm(int m1, int m2, int m3)
00838 {
00839     double msq, m1s, m2s, m3s, tmp1, tmp2, bsp_Cm;
00840     
00841     if((m1==0)&&(m2==0)&&(m3==0))
00842     {
00843         bsp_Cm=0.0;
00844     }
00845     else
00846     {
00847         if(m1>(PME_K1/2)){m1=m1-PME_K1;}
00848         if(m2>(PME_K2/2)){m2=m2-PME_K2;}
00849         if(m3>(PME_K3/2)){m3=m3-PME_K3;}
00850         m1s=((double)m1)/Ewald_H_L1; // m1/L1 => double
00851         m2s=((double)m2)/Ewald_H_L2; // m2/L2 => double
00852         m3s=((double)m3)/Ewald_H_L3; // m3/L3 => double
00853         msq=SQR(m1s)+SQR(m2s)+SQR(m3s);
00854         tmp1=-SQR(M_PI)*msq/(SQR(Ewald_Alpha));
00855         tmp2=1.0/M_PI/Ewald_cell_V/msq;
00856         bsp_Cm=tmp2*exp(tmp1);
00857     }
00858     return bsp_Cm;
00859 }
00860 
00861 //-------------------------------------------------------------------
00862 //  Particle Mesh Ewald Summation
00863 //-------------------------------------------------------------------
00864 void MDFrame::PME_Rec()
00865 {
00866     int i, j, i1, i2, i3, m1, m2, m3, pm1, pm2, pm3;
00867     int tmp_ind_i, tmp_ind_1, tmp_ind_2, tmp_ind_3, tmp_index;
00868     double tmp_double, tmp_charge, m1s, m2s, m3s, msq, pfac, UE;
00869 
00870      // start time of reciprocal part
00871     PME_sREC=clock();
00872 
00873     // recalculate system size and volume (in case the box is changing)
00874     if(conj_fixbox==0)
00875       {
00876         Ewald_H_L1=sqrt(SQR(_H[0][0])+SQR(_H[1][0])+SQR(_H[2][0]));
00877         Ewald_H_L2=sqrt(SQR(_H[0][1])+SQR(_H[1][1])+SQR(_H[2][1]));
00878         Ewald_H_L3=sqrt(SQR(_H[0][2])+SQR(_H[1][2])+SQR(_H[2][2]));
00879         Ewald_cell_V=_H.det();
00880       }
00881 
00882     // initialize Q and recalculate C and BC if neccessary
00883     for (i1=0;i1<PME_K1;i1++) // from 0 to (K1-1)
00884       {
00885         for (i2=0;i2<PME_K2;i2++) // from 0 to (K2-1)
00886           {
00887             for (i3=0;i3<PME_K3;i3++) // from 0 to (K3-1)
00888               {
00889                 tmp_index=i3+i2*PME_K3+i1*PME_K23;
00890         
00891                 if(conj_fixbox==0)
00892                   {
00893                     // calculate each B and C component
00894                     tmp_double=PME_cal_bsp_Cm(i1,i2,i3);
00895                 
00896                     // compute BC matrix entries
00897                     PME_C[tmp_index]=tmp_double;
00898                     PME_BC[tmp_index]=PME_B[tmp_index]*tmp_double;
00899                   }
00900 
00901                 // initialize Q matrix
00902                 PME_Q[tmp_index]=0.0;
00903               }
00904           }
00905       }
00906 
00907     // initialize force, potential and virials
00908     _EPOT_Ewald_Rec = 0.0;
00909     _VIRIAL_Ewald_Rec.clear();
00910     for(i=0;i<_NP;i++)
00911       {
00912         _F_Ewald_Rec[i].clear();
00913         _EPOT_IND_Ewald_Rec[i]=0.0;
00914       }
00915 
00916     // get the scaled and manipulated coordinates UR
00917     for (i=0;i<_NP;i++)
00918       {
00919         _UR[i].x=(_SR[i].x+0.5)*PME_K1d;
00920         _UR[i].y=(_SR[i].y+0.5)*PME_K2d;
00921         _UR[i].z=(_SR[i].z+0.5)*PME_K3d;
00922       }
00923     
00924     // calculate Q matrix
00925     for (i=0;i<_NP;i++) // for all atoms
00926       {
00927         /* if fixed == -1, simply ignore this atom */
00928         if(fixed[i]==-1) continue;
00929         
00930         tmp_ind_i=PME_bsp_n*i;
00931         tmp_charge=P_SQRT_CLMB*_ATOMCHARGE[species[i]];
00932         
00933         // since function bsp_Mu has bounded support, we can find the initial and
00934         // final k-space indices of each atom where Q matrix entries are non-zero
00935         // [ start index : floor(u)-n+1    end index : floor(u) ]
00936         
00937         PME_m1s[i]=((int)floor(_UR[i].x))-PME_bsp_n+1;
00938         PME_m2s[i]=((int)floor(_UR[i].y))-PME_bsp_n+1;
00939         PME_m3s[i]=((int)floor(_UR[i].z))-PME_bsp_n+1;
00940         
00941         for (j=0;j<PME_bsp_n;j++) // from 0 to (n-1)
00942           {
00943             // calculate temporary index
00944             tmp_index =j+tmp_ind_i;
00945             
00946             // goes through effective range of k-space for atoms i
00947             PME_m1[j]=PME_m1s[i]+j; // in x-direction
00948             PME_m2[j]=PME_m2s[i]+j; // in y-direction
00949             PME_m3[j]=PME_m3s[i]+j; // in z-direction
00950             
00951             // apply periodic boundary condition
00952             PME_m1mod[tmp_index]=PME_int_mod(PME_m1[j],PME_K1);
00953             PME_m2mod[tmp_index]=PME_int_mod(PME_m2[j],PME_K2);
00954             PME_m3mod[tmp_index]=PME_int_mod(PME_m3[j],PME_K3);
00955             
00956             // calculate (u-k)
00957             PME_x1[tmp_index]=_UR[i].x-((double)PME_m1[j]);
00958             PME_x2[tmp_index]=_UR[i].y-((double)PME_m2[j]);
00959             PME_x3[tmp_index]=_UR[i].z-((double)PME_m3[j]);
00960             
00961             // calculate Mn(u) coefficients
00962             PME_MU1[tmp_index]=PME_cal_bsp_Mu(PME_bsp_n,PME_x1[tmp_index]);
00963             PME_MU2[tmp_index]=PME_cal_bsp_Mu(PME_bsp_n,PME_x2[tmp_index]);
00964             PME_MU3[tmp_index]=PME_cal_bsp_Mu(PME_bsp_n,PME_x3[tmp_index]);
00965           }
00966         
00967         // scan the atomwise small cubic which has nonzero entries
00968         for (i1=0;i1<PME_bsp_n;i1++) // from 0 to (n-1)
00969           {
00970             for (i2=0;i2<PME_bsp_n;i2++) // from 0 to (n-1)
00971               {
00972                 for (i3=0;i3<PME_bsp_n;i3++) // from 0 to (n-1)
00973                   {
00974                     tmp_ind_1=i1+tmp_ind_i;
00975                     tmp_ind_2=i2+tmp_ind_i;
00976                     tmp_ind_3=i3+tmp_ind_i;
00977                     tmp_index=PME_m3mod[tmp_ind_3]+PME_m2mod[tmp_ind_2]*PME_K3+PME_m1mod[tmp_ind_1]*PME_K23;
00978                     PME_Q[tmp_index]+=tmp_charge*PME_MU1[tmp_ind_1]*PME_MU2[tmp_ind_2]*PME_MU3[tmp_ind_3];
00979                   }
00980               }
00981           }
00982       }
00983     
00984     // out-of-place backward fft of matrix PME_Q (backward fft in the paper => forward fft for the fftw)
00985 #ifdef _USEFFTW
00986     fftw_execute(PME_fft_plan1); // real to complex
00987 #else
00988     INFO("FFTW not installed! Cannot use PME!");
00989 #endif
00990 
00991     // calculate energy and stress and
00992     // matrix multiplication of PME_BC and PME_IQ -> (IQ*BC)
00993     for (m1=0;m1<PME_K1;m1++) // from 0 to (PME_K1-1)
00994       {
00995         for (m2=0;m2<PME_K2;m2++) // from 0 to (PME_K2-1)
00996           {
00997             for (m3=0;m3<PME_K3S;m3++) // from 0 to (PME_K3S-1)
00998               {
00999                 tmp_ind_1=m3+m2*PME_K3S+m1*PME_K23S;
01000                 tmp_index=m3+m2*PME_K3+m1*PME_K23;
01001                 tmp_double=PME_BC[tmp_index]/PME_K123d;
01002                 UE=0.0;
01003                 
01004                 // calculate the reciprocal energy
01005                 if((m3==(PME_K3S-1))||(m3==0))
01006                   {UE=0.5*PME_BC[tmp_index]*(SQR(PME_IQ[tmp_ind_1][0])+SQR(PME_IQ[tmp_ind_1][1]));}
01007                 else
01008                   {UE=PME_BC[tmp_index]*(SQR(PME_IQ[tmp_ind_1][0])+SQR(PME_IQ[tmp_ind_1][1]));}
01009                 _EPOT_Ewald_Rec+=UE;
01010                 
01011                 // calculate BC*IQ
01012                 PME_IQ[tmp_ind_1][0]*=tmp_double;
01013                 PME_IQ[tmp_ind_1][1]*=tmp_double;
01014                 
01015                 // apply periodic condition
01016                 if(m1>(PME_K1/2)){pm1=m1-PME_K1;}
01017                 else{pm1=m1;}
01018                 if(m2>(PME_K2/2)){pm2=m2-PME_K2;}
01019                 else{pm2=m2;}
01020                 if(m3>(PME_K3/2)){pm3=m3-PME_K3;}
01021                 else{pm3=m3;}
01022 
01023                 //calculate the pre-factor
01024                 m1s=((double)pm1)/Ewald_H_L1; // m1/L1 => double
01025                 m2s=((double)pm2)/Ewald_H_L2; // m2/L2 => double
01026                 m3s=((double)pm3)/Ewald_H_L3; // m3/L3 => double
01027                 msq=SQR(m1s)+SQR(m2s)+SQR(m3s);
01028                 pfac=(1.0/msq+(SQR(M_PI))/(SQR(Ewald_Alpha)));
01029                 
01030                 // Virial Stress Calculation
01031                 if(!((m1==0)&&(m2==0)&&(m3==0)))
01032                   {
01033                     _VIRIAL_Ewald_Rec[0][0]+=(UE*(1.0-2.0*((m1s*m1s))*pfac)); //xx
01034                     _VIRIAL_Ewald_Rec[1][1]+=(UE*(1.0-2.0*((m2s*m2s))*pfac)); //yy
01035                     _VIRIAL_Ewald_Rec[2][2]+=(UE*(1.0-2.0*((m3s*m3s))*pfac)); //zz
01036                     _VIRIAL_Ewald_Rec[0][1]+=(UE*(-2.0*((m1s*m2s))*pfac));  //xy
01037                     _VIRIAL_Ewald_Rec[1][0]+=(UE*(-2.0*((m1s*m2s))*pfac));  //yx
01038                     _VIRIAL_Ewald_Rec[0][2]+=(UE*(-2.0*((m1s*m3s))*pfac));  //xz
01039                     _VIRIAL_Ewald_Rec[2][0]+=(UE*(-2.0*((m1s*m3s))*pfac));  //zx
01040                     _VIRIAL_Ewald_Rec[1][2]+=(UE*(-2.0*((m2s*m3s))*pfac));  //yz
01041                     _VIRIAL_Ewald_Rec[2][1]+=(UE*(-2.0*((m2s*m3s))*pfac));  //zy
01042                   }
01043               }
01044           }
01045       }
01046     
01047     // in-place forward fft of matrix PME_CONV (forward fft in the paper => backward fft for the fftw)
01048 #ifdef _USEFFTW
01049     fftw_execute(PME_fft_plan2); // complex to real
01050 #else
01051     INFO("FFTW not installed! Cannot use PME!");
01052 #endif
01053 
01054     /*
01055     // matrix multiplication of PME_Q and PME_CONV and calculate the reciprocal energy simultaneously
01056     for (m1=0;m1<PME_K1;m1++) // from 0 to (PME_K1-1)
01057     {
01058         for (m2=0;m2<PME_K2;m2++)  // from 0 to (PME_K2-1)
01059         {
01060             for (m3=0;m3<PME_K3;m3++) // from 0 to (PME_K3-1)
01061             {
01062                 tmp_index=m3+m2*PME_K3+m1*PME_K23;
01063                 _EPOT_Ewald_Rec+=(PME_Q[tmp_index]*PME_CONV[tmp_index]);
01064             }
01065         }
01066     }
01067     _EPOT_Ewald_Rec*=(PME_K123d/(2.0));
01068     //INFO("[PME] Reciprocal Energy Calculated by PME = "<<_EPOT_Ewald_Rec);
01069     */
01070     
01071     // force calculation - Calculate dQ/dr matrix
01072     for (i=0;i<_NP;i++) // for all atoms
01073       {
01074         /* if fixed == -1, simply ignore this atom */
01075         if(fixed[i]==-1) continue;
01076         
01077         tmp_ind_i=PME_bsp_n*i;
01078         for (j=0;j<PME_bsp_n;j++) // from 0 to (n-1)
01079           {
01080             // calculate temporary index
01081             tmp_index =j+tmp_ind_i;
01082             
01083             PME_d1MU1[j]=PME_cal_bsp_Mu((PME_bsp_n-1),PME_x1[tmp_index]); // Mn-1(u1-k1)
01084             PME_d1MU2[j]=PME_cal_bsp_Mu((PME_bsp_n-1),PME_x2[tmp_index]); // Mn-1(u2-k2)
01085             PME_d1MU3[j]=PME_cal_bsp_Mu((PME_bsp_n-1),PME_x3[tmp_index]); // Mn-1(u3-k3)
01086             
01087             PME_d2MU1[j]=PME_cal_bsp_Mu((PME_bsp_n-1),(PME_x1[tmp_index]-1.0)); // Mn-1(u1-k1-1)
01088             PME_d2MU2[j]=PME_cal_bsp_Mu((PME_bsp_n-1),(PME_x2[tmp_index]-1.0)); // Mn-1(u2-k2-1)
01089             PME_d2MU3[j]=PME_cal_bsp_Mu((PME_bsp_n-1),(PME_x3[tmp_index]-1.0)); // Mn-1(u3-k3-1)
01090           }
01091         
01092         // scan the atomwise small cubic which has nonzero entries
01093         for (i1=0;i1<PME_bsp_n;i1++) // from 0 to (n-1)
01094           {
01095             for (i2=0;i2<PME_bsp_n;i2++) // from 0 to (n-1)
01096               {
01097                 for (i3=0;i3<PME_bsp_n;i3++) // from 0 to (n-1)
01098                   {
01099                     tmp_ind_1=i1+tmp_ind_i;
01100                     tmp_ind_2=i2+tmp_ind_i;
01101                     tmp_ind_3=i3+tmp_ind_i;
01102                     tmp_index=PME_m3mod[tmp_ind_3]+PME_m2mod[tmp_ind_2]*PME_K3+PME_m1mod[tmp_ind_1]*PME_K23;
01103                     _F_Ewald_Rec[i].x+=((PME_d2MU1[i1]-PME_d1MU1[i1])*PME_MU2[tmp_ind_2]*PME_MU3[tmp_ind_3]*PME_CONV[tmp_index]);
01104                     _F_Ewald_Rec[i].y+=((PME_d2MU2[i2]-PME_d1MU2[i2])*PME_MU1[tmp_ind_1]*PME_MU3[tmp_ind_3]*PME_CONV[tmp_index]);
01105                     _F_Ewald_Rec[i].z+=((PME_d2MU3[i3]-PME_d1MU3[i3])*PME_MU1[tmp_ind_1]*PME_MU2[tmp_ind_2]*PME_CONV[tmp_index]);
01106                   }
01107               }
01108           }
01109         _F_Ewald_Rec[i].x*=(P_SQRT_CLMB*_ATOMCHARGE[species[i]]*PME_K123d*PME_K1d/Ewald_H_L1);
01110         _F_Ewald_Rec[i].y*=(P_SQRT_CLMB*_ATOMCHARGE[species[i]]*PME_K123d*PME_K2d/Ewald_H_L2);
01111         _F_Ewald_Rec[i].z*=(P_SQRT_CLMB*_ATOMCHARGE[species[i]]*PME_K123d*PME_K3d/Ewald_H_L3);
01112         //_F_Ewald_Rec[i].x*=PME_C1[species[i]];
01113         //_F_Ewald_Rec[i].y*=PME_C2[species[i]];
01114         //_F_Ewald_Rec[i].z*=PME_C3[species[i]];
01115       }
01116     
01117     // end time of reciprocal part
01118     PME_eREC=clock();
01119 }
01120 
01121 //-------------------------------------------------------------------
01122 //  Clear Particle Mesh Ewald
01123 //-------------------------------------------------------------------
01124 void MDFrame::PME_clear()
01125 {
01126     // free variables used in the Particle Mesh Ewald summation
01127     free(_UR); free(PME_bsp_Mk); free(PME_bsp_fac);
01128     
01129     free(PME_m1s); free(PME_m2s); free(PME_m3s);
01130     free(PME_m1); free(PME_m2); free(PME_m3);
01131     free(PME_m1mod); free(PME_m2mod); free(PME_m3mod);
01132     free(PME_x1); free(PME_x2); free(PME_x3);
01133     
01134     free(PME_MU1); free(PME_MU2); free(PME_MU3);
01135     free(PME_d1MU1); free(PME_d1MU2); free(PME_d1MU3);
01136     free(PME_d2MU1); free(PME_d2MU2); free(PME_d2MU3);
01137     
01138     free(PME_Q); free(PME_IQ);
01139     free(PME_BC); free(PME_CONV);
01140 
01141 #ifdef _USEFFTW    
01142     fftw_destroy_plan(PME_fft_plan1);
01143     fftw_destroy_plan(PME_fft_plan2);
01144 #endif
01145 }
01146