Change indexing of all MEAM element arrays

- arrays in MEAM class
- eltind setting
- remove fmap translation

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 Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];
-  double Omega_meam[maxelt + 1], Z_meam[maxelt + 1];
+  double Omega_meam[maxelt], Z_meam[maxelt];
-  double A_meam[maxelt + 1], alpha_meam[maxelt + 1][maxelt + 1],
+  double A_meam[maxelt], alpha_meam[maxelt][maxelt],
-    rho0_meam[maxelt + 1];
+    rho0_meam[maxelt];
-  double delta_meam[maxelt + 1][maxelt + 1];
+  double delta_meam[maxelt][maxelt];
-  double beta0_meam[maxelt + 1], beta1_meam[maxelt + 1];
+  double beta0_meam[maxelt], beta1_meam[maxelt];
-  double beta2_meam[maxelt + 1], beta3_meam[maxelt + 1];
+  double beta2_meam[maxelt], beta3_meam[maxelt];
-  double t0_meam[maxelt + 1], t1_meam[maxelt + 1];
+  double t0_meam[maxelt], t1_meam[maxelt];
-  double t2_meam[maxelt + 1], t3_meam[maxelt + 1];
+  double t2_meam[maxelt], t3_meam[maxelt];
-  double rho_ref_meam[maxelt + 1];
+  double rho_ref_meam[maxelt];
-  int ibar_meam[maxelt + 1], ielt_meam[maxelt + 1];
+  int ibar_meam[maxelt], ielt_meam[maxelt];
-  lattice_t lattce_meam[maxelt + 1][maxelt + 1];
+  lattice_t lattce_meam[maxelt][maxelt];
-  int nn2_meam[maxelt + 1][maxelt + 1];
+  int nn2_meam[maxelt][maxelt];
-  int zbl_meam[maxelt + 1][maxelt + 1];
+  int zbl_meam[maxelt][maxelt];
-  int eltind[maxelt + 1][maxelt + 1];
+  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],
+  double attrac_meam[maxelt][maxelt],
-    repuls_meam[maxelt + 1][maxelt + 1];
+    repuls_meam[maxelt][maxelt];
 
