From cd6d2c55d171e619b163f0d77fe3e919842daf52 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Thu, 22 Aug 2019 16:52:43 +0200
Subject: [PATCH 1/9] MEAM/C: helper function for x=1

---
 src/USER-MEAMC/meam.h             | 4 ++++
 src/USER-MEAMC/meam_dens_init.cpp | 2 +-
 src/USER-MEAMC/meam_force.cpp     | 4 ++--
 src/USER-MEAMC/pair_meamc.cpp     | 2 +-
 4 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
index 4e5b298ffb..9be4628022 100644
--- a/src/USER-MEAMC/meam.h
+++ b/src/USER-MEAMC/meam.h
@@ -234,6 +234,10 @@ static inline bool iszero(const double f) {
   return fabs(f) < 1e-20;
 }
 
+static inline bool isone(const double f) {
+  return fabs(f-1.0) < 1e-20;
+}
+
 // Helper functions
 
 static inline double fdiv_zero(const double n, const double d) {
diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp
index 155941422c..39e3ff9467 100644
--- a/src/USER-MEAMC/meam_dens_init.cpp
+++ b/src/USER-MEAMC/meam_dens_init.cpp
@@ -197,7 +197,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
     //     Now compute derivatives
     dscrfcn[jn] = 0.0;
     sfcij = sij * fcij;
-    if (!iszero(sfcij) && !iszero(sfcij - 1.0)) {
+    if (!iszero(sfcij) && !isone(sfcij)) {
       for (kn = 0; kn < numneigh_full; kn++) {
         k = firstneigh_full[kn];
         if (k == j) continue;
diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp
index 09aad90111..8c50296358 100644
--- a/src/USER-MEAMC/meam_force.cpp
+++ b/src/USER-MEAMC/meam_force.cpp
@@ -410,7 +410,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
 
         //     Now compute forces on other atoms k due to change in sij
 
-        if (iszero(sij) || iszero(sij - 1.0)) continue; //: cont jn loop
+        if (iszero(sij) || isone(sij)) continue; //: cont jn loop
 
         double dxik(0), dyik(0), dzik(0);
         double dxjk(0), dyjk(0), dzjk(0);
@@ -429,7 +429,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
 
             dsij1 = 0.0;
             dsij2 = 0.0;
-            if (!iszero(sij) && !iszero(sij - 1.0)) {
+            if (!iszero(sij) && !isone(sij)) {
               const double rbound = rij2 * this->ebound_meam[elti][eltj];
               delc = Cmax - Cmin;
               dxjk = x[k][0] - x[j][0];
diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp
index a9034a1af3..dc2e194dea 100644
--- a/src/USER-MEAMC/pair_meamc.cpp
+++ b/src/USER-MEAMC/pair_meamc.cpp
@@ -458,7 +458,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
     rozero[i] = atof(words[17]);
     ibar[i] = atoi(words[18]);
 
-    if (!iszero(t0[i]-1.0))
+    if (!isone(t0[i]))
       error->all(FLERR,"Unsupported parameter in MEAM potential file");
 
     nset++;

From 204529bcafcf76326edd2e25196b2cdfa9673ff4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Wed, 21 Aug 2019 18:38:45 +0200
Subject: [PATCH 2/9] MEAM/C: remove unused vars, refactoring for extensibility

---
 src/USER-MEAMC/meam.h                |  29 ++++-
 src/USER-MEAMC/meam_impl.cpp         |   2 +-
 src/USER-MEAMC/meam_setup_done.cpp   | 184 ++++++++++++++-------------
 src/USER-MEAMC/meam_setup_global.cpp |  28 ++--
 src/USER-MEAMC/pair_meamc.cpp        |  24 +---
 5 files changed, 151 insertions(+), 116 deletions(-)

diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
index 9be4628022..698f620694 100644
--- a/src/USER-MEAMC/meam.h
+++ b/src/USER-MEAMC/meam.h
@@ -2,6 +2,7 @@
 #define LMP_MEAM_H
 
 #include <cmath>
+#include <cstring>
 
 #define maxelt 5
 
@@ -26,7 +27,6 @@ private:
 
   // 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
@@ -65,7 +65,7 @@ private:
   // nrar,rdrar = spline coeff array parameters
 
   double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];
-  double Omega_meam[maxelt], Z_meam[maxelt];
+  double 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];
@@ -212,6 +212,31 @@ protected:
   void interpolate_meam(int);
 
 public:
+  //-----------------------------------------------------------------------------
+  // convert lattice spec to lattice_t
+  // only use single-element lattices if single=true
+  // return false on failure
+  // return true and set lat on success
+  static bool str_to_lat(const char* str, bool single, lattice_t& lat)
+  {
+    if (strcmp(str,"fcc") == 0) lat = FCC;
+    else if (strcmp(str,"bcc") == 0) lat = BCC;
+    else if (strcmp(str,"hcp") == 0) lat = HCP;
+    else if (strcmp(str,"dim") == 0) lat = DIM;
+    else if (strcmp(str,"dia") == 0) lat = DIA;
+    else {
+      if (single)
+        return false;
+      
+      if (strcmp(str,"b1")  == 0) lat = B1;
+      else if (strcmp(str,"c11") == 0) lat = C11;
+      else if (strcmp(str,"l12") == 0) lat = L12;
+      else if (strcmp(str,"b2")  == 0) lat = B2;
+      else return false;
+    }
+    return true;
+  }
+  
   void 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,
diff --git a/src/USER-MEAMC/meam_impl.cpp b/src/USER-MEAMC/meam_impl.cpp
index 74e8af1cde..96ef364fa3 100644
--- a/src/USER-MEAMC/meam_impl.cpp
+++ b/src/USER-MEAMC/meam_impl.cpp
@@ -38,7 +38,7 @@ MEAM::MEAM(Memory* mem)
 
   neltypes = 0;
   for (int i = 0; i < maxelt; i++) {
-    Omega_meam[i] = Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] =
+    Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] =
       beta1_meam[i]= beta2_meam[i] = beta3_meam[i] =
       t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] =
       rho_ref_meam[i] = ibar_meam[i] = ielt_meam[i] = 0.0;
diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp
index 37bfce5873..1c0fc3a53c 100644
--- a/src/USER-MEAMC/meam_setup_done.cpp
+++ b/src/USER-MEAMC/meam_setup_done.cpp
@@ -161,10 +161,10 @@ void
 MEAM::compute_pair_meam(void)
 {
 
-  double r /*ununsed:, temp*/;
+  double r;
   int j, a, b, nv2;
   double astar, frac, phizbl;
-  int n, nmax, Z1, Z2;
+  int n, Z1, Z2;
   double arat, rarat, scrn, scrn2;
   double phiaa, phibb /*unused:,phitmp*/;
   double C, s111, s112, s221, S11, S22;
@@ -230,9 +230,8 @@ MEAM::compute_pair_meam(void)
             Z1 = get_Zij(this->lattce_meam[a][a]);
             Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
                      this->Cmax_meam[a][a][a], arat, scrn);
-            nmax = 10;
             if (scrn > 0.0) {
-              for (n = 1; n <= nmax; n++) {
+              for (n = 1; n <= 10; n++) {
                 phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a);
               }
             }
@@ -242,9 +241,8 @@ MEAM::compute_pair_meam(void)
             Z1 = get_Zij(this->lattce_meam[b][b]);
             Z2 = get_Zij2(this->lattce_meam[b][b], this->Cmin_meam[b][b][b],
                      this->Cmax_meam[b][b][b], arat, scrn);
-            nmax = 10;
             if (scrn > 0.0) {
-              for (n = 1; n <= nmax; n++) {
+              for (n = 1; n <= 10; n++) {
                 phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b);
               }
             }
@@ -279,8 +277,7 @@ MEAM::compute_pair_meam(void)
             }
 
           } else {
-            nmax = 10;
-            for (n = 1; n <= nmax; n++) {
+            for (n = 1; n <= 10; n++) {
               this->phir[nv2][j] =
                 this->phir[nv2][j] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b);
             }
@@ -330,7 +327,7 @@ MEAM::phi_meam(double r, int a, int b)
   double Eu;
   double arat, scrn /*unused:,scrn2*/;
   int Z12, errorflag;
-  int n, nmax, Z1nn, Z2nn;
+  int n, Z1nn, Z2nn;
   lattice_t latta /*unused:,lattb*/;
   double rho_bkgd1, rho_bkgd2;
 
@@ -476,9 +473,8 @@ MEAM::phi_meam(double r, int a, int b)
     Z1nn = get_Zij(this->lattce_meam[a][a]);
     Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
              this->Cmax_meam[a][a][a], arat, scrn);
