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

00001 /*
00002   eam.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Wed Oct  8 22:42:59 2008
00005 
00006   FUNCTION  : MD++ implementation of EAM/MEAM potential
00007 */
00008 
00009 #include "eam.h"
00010 
00011 #define _CUBICSPLINE
00012 
00013 inline double spline(double spline_coeff[NGRID][4],int ind, double qq)
00014 {
00015   double a, b, c, d, qq2, qq3, f;
00016   a = spline_coeff[ind][0];
00017   b = spline_coeff[ind][1];
00018   c = spline_coeff[ind][2];
00019   d = spline_coeff[ind][3];
00020   qq2=qq*qq; qq3=qq2*qq;
00021   f=a+b*qq+c*qq2+d*qq3;
00022   return f;
00023 }  
00024 
00025 inline double spline1(double spline_coeff[NGRID][4],int ind, double qq)
00026 {
00027   double a, b, c, d, qq2, qq3, fp;
00028   a = spline_coeff[ind][0];
00029   b = spline_coeff[ind][1];
00030   c = spline_coeff[ind][2];
00031   d = spline_coeff[ind][3];
00032   qq2=qq*qq; qq3=qq2*qq;
00033   fp=b+2*c*qq+3*d*qq2;
00034   return fp;
00035 }
00036 
00037 
00038 void compute_spline_coeff(int ngrid, double func[],double deriv[],double dr,double spline_coeff[NGRID][4])
00039 {
00040   int ind;
00041   double a, b, c, d, f1, p1, A1, A2, dr2, dr3;
00042   
00043   dr2=dr*dr; dr3=dr2*dr;
00044 
00045   for(ind=0;ind<ngrid;ind++)
00046   {
00047     a = func[ind];
00048     b = deriv[ind];
00049     f1 = func[ind+1];
00050     p1 = deriv[ind+1];
00051     A1 = f1-a-b*dr;
00052     A2 = (p1-b)*dr;
00053     d = (A2-2*A1)/dr3;
00054     c = (3*A1-A2)/dr2;
00055     spline_coeff[ind][0] = a;
00056     spline_coeff[ind][1] = b;
00057     spline_coeff[ind][2] = c;
00058     spline_coeff[ind][3] = d;
00059   }
00060 }
00061 
00062 int EAMFrame::readeam()
00063 {
00064     FILE *fp;
00065     char string[500];
00066     int i;
00067     
00068     LFile::SubHomeDir(potfile,potfile);
00069     
00070     fp=fopen(potfile,"r");
00071     if(fp==NULL)
00072     {
00073         FATAL("EAMFrame::readeam file ("<<potfile<<") not found!");
00074         return -1;
00075     }
00076     fgets(title3,100,fp);
00077     INFO_Printf("%s\n",title3);
00078 
00079     fgets(string,500,fp);
00080     sscanf(string,"%d %lf %lf",&ntyp,&rmass,&rlatt);
00081     fgets(string,500,fp);
00082     sscanf(string,"%le %le %le %le",&drar,&drhoar,&actual,&rmin);
00083 
00084     /*debug */
00085     actual*=0.995;
00086     
00087     actual2 = actual*actual;
00088 
00089     INFO_Printf("ntyp=%d rmass=%f rlatt=%f\n",ntyp,rmass,rlatt);
00090     INFO_Printf("dr=%e drho=%e rcut=%e rmin=%e\n",drar,drhoar,actual,rmin);
00091 
00092     _RLIST = 1.1 * actual;
00093     _SKIN = _RLIST - actual;
00094 
00095     for(i=0;i<eamgrid;i++)
00096     {
00097       fgets(string,500,fp);
00098       sscanf(string,"%le %le %le %le",&(rho[0][i]),&(rhop[0][i]),&(rho[1][i]),&(rhop[1][i]));
00099     }
00100     for(i=0;i<eamgrid;i++)
00101     {
00102       fgets(string,500,fp);
00103       sscanf(string,"%le %le %le %le",&(phi[0][i]),&(phip[0][i]),&(phi[1][i]),&(phip[1][i]));
00104     }
00105     if(eamfiletype==2) 
00106     {
00107         for(i=0;i<eamgrid;i++) 
00108         {
00109           fgets(string,500,fp);
00110           sscanf(string,"%le %le",&(phix[i]),&(phipx[i]));
00111         }
00112     }
00113     for(i=0;i<eamgrid;i++)
00114     {
00115       fgets(string,500,fp);
00116       sscanf(string,"%le %le %le %le",&(frho[0][i]),&(frhop[0][i]),&(frho[1][i]),&(frhop[1][i]));
00117     }  
00118     INFO("Finished reading EAM potential");
00119     
00120     /*   Initialize linear grid points reduntantly on each node to save
00121      *   transfer time.
00122      */
00123     for(i=0;i<eamgrid;i++)
00124     {
00125         rval[i] = i*drar;
00126         rhoval[i] = i*drhoar;
00127     }
00128     pottype = 1; /* set potential to EAM */
00129 
00130 #ifdef _CUBICSPLINE
00131     /* compute spline coefficients */
00132     compute_spline_coeff(eamgrid,rho[0],rhop[0],drar,rho_spline[0]);
00133     compute_spline_coeff(eamgrid,phi[0],phip[0],drar,phi_spline[0]);
00134     compute_spline_coeff(eamgrid,frho[0],frhop[0],drhoar,frho_spline[0]);
00135 
00136     if(eamfiletype==2)
00137     {
00138       compute_spline_coeff(eamgrid,rho[1],rhop[1],drar,rho_spline[1]);
00139       compute_spline_coeff(eamgrid,phi[1],phip[1],drar,phi_spline[1]);
00140       compute_spline_coeff(eamgrid,phix,  phipx,  drar,phix_spline);
00141       compute_spline_coeff(eamgrid,frho[1],frhop[1],drhoar,frho_spline[1]);
00142     }
00143 #endif
00144     return 0;
00145 }
00146 
00147 
00148 #ifdef _TORSION_OR_BENDING
00149 #include "../cookies/src/eam_torsion_bending.cpp"
00150 #else
00151 
00152 void EAMFrame::rhoeam()
00153 {
00154     int i, j, jpt, idx, jdx, ind;
00155     Vector3 sij,rij;
00156     double r2ij, rmagg, qq, qr;
00157     double rhoi, rhoj;
00158 
00159     /*INFO("rhoeam");*/
00160     petrip = 0;
00161     rhocon = 1e-10;
00162     rhomax = eamgrid*drhoar;
00163 
00164     for(i=0;i<_NP;i++)
00165     {
00166         atpe3b[i]=0;
00167         rhotot[i]=0;
00168         nbst[i]=0;
00169     }
00170     for(i=0;i<_NP;i++)
00171     {
00172         _F[i].clear(); _EPOT_IND[i]=0;
00173     }
00174     _EPOT=0;
00175     _VIRIAL.clear();
00176     
00177     /*  start calculation
00178      *  calculate pair contribution to total energy and forces
00179      */
00180 
00181     /* do on i-particles */
00182     for(i=0;i<_NP;i++)
00183     {
00184         /* modify here for binary systems (0/1) */
00185         idx = species[i]; /* type of atom (i) */
00186         /* do on j-particles */
00187         for(j=0;j<nn[i];j++)
00188         {
00189             /* modify here for binary systems (0/1) */
00190             jpt=nindex[i][j];
00191             jdx = species[jpt]; /* type of atom (j) */
00192             if(i>=jpt) continue;
00193             sij=_SR[jpt]-_SR[i];
00194             sij.subint();
00195             rij=_H*sij;
00196             r2ij=rij.norm2();
00197             if(r2ij>actual2) continue;
00198             rmagg=sqrt(r2ij)-rmin;
00199             
00200             nbst[i]++;
00201             ind = (int)(rmagg/drar);
00202 
00203             if(ind>=eamgrid)
00204             {
00205                 ind=eamgrid-1;
00206                 INFO_Printf("ind = %d r=%f in RHOEAM\n",ind, rmagg+rmin);
00207             }
00208             else if(ind<0)
00209             {
00210                 ind=0;
00211                 INFO_Printf("ind = %d r=%f in RHOEAM\n",ind, rmagg+rmin);
00212             }
00213             qq=rmagg-rval[ind];
00214 
00215             if((idx==jdx))
00216             {
00217 #ifndef _CUBICSPLINE
00218             rhoi=rho[jdx][ind]+ qq*rhop[jdx][ind];
00219 #else
00220             //rhoi = interp(rho[jdx],rhop[jdx],drar,ind,qq);
00221             rhoi = spline(rho_spline[jdx],ind,qq);
00222 #endif
00223             rhotot[i]+=rhoi;
00224             rhotot[jpt]+=rhoi;
00225             }
00226             else
00227             {
00228 #ifndef _CUBICSPLINE
00229               rhoi=rho[jdx][ind]+ qq*rhop[jdx][ind];
00230               rhoj=rho[idx][ind]+ qq*rhop[idx][ind];
00231 #else
00232               //rhoi = interp(rho[jdx],rhop[jdx],drar,ind,qq);
00233               //rhoj = interp(rho[idx],rhop[idx],drar,ind,qq);
00234               rhoi = spline(rho_spline[jdx],ind,qq);
00235               rhoj = spline(rho_spline[idx],ind,qq);
00236 #endif
00237               rhotot[i]+=rhoi;
00238               rhotot[jpt]+=rhoj;
00239             }
00240         }
00241     }
00242 
00243     /* debug */
00244 #if 0
00245     for(i=0;i<_NP;i++)
00246     {
00247         INFO_Printf("atom[%d] rho=%e\n",i,rhotot[i]);
00248     }
00249 #endif
00250     
00251     for(i=0;i<_NP;i++)
00252     {
00253         /* modify here for binary systems (0/1) */
00254         idx = species[i]; /* type of atom (i) */
00255         if(rhotot[i]<rhocon)
00256         {
00257             rhotot[i]=rhocon;
00258         }
00259         if(rhotot[i]>rhomax)
00260         {
00261             rhotot[i]=rhomax;
00262         }
00263         ind = (int)(rhotot[i]/drhoar);
00264         if(ind>=eamgrid-1) ind=eamgrid-1;
00265         qr = rhotot[i] - rhoval[ind];
00266 
00267 #ifndef _CUBICSPLINE
00268         embf[i] = frho[idx][ind] + qr*frhop[idx][ind];
00269         embfp[i] = frhop[idx][ind] +
00270           qr*(frhop[idx][ind+1]-frhop[idx][ind])/drhoar;
00271 #else
00272         //embf[i] = interp(frho[idx],frhop[idx],drhoar,ind,qr);
00273         //embfp[i] = interp1(frho[idx],frhop[idx],drhoar,ind,qr);
00274         embf[i] = spline(frho_spline[idx],ind,qr);
00275         embfp[i] = spline1(frho_spline[idx],ind,qr);
00276 #endif
00277         atpe3b[i] = embf[i];
00278         _EPOT_IND[i]+=atpe3b[i];
00279         _EPOT+=atpe3b[i];
00280     }    
00281 }
00282 
00283 void EAMFrame::kraeam()
00284 {
00285     int i, j, jpt, idx, jdx, ind;
00286     Vector3 sij,rij,fij;
00287     double r2ij, rmagg, pp, qq, fcp, fpp, fp, denspi, denspj;
00288     
00289     /*  start calculation
00290      *  calculate pair contribution to total energy and forces
00291      */
00292     
00293     /* do on i-particles */
00294     for(i=0;i<_NP;i++)
00295     {
00296         /* modify here for binary systems (0/1) */
00297         idx = species[i]; /* type of atom (i) */
00298         /* do on j-particles */
00299         for(j=0;j<nn[i];j++)
00300         {
00301             /* modify here for binary systems (0/1) */
00302             jpt=nindex[i][j];
00303             jdx = species[jpt]; /* type of atom (j) */
00304             if(i>=jpt) continue;
00305             sij=_SR[jpt]-_SR[i];
00306             sij.subint();
00307             rij=_H*sij;
00308             r2ij=rij.norm2();
00309             if(r2ij>actual2) continue;
00310             rmagg=sqrt(r2ij)-rmin;
00311 
00312             ind = int(rmagg/drar);
00313             
00314             if(ind>=eamgrid)
00315             {
00316                 ind=eamgrid-1;
00317                 INFO_Printf("ind = %d in RHOEAM\n",ind);
00318             }
00319             else if(ind<0)
00320             {
00321                 ind=0;
00322                 INFO_Printf("ind = %d in RHOEAM\n",ind);
00323             }
00324             qq=rmagg-rval[ind];
00325 
00326             if(idx==jdx)
00327             {
00328 #ifndef _CUBICSPLINE
00329               pp = phi[idx][ind] + qq*phip[idx][ind];
00330               fpp = phip[idx][ind] +
00331                 qq*(phip[idx][ind+1]-phip[idx][ind])/drar;
00332 #else
00333               //pp = interp(phi[idx],phip[idx],drar,ind,qq);
00334               //fpp = interp1(phi[idx],phip[idx],drar,ind,qq);
00335               pp = spline(phi_spline[idx],ind,qq);
00336               fpp = spline1(phi_spline[idx],ind,qq);
00337 #endif
00338             }
00339             else
00340             {
00341 #ifndef _CUBICSPLINE
00342               pp = phix[ind] + qq*phipx[ind];
00343               fpp = phipx[ind] + qq*(phipx[ind+1]-phipx[ind])/drar;
00344 #else
00345               //pp = interp(phix,phipx,drar,ind,qq);
00346               //fpp = interp1(phix,phipx,drar,ind,qq);
00347               pp = spline(phix_spline,ind,qq);
00348               fpp = spline1(phix_spline,ind,qq);
00349 #endif
00350             }
00351 #ifndef _CUBICSPLINE
00352             denspi = rhop[idx][ind] +
00353                 qq*(rhop[idx][ind+1]-rhop[idx][ind])/drar ;
00354             denspj = rhop[jdx][ind] +
00355                 qq*(rhop[jdx][ind+1]-rhop[idx][ind])/drar ;
00356 #else
00357             //denspi = interp1(rho[idx],rhop[idx],drar,ind,qq);
00358             //denspj = interp1(rho[jdx],rhop[jdx],drar,ind,qq);
00359             denspi = spline1(rho_spline[idx],ind,qq);
00360             denspj = spline1(rho_spline[jdx],ind,qq);
00361 #endif
00362 
00363             fcp = denspj * embfp[i] + denspi * embfp[jpt];
00364             fp = (fpp + fcp) / rmagg;
00365             _EPOT_IND[i]+= 0.5*pp;
00366             _EPOT_IND[jpt]+= 0.5*pp;
00367             _EPOT+=pp;
00368             
00369             fij=rij*fp;
00370             _F[i]+=fij;
00371             _F[jpt]-=fij;
00372 #if 0
00373             if(!(fixed[i]||fixed[jpt]))
00374                 _VIRIAL.addnvv(-fp,rij,rij);
00375             else if(!(fixed[i]&&fixed[jpt]))
00376                 _VIRIAL.addnvv(-0.5*fp,rij,rij);
00377 #else
00378                 _VIRIAL.addnvv(-fp,rij,rij);
00379 #endif
00380         }
00381     }
00382 }
00383 #endif
00384 
00385 
00386 double EAMFrame::rhoeam(int iatom)
00387 {
00388     /* returning the electron density at iatom */
00389     WARNING("EAMFrame::rhoeam(iatom) not implemented yet!");
00390     return 0;
00391 }
00392 
00393 double EAMFrame::kraeam_energyonly(int iatom)
00394 {
00395     /* returning the potential energy contribution involving iatom */
00396     WARNING("EAMFrame::kraeam_energyonly(iatom) not implemented yet!");
00397     return 0;
00398 }
00399 
00400 int EAMFrame::readMEAM()
00401 {
00402     pottype = 2; /* set potential to MEAM */
00403     return -1;
00404 };
00405 
00406 void EAMFrame::rhoMEAM()
00407 {
00408 
00409 }
00410 
00411 void EAMFrame::kraMEAM()
00412 {
00413 
00414 }
00415 
00416 void EAMFrame::Alloc()
00417 {
00418     MDPARALLELFrame::Alloc();
00419     int size;
00420     size=_NP*allocmultiple;
00421     
00422     /* Shared Memory */
00423     Realloc(atpe3b,double,size);
00424     Realloc(rhotot,double,size);
00425     Realloc(embf,double,size);
00426     Realloc(embfp,double,size);
00427     Realloc(nbst,int,size);
00428 }    
00429 
00430 void EAMFrame::potential()
00431 {
00432     refreshneighborlist(); // 2007.Jul.10 Chris Weinberger & Wei Cai
00433     if(pottype==1)
00434     {
00435         rhoeam();
00436         kraeam();
00437     }
00438     else if(pottype==2)
00439     {
00440         rhoMEAM();
00441         kraMEAM();
00442     }
00443     else
00444     {
00445         WARNING("Unknown pottype ("<<pottype<<")"
00446                 " no calculation performed");
00447     }
00448 }
00449 
00450 void EAMFrame::potential_energyonly()
00451 {
00452     ERROR("EAMFrame::potential_energyonly() not implemented yet!");
00453 }
00454 
00455 double EAMFrame::potential_energyonly(int iatom)
00456 {
00457     rhoeam(iatom);
00458     return kraeam_energyonly(iatom);
00459 }
00460 
00461 
00462 void EAMFrame::initvars()
00463 {
00464     MDPARALLELFrame::initvars();
00465    
00466     strcpy(potfile,"eamdata");
00467 }
00468 void EAMFrame::initparser()
00469 {
00470     MDPARALLELFrame::initparser();
00471     bindvar("eamgrid",&eamgrid,INT);
00472     bindvar("pottype",&pottype,INT);
00473     bindvar("eamfiletype",&eamfiletype,INT);
00474 }
00475 
00476 int EAMFrame::exec(char *name)
00477 {
00478     if(MDPARALLELFrame::exec(name)==0) return 0;
00479     bindcommand(name,"readeam",readeam());
00480     bindcommand(name,"readMEAM",readMEAM());
00481 
00482 #ifdef _PARALLEL
00483     bindcommand(name,"Broadcast_EAM_Param",Broadcast_EAM_Param());
00484 #endif
00485     return -1;
00486 }
00487    
00488 #ifdef _PARALLEL
00489 void EAMFrame::Broadcast_EAM_Param()
00490 {
00491     double *buf_double;
00492     int nparam, i, j, noffset1, noffset2;
00493 
00494     /* master asking slaves to call the same function */
00495     if(myDomain==0) Master_to_Slave("Broadcast_EAM_Param");    
00496 
00497     /* broadcasting */
00498     MPI_Bcast(&pottype,1,MPI_INT,0,MPI_COMM_WORLD);
00499     MPI_Bcast(&eamgrid,1,MPI_INT,0,MPI_COMM_WORLD);
00500 
00501     /* packing parameters */
00502     noffset1 = 12;
00503     noffset2 = noffset1 + 14*eamgrid;
00504     nparam = noffset2 + 7*4*eamgrid;
00505 
00506     buf_double = (double *) malloc(nparam*sizeof(double));
00507     if(myDomain==0)
00508     {
00509     buf_double[0]  = (double) ntyp;
00510     buf_double[1]  = rmass;
00511     buf_double[2]  = rlatt;
00512     buf_double[3]  = drar;
00513     buf_double[4]  = drhoar;
00514     buf_double[5]  = actual;
00515     buf_double[6]  = actual2;
00516     buf_double[7]  = rmin;
00517     buf_double[8]  = _RLIST;
00518     buf_double[9]  = _SKIN;
00519     buf_double[10] = (double) _NNM;
00520     buf_double[11] = (double) _NIC;
00521 
00522     for(i=0;i<eamgrid;i++)
00523     {
00524       buf_double[i*14 + noffset1 + 0] = rho[0][i];
00525       buf_double[i*14 + noffset1 + 1] = rhop[0][i];
00526 
00527       buf_double[i*14 + noffset1 + 2] = rho[1][i];
00528       buf_double[i*14 + noffset1 + 3] = rhop[1][i];
00529 
00530       buf_double[i*14 + noffset1 + 4] = phi[0][i];
00531       buf_double[i*14 + noffset1 + 5] = phip[0][i];
00532 
00533       buf_double[i*14 + noffset1 + 6] = phi[1][i];
00534       buf_double[i*14 + noffset1 + 7] = phip[1][i];
00535 
00536       buf_double[i*14 + noffset1 + 8] = frho[0][i];
00537       buf_double[i*14 + noffset1 + 9] = frhop[0][i];
00538 
00539       buf_double[i*14 + noffset1 +10] = frho[1][i];
00540       buf_double[i*14 + noffset1 +11] = frhop[1][i];
00541 
00542       buf_double[i*14 + noffset1 +12] = phix[i];
00543       buf_double[i*14 + noffset1 +13] = phipx[i];
00544     }
00545     
00546     for(i=0;i<eamgrid;i++)
00547        for(j=0;j<4;j++)
00548        {
00549           buf_double[(i*4+j)*7 + noffset2 + 0] = rho_spline[0][i][j];
00550           buf_double[(i*4+j)*7 + noffset2 + 1] = rho_spline[1][i][j];
00551           buf_double[(i*4+j)*7 + noffset2 + 2] = phi_spline[0][i][j];
00552           buf_double[(i*4+j)*7 + noffset2 + 3] = phi_spline[1][i][j];
00553           buf_double[(i*4+j)*7 + noffset2 + 4] = frho_spline[0][i][j];
00554           buf_double[(i*4+j)*7 + noffset2 + 5] = frho_spline[1][i][j];
00555           buf_double[(i*4+j)*7 + noffset2 + 6] = phix_spline[i][j];
00556        }
00557     }
00558     /* broadcasting */
00559     MPI_Bcast(buf_double,nparam,MPI_DOUBLE,0,MPI_COMM_WORLD);
00560 
00561     if(myDomain!=0)
00562     {
00563     /* unpacking parameters
00564      * the following lines can be generated from the above lines by
00565      * replace regexp:
00566      *  \(buf_double\[[0-9]+\]\)[ ]+= \([^;]+\);
00567      *  \2 = \1;
00568     */
00569     ntyp   = (int) buf_double[0];
00570     rmass  = buf_double[1];
00571     rlatt  = buf_double[2];
00572     drar   = buf_double[3];
00573     drhoar = buf_double[4];
00574     actual = buf_double[5];
00575     actual2= buf_double[6];
00576     rmin   = buf_double[7];
00577     _RLIST = buf_double[8];
00578     _SKIN  = buf_double[9]; 
00579     _NNM = (int) buf_double[10];
00580     _NIC = (int) buf_double[11];
00581 
00582     for(i=0;i<eamgrid;i++)
00583     {
00584       rho[0][i]  = buf_double[i*14 + noffset1 + 0];
00585       rhop[0][i] = buf_double[i*14 + noffset1 + 1];
00586 
00587       rho[1][i]  = buf_double[i*14 + noffset1 + 2];
00588       rhop[1][i] = buf_double[i*14 + noffset1 + 3];
00589 
00590       phi[0][i]  = buf_double[i*14 + noffset1 + 4];
00591       phip[0][i] = buf_double[i*14 + noffset1 + 5];
00592 
00593       phi[1][i]  = buf_double[i*14 + noffset1 + 6];
00594       phip[1][i] = buf_double[i*14 + noffset1 + 7];
00595 
00596       frho[0][i] = buf_double[i*14 + noffset1 + 8];
00597       frhop[0][i]= buf_double[i*14 + noffset1 + 9];
00598 
00599       frho[1][i] = buf_double[i*14 + noffset1 +10];
00600       frhop[1][i]= buf_double[i*14 + noffset1 +11];
00601 
00602       phix[i]    = buf_double[i*14 + noffset1 +12];
00603       phipx[i]   = buf_double[i*14 + noffset1 +13];
00604     }
00605 
00606     for(i=0;i<eamgrid;i++)
00607       for(j=0;j<4;j++)
00608       {
00609          rho_spline[0][i][j]  = buf_double[(i*4+j)*7 + noffset2 + 0];
00610          rho_spline[1][i][j]  = buf_double[(i*4+j)*7 + noffset2 + 1];
00611          phi_spline[0][i][j]  = buf_double[(i*4+j)*7 + noffset2 + 2];
00612          phi_spline[1][i][j]  = buf_double[(i*4+j)*7 + noffset2 + 3];
00613          frho_spline[0][i][j] = buf_double[(i*4+j)*7 + noffset2 + 4];
00614          frho_spline[1][i][j] = buf_double[(i*4+j)*7 + noffset2 + 5];
00615          phix_spline[i][j]    = buf_double[(i*4+j)*7 + noffset2 + 6];
00616       }
00617 
00618     }
00619     
00620     free(buf_double);
00621 
00622     for(i=0;i<eamgrid;i++)
00623     {
00624         rval[i] = i*drar;
00625         rhoval[i] = i*drhoar;
00626     }
00627 
00628 }
00629 #endif//_PARALLEL
00630 
00631 #ifdef _TEST
00632 
00633 /* Main Program Begins */
00634 class EAMFrame sim;
00635 
00636 /* The main program is defined here */
00637 #include "main.cpp"
00638 
00639 #endif//_TEST
00640 
00641 
00642 
00643 
00644 
00645 
00646 
00647 
00648 
00649 
00650 
00651 
00652 
00653 
00654 
00655 
00656 
00657 
00658 
00659 
00660 
00661 
00662 
00663 
00664 
00665 
00666 
00667 
00668 
00669 
00670 
00671 
00672 
00673 
00674 
00675 
00676 
00677 
00678 
00679 
00680 
00681 
00682 
00683 
00684 
00685 
00686 
00687 
00688 
00689 
00690 
00691 
00692 
00693 
00694 
00695 
00696 
00697 
00698 
00699 
00700 
00701 /* old code */
00702 #if 0
00703 
00704 inline double interp(double func[],double deriv[],double dr,int ind,double qq)
00705 {
00706     double f, a, b, c, d, f1, p1, A1, A2, dr2, dr3, qq2, qq3;
00707 //    f = func[ind] + qq*deriv[ind];
00708     dr2=dr*dr; dr3=dr2*dr;
00709     qq2=qq*qq; qq3=qq2*qq;
00710     a = func[ind];
00711     b = deriv[ind];
00712     f1 = func[ind+1];
00713     p1 = deriv[ind+1];
00714     A1 = f1-a-b*dr;
00715     A2 = (p1-b)*dr;
00716     d = (A2-2*A1)/dr3;
00717     c = (3*A1-A2)/dr2;
00718     f=a+b*qq+c*qq2+d*qq3;
00719     return f;
00720 }
00721     
00722     
00723 inline double interp1(double func[],double deriv[],double dr,int ind,double qq)
00724 {
00725     double fp, a, b, c, d, f1, p1, A1, A2, dr2, dr3, qq2, qq3;
00726 //    fp = deriv[ind] + qq1/dr*(deriv[ind+1]-deriv[ind]);
00727     dr2=dr*dr; dr3=dr2*dr;
00728     qq2=qq*qq; qq3=qq2*qq;
00729     a = func[ind];
00730     b = deriv[ind];
00731     f1 = func[ind+1];
00732     p1 = deriv[ind+1];
00733     A1 = f1-a-b*dr;
00734     A2 = (p1-b)*dr;
00735     d = (A2-2*A1)/dr3;
00736     c = (3*A1-A2)/dr2;
00737     fp=b+2*c*qq+3*d*qq2;
00738     return fp;
00739 }
00740 #endif