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

00001 /*
00002   sw.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Tue Jul 31 14:22:32 2007
00005 
00006   FUNCTION  :  MD simulation package of Si using Stillinger-Weber potential
00007 */
00008 
00009 #include "sw.h"
00010 
00011 void SWFrame::initparser()
00012 {
00013     MDPARALLELFrame::initparser();
00014 
00015     /* input */
00016     bindvar("sw3b_multiply",&_SW3B_MUL,DOUBLE);
00017     bindvar("sw2b_multiply",&_SW2B_MUL,DOUBLE);
00018     bindvar("tote2",&tote2,DOUBLE);
00019     bindvar("tote3",&tote3,DOUBLE);
00020 }
00021 
00022 void SWFrame::initvars()
00023 {
00024     /* convresion between original PRB 31, 5262 (1985) and modified PRB 46,2250 (1992)
00025        Rc = acut
00026        mu = pss
00027        A1 = aa*bb
00028        A2 = aa
00029        lambda1 = rho
00030        lambda2 = 0
00031        Z = plam
00032        alpha = pgam/pss
00033        theta0 = acos(1./3.)
00034     */
00035 
00036 #ifdef _SW_Si
00037 #ifdef _SW_ORIG    
00038     /* original Si version PRB 31, 5262 (1985) */
00039        aa=15.27991323; bb=11.60319228; plam=45.51575;
00040        pgam=2.51412; acut=3.77118; pss=2.0951; rho=4.0;
00041 #else   
00042     /* modified Si parameters PRB 46, 2250 (1992) */
00043        aa=16.31972277; bb=11.60319228; plam=48.61499998;
00044        pgam=2.51412;  acut=3.77118; pss=2.0951; rho=4.;
00045 #endif    
00046 #endif
00047 
00048 #ifdef _SW_Ge       
00049    /* Ge version Ding and Anderson, PRB 34, 6987 (1986) */
00050        aa=13.60564361; bb=13.62639971; plam=59.83;
00051        pgam=2.6172; acut=3.9258; pss=2.181; rho=4.0; 
00052 #endif
00053 
00054        _RLIST=acut*1.1;
00055        /*_RLIST=3.8;*/
00056        /*_RLIST=4.2;*/   /* 1.1*acut */
00057        /*_RLIST=5.6;*/   /* large skin to avoid reconstruction in Free energy calculation */
00058        _SKIN=_RLIST-acut;
00059        rho1=rho+1; acutsq=acut*acut;
00060        
00061        strcpy(incnfile,"../si.cn");
00062        DUMP(HIG"SWFrame initvars"NOR);
00063        MDPARALLELFrame::initvars();
00064 }
00065 
00066 
00067 #ifdef _TORSION_OR_BENDING
00068 #include "../cookies/src/sw_torsion_bending.cpp"
00069 #else
00070 
00071 void SWFrame::stillinger_weber()
00072 {
00073 #ifndef NNM /* need to be taken care of later */
00074 #define NNM 60
00075 #endif
00076     /* Multi process function */
00077     int i,j,k,ipt,jpt,kpt;
00078     int neg;
00079     Vector3 sij,rij,f3i,f3j,f3k,f2;
00080     Vector3 rg[NNM],rt[NNM];
00081     double rr[NNM],ri[NNM],rdi[NNM],rdi2[NNM],xpon[NNM],xpon2[NNM];
00082     int list2[NNM];
00083     double r2ij,rrij;
00084     double ang,tm1,tm2,tm3,dhij,dhik,dhmu,tmij,tmik,tm2ij,tm2ik,
00085         tm3ik,tm3ij,eterm,au,pplm,ff2,eterm2;
00086     int n0, n1;
00087 
00088     DUMP("SW");
00089 
00090     _SKIN=_RLIST-acut;
00091     refreshneighborlist();
00092 
00093     /* _EPOT, _F, _EPOT_IND, _VIRIAL all local */
00094     _EPOT=0;tote2=tote3=0;
00095     for(i=0;i<_NP;i++)
00096     { _F[i].clear(); _EPOT_IND[i]=0; _EPOT_RMV[i]=0; _VIRIAL_IND[i].clear(); }
00097     _VIRIAL.clear();
00098     
00099 #ifdef _TORSION_OR_BENDING
00100     //if (_NIMAGES>0)
00101     {
00102         if (_TORSIONSIM) _TORQUE = 0;
00103         if (_BENDSIM)    _BENDMOMENT = 0;
00104     }
00105 #endif
00106     
00107     n0=0;
00108     n1=_NP;
00109     
00110     for(ipt=n0;ipt<n1;ipt++)
00111     {
00112         /* if fixed == -1, simply ignore this atom */
00113         if(fixed[ipt]==-1) continue;
00114         neg=0;
00115         for(j=0;j<nn[ipt];j++)
00116         {
00117             jpt=nindex[ipt][j];
00118             /* if fixed == -1, simply ignore this atom */
00119             if(fixed[jpt]==-1) continue;
00120             
00121             /*sij=_SR[jpt]-_SR[ipt];*/
00122             sij.subtract(_SR[jpt],_SR[ipt]);
00123             sij.subint();
00124             /*rij=_H*sij;*/
00125             _H.multiply(sij,rij);
00126             
00127             r2ij=rij.norm2();
00128             rrij=sqrt(r2ij);
00129             if(r2ij>acutsq) continue;
00130             rg[neg]=rij;
00131             /*rt[neg]=rij/rrij;*/
00132             rt[neg]=rij; rt[neg]/=rrij;
00133             rr[neg]=rrij;
00134             ri[neg]=1./rrij;
00135             rdi[neg]=1./(rrij-acut);
00136             rdi2[neg]=rdi[neg]*rdi[neg];
00137             if(fabs(pgam*rdi[neg])>30) xpon[neg]=0;/* avoid denormalization*/
00138             else xpon[neg]=exp(pgam*rdi[neg]);
00139             if(fabs(pss*rdi[neg])>30) xpon2[neg]=0;/* avoid denormalization*/
00140             else xpon2[neg]=exp(pss*rdi[neg]);
00141             list2[neg]=jpt;
00142             neg++;
00143         }
00144 
00145         /*
00146           for(j=0;j<neg;j++)
00147           INFO_Printf("xpon[%d]=%e rdi=%e pgam=%e\n",j,xpon[j],rdi[j],pgam);
00148         */
00149         
00150         /* second inner loop */
00151         for(j=0;j<neg;j++)
00152         {
00153             /* three body loop */
00154             jpt=list2[j];
00155             /* if fixed == -1, simply ignore this atom */
00156             if(fixed[jpt]==-1) continue;
00157             
00158             if(_SW3B_MUL!=0)
00159             {
00160                 for(k=j+1;k<neg;k++)
00161                 {
00162                     kpt=list2[k];
00163                     /* if fixed == -1, simply ignore this atom */
00164                     if(fixed[kpt]==-1) continue;
00165                     ang=rg[j]*rg[k]/(rr[j]*rr[k]);
00166                     tm1=ang+1./3.;
00167                     tm2=plam*pgam*tm1*tm1;
00168                     tm3=xpon[j]*xpon[k];
00169                     
00170                     dhij=-1.*tm2*rdi[j]*rdi[j]*tm3;
00171                     dhik=-1.*tm2*rdi[k]*rdi[k]*tm3;
00172                     dhmu=2.*plam*tm3*tm1;
00173                     
00174                     tmij=dhij+dhmu*(ri[k]-ang*ri[j]);
00175                     tmik=dhik+dhmu*(ri[j]-ang*ri[k]);
00176                     tm2ij=-dhij+dhmu*ang*ri[j];
00177                     tm2ik=-dhmu*ri[j];
00178                     tm3ik=-dhik+dhmu*ang*ri[k]; 
00179                     tm3ij=-dhmu*ri[k];
00180                     
00181                     /*f3i=rt[j]*tmij+rt[k]*tmik;*/
00182                     f3i.clear();f3i.addnv(tmij,rt[j]);f3i.addnv(tmik,rt[k]);
00183                     f3i *= _SW3B_MUL; /* multiplication factor for 3-body term */
00184                     _F[ipt]+=f3i;
00185                     
00186                     /*f3j=rt[j]*tm2ij+rt[k]*tm2ik;*/
00187                     f3j.clear();f3j.addnv(tm2ij,rt[j]);f3j.addnv(tm2ik,rt[k]);
00188                     f3j *= _SW3B_MUL; /* multiplication factor for 3-body term */
00189                     _F[jpt]+=f3j;
00190                     
00191                     /*f3k=rt[k]*tm3ik+rt[j]*tm3ij;*/
00192                     f3k.clear();f3k.addnv(tm3ik,rt[k]);f3k.addnv(tm3ij,rt[j]);
00193                     f3k *= _SW3B_MUL; /* multiplication factor for 3-body term */
00194                     _F[kpt]+=f3k;
00195                     
00196                     if(fixedatomenergypartition==0)
00197                     {
00198                         _VIRIAL.addnvv(1.,f3j,rg[j]);
00199                         _VIRIAL.addnvv(1.,f3k,rg[k]);
00200                         
00201                         _VIRIAL_IND[ipt].addnvv(0.5,f3j,rg[j]);
00202                         _VIRIAL_IND[jpt].addnvv(0.5,f3j,rg[j]);
00203                         _VIRIAL_IND[ipt].addnvv(0.5,f3k,rg[k]);
00204                         _VIRIAL_IND[kpt].addnvv(0.5,f3k,rg[k]);
00205                         
00206                     }
00207                     else
00208                     {
00209                         if(!(fixed[ipt]||fixed[jpt]))
00210                             _VIRIAL.addnvv(1.,f3j,rg[j]);
00211                         else if(!(fixed[ipt]&&fixed[jpt]))
00212                             _VIRIAL.addnvv(0.5,f3j,rg[j]);
00213                         if(!(fixed[ipt]||fixed[kpt]))
00214                             _VIRIAL.addnvv(1.,f3k,rg[k]);
00215                         else if(!(fixed[ipt]&&fixed[kpt]))
00216                             _VIRIAL.addnvv(0.5,f3k,rg[k]);
00217                     }
00218                     /*_VIRIAL.addnvv(1.,f3k,rg[k]);*/
00219                     
00220                     eterm=tm2*tm3/pgam;
00221                     eterm *= _SW3B_MUL; /* multiplication factor for 3-body term */
00222                     
00223                     if(fixedatomenergypartition==0)
00224                     {
00225                         tote3+=eterm;
00226                     }
00227                     else
00228                     {
00229                         if(!(fixed[ipt]&&fixed[jpt]&&fixed[kpt]))
00230                             tote3+=eterm;
00231                     }
00232                     
00233                     _EPOT_IND[ipt]+=eterm/3;
00234                     _EPOT_IND[jpt]+=eterm/3;
00235                     _EPOT_IND[kpt]+=eterm/3;
00236                     
00237                     _EPOT_RMV[ipt]+=eterm;
00238                     _EPOT_RMV[jpt]+=eterm;
00239                     _EPOT_RMV[kpt]+=eterm;                
00240                 }
00241             }
00242             /* two body terms */
00243             if(!Bond(ipt,jpt)) continue;
00244             au=aa*xpon2[j];
00245             pplm=pss*(bb*pow(ri[j],rho)-1.);
00246             ff2=au*(rho*bb*pow(ri[j],rho1)+pplm*rdi2[j]);
00247             /*f2=rt[j]*ff2;*/
00248             f2.clear();f2.addnv(ff2,rt[j]);
00249             f2 *= _SW2B_MUL; /* multiplication factor for 2-body term */                
00250             
00251             _F[ipt]-=f2;
00252             _F[jpt]+=f2;
00253             
00254             eterm2=aa*(bb*pow(ri[j],rho)-1.)*xpon2[j];
00255             eterm2 *= _SW2B_MUL; /* multiplication factor for 2-body term */
00256             
00257             if(fixedatomenergypartition==0)
00258             {
00259                 tote2+=eterm2;
00260             }
00261             else
00262             {
00263                 if(!(fixed[ipt]&&fixed[jpt]))
00264                     tote2+=eterm2;
00265             }
00266             _EPOT_IND[ipt]+=eterm2/2;
00267             _EPOT_IND[jpt]+=eterm2/2;
00268             
00269             _EPOT_RMV[ipt]+=eterm2;
00270             _EPOT_RMV[jpt]+=eterm2;            
00271             
00272             if(fixedatomenergypartition==0)
00273             {
00274                 _VIRIAL.addnvv(1.,f2,rg[j]);
00275                 _VIRIAL_IND[ipt].addnvv(0.5,f2,rg[j]);
00276                 _VIRIAL_IND[jpt].addnvv(0.5,f2,rg[j]);
00277             }
00278             else
00279             {
00280                 if(!(fixed[ipt]||fixed[jpt]))
00281                     _VIRIAL.addnvv(1.,f2,rg[j]);
00282                 else if(!(fixed[ipt]&&fixed[jpt]))
00283                     _VIRIAL.addnvv(0.5,f2,rg[j]);
00284             }
00285         }  
00286     }
00287     _EPOT=tote2+tote3;
00288     /* zero out the forces on fixed atoms */
00289     for(i=0;i<_NP;i++)
00290     {
00291         if(fixed[i])
00292         {
00293             _F0[i]=_F[i];
00294             if(conj_fixdir==0)
00295                 _F[i].clear();
00296             else if(conj_fixdir==1)
00297                 _F[i].x=0;
00298             else if(conj_fixdir==2)
00299                 _F[i].y=0;
00300             else if(conj_fixdir==3)
00301                 _F[i].z=0;
00302         }
00303     }
00304     
00305 //#define CONSTRAIN_SUBLATTICE
00306 #ifdef CONSTRAIN_SUBLATTICE
00307     /* apply constraints for ideal shear strength calculation*/
00308     Vector3 favg0, favg1;
00309     n0 = 0; n1 = 0; favg0.clear(); favg1.clear();
00310     
00311     for(i=0;i<_NP;i++)
00312     {
00313         if(species[i]==0)
00314         {
00315             favg0+=_F[i]; n0++;
00316         }
00317         else
00318         {
00319             favg1+=_F[i]; n1++;
00320         }
00321     }
00322     favg0/=n0; favg1/=n1;
00323     for(i=0;i<_NP;i++)
00324     {
00325         if(species[i]==0)
00326             _F[i]=favg0;
00327         else
00328             _F[i]=favg1;
00329     }
00330 #endif
00331     
00332 //#define TESTVIRIAL
00333 #ifdef TESTVIRIAL
00334     SHtoR();
00335     for(i=0;i<_NP;i++)
00336     {
00337         _TOPOL[i] = dot(_R[i],_F[i]);
00338     }
00339 #endif
00340     
00341     DUMP("SW complete");
00342 }
00343 #endif
00344 
00345 
00346 void SWFrame::SWITCHpotential_user(double lambda)
00347 {
00348     stillinger_weber_switch_3body(lambda);
00349 }
00350 
00351 void SWFrame::stillinger_weber_switch_3body(double lambda)
00352 {
00353     _SW3B_MUL = lambda;
00354     stillinger_weber();
00355     dEdlambda = tote3;
00356     _SW3B_MUL = 1;
00357 }
00358 
00359 
00360 void SWFrame::stillinger_weber_energyonly()
00361 {             
00362     /* Multi process function */
00363     int i,j,k,ipt,jpt,kpt;
00364     int neg;
00365     Vector3 sij,rij;
00366     Vector3 rg[NNM],rt[NNM];
00367     double rr[NNM],ri[NNM],rdi[NNM],rdi2[NNM],xpon[NNM],xpon2[NNM];
00368     int list2[NNM];
00369     double r2ij,rrij;
00370     double ang,tm1,tm2,tm3,eterm,eterm2;
00371     int n0, n1;
00372     
00373     DUMP("SW");
00374     
00375     refreshneighborlist();
00376     
00377     /* _EPOT, _F, _EPOT_IND, _VIRIAL all local */
00378     _EPOT=0;tote2=tote3=0;
00379     for(i=0;i<_NP;i++)
00380     { _F[i].clear(); _EPOT_IND[i]=0; _EPOT_RMV[i]=0; }
00381     _VIRIAL.clear();
00382 
00383     n0=0;
00384     n1=_NP;
00385     
00386     for(ipt=n0;ipt<n1;ipt++)
00387     {
00388         /* if fixed == -1, simply ignore this atom */
00389         if(fixed[ipt]==-1) continue;
00390         neg=0;
00391         for(j=0;j<nn[ipt];j++)
00392         {
00393             jpt=nindex[ipt][j];
00394             /* if fixed == -1, simply ignore this atom */
00395             if(fixed[jpt]==-1) continue;
00396             
00397             /*sij=_SR[jpt]-_SR[ipt];*/
00398             sij.subtract(_SR[jpt],_SR[ipt]);
00399             sij.subint();
00400             /*rij=_H*sij;*/
00401             _H.multiply(sij,rij);
00402             
00403             r2ij=rij.norm2();
00404             rrij=sqrt(r2ij);
00405             if(r2ij>acutsq) continue;
00406             rg[neg]=rij;
00407             /*rt[neg]=rij/rrij;*/
00408             rt[neg]=rij; rt[neg]/=rrij;
00409             rr[neg]=rrij;
00410             ri[neg]=1./rrij;
00411             rdi[neg]=1./(rrij-acut);
00412             rdi2[neg]=rdi[neg]*rdi[neg];
00413             if(fabs(pgam*rdi[neg])>30) xpon[neg]=0;/* avoid denormalization*/
00414             else xpon[neg]=exp(pgam*rdi[neg]);
00415             if(fabs(pss*rdi[neg])>30) xpon2[neg]=0;/* avoid denormalization*/
00416             else xpon2[neg]=exp(pss*rdi[neg]);
00417             list2[neg]=jpt;
00418             neg++;
00419         }
00420         /* second inner loop */
00421         for(j=0;j<neg;j++)
00422         {
00423             /* three body loop */
00424             jpt=list2[j];
00425             /* if fixed == -1, simply ignore this atom */
00426             if(fixed[jpt]==-1) continue;
00427             
00428             if(_SW3B_MUL!=0)
00429             {
00430                 for(k=j+1;k<neg;k++)
00431                 {
00432                     kpt=list2[k];
00433                     /* if fixed == -1, simply ignore this atom */
00434                     if(fixed[kpt]==-1) continue;
00435                     ang=rg[j]*rg[k]/(rr[j]*rr[k]);
00436                     tm1=ang+1./3.;
00437                     tm2=plam*pgam*tm1*tm1;
00438                     tm3=xpon[j]*xpon[k];
00439                     
00440                     eterm=tm2*tm3/pgam;
00441                     eterm *= _SW3B_MUL; /* multiplication factor for 3-body term */
00442                     
00443                     if(fixedatomenergypartition==0)
00444                     {
00445                         tote3+=eterm;
00446                     }
00447                     else
00448                     {
00449                         if(!(fixed[ipt]&&fixed[jpt]&&fixed[kpt]))
00450                             tote3+=eterm;
00451                     }
00452                     _EPOT_IND[ipt]+=eterm/3;
00453                     _EPOT_IND[jpt]+=eterm/3;
00454                     _EPOT_IND[kpt]+=eterm/3;
00455                     
00456                     _EPOT_RMV[ipt]+=eterm;
00457                     _EPOT_RMV[jpt]+=eterm;
00458                     _EPOT_RMV[kpt]+=eterm;
00459                     
00460                 }
00461             }
00462             /* two body terms */
00463             if(!Bond(ipt,jpt)) continue;
00464             
00465             eterm2=aa*(bb*pow(ri[j],rho)-1.)*xpon2[j];
00466             
00467             
00468             if(fixedatomenergypartition==0)
00469             {
00470                 tote2+=eterm2;
00471             }
00472             else
00473             {
00474                 if(!(fixed[ipt]&&fixed[jpt]))
00475                     tote2+=eterm2;
00476             }
00477             _EPOT_IND[ipt]+=eterm2/2;
00478             _EPOT_IND[jpt]+=eterm2/2;
00479             
00480             _EPOT_RMV[ipt]+=eterm2;
00481             _EPOT_RMV[jpt]+=eterm2;
00482         }
00483     }          
00484     _EPOT=tote2+tote3;
00485     
00486 }
00487 
00488 double SWFrame::stillinger_weber_energyonly(int iatom)
00489 {             
00490     /* Multi process function */
00491     int n0,n1,j,k,ipt,jpt,kpt;
00492     int neg;
00493     Vector3 sij,rij;
00494     Vector3 rg[NNM],rt[NNM];
00495     double rr[NNM],ri[NNM],rdi[NNM],rdi2[NNM],xpon[NNM],xpon2[NNM];
00496     int list2[NNM];
00497     double r2ij,rrij;
00498     double ang,tm1,tm2,tm3,eterm,eterm2;
00499     double Eatom; //, maxd;
00500     int relevant;
00501     
00502     if(iatom<_NP)
00503     {
00504         refreshneighborlist();
00505     }
00506 
00507     Eatom=0; tote2=tote3=0; _EPOT_IND[iatom]=0; _EPOT_RMV[iatom]=0;
00508 
00509     n0=0;
00510     n1=_NP;
00511 
00512     if(iatom>=_NP) n1=_NP+1;     /* test particle */
00513     
00514     for(ipt=n0;ipt<n1;ipt++)
00515     {
00516         if(ipt>=_NP) ipt=iatom; /* test particle */
00517         
00518         /* if fixed == -1, simply ignore this atom */
00519         if(fixed[ipt]==-1) continue;
00520         
00521         relevant=0;
00522         if(ipt==iatom) relevant=1;
00523 #if 1 /* also need to consider 3-body terms when iatom is not at the center */
00524         if(!relevant)
00525         {
00526             for(j=0;j<nn[ipt];j++)
00527             {
00528                 jpt=nindex[ipt][j];
00529                 if(jpt==iatom) relevant=1;
00530             }
00531         }
00532 #endif
00533         if(!relevant) continue;
00534         
00535         neg=0;
00536         for(j=0;j<nn[ipt];j++)
00537         {
00538             jpt=nindex[ipt][j];
00539             /* if fixed == -1, simply ignore this atom */
00540             if(fixed[jpt]==-1) continue;
00541             
00542             /*sij=_SR[jpt]-_SR[ipt];*/
00543             sij.subtract(_SR[jpt],_SR[ipt]);
00544             sij.subint();
00545             /*rij=_H*sij;*/
00546             _H.multiply(sij,rij);
00547             
00548             r2ij=rij.norm2();
00549             rrij=sqrt(r2ij);
00550             if(r2ij>acutsq) continue;
00551             rg[neg]=rij;
00552             /*rt[neg]=rij/rrij;*/
00553             rt[neg]=rij; rt[neg]/=rrij;
00554             rr[neg]=rrij;
00555             ri[neg]=1./rrij;
00556             rdi[neg]=1./(rrij-acut);
00557             rdi2[neg]=rdi[neg]*rdi[neg];
00558             if(fabs(pgam*rdi[neg])>30) xpon[neg]=0;/* avoid denormalization*/
00559             else xpon[neg]=exp(pgam*rdi[neg]);
00560             if(fabs(pss*rdi[neg])>30) xpon2[neg]=0;/* avoid denormalization*/
00561             else xpon2[neg]=exp(pss*rdi[neg]);
00562             list2[neg]=jpt;
00563             neg++;
00564         }
00565         /* second inner loop */
00566         for(j=0;j<neg;j++)
00567         {
00568             /* three body loop */
00569             jpt=list2[j];
00570             /* if fixed == -1, simply ignore this atom */
00571             if(fixed[jpt]==-1) continue;
00572             
00573             for(k=j+1;k<neg;k++)
00574             {
00575                 kpt=list2[k];
00576                 /* if fixed == -1, simply ignore this atom */
00577                 if(fixed[kpt]==-1) continue;
00578                 ang=rg[j]*rg[k]/(rr[j]*rr[k]);
00579                 tm1=ang+1./3.;
00580                 tm2=plam*pgam*tm1*tm1;
00581                 tm3=xpon[j]*xpon[k];
00582                 
00583                 eterm=tm2*tm3/pgam;
00584                 eterm *= _SW3B_MUL; /* multiplication factor for 3-body term */
00585                 
00586                 if(fixedatomenergypartition==0)
00587                 {
00588                     tote3+=eterm;
00589                 }
00590                 else
00591                 {
00592                     if(!(fixed[ipt]&&fixed[jpt]&&fixed[kpt]))
00593                         tote3+=eterm;
00594                 }
00595                 if(ipt==iatom) _EPOT_IND[ipt]+=eterm/3;
00596                 if(jpt==iatom) _EPOT_IND[jpt]+=eterm/3;
00597                 if(kpt==iatom) _EPOT_IND[kpt]+=eterm/3;
00598                 
00599                 if(ipt==iatom) _EPOT_RMV[ipt]+=eterm;
00600                 if(jpt==iatom) _EPOT_RMV[jpt]+=eterm;
00601                 if(kpt==iatom) _EPOT_RMV[kpt]+=eterm;
00602                 
00603             }
00604             /* two body terms */
00605             if(!Bond(ipt,jpt)) continue;
00606             
00607             eterm2=aa*(bb*pow(ri[j],rho)-1.)*xpon2[j];
00608             
00609             
00610             if(fixedatomenergypartition==0)
00611             {
00612                 tote2+=eterm2;
00613             }
00614             else
00615             {
00616                 if(!(fixed[ipt]&&fixed[jpt]))
00617                     tote2+=eterm2;
00618             }
00619             if(ipt==iatom) _EPOT_IND[ipt]+=eterm2/2;
00620             if(jpt==iatom) _EPOT_IND[jpt]+=eterm2/2;
00621             
00622             if(ipt==iatom) _EPOT_RMV[ipt]+=eterm2;
00623             if(jpt==iatom) _EPOT_RMV[jpt]+=eterm2;
00624             
00625         }
00626     }          
00627     Eatom=tote2+tote3;
00628     
00629     return Eatom;
00630 }
00631 
00632 void SWFrame::potential()
00633 {
00634     stillinger_weber();
00635 }
00636 
00637 void SWFrame::potential_energyonly()
00638 {
00639     stillinger_weber_energyonly();
00640 }
00641 
00642 double SWFrame::potential_energyonly(int iatom)
00643 {
00644     return stillinger_weber_energyonly(iatom);
00645 }
00646 
00647 #ifdef _TEST
00648 
00649 /* Main Program Begins */
00650 class SWFrame sim;
00651 
00652 /* The main program is defined here */
00653 #include "main.cpp"
00654 
00655 #endif//_TEST
00656