Create package USER-MEAMC

Step 1: very literal translation of lib/meam

examples/meam/in.meamc

@@ -0,0 +1,30 @@
+# Test of MEAM potential for SiC system
+
+units		metal
+boundary	p p p
+
+atom_style	atomic
+
+read_data	data.meam
+
+pair_style	meam/c
+pair_coeff	* * library.meam Si C SiC.meam Si C
+
+neighbor	0.3 bin
+neigh_modify	delay 10
+
+fix		1 all nve
+thermo		10
+timestep	0.001
+
+#dump		1 all atom 50 dump.meam
+
+#dump		2 all image 10 image.*.jpg element element &
+#		axes yes 0.8 0.02 view 60 -30
+#dump_modify	2 pad 3 element Si C
+
+#dump		3 all movie 10 movie.mpg element element &
+#		axes yes 0.8 0.02 view 60 -30
+#dump_modify	3 pad 3 element Si C
+
+run		100

examples/meam/in.meamc.shear

@@ -0,0 +1,79 @@
+# 3d metal shear simulation
+
+units		metal
+boundary	s s p
+
+atom_style	atomic
+lattice		fcc 3.52
+region		box block 0 16.0 0 10.0 0 2.828427
+create_box	3 box
+
+lattice		fcc 3.52 orient	x 1 0 0 orient y 0 1 1 orient z 0 -1 1 &
+		origin 0.5 0 0 
+create_atoms	1 box
+
+pair_style	meam/c
+pair_coeff	* * library.meam Ni4 Ni.meam Ni4 Ni4 Ni4
+
+neighbor	0.3 bin
+neigh_modify	delay 5
+
+region		lower block INF INF INF 0.9 INF INF
+region		upper block INF INF 6.1 INF INF INF
+group		lower region lower
+group		upper region upper
+group		boundary union lower upper
+group		mobile subtract all boundary
+
+set		group lower type 2
+set		group upper type 3
+
+# void
+
+#region		void cylinder z 8 5 2.5 INF INF
+#delete_atoms	region void
+
+# temp controllers
+
+compute		new3d mobile temp
+compute		new2d mobile temp/partial 0 1 1
+
+# equilibrate
+
+velocity	mobile create 300.0 5812775 temp new3d
+fix		1 all nve
+fix		2 boundary setforce 0.0 0.0 0.0
+
+fix		3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
+fix_modify	3 temp new3d
+
+thermo		25
+thermo_modify	temp new3d
+
+timestep	0.001
+run		100
+
+# shear
+
+velocity	upper set 1.0 0 0
+velocity	mobile ramp vx 0.0 1.0 y 1.4 8.6 sum yes
+
+unfix		3
+fix		3 mobile temp/rescale 10 300.0 300.0 10.0 1.0
+fix_modify	3 temp new2d
+
+#dump		1 all atom 500 dump.meam.shear
+
+#dump		2 all image 100 image.*.jpg type type &
+#		axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
+#dump_modify	2 pad 4
+
+#dump		3 all movie 100 movie.mpg type type &
+#		axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0
+#dump_modify	3 pad 4
+
+thermo		100
+thermo_modify	temp new2d
+
+reset_timestep	0
+run		3000

src/Makefile

@@ -59,7 +59,7 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \
 
 PACKUSER = user-atc user-awpmd user-cgdna user-cgsdk user-colvars \
 	   user-diffraction user-dpd user-drude user-eff user-fep user-h5md \
-	   user-intel user-lb user-manifold user-mgpt user-misc user-molfile \
+	   user-intel user-lb user-manifold user-meamc user-mgpt user-misc user-molfile \
 	   user-netcdf user-omp user-phonon user-qmmm user-qtb \
 	   user-quip user-reaxc user-smd user-smtbq user-sph user-tally \
 	   user-vtk

src/USER-MEAMC/Install.sh

@@ -0,0 +1,47 @@
+# Install/unInstall package files in LAMMPS
+# mode = 0/1/2 for uninstall/install/update
+
+# this is default Install.sh for all packages
+# if package has an auxiliary library or a file with a dependency,
+# then package dir has its own customized Install.sh
+
+mode=$1
+
+# enforce using portable C locale
+LC_ALL=C
+export LC_ALL
+
+# arg1 = file, arg2 = file it depends on
+
+action () {
+  if (test $mode = 0) then
+    rm -f ../$1
+  elif (! cmp -s $1 ../$1) then
+    if (test -z "$2" || test -e ../$2) then
+      cp $1 ..
+      if (test $mode = 2) then
+        echo "  updating src/$1"
+      fi
+    fi
+  elif (test -n "$2") then
+    if (test ! -e ../$2) then
+      rm -f ../$1
+    fi
+  fi
+}
+
+# all package files with no dependencies
+
+for file in *.cpp *.h; do
+  test -f ${file} && action $file
+done
+
+
+# additional files
+
+for file in meam_*.c; do
+  test -f ${file} && action ${file}
+done
+
+action fm_exp.c
+action meam.h

src/USER-MEAMC/README

@@ -0,0 +1,26 @@
+This package implements the MEAM/C potential as a LAMMPS pair style.
+
++==============================================================================+
+
+This package is a translation of the MEAM package to native C. See
+that package as well as the Fortran code distributed in lib/meam for
+the original sources and information.
+
+
+Translation by
+ Sebastian Hütter, sebastian.huetter@ovgu.de
+ Institute of Materials and Joining Technology
+ Otto-von-Guericke University Magdeburg, Germany
+
+The original Fortran implementation was created by
+  Greg Wagner (while at Sandia, now at Northwestern U).
+
++==============================================================================+
+
+Use "make yes-user-meamc" to enable this package when building LAMMPS.
+
+In your LAMMPS input script, specifiy
+  pair_style meam/c
+to enable the use of this implementation. All parameters, input files and
+outputs are exactly identical to these used with pair_style meam.
+

src/USER-MEAMC/fm_exp.c

@@ -0,0 +1,133 @@
+/*
+   Copyright (c) 2012,2013   Axel Kohlmeyer <akohlmey@gmail.com>
+   All rights reserved.
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   * Redistributions of source code must retain the above copyright
+     notice, this list of conditions and the following disclaimer.
+   * Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+   * Neither the name of the <organization> nor the
+     names of its contributors may be used to endorse or promote products
+     derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+/* faster versions of 2**x, e**x, and 10**x in single and double precision.
+ *
+ * Based on the Cephes math library 2.8
+ */
+
+#include <math.h>
+#include <stdint.h>
+
+/* internal definitions for the fastermath library */
+
+/* IEEE 754 double precision floating point data manipulation */
+typedef union
+{
+    double   f;
+    uint64_t u;
+    struct {int32_t  i0,i1;};
+}  udi_t;
+#define FM_DOUBLE_BIAS 1023
+#define FM_DOUBLE_EMASK 2146435072
+#define FM_DOUBLE_MBITS 20
+#define FM_DOUBLE_MMASK 1048575
+#define FM_DOUBLE_EZERO 1072693248
+
+/* generate 2**num in floating point by bitshifting */
+#define FM_DOUBLE_INIT_EXP(var,num)                 \
+    var.i0 = 0;                                     \
+    var.i1 = (((int) num) + FM_DOUBLE_BIAS) << 20
+
+/* double precision constants */
+#define FM_DOUBLE_LOG2OFE  1.4426950408889634074
+#define FM_DOUBLE_LOGEOF2  6.9314718055994530942e-1
+#define FM_DOUBLE_LOG2OF10 3.32192809488736234789
+#define FM_DOUBLE_LOG10OF2 3.0102999566398119521e-1
+#define FM_DOUBLE_LOG10OFE 4.3429448190325182765e-1
+#define FM_DOUBLE_SQRT2    1.41421356237309504880
+#define FM_DOUBLE_SQRTH    0.70710678118654752440
+
+/* optimizer friendly implementation of exp2(x).
+ *
+ * strategy:
+ *
+ * split argument into an integer part and a fraction:
+ * ipart = floor(x+0.5);
+ * fpart = x - ipart;
+ *
+ * compute exp2(ipart) from setting the ieee754 exponent
+ * compute exp2(fpart) using a pade' approximation for x in [-0.5;0.5[
+ *
+ * the result becomes: exp2(x) = exp2(ipart) * exp2(fpart)
+ */
+
+static const double fm_exp2_q[] = {
+/*  1.00000000000000000000e0, */
+    2.33184211722314911771e2,
+    4.36821166879210612817e3
+};
+static const double fm_exp2_p[] = {
+    2.30933477057345225087e-2,
+    2.02020656693165307700e1,
+    1.51390680115615096133e3
+};
+
+static double fm_exp2(double x)
+{
+    double   ipart, fpart, px, qx;
+    udi_t    epart;
+
+    ipart = floor(x+0.5);
+    fpart = x - ipart;
+    FM_DOUBLE_INIT_EXP(epart,ipart);
+
+    x = fpart*fpart;
+
+    px =        fm_exp2_p[0];
+    px = px*x + fm_exp2_p[1];
+    qx =    x + fm_exp2_q[0];
+    px = px*x + fm_exp2_p[2];
+    qx = qx*x + fm_exp2_q[1];
+
+    px = px * fpart;
+
+    x = 1.0 + 2.0*(px/(qx-px));
+    return epart.f*x;
+}
+
+double fm_exp(double x)
+{
+#if defined(__BYTE_ORDER__)
+#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
+    return fm_exp2(FM_DOUBLE_LOG2OFE * x);
+#endif
+#endif
+    return exp(x);
+}
+
+/*
+ * Local Variables:
+ * mode: c
+ * compile-command: "make -C .."
+ * c-basic-offset: 4
+ * fill-column: 76
+ * indent-tabs-mode: nil
+ * End:
+ */

src/USER-MEAMC/meam.h