-    nmax = 10;
     if (scrn > 0.0) {
-      for (n = 1; n <= nmax; n++) {
+      for (n = 1; n <= 10; n++) {
         phiaa = phiaa + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a);
       }
     }
@@ -555,30 +551,38 @@ MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, dou
     *t12av = t12;
     *t22av = t22;
     *t32av = t32;
-  } else if (latt == FCC || latt == BCC || latt == DIA || latt == HCP || latt == B1 || latt == DIM || latt == B2) {
-    //     all neighbors are of the opposite type
-    *t11av = t12;
-    *t21av = t22;
-    *t31av = t32;
-    *t12av = t11;
-    *t22av = t21;
-    *t32av = t31;
-  } else {
-    a1 = r / this->re_meam[a][a] - 1.0;
-    a2 = r / this->re_meam[b][b] - 1.0;
-    rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
-    rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
-    if (latt == L12) {
-      rho01 = 8 * rhoa01 + 4 * rhoa02;
-      *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01;
+  } else switch (latt)  {
+    case FCC:
+    case BCC:
+    case DIA:
+    case HCP:
+    case B1:
+    case DIM:
+    case B2:
+      //     all neighbors are of the opposite type
+      *t11av = t12;
+      *t21av = t22;
+      *t31av = t32;
       *t12av = t11;
-      *t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01;
       *t22av = t21;
-      *t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01;
       *t32av = t31;
-    } else {
-      //      call error('Lattice not defined in get_tavref.')
-    }
+      break;
+    default:
+      a1 = r / this->re_meam[a][a] - 1.0;
+      a2 = r / this->re_meam[b][b] - 1.0;
+      rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
+      rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2);
+      if (latt == L12) {
+        rho01 = 8 * rhoa01 + 4 * rhoa02;
+        *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01;
+        *t12av = t11;
+        *t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01;
+        *t22av = t21;
+        *t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01;
+        *t32av = t31;
+      } else {
+        //      call error('Lattice not defined in get_tavref.')
+      }
   }
 }
 
@@ -627,59 +631,69 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
   *rho12 = 0.0;
   *rho22 = 0.0;
   *rho32 = 0.0;
-
-  if (lat == FCC) {
-    *rho01 = 12.0 * rhoa02;
-    *rho02 = 12.0 * rhoa01;
-  } else if (lat == BCC) {
-    *rho01 = 8.0 * rhoa02;
-    *rho02 = 8.0 * rhoa01;
-  } else if (lat == B1) {
-    *rho01 = 6.0 * rhoa02;
-    *rho02 = 6.0 * rhoa01;
-  } else if (lat == DIA) {
-    *rho01 = 4.0 * rhoa02;
-    *rho02 = 4.0 * rhoa01;
-    *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
-    *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
-  } else if (lat == HCP) {
-    *rho01 = 12 * rhoa02;
-    *rho02 = 12 * rhoa01;
-    *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
-    *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
-  } else if (lat == DIM) {
-    get_shpfcn(DIM, s);
-    *rho01 = rhoa02;
-    *rho02 = rhoa01;
-    *rho11 = s[0] * rhoa12 * rhoa12;
-    *rho12 = s[0] * rhoa11 * rhoa11;
-    *rho21 = s[1] * rhoa22 * rhoa22;
-    *rho22 = s[1] * rhoa21 * rhoa21;
-    *rho31 = s[2] * rhoa32 * rhoa32;
-    *rho32 = s[2] * rhoa31 * rhoa31;
-  } else if (lat == C11) {
-    *rho01 = rhoa01;
-    *rho02 = rhoa02;
-    *rho11 = rhoa11;
-    *rho12 = rhoa12;
-    *rho21 = rhoa21;
-    *rho22 = rhoa22;
-    *rho31 = rhoa31;
-    *rho32 = rhoa32;
-  } else if (lat == L12) {
-    *rho01 = 8 * rhoa01 + 4 * rhoa02;
-    *rho02 = 12 * rhoa01;
-    if (this->ialloy == 1) {
-      *rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]);
-      denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]);
-      if (denom > 0.)
-        *rho21 = *rho21 / denom * *rho01;
-    } else
-      *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
-  } else if (lat == B2) {
-    *rho01 = 8.0 * rhoa02;
-    *rho02 = 8.0 * rhoa01;
-  } else {
+  
+  switch (lat) {
+    case FCC:
+      *rho01 = 12.0 * rhoa02;
+      *rho02 = 12.0 * rhoa01;
+      break;
+    case BCC:
+      *rho01 = 8.0 * rhoa02;
+      *rho02 = 8.0 * rhoa01;
+      break;
+    case B1:
+      *rho01 = 6.0 * rhoa02;
+      *rho02 = 6.0 * rhoa01;
+      break;
+    case DIA:
+      *rho01 = 4.0 * rhoa02;
+      *rho02 = 4.0 * rhoa01;
+      *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
+      *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
+      break;
+    case HCP:
+      *rho01 = 12 * rhoa02;
+      *rho02 = 12 * rhoa01;
+      *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
+      *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
+      break;
+    case DIM:
+      get_shpfcn(DIM, s);
+      *rho01 = rhoa02;
+      *rho02 = rhoa01;
+      *rho11 = s[0] * rhoa12 * rhoa12;
+      *rho12 = s[0] * rhoa11 * rhoa11;
+      *rho21 = s[1] * rhoa22 * rhoa22;
+      *rho22 = s[1] * rhoa21 * rhoa21;
+      *rho31 = s[2] * rhoa32 * rhoa32;
+      *rho32 = s[2] * rhoa31 * rhoa31;
+      break;
+    case C11:
+      *rho01 = rhoa01;
+      *rho02 = rhoa02;
+      *rho11 = rhoa11;
+      *rho12 = rhoa12;
+      *rho21 = rhoa21;
+      *rho22 = rhoa22;
+      *rho31 = rhoa31;
+      *rho32 = rhoa32;
+      break;
+    case L12:
+      *rho01 = 8 * rhoa01 + 4 * rhoa02;
+      *rho02 = 12 * rhoa01;
+      if (this->ialloy == 1) {
+        *rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]);
+        denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]);
+        if (denom > 0.)
+          *rho21 = *rho21 / denom * *rho01;
+      } else
+        *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
+      break;
+    case B2:
+      *rho01 = 8.0 * rhoa02;
+      *rho02 = 8.0 * rhoa01;
+      break;
+    // default:
     //        call error('Lattice not defined in get_densref.')
   }
 
diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp
index b66aafa06c..2306ac53c6 100644
--- a/src/USER-MEAMC/meam_setup_global.cpp
+++ b/src/USER-MEAMC/meam_setup_global.cpp
@@ -49,17 +49,23 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub
     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);
-    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 {
+    switch(this->lattce_meam[i][i]) {
+      case FCC:
+        this->re_meam[i][i] = tmplat[i] / sqrt(2.0);
+        break;
+      case BCC:
+        this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
+        break;
+      case HCP:
+        this->re_meam[i][i] = tmplat[i];
+        break;
+      case DIM:
+        this->re_meam[i][i] = tmplat[i];
+        break;
+      case DIA:
+        this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
+        break;
+      //default:
       //           error
     }
   }
diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp
index dc2e194dea..bc366854b6 100644
--- a/src/USER-MEAMC/pair_meamc.cpp
+++ b/src/USER-MEAMC/pair_meamc.cpp
@@ -430,13 +430,9 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
     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");
+    
+    if (!MEAM::str_to_lat(words[1], true, lat[i]))
+      error->all(FLERR,"Unrecognized lattice type in MEAM file 1");
 
     // store parameters
 
@@ -523,6 +519,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
 
   int which;
   double value;
+  lattice_t latt;
   int nindex,index[3];
   int maxparams = 6;
   char **params = new char*[maxparams];
@@ -571,16 +568,9 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
     // map lattce_meam value to an integer
 
     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");
+      if (!MEAM::str_to_lat(params[nparams-1], false, latt))
+        error->all(FLERR,"Unrecognized lattice type in MEAM file 2");
+      value = latt;
     }
     else value = atof(params[nparams-1]);
 

From 2c65659421f935b4cceb874cc938241d87822b50 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Sat, 6 Jul 2019 18:23:53 +0200
Subject: [PATCH 3/9] MEAM/C: implement scaling factor for reversible scaling
 calculations

