forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@103 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
parent
ed0d30f88d
commit
290625c0cf
|
@ -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<int> (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<int> (rmax/dr + 0.5);
|
||||
nrho = static_cast<int> (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<int> (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<int> (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<int> (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<int> (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<int> (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<int> (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];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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<int> (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<int> (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<int> (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<int> (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;
|
||||
}
|
||||
}
|
||||
|
|
|
@ -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
|
||||
|
|
Loading…
Reference in New Issue