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

00001 /*
00002  * Last Modified : Wed Mar 12 22:40:16 2008
00003  *
00004  * Purpose:
00005  *   Numberical integrators for MD++, included in md.cpp
00006  *   different ensemble_type's and integrator_type's
00007  *
00008  * To Do:
00009  *
00010  *   1. HELMP not conserved in NPH Gear6, HELMP oscilliates in NPH VVerlet
00011  *
00012  *   2. NPT VVerlet need to be implemented
00013  */
00014 
00015 
00016 /*******************************************************
00017  *
00018  *  NVE ensemble
00019  *  Gear6 predictor-corrector integrator
00020  *
00021  ******************************************************/
00022 #define F02  (3./16.)
00023 #define F12  (251./360.)
00024 #define F32  (11./18.)
00025 #define F42  (1./6.)
00026 #define F52  (1./60.)
00027 
00028 void MDFrame::NVE_Gear6()
00029 {
00030     int i;
00031     double mass, tmp;
00032     Matrix33 hinv; Vector3 serr;
00033 
00034     /* initializing memory */    
00035     if(npold!=_NP) Gear6_init();    
00036 
00037     /* predictor */
00038     for(i=0;i<_NP;i++)
00039     {
00040         if(fixed[i]) continue;
00041         _SR[i]+=_VSR[i]+s2[i]+s3[i]+s4[i]+s5[i];
00042         _VSR[i]+=s2[i]*2.+s3[i]*3.+s4[i]*4.+s5[i]*5.;
00043         s2[i]+=s3[i]*3.+s4[i]*6.+s5[i]*10.;
00044         s3[i]+=s4[i]*4.+s5[i]*10.;
00045         s4[i]+=s5[i]*5.;
00046     }
00047 
00048     /* evaluator */
00049     call_potential();
00050     calcprop();
00051     
00052     hinv=_H.inv();
00053 
00054     // tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24 /(_ATOMMASS[0]*1e-3/AVO)*EV/1e-20;
00055     
00056     /* ATOMMASS has unit of g/mol
00057      * mass     has unit of eV ps^2 / A^2 
00058      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
00059     mass = _ATOMMASS[0]*MASSCONVERT;
00060     tmp  = 0.5 / mass * _TIMESTEP*_TIMESTEP;
00061         
00062     for(i=0;i<_NP;i++)
00063     {
00064         if(fixed[i]) continue;
00065         if(species[i]==0)
00066             s2real[i]=hinv*_F[i]*tmp;
00067         else
00068             s2real[i]=hinv*_F[i]*tmp *_ATOMMASS[0]/_ATOMMASS[species[i]];
00069     }
00070 
00071     /* corrector */
00072     for(i=0;i<_NP;i++)
00073     {
00074         if(fixed[i]) continue;
00075         serr=s2[i]-s2real[i];
00076         _SR[i]-=serr*F02;
00077         _VSR[i]-=serr*F12;
00078         s2[i]-=serr;
00079         s3[i]-=serr*F32;
00080         s4[i]-=serr*F42;
00081         s5[i]-=serr*F52;
00082     }
00083 }
00084 
00085 /*******************************************************
00086  *
00087  *  NVT ensemble
00088  *  Gear6 predictor-corrector integrator
00089  *
00090  ******************************************************/
00091 void MDFrame::NVT_Gear6()
00092 {
00093     int i;
00094     double mass, tmp, zetaerr;
00095     Matrix33 hinv; Vector3 serr;    
00096 
00097     /* initializing memory */    
00098     if(npold!=_NP) Gear6_init();    
00099     
00100      /* predictor */
00101     for(i=0;i<_NP;i++)
00102     {
00103         if(fixed[i]) continue;
00104         _SR[i]+=_VSR[i]+s2[i]+s3[i]+s4[i]+s5[i];
00105         _VSR[i]+=s2[i]*2.+s3[i]*3.+s4[i]*4.+s5[i]*5.;
00106         s2[i]+=s3[i]*3.+s4[i]*6.+s5[i]*10.;
00107         s3[i]+=s4[i]*4.+s5[i]*10.;
00108         s4[i]+=s5[i]*5.;
00109     }
00110 
00111     zeta+=zetav+zeta2+zeta3+zeta4+zeta5;
00112     zetav+=zeta2*2.+zeta3*3.+zeta4*4.+zeta5*5.;
00113     zeta2+=zeta3*3.+zeta4*6.+zeta5*10.;
00114     zeta3+=zeta4*4.+zeta5*10.;
00115     zeta4+=zeta5*5.;
00116     
00117     /* evaluator */
00118     call_potential();
00119     calcprop();
00120     
00121     hinv=_H.inv();
00122 
00123     // tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24 /(_ATOMMASS[0]*1e-3/AVO)*EV/1e-20;
00124     
00125     /* ATOMMASS has unit of g/mol
00126      * mass     has unit of eV ps^2 / A^2 
00127      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
00128     mass = _ATOMMASS[0]*MASSCONVERT;
00129     tmp  = 0.5 / mass * _TIMESTEP*_TIMESTEP;
00130     
00131     for(i=0;i<_NP;i++)
00132     {
00133         if(fixed[i]) continue;
00134         if(species[i]==0)
00135             s2real[i]=hinv*_F[i]*tmp;
00136         else
00137             s2real[i]=hinv*_F[i]*tmp *_ATOMMASS[0]/_ATOMMASS[species[i]];
00138     }
00139 
00140     if(NHMass[0]==0)
00141     { /* old input format */
00142         zetaa=vt2*(_T/_TDES-1)*_TIMESTEP*_TIMESTEP*1e-24*0.5;
00143     }
00144     else
00145     { /* new input format: Q = NHMass = (3NKT)*1e24/vt2 (in unit of eV*ps^2)
00146        *  T = 300K KT = 0.025 N = 1000 3NKT = 75 --> vt2=1e28 <=> NHMass=7.5e-3)
00147        */
00148         zetaa = (_KATOM*2 - 3*_NPfree*KB*_TDES) / NHMass[0] * (0.5*_TIMESTEP*_TIMESTEP);
00149     }
00150         
00151 
00152     for(i=0;i<_NP;i++)
00153     {
00154         if(fixed[i]) continue;
00155         s2real[i]-=_VSR[i]*zetav*0.5;
00156     }
00157 
00158     /* corrector */
00159     for(i=0;i<_NP;i++)
00160     {
00161         if(fixed[i]) continue;
00162         serr=s2[i]-s2real[i];
00163         _SR[i]-=serr*F02;
00164         _VSR[i]-=serr*F12;
00165         s2[i]-=serr;
00166         s3[i]-=serr*F32;
00167         s4[i]-=serr*F42;
00168         s5[i]-=serr*F52;
00169     }
00170     zetaerr=zeta2-zetaa;
00171     zeta-=zetaerr*F02;
00172     zetav-=zetaerr*F12;
00173     zeta2-=zetaerr;
00174     zeta3-=zetaerr*F32;
00175     zeta4-=zetaerr*F42;
00176     zeta5-=zetaerr*F52;
00177     
00178 }
00179 
00180 /*******************************************************
00181  *
00182  *  NVT ensemble (Nose Hoover Chain)
00183  *  Gear6 predictor-corrector integrator
00184  *
00185  ******************************************************/
00186 void MDFrame::NVTC_Gear6()
00187 {
00188     int i;
00189     double mass, tmp, zetaNHCerr[MAXNHCLEN];
00190     Matrix33 hinv; Vector3 serr;    
00191 
00192     //INFO_Printf("NVHC_Gear6: Hello!!!!  NHChainLen = %d\n",NHChainLen);
00193     
00194     /* initializing memory */
00195     if(npold!=_NP) Gear6_init();
00196     
00197      /* predictor */
00198     for(i=0;i<_NP;i++)
00199     {
00200         if(fixed[i]) continue;
00201         _SR[i]+=_VSR[i]+s2[i]+s3[i]+s4[i]+s5[i];
00202         _VSR[i]+=s2[i]*2.+s3[i]*3.+s4[i]*4.+s5[i]*5.;
00203         s2[i]+=s3[i]*3.+s4[i]*6.+s5[i]*10.;
00204         s3[i]+=s4[i]*4.+s5[i]*10.;
00205         s4[i]+=s5[i]*5.;
00206     }
00207 
00208     for(i=0;i<NHChainLen;i++)
00209     {
00210         zetaNHC[i] +=zetaNHCv[i]+zetaNHC2[i]+zetaNHC3[i]+zetaNHC4[i]+zetaNHC5[i];
00211         zetaNHCv[i]+=zetaNHC2[i]*2.+zetaNHC3[i]*3.+zetaNHC4[i]*4.+zetaNHC5[i]*5.;
00212         zetaNHC2[i]+=zetaNHC3[i]*3.+zetaNHC4[i]*6.+zetaNHC5[i]*10.;
00213         zetaNHC3[i]+=zetaNHC4[i]*4.+zetaNHC5[i]*10.;
00214         zetaNHC4[i]+=zetaNHC5[i]*5.;
00215     }
00216 
00217     /* evaluator */
00218     call_potential();
00219     calcprop();
00220     
00221     hinv=_H.inv();
00222 
00223     // tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24 /(_ATOMMASS[0]*1e-3/AVO)*EV/1e-20;     //(A^2)
00224 
00225     /* ATOMMASS has unit of g/mol
00226      * mass     has unit of eV ps^2 / A^2 
00227      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
00228     mass = _ATOMMASS[0]*MASSCONVERT;
00229     tmp  = 0.5 / mass * _TIMESTEP*_TIMESTEP;
00230     
00231     
00232     for(i=0;i<_NP;i++)
00233     {
00234         if(fixed[i]) continue;
00235         if(species[i]==0)
00236             s2real[i]=hinv*_F[i]*tmp;
00237         else
00238             s2real[i]=hinv*_F[i]*tmp *_ATOMMASS[0]/_ATOMMASS[species[i]];
00239     }
00240 
00241     //INFO_Printf("NVHC: NHChainLen = %d\n",NHChainLen);
00242     
00243     if(NHChainLen<=1)
00244     {
00245         /* new input format: Q = NHMass = (3NKT)*1e24/vt2 (in unit of eV*ps^2)
00246          *  T = 300K KT = 0.025 N = 1000 3NKT = 75 --> vt2=1e28 <=> NHMass=7.5e-3)
00247          */
00248         zetaNHCa[0] = (_KATOM*2 - 3*_NPfree*KB*_TDES) / NHMass[0] * (0.5*_TIMESTEP*_TIMESTEP);
00249     }
00250     else
00251     {    
00252         /* new input format: Q = NHMass = (3NKT)*1e24/vt2 (in unit of eV*ps^2)
00253          *  T = 300K KT = 0.025 N = 1000 3NKT = 75 --> vt2=1e28 <=> NHMass=7.5e-3)
00254          */
00255         zetaNHCa[0] = (_KATOM*2 - 3*_NPfree*KB*_TDES) / NHMass[0] * (0.5*_TIMESTEP*_TIMESTEP)
00256                       - zetaNHCv[0]*zetaNHCv[1] * 0.5;
00257                          
00258         for(i=1;i<NHChainLen-1;i++)
00259         {
00260             zetaNHCa[i] = (NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES) / NHMass[i] * (0.5*_TIMESTEP*_TIMESTEP)
00261                           -  zetaNHCv[i]*zetaNHCv[i+1] * 0.5;
00262         }
00263         i = NHChainLen-1;
00264         zetaNHCa[i] = (NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES) / NHMass[i] * (0.5*_TIMESTEP*_TIMESTEP);       
00265     }
00266     
00267     for(i=0;i<_NP;i++)
00268     {
00269         if(fixed[i]) continue;
00270         s2real[i]-=_VSR[i]*zetaNHCv[0]*0.5;
00271     }
00272 
00273     /* corrector */
00274     for(i=0;i<_NP;i++)
00275     {
00276         if(fixed[i]) continue;
00277         serr=s2[i]-s2real[i];
00278         _SR[i]-=serr*F02;
00279         _VSR[i]-=serr*F12;
00280         s2[i]-=serr;
00281         s3[i]-=serr*F32;
00282         s4[i]-=serr*F42;
00283         s5[i]-=serr*F52;
00284     }
00285 
00286 
00287     for(i=0;i<NHChainLen;i++)
00288     {
00289         zetaNHCerr[i] = zetaNHC2[i]-zetaNHCa[i];
00290     }
00291 
00292     for(i=0;i<NHChainLen;i++)
00293     {
00294         zetaNHC[i]  -= zetaNHCerr[i]*F02;
00295         zetaNHCv[i] -= zetaNHCerr[i]*F12;
00296         zetaNHC2[i] -= zetaNHCerr[i];
00297         zetaNHC3[i] -= zetaNHCerr[i]*F32;
00298         zetaNHC4[i] -= zetaNHCerr[i]*F42;
00299         zetaNHC5[i] -= zetaNHCerr[i]*F52;        
00300     }
00301 
00302 }
00303 
00304 /*******************************************************
00305  *
00306  *  NPH ensemble
00307  *  Gear6 predictor-corrector integrator
00308  *
00309  ******************************************************/
00310 void MDFrame::NPH_Gear6()
00311 {
00312     int i, j;
00313     double mass, tmp;
00314     Matrix33 hinv, htran, vhtran, G, Ginv, VG, VG1, VG2, GiVG, herr;
00315     Vector3 serr;    
00316 
00317     /* initializing memory */    
00318     if(npold!=_NP) Gear6_init();    
00319     
00320      /* predictor */
00321     for(i=0;i<_NP;i++)
00322     {
00323         if(fixed[i]) continue;
00324         _SR[i]+=_VSR[i]+s2[i]+s3[i]+s4[i]+s5[i];
00325         _VSR[i]+=s2[i]*2.+s3[i]*3.+s4[i]*4.+s5[i]*5.;
00326         s2[i]+=s3[i]*3.+s4[i]*6.+s5[i]*10.;
00327         s3[i]+=s4[i]*4.+s5[i]*10.;
00328         s4[i]+=s5[i]*5.;
00329     }
00330     for(i=0;i<3;i++)
00331         for(j=0;j<3;j++)
00332             if(conj_fixboxvec[i][j]==1)
00333                 _VH[i][j]=h2[i][j]=h3[i][j]=h4[i][j]=h5[i][j]=0;
00334 
00335     _H+=_VH+h2+h3+h4+h5;
00336     _VH+= h2*2.+ h3*3.+ h4*4.+ h5*5.;
00337     h2+= h3*3.+ h4*6.+ h5*10.;
00338     h3+= h4*4.+ h5*10.;
00339     h4+= h5*5.;
00340 
00341     /* evaluator */
00342     call_potential();
00343     calcprop();
00344     
00345     hinv=_H.inv();
00346 
00347     // tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24 /(_ATOMMASS[0]*1e-3/AVO)*EV/1e-20;     //(A^2)
00348 
00349     /* ATOMMASS has unit of g/mol
00350      * mass     has unit of eV ps^2 / A^2 
00351      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
00352     mass = _ATOMMASS[0]*MASSCONVERT;
00353     tmp  = 0.5 / mass * _TIMESTEP*_TIMESTEP;
00354     
00355     
00356     for(i=0;i<_NP;i++)
00357     {
00358         if(fixed[i]) continue;
00359         if(species[i]==0)
00360             s2real[i]=hinv*_F[i]*tmp;
00361         else
00362             s2real[i]=hinv*_F[i]*tmp *_ATOMMASS[0]/_ATOMMASS[species[i]];
00363     }
00364     h2real=_GH* tmp * _ATOMMASS[0]/_WALLMASS;
00365 
00366     if(_BOXDAMP!=0)
00367         for(i=0;i<3;i++)
00368             for(j=0;j<3;j++)
00369                 if(conj_fixboxvec[i][j]==0)
00370                     h2real[i][j]-=_VH[i][j]*_BOXDAMP;
00371 
00372     htran=_H.tran(); G=htran*_H; Ginv=G.inv();
00373     vhtran=_VH.tran(); VG1=vhtran*_H; VG2=htran*_VH;
00374     VG=VG1+VG2; GiVG=Ginv*VG;
00375 
00376     for(i=0;i<_NP;i++)
00377         s2real[i]-=GiVG*_VSR[i];
00378 
00379     /* corrector */
00380     for(i=0;i<_NP;i++)
00381     {
00382         if(fixed[i]) continue;
00383         serr=s2[i]-s2real[i];
00384         _SR[i]-=serr*F02;
00385         _VSR[i]-=serr*F12;
00386         s2[i]-=serr;
00387         s3[i]-=serr*F32;
00388         s4[i]-=serr*F42;
00389         s5[i]-=serr*F52;
00390     }    
00391 
00392     /* Parrinello-Rahman boundary condition */
00393     herr=h2-h2real;
00394     for(i=0;i<3;i++)
00395         for(j=0;j<3;j++)
00396             if(conj_fixboxvec[i][j]==0)
00397             {
00398                  _H[i][j] -= herr[i][j]*F02;
00399                 _VH[i][j] -= herr[i][j]*F12;
00400                  h2[i][j] -= herr[i][j];
00401                  h3[i][j] -= herr[i][j]*F32;
00402                  h4[i][j] -= herr[i][j]*F42;
00403                  h5[i][j] -= herr[i][j]*F52;
00404             }
00405 
00406 }
00407 
00408 
00409 /*******************************************************
00410  *
00411  *  NPT ensemble
00412  *  Gear6 predictor-corrector integrator
00413  *
00414  ******************************************************/
00415 void MDFrame::NPT_Gear6()
00416 {
00417     int i, j;
00418     double mass, tmp, zetaerr;
00419     Matrix33 hinv, htran, vhtran, G, Ginv, VG, VG1, VG2, GiVG, herr;
00420     Vector3 serr;    
00421 
00422     /* initializing memory */    
00423     if(npold!=_NP) Gear6_init();    
00424     
00425      /* predictor */
00426     for(i=0;i<_NP;i++)
00427     {
00428         if(fixed[i]) continue;
00429         _SR[i]+=_VSR[i]+s2[i]+s3[i]+s4[i]+s5[i];
00430         _VSR[i]+=s2[i]*2.+s3[i]*3.+s4[i]*4.+s5[i]*5.;
00431         s2[i]+=s3[i]*3.+s4[i]*6.+s5[i]*10.;
00432         s3[i]+=s4[i]*4.+s5[i]*10.;
00433         s4[i]+=s5[i]*5.;
00434     }
00435     for(i=0;i<3;i++)
00436         for(j=0;j<3;j++)
00437             if(conj_fixboxvec[i][j]==1)
00438                 _VH[i][j]=h2[i][j]=h3[i][j]=h4[i][j]=h5[i][j]=0;
00439 
00440     zeta+=zetav+zeta2+zeta3+zeta4+zeta5;
00441     zetav+=zeta2*2.+zeta3*3.+zeta4*4.+zeta5*5.;
00442     zeta2+=zeta3*3.+zeta4*6.+zeta5*10.;
00443     zeta3+=zeta4*4.+zeta5*10.;
00444     zeta4+=zeta5*5.;
00445     
00446     _H+=_VH+h2+h3+h4+h5;
00447     _VH+= h2*2.+ h3*3.+ h4*4.+ h5*5.;
00448     h2+= h3*3.+ h4*6.+ h5*10.;
00449     h3+= h4*4.+ h5*10.;
00450     h4+= h5*5.;
00451     
00452     /* evaluator */
00453     call_potential();
00454     calcprop();
00455     
00456     hinv=_H.inv();
00457 
00458     // tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24 /(_ATOMMASS[0]*1e-3/AVO)*EV/1e-20;     //(A^2)
00459 
00460     /* ATOMMASS has unit of g/mol
00461      * mass     has unit of eV ps^2 / A^2 
00462      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
00463     mass = _ATOMMASS[0]*MASSCONVERT;
00464     tmp  = 0.5 / mass * _TIMESTEP*_TIMESTEP;
00465     
00466     
00467     for(i=0;i<_NP;i++)
00468     {
00469         if(fixed[i]) continue;
00470         if(species[i]==0)
00471             s2real[i]=hinv*_F[i]*tmp;
00472         else
00473             s2real[i]=hinv*_F[i]*tmp *_ATOMMASS[0]/_ATOMMASS[species[i]];
00474     }
00475 
00476     zetaa=vt2*(_T/_TDES-1)*_TIMESTEP*_TIMESTEP*1e-24*0.5;
00477 
00478     for(i=0;i<_NP;i++)
00479     {
00480         if(fixed[i]) continue;
00481         s2real[i]-=_VSR[i]*zetav*0.5;
00482     }
00483 
00484     h2real=_GH* tmp * _ATOMMASS[0]/_WALLMASS;
00485     
00486     if(_BOXDAMP!=0)
00487         for(i=0;i<3;i++)
00488             for(j=0;j<3;j++)
00489                 if(conj_fixboxvec[i][j]==0)
00490                     h2real[i][j]-=_VH[i][j]*_BOXDAMP;
00491 
00492     htran=_H.tran(); G=htran*_H; Ginv=G.inv();
00493     vhtran=_VH.tran(); VG1=vhtran*_H; VG2=htran*_VH;
00494     VG=VG1+VG2; GiVG=Ginv*VG;
00495 
00496     for(i=0;i<_NP;i++)
00497         s2real[i]-=GiVG*_VSR[i];
00498     
00499     /* corrector */
00500     for(i=0;i<_NP;i++)
00501     {
00502         if(fixed[i]) continue;
00503         serr=s2[i]-s2real[i];
00504         _SR[i]-=serr*F02;
00505         _VSR[i]-=serr*F12;
00506         s2[i]-=serr;
00507         s3[i]-=serr*F32;
00508         s4[i]-=serr*F42;
00509         s5[i]-=serr*F52;
00510     }
00511     zetaerr=zeta2-zetaa;
00512     zeta-=zetaerr*F02;
00513     zetav-=zetaerr*F12;
00514     zeta2-=zetaerr;
00515     zeta3-=zetaerr*F32;
00516     zeta4-=zetaerr*F42;
00517     zeta5-=zetaerr*F52;
00518 
00519     herr=h2-h2real;
00520     for(i=0;i<3;i++)
00521         for(j=0;j<3;j++)
00522             if(conj_fixboxvec[i][j]==0)
00523             {
00524                  _H[i][j] -= herr[i][j]*F02;
00525                 _VH[i][j] -= herr[i][j]*F12;
00526                  h2[i][j] -= herr[i][j];
00527                  h3[i][j] -= herr[i][j]*F32;
00528                  h4[i][j] -= herr[i][j]*F42;
00529                  h5[i][j] -= herr[i][j]*F52;
00530             }        
00531 }
00532 
00533 /*******************************************************
00534  *
00535  *  NPT ensemble (Nose Hoover Chain)
00536  *  Gear6 predictor-corrector integrator
00537  *
00538  ******************************************************/
00539 void MDFrame::NPTC_Gear6()
00540 {
00541     int i, j; //,boxdof;
00542     double mass, tmp, zetaNHCerr[MAXNHCLEN]; //, zetaBNHCerr[MAXNHCLEN];
00543     Matrix33 hinv, htran, vhtran, G, Ginv, VG, VG1, VG2, GiVG, herr;
00544     Vector3 serr;    
00545 
00546     /* initializing memory */    
00547     if(npold!=_NP) Gear6_init();    
00548     
00549      /* predictor */
00550     for(i=0;i<_NP;i++)
00551     {
00552         if(fixed[i]) continue;
00553         _SR[i]+=_VSR[i]+s2[i]+s3[i]+s4[i]+s5[i];
00554         _VSR[i]+=s2[i]*2.+s3[i]*3.+s4[i]*4.+s5[i]*5.;
00555         s2[i]+=s3[i]*3.+s4[i]*6.+s5[i]*10.;
00556         s3[i]+=s4[i]*4.+s5[i]*10.;
00557         s4[i]+=s5[i]*5.;
00558     }
00559     for(i=0;i<3;i++)
00560         for(j=0;j<3;j++)
00561             if(conj_fixboxvec[i][j]==1)
00562                 _VH[i][j]=h2[i][j]=h3[i][j]=h4[i][j]=h5[i][j]=0;
00563 
00564     /* Nose Hoover Chain */
00565     for(i=0;i<NHChainLen;i++)
00566     {
00567         zetaNHC[i] +=zetaNHCv[i]+zetaNHC2[i]+zetaNHC3[i]+zetaNHC4[i]+zetaNHC5[i];
00568         zetaNHCv[i]+=zetaNHC2[i]*2.+zetaNHC3[i]*3.+zetaNHC4[i]*4.+zetaNHC5[i]*5.;
00569         zetaNHC2[i]+=zetaNHC3[i]*3.+zetaNHC4[i]*6.+zetaNHC5[i]*10.;
00570         zetaNHC3[i]+=zetaNHC4[i]*4.+zetaNHC5[i]*10.;
00571         zetaNHC4[i]+=zetaNHC5[i]*5.;
00572     }
00573 
00574     /* Nose Hoover Chain for Box */
00575 //    for(i=0;i<NHChainLen;i++)
00576 //    {
00577 //        zetaBNHC[i] +=zetaBNHCv[i]+zetaBNHC2[i]+zetaBNHC3[i]+zetaBNHC4[i]+zetaBNHC5[i];
00578 //        zetaBNHCv[i]+=zetaBNHC2[i]*2.+zetaBNHC3[i]*3.+zetaBNHC4[i]*4.+zetaBNHC5[i]*5.;
00579 //        zetaBNHC2[i]+=zetaBNHC3[i]*3.+zetaBNHC4[i]*6.+zetaBNHC5[i]*10.;
00580 //        zetaBNHC3[i]+=zetaBNHC4[i]*4.+zetaBNHC5[i]*10.;
00581 //        zetaBNHC4[i]+=zetaBNHC5[i]*5.;
00582 //    }
00583     
00584     
00585     _H+=_VH+h2+h3+h4+h5;
00586     _VH+= h2*2.+ h3*3.+ h4*4.+ h5*5.;
00587     h2+= h3*3.+ h4*6.+ h5*10.;
00588     h3+= h4*4.+ h5*10.;
00589     h4+= h5*5.;
00590     
00591     /* evaluator */
00592     call_potential();
00593     calcprop();
00594     
00595     hinv=_H.inv();
00596 
00597     // tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24 /(_ATOMMASS[0]*1e-3/AVO)*EV/1e-20;     //(A^2)
00598 
00599     /* ATOMMASS has unit of g/mol
00600      * mass     has unit of eV ps^2 / A^2 
00601      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
00602     mass = _ATOMMASS[0]*MASSCONVERT;
00603     tmp  = 0.5 / mass * _TIMESTEP*_TIMESTEP;
00604     
00605     
00606     for(i=0;i<_NP;i++)
00607     {
00608         if(fixed[i]) continue;
00609         if(species[i]==0)
00610             s2real[i]=hinv*_F[i]*tmp;
00611         else
00612             s2real[i]=hinv*_F[i]*tmp *_ATOMMASS[0]/_ATOMMASS[species[i]];
00613     }
00614 
00615     if(NHChainLen<=1)
00616     {
00617         /* new input format: Q = NHMass = (3NKT)*1e24/vt2 (in unit of eV*ps^2)
00618          *  T = 300K KT = 0.025 N = 1000 3NKT = 75 --> vt2=1e28 <=> NHMass=7.5e-3)
00619          */
00620         zetaNHCa[0] = (_KATOM*2 - 3*_NPfree*KB*_TDES) / NHMass[0] * (0.5*_TIMESTEP*_TIMESTEP);
00621 
00622 //        boxdof = 9; for(i=0;i<3;i++) for(j=0;j<3;j++) if(conj_fixboxvec[i][j]) boxdof--;
00623 //        zetaBNHCa[0] = (_KBOX*2 - boxdof*KB*_TDES) / BNHMass[0] * (0.5*_TIMESTEP*_TIMESTEP);        
00624     }
00625     else
00626     {    
00627         /* new input format: Q = NHMass = (3NKT)*1e24/vt2 (in unit of eV*ps^2)
00628          *  T = 300K KT = 0.025 N = 1000 3NKT = 75 --> vt2=1e28 <=> NHMass=7.5e-3)
00629          */
00630         zetaNHCa[0] = (_KATOM*2 - 3*_NPfree*KB*_TDES) / NHMass[0] * (0.5*_TIMESTEP*_TIMESTEP)
00631                       - zetaNHCv[0]*zetaNHCv[1] * 0.5;
00632 
00633 //        INFO_Printf("KBOX = %20.12e\n",_KBOX);
00634 //        boxdof = 9; for(i=0;i<3;i++) for(j=0;j<3;j++) if(conj_fixboxvec[i][j]) boxdof--;
00635 //        zetaBNHCa[0] = (_KBOX*2 - boxdof*KB*_TDES) / BNHMass[0] * (0.5*_TIMESTEP*_TIMESTEP)
00636 //                      - zetaBNHCv[0]*zetaBNHCv[1] * 0.5;
00637                 
00638         
00639         for(i=1;i<NHChainLen-1;i++)
00640         {
00641             zetaNHCa[i] = (NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES) / NHMass[i] * (0.5*_TIMESTEP*_TIMESTEP)
00642                           -  zetaNHCv[i]*zetaNHCv[i+1] * 0.5;
00643 //            zetaBNHCa[i] = (BNHMass[i-1]*SQR(zetaBNHCv[i-1]/_TIMESTEP)-KB*_TDES) / BNHMass[i] * (0.5*_TIMESTEP*_TIMESTEP)
00644 //                          -  zetaBNHCv[i]*zetaBNHCv[i+1] * 0.5;            
00645         }
00646         i = NHChainLen-1;
00647         zetaNHCa[i] = (NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES) / NHMass[i] * (0.5*_TIMESTEP*_TIMESTEP);       
00648 //        zetaBNHCa[i] = (BNHMass[i-1]*SQR(zetaBNHCv[i-1]/_TIMESTEP)-KB*_TDES) / BNHMass[i] * (0.5*_TIMESTEP*_TIMESTEP);       
00649     }
00650 
00651     for(i=0;i<_NP;i++)
00652     {
00653         if(fixed[i]) continue;
00654         s2real[i]-=_VSR[i]*zetaNHCv[0]*0.5;
00655     }
00656 
00657     h2real=_GH * tmp * _ATOMMASS[0]/_WALLMASS;
00658 
00659 //    for(i=0;i<3;i++)
00660 //        for(j=0;j<3;j++)
00661 //            if(conj_fixboxvec[i][j]==0)
00662 //                h2real[i][j]-=_VH[i][j]*zetaBNHCv[0];
00663     
00664     if(_BOXDAMP!=0)
00665         for(i=0;i<3;i++)
00666             for(j=0;j<3;j++)
00667                 if(conj_fixboxvec[i][j]==0)
00668                     h2real[i][j]-=_VH[i][j]*_BOXDAMP;
00669 
00670     htran=_H.tran(); G=htran*_H; Ginv=G.inv();
00671     vhtran=_VH.tran(); VG1=vhtran*_H; VG2=htran*_VH;
00672     VG=VG1+VG2; GiVG=Ginv*VG;
00673 
00674     for(i=0;i<_NP;i++)
00675         s2real[i]-=GiVG*_VSR[i];
00676     
00677     /* corrector */
00678     for(i=0;i<_NP;i++)
00679     {
00680         if(fixed[i]) continue;
00681         serr=s2[i]-s2real[i];
00682         _SR[i]-=serr*F02;
00683         _VSR[i]-=serr*F12;
00684         s2[i]-=serr;
00685         s3[i]-=serr*F32;
00686         s4[i]-=serr*F42;
00687         s5[i]-=serr*F52;
00688     }
00689 
00690     for(i=0;i<NHChainLen;i++)
00691     {
00692         zetaNHCerr[i] = zetaNHC2[i]-zetaNHCa[i];
00693 //        zetaBNHCerr[i] = zetaBNHC2[i]-zetaBNHCa[i];
00694     }
00695 
00696     for(i=0;i<NHChainLen;i++)
00697     {
00698         zetaNHC[i]  -= zetaNHCerr[i]*F02;
00699         zetaNHCv[i] -= zetaNHCerr[i]*F12;
00700         zetaNHC2[i] -= zetaNHCerr[i];
00701         zetaNHC3[i] -= zetaNHCerr[i]*F32;
00702         zetaNHC4[i] -= zetaNHCerr[i]*F42;
00703         zetaNHC5[i] -= zetaNHCerr[i]*F52;
00704 
00705 //        zetaBNHC[i]  -= zetaBNHCerr[i]*F02;
00706 //        zetaBNHCv[i] -= zetaBNHCerr[i]*F12;
00707 //        zetaBNHC2[i] -= zetaBNHCerr[i];
00708 //        zetaBNHC3[i] -= zetaBNHCerr[i]*F32;
00709 //        zetaBNHC4[i] -= zetaBNHCerr[i]*F42;
00710 //        zetaBNHC5[i] -= zetaBNHCerr[i]*F52;                
00711     }
00712     
00713 
00714     herr=h2-h2real;
00715     for(i=0;i<3;i++)
00716         for(j=0;j<3;j++)
00717             if(conj_fixboxvec[i][j]==0)
00718             {
00719                  _H[i][j] -= herr[i][j]*F02;
00720                 _VH[i][j] -= herr[i][j]*F12;
00721                  h2[i][j] -= herr[i][j];
00722                  h3[i][j] -= herr[i][j]*F32;
00723                  h4[i][j] -= herr[i][j]*F42;
00724                  h5[i][j] -= herr[i][j]*F52;
00725             }        
00726 }
00727 
00728 
00729 
00730 /* Gear6 initializing memory */   
00731 void MDFrame::Gear6_init()
00732 {
00733     Realloc(s2,Vector3,_NP);
00734     Realloc(s3,Vector3,_NP);
00735     Realloc(s4,Vector3,_NP);
00736     Realloc(s5,Vector3,_NP);
00737     Realloc(s2real,Vector3,_NP);
00738     
00739     memset(s2,0,sizeof(Vector3)*_NP);
00740     memset(s3,0,sizeof(Vector3)*_NP);
00741     memset(s4,0,sizeof(Vector3)*_NP);
00742     memset(s5,0,sizeof(Vector3)*_NP);
00743     memset(s2real,0,sizeof(Vector3)*_NP);
00744 
00745     memset(zetaNHC,0,sizeof(double)*MAXNHCLEN);
00746     memset(zetaNHCv,0,sizeof(double)*MAXNHCLEN);
00747     memset(zetaNHCa,0,sizeof(double)*MAXNHCLEN);
00748     memset(zetaNHC2,0,sizeof(double)*MAXNHCLEN);
00749     memset(zetaNHC3,0,sizeof(double)*MAXNHCLEN);
00750     memset(zetaNHC4,0,sizeof(double)*MAXNHCLEN);
00751     memset(zetaNHC5,0,sizeof(double)*MAXNHCLEN);
00752     
00753 //    memset(zetaBNHC,0,sizeof(double)*MAXNHCLEN);
00754 //    memset(zetaBNHCv,0,sizeof(double)*MAXNHCLEN);
00755 //    memset(zetaBNHCa,0,sizeof(double)*MAXNHCLEN);
00756 //    memset(zetaBNHC2,0,sizeof(double)*MAXNHCLEN);
00757 //    memset(zetaBNHC3,0,sizeof(double)*MAXNHCLEN);
00758 //    memset(zetaBNHC4,0,sizeof(double)*MAXNHCLEN);
00759 //    memset(zetaBNHC5,0,sizeof(double)*MAXNHCLEN);
00760 
00761     npold=_NP;
00762 }
00763 
00764 
00765 
00766 
00767 
00768 /*******************************************************
00769  *
00770  *  NVE ensemble
00771  *  Velocity Verlet integrator
00772  *
00773  ******************************************************/
00774 void MDFrame::NVE_VVerlet(int nstep)
00775 {
00776     int i;
00777     
00778     if(npold!=_NP)
00779     {
00780         Realloc(s2,Vector3,_NP);
00781         memset(s2,0,sizeof(Vector3)*_NP);
00782         npold=_NP;
00783     }
00784 
00785     /* Jan 22 2007 Keonwook Kang */
00786     if(nstep == step0) VVerlet_Get_s2();
00787     
00788     /* s2[i] is the acceleration in previous time step (n) */
00789     for(i=0;i<_NP;i++)
00790     {
00791         if(fixed[i]) continue;
00792         _VSR[i]+=s2[i];        /* p(n+1/2) */
00793         _SR[i]+=_VSR[i];       /* q(n+1) */
00794     }
00795     /* _SR[i] is the position at current step,  q(n+1) */
00796     /* _VSR[i] is the velocity at half-step,    p(n+1/2) */
00797     
00798     VVerlet_Get_s2();
00799         
00800     /* s2[i] is the acceleration at current step, 0.5*a(n+1)*(dt)^2 */
00801     
00802     for(i=0;i<_NP;i++)
00803     {
00804         if(fixed[i]) continue;
00805         _VSR[i]+=s2[i];
00806     }
00807     /* _VSR[i] is the velocity at current step, p(n+1) */
00808 }
00809 
00810 
00811 /*******************************************************
00812  *
00813  *  NVT ensemble
00814  *  Velocity Verlet integrator, implicit algorithm
00815  *
00816  ******************************************************/
00817 void MDFrame::NVT_VVerlet_Implicit(int nstep)
00818 {
00819     int i, maxiter;
00820     double tmp, xinew, B, C;
00821     Vector3 dvs;
00822     
00823     if(npold!=_NP)
00824     {
00825         Realloc(s2,Vector3,_NP);
00826         memset(s2,0,sizeof(Vector3)*_NP);
00827         npold=_NP;
00828     }
00829 
00830     /* Jan 22 2007 Keonwook Kang */
00831     if(nstep == step0)  VVerlet_Get_s2();
00832 
00833     zetaa=vt2* (_T/_TDES-1) *_TIMESTEP*_TIMESTEP*1e-24*0.5;
00834     /* zetaa = xi'(n) *(dt)^2 *0.5 */
00835 
00836     /* s2[i] is the acceleration in previous time step (n) */
00837     for(i=0;i<_NP;i++)
00838     {
00839         if(fixed[i]) continue;                
00840         dvs=_VSR[i]; dvs*=zetav*0.5; /* viscous force */
00841         _VSR[i]+=s2[i]; _VSR[i]-=dvs; /* p(n+1/2) */
00842         _SR[i]+=_VSR[i];              /* q(n+1) */
00843     }
00844     zetav+=zetaa;  /* zetav = xi(n+1/2) *(dt) */
00845     zeta +=zetav;  /* zeta  = eta(n+1) */
00846     /* _SR[i] is the position at current step,  q(n+1) */
00847     /* _VSR[i] is the velocity at half-step,    p(n+1/2) */
00848     
00849     VVerlet_Get_s2();
00850     
00851     /* s2[i] is the acceleration at current step, 0.5*a(n+1)*(dt)^2 */
00852     
00853     for(i=0;i<_NP;i++)
00854     {
00855         if(fixed[i]) continue;
00856         _VSR[i]+=s2[i];
00857     }
00858     /* _VSR[i] is the velocity at current step with an unknown factor */
00859     /*    p(n+1) * (1 + xi(n+1)*dt/2) */
00860     /* need to solve for xi(n+1) */
00861         
00862     /* first compute instananeous temperature */
00863     VVerlet_Get_Tinst();
00864         
00865     /* xi(n+1) satisfies the equation:
00866      *
00867      *   xi(n+1) = xi(n+1/2) + vt2*(_T/_TDES/(1+0.5*xi(n+1))^2 - 1)*0.5
00868      *           = xi(n+1/2) + (B/(1+0.5*xi(n+1))^2 - 1)*C
00869      */
00870     maxiter = 1000;
00871     B = _T/_TDES;
00872     C = vt2*_TIMESTEP*_TIMESTEP*1e-24*0.5;
00873     xinew = zetav;
00874     for(i=0;i<maxiter;i++)
00875     {
00876         tmp = zetav + (B/(1+0.5*xinew)/(1+0.5*xinew) - 1)*C;
00877         if(fabs(tmp-xinew)<1e-14)
00878             break;
00879         else
00880             xinew = tmp;
00881     }
00882         
00883     if(i>=maxiter)
00884     {
00885         INFO_Printf("maxiter = %d xinew = %e %e\n",maxiter,xinew,tmp);
00886         ERROR("maximum iteration exceeded");
00887         zetav = xinew;        
00888     }
00889     else
00890     {
00891         //INFO_Printf("zetav_old = %e zetav = %e iter = %d\n",zetav,xinew,i);
00892         zetav = xinew;        
00893     }
00894         
00895     tmp = 1.0/(1+zetav*0.5);
00896     for(i=0;i<_NP;i++)
00897     {
00898         if(fixed[i]) continue;
00899         _VSR[i]*=tmp;
00900     }
00901     tmp = tmp*tmp; _KATOM *= tmp; _T *= tmp;
00902     
00903 }
00904 
00905 /*******************************************************
00906  *
00907  *  NVT ensemble
00908  *  Velocity Verlet integrator, explicit algorithm 1
00909  *
00910  ******************************************************/
00911 void MDFrame::NVT_VVerlet_Explicit_1(int nstep)
00912 {
00913     int i;
00914     Vector3 dvs;
00915     
00916     if(npold!=_NP)
00917     {
00918         Realloc(s2,Vector3,_NP);
00919         memset(s2,0,sizeof(Vector3)*_NP);
00920         npold=_NP;
00921     }
00922 
00923     /* Jan 22 2007 Keonwook Kang */    
00924     if(nstep == step0) VVerlet_Get_s2();
00925 
00926     
00927     /****************************************************/
00928     /* s2[i] is the acceleration in previous time step (n) */
00929     for(i=0;i<_NP;i++)
00930     {
00931         if(fixed[i]) continue;                
00932         _VSR[i]+=s2[i]; _VSR[i]/=(1+zetav*.5); /* Keonwook Kang Feb 12 2007*/
00933         //dvs=_VSR[i]; dvs*=zetav*0.5; /* viscous force */
00934         //_VSR[i]+=s2[i]; _VSR[i]-=dvs; /* p(n+1/2) */
00935         _SR[i]+=_VSR[i];              /* q(n+1) */
00936     }
00937 
00938     /* compute instananeous temperature at half step*/
00939     VVerlet_Get_Tinst();
00940 
00941     zetaa=vt2* (_T/_TDES-1) *_TIMESTEP*_TIMESTEP*1e-24*0.5;
00942     /* zetaa = xi'(n) *(dt)^2 *0.5 */
00943 
00944     zeta +=zetav+zetaa;/* zeta  = eta(n+1) */
00945     zetav+=zetaa*2.0;  /* zetav = xi(n+1/2) *(dt) */
00946 
00947     /* Above two lines are originally given in the following
00948        three steps.                                     */
00949     /* Keonwook Kang Feb 13 2007*/
00950 //    zeta +=zetav*.5;
00951 //    zetav+=zetaa*2.0;  /* zetav = xi(n+1/2) *(dt) */
00952 //    zeta +=zetav*.5;   /* zeta  = eta(n+1) */
00953 
00954     /****************************************************/
00955     /* compute acceleration a(n+1),  s2[i]=0.5*a*(dt)^2 */
00956     VVerlet_Get_s2();
00957 
00958     /****************************************************/    
00959     for(i=0;i<_NP;i++)
00960     {
00961         if(fixed[i]) continue;
00962         dvs=_VSR[i]; dvs*=zetav*0.5; /* viscous force */
00963         _VSR[i]+=s2[i]; _VSR[i]-=dvs; /* p(n+1/2) */
00964         }
00965     /* _VSR[i] is the velocity at current step,    p(n+1) */
00966 
00967 }
00968 
00969 /*******************************************************
00970  *
00971  *  NVT ensemble
00972  *  Velocity Verlet integrator, explicit algorithm 2
00973  *  Jan 23 2007 Keonwook Kang
00974  ******************************************************/
00975 void MDFrame::NVT_VVerlet_Explicit_2(int nstep)
00976 {
00977     int i;
00978     Vector3 dvs;
00979 
00980     
00981     if(npold!=_NP)
00982     {
00983         Realloc(s2,Vector3,_NP);
00984         memset(s2,0,sizeof(Vector3)*_NP);
00985         npold=_NP;
00986     }
00987 
00988     /* Jan 22 2007 Keonwook Kang */    
00989     if(nstep == step0) VVerlet_Get_s2();
00990 
00991     /************************************************/
00992     /* compute instananeous temperature at half step*/
00993     VVerlet_Get_Tinst();
00994 
00995     zetaa=vt2* (_T/_TDES-1) *_TIMESTEP*_TIMESTEP*1e-24*0.5;
00996     /* zetaa = xi'(n) *(dt)^2 *0.5 */
00997 
00998     zetav+=zetaa;  /* zetav = xi(n+1/2) *(dt) */        
00999     zeta +=zetav;  /* zeta  = eta(n+1) */
01000     /* Originally two steps: zeta += zetav*.5; zeta += zetav*.5 */
01001     
01002     /* s2[i] is the acceleration in previous time step a(n)*(dt)^2*0.5 */
01003     for(i=0;i<_NP;i++)
01004     {
01005         if(fixed[i]) continue;
01006         dvs=_VSR[i]; dvs*=exp(-zetav*0.5); /* viscous force */
01007         _VSR[i]=dvs + s2[i]; /* p(n+1/2) */            
01008         _SR[i]+=_VSR[i];              /* q(n+1) */
01009     }
01010     
01011     VVerlet_Get_s2();
01012     
01013     for(i=0;i<_NP;i++)
01014     {
01015         if(fixed[i]) continue;
01016         _VSR[i] += s2[i];
01017         _VSR[i] *= exp(-zetav*0.5); /* viscous force */
01018     }
01019 
01020 //        zeta +=zetav*0.5;/* zeta  = eta(n+1) */                
01021     /* compute instananeous temperature at half step*/
01022     VVerlet_Get_Tinst();
01023 
01024     zetaa=vt2* (_T/_TDES-1) *_TIMESTEP*_TIMESTEP*1e-24*0.5;
01025     /* zetaa = xi'(n) *(dt)^2 *0.5 */
01026 
01027     zetav+=zetaa;  /* zetav = xi(n+1/2) *(dt) */
01028 }
01029 
01030 /*******************************************************
01031  *
01032  *  NVT ensemble Nose-Hoover Chain
01033  *  Velocity Verlet integrator, 
01034  *  Mar 20 2007 implemented by Keonwook Kang
01035  *  Ref. G. J. Martyna, Molecular Physics (1996) 87 1117-1157
01036  *
01037  ******************************************************/
01038 void MDFrame::NVTC_VVerlet_Explicit(int nstep)
01039 {
01040     int i;
01041     Vector3 dvs;
01042     
01043     if(npold!=_NP) /* initialization */
01044     {
01045         Realloc(s2,Vector3,_NP);
01046         memset(s2,0,sizeof(Vector3)*_NP);
01047         memset(zetaNHC,0,sizeof(double)*MAXNHCLEN);
01048         memset(zetaNHCv,0,sizeof(double)*MAXNHCLEN);
01049         memset(zetaNHCa,0,sizeof(double)*MAXNHCLEN);
01050         npold=_NP;
01051     }
01052 
01053     VVerlet_NHCINT(); /* Nose-Hoover Chain integrator */
01054 
01055     /* Jan 22 2007 Keonwook Kang */    
01056     if(nstep == step0) VVerlet_Get_s2();
01057     /* s2[i] is the acceleration in previous time step a(n)*(dt)^2*0.5 */
01058     for(i=0;i<_NP;i++)
01059     {
01060         if(fixed[i]) continue;
01061         _VSR[i] += s2[i]; /* p(n+1/2) */            
01062         _SR[i] += _VSR[i];              /* q(n+1) */
01063     }
01064     
01065     VVerlet_Get_s2();
01066     
01067     for(i=0;i<_NP;i++)
01068     {
01069         if(fixed[i]) continue;
01070         _VSR[i] += s2[i];
01071     }
01072 
01073     VVerlet_NHCINT(); /* Nose-Hoover Chain integrator */
01074 }
01075 
01076 /*******************************************************
01077  *
01078  *  NPH ensemble
01079  *  Velocity Verlet integrator
01080  *
01081  ******************************************************/
01082 void MDFrame::NPH_VVerlet(int nstep)
01083 {
01084     int i, j;
01085     Matrix33 hinv;
01086     
01087     if(npold!=_NP)
01088     {
01089         Realloc(s2,Vector3,_NP);
01090         memset(s2,0,sizeof(Vector3)*_NP);
01091         npold=_NP;
01092     }
01093 
01094     /* Jan 22 2007 Keonwook Kang */    
01095     if(nstep == step0)
01096     {
01097         VVerlet_Get_s2();
01098         VVerlet_Get_h2();
01099     }
01100 
01101 //    for(i=0;i<_NP;i++)
01102 //        s2[i]-=GiVG*_VSR[i];
01103     
01104     /* s2[i] is the acceleration in previous time step (n) */
01105     for(i=0;i<_NP;i++)
01106     {
01107         if(fixed[i]) continue;
01108         _VSR[i]+=s2[i];        /* p(n+1/2) */
01109         _SR[i]+=_VSR[i];       /* q(n+1) */
01110     }
01111 
01112     for(i=0;i<3;i++)
01113         for(j=0;j<3;j++)
01114             if(conj_fixboxvec[i][j]==1)
01115                 _VH[i][j]=h2[i][j]=0;
01116     _VH+=h2;
01117     _H+=_VH;
01118     /* _SR[i] is the position at current step,  q(n+1) */
01119     /* _VSR[i] is the velocity at half-step,    p(n+1/2) */
01120 
01121     VVerlet_Get_s2();
01122     
01123     /* s2[i] is the acceleration at current step, 0.5*a(n+1)*(dt)^2 */
01124     for(i=0;i<_NP;i++)
01125     {
01126         if(fixed[i]) continue;
01127         _VSR[i]+=s2[i];
01128     }
01129     for(i=0;i<3;i++)
01130         for(j=0;j<3;j++)
01131             if(conj_fixboxvec[i][j]==1)
01132                 _VH[i][j]=h2[i][j]=0;
01133     _VH+=h2;
01134     /* _VSR[i] is the velocity at current step, p(n+1) */
01135 }
01136 
01137 
01138 /*******************************************************
01139  *
01140  *  NPH ensemble
01141  *  Velocity Verlet integrator
01142  *  Feb 15 2007 Keonwook Kang
01143  *  Finished Feb 19 2007
01144  ******************************************************/
01145 void MDFrame::NPH_VVerlet_Explicit_1(int nstep)
01146 {
01147     int i, j, n;
01148     Matrix33 htran, G, Ginv, vhtran, VG1, VG2, VG;
01149     Matrix33 GiVG, GiVGoffDiag, I, Iinv;
01150     Vector3 dvs;
01151 
01152     I.eye();
01153     
01154     if(npold!=_NP)
01155     {
01156         Realloc(s2,Vector3,_NP);
01157         memset(s2,0,sizeof(Vector3)*_NP);
01158         npold=_NP;
01159     }
01160 
01161     if(nstep == step0) VVerlet_Get_h2();
01162 
01163     for(i=0;i<3;i++)
01164         for(j=0;j<3;j++)
01165             if(conj_fixboxvec[i][j]==1)
01166                 _VH[i][j]=h2[i][j]=0;
01167     _VH+=h2;
01168     _H+=_VH*.5;
01169 
01170     VVerlet_Get_s2();   /* eval s2 due to box change */
01171     
01172     htran = _H.tran(); G=htran*_H; Ginv = G.inv();
01173     vhtran=_VH.tran(); VG1=vhtran*_H; VG2=htran*_VH;
01174     VG=VG1+VG2; GiVG = Ginv*VG;
01175     GiVGoffDiag.set(0, GiVG[0][1], GiVG[0][2],
01176                     GiVG[1][0], 0, GiVG[1][2],
01177                     GiVG[2][0], GiVG[2][1], 0);
01178 
01179     /* s2[i] is the acceleration in previous time step (n) */
01180     for(i=0;i<_NP;i++)
01181     {
01182         if(fixed[i]) continue;
01183         s2[i] -= GiVGoffDiag*_VSR[i]*.5;
01184         _VSR[i] += s2[i];
01185         dvs = _VSR[i]; 
01186         for(n=0;n<3;n++)
01187             dvs[n] = dvs[n]*exp((-.5)*GiVG[n][n]);        
01188         _VSR[i] = dvs;
01189         _SR[i] += _VSR[i];        
01190     }
01191 
01192     /* _SR[i] is the position at full step,  q(n+1) */
01193     /* _VSR[i] is the velocity at half-step,    p(n+1/2) */
01194 
01195     VVerlet_Get_s2(); /* eval s2 due to position change */
01196 
01197 //    INFO("H in NPH=["<<_H<<"]");
01198 
01199     I += GiVGoffDiag*.5; Iinv = I.inv();
01200     for(i=0;i<_NP;i++)
01201     {
01202         if(fixed[i]) continue;
01203         dvs = _VSR[i];
01204         for(n=0;n<3;n++)
01205             dvs[n] = dvs[n]*exp((-.5)*GiVG[n][n]);
01206         _VSR[i] = dvs + s2[i];
01207         _VSR[i] = Iinv*_VSR[i];
01208     }
01209     /* _VSR[i] is the velocity at current step, p(n+1) */
01210     
01211     _H+=_VH*.5;
01212     VVerlet_Get_s2(); /* eval s2 due to box change */     
01213     VVerlet_Get_h2(); /* Need to eval h2 due to Virial stress change */
01214     for(i=0;i<3;i++)
01215         for(j=0;j<3;j++)
01216             if(conj_fixboxvec[i][j]==1)
01217                 _VH[i][j]=h2[i][j]=0;
01218     _VH+=h2;
01219 
01220 }
01221 
01222 /*******************************************************
01223  *
01224  *  NPT ensemble
01225  *  Velocity Verlet integrator
01226  *
01227  ******************************************************/
01228 void MDFrame::NPT_VVerlet_Implicit(int nstep)
01229 {
01230     FATAL("NPT_VVerlet not implemented yet.  Come back soon!");
01231 }
01232 
01233 void MDFrame::NPT_VVerlet_Explicit_1(int nstep)
01234 {
01235     FATAL("NPT_VVerlet not implemented yet.  Come back soon!");
01236 }
01237 
01238 void MDFrame::NPT_VVerlet_Explicit_2(int nstep)
01239 {
01240     FATAL("NPT_VVerlet not implemented yet.  Come back soon!");
01241 }
01242 
01243 
01244 /* VVerlet Get s2[i] for the 1st run */
01245 /* Jan 22 Keonwook Kang              */
01246 void MDFrame::VVerlet_Get_s2()
01247 {
01248     int i;
01249     double tmp;
01250     Matrix33 hinv;
01251     
01252     /* compute acceleration a(n+1),  s2[i]=0.5*a*(dt)^2 */
01253     call_potential();
01254     calcprop();
01255 
01256     hinv=_H.inv();
01257 
01258     /* fix me! need to change the unit of velocity to A/ps */
01259     tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24  //(ps^2)
01260         /(_ATOMMASS[0]*1e-3/AVO)*EV/1e-20; //(A^2)
01261 
01262     for(i=0;i<_NP;i++)
01263     {
01264         if(fixed[i]) continue;
01265         if(species[i]==0)
01266             s2[i]=hinv*_F[i]*tmp;
01267         else
01268             s2[i]=hinv*_F[i]*tmp *_ATOMMASS[0]/_ATOMMASS[species[i]];
01269     }    
01270 }
01271 
01272 void MDFrame::VVerlet_Get_h2()
01273 {
01274     int i, j;    
01275     double mass, tmp;
01276 
01277 //    call_potential();
01278     calcprop();     /* To get GH */
01279     
01280     /* fix me! need to change the unit of velocity to A/ps */
01281     /* Feb 15 2007 Keonwook Kang
01282        Notice no "1e-20" in the denominatro!!! */
01283 //    tmp=0.5*(_TIMESTEP)*(_TIMESTEP)*1e-24  //(ps^2)
01284 //        /(_ATOMMASS[0]*1e-3/AVO)*EV; //
01285     
01286     /* ATOMMASS has unit of g/mol
01287      * mass     has unit of eV ps^2 / A^2 
01288      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
01289 //    mass = (_ATOMMASS[0]*1e-3/AVO/EV/1e20)*1e24;
01290     mass = _ATOMMASS[0]*MASSCONVERT;
01291     tmp  = 0.5 / mass * _TIMESTEP*_TIMESTEP;
01292     
01293 //  h2=_GH*(-tmp)* _ATOMMASS[0]/_WALLMASS;
01294     h2=_GH*tmp*_ATOMMASS[0]/_WALLMASS;  /* unit of angstrom */
01295 // change its sign after commenting out _GH*=-1; in calcprop() of md.cpp
01296 
01297     if(_BOXDAMP!=0)
01298         for(i=0;i<3;i++)
01299             for(j=0;j<3;j++)
01300                 if(conj_fixboxvec[i][j]==0)
01301                     h2[i][j]-=_VH[i][j]*_BOXDAMP;
01302 
01303 }
01304 
01305 
01306 void MDFrame::VVerlet_Get_Tinst()
01307 {
01308     int i;
01309     double mass, tmp;
01310     /* ATOMMASS has unit of g/mol
01311      * mass     has unit of eV ps^2 / A^2 
01312      * tmp = 0.5 * dt^2 / m  has unit of (A^2/eV) */
01313     mass = _ATOMMASS[0]*MASSCONVERT;
01314     tmp  = 0.5 * mass / (_TIMESTEP*_TIMESTEP);    
01315     
01316     _KATOM=0; _NPfree=0;
01317     for(i=0;i<_NP;i++)
01318     {
01319         if(fixed[i]) continue;
01320         _NPfree++;
01321         _VR[i]=_H*_VSR[i];
01322         if(species[i]==0)
01323             _KATOM+=_VR[i].norm2(); // in the unit of Angstrom^2
01324         else
01325             _KATOM+=_VR[i].norm2()*_ATOMMASS[species[i]]/_ATOMMASS[0];
01326     }
01327     _KATOM *= tmp; // 1/2 times (unit conversion from Angstrom^2 to eV )
01328 
01329     _T=_KATOM/(1.5*KB*_NPfree);    
01330 }
01331     
01332 void MDFrame::VVerlet_NHCINT()
01333 { /* Nose-Hoover Chain integrator from t=0 to t=_TIMESTEP/2 */
01334     int i;
01335     double expterm;
01336 
01337     VVerlet_Get_Tinst();
01338     
01339     i=NHChainLen-1;
01340 
01341     if(NHChainLen==1)
01342         zetaNHCa[i]=(_KATOM*2 - 3*_NPfree*KB*_TDES)/NHMass[i]
01343                     *(0.5*SQR(_TIMESTEP));
01344     else
01345         zetaNHCa[i]=(NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES)/NHMass[i]
01346                 *(0.5*SQR(_TIMESTEP));
01347     zetaNHCv[i]+=.5*zetaNHCa[i];
01348     
01349     for(i=NHChainLen-2;i>=0;i--)
01350     {
01351         expterm = exp(zetaNHCv[i+1]/(-8.0));
01352         if (i!=0)
01353         {
01354             zetaNHCa[i]=(NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES)/NHMass[i]
01355                         *(0.5*SQR(_TIMESTEP));
01356         }
01357         else
01358         {
01359             zetaNHCa[i]=(_KATOM*2 - 3*_NPfree*KB*_TDES)/NHMass[i]
01360                         *(0.5*SQR(_TIMESTEP));
01361         }
01362         zetaNHCv[i] *= SQR(expterm);
01363         zetaNHCv[i] += 0.5*zetaNHCa[i]*expterm;
01364     }
01365     
01366     expterm = exp(-0.5*zetaNHCv[0]);
01367     for(i=0;i<_NP;i++)  _VSR[i]*=expterm;
01368     for(i=0;i<=NHChainLen-1;i++)
01369         zetaNHC[i]+=0.5*zetaNHCv[i];
01370 
01371     VVerlet_Get_Tinst();
01372 
01373     for(i=0;i<NHChainLen-1;i++)
01374     {
01375         expterm = exp(zetaNHCv[i+1]/(-8.0));
01376         if (i!=0)
01377         {
01378             zetaNHCa[i]=(NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES)/NHMass[i]
01379                         *(0.5*SQR(_TIMESTEP));
01380         }
01381         else
01382         {
01383             zetaNHCa[i]=(_KATOM*2 - 3*_NPfree*KB*_TDES)/NHMass[i]
01384                         *(0.5*SQR(_TIMESTEP));
01385         }
01386         zetaNHCv[i] *= SQR(expterm);
01387         zetaNHCv[i] += 0.5*zetaNHCa[i]*expterm;
01388     }
01389     i=NHChainLen-1;
01390     if(NHChainLen==1)
01391         zetaNHCa[i]=(_KATOM*2 - 3*_NPfree*KB*_TDES)/NHMass[i]
01392                     *(0.5*SQR(_TIMESTEP));
01393     else
01394         zetaNHCa[i]=(NHMass[i-1]*SQR(zetaNHCv[i-1]/_TIMESTEP)-KB*_TDES)/NHMass[i]
01395                 *(0.5*SQR(_TIMESTEP));
01396     zetaNHCv[i]+=.5*zetaNHCa[i];    
01397 }