---
 src/USER-MEAMC/meam.h              |  4 ++--
 src/USER-MEAMC/meam_dens_final.cpp |  5 ++++-
 src/USER-MEAMC/meam_force.cpp      | 18 ++++++++++++++----
 src/USER-MEAMC/pair_meamc.cpp      | 27 ++++++++++++++++++++++-----
 src/USER-MEAMC/pair_meamc.h        |  2 ++
 5 files changed, 44 insertions(+), 12 deletions(-)

diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
index 698f620694..6e7dedfe2d 100644
--- a/src/USER-MEAMC/meam.h
+++ b/src/USER-MEAMC/meam.h
@@ -247,9 +247,9 @@ public:
   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);
   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);
+                       double* eatom, int ntype, int* type, int* fmap, double** scale, int& errorflag);
   void meam_force(int i, 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,
+                  double* eatom, int ntype, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh,
                   int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom);
 };
 
diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp
index 867106df88..f59ad1b08a 100644
--- a/src/USER-MEAMC/meam_dens_final.cpp
+++ b/src/USER-MEAMC/meam_dens_final.cpp
@@ -4,18 +4,20 @@ using namespace LAMMPS_NS;
 
 void
 MEAM::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)
+                      double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, int& errorflag)
 {
   int i, elti;
   int m;
   double rhob, G, dG, Gbar, dGbar, gam, shp[3], Z;
   double denom, rho_bkgd, Fl;
+  double scaleii;
 
   //     Complete the calculation of density
 
   for (i = 0; i < nlocal; i++) {
     elti = fmap[type[i]];
     if (elti >= 0) {
+      scaleii = scale[type[i]][type[i]];
       rho1[i] = 0.0;
       rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i];
       rho3[i] = 0.0;
@@ -113,6 +115,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
       Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]);
 
       if (eflag_either != 0) {
+        Fl *= scaleii;
         if (eflag_global != 0) {
           *eng_vdwl = *eng_vdwl + Fl;
         }
diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp
index 8c50296358..924d88f6ff 100644
--- a/src/USER-MEAMC/meam_force.cpp
+++ b/src/USER-MEAMC/meam_force.cpp
@@ -8,7 +8,7 @@ using namespace LAMMPS_NS;
 
 void
 MEAM::meam_force(int i, 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,
+                 double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh,
                  int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom)
 {
   int j, jn, k, kn, kk, m, n, p, q;
@@ -42,6 +42,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
   double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3;
   double dsij1, dsij2, force1, force2;
   double t1i, t2i, t3i, t1j, t2j, t3j;
+  double scaleij;
 
   third = 1.0 / 3.0;
   sixth = 1.0 / 6.0;
@@ -59,6 +60,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
   for (jn = 0; jn < numneigh; jn++) {
     j = firstneigh[jn];
     eltj = fmap[type[j]];
+    scaleij = scale[type[i]][type[j]];
 
     if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) {
 
@@ -84,11 +86,12 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
         recip = 1.0 / r;
 
         if (eflag_either != 0) {
+          double phi_sc = phi * scaleij;
           if (eflag_global != 0)
-            *eng_vdwl = *eng_vdwl + phi * sij;
+            *eng_vdwl = *eng_vdwl + phi_sc * sij;
           if (eflag_atom != 0) {
-            eatom[i] = eatom[i] + 0.5 * phi * sij;
-            eatom[j] = eatom[j] + 0.5 * phi * sij;
+            eatom[i] = eatom[i] + 0.5 * phi_sc * sij;
+            eatom[j] = eatom[j] + 0.5 * phi_sc * sij;
           }
         }
 
@@ -379,6 +382,13 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
         for (m = 0; m < 3; m++) {
           dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];
         }
+        if (!isone(scaleij)) {
+          dUdrij *= scaleij;
+          dUdsij *= scaleij;
+          dUdrijm[0] *= scaleij;
+          dUdrijm[1] *= scaleij;
+          dUdrijm[2] *= scaleij;
+        }
 
         //     Add the part of the force due to dUdrij and dUdsij
 
diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp
index bc366854b6..2a2a168fd4 100644
--- a/src/USER-MEAMC/pair_meamc.cpp
+++ b/src/USER-MEAMC/pair_meamc.cpp
@@ -56,6 +56,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp)
   elements = NULL;
   mass = NULL;
   meam_inst = new MEAM(memory);
+  scale = NULL;
 
   // set comm size needed by this Pair
 
@@ -79,6 +80,7 @@ PairMEAMC::~PairMEAMC()
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
+    memory->destroy(scale);
     delete [] map;
   }
 }
@@ -143,7 +145,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,map,errorflag);
+                   &eng_vdwl,eatom,ntype,type,map,scale,errorflag);
   if (errorflag) {
     char str[128];
     sprintf(str,"MEAM library error %d",errorflag);
@@ -164,7 +166,7 @@ void PairMEAMC::compute(int eflag, int vflag)
   for (ii = 0; ii < inum_half; ii++) {
     i = ilist_half[ii];
     meam_inst->meam_force(i,eflag_either,eflag_global,eflag_atom,
-                vflag_atom,&eng_vdwl,eatom,ntype,type,map,x,
+                vflag_atom,&eng_vdwl,eatom,ntype,type,map,scale,x,
                 numneigh_half[i],firstneigh_half[i],
                 numneigh_full[i],firstneigh_full[i],
                 offset,f,vptr);
@@ -183,6 +185,7 @@ void PairMEAMC::allocate()
 
   memory->create(setflag,n+1,n+1,"pair:setflag");
   memory->create(cutsq,n+1,n+1,"pair:cutsq");
+  memory->create(scale,n+1,n+1,"pair:scale");
 
   map = new int[n+1];
 }
@@ -267,13 +270,16 @@ void PairMEAMC::coeff(int narg, char **arg)
   // 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++)
+  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++;
       }
+      scale[i][j] = 1.0;
+    }
+  }
 
   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
 }
@@ -312,8 +318,10 @@ void PairMEAMC::init_list(int id, NeighList *ptr)
    init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */
 
-double PairMEAMC::init_one(int /*i*/, int /*j*/)
+double PairMEAMC::init_one(int i, int j)
 {
+  if (setflag[i][j] == 0) scale[i][j] = 1.0;
+  scale[j][i] = scale[i][j];
   return cutmax;
 }
 
@@ -773,3 +781,12 @@ void PairMEAMC::neigh_strip(int inum, int *ilist,
     for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;
   }
 }
+
+/* ---------------------------------------------------------------------- */
+
+void *PairMEAMC::extract(const char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"scale") == 0) return (void *) scale;
+  return NULL;
+}
diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h
index 31dd8ba022..102985f0ca 100644
--- a/src/USER-MEAMC/pair_meamc.h
+++ b/src/USER-MEAMC/pair_meamc.h
@@ -35,6 +35,7 @@ class PairMEAMC : public Pair {
   void init_style();
   void init_list(int, class NeighList *);
   double init_one(int, int);
+  virtual void *extract(const char *, int &);
 
   int pack_forward_comm(int, int *, double *, int, int *);
   void unpack_forward_comm(int, int, double *);
@@ -50,6 +51,7 @@ class PairMEAMC : public Pair {
   double *mass;                 // mass of each element
 
   int *map;                     // mapping from atom types (1-indexed) to elements (1-indexed)
+  double **scale;               // scaling factor for adapt
 
   void allocate();
   void read_files(char *, char *);

From 7e14dda789c071b27a26546beff3fb0bfc313cd5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Thu, 22 Aug 2019 11:36:25 +0200
Subject: [PATCH 4/9] MEAM/C: warn if z given and expected by lattice do not
 agree

---
 src/USER-MEAMC/meam.h         | 2 +-
 src/USER-MEAMC/pair_meamc.cpp | 5 ++++-
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
index 6e7dedfe2d..d19a653ee2 100644
--- a/src/USER-MEAMC/meam.h
+++ b/src/USER-MEAMC/meam.h
@@ -192,7 +192,6 @@ protected:
   static double erose(const double r, const double re, const double alpha, const double Ec, const double repuls, const double attrac, const int form);
 
   static void get_shpfcn(const lattice_t latt, double (&s)[3]);
-  static int get_Zij(const lattice_t latt);
   static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, double &a, double &S);
 protected:
   void meam_checkindex(int, int, int, int*, int*);
@@ -237,6 +236,7 @@ public:
     return true;
   }
   
+  static int get_Zij(const lattice_t latt);
   void 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,
diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp
index 2a2a168fd4..37e3fdd8c8 100644
--- a/src/USER-MEAMC/pair_meamc.cpp
+++ b/src/USER-MEAMC/pair_meamc.cpp
@@ -463,7 +463,10 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
     ibar[i] = atoi(words[18]);
 
     if (!isone(t0[i]))
