From 51743655703cbfa3e31c41894316df5f451343bb Mon Sep 17 00:00:00 2001
From: ares20105 <ares201005@gmail.com>
Date: Wed, 20 Nov 2019 13:21:41 -0700
Subject: [PATCH 01/22] add force modifications. Previous code does not call
 force modify, thus the dynamical matrix calculation does not work for other
 potentials defined via modify

---
 src/USER-PHONON/dynamical_matrix.cpp | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/USER-PHONON/dynamical_matrix.cpp b/src/USER-PHONON/dynamical_matrix.cpp
index 1495219124..6bb843c16e 100644
--- a/src/USER-PHONON/dynamical_matrix.cpp
+++ b/src/USER-PHONON/dynamical_matrix.cpp
@@ -19,6 +19,7 @@
 #include "improper.h"
 #include "kspace.h"
 #include "update.h"
+#include "modify.h"
 #include "neighbor.h"
 #include "pair.h"
 #include "timer.h"
@@ -385,6 +386,7 @@ void DynamicalMatrix::displace_atom(int local_idx, int direction, int magnitude)
 void DynamicalMatrix::update_force()
 {
     force_clear();
+    int n_post_force = modify->n_post_force;
 
     if (pair_compute_flag) {
         force->pair->compute(eflag,vflag);
@@ -405,6 +407,11 @@ void DynamicalMatrix::update_force()
         comm->reverse_comm();
         timer->stamp(Timer::COMM);
     }
+
+    // force modifications, 
+    if (n_post_force) modify->post_force(vflag);
+    timer->stamp(Timer::MODIFY);
+
     ++ update->nsteps;
 }
 

From 36e102516fcc36b7cd6741dd8987263260087e9f Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Mon, 2 Dec 2019 10:34:03 -0500
Subject: [PATCH 02/22] angle constraint bugfix

ghost atom fix
---
 src/USER-MISC/fix_bond_react.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
index 33aab5b4ad..3b5719156a 100644
--- a/src/USER-MISC/fix_bond_react.cpp
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -1675,6 +1675,7 @@ int FixBondReact::check_constraints()
         delx1 = x[atom1][0] - x[atom2][0];
         dely1 = x[atom1][1] - x[atom2][1];
         delz1 = x[atom1][2] - x[atom2][2];
+        domain->minimum_image(delx1,dely1,delz1); // ghost location fix
         rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
         r1 = sqrt(rsq1);
 
@@ -1682,6 +1683,7 @@ int FixBondReact::check_constraints()
         delx2 = x[atom3][0] - x[atom2][0];
         dely2 = x[atom3][1] - x[atom2][1];
         delz2 = x[atom3][2] - x[atom2][2];
+        domain->minimum_image(delx2,dely2,delz2); // ghost location fix
         rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
         r2 = sqrt(rsq2);
 

From 819fe9ec5668319fb17c465e8998c3cf21f9379c Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Mon, 2 Dec 2019 12:27:57 -0500
Subject: [PATCH 03/22] add option to enforce atom chirality

---
 src/USER-MISC/fix_bond_react.cpp | 94 +++++++++++++++++++++++++++++++-
 src/USER-MISC/fix_bond_react.h   | 17 +++++-
 2 files changed, 109 insertions(+), 2 deletions(-)

diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
index 3b5719156a..dbdb8ef0a7 100644
--- a/src/USER-MISC/fix_bond_react.cpp
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -36,6 +36,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
 #include "group.h"
 #include "citeme.h"
 #include "math_const.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "error.h"
 
@@ -297,6 +298,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms");
   memory->create(custom_edges,max_natoms,nreacts,"bond/react:custom_edges");
   memory->create(delete_atoms,max_natoms,nreacts,"bond/react:delete_atoms");
+  memory->create(chiral_atoms,max_natoms,6,nreacts,"bond/react:chiral_atoms");
 
   for (int j = 0; j < nreacts; j++)
     for (int i = 0; i < max_natoms; i++) {
@@ -304,6 +306,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
       if (update_edges_flag[j] == 1) custom_edges[i][j] = 1;
       else custom_edges[i][j] = 0;
       delete_atoms[i][j] = 0;
+      for (int k = 0; k < 6; k++) {
+        chiral_atoms[i][k][j] = 0;
+      }
     }
 
   // read all map files afterward
