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

00001 /*
00002   swlj.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Mon Jan  1 21:30:24 2007
00005 
00006   FUNCTION  :  MD simulation package of Si using Stillinger-Weber potential
00007                and Lennard-Jones potential (depending on group ID)
00008                group 0 and group 0 interact with SW potential
00009                group 1 and group 1 interact with SW potential
00010                group 0 and group 1 interact with LJ potential
00011 
00012                This model is developed for Nano-indentation simulation.
00013 */
00014 
00015 #include "swlj.h"
00016 
00017 void SWFrame::initparser()
00018 {
00019     MDPARALLELFrame::initparser();
00020 
00021     /* input */
00022     bindvar("C12",&_ALJ,DOUBLE);
00023     bindvar("C6",&_BLJ,DOUBLE);
00024 
00025     /* output */
00026     bindvar("tipforce",&(tipforce[0]),DOUBLE);
00027 }
00028 
00029 int SWFrame::exec(char *name)
00030 {
00031     if(MDPARALLELFrame::exec(name)==0) return 0;
00032     bindcommand(name,"initLJ",initLJ());
00033     return -1;
00034 }
00035 
00036 void SWFrame::initvars()
00037 {
00038     /* modified Si parameters PRB 46, 2250 (1992) */
00039     aa=16.31972277; bb=11.60319228; plam=48.61499998;
00040     pgam=2.51412;  acut=3.77118; pss=2.0951; rho=4.;
00041 
00042     
00043     _RLIST=acut*1.1;
00044     _SKIN=_RLIST-acut;
00045     rho1=rho+1; acutsq=acut*acut;
00046 
00047     strcpy(incnfile,"../si.cn");
00048     DUMP(HIG"SWFrame initvars"NOR);
00049     MDPARALLELFrame::initvars();
00050 }
00051 
00052 void SWFrame::initLJ()
00053 {
00054     ALJ=_ALJ*LJ_ENERGY*POW12(LJ_LENGTH);
00055     BLJ=_BLJ*LJ_ENERGY*POW6(LJ_LENGTH);
00056     Uc=(ALJ*POW6(1/LJ_RC)-BLJ)*POW6(1/LJ_RC);
00057     DUDRc=-(12*ALJ/POW6(LJ_RC)-6*BLJ)/POW6(LJ_RC)/LJ_RC;
00058 }
00059 
00060 void SWFrame::calcprop()
00061 {
00062     MDPARALLELFrame::calcprop();
00063 }
00064 
00065 void SWFrame::stillinger_weber()
00066 {
00067 #ifndef NNM /* need to be taken care of later */
00068 #define NNM 60
00069 #endif
00070     /* Multi process function */
00071     int i,j,k,ipt,jpt,kpt;
00072     int neg;
00073     Vector3 sij,rij,f3i,f3j,f3k,f2;
00074     Vector3 rg[NNM],rt[NNM];
00075     double rr[NNM],ri[NNM],rdi[NNM],rdi2[NNM],xpon[NNM],xpon2[NNM];
00076     double rri6[NNM];
00077     int list2[NNM];
00078     double r2ij,rrij;
00079     double ang,tm1,tm2,tm3,dhij,dhik,dhmu,tmij,tmik,tm2ij,tm2ik,
00080         tm3ik,tm3ij,eterm,tote3,au,pplm,ff2,eterm2,tote2;
00081     int n0, n1;
00082     int group_01;
00083     
00084     DUMP("SW");
00085 
00086     _SKIN=_RLIST-acut;
00087     refreshneighborlist();
00088 
00089     /* _EPOT, _F, _EPOT_IND, _VIRIAL all local */
00090     _EPOT=0;tote2=tote3=0;
00091     for(i=0;i<_NP;i++)
00092     { _F[i].clear(); _EPOT_IND[i]=0;}
00093     _VIRIAL.clear();
00094 
00095     n0=0;
00096     n1=_NP;
00097 
00098     for(ipt=n0;ipt<n1;ipt++)
00099     {        
00100         /* if fixed == -1, simply ignore this atom */
00101         if(fixed[ipt]==-1) continue;
00102         neg=0;
00103         for(j=0;j<nn[ipt];j++)
00104         {
00105             jpt=nindex[ipt][j];
00106             /* if fixed == -1, simply ignore this atom */
00107             if(fixed[jpt]==-1) continue;
00108             
00109             /*sij=_SR[jpt]-_SR[ipt];*/
00110             sij.subtract(_SR[jpt],_SR[ipt]);
00111             sij.subint();
00112             /*rij=_H*sij;*/
00113             _H.multiply(sij,rij);
00114             
00115             r2ij=rij.norm2();
00116             rrij=sqrt(r2ij);
00117             if(r2ij>acutsq) continue;
00118             rg[neg]=rij;
00119             /*rt[neg]=rij/rrij;*/
00120             rt[neg]=rij; rt[neg]/=rrij;
00121             rr[neg]=rrij;
00122             ri[neg]=1./rrij;
00123             rdi[neg]=1./(rrij-acut);
00124             rdi2[neg]=rdi[neg]*rdi[neg];
00125 
00126             rri6[neg]=1./(r2ij*r2ij*r2ij);
00127             
00128             if(fabs(pgam*rdi[neg])>30) xpon[neg]=0;/* avoid denormalization*/
00129             else xpon[neg]=exp(pgam*rdi[neg]);
00130             if(fabs(pss*rdi[neg])>30) xpon2[neg]=0;/* avoid denormalization*/
00131             else xpon2[neg]=exp(pss*rdi[neg]);
00132             list2[neg]=jpt;
00133             neg++;
00134         }
00135         /* second inner loop */
00136         for(j=0;j<neg;j++)
00137         {
00138             /* three body loop */
00139             jpt=list2[j];
00140             /* if fixed == -1, simply ignore this atom */
00141             if(fixed[jpt]==-1) continue;
00142 
00143             /* between group 0 and 1, no three-body interaction*/
00144             if(((group[ipt]==0) && (group[jpt]==1))||
00145                ((group[ipt]==1) && (group[jpt]==0)))
00146                 group_01 = 1;
00147             else
00148                 group_01 = 0;
00149 
00150 
00151             if(!group_01)
00152             for(k=j+1;k<neg;k++)
00153             {
00154                 kpt=list2[k];
00155                 /* if fixed == -1, simply ignore this atom */
00156                 if(fixed[kpt]==-1) continue;
00157 
00158                 /* between group 0 and 1, no three-body interaction*/
00159                 if(((group[ipt]==0) && (group[kpt]==1))||
00160                    ((group[ipt]==1) && (group[kpt]==0)))
00161                     continue;
00162                 
00163                 ang=rg[j]*rg[k]/(rr[j]*rr[k]);
00164                 tm1=ang+1./3.;
00165                 tm2=plam*pgam*tm1*tm1;
00166                 tm3=xpon[j]*xpon[k];
00167 
00168                 dhij=-1.*tm2*rdi[j]*rdi[j]*tm3;
00169                 dhik=-1.*tm2*rdi[k]*rdi[k]*tm3;
00170                 dhmu=2.*plam*tm3*tm1;
00171                     
00172                 tmij=dhij+dhmu*(ri[k]-ang*ri[j]);
00173                 tmik=dhik+dhmu*(ri[j]-ang*ri[k]);
00174                 tm2ij=-dhij+dhmu*ang*ri[j];
00175                 tm2ik=-dhmu*ri[j];
00176                 tm3ik=-dhik+dhmu*ang*ri[k]; 
00177                 tm3ij=-dhmu*ri[k];
00178 
00179                 /*f3i=rt[j]*tmij+rt[k]*tmik;*/
00180                 f3i.clear();f3i.addnv(tmij,rt[j]);f3i.addnv(tmik,rt[k]);
00181                 _F[ipt]+=f3i;
00182                     
00183                 /*f3j=rt[j]*tm2ij+rt[k]*tm2ik;*/
00184                 f3j.clear();f3j.addnv(tm2ij,rt[j]);f3j.addnv(tm2ik,rt[k]);
00185                 _F[jpt]+=f3j;
00186 
00187                 /*f3k=rt[k]*tm3ik+rt[j]*tm3ij;*/
00188                 f3k.clear();f3k.addnv(tm3ik,rt[k]);f3k.addnv(tm3ij,rt[j]);
00189                 _F[kpt]+=f3k;
00190                 if(fixedatomenergypartition==0)
00191                 {
00192                     _VIRIAL.addnvv(1.,f3j,rg[j]);
00193                     _VIRIAL.addnvv(1.,f3k,rg[k]);
00194                 }
00195                 else
00196                 {
00197                     if(!(fixed[ipt]||fixed[jpt]))
00198                         _VIRIAL.addnvv(1.,f3j,rg[j]);
00199                     else if(!(fixed[ipt]&&fixed[jpt]))
00200                         _VIRIAL.addnvv(0.5,f3j,rg[j]);
00201                     if(!(fixed[ipt]||fixed[kpt]))
00202                         _VIRIAL.addnvv(1.,f3k,rg[k]);
00203                     else if(!(fixed[ipt]&&fixed[kpt]))
00204                         _VIRIAL.addnvv(0.5,f3k,rg[k]);
00205                 }
00206                 /*_VIRIAL.addnvv(1.,f3k,rg[k]);*/
00207 
00208                 eterm=tm2*tm3/pgam;
00209 
00210                 if(fixedatomenergypartition==0)
00211                 {
00212                     tote3+=eterm;
00213                 }
00214                 else
00215                 {
00216                     if(!(fixed[ipt]&&fixed[jpt]&&fixed[kpt]))
00217                         tote3+=eterm;
00218 //                    tote3+=eterm*(1.0-(fixed[ipt]+fixed[jpt]+fixed[kpt])/3.0);
00219                     /*tote3+=eterm;*/
00220                 }
00221                 _EPOT_IND[ipt]+=eterm/3;
00222                 _EPOT_IND[jpt]+=eterm/3;
00223                 _EPOT_IND[kpt]+=eterm/3;
00224             }
00225 
00226             /* two body terms */            
00227             if(!Bond(ipt,jpt)) continue;
00228 
00229             if(group_01)
00230             {
00231                 eterm2=(ALJ*rri6[j]-BLJ)*rri6[j]-rr[j]*DUDRc+(LJ_RC*DUDRc-Uc);
00232                 f2.clear();   
00233                 f2=rt[j]*(((12*ALJ*rri6[j]-6*BLJ)*rri6[j])*ri[j]+DUDRc);
00234             }
00235             else
00236             {
00237                 au=aa*xpon2[j];
00238                 pplm=pss*(bb*pow(ri[j],rho)-1.);
00239                 eterm2=aa*(bb*pow(ri[j],rho)-1.)*xpon2[j];
00240                 ff2=au*(rho*bb*pow(ri[j],rho1)+pplm*rdi2[j]);
00241                 f2.clear();f2.addnv(ff2,rt[j]);
00242             }
00243             
00244             _F[ipt]-=f2;
00245             _F[jpt]+=f2;
00246 
00247 
00248             
00249             if(fixedatomenergypartition==0)
00250             {
00251                 tote2+=eterm2;
00252             }
00253             else
00254             {
00255                 if(!(fixed[ipt]&&fixed[jpt]))
00256                     tote2+=eterm2;
00257             }
00258 //            tote2+=eterm2*(1.0-(fixed[ipt]+fixed[jpt])/2.0);
00259             /*tote2+=eterm2;*/
00260             _EPOT_IND[ipt]+=eterm2/2;
00261             _EPOT_IND[jpt]+=eterm2/2;
00262             
00263             if(fixedatomenergypartition==0)
00264             {
00265                 _VIRIAL.addnvv(1.,f2,rg[j]);
00266             }
00267             else
00268             {
00269                 if(!(fixed[ipt]||fixed[jpt]))
00270                     _VIRIAL.addnvv(1.,f2,rg[j]);
00271                 else if(!(fixed[ipt]&&fixed[jpt]))
00272                     _VIRIAL.addnvv(0.5,f2,rg[j]);
00273                 /*_VIRIAL.addnvv(1.,f2,rg[j]);*/
00274             }
00275         }
00276     }          
00277     _EPOT=tote2+tote3;
00278 
00279 //    _EPOT_TIP=tipe2+tipe3;
00280 //    _EPOT_SUB=sube2+sube3;
00281     
00282     tipforce.clear();    
00283     for(i=0;i<_NP;i++)
00284     {
00285         if((group[i]==0)||(group[i]==4))
00286             tipforce+=_F[i];
00287     }
00288 
00289     
00290     /* zero out the forces on fixed atoms */
00291     for(i=0;i<_NP;i++)
00292     {
00293         if(fixed[i])
00294             _F[i].clear();
00295     }
00296 
00297     DUMP("SW complete");
00298 }
00299 
00300 void SWFrame::stillinger_weber_energyonly()
00301 {
00302     ERROR("SWFrame: SWLJ energy only is not implemented!");
00303 }
00304 
00305 double SWFrame::stillinger_weber_energyonly(int iatom)
00306 {             
00307     ERROR("SWFrame: SWLJ energy only is not implemented!");
00308     return 0;
00309 }
00310 
00311 void SWFrame::potential()
00312 {
00313     stillinger_weber();
00314 }
00315 
00316 void SWFrame::potential_energyonly()
00317 {
00318     stillinger_weber_energyonly();
00319 }
00320 
00321 double SWFrame::potential_energyonly(int iatom)
00322 {
00323     return stillinger_weber_energyonly(iatom);
00324 }
00325 
00326 #ifdef _TEST
00327 
00328 /* Main Program Begins */
00329 class SWFrame sim;
00330 
00331 /* The main program is defined here */
00332 #include "main.cpp"
00333 
00334 #endif//_TEST
00335