diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 3f998d4746..789e21e6fe 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -62,21 +62,21 @@ class MEAM { // 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]; + double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt]; + double Omega_meam[maxelt], Z_meam[maxelt]; + double A_meam[maxelt], alpha_meam[maxelt][maxelt], + rho0_meam[maxelt]; + double delta_meam[maxelt][maxelt]; + double beta0_meam[maxelt], beta1_meam[maxelt]; + double beta2_meam[maxelt], beta3_meam[maxelt]; + double t0_meam[maxelt], t1_meam[maxelt]; + double t2_meam[maxelt], t3_meam[maxelt]; + double rho_ref_meam[maxelt]; + int ibar_meam[maxelt], ielt_meam[maxelt]; + lattice_t lattce_meam[maxelt][maxelt]; + int nn2_meam[maxelt][maxelt]; + int zbl_meam[maxelt][maxelt]; + int eltind[maxelt][maxelt]; int neltypes; double **phir; @@ -84,12 +84,12 @@ class MEAM { double **phirar, **phirar1, **phirar2, **phirar3, **phirar4, **phirar5, **phirar6; - double attrac_meam[maxelt + 1][maxelt + 1], - repuls_meam[maxelt + 1][maxelt + 1]; + double attrac_meam[maxelt][maxelt], + repuls_meam[maxelt][maxelt]; - 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]; + double Cmin_meam[maxelt][maxelt][maxelt]; + double Cmax_meam[maxelt][maxelt][maxelt]; + double rc_meam, delr_meam, ebound_meam[maxelt][maxelt]; int augt1, ialloy, mix_ref_t, erose_form; int emb_lin_neg, bkgd_dyn; double gsmooth_factor; @@ -165,15 +165,15 @@ public: #define setall2d(arr, v) \ { \ - for (int __i = 1; __i <= maxelt; __i++) \ - for (int __j = 1; __j <= maxelt; __j++) \ + for (int __i = 0; __i < maxelt; __i++) \ + for (int __j = 0; __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++) \ + for (int __i = 0; __i < maxelt; __i++) \ + for (int __j = 0; __j < maxelt; __j++) \ + for (int __k = 0; __k < maxelt; __k++) \ arr[__i][__j][__k] = v; \ } diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index a90e769730..347add205b 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -33,8 +33,8 @@ MEAM::meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, // Complete the calculation of density for (i = 1; i <= *nlocal; i++) { - elti = arr1v(fmap, arr1v(type, i)); - if (elti > 0) { + elti = fmap[arr1v(type, i)]; + if (elti >= 0) { arr1v(rho1, i) = 0.0; arr1v(rho2, i) = -1.0 / 3.0 * arr1v(arho2b, i) * arr1v(arho2b, i); arr1v(rho3, i) = 0.0; diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index bfa501a4d2..c265ed5dca 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -135,9 +135,9 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double rnorm, fc, dfc, drinv; drinv = 1.0 / this->delr_meam; - elti = arr1v(fmap, arr1v(type, i)); + elti = fmap[arr1v(type, i)]; - if (elti > 0) { + if (elti >= 0) { xitmp = arr2v(x, 1, i); yitmp = arr2v(x, 2, i); @@ -146,8 +146,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, for (jn = 1; jn <= numneigh; jn++) { j = arr1v(firstneigh, jn); - eltj = arr1v(fmap, arr1v(type, j)); - if (eltj > 0) { + eltj = fmap[arr1v(type, j)]; + if (eltj >= 0) { // First compute screening function itself, sij xjtmp = arr2v(x, 1, j); @@ -180,8 +180,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, k = arr1v(firstneigh_full, kn); if (k == j) continue; - eltk = arr1v(fmap, arr1v(type, k)); - if (eltk == 0) + eltk = fmap[arr1v(type, k)]; + if (eltk < 0) continue; xktmp = arr2v(x, 1, k); yktmp = arr2v(x, 2, k); @@ -252,7 +252,7 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, double ro0i, ro0j; double rhoa0i, rhoa1i, rhoa2i, rhoa3i, A1i, A2i, A3i; - elti = arr1v(fmap, arr1v(type, i)); + elti = fmap[arr1v(type, i)]; xtmp = arr2v(x, 1, i); ytmp = arr2v(x, 2, i); ztmp = arr2v(x, 3, i); @@ -265,7 +265,7 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, delij[3] = arr2v(x, 3, j) - ztmp; rij2 = delij[1] * delij[1] + delij[2] * delij[2] + delij[3] * delij[3]; if (rij2 < this->cutforcesq) { - eltj = arr1v(fmap, arr1v(type, j)); + eltj = fmap[arr1v(type, j)]; rij = sqrt(rij2); ai = rij / this->re_meam[elti][elti] - 1.0; aj = rij / this->re_meam[eltj][eltj] - 1.0; @@ -379,15 +379,15 @@ MEAM::screen(int i, int j, double** x, double rijsq, double* sij, double Cmax, Cmin, rbound; *sij = 1.0; - eltj = arr1v(fmap, arr1v(type, j)); - elti = arr1v(fmap, arr1v(type, j)); + eltj = fmap[arr1v(type, j)]; + elti = fmap[arr1v(type, j)]; // if rjksq > ebound*rijsq, atom k is definitely outside the ellipse rbound = this->ebound_meam[elti][eltj] * rijsq; for (nk = 1; nk <= numneigh_full; nk++) { k = arr1v(firstneigh_full, nk); - eltk = arr1v(fmap, arr1v(type, k)); + eltk = fmap[arr1v(type, k)]; if (k == j) continue; delxjk = arr2v(x, 1, k) - arr2v(x, 1, j); @@ -451,9 +451,9 @@ MEAM::dsij(int i, int j, int k, int jn, int numneigh, double rij2, double Cmax, Cmin, dCikj1, dCikj2; sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn); - elti = arr1v(fmap, arr1v(type, i)); - eltj = arr1v(fmap, arr1v(type, j)); - eltk = arr1v(fmap, arr1v(type, k)); + elti = fmap[arr1v(type, i)]; + eltj = fmap[arr1v(type, j)]; + eltk = fmap[arr1v(type, k)]; Cmax = this->Cmax_meam[elti][eltj][eltk]; Cmin = this->Cmin_meam[elti][eltj][eltk]; diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index 9bfeb571f5..5a5c018435 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -74,9 +74,9 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, // Compute forces atom i - elti = arr1v(fmap, arr1v(type, i)); + elti = fmap[arr1v(type, i)]; - if (elti > 0) { + if (elti >= 0) { xitmp = arr2v(x, 1, i); yitmp = arr2v(x, 2, i); zitmp = arr2v(x, 3, i); @@ -84,9 +84,9 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, // Treat each pair for (jn = 1; jn <= *numneigh; jn++) { j = arr1v(firstneigh, jn); - eltj = arr1v(fmap, arr1v(type, j)); + eltj = fmap[arr1v(type, j)]; - if (!iszero(arr1v(scrfcn, fnoffset + jn)) && eltj > 0) { + if (!iszero(arr1v(scrfcn, fnoffset + jn)) && eltj >= 0) { sij = arr1v(scrfcn, fnoffset + jn) * arr1v(fcpair, fnoffset + jn); delij[1] = arr2v(x, 1, j) - xitmp; @@ -98,7 +98,7 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, r = rij; // Compute phi and phip - ind = this->eltind[elti][eltj] - 1; //: TODO Remove -1 when reindexing eltind + ind = this->eltind[elti][eltj]; pp = rij * this->rdrar; kk = (int)pp; kk = std::min(kk, this->nrar - 2); @@ -514,8 +514,8 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, continue; //: cont jn loop for (kn = 1; kn <= *numneigh_full; kn++) { k = arr1v(firstneigh_full, kn); - eltk = arr1v(fmap, arr1v(type, k)); - if (k != j && eltk > 0) { + eltk = fmap[arr1v(type, k)]; + if (k != j && eltk >= 0) { dsij(i, j, k, jn, *numneigh, rij2, &dsij1, &dsij2, *ntype, type, fmap, x, &scrfcn[fnoffset], &fcpair[fnoffset]); if (!iszero(dsij1) || !iszero(dsij2)) { @@ -579,6 +579,6 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, // end of j loop } - // else if elti=0, this is not a meam atom + // else if elti<0, this is not a meam atom } } diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp index 1f701061a3..60481bf661 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -26,7 +26,7 @@ MEAM::meam_setup_done(double* cutmax) *cutmax = this->cutforce; // Augment t1 term - for (int i = 1; i <= maxelt; i++) + for (int i = 0; i < maxelt; i++) this->t1_meam[i] = this->t1_meam[i] + this->augt1 * 3.0 / 5.0 * this->t3_meam[i]; @@ -71,9 +71,9 @@ MEAM::meam_setup_done(double* cutmax) this->v3D[9] = 3; this->v3D[10] = 1; - nv2 = 1; - for (m = 1; m <= this->neltypes; m++) { - for (n = m; n <= this->neltypes; n++) { + nv2 = 0; + for (m = 0; m < this->neltypes; m++) { + for (n = m; n < this->neltypes; n++) { this->eltind[m][n] = nv2; this->eltind[n][m] = nv2; nv2 = nv2 + 1; @@ -99,8 +99,8 @@ MEAM::alloyparams(void) double eb; // Loop over pairs - for (i = 1; i <= this->neltypes; i++) { - for (j = 1; i <= this->neltypes; i++) { + for (i = 0; i < this->neltypes; i++) { + for (j = 0; i < this->neltypes; i++) { // Treat off-diagonal pairs // If i>j, set all equal to ij, set equal to the ineltypes; i++) { - for (j = 1; j <= i - 1; j++) { - for (k = 1; k <= this->neltypes; k++) { + for (i = 1; i < this->neltypes; i++) { + for (j = 0; j < i; j++) { + for (k = 0; k < this->neltypes; k++) { this->Cmin_meam[i][j][k] = this->Cmin_meam[j][i][k]; this->Cmax_meam[i][j][k] = this->Cmax_meam[j][i][k]; } @@ -157,9 +157,9 @@ MEAM::alloyparams(void) // atom k definitely lies outside the screening function ellipse (so // there is no need to calculate its effects). Here, compute it for all // triplets [i][j][k] so that ebound[i][j] is the maximized over k - for (i = 2; i <= this->neltypes; i++) { - for (j = 1; j <= this->neltypes; j++) { - for (k = 1; k <= this->neltypes; k++) { + for (i = 0; i < this->neltypes; i++) { + for (j = 0; j < this->neltypes; j++) { + for (k = 0; k < this->neltypes; k++) { eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) / (4.0 * (this->Cmax_meam[i][j][k] - 1.0)); this->ebound_meam[i][j] = std::max(this->ebound_meam[i][j], eb); @@ -240,9 +240,9 @@ MEAM::compute_pair_meam(void) "pair:phirar6"); // loop over pairs of element types - nv2 = 0; - for (a = 1; a <= this->neltypes; a++) { - for (b = a; b <= this->neltypes; b++) { + nv2 = 0; + for (a = 0; a < this->neltypes; a++) { + for (b = a; b < this->neltypes; b++) { // loop over r values and compute for (j = 0; j < this->nr; j++) { r = j * this->dr; @@ -593,7 +593,7 @@ MEAM::compute_reference_density(void) double rho0, rho0_2nn, arat, scrn; // loop over element types - for (a = 1; a <= this->neltypes; a++) { + for (a = 0; a < this->neltypes; a++) { Z = (int)this->Z_meam[a]; if (this->ibar_meam[a] <= 0) Gbar = 1.0; diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index 81dc619316..b4ef8e4eb8 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -30,25 +30,25 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub this->neltypes = nelt; - for (i = 1; i <= nelt; i++) { - this->lattce_meam[i][i] = arr1v(lat, i); + for (i = 0; i < nelt; i++) { + this->lattce_meam[i][i] = 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); + this->Z_meam[i] = z[i]; + this->ielt_meam[i] = ielement[i]; + this->alpha_meam[i][i] = alpha[i]; + this->beta0_meam[i] = b0[i]; + this->beta1_meam[i] = b1[i]; + this->beta2_meam[i] = b2[i]; + this->beta3_meam[i] = b3[i]; + tmplat[i] = alat[i]; + this->Ec_meam[i][i] = esub[i]; + this->A_meam[i] = asub[i]; + this->t0_meam[i] = t0[i]; + this->t1_meam[i] = t1[i]; + this->t2_meam[i] = t2[i]; + this->t3_meam[i] = t3[i]; + this->rho0_meam[i] = rozero[i]; + this->ibar_meam[i] = ibar[i]; if (this->lattce_meam[i][i] == FCC) this->re_meam[i][i] = tmplat[i] / sqrt(2.0); diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index ca12a96093..efcf26b104 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -15,7 +15,7 @@ MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr } for (int i = 0; i < num; i++) { - if ((idx[i] < 1) || (idx[i] > lim)) { + if ((idx[i] < 0) || (idx[i] >= lim)) { *ierr = 3; return; } diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index 2d9c0b737b..7ffeec360e 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -82,7 +82,6 @@ PairMEAMC::~PairMEAMC() memory->destroy(setflag); memory->destroy(cutsq); delete [] map; - delete [] fmap; } } @@ -145,7 +144,7 @@ void PairMEAMC::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_inst->meam_dens_init(&ifort,&ntype,type,fmap,x, + meam_inst->meam_dens_init(&ifort,&ntype,type,map,x, &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], offset, @@ -161,7 +160,7 @@ void PairMEAMC::compute(int eflag, int vflag) comm->reverse_comm_pair(this); meam_inst->meam_dens_final(&nlocal,&eflag_either,&eflag_global,&eflag_atom, - &eng_vdwl,eatom,&ntype,type,fmap, + &eng_vdwl,eatom,&ntype,type,map, &errorflag); if (errorflag) { char str[128]; @@ -184,7 +183,7 @@ void PairMEAMC::compute(int eflag, int vflag) 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, + &vflag_atom,&eng_vdwl,eatom,&ntype,type,map,x, &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], offset, @@ -216,7 +215,6 @@ void PairMEAMC::allocate() memory->create(cutsq,n+1,n+1,"pair:cutsq"); map = new int[n+1]; - fmap = new int[n]; } /* ---------------------------------------------------------------------- @@ -323,12 +321,6 @@ void PairMEAMC::init_style() 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; } /* ---------------------------------------------------------------------- @@ -597,7 +589,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) error->all(FLERR,str); } nindex = nparams - 2; - for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]); + for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]) - 1; // map lattce_meam value to an integer diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h index 0431c678a8..fc31cd55c6 100644 --- a/src/USER-MEAMC/pair_meamc.h +++ b/src/USER-MEAMC/pair_meamc.h @@ -49,8 +49,7 @@ class PairMEAMC : public Pair { char **elements; // names of unique elements double *mass; // mass of each element - int *map; // mapping from atom types to elements - int *fmap; // Fortran version of map array for MEAM lib + int *map; // mapping from atom types (1-indexed) to elements (1-indexed) void allocate(); void read_files(char *, char *);