@@ -430,6 +435,7 @@ FixBondReact::~FixBondReact()
   memory->destroy(landlocked_atoms);
   memory->destroy(custom_edges);
   memory->destroy(delete_atoms);
+  memory->destroy(chiral_atoms);
 
   memory->destroy(nevery);
   memory->destroy(cutsq);
@@ -1649,10 +1655,11 @@ evaluate constraints: return 0 if any aren't satisfied
 int FixBondReact::check_constraints()
 {
   tagint atom1,atom2,atom3;
-  double delx,dely,delz,rsq;
+  double unwrap[3],delx,dely,delz,rsq;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double rsq1,rsq2,r1,r2,c,t,prrhob;
 
+  imageint *image = atom->image;
   double **x = atom->x;
 
   for (int i = 0; i < nconstraints; i++) {
@@ -1701,6 +1708,30 @@ int FixBondReact::check_constraints()
       }
     }
   }
+
+  // let's also check chirality within 'check_constraint'
+  for (int i = 0; i < onemol->natoms; i++) {
+    if (chiral_atoms[i][0][rxnID] == 1) {
+      double my4coords[12];
+      // already ensured, by transitive property, that chiral simulation atom has four neighs
+      for (int j = 0; j < 4; j++) {
+        atom1 = atom->map(glove[i][1]);
+        // loop over known types involved in chiral center
+        for (int jj = 0; jj < 4; jj++) {
+          if (atom->type[atom->map(xspecial[atom1][j])] == chiral_atoms[i][jj+2][rxnID]) {
+            atom2 = atom->map(xspecial[atom1][j]);
+            domain->unmap(x[atom2],image[atom2],unwrap);
+            for (int k = 0; k < 3; k++) {
+              my4coords[3*jj+k] = unwrap[k];
+            }
+            break;
+          }
+        }
+      }
+      if (get_chirality(my4coords) != chiral_atoms[i][1][rxnID]) return 0;
+    }
+  }
+
   return 1;
 }
 
@@ -1741,6 +1772,33 @@ double FixBondReact::get_temperature()
   return t;
 }
 
+/* ----------------------------------------------------------------------
+return handedness (1 or -1) of a chiral center, given ordered set of coordinates
+------------------------------------------------------------------------- */
+
+int FixBondReact::get_chirality(double four_coords[12])
+{
+  // define oriented plane with first three coordinates
+  double vec1[3],vec2[3],vec3[3],vec4[3],mean3[3],dot;
+
+  for (int i = 0; i < 3; i++) {
+    vec1[i] = four_coords[i]-four_coords[i+3];
+    vec2[i] = four_coords[i+3]-four_coords[i+6];
+  }
+
+  MathExtra::cross3(vec1,vec2,vec3);
+
+  for (int i = 0; i < 3; i++) {
+    mean3[i] = (four_coords[i] + four_coords[i+3] +
+                                 four_coords[i+6])/3;
+    vec4[i] = four_coords[i+9] - mean3[i];
+  }
+
+  dot = MathExtra::dot3(vec3,vec4);
+  dot = dot/fabs(dot);
+  return (int) dot;
+}
+
 /* ----------------------------------------------------------------------
   Get xspecials for current molecule templates
 ------------------------------------------------------------------------- */
@@ -2869,6 +2927,7 @@ void FixBondReact::read(int myrxn)
     }
     else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);
     else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);
+    else if (strstr(line,"chiralIDs")) sscanf(line,"%d",&nchiral);
     else if (strstr(line,"constraints")) {
       sscanf(line,"%d",&nconstr);
       memory->grow(constraints,nconstraints+nconstr,MAXCONARGS,"bond/react:constraints");
@@ -2900,6 +2959,8 @@ void FixBondReact::read(int myrxn)
       CustomEdges(line, myrxn);
     } else if (strcmp(keyword,"DeleteIDs") == 0) {
       DeleteAtoms(line, myrxn);
+    } else if (strcmp(keyword,"ChiralIDs") == 0) {
+      ChiralCenters(line, myrxn);
     } else if (strcmp(keyword,"Constraints") == 0) {
       Constraints(line, myrxn);
     } else error->one(FLERR,"Bond/react: Unknown section in map file");