-  double Cmin_meam[maxelt + 1][maxelt + 1][maxelt + 1];
+  double Cmin_meam[maxelt][maxelt][maxelt];
-  double Cmax_meam[maxelt + 1][maxelt + 1][maxelt + 1];
+  double Cmax_meam[maxelt][maxelt][maxelt];
-  double rc_meam, delr_meam, ebound_meam[maxelt + 1][maxelt + 1];
+  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 __i = 0; __i < maxelt; __i++)                                     \
-      for (int __j = 1; __j <= maxelt; __j++)                                  \
+      for (int __j = 0; __j < maxelt; __j++)                                   \
         arr[__i][__j] = v;                                                     \
   }
 #define setall3d(arr, v)                                                       \
   {                                                                            \
-    for (int __i = 1; __i <= maxelt; __i++)                                    \
+    for (int __i = 0; __i < maxelt; __i++)                                     \
-      for (int __j = 1; __j <= maxelt; __j++)                                  \
+      for (int __j = 0; __j < maxelt; __j++)                                   \
-        for (int __k = 1; __k <= maxelt; __k++)                                \
+        for (int __k = 0; __k < maxelt; __k++)                                 \
           arr[__i][__j][__k] = v;                                              \
   }
 

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));
+    elti = fmap[arr1v(type, i)];
-    if (elti > 0) {
+    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;

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));
+      eltj = fmap[arr1v(type, j)];
-      if (eltj > 0) {
+      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));
+          eltk = fmap[arr1v(type, k)];
-          if (eltk == 0)
+          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));
+  eltj = fmap[arr1v(type, j)];
-  elti = arr1v(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));
+  elti = fmap[arr1v(type, i)];
-  eltj = arr1v(fmap, arr1v(type, j));
+  eltj = fmap[arr1v(type, j)];
-  eltk = arr1v(fmap, arr1v(type, k));
+  eltk = fmap[arr1v(type, k)];
   Cmax = this->Cmax_meam[elti][eltj][eltk];
   Cmin = this->Cmin_meam[elti][eltj][eltk];
 

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));
+            eltk = fmap[arr1v(type, k)];
-            if (k != j && eltk > 0) {
+            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
   }
 }

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;
+  nv2 = 0;
-  for (m = 1; m <= this->neltypes; m++) {
+  for (m = 0; m < this->neltypes; m++) {
-    for (n = m; n <= this->neltypes; n++) {
+    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 (i = 0; i < this->neltypes; i++) {
-    for (j = 1; i <= this->neltypes; i++) {
+    for (j = 0; i < this->neltypes; i++) {
       // Treat off-diagonal pairs
       // If i>j, set all equal to i<j case (which has aready been set,
       // here or in the input file)
@@ -144,9 +144,9 @@ MEAM::alloyparams(void)
 
   // Cmin[i][k][j] is symmetric in i-j, but not k.  For all triplets
   // where i>j, set equal to the i<j element.  Likewise for Cmax.
-  for (i = 2; i <= this->neltypes; i++) {
+  for (i = 1; i < this->neltypes; i++) {
-    for (j = 1; j <= i - 1; j++) {
+    for (j = 0; j < i; j++) {
-      for (k = 1; k <= this->neltypes; k++) {
+      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 (i = 0; i < this->neltypes; i++) {
-    for (j = 1; j <= this->neltypes; j++) {
+    for (j = 0; j < this->neltypes; j++) {
-      for (k = 1; k <= this->neltypes; k++) {
+      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;
+  nv2 = 0;
-  for (a = 1; a <= this->neltypes; a++) {
+  for (a = 0; a < this->neltypes; a++) {
-    for (b = a; b <= this->neltypes; b++) {
+    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;

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++) {
+  for (i = 0; i < nelt; i++) {
-    this->lattce_meam[i][i] = arr1v(lat, i);
+    this->lattce_meam[i][i] = lat[i];
 
-    this->Z_meam[i] = arr1v(z, i);
+    this->Z_meam[i] = z[i];
-    this->ielt_meam[i] = arr1v(ielement, i);
+    this->ielt_meam[i] = ielement[i];
-    this->alpha_meam[i][i] = arr1v(alpha, i);
+    this->alpha_meam[i][i] = alpha[i];
-    this->beta0_meam[i] = arr1v(b0, i);
+    this->beta0_meam[i] = b0[i];
-    this->beta1_meam[i] = arr1v(b1, i);
+    this->beta1_meam[i] = b1[i];
-    this->beta2_meam[i] = arr1v(b2, i);
+    this->beta2_meam[i] = b2[i];
-    this->beta3_meam[i] = arr1v(b3, i);
+    this->beta3_meam[i] = b3[i];
-    tmplat[i] = arr1v(alat, i);
+    tmplat[i] = alat[i];
-    this->Ec_meam[i][i] = arr1v(esub, i);
+    this->Ec_meam[i][i] = esub[i];
-    this->A_meam[i] = arr1v(asub, i);
+    this->A_meam[i] = asub[i];
-    this->t0_meam[i] = arr1v(t0, i);
+    this->t0_meam[i] = t0[i];
-    this->t1_meam[i] = arr1v(t1, i);
+    this->t1_meam[i] = t1[i];
-    this->t2_meam[i] = arr1v(t2, i);
+    this->t2_meam[i] = t2[i];
-    this->t3_meam[i] = arr1v(t3, i);
+    this->t3_meam[i] = t3[i];
-    this->rho0_meam[i] = arr1v(rozero, i);
+    this->rho0_meam[i] = rozero[i];
-    this->ibar_meam[i] = arr1v(ibar, i);
+    this->ibar_meam[i] = ibar[i];
 
     if (this->lattce_meam[i][i] == FCC)
       this->re_meam[i][i] = tmplat[i] / sqrt(2.0);

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;
     }

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
 

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 *map;                     // mapping from atom types (1-indexed) to elements (1-indexed)
-  int *fmap;                    // Fortran version of map array for MEAM lib
 
   void allocate();
   void read_files(char *, char *);