-      error->all(FLERR,"Unsupported parameter in MEAM potential file");
+      error->all(FLERR,"Unsupported parameter in MEAM potential file: t0!=1");
+    
+    if (z[i] != MEAM::get_Zij(lat[i]))
+      error->warning(FLERR,"Unsupported parameter in MEAM potential file: z!=lat", true);
 
     nset++;
   }

From ce05ed4cca9f02250ea75c36f1c3465a9c612566 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Thu, 22 Aug 2019 13:21:56 +0200
Subject: [PATCH 5/9] MEAM/C: infer z parameter from lattice structure,
 eliminates possible user mistakes

---
 doc/txt/pair_meamc.txt               |  2 ++
 src/USER-MEAMC/meam.h                |  4 +---
 src/USER-MEAMC/meam_dens_final.cpp   |  2 +-
 src/USER-MEAMC/meam_funcs.cpp        |  6 +++---
 src/USER-MEAMC/meam_impl.cpp         |  2 +-
 src/USER-MEAMC/meam_setup_done.cpp   | 17 ++++++++---------
 src/USER-MEAMC/meam_setup_global.cpp |  3 +--
 src/USER-MEAMC/pair_meamc.cpp        |  5 +++--
 8 files changed, 20 insertions(+), 21 deletions(-)

diff --git a/doc/txt/pair_meamc.txt b/doc/txt/pair_meamc.txt
index 7c42e9d2f2..2831600e08 100644
--- a/doc/txt/pair_meamc.txt
+++ b/doc/txt/pair_meamc.txt
@@ -139,6 +139,8 @@ potential parameters:
 
 lat      = lattice structure of reference configuration
 z        = number of nearest neighbors in the reference structure
+             This field is only read for compatibility, the correct
+             value is inferred from the lattice structure
 ielement = atomic number
 atwt     = atomic weight
 alat     = lattice constant of reference structure
diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
index d19a653ee2..6bf3b93652 100644
--- a/src/USER-MEAMC/meam.h
+++ b/src/USER-MEAMC/meam.h
@@ -28,7 +28,6 @@ private:
   // Ec_meam = cohesive energy
   // re_meam = nearest-neighbor distance
   // 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)
@@ -65,7 +64,6 @@ private:
   // nrar,rdrar = spline coeff array parameters
 
   double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];
-  double 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];
@@ -237,7 +235,7 @@ public:
   }
   
   static int get_Zij(const lattice_t latt);
-  void meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha,
+  void meam_setup_global(int nelt, lattice_t* lat, 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);
diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp
index f59ad1b08a..eac8af44dc 100644
--- a/src/USER-MEAMC/meam_dens_final.cpp
+++ b/src/USER-MEAMC/meam_dens_final.cpp
@@ -54,7 +54,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
         gamma[i] = gamma[i] / (rho0[i] * rho0[i]);
       }
 
-      Z = this->Z_meam[elti];
+      Z = get_Zij(this->lattce_meam[elti][elti]);
 
       G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);
       if (errorflag != 0)
diff --git a/src/USER-MEAMC/meam_funcs.cpp b/src/USER-MEAMC/meam_funcs.cpp
index 312caff686..15e06ce78f 100644
--- a/src/USER-MEAMC/meam_funcs.cpp
+++ b/src/USER-MEAMC/meam_funcs.cpp
@@ -232,7 +232,7 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
 }
 
 //-----------------------------------------------------------------------------
-// Number of neighbors for the reference structure
+// Number of first neighbors for reference structure
 //
 int
 MEAM::get_Zij(const lattice_t latt)
@@ -244,12 +244,12 @@ MEAM::get_Zij(const lattice_t latt)
       return 8;
     case HCP:
       return 12;
-    case B1:
-      return 6;
     case DIA:
       return 4;
     case DIM:
       return 1;
+    case B1:
+      return 6;
     case C11:
       return 10;
     case L12:
diff --git a/src/USER-MEAMC/meam_impl.cpp b/src/USER-MEAMC/meam_impl.cpp
index 96ef364fa3..0a0b95e14a 100644
--- a/src/USER-MEAMC/meam_impl.cpp
+++ b/src/USER-MEAMC/meam_impl.cpp
@@ -38,7 +38,7 @@ MEAM::MEAM(Memory* mem)
 
   neltypes = 0;
   for (int i = 0; i < maxelt; i++) {
-    Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] =
+    A_meam[i] = rho0_meam[i] = beta0_meam[i] =
       beta1_meam[i]= beta2_meam[i] = beta3_meam[i] =
       t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] =
       rho_ref_meam[i] = ibar_meam[i] = ielt_meam[i] = 0.0;
diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp
index 1c0fc3a53c..8232c2aed5 100644
--- a/src/USER-MEAMC/meam_setup_done.cpp
+++ b/src/USER-MEAMC/meam_setup_done.cpp
@@ -327,7 +327,7 @@ MEAM::phi_meam(double r, int a, int b)
   double Eu;
   double arat, scrn /*unused:,scrn2*/;
   int Z12, errorflag;
-  int n, Z1nn, Z2nn;
+  int n, Z1nn;
   lattice_t latta /*unused:,lattb*/;
   double rho_bkgd1, rho_bkgd2;
 
@@ -338,6 +338,8 @@ MEAM::phi_meam(double r, int a, int b)
 
   // get number of neighbors in the reference structure
   //   Nref[i][j] = # of i's neighbors of type j
+  Z1 = get_Zij(this->lattce_meam[a][a]);
+  Z2 = get_Zij(this->lattce_meam[b][b]);
   Z12 = get_Zij(this->lattce_meam[a][b]);
 
   get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32);