@@ -2977,6 +3038,37 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn)
   }
 }
 
+void FixBondReact::ChiralCenters(char *line, int myrxn)
+{
+  int tmp;
+  for (int i = 0; i < nchiral; i++) {
+    readline(line);
+    sscanf(line,"%d",&tmp);
+    chiral_atoms[tmp-1][0][myrxn] = 1;
+    if (onemol->xflag == 0)
+      error->one(FLERR,"Bond/react: Molecule template 'Coords' section required for chiralIDs keyword");
+    if ((int) onemol_nxspecial[tmp-1][0] != 4)
+      error->one(FLERR,"Bond/react: Chiral atoms must have exactly four first neighbors");
+    for (int j = 0; j < 4; j++) {
+      for (int k = j+1; k < 4; k++) {
+        if (onemol->type[onemol_xspecial[tmp-1][j]-1] ==
+            onemol->type[onemol_xspecial[tmp-1][k]-1])
+          error->one(FLERR,"Bond/react: First neighbors of chiral atoms must be of mutually different types");
+      }
+    }
+    // record order of atom types, and coords
+    double my4coords[12];
+    for (int j = 0; j < 4; j++) {
+      chiral_atoms[tmp-1][j+2][myrxn] = onemol->type[onemol_xspecial[tmp-1][j]-1];
+      for (int k = 0; k < 3; k++) {
+        my4coords[3*j+k] = onemol->x[onemol_xspecial[tmp-1][j]-1][k];
+      }
+    }
+    // get orientation
+    chiral_atoms[tmp-1][1][myrxn] = get_chirality(my4coords);
+  }
+}
+
 void FixBondReact::Constraints(char *line, int myrxn)
 {
   double tmp[MAXCONARGS];
diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h
index f59e051ed1..a5bf2bc587 100644
--- a/src/USER-MISC/fix_bond_react.h
+++ b/src/USER-MISC/fix_bond_react.h
@@ -110,7 +110,7 @@ class FixBondReact : public Fix {
 
   int *ibonding,*jbonding;
   int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
-  int nedge,nequivalent,ncustom,ndelete,nconstr; // # edge, equivalent, custom atoms in mapping file
+  int nedge,nequivalent,ncustom,ndelete,nchiral,nconstr; // # edge, equivalent, custom atoms in mapping file
   int attempted_rxn; // there was an attempt!
   int *local_rxn_count;
   int *ghostly_rxn_count;
@@ -126,6 +126,7 @@ class FixBondReact : public Fix {
   int **landlocked_atoms; // all atoms at least three bonds away from edge atoms
   int **custom_edges; // atoms in molecule templates with incorrect valences
   int **delete_atoms; // atoms in pre-reacted templates to delete
+  int ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types
 
   int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors
   tagint **xspecial,**onemol_xspecial,**twomol_xspecial; // full 1-4 neighbor list
@@ -149,6 +150,7 @@ class FixBondReact : public Fix {
   void Equivalences(char *,int);
   void CustomEdges(char *,int);
   void DeleteAtoms(char *,int);
+  void ChiralCenters(char *,int);
   void Constraints(char *,int);
 
   void make_a_guess ();
@@ -159,6 +161,7 @@ class FixBondReact : public Fix {
   void ring_check();
   int check_constraints();
   double get_temperature();
+  int get_chirality(double[12]); // get handedness given an ordered set of coordinates
 
   void open(char *);
   void readline(char *);
@@ -249,6 +252,18 @@ E: Bond/react: A deleted atom cannot remain bonded to an atom that is not delete
 
 Self-explanatory.
 
+E: Bond/react: First neighbors of chiral atoms must be of mutually different types
+
+Self-explanatory.
+
+E: Bond/react: Chiral atoms must have exactly four first neighbors
+
+Self-explanatory.
+
+E: Bond/react: Molecule template 'Coords' section required for chiralIDs keyword
+
+The coordinates of atoms in the pre-reacted template are used to determine chirality.
+
 E: Bond/react special bond generation overflow
 
 The number of special bonds per-atom created by a reaction exceeds the

From 4cb797e63dc3dd1bdab4ec3faf7002206ed8d241 Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Mon, 2 Dec 2019 14:25:58 -0500
Subject: [PATCH 04/22] correct image atom mistake

---
 src/USER-MISC/fix_bond_react.cpp | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
index dbdb8ef0a7..2995a596a5 100644
--- a/src/USER-MISC/fix_bond_react.cpp
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -1655,11 +1655,10 @@ evaluate constraints: return 0 if any aren't satisfied
 int FixBondReact::check_constraints()
 {
   tagint atom1,atom2,atom3;
-  double unwrap[3],delx,dely,delz,rsq;
+  double delx,dely,delz,rsq;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double rsq1,rsq2,r1,r2,c,t,prrhob;
 
-  imageint *image = atom->image;
   double **x = atom->x;
 
   for (int i = 0; i < nconstraints; i++) {
@@ -1720,9 +1719,9 @@ int FixBondReact::check_constraints()
         for (int jj = 0; jj < 4; jj++) {
           if (atom->type[atom->map(xspecial[atom1][j])] == chiral_atoms[i][jj+2][rxnID]) {
             atom2 = atom->map(xspecial[atom1][j]);
-            domain->unmap(x[atom2],image[atom2],unwrap);
+            atom2 = domain->closest_image(atom1,atom2);
             for (int k = 0; k < 3; k++) {
-              my4coords[3*jj+k] = unwrap[k];
+              my4coords[3*jj+k] = x[atom2][k];
             }
             break;
           }

From 28fda045260a885f1c2bc9a2a4bcf4d552c9832c Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Mon, 2 Dec 2019 15:11:59 -0500
Subject: [PATCH 05/22] chiral centers docs

---
 doc/src/fix_bond_react.rst | 26 +++++++++++++++++++-------
 1 file changed, 19 insertions(+), 7 deletions(-)

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index ab63aa232c..7f72d1f428 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -20,9 +20,9 @@ Syntax
 * the common keyword/values may be appended directly after 'bond/react'
 * this applies to all reaction specifications (below)
 * common\_keyword = *stabilization*
-  
+
   .. parsed-literal::
-  
+
        *stabilization* values = *no* or *yes* *group-ID* *xmax*
          *no* = no reaction site stabilization
          *yes* = perform reaction site stabilization
@@ -40,9 +40,9 @@ Syntax
 * map\_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates
 * zero or more individual keyword/value pairs may be appended to each react argument
 * individual\_keyword = *prob* or *max\_rxn* or *stabilize\_steps* or *update\_edges*
-  
+
   .. parsed-literal::
-  
+
          *prob* values = fraction seed
            fraction = initiate reaction with this probability if otherwise eligible
            seed = random number seed (positive integer)
@@ -269,10 +269,11 @@ The optional keywords are 'edgeIDs', 'deleteIDs', 'customIDs' and
 
    N *edgeIDs* = # of edge atoms N in the pre-reacted molecule template
    N *deleteIDs* = # of atoms N that are specified for deletion
+   N *chiralIDs* = # of specified chiral centers N
    N *customIDs* = # of atoms N that are specified for a custom update
    N *constraints* = # of specified reaction constraints N
 
-The body of the map file contains two mandatory sections and four
+The body of the map file contains two mandatory sections and five
 optional sections. The first mandatory section begins with the keyword
 'BondingIDs' and lists the atom IDs of the bonding atom pair in the
 pre-reacted molecule template. The second mandatory section begins
@@ -284,12 +285,14 @@ molecule template. The first optional section begins with the keyword
 'EdgeIDs' and lists the atom IDs of edge atoms in the pre-reacted
 molecule template. The second optional section begins with the keyword
 'DeleteIDs' and lists the atom IDs of pre-reaction template atoms to
-delete. The third optional section begins with the keyword 'Custom
+delete. The third optional section begins with the keyword 'ChiralIDs'
+lists the atom IDs of chiral atoms whose handedness should be
+enforced. The fourth optional section begins with the keyword 'Custom
 Edges' and allows for forcing the update of a specific atom's atomic
 charge. The first column is the ID of an atom near the edge of the
 pre-reacted molecule template, and the value of the second column is
 either 'none' or 'charges.' Further details are provided in the
-discussion of the 'update\_edges' keyword. The fourth optional section
+discussion of the 'update\_edges' keyword. The fifth optional section
 begins with the keyword 'Constraints' and lists additional criteria
 that must be satisfied in order for the reaction to occur. Currently,
 there are three types of constraints available, as discussed below.
@@ -332,6 +335,15 @@ A sample map file is given below:
 ----------
 
 
+The handedness of atoms that are chiral centers can be enforced by
+listing their IDs in the ChiralIDs section. A chiral atom must be
+bonded to four atoms with mutually different atom types. This feature
+uses the coordinates and types of the involved atoms in the
+pre-reaction template to determine handedness. Three atoms bonded to
+the chiral center are arbitrarily chosen, to define an oriented plane,
+and the relative position of the fourth bonded atom determines the
+chiral center's handedness.
+
 Any number of additional constraints may be specified in the
 Constraints section of the map file. The constraint of type 'distance'
 has syntax as follows:

From e69e96ffbe900d48e6b7d7e592af004903d86016 Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Mon, 2 Dec 2019 15:20:01 -0500
Subject: [PATCH 06/22] Update Errors_messages.rst

---
 doc/src/Errors_messages.rst | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst
index 71de83afa9..a1795a6776 100644
--- a/doc/src/Errors_messages.rst
+++ b/doc/src/Errors_messages.rst
@@ -518,6 +518,16 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
 *Bond/react: A deleted atom cannot remain bonded to an atom that is not deleted*
    Self-explanatory.
 
+*Bond/react: First neighbors of chiral atoms must be of mutually different types*
+   Self-explanatory.
+
+*Bond/react: Chiral atoms must have exactly four first neighbors*
+   Self-explanatory.
+
+*Bond/react: Molecule template 'Coords' section required for chiralIDs keyword*
+   The coordinates of atoms in the pre-reacted template are used to determine
+   chirality.
+
 *Bond/react special bond generation overflow*
    The number of special bonds per-atom created by a reaction exceeds the
    system setting. See the read\_data or create\_box command for how to

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 07/22] 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 08/22] 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 09/22] 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 10/22] 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 11/22] 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 12/22] * 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 99ba15bf6afb3dc6843e2bf18cb939ce3ade276e Mon Sep 17 00:00:00 2001
From: Trung Nguyen <ndactrung@gmail.com>
Date: Wed, 11 Dec 2019 16:29:42 -0600
Subject: [PATCH 13/22] Updated README and added -L$(CUDA_HOME)/lib64/stubs to
 the Makefile's

---
 lib/gpu/Makefile.fermi           | 2 +-
 lib/gpu/Makefile.lammps.standard | 2 +-
 lib/gpu/Makefile.linux           | 2 +-
 lib/gpu/Makefile.linux.double    | 2 +-
 lib/gpu/Makefile.linux.mixed     | 2 +-
 lib/gpu/Makefile.linux.single    | 2 +-
 lib/gpu/Makefile.linux_multi     | 2 +-
 lib/gpu/Makefile.mac             | 2 +-
 lib/gpu/Makefile.serial          | 2 +-
 lib/gpu/Makefile.shannon         | 2 +-
 lib/gpu/Makefile.xk7             | 2 +-
 lib/gpu/README                   | 5 +++--
 12 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/lib/gpu/Makefile.fermi b/lib/gpu/Makefile.fermi
index ce5ccaaf78..81ab0bcde5 100644
--- a/lib/gpu/Makefile.fermi
+++ b/lib/gpu/Makefile.fermi
@@ -4,7 +4,7 @@ EXTRAMAKE = Makefile.lammps.standard
 CUDA_ARCH = -arch=sm_35
 CUDA_PRECISION = -D_SINGLE_DOUBLE
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64 -Xlinker -rpath -Xlinker $(CUDA_HOME)/lib64 -lcudart
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -Xlinker -rpath -Xlinker $(CUDA_HOME)/lib64/stubs -lcudart
 CUDA_OPTS = -DUNIX -O3 --use_fast_math --ftz=true
 
 CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -I$(CUDA_HOME)/include
diff --git a/lib/gpu/Makefile.lammps.standard b/lib/gpu/Makefile.lammps.standard
index 1ceb93b342..9526e8e373 100644
--- a/lib/gpu/Makefile.lammps.standard
+++ b/lib/gpu/Makefile.lammps.standard
@@ -7,5 +7,5 @@ endif
 
 gpu_SYSINC =
 gpu_SYSLIB =  -lcudart -lcuda
-gpu_SYSPATH = -L$(CUDA_HOME)/lib64
+gpu_SYSPATH = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 
diff --git a/lib/gpu/Makefile.linux b/lib/gpu/Makefile.linux
index 7001c6d8b9..a26fbe114c 100644
--- a/lib/gpu/Makefile.linux
+++ b/lib/gpu/Makefile.linux
@@ -54,7 +54,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
 CUDA_PRECISION = -D_SINGLE_DOUBLE
 
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math $(LMP_INC) -Xcompiler -fPIC
 
 CUDR_CPP = mpicxx -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK -DOMPI_SKIP_MPICXX=1 -fPIC
diff --git a/lib/gpu/Makefile.linux.double b/lib/gpu/Makefile.linux.double
index e65647f160..05f083697d 100644
--- a/lib/gpu/Makefile.linux.double
+++ b/lib/gpu/Makefile.linux.double
@@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
 CUDA_PRECISION = -D_DOUBLE_DOUBLE
 
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math
 
 CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK
diff --git a/lib/gpu/Makefile.linux.mixed b/lib/gpu/Makefile.linux.mixed
index a036b84ee3..ca414f1fc1 100644
--- a/lib/gpu/Makefile.linux.mixed
+++ b/lib/gpu/Makefile.linux.mixed
@@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
 CUDA_PRECISION = -D_SINGLE_DOUBLE
 
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math
 
 CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK
diff --git a/lib/gpu/Makefile.linux.single b/lib/gpu/Makefile.linux.single
index 808647cea7..1b349faac2 100644
--- a/lib/gpu/Makefile.linux.single
+++ b/lib/gpu/Makefile.linux.single
@@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
 CUDA_PRECISION = -D_SINGLE_SINGLE
 
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math
 
 CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK
diff --git a/lib/gpu/Makefile.linux_multi b/lib/gpu/Makefile.linux_multi
index 2f75ca0e2b..ba50170f39 100644
--- a/lib/gpu/Makefile.linux_multi
+++ b/lib/gpu/Makefile.linux_multi
@@ -40,7 +40,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
 CUDA_PRECISION = -D_SINGLE_DOUBLE
 
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math $(LMP_INC) -Xcompiler "-fPIC -std=c++98"
 
 CUDR_CPP = mpicxx -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK -DOMPI_SKIP_MPICXX=1 -fPIC
diff --git a/lib/gpu/Makefile.mac b/lib/gpu/Makefile.mac
index f9f8d5179a..b571e8d85a 100644
--- a/lib/gpu/Makefile.mac
+++ b/lib/gpu/Makefile.mac
@@ -14,7 +14,7 @@ NVCC = nvcc -m64
 CUDA_ARCH = -arch=sm_11
 CUDA_PRECISION = -D_SINGLE_SINGLE
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib
+CUDA_LIB = -L$(CUDA_HOME)/lib -L$(CUDA_HOME)/lib/stubs
 CUDA_OPTS = -DUNIX -DUCL_NO_EXIT -O3 --use_fast_math
 
 CUDR_CPP = mpic++ -m64
diff --git a/lib/gpu/Makefile.serial b/lib/gpu/Makefile.serial
index 99153fc471..b0cfb3c86b 100644
--- a/lib/gpu/Makefile.serial
+++ b/lib/gpu/Makefile.serial
@@ -34,7 +34,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
 CUDA_PRECISION = -D_SINGLE_DOUBLE
 
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64 -L../../src/STUBS -lmpi_stubs
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs -L../../src/STUBS -lmpi_stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math $(LMP_INC)
 
 CUDR_CPP = g++ -DMPI_GERYON -DUCL_NO_EXIT -fPIC -I../../src/STUBS
diff --git a/lib/gpu/Makefile.shannon b/lib/gpu/Makefile.shannon
index 22c2dc89d7..b14bd01ec2 100644
--- a/lib/gpu/Makefile.shannon
+++ b/lib/gpu/Makefile.shannon
@@ -32,7 +32,7 @@ LMP_INC = -DLAMMPS_SMALLBIG
 CUDA_PRECISION = -D_DOUBLE_DOUBLE
 
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math
 
 CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK
diff --git a/lib/gpu/Makefile.xk7 b/lib/gpu/Makefile.xk7
index 0b9f029399..f3050575e1 100644
--- a/lib/gpu/Makefile.xk7
+++ b/lib/gpu/Makefile.xk7
@@ -14,7 +14,7 @@ CUDA_ARCH = -arch=sm_35
 
 CUDA_PRECISION = -D_SINGLE_DOUBLE
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib64/stubs
 CUDA_OPTS = -DUNIX -O3 --use_fast_math
 
 CUDR_CPP = CC -DCUDA_PROXY -DMPI_GERYON -DUCL_NO_EXIT -DMPICH_IGNORE_CXX_SEEK
diff --git a/lib/gpu/README b/lib/gpu/README
index 6cc386270d..92be2d408c 100644
--- a/lib/gpu/README
+++ b/lib/gpu/README
@@ -3,7 +3,7 @@
                   --------------------------------
                      
                        W. Michael Brown (ORNL)
-                        Trung Dac Nguyen (ORNL)
+                        Trung Dac Nguyen (ORNL/Northwestern)
                           Peng Wang (NVIDIA)
                        Axel Kohlmeyer (Temple)
                          Steve Plimpton (SNL)
@@ -129,7 +129,8 @@ that ships with the CUDA toolkit, but also with the CUDA driver library
 on the head node of a GPU cluster, this library may not be installed,
 so you may need to copy it over from one of the compute nodes (best into
 this directory). Recent CUDA toolkits starting from CUDA 9 provide a dummy
-libcuda.so library, that can be used for linking (but not for running).
+libcuda.so library (typically under $(CUDA_HOME)/lib64/stubs), that can be used for
+linking (but not for running).
 
 The gpu library supports 3 precision modes as determined by 
 the CUDA_PRECISION variable:

From bde8b57f0bae7fb8876fbf63a030e7e1c8347372 Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Wed, 11 Dec 2019 23:10:11 -0700
Subject: [PATCH 14/22] Update fix_bond_react.rst

---
 doc/src/fix_bond_react.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 7f72d1f428..b2866eb9c7 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -253,7 +253,7 @@ A discussion of correctly handling this is also provided on the
 The map file is a text document with the following format:
 
 A map file has a header and a body. The header of map file the
-contains one mandatory keyword and four optional keywords. The
+contains one mandatory keyword and five optional keywords. The
 mandatory keyword is 'equivalences':
 
 

From 7e8a04d9857c9ed3685d7d933bf2f5a6981b3bb5 Mon Sep 17 00:00:00 2001
From: Trung Nguyen <ndactrung@gmail.com>
Date: Thu, 12 Dec 2019 09:57:49 -0600
Subject: [PATCH 15/22] More minor update to README

---
 lib/gpu/README | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/gpu/README b/lib/gpu/README
index 92be2d408c..2d98749a40 100644
--- a/lib/gpu/README
+++ b/lib/gpu/README
@@ -130,7 +130,7 @@ on the head node of a GPU cluster, this library may not be installed,
 so you may need to copy it over from one of the compute nodes (best into
 this directory). Recent CUDA toolkits starting from CUDA 9 provide a dummy
 libcuda.so library (typically under $(CUDA_HOME)/lib64/stubs), that can be used for
-linking (but not for running).
+linking.
 
 The gpu library supports 3 precision modes as determined by 
 the CUDA_PRECISION variable:

From 4c3ec145f3c1f2f456e4be1ed90fc36394da2114 Mon Sep 17 00:00:00 2001
From: Trung Nguyen <ndactrung@gmail.com>
Date: Tue, 17 Dec 2019 14:05:39 -0600
Subject: [PATCH 16/22] Corrected the wrong use of the stubs path in -rpath

---
 lib/gpu/Makefile.fermi | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/gpu/Makefile.fermi b/lib/gpu/Makefile.fermi
index 81ab0bcde5..ce5ccaaf78 100644
--- a/lib/gpu/Makefile.fermi
+++ b/lib/gpu/Makefile.fermi
@@ -4,7 +4,7 @@ EXTRAMAKE = Makefile.lammps.standard
 CUDA_ARCH = -arch=sm_35
 CUDA_PRECISION = -D_SINGLE_DOUBLE
 CUDA_INCLUDE = -I$(CUDA_HOME)/include
-CUDA_LIB = -L$(CUDA_HOME)/lib64 -Xlinker -rpath -Xlinker $(CUDA_HOME)/lib64/stubs -lcudart
+CUDA_LIB = -L$(CUDA_HOME)/lib64 -Xlinker -rpath -Xlinker $(CUDA_HOME)/lib64 -lcudart
 CUDA_OPTS = -DUNIX -O3 --use_fast_math --ftz=true
 
 CUDR_CPP = mpic++ -DMPI_GERYON -DUCL_NO_EXIT -I$(CUDA_HOME)/include

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 17/22] 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 18/22] 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 19/22] 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);