@@ -0,0 +1,157 @@
+#include <stdlib.h>
+
+      #define maxelt 5
+      double fm_exp(double);
+	  
+	  typedef enum {FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2} lattice_t;
+
+typedef struct {
+  int dim1, dim2;
+  double *ptr;
+} allocatable_double_2d;
+    
+typedef struct {
+// cutforce = force cutoff
+// cutforcesq = force cutoff squared
+
+      double cutforce,cutforcesq;
+
+// Ec_meam = cohesive energy
+// re_meam = nearest-neighbor distance
+// Omega_meam = atomic volume
+// B_meam = bulk modulus
+// Z_meam = number of first neighbors for reference structure
+// ielt_meam = atomic number of element
+// A_meam = adjustable parameter
+// alpha_meam = sqrt(9*Omega*B/Ec)
+// rho0_meam = density scaling parameter
+// delta_meam = heat of formation for alloys
+// beta[0-3]_meam = electron density constants
+// t[0-3]_meam = coefficients on densities in Gamma computation
+// rho_ref_meam = background density for reference structure
+// ibar_meam(i) = selection parameter for Gamma function for elt i,
+// lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j)
+// neltypes = maximum number of element type defined
+// eltind = index number of pair (similar to Voigt notation; ij = ji)
+// phir = pair potential function array
+// phirar[1-6] = spline coeffs
+// attrac_meam = attraction parameter in Rose energy
+// repuls_meam = repulsion parameter in Rose energy
+// nn2_meam = 1 if second nearest neighbors are to be computed, else 0
+// zbl_meam = 1 if zbl potential for small r to be use, else 0
+// emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0
+// bkgd_dyn = 1 if reference densities follows Dynamo, else 0
+// Cmin_meam, Cmax_meam = min and max values in screening cutoff
+// rc_meam = cutoff distance for meam
+// delr_meam = cutoff region for meam
+// ebound_meam = factor giving maximum boundary of sceen fcn ellipse
+// augt1 = flag for whether t1 coefficient should be augmented
+// ialloy = flag for newer alloy formulation (as in dynamo code)
+// mix_ref_t = flag to recover "old" way of computing t in reference config
+// erose_form = selection parameter for form of E_rose function
+// gsmooth_factor = factor determining length of G smoothing region
+// vind[23]D = Voight notation index maps for 2 and 3D
+// v2D,v3D = array of factors to apply for Voight notation
+
+// nr,dr = pair function discretization parameters
+// nrar,rdrar = spline coeff array parameters
+
+      double Ec_meam[maxelt+1][maxelt+1],re_meam[maxelt+1][maxelt+1];
+      double Omega_meam[maxelt+1],Z_meam[maxelt+1];
+      double A_meam[maxelt+1],alpha_meam[maxelt+1][maxelt+1],rho0_meam[maxelt+1];
+      double delta_meam[maxelt+1][maxelt+1];
+      double beta0_meam[maxelt+1],beta1_meam[maxelt+1];
+      double beta2_meam[maxelt+1],beta3_meam[maxelt+1];
+      double t0_meam[maxelt+1],t1_meam[maxelt+1];
+      double t2_meam[maxelt+1],t3_meam[maxelt+1];
+      double rho_ref_meam[maxelt+1];
+      int ibar_meam[maxelt+1],ielt_meam[maxelt+1];
+	  lattice_t lattce_meam[maxelt+1][maxelt+1];
+      int nn2_meam[maxelt+1][maxelt+1];
+      int zbl_meam[maxelt+1][maxelt+1];
+      int eltind[maxelt+1][maxelt+1];
+      int neltypes;
+
+      allocatable_double_2d phir; // [:,:]
+
+      allocatable_double_2d phirar,phirar1,phirar2,phirar3,phirar4,phirar5,phirar6;  // [:,:]
+
+      double attrac_meam[maxelt+1][maxelt+1],repuls_meam[maxelt+1][maxelt+1];
+
+      double Cmin_meam[maxelt+1][maxelt+1][maxelt+1];
+      double Cmax_meam[maxelt+1][maxelt+1][maxelt+1];
+      double rc_meam,delr_meam,ebound_meam[maxelt+1][maxelt+1];
+      int augt1, ialloy, mix_ref_t, erose_form;
+      int emb_lin_neg, bkgd_dyn;
+      double  gsmooth_factor;
+      int vind2D[3+1][3+1],vind3D[3+1][3+1][3+1];
+      int v2D[6+1],v3D[10+1];
+
+      int nr,nrar;
+      double dr,rdrar;
+} meam_data_t;
+
+meam_data_t meam_data;
+
+// Functions we need for compat
+#ifndef max
+#define max(a,b) \
+   ({ __typeof__ (a) _a = (a); \
+       __typeof__ (b) _b = (b); \
+     _a > _b ? _a : _b; })
+#endif
+
+#ifndef min
+#define min(a,b) \
+   ({ __typeof__ (a) _a = (a); \
+       __typeof__ (b) _b = (b); \
+     _a < _b ? _a : _b; })
+#endif
+
+#define iszero(f) (fabs(f)<1e-20)
+
+#define setall2d(arr, v) {for(int __i=1;__i<=maxelt;__i++) for(int __j=1;__j<=maxelt;__j++) arr[__i][__j] = v;}
+#define setall3d(arr, v) {for(int __i=1;__i<=maxelt;__i++) for(int __j=1;__j<=maxelt;__j++) for(int __k=1;__k<=maxelt;__k++) arr[__i][__j][__k] = v;}
+
+/*
+Fortran Array Semantics in C.
+ - Stack-Allocated and global arrays are 1-based, declared as foo[N+1] and simply ignoring the first element
+    - Multi-Dimensional MUST be declared in reverse, so that the order when accessing is the same as in Fortran
+ - arrays that are passed externally via the meam_* functions must use the arr*v() functions below
+   (or be used with 0-based indexing)
+ - allocatable arrays (only global phir*) are actually a struct, and must be accessed with arr2()
+*/
+
+// we receive a pointer to the first element, and F dimensions is ptr(a,b,c)
+// we know c data structure is ptr[c][b][a]
+#define arrdim2v(ptr,a,b) \
+  const int DIM1__##ptr = a; const int DIM2__##ptr = b;\
+  (void)(DIM1__##ptr); (void)(DIM2__##ptr);
+#define arrdim3v(ptr,a,b,c) \
+  const int DIM1__##ptr = a; const int DIM2__##ptr = b; const int DIM3__##ptr = c;\
+  (void)(DIM1__##ptr); (void)(DIM2__##ptr; (void)(DIM3__##ptr);
+  
+
+// access data with same index as used in fortran (1-based)
+#define arr1v(ptr,i) \
+  ptr[i-1]
+#define arr2v(ptr,i,j) \
+  ptr[(DIM1__##ptr)*(j-1) + (i-1)]
+#define arr3v(ptr,i,j,k) \
+  ptr[(i-1) + (j-1)*(DIM1__##ptr) + (k-1)*(DIM1__##ptr)*(DIM2__##ptr)]
+
+// allocatable arrays via special type
+#define allocate_2d(arr,cols,rows) \
+  arr.dim1 = cols; arr.dim2=rows; arr.ptr=(double*)calloc((size_t)(cols)*(size_t)(rows),sizeof(double));
+#define allocated(a) \
+  (a.ptr!=NULL)
+#define deallocate(a) \
+  ({free(a.ptr); a.ptr=NULL;})
+// access data with same index as used in fortran (1-based)
+#define arr2(arr,i,j) \
+  arr.ptr[(arr.dim1)*(j-1) + (i-1)]
+
+
+
+
+

src/USER-MEAMC/meam_cleanup.c

@@ -0,0 +1,14 @@
+#include "meam.h"
+
+void meam_cleanup_(void) {
+
+    deallocate(meam_data.phirar6);
+    deallocate(meam_data.phirar5);
+    deallocate(meam_data.phirar4);
+    deallocate(meam_data.phirar3);
+    deallocate(meam_data.phirar2);
+    deallocate(meam_data.phirar1);
+    deallocate(meam_data.phirar);
+    deallocate(meam_data.phir);
+
+}

src/USER-MEAMC/meam_data.c

@@ -0,0 +1,3 @@
+#include "meam.h"
+
+meam_data_t meam_data = {};

src/USER-MEAMC/meam_dens_final.c

@@ -0,0 +1,262 @@
+#include "meam.h"
+#include <math.h>
+
+      void G_gam(double Gamma, int ibar, double gsmooth_factor, double *G, int *errorflag);
+      void dG_gam(double Gamma, int ibar, double gsmooth_factor,double *G, double *dG);
+
+// in meam_setup_done
+      void get_shpfcn(double *s /* s(3) */, lattice_t latt);
+
+// Extern "C" declaration has the form:
+//
+//  void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
+//                int *, int *, int *,
+//		 double *, double *, double *, double *, double *, double *,
+//		 double *, double *, double *, double *, double *, double *,
+//		 double *, double *, double *, double *, double *, int *);
+//
+// Call from pair_meam.cpp has the form:
+//
+//  meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
+//            &eng_vdwl,eatom,ntype,type,fmap,
+//	     &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
+//	     &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
+//	     dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
+//
+
+      void meam_dens_final_(int *nlocal, int *nmax,
+           int *eflag_either, int *eflag_global, int *eflag_atom, double *eng_vdwl, double *eatom,
+           int *ntype, int *type, int *fmap,
+           double *Arho1, double *Arho2, double *Arho2b, double *Arho3, double *Arho3b, double *t_ave, double *tsq_ave,
+           double *Gamma, double *dGamma1, double *dGamma2, double *dGamma3,
+           double *rho, double *rho0, double *rho1, double *rho2, double *rho3, double *fp, int *errorflag)
+      {
+        
+      arrdim2v(Arho1,3,*nmax)
+      arrdim2v(Arho2,6,*nmax)
+      arrdim2v(Arho3,10,*nmax)
+      arrdim2v(Arho3b,3,*nmax)
+      arrdim2v(t_ave,3,*nmax)
+      arrdim2v(tsq_ave,3,*nmax)
+
+      int i, elti;
+      int m;
+      double  rhob, G, dG, Gbar, dGbar, gam, shp[3+1], Z;
+      double  B, denom, rho_bkgd;
+
+//     Complete the calculation of density
+
+      for(i=1; i<=*nlocal; i++) {
+        elti = arr1v(fmap,arr1v(type,i));
+        if (elti > 0) {
+          arr1v(rho1,i) = 0.0;
+          arr1v(rho2,i) = -1.0/3.0*arr1v(Arho2b,i)*arr1v(Arho2b,i);
+          arr1v(rho3,i) = 0.0;
+          for(m=1; m<=3; m++) {
+            arr1v(rho1,i) = arr1v(rho1,i) + arr2v(Arho1,m,i)*arr2v(Arho1,m,i);
+            arr1v(rho3,i) = arr1v(rho3,i) - 3.0/5.0*arr2v(Arho3b,m,i)*arr2v(Arho3b,m,i);
+          }
+          for(m=1; m<=6; m++) {
+            arr1v(rho2,i) = arr1v(rho2,i) + meam_data.v2D[m]*arr2v(Arho2,m,i)*arr2v(Arho2,m,i);
+          }
+          for(m=1; m<=10; m++) {
+            arr1v(rho3,i) = arr1v(rho3,i) + meam_data.v3D[m]*arr2v(Arho3,m,i)*arr2v(Arho3,m,i);
+          }
+
+          if( arr1v(rho0,i)  >  0.0 ) {
+            if (meam_data.ialloy==1) {
+              arr2v(t_ave,1,i) = arr2v(t_ave,1,i)/arr2v(tsq_ave,1,i);
+              arr2v(t_ave,2,i) = arr2v(t_ave,2,i)/arr2v(tsq_ave,2,i);
+              arr2v(t_ave,3,i) = arr2v(t_ave,3,i)/arr2v(tsq_ave,3,i);
+            } else if (meam_data.ialloy==2) {
+              arr2v(t_ave,1,i) = meam_data.t1_meam[elti];
+              arr2v(t_ave,2,i) = meam_data.t2_meam[elti];
+              arr2v(t_ave,3,i) = meam_data.t3_meam[elti];
+            } else {
+              arr2v(t_ave,1,i) = arr2v(t_ave,1,i)/arr1v(rho0,i);
+              arr2v(t_ave,2,i) = arr2v(t_ave,2,i)/arr1v(rho0,i);
+              arr2v(t_ave,3,i) = arr2v(t_ave,3,i)/arr1v(rho0,i);
+            }
+          }
+
+          arr1v(Gamma,i) = arr2v(t_ave,1,i)*arr1v(rho1,i) + arr2v(t_ave,2,i)*arr1v(rho2,i) + arr2v(t_ave,3,i)*arr1v(rho3,i);
+
+          if( arr1v(rho0,i)  >  0.0 ) {
+            arr1v(Gamma,i) = arr1v(Gamma,i)/(arr1v(rho0,i)*arr1v(rho0,i));
+          }
+
+          Z = meam_data.Z_meam[elti];
+
+          G_gam(arr1v(Gamma,i),meam_data.ibar_meam[elti], meam_data.gsmooth_factor,&G,errorflag);
+          if (*errorflag!=0) return;
+          get_shpfcn(shp,meam_data.lattce_meam[elti][elti]);
+          if (meam_data.ibar_meam[elti]<=0) {
+            Gbar = 1.0;
+            dGbar = 0.0;
+          } else {
+            if (meam_data.mix_ref_t==1) {
+              gam = (arr2v(t_ave,1,i)*shp[1]+arr2v(t_ave,2,i)*shp[2] +arr2v(t_ave,3,i)*shp[3])/(Z*Z);
+            } else {
+              gam = (meam_data.t1_meam[elti]*shp[1]+meam_data.t2_meam[elti]*shp[2] +meam_data.t3_meam[elti]*shp[3])/(Z*Z);
+            }
+            G_gam(gam,meam_data.ibar_meam[elti],meam_data.gsmooth_factor,&Gbar,errorflag);
+          }
+          arr1v(rho,i) = arr1v(rho0,i) * G;
+
+          if (meam_data.mix_ref_t==1) {
+            if (meam_data.ibar_meam[elti]<=0) {
+              Gbar = 1.0;
+              dGbar = 0.0;
+            } else {
+              gam = (arr2v(t_ave,1,i)*shp[1]+arr2v(t_ave,2,i)*shp[2] +arr2v(t_ave,3,i)*shp[3])/(Z*Z);
+              dG_gam(gam,meam_data.ibar_meam[elti],meam_data.gsmooth_factor, &Gbar,&dGbar);
+            }
+            rho_bkgd = meam_data.rho0_meam[elti]*Z*Gbar;
+          } else {
+            if (meam_data.bkgd_dyn==1) {
+              rho_bkgd = meam_data.rho0_meam[elti]*Z;
+            } else {
+              rho_bkgd = meam_data.rho_ref_meam[elti];
+            }
+          }
+          rhob = arr1v(rho,i)/rho_bkgd;
+          denom = 1.0/rho_bkgd;
+
+          dG_gam(arr1v(Gamma,i),meam_data.ibar_meam[elti],meam_data.gsmooth_factor,&G,&dG);
+
+          arr1v(dGamma1,i) = (G - 2*dG*arr1v(Gamma,i))*denom;
+
+          if( !iszero(arr1v(rho0,i)) ) {
+            arr1v(dGamma2,i) = (dG/arr1v(rho0,i))*denom;
+          } else {
+            arr1v(dGamma2,i) = 0.0;
+          }
+
+//     dGamma3 is nonzero only if we are using the "mixed" rule for
+//     computing t in the reference system (which is not correct, but
+//     included for backward compatibility
+          if (meam_data.mix_ref_t==1) {
+            arr1v(dGamma3,i) = arr1v(rho0,i)*G*dGbar/(Gbar*Z*Z)*denom;
+          } else {
+            arr1v(dGamma3,i) = 0.0;
+          }
+
+          B = meam_data.A_meam[elti]*meam_data.Ec_meam[elti][elti];
+
+          if( !iszero(rhob) ) {
+            if (meam_data.emb_lin_neg==1 && rhob<=0) {
+              arr1v(fp,i) = -B;
+            } else {
+              arr1v(fp,i) = B*(log(rhob)+1.0);
+            }
+            if (*eflag_either!=0) {
+              if (*eflag_global!=0) {
+                if (meam_data.emb_lin_neg==1 && rhob<=0) {
+                  *eng_vdwl = *eng_vdwl - B*rhob;
+                } else {
+                  *eng_vdwl = *eng_vdwl + B*rhob*log(rhob);
+                }
+              }
+              if (*eflag_atom!=0) {
+                if (meam_data.emb_lin_neg==1 && rhob<=0) {
+                  arr1v(eatom,i) = arr1v(eatom,i) - B*rhob;
+                } else {
+                  arr1v(eatom,i) = arr1v(eatom,i) + B*rhob*log(rhob);
+                }
+              }
+            }
+          } else {
+            if (meam_data.emb_lin_neg==1) {
+              arr1v(fp,i) = -B;
+            } else {
+              arr1v(fp,i) = B;
+            }
+          }
+        }
+      }
+
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void G_gam(double Gamma, int ibar, double gsmooth_factor, double *G, int *errorflag)
+      {
+//     Compute G(Gamma) based on selection flag ibar:
+//   0 => G = sqrt(1+Gamma)
+//   1 => G = exp(Gamma/2)
+//   2 => not implemented
+//   3 => G = 2/(1+exp(-Gamma))
+//   4 => G = sqrt(1+Gamma)
+//  -5 => G = +-sqrt(abs(1+Gamma))
+      
+      double gsmooth_switchpoint;
+      if (ibar==0 || ibar==4) {
+        gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1);
+        if (Gamma < gsmooth_switchpoint) {
+//         e.g. gsmooth_factor is 99, {:
+//         gsmooth_switchpoint = -0.99
+//         G = 0.01*(-0.99/Gamma)**99
+          *G = 1/(gsmooth_factor+1) * pow((gsmooth_switchpoint/Gamma),gsmooth_factor);
+          *G = sqrt(*G);
+        } else {
+          *G = sqrt(1.0+Gamma);
+        }
+      } else if (ibar==1) {
+        *G = fm_exp(Gamma/2.0);
+      } else if (ibar==3) {
+        *G = 2.0/(1.0+exp(-Gamma));
+      } else if (ibar==-5) {
+        if ((1.0+Gamma)>=0) {
+           *G = sqrt(1.0+Gamma);
+        } else {
+           *G = -sqrt(-1.0-Gamma);
+        }
+      } else {
+         *errorflag = 1;
+      }
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void dG_gam(double Gamma, int ibar, double gsmooth_factor,double *G, double *dG)
+      {
+// Compute G(Gamma) and dG(gamma) based on selection flag ibar:
+//   0 => G = sqrt(1+Gamma)
+//   1 => G = fm_exp(Gamma/2)
+//   2 => not implemented
+//   3 => G = 2/(1+fm_exp(-Gamma))
+//   4 => G = sqrt(1+Gamma)
+//  -5 => G = +-sqrt(abs(1+Gamma))
+      double gsmooth_switchpoint;
+      
+      if (ibar==0 || ibar==4) {
+        gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1);
+        if (Gamma < gsmooth_switchpoint) {
+//         e.g. gsmooth_factor is 99, {:
+//         gsmooth_switchpoint = -0.99
+//         G = 0.01*(-0.99/Gamma)**99
+          *G = 1/(gsmooth_factor+1) * pow((gsmooth_switchpoint/Gamma), gsmooth_factor);
+          *G = sqrt(*G);
+          *dG = -gsmooth_factor* *G/(2.0*Gamma);
+        } else {
+          *G = sqrt(1.0+Gamma);
+          *dG = 1.0/(2.0* *G);
+        }
+      } else if (ibar==1) {
+        *G = fm_exp(Gamma/2.0);
+        *dG = *G/2.0;
+      } else if (ibar==3) {
+        *G = 2.0/(1.0+fm_exp(-Gamma));
+        *dG = *G*(2.0-*G)/2;
+      } else if (ibar==-5) {
+        if ((1.0+Gamma)>=0) {
+           *G = sqrt(1.0+Gamma);
+           *dG = 1.0/(2.0* *G);
+        } else {
+           *G = -sqrt(-1.0-Gamma);
+           *dG = -1.0/(2.0* *G);
+        }
+      }
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

src/USER-MEAMC/meam_dens_init.c

@@ -0,0 +1,504 @@
+#include "meam.h"
+#include <math.h>
+
+      void getscreen(int i, int nmax, double *scrfcn, double *dscrfcn, double *fcpair, double *x,
+          int numneigh, int *firstneigh,
+          int numneigh_full, int *firstneigh_full,
+          int ntype, int *type, int *fmap);
+      void calc_rho1(int i, int nmax, int ntype, int *type, int *fmap, double *x,
+          int numneigh, int *firstneigh,
+          double *scrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b,
+          double *arho3, double *arho3b, double *t_ave, double *tsq_ave);
+          
+      void screen(int i, int j, int nmax, double *x, double rijsq, double *sij, int numneigh_full, int *firstneigh_full, int ntype, int *type, int *fmap);
+      void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair);
+      void fcut(double xi, double * fc);
+      void dfcut(double xi, double *fc, double *dfc);
+      void dCfunc(double rij2,double rik2,double rjk2,double *dCikj);
+      void dCfunc2(double rij2,double rik2,double rjk2,double *dCikj1,double *dCikj2);
+
+//     Extern "C" declaration has the form:
+//
+//  void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *,
+//		 int *, int *, int *, int *,
+//		 double *, double *, double *, double *, double *, double *,
+//		 double *, double *, double *, double *, double *, int *);
+//
+//
+//     Call from pair_meam.cpp has the form:
+//
+//    meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0],
+//	       &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
+//	       &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
+//	       rho0,&arho1[0][0],&arho2[0][0],arho2b,
+//	       &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag);
+//
+
+      void meam_dens_init_(int *i, int *nmax,
+          int *ntype, int *type, int *fmap, double *x,
+          int *numneigh, int *firstneigh,
+          int *numneigh_full, int *firstneigh_full,
+          double *scrfcn, double *dscrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b,
+          double *arho3, double *arho3b, double *t_ave, double *tsq_ave, int *errorflag)
+      {
+      *errorflag = 0;
+
+//     Compute screening function and derivatives
+      getscreen(*i, *nmax, scrfcn, dscrfcn, fcpair, x,
+          *numneigh, firstneigh,
+          *numneigh_full, firstneigh_full,
+          *ntype, type, fmap);
+
+//     Calculate intermediate density terms to be communicated
+      calc_rho1(*i, *nmax, *ntype, type, fmap, x,
+          *numneigh, firstneigh,
+          scrfcn, fcpair, rho0, arho1, arho2, arho2b,
+          arho3, arho3b, t_ave, tsq_ave);
+
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void getscreen(int i, int nmax, double *scrfcn, double *dscrfcn, double *fcpair, double *x,
+          int numneigh, int *firstneigh,
+          int numneigh_full, int *firstneigh_full,
+          int ntype, int *type, int *fmap)
+      {
+        
+        arrdim2v(x,3,nmax)
+
+      int jn,j,kn,k;
+      int elti,eltj,eltk;
+      double xitmp,yitmp,zitmp,delxij,delyij,delzij,rij2,rij;
+      double xjtmp,yjtmp,zjtmp,delxik,delyik,delzik,rik2/*,rik*/;
+      double xktmp,yktmp,zktmp,delxjk,delyjk,delzjk,rjk2/*,rjk*/;
+      double xik,xjk,sij,fcij,sfcij,dfcij,sikj,dfikj,cikj;
+      double Cmin,Cmax,delc,/*ebound,*/rbound,a,coef1,coef2;
+      //double coef1a,coef1b,coef2a,coef2b;
+      double dCikj;
+      /*unused:double dC1a,dC1b,dC2a,dC2b;*/
+      double rnorm,fc,dfc,drinv;
+
+      drinv = 1.0/meam_data.delr_meam;
+      elti = arr1v(fmap,arr1v(type,i));
+
+      if (elti > 0) {
+
+        xitmp = arr2v(x,1,i);
+        yitmp = arr2v(x,2,i);
+        zitmp = arr2v(x,3,i);
+
+        for(jn=1; jn<=numneigh; jn++) {
+          j = arr1v(firstneigh,jn);
+
+          eltj = arr1v(fmap,arr1v(type,j));
+          if (eltj > 0) {
+
+//     First compute screening function itself, sij
+            xjtmp = arr2v(x,1,j);
+            yjtmp = arr2v(x,2,j);
+            zjtmp = arr2v(x,3,j);
+            delxij = xjtmp - xitmp;
+            delyij = yjtmp - yitmp;
+            delzij = zjtmp - zitmp;
+            rij2 = delxij*delxij + delyij*delyij + delzij*delzij;
+            rij = sqrt(rij2);
+            if (rij > meam_data.rc_meam) {
+              fcij = 0.0;
+              dfcij = 0.0;
+              sij = 0.0;
+            } else {
+              rnorm = (meam_data.rc_meam-rij)*drinv;
+              screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full, ntype, type, fmap);
+              dfcut(rnorm,&fc,&dfc);
+              fcij = fc;
+              dfcij = dfc*drinv;
+            }
+
+//     Now compute derivatives
+            arr1v(dscrfcn,jn) = 0.0;
+            sfcij = sij*fcij;
+            if (iszero(sfcij) || iszero(sfcij-1.0)) goto LABEL_100;
+            rbound = meam_data.ebound_meam[elti][eltj] * rij2;
+            for(kn=1; kn<=numneigh_full; kn++) {
+              k = arr1v(firstneigh_full,kn);
+              if (k==j) continue;
+              eltk = arr1v(fmap,arr1v(type,k));
+              if (eltk==0) continue;
+              xktmp = arr2v(x,1,k);
+              yktmp = arr2v(x,2,k);
+              zktmp = arr2v(x,3,k);
+              delxjk = xktmp - xjtmp;
+              delyjk = yktmp - yjtmp;
+              delzjk = zktmp - zjtmp;
+              rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk;
+              if (rjk2 > rbound) continue;
+              delxik = xktmp - xitmp;
+              delyik = yktmp - yitmp;
+              delzik = zktmp - zitmp;
+              rik2 = delxik*delxik + delyik*delyik + delzik*delzik;
+              if (rik2 > rbound) continue;
+              xik = rik2/rij2;
+              xjk = rjk2/rij2;
+              a =  1 - (xik-xjk)*(xik-xjk);
+//     if a < 0, then ellipse equation doesn't describe this case and
+//     atom k can't possibly screen i-j
+              if (a<=0.0) continue;
+              cikj = (2.0*(xik+xjk) + a - 2.0)/a;
+              Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
+              Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
+              if (cikj>=Cmax) {
+                continue;
+//     Note that cikj may be slightly negative (within numerical
+//     tolerance) if atoms are colinear, so don't reject that case here
+//     (other negative cikj cases were handled by the test on "a" above)
+//     Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
+              } else {
+                delc = Cmax - Cmin;
+                cikj = (cikj-Cmin)/delc;
+                dfcut(cikj,&sikj,&dfikj);
+                coef1 = dfikj/(delc*sikj);
+                dCfunc(rij2,rik2,rjk2,&dCikj);
+                arr1v(dscrfcn,jn) = arr1v(dscrfcn,jn) + coef1*dCikj;
+              }
+            }
+            coef1 = sfcij;
+            coef2 = sij*dfcij/rij;
+            arr1v(dscrfcn,jn) = arr1v(dscrfcn,jn)*coef1 - coef2;
+
+ LABEL_100:
+            arr1v(scrfcn,jn) = sij;
+            arr1v(fcpair,jn) = fcij;
+
+          }
+
+        }
+
+      }
+      }
+
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void calc_rho1(int i, int nmax, int ntype, int *type, int *fmap, double *x,
+          int numneigh, int *firstneigh,
+          double *scrfcn, double *fcpair, double *rho0, double *arho1, double *arho2, double *arho2b,
+          double *arho3, double *arho3b, double *t_ave, double *tsq_ave)
+      {
+        arrdim2v(x,3,nmax)
+        arrdim2v(arho1,3,nmax)
+        arrdim2v(arho2,6,nmax)
+        arrdim2v(arho3,10,nmax)
+        arrdim2v(arho3b,3,nmax)
+        arrdim2v(t_ave,3,nmax)
+        arrdim2v(tsq_ave,3,nmax)
+
+      int jn,j,m,n,p,elti,eltj;
+      int nv2,nv3;
+      double xtmp,ytmp,ztmp,delij[3+1],rij2,rij,sij;
+      double ai,aj,rhoa0j,rhoa1j,rhoa2j,rhoa3j,A1j,A2j,A3j;
+      //double G,Gbar,gam,shp[3+1];
+      double ro0i,ro0j;
+      double rhoa0i,rhoa1i,rhoa2i,rhoa3i,A1i,A2i,A3i;
+
+      elti = arr1v(fmap,arr1v(type,i));
+      xtmp = arr2v(x,1,i);
+      ytmp = arr2v(x,2,i);
+      ztmp = arr2v(x,3,i);
+      for(jn=1; jn<=numneigh; jn++) {
+        if (!iszero(arr1v(scrfcn,jn))) {
+          j = arr1v(firstneigh,jn);
+          sij = arr1v(scrfcn,jn)*arr1v(fcpair,jn);
+          delij[1] = arr2v(x,1,j) - xtmp;
+          delij[2] = arr2v(x,2,j) - ytmp;
+          delij[3] = arr2v(x,3,j) - ztmp;
+          rij2 = delij[1]*delij[1] + delij[2]*delij[2] + delij[3]*delij[3];
+          if (rij2 < meam_data.cutforcesq) {
+            eltj = arr1v(fmap,arr1v(type,j));
+            rij = sqrt(rij2);
+            ai = rij/meam_data.re_meam[elti][elti] - 1.0;
+            aj = rij/meam_data.re_meam[eltj][eltj] - 1.0;
+            ro0i = meam_data.rho0_meam[elti];
+            ro0j = meam_data.rho0_meam[eltj];
+            rhoa0j = ro0j*fm_exp(-meam_data.beta0_meam[eltj]*aj)*sij;
+            rhoa1j = ro0j*fm_exp(-meam_data.beta1_meam[eltj]*aj)*sij;
+            rhoa2j = ro0j*fm_exp(-meam_data.beta2_meam[eltj]*aj)*sij;
+            rhoa3j = ro0j*fm_exp(-meam_data.beta3_meam[eltj]*aj)*sij;
+            rhoa0i = ro0i*fm_exp(-meam_data.beta0_meam[elti]*ai)*sij;
+            rhoa1i = ro0i*fm_exp(-meam_data.beta1_meam[elti]*ai)*sij;
+            rhoa2i = ro0i*fm_exp(-meam_data.beta2_meam[elti]*ai)*sij;
+            rhoa3i = ro0i*fm_exp(-meam_data.beta3_meam[elti]*ai)*sij;
+            if (meam_data.ialloy==1) {
+              rhoa1j = rhoa1j * meam_data.t1_meam[eltj];
+              rhoa2j = rhoa2j * meam_data.t2_meam[eltj];
+              rhoa3j = rhoa3j * meam_data.t3_meam[eltj];
+              rhoa1i = rhoa1i * meam_data.t1_meam[elti];
+              rhoa2i = rhoa2i * meam_data.t2_meam[elti];
+              rhoa3i = rhoa3i * meam_data.t3_meam[elti];
+            }
+            arr1v(rho0,i) = arr1v(rho0,i) + rhoa0j;
+            arr1v(rho0,j) = arr1v(rho0,j) + rhoa0i;
+// For ialloy = 2, use single-element value (not average)
+            if (meam_data.ialloy!=2) {
+              arr2v(t_ave,1,i) = arr2v(t_ave,1,i) + meam_data.t1_meam[eltj]*rhoa0j;
+              arr2v(t_ave,2,i) = arr2v(t_ave,2,i) + meam_data.t2_meam[eltj]*rhoa0j;
+              arr2v(t_ave,3,i) = arr2v(t_ave,3,i) + meam_data.t3_meam[eltj]*rhoa0j;
+              arr2v(t_ave,1,j) = arr2v(t_ave,1,j) + meam_data.t1_meam[elti]*rhoa0i;
+              arr2v(t_ave,2,j) = arr2v(t_ave,2,j) + meam_data.t2_meam[elti]*rhoa0i;
+              arr2v(t_ave,3,j) = arr2v(t_ave,3,j) + meam_data.t3_meam[elti]*rhoa0i;
+            }
+            if (meam_data.ialloy==1) {
+              arr2v(tsq_ave,1,i) = arr2v(tsq_ave,1,i) + meam_data.t1_meam[eltj]*meam_data.t1_meam[eltj]*rhoa0j;
+              arr2v(tsq_ave,2,i) = arr2v(tsq_ave,2,i) + meam_data.t2_meam[eltj]*meam_data.t2_meam[eltj]*rhoa0j;
+              arr2v(tsq_ave,3,i) = arr2v(tsq_ave,3,i) + meam_data.t3_meam[eltj]*meam_data.t3_meam[eltj]*rhoa0j;
+              arr2v(tsq_ave,1,j) = arr2v(tsq_ave,1,j) + meam_data.t1_meam[elti]*meam_data.t1_meam[elti]*rhoa0i;
+              arr2v(tsq_ave,2,j) = arr2v(tsq_ave,2,j) + meam_data.t2_meam[elti]*meam_data.t2_meam[elti]*rhoa0i;
+              arr2v(tsq_ave,3,j) = arr2v(tsq_ave,3,j) + meam_data.t3_meam[elti]*meam_data.t3_meam[elti]*rhoa0i;
+            }
+            arr1v(arho2b,i) = arr1v(arho2b,i) + rhoa2j;
+            arr1v(arho2b,j) = arr1v(arho2b,j) + rhoa2i;
+
+            A1j = rhoa1j/rij;
+            A2j = rhoa2j/rij2;
+            A3j = rhoa3j/(rij2*rij);
+            A1i = rhoa1i/rij;
+            A2i = rhoa2i/rij2;
+            A3i = rhoa3i/(rij2*rij);
+            nv2 = 1;
+            nv3 = 1;
+            for(m=1; m<=3; m++) {
+              arr2v(arho1,m,i)  = arr2v(arho1,m,i) + A1j*delij[m];
+              arr2v(arho1,m,j)  = arr2v(arho1,m,j) - A1i*delij[m];
+              arr2v(arho3b,m,i) = arr2v(arho3b,m,i) + rhoa3j*delij[m]/rij;
+              arr2v(arho3b,m,j) = arr2v(arho3b,m,j) - rhoa3i*delij[m]/rij;
+              for(n=m; n<=3; n++) {
+                arr2v(arho2,nv2,i) = arr2v(arho2,nv2,i) + A2j*delij[m]*delij[n];
+                arr2v(arho2,nv2,j) = arr2v(arho2,nv2,j) + A2i*delij[m]*delij[n];
+                nv2 = nv2+1;
+                for(p=n; p<=3; p++) {
+                  arr2v(arho3,nv3,i) = arr2v(arho3,nv3,i) + A3j*delij[m]*delij[n]*delij[p];
+                  arr2v(arho3,nv3,j) = arr2v(arho3,nv3,j) - A3i*delij[m]*delij[n]*delij[p];
+                  nv3 = nv3+1;
+                }
+              }
+            }
+          }
+        }
+      }
+
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void screen(int i, int j, int nmax, double *x, double rijsq, double *sij, int numneigh_full, int *firstneigh_full, int ntype, int *type, int *fmap)
+//     Screening function
+//     Inputs:  i = atom 1 id (integer)
+//     j = atom 2 id (integer)
+//     rijsq = squared distance between i and j
+//     Outputs: sij = screening function
+      {
+        
+        arrdim2v(x,3,nmax)
+
+      int k,nk/*,m*/;
+      int elti,eltj,eltk;
+      double delxik,delyik,delzik;
+      double delxjk,delyjk,delzjk;
+      double riksq,rjksq,xik,xjk,cikj,a,delc,sikj/*,fcij,rij*/;
+      double Cmax,Cmin,rbound;
+
+      *sij = 1.0;
+      eltj = arr1v(fmap,arr1v(type,j));
+      elti = arr1v(fmap,arr1v(type,j));
+
+//     if rjksq > ebound*rijsq, atom k is definitely outside the ellipse
+      rbound = meam_data.ebound_meam[elti][eltj]*rijsq;
+
+      for(nk=1; nk<=numneigh_full; nk++) {
+        k = arr1v(firstneigh_full,nk);
+        eltk = arr1v(fmap,arr1v(type,k));
+        if (k==j) continue;
+        delxjk = arr2v(x,1,k) - arr2v(x,1,j);
+        delyjk = arr2v(x,2,k) - arr2v(x,2,j);
+        delzjk = arr2v(x,3,k) - arr2v(x,3,j);
+        rjksq = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk;
+        if (rjksq > rbound) continue;
+        delxik = arr2v(x,1,k) - arr2v(x,1,i);
+        delyik = arr2v(x,2,k) - arr2v(x,2,i);
+        delzik = arr2v(x,3,k) - arr2v(x,3,i);
+        riksq = delxik*delxik + delyik*delyik + delzik*delzik;
+        if (riksq > rbound) continue;
+        xik = riksq/rijsq;
+        xjk = rjksq/rijsq;
+        a = 1 - (xik-xjk)*(xik-xjk);
+//     if a < 0, then ellipse equation doesn't describe this case and
+//     atom k can't possibly screen i-j
+        if (a<=0.0) continue;
+        cikj = (2.0*(xik+xjk) + a - 2.0)/a;
+        Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
+        Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
+        if (cikj>=Cmax)
+          continue;
+//     note that cikj may be slightly negative (within numerical
+//     tolerance) if atoms are colinear, so don't reject that case here
+//     (other negative cikj cases were handled by the test on "a" above)
+        else if (cikj<=Cmin) {
+          *sij = 0.0;
+          break;
+        } else {
+          delc = Cmax - Cmin;
+          cikj = (cikj-Cmin)/delc;
+          fcut(cikj,&sikj);
+        }
+        *sij = *sij * sikj;
+      }
+
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair)
+      {
+//     Inputs: i,j,k = id's of 3 atom triplet
+//     jn = id of i-j pair
+//     rij2 = squared distance between i and j
+//     Outputs: dsij1 = deriv. of sij w.r.t. rik
+//     dsij2 = deriv. of sij w.r.t. rjk
+
+        arrdim2v(x,3,nmax)
+      int elti,eltj,eltk;
+      double rik2,rjk2;
+
+      double dxik,dyik,dzik;
+      double dxjk,dyjk,dzjk;
+      double rbound,delc,sij,xik,xjk,cikj,sikj,dfc,a;
+      double Cmax,Cmin,dCikj1,dCikj2;
+
+      sij = arr1v(scrfcn,jn)*arr1v(fcpair,jn);
+      elti = arr1v(fmap,arr1v(type,i));
+      eltj = arr1v(fmap,arr1v(type,j));
+      eltk = arr1v(fmap,arr1v(type,k));
+      Cmax = meam_data.Cmax_meam[elti][eltj][eltk];
+      Cmin = meam_data.Cmin_meam[elti][eltj][eltk];
+
+      *dsij1 = 0.0;
+      *dsij2 = 0.0;
+      if (!iszero(sij) && !iszero(sij-1.0)) {
+        rbound = rij2*meam_data.ebound_meam[elti][eltj];
+        delc = Cmax-Cmin;
+        dxjk = arr2v(x,1,k) - arr2v(x,1,j);
+        dyjk = arr2v(x,2,k) - arr2v(x,2,j);
+        dzjk = arr2v(x,3,k) - arr2v(x,3,j);
+        rjk2 = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk;
+        if (rjk2<=rbound) {
+          dxik = arr2v(x,1,k) - arr2v(x,1,i);
+          dyik = arr2v(x,2,k) - arr2v(x,2,i);
+          dzik = arr2v(x,3,k) - arr2v(x,3,i);
+          rik2 = dxik*dxik + dyik*dyik + dzik*dzik;
+          if (rik2<=rbound) {
+            xik = rik2/rij2;
+            xjk = rjk2/rij2;
+            a = 1 - (xik-xjk)*(xik-xjk);
+            if (!iszero(a)) {
+              cikj = (2.0*(xik+xjk) + a - 2.0)/a;
+              if (cikj>=Cmin && cikj<=Cmax) {
+                cikj = (cikj-Cmin)/delc;
+                dfcut(cikj,&sikj,&dfc);
+                dCfunc2(rij2,rik2,rjk2,&dCikj1,&dCikj2);
+                a = sij/delc*dfc/sikj;
+                *dsij1 = a*dCikj1;
+                *dsij2 = a*dCikj2;
+              }
+            }
+          }
+        }
+      }
+      }
+
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void fcut(double xi, double * fc)
+      {
+//     cutoff function
+      double a;
+      if (xi>=1.0)
+        *fc = 1.0;
+      else if (xi<=0.0)
+        *fc = 0.0;
+      else {
+        a = 1.0-xi;
+        a = a*a;
+        a = a*a;
+        a = 1.0-a;
+        *fc = a*a;
+//     fc = xi
+      }
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void dfcut(double xi, double *fc, double *dfc)
+      {
+//     cutoff function and its derivative
+      double a,a3,a4;
+      if (xi>=1.0) {
+        *fc = 1.0;
+        *dfc = 0.0;
+      } else if (xi<=0.0) {
+        *fc = 0.0;
+        *dfc = 0.0;
+      } else {
+        a = 1.0-xi;
+        a3 = a*a*a;
+        a4 = a*a3;
+        *fc = pow((1.0-a4), 2);
+        *dfc = 8*(1.0-a4)*a3;
+//     fc = xi
+//     dfc = 1.d0
+      }
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void dCfunc(double rij2,double rik2,double rjk2,double *dCikj)
+      {
+//     Inputs: rij,rij2,rik2,rjk2
+//     Outputs: dCikj = derivative of Cikj w.r.t. rij
+      double rij4,a,b,denom;
+
+      rij4 = rij2*rij2;
+      a = rik2-rjk2;
+      b = rik2+rjk2;
+      denom = rij4 - a*a;
+      denom = denom*denom;
+      *dCikj = -4*(-2*rij2*a*a + rij4*b + a*a*b)/denom;
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      void dCfunc2(double rij2,double rik2,double rjk2,double *dCikj1,double *dCikj2)
+      {
+//     Inputs: rij,rij2,rik2,rjk2
+//     Outputs: dCikj1 = derivative of Cikj w.r.t. rik
+//     dCikj2 = derivative of Cikj w.r.t. rjk
+      double rij4,rik4,rjk4,a,b,denom;
+
+      rij4 = rij2*rij2;
+      rik4 = rik2*rik2;
+      rjk4 = rjk2*rjk2;
+      a = rik2-rjk2;
+      b = rik2+rjk2;
+      denom = rij4 - a*a;
+      denom = denom*denom;
+      *dCikj1 = 4*rij2*(rij4 + rik4 + 2*rik2*rjk2 - 3*rjk4 - 2*rij2*a)/ denom;
+      *dCikj2 = 4*rij2*(rij4 - 3*rik4 + 2*rik2*rjk2 + rjk4 + 2*rij2*a)/ denom;
+      
+      (void)(b);
+      }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+
+
+

src/USER-MEAMC/meam_force.c

@@ -0,0 +1,551 @@
+#include "meam.h"
+#include <math.h>
+
+// in meam_setup_done
+void get_shpfcn(double *s /* s(3) */, lattice_t latt);
+
+// in meam_dens_init
+void dsij(int i,int j,int k,int jn,int nmax,int numneigh,double rij2,double *dsij1,double *dsij2, int ntype, int *type, int *fmap,double *x,double *scrfcn, double*fcpair);
+
+
+//     Extern "C" declaration has the form:
+//
+//  void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *,
+//		 int *, int *, int *, int *, double *, double *,
+//		 double *, double *, double *, double *, double *, double *,
+//		 double *, double *, double *, double *, double *, double *,
+//		 double *, double *, double *, double *, double *, double *, int *);
+//
+//     Call from pair_meam.cpp has the form:
+//
+//    meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom,
+//              &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
+//	       &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
+//	       &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
+//	       dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
+//	       &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
+//	       &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag);
+//
+
+      void meam_force_(int *iptr, int *nmax,
+          int *eflag_either, int *eflag_global, int *eflag_atom, int *vflag_atom,
+          double *eng_vdwl, double *eatom, int *ntype, int *type, int *fmap, double *x,
+          int *numneigh, int *firstneigh, int *numneigh_full, int *firstneigh_full,
+          double *scrfcn, double *dscrfcn, double *fcpair,
+          double *dGamma1, double *dGamma2, double *dGamma3, double *rho0, double *rho1, double *rho2, double *rho3, double *fp,
+          double *Arho1, double *Arho2, double *Arho2b, double *Arho3, double *Arho3b, double *t_ave, double *tsq_ave, double *f,
+          double *vatom, int *errorflag)
+      {
+        arrdim2v(x,3,*nmax)
+        arrdim2v(Arho1,3,*nmax)
+        arrdim2v(Arho2,6,*nmax)
+        arrdim2v(Arho3,10,*nmax)
+        arrdim2v(Arho3b,3,*nmax)
+        arrdim2v(t_ave,3,*nmax)
+        arrdim2v(tsq_ave,3,*nmax)
+        arrdim2v(f,3,*nmax)
+        arrdim2v(vatom,6,*nmax)
+
+      int i,j,jn,k,kn,kk,m,n,p,q;
+      int nv2,nv3,elti,eltj,eltk,ind;
+      double xitmp,yitmp,zitmp,delij[3+1],rij2,rij,rij3;
+      double delik[3+1],deljk[3+1],v[6+1],fi[3+1],fj[3+1];
+      double third,sixth;
+      double pp,dUdrij,dUdsij,dUdrijm[3+1],force,forcem;
+      double r,recip,phi,phip;
+      double sij;
+      double a1,a1i,a1j,a2,a2i,a2j;
+      double a3i,a3j;
+      double shpi[3+1],shpj[3+1];
+      double ai,aj,ro0i,ro0j,invrei,invrej;
+      double rhoa0j,drhoa0j,rhoa0i,drhoa0i;
+      double rhoa1j,drhoa1j,rhoa1i,drhoa1i;
+      double rhoa2j,drhoa2j,rhoa2i,drhoa2i;
+      double a3,a3a,rhoa3j,drhoa3j,rhoa3i,drhoa3i;
+      double drho0dr1,drho0dr2,drho0ds1,drho0ds2;
+      double drho1dr1,drho1dr2,drho1ds1,drho1ds2;
+      double drho1drm1[3+1],drho1drm2[3+1];
+      double drho2dr1,drho2dr2,drho2ds1,drho2ds2;
+      double drho2drm1[3+1],drho2drm2[3+1];
+      double drho3dr1,drho3dr2,drho3ds1,drho3ds2;
+      double drho3drm1[3+1],drho3drm2[3+1];
+      double dt1dr1,dt1dr2,dt1ds1,dt1ds2;
+      double dt2dr1,dt2dr2,dt2ds1,dt2ds2;
+      double dt3dr1,dt3dr2,dt3ds1,dt3ds2;
+      double drhodr1,drhodr2,drhods1,drhods2,drhodrm1[3+1],drhodrm2[3+1];
+      double arg;
+      double arg1i1,arg1j1,arg1i2,arg1j2,arg1i3,arg1j3,arg3i3,arg3j3;
+      double dsij1,dsij2,force1,force2;
+      double t1i,t2i,t3i,t1j,t2j,t3j;
+
+      *errorflag = 0;
+      third = 1.0/3.0;
+      sixth = 1.0/6.0;
+      
+      //: aliased
+      i = *iptr;
+
+//     Compute forces atom i
+
+      elti = arr1v(fmap,arr1v(type,i));
+
+      if (elti > 0) {
+        xitmp = arr2v(x,1,i);
+        yitmp = arr2v(x,2,i);
+        zitmp = arr2v(x,3,i);
+
+//     Treat each pair
+        for(jn=1; jn<=*numneigh; jn++) {
+          j = arr1v(firstneigh,jn);
+          eltj = arr1v(fmap,arr1v(type,j));
+
+          if (!iszero(arr1v(scrfcn,jn)) && eltj > 0) {
+
+            sij = arr1v(scrfcn,jn)*arr1v(fcpair,jn);
+            delij[1] = arr2v(x,1,j) - xitmp;
+            delij[2] = arr2v(x,2,j) - yitmp;
+            delij[3] = arr2v(x,3,j) - zitmp;
+            rij2 = delij[1]*delij[1] + delij[2]*delij[2] + delij[3]*delij[3];
+            if (rij2 < meam_data.cutforcesq) {
+              rij = sqrt(rij2);
+              r = rij;
+
+//     Compute phi and phip
+              ind = meam_data.eltind[elti][eltj];
+              pp = rij*meam_data.rdrar + 1.0;
+              kk = (int)pp;
+              kk = min(kk,meam_data.nrar-1);
+              pp = pp - kk;
+              pp = min(pp,1.0);
+              phi = ((arr2(meam_data.phirar3,kk,ind)*pp + arr2(meam_data.phirar2,kk,ind))*pp + arr2(meam_data.phirar1,kk,ind))*pp + arr2(meam_data.phirar,kk,ind);
+              phip = (arr2(meam_data.phirar6,kk,ind)*pp + arr2(meam_data.phirar5,kk,ind))*pp + arr2(meam_data.phirar4,kk,ind);
+              recip = 1.0/r;
+
+              if (*eflag_either!=0) {
+                if (*eflag_global!=0) *eng_vdwl = *eng_vdwl + phi*sij;
+                if (*eflag_atom!=0) {
+                  arr1v(eatom,i) = arr1v(eatom,i) + 0.5*phi*sij;
+                  arr1v(eatom,j) = arr1v(eatom,j) + 0.5*phi*sij;
+                }
+              }
+
+//     write(1,*) "force_meamf: phi: ",phi
+//     write(1,*) "force_meamf: phip: ",phip
+
+//     Compute pair densities and derivatives
+              invrei = 1.0/meam_data.re_meam[elti][elti];
+              ai = rij*invrei - 1.0;
+              ro0i = meam_data.rho0_meam[elti];
+              rhoa0i = ro0i*fm_exp(-meam_data.beta0_meam[elti]*ai);
+              drhoa0i = -meam_data.beta0_meam[elti]*invrei*rhoa0i;
+              rhoa1i = ro0i*fm_exp(-meam_data.beta1_meam[elti]*ai);
+              drhoa1i = -meam_data.beta1_meam[elti]*invrei*rhoa1i;
+              rhoa2i = ro0i*fm_exp(-meam_data.beta2_meam[elti]*ai);
+              drhoa2i = -meam_data.beta2_meam[elti]*invrei*rhoa2i;
+              rhoa3i = ro0i*fm_exp(-meam_data.beta3_meam[elti]*ai);
+              drhoa3i = -meam_data.beta3_meam[elti]*invrei*rhoa3i;
+
+              if (elti!=eltj) {
+                invrej = 1.0/meam_data.re_meam[eltj][eltj];
+                aj = rij*invrej - 1.0;
+                ro0j = meam_data.rho0_meam[eltj];
+                rhoa0j = ro0j*fm_exp(-meam_data.beta0_meam[eltj]*aj);
+                drhoa0j = -meam_data.beta0_meam[eltj]*invrej*rhoa0j;
+                rhoa1j = ro0j*fm_exp(-meam_data.beta1_meam[eltj]*aj);
+                drhoa1j = -meam_data.beta1_meam[eltj]*invrej*rhoa1j;
+                rhoa2j = ro0j*fm_exp(-meam_data.beta2_meam[eltj]*aj);
+                drhoa2j = -meam_data.beta2_meam[eltj]*invrej*rhoa2j;
+                rhoa3j = ro0j*fm_exp(-meam_data.beta3_meam[eltj]*aj);
+                drhoa3j = -meam_data.beta3_meam[eltj]*invrej*rhoa3j;
+              } else {
+                rhoa0j = rhoa0i;
+                drhoa0j = drhoa0i;
+                rhoa1j = rhoa1i;
+                drhoa1j = drhoa1i;
+                rhoa2j = rhoa2i;
+                drhoa2j = drhoa2i;
+                rhoa3j = rhoa3i;
+                drhoa3j = drhoa3i;
+              }
+
+              if (meam_data.ialloy==1) {
+                rhoa1j = rhoa1j * meam_data.t1_meam[eltj];
+                rhoa2j = rhoa2j * meam_data.t2_meam[eltj];
+                rhoa3j = rhoa3j * meam_data.t3_meam[eltj];
+                rhoa1i = rhoa1i * meam_data.t1_meam[elti];
+                rhoa2i = rhoa2i * meam_data.t2_meam[elti];
+                rhoa3i = rhoa3i * meam_data.t3_meam[elti];
+                drhoa1j = drhoa1j * meam_data.t1_meam[eltj];
+                drhoa2j = drhoa2j * meam_data.t2_meam[eltj];
+                drhoa3j = drhoa3j * meam_data.t3_meam[eltj];
+                drhoa1i = drhoa1i * meam_data.t1_meam[elti];
+                drhoa2i = drhoa2i * meam_data.t2_meam[elti];
+                drhoa3i = drhoa3i * meam_data.t3_meam[elti];
+              }
+
+              nv2 = 1;
+              nv3 = 1;
+              arg1i1 = 0.0;
+              arg1j1 = 0.0;
+              arg1i2 = 0.0;
+              arg1j2 = 0.0;
+              arg1i3 = 0.0;
+              arg1j3 = 0.0;
+              arg3i3 = 0.0;
+              arg3j3 = 0.0;
+              for(n=1; n<=3; n++) {
+                for(p=n; p<=3; p++) {
+                  for(q=p; q<=3; q++) {
+                    arg = delij[n]*delij[p]*delij[q]*meam_data.v3D[nv3];
+                    arg1i3 = arg1i3 + arr2v(Arho3,nv3,i)*arg;
+                    arg1j3 = arg1j3 - arr2v(Arho3,nv3,j)*arg;
+                    nv3 = nv3+1;
+                  }
+                  arg = delij[n]*delij[p]*meam_data.v2D[nv2];
+                  arg1i2 = arg1i2 + arr2v(Arho2,nv2,i)*arg;
+                  arg1j2 = arg1j2 + arr2v(Arho2,nv2,j)*arg;
+                  nv2 = nv2+1;
+                }
+                arg1i1 = arg1i1 + arr2v(Arho1,n,i)*delij[n];
+                arg1j1 = arg1j1 - arr2v(Arho1,n,j)*delij[n];
+                arg3i3 = arg3i3 + arr2v(Arho3b,n,i)*delij[n];
+                arg3j3 = arg3j3 - arr2v(Arho3b,n,j)*delij[n];
+              }
+
+//     rho0 terms
+              drho0dr1 = drhoa0j * sij;
+              drho0dr2 = drhoa0i * sij;
+
+//     rho1 terms
+              a1 = 2*sij/rij;
+              drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1;
+              drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1;
+              a1 = 2.0*sij/rij;
+              for(m=1; m<=3; m++) {
+                drho1drm1[m] =  a1*rhoa1j*arr2v(Arho1,m,i);
+                drho1drm2[m] = -a1*rhoa1i*arr2v(Arho1,m,j);
+              }
+
+//     rho2 terms
+              a2 = 2*sij/rij2;
+              drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2            - 2.0/3.0*arr1v(Arho2b,i)*drhoa2j*sij;
+              drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2            - 2.0/3.0*arr1v(Arho2b,j)*drhoa2i*sij;
+              a2 = 4*sij/rij2;
+              for(m=1; m<=3; m++) {
+                drho2drm1[m] = 0.0;
+                drho2drm2[m] = 0.0;
+                for(n=1; n<=3; n++) {
+                  drho2drm1[m] = drho2drm1[m]               + arr2v(Arho2,meam_data.vind2D[m][n],i)*delij[n];
+                  drho2drm2[m] = drho2drm2[m]               - arr2v(Arho2,meam_data.vind2D[m][n],j)*delij[n];
+                }
+                drho2drm1[m] =  a2*rhoa2j*drho2drm1[m];
+                drho2drm2[m] = -a2*rhoa2i*drho2drm2[m];
+              }
+
+//     rho3 terms
+              rij3 = rij*rij2;
+              a3 = 2*sij/rij3;
+              a3a = 6.0/5.0*sij/rij;
+              drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3  - a3a*(drhoa3j - rhoa3j/rij)*arg3i3;
+              drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3  - a3a*(drhoa3i - rhoa3i/rij)*arg3j3;
+              a3 = 6*sij/rij3;
+              a3a = 6*sij/(5*rij);
+              for(m=1; m<=3; m++) {
+                drho3drm1[m] = 0.0;
+                drho3drm2[m] = 0.0;
+                nv2 = 1;
+                for(n=1; n<=3; n++) {
+                  for(p=n; p<=3; p++) {
+                    arg = delij[n]*delij[p]*meam_data.v2D[nv2];
+                    drho3drm1[m] = drho3drm1[m]     + arr2v(Arho3,meam_data.vind3D[m][n][p],i)*arg;
+                    drho3drm2[m] = drho3drm2[m]     + arr2v(Arho3,meam_data.vind3D[m][n][p],j)*arg;
+                    nv2 = nv2 + 1;
+                  }
+                }
+                drho3drm1[m] = ( a3*drho3drm1[m] - a3a*arr2v(Arho3b,m,i))  *rhoa3j;
+                drho3drm2[m] = (-a3*drho3drm2[m] + a3a*arr2v(Arho3b,m,j))  *rhoa3i;
+              }
+
+//     Compute derivatives of weighting functions t wrt rij
+              t1i = arr2v(t_ave,1,i);
+              t2i = arr2v(t_ave,2,i);
+              t3i = arr2v(t_ave,3,i);
+              t1j = arr2v(t_ave,1,j);
+              t2j = arr2v(t_ave,2,j);
+              t3j = arr2v(t_ave,3,j);
+
+              if (meam_data.ialloy==1) {
+
+                a1i = 0.0;
+                a1j = 0.0;
+                a2i = 0.0;
+                a2j = 0.0;
+                a3i = 0.0;
+                a3j = 0.0;
+                if ( !iszero(arr2v(tsq_ave,1,i)) )  a1i = drhoa0j*sij/arr2v(tsq_ave,1,i);
+                if ( !iszero(arr2v(tsq_ave,1,j)) )  a1j = drhoa0i*sij/arr2v(tsq_ave,1,j);
+                if ( !iszero(arr2v(tsq_ave,2,i)) )  a2i = drhoa0j*sij/arr2v(tsq_ave,2,i);
+                if ( !iszero(arr2v(tsq_ave,2,j)) )  a2j = drhoa0i*sij/arr2v(tsq_ave,2,j);
+                if ( !iszero(arr2v(tsq_ave,3,i)) )  a3i = drhoa0j*sij/arr2v(tsq_ave,3,i);
+                if ( !iszero(arr2v(tsq_ave,3,j)) )  a3j = drhoa0i*sij/arr2v(tsq_ave,3,j);
+
+                dt1dr1 = a1i*(meam_data.t1_meam[eltj]-t1i*pow(meam_data.t1_meam[eltj],2));
+                dt1dr2 = a1j*(meam_data.t1_meam[elti]-t1j*pow(meam_data.t1_meam[elti],2));
+                dt2dr1 = a2i*(meam_data.t2_meam[eltj]-t2i*pow(meam_data.t2_meam[eltj],2));
+                dt2dr2 = a2j*(meam_data.t2_meam[elti]-t2j*pow(meam_data.t2_meam[elti],2));
+                dt3dr1 = a3i*(meam_data.t3_meam[eltj]-t3i*pow(meam_data.t3_meam[eltj],2));
+                dt3dr2 = a3j*(meam_data.t3_meam[elti]-t3j*pow(meam_data.t3_meam[elti],2));
+
+              } else if (meam_data.ialloy==2) {
+
+                dt1dr1 = 0.0;
+                dt1dr2 = 0.0;
+                dt2dr1 = 0.0;
+                dt2dr2 = 0.0;
+                dt3dr1 = 0.0;
+                dt3dr2 = 0.0;
+
+              } else {
+
+                ai = 0.0;
+                if( !iszero(arr1v(rho0,i)) )
+                  ai = drhoa0j*sij/arr1v(rho0,i);
+                aj = 0.0;
+                if( !iszero(arr1v(rho0,j)) )
+                  aj = drhoa0i*sij/arr1v(rho0,j);
+
+                dt1dr1 = ai*(meam_data.t1_meam[elti]-t1i);
+                dt1dr2 = aj*(meam_data.t1_meam[elti]-t1j);
+                dt2dr1 = ai*(meam_data.t2_meam[elti]-t2i);
+                dt2dr2 = aj*(meam_data.t2_meam[elti]-t2j);
+                dt3dr1 = ai*(meam_data.t3_meam[elti]-t3i);
+                dt3dr2 = aj*(meam_data.t3_meam[elti]-t3j);
+
+              }
+
+//     Compute derivatives of total density wrt rij, sij and rij(3)
+              get_shpfcn(shpi,meam_data.lattce_meam[elti][elti]);
+              get_shpfcn(shpj,meam_data.lattce_meam[eltj][eltj]);
+              drhodr1 = arr1v(dGamma1,i)*drho0dr1
+                  + arr1v(dGamma2,i)*
+                  (dt1dr1*arr1v(rho1,i)+t1i*drho1dr1
+                  + dt2dr1*arr1v(rho2,i)+t2i*drho2dr1
+                  + dt3dr1*arr1v(rho3,i)+t3i*drho3dr1)
+                  - arr1v(dGamma3,i)*
+                  (shpi[1]*dt1dr1+shpi[2]*dt2dr1+shpi[3]*dt3dr1);
+              drhodr2 = arr1v(dGamma1,j)*drho0dr2
+                  + arr1v(dGamma2,j)*
+                  (dt1dr2*arr1v(rho1,j)+t1j*drho1dr2
+                  + dt2dr2*arr1v(rho2,j)+t2j*drho2dr2
+                  + dt3dr2*arr1v(rho3,j)+t3j*drho3dr2)
+                  - arr1v(dGamma3,j)*
+                  (shpj[1]*dt1dr2+shpj[2]*dt2dr2+shpj[3]*dt3dr2);
+              for(m=1; m<=3; m++) {
+                drhodrm1[m] = 0.0;
+                drhodrm2[m] = 0.0;
+                drhodrm1[m] = arr1v(dGamma2,i)*
+                     (t1i*drho1drm1[m]
+                    + t2i*drho2drm1[m]
+                    + t3i*drho3drm1[m]);
+                drhodrm2[m] = arr1v(dGamma2,j)*
+                     (t1j*drho1drm2[m]
+                    + t2j*drho2drm2[m]
+                    + t3j*drho3drm2[m]);
+              }
+
+//     Compute derivatives wrt sij, but only if necessary
+              if ( !iszero(arr1v(dscrfcn,jn))) {
+                drho0ds1 = rhoa0j;
+                drho0ds2 = rhoa0i;
+                a1 = 2.0/rij;
+                drho1ds1 = a1*rhoa1j*arg1i1;
+                drho1ds2 = a1*rhoa1i*arg1j1;
+                a2 = 2.0/rij2;
+                drho2ds1 = a2*rhoa2j*arg1i2  - 2.0/3.0*arr1v(Arho2b,i)*rhoa2j;
+                drho2ds2 = a2*rhoa2i*arg1j2  - 2.0/3.0*arr1v(Arho2b,j)*rhoa2i;
+                a3 = 2.0/rij3;
+                a3a = 6.0/(5.0*rij);
+                drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3;
+                drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3;
+
+                if (meam_data.ialloy==1) {
+
+                  a1i = 0.0;
+                  a1j = 0.0;
+                  a2i = 0.0;
+                  a2j = 0.0;
+                  a3i = 0.0;
+                  a3j = 0.0;
+                  if ( !iszero(arr2v(tsq_ave,1,i)) ) a1i = rhoa0j/arr2v(tsq_ave,1,i);
+                  if ( !iszero(arr2v(tsq_ave,1,j)) ) a1j = rhoa0i/arr2v(tsq_ave,1,j);
+                  if ( !iszero(arr2v(tsq_ave,2,i)) ) a2i = rhoa0j/arr2v(tsq_ave,2,i);
+                  if ( !iszero(arr2v(tsq_ave,2,j)) ) a2j = rhoa0i/arr2v(tsq_ave,2,j);
+                  if ( !iszero(arr2v(tsq_ave,3,i)) ) a3i = rhoa0j/arr2v(tsq_ave,3,i);
+                  if ( !iszero(arr2v(tsq_ave,3,j)) ) a3j = rhoa0i/arr2v(tsq_ave,3,j);
+
+                  dt1ds1 = a1i*(meam_data.t1_meam[eltj]-t1i*pow(meam_data.t1_meam[eltj],2));
+                  dt1ds2 = a1j*(meam_data.t1_meam[elti]-t1j*pow(meam_data.t1_meam[elti],2));
+                  dt2ds1 = a2i*(meam_data.t2_meam[eltj]-t2i*pow(meam_data.t2_meam[eltj],2));
+                  dt2ds2 = a2j*(meam_data.t2_meam[elti]-t2j*pow(meam_data.t2_meam[elti],2));
+                  dt3ds1 = a3i*(meam_data.t3_meam[eltj]-t3i*pow(meam_data.t3_meam[eltj],2));
+                  dt3ds2 = a3j*(meam_data.t3_meam[elti]-t3j*pow(meam_data.t3_meam[elti],2));
+
+                } else if (meam_data.ialloy==2) {
+
+                  dt1ds1 = 0.0;
+                  dt1ds2 = 0.0;
+                  dt2ds1 = 0.0;
+                  dt2ds2 = 0.0;
+                  dt3ds1 = 0.0;
+                  dt3ds2 = 0.0;
+
+                } else {
+
+                  ai = 0.0;
+                  if( !iszero(arr1v(rho0,i)) )
+                    ai = rhoa0j/arr1v(rho0,i);
+                  aj = 0.0;
+                  if( !iszero(arr1v(rho0,j)) )
+                    aj = rhoa0i/arr1v(rho0,j);
+
+                  dt1ds1 = ai*(meam_data.t1_meam[eltj]-t1i);
+                  dt1ds2 = aj*(meam_data.t1_meam[elti]-t1j);
+                  dt2ds1 = ai*(meam_data.t2_meam[eltj]-t2i);
+                  dt2ds2 = aj*(meam_data.t2_meam[elti]-t2j);
+                  dt3ds1 = ai*(meam_data.t3_meam[eltj]-t3i);
+                  dt3ds2 = aj*(meam_data.t3_meam[elti]-t3j);
+
+                }
+
+                drhods1 = arr1v(dGamma1,i)*drho0ds1
+                    + arr1v(dGamma2,i)*
+                    (dt1ds1*arr1v(rho1,i)+t1i*drho1ds1
+                    + dt2ds1*arr1v(rho2,i)+t2i*drho2ds1
+                    + dt3ds1*arr1v(rho3,i)+t3i*drho3ds1)
+                    - arr1v(dGamma3,i)*
+                    (shpi[1]*dt1ds1+shpi[2]*dt2ds1+shpi[3]*dt3ds1);
+                drhods2 = arr1v(dGamma1,j)*drho0ds2
+                    + arr1v(dGamma2,j)*
+                    (dt1ds2*arr1v(rho1,j)+t1j*drho1ds2
+                    + dt2ds2*arr1v(rho2,j)+t2j*drho2ds2
+                    + dt3ds2*arr1v(rho3,j)+t3j*drho3ds2)
+                    - arr1v(dGamma3,j)*
+                    (shpj[1]*dt1ds2+shpj[2]*dt2ds2+shpj[3]*dt3ds2);
+              }
+
+//     Compute derivatives of energy wrt rij, sij and rij[3]
+              dUdrij = phip*sij           + arr1v(fp,i)*drhodr1 + arr1v(fp,j)*drhodr2;
+              dUdsij = 0.0;
+              if ( !iszero(arr1v(dscrfcn,jn)) ) {
+                dUdsij = phi              + arr1v(fp,i)*drhods1 + arr1v(fp,j)*drhods2;
+              }
+              for(m=1; m<=3; m++) {
+                dUdrijm[m] = arr1v(fp,i)*drhodrm1[m] + arr1v(fp,j)*drhodrm2[m];
+              }
+
+//     Add the part of the force due to dUdrij and dUdsij
+
+              force = dUdrij*recip + dUdsij*arr1v(dscrfcn,jn);
+              for(m=1; m<=3; m++) {
+                forcem = delij[m]*force + dUdrijm[m];
+                arr2v(f,m,i) = arr2v(f,m,i) + forcem;
+                arr2v(f,m,j) = arr2v(f,m,j) - forcem;
+              }
+
+//     Tabulate per-atom virial as symmetrized stress tensor
+
+              if (*vflag_atom!=0) {
+                fi[1] = delij[1]*force + dUdrijm[1];
+                fi[2] = delij[2]*force + dUdrijm[2];
+                fi[3] = delij[3]*force + dUdrijm[3];
+                v[1] = -0.5 * (delij[1] * fi[1]);
+                v[2] = -0.5 * (delij[2] * fi[2]);
+                v[3] = -0.5 * (delij[3] * fi[3]);
+                v[4] = -0.25 * (delij[1]*fi[2] + delij[2]*fi[1]);
+                v[5] = -0.25 * (delij[1]*fi[3] + delij[3]*fi[1]);
+                v[6] = -0.25 * (delij[2]*fi[3] + delij[3]*fi[2]);
+
+                arr2v(vatom,1,i) = arr2v(vatom,1,i) + v[1];
+                arr2v(vatom,2,i) = arr2v(vatom,2,i) + v[2];
+                arr2v(vatom,3,i) = arr2v(vatom,3,i) + v[3];
+                arr2v(vatom,4,i) = arr2v(vatom,4,i) + v[4];
+                arr2v(vatom,5,i) = arr2v(vatom,5,i) + v[5];
+                arr2v(vatom,6,i) = arr2v(vatom,6,i) + v[6];
+                arr2v(vatom,1,j) = arr2v(vatom,1,j) + v[1];
+                arr2v(vatom,2,j) = arr2v(vatom,2,j) + v[2];
+                arr2v(vatom,3,j) = arr2v(vatom,3,j) + v[3];
+                arr2v(vatom,4,j) = arr2v(vatom,4,j) + v[4];
+                arr2v(vatom,5,j) = arr2v(vatom,5,j) + v[5];
+                arr2v(vatom,6,j) = arr2v(vatom,6,j) + v[6];
+              }
+
+//     Now compute forces on other atoms k due to change in sij
+
+              if (iszero(sij) || iszero(sij-1.0)) continue; //: cont jn loop
+              for(kn=1; kn<=*numneigh_full; kn++) {
+                k = arr1v(firstneigh_full,kn);
+                eltk = arr1v(fmap,arr1v(type,k));
+                if (k!=j && eltk > 0) {
+                  dsij(i,j,k,jn,*nmax,*numneigh,rij2,&dsij1,&dsij2,*ntype,type,fmap,x,scrfcn,fcpair);
+                  if (!iszero(dsij1) || !iszero(dsij2)) {
+                    force1 = dUdsij*dsij1;
+                    force2 = dUdsij*dsij2;
+                    for(m=1; m<=3; m++) {
+                      delik[m] = arr2v(x,m,k) - arr2v(x,m,i);
+                      deljk[m] = arr2v(x,m,k) - arr2v(x,m,j);
+                    }
+                    for(m=1; m<=3; m++) {
+                      arr2v(f,m,i) = arr2v(f,m,i) + force1*delik[m];
+                      arr2v(f,m,j) = arr2v(f,m,j) + force2*deljk[m];
+                      arr2v(f,m,k) = arr2v(f,m,k) - force1*delik[m] - force2*deljk[m];
+                    }
+
+//     Tabulate per-atom virial as symmetrized stress tensor
+
+                    if (*vflag_atom!=0) {
+                      fi[1] = force1*delik[1];
+                      fi[2] = force1*delik[2];
+                      fi[3] = force1*delik[3];
+                      fj[1] = force2*deljk[1];
+                      fj[2] = force2*deljk[2];
+                      fj[3] = force2*deljk[3];
+                      v[1] = -third * (delik[1]*fi[1] + deljk[1]*fj[1]);
+                      v[2] = -third * (delik[2]*fi[2] + deljk[2]*fj[2]);
+                      v[3] = -third * (delik[3]*fi[3] + deljk[3]*fj[3]);
+                      v[4] = -sixth * (delik[1]*fi[2] + deljk[1]*fj[2] +                    delik[2]*fi[1] + deljk[2]*fj[1]);
+                      v[5] = -sixth * (delik[1]*fi[3] + deljk[1]*fj[3] +                    delik[3]*fi[1] + deljk[3]*fj[1]);
+                      v[6] = -sixth * (delik[2]*fi[3] + deljk[2]*fj[3] +                    delik[3]*fi[2] + deljk[3]*fj[2]);
+
+                      arr2v(vatom,1,i) = arr2v(vatom,1,i) + v[1];
+                      arr2v(vatom,2,i) = arr2v(vatom,2,i) + v[2];
+                      arr2v(vatom,3,i) = arr2v(vatom,3,i) + v[3];
+                      arr2v(vatom,4,i) = arr2v(vatom,4,i) + v[4];
+                      arr2v(vatom,5,i) = arr2v(vatom,5,i) + v[5];
+                      arr2v(vatom,6,i) = arr2v(vatom,6,i) + v[6];
+                      arr2v(vatom,1,j) = arr2v(vatom,1,j) + v[1];
+                      arr2v(vatom,2,j) = arr2v(vatom,2,j) + v[2];
+                      arr2v(vatom,3,j) = arr2v(vatom,3,j) + v[3];
+                      arr2v(vatom,4,j) = arr2v(vatom,4,j) + v[4];
+                      arr2v(vatom,5,j) = arr2v(vatom,5,j) + v[5];
+                      arr2v(vatom,6,j) = arr2v(vatom,6,j) + v[6];
+                      arr2v(vatom,1,k) = arr2v(vatom,1,k) + v[1];
+                      arr2v(vatom,2,k) = arr2v(vatom,2,k) + v[2];
+                      arr2v(vatom,3,k) = arr2v(vatom,3,k) + v[3];
+                      arr2v(vatom,4,k) = arr2v(vatom,4,k) + v[4];
+                      arr2v(vatom,5,k) = arr2v(vatom,5,k) + v[5];
+                      arr2v(vatom,6,k) = arr2v(vatom,6,k) + v[6];
+                    }
+
+                  }
+                }
+//     end of k loop
+              }
+            }
+          }
+//     end of j loop
+        }
+
+//     else if elti=0, this is not a meam atom
+      }
+
+      }

src/USER-MEAMC/meam_setup_done.c

@@ -0,0 +1,945 @@
+#include "meam.h"
+#include <math.h>
+
+void alloyparams();
+void compute_pair_meam();
+double phi_meam(double r,int a, int b);
+void compute_reference_density();
+void get_shpfcn(double *s /* s(3) */, lattice_t latt);
+void get_tavref(double *t11av,double *t21av,double *t31av,double *t12av,double *t22av,double *t32av, double t11,double t21,double t31,double t12,double t22,double t32, double r,int a,int b,lattice_t latt);
+void get_Zij(int *Zij, lattice_t latt);
+void get_Zij2(int *Zij2, double *a, double*S, lattice_t latt,double cmin,double cmax);
+void get_sijk(double C,int i,int j,int k, double *sijk);
+void get_densref(double r,int a,int b,double *rho01,double *rho11,double *rho21,double *rho31, double *rho02,double *rho12,double *rho22,double *rho32);
+double zbl(double r, int z1, int z2);
+double erose(double r,double re,double alpha,double Ec,double repuls,double attrac,int form);
+void interpolate_meam(int ind);
+double compute_phi(double rij, int elti, int eltj);
+
+// in meam_dens_init
+void fcut(double xi, double *fc);
+// in meam_dens_final
+void G_gam(double Gamma,int ibar,double gsmooth_factor, double *G, int *errorflag);
+
+// Declaration in pair_meam.h:
+//
+// void meam_setup_done(double *)
+//
+// Call from pair_meam.cpp:
+//
+// meam_setup_done(&cutmax)
+//
+
+      void meam_setup_done_(double *cutmax)
+      {
+      int nv2, nv3, m, n, p;
+
+//     Force cutoff
+      meam_data.cutforce = meam_data.rc_meam;
+      meam_data.cutforcesq = meam_data.cutforce*meam_data.cutforce;
+
+//     Pass cutoff back to calling program
+      *cutmax = meam_data.cutforce;
+
+//     Augment t1 term
+      for (int i=1; i<=maxelt; i++)
+        meam_data.t1_meam[i] = meam_data.t1_meam[i] + meam_data.augt1 * 3.0/5.0 * meam_data.t3_meam[i];
+
+//     Compute off-diagonal alloy parameters
+      alloyparams();
+
+// indices and factors for Voight notation
+      nv2 = 1;
+      nv3 = 1;
+      for(m = 1; m<=3; m++) {
+        for(n = m; n<=3; n++) {
+          meam_data.vind2D[m][n] = nv2;
+          meam_data.vind2D[n][m] = nv2;
+          nv2 = nv2+1;
+          for (p = n; p<=3; p++) {
+            meam_data.vind3D[m][n][p] = nv3;
+            meam_data.vind3D[m][p][n] = nv3;
+            meam_data.vind3D[n][m][p] = nv3;
+            meam_data.vind3D[n][p][m] = nv3;
+            meam_data.vind3D[p][m][n] = nv3;
+            meam_data.vind3D[p][n][m] = nv3;
+            nv3 = nv3+1;
+          }
+        }
+      }
+
+      meam_data.v2D[1] = 1;
+      meam_data.v2D[2] = 2;
+      meam_data.v2D[3] = 2;
+      meam_data.v2D[4] = 1;
+      meam_data.v2D[5] = 2;
+      meam_data.v2D[6] = 1;
+
+      meam_data.v3D[1] = 1;
+      meam_data.v3D[2] = 3;
+      meam_data.v3D[3] = 3;
+      meam_data.v3D[4] = 3;
+      meam_data.v3D[5] = 6;
+      meam_data.v3D[6] = 3;
+      meam_data.v3D[7] = 1;
+      meam_data.v3D[8] = 3;
+      meam_data.v3D[9] = 3;
+      meam_data.v3D[10] = 1;
+
+      nv2 = 1;
+      for(m = 1; m<=meam_data.neltypes; m++) {
+        for(n = m; n<=meam_data.neltypes; n++) {
+          meam_data.eltind[m][n] = nv2;
+          meam_data.eltind[n][m] = nv2;
+          nv2 = nv2+1;
+        }
+      }
+
+//     Compute background densities for reference structure
+      compute_reference_density();
+
+//     Compute pair potentials and setup arrays for interpolation
+      meam_data.nr = 1000;
+      meam_data.dr = 1.1*meam_data.rc_meam/meam_data.nr;
+      compute_pair_meam();
+    }
+
+//ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+// Fill off-diagonal alloy parameters
+      void alloyparams(void)
+      {
+      
+      int i,j,k;
+      double eb;
+
+// Loop over pairs
+      for(i = 1; i<=meam_data.neltypes; i++) {
+        for(j = 1;i<=meam_data.neltypes; i++) {
+// Treat off-diagonal pairs
+// If i>j, set all equal to i<j case (which has aready been set,
+// here or in the input file)
+          if (i > j) {
+            meam_data.re_meam[i][j] = meam_data.re_meam[j][i];
+            meam_data.Ec_meam[i][j] = meam_data.Ec_meam[j][i];
+            meam_data.alpha_meam[i][j] = meam_data.alpha_meam[j][i];
+            meam_data.lattce_meam[i][j] = meam_data.lattce_meam[j][i];
+            meam_data.nn2_meam[i][j] = meam_data.nn2_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) {
+            if (iszero(meam_data.Ec_meam[i][j])) {
+              if (meam_data.lattce_meam[i][j]==L12)
+                meam_data.Ec_meam[i][j] = (3*meam_data.Ec_meam[i][i]+meam_data.Ec_meam[j][j])/4.0 - meam_data.delta_meam[i][j];
+              else if (meam_data.lattce_meam[i][j]==C11) {
+                if (meam_data.lattce_meam[i][i]==DIA)
+                  meam_data.Ec_meam[i][j] = (2*meam_data.Ec_meam[i][i]+meam_data.Ec_meam[j][j])/3.0 - meam_data.delta_meam[i][j];
+                else
+                  meam_data.Ec_meam[i][j] = (meam_data.Ec_meam[i][i]+2*meam_data.Ec_meam[j][j])/3.0 - meam_data.delta_meam[i][j];
+              } else
+                meam_data.Ec_meam[i][j] = (meam_data.Ec_meam[i][i]+meam_data.Ec_meam[j][j])/2.0 - meam_data.delta_meam[i][j];
+            }
+            if (iszero(meam_data.alpha_meam[i][j]))
+              meam_data.alpha_meam[i][j] = (meam_data.alpha_meam[i][i]+meam_data.alpha_meam[j][j])/2.0;
+            if (iszero(meam_data.re_meam[i][j]))
+              meam_data.re_meam[i][j] = (meam_data.re_meam[i][i]+meam_data.re_meam[j][j])/2.0;
+          }
+        }
+      }
+
+// Cmin[i][k][j] is symmetric in i-j, but not k.  For all triplets
+// where i>j, set equal to the i<j element.  Likewise for Cmax.
+      for(i = 2;i<=meam_data.neltypes;i++){
+        for(j = 1;j<=i-1;j++){
+          for(k = 1;k<=meam_data.neltypes;k++){
+          meam_data.Cmin_meam[i][j][k] = meam_data.Cmin_meam[j][i][k];
+          meam_data.Cmax_meam[i][j][k] = meam_data.Cmax_meam[j][i][k];
+          }
+        }
+      }
+        
+
+// ebound gives the squared distance such that, for rik2 or rjk2>ebound,
+// atom k definitely lies outside the screening function ellipse (so
+// there is no need to calculate its effects).  Here, compute it for all
+// triplets [i][j][k] so that ebound[i][j] is the maximized over k
+      for(i = 2;i<=meam_data.neltypes;i++){
+        for(j = 1;j<=meam_data.neltypes;j++){
+          for(k = 1;k<=meam_data.neltypes;k++){
+            eb = (meam_data.Cmax_meam[i][j][k]*meam_data.Cmax_meam[i][j][k]) / (4.0*(meam_data.Cmax_meam[i][j][k]-1.0));
+            meam_data.ebound_meam[i][j] = max(meam_data.ebound_meam[i][j],eb);
+          }
+        }
+      }
+    }
+
+//-----------------------------------------------------------------------
+// compute MEAM pair potential for each pair of element types
+//
+
+      void compute_pair_meam(void)
+      {
+
+      double r/*ununsed:, temp*/;
+      int j,a,b,nv2;
+      double astar,frac,phizbl;
+      int n,nmax,Z1,Z2;
+      double arat,rarat,scrn,scrn2;
+      double phiaa,phibb/*unused:,phitmp*/;
+      double C,s111,s112,s221,S11,S22;
+
+// check for previously allocated arrays and free them
+      if(allocated(meam_data.phir)) deallocate(meam_data.phir);
+      if(allocated(meam_data.phirar)) deallocate(meam_data.phirar);
+      if(allocated(meam_data.phirar1)) deallocate(meam_data.phirar1);
+      if(allocated(meam_data.phirar2)) deallocate(meam_data.phirar2);
+      if(allocated(meam_data.phirar3)) deallocate(meam_data.phirar3);
+      if(allocated(meam_data.phirar4)) deallocate(meam_data.phirar4);
+      if(allocated(meam_data.phirar5)) deallocate(meam_data.phirar5);
+      if(allocated(meam_data.phirar6)) deallocate(meam_data.phirar6);
+
+// allocate memory for array that defines the potential
+      allocate_2d(meam_data.phir,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+
+// allocate coeff memory
+
+      allocate_2d(meam_data.phirar,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+      allocate_2d(meam_data.phirar1,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+      allocate_2d(meam_data.phirar2,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+      allocate_2d(meam_data.phirar3,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+      allocate_2d(meam_data.phirar4,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+      allocate_2d(meam_data.phirar5,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+      allocate_2d(meam_data.phirar6,meam_data.nr,(meam_data.neltypes*(meam_data.neltypes+1))/2);
+
+// loop over pairs of element types
+      nv2 = 0;
+      for(a=1; a<=meam_data.neltypes; a++) {
+        for(b=a; b<=meam_data.neltypes; b++) {
+          nv2 = nv2 + 1;
+
+// loop over r values and compute
+          for(j=1; j<=meam_data.nr; j++) {
+            r = (j-1)*meam_data.dr;
+
+            arr2(meam_data.phir,j,nv2) = phi_meam(r,a,b);
+
+// if using second-nearest neighbor, solve recursive problem
+// (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
+            if (meam_data.nn2_meam[a][b]==1) {
+              get_Zij(&Z1,meam_data.lattce_meam[a][b]);
+              get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][b],meam_data.Cmin_meam[a][a][b],meam_data.Cmax_meam[a][a][b]);
+
+//     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
+//     pairs, but need to include NN2 contributions to those pairs as
+//     well.
+              if (meam_data.lattce_meam[a][b]==B1 || meam_data.lattce_meam[a][b]==B2 || meam_data.lattce_meam[a][b]==L12) {
+                rarat = r*arat;
+
+//               phi_aa
+                phiaa = phi_meam(rarat,a,a);
+                get_Zij(&Z1,meam_data.lattce_meam[a][a]);
+                get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][a], meam_data.Cmin_meam[a][a][a],meam_data.Cmax_meam[a][a][a]);
+                nmax = 10;
+                if (scrn > 0.0) {
+                  for(n=1; n<=nmax; n++) {
+                    phiaa = phiaa + pow((-Z2*scrn/Z1),n) * phi_meam(rarat*pow(arat,n),a,a);
+                  }
+                }
+
+//               phi_bb
+                phibb = phi_meam(rarat,b,b);
+                get_Zij(&Z1,meam_data.lattce_meam[b][b]);
+                get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[b][b], meam_data.Cmin_meam[b][b][b],meam_data.Cmax_meam[b][b][b]);
+                nmax = 10;
+                if (scrn > 0.0) {
+                  for(n=1; n<=nmax; n++) {
+                    phibb = phibb + pow((-Z2*scrn/Z1),n) * phi_meam(rarat*pow(arat,n),b,b);
+                  }
+                }
+
+                if (meam_data.lattce_meam[a][b]==B1 || meam_data.lattce_meam[a][b]==B2) {
+//     Add contributions to the B1 or B2 potential
+                  get_Zij(&Z1,meam_data.lattce_meam[a][b]);
+                  get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][b], meam_data.Cmin_meam[a][a][b],meam_data.Cmax_meam[a][a][b]);
+                  arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) - Z2*scrn/(2*Z1) * phiaa;
+                  get_Zij2(&Z2,&arat,&scrn2,meam_data.lattce_meam[a][b], meam_data.Cmin_meam[b][b][a],meam_data.Cmax_meam[b][b][a]);
+                  arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) - Z2*scrn2/(2*Z1) * phibb;
+
+                } else if (meam_data.lattce_meam[a][b]==L12) {
+//     The L12 case has one last trick; we have to be careful to compute
+//     the correct screening between 2nd-neighbor pairs.  1-1
+//     second-neighbor pairs are screened by 2 type 1 atoms and two type
+//     2 atoms.  2-2 second-neighbor pairs are screened by 4 type 1
+//     atoms.
+                  C = 1.0;
+                  get_sijk(C,a,a,a,&s111);
+                  get_sijk(C,a,a,b,&s112);
+                  get_sijk(C,b,b,a,&s221);
+                  S11 = s111 * s111 * s112 * s112;
+                  S22 = pow(s221,4);
+                  arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) - 0.75*S11*phiaa - 0.25*S22*phibb;
+
+                }
+
+              } else {
+                nmax = 10;
+                for(n=1; n<=nmax; n++) {
+                  arr2(meam_data.phir,j,nv2) = arr2(meam_data.phir,j,nv2) + pow((-Z2*scrn/Z1),n) * phi_meam(r*pow(arat,n),a,b);
+                }
+              }
+
+            }
+
+// For Zbl potential:
+// if astar <= -3
+//   potential is zbl potential
+// else if -3 < astar < -1
+//   potential is linear combination with zbl potential
+// endif
+            if (meam_data.zbl_meam[a][b]==1) {
+              astar = meam_data.alpha_meam[a][b] * (r/meam_data.re_meam[a][b] - 1.0);
+              if (astar <= -3.0)
+                arr2(meam_data.phir,j,nv2) = zbl(r,meam_data.ielt_meam[a],meam_data.ielt_meam[b]);
+              else if (astar > -3.0 && astar < -1.0) {
+                fcut(1-(astar+1.0)/(-3.0+1.0),&frac);
+                phizbl = zbl(r,meam_data.ielt_meam[a],meam_data.ielt_meam[b]);
+                arr2(meam_data.phir,j,nv2) = frac*arr2(meam_data.phir,j,nv2) + (1-frac)*phizbl;
+              }
+            }
+
+          }
+
+// call interpolation
+          interpolate_meam(nv2);
+
+        }
+      }
+
+      }
+
+
+//----------------------------------------------------------------------c
+// Compute MEAM pair potential for distance r, element types a and b
+//
+      double phi_meam(double r,int a, int b)
+      {
+      /*unused:double a1,a2,a12;*/
+      double t11av,t21av,t31av,t12av,t22av,t32av;
+      double G1,G2,s1[3+1],s2[3+1]/*,s12[3+1]*/,rho0_1,rho0_2;
+      double Gam1,Gam2,Z1,Z2;
+      double rhobar1,rhobar2,F1,F2;
+      double rho01,rho11,rho21,rho31;
+      double rho02,rho12,rho22,rho32;
+      double scalfac,phiaa,phibb;
+      double Eu;
+      double arat,scrn/*unused:,scrn2*/;
+      int Z12, errorflag;
+      int n,nmax,Z1nn,Z2nn;
+      lattice_t latta/*unused:,lattb*/;
+      double rho_bkgd1, rho_bkgd2;
+      
+      double phi_m = 0.0;
+
+// Equation numbers below refer to:
+//   I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615
+
+// get number of neighbors in the reference structure
+//   Nref[i][j] = # of i's neighbors of type j
+      get_Zij(&Z12,meam_data.lattce_meam[a][b]);
+
+      get_densref(r,a,b,&rho01,&rho11,&rho21,&rho31,&rho02,&rho12,&rho22,&rho32);
+
+// if densities are too small, numerical problems may result; just return zero
+      if (rho01<=1e-14 && rho02<=1e-14)
+        return 0.0;
+
+// calculate average weighting factors for the reference structure
+      if (meam_data.lattce_meam[a][b]==C11) {
+        if (meam_data.ialloy==2) {
+          t11av = meam_data.t1_meam[a];
+          t12av = meam_data.t1_meam[b];
+          t21av = meam_data.t2_meam[a];
+          t22av = meam_data.t2_meam[b];
+          t31av = meam_data.t3_meam[a];
+          t32av = meam_data.t3_meam[b];
+        } else {
+          scalfac = 1.0/(rho01+rho02);
+          t11av = scalfac*(meam_data.t1_meam[a]*rho01 + meam_data.t1_meam[b]*rho02);
+          t12av = t11av;
+          t21av = scalfac*(meam_data.t2_meam[a]*rho01 + meam_data.t2_meam[b]*rho02);
+          t22av = t21av;
+          t31av = scalfac*(meam_data.t3_meam[a]*rho01 + meam_data.t3_meam[b]*rho02);
+          t32av = t31av;
+        }
+      } else {
+// average weighting factors for the reference structure, eqn. I.8
+         get_tavref(&t11av,&t21av,&t31av,&t12av,&t22av,&t32av, meam_data.t1_meam[a],meam_data.t2_meam[a],meam_data.t3_meam[a], meam_data.t1_meam[b],meam_data.t2_meam[b],meam_data.t3_meam[b], r,a,b,meam_data.lattce_meam[a][b]);
+      }
+
+// for c11b structure, calculate background electron densities
+      if (meam_data.lattce_meam[a][b]==C11) {
+         latta = meam_data.lattce_meam[a][a];
+         if (latta==DIA) {
+            rhobar1 = pow(((Z12/2)*(rho02+rho01)),2) + t11av*pow((rho12-rho11),2) + t21av/6.0*pow(rho22+rho21,2)+ 121.0/40.0*t31av*pow((rho32-rho31),2);
+            rhobar1 = sqrt(rhobar1);
+            rhobar2 = pow(Z12*rho01,2) + 2.0/3.0*t21av*pow(rho21,2);
+            rhobar2 = sqrt(rhobar2);
+         } else {
+            rhobar2 = pow(((Z12/2)*(rho01+rho02)),2) + t12av*pow((rho11-rho12),2) + t22av/6.0*pow(rho21+rho22,2) + 121.0/40.0*t32av*pow((rho31-rho32),2);
+            rhobar2 = sqrt(rhobar2);
+            rhobar1 = pow(Z12*rho02,2) + 2.0/3.0*t22av*pow(rho22,2);
+            rhobar1 = sqrt(rhobar1);
+         }
+      } else {
+// for other structures, use formalism developed in Huang's paper
+//
+//     composition-dependent scaling, equation I.7
+//     If using mixing rule for t, apply to reference structure; else
+//     use precomputed values
+        if (meam_data.mix_ref_t==1) {
+          Z1 = meam_data.Z_meam[a];
+          Z2 = meam_data.Z_meam[b];
+          if (meam_data.ibar_meam[a]<=0)
+            G1 = 1.0;
+          else {
+            get_shpfcn(s1,meam_data.lattce_meam[a][a]);
+            Gam1 = (s1[1]*t11av+s1[2]*t21av+s1[3]*t31av)/(Z1*Z1);
+            G_gam(Gam1,meam_data.ibar_meam[a],meam_data.gsmooth_factor,&G1,&errorflag);
+          }
+          if (meam_data.ibar_meam[b]<=0)
+            G2 = 1.0;
+          else {
+            get_shpfcn(s2,meam_data.lattce_meam[b][b]);
+            Gam2 = (s2[1]*t12av+s2[2]*t22av+s2[3]*t32av)/(Z2*Z2);
+            G_gam(Gam2,meam_data.ibar_meam[b],meam_data.gsmooth_factor,&G2,&errorflag);
+          }
+          rho0_1 = meam_data.rho0_meam[a]*Z1*G1;
+          rho0_2 = meam_data.rho0_meam[b]*Z2*G2;
+        }
+        Gam1 = (t11av*rho11+t21av*rho21+t31av*rho31);
+        if (rho01 < 1.0e-14)
+          Gam1 = 0.0;
+        else
+          Gam1 = Gam1/(rho01*rho01);
+
+        Gam2 = (t12av*rho12+t22av*rho22+t32av*rho32);
+        if (rho02 < 1.0e-14)
+          Gam2 = 0.0;
+        else
+          Gam2 = Gam2/(rho02*rho02);
+
+        G_gam(Gam1,meam_data.ibar_meam[a],meam_data.gsmooth_factor,&G1,&errorflag);
+        G_gam(Gam2,meam_data.ibar_meam[b],meam_data.gsmooth_factor,&G2,&errorflag);
+        if (meam_data.mix_ref_t==1) {
+          rho_bkgd1 = rho0_1;
+          rho_bkgd2 = rho0_2;
+        } else {
+          if (meam_data.bkgd_dyn==1) {
+            rho_bkgd1 = meam_data.rho0_meam[a]*meam_data.Z_meam[a];
+            rho_bkgd2 = meam_data.rho0_meam[b]*meam_data.Z_meam[b];
+          } else {
+            rho_bkgd1 = meam_data.rho_ref_meam[a];
+            rho_bkgd2 = meam_data.rho_ref_meam[b];
+          }
+        }
+        rhobar1 = rho01/rho_bkgd1*G1;
+        rhobar2 = rho02/rho_bkgd2*G2;
+
+      }
+
+// compute embedding functions, eqn I.5
+      if (iszero(rhobar1))
+        F1 = 0.0;
+      else {
+        if (meam_data.emb_lin_neg==1 && rhobar1<=0)
+          F1 = -meam_data.A_meam[a]*meam_data.Ec_meam[a][a]*rhobar1;
+        else
+          F1 = meam_data.A_meam[a]*meam_data.Ec_meam[a][a]*rhobar1*log(rhobar1);
+      }
+      if (iszero(rhobar2))
+        F2 = 0.0;
+      else {
+        if (meam_data.emb_lin_neg==1  &&  rhobar2<=0)
+          F2 = -meam_data.A_meam[b]*meam_data.Ec_meam[b][b]*rhobar2;
+        else
+          F2 = meam_data.A_meam[b]*meam_data.Ec_meam[b][b]*rhobar2*log(rhobar2);
+      }
+
+// compute Rose function, I.16
+      Eu = erose(r,meam_data.re_meam[a][b],meam_data.alpha_meam[a][b], meam_data.Ec_meam[a][b],meam_data.repuls_meam[a][b],meam_data.attrac_meam[a][b],meam_data.erose_form);
+
+// calculate the pair energy
+      if (meam_data.lattce_meam[a][b]==C11) {
+        latta = meam_data.lattce_meam[a][a];
+        if (latta==DIA) {
+          phiaa = phi_meam(r,a,a);
+          phi_m = (3*Eu - F2 - 2*F1 - 5*phiaa)/Z12;
+        } else {
+          phibb = phi_meam(r,b,b);
+          phi_m = (3*Eu - F1 - 2*F2 - 5*phibb)/Z12;
+        }
+      } else if (meam_data.lattce_meam[a][b]==L12) {
+        phiaa = phi_meam(r,a,a);
+//       account for second neighbor a-a potential here...
+        get_Zij(&Z1nn,meam_data.lattce_meam[a][a]);
+        get_Zij2(&Z2nn,&arat,&scrn,meam_data.lattce_meam[a][a],meam_data.Cmin_meam[a][a][a],meam_data.Cmax_meam[a][a][a]);
+        nmax = 10;
+        if (scrn > 0.0) {
+          for(n=1; n<=nmax; n++) {
+            phiaa = phiaa + pow((-Z2nn*scrn/Z1nn),n) * phi_meam(r*pow(arat,n),a,a);
+          }
+        }
+        phi_m = Eu/3.0 - F1/4.0 - F2/12.0 - phiaa;
+         
+      } else {
+//
+// potential is computed from Rose function and embedding energy
+         phi_m = (2*Eu - F1 - F2)/Z12;
+//
+      }
+
+// if r = 0, just return 0
+      if (iszero(r)) {
+        phi_m = 0.0;
+      }
+      
+      return phi_m;
+    }
+
+//----------------------------------------------------------------------c
+// Compute background density for reference structure of each element
+      void compute_reference_density(void)
+      {
+      int a,Z,Z2,errorflag;
+      double gam,Gbar,shp[3+1];
+      double rho0,rho0_2nn,arat,scrn;
+
+// loop over element types
+      for(a=1; a<=meam_data.neltypes; a++) {
+        Z = (int)meam_data.Z_meam[a];
+        if (meam_data.ibar_meam[a]<=0)
+          Gbar = 1.0;
+        else {
+          get_shpfcn(shp,meam_data.lattce_meam[a][a]);
+          gam = (meam_data.t1_meam[a]*shp[1]+meam_data.t2_meam[a]*shp[2]+meam_data.t3_meam[a]*shp[3])/(Z*Z);
+          G_gam(gam,meam_data.ibar_meam[a],meam_data.gsmooth_factor,&Gbar,&errorflag);
+        }
+
+//     The zeroth order density in the reference structure, with
+//     equilibrium spacing, is just the number of first neighbors times
+//     the rho0_meam coefficient...
+        rho0 = meam_data.rho0_meam[a]*Z;
+
+//     ...unless we have unscreened second neighbors, in which case we
+//     add on the contribution from those (accounting for partial
+//     screening)
+        if (meam_data.nn2_meam[a][a]==1) {
+          get_Zij2(&Z2,&arat,&scrn,meam_data.lattce_meam[a][a],meam_data.Cmin_meam[a][a][a],meam_data.Cmax_meam[a][a][a]);
+          rho0_2nn = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*(arat-1));
+          rho0 = rho0 + Z2*rho0_2nn*scrn;
+        }
+
+        meam_data.rho_ref_meam[a] = rho0*Gbar;
+      }
+      }
+
+//----------------------------------------------------------------------c
+// Shape factors for various configurations
+      void get_shpfcn(double *s /* s(3) */, lattice_t latt)
+      {
+      if (latt==FCC || latt==BCC || latt==B1 || latt==B2) {
+        s[1] = 0.0;
+        s[2] = 0.0;
+        s[3] = 0.0;
+      } else if (latt==HCP) {
+        s[1] = 0.0;
+        s[2] = 0.0;
+        s[3] = 1.0/3.0;
+      } else if (latt==DIA) {
+        s[1] = 0.0;
+        s[2] = 0.0;
+        s[3] = 32.0/9.0;
+      } else if (latt==DIM) {
+        s[1] = 1.0;
+        s[2] = 2.0/3.0;
+//        s(3) = 1.d0
+        s[3] = 0.40;
+      } else {
+        s[1] = 0.0;
+//        call error('Lattice not defined in get_shpfcn.')
+      }
+      }
+      
+//------------------------------------------------------------------------------c
+// Average weighting factors for the reference structure
+      void get_tavref(double *t11av,double *t21av,double *t31av,double *t12av,double *t22av,double *t32av, double t11,double t21,double t31,double t12,double t22,double t32, double r,int a,int b,lattice_t latt)
+      {
+      double rhoa01,rhoa02,a1,a2,rho01/*,rho02*/;
+
+//     For ialloy = 2, no averaging is done
+      if (meam_data.ialloy==2) {
+          *t11av = t11;
+          *t21av = t21;
+          *t31av = t31;
+          *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/meam_data.re_meam[a][a] - 1.0;
+          a2  = r/meam_data.re_meam[b][b] - 1.0;
+          rhoa01 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*a1);
+          rhoa02 = meam_data.rho0_meam[b]*fm_exp(-meam_data.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.')
+          }
+        }
+      }
+      }
+      
+//------------------------------------------------------------------------------c
+// Number of neighbors for the reference structure
+      void get_Zij(int *Zij, lattice_t latt)
+      {
+      if (latt==FCC)
+        *Zij = 12;
+      else if (latt==BCC)
+        *Zij = 8;
+      else if (latt==HCP)
+        *Zij = 12;
+      else if (latt==B1)
+        *Zij = 6;
+      else if (latt==DIA)
+        *Zij = 4;
+      else if (latt==DIM)
+        *Zij = 1;
+      else if (latt==C11)
+        *Zij = 10;
+      else if (latt==L12)
+        *Zij = 12;
+      else if (latt==B2)
+        *Zij = 8;
+      else {
+//        call error('Lattice not defined in get_Zij.')
+      }
+      }
+     
+//------------------------------------------------------------------------------c
+// Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second
+// neighbor screening function for lattice type "latt"
+
+      void get_Zij2(int *Zij2, double *a, double*S, lattice_t latt,double cmin,double cmax)
+      {
+        
+      double /*rratio,*/C,x,sijk;
+      int numscr = 0;
+
+      if (latt==BCC) {
+        *Zij2 = 6;
+        *a = 2.0/sqrt(3.0);
+        numscr = 4;
+      } else if (latt==FCC) {
+        *Zij2 = 6;
+        *a = sqrt(2.0);
+        numscr = 4;
+      } else if (latt==DIA) {
+        *Zij2 = 0;
+        *a = sqrt(8.0/3.0);
+        numscr = 4;
+        if (cmin < 0.500001) {
+//          call error('can not do 2NN MEAM for dia')
+        }
+      } else if (latt==HCP) {
+        *Zij2 = 6;
+        *a = sqrt(2.0);
+        numscr = 4;
+      } else if (latt==B1) {
+        *Zij2 = 12;
+        *a = sqrt(2.0);
+        numscr = 2;
+      } else if (latt==L12) {
+        *Zij2 = 6;
+        *a = sqrt(2.0);
+        numscr = 4;
+      } else if (latt==B2) {
+        *Zij2 = 6;
+        *a = 2.0/sqrt(3.0);
+        numscr = 4;
+      } else if (latt==DIM) {
+//        this really shouldn't be allowed; make sure screening is zero
+         *Zij2 = 0;
+         *a = 1;
+         *S = 0;
+         return;
+      } else {
+//        call error('Lattice not defined in get_Zij2.')
+      }
+
+// Compute screening for each first neighbor
+      C = 4.0/(*a * *a) - 1.0;
+      x = (C-cmin)/(cmax-cmin);
+      fcut(x,&sijk);
+// There are numscr first neighbors screening the second neighbors
+      *S = pow(sijk,numscr);
+      
+      }
+      
+
+
+//------------------------------------------------------------------------------c
+      void get_sijk(double C,int i,int j,int k, double *sijk)
+      {
+      double x;
+      x = (C-meam_data.Cmin_meam[i][j][k])/(meam_data.Cmax_meam[i][j][k]-meam_data.Cmin_meam[i][j][k]);
+      fcut(x,sijk);
+      }
+
+//------------------------------------------------------------------------------c
+// Calculate density functions, assuming reference configuration
+      void get_densref(double r,int a,int b,double *rho01,double *rho11,double *rho21,double *rho31, double *rho02,double *rho12,double *rho22,double *rho32)
+      {
+      double a1,a2;
+      double s[3+1];
+      lattice_t lat;
+      int  Zij1nn,Zij2nn;
+      double rhoa01nn,rhoa02nn;
+      double rhoa01,rhoa11,rhoa21,rhoa31;
+      double rhoa02,rhoa12,rhoa22,rhoa32;
+      double arat,scrn,denom;
+      double C,s111,s112,s221,S11,S22;
+
+      a1  = r/meam_data.re_meam[a][a] - 1.0;
+      a2  = r/meam_data.re_meam[b][b] - 1.0;
+
+      rhoa01 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*a1);
+      rhoa11 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta1_meam[a]*a1);
+      rhoa21 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta2_meam[a]*a1);
+      rhoa31 = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta3_meam[a]*a1);
+      rhoa02 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta0_meam[b]*a2);
+      rhoa12 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta1_meam[b]*a2);
+      rhoa22 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta2_meam[b]*a2);
+      rhoa32 = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta3_meam[b]*a2);
+
+      lat = meam_data.lattce_meam[a][b];
+
+      *rho11 = 0.0;
+      *rho21 = 0.0;
+      *rho31 = 0.0;
+      *rho12 = 0.0;
+      *rho22 = 0.0;
+      *rho32 = 0.0;
+
+      get_Zij(&Zij1nn,lat);
+
+      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(s,DIM);
+        *rho01 = rhoa02;
+        *rho02 = rhoa01;
+        *rho11 = s[1]*rhoa12*rhoa12;
+        *rho12 = s[1]*rhoa11*rhoa11;
+        *rho21 = s[2]*rhoa22*rhoa22;
+        *rho22 = s[2]*rhoa21*rhoa21;
+        *rho31 = s[3]*rhoa32*rhoa32;
+        *rho32 = s[3]*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 (meam_data.ialloy==1) {
+          *rho21 = 8./3.*pow(rhoa21*meam_data.t2_meam[a]-rhoa22*meam_data.t2_meam[b],2);
+          denom = 8*rhoa01*pow(meam_data.t2_meam[a],2) + 4*rhoa02*pow(meam_data.t2_meam[b],2);
+          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 {
+//        call error('Lattice not defined in get_densref.')
+      }
+
+      if (meam_data.nn2_meam[a][b]==1) {
+
+        get_Zij2(&Zij2nn,&arat,&scrn,lat,meam_data.Cmin_meam[a][a][b],meam_data.Cmax_meam[a][a][b]);
+
+        a1 = arat*r/meam_data.re_meam[a][a] - 1.0;
+        a2 = arat*r/meam_data.re_meam[b][b] - 1.0;
+
+        rhoa01nn = meam_data.rho0_meam[a]*fm_exp(-meam_data.beta0_meam[a]*a1);
+        rhoa02nn = meam_data.rho0_meam[b]*fm_exp(-meam_data.beta0_meam[b]*a2);
+
+        if (lat==L12) {
+//     As usual, L12 thinks it's special; we need to be careful computing
+//     the screening functions
+          C = 1.0;
+          get_sijk(C,a,a,a,&s111);
+          get_sijk(C,a,a,b,&s112);
+          get_sijk(C,b,b,a,&s221);
+          S11 = s111 * s111 * s112 * s112;
+          S22 = pow(s221,4);
+          *rho01 = *rho01 + 6*S11*rhoa01nn;
+          *rho02 = *rho02 + 6*S22*rhoa02nn;
+
+        } else {
+//     For other cases, assume that second neighbor is of same type,
+//     first neighbor may be of different type
+
+          *rho01 = *rho01 + Zij2nn*scrn*rhoa01nn;
+
+//     Assume Zij2nn and arat don't depend on order, but scrn might
+          get_Zij2(&Zij2nn,&arat,&scrn,lat,meam_data.Cmin_meam[b][b][a],meam_data.Cmax_meam[b][b][a]);
+          *rho02 = *rho02 + Zij2nn*scrn*rhoa02nn;
+        }
+      }
+      }
+
+//---------------------------------------------------------------------
+// Compute ZBL potential
+//
+      double zbl(double r, int z1, int z2)
+      {
+      int i;
+      const double c[]={0.028171,0.28022,0.50986,0.18175};
+      const double d[]={0.20162,0.40290,0.94229,3.1998};
+      const double azero=0.4685;
+      const double cc=14.3997;
+      double a,x;
+// azero = (9pi^2/128)^1/3 (0.529) Angstroms
+      a = azero/(pow(z1,0.23)+pow(z2,0.23));
+      double result = 0.0;
+      x = r/a;
+      for(i=0; i<=3; i++) {
+        result = result + c[i]*fm_exp(-d[i]*x);
+      }
+      if (r > 0.0) result = result*z1*z2/r*cc;
+      return result;
+      }
+
+//---------------------------------------------------------------------
+// Compute Rose energy function, I.16
+//
+      double erose(double r,double re,double alpha,double Ec,double repuls,double attrac,int form)
+      {
+      double astar,a3;
+      double result = 0.0;
+
+      if (r > 0.0) {
+        astar = alpha * (r/re - 1.0);
+        a3 = 0.0;
+        if (astar>=0)
+          a3 = attrac;
+        else if (astar < 0)
+          a3 = repuls;
+
+        if (form==1)
+          result = -Ec*(1+astar+(-attrac+repuls/r)* pow(astar,3))*fm_exp(-astar);
+        else if (form==2)
+          result = -Ec * (1 +astar + a3*pow(astar,3))*fm_exp(-astar);
+        else
+          result = -Ec * (1+ astar + a3*pow(astar,3)/(r/re))*fm_exp(-astar);
+       
+      }
+      return result;
+      }
+
+// -----------------------------------------------------------------------
+
+      void interpolate_meam(int ind)
+      {
+      int j;
+      double drar;
+
+// map to coefficient space
+
+      meam_data.nrar = meam_data.nr;
+      drar = meam_data.dr;
+      meam_data.rdrar = 1.0/drar;
+
+// phir interp
+      for(j=1; j<=meam_data.nrar; j++) {
+        arr2(meam_data.phirar,j,ind) = arr2(meam_data.phir,j,ind);
+      }
+      arr2(meam_data.phirar1,1,ind) = arr2(meam_data.phirar,2,ind)-arr2(meam_data.phirar,1,ind);
+      arr2(meam_data.phirar1,2,ind) = 0.5*(arr2(meam_data.phirar,3,ind)-arr2(meam_data.phirar,1,ind));
+      arr2(meam_data.phirar1,meam_data.nrar-1,ind) = 0.5*(arr2(meam_data.phirar,meam_data.nrar,ind)      -arr2(meam_data.phirar,meam_data.nrar-2,ind));
+      arr2(meam_data.phirar1,meam_data.nrar,ind) = 0.0;
+      for(j=3; j<=meam_data.nrar-2; j++) {
+        arr2(meam_data.phirar1,j,ind) = ((arr2(meam_data.phirar,j-2,ind)-arr2(meam_data.phirar,j+2,ind)) +        8.0*(arr2(meam_data.phirar,j+1,ind)-arr2(meam_data.phirar,j-1,ind)))/12.;
+      }
+
+      for(j=1; j<=meam_data.nrar-1; j++) {
+        arr2(meam_data.phirar2,j,ind) = 3.0*(arr2(meam_data.phirar,j+1,ind)-arr2(meam_data.phirar,j,ind)) -        2.0*arr2(meam_data.phirar1,j,ind) - arr2(meam_data.phirar1,j+1,ind);
+        arr2(meam_data.phirar3,j,ind) = arr2(meam_data.phirar1,j,ind) + arr2(meam_data.phirar1,j+1,ind) -        2.0*(arr2(meam_data.phirar,j+1,ind)-arr2(meam_data.phirar,j,ind));
+      }
+      arr2(meam_data.phirar2,meam_data.nrar,ind) = 0.0;
+      arr2(meam_data.phirar3,meam_data.nrar,ind) = 0.0;
+
+      for(j=1; j<=meam_data.nrar; j++) {
+        arr2(meam_data.phirar4,j,ind) = arr2(meam_data.phirar1,j,ind)/drar;
+        arr2(meam_data.phirar5,j,ind) = 2.0*arr2(meam_data.phirar2,j,ind)/drar;
+        arr2(meam_data.phirar6,j,ind) = 3.0*arr2(meam_data.phirar3,j,ind)/drar;
+      }  
+           
+      }
+
+//---------------------------------------------------------------------
+// Compute Rose energy function, I.16
+//
+      double compute_phi(double rij, int elti, int eltj)
+      {
+      double pp;
+      int ind, kk;
+
+      ind = meam_data.eltind[elti][eltj];
+      pp = rij*meam_data.rdrar + 1.0;
+      kk = (int)pp;
+      kk = min(kk,meam_data.nrar-1);
+      pp = pp - kk;
+      pp = min(pp,1.0);
+      double result = ((arr2(meam_data.phirar3,kk,ind)*pp + arr2(meam_data.phirar2,kk,ind))*pp + arr2(meam_data.phirar1,kk,ind))*pp + arr2(meam_data.phirar,kk,ind);
+
+      return result;
+      }