@@ -398,8 +400,6 @@ MEAM::phi_meam(double r, int a, int b)
     //     If using mixing rule for t, apply to reference structure; else
     //     use precomputed values
     if (this->mix_ref_t == 1) {
-      Z1 = this->Z_meam[a];
-      Z2 = this->Z_meam[b];
       if (this->ibar_meam[a] <= 0)
         G1 = 1.0;
       else {
@@ -436,8 +436,8 @@ MEAM::phi_meam(double r, int a, int b)
       rho_bkgd2 = rho0_2;
     } else {
       if (this->bkgd_dyn == 1) {
-        rho_bkgd1 = this->rho0_meam[a] * this->Z_meam[a];
-        rho_bkgd2 = this->rho0_meam[b] * this->Z_meam[b];
+        rho_bkgd1 = this->rho0_meam[a] * Z1;
+        rho_bkgd2 = this->rho0_meam[b] * Z2;
       } else {
         rho_bkgd1 = this->rho_ref_meam[a];
         rho_bkgd2 = this->rho_ref_meam[b];
@@ -470,12 +470,11 @@ MEAM::phi_meam(double r, int a, int b)
   } else if (this->lattce_meam[a][b] == L12) {
     phiaa = phi_meam(r, a, a);
     //       account for second neighbor a-a potential here...
-    Z1nn = get_Zij(this->lattce_meam[a][a]);
-    Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
+    Z1nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
              this->Cmax_meam[a][a][a], arat, scrn);
     if (scrn > 0.0) {
       for (n = 1; n <= 10; n++) {
-        phiaa = phiaa + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a);
+        phiaa = phiaa + pow((-Z1nn * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, a);
       }
     }
     phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;
@@ -506,7 +505,7 @@ MEAM::compute_reference_density(void)
 
   // loop over element types
   for (a = 0; a < this->neltypes; a++) {
-    Z = (int)this->Z_meam[a];
+    Z = get_Zij(this->lattce_meam[a][a]);
     if (this->ibar_meam[a] <= 0)
       Gbar = 1.0;
     else {
diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp
index 2306ac53c6..56bac1d5af 100644
--- a/src/USER-MEAMC/meam_setup_global.cpp
+++ b/src/USER-MEAMC/meam_setup_global.cpp
@@ -18,7 +18,7 @@ static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) {
 }
 
 void
-MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* /*atwt*/, double* alpha,
+MEAM::meam_setup_global(int nelt, lattice_t* lat, 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)
@@ -32,7 +32,6 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub
   for (i = 0; i < nelt; i++) {
     this->lattce_meam[i][i] = lat[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];
diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp
index 37e3fdd8c8..77a5e9a461 100644
--- a/src/USER-MEAMC/pair_meamc.cpp
+++ b/src/USER-MEAMC/pair_meamc.cpp
@@ -465,8 +465,9 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
     if (!isone(t0[i]))
       error->all(FLERR,"Unsupported parameter in MEAM potential file: t0!=1");
     
+    // z given is ignored: if this is mismatched, we definitely won't do what the user said -> fatal error
     if (z[i] != MEAM::get_Zij(lat[i]))
-      error->warning(FLERR,"Unsupported parameter in MEAM potential file: z!=lat", true);
+      error->all(FLERR,"Mismatched parameter in MEAM potential file: z!=lat");
 
     nset++;
   }
@@ -478,7 +479,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile)
 
   // 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,ielement,atwt,alpha,b0,b1,b2,b3,
                        alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
 
   // set element masses

From 81fb0d613f97084f4d60751795204fe266c516d7 Mon Sep 17 00:00:00 2001
From: Sungkwang Mun <munsungkwang@gmail.com>
Date: Fri, 30 Aug 2019 17:54:43 +0200
Subject: [PATCH 6/9] * This commit includes the addition of new reference
 structures such as ch4: methane-like structure only for binary system. dia3:
 diamond structure with primary 1NN and secondary 3NN inteation tri: H2O-like
 structure that has an angle zig: zigzag structure with a uniform angle lin:
 linear structure (180 degree angle)

** tri, zig, and lin reference structures require angle information (in degree)
such as the following.
   theta = 109.5
---
 src/USER-MEAMC/meam.h                |  32 ++++-
 src/USER-MEAMC/meam_dens_final.cpp   |   4 +-
 src/USER-MEAMC/meam_force.cpp        |   5 +-
 src/USER-MEAMC/meam_funcs.cpp        |  90 +++++++++++--
 src/USER-MEAMC/meam_setup_done.cpp   | 192 +++++++++++++++++++++------
 src/USER-MEAMC/meam_setup_global.cpp |  10 +-
 src/USER-MEAMC/meam_setup_param.cpp  |  17 +++
 src/USER-MEAMC/pair_meamc.cpp        |   4 +-
 8 files changed, 294 insertions(+), 60 deletions(-)

diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
index 6bf3b93652..719e1af6f0 100644
--- a/src/USER-MEAMC/meam.h
+++ b/src/USER-MEAMC/meam.h
@@ -3,13 +3,14 @@
 
 #include <cmath>
 #include <cstring>
+#include "math_const.h"
 
 #define maxelt 5
 
 namespace LAMMPS_NS {
 class Memory;
 
-typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t;
+typedef enum { FCC, BCC, HCP, DIM, DIA, DIA3, B1, C11, L12, B2, CH4, LIN, ZIG, TRI } lattice_t;
 
 class MEAM
 {
@@ -63,6 +64,10 @@ private:
   // nr,dr = pair function discretization parameters
   // nrar,rdrar = spline coeff array parameters
 
+  // theta = angle between three atoms in line, zigzag, and trimer reference structures
+  // stheta_meam = sin(theta/2) in radian used in line, zigzag, and trimer reference structures
+  // ctheta_meam = cos(theta/2) in radian used in line, zigzag, and trimer reference structures
+
   double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];
   double A_meam[maxelt], alpha_meam[maxelt][maxelt], rho0_meam[maxelt];
   double delta_meam[maxelt][maxelt];
@@ -106,6 +111,10 @@ public:
   int maxneigh;
   double *scrfcn, *dscrfcn, *fcpair;
 
+  //angle for trimer, zigzag, line reference structures
+  double stheta_meam[maxelt][maxelt];
+  double ctheta_meam[maxelt][maxelt];
+
 protected:
   // meam_funcs.cpp
 
@@ -189,8 +198,12 @@ protected:
   double embedding(const double A, const double Ec, const double rhobar, double& dF) const;
   static double erose(const double r, const double re, const double alpha, const double Ec, const double repuls, const double attrac, const int form);
 
-  static void get_shpfcn(const lattice_t latt, double (&s)[3]);
-  static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, double &a, double &S);
+  static void get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3]);
+
+  static int get_Zij2(const lattice_t latt, const double cmin, const double cmax,
+                      const double sthe, double &a, double &S);
+  static int get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S);
+
 protected:
   void meam_checkindex(int, int, int, int*, int*);
   void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,
@@ -201,6 +214,7 @@ protected:
   void alloyparams();
   void compute_pair_meam();
   double phi_meam(double, int, int);
+  const double phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat);
   void compute_reference_density();
   void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double,
                   double, double, double, int, int, lattice_t);
@@ -221,19 +235,27 @@ public:
     else if (strcmp(str,"hcp") == 0) lat = HCP;
     else if (strcmp(str,"dim") == 0) lat = DIM;
     else if (strcmp(str,"dia") == 0) lat = DIA;
+    else if (strcmp(str,"dia3") == 0) lat = DIA3;
+    else if (strcmp(str,"lin") == 0) lat = LIN;
+    else if (strcmp(str,"zig") == 0) lat = ZIG;
+    else if (strcmp(str,"tri") == 0) lat = TRI;
     else {
       if (single)
         return false;
-      
+
       if (strcmp(str,"b1")  == 0) lat = B1;
       else if (strcmp(str,"c11") == 0) lat = C11;
       else if (strcmp(str,"l12") == 0) lat = L12;
       else if (strcmp(str,"b2")  == 0) lat = B2;
+      else if (strcmp(str,"ch4")  == 0) lat = CH4;
+      else if (strcmp(str,"lin")  == 0) lat =LIN;
+      else if (strcmp(str,"zig")  == 0) lat = ZIG;
+      else if (strcmp(str,"tri")  == 0) lat = TRI;
       else return false;
     }
     return true;
   }
-  
+
   static int get_Zij(const lattice_t latt);
   void meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* atwt, double* alpha,
                          double* b0, double* b1, double* b2, double* b3, double* alat, double* esub,
diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp
index eac8af44dc..fe9e74ca7c 100644
--- a/src/USER-MEAMC/meam_dens_final.cpp
+++ b/src/USER-MEAMC/meam_dens_final.cpp
@@ -59,7 +59,9 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
       G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);
       if (errorflag != 0)
         return;
-      get_shpfcn(this->lattce_meam[elti][elti], shp);
+
+      get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shp);
+
       if (this->ibar_meam[elti] <= 0) {
         Gbar = 1.0;
         dGbar = 0.0;
diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp
index 924d88f6ff..5a2f544d8b 100644
--- a/src/USER-MEAMC/meam_force.cpp
+++ b/src/USER-MEAMC/meam_force.cpp
@@ -290,8 +290,9 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
         }
 
         //     Compute derivatives of total density wrt rij, sij and rij(3)
-        get_shpfcn(this->lattce_meam[elti][elti], shpi);
-        get_shpfcn(this->lattce_meam[eltj][eltj], shpj);
+        get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpi);
+        get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpj);
+
         drhodr1 = dgamma1[i] * drho0dr1 +
           dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 +
                         dt3dr1 * rho3[i] + t3i * drho3dr1) -
diff --git a/src/USER-MEAMC/meam_funcs.cpp b/src/USER-MEAMC/meam_funcs.cpp
index 15e06ce78f..da984cfeb7 100644
--- a/src/USER-MEAMC/meam_funcs.cpp
+++ b/src/USER-MEAMC/meam_funcs.cpp
@@ -198,7 +198,7 @@ MEAM::erose(const double r, const double re, const double alpha, const double Ec
 // Shape factors for various configurations
 //
 void
-MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
+MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3])
 {
   switch (latt) {
     case FCC:
@@ -214,7 +214,9 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
       s[1] = 0.0;
       s[2] = 1.0 / 3.0;
       break;
+    case CH4: // CH4 actually needs shape factor for diamond for C, dimer for H
     case DIA:
+    case DIA3:
       s[0] = 0.0;
       s[1] = 0.0;
       s[2] = 32.0 / 9.0;
@@ -222,8 +224,20 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])
     case DIM:
       s[0] = 1.0;
       s[1] = 2.0 / 3.0;
-      //        s(3) = 1.d0
-      s[2] = 0.40;
+      //        s(3) = 1.d0 // this should be 0.4 unless (1-legendre) is multiplied in the density calc.
+      s[2] = 0.40; // this is (1-legendre) where legendre = 0.6 in dynamo is accounted.
+      break;
+    case LIN: //linear, theta being 180
+      s[0] = 0.0;
+      s[1] = 8.0 / 3.0; // 4*(co**4 + si**4 - 1.0/3.0) in zig become 4*(1-1/3)
+      s[2] = 0.0;
+      break;
+    case ZIG: //zig-zag
+    case TRI: //trimer e.g. H2O
+      s[0] = 4.0*pow(cthe,2);
+      s[1] = 4.0*(pow(cthe,4) + pow(sthe,4) - 1.0/3.0);
+      s[2] = 4.0*(pow(cthe,2) * (3*pow(sthe,4) + pow(cthe,4)));
+      s[2] = s[2] - 0.6*s[0]; //legend in dyn, 0.6 is default value.
       break;
     default:
       s[0] = 0.0;
