From 290625c0cf867dbcb208518ca1af59d6cf4432d9 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 1 Nov 2006 23:00:53 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@103 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MANYBODY/pair_eam.cpp | 777 ++++++++++++++--------------- src/MANYBODY/pair_eam_alloy.cpp | 538 +++++++------------- src/MANYBODY/pair_eam_alloy.h | 8 +- src/MANYBODY/pair_eam_fs.cpp | 840 ++++++++------------------------ src/MANYBODY/pair_eam_fs.h | 16 +- 5 files changed, 749 insertions(+), 1430 deletions(-) diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 651d32dde8..ac1e7e04f9 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -24,8 +24,8 @@ #include "force.h" #include "update.h" #include "comm.h" -#include "memory.h" #include "neighbor.h" +#include "memory.h" #include "error.h" #define MIN(a,b) ((a) < (b) ? (a) : (b)) @@ -40,16 +40,20 @@ PairEAM::PairEAM() nmax = 0; rho = NULL; fp = NULL; - - ntables = 0; - tables = NULL; + + nfuncfl = 0; + funcfl = NULL; + + setfl = NULL; + fs = NULL; + frho = NULL; - frho_0 = NULL; - - // set rhor to NULL so memory deallocation will work - // even from derived classes that don't use rhor - rhor = NULL; + z2r = NULL; + + frho_spline = NULL; + rhor_spline = NULL; + z2r_spline = NULL; } /* ---------------------------------------------------------------------- @@ -65,26 +69,49 @@ PairEAM::~PairEAM() if (allocated) { memory->destroy_2d_int_array(setflag); memory->destroy_2d_double_array(cutsq); - memory->destroy_2d_int_array(tabindex); + delete [] map; + delete [] type2frho; + memory->destroy_2d_int_array(type2rhor); + memory->destroy_2d_int_array(type2z2r); } - for (int m = 0; m < ntables; m++) { - delete [] tables[m].filename; - delete [] tables[m].frho; - delete [] tables[m].rhor; - delete [] tables[m].zr; - delete [] tables[m].z2r; - } - memory->sfree(tables); - - if (frho) { - memory->destroy_2d_double_array(frho); - memory->destroy_2d_double_array(rhor); - memory->destroy_2d_double_array(zrtmp); - memory->destroy_3d_double_array(z2r); + if (funcfl) { + for (int i = 0; i < nfuncfl; i++) { + delete [] funcfl[i].file; + memory->sfree(funcfl[i].frho); + memory->sfree(funcfl[i].rhor); + memory->sfree(funcfl[i].zr); + } + memory->sfree(funcfl); } - if (frho_0) interpolate_deallocate(); + if (setfl) { + for (int i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i]; + delete [] setfl->elements; + delete [] setfl->mass; + memory->destroy_2d_double_array(setfl->frho); + memory->destroy_2d_double_array(setfl->rhor); + memory->destroy_3d_double_array(setfl->z2r); + delete setfl; + } + + if (fs) { + for (int i = 0; i < fs->nelements; i++) delete [] fs->elements[i]; + delete [] fs->elements; + delete [] fs->mass; + memory->destroy_2d_double_array(fs->frho); + memory->destroy_3d_double_array(fs->rhor); + memory->destroy_3d_double_array(fs->z2r); + delete fs; + } + + memory->destroy_2d_double_array(frho); + memory->destroy_2d_double_array(rhor); + memory->destroy_2d_double_array(z2r); + + memory->destroy_3d_double_array(frho_spline); + memory->destroy_3d_double_array(rhor_spline); + memory->destroy_3d_double_array(z2r_spline); } /* ---------------------------------------------------------------------- */ @@ -94,6 +121,7 @@ void PairEAM::compute(int eflag, int vflag) int i,j,k,m,numneigh,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r,p,fforce,rhoip,rhojp,z2,z2p,recip,phi,phip,psip; + double *coeff; int *neighs; double **f; @@ -103,8 +131,8 @@ void PairEAM::compute(int eflag, int vflag) memory->sfree(rho); memory->sfree(fp); nmax = atom->nmax; - rho = (double *) memory->smalloc(nmax*sizeof(double),"eam:rho"); - fp = (double *) memory->smalloc(nmax*sizeof(double),"eam:fp"); + rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho"); + fp = (double *) memory->smalloc(nmax*sizeof(double),"pair:fp"); } eng_vdwl = 0.0; @@ -150,11 +178,12 @@ void PairEAM::compute(int eflag, int vflag) m = MIN(m,nr-1); p -= m; p = MIN(p,1.0); - rho[i] += ((rhor_3[jtype][m]*p + rhor_2[jtype][m])*p + - rhor_1[jtype][m])*p + rhor_0[jtype][m]; - if (newton_pair || j < nlocal) - rho[j] += ((rhor_3[itype][m]*p + rhor_2[itype][m])*p + - rhor_1[itype][m])*p + rhor_0[itype][m]; + coeff = rhor_spline[type2rhor[jtype][itype]][m]; + rho[i] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + if (newton_pair || j < nlocal) { + coeff = rhor_spline[type2rhor[itype][jtype]][m]; + rho[j] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + } } } } @@ -167,18 +196,14 @@ void PairEAM::compute(int eflag, int vflag) // phi = embedding energy at each atom for (i = 0; i < nlocal; i++) { - itype = type[i]; p = rho[i]*rdrho + 1.0; m = static_cast (p); m = MAX(1,MIN(m,nrho-1)); p -= m; p = MIN(p,1.0); - fp[i] = (frho_6[itype][m]*p + frho_5[itype][m])*p + frho_4[itype][m]; - if (eflag) { - phi = ((frho_3[itype][m]*p + frho_2[itype][m])*p + - frho_1[itype][m])*p + frho_0[itype][m]; - eng_vdwl += phi; - } + coeff = frho_spline[type2frho[type[i]]][m]; + fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2]; + if (eflag) eng_vdwl += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; } // communicate derivative of embedding function @@ -223,20 +248,20 @@ void PairEAM::compute(int eflag, int vflag) // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji) // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip - rhoip = (rhor_6[itype][m]*p + rhor_5[itype][m])*p + - rhor_4[itype][m]; - rhojp = (rhor_6[jtype][m]*p + rhor_5[jtype][m])*p + - rhor_4[jtype][m]; - z2 = ((z2r_3[itype][jtype][m]*p + z2r_2[itype][jtype][m])*p + - z2r_1[itype][jtype][m])*p + z2r_0[itype][jtype][m]; - z2p = (z2r_6[itype][jtype][m]*p + z2r_5[itype][jtype][m])*p + - z2r_4[itype][jtype][m]; - + coeff = rhor_spline[type2rhor[itype][jtype]][m]; + rhoip = (coeff[0]*p + coeff[1])*p + coeff[2]; + coeff = rhor_spline[type2rhor[jtype][itype]][m]; + rhojp = (coeff[0]*p + coeff[1])*p + coeff[2]; + coeff = z2r_spline[type2z2r[itype][jtype]][m]; + z2p = (coeff[0]*p + coeff[1])*p + coeff[2]; + z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; + recip = 1.0/r; phi = z2*recip; phip = z2p*recip - phi*recip; psip = fp[i]*rhojp + fp[j]*rhoip + phip; fforce = psip*recip; + f[i][0] -= delx*fforce; f[i][1] -= dely*fforce; f[i][2] -= delz*fforce; @@ -289,7 +314,13 @@ void PairEAM::allocate() setflag[i][j] = 0; cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); - tabindex = memory->create_2d_int_array(n+1,n+1,"pair:tabindex"); + + map = new int[n+1]; + for (int i = 1; i <= n; i++) map[n] = -1; + + type2frho = new int[n+1]; + type2rhor = memory->create_2d_int_array(n+1,n+1,"pair:type2rhor"); + type2z2r = memory->create_2d_int_array(n+1,n+1,"pair:type2z2r"); } /* ---------------------------------------------------------------------- @@ -303,7 +334,7 @@ void PairEAM::settings(int narg, char **arg) /* ---------------------------------------------------------------------- set coeffs for one or more type pairs - reading multiple funcfl files defines a funcfl alloy simulation + read DYNAMO funcfl file ------------------------------------------------------------------------- */ void PairEAM::coeff(int narg, char **arg) @@ -318,20 +349,33 @@ void PairEAM::coeff(int narg, char **arg) force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); - // read funcfl file only for i,i pairs - // only setflag i,i will be set - // set mass of each atom type + // read funcfl file if hasn't already been read + // store filename in Funcfl data struct - int itable; + int ifuncfl; + for (ifuncfl = 0; ifuncfl < nfuncfl; ifuncfl++) + if (strcmp(arg[2],funcfl[ifuncfl].file) == 0) break; + + if (ifuncfl == nfuncfl) { + nfuncfl++; + funcfl = (Funcfl *) + memory->srealloc(funcfl,nfuncfl*sizeof(Funcfl),"pair:funcfl"); + read_file(arg[2]); + int n = strlen(arg[2]) + 1; + funcfl[ifuncfl].file = new char[n]; + strcpy(funcfl[ifuncfl].file,arg[2]); + } + + // set setflag and map only for i,i type pairs + // set mass of atom type if i = j int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { if (i == j) { - itable = read_funcfl(arg[2]); - atom->set_mass(i,tables[itable].mass); - tabindex[i][i] = itable; setflag[i][i] = 1; + map[i] = ifuncfl; + atom->set_mass(i,funcfl[ifuncfl].mass); count++; } } @@ -346,23 +390,16 @@ void PairEAM::coeff(int narg, char **arg) double PairEAM::init_one(int i, int j) { - // only setflag I,I was set by coeff - // mixing will occur in init_style if both I,I and J,J were set + // single global cutoff = max of cut from all files read in + // for funcfl could be multiple files + // for setfl or fs, just one file - if (setflag[i][i] == 0 || setflag[j][j] == 0) - error->all("All EAM pair coeffs are not set"); - - // EAM has only one cutoff = max of all pairwise cutoffs - // determine max by checking table assigned to all type pairs - // only setflag[i][j] = 1 is relevant (if hybrid, some may not be set) - - cutmax = 0.0; - for (int ii = 1; ii <= atom->ntypes; ii++) { - for (int jj = ii; jj <= atom->ntypes; jj++) { - if (setflag[ii][jj] == 0) continue; - cutmax = MAX(cutmax,tables[tabindex[ii][jj]].cut); - } - } + if (funcfl) { + cutmax = 0.0; + for (int m = 0; m < nfuncfl; m++) + cutmax = MAX(cutmax,funcfl[m].cut); + } else if (setfl) cutmax = setfl->cut; + else if (fs) cutmax = fs->cut; return cutmax; } @@ -378,120 +415,98 @@ void PairEAM::init_style() comm->maxforward_pair = MAX(comm->maxforward_pair,1); comm->maxreverse_pair = MAX(comm->maxreverse_pair,1); - // convert read-in funcfl tables to multi-type setfl format and mix I,J - // interpolate final spline coeffs - - convert_funcfl(); - interpolate(); - + // convert read-in file(s) to arrays and spline them + + file2array(); + array2spline(); + cutforcesq = cutmax*cutmax; } /* ---------------------------------------------------------------------- - read potential values from a single element EAM file - read values into table and bcast values + read potential values from a DYNAMO single element funcfl file ------------------------------------------------------------------------- */ -int PairEAM::read_funcfl(char *file) +void PairEAM::read_file(char *filename) { - // check if same file has already been read - // if yes, return index of table entry - // if no, extend table list - - for (int i = 0; i < ntables; i++) - if (strcmp(file,tables->filename) == 0) return i; - - tables = (Table *) - memory->srealloc(tables,(ntables+1)*sizeof(Table),"pair:tables"); - - Table *tb = &tables[ntables]; - int n = strlen(file) + 1; - tb->filename = new char[n]; - strcpy(tb->filename,file); - tb->ith = tb->jth = 0; - - // open potential file + Funcfl *file = &funcfl[nfuncfl-1]; int me = comm->me; FILE *fp; char line[MAXLINE]; if (me == 0) { - fp = fopen(file,"r"); + fp = fopen(filename,"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open EAM potential file %s",file); + sprintf(str,"Cannot open EAM potential file %s",filename); error->one(str); } } - // read and broadcast header - int tmp; if (me == 0) { fgets(line,MAXLINE,fp); fgets(line,MAXLINE,fp); - sscanf(line,"%d %lg",&tmp,&tb->mass); + sscanf(line,"%d %lg",&tmp,&file->mass); fgets(line,MAXLINE,fp); sscanf(line,"%d %lg %d %lg %lg", - &tb->nrho,&tb->drho,&tb->nr,&tb->dr,&tb->cut); + &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut); } - MPI_Bcast(&tb->mass,1,MPI_DOUBLE,0,world); - MPI_Bcast(&tb->nrho,1,MPI_INT,0,world); - MPI_Bcast(&tb->drho,1,MPI_DOUBLE,0,world); - MPI_Bcast(&tb->nr,1,MPI_INT,0,world); - MPI_Bcast(&tb->dr,1,MPI_DOUBLE,0,world); - MPI_Bcast(&tb->cut,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->mass,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->nrho,1,MPI_INT,0,world); + MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->nr,1,MPI_INT,0,world); + MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world); - // allocate potential arrays and read/bcast them - // set z2r to NULL (setfl array) so it can be deallocated + file->frho = (double *) memory->smalloc((file->nrho+1)*sizeof(double), + "pair:frho"); + file->rhor = (double *) memory->smalloc((file->nr+1)*sizeof(double), + "pair:rhor"); + file->zr = (double *) memory->smalloc((file->nr+1)*sizeof(double), + "pair:zr"); - tb->frho = new double[tb->nrho+1]; - tb->zr = new double[tb->nr+1]; - tb->rhor = new double[tb->nr+1]; - tb->z2r = NULL; + if (me == 0) grab(fp,file->nrho,&file->frho[1]); + MPI_Bcast(&file->frho[1],file->nrho,MPI_DOUBLE,0,world); - if (me == 0) grab(fp,tb->nrho,&tb->frho[1]); - MPI_Bcast(&tb->frho[1],tb->nrho,MPI_DOUBLE,0,world); + if (me == 0) grab(fp,file->nr,&file->zr[1]); + MPI_Bcast(&file->zr[1],file->nr,MPI_DOUBLE,0,world); - if (me == 0) grab(fp,tb->nr,&tb->zr[1]); - MPI_Bcast(&tb->zr[1],tb->nr,MPI_DOUBLE,0,world); - - if (me == 0) grab(fp,tb->nr,&tb->rhor[1]); - MPI_Bcast(&tb->rhor[1],tb->nr,MPI_DOUBLE,0,world); - - // close the potential file + if (me == 0) grab(fp,file->nr,&file->rhor[1]); + MPI_Bcast(&file->rhor[1],file->nr,MPI_DOUBLE,0,world); if (me == 0) fclose(fp); - - ntables++; - return ntables-1; } /* ---------------------------------------------------------------------- - convert read-in funcfl potentials to multi-type setfl format + convert read-in funcfl potential(s) to standard array format + interpolate all file values to a single grid and cutoff ------------------------------------------------------------------------- */ -void PairEAM::convert_funcfl() +void PairEAM::file2array() { - int i,j,k,m; - + int i,j,k,m,n; int ntypes = atom->ntypes; - // determine max values for all i,i type pairs - // skip if setflag = 0 (if hybrid, some may not be set) + // determine max function params from all active funcfl files + // active means some element is pointing at it via map + int active; double rmax,rhomax; dr = drho = rmax = rhomax = 0.0; - for (int i = 1; i <= ntypes; i++) { - if (setflag[i][i] == 0) continue; - Table *tb = &tables[tabindex[i][i]]; - dr = MAX(dr,tb->dr); - drho = MAX(drho,tb->drho); - rmax = MAX(rmax,(tb->nr-1) * tb->dr); - rhomax = MAX(rhomax,(tb->nrho-1) * tb->drho); + for (int i = 0; i < nfuncfl; i++) { + active = 0; + for (j = 1; j <= ntypes; j++) + if (map[j] == i) active = 1; + if (active == 0) continue; + Funcfl *file = &funcfl[i]; + dr = MAX(dr,file->dr); + drho = MAX(drho,file->drho); + rmax = MAX(rmax,(file->nr-1) * file->dr); + rhomax = MAX(rhomax,(file->nrho-1) * file->drho); } // set nr,nrho from cutoff and spacings @@ -500,38 +515,29 @@ void PairEAM::convert_funcfl() nr = static_cast (rmax/dr + 0.5); nrho = static_cast (rhomax/drho + 0.5); - // allocate multi-type setfl arrays + // ------------------------------------------------------------------ + // setup frho arrays + // ------------------------------------------------------------------ - if (frho) { - memory->destroy_2d_double_array(frho); - memory->destroy_2d_double_array(rhor); - memory->destroy_2d_double_array(zrtmp); - memory->destroy_3d_double_array(z2r); - } + // allocate frho arrays + // nfrho = # of funcfl files + 1 for zero array + + nfrho = nfuncfl + 1; + memory->destroy_2d_double_array(frho); + frho = (double **) memory->create_2d_double_array(nfrho,nrho+1,"pair:frho"); - frho = (double **) - memory->create_2d_double_array(ntypes+1,nrho+1,"eam:frho"); - rhor = (double **) - memory->create_2d_double_array(ntypes+1,nr+1,"eam:rhor"); - zrtmp = (double **) - memory->create_2d_double_array(ntypes+1,nr+1,"eam:zrtmp"); - z2r = (double ***) - memory->create_3d_double_array(ntypes+1,ntypes+1,nr+1,"eam:frho"); - - // interpolate all potentials to a single grid and cutoff for all atom types - // frho,rhor are 1:ntypes, z2r is 1:ntypes,1:ntypes - // skip if setflag i,i or j,j = 0 (if hybrid, some may not be set) + // interpolate each file's frho to a single grid and cutoff double r,p,cof1,cof2,cof3,cof4; - for (i = 1; i <= ntypes; i++) { - if (setflag[i][i] == 0) continue; - Table *tb = &tables[tabindex[i][i]]; + n = 0; + for (i = 0; i < nfuncfl; i++) { + Funcfl *file = &funcfl[i]; for (m = 1; m <= nrho; m++) { r = (m-1)*drho; - p = r/tb->drho + 1.0; + p = r/file->drho + 1.0; k = static_cast (p); - k = MIN(k,tb->nrho-2); + k = MIN(k,file->nrho-2); k = MAX(k,2); p -= k; p = MIN(p,2.0); @@ -539,207 +545,203 @@ void PairEAM::convert_funcfl() cof2 = 0.5*(p*p-1.0)*(p-2.0); cof3 = -0.5*p*(p+1.0)*(p-2.0); cof4 = 0.166666667*p*(p*p-1.0); - frho[i][m] = cof1*tb->frho[k-1] + cof2*tb->frho[k] + - cof3*tb->frho[k+1] + cof4*tb->frho[k+2]; + frho[n][m] = cof1*file->frho[k-1] + cof2*file->frho[k] + + cof3*file->frho[k+1] + cof4*file->frho[k+2]; } + n++; } - for (i = 1; i <= ntypes; i++) { - if (setflag[i][i] == 0) continue; - Table *tb = &tables[tabindex[i][i]]; - for (m = 1; m <= nr; m++) { - r = (m-1)*dr; - p = r/tb->dr + 1.0; - k = static_cast (p); - k = MIN(k,tb->nr-2); - k = MAX(k,2); - p -= k; - p = MIN(p,2.0); - cof1 = -0.166666667*p*(p-1.0)*(p-2.0); - cof2 = 0.5*(p*p-1.0)*(p-2.0); - cof3 = -0.5*p*(p+1.0)*(p-2.0); - cof4 = 0.166666667*p*(p*p-1.0); - rhor[i][m] = cof1*tb->rhor[k-1] + cof2*tb->rhor[k] + - cof3*tb->rhor[k+1] + cof4*tb->rhor[k+2]; - zrtmp[i][m] = cof1*tb->zr[k-1] + cof2*tb->zr[k] + - cof3*tb->zr[k+1] + cof4*tb->zr[k+2]; - } - } + // add extra frho of zeroes for non-EAM types to point to (pair hybrid) + // this is necessary b/c fp is still computed for non-EAM atoms + + for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0; + + // type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to + // if atom type doesn't point to file (non-EAM atom in pair hybrid) + // then map it to last frho array of zeroes for (i = 1; i <= ntypes; i++) - for (j = i; j <= ntypes; j++) { - if (setflag[i][i] == 0 || setflag[j][j] == 0) continue; - for (m = 1; m <= nr; m++) - z2r[i][j][m] = 27.2*0.529 * zrtmp[i][m]*zrtmp[j][m]; + if (map[i] >= 0) type2frho[i] = map[i]; + else type2frho[i] = nfrho-1; + + // ------------------------------------------------------------------ + // setup rhor arrays + // ------------------------------------------------------------------ + + // allocate rhor arrays + // nrhor = # of funcfl files + + nrhor = nfuncfl; + memory->destroy_2d_double_array(rhor); + rhor = (double **) memory->create_2d_double_array(nrhor,nr+1,"pair:rhor"); + + // interpolate each file's rhor to a single grid and cutoff + + n = 0; + for (i = 0; i < nfuncfl; i++) { + Funcfl *file = &funcfl[i]; + for (m = 1; m <= nr; m++) { + r = (m-1)*dr; + p = r/file->dr + 1.0; + k = static_cast (p); + k = MIN(k,file->nr-2); + k = MAX(k,2); + p -= k; + p = MIN(p,2.0); + cof1 = -0.166666667*p*(p-1.0)*(p-2.0); + cof2 = 0.5*(p*p-1.0)*(p-2.0); + cof3 = -0.5*p*(p+1.0)*(p-2.0); + cof4 = 0.166666667*p*(p*p-1.0); + rhor[n][m] = cof1*file->rhor[k-1] + cof2*file->rhor[k] + + cof3*file->rhor[k+1] + cof4*file->rhor[k+2]; } + n++; + } + + // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to + // for funcfl files, I,J mapping only depends on I + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used + + for (i = 1; i <= ntypes; i++) + for (j = 1; j <= ntypes; j++) + type2rhor[i][j] = map[i]; + + // ------------------------------------------------------------------ + // setup z2r arrays + // ------------------------------------------------------------------ + + // allocate z2r arrays + // nz2r = N*(N+1)/2 where N = # of funcfl files + + nz2r = nfuncfl*(nfuncfl+1)/2; + memory->destroy_2d_double_array(z2r); + z2r = (double **) memory->create_2d_double_array(nz2r,nr+1,"pair:z2r"); + + // create a z2r array for each file against other files, only for I >= J + // interpolate zri and zrj to a single grid and cutoff + + double zri,zrj; + + n = 0; + for (i = 0; i < nfuncfl; i++) { + Funcfl *ifile = &funcfl[i]; + for (j = 0; j <= i; j++) { + Funcfl *jfile = &funcfl[j]; + + for (m = 1; m <= nr; m++) { + r = (m-1)*dr; + + p = r/ifile->dr + 1.0; + k = static_cast (p); + k = MIN(k,ifile->nr-2); + k = MAX(k,2); + p -= k; + p = MIN(p,2.0); + cof1 = -0.166666667*p*(p-1.0)*(p-2.0); + cof2 = 0.5*(p*p-1.0)*(p-2.0); + cof3 = -0.5*p*(p+1.0)*(p-2.0); + cof4 = 0.166666667*p*(p*p-1.0); + zri = cof1*ifile->zr[k-1] + cof2*ifile->zr[k] + + cof3*ifile->zr[k+1] + cof4*ifile->zr[k+2]; + + p = r/jfile->dr + 1.0; + k = static_cast (p); + k = MIN(k,jfile->nr-2); + k = MAX(k,2); + p -= k; + p = MIN(p,2.0); + cof1 = -0.166666667*p*(p-1.0)*(p-2.0); + cof2 = 0.5*(p*p-1.0)*(p-2.0); + cof3 = -0.5*p*(p+1.0)*(p-2.0); + cof4 = 0.166666667*p*(p*p-1.0); + zrj = cof1*jfile->zr[k-1] + cof2*jfile->zr[k] + + cof3*jfile->zr[k+1] + cof4*jfile->zr[k+2]; + + z2r[n][m] = 27.2*0.529 * zri*zrj; + } + n++; + } + } + + // type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to + // set of z2r arrays only fill lower triangular Nelement matrix + // value = n = sum over rows of lower-triangular matrix until reach irow,icol + // swap indices when irow < icol to stay lower triangular + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2z2r not used + + int irow,icol; + for (i = 1; i <= ntypes; i++) { + irow = map[i]; + if (irow == -1) continue; + for (j = 1; j <= ntypes; j++) { + icol = map[j]; + if (icol == -1) continue; + if (irow < icol) { + irow = map[j]; + icol = map[i]; + } + n = 0; + for (m = 0; m < irow; m++) n += m + 1; + n += icol; + type2z2r[i][j] = n; + } + } } -/* ---------------------------------------------------------------------- - interpolate EAM potentials -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ -void PairEAM::interpolate() +void PairEAM::array2spline() { - // free memory from previous interpolation - - if (frho_0) interpolate_deallocate(); - - // interpolation spacings - rdr = 1.0/dr; rdrho = 1.0/drho; - // allocate coeff arrays + memory->destroy_3d_double_array(frho_spline); + memory->destroy_3d_double_array(rhor_spline); + memory->destroy_3d_double_array(z2r_spline); - int n = atom->ntypes; + frho_spline = memory->create_3d_double_array(nfrho,nrho+1,7,"pair:frho"); + rhor_spline = memory->create_3d_double_array(nrhor,nr+1,7,"pair:rhor"); + z2r_spline = memory->create_3d_double_array(nz2r,nr+1,7,"pair:z2r"); - frho_0 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_0"); - frho_1 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_1"); - frho_2 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_2"); - frho_3 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_3"); - frho_4 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_4"); - frho_5 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_5"); - frho_6 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_6"); + for (int i = 0; i < nfrho; i++) + interpolate(nrho,drho,frho[i],frho_spline[i]); - rhor_0 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_0"); - rhor_1 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_1"); - rhor_2 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_2"); - rhor_3 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_3"); - rhor_4 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_4"); - rhor_5 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_5"); - rhor_6 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_6"); + for (int i = 0; i < nrhor; i++) + interpolate(nr,dr,rhor[i],rhor_spline[i]); - z2r_0 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_0"); - z2r_1 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_1"); - z2r_2 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_2"); - z2r_3 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_3"); - z2r_4 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_4"); - z2r_5 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_5"); - z2r_6 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_6"); + for (int i = 0; i < nz2r; i++) + interpolate(nr,dr,z2r[i],z2r_spline[i]); +} - // frho interpolation for 1:ntypes - // skip if setflag = 0 (if hybrid, some may not be set) - // if skip, set frho arrays to 0.0, since they will still be accessed - // for non-EAM atoms when compute() calculates embedding function +/* ---------------------------------------------------------------------- */ - int i,j,m; +void PairEAM::interpolate(int n, double delta, double *f, double **spline) +{ + for (int m = 1; m <= n; m++) spline[m][6] = f[m]; - for (i = 1; i <= atom->ntypes; i++) { - if (setflag[i][i] == 0) { - for (j = 1; j <= n; j++) - for (m = 1; m <= nrho; m++) - frho_0[j][m] = frho_1[j][m] = frho_2[j][m] = frho_3[j][m] = - frho_4[j][m] = frho_5[j][m] = frho_6[j][m] = 0.0; - continue; - } - - for (m = 1; m <= nrho; m++) frho_0[i][m] = frho[i][m]; - - frho_1[i][1] = frho_0[i][2]-frho_0[i][1]; - frho_1[i][2] = 0.5*(frho_0[i][3]-frho_0[i][1]); - frho_1[i][nrho-1] = 0.5*(frho_0[i][nrho]-frho_0[i][nrho-2]); - frho_1[i][nrho] = frho_0[i][nrho]-frho_0[i][nrho-1]; - - for (m = 3; m <= nrho-2; m++) - frho_1[i][m] = ((frho_0[i][m-2]-frho_0[i][m+2]) + - 8.0*(frho_0[i][m+1]-frho_0[i][m-1]))/12.0; - - for (m = 1; m <= nrho-1; m++) { - frho_2[i][m] = 3.*(frho_0[i][m+1]-frho_0[i][m]) - - 2.0*frho_1[i][m] - frho_1[i][m+1]; - frho_3[i][m] = frho_1[i][m] + frho_1[i][m+1] - - 2.0*(frho_0[i][m+1]-frho_0[i][m]); - } - - frho_2[i][nrho] = 0.0; - frho_3[i][nrho] = 0.0; - - for (m = 1; m <= nrho; m++) { - frho_4[i][m] = frho_1[i][m]/drho; - frho_5[i][m] = 2.0*frho_2[i][m]/drho; - frho_6[i][m] = 3.0*frho_3[i][m]/drho; - } + spline[1][5] = spline[2][6] - spline[1][6]; + spline[2][5] = 0.5 * (spline[3][6]-spline[1][6]); + spline[n-1][5] = 0.5 * (spline[n][6]-spline[n-2][6]); + spline[n][5] = spline[n][6] - spline[n-1][6]; + + for (int m = 3; m <= n-2; m++) + spline[m][5] = ((spline[m-2][6]-spline[m+2][6]) + + 8.0*(spline[m+1][6]-spline[m-1][6])) / 12.0; + + for (int m = 1; m <= n-1; m++) { + spline[m][4] = 3.0*(spline[m+1][6]-spline[m][6]) - + 2.0*spline[m][5] - spline[m+1][5]; + spline[m][3] = spline[m][5] + spline[m+1][5] - + 2.0*(spline[m+1][6]-spline[m][6]); } - - // rhor interpolation for 1:ntypes - // skip if setflag = 0 (if hybrid, some may not be set) - - for (i = 1; i <= atom->ntypes; i++) { - if (setflag[i][i] == 0) continue; - - for (m = 1; m <= nr; m++) rhor_0[i][m] = rhor[i][m]; - - rhor_1[i][1] = rhor_0[i][2]-rhor_0[i][1]; - rhor_1[i][2] = 0.5*(rhor_0[i][3]-rhor_0[i][1]); - rhor_1[i][nr-1] = 0.5*(rhor_0[i][nr]-rhor_0[i][nr-2]); - rhor_1[i][nr] = 0.0; - - for (m = 3; m <= nr-2; m++) - rhor_1[i][m] = ((rhor_0[i][m-2]-rhor_0[i][m+2]) + - 8.0*(rhor_0[i][m+1]-rhor_0[i][m-1]))/12.; - - for (m = 1; m <= nr-1; m++) { - rhor_2[i][m] = 3.0*(rhor_0[i][m+1]-rhor_0[i][m]) - - 2.0*rhor_1[i][m] - rhor_1[i][m+1]; - rhor_3[i][m] = rhor_1[i][m] + rhor_1[i][m+1] - - 2.0*(rhor_0[i][m+1]-rhor_0[i][m]); - } - - rhor_2[i][nr] = 0.0; - rhor_3[i][nr] = 0.0; - - for (m = 1; m <= nr; m++) { - rhor_4[i][m] = rhor_1[i][m]/dr; - rhor_5[i][m] = 2.0*rhor_2[i][m]/dr; - rhor_6[i][m] = 3.0*rhor_3[i][m]/dr; - } - } - - // z2r interpolation for 1:ntypes,1:ntypes - // skip if setflag i,i or j,j = 0 (if hybrid, some may not be set) - // set j,i coeffs = i,j coeffs - - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - if (setflag[i][i] == 0 || setflag[j][j] == 0) continue; - - for (m = 1; m <= nr; m++) z2r_0[i][j][m] = z2r[i][j][m]; - - z2r_1[i][j][1] = z2r_0[i][j][2]-z2r_0[i][j][1]; - z2r_1[i][j][2] = 0.5*(z2r_0[i][j][3]-z2r_0[i][j][1]); - z2r_1[i][j][nr-1] = 0.5*(z2r_0[i][j][nr]-z2r_0[i][j][nr-2]); - z2r_1[i][j][nr] = 0.0; - - for (m = 3; m <= nr-2; m++) - z2r_1[i][j][m] = ((z2r_0[i][j][m-2]-z2r_0[i][j][m+2]) + - 8.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m-1]))/12.; - - for (m = 1; m <= nr-1; m++) { - z2r_2[i][j][m] = 3.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m]) - - 2.0*z2r_1[i][j][m] - z2r_1[i][j][m+1]; - z2r_3[i][j][m] = z2r_1[i][j][m] + z2r_1[i][j][m+1] - - 2.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m]); - } - - z2r_2[i][j][nr] = 0.0; - z2r_3[i][j][nr] = 0.0; - - for (m = 1; m <= nr; m++) { - z2r_4[i][j][m] = z2r_1[i][j][m]/dr; - z2r_5[i][j][m] = 2.0*z2r_2[i][j][m]/dr; - z2r_6[i][j][m] = 3.0*z2r_3[i][j][m]/dr; - } - - for (m = 1; m <= nr; m++) { - z2r_0[j][i][m] = z2r_0[i][j][m]; - z2r_1[j][i][m] = z2r_1[i][j][m]; - z2r_2[j][i][m] = z2r_2[i][j][m]; - z2r_3[j][i][m] = z2r_3[i][j][m]; - z2r_4[j][i][m] = z2r_4[i][j][m]; - z2r_5[j][i][m] = z2r_5[i][j][m]; - z2r_6[j][i][m] = z2r_6[i][j][m]; - } - } + + spline[n][4] = 0.0; + spline[n][3] = 0.0; + + for (int m = 1; m <= n; m++) { + spline[m][2] = spline[m][5]/delta; + spline[m][1] = 2.0*spline[m][4]/delta; + spline[m][0] = 3.0*spline[m][3]/delta; } } @@ -763,64 +765,15 @@ void PairEAM::grab(FILE *fp, int n, double *list) } } -/* ---------------------------------------------------------------------- - skip n values from file fp - values can be several to a line - only called by proc 0 -------------------------------------------------------------------------- */ - -void PairEAM::skip(FILE *fp, int n) -{ - char line[MAXLINE]; - - int i = 0; - while (i < n) { - fgets(line,MAXLINE,fp); - strtok(line," \t\n\r\f"); - i++; - while (strtok(NULL," \t\n\r\f")) i++; - } -} - -/* ---------------------------------------------------------------------- - deallocate spline interpolation arrays -------------------------------------------------------------------------- */ - -void PairEAM::interpolate_deallocate() -{ - memory->destroy_2d_double_array(frho_0); - memory->destroy_2d_double_array(frho_1); - memory->destroy_2d_double_array(frho_2); - memory->destroy_2d_double_array(frho_3); - memory->destroy_2d_double_array(frho_4); - memory->destroy_2d_double_array(frho_5); - memory->destroy_2d_double_array(frho_6); - - memory->destroy_2d_double_array(rhor_0); - memory->destroy_2d_double_array(rhor_1); - memory->destroy_2d_double_array(rhor_2); - memory->destroy_2d_double_array(rhor_3); - memory->destroy_2d_double_array(rhor_4); - memory->destroy_2d_double_array(rhor_5); - memory->destroy_2d_double_array(rhor_6); - - memory->destroy_3d_double_array(z2r_0); - memory->destroy_3d_double_array(z2r_1); - memory->destroy_3d_double_array(z2r_2); - memory->destroy_3d_double_array(z2r_3); - memory->destroy_3d_double_array(z2r_4); - memory->destroy_3d_double_array(z2r_5); - memory->destroy_3d_double_array(z2r_6); -} - /* ---------------------------------------------------------------------- */ void PairEAM::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_lj, - int eflag, One &one) + double rsq, double factor_coul, double factor_lj, + int eflag, One &one) { - double r,p,rhoip,rhojp,z2,z2p,recip,phi,phip,psip; int m; + double r,p,rhoip,rhojp,z2,z2p,recip,phi,phip,psip; + double *coeff; r = sqrt(rsq); p = r*rdr + 1.0; @@ -828,15 +781,14 @@ void PairEAM::single(int i, int j, int itype, int jtype, m = MIN(m,nr-1); p -= m; p = MIN(p,1.0); - - rhoip = (rhor_6[itype][m]*p + rhor_5[itype][m])*p + - rhor_4[itype][m]; - rhojp = (rhor_6[jtype][m]*p + rhor_5[jtype][m])*p + - rhor_4[jtype][m]; - z2 = ((z2r_3[itype][jtype][m]*p + z2r_2[itype][jtype][m])*p + - z2r_1[itype][jtype][m])*p + z2r_0[itype][jtype][m]; - z2p = (z2r_6[itype][jtype][m]*p + z2r_5[itype][jtype][m])*p + - z2r_4[itype][jtype][m]; + + coeff = rhor_spline[type2rhor[itype][jtype]][m]; + rhoip = (coeff[0]*p + coeff[1])*p + coeff[2]; + coeff = rhor_spline[type2rhor[jtype][itype]][m]; + rhojp = (coeff[0]*p + coeff[1])*p + coeff[2]; + coeff = z2r_spline[type2z2r[itype][jtype]][m]; + z2p = (coeff[0]*p + coeff[1])*p + coeff[2]; + z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; recip = 1.0/r; phi = z2*recip; @@ -859,11 +811,10 @@ void PairEAM::single_embed(int i, int itype, double &fpi, int m = static_cast (p); m = MAX(1,MIN(m,nrho-1)); p -= m; - - fpi = (frho_6[itype][m]*p + frho_5[itype][m])*p + frho_4[itype][m]; - if (eflag) - phi = ((frho_3[itype][m]*p + frho_2[itype][m])*p + - frho_1[itype][m])*p + frho_0[itype][m]; + + double *coeff = frho_spline[type2frho[itype]][m]; + fpi = (coeff[0]*p + coeff[1])*p + coeff[2]; + if (eflag) phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; } /* ---------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_eam_alloy.cpp b/src/MANYBODY/pair_eam_alloy.cpp index ccd832c6ff..f43502a163 100644 --- a/src/MANYBODY/pair_eam_alloy.cpp +++ b/src/MANYBODY/pair_eam_alloy.cpp @@ -15,36 +15,26 @@ Contributing authors: Stephen Foiles (SNL), Murray Daw (SNL) ------------------------------------------------------------------------- */ -#include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_eam_alloy.h" #include "atom.h" -#include "force.h" #include "comm.h" #include "memory.h" -#include "neighbor.h" #include "error.h" -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - #define MAXLINE 1024 -/* ---------------------------------------------------------------------- */ - -PairEAMAlloy::PairEAMAlloy() -{ - one_coeff = 1; -} - /* ---------------------------------------------------------------------- set coeffs for one or more type pairs + read DYNAMO setfl file ------------------------------------------------------------------------- */ void PairEAMAlloy::coeff(int narg, char **arg) { + int i,j; + if (!allocated) allocate(); if (narg != 3 + atom->ntypes) @@ -55,32 +45,50 @@ void PairEAMAlloy::coeff(int narg, char **arg) if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->all("Incorrect args for pair coefficients"); - // read EAM setfl file, possibly multiple times - // first clear setflag since are doing this once for I,J = *,* - // read for all i,j pairs where ith,jth mapping is non-zero - // set setflag i,j for non-zero pairs - // set mass of atom type if i = j + // read EAM setfl file + + if (setfl) { + for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i]; + delete [] setfl->elements; + delete [] setfl->mass; + memory->destroy_2d_double_array(setfl->frho); + memory->destroy_2d_double_array(setfl->rhor); + memory->destroy_3d_double_array(setfl->z2r); + memory->sfree(setfl); + } + setfl = new Setfl(); + read_file(arg[2]); + + // read args that map atom types to elements in potential file + // map[i] = which element the Ith atom type is, -1 if NULL + + for (i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } + for (j = 0; j < setfl->nelements; j++) + if (strcmp(arg[i],setfl->elements[j]) == 0) break; + if (j < setfl->nelements) map[i-2] = j; + else error->all("No matching element in EAM potential file"); + } + + // clear setflag since coeff() called once with I,J = * * int n = atom->ntypes; - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) + for (i = 1; i <= n; i++) + for (j = i; j <= n; j++) setflag[i][j] = 0; - int itable,ith,jth; - int ilo,ihi,jlo,jhi; - ilo = jlo = 1; - ihi = jhi = n; + // set setflag i,j for type pairs where both are mapped to elements + // set mass of atom type if i = j int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - ith = atoi(arg[2+i]); - jth = atoi(arg[2+j]); - if (ith > 0 && jth > 0) { - itable = read_setfl(arg[2],ith,jth); - if (i == j) atom->set_mass(i,tables[itable].mass); - tabindex[i][j] = itable; + for (i = 1; i <= n; i++) { + for (j = i; j <= n; j++) { + if (map[i] >= 0 && map[j] >= 0) { setflag[i][j] = 1; + if (i == j) atom->set_mass(i,setfl->mass[map[i]]); count++; } } @@ -90,72 +98,12 @@ void PairEAMAlloy::coeff(int narg, char **arg) } /* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i + read a multi-element DYNAMO setfl file ------------------------------------------------------------------------- */ -double PairEAMAlloy::init_one(int i, int j) +void PairEAMAlloy::read_file(char *filename) { - if (setflag[i][j] == 0) - error->all("All EAM pair coeffs are not set"); - - // EAM has only one cutoff = max of all pairwise cutoffs - // determine max by checking table assigned to all type pairs - // only setflag[i][j] = 1 is relevant (if hybrid, some may not be set) - - cutmax = 0.0; - for (int ii = 1; ii <= atom->ntypes; ii++) { - for (int jj = ii; jj <= atom->ntypes; jj++) { - if (setflag[ii][jj] == 0) continue; - cutmax = MAX(cutmax,tables[tabindex[ii][jj]].cut); - } - } - - return cutmax; -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairEAMAlloy::init_style() -{ - // set communication sizes in comm class - - comm->maxforward_pair = MAX(comm->maxforward_pair,1); - comm->maxreverse_pair = MAX(comm->maxreverse_pair,1); - - // copy read-in-tables to multi-type setfl format - // interpolate final spline coeffs - - store_setfl(); - interpolate(); - - cutforcesq = cutmax*cutmax; -} - -/* ---------------------------------------------------------------------- - read ith,jth potential values from a multi-element alloy EAM file - read values into table and bcast values -------------------------------------------------------------------------- */ - -int PairEAMAlloy::read_setfl(char *file, int ith, int jth) -{ - // check if ith,jth portion of same file has already been read - // if yes, return index of table entry - // if no, extend table list - - for (int i = 0; i < ntables; i++) - if (ith == tables[i].ith && jth == tables[i].jth) return i; - - tables = (Table *) - memory->srealloc(tables,(ntables+1)*sizeof(Table),"pair:tables"); - - Table *tb = &tables[ntables]; - int n = strlen(file) + 1; - tb->filename = new char[n]; - strcpy(tb->filename,file); - tb->ith = ith; - tb->jth = jth; + Setfl *file = setfl; // open potential file @@ -164,309 +112,199 @@ int PairEAMAlloy::read_setfl(char *file, int ith, int jth) char line[MAXLINE]; if (me == 0) { - fp = fopen(file,"r"); + fp = fopen(filename,"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open EAM potential file %s",file); + sprintf(str,"Cannot open EAM potential file %s",filename); error->one(str); } } // read and broadcast header + // extract element names from nelements line - int ntypes; + int n; if (me == 0) { fgets(line,MAXLINE,fp); fgets(line,MAXLINE,fp); fgets(line,MAXLINE,fp); fgets(line,MAXLINE,fp); - sscanf(line,"%d",&ntypes); + n = strlen(line) + 1; + } + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + + sscanf(line,"%d",&file->nelements); + int nwords = atom->count_words(line); + if (nwords != file->nelements + 1) + error->all("Incorrect element names in EAM potential file"); + + char *words[file->nelements]; + nwords = 0; + char *first = strtok(line," \t\n\r\f"); + while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; + + file->elements = new char*[file->nelements]; + for (int i = 0; i < file->nelements; i++) { + n = strlen(words[i]) + 1; + file->elements[i] = new char[n]; + strcpy(file->elements[i],words[i]); + } + + if (me == 0) { fgets(line,MAXLINE,fp); sscanf(line,"%d %lg %d %lg %lg", - &tb->nrho,&tb->drho,&tb->nr,&tb->dr,&tb->cut); + &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut); } - MPI_Bcast(&ntypes,1,MPI_INT,0,world); - MPI_Bcast(&tb->nrho,1,MPI_INT,0,world); - MPI_Bcast(&tb->drho,1,MPI_DOUBLE,0,world); - MPI_Bcast(&tb->nr,1,MPI_INT,0,world); - MPI_Bcast(&tb->dr,1,MPI_DOUBLE,0,world); - MPI_Bcast(&tb->cut,1,MPI_DOUBLE,0,world); - - // check if ith,jth are consistent with ntypes - - if (ith > ntypes || jth > ntypes) - error->all("Requested atom types in EAM setfl file do not exist"); - - // allocate potential arrays and read/bcast them - // skip sections of file that don't correspond to ith,jth - // extract mass, frho, rhor for i,i from ith element section - // extract z2r for i,j from ith,jth array of z2r section - // note that ith can be < or > than jth - // set zr to NULL (funcl array) so it can be deallocated - - tb->frho = new double[tb->nrho+1]; - tb->rhor = new double[tb->nr+1]; - tb->z2r = new double[tb->nr+1]; - tb->zr = NULL; + MPI_Bcast(&file->nrho,1,MPI_INT,0,world); + MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->nr,1,MPI_INT,0,world); + MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world); + file->mass = new double[file->nelements]; + file->frho = memory->create_2d_double_array(file->nelements,file->nrho+1, + "pair:frho"); + file->rhor = memory->create_2d_double_array(file->nelements,file->nr+1, + "pair:rhor"); + file->z2r = memory->create_3d_double_array(file->nelements,file->nelements, + file->nr+1,"pair:z2r"); int i,j,tmp; - double mass; - - for (i = 1; i <= ntypes; i++) { + for (i = 0; i < file->nelements; i++) { if (me == 0) { fgets(line,MAXLINE,fp); - sscanf(line,"%d %lg",&tmp,&mass); + sscanf(line,"%d %lg",&tmp,&file->mass[i]); } - MPI_Bcast(&mass,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world); - if (i == ith && ith == jth) { - tb->mass = mass; - if (me == 0) grab(fp,tb->nrho,&tb->frho[1]); - MPI_Bcast(&tb->frho[1],tb->nrho,MPI_DOUBLE,0,world); - if (me == 0) grab(fp,tb->nr,&tb->rhor[1]); - MPI_Bcast(&tb->rhor[1],tb->nr,MPI_DOUBLE,0,world); - } else { - if (me == 0) skip(fp,tb->nrho); - if (me == 0) skip(fp,tb->nr); - } + if (me == 0) grab(fp,file->nrho,&file->frho[i][1]); + MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world); + if (me == 0) grab(fp,file->nr,&file->rhor[i][1]); + MPI_Bcast(&file->rhor[i][1],file->nr,MPI_DOUBLE,0,world); } - for (i = 1; i <= ntypes; i++) { - for (j = 1; j <= i; j++) { - if ((i == ith && j == jth) || (j == ith && i == jth)) { - if (me == 0) grab(fp,tb->nr,&tb->z2r[1]); - MPI_Bcast(&tb->z2r[1],tb->nr,MPI_DOUBLE,0,world); - } else if (me == 0) skip(fp,tb->nr); + for (i = 0; i < file->nelements; i++) + for (j = 0; j <= i; j++) { + if (me == 0) grab(fp,file->nr,&file->z2r[i][j][1]); + MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world); } - } // close the potential file if (me == 0) fclose(fp); - - ntables++; - return ntables-1; } /* ---------------------------------------------------------------------- - store read-in setfl values in multi-type setfl format + copy read-in setfl potential to standard array format ------------------------------------------------------------------------- */ -void PairEAMAlloy::store_setfl() +void PairEAMAlloy::file2array() { - int i,j,m; - + int i,j,m,n; int ntypes = atom->ntypes; + + // set function params directly from setfl file + + nrho = setfl->nrho; + nr = setfl->nr; + drho = setfl->drho; + dr = setfl->dr; + + // ------------------------------------------------------------------ + // setup frho arrays + // ------------------------------------------------------------------ + + // allocate frho arrays + // nfrho = # of setfl elements + 1 for zero array - // set nrho,nr,drho,dr from any i,i table entry since all the same + nfrho = setfl->nelements + 1; + memory->destroy_2d_double_array(frho); + frho = (double **) memory->create_2d_double_array(nfrho,nrho+1,"pair:frho"); + + // copy each element's frho to global frho + + for (i = 0; i < setfl->nelements; i++) + for (m = 1; m <= nrho; m++) frho[i][m] = setfl->frho[i][m]; + + // add extra frho of zeroes for non-EAM types to point to (pair hybrid) + // this is necessary b/c fp is still computed for non-EAM atoms + + for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0; + + // type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to + // if atom type doesn't point to element (non-EAM atom in pair hybrid) + // then map it to last frho array of zeroes for (i = 1; i <= ntypes; i++) - if (setflag[i][i]) break; + if (map[i] >= 0) type2frho[i] = map[i]; + else type2frho[i] = nfrho-1; - nrho = tables[tabindex[i][i]].nrho; - nr = tables[tabindex[i][i]].nr; - drho = tables[tabindex[i][i]].drho; - dr = tables[tabindex[i][i]].dr; + // ------------------------------------------------------------------ + // setup rhor arrays + // ------------------------------------------------------------------ - // allocate multi-type setfl arrays + // allocate rhor arrays + // nrhor = # of setfl elements - if (frho) { - memory->destroy_2d_double_array(frho); - memory->destroy_2d_double_array(rhor); - memory->destroy_2d_double_array(zrtmp); - memory->destroy_3d_double_array(z2r); - } + nrhor = setfl->nelements; + memory->destroy_2d_double_array(rhor); + rhor = (double **) memory->create_2d_double_array(nrhor,nr+1,"pair:rhor"); - frho = (double **) - memory->create_2d_double_array(ntypes+1,nrho+1,"eam:frho"); - rhor = (double **) - memory->create_2d_double_array(ntypes+1,nr+1,"eam:rhor"); - zrtmp = (double **) - memory->create_2d_double_array(ntypes+1,nr+1,"eam:zrtmp"); - z2r = (double ***) - memory->create_3d_double_array(ntypes+1,ntypes+1,nr+1,"eam:frho"); + // copy each element's rhor to global rhor - // copy from read-in tables to multi-type setfl arrays - // frho,rhor are 1:ntypes, z2r is 1:ntypes,1:ntypes - // skip if setflag i,j = 0 (if hybrid, some may not be set) + for (i = 0; i < setfl->nelements; i++) + for (m = 1; m <= nr; m++) rhor[i][m] = setfl->rhor[i][m]; + + // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to + // for setfl files, I,J mapping only depends on I + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used for (i = 1; i <= ntypes; i++) - for (j = i; j <= ntypes; j++) { - if (setflag[i][j] == 0) continue; - Table *tb = &tables[tabindex[i][j]]; - if (i == j) { - for (m = 1; m <= nrho; m++) frho[i][m] = tb->frho[m]; - for (m = 1; m <= nr; m++) rhor[i][m] = tb->rhor[m]; - } - for (m = 1; m <= nr; m++) z2r[i][j][m] = tb->z2r[m]; - } -} - -/* ---------------------------------------------------------------------- - interpolate EAM potentials -------------------------------------------------------------------------- */ - -void PairEAMAlloy::interpolate() -{ - // free memory from previous interpolation - - if (frho_0) interpolate_deallocate(); - - // interpolation spacings - - rdr = 1.0/dr; - rdrho = 1.0/drho; - - // allocate coeff arrays - - int n = atom->ntypes; - - frho_0 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_0"); - frho_1 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_1"); - frho_2 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_2"); - frho_3 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_3"); - frho_4 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_4"); - frho_5 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_5"); - frho_6 = memory->create_2d_double_array(n+1,nrho+1,"eam:frho_6"); - - rhor_0 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_0"); - rhor_1 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_1"); - rhor_2 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_2"); - rhor_3 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_3"); - rhor_4 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_4"); - rhor_5 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_5"); - rhor_6 = memory->create_2d_double_array(n+1,nr+1,"eam:rhor_6"); - - z2r_0 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_0"); - z2r_1 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_1"); - z2r_2 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_2"); - z2r_3 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_3"); - z2r_4 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_4"); - z2r_5 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_5"); - z2r_6 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam:z2r_6"); - - // frho interpolation for 1:ntypes - // skip if setflag = 0 (if hybrid, some may not be set) - // if skip, set frho arrays to 0.0, since they will still be accessed - // for non-EAM atoms when compute() calculates embedding function - - int i,j,m; - - for (i = 1; i <= atom->ntypes; i++) { - if (setflag[i][i] == 0) { - for (m = 1; m <= nrho; m++) - frho_0[i][m] = frho_1[i][m] = frho_2[i][m] = frho_3[i][m] = - frho_4[i][m] = frho_5[i][m] = frho_6[i][m] = 0.0; - continue; - } - - for (m = 1; m <= nrho; m++) frho_0[i][m] = frho[i][m]; - - frho_1[i][1] = frho_0[i][2]-frho_0[i][1]; - frho_1[i][2] = 0.5*(frho_0[i][3]-frho_0[i][1]); - frho_1[i][nrho-1] = 0.5*(frho_0[i][nrho]-frho_0[i][nrho-2]); - frho_1[i][nrho] = frho_0[i][nrho]-frho_0[i][nrho-1]; - - for (m = 3; m <= nrho-2; m++) - frho_1[i][m] = ((frho_0[i][m-2]-frho_0[i][m+2]) + - 8.0*(frho_0[i][m+1]-frho_0[i][m-1]))/12.0; - - for (m = 1; m <= nrho-1; m++) { - frho_2[i][m] = 3.*(frho_0[i][m+1]-frho_0[i][m]) - - 2.0*frho_1[i][m] - frho_1[i][m+1]; - frho_3[i][m] = frho_1[i][m] + frho_1[i][m+1] - - 2.0*(frho_0[i][m+1]-frho_0[i][m]); - } - - frho_2[i][nrho] = 0.0; - frho_3[i][nrho] = 0.0; - - for (m = 1; m <= nrho; m++) { - frho_4[i][m] = frho_1[i][m]/drho; - frho_5[i][m] = 2.0*frho_2[i][m]/drho; - frho_6[i][m] = 3.0*frho_3[i][m]/drho; - } - } - - // rhor interpolation for 1:ntypes - // skip if setflag = 0 (if hybrid, some may not be set) - - for (i = 1; i <= atom->ntypes; i++) { - if (setflag[i][i] == 0) continue; - - for (m = 1; m <= nr; m++) rhor_0[i][m] = rhor[i][m]; - - rhor_1[i][1] = rhor_0[i][2]-rhor_0[i][1]; - rhor_1[i][2] = 0.5*(rhor_0[i][3]-rhor_0[i][1]); - rhor_1[i][nr-1] = 0.5*(rhor_0[i][nr]-rhor_0[i][nr-2]); - rhor_1[i][nr] = 0.0; - - for (m = 3; m <= nr-2; m++) - rhor_1[i][m] = ((rhor_0[i][m-2]-rhor_0[i][m+2]) + - 8.0*(rhor_0[i][m+1]-rhor_0[i][m-1]))/12.; - - for (m = 1; m <= nr-1; m++) { - rhor_2[i][m] = 3.0*(rhor_0[i][m+1]-rhor_0[i][m]) - - 2.0*rhor_1[i][m] - rhor_1[i][m+1]; - rhor_3[i][m] = rhor_1[i][m] + rhor_1[i][m+1] - - 2.0*(rhor_0[i][m+1]-rhor_0[i][m]); - } - - rhor_2[i][nr] = 0.0; - rhor_3[i][nr] = 0.0; - - for (m = 1; m <= nr; m++) { - rhor_4[i][m] = rhor_1[i][m]/dr; - rhor_5[i][m] = 2.0*rhor_2[i][m]/dr; - rhor_6[i][m] = 3.0*rhor_3[i][m]/dr; - } - } - - // z2r interpolation for 1:ntypes,1:ntypes - // skip if setflag i,j = 0 (if hybrid, some may not be set) - // set j,i coeffs = i,j coeffs - - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - if (setflag[i][j] == 0) continue; - - for (m = 1; m <= nr; m++) z2r_0[i][j][m] = z2r[i][j][m]; - - z2r_1[i][j][1] = z2r_0[i][j][2]-z2r_0[i][j][1]; - z2r_1[i][j][2] = 0.5*(z2r_0[i][j][3]-z2r_0[i][j][1]); - z2r_1[i][j][nr-1] = 0.5*(z2r_0[i][j][nr]-z2r_0[i][j][nr-2]); - z2r_1[i][j][nr] = 0.0; - - for (m = 3; m <= nr-2; m++) - z2r_1[i][j][m] = ((z2r_0[i][j][m-2]-z2r_0[i][j][m+2]) + - 8.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m-1]))/12.; - - for (m = 1; m <= nr-1; m++) { - z2r_2[i][j][m] = 3.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m]) - - 2.0*z2r_1[i][j][m] - z2r_1[i][j][m+1]; - z2r_3[i][j][m] = z2r_1[i][j][m] + z2r_1[i][j][m+1] - - 2.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m]); - } - - z2r_2[i][j][nr] = 0.0; - z2r_3[i][j][nr] = 0.0; - - for (m = 1; m <= nr; m++) { - z2r_4[i][j][m] = z2r_1[i][j][m]/dr; - z2r_5[i][j][m] = 2.0*z2r_2[i][j][m]/dr; - z2r_6[i][j][m] = 3.0*z2r_3[i][j][m]/dr; - } - - for (m = 1; m <= nr; m++) { - z2r_0[j][i][m] = z2r_0[i][j][m]; - z2r_1[j][i][m] = z2r_1[i][j][m]; - z2r_2[j][i][m] = z2r_2[i][j][m]; - z2r_3[j][i][m] = z2r_3[i][j][m]; - z2r_4[j][i][m] = z2r_4[i][j][m]; - z2r_5[j][i][m] = z2r_5[i][j][m]; - z2r_6[j][i][m] = z2r_6[i][j][m]; + for (j = 1; j <= ntypes; j++) + type2rhor[i][j] = map[i]; + + // ------------------------------------------------------------------ + // setup z2r arrays + // ------------------------------------------------------------------ + + // allocate z2r arrays + // nz2r = N*(N+1)/2 where N = # of setfl elements + + nz2r = setfl->nelements * (setfl->nelements+1) / 2; + memory->destroy_2d_double_array(z2r); + z2r = (double **) memory->create_2d_double_array(nz2r,nr+1,"pair:z2r"); + + // copy each element pair z2r to global z2r, only for I >= J + + n = 0; + for (i = 0; i < setfl->nelements; i++) + for (j = 0; j <= i; j++) { + for (m = 1; m <= nr; m++) z2r[n][m] = setfl->z2r[i][j][m]; + n++; + } + + // type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to + // set of z2r arrays only fill lower triangular Nelement matrix + // value = n = sum over rows of lower-triangular matrix until reach irow,icol + // swap indices when irow < icol to stay lower triangular + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2z2r not used + + int irow,icol; + for (i = 1; i <= ntypes; i++) { + irow = map[i]; + if (irow == -1) continue; + for (j = 1; j <= ntypes; j++) { + icol = map[j]; + if (icol == -1) continue; + if (irow < icol) { + irow = map[j]; + icol = map[i]; } + n = 0; + for (m = 0; m < irow; m++) n += m + 1; + n += icol; + type2z2r[i][j] = n; } } } diff --git a/src/MANYBODY/pair_eam_alloy.h b/src/MANYBODY/pair_eam_alloy.h index 9296baff48..7231de3db6 100644 --- a/src/MANYBODY/pair_eam_alloy.h +++ b/src/MANYBODY/pair_eam_alloy.h @@ -18,15 +18,11 @@ class PairEAMAlloy : public PairEAM { public: - PairEAMAlloy(); void coeff(int, char **); - double init_one(int, int); - void init_style(); private: - int read_setfl(char *, int, int); - void store_setfl(); - void interpolate(); + void read_file(char *); + void file2array(); }; #endif diff --git a/src/MANYBODY/pair_eam_fs.cpp b/src/MANYBODY/pair_eam_fs.cpp index e1cf0cc78b..ddd61b6ae7 100644 --- a/src/MANYBODY/pair_eam_fs.cpp +++ b/src/MANYBODY/pair_eam_fs.cpp @@ -12,286 +12,83 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Tim Lau (MIT) + Contributing authors: Tim Lau (MIT) ------------------------------------------------------------------------- */ -#include "math.h" #include "stdio.h" #include "stdlib.h" #include "string.h" #include "pair_eam_fs.h" #include "atom.h" -#include "force.h" -#include "update.h" #include "comm.h" #include "memory.h" -#include "neighbor.h" #include "error.h" -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - #define MAXLINE 1024 -/* ---------------------------------------------------------------------- */ - -PairEAMFS::PairEAMFS() -{ - one_coeff = 1; -} - -/* ---------------------------------------------------------------------- */ - -PairEAMFS::~PairEAMFS() { - - // deallocate array unique to derived class, parent will deallocate the rest - - if (frho) memory->destroy_3d_double_array(rhor_fs); - - // insure derived-class's deallocate is called - // set ptr to NULL to prevent parent class from calling it's deallocate - - if (frho_0) interpolate_deallocate(); - frho_0 = NULL; -} - -/* ---------------------------------------------------------------------- */ - -void PairEAMFS::compute(int eflag, int vflag) -{ - int i,j,k,m,numneigh,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz; - double rsq,r,p,fforce,rhoip,rhojp,z2,z2p,recip,phi,phip,psip; - int *neighs; - double **f; - - // grow energy array if necessary - - if (atom->nmax > nmax) { - memory->sfree(rho); - memory->sfree(fp); - nmax = atom->nmax; - rho = (double *) memory->smalloc(nmax*sizeof(double),"eam:rho"); - fp = (double *) memory->smalloc(nmax*sizeof(double),"eam:fp"); - } - - eng_vdwl = 0.0; - if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - - if (vflag == 2) f = update->f_pair; - else f = atom->f; - double **x = atom->x; - int *type = atom->type; - int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; - - // zero out density - - if (newton_pair) { - m = nlocal + atom->nghost; - for (i = 0; i < m; i++) rho[i] = 0.0; - } else for (i = 0; i < nlocal; i++) rho[i] = 0.0; - - // rho = density at each atom - // loop over neighbors of my atoms - // FS has type-specific rho functional - - for (i = 0; i < nlocal; i++) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; - - for (k = 0; k < numneigh; k++) { - j = neighs[k]; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if (rsq < cutforcesq) { - jtype = type[j]; - p = sqrt(rsq)*rdr + 1.0; - m = static_cast (p); - m = MIN(m,nr-1); - p -= m; - p = MIN(p,1.0); - rho[i] += - ((rhor_fs_3[jtype][itype][m]*p + rhor_fs_2[jtype][itype][m])*p + - rhor_fs_1[jtype][itype][m])*p + rhor_fs_0[jtype][itype][m]; - if (newton_pair || j < nlocal) - rho[j] += - ((rhor_fs_3[itype][jtype][m]*p + rhor_fs_2[itype][jtype][m])*p + - rhor_fs_1[itype][jtype][m])*p + rhor_fs_0[itype][jtype][m]; - } - } - } - - // communicate and sum densities - - if (newton_pair) comm->reverse_comm_pair(this); - - // fp = derivative of embedding energy at each atom - // phi = embedding energy at each atom - // FS is same as standard EAM - - for (i = 0; i < nlocal; i++) { - itype = type[i]; - p = rho[i]*rdrho + 1.0; - m = static_cast (p); - m = MAX(1,MIN(m,nrho-1)); - p -= m; - p = MIN(p,1.0); - fp[i] = (frho_6[itype][m]*p + frho_5[itype][m])*p + frho_4[itype][m]; - if (eflag) { - phi = ((frho_3[itype][m]*p + frho_2[itype][m])*p + - frho_1[itype][m])*p + frho_0[itype][m]; - eng_vdwl += phi; - } - } - - // communicate derivative of embedding function - - comm->comm_pair(this); - - // compute forces on each atom - // loop over neighbors of my atoms - // FS has type-specific rhoip and rhojp - - for (i = 0; i < nlocal; i++) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; - - for (k = 0; k < numneigh; k++) { - j = neighs[k]; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if (rsq < cutforcesq) { - jtype = type[j]; - r = sqrt(rsq); - p = r*rdr + 1.0; - m = static_cast (p); - m = MIN(m,nr-1); - p -= m; - p = MIN(p,1.0); - - // rhoip = derivative of (density at atom j due to atom i) - // rhojp = derivative of (density at atom i due to atom j) - // phi = pair potential energy - // phip = phi' - // z2 = phi * r - // z2p = (phi * r)' = (phi' r) + phi - // psip needs both fp[i] and fp[j] terms since r_ij appears in two - // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji) - // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip - - rhoip = - (rhor_fs_6[itype][jtype][m]*p + rhor_fs_5[itype][jtype][m])*p + - rhor_fs_4[itype][jtype][m]; - rhojp = - (rhor_fs_6[jtype][itype][m]*p + rhor_fs_5[jtype][itype][m])*p + - rhor_fs_4[jtype][itype][m]; - z2 = ((z2r_3[itype][jtype][m]*p + z2r_2[itype][jtype][m])*p + - z2r_1[itype][jtype][m])*p + z2r_0[itype][jtype][m]; - z2p = (z2r_6[itype][jtype][m]*p + z2r_5[itype][jtype][m])*p + - z2r_4[itype][jtype][m]; - - recip = 1.0/r; - phi = z2*recip; - phip = z2p*recip - phi*recip; - psip = fp[i]*rhojp + fp[j]*rhoip + phip; - fforce = psip*recip; - f[i][0] -= delx*fforce; - f[i][1] -= dely*fforce; - f[i][2] -= delz*fforce; - if (newton_pair || j < nlocal) { - f[j][0] += delx*fforce; - f[j][1] += dely*fforce; - f[j][2] += delz*fforce; - } - - if (eflag) { - if (newton_pair || j < nlocal) eng_vdwl += phi; - else eng_vdwl += 0.5*phi; - } - - if (vflag == 1) { - if (newton_pair || j < nlocal) { - virial[0] -= delx*delx*fforce; - virial[1] -= dely*dely*fforce; - virial[2] -= delz*delz*fforce; - virial[3] -= delx*dely*fforce; - virial[4] -= delx*delz*fforce; - virial[5] -= dely*delz*fforce; - } else { - virial[0] -= 0.5*delx*delx*fforce; - virial[1] -= 0.5*dely*dely*fforce; - virial[2] -= 0.5*delz*delz*fforce; - virial[3] -= 0.5*delx*dely*fforce; - virial[4] -= 0.5*delx*delz*fforce; - virial[5] -= 0.5*dely*delz*fforce; - } - } - } - } - } - if (vflag == 2) virial_compute(); -} - /* ---------------------------------------------------------------------- set coeffs for one or more type pairs + read EAM Finnis-Sinclair file ------------------------------------------------------------------------- */ void PairEAMFS::coeff(int narg, char **arg) { + int i,j; + if (!allocated) allocate(); if (narg != 3 + atom->ntypes) error->all("Incorrect args for pair coefficients"); // insure I,J args are * * - // can only have one file, so clear all setflag settings if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->all("Incorrect args for pair coefficients"); - // read Finnis/Sinclair modified setfl file - // first clear setflag since are doing this once for I,J = *,* - // read for all i,j pairs where ith,jth mapping is non-zero - // set setflag i,j for non-zero pairs - // set mass of atom type if i = j + // read EAM Finnis-Sinclair file + + if (fs) { + for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i]; + delete [] fs->elements; + delete [] fs->mass; + memory->destroy_2d_double_array(fs->frho); + memory->destroy_3d_double_array(fs->rhor); + memory->destroy_3d_double_array(fs->z2r); + memory->sfree(fs); + } + fs = new Fs(); + read_file(arg[2]); + + // read args that map atom types to elements in potential file + // map[i] = which element the Ith atom type is, -1 if NULL + + for (i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } + for (j = 0; j < fs->nelements; j++) + if (strcmp(arg[i],fs->elements[j]) == 0) break; + if (j < fs->nelements) map[i-2] = j; + else error->all("No matching element in EAM potential file"); + } + + // clear setflag since coeff() called once with I,J = * * int n = atom->ntypes; - for (int i = 1; i <= n; i++) - for (int j = 1; j <= n; j++) + for (i = 1; i <= n; i++) + for (j = i; j <= n; j++) setflag[i][j] = 0; - int itable,ith,jth; - int ilo,ihi,jlo,jhi; - ilo = jlo = 1; - ihi = jhi = n; + // set setflag i,j for type pairs where both are mapped to elements + // set mass of atom type if i = j int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = jlo; j <= jhi; j++) { - ith = atoi(arg[2+i]); - jth = atoi(arg[2+j]); - if (ith > 0 && jth > 0) { - itable = read_setfl(arg[2],ith,jth); - if (i == j) atom->set_mass(i,tables[itable].mass); - tabindex[i][j] = itable; + for (i = 1; i <= n; i++) { + for (j = i; j <= n; j++) { + if (map[i] >= 0 && map[j] >= 0) { setflag[i][j] = 1; + if (i == j) atom->set_mass(i,fs->mass[map[i]]); count++; } } @@ -301,73 +98,12 @@ void PairEAMFS::coeff(int narg, char **arg) } /* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i + read a multi-element DYNAMO setfl file ------------------------------------------------------------------------- */ -double PairEAMFS::init_one(int i, int j) +void PairEAMFS::read_file(char *filename) { - if (setflag[i][j] == 0) - error->all("All EAM pair coeffs are not set"); - - // EAM has only one cutoff = max of all pairwise cutoffs - // determine max by checking table assigned to all type pairs - // only setflag[i][j] = 1 is relevant (if hybrid, some may not be set) - - cutmax = 0.0; - for (int ii = 1; ii <= atom->ntypes; ii++) { - for (int jj = ii; jj <= atom->ntypes; jj++) { - if (setflag[ii][jj] == 0) continue; - cutmax = MAX(cutmax,tables[tabindex[ii][jj]].cut); - } - } - - return cutmax; -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairEAMFS::init_style() -{ - // set communication sizes in comm class - - comm->maxforward_pair = MAX(comm->maxforward_pair,1); - comm->maxreverse_pair = MAX(comm->maxreverse_pair,1); - - // copy read-in-tables to multi-type setfl format - // interpolate final spline coeffs - - store_setfl(); - interpolate(); - - cutforcesq = cutmax*cutmax; -} - -/* ---------------------------------------------------------------------- - read ith,jth potential values from a multi-element alloy EAM/FS file - read values into table and bcast values - this file has different format than standard EAM setfl file -------------------------------------------------------------------------- */ - -int PairEAMFS::read_setfl(char *file, int ith, int jth) -{ - // check if ith,jth portion of same file has already been read - // if yes, return index of table entry - // if no, extend table list - - for (int i = 0; i < ntables; i++) - if (ith == tables[i].ith && jth == tables[i].jth) return i; - - tables = (Table *) - memory->srealloc(tables,(ntables+1)*sizeof(Table),"pair:tables"); - - Table *tb = &tables[ntables]; - int n = strlen(file) + 1; - tb->filename = new char[n]; - strcpy(tb->filename,file); - tb->ith = ith; - tb->jth = jth; + Fs *file = fs; // open potential file @@ -376,396 +112,206 @@ int PairEAMFS::read_setfl(char *file, int ith, int jth) char line[MAXLINE]; if (me == 0) { - fp = fopen(file,"r"); + fp = fopen(filename,"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open EAM potential file %s",file); + sprintf(str,"Cannot open EAM potential file %s",filename); error->one(str); } } // read and broadcast header + // extract element names from nelements line - int ntypes; + int n; if (me == 0) { fgets(line,MAXLINE,fp); fgets(line,MAXLINE,fp); fgets(line,MAXLINE,fp); fgets(line,MAXLINE,fp); - sscanf(line,"%d",&ntypes); + n = strlen(line) + 1; + } + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + + sscanf(line,"%d",&file->nelements); + int nwords = atom->count_words(line); + if (nwords != file->nelements + 1) + error->all("Incorrect element names in EAM potential file"); + + char *words[file->nelements]; + nwords = 0; + char *first = strtok(line," \t\n\r\f"); + while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; + + file->elements = new char*[file->nelements]; + for (int i = 0; i < file->nelements; i++) { + n = strlen(words[i]) + 1; + file->elements[i] = new char[n]; + strcpy(file->elements[i],words[i]); + } + + if (me == 0) { fgets(line,MAXLINE,fp); sscanf(line,"%d %lg %d %lg %lg", - &tb->nrho,&tb->drho,&tb->nr,&tb->dr,&tb->cut); + &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut); } - MPI_Bcast(&ntypes,1,MPI_INT,0,world); - MPI_Bcast(&tb->nrho,1,MPI_INT,0,world); - MPI_Bcast(&tb->drho,1,MPI_DOUBLE,0,world); - MPI_Bcast(&tb->nr,1,MPI_INT,0,world); - MPI_Bcast(&tb->dr,1,MPI_DOUBLE,0,world); - MPI_Bcast(&tb->cut,1,MPI_DOUBLE,0,world); - - // check if ith,jth are consistent with ntypes - - if (ith > ntypes || jth > ntypes) - error->all("Requested atom types in EAM setfl file do not exist"); - - // allocate potential arrays and read/bcast them - // skip sections of file that don't correspond to ith,jth - // extract mass, frho for i,i from ith element section - // extract rhor for i,j from jth array of ith element section - // extract z2r for i,j from ith,jth array of z2r section - // note that ith can be < or > than jth - // set zr to NULL (funcl array) so it can be deallocated - - tb->frho = new double[tb->nrho+1]; - tb->rhor = new double[tb->nr+1]; - tb->z2r = new double[tb->nr+1]; - tb->zr = NULL; + MPI_Bcast(&file->nrho,1,MPI_INT,0,world); + MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->nr,1,MPI_INT,0,world); + MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world); + file->mass = new double[file->nelements]; + file->frho = memory->create_2d_double_array(file->nelements,file->nrho+1, + "pair:frho"); + file->rhor = memory->create_3d_double_array(file->nelements,file->nelements, + file->nr+1,"pair:rhor"); + file->z2r = memory->create_3d_double_array(file->nelements,file->nelements, + file->nr+1,"pair:z2r"); int i,j,tmp; - double mass; - - for (i = 1; i <= ntypes; i++) { + for (i = 0; i < file->nelements; i++) { if (me == 0) { fgets(line,MAXLINE,fp); - sscanf(line,"%d %lg",&tmp,&mass); + sscanf(line,"%d %lg",&tmp,&file->mass[i]); } - MPI_Bcast(&mass,1,MPI_DOUBLE,0,world); + MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world); - if (i == ith && ith == jth) { - tb->mass = mass; - if (me == 0) grab(fp,tb->nrho,&tb->frho[1]); - MPI_Bcast(&tb->frho[1],tb->nrho,MPI_DOUBLE,0,world); + if (me == 0) grab(fp,file->nrho,&file->frho[i][1]); + MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world); - for (j = 1; j <= ntypes; j++) { - if (j == jth) { - if (me == 0) grab(fp,tb->nr,&tb->rhor[1]); - MPI_Bcast(&tb->rhor[1],tb->nr,MPI_DOUBLE,0,world); - } else if (me == 0) skip(fp,tb->nr); - } - - } else if (i == ith && ith != jth) { - if (me == 0) skip(fp,tb->nrho); - for (j = 1; j <= ntypes; j++) { - if (j == jth) { - if (me == 0) grab(fp,tb->nr,&tb->rhor[1]); - MPI_Bcast(&tb->rhor[1],tb->nr,MPI_DOUBLE,0,world); - } else if (me == 0) skip(fp,tb->nr); - } - - } else { - if (me == 0) skip(fp,tb->nrho); - if (me == 0) for (j = 1; j <= ntypes; j = j + 1) skip(fp,tb->nr); + for (j = 0; j < file->nelements; j++) { + if (me == 0) grab(fp,file->nr,&file->rhor[i][j][1]); + MPI_Bcast(&file->rhor[i][j][1],file->nr,MPI_DOUBLE,0,world); } } - for (i = 1; i <= ntypes; i++) { - for (j = 1; j <= i; j++) { - if ((i == ith && j == jth) || (j == ith && i == jth)) { - if (me == 0) grab(fp,tb->nr,&tb->z2r[1]); - MPI_Bcast(&tb->z2r[1],tb->nr,MPI_DOUBLE,0,world); - } else if (me == 0) skip(fp,tb->nr); + for (i = 0; i < file->nelements; i++) + for (j = 0; j <= i; j++) { + if (me == 0) grab(fp,file->nr,&file->z2r[i][j][1]); + MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world); } - } // close the potential file if (me == 0) fclose(fp); - - ntables++; - return ntables-1; } /* ---------------------------------------------------------------------- - store read-in setfl values in multi-type setfl format + copy read-in setfl potential to standard array format ------------------------------------------------------------------------- */ -void PairEAMFS::store_setfl() +void PairEAMFS::file2array() { - int i,j,m; - + int i,j,m,n; int ntypes = atom->ntypes; + + // set function params directly from fs file + + nrho = fs->nrho; + nr = fs->nr; + drho = fs->drho; + dr = fs->dr; + + // ------------------------------------------------------------------ + // setup frho arrays + // ------------------------------------------------------------------ + + // allocate frho arrays + // nfrho = # of fs elements + 1 for zero array - // set nrho,nr,drho,dr from any i,i table entry since all the same + nfrho = fs->nelements + 1; + memory->destroy_2d_double_array(frho); + frho = (double **) memory->create_2d_double_array(nfrho,nrho+1,"pair:frho"); + + // copy each element's frho to global frho + + for (i = 0; i < fs->nelements; i++) + for (m = 1; m <= nrho; m++) frho[i][m] = fs->frho[i][m]; + + // add extra frho of zeroes for non-EAM types to point to (pair hybrid) + // this is necessary b/c fp is still computed for non-EAM atoms + + for (m = 1; m <= nrho; m++) frho[nfrho-1][m] = 0.0; + + // type2frho[i] = which frho array (0 to nfrho-1) each atom type maps to + // if atom type doesn't point to element (non-EAM atom in pair hybrid) + // then map it to last frho array of zeroes for (i = 1; i <= ntypes; i++) - if (setflag[i][i]) break; + if (map[i] >= 0) type2frho[i] = map[i]; + else type2frho[i] = nfrho-1; - nrho = tables[tabindex[i][i]].nrho; - nr = tables[tabindex[i][i]].nr; - drho = tables[tabindex[i][i]].drho; - dr = tables[tabindex[i][i]].dr; + // ------------------------------------------------------------------ + // setup rhor arrays + // ------------------------------------------------------------------ - // allocate multi-type setfl arrays + // allocate rhor arrays + // nrhor = square of # of fs elements - if (frho) { - memory->destroy_2d_double_array(frho); - memory->destroy_3d_double_array(rhor_fs); - memory->destroy_2d_double_array(zrtmp); - memory->destroy_3d_double_array(z2r); - } + nrhor = fs->nelements * fs->nelements; + memory->destroy_2d_double_array(rhor); + rhor = (double **) memory->create_2d_double_array(nrhor,nr+1,"pair:rhor"); - frho = (double **) - memory->create_2d_double_array(ntypes+1,nrho+1,"eam/fs:frho"); - rhor_fs = (double ***) - memory->create_3d_double_array(ntypes+1,ntypes+1,nr+1,"eam/fs:rhor_fs"); - zrtmp = (double **) - memory->create_2d_double_array(ntypes+1,nr+1,"eam/fs:zrtmp"); - z2r = (double ***) - memory->create_3d_double_array(ntypes+1,ntypes+1,nr+1,"eam/fs:frho"); + // copy each element pair rhor to global rhor - // copy from read-in tables to multi-type setfl arrays - // frho,rhor are 1:ntypes, z2r is 1:ntypes,1:ntypes - // skip if setflag i,j = 0 (if hybrid, some may not be set) + n = 0; + for (i = 0; i < fs->nelements; i++) + for (j = 0; j < fs->nelements; j++) { + for (m = 1; m <= nr; m++) rhor[n][m] = fs->rhor[i][j][m]; + n++; + } + + // type2rhor[i][j] = which rhor array (0 to nrhor-1) each type pair maps to + // for fs files, there is a full NxN set of rhor arrays + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2rhor not used for (i = 1; i <= ntypes; i++) + for (j = 1; j <= ntypes; j++) + type2rhor[i][j] = map[i] * fs->nelements + map[j]; + + // ------------------------------------------------------------------ + // setup z2r arrays + // ------------------------------------------------------------------ + + // allocate z2r arrays + // nz2r = N*(N+1)/2 where N = # of fs elements + + nz2r = fs->nelements * (fs->nelements+1) / 2; + memory->destroy_2d_double_array(z2r); + z2r = (double **) memory->create_2d_double_array(nz2r,nr+1,"pair:z2r"); + + // copy each element pair z2r to global z2r, only for I >= J + + n = 0; + for (i = 0; i < fs->nelements; i++) + for (j = 0; j <= i; j++) { + for (m = 1; m <= nr; m++) z2r[n][m] = fs->z2r[i][j][m]; + n++; + } + + // type2z2r[i][j] = which z2r array (0 to nz2r-1) each type pair maps to + // set of z2r arrays only fill lower triangular Nelement matrix + // value = n = sum over rows of lower-triangular matrix until reach irow,icol + // swap indices when irow < icol to stay lower triangular + // OK if map = -1 (non-EAM atom in pair hybrid) b/c type2z2r not used + + int irow,icol; + for (i = 1; i <= ntypes; i++) { + irow = map[i]; + if (irow == -1) continue; for (j = 1; j <= ntypes; j++) { - if (setflag[i][j] == 0) continue; - Table *tb = &tables[tabindex[i][j]]; - if (i == j) for (m = 1; m <= nrho; m++) frho[i][m] = tb->frho[m]; - for (m = 1; m <= nr; m++) { - rhor_fs[i][j][m] = tb->rhor[m]; - z2r[i][j][m] = tb->z2r[m]; - } - } -} - -/* ---------------------------------------------------------------------- - interpolate EAM potentials -------------------------------------------------------------------------- */ - -void PairEAMFS::interpolate() -{ - // free memory from previous interpolation - - if (frho_0) interpolate_deallocate(); - - // interpolation spacings - - rdr = 1.0/dr; - rdrho = 1.0/drho; - - // allocate coeff arrays - - int n = atom->ntypes; - - frho_0 = memory->create_2d_double_array(n+1,nrho+1,"eam/fs:frho_0"); - frho_1 = memory->create_2d_double_array(n+1,nrho+1,"eam/fs:frho_1"); - frho_2 = memory->create_2d_double_array(n+1,nrho+1,"eam/fs:frho_2"); - frho_3 = memory->create_2d_double_array(n+1,nrho+1,"eam/fs:frho_3"); - frho_4 = memory->create_2d_double_array(n+1,nrho+1,"eam/fs:frho_4"); - frho_5 = memory->create_2d_double_array(n+1,nrho+1,"eam/fs:frho_5"); - frho_6 = memory->create_2d_double_array(n+1,nrho+1,"eam/fs:frho_6"); - - rhor_fs_0 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:rhor_fs_0"); - rhor_fs_1 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:rhor_fs_1"); - rhor_fs_2 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:rhor_fs_2"); - rhor_fs_3 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:rhor_fs_3"); - rhor_fs_4 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:rhor_fs_4"); - rhor_fs_5 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:rhor_fs_5"); - rhor_fs_6 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:rhor_fs_6"); - - z2r_0 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:z2r_0"); - z2r_1 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:z2r_1"); - z2r_2 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:z2r_2"); - z2r_3 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:z2r_3"); - z2r_4 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:z2r_4"); - z2r_5 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:z2r_5"); - z2r_6 = memory->create_3d_double_array(n+1,n+1,nr+1,"eam/fs:z2r_6"); - - // frho interpolation for 1:ntypes - // skip if setflag = 0 (if hybrid, some may not be set) - // if skip, set frho arrays to 0.0, since they will still be accessed - // for non-EAM atoms when compute() calculates embedding function - - int i,j,m; - - for (i = 1; i <= atom->ntypes; i++) { - if (setflag[i][i] == 0) { - for (j = 1; j <= n; j++) - for (m = 1; m <= nrho; m++) - frho_0[j][m] = frho_1[j][m] = frho_2[j][m] = frho_3[j][m] = - frho_4[j][m] = frho_5[j][m] = frho_6[j][m] = 0.0; - continue; - } - - for (m = 1; m <= nrho; m++) frho_0[i][m] = frho[i][m]; - - frho_1[i][1] = frho_0[i][2]-frho_0[i][1]; - frho_1[i][2] = 0.5*(frho_0[i][3]-frho_0[i][1]); - frho_1[i][nrho-1] = 0.5*(frho_0[i][nrho]-frho_0[i][nrho-2]); - frho_1[i][nrho] = frho_0[i][nrho]-frho_0[i][nrho-1]; - - for (m = 3; m <= nrho-2; m++) - frho_1[i][m] = ((frho_0[i][m-2]-frho_0[i][m+2]) + - 8.0*(frho_0[i][m+1]-frho_0[i][m-1]))/12.0; - - for (m = 1; m <= nrho-1; m++) { - frho_2[i][m] = 3.*(frho_0[i][m+1]-frho_0[i][m]) - - 2.0*frho_1[i][m] - frho_1[i][m+1]; - frho_3[i][m] = frho_1[i][m] + frho_1[i][m+1] - - 2.0*(frho_0[i][m+1]-frho_0[i][m]); - } - - frho_2[i][nrho] = 0.0; - frho_3[i][nrho] = 0.0; - - for (m = 1; m <= nrho; m++) { - frho_4[i][m] = frho_1[i][m]/drho; - frho_5[i][m] = 2.0*frho_2[i][m]/drho; - frho_6[i][m] = 3.0*frho_3[i][m]/drho; - } - } - - // rhor interpolation for 1:ntypes - // skip if setflag = 0 (if hybrid, some may not be set) - - for (i = 1; i <= atom->ntypes; i++) { - for (j = 1; j <= atom->ntypes; j++) { - if (setflag[i][j] == 0) continue; - - for (m = 1; m <= nr; m++) rhor_fs_0[i][j][m] = rhor_fs[i][j][m]; - - rhor_fs_1[i][j][1] = rhor_fs_0[i][j][2]-rhor_fs_0[i][j][1]; - rhor_fs_1[i][j][2] = 0.5*(rhor_fs_0[i][j][3]-rhor_fs_0[i][j][1]); - rhor_fs_1[i][j][nr-1] = 0.5*(rhor_fs_0[i][j][nr]-rhor_fs_0[i][j][nr-2]); - rhor_fs_1[i][j][nr] = 0.0; - - for (m = 3; m <= nr-2; m++) - rhor_fs_1[i][j][m] = ((rhor_fs_0[i][j][m-2]-rhor_fs_0[i][j][m+2]) + - 8.0*(rhor_fs_0[i][j][m+1]-rhor_fs_0[i][j][m-1]))/12.; - - for (m = 1; m <= nr-1; m++) { - rhor_fs_2[i][j][m] = 3.0*(rhor_fs_0[i][j][m+1]-rhor_fs_0[i][j][m]) - - 2.0*rhor_fs_1[i][j][m] - rhor_fs_1[i][j][m+1]; - rhor_fs_3[i][j][m] = rhor_fs_1[i][j][m] + rhor_fs_1[i][j][m+1] - - 2.0*(rhor_fs_0[i][j][m+1]-rhor_fs_0[i][j][m]); - } - - rhor_fs_2[i][j][nr] = 0.0; - rhor_fs_3[i][j][nr] = 0.0; - - for (m = 1; m <= nr; m++) { - rhor_fs_4[i][j][m] = rhor_fs_1[i][j][m]/dr; - rhor_fs_5[i][j][m] = 2.0*rhor_fs_2[i][j][m]/dr; - rhor_fs_6[i][j][m] = 3.0*rhor_fs_3[i][j][m]/dr; - } - } - } - - // z2r interpolation for 1:ntypes,1:ntypes - // skip if setflag i,j = 0 (if hybrid, some may not be set) - // set j,i coeffs = i,j coeffs - - for (i = 1; i <= atom->ntypes; i++) { - for (j = i; j <= atom->ntypes; j++) { - if (setflag[i][j] == 0) continue; - - for (m = 1; m <= nr; m++) z2r_0[i][j][m] = z2r[i][j][m]; - - z2r_1[i][j][1] = z2r_0[i][j][2]-z2r_0[i][j][1]; - z2r_1[i][j][2] = 0.5*(z2r_0[i][j][3]-z2r_0[i][j][1]); - z2r_1[i][j][nr-1] = 0.5*(z2r_0[i][j][nr]-z2r_0[i][j][nr-2]); - z2r_1[i][j][nr] = 0.0; - - for (m = 3; m <= nr-2; m++) - z2r_1[i][j][m] = ((z2r_0[i][j][m-2]-z2r_0[i][j][m+2]) + - 8.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m-1]))/12.; - - for (m = 1; m <= nr-1; m++) { - z2r_2[i][j][m] = 3.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m]) - - 2.0*z2r_1[i][j][m] - z2r_1[i][j][m+1]; - z2r_3[i][j][m] = z2r_1[i][j][m] + z2r_1[i][j][m+1] - - 2.0*(z2r_0[i][j][m+1]-z2r_0[i][j][m]); - } - - z2r_2[i][j][nr] = 0.0; - z2r_3[i][j][nr] = 0.0; - - for (m = 1; m <= nr; m++) { - z2r_4[i][j][m] = z2r_1[i][j][m]/dr; - z2r_5[i][j][m] = 2.0*z2r_2[i][j][m]/dr; - z2r_6[i][j][m] = 3.0*z2r_3[i][j][m]/dr; - } - - for (m = 1; m <= nr; m++) { - z2r_0[j][i][m] = z2r_0[i][j][m]; - z2r_1[j][i][m] = z2r_1[i][j][m]; - z2r_2[j][i][m] = z2r_2[i][j][m]; - z2r_3[j][i][m] = z2r_3[i][j][m]; - z2r_4[j][i][m] = z2r_4[i][j][m]; - z2r_5[j][i][m] = z2r_5[i][j][m]; - z2r_6[j][i][m] = z2r_6[i][j][m]; + icol = map[j]; + if (icol == -1) continue; + if (irow < icol) { + irow = map[j]; + icol = map[i]; } + n = 0; + for (m = 0; m < irow; m++) n += m + 1; + n += icol; + type2z2r[i][j] = n; } } } - -/* ---------------------------------------------------------------------- - deallocate spline interpolation arrays -------------------------------------------------------------------------- */ - -void PairEAMFS::interpolate_deallocate() -{ - memory->destroy_2d_double_array(frho_0); - memory->destroy_2d_double_array(frho_1); - memory->destroy_2d_double_array(frho_2); - memory->destroy_2d_double_array(frho_3); - memory->destroy_2d_double_array(frho_4); - memory->destroy_2d_double_array(frho_5); - memory->destroy_2d_double_array(frho_6); - - memory->destroy_3d_double_array(rhor_fs_0); - memory->destroy_3d_double_array(rhor_fs_1); - memory->destroy_3d_double_array(rhor_fs_2); - memory->destroy_3d_double_array(rhor_fs_3); - memory->destroy_3d_double_array(rhor_fs_4); - memory->destroy_3d_double_array(rhor_fs_5); - memory->destroy_3d_double_array(rhor_fs_6); - - memory->destroy_3d_double_array(z2r_0); - memory->destroy_3d_double_array(z2r_1); - memory->destroy_3d_double_array(z2r_2); - memory->destroy_3d_double_array(z2r_3); - memory->destroy_3d_double_array(z2r_4); - memory->destroy_3d_double_array(z2r_5); - memory->destroy_3d_double_array(z2r_6); -} - -/* ---------------------------------------------------------------------- */ - -void PairEAMFS::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_lj, - int eflag, One &one) -{ - double r,p,rhoip,rhojp,z2,z2p,recip,phi,phip,psip; - int m; - - r = sqrt(rsq); - p = r*rdr + 1.0; - m = static_cast (p); - m = MIN(m,nr-1); - p -= m; - p = MIN(p,1.0); - - rhoip = (rhor_fs_6[itype][jtype][m]*p + rhor_fs_5[itype][jtype][m])*p + - rhor_fs_4[itype][jtype][m]; - rhojp = (rhor_fs_6[jtype][itype][m]*p + rhor_fs_5[jtype][itype][m])*p + - rhor_fs_4[jtype][itype][m]; - z2 = ((z2r_3[itype][jtype][m]*p + z2r_2[itype][jtype][m])*p + - z2r_1[itype][jtype][m])*p + z2r_0[itype][jtype][m]; - z2p = (z2r_6[itype][jtype][m]*p + z2r_5[itype][jtype][m])*p + - z2r_4[itype][jtype][m]; - - recip = 1.0/r; - phi = z2*recip; - phip = z2p*recip - phi*recip; - psip = fp[i]*rhojp + fp[j]*rhoip + phip; - one.fforce = -psip*recip; - - if (eflag) { - one.eng_vdwl = phi; - one.eng_coul = 0.0; - } -} diff --git a/src/MANYBODY/pair_eam_fs.h b/src/MANYBODY/pair_eam_fs.h index 1787d1beab..baa977b2da 100644 --- a/src/MANYBODY/pair_eam_fs.h +++ b/src/MANYBODY/pair_eam_fs.h @@ -18,23 +18,11 @@ class PairEAMFS : public PairEAM { public: - PairEAMFS(); - ~PairEAMFS(); - void compute(int, int); void coeff(int, char **); - double init_one(int, int); - void init_style(); - void single(int, int, int, int, double, double, double, int, One &); private: - double ***rhor_fs; - double ***rhor_fs_0,***rhor_fs_1,***rhor_fs_2,***rhor_fs_3; - double ***rhor_fs_4,***rhor_fs_5,***rhor_fs_6; - - int read_setfl(char *, int, int); - void store_setfl(); - void interpolate(); - void interpolate_deallocate(); + void read_file(char *); + void file2array(); }; #endif