src/USER-MEAMC/meam_setup_global.c

@@ -0,0 +1,99 @@
+#include "meam.h"
+#include <math.h>
+
+//
+// declaration in pair_meam.h:
+//
+//  void meam_setup_global(int *, int *, double *, int *, double *, double *,
+//			 double *, double *, double *, double *, double *,
+//			 double *, double *, double *, double *, double *,
+//			 double *, double *, int *);
+//
+// call in pair_meam.cpp:
+//
+//  meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
+//		    alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
+//
+//
+
+      void meam_setup_global_(int *nelt, int *lat, double *z, int *ielement, double *atwt, double *alpha,
+          double *b0, double *b1, double *b2, double *b3, double *alat, double *esub, double *asub,
+          double *t0, double *t1, double *t2, double *t3, double *rozero, int *ibar)
+      {
+      
+      int i;
+      double tmplat[maxelt];
+
+      meam_data.neltypes = *nelt;
+
+      for(i=1; i<=*nelt; i++) {
+         if (arr1v(lat,i)==0)
+            meam_data.lattce_meam[i][i] = FCC;
+         else if (arr1v(lat,i)==1)
+            meam_data.lattce_meam[i][i] = BCC;
+         else if (arr1v(lat,i)==2)
+            meam_data.lattce_meam[i][i] = HCP;
+         else if (arr1v(lat,i)==3)
+            meam_data.lattce_meam[i][i] = DIM;
+         else if (arr1v(lat,i)==4)
+            meam_data.lattce_meam[i][i] = DIA;
+         else {
+//           unknown
+         }
+
+         meam_data.Z_meam[i] = arr1v(z,i);
+         meam_data.ielt_meam[i] = arr1v(ielement,i);
+         meam_data.alpha_meam[i][i] = arr1v(alpha,i);
+         meam_data.beta0_meam[i] = arr1v(b0,i);
+         meam_data.beta1_meam[i] = arr1v(b1,i);
+         meam_data.beta2_meam[i] = arr1v(b2,i);
+         meam_data.beta3_meam[i] = arr1v(b3,i);
+         tmplat[i] = arr1v(alat,i);
+         meam_data.Ec_meam[i][i] = arr1v(esub,i);
+         meam_data.A_meam[i] = arr1v(asub,i);
+         meam_data.t0_meam[i] = arr1v(t0,i);
+         meam_data.t1_meam[i] = arr1v(t1,i);
+         meam_data.t2_meam[i] = arr1v(t2,i);
+         meam_data.t3_meam[i] = arr1v(t3,i);
+         meam_data.rho0_meam[i] = arr1v(rozero,i);
+         meam_data.ibar_meam[i] = arr1v(ibar,i);
+
+         if (meam_data.lattce_meam[i][i]==FCC)
+            meam_data.re_meam[i][i] = tmplat[i]/sqrt(2.0);
+         else if (meam_data.lattce_meam[i][i]==BCC)
+            meam_data.re_meam[i][i] = tmplat[i]*sqrt(3.0)/2.0;
+         else if (meam_data.lattce_meam[i][i]==HCP)
+            meam_data.re_meam[i][i] = tmplat[i];
+         else if (meam_data.lattce_meam[i][i]==DIM)
+            meam_data.re_meam[i][i] = tmplat[i];
+         else if (meam_data.lattce_meam[i][i]==DIA)
+            meam_data.re_meam[i][i] = tmplat[i]*sqrt(3.0)/4.0;
+         else {
+//           error
+         }
+
+      }
+
+
+// Set some defaults
+      meam_data.rc_meam = 4.0;
+      meam_data.delr_meam = 0.1;
+      setall2d(meam_data.attrac_meam, 0.0);
+      setall2d(meam_data.repuls_meam, 0.0);
+      setall3d(meam_data.Cmax_meam, 2.8);
+      setall3d(meam_data.Cmin_meam, 2.0);
+      setall2d(meam_data.ebound_meam, pow(2.8,2)/(4.0*(2.8-1.0)));
+      setall2d(meam_data.delta_meam, 0.0);
+      setall2d(meam_data.nn2_meam, 0);
+      setall2d(meam_data.zbl_meam, 1);
+      meam_data.gsmooth_factor = 99.0;
+      meam_data.augt1 = 1;
+      meam_data.ialloy = 0;
+      meam_data.mix_ref_t = 0;
+      meam_data.emb_lin_neg = 0;
+      meam_data.bkgd_dyn = 0;
+      meam_data.erose_form = 0;
+
+      }
+
+