From 3be04e4671734ce6dab12bda8d5af5351431e3ac Mon Sep 17 00:00:00 2001
From: Michael Brown <michael.w.brown@intel.com>
Date: Wed, 18 Dec 2019 01:27:39 -0800
Subject: [PATCH 20/22] Bug fix for USER-INTEL package with triclinic neighbor
 builds.

---
 src/USER-INTEL/npair_intel.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/USER-INTEL/npair_intel.cpp b/src/USER-INTEL/npair_intel.cpp
index a82d3f29e5..4256e03b3c 100644
--- a/src/USER-INTEL/npair_intel.cpp
+++ b/src/USER-INTEL/npair_intel.cpp
@@ -360,7 +360,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
             if (THREE) ttag[u] = tag[j];
           }
 
-          if (FULL == 0 || TRI == 1) {
+          if (FULL == 0 && TRI != 1) {
             icount = 0;
             istart = ncount;
             IP_PRE_edge_align(istart, sizeof(int));
@@ -392,7 +392,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list,
         // ---------------------- Loop over i bin
 
         int n = 0;
-        if (FULL == 0 || TRI == 1) {
+        if (FULL == 0 && TRI != 1) {
           #if defined(LMP_SIMD_COMPILER)
           #pragma vector aligned
           #pragma ivdep

From d05f32d15228c76a216d2097cb68e2bc14d069f2 Mon Sep 17 00:00:00 2001
From: Oliver Henrich <ohenrich@users.noreply.github.com>
Date: Wed, 18 Dec 2019 11:31:02 +0000
Subject: [PATCH 21/22] Added oxRNA2 files to .gitignore

---
 src/.gitignore | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/.gitignore b/src/.gitignore
index 5848874d94..1b7dc7e6b2 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -280,6 +280,8 @@
 /bond_oxdna_fene.h
 /bond_oxdna2_fene.cpp
 /bond_oxdna2_fene.h
+/bond_oxrna2_fene.cpp
+/bond_oxrna2_fene.h
 /bond_quartic.cpp
 /bond_quartic.h
 /bond_table.cpp
@@ -982,6 +984,8 @@
 /pair_oxdna_*.h
 /pair_oxdna2_*.cpp
 /pair_oxdna2_*.h
+/pair_oxrna2_*.cpp
+/pair_oxrna2_*.h
 /mf_oxdna.h
 /pair_peri_eps.cpp
 /pair_peri_eps.h

From 7d16783366c33dc19e8cf1afb67d005dfa5f567f Mon Sep 17 00:00:00 2001
From: Stan Moore <stamoor@sandia.gov>
Date: Thu, 19 Dec 2019 09:31:30 -0700
Subject: [PATCH 22/22] Fix issue in Kokkos minimize

---
 src/min.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/min.h b/src/min.h
index fe7f0be8c6..ba47f95410 100644
--- a/src/min.h
+++ b/src/min.h
@@ -39,9 +39,9 @@ class Min : protected Pointers {
   virtual bigint memory_usage() {return 0;}
   void modify_params(int, char **);
   virtual int modify_param(int, char **) {return 0;}
-  double fnorm_sqr();
-  double fnorm_inf();
-  double fnorm_max();
+  virtual double fnorm_sqr();
+  virtual double fnorm_inf();
+  virtual double fnorm_max();
 
   enum{TWO,MAX,INF};