@@ -245,6 +259,7 @@ MEAM::get_Zij(const lattice_t latt)
     case HCP:
       return 12;
     case DIA:
+    case DIA3:
       return 4;
     case DIM:
       return 1;
@@ -256,6 +271,12 @@ MEAM::get_Zij(const lattice_t latt)
       return 12;
     case B2:
       return 8;
+    case CH4: // DYNAMO currenly implemented this way while it needs two Z values, 4 and 1
+      return 4;
+    case LIN:
+    case ZIG:
+    case TRI:
+      return 2;
       //        call error('Lattice not defined in get_Zij.')
   }
   return 0;
@@ -263,11 +284,13 @@ MEAM::get_Zij(const lattice_t latt)
 
 //-----------------------------------------------------------------------------
 // Number of second neighbors for the reference structure
-//   a = distance ratio R1/R2
-//   S = second neighbor screening function
+//   a = distance ratio R1/R2 (a2nn in dynamo)
+//   numscr = number of atoms that screen the 2NN bond
+//   S = second neighbor screening function (xfac, a part of b2nn in dynamo)
 //
 int
-MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, double& a, double& S)
+MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax,
+                const double stheta, double& a, double& S)
 {
 
   double C, x, sijk;
@@ -299,7 +322,7 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl
     numscr = 2;
     break;
 
-  case DIA:
+  case DIA: // 2NN
     Zij2 = 12;
     a = sqrt(8.0 / 3.0);
     numscr = 1;
@@ -308,12 +331,35 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl
     }
     break;
 
+  case DIA3: // 3NN
+    Zij2 = 12;
+    a = sqrt(11.0 / 3.0);
+    numscr = 4;
+    if (cmin < 0.500001) {
+        //          call error('can not do 2NN MEAM for dia')
+    }
+    break;
+
+  case CH4: //does not have 2nn structure so it returns 0
+  case LIN: //line
   case DIM:
     //        this really shouldn't be allowed; make sure screening is zero
     a = 1.0;
     S = 0.0;
     return 0;
 
+  case TRI: //TRI
+    Zij2 = 1;
+    a = 2.0*stheta;
+    numscr=2;
+    break;
+
+  case ZIG: //zig-zag
+    Zij2 = 2;
+    a = 2.0*stheta;
+    numscr=1;
+    break;
+
   case L12:
     Zij2 = 6;
     a = sqrt(2.0);
@@ -335,10 +381,38 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl
   }
 
   // Compute screening for each first neighbor
-  C = 4.0 / (a * a) - 1.0;
+  if (latt==DIA3){ // special case for 3NN diamond structure
+	C = 1.0;
+  } else {
+	C = 4.0 / (a * a) - 1.0;
+  }
   x = (C - cmin) / (cmax - cmin);
   sijk = fcut(x);
   // There are numscr first neighbors screening the second neighbors
   S = MathSpecial::powint(sijk, numscr);
   return Zij2;
 }
+
+int
+MEAM::get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S)
+{
+
+  double x, sijk, C;
+  int numscr = 0, Zij2 = 0;
+  switch (latt) {
+  case ZIG: //zig-zag for b11s and b22s term
+  case TRI: //trimer for b11s
+    Zij2 = 2;
+    numscr=1;
+    break;
+  default:
+    // unknown lattic flag in get Zij
+    //        call error('Lattice not defined in get_Zij.')
+    break;
+  }
+  C = 1.0;
+  x = (C - cmin) / (cmax - cmin);
+  sijk = fcut(x);
+  S = MathSpecial::powint(sijk, numscr);
+  return Zij2;
+}
\ No newline at end of file
diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp
index 8232c2aed5..99b330a83a 100644
--- a/src/USER-MEAMC/meam_setup_done.cpp
+++ b/src/USER-MEAMC/meam_setup_done.cpp
@@ -103,6 +103,9 @@ MEAM::alloyparams(void)
         this->alpha_meam[i][j] = this->alpha_meam[j][i];
         this->lattce_meam[i][j] = this->lattce_meam[j][i];
         this->nn2_meam[i][j] = this->nn2_meam[j][i];
+	    // theta for lin,tri,zig references
+        this->stheta_meam[i][j] = this->stheta_meam[j][i];
+        this->ctheta_meam[i][j] = this->ctheta_meam[j][i];
         // If i<j and term is unset, use default values (e.g. mean of i-i and
         // j-j)
       } else if (j > i) {
@@ -161,7 +164,7 @@ void
 MEAM::compute_pair_meam(void)
 {
 
-  double r;
+  double r, b2nn, phi_val;
   int j, a, b, nv2;
   double astar, frac, phizbl;
   int n, Z1, Z2;
@@ -215,7 +218,7 @@ MEAM::compute_pair_meam(void)
         if (this->nn2_meam[a][b] == 1) {
           Z1 = get_Zij(this->lattce_meam[a][b]);
           Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
-                   this->Cmax_meam[a][a][b], arat, scrn);
+                     this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn);
 
           //     The B1, B2,  and L12 cases with NN2 have a trick to them; we need to
           //     compute the contributions from second nearest neighbors, like a-a
@@ -229,33 +232,26 @@ MEAM::compute_pair_meam(void)
             phiaa = phi_meam(rarat, a, a);
             Z1 = get_Zij(this->lattce_meam[a][a]);
             Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
-                     this->Cmax_meam[a][a][a], arat, scrn);
-            if (scrn > 0.0) {
-              for (n = 1; n <= 10; n++) {
-                phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a);
-              }
-            }
+                     this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);
+            phiaa+= phi_meam_series(scrn, Z1, Z2, a, a, rarat, arat);
 
             //               phi_bb
             phibb = phi_meam(rarat, b, b);
             Z1 = get_Zij(this->lattce_meam[b][b]);
             Z2 = get_Zij2(this->lattce_meam[b][b], this->Cmin_meam[b][b][b],
-                     this->Cmax_meam[b][b][b], arat, scrn);
-            if (scrn > 0.0) {
-              for (n = 1; n <= 10; n++) {
-                phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b);
-              }
-            }
+                     this->Cmax_meam[b][b][b], this->stheta_meam[b][b], arat, scrn);
+            phibb+= phi_meam_series(scrn, Z1, Z2, b, b, rarat, arat);
 
             if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||
                 this->lattce_meam[a][b] == DIA) {
               //     Add contributions to the B1 or B2 potential
               Z1 = get_Zij(this->lattce_meam[a][b]);
               Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
-                       this->Cmax_meam[a][a][b], arat, scrn);
+                       this->Cmax_meam[a][a][b], this->stheta_meam[a][b],  arat, scrn);
               this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;
               Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[b][b][a],
-                       this->Cmax_meam[b][b][a], arat, scrn2);
+                       this->Cmax_meam[b][b][a], this->stheta_meam[a][b], arat, scrn2);
+
               this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
 
             } else if (this->lattce_meam[a][b] == L12) {
@@ -277,10 +273,7 @@ MEAM::compute_pair_meam(void)
             }
 
           } else {
-            for (n = 1; n <= 10; n++) {
-              this->phir[nv2][j] =
-                this->phir[nv2][j] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b);
-            }
+            this->phir[nv2][j]+= phi_meam_series(scrn, Z1, Z2, a, b, r, arat);
           }
         }
 
@@ -325,11 +318,12 @@ MEAM::phi_meam(double r, int a, int b)
   double rho02, rho12, rho22, rho32;
   double scalfac, phiaa, phibb;
   double Eu;
-  double arat, scrn /*unused:,scrn2*/;
+  double arat, scrn, scrn2;
   int Z12, errorflag;
-  int n, Z1nn;
+  int n, Z1nn, Z2nn;
   lattice_t latta /*unused:,lattb*/;
   double rho_bkgd1, rho_bkgd2;
+  double b11s, b22s;
 
   double phi_m = 0.0;
 
