diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 253edb343..3f998d474 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -1,194 +1,194 @@ #ifndef LMP_MEAM_H #define LMP_MEAM_H #include #include "memory.h" #define maxelt 5 namespace LAMMPS_NS { typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t; class MEAM { public: MEAM(Memory *mem); ~MEAM(); private: Memory *memory; // cutforce = force cutoff // cutforcesq = force cutoff squared double cutforce, cutforcesq; // Ec_meam = cohesive energy // re_meam = nearest-neighbor distance // Omega_meam = atomic volume // B_meam = bulk modulus // Z_meam = number of first neighbors for reference structure // ielt_meam = atomic number of element // A_meam = adjustable parameter // alpha_meam = sqrt(9*Omega*B/Ec) // rho0_meam = density scaling parameter // delta_meam = heat of formation for alloys // beta[0-3]_meam = electron density constants // t[0-3]_meam = coefficients on densities in Gamma computation // rho_ref_meam = background density for reference structure // ibar_meam(i) = selection parameter for Gamma function for elt i, // lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j) // neltypes = maximum number of element type defined // eltind = index number of pair (similar to Voigt notation; ij = ji) // phir = pair potential function array // phirar[1-6] = spline coeffs // attrac_meam = attraction parameter in Rose energy // repuls_meam = repulsion parameter in Rose energy // nn2_meam = 1 if second nearest neighbors are to be computed, else 0 // zbl_meam = 1 if zbl potential for small r to be use, else 0 // emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0 // bkgd_dyn = 1 if reference densities follows Dynamo, else 0 // Cmin_meam, Cmax_meam = min and max values in screening cutoff // rc_meam = cutoff distance for meam // delr_meam = cutoff region for meam // ebound_meam = factor giving maximum boundary of sceen fcn ellipse // augt1 = flag for whether t1 coefficient should be augmented // ialloy = flag for newer alloy formulation (as in dynamo code) // mix_ref_t = flag to recover "old" way of computing t in reference config // erose_form = selection parameter for form of E_rose function // gsmooth_factor = factor determining length of G smoothing region // vind[23]D = Voight notation index maps for 2 and 3D // v2D,v3D = array of factors to apply for Voight notation // nr,dr = pair function discretization parameters // nrar,rdrar = spline coeff array parameters double Ec_meam[maxelt + 1][maxelt + 1], re_meam[maxelt + 1][maxelt + 1]; double Omega_meam[maxelt + 1], Z_meam[maxelt + 1]; double A_meam[maxelt + 1], alpha_meam[maxelt + 1][maxelt + 1], rho0_meam[maxelt + 1]; double delta_meam[maxelt + 1][maxelt + 1]; double beta0_meam[maxelt + 1], beta1_meam[maxelt + 1]; double beta2_meam[maxelt + 1], beta3_meam[maxelt + 1]; double t0_meam[maxelt + 1], t1_meam[maxelt + 1]; double t2_meam[maxelt + 1], t3_meam[maxelt + 1]; double rho_ref_meam[maxelt + 1]; int ibar_meam[maxelt + 1], ielt_meam[maxelt + 1]; lattice_t lattce_meam[maxelt + 1][maxelt + 1]; int nn2_meam[maxelt + 1][maxelt + 1]; int zbl_meam[maxelt + 1][maxelt + 1]; int eltind[maxelt + 1][maxelt + 1]; int neltypes; double **phir; double **phirar, **phirar1, **phirar2, **phirar3, **phirar4, **phirar5, **phirar6; double attrac_meam[maxelt + 1][maxelt + 1], repuls_meam[maxelt + 1][maxelt + 1]; double Cmin_meam[maxelt + 1][maxelt + 1][maxelt + 1]; double Cmax_meam[maxelt + 1][maxelt + 1][maxelt + 1]; double rc_meam, delr_meam, ebound_meam[maxelt + 1][maxelt + 1]; int augt1, ialloy, mix_ref_t, erose_form; int emb_lin_neg, bkgd_dyn; double gsmooth_factor; int vind2D[3 + 1][3 + 1], vind3D[3 + 1][3 + 1][3 + 1]; int v2D[6 + 1], v3D[10 + 1]; int nr, nrar; double dr, rdrar; public: int nmax; double *rho,*rho0,*rho1,*rho2,*rho3,*frhop; double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b; double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave; int maxneigh; double *scrfcn,*dscrfcn,*fcpair; protected: void meam_checkindex(int, int, int, int*, int*); void G_gam(double, int, double, double*, int*); void dG_gam(double, int, double, double*, double*); void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh, int* firstneigh, int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap); void screen(int i, int j, double** x, double rijsq, double* sij, int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap); void calc_rho1(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh, double* scrfcn, double* fcpair); void dsij(int i, int j, int k, int jn, int numneigh, double rij2, double* dsij1, double* dsij2, int ntype, int* type, int* fmap, double** x, double* scrfcn, double* fcpair); void fcut(double, double*); void dfcut(double, double*, double*); void dCfunc(double, double, double, double*); void dCfunc2(double, double, double, double*, double*); void alloyparams(); void compute_pair_meam(); double phi_meam(double, int, int); void compute_reference_density(); void get_shpfcn(double*, lattice_t); void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double, double, double, double, int, int, lattice_t); void get_Zij(int*, lattice_t); void get_Zij2(int*, double*, double*, lattice_t, double, double); void get_sijk(double, int, int, int, double*); void get_densref(double, int, int, double*, double*, double*, double*, double*, double*, double*, double*); double zbl(double, int, int); double erose(double, double, double, double, double, double, int); void interpolate_meam(int); double compute_phi(double, int, int); public: - void meam_setup_global(int*, int*, double*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); - void meam_setup_param(int*, double*, int*, int*, int*); + void meam_setup_global(int, lattice_t*, double*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); + void meam_setup_param(int, double, int, int*, int*); void meam_setup_done(double*); void meam_dens_setup(int, int, int); void meam_dens_init(int* i, int* ntype, int* type, int* fmap, double** x, int* numneigh, int* firstneigh, int* numneigh_full, int* firstneigh_full, int fnoffset, int* errorflag); void meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, int* eflag_atom, double* eng_vdwl, double* eatom, int* ntype, int* type, int* fmap, int* errorflag); void meam_force(int* iptr, int* eflag_either, int* eflag_global, int* eflag_atom, int* vflag_atom, double* eng_vdwl, double* eatom, int* ntype, int* type, int* fmap, double** x, int* numneigh, int* firstneigh, int* numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom, int* errorflag); }; // Functions we need for compat #define iszero(f) (fabs(f) < 1e-20) #define setall2d(arr, v) \ { \ for (int __i = 1; __i <= maxelt; __i++) \ for (int __j = 1; __j <= maxelt; __j++) \ arr[__i][__j] = v; \ } #define setall3d(arr, v) \ { \ for (int __i = 1; __i <= maxelt; __i++) \ for (int __j = 1; __j <= maxelt; __j++) \ for (int __k = 1; __k <= maxelt; __k++) \ arr[__i][__j][__k] = v; \ } /* Fortran Array Semantics in C. - Stack-Allocated and global arrays are 1-based, declared as foo[N+1] and simply ignoring the first element - Multi-Dimensional MUST be declared in reverse, so that the order when accessing is the same as in Fortran - arrays that are passed externally via the meam_* functions must use the arr*v() functions below (or be used with 0-based indexing) */ // access data with same index as used in fortran (1-based) #define arr1v(ptr, i) ptr[i - 1] #define arr2v(ptr, i, j) ptr[j - 1][i - 1] }; #endif diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index 2df434d60..81dc61931 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -1,98 +1,86 @@ #include "meam.h" #include using namespace LAMMPS_NS; // // declaration in pair_meam.h: // // void meam_setup_global(int *, int *, double *, int *, double *, double *, // double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, // double *, double *, int *); // // call in pair_meam.cpp: // // meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, // alat,esub,asub,t0,t1,t2,t3,rozero,ibar); // // void -MEAM::meam_setup_global(int* nelt, int* lat, double* z, int* ielement, double* atwt, +MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha, double* b0, double* b1, double* b2, double* b3, double* alat, double* esub, double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero, int* ibar) { int i; double tmplat[maxelt]; - this->neltypes = *nelt; + this->neltypes = nelt; - for (i = 1; i <= *nelt; i++) { - if (arr1v(lat, i) == 0) - this->lattce_meam[i][i] = FCC; - else if (arr1v(lat, i) == 1) - this->lattce_meam[i][i] = BCC; - else if (arr1v(lat, i) == 2) - this->lattce_meam[i][i] = HCP; - else if (arr1v(lat, i) == 3) - this->lattce_meam[i][i] = DIM; - else if (arr1v(lat, i) == 4) - this->lattce_meam[i][i] = DIA; - else { - // unknown - } + for (i = 1; i <= nelt; i++) { + this->lattce_meam[i][i] = arr1v(lat, i); this->Z_meam[i] = arr1v(z, i); this->ielt_meam[i] = arr1v(ielement, i); this->alpha_meam[i][i] = arr1v(alpha, i); this->beta0_meam[i] = arr1v(b0, i); this->beta1_meam[i] = arr1v(b1, i); this->beta2_meam[i] = arr1v(b2, i); this->beta3_meam[i] = arr1v(b3, i); tmplat[i] = arr1v(alat, i); this->Ec_meam[i][i] = arr1v(esub, i); this->A_meam[i] = arr1v(asub, i); this->t0_meam[i] = arr1v(t0, i); this->t1_meam[i] = arr1v(t1, i); this->t2_meam[i] = arr1v(t2, i); this->t3_meam[i] = arr1v(t3, i); this->rho0_meam[i] = arr1v(rozero, i); this->ibar_meam[i] = arr1v(ibar, i); if (this->lattce_meam[i][i] == FCC) this->re_meam[i][i] = tmplat[i] / sqrt(2.0); else if (this->lattce_meam[i][i] == BCC) this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0; else if (this->lattce_meam[i][i] == HCP) this->re_meam[i][i] = tmplat[i]; else if (this->lattce_meam[i][i] == DIM) this->re_meam[i][i] = tmplat[i]; else if (this->lattce_meam[i][i] == DIA) this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0; else { // error } } // Set some defaults this->rc_meam = 4.0; this->delr_meam = 0.1; setall2d(this->attrac_meam, 0.0); setall2d(this->repuls_meam, 0.0); setall3d(this->Cmax_meam, 2.8); setall3d(this->Cmin_meam, 2.0); setall2d(this->ebound_meam, pow(2.8, 2) / (4.0 * (2.8 - 1.0))); setall2d(this->delta_meam, 0.0); setall2d(this->nn2_meam, 0); setall2d(this->zbl_meam, 1); this->gsmooth_factor = 99.0; this->augt1 = 1; this->ialloy = 0; this->mix_ref_t = 0; this->emb_lin_neg = 0; this->bkgd_dyn = 0; this->erose_form = 0; } diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index 0af0cbb45..ca12a9609 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -1,240 +1,220 @@ #include "meam.h" #include using namespace LAMMPS_NS; // // do a sanity check on index parameters void MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr) { //: idx[0..2] *ierr = 0; if (nidx < num) { *ierr = 2; return; } for (int i = 0; i < num; i++) { if ((idx[i] < 1) || (idx[i] > lim)) { *ierr = 3; return; } } } // // Declaration in pair_meam.h: // // void meam_setup_param(int *, double *, int *, int *, int *); // // in pair_meam.cpp // // meam_setup_param(&which,&value,&nindex,index,&errorflag); // // // // The "which" argument corresponds to the index of the "keyword" array // in pair_meam.cpp: // // 0 = Ec_meam // 1 = alpha_meam // 2 = rho0_meam // 3 = delta_meam // 4 = lattce_meam // 5 = attrac_meam // 6 = repuls_meam // 7 = nn2_meam // 8 = Cmin_meam // 9 = Cmax_meam // 10 = rc_meam // 11 = delr_meam // 12 = augt1 // 13 = gsmooth_factor // 14 = re_meam // 15 = ialloy // 16 = mixture_ref_t // 17 = erose_form // 18 = zbl_meam // 19 = emb_lin_neg // 20 = bkgd_dyn void -MEAM::meam_setup_param(int* which_p, double* value_p, int* nindex_p, - int* index /*index(3)*/, int* errorflag) +MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag) { //: index[0..2] - int i1, i2, val; + int i1, i2; + lattice_t vlat; *errorflag = 0; - int which = *which_p; - double value = *value_p; - int nindex = *nindex_p; switch (which) { // 0 = Ec_meam case 0: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->Ec_meam[index[0]][index[1]] = value; break; // 1 = alpha_meam case 1: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->alpha_meam[index[0]][index[1]] = value; break; // 2 = rho0_meam case 2: - meam_checkindex(1, maxelt, nindex, index, errorflag); + meam_checkindex(1, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->rho0_meam[index[0]] = value; break; // 3 = delta_meam case 3: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->delta_meam[index[0]][index[1]] = value; break; // 4 = lattce_meam case 4: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; - val = (int)value; - - if (val == 0) - this->lattce_meam[index[0]][index[1]] = FCC; - else if (val == 1) - this->lattce_meam[index[0]][index[1]] = BCC; - else if (val == 2) - this->lattce_meam[index[0]][index[1]] = HCP; - else if (val == 3) - this->lattce_meam[index[0]][index[1]] = DIM; - else if (val == 4) - this->lattce_meam[index[0]][index[1]] = DIA; - else if (val == 5) - this->lattce_meam[index[0]][index[1]] = B1; - else if (val == 6) - this->lattce_meam[index[0]][index[1]] = C11; - else if (val == 7) - this->lattce_meam[index[0]][index[1]] = L12; - else if (val == 8) - this->lattce_meam[index[0]][index[1]] = B2; + vlat = (lattice_t)value; + + this->lattce_meam[index[0]][index[1]] = vlat; break; // 5 = attrac_meam case 5: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->attrac_meam[index[0]][index[1]] = value; break; // 6 = repuls_meam case 6: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->repuls_meam[index[0]][index[1]] = value; break; // 7 = nn2_meam case 7: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; i1 = std::min(index[0], index[1]); i2 = std::max(index[0], index[1]); this->nn2_meam[i1][i2] = (int)value; break; // 8 = Cmin_meam case 8: - meam_checkindex(3, maxelt, nindex, index, errorflag); + meam_checkindex(3, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->Cmin_meam[index[0]][index[1]][index[2]] = value; break; // 9 = Cmax_meam case 9: - meam_checkindex(3, maxelt, nindex, index, errorflag); + meam_checkindex(3, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->Cmax_meam[index[0]][index[1]][index[2]] = value; break; // 10 = rc_meam case 10: this->rc_meam = value; break; // 11 = delr_meam case 11: this->delr_meam = value; break; // 12 = augt1 case 12: this->augt1 = (int)value; break; // 13 = gsmooth case 13: this->gsmooth_factor = value; break; // 14 = re_meam case 14: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; this->re_meam[index[0]][index[1]] = value; break; // 15 = ialloy case 15: this->ialloy = (int)value; break; // 16 = mixture_ref_t case 16: this->mix_ref_t = (int)value; break; // 17 = erose_form case 17: this->erose_form = (int)value; break; // 18 = zbl_meam case 18: - meam_checkindex(2, maxelt, nindex, index, errorflag); + meam_checkindex(2, neltypes, nindex, index, errorflag); if (*errorflag != 0) return; i1 = std::min(index[0], index[1]); i2 = std::max(index[0], index[1]); this->zbl_meam[i1][i2] = (int)value; break; // 19 = emb_lin_neg case 19: this->emb_lin_neg = (int)value; break; // 20 = bkgd_dyn case 20: this->bkgd_dyn = (int)value; break; default: *errorflag = 1; } } diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index 6eeb89050..2d9c0b737 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -1,847 +1,847 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Greg Wagner (SNL) ------------------------------------------------------------------------- */ #include #include #include #include #include "meam.h" #include "pair_meamc.h" #include "atom.h" #include "force.h" #include "comm.h" #include "memory.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define MAXLINE 1024 static const int nkeywords = 21; static const char *keywords[] = { "Ec","alpha","rho0","delta","lattce", "attrac","repuls","nn2","Cmin","Cmax","rc","delr", "augt1","gsmooth_factor","re","ialloy", "mixture_ref_t","erose_form","zbl", "emb_lin_neg","bkgd_dyn"}; /* ---------------------------------------------------------------------- */ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; restartinfo = 0; one_coeff = 1; manybody_flag = 1; allocated = 0; nelements = 0; elements = NULL; mass = NULL; meam_inst = new MEAM(memory); // set comm size needed by this Pair comm_forward = 38; comm_reverse = 30; } /* ---------------------------------------------------------------------- free all arrays check if allocated, since class can be destructed when incomplete ------------------------------------------------------------------------- */ PairMEAMC::~PairMEAMC() { delete meam_inst; for (int i = 0; i < nelements; i++) delete [] elements[i]; delete [] elements; delete [] mass; if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); delete [] map; delete [] fmap; } } /* ---------------------------------------------------------------------- */ void PairMEAMC::compute(int eflag, int vflag) { int i,j,ii,n,inum_half,errorflag; int *ilist_half,*numneigh_half,**firstneigh_half; int *numneigh_full,**firstneigh_full; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = eflag_global = vflag_global = eflag_atom = vflag_atom = 0; // neighbor list info inum_half = listhalf->inum; ilist_half = listhalf->ilist; numneigh_half = listhalf->numneigh; firstneigh_half = listhalf->firstneigh; numneigh_full = listfull->numneigh; firstneigh_full = listfull->firstneigh; // strip neighbor lists of any special bond flags before using with MEAM // necessary before doing neigh_f2c and neigh_c2f conversions each step if (neighbor->ago == 0) { neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full); } // check size of scrfcn based on half neighbor list int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; n = 0; for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]]; meam_inst->meam_dens_setup(atom->nmax, nall, n); double **x = atom->x; double **f = atom->f; int *type = atom->type; int ntype = atom->ntypes; // change neighbor list indices to Fortran indexing neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full); // 3 stages of MEAM calculation // loop over my atoms followed by communication int ifort; int offset = 0; errorflag = 0; for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; meam_inst->meam_dens_init(&ifort,&ntype,type,fmap,x, &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], offset, &errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); error->one(FLERR,str); } offset += numneigh_half[i]; } comm->reverse_comm_pair(this); meam_inst->meam_dens_final(&nlocal,&eflag_either,&eflag_global,&eflag_atom, &eng_vdwl,eatom,&ntype,type,fmap, &errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); error->one(FLERR,str); } comm->forward_comm_pair(this); offset = 0; // vptr is first value in vatom if it will be used by meam_force() // else vatom may not exist, so pass dummy ptr double **vptr; if (vflag_atom) vptr = vatom; else vptr = NULL; for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; meam_inst->meam_force(&ifort,&eflag_either,&eflag_global,&eflag_atom, &vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,x, &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], offset, f,vptr,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); error->one(FLERR,str); } offset += numneigh_half[i]; } // change neighbor list indices back to C indexing neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- */ void PairMEAMC::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); map = new int[n+1]; fmap = new int[n]; } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairMEAMC::settings(int narg, char **arg) { if (narg != 0) error->all(FLERR,"Illegal pair_style command"); } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairMEAMC::coeff(int narg, char **arg) { int i,j,m,n; if (!allocated) allocate(); if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients"); // insure I,J args are * * if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->all(FLERR,"Incorrect args for pair coefficients"); // read MEAM element names between 2 filenames // nelements = # of MEAM elements // elements = list of unique element names if (nelements) { for (i = 0; i < nelements; i++) delete [] elements[i]; delete [] elements; delete [] mass; } nelements = narg - 4 - atom->ntypes; if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients"); elements = new char*[nelements]; mass = new double[nelements]; for (i = 0; i < nelements; i++) { n = strlen(arg[i+3]) + 1; elements[i] = new char[n]; strcpy(elements[i],arg[i+3]); } // read MEAM library and parameter files // pass all parameters to MEAM package // tell MEAM package that setup is done read_files(arg[2],arg[2+nelements+1]); meam_inst->meam_setup_done(&cutmax); // read args that map atom types to MEAM elements // map[i] = which element the Ith atom type is, -1 if not mapped for (i = 4 + nelements; i < narg; i++) { m = i - (4+nelements) + 1; for (j = 0; j < nelements; j++) if (strcmp(arg[i],elements[j]) == 0) break; if (j < nelements) map[m] = j; else if (strcmp(arg[i],"NULL") == 0) map[m] = -1; else error->all(FLERR,"Incorrect args for pair coefficients"); } // clear setflag since coeff() called once with I,J = * * n = atom->ntypes; for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; // set setflag i,j for type pairs where both are mapped to elements // set mass for i,i in atom class int count = 0; for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) if (map[i] >= 0 && map[j] >= 0) { setflag[i][j] = 1; if (i == j) atom->set_mass(FLERR,i,mass[map[i]]); count++; } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairMEAMC::init_style() { if (force->newton_pair == 0) error->all(FLERR,"Pair style MEAM requires newton pair on"); // need full and half neighbor list int irequest_full = neighbor->request(this,instance_me); neighbor->requests[irequest_full]->id = 1; neighbor->requests[irequest_full]->half = 0; neighbor->requests[irequest_full]->full = 1; int irequest_half = neighbor->request(this,instance_me); neighbor->requests[irequest_half]->id = 2; // setup Fortran-style mapping array needed by MEAM package // fmap is indexed from 1:ntypes by Fortran and stores a Fortran index // if type I is not a MEAM atom, fmap stores a 0 for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1; } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use half or full ------------------------------------------------------------------------- */ void PairMEAMC::init_list(int id, NeighList *ptr) { if (id == 1) listfull = ptr; else if (id == 2) listhalf = ptr; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairMEAMC::init_one(int i, int j) { return cutmax; } /* ---------------------------------------------------------------------- */ void PairMEAMC::read_files(char *globalfile, char *userfile) { // open global meamf file on proc 0 FILE *fp; if (comm->me == 0) { fp = force->open_potential(globalfile); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open MEAM potential file %s",globalfile); error->one(FLERR,str); } } // allocate parameter arrays int params_per_line = 19; - int *lat = new int[nelements]; + lattice_t *lat = new lattice_t[nelements]; int *ielement = new int[nelements]; int *ibar = new int[nelements]; double *z = new double[nelements]; double *atwt = new double[nelements]; double *alpha = new double[nelements]; double *b0 = new double[nelements]; double *b1 = new double[nelements]; double *b2 = new double[nelements]; double *b3 = new double[nelements]; double *alat = new double[nelements]; double *esub = new double[nelements]; double *asub = new double[nelements]; double *t0 = new double[nelements]; double *t1 = new double[nelements]; double *t2 = new double[nelements]; double *t3 = new double[nelements]; double *rozero = new double[nelements]; bool *found = new bool[nelements]; for (int i = 0; i < nelements; i++) found[i] = false; // read each set of params from global MEAM file // one set of params can span multiple lines // store params if element name is in element list // if element name appears multiple times, only store 1st entry int i,n,nwords; char **words = new char*[params_per_line+1]; char line[MAXLINE],*ptr; int eof = 0; int nset = 0; while (1) { if (comm->me == 0) { ptr = fgets(line,MAXLINE,fp); if (ptr == NULL) { eof = 1; fclose(fp); } else n = strlen(line) + 1; } MPI_Bcast(&eof,1,MPI_INT,0,world); if (eof) break; MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); // strip comment, skip line if blank if ((ptr = strchr(line,'#'))) *ptr = '\0'; nwords = atom->count_words(line); if (nwords == 0) continue; // concatenate additional lines until have params_per_line words while (nwords < params_per_line) { n = strlen(line); if (comm->me == 0) { ptr = fgets(&line[n],MAXLINE-n,fp); if (ptr == NULL) { eof = 1; fclose(fp); } else n = strlen(line) + 1; } MPI_Bcast(&eof,1,MPI_INT,0,world); if (eof) break; MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); if ((ptr = strchr(line,'#'))) *ptr = '\0'; nwords = atom->count_words(line); } if (nwords != params_per_line) error->all(FLERR,"Incorrect format in MEAM potential file"); // words = ptrs to all words in line // strip single and double quotes from words nwords = 0; words[nwords++] = strtok(line,"' \t\n\r\f"); while ((words[nwords++] = strtok(NULL,"' \t\n\r\f"))) continue; // skip if element name isn't in element list for (i = 0; i < nelements; i++) if (strcmp(words[0],elements[i]) == 0) break; if (i >= nelements) continue; // skip if element already appeared if (found[i] == true) continue; found[i] = true; // map lat string to an integer if (strcmp(words[1],"fcc") == 0) lat[i] = FCC; else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC; else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP; else if (strcmp(words[1],"dim") == 0) lat[i] = DIM; else if (strcmp(words[1],"dia") == 0) lat[i] = DIA; else error->all(FLERR,"Unrecognized lattice type in MEAM file 1"); // store parameters z[i] = atof(words[2]); ielement[i] = atoi(words[3]); atwt[i] = atof(words[4]); alpha[i] = atof(words[5]); b0[i] = atof(words[6]); b1[i] = atof(words[7]); b2[i] = atof(words[8]); b3[i] = atof(words[9]); alat[i] = atof(words[10]); esub[i] = atof(words[11]); asub[i] = atof(words[12]); t0[i] = atof(words[13]); t1[i] = atof(words[14]); t2[i] = atof(words[15]); t3[i] = atof(words[16]); rozero[i] = atof(words[17]); ibar[i] = atoi(words[18]); nset++; } // error if didn't find all elements in file if (nset != nelements) error->all(FLERR,"Did not find all elements in MEAM library file"); // pass element parameters to MEAM package - meam_inst->meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, + meam_inst->meam_setup_global(nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, alat,esub,asub,t0,t1,t2,t3,rozero,ibar); // set element masses for (i = 0; i < nelements; i++) mass[i] = atwt[i]; // clean-up memory delete [] words; delete [] lat; delete [] ielement; delete [] ibar; delete [] z; delete [] atwt; delete [] alpha; delete [] b0; delete [] b1; delete [] b2; delete [] b3; delete [] alat; delete [] esub; delete [] asub; delete [] t0; delete [] t1; delete [] t2; delete [] t3; delete [] rozero; delete [] found; // done if user param file is NULL if (strcmp(userfile,"NULL") == 0) return; // open user param file on proc 0 if (comm->me == 0) { fp = force->open_potential(userfile); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open MEAM potential file %s",userfile); error->one(FLERR,str); } } // read settings // pass them one at a time to MEAM package // match strings to list of corresponding ints int which; double value; int nindex,index[3]; int maxparams = 6; char **params = new char*[maxparams]; int nparams; eof = 0; while (1) { if (comm->me == 0) { ptr = fgets(line,MAXLINE,fp); if (ptr == NULL) { eof = 1; fclose(fp); } else n = strlen(line) + 1; } MPI_Bcast(&eof,1,MPI_INT,0,world); if (eof) break; MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); // strip comment, skip line if blank if ((ptr = strchr(line,'#'))) *ptr = '\0'; nparams = atom->count_words(line); if (nparams == 0) continue; // words = ptrs to all words in line nparams = 0; params[nparams++] = strtok(line,"=(), '\t\n\r\f"); while (nparams < maxparams && (params[nparams++] = strtok(NULL,"=(), '\t\n\r\f"))) continue; nparams--; for (which = 0; which < nkeywords; which++) if (strcmp(params[0],keywords[which]) == 0) break; if (which == nkeywords) { char str[128]; sprintf(str,"Keyword %s in MEAM parameter file not recognized", params[0]); error->all(FLERR,str); } nindex = nparams - 2; for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]); // map lattce_meam value to an integer - if (which == 4) { + if (which == 4) { if (strcmp(params[nparams-1],"fcc") == 0) value = FCC; else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC; else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP; else if (strcmp(params[nparams-1],"dim") == 0) value = DIM; else if (strcmp(params[nparams-1],"dia") == 0) value = DIA; else if (strcmp(params[nparams-1],"b1") == 0) value = B1; else if (strcmp(params[nparams-1],"c11") == 0) value = C11; else if (strcmp(params[nparams-1],"l12") == 0) value = L12; else if (strcmp(params[nparams-1],"b2") == 0) value = B2; else error->all(FLERR,"Unrecognized lattice type in MEAM file 2"); } else value = atof(params[nparams-1]); // pass single setting to MEAM package int errorflag = 0; - meam_inst->meam_setup_param(&which,&value,&nindex,index,&errorflag); + meam_inst->meam_setup_param(which,value,nindex,index,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); error->all(FLERR,str); } } delete [] params; } /* ---------------------------------------------------------------------- */ int PairMEAMC::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,k,m; m = 0; for (i = 0; i < n; i++) { j = list[i]; buf[m++] = meam_inst->rho0[j]; buf[m++] = meam_inst->rho1[j]; buf[m++] = meam_inst->rho2[j]; buf[m++] = meam_inst->rho3[j]; buf[m++] = meam_inst->frhop[j]; buf[m++] = meam_inst->gamma[j]; buf[m++] = meam_inst->dgamma1[j]; buf[m++] = meam_inst->dgamma2[j]; buf[m++] = meam_inst->dgamma3[j]; buf[m++] = meam_inst->arho2b[j]; buf[m++] = meam_inst->arho1[j][0]; buf[m++] = meam_inst->arho1[j][1]; buf[m++] = meam_inst->arho1[j][2]; buf[m++] = meam_inst->arho2[j][0]; buf[m++] = meam_inst->arho2[j][1]; buf[m++] = meam_inst->arho2[j][2]; buf[m++] = meam_inst->arho2[j][3]; buf[m++] = meam_inst->arho2[j][4]; buf[m++] = meam_inst->arho2[j][5]; for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[j][k]; buf[m++] = meam_inst->arho3b[j][0]; buf[m++] = meam_inst->arho3b[j][1]; buf[m++] = meam_inst->arho3b[j][2]; buf[m++] = meam_inst->t_ave[j][0]; buf[m++] = meam_inst->t_ave[j][1]; buf[m++] = meam_inst->t_ave[j][2]; buf[m++] = meam_inst->tsq_ave[j][0]; buf[m++] = meam_inst->tsq_ave[j][1]; buf[m++] = meam_inst->tsq_ave[j][2]; } return m; } /* ---------------------------------------------------------------------- */ void PairMEAMC::unpack_forward_comm(int n, int first, double *buf) { int i,k,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { meam_inst->rho0[i] = buf[m++]; meam_inst->rho1[i] = buf[m++]; meam_inst->rho2[i] = buf[m++]; meam_inst->rho3[i] = buf[m++]; meam_inst->frhop[i] = buf[m++]; meam_inst->gamma[i] = buf[m++]; meam_inst->dgamma1[i] = buf[m++]; meam_inst->dgamma2[i] = buf[m++]; meam_inst->dgamma3[i] = buf[m++]; meam_inst->arho2b[i] = buf[m++]; meam_inst->arho1[i][0] = buf[m++]; meam_inst->arho1[i][1] = buf[m++]; meam_inst->arho1[i][2] = buf[m++]; meam_inst->arho2[i][0] = buf[m++]; meam_inst->arho2[i][1] = buf[m++]; meam_inst->arho2[i][2] = buf[m++]; meam_inst->arho2[i][3] = buf[m++]; meam_inst->arho2[i][4] = buf[m++]; meam_inst->arho2[i][5] = buf[m++]; for (k = 0; k < 10; k++) meam_inst->arho3[i][k] = buf[m++]; meam_inst->arho3b[i][0] = buf[m++]; meam_inst->arho3b[i][1] = buf[m++]; meam_inst->arho3b[i][2] = buf[m++]; meam_inst->t_ave[i][0] = buf[m++]; meam_inst->t_ave[i][1] = buf[m++]; meam_inst->t_ave[i][2] = buf[m++]; meam_inst->tsq_ave[i][0] = buf[m++]; meam_inst->tsq_ave[i][1] = buf[m++]; meam_inst->tsq_ave[i][2] = buf[m++]; } } /* ---------------------------------------------------------------------- */ int PairMEAMC::pack_reverse_comm(int n, int first, double *buf) { int i,k,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { buf[m++] = meam_inst->rho0[i]; buf[m++] = meam_inst->arho2b[i]; buf[m++] = meam_inst->arho1[i][0]; buf[m++] = meam_inst->arho1[i][1]; buf[m++] = meam_inst->arho1[i][2]; buf[m++] = meam_inst->arho2[i][0]; buf[m++] = meam_inst->arho2[i][1]; buf[m++] = meam_inst->arho2[i][2]; buf[m++] = meam_inst->arho2[i][3]; buf[m++] = meam_inst->arho2[i][4]; buf[m++] = meam_inst->arho2[i][5]; for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[i][k]; buf[m++] = meam_inst->arho3b[i][0]; buf[m++] = meam_inst->arho3b[i][1]; buf[m++] = meam_inst->arho3b[i][2]; buf[m++] = meam_inst->t_ave[i][0]; buf[m++] = meam_inst->t_ave[i][1]; buf[m++] = meam_inst->t_ave[i][2]; buf[m++] = meam_inst->tsq_ave[i][0]; buf[m++] = meam_inst->tsq_ave[i][1]; buf[m++] = meam_inst->tsq_ave[i][2]; } return m; } /* ---------------------------------------------------------------------- */ void PairMEAMC::unpack_reverse_comm(int n, int *list, double *buf) { int i,j,k,m; m = 0; for (i = 0; i < n; i++) { j = list[i]; meam_inst->rho0[j] += buf[m++]; meam_inst->arho2b[j] += buf[m++]; meam_inst->arho1[j][0] += buf[m++]; meam_inst->arho1[j][1] += buf[m++]; meam_inst->arho1[j][2] += buf[m++]; meam_inst->arho2[j][0] += buf[m++]; meam_inst->arho2[j][1] += buf[m++]; meam_inst->arho2[j][2] += buf[m++]; meam_inst->arho2[j][3] += buf[m++]; meam_inst->arho2[j][4] += buf[m++]; meam_inst->arho2[j][5] += buf[m++]; for (k = 0; k < 10; k++) meam_inst->arho3[j][k] += buf[m++]; meam_inst->arho3b[j][0] += buf[m++]; meam_inst->arho3b[j][1] += buf[m++]; meam_inst->arho3b[j][2] += buf[m++]; meam_inst->t_ave[j][0] += buf[m++]; meam_inst->t_ave[j][1] += buf[m++]; meam_inst->t_ave[j][2] += buf[m++]; meam_inst->tsq_ave[j][0] += buf[m++]; meam_inst->tsq_ave[j][1] += buf[m++]; meam_inst->tsq_ave[j][2] += buf[m++]; } } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double PairMEAMC::memory_usage() { double bytes = 11 * meam_inst->nmax * sizeof(double); bytes += (3 + 6 + 10 + 3 + 3 + 3) * meam_inst->nmax * sizeof(double); bytes += 3 * meam_inst->maxneigh * sizeof(double); return bytes; } /* ---------------------------------------------------------------------- strip special bond flags from neighbor list entries are not used with MEAM need to do here so Fortran lib doesn't see them done once per reneighbor so that neigh_f2c and neigh_c2f don't see them ------------------------------------------------------------------------- */ void PairMEAMC::neigh_strip(int inum, int *ilist, int *numneigh, int **firstneigh) { int i,j,ii,jnum; int *jlist; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK; } } /* ---------------------------------------------------------------------- toggle neighbor list indices between zero- and one-based values needed for access by MEAM Fortran library ------------------------------------------------------------------------- */ void PairMEAMC::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh) { int i,j,ii,jnum; int *jlist; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; for (j = 0; j < jnum; j++) jlist[j]--; } } void PairMEAMC::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh) { int i,j,ii,jnum; int *jlist; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; for (j = 0; j < jnum; j++) jlist[j]++; } }