src/USER-MEAMC/meam_setup_param.c

@@ -0,0 +1,223 @@
+#include "meam.h"
+
+//
+//     do a sanity check on index parameters
+      void meam_checkindex(int num, int lim, int nidx, int *idx /*idx(3)*/, int *ierr)
+      {
+        //: idx[0..2]
+        *ierr = 0;
+        if (nidx < num) {
+           *ierr = 2;
+           return;
+        }
+
+        for (int i=0; i<num; i++) {
+          if ((idx[i] < 1) || (idx[i] > lim)) {
+              *ierr = 3;
+              return;
+          }
+        }
+      }
+
+//
+//     Declaration in pair_meam.h:
+//
+//     void meam_setup_param(int *, double *, int *, int *, int *);
+//
+//     in pair_meam.cpp
+//
+//     meam_setup_param(&which,&value,&nindex,index,&errorflag);
+//
+//
+//
+//     The "which" argument corresponds to the index of the "keyword" array
+//     in pair_meam.cpp:
+//
+//     0 = Ec_meam
+//     1 = alpha_meam
+//     2 = rho0_meam
+//     3 = delta_meam
+//     4 = lattce_meam
+//     5 = attrac_meam
+//     6 = repuls_meam
+//     7 = nn2_meam
+//     8 = Cmin_meam
+//     9 = Cmax_meam
+//     10 = rc_meam
+//     11 = delr_meam
+//     12 = augt1
+//     13 = gsmooth_factor
+//     14 = re_meam
+//     15 = ialloy
+//     16 = mixture_ref_t
+//     17 = erose_form
+//     18 = zbl_meam
+//     19 = emb_lin_neg
+//     20 = bkgd_dyn
+
+      void meam_setup_param_(int *which_p, double *value_p, int *nindex_p, int *index /*index(3)*/, int *errorflag)
+      {
+        //: index[0..2]
+        int i1, i2;
+        *errorflag = 0;
+        int which = *which_p;
+        double value = *value_p;
+        int nindex = *nindex_p;
+
+      switch(which) {
+  //     0 = Ec_meam
+        case 0:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.Ec_meam[index[0]][index[1]] = value;
+          break;
+
+  //     1 = alpha_meam
+        case 1:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.alpha_meam[index[0]][index[1]] = value;
+          break;
+
+  //     2 = rho0_meam
+        case 2:
+          meam_checkindex(1,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.rho0_meam[index[0]] = value;
+          break;
+
+  //     3 = delta_meam
+        case 3:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.delta_meam[index[0]][index[1]] = value;
+          break;
+
+  //     4 = lattce_meam
+        case 4:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          int val = (int)value;
+
+          if (val==0)
+            meam_data.lattce_meam[index[0]][index[1]] = FCC;
+          else if (val==1)
+            meam_data.lattce_meam[index[0]][index[1]] = BCC;
+          else if (val==2)
+            meam_data.lattce_meam[index[0]][index[1]] = HCP;
+          else if (val==3)
+            meam_data.lattce_meam[index[0]][index[1]] = DIM;
+          else if (val==4)
+            meam_data.lattce_meam[index[0]][index[1]] = DIA;
+          else if (val==5)
+            meam_data.lattce_meam[index[0]][index[1]] = B1;
+          else if (val==6)
+            meam_data.lattce_meam[index[0]][index[1]] = C11;
+          else if (val==7)
+            meam_data.lattce_meam[index[0]][index[1]] = L12;
+          else if (val==8)
+            meam_data.lattce_meam[index[0]][index[1]] = B2;
+          break;
+
+  //     5 = attrac_meam
+        case 5:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.attrac_meam[index[0]][index[1]] = value;
+          break;
+
+  //     6 = repuls_meam
+        case 6:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.repuls_meam[index[0]][index[1]] = value;
+          break;
+
+  //     7 = nn2_meam
+        case 7:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          i1 = min(index[0],index[1]);
+          i2 = max(index[0],index[1]);
+          meam_data.nn2_meam[i1][i2] = (int)value;
+          break;
+
+  //     8 = Cmin_meam
+        case 8:
+          meam_checkindex(3,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.Cmin_meam[index[0]][index[1]][index[2]] = value;
+          break;
+
+  //     9 = Cmax_meam
+        case 9:
+          meam_checkindex(3,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.Cmax_meam[index[0]][index[1]][index[2]] = value;
+          break;
+
+  //     10 = rc_meam
+        case 10:
+          meam_data.rc_meam = value;
+          break;
+
+  //     11 = delr_meam
+        case 11:
+          meam_data.delr_meam = value;
+          break;
+
+  //     12 = augt1
+        case 12:
+          meam_data.augt1 = (int)value;
+          break;
+
+  //     13 = gsmooth
+        case 13:
+          meam_data.gsmooth_factor = value;
+          break;
+
+  //     14 = re_meam
+        case 14:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          meam_data.re_meam[index[0]][index[1]] = value;
+          break;
+
+  //     15 = ialloy
+        case 15:
+          meam_data.ialloy = (int)value;
+          break;
+
+  //     16 = mixture_ref_t
+        case 16:
+          meam_data.mix_ref_t = (int)value;
+          break;
+
+  //     17 = erose_form
+        case 17:
+          meam_data.erose_form = (int)value;
+          break;
+
+  //     18 = zbl_meam
+        case 18:
+          meam_checkindex(2,maxelt,nindex,index,errorflag);
+          if (*errorflag!=0) return;
+          i1 = min(index[0],index[1]);
+          i2 = max(index[0],index[1]);
+          meam_data.zbl_meam[i1][i2] = (int)value;
+          break;
+
+  //     19 = emb_lin_neg
+        case 19:
+          meam_data.emb_lin_neg = (int)value;
+          break;
+
+  //     20 = bkgd_dyn
+        case 20:
+          meam_data.bkgd_dyn = (int)value;
+          break;
+
+        default:
+          *errorflag = 1;
+      }
+    }

src/USER-MEAMC/pair_meamc.cpp

@@ -0,0 +1,945 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Greg Wagner (SNL)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "pair_meamc.h"
+#include "atom.h"
+#include "force.h"
+#include "comm.h"
+#include "memory.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MAXLINE 1024
+
+enum{FCC,BCC,HCP,DIM,DIAMOND,B1,C11,L12,B2};
+static const int nkeywords = 21;
+static const char *keywords[] = {
+  "Ec","alpha","rho0","delta","lattce",
+  "attrac","repuls","nn2","Cmin","Cmax","rc","delr",
+  "augt1","gsmooth_factor","re","ialloy",
+  "mixture_ref_t","erose_form","zbl",
+  "emb_lin_neg","bkgd_dyn"};
+
+/* ---------------------------------------------------------------------- */
+
+PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp)
+{
+  single_enable = 0;
+  restartinfo = 0;
+  one_coeff = 1;
+  manybody_flag = 1;
+
+  nmax = 0;
+  rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
+  gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
+  arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
+
+  maxneigh = 0;
+  allocated = 0;
+  scrfcn = dscrfcn = fcpair = NULL;
+
+  nelements = 0;
+  elements = NULL;
+  mass = NULL;
+
+  // set comm size needed by this Pair
+
+  comm_forward = 38;
+  comm_reverse = 30;
+}
+
+/* ----------------------------------------------------------------------
+   free all arrays
+   check if allocated, since class can be destructed when incomplete
+------------------------------------------------------------------------- */
+
+PairMEAMC::~PairMEAMC()
+{
+  meam_cleanup_();
+
+  memory->destroy(rho);
+  memory->destroy(rho0);
+  memory->destroy(rho1);
+  memory->destroy(rho2);
+  memory->destroy(rho3);
+  memory->destroy(frhop);
+  memory->destroy(gamma);
+  memory->destroy(dgamma1);
+  memory->destroy(dgamma2);
+  memory->destroy(dgamma3);
+  memory->destroy(arho2b);
+
+  memory->destroy(arho1);
+  memory->destroy(arho2);
+  memory->destroy(arho3);
+  memory->destroy(arho3b);
+  memory->destroy(t_ave);
+  memory->destroy(tsq_ave);
+
+  memory->destroy(scrfcn);
+  memory->destroy(dscrfcn);
+  memory->destroy(fcpair);
+
+  for (int i = 0; i < nelements; i++) delete [] elements[i];
+  delete [] elements;
+  delete [] mass;
+
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+    delete [] map;
+    delete [] fmap;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAMC::compute(int eflag, int vflag)
+{
+  int i,j,ii,n,inum_half,errorflag;
+  int *ilist_half,*numneigh_half,**firstneigh_half;
+  int *numneigh_full,**firstneigh_full;
+
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = eflag_global = vflag_global =
+         eflag_atom = vflag_atom = 0;
+
+  // grow local arrays if necessary
+
+  if (atom->nmax > nmax) {
+    memory->destroy(rho);
+    memory->destroy(rho0);
+    memory->destroy(rho1);
+    memory->destroy(rho2);
+    memory->destroy(rho3);
+    memory->destroy(frhop);
+    memory->destroy(gamma);
+    memory->destroy(dgamma1);
+    memory->destroy(dgamma2);
+    memory->destroy(dgamma3);
+    memory->destroy(arho2b);
+    memory->destroy(arho1);
+    memory->destroy(arho2);
+    memory->destroy(arho3);
+    memory->destroy(arho3b);
+    memory->destroy(t_ave);
+    memory->destroy(tsq_ave);
+
+    nmax = atom->nmax;
+
+    memory->create(rho,nmax,"pair:rho");
+    memory->create(rho0,nmax,"pair:rho0");
+    memory->create(rho1,nmax,"pair:rho1");
+    memory->create(rho2,nmax,"pair:rho2");
+    memory->create(rho3,nmax,"pair:rho3");
+    memory->create(frhop,nmax,"pair:frhop");
+    memory->create(gamma,nmax,"pair:gamma");
+    memory->create(dgamma1,nmax,"pair:dgamma1");
+    memory->create(dgamma2,nmax,"pair:dgamma2");
+    memory->create(dgamma3,nmax,"pair:dgamma3");
+    memory->create(arho2b,nmax,"pair:arho2b");
+    memory->create(arho1,nmax,3,"pair:arho1");
+    memory->create(arho2,nmax,6,"pair:arho2");
+    memory->create(arho3,nmax,10,"pair:arho3");
+    memory->create(arho3b,nmax,3,"pair:arho3b");
+    memory->create(t_ave,nmax,3,"pair:t_ave");
+    memory->create(tsq_ave,nmax,3,"pair:tsq_ave");
+  }
+
+  // neighbor list info
+
+  inum_half = listhalf->inum;
+  ilist_half = listhalf->ilist;
+  numneigh_half = listhalf->numneigh;
+  firstneigh_half = listhalf->firstneigh;
+  numneigh_full = listfull->numneigh;
+  firstneigh_full = listfull->firstneigh;
+
+  // strip neighbor lists of any special bond flags before using with MEAM
+  // necessary before doing neigh_f2c and neigh_c2f conversions each step
+
+  if (neighbor->ago == 0) {
+    neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half);
+    neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full);
+  }
+
+  // check size of scrfcn based on half neighbor list
+
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+
+  n = 0;
+  for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
+
+  if (n > maxneigh) {
+    memory->destroy(scrfcn);
+    memory->destroy(dscrfcn);
+    memory->destroy(fcpair);
+    maxneigh = n;
+    memory->create(scrfcn,maxneigh,"pair:scrfcn");
+    memory->create(dscrfcn,maxneigh,"pair:dscrfcn");
+    memory->create(fcpair,maxneigh,"pair:fcpair");
+  }
+
+  // zero out local arrays
+
+  for (i = 0; i < nall; i++) {
+    rho0[i] = 0.0;
+    arho2b[i] = 0.0;
+    arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
+    for (j = 0; j < 6; j++) arho2[i][j] = 0.0;
+    for (j = 0; j < 10; j++) arho3[i][j] = 0.0;
+    arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
+    t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
+    tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
+  }
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int ntype = atom->ntypes;
+
+  // change neighbor list indices to Fortran indexing
+
+  neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half);
+  neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full);
+
+  // 3 stages of MEAM calculation
+  // loop over my atoms followed by communication
+
+  int ifort;
+  int offset = 0;
+  errorflag = 0;
+
+  for (ii = 0; ii < inum_half; ii++) {
+    i = ilist_half[ii];
+    ifort = i+1;
+    meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0],
+                    &numneigh_half[i],firstneigh_half[i],
+                    &numneigh_full[i],firstneigh_full[i],
+                    &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
+                    rho0,&arho1[0][0],&arho2[0][0],arho2b,
+                    &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],
+                    &errorflag);
+    if (errorflag) {
+      char str[128];
+      sprintf(str,"MEAM library error %d",errorflag);
+      error->one(FLERR,str);
+    }
+    offset += numneigh_half[i];
+  }
+
+  comm->reverse_comm_pair(this);
+
+  meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
+                   &eng_vdwl,eatom,&ntype,type,fmap,
+                   &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
+                   &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
+                   dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
+  if (errorflag) {
+    char str[128];
+    sprintf(str,"MEAM library error %d",errorflag);
+    error->one(FLERR,str);
+  }
+
+  comm->forward_comm_pair(this);
+
+  offset = 0;
+
+  // vptr is first value in vatom if it will be used by meam_force()
+  // else vatom may not exist, so pass dummy ptr
+
+  double *vptr;
+  if (vflag_atom) vptr = &vatom[0][0];
+  else vptr = &cutmax;
+
+  for (ii = 0; ii < inum_half; ii++) {
+    i = ilist_half[ii];
+    ifort = i+1;
+    meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom,
+                &vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
+                &numneigh_half[i],firstneigh_half[i],
+                &numneigh_full[i],firstneigh_full[i],
+                &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
+                dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
+                &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
+                &t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag);
+    if (errorflag) {
+      char str[128];
+      sprintf(str,"MEAM library error %d",errorflag);
+      error->one(FLERR,str);
+    }
+    offset += numneigh_half[i];
+  }
+
+  // change neighbor list indices back to C indexing
+
+  neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half);
+  neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full);
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAMC::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+  map = new int[n+1];
+  fmap = new int[n];
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairMEAMC::settings(int narg, char **arg)
+{
+  if (narg != 0) error->all(FLERR,"Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairMEAMC::coeff(int narg, char **arg)
+{
+  int i,j,m,n;
+
+  if (!allocated) allocate();
+
+  if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients");
+
+  // insure I,J args are * *
+
+  if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
+    error->all(FLERR,"Incorrect args for pair coefficients");
+
+  // read MEAM element names between 2 filenames
+  // nelements = # of MEAM elements
+  // elements = list of unique element names
+
+  if (nelements) {
+    for (i = 0; i < nelements; i++) delete [] elements[i];
+    delete [] elements;
+    delete [] mass;
+  }
+  nelements = narg - 4 - atom->ntypes;
+  if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
+  elements = new char*[nelements];
+  mass = new double[nelements];
+
+  for (i = 0; i < nelements; i++) {
+    n = strlen(arg[i+3]) + 1;
+    elements[i] = new char[n];
+    strcpy(elements[i],arg[i+3]);
+  }
+
+  // read MEAM library and parameter files
+  // pass all parameters to MEAM package
+  // tell MEAM package that setup is done
+
+  read_files(arg[2],arg[2+nelements+1]);
+  meam_setup_done_(&cutmax);
+
+  // read args that map atom types to MEAM elements
+  // map[i] = which element the Ith atom type is, -1 if not mapped
+
+  for (i = 4 + nelements; i < narg; i++) {
+    m = i - (4+nelements) + 1;
+    for (j = 0; j < nelements; j++)
+      if (strcmp(arg[i],elements[j]) == 0) break;
+    if (j < nelements) map[m] = j;
+    else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
+    else error->all(FLERR,"Incorrect args for pair coefficients");
+  }
+
+  // clear setflag since coeff() called once with I,J = * *
+
+  n = atom->ntypes;
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  // set setflag i,j for type pairs where both are mapped to elements
+  // set mass for i,i in atom class
+
+  int count = 0;
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      if (map[i] >= 0 && map[j] >= 0) {
+        setflag[i][j] = 1;
+        if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);
+        count++;
+      }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairMEAMC::init_style()
+{
+  if (force->newton_pair == 0)
+    error->all(FLERR,"Pair style MEAM requires newton pair on");
+
+  // need full and half neighbor list
+
+  int irequest_full = neighbor->request(this,instance_me);
+  neighbor->requests[irequest_full]->id = 1;
+  neighbor->requests[irequest_full]->half = 0;
+  neighbor->requests[irequest_full]->full = 1;
+  int irequest_half = neighbor->request(this,instance_me);
+  neighbor->requests[irequest_half]->id = 2;
+
+  // setup Fortran-style mapping array needed by MEAM package
+  // fmap is indexed from 1:ntypes by Fortran and stores a Fortran index
+  // if type I is not a MEAM atom, fmap stores a 0
+
+  for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1;
+}
+
+/* ----------------------------------------------------------------------
+   neighbor callback to inform pair style of neighbor list to use
+   half or full
+------------------------------------------------------------------------- */
+
+void PairMEAMC::init_list(int id, NeighList *ptr)
+{
+  if (id == 1) listfull = ptr;
+  else if (id == 2) listhalf = ptr;
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairMEAMC::init_one(int i, int j)
+{
+  return cutmax;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAMC::read_files(char *globalfile, char *userfile)
+{
+  // open global meamf file on proc 0
+
+  FILE *fp;
+  if (comm->me == 0) {
+    fp = force->open_potential(globalfile);
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open MEAM potential file %s",globalfile);
+      error->one(FLERR,str);
+    }
+  }
+
+  // allocate parameter arrays
+
+  int params_per_line = 19;
+
+  int *lat = new int[nelements];
+  int *ielement = new int[nelements];
+  int *ibar = new int[nelements];
+  double *z = new double[nelements];
+  double *atwt = new double[nelements];
+  double *alpha = new double[nelements];
+  double *b0 = new double[nelements];
+  double *b1 = new double[nelements];
+  double *b2 = new double[nelements];
+  double *b3 = new double[nelements];
+  double *alat = new double[nelements];
+  double *esub = new double[nelements];
+  double *asub = new double[nelements];
+  double *t0 = new double[nelements];
+  double *t1 = new double[nelements];
+  double *t2 = new double[nelements];
+  double *t3 = new double[nelements];
+  double *rozero = new double[nelements];
+
+  bool *found = new bool[nelements];
+  for (int i = 0; i < nelements; i++) found[i] = false;
+
+  // read each set of params from global MEAM file
+  // one set of params can span multiple lines
+  // store params if element name is in element list
+  // if element name appears multiple times, only store 1st entry
+
+  int i,n,nwords;
+  char **words = new char*[params_per_line+1];
+  char line[MAXLINE],*ptr;
+  int eof = 0;
+
+  int nset = 0;
+  while (1) {
+    if (comm->me == 0) {
+      ptr = fgets(line,MAXLINE,fp);
+      if (ptr == NULL) {
+        eof = 1;
+        fclose(fp);
+      } else n = strlen(line) + 1;
+    }
+    MPI_Bcast(&eof,1,MPI_INT,0,world);
+    if (eof) break;
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    MPI_Bcast(line,n,MPI_CHAR,0,world);
+
+    // strip comment, skip line if blank
+
+    if ((ptr = strchr(line,'#'))) *ptr = '\0';
+    nwords = atom->count_words(line);
+    if (nwords == 0) continue;
+
+    // concatenate additional lines until have params_per_line words
+
+    while (nwords < params_per_line) {
+      n = strlen(line);
+      if (comm->me == 0) {
+        ptr = fgets(&line[n],MAXLINE-n,fp);
+        if (ptr == NULL) {
+          eof = 1;
+          fclose(fp);
+        } else n = strlen(line) + 1;
+      }
+      MPI_Bcast(&eof,1,MPI_INT,0,world);
+      if (eof) break;
+      MPI_Bcast(&n,1,MPI_INT,0,world);
+      MPI_Bcast(line,n,MPI_CHAR,0,world);
+      if ((ptr = strchr(line,'#'))) *ptr = '\0';
+      nwords = atom->count_words(line);
+    }
+
+    if (nwords != params_per_line)
+      error->all(FLERR,"Incorrect format in MEAM potential file");
+
+    // words = ptrs to all words in line
+    // strip single and double quotes from words
+
+    nwords = 0;
+    words[nwords++] = strtok(line,"' \t\n\r\f");
+    while ((words[nwords++] = strtok(NULL,"' \t\n\r\f"))) continue;
+
+    // skip if element name isn't in element list
+
+    for (i = 0; i < nelements; i++)
+      if (strcmp(words[0],elements[i]) == 0) break;
+    if (i >= nelements) continue;
+
+    // skip if element already appeared
+
+    if (found[i] == true) continue;
+    found[i] = true;
+
+    // map lat string to an integer
+
+    if (strcmp(words[1],"fcc") == 0) lat[i] = FCC;
+    else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC;
+    else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP;
+    else if (strcmp(words[1],"dim") == 0) lat[i] = DIM;
+    else if (strcmp(words[1],"dia") == 0) lat[i] = DIAMOND;
+    else error->all(FLERR,"Unrecognized lattice type in MEAM file 1");
+
+    // store parameters
+
+    z[i] = atof(words[2]);
+    ielement[i] = atoi(words[3]);
+    atwt[i] = atof(words[4]);
+    alpha[i] = atof(words[5]);
+    b0[i] = atof(words[6]);
+    b1[i] = atof(words[7]);
+    b2[i] = atof(words[8]);
+    b3[i] = atof(words[9]);
+    alat[i] = atof(words[10]);
+    esub[i] = atof(words[11]);
+    asub[i] = atof(words[12]);
+    t0[i] = atof(words[13]);
+    t1[i] = atof(words[14]);
+    t2[i] = atof(words[15]);
+    t3[i] = atof(words[16]);
+    rozero[i] = atof(words[17]);
+    ibar[i] = atoi(words[18]);
+
+    nset++;
+  }
+
+  // error if didn't find all elements in file
+
+  if (nset != nelements)
+    error->all(FLERR,"Did not find all elements in MEAM library file");
+
+  // pass element parameters to MEAM package
+
+  meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
+                       alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
+
+  // set element masses
+
+  for (i = 0; i < nelements; i++) mass[i] = atwt[i];
+
+  // clean-up memory
+
+  delete [] words;
+
+  delete [] lat;
+  delete [] ielement;
+  delete [] ibar;
+  delete [] z;
+  delete [] atwt;
+  delete [] alpha;
+  delete [] b0;
+  delete [] b1;
+  delete [] b2;
+  delete [] b3;
+  delete [] alat;
+  delete [] esub;
+  delete [] asub;
+  delete [] t0;
+  delete [] t1;
+  delete [] t2;
+  delete [] t3;
+  delete [] rozero;
+  delete [] found;
+
+  // done if user param file is NULL
+
+  if (strcmp(userfile,"NULL") == 0) return;
+
+  // open user param file on proc 0
+
+  if (comm->me == 0) {
+    fp = force->open_potential(userfile);
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open MEAM potential file %s",userfile);
+      error->one(FLERR,str);
+    }
+  }
+
+  // read settings
+  // pass them one at a time to MEAM package
+  // match strings to list of corresponding ints
+
+  int which;
+  double value;
+  int nindex,index[3];
+  int maxparams = 6;
+  char **params = new char*[maxparams];
+  int nparams;
+
+  eof = 0;
+  while (1) {
+    if (comm->me == 0) {
+      ptr = fgets(line,MAXLINE,fp);
+      if (ptr == NULL) {
+        eof = 1;
+        fclose(fp);
+      } else n = strlen(line) + 1;
+    }
+    MPI_Bcast(&eof,1,MPI_INT,0,world);
+    if (eof) break;
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    MPI_Bcast(line,n,MPI_CHAR,0,world);
+
+    // strip comment, skip line if blank
+
+    if ((ptr = strchr(line,'#'))) *ptr = '\0';
+    nparams = atom->count_words(line);
+    if (nparams == 0) continue;
+
+    // words = ptrs to all words in line
+
+    nparams = 0;
+    params[nparams++] = strtok(line,"=(), '\t\n\r\f");
+    while (nparams < maxparams &&
+           (params[nparams++] = strtok(NULL,"=(), '\t\n\r\f")))
+      continue;
+    nparams--;
+
+    for (which = 0; which < nkeywords; which++)
+      if (strcmp(params[0],keywords[which]) == 0) break;
+    if (which == nkeywords) {
+      char str[128];
+      sprintf(str,"Keyword %s in MEAM parameter file not recognized",
+              params[0]);
+      error->all(FLERR,str);
+    }
+    nindex = nparams - 2;
+    for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]);
+
+    // map lattce_meam value to an integer
+
+    if (which == 4) {
+      if (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 = DIAMOND;
+      else if (strcmp(params[nparams-1],"b1")  == 0) value = B1;
+      else if (strcmp(params[nparams-1],"c11") == 0) value = C11;
+      else if (strcmp(params[nparams-1],"l12") == 0) value = L12;
+      else if (strcmp(params[nparams-1],"b2")  == 0) value = B2;
+      else error->all(FLERR,"Unrecognized lattice type in MEAM file 2");
+    }
+    else value = atof(params[nparams-1]);
+
+    // pass single setting to MEAM package
+
+    int errorflag = 0;
+    meam_setup_param_(&which,&value,&nindex,index,&errorflag);
+    if (errorflag) {
+      char str[128];
+      sprintf(str,"MEAM library error %d",errorflag);
+      error->all(FLERR,str);
+    }
+  }
+
+  delete [] params;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairMEAMC::pack_forward_comm(int n, int *list, double *buf,
+                                int pbc_flag, int *pbc)
+{
+  int i,j,k,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    buf[m++] = rho0[j];
+    buf[m++] = rho1[j];
+    buf[m++] = rho2[j];
+    buf[m++] = rho3[j];
+    buf[m++] = frhop[j];
+    buf[m++] = gamma[j];
+    buf[m++] = dgamma1[j];
+    buf[m++] = dgamma2[j];
+    buf[m++] = dgamma3[j];
+    buf[m++] = arho2b[j];
+    buf[m++] = arho1[j][0];
+    buf[m++] = arho1[j][1];
+    buf[m++] = arho1[j][2];
+    buf[m++] = arho2[j][0];
+    buf[m++] = arho2[j][1];
+    buf[m++] = arho2[j][2];
+    buf[m++] = arho2[j][3];
+    buf[m++] = arho2[j][4];
+    buf[m++] = arho2[j][5];
+    for (k = 0; k < 10; k++) buf[m++] = arho3[j][k];
+    buf[m++] = arho3b[j][0];
+    buf[m++] = arho3b[j][1];
+    buf[m++] = arho3b[j][2];
+    buf[m++] = t_ave[j][0];
+    buf[m++] = t_ave[j][1];
+    buf[m++] = t_ave[j][2];
+    buf[m++] = tsq_ave[j][0];
+    buf[m++] = tsq_ave[j][1];
+    buf[m++] = tsq_ave[j][2];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAMC::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i,k,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    rho0[i] = buf[m++];
+    rho1[i] = buf[m++];
+    rho2[i] = buf[m++];
+    rho3[i] = buf[m++];
+    frhop[i] = buf[m++];
+    gamma[i] = buf[m++];
+    dgamma1[i] = buf[m++];
+    dgamma2[i] = buf[m++];
+    dgamma3[i] = buf[m++];
+    arho2b[i] = buf[m++];
+    arho1[i][0] = buf[m++];
+    arho1[i][1] = buf[m++];
+    arho1[i][2] = buf[m++];
+    arho2[i][0] = buf[m++];
+    arho2[i][1] = buf[m++];
+    arho2[i][2] = buf[m++];
+    arho2[i][3] = buf[m++];
+    arho2[i][4] = buf[m++];
+    arho2[i][5] = buf[m++];
+    for (k = 0; k < 10; k++) arho3[i][k] = buf[m++];
+    arho3b[i][0] = buf[m++];
+    arho3b[i][1] = buf[m++];
+    arho3b[i][2] = buf[m++];
+    t_ave[i][0] = buf[m++];
+    t_ave[i][1] = buf[m++];
+    t_ave[i][2] = buf[m++];
+    tsq_ave[i][0] = buf[m++];
+    tsq_ave[i][1] = buf[m++];
+    tsq_ave[i][2] = buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairMEAMC::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,k,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    buf[m++] = rho0[i];
+    buf[m++] = arho2b[i];
+    buf[m++] = arho1[i][0];
+    buf[m++] = arho1[i][1];
+    buf[m++] = arho1[i][2];
+    buf[m++] = arho2[i][0];
+    buf[m++] = arho2[i][1];
+    buf[m++] = arho2[i][2];
+    buf[m++] = arho2[i][3];
+    buf[m++] = arho2[i][4];
+    buf[m++] = arho2[i][5];
+    for (k = 0; k < 10; k++) buf[m++] = arho3[i][k];
+    buf[m++] = arho3b[i][0];
+    buf[m++] = arho3b[i][1];
+    buf[m++] = arho3b[i][2];
+    buf[m++] = t_ave[i][0];
+    buf[m++] = t_ave[i][1];
+    buf[m++] = t_ave[i][2];
+    buf[m++] = tsq_ave[i][0];
+    buf[m++] = tsq_ave[i][1];
+    buf[m++] = tsq_ave[i][2];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAMC::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,j,k,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    rho0[j] += buf[m++];
+    arho2b[j] += buf[m++];
+    arho1[j][0] += buf[m++];
+    arho1[j][1] += buf[m++];
+    arho1[j][2] += buf[m++];
+    arho2[j][0] += buf[m++];
+    arho2[j][1] += buf[m++];
+    arho2[j][2] += buf[m++];
+    arho2[j][3] += buf[m++];
+    arho2[j][4] += buf[m++];
+    arho2[j][5] += buf[m++];
+    for (k = 0; k < 10; k++) arho3[j][k] += buf[m++];
+    arho3b[j][0] += buf[m++];
+    arho3b[j][1] += buf[m++];
+    arho3b[j][2] += buf[m++];
+    t_ave[j][0] += buf[m++];
+    t_ave[j][1] += buf[m++];
+    t_ave[j][2] += buf[m++];
+    tsq_ave[j][0] += buf[m++];
+    tsq_ave[j][1] += buf[m++];
+    tsq_ave[j][2] += buf[m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double PairMEAMC::memory_usage()
+{
+  double bytes = 11 * nmax * sizeof(double);
+  bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double);
+  bytes += 3 * maxneigh * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   strip special bond flags from neighbor list entries
+   are not used with MEAM
+   need to do here so Fortran lib doesn't see them
+   done once per reneighbor so that neigh_f2c and neigh_c2f don't see them
+------------------------------------------------------------------------- */
+
+void PairMEAMC::neigh_strip(int inum, int *ilist,
+                           int *numneigh, int **firstneigh)
+{
+  int i,j,ii,jnum;
+  int *jlist;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   toggle neighbor list indices between zero- and one-based values
+   needed for access by MEAM Fortran library
+------------------------------------------------------------------------- */
+
+void PairMEAMC::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh)
+{
+  int i,j,ii,jnum;
+  int *jlist;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    for (j = 0; j < jnum; j++) jlist[j]--;
+  }
+}
+
+void PairMEAMC::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh)
+{
+  int i,j,ii,jnum;
+  int *jlist;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    for (j = 0; j < jnum; j++) jlist[j]++;
+  }
+}

src/USER-MEAMC/pair_meamc.h

@@ -0,0 +1,155 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(meam/c,PairMEAMC)
+
+#else
+
+#ifndef LMP_PAIR_MEAMC_H
+#define LMP_PAIR_MEAMC_H
+
+extern "C" {
+  void meam_setup_global_(int *, int *, double *, int *, double *, double *,
+                         double *, double *, double *, double *, double *,
+                         double *, double *, double *, double *, double *,
+                         double *, double *, int *);
+  void meam_setup_param_(int *, double *, int *, int *, int *);
+  void meam_setup_done_(double *);
+
+  void meam_dens_init_(int *, int *, int *, int *, int *,
+                       double *, int *, int *, int *, int *,
+                       double *, double *, double *, double *,
+                       double *, double *,
+                       double *, double *, double *, double *, double *,
+                       int *);
+
+  void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
+                        int *, int *, int *,
+                        double *, double *, double *, double *,
+                        double *, double *, double *,
+                        double *, double *, double *, double *,
+                        double *, double *,
+                        double *, double *, double *, double *, int *);
+
+  void meam_force_(int *, int *, int *, int *, int *, int *,
+                   double *, double *, int *, int *, int *,
+                   double *, int *, int *, int *, int *, double *, double *,
+                   double *, double *, double *, double *, double *, double *,
+                   double *, double *, double *, double *, double *, double *,
+                   double *, double *, double *, double *, double *, double *, int *);
+
+  void meam_cleanup_();
+}
+
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairMEAMC : public Pair {
+ public:
+  PairMEAMC(class LAMMPS *);
+  ~PairMEAMC();
+  void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  double init_one(int, int);
+
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+  double memory_usage();
+
+ private:
+  double cutmax;                // max cutoff for all elements
+  int nelements;                // # of unique elements
+  char **elements;              // names of unique elements
+  double *mass;                 // mass of each element
+
+  int *map;                     // mapping from atom types to elements
+  int *fmap;                    // Fortran version of map array for MEAM lib
+
+  int maxneigh;
+  double *scrfcn,*dscrfcn,*fcpair;
+
+  int nmax;
+  double *rho,*rho0,*rho1,*rho2,*rho3,*frhop;
+  double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b;
+  double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave;
+
+  void allocate();
+  void read_files(char *, char *);
+  void neigh_strip(int, int *, int *, int **);
+  void neigh_f2c(int, int *, int *, int **);
+  void neigh_c2f(int, int *, int *, int **);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: MEAM library error %d
+
+A call to the MEAM Fortran library returned an error.
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Pair style MEAM requires newton pair on
+
+See the newton command.  This is a restriction to use the MEAM
+potential.
+
+E: Cannot open MEAM potential file %s
+
+The specified MEAM potential file cannot be opened.  Check that the
+path and name are correct.
+
+E: Incorrect format in MEAM potential file
+
+Incorrect number of words per line in the potential file.
+
+E: Unrecognized lattice type in MEAM file 1
+
+The lattice type in an entry of the MEAM library file is not
+valid.
+
+E: Did not find all elements in MEAM library file
+
+The requested elements were not found in the MEAM file.
+
+E: Keyword %s in MEAM parameter file not recognized
+
+Self-explanatory.
+
+E: Unrecognized lattice type in MEAM file 2
+
+The lattice type in an entry of the MEAM parameter file is not
+valid.
+
+*/