From bab292b551ba7b65f205c52c7570e7c2948cab93 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= <sebastian.huetter@ovgu.de>
Date: Fri, 9 Jun 2017 17:11:29 +0200
Subject: [PATCH] Create package USER-MEAMC
Step 1: very literal translation of lib/meam
---
examples/meam/in.meamc | 30 +
examples/meam/in.meamc.shear | 79 +++
src/Makefile | 2 +-
src/USER-MEAMC/Install.sh | 47 ++
src/USER-MEAMC/README | 26 +
src/USER-MEAMC/fm_exp.c | 133 ++++
src/USER-MEAMC/meam.h | 157 +++++
src/USER-MEAMC/meam_cleanup.c | 14 +
src/USER-MEAMC/meam_data.c | 3 +
src/USER-MEAMC/meam_dens_final.c | 262 ++++++++
src/USER-MEAMC/meam_dens_init.c | 504 +++++++++++++++
src/USER-MEAMC/meam_force.c | 551 +++++++++++++++++
src/USER-MEAMC/meam_setup_done.c | 945 +++++++++++++++++++++++++++++
src/USER-MEAMC/meam_setup_global.c | 99 +++
src/USER-MEAMC/meam_setup_param.c | 223 +++++++
src/USER-MEAMC/pair_meamc.cpp | 945 +++++++++++++++++++++++++++++
src/USER-MEAMC/pair_meamc.h | 155 +++++
17 files changed, 4174 insertions(+), 1 deletion(-)
create mode 100644 examples/meam/in.meamc
create mode 100644 examples/meam/in.meamc.shear
create mode 100644 src/USER-MEAMC/Install.sh
create mode 100644 src/USER-MEAMC/README
create mode 100644 src/USER-MEAMC/fm_exp.c
create mode 100644 src/USER-MEAMC/meam.h
create mode 100644 src/USER-MEAMC/meam_cleanup.c
create mode 100644 src/USER-MEAMC/meam_data.c
create mode 100644 src/USER-MEAMC/meam_dens_final.c
create mode 100644 src/USER-MEAMC/meam_dens_init.c
create mode 100644 src/USER-MEAMC/meam_force.c
create mode 100644 src/USER-MEAMC/meam_setup_done.c
create mode 100644 src/USER-MEAMC/meam_setup_global.c
create mode 100644 src/USER-MEAMC/meam_setup_param.c
create mode 100644 src/USER-MEAMC/pair_meamc.cpp
create mode 100644 src/USER-MEAMC/pair_meamc.h
diff --git a/examples/meam/in.meamc b/examples/meam/in.meamc
new file mode 100644
index 0000000000..f6815cd7d4
--- /dev/null
+++ b/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
diff --git a/examples/meam/in.meamc.shear b/examples/meam/in.meamc.shear
new file mode 100644
index 0000000000..e4584d9744
--- /dev/null
+++ b/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
diff --git a/src/Makefile b/src/Makefile
index 92a430a747..c7b20dcb13 100644
--- a/src/Makefile
+++ b/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
diff --git a/src/USER-MEAMC/Install.sh b/src/USER-MEAMC/Install.sh
new file mode 100644
index 0000000000..2b4098a473
--- /dev/null
+++ b/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
diff --git a/src/USER-MEAMC/README b/src/USER-MEAMC/README
new file mode 100644
index 0000000000..8d6bd6360c
--- /dev/null
+++ b/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.
+
diff --git a/src/USER-MEAMC/fm_exp.c b/src/USER-MEAMC/fm_exp.c
new file mode 100644
index 0000000000..00b957b27c
--- /dev/null
+++ b/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:
+ */
diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h
new file mode 100644
index 0000000000..ab2c756dc3
--- /dev/null
+++ b/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)]
+
+
+
+
+
diff --git a/src/USER-MEAMC/meam_cleanup.c b/src/USER-MEAMC/meam_cleanup.c
new file mode 100644
index 0000000000..5c81fc3b49
--- /dev/null
+++ b/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);
+
+}
diff --git a/src/USER-MEAMC/meam_data.c b/src/USER-MEAMC/meam_data.c
new file mode 100644
index 0000000000..54ddf81b64
--- /dev/null
+++ b/src/USER-MEAMC/meam_data.c
@@ -0,0 +1,3 @@
+#include "meam.h"
+
+meam_data_t meam_data = {};
\ No newline at end of file
diff --git a/src/USER-MEAMC/meam_dens_final.c b/src/USER-MEAMC/meam_dens_final.c
new file mode 100644
index 0000000000..1e553866e4
--- /dev/null
+++ b/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
diff --git a/src/USER-MEAMC/meam_dens_init.c b/src/USER-MEAMC/meam_dens_init.c
new file mode 100644
index 0000000000..8c1d84f6d4
--- /dev/null
+++ b/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
+
+
+
+
diff --git a/src/USER-MEAMC/meam_force.c b/src/USER-MEAMC/meam_force.c
new file mode 100644
index 0000000000..d27d8ca7eb
--- /dev/null
+++ b/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
+ }
+
+ }
diff --git a/src/USER-MEAMC/meam_setup_done.c b/src/USER-MEAMC/meam_setup_done.c
new file mode 100644
index 0000000000..e777025067
--- /dev/null
+++ b/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;
+ }
diff --git a/src/USER-MEAMC/meam_setup_global.c b/src/USER-MEAMC/meam_setup_global.c
new file mode 100644
index 0000000000..c8c4ef8468
--- /dev/null
+++ b/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;
+
+ }
+
+
diff --git a/src/USER-MEAMC/meam_setup_param.c b/src/USER-MEAMC/meam_setup_param.c
new file mode 100644
index 0000000000..2c49481a28
--- /dev/null
+++ b/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;
+ }
+ }
diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp
new file mode 100644
index 0000000000..f2c7da9239
--- /dev/null
+++ b/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]++;
+ }
+}
diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h
new file mode 100644
index 0000000000..c4e03aef55
--- /dev/null
+++ b/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.
+
+*/
--
GitLab