@@ -403,14 +397,14 @@ MEAM::phi_meam(double r, int a, int b)
       if (this->ibar_meam[a] <= 0)
         G1 = 1.0;
       else {
-        get_shpfcn(this->lattce_meam[a][a], s1);
+        get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], s1);
         Gam1 = (s1[0] * t11av + s1[1] * t21av + s1[2] * t31av) / (Z1 * Z1);
         G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);
       }
       if (this->ibar_meam[b] <= 0)
         G2 = 1.0;
       else {
-        get_shpfcn(this->lattce_meam[b][b], s2);
+        get_shpfcn(this->lattce_meam[b][b], this->stheta_meam[b][b], this->ctheta_meam[b][b],  s2);
         Gam2 = (s2[0] * t12av + s2[1] * t22av + s2[2] * t32av) / (Z2 * Z2);
         G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);
       }
@@ -470,20 +464,46 @@ MEAM::phi_meam(double r, int a, int b)
   } else if (this->lattce_meam[a][b] == L12) {
     phiaa = phi_meam(r, a, a);
     //       account for second neighbor a-a potential here...
-    Z1nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
-             this->Cmax_meam[a][a][a], arat, scrn);
-    if (scrn > 0.0) {
-      for (n = 1; n <= 10; n++) {
-        phiaa = phiaa + pow((-Z1nn * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, a);
-      }
-    }
+    Z1nn = get_Zij(this->lattce_meam[a][a]);
+    Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
+             this->Cmax_meam[a][a][a], this->stheta_meam[a][b], arat, scrn);
+
+
+    phiaa += phi_meam_series(scrn, Z1nn, Z2nn, a, a, r, arat);
     phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;
 
+  } else if (this->lattce_meam[a][b] == CH4) {
+    phi_m = (5 * Eu - F1 - 4*F2)/4;
+
+  } else if (this->lattce_meam[a][b] == ZIG){
+      if (a==b){
+        phi_m = (2 * Eu - F1 - F2) / Z12;
+      } else{
+        Z1 = get_Zij(this->lattce_meam[a][b]);
+        Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);
+        b11s = -Z2/Z1*scrn;
+        Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], scrn2);
+        b22s = -Z2/Z1*scrn2;
+
+        phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);
+        phibb = phi_meam(2.0*this->stheta_meam[a][b]*r, b, b);
+        phi_m = (2.0*Eu - F1 - F2 + phiaa*b11s + phibb*b22s) / Z12;
+      }
+
+  } else if (this->lattce_meam[a][b] == TRI) {
+      if (a==b){
+        phi_m = (3.0*Eu - 2.0*F1 - F2) / Z12;
+     } else {
+        Z1 = get_Zij(this->lattce_meam[a][b]);
+        Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);
+        b11s = -Z2/Z1*scrn;
+        phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);
+        phi_m = (3.0*Eu - 2.0*F1 - F2 + phiaa*b11s) / Z12;
+      }
+
   } else {
-    //
     // potential is computed from Rose function and embedding energy
     phi_m = (2 * Eu - F1 - F2) / Z12;
-    //
   }
 
   // if r = 0, just return 0
@@ -494,6 +514,31 @@ MEAM::phi_meam(double r, int a, int b)
   return phi_m;
 }
 
+//----------------------------------------------------------------------c
+// Compute 2NN series terms for phi
+//   To avoid nan values of phir due to rapid decrease of b2nn^n or/and
+//   argument of phi_meam, i.e. r*arat^n, in some cases (3NN dia with low Cmin value)
+//
+const double
+MEAM::phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat)
+{
+  double phi_sum = 0.0;
+  double b2nn, phi_val;
+  if (scrn > 0.0) {
+    b2nn = -Z2*scrn/Z1;
+    for (int n = 1; n <= 10; n++) {
+      phi_val = MathSpecial::powint(b2nn,n) * phi_meam(r * MathSpecial::powint(arat, n), a, b);
+      if (iszero(phi_val)) {
+        // once either term becomes zero at some point, all folliwng will also be zero
+        // necessary to avoid numerical error (nan or infty) due to exponential decay in phi_meam
+        break;
+      }
+      phi_sum += phi_val;
+    }
+  }
+  return phi_sum;
+}
+
 //----------------------------------------------------------------------c
 // Compute background density for reference structure of each element
 void
@@ -509,7 +554,7 @@ MEAM::compute_reference_density(void)
     if (this->ibar_meam[a] <= 0)
       Gbar = 1.0;
     else {
-      get_shpfcn(this->lattce_meam[a][a], shp);
+      get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], shp);
       gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] + this->t3_meam[a] * shp[2]) / (Z * Z);
       Gbar = G_gam(gam, this->ibar_meam[a], errorflag);
     }
@@ -524,7 +569,7 @@ MEAM::compute_reference_density(void)
     //     screening)
     if (this->nn2_meam[a][a] == 1) {
       Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
-               this->Cmax_meam[a][a][a], arat, scrn);
+               this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);
       rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));
       rho0 = rho0 + Z2 * rho0_2nn * scrn;
     }
@@ -554,10 +599,15 @@ MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, dou
     case FCC:
     case BCC:
     case DIA:
+    case DIA3:
     case HCP:
     case B1:
     case DIM:
     case B2:
+    case CH4:
+    case LIN:
+    case ZIG:
+    case TRI:
       //     all neighbors are of the opposite type
       *t11av = t12;
       *t21av = t22;
@@ -603,7 +653,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
   double a1, a2;
   double s[3];
   lattice_t lat;
-  int Zij2nn;
+  int Zij,Zij2nn;
   double rhoa01nn, rhoa02nn;
   double rhoa01, rhoa11, rhoa21, rhoa31;
   double rhoa02, rhoa12, rhoa22, rhoa32;
@@ -624,13 +674,15 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
 
   lat = this->lattce_meam[a][b];
 
+  Zij = get_Zij(lat);
+
   *rho11 = 0.0;
   *rho21 = 0.0;
   *rho31 = 0.0;
   *rho12 = 0.0;
   *rho22 = 0.0;
   *rho32 = 0.0;
-  
+
   switch (lat) {
     case FCC:
       *rho01 = 12.0 * rhoa02;
@@ -645,6 +697,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
       *rho02 = 6.0 * rhoa01;
       break;
     case DIA:
+    case DIA3:
       *rho01 = 4.0 * rhoa02;
       *rho02 = 4.0 * rhoa01;
       *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
@@ -657,7 +710,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
       *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
       break;
     case DIM:
-      get_shpfcn(DIM, s);
+      get_shpfcn(DIM, 0, 0, s);
       *rho01 = rhoa02;
       *rho02 = rhoa01;
       *rho11 = s[0] * rhoa12 * rhoa12;
@@ -692,13 +745,71 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
       *rho01 = 8.0 * rhoa02;
       *rho02 = 8.0 * rhoa01;
       break;
+    case CH4:
+      *rho01 = 4.0 * rhoa02; //in assumption that 'a' represent carbon
+      *rho02 = rhoa01;	//in assumption that 'b' represent hydrogen
+
+      get_shpfcn(DIM, 0, 0, s);	//H
+      *rho12 = s[0] * rhoa11 * rhoa11;
+      *rho22 = s[1] * rhoa21 * rhoa21;
+      *rho32 = s[2] * rhoa31 * rhoa31;
+
+      get_shpfcn(CH4, 0, 0, s); //C
+      *rho11 = s[0] * rhoa12 * rhoa12;
+      *rho21 = s[1] * rhoa22 * rhoa22;
+      *rho31 = s[2] * rhoa32 * rhoa32;
+      break;
+    case LIN:
+      *rho01 = rhoa02*Zij;
+      *rho02 = rhoa01*Zij;
+
+      get_shpfcn(LIN, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
+      *rho12 = s[0] * rhoa11 * rhoa11;
+      *rho22 = s[1] * rhoa21 * rhoa21;
+      *rho32 = s[2] * rhoa31 * rhoa31;
+      *rho11 = s[0] * rhoa12 * rhoa12;
+      *rho21 = s[1] * rhoa22 * rhoa22;
+      *rho31 = s[2] * rhoa32 * rhoa32;
+      break;
+    case ZIG:
+      *rho01 = rhoa02*Zij;
+      *rho02 = rhoa01*Zij;
+
+      get_shpfcn(ZIG, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
+      *rho12 = s[0] * rhoa11 * rhoa11;
+      *rho22 = s[1] * rhoa21 * rhoa21;
+      *rho32 = s[2] * rhoa31 * rhoa31;
+      *rho11 = s[0] * rhoa12 * rhoa12;
+      *rho21 = s[1] * rhoa22 * rhoa22;
+      *rho31 = s[2] * rhoa32 * rhoa32;
+      break;
+    case TRI:
+      *rho01 = rhoa02;
+      *rho02 = rhoa01*Zij;
+
+      get_shpfcn(TRI, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);
+      *rho12 = s[0] * rhoa11 * rhoa11;
+      *rho22 = s[1] * rhoa21 * rhoa21;
+      *rho32 = s[2] * rhoa31 * rhoa31;
+      s[0] = 1.0;
+      s[1] = 2.0/3.0;
+      s[2] = 1.0 - 0.6*s[0];
+
+      *rho11 = s[0] * rhoa12 * rhoa12;
+      *rho21 = s[1] * rhoa22 * rhoa22;
+      *rho31 = s[2] * rhoa32 * rhoa32;
+      break;
+
+
     // default:
     //        call error('Lattice not defined in get_densref.')
   }
 
   if (this->nn2_meam[a][b] == 1) {
 
-    Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], arat, scrn);
+
+    Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b],
+                      this->stheta_meam[a][b], arat, scrn);
 
     a1 = arat * r / this->re_meam[a][a] - 1.0;
     a2 = arat * r / this->re_meam[b][b] - 1.0;
