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

00001 /*
00002   siedip.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Thu Jul 19 18:55:46 2007
00005 
00006   FUNCTION  :  MD simulation package of Si using EDIP potential
00007 */
00008 
00009 #include "siedip.h"
00010 
00011 void SIEDIPFrame::initvars()
00012 {
00013     A = 5.6714030;
00014     B = 2.0002804;
00015     rh = 1.2085196;
00016     a = 3.1213820;
00017     sig = 0.5774108;
00018     lam = 1.4533108;
00019     gam = 1.1247945;
00020     b = 3.1213820;
00021     c = 2.5609104;
00022     delta = 78.7590539;
00023     mu = 0.6966326;
00024     Qo = 312.1341346;
00025     palp = 1.4074424;
00026     bet = 0.0070975;
00027     alp = 3.1083847;
00028     
00029     u1 = -0.165799;
00030     u2 = 32.557;
00031     u3 = 0.286198;
00032     u4 = 0.66; 
00033     
00034     INFO_Printf("\n EDIP-Si Parameters: \n");
00035     INFO_Printf("%e %e %e %e %e \n",A,B,rh,a,sig);
00036     INFO_Printf("%e %e %e %e %e \n",lam,gam,b,c,delta);
00037     INFO_Printf("%e %e %e %e %e \n",mu,Qo,palp,bet,alp);
00038     bg=a;
00039     eta = delta/Qo;
00040     
00041     _RLIST=a*1.1;
00042     _SKIN=_RLIST-a;
00043 
00044     strcpy(incnfile,"../siedip.cn");
00045     DUMP(HIG"SIEDIPFrame initvars"NOR);
00046     MDFrame::initvars();
00047 }
00048 
00049 void SIEDIPFrame::edip()
00050 {
00051     /*------------------------- VARIABLE DECLARATIONS -------------------------*/
00052     
00053     int i,j;
00054     double dx,dy,dz,r,rsqr,asqr;
00055     double rinv,rmainv,xinv,xinv3,den,Z,fZ;
00056     double dV2j,dV2ijx,dV2ijy,dV2ijz,pZ,dp;
00057     double temp0,temp1;
00058     double Qort,muhalf,u5;
00059     double rmbinv,winv,dwinv,tau,dtau,lcos,x,H,dHdx,dhdl;
00060     double dV3rij,dV3rijx,dV3rijy,dV3rijz;
00061     double dV3rik,dV3rikx,dV3riky,dV3rikz;
00062     double dV3l,dV3ljx,dV3ljy,dV3ljz,dV3lkx,dV3lky,dV3lkz;
00063     double dV2dZ,dxdZ,dV3dZ;
00064     double dEdrl,dEdrlx,dEdrly,dEdrlz;
00065     double bmc,cmbinv;
00066     double fjx,fjy,fjz,fkx,fky,fkz;
00067     
00068 #ifndef NNM /* need to be taken care of later */
00069 #define NNM 60
00070 #endif
00071     typedef struct{      
00072         double t0,t1,t2,t3;  /* various V2 functions and derivatives */
00073         double dx,dy,dz;     /* unit separation vector */
00074         double r;                /* bond length (only needed for virial) */
00075     } store2;
00076     store2 s2[NNM]; /* two-body interaction quantities, r<a */  
00077     int n2;                /* size of s2[] */
00078     int num2[NNM];  /* atom ID numbers for s2[] */
00079     
00080     typedef struct{
00081         double g,dg;         /* 3-body radial function and its derivative */
00082         double rinv;         /* 1/r */
00083         double dx,dy,dz;     /* unit separation vector */
00084         double r;                /* bond length (only needed for virial) */
00085     } store3;
00086     store3 s3[NNM]; /* three-body interaction quantities, r<b */
00087     int n3;                /* size of s3[] */
00088     int num3[NNM];  /* atom ID numbers for s3[] */
00089     
00090     typedef struct{
00091         double df;           /* derivative of neighbor function f'(r) */
00092         double sum;              /* array to accumulate coordination force prefactors */
00093         double dx,dy,dz;     /* unit separation vector */
00094         double r;                /* bond length (only needed for virial) */
00095     } storez;
00096     storez sz[NNM]; /* coordination number stuff, c<r<b */
00097     int nz;                /* size of sz[] */
00098     int numz[NNM];  /* atom ID numbers for sz[] */
00099     
00100     int nj,nk,nl;         /* indices for the store arrays */
00101 
00102     int n0, n1;
00103     int ipt, jpt, kpt, lpt;
00104     double V2, V3, virial; /* virial is the diagonal part of _VIRIAL */
00105     double eterm2, eterm3;
00106     Vector3 sij, rij, fgj, rgj, fgk, rgk;
00107     
00108     DUMP("EDIP");
00109     
00110     refreshneighborlist();
00111 
00112 //    n0=(int)rint(MP_cpu()*_NP/MP_ncpu());
00113 //    n1=(int)rint((MP_cpu()+1)*_NP/MP_ncpu());
00114     n0=0;
00115     n1=_NP;
00116 
00117     /* _EPOT, _F, _EPOT_IND, _VIRIAL all local */
00118     for(i=n0;i<n1;i++)
00119     { _F[i].clear(); _EPOT_IND[i]=0;}
00120     _VIRIAL.clear();
00121 
00122     _EPOT=V2=V3=0; /* _EPOT = V2 + V3 */
00123 
00124     /* COMBINE COEFFICIENTS */
00125     
00126     asqr = a*a;
00127     Qort = sqrt(Qo);
00128     muhalf = mu*0.5;
00129     u5 = u2*u4;    
00130     bmc = b-c;
00131     cmbinv = 1.0/(c-b);
00132 
00133     /* strange, why needed? */
00134     dx=dy=dz=0;
00135     
00136     /*--- LEVEL 1: OUTER LOOP OVER ATOMS ---*/
00137     for(ipt=n0; ipt<n1; ipt++)
00138     {   /* not parallelized yet */
00139         /* RESET COORDINATION AND NEIGHBOR NUMBERS */
00140 
00141         Z = 0.0;
00142         n2 = 0;  
00143         n3 = 0;
00144         nz = 0;
00145 
00146         /*--- LEVEL 2: LOOP PREPASS OVER PAIRS ---*/
00147 
00148         for(j=0; j<nn[ipt]; j++)
00149         {
00150             jpt=nindex[ipt][j];
00151             sij=_SR[jpt]-_SR[ipt];
00152             sij.subint();
00153             rij=_H*sij;
00154             dx=rij.x; dy=rij.y; dz=rij.z;
00155 
00156             /* TEST IF WITHIN OUTER CUTOFF */
00157             rsqr=rij.norm2();
00158             if(rsqr>=asqr) continue;
00159             r=sqrt(rsqr);
00160             
00161             /* PARTS OF TWO-BODY INTERACTION r<a */
00162  
00163             num2[n2] = jpt;
00164             rinv = 1.0/r;
00165             dx *= rinv;
00166             dy *= rinv;
00167             dz *= rinv;
00168             rmainv = 1.0/(r-a);
00169             s2[n2].t0 = A*exp(sig*rmainv);
00170             s2[n2].t1 = pow(B*rinv,rh);
00171             s2[n2].t2 = rh*rinv;
00172             s2[n2].t3 = sig*rmainv*rmainv;
00173             s2[n2].dx = dx;
00174             s2[n2].dy = dy;
00175             s2[n2].dz = dz;
00176             s2[n2].r = r;
00177             n2++;
00178 
00179             /* RADIAL PARTS OF THREE-BODY INTERACTION r<b */
00180             
00181             if(r < bg) {
00182                 
00183                 num3[n3] = jpt;
00184                 rmbinv = 1.0/(r-bg);
00185                 temp1 = gam*rmbinv;
00186                 temp0 = exp(temp1);
00187 #if V3g_on
00188                 s3[n3].g = temp0;
00189                 s3[n3].dg = -rmbinv*temp1*temp0;
00190 #else
00191                 s3[n3].g = 1;
00192                 s3[n3].dg = 0;
00193 #endif
00194                 s3[n3].dx = dx;
00195                 s3[n3].dy = dy;
00196                 s3[n3].dz = dz;
00197                 s3[n3].rinv = rinv;
00198                 s3[n3].r = r;
00199                 n3++;
00200                 
00201                 /* COORDINATION AND NEIGHBOR FUNCTION c<r<b */
00202 
00203                 if(r<b) {
00204                     if(r < c) 
00205                         Z += 1.0;
00206                     else {
00207                         xinv = bmc/(r-c);
00208                         xinv3 = xinv*xinv*xinv;
00209                         den = 1.0/(1 - xinv3);
00210                         temp1 = alp*den;
00211                         fZ = exp(temp1);
00212                         Z += fZ;
00213                         numz[nz] = jpt;
00214                         sz[nz].df = fZ*temp1*den*3.0*xinv3*xinv*cmbinv;   /* df/dr */
00215                         sz[nz].dx = dx;
00216                         sz[nz].dy = dy;
00217                         sz[nz].dz = dz;
00218                         sz[nz].r = r;
00219                         nz++;
00220                     }      
00221                 }
00222             }
00223         }
00224         //local[ipt].z=Z;
00225 
00226         //if(measure) coord_total += Z;
00227 
00228         /* ZERO ACCUMULATION ARRAY FOR ENVIRONMENT FORCES */
00229 
00230         for(nl=0; nl<nz; nl++) sz[nl].sum=0.0;
00231 
00232         /* ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION */
00233 
00234 #if V2_on
00235         
00236 #if V2Z_on
00237         temp0 = bet*Z;
00238         pZ = palp*exp(-temp0*Z);         /* bond order */
00239         dp = -2.0*temp0*pZ;         /* derivative of bond order */
00240 #else
00241         pZ = palp*exp(-bet*16);
00242         dp=0.0;
00243 #endif
00244 
00245         /*--- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---*/
00246 
00247         for(nj=0; nj<n2; nj++)
00248         {
00249 
00250             temp0 = s2[nj].t1 - pZ;
00251 
00252             /* two-body energy V2(rij,Z) */
00253             eterm2 = temp0*s2[nj].t0; 
00254             V2 += eterm2;
00255             jpt = num2[nj];
00256 
00257             _EPOT_IND[ipt]+=eterm2*0.5;
00258             _EPOT_IND[jpt]+=eterm2*0.5;
00259             
00260             /* two-body forces */
00261 
00262             dV2j = - (s2[nj].t0) * ((s2[nj].t1)*(s2[nj].t2) 
00263                                     + temp0 * (s2[nj].t3));   /* dV2/dr */
00264             dV2ijx = dV2j * s2[nj].dx;
00265             dV2ijy = dV2j * s2[nj].dy;
00266             dV2ijz = dV2j * s2[nj].dz;
00267             fgj.set(dV2ijx, dV2ijy, dV2ijz);
00268             _F[ipt] += fgj;
00269             _F[jpt] -= fgj;
00270             
00271             /* dV2/dr contribution to virial */
00272             virial -= s2[nj].r * (dV2ijx*s2[nj].dx 
00273                               + dV2ijy*s2[nj].dy + dV2ijz*s2[nj].dz);
00274 
00275             rgj.set(s2[nj].dx, s2[nj].dy, s2[nj].dz);
00276             _VIRIAL.addnvv(-s2[nj].r,fgj,rgj);
00277             
00278             /*--- LEVEL 3: LOOP FOR PAIR COORDINATION FORCES ---*/
00279             
00280             dV2dZ = - dp * s2[nj].t0;
00281             for(nl=0; nl<nz; nl++) sz[nl].sum += dV2dZ;
00282             
00283         }
00284 #endif
00285 
00286         /* COORDINATION-DEPENDENCE OF THREE-BODY INTERACTION */
00287         
00288 #if V3_on
00289         
00290 #if V3Z_on
00291         winv = Qort*exp(-muhalf*Z); /* inverse width of angular function */    
00292         dwinv = -muhalf*winv;       /* its derivative */
00293         temp0 = exp(-u4*Z);
00294         tau = u1+u2*temp0*(u3-temp0); /* -cosine of angular minimum */
00295         dtau = u5*temp0*(2*temp0-u3); /* its derivative */
00296 #else
00297         winv =  Qort*exp(-muhalf*4);
00298         dwinv = 0.0;
00299         tau = 1.0/3.0;
00300         dtau=0.0;
00301 #endif
00302 
00303         /*--- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---*/
00304         
00305         for(nj=0; nj<(n3-1); nj++)
00306         {
00307             jpt=num3[nj];
00308             
00309             /*--- LEVEL 3: SECOND LOOP FOR THREE-BODY INTERACTIONS ---*/
00310             for(nk=nj+1; nk<n3; nk++)
00311             {
00312                 kpt=num3[nk];
00313 
00314                 /* angular function h(l,Z) */
00315                 
00316                 lcos = s3[nj].dx * s3[nk].dx + s3[nj].dy * s3[nk].dy
00317                     + s3[nj].dz * s3[nk].dz;
00318                 x = (lcos + tau)*winv;
00319                 temp0 = exp(-x*x);
00320 #if V3h_on
00321                 H = lam*(1 - temp0 + eta*x*x);
00322                 dHdx = 2*lam*x*(temp0 + eta);
00323                 dhdl = dHdx*winv;
00324 #else
00325                 H=1.0;
00326                 dhdl=0.0;
00327 #endif
00328 
00329                 /* three-body energy */
00330                 
00331                 temp1 = s3[nj].g * s3[nk].g;
00332                 eterm3 = temp1*H;
00333                 
00334                 V3 += eterm3;
00335                 _EPOT_IND[ipt] += eterm3/3;
00336                 _EPOT_IND[jpt] += eterm3/3;
00337                 _EPOT_IND[kpt] += eterm3/3;
00338                 
00339                 /* (-) radial force on atom j */
00340                 dV3rij = s3[nj].dg * s3[nk].g * H;
00341                 dV3rijx = dV3rij * s3[nj].dx;
00342                 dV3rijy = dV3rij * s3[nj].dy;
00343                 dV3rijz = dV3rij * s3[nj].dz;
00344                 fjx = dV3rijx;
00345                 fjy = dV3rijy;
00346                 fjz = dV3rijz;
00347                 
00348                 /* (-) radial force on atom k */
00349                 
00350                 dV3rik = s3[nj].g * s3[nk].dg * H;
00351                 dV3rikx = dV3rik * s3[nk].dx;
00352                 dV3riky = dV3rik * s3[nk].dy;
00353                 dV3rikz = dV3rik * s3[nk].dz;
00354                 fkx = dV3rikx;
00355                 fky = dV3riky;
00356                 fkz = dV3rikz;
00357                 
00358                 /* (-) angular force on j */
00359                 
00360                 dV3l = temp1*dhdl;
00361                 dV3ljx = dV3l * (s3[nk].dx - lcos * s3[nj].dx) * s3[nj].rinv;
00362                 dV3ljy = dV3l * (s3[nk].dy - lcos * s3[nj].dy) * s3[nj].rinv;
00363                 dV3ljz = dV3l * (s3[nk].dz - lcos * s3[nj].dz) * s3[nj].rinv;
00364                 fjx += dV3ljx;
00365                 fjy += dV3ljy;
00366                 fjz += dV3ljz;
00367                 
00368                 /* (-) angular force on k */
00369                 
00370                 dV3lkx = dV3l * (s3[nj].dx - lcos * s3[nk].dx) * s3[nk].rinv;
00371                 dV3lky = dV3l * (s3[nj].dy - lcos * s3[nk].dy) * s3[nk].rinv;
00372                 dV3lkz = dV3l * (s3[nj].dz - lcos * s3[nk].dz) * s3[nk].rinv;
00373                 fkx += dV3lkx;
00374                 fky += dV3lky;
00375                 fkz += dV3lkz;
00376                 
00377                 /* apply radial + angular forces to i, j, k */
00378 
00379                 fgj.set(fjx, fjy, fjz);
00380                 fgk.set(fkx, fky, fkz);
00381 
00382                 _F[jpt] -= fgj;
00383                 _F[kpt] -= fgk;
00384                 _F[ipt] += fgj; _F[ipt] += fgk;
00385 
00386                 /* dV3/dR contributions to virial */
00387                 
00388                 virial -= s3[nj].r * (fjx*s3[nj].dx 
00389                                       + fjy*s3[nj].dy + fjz*s3[nj].dz); 
00390                 virial -= s3[nk].r * (fkx*s3[nk].dx 
00391                                       + fky*s3[nk].dy + fkz*s3[nk].dz);
00392                 
00393                 rgj.set(s3[nj].dx,s3[nj].dy,s3[nj].dz); 
00394                 rgk.set(s3[nk].dx,s3[nk].dy,s3[nk].dz);
00395                 _VIRIAL.addnvv(-s3[nj].r,fgj,rgj);
00396                 _VIRIAL.addnvv(-s3[nk].r,fgk,rgk);
00397                 
00398                 /* prefactor for 4-body forces from coordination */
00399 #if V3Z_on
00400                 dxdZ = dwinv*(lcos + tau) + winv*dtau;
00401                 dV3dZ = temp1*dHdx*dxdZ;
00402                 
00403                 /*--- LEVEL 4: LOOP FOR THREE-BODY COORDINATION FORCES ---*/
00404                 
00405                 for(nl=0; nl<nz; nl++) sz[nl].sum += dV3dZ; 
00406 #endif
00407             }
00408         }
00409 #endif
00410 
00411         /*--- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---*/
00412         
00413         for(nl=0; nl<nz; nl++)
00414         {
00415             dEdrl = sz[nl].sum * sz[nl].df;
00416             dEdrlx = dEdrl * sz[nl].dx;
00417             dEdrly = dEdrl * sz[nl].dy;
00418             dEdrlz = dEdrl * sz[nl].dz;
00419             _F[ipt].x += dEdrlx;
00420             _F[ipt].y += dEdrly;
00421             _F[ipt].z += dEdrlz;
00422             lpt = numz[nl];
00423             _F[lpt].x -= dEdrlx;
00424             _F[lpt].y -= dEdrly;
00425             _F[lpt].z -= dEdrlz;
00426 
00427             /* dE/dZ*dZ/dr contribution to virial */
00428             
00429             virial -= sz[nl].r * (dEdrlx*sz[nl].dx 
00430                     + dEdrly*sz[nl].dy + dEdrlz*sz[nl].dz); 
00431         }
00432     }
00433 
00434     _EPOT = V2 + V3;
00435     virial /= 3.0;
00436 
00437     /* zero out the forces on fixed atoms */
00438     for(i=0;i<_NP;i++)
00439     {
00440         if(fixed[i])
00441             _F[i].clear();
00442     }
00443 
00444 #if 0
00445     DUMP("EDIP Collect");
00446     /* _EPOT, _F, _EPOT_IND, _VIRIAL will be summed */
00447     MP_AccumulateTo(_shmEPOTptr,&(_EPOT),1);
00448     MP_AccumulateTo((double *)_shmF,(double *)_F,3*_NP);
00449     MP_AccumulateTo(_shmEPOT_IND,_EPOT_IND,_NP);
00450     MP_AccumulateTo((double *)((*_shmVIRIALptr)[0]),
00451                         (double *)(_VIRIAL[0]),9);
00452     
00453     StrMemCpy(&(_EPOT),_shmEPOTptr,sizeof(double));
00454     StrMemCpy((double *)_F,(double *)_shmF,3*_NP*sizeof(double));
00455     StrMemCpy(_EPOT_IND,_shmEPOT_IND,_NP*sizeof(double));
00456     StrMemCpy((double *)(_VIRIAL[0]),
00457               (double *)((*_shmVIRIALptr)[0]),9*sizeof(double));
00458     
00459     DUMP("EDIP complete");
00460 #endif
00461 }
00462     
00463 void SIEDIPFrame::potential()
00464 {
00465     edip();
00466 }
00467 
00468 
00469 
00470 
00471 #ifdef _TEST
00472 
00473 /* Main Program Begins */
00474 class SIEDIPFrame sim;
00475 
00476 /* The main program is defined here */
00477 #include "main.cpp"
00478 
00479 #endif//_TEST
00480