@@ -725,7 +836,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
       *rho01 = *rho01 + Zij2nn * scrn * rhoa01nn;
 
       //     Assume Zij2nn and arat don't depend on order, but scrn might
-      Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], arat, scrn);
+      Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a],
+                        this->stheta_meam[a][b], arat, scrn);
       *rho02 = *rho02 + Zij2nn * scrn * rhoa02nn;
     }
   }
diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp
index 56bac1d5af..ac5718914c 100644
--- a/src/USER-MEAMC/meam_setup_global.cpp
+++ b/src/USER-MEAMC/meam_setup_global.cpp
@@ -56,12 +56,15 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt*
         this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0;
         break;
       case HCP:
-        this->re_meam[i][i] = tmplat[i];
-        break;
       case DIM:
+      case CH4:
+      case LIN:
+      case ZIG:
+      case TRI:
         this->re_meam[i][i] = tmplat[i];
         break;
       case DIA:
+      case DIA3:
         this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0;
         break;
       //default:
@@ -87,4 +90,7 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt*
   this->emb_lin_neg = 0;
   this->bkgd_dyn = 0;
   this->erose_form = 0;
+  // for trimer, zigzag, line refernece structure, sungkwang
+  setall2d(this->stheta_meam, 1.0); // stheta = sin(theta/2*pi/180) where theta is 180, so 1.0
+  setall2d(this->ctheta_meam, 0.0); // stheta = cos(theta/2*pi/180) where theta is 180, so 0
 }
diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp
index 585b5b5712..70d8b63598 100644
--- a/src/USER-MEAMC/meam_setup_param.cpp
+++ b/src/USER-MEAMC/meam_setup_param.cpp
@@ -1,6 +1,7 @@
 #include "meam.h"
 #include <algorithm>
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 //
 //     do a sanity check on index parameters
@@ -46,6 +47,8 @@ MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr
 //     18 = zbl_meam
 //     19 = emb_lin_neg
 //     20 = bkgd_dyn
+//     21 = theta
+
 
 void
 MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag)
@@ -203,6 +206,20 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3
       this->bkgd_dyn = (int)value;
       break;
 
+    //     21 = theta
+    // see alloyparams(void) in meam_setup_done.cpp
+    case 21:
+      // const double PI = 3.141592653589793238463;
+      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]);
+      // we don't use theta, instead stheta and ctheta
+      this->stheta_meam[i1][i2] = sin(value/2*MY_PI/180.0);
+      this->ctheta_meam[i1][i2] = cos(value/2*MY_PI/180.0);
+      break;
+
     default:
       *errorflag = 1;
   }
diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp
index 77a5e9a461..9733ee2dad 100644
--- a/src/USER-MEAMC/pair_meamc.cpp
+++ b/src/USER-MEAMC/pair_meamc.cpp
@@ -33,13 +33,13 @@ using namespace LAMMPS_NS;
 
 #define MAXLINE 1024
 
-static const int nkeywords = 21;
+static const int nkeywords = 22;
 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"};
+  "emb_lin_neg","bkgd_dyn", "theta"};
 
 /* ---------------------------------------------------------------------- */
 

From 1be9364a895ceb67cbc352ad114fc4924d02c563 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Tue, 3 Dec 2019 19:47:23 +0100
Subject: [PATCH 7/9] MEAM/C: document new reference structures

---
 doc/src/pair_meamc.rst | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/doc/src/pair_meamc.rst b/doc/src/pair_meamc.rst
index daa10906c7..ac858ce24b 100644
--- a/doc/src/pair_meamc.rst
+++ b/doc/src/pair_meamc.rst
@@ -233,15 +233,20 @@ where
    Cmin(I,J,K) = Cmin screening parameter when I-J pair is screened
                  by K (I<=J); default = 2.0
    lattce(I,J) = lattice structure of I-J reference structure:
-                   dia = diamond (interlaced fcc for alloy)
                    fcc = face centered cubic
                    bcc = body centered cubic
-                   dim = dimer
-                   b1  = rock salt (NaCl structure)
                    hcp = hexagonal close-packed
+                   dim = dimer
+                   dia = diamond (interlaced fcc for alloy)
+                   dia3= diamond structure with primary 1NN and secondary 3NN interation
+                   b1  = rock salt (NaCl structure)
                    c11 = MoSi2 structure
                    l12 = Cu3Au structure (lower case L, followed by 12)
                    b2  = CsCl structure (interpenetrating simple cubic)
+                   ch4 = methane-like structure, only for binary system
+                   lin = linear structure (180 degree angle)
+                   zig = zigzag structure with a uniform angle
+                   tri = H2O-like structure that has an angle
    nn2(I,J)    = turn on second-nearest neighbor MEAM formulation for
                  I-J pair (see for example :ref:`(Lee) <Lee>`).
                    0 = second-nearest neighbor formulation off
@@ -254,6 +259,8 @@ where
    zbl(I,J)    = blend the MEAM I-J pair potential with the ZBL potential for small
                  atom separations :ref:`(ZBL) <ZBL>`
                    default = 1
+   theta(I,J)  = angle between three atoms in line, zigzag, and trimer reference structures in degrees
+                   default = 180
    gsmooth_factor  = factor determining the length of the G-function smoothing
                      region; only significant for ibar=0 or ibar=4.
                          99.0 = short smoothing region, sharp step

From a231197c10a246448720f7e99cca61256a6c1e73 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Tue, 17 Dec 2019 21:47:15 +0100
Subject: [PATCH 8/9] MEAM/C: remove unused variable alias

---
 src/USER-MEAMC/meam_force.cpp | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp
index 5a2f544d8b..2b6832e155 100644
--- a/src/USER-MEAMC/meam_force.cpp
+++ b/src/USER-MEAMC/meam_force.cpp
@@ -17,7 +17,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
   double v[6], fi[3], fj[3];
   double third, sixth;
   double pp, dUdrij, dUdsij, dUdrijm[3], force, forcem;
-  double r, recip, phi, phip;
+  double recip, phi, phip;
   double sij;
   double a1, a1i, a1j, a2, a2i, a2j;
   double a3i, a3j;
@@ -71,7 +71,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
       rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
       if (rij2 < this->cutforcesq) {
         rij = sqrt(rij2);
-        r = rij;
+        recip = 1.0 / rij;
 
         //     Compute phi and phip
         ind = this->eltind[elti][eltj];
@@ -80,10 +80,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
         kk = std::min(kk, this->nrar - 2);
         pp = pp - kk;
         pp = std::min(pp, 1.0);
-        phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp +
-          this->phirar[ind][kk];
+        phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + this->phirar[ind][kk];
         phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk];
-        recip = 1.0 / r;
 
         if (eflag_either != 0) {
           double phi_sc = phi * scaleij;

From cb20cb9f25c9f18cf9891dcfef8795a2a68b6006 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Tue, 17 Dec 2019 22:23:19 +0100
Subject: [PATCH 9/9] Add range checks for MathSpecial::fm_exp

---
 src/math_special.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/math_special.cpp b/src/math_special.cpp
index bf11a1ad45..d4abc36f25 100644
--- a/src/math_special.cpp
+++ b/src/math_special.cpp
@@ -538,6 +538,8 @@ double MathSpecial::exp2_x86(double x)
 double MathSpecial::fm_exp(double x)
 {
 #if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__)
+    if (x < -1022.0/FM_DOUBLE_LOG2OFE) return 0;
+    if (x > 1023.0/FM_DOUBLE_LOG2OFE) return INFINITY;
     return exp2_x86(FM_DOUBLE_LOG2OFE * x);
 #else
     return ::exp(x);