From 03c93b31d61a961de1385c63f314e59683e3d948 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= Date: Fri, 9 Jun 2017 19:44:59 +0200 Subject: [PATCH] Convert to C++, allow multiple instances --- src/USER-MEAMC/README | 2 +- src/USER-MEAMC/meam.h | 288 +++++++----- src/USER-MEAMC/meam_cleanup.cpp | 20 +- src/USER-MEAMC/meam_data.cpp | 5 - src/USER-MEAMC/meam_dens_final.cpp | 80 ++-- src/USER-MEAMC/meam_dens_init.cpp | 132 +++--- src/USER-MEAMC/meam_force.cpp | 198 ++++----- src/USER-MEAMC/meam_setup_done.cpp | 634 +++++++++++++-------------- src/USER-MEAMC/meam_setup_global.cpp | 100 +++-- src/USER-MEAMC/meam_setup_param.cpp | 64 ++- src/USER-MEAMC/pair_meamc.cpp | 19 +- src/USER-MEAMC/pair_meamc.h | 36 +- 12 files changed, 762 insertions(+), 816 deletions(-) delete mode 100644 src/USER-MEAMC/meam_data.cpp diff --git a/src/USER-MEAMC/README b/src/USER-MEAMC/README index 8d6bd6360c..19756d98a3 100644 --- a/src/USER-MEAMC/README +++ b/src/USER-MEAMC/README @@ -2,7 +2,7 @@ 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 +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. diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index c8e46903c5..3ef6bf8a42 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -1,117 +1,176 @@ +#ifndef LMP_MEAM_H +#define LMP_MEAM_H + #include - #define maxelt 5 - double fm_exp(double); - - typedef enum {FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2} lattice_t; +#define maxelt 5 -typedef struct { +extern "C" { + double fm_exp(double); +} + +typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t; + +typedef struct +{ int dim1, dim2; - double *ptr; + double* ptr; } allocatable_double_2d; - -typedef struct { -// cutforce = force cutoff -// cutforcesq = force cutoff squared - double cutforce,cutforcesq; +class MEAM { +private: + // cutforce = force cutoff + // cutforcesq = force cutoff squared -// 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 + double cutforce, cutforcesq; -// nr,dr = pair function discretization parameters -// nrar,rdrar = spline coeff array parameters + // 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 - 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; + // nr,dr = pair function discretization parameters + // nrar,rdrar = spline coeff array parameters - allocatable_double_2d phir; // [:,:] + 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 phirar,phirar1,phirar2,phirar3,phirar4,phirar5,phirar6; // [:,:] + allocatable_double_2d phir; // [:,:] - double attrac_meam[maxelt+1][maxelt+1],repuls_meam[maxelt+1][maxelt+1]; + allocatable_double_2d phirar, phirar1, phirar2, phirar3, phirar4, phirar5, + phirar6; // [:,:] - 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]; + double attrac_meam[maxelt + 1][maxelt + 1], + repuls_meam[maxelt + 1][maxelt + 1]; - int nr,nrar; - double dr,rdrar; -} meam_data_t; + 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]; -extern meam_data_t meam_data; + int nr, nrar; + double dr, rdrar; +public: +void meam_checkindex(int, int, int, int*, int*); +void meam_setup_param_(int*, double*, int*, int*, 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 G_gam(double, int, double, double*, int*); +void dG_gam(double, int, double, double*, 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 getscreen(int, int, double*, double*, double*, double*, int, int*, int, int*, int, int*, int*); +void screen(int, int, int, double*, double, double*, int, int*, int, int*, int*); +void calc_rho1(int, int, int, int*, int*, double*, int, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*); +void dsij(int, int, int, int, int, int, double, double*, double*, int, int*, int*, double*, double*, double*); +void fcut(double, double*); +void dfcut(double, double*, double*); +void dCfunc(double, double, double, double*); +void dCfunc2(double, double, double, double*, double*); + +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_(); +void meam_setup_done_(double*); +void alloyparams(); +void compute_pair_meam(); +double phi_meam(double, int, int); +void compute_reference_density(); +void get_shpfcn(double*, lattice_t); +void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double, double, double, double, int, int, lattice_t); +void get_Zij(int*, lattice_t); +void get_Zij2(int*, double*, double*, lattice_t, double, double); +void get_sijk(double, int, int, int, double*); +void get_densref(double, int, int, double*, double*, double*, double*, double*, double*, double*, double*); +double zbl(double, int, int); +double erose(double, double, double, double, double, double, int); +void interpolate_meam(int); +double compute_phi(double, int, int); +void meam_setup_global_(int*, int*, double*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); + +}; // Functions we need for compat #ifndef max -#define max(a,b) \ - ({ __typeof__ (a) _a = (a); \ - __typeof__ (b) _b = (b); \ - _a > _b ? _a : _b; }) +#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; }) +#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 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;} +#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. @@ -124,34 +183,37 @@ Fortran Array Semantics in C. // 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;\ +#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)] +#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;}) +#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 meam_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)] - - - +#define arr2(arr, i, j) arr.ptr[(arr.dim1) * (j - 1) + (i - 1)] +#endif \ No newline at end of file diff --git a/src/USER-MEAMC/meam_cleanup.cpp b/src/USER-MEAMC/meam_cleanup.cpp index 3b1edc8cd9..b0ee013c4e 100644 --- a/src/USER-MEAMC/meam_cleanup.cpp +++ b/src/USER-MEAMC/meam_cleanup.cpp @@ -1,17 +1,15 @@ -extern "C" { #include "meam.h" void -meam_cleanup_(void) +MEAM::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); -} + meam_deallocate(this->phirar6); + meam_deallocate(this->phirar5); + meam_deallocate(this->phirar4); + meam_deallocate(this->phirar3); + meam_deallocate(this->phirar2); + meam_deallocate(this->phirar1); + meam_deallocate(this->phirar); + meam_deallocate(this->phir); } \ No newline at end of file diff --git a/src/USER-MEAMC/meam_data.cpp b/src/USER-MEAMC/meam_data.cpp deleted file mode 100644 index c5fd152147..0000000000 --- a/src/USER-MEAMC/meam_data.cpp +++ /dev/null @@ -1,5 +0,0 @@ -extern "C" { -#include "meam.h" - -meam_data_t meam_data = {}; -} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index 587d3bbb5f..8a63fbd820 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -1,15 +1,6 @@ -extern "C" { #include "meam.h" #include -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 *, @@ -28,7 +19,7 @@ void get_shpfcn(double* s /* s(3) */, lattice_t latt); // void -meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, +MEAM::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, @@ -67,23 +58,23 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, for (m = 1; m <= 6; m++) { arr1v(rho2, i) = arr1v(rho2, i) + - meam_data.v2D[m] * arr2v(Arho2, m, i) * arr2v(Arho2, m, i); + this->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); + this->v3D[m] * arr2v(Arho3, m, i) * arr2v(Arho3, m, i); } if (arr1v(rho0, i) > 0.0) { - if (meam_data.ialloy == 1) { + if (this->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 if (this->ialloy == 2) { + arr2v(t_ave, 1, i) = this->t1_meam[elti]; + arr2v(t_ave, 2, i) = this->t2_meam[elti]; + arr2v(t_ave, 3, i) = this->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); @@ -99,56 +90,56 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, arr1v(Gamma, i) = arr1v(Gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i)); } - Z = meam_data.Z_meam[elti]; + Z = this->Z_meam[elti]; - G_gam(arr1v(Gamma, i), meam_data.ibar_meam[elti], - meam_data.gsmooth_factor, &G, errorflag); + G_gam(arr1v(Gamma, i), this->ibar_meam[elti], + this->gsmooth_factor, &G, errorflag); if (*errorflag != 0) return; - get_shpfcn(shp, meam_data.lattce_meam[elti][elti]); - if (meam_data.ibar_meam[elti] <= 0) { + get_shpfcn(shp, this->lattce_meam[elti][elti]); + if (this->ibar_meam[elti] <= 0) { Gbar = 1.0; dGbar = 0.0; } else { - if (meam_data.mix_ref_t == 1) { + if (this->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]) / + gam = (this->t1_meam[elti] * shp[1] + + this->t2_meam[elti] * shp[2] + + this->t3_meam[elti] * shp[3]) / (Z * Z); } - G_gam(gam, meam_data.ibar_meam[elti], meam_data.gsmooth_factor, &Gbar, + G_gam(gam, this->ibar_meam[elti], this->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) { + if (this->mix_ref_t == 1) { + if (this->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, + dG_gam(gam, this->ibar_meam[elti], this->gsmooth_factor, &Gbar, &dGbar); } - rho_bkgd = meam_data.rho0_meam[elti] * Z * Gbar; + rho_bkgd = this->rho0_meam[elti] * Z * Gbar; } else { - if (meam_data.bkgd_dyn == 1) { - rho_bkgd = meam_data.rho0_meam[elti] * Z; + if (this->bkgd_dyn == 1) { + rho_bkgd = this->rho0_meam[elti] * Z; } else { - rho_bkgd = meam_data.rho_ref_meam[elti]; + rho_bkgd = this->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); + dG_gam(arr1v(Gamma, i), this->ibar_meam[elti], + this->gsmooth_factor, &G, &dG); arr1v(dGamma1, i) = (G - 2 * dG * arr1v(Gamma, i)) * denom; @@ -161,30 +152,30 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, // 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) { + if (this->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]; + B = this->A_meam[elti] * this->Ec_meam[elti][elti]; if (!iszero(rhob)) { - if (meam_data.emb_lin_neg == 1 && rhob <= 0) { + if (this->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) { + if (this->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) { + if (this->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); @@ -192,7 +183,7 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, } } } else { - if (meam_data.emb_lin_neg == 1) { + if (this->emb_lin_neg == 1) { arr1v(fp, i) = -B; } else { arr1v(fp, i) = B; @@ -205,7 +196,7 @@ meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* errorflag) +MEAM::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) @@ -246,7 +237,7 @@ G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* errorflag) // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* dG) +MEAM::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) @@ -289,4 +280,3 @@ dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* dG) } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index 4ea56ded9c..b94eec65cb 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -1,27 +1,6 @@ -extern "C" { #include "meam.h" #include -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 *, @@ -41,7 +20,7 @@ void dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, // void -meam_dens_init_(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x, +MEAM::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, @@ -62,7 +41,7 @@ meam_dens_init_(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x, // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, +MEAM::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) { @@ -81,7 +60,7 @@ getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, /*unused:double dC1a,dC1b,dC2a,dC2b;*/ double rnorm, fc, dfc, drinv; - drinv = 1.0 / meam_data.delr_meam; + drinv = 1.0 / this->delr_meam; elti = arr1v(fmap, arr1v(type, i)); if (elti > 0) { @@ -105,12 +84,12 @@ getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, delzij = zjtmp - zitmp; rij2 = delxij * delxij + delyij * delyij + delzij * delzij; rij = sqrt(rij2); - if (rij > meam_data.rc_meam) { + if (rij > this->rc_meam) { fcij = 0.0; dfcij = 0.0; sij = 0.0; } else { - rnorm = (meam_data.rc_meam - rij) * drinv; + rnorm = (this->rc_meam - rij) * drinv; screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full, ntype, type, fmap); dfcut(rnorm, &fc, &dfc); @@ -123,7 +102,7 @@ getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, sfcij = sij * fcij; if (iszero(sfcij) || iszero(sfcij - 1.0)) goto LABEL_100; - rbound = meam_data.ebound_meam[elti][eltj] * rij2; + rbound = this->ebound_meam[elti][eltj] * rij2; for (kn = 1; kn <= numneigh_full; kn++) { k = arr1v(firstneigh_full, kn); if (k == j) @@ -154,8 +133,8 @@ getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, 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]; + Cmax = this->Cmax_meam[elti][eltj][eltk]; + Cmin = this->Cmin_meam[elti][eltj][eltk]; if (cikj >= Cmax) { continue; // Note that cikj may be slightly negative (within numerical @@ -189,7 +168,7 @@ getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair, // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x, +MEAM::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) @@ -222,65 +201,65 @@ calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x, 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) { + if (rij2 < this->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]; + ai = rij / this->re_meam[elti][elti] - 1.0; + aj = rij / this->re_meam[eltj][eltj] - 1.0; + ro0i = this->rho0_meam[elti]; + ro0j = this->rho0_meam[eltj]; + rhoa0j = ro0j * fm_exp(-this->beta0_meam[eltj] * aj) * sij; + rhoa1j = ro0j * fm_exp(-this->beta1_meam[eltj] * aj) * sij; + rhoa2j = ro0j * fm_exp(-this->beta2_meam[eltj] * aj) * sij; + rhoa3j = ro0j * fm_exp(-this->beta3_meam[eltj] * aj) * sij; + rhoa0i = ro0i * fm_exp(-this->beta0_meam[elti] * ai) * sij; + rhoa1i = ro0i * fm_exp(-this->beta1_meam[elti] * ai) * sij; + rhoa2i = ro0i * fm_exp(-this->beta2_meam[elti] * ai) * sij; + rhoa3i = ro0i * fm_exp(-this->beta3_meam[elti] * ai) * sij; + if (this->ialloy == 1) { + rhoa1j = rhoa1j * this->t1_meam[eltj]; + rhoa2j = rhoa2j * this->t2_meam[eltj]; + rhoa3j = rhoa3j * this->t3_meam[eltj]; + rhoa1i = rhoa1i * this->t1_meam[elti]; + rhoa2i = rhoa2i * this->t2_meam[elti]; + rhoa3i = rhoa3i * this->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) { + if (this->ialloy != 2) { arr2v(t_ave, 1, i) = - arr2v(t_ave, 1, i) + meam_data.t1_meam[eltj] * rhoa0j; + arr2v(t_ave, 1, i) + this->t1_meam[eltj] * rhoa0j; arr2v(t_ave, 2, i) = - arr2v(t_ave, 2, i) + meam_data.t2_meam[eltj] * rhoa0j; + arr2v(t_ave, 2, i) + this->t2_meam[eltj] * rhoa0j; arr2v(t_ave, 3, i) = - arr2v(t_ave, 3, i) + meam_data.t3_meam[eltj] * rhoa0j; + arr2v(t_ave, 3, i) + this->t3_meam[eltj] * rhoa0j; arr2v(t_ave, 1, j) = - arr2v(t_ave, 1, j) + meam_data.t1_meam[elti] * rhoa0i; + arr2v(t_ave, 1, j) + this->t1_meam[elti] * rhoa0i; arr2v(t_ave, 2, j) = - arr2v(t_ave, 2, j) + meam_data.t2_meam[elti] * rhoa0i; + arr2v(t_ave, 2, j) + this->t2_meam[elti] * rhoa0i; arr2v(t_ave, 3, j) = - arr2v(t_ave, 3, j) + meam_data.t3_meam[elti] * rhoa0i; + arr2v(t_ave, 3, j) + this->t3_meam[elti] * rhoa0i; } - if (meam_data.ialloy == 1) { + if (this->ialloy == 1) { arr2v(tsq_ave, 1, i) = arr2v(tsq_ave, 1, i) + - meam_data.t1_meam[eltj] * meam_data.t1_meam[eltj] * rhoa0j; + this->t1_meam[eltj] * this->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; + this->t2_meam[eltj] * this->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; + this->t3_meam[eltj] * this->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; + this->t1_meam[elti] * this->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; + this->t2_meam[elti] * this->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; + this->t3_meam[elti] * this->t3_meam[elti] * rhoa0i; } arr1v(arho2b, i) = arr1v(arho2b, i) + rhoa2j; arr1v(arho2b, j) = arr1v(arho2b, j) + rhoa2i; @@ -321,7 +300,7 @@ calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x, // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -screen(int i, int j, int nmax, double* x, double rijsq, double* sij, +MEAM::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) @@ -345,7 +324,7 @@ screen(int i, int j, int nmax, double* x, double rijsq, double* sij, 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; + rbound = this->ebound_meam[elti][eltj] * rijsq; for (nk = 1; nk <= numneigh_full; nk++) { k = arr1v(firstneigh_full, nk); @@ -372,8 +351,8 @@ screen(int i, int j, int nmax, double* x, double rijsq, double* sij, 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]; + Cmax = this->Cmax_meam[elti][eltj][eltk]; + Cmin = this->Cmin_meam[elti][eltj][eltk]; if (cikj >= Cmax) continue; // note that cikj may be slightly negative (within numerical @@ -394,7 +373,7 @@ screen(int i, int j, int nmax, double* x, double rijsq, double* sij, // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2, +MEAM::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) { @@ -416,13 +395,13 @@ dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2, 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]; + Cmax = this->Cmax_meam[elti][eltj][eltk]; + Cmin = this->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]; + rbound = rij2 * this->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); @@ -456,7 +435,7 @@ dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2, // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -fcut(double xi, double* fc) +MEAM::fcut(double xi, double* fc) { // cutoff function double a; @@ -477,7 +456,7 @@ fcut(double xi, double* fc) // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -dfcut(double xi, double* fc, double* dfc) +MEAM::dfcut(double xi, double* fc, double* dfc) { // cutoff function and its derivative double a, a3, a4; @@ -501,7 +480,7 @@ dfcut(double xi, double* fc, double* dfc) // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -dCfunc(double rij2, double rik2, double rjk2, double* dCikj) +MEAM::dCfunc(double rij2, double rik2, double rjk2, double* dCikj) { // Inputs: rij,rij2,rik2,rjk2 // Outputs: dCikj = derivative of Cikj w.r.t. rij @@ -518,7 +497,7 @@ dCfunc(double rij2, double rik2, double rjk2, double* dCikj) // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void -dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, double* dCikj2) +MEAM::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 @@ -541,4 +520,3 @@ dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, double* dCikj2) } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index 8f2efa2daa..80e34c66d6 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -1,15 +1,6 @@ -extern "C" { #include "meam.h" #include -// 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 @@ -32,7 +23,7 @@ void dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2, // void -meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, +MEAM::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, @@ -113,27 +104,27 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, 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) { + if (rij2 < this->cutforcesq) { rij = sqrt(rij2); r = rij; // Compute phi and phip - ind = meam_data.eltind[elti][eltj]; - pp = rij * meam_data.rdrar + 1.0; + ind = this->eltind[elti][eltj]; + pp = rij * this->rdrar + 1.0; kk = (int)pp; - kk = min(kk, meam_data.nrar - 1); + kk = min(kk, this->nrar - 1); pp = pp - kk; pp = min(pp, 1.0); - phi = ((arr2(meam_data.phirar3, kk, ind) * pp + - arr2(meam_data.phirar2, kk, ind)) * + phi = ((arr2(this->phirar3, kk, ind) * pp + + arr2(this->phirar2, kk, ind)) * pp + - arr2(meam_data.phirar1, kk, ind)) * + arr2(this->phirar1, kk, ind)) * pp + - arr2(meam_data.phirar, kk, ind); - phip = (arr2(meam_data.phirar6, kk, ind) * pp + - arr2(meam_data.phirar5, kk, ind)) * + arr2(this->phirar, kk, ind); + phip = (arr2(this->phirar6, kk, ind) * pp + + arr2(this->phirar5, kk, ind)) * pp + - arr2(meam_data.phirar4, kk, ind); + arr2(this->phirar4, kk, ind); recip = 1.0 / r; if (*eflag_either != 0) { @@ -149,30 +140,30 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, // write(1,*) "force_meamf: phip: ",phip // Compute pair densities and derivatives - invrei = 1.0 / meam_data.re_meam[elti][elti]; + invrei = 1.0 / this->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; + ro0i = this->rho0_meam[elti]; + rhoa0i = ro0i * fm_exp(-this->beta0_meam[elti] * ai); + drhoa0i = -this->beta0_meam[elti] * invrei * rhoa0i; + rhoa1i = ro0i * fm_exp(-this->beta1_meam[elti] * ai); + drhoa1i = -this->beta1_meam[elti] * invrei * rhoa1i; + rhoa2i = ro0i * fm_exp(-this->beta2_meam[elti] * ai); + drhoa2i = -this->beta2_meam[elti] * invrei * rhoa2i; + rhoa3i = ro0i * fm_exp(-this->beta3_meam[elti] * ai); + drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i; if (elti != eltj) { - invrej = 1.0 / meam_data.re_meam[eltj][eltj]; + invrej = 1.0 / this->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; + ro0j = this->rho0_meam[eltj]; + rhoa0j = ro0j * fm_exp(-this->beta0_meam[eltj] * aj); + drhoa0j = -this->beta0_meam[eltj] * invrej * rhoa0j; + rhoa1j = ro0j * fm_exp(-this->beta1_meam[eltj] * aj); + drhoa1j = -this->beta1_meam[eltj] * invrej * rhoa1j; + rhoa2j = ro0j * fm_exp(-this->beta2_meam[eltj] * aj); + drhoa2j = -this->beta2_meam[eltj] * invrej * rhoa2j; + rhoa3j = ro0j * fm_exp(-this->beta3_meam[eltj] * aj); + drhoa3j = -this->beta3_meam[eltj] * invrej * rhoa3j; } else { rhoa0j = rhoa0i; drhoa0j = drhoa0i; @@ -184,19 +175,19 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, 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]; + if (this->ialloy == 1) { + rhoa1j = rhoa1j * this->t1_meam[eltj]; + rhoa2j = rhoa2j * this->t2_meam[eltj]; + rhoa3j = rhoa3j * this->t3_meam[eltj]; + rhoa1i = rhoa1i * this->t1_meam[elti]; + rhoa2i = rhoa2i * this->t2_meam[elti]; + rhoa3i = rhoa3i * this->t3_meam[elti]; + drhoa1j = drhoa1j * this->t1_meam[eltj]; + drhoa2j = drhoa2j * this->t2_meam[eltj]; + drhoa3j = drhoa3j * this->t3_meam[eltj]; + drhoa1i = drhoa1i * this->t1_meam[elti]; + drhoa2i = drhoa2i * this->t2_meam[elti]; + drhoa3i = drhoa3i * this->t3_meam[elti]; } nv2 = 1; @@ -212,12 +203,12 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, 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]; + arg = delij[n] * delij[p] * delij[q] * this->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]; + arg = delij[n] * delij[p] * this->v2D[nv2]; arg1i2 = arg1i2 + arr2v(Arho2, nv2, i) * arg; arg1j2 = arg1j2 + arr2v(Arho2, nv2, j) * arg; nv2 = nv2 + 1; @@ -254,9 +245,9 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, drho2drm2[m] = 0.0; for (n = 1; n <= 3; n++) { drho2drm1[m] = drho2drm1[m] + - arr2v(Arho2, meam_data.vind2D[m][n], i) * delij[n]; + arr2v(Arho2, this->vind2D[m][n], i) * delij[n]; drho2drm2[m] = drho2drm2[m] - - arr2v(Arho2, meam_data.vind2D[m][n], j) * delij[n]; + arr2v(Arho2, this->vind2D[m][n], j) * delij[n]; } drho2drm1[m] = a2 * rhoa2j * drho2drm1[m]; drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m]; @@ -278,11 +269,11 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, nv2 = 1; for (n = 1; n <= 3; n++) { for (p = n; p <= 3; p++) { - arg = delij[n] * delij[p] * meam_data.v2D[nv2]; + arg = delij[n] * delij[p] * this->v2D[nv2]; drho3drm1[m] = drho3drm1[m] + - arr2v(Arho3, meam_data.vind3D[m][n][p], i) * arg; + arr2v(Arho3, this->vind3D[m][n][p], i) * arg; drho3drm2[m] = drho3drm2[m] + - arr2v(Arho3, meam_data.vind3D[m][n][p], j) * arg; + arr2v(Arho3, this->vind3D[m][n][p], j) * arg; nv2 = nv2 + 1; } } @@ -300,7 +291,7 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, t2j = arr2v(t_ave, 2, j); t3j = arr2v(t_ave, 3, j); - if (meam_data.ialloy == 1) { + if (this->ialloy == 1) { a1i = 0.0; a1j = 0.0; @@ -321,20 +312,20 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, 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)); + dt1dr1 = a1i * (this->t1_meam[eltj] - + t1i * pow(this->t1_meam[eltj], 2)); + dt1dr2 = a1j * (this->t1_meam[elti] - + t1j * pow(this->t1_meam[elti], 2)); + dt2dr1 = a2i * (this->t2_meam[eltj] - + t2i * pow(this->t2_meam[eltj], 2)); + dt2dr2 = a2j * (this->t2_meam[elti] - + t2j * pow(this->t2_meam[elti], 2)); + dt3dr1 = a3i * (this->t3_meam[eltj] - + t3i * pow(this->t3_meam[eltj], 2)); + dt3dr2 = a3j * (this->t3_meam[elti] - + t3j * pow(this->t3_meam[elti], 2)); - } else if (meam_data.ialloy == 2) { + } else if (this->ialloy == 2) { dt1dr1 = 0.0; dt1dr2 = 0.0; @@ -352,17 +343,17 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, 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); + dt1dr1 = ai * (this->t1_meam[elti] - t1i); + dt1dr2 = aj * (this->t1_meam[elti] - t1j); + dt2dr1 = ai * (this->t2_meam[elti] - t2i); + dt2dr2 = aj * (this->t2_meam[elti] - t2j); + dt3dr1 = ai * (this->t3_meam[elti] - t3i); + dt3dr2 = aj * (this->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]); + get_shpfcn(shpi, this->lattce_meam[elti][elti]); + get_shpfcn(shpj, this->lattce_meam[eltj][eltj]); drhodr1 = arr1v(dGamma1, i) * drho0dr1 + arr1v(dGamma2, i) * (dt1dr1 * arr1v(rho1, i) + t1i * drho1dr1 + @@ -405,7 +396,7 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3; drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3; - if (meam_data.ialloy == 1) { + if (this->ialloy == 1) { a1i = 0.0; a1j = 0.0; @@ -426,20 +417,20 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, 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)); + dt1ds1 = a1i * (this->t1_meam[eltj] - + t1i * pow(this->t1_meam[eltj], 2)); + dt1ds2 = a1j * (this->t1_meam[elti] - + t1j * pow(this->t1_meam[elti], 2)); + dt2ds1 = a2i * (this->t2_meam[eltj] - + t2i * pow(this->t2_meam[eltj], 2)); + dt2ds2 = a2j * (this->t2_meam[elti] - + t2j * pow(this->t2_meam[elti], 2)); + dt3ds1 = a3i * (this->t3_meam[eltj] - + t3i * pow(this->t3_meam[eltj], 2)); + dt3ds2 = a3j * (this->t3_meam[elti] - + t3j * pow(this->t3_meam[elti], 2)); - } else if (meam_data.ialloy == 2) { + } else if (this->ialloy == 2) { dt1ds1 = 0.0; dt1ds2 = 0.0; @@ -457,12 +448,12 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, 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); + dt1ds1 = ai * (this->t1_meam[eltj] - t1i); + dt1ds2 = aj * (this->t1_meam[elti] - t1j); + dt2ds1 = ai * (this->t2_meam[eltj] - t2i); + dt2ds2 = aj * (this->t2_meam[elti] - t2j); + dt3ds1 = ai * (this->t3_meam[eltj] - t3i); + dt3ds2 = aj * (this->t3_meam[elti] - t3j); } drhods1 = @@ -602,4 +593,3 @@ meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, // else if elti=0, this is not a meam atom } } -} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp index 0b20827420..b05e3824b5 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -1,35 +1,6 @@ -extern "C" { #include "meam.h" #include -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 *) @@ -40,21 +11,21 @@ void G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, // void -meam_setup_done_(double* cutmax) +MEAM::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; + this->cutforce = this->rc_meam; + this->cutforcesq = this->cutforce * this->cutforce; // Pass cutoff back to calling program - *cutmax = meam_data.cutforce; + *cutmax = this->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]; + this->t1_meam[i] = + this->t1_meam[i] + this->augt1 * 3.0 / 5.0 * this->t3_meam[i]; // Compute off-diagonal alloy parameters alloyparams(); @@ -64,44 +35,44 @@ meam_setup_done_(double* cutmax) 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; + this->vind2D[m][n] = nv2; + this->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; + this->vind3D[m][n][p] = nv3; + this->vind3D[m][p][n] = nv3; + this->vind3D[n][m][p] = nv3; + this->vind3D[n][p][m] = nv3; + this->vind3D[p][m][n] = nv3; + this->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; + this->v2D[1] = 1; + this->v2D[2] = 2; + this->v2D[3] = 2; + this->v2D[4] = 1; + this->v2D[5] = 2; + this->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; + this->v3D[1] = 1; + this->v3D[2] = 3; + this->v3D[3] = 3; + this->v3D[4] = 3; + this->v3D[5] = 6; + this->v3D[6] = 3; + this->v3D[7] = 1; + this->v3D[8] = 3; + this->v3D[9] = 3; + this->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; + for (m = 1; m <= this->neltypes; m++) { + for (n = m; n <= this->neltypes; n++) { + this->eltind[m][n] = nv2; + this->eltind[n][m] = nv2; nv2 = nv2 + 1; } } @@ -110,71 +81,71 @@ meam_setup_done_(double* cutmax) 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; + this->nr = 1000; + this->dr = 1.1 * this->rc_meam / this->nr; compute_pair_meam(); } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc // Fill off-diagonal alloy parameters void -alloyparams(void) +MEAM::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++) { + for (i = 1; i <= this->neltypes; i++) { + for (j = 1; i <= this->neltypes; i++) { // Treat off-diagonal pairs // If i>j, set all equal to 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]; + this->re_meam[i][j] = this->re_meam[j][i]; + this->Ec_meam[i][j] = this->Ec_meam[j][i]; + this->alpha_meam[i][j] = this->alpha_meam[j][i]; + this->lattce_meam[i][j] = this->lattce_meam[j][i]; + this->nn2_meam[i][j] = this->nn2_meam[j][i]; // If i 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]; + if (iszero(this->Ec_meam[i][j])) { + if (this->lattce_meam[i][j] == L12) + this->Ec_meam[i][j] = + (3 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 4.0 - + this->delta_meam[i][j]; + else if (this->lattce_meam[i][j] == C11) { + if (this->lattce_meam[i][i] == DIA) + this->Ec_meam[i][j] = + (2 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 3.0 - + this->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]; + this->Ec_meam[i][j] = + (this->Ec_meam[i][i] + 2 * this->Ec_meam[j][j]) / 3.0 - + this->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]; + this->Ec_meam[i][j] = + (this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 2.0 - + this->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; + if (iszero(this->alpha_meam[i][j])) + this->alpha_meam[i][j] = + (this->alpha_meam[i][i] + this->alpha_meam[j][j]) / 2.0; + if (iszero(this->re_meam[i][j])) + this->re_meam[i][j] = + (this->re_meam[i][i] + this->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 ineltypes; 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]; + for (k = 1; k <= this->neltypes; k++) { + this->Cmin_meam[i][j][k] = this->Cmin_meam[j][i][k]; + this->Cmax_meam[i][j][k] = this->Cmax_meam[j][i][k]; } } } @@ -183,12 +154,12 @@ alloyparams(void) // atom k definitely lies outside the screening function ellipse (so // there is no need to calculate its effects). Here, compute it for all // triplets [i][j][k] so that ebound[i][j] is the maximized over k - for (i = 2; i <= 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); + for (i = 2; i <= this->neltypes; i++) { + for (j = 1; j <= this->neltypes; j++) { + for (k = 1; k <= this->neltypes; k++) { + eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) / + (4.0 * (this->Cmax_meam[i][j][k] - 1.0)); + this->ebound_meam[i][j] = max(this->ebound_meam[i][j], eb); } } } @@ -199,7 +170,7 @@ alloyparams(void) // void -compute_pair_meam(void) +MEAM::compute_pair_meam(void) { double r /*ununsed:, temp*/; @@ -211,62 +182,62 @@ compute_pair_meam(void) 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); + if (allocated(this->phir)) + meam_deallocate(this->phir); + if (allocated(this->phirar)) + meam_deallocate(this->phirar); + if (allocated(this->phirar1)) + meam_deallocate(this->phirar1); + if (allocated(this->phirar2)) + meam_deallocate(this->phirar2); + if (allocated(this->phirar3)) + meam_deallocate(this->phirar3); + if (allocated(this->phirar4)) + meam_deallocate(this->phirar4); + if (allocated(this->phirar5)) + meam_deallocate(this->phirar5); + if (allocated(this->phirar6)) + meam_deallocate(this->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_2d(this->phir, this->nr, + (this->neltypes * (this->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); + allocate_2d(this->phirar, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2); + allocate_2d(this->phirar1, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2); + allocate_2d(this->phirar2, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2); + allocate_2d(this->phirar3, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2); + allocate_2d(this->phirar4, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2); + allocate_2d(this->phirar5, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2); + allocate_2d(this->phirar6, this->nr, + (this->neltypes * (this->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++) { + for (a = 1; a <= this->neltypes; a++) { + for (b = a; b <= this->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; + for (j = 1; j <= this->nr; j++) { + r = (j - 1) * this->dr; - arr2(meam_data.phir, j, nv2) = phi_meam(r, a, b); + arr2(this->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]); + if (this->nn2_meam[a][b] == 1) { + get_Zij(&Z1, this->lattce_meam[a][b]); + get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][b], + this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b]); // The B1, B2, and L12 cases with NN2 have a trick to them; we // need to @@ -274,17 +245,17 @@ compute_pair_meam(void) // 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) { + if (this->lattce_meam[a][b] == B1 || + this->lattce_meam[a][b] == B2 || + this->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]); + get_Zij(&Z1, this->lattce_meam[a][a]); + get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][a], + this->Cmin_meam[a][a][a], + this->Cmax_meam[a][a][a]); nmax = 10; if (scrn > 0.0) { for (n = 1; n <= nmax; n++) { @@ -296,10 +267,10 @@ compute_pair_meam(void) // 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]); + get_Zij(&Z1, this->lattce_meam[b][b]); + get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[b][b], + this->Cmin_meam[b][b][b], + this->Cmax_meam[b][b][b]); nmax = 10; if (scrn > 0.0) { for (n = 1; n <= nmax; n++) { @@ -309,22 +280,22 @@ compute_pair_meam(void) } } - if (meam_data.lattce_meam[a][b] == B1 || - meam_data.lattce_meam[a][b] == B2) { + if (this->lattce_meam[a][b] == B1 || + this->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; + get_Zij(&Z1, this->lattce_meam[a][b]); + get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][b], + this->Cmin_meam[a][a][b], + this->Cmax_meam[a][a][b]); + arr2(this->phir, j, nv2) = + arr2(this->phir, j, nv2) - Z2 * scrn / (2 * Z1) * phiaa; + get_Zij2(&Z2, &arat, &scrn2, this->lattce_meam[a][b], + this->Cmin_meam[b][b][a], + this->Cmax_meam[b][b][a]); + arr2(this->phir, j, nv2) = + arr2(this->phir, j, nv2) - Z2 * scrn2 / (2 * Z1) * phibb; - } else if (meam_data.lattce_meam[a][b] == L12) { + } else if (this->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 @@ -339,7 +310,7 @@ compute_pair_meam(void) 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) - + arr2(this->phir, j, nv2) = arr2(this->phir, j, nv2) - 0.75 * S11 * phiaa - 0.25 * S22 * phibb; } @@ -347,8 +318,8 @@ compute_pair_meam(void) } else { nmax = 10; for (n = 1; n <= nmax; n++) { - arr2(meam_data.phir, j, nv2) = - arr2(meam_data.phir, j, nv2) + + arr2(this->phir, j, nv2) = + arr2(this->phir, j, nv2) + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b); } } @@ -360,17 +331,17 @@ compute_pair_meam(void) // else if -3 < astar < -1 // potential is linear combination with zbl potential // endif - if (meam_data.zbl_meam[a][b] == 1) { + if (this->zbl_meam[a][b] == 1) { astar = - meam_data.alpha_meam[a][b] * (r / meam_data.re_meam[a][b] - 1.0); + this->alpha_meam[a][b] * (r / this->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]); + arr2(this->phir, j, nv2) = + zbl(r, this->ielt_meam[a], this->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; + phizbl = zbl(r, this->ielt_meam[a], this->ielt_meam[b]); + arr2(this->phir, j, nv2) = + frac * arr2(this->phir, j, nv2) + (1 - frac) * phizbl; } } } @@ -385,7 +356,7 @@ compute_pair_meam(void) // Compute MEAM pair potential for distance r, element types a and b // double -phi_meam(double r, int a, int b) +MEAM::phi_meam(double r, int a, int b) { /*unused:double a1,a2,a12;*/ double t11av, t21av, t31av, t12av, t22av, t32av; @@ -409,7 +380,7 @@ phi_meam(double r, int a, int b) // get number of neighbors in the reference structure // Nref[i][j] = # of i's neighbors of type j - get_Zij(&Z12, meam_data.lattce_meam[a][b]); + get_Zij(&Z12, this->lattce_meam[a][b]); get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32); @@ -419,37 +390,37 @@ phi_meam(double r, int a, int b) 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]; + if (this->lattce_meam[a][b] == C11) { + if (this->ialloy == 2) { + t11av = this->t1_meam[a]; + t12av = this->t1_meam[b]; + t21av = this->t2_meam[a]; + t22av = this->t2_meam[b]; + t31av = this->t3_meam[a]; + t32av = this->t3_meam[b]; } else { scalfac = 1.0 / (rho01 + rho02); t11av = - scalfac * (meam_data.t1_meam[a] * rho01 + meam_data.t1_meam[b] * rho02); + scalfac * (this->t1_meam[a] * rho01 + this->t1_meam[b] * rho02); t12av = t11av; t21av = - scalfac * (meam_data.t2_meam[a] * rho01 + meam_data.t2_meam[b] * rho02); + scalfac * (this->t2_meam[a] * rho01 + this->t2_meam[b] * rho02); t22av = t21av; t31av = - scalfac * (meam_data.t3_meam[a] * rho01 + meam_data.t3_meam[b] * rho02); + scalfac * (this->t3_meam[a] * rho01 + this->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]); + this->t1_meam[a], this->t2_meam[a], this->t3_meam[a], + this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], + r, a, b, this->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 (this->lattce_meam[a][b] == C11) { + latta = this->lattce_meam[a][a]; if (latta == DIA) { rhobar1 = pow(((Z12 / 2) * (rho02 + rho01)), 2) + t11av * pow((rho12 - rho11), 2) + @@ -473,27 +444,27 @@ phi_meam(double r, int a, int b) // 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) + if (this->mix_ref_t == 1) { + Z1 = this->Z_meam[a]; + Z2 = this->Z_meam[b]; + if (this->ibar_meam[a] <= 0) G1 = 1.0; else { - get_shpfcn(s1, meam_data.lattce_meam[a][a]); + get_shpfcn(s1, this->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, + G_gam(Gam1, this->ibar_meam[a], this->gsmooth_factor, &G1, &errorflag); } - if (meam_data.ibar_meam[b] <= 0) + if (this->ibar_meam[b] <= 0) G2 = 1.0; else { - get_shpfcn(s2, meam_data.lattce_meam[b][b]); + get_shpfcn(s2, this->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, + G_gam(Gam2, this->ibar_meam[b], this->gsmooth_factor, &G2, &errorflag); } - rho0_1 = meam_data.rho0_meam[a] * Z1 * G1; - rho0_2 = meam_data.rho0_meam[b] * Z2 * G2; + rho0_1 = this->rho0_meam[a] * Z1 * G1; + rho0_2 = this->rho0_meam[b] * Z2 * G2; } Gam1 = (t11av * rho11 + t21av * rho21 + t31av * rho31); if (rho01 < 1.0e-14) @@ -507,20 +478,20 @@ phi_meam(double r, int a, int b) else Gam2 = Gam2 / (rho02 * rho02); - G_gam(Gam1, meam_data.ibar_meam[a], meam_data.gsmooth_factor, &G1, + G_gam(Gam1, this->ibar_meam[a], this->gsmooth_factor, &G1, &errorflag); - G_gam(Gam2, meam_data.ibar_meam[b], meam_data.gsmooth_factor, &G2, + G_gam(Gam2, this->ibar_meam[b], this->gsmooth_factor, &G2, &errorflag); - if (meam_data.mix_ref_t == 1) { + if (this->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]; + if (this->bkgd_dyn == 1) { + rho_bkgd1 = this->rho0_meam[a] * this->Z_meam[a]; + rho_bkgd2 = this->rho0_meam[b] * this->Z_meam[b]; } else { - rho_bkgd1 = meam_data.rho_ref_meam[a]; - rho_bkgd2 = meam_data.rho_ref_meam[b]; + rho_bkgd1 = this->rho_ref_meam[a]; + rho_bkgd2 = this->rho_ref_meam[b]; } } rhobar1 = rho01 / rho_bkgd1 * G1; @@ -531,30 +502,30 @@ phi_meam(double r, int a, int b) 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; + if (this->emb_lin_neg == 1 && rhobar1 <= 0) + F1 = -this->A_meam[a] * this->Ec_meam[a][a] * rhobar1; else F1 = - meam_data.A_meam[a] * meam_data.Ec_meam[a][a] * rhobar1 * log(rhobar1); + this->A_meam[a] * this->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; + if (this->emb_lin_neg == 1 && rhobar2 <= 0) + F2 = -this->A_meam[b] * this->Ec_meam[b][b] * rhobar2; else F2 = - meam_data.A_meam[b] * meam_data.Ec_meam[b][b] * rhobar2 * log(rhobar2); + this->A_meam[b] * this->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); + Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b], + this->Ec_meam[a][b], this->repuls_meam[a][b], + this->attrac_meam[a][b], this->erose_form); // calculate the pair energy - if (meam_data.lattce_meam[a][b] == C11) { - latta = meam_data.lattce_meam[a][a]; + if (this->lattce_meam[a][b] == C11) { + latta = this->lattce_meam[a][a]; if (latta == DIA) { phiaa = phi_meam(r, a, a); phi_m = (3 * Eu - F2 - 2 * F1 - 5 * phiaa) / Z12; @@ -562,12 +533,12 @@ phi_meam(double r, int a, int b) 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) { + } else if (this->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]); + get_Zij(&Z1nn, this->lattce_meam[a][a]); + get_Zij2(&Z2nn, &arat, &scrn, this->lattce_meam[a][a], + this->Cmin_meam[a][a][a], this->Cmax_meam[a][a][a]); nmax = 10; if (scrn > 0.0) { for (n = 1; n <= nmax; n++) { @@ -596,50 +567,50 @@ phi_meam(double r, int a, int b) //----------------------------------------------------------------------c // Compute background density for reference structure of each element void -compute_reference_density(void) +MEAM::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) + for (a = 1; a <= this->neltypes; a++) { + Z = (int)this->Z_meam[a]; + if (this->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]) / + get_shpfcn(shp, this->lattce_meam[a][a]); + gam = (this->t1_meam[a] * shp[1] + this->t2_meam[a] * shp[2] + + this->t3_meam[a] * shp[3]) / (Z * Z); - G_gam(gam, meam_data.ibar_meam[a], meam_data.gsmooth_factor, &Gbar, + G_gam(gam, this->ibar_meam[a], this->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; + rho0 = this->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]); + if (this->nn2_meam[a][a] == 1) { + get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][a], + this->Cmin_meam[a][a][a], this->Cmax_meam[a][a][a]); rho0_2nn = - meam_data.rho0_meam[a] * fm_exp(-meam_data.beta0_meam[a] * (arat - 1)); + this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * (arat - 1)); rho0 = rho0 + Z2 * rho0_2nn * scrn; } - meam_data.rho_ref_meam[a] = rho0 * Gbar; + this->rho_ref_meam[a] = rho0 * Gbar; } } //----------------------------------------------------------------------c // Shape factors for various configurations void -get_shpfcn(double* s /* s(3) */, lattice_t latt) +MEAM::get_shpfcn(double* s /* s(3) */, lattice_t latt) { if (latt == FCC || latt == BCC || latt == B1 || latt == B2) { s[1] = 0.0; @@ -667,7 +638,7 @@ get_shpfcn(double* s /* s(3) */, lattice_t latt) //------------------------------------------------------------------------------c // Average weighting factors for the reference structure void -get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, +MEAM::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) @@ -675,7 +646,7 @@ get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, double rhoa01, rhoa02, a1, a2, rho01 /*,rho02*/; // For ialloy = 2, no averaging is done - if (meam_data.ialloy == 2) { + if (this->ialloy == 2) { *t11av = t11; *t21av = t21; *t31av = t31; @@ -693,10 +664,10 @@ get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, *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); + a1 = r / this->re_meam[a][a] - 1.0; + a2 = r / this->re_meam[b][b] - 1.0; + rhoa01 = this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * a1); + rhoa02 = this->rho0_meam[b] * fm_exp(-this->beta0_meam[b] * a2); if (latt == L12) { rho01 = 8 * rhoa01 + 4 * rhoa02; *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01; @@ -715,7 +686,7 @@ get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, //------------------------------------------------------------------------------c // Number of neighbors for the reference structure void -get_Zij(int* Zij, lattice_t latt) +MEAM::get_Zij(int* Zij, lattice_t latt) { if (latt == FCC) *Zij = 12; @@ -745,7 +716,7 @@ get_Zij(int* Zij, lattice_t latt) // neighbor screening function for lattice type "latt" void -get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin, +MEAM::get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin, double cmax) { @@ -803,18 +774,18 @@ get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin, //------------------------------------------------------------------------------c void -get_sijk(double C, int i, int j, int k, double* sijk) +MEAM::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]); + x = (C - this->Cmin_meam[i][j][k]) / + (this->Cmax_meam[i][j][k] - this->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, +MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, double* rho31, double* rho02, double* rho12, double* rho22, double* rho32) { @@ -828,19 +799,19 @@ get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, 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; + a1 = r / this->re_meam[a][a] - 1.0; + a2 = r / this->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); + rhoa01 = this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * a1); + rhoa11 = this->rho0_meam[a] * fm_exp(-this->beta1_meam[a] * a1); + rhoa21 = this->rho0_meam[a] * fm_exp(-this->beta2_meam[a] * a1); + rhoa31 = this->rho0_meam[a] * fm_exp(-this->beta3_meam[a] * a1); + rhoa02 = this->rho0_meam[b] * fm_exp(-this->beta0_meam[b] * a2); + rhoa12 = this->rho0_meam[b] * fm_exp(-this->beta1_meam[b] * a2); + rhoa22 = this->rho0_meam[b] * fm_exp(-this->beta2_meam[b] * a2); + rhoa32 = this->rho0_meam[b] * fm_exp(-this->beta3_meam[b] * a2); - lat = meam_data.lattce_meam[a][b]; + lat = this->lattce_meam[a][b]; *rho11 = 0.0; *rho21 = 0.0; @@ -892,12 +863,12 @@ get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, } else if (lat == L12) { *rho01 = 8 * rhoa01 + 4 * rhoa02; *rho02 = 12 * rhoa01; - if (meam_data.ialloy == 1) { + if (this->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); + pow(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b], 2); + denom = 8 * rhoa01 * pow(this->t2_meam[a], 2) + + 4 * rhoa02 * pow(this->t2_meam[b], 2); if (denom > 0.) *rho21 = *rho21 / denom * *rho01; } else @@ -909,16 +880,16 @@ get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, // call error('Lattice not defined in get_densref.') } - if (meam_data.nn2_meam[a][b] == 1) { + if (this->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]); + get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[a][a][b], + this->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; + a1 = arat * r / this->re_meam[a][a] - 1.0; + a2 = arat * r / this->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); + rhoa01nn = this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * a1); + rhoa02nn = this->rho0_meam[b] * fm_exp(-this->beta0_meam[b] * a2); if (lat == L12) { // As usual, L12 thinks it's special; we need to be careful computing @@ -939,8 +910,8 @@ get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, *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]); + get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[b][b][a], + this->Cmax_meam[b][b][a]); *rho02 = *rho02 + Zij2nn * scrn * rhoa02nn; } } @@ -950,7 +921,7 @@ get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, // Compute ZBL potential // double -zbl(double r, int z1, int z2) +MEAM::zbl(double r, int z1, int z2) { int i; const double c[] = { 0.028171, 0.28022, 0.50986, 0.18175 }; @@ -974,7 +945,7 @@ zbl(double r, int z1, int z2) // Compute Rose energy function, I.16 // double -erose(double r, double re, double alpha, double Ec, double repuls, +MEAM::erose(double r, double re, double alpha, double Ec, double repuls, double attrac, int form) { double astar, a3; @@ -1003,58 +974,58 @@ erose(double r, double re, double alpha, double Ec, double repuls, // ----------------------------------------------------------------------- void -interpolate_meam(int ind) +MEAM::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; + this->nrar = this->nr; + drar = this->dr; + this->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); + for (j = 1; j <= this->nrar; j++) { + arr2(this->phirar, j, ind) = arr2(this->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))) / + arr2(this->phirar1, 1, ind) = + arr2(this->phirar, 2, ind) - arr2(this->phirar, 1, ind); + arr2(this->phirar1, 2, ind) = + 0.5 * (arr2(this->phirar, 3, ind) - arr2(this->phirar, 1, ind)); + arr2(this->phirar1, this->nrar - 1, ind) = + 0.5 * (arr2(this->phirar, this->nrar, ind) - + arr2(this->phirar, this->nrar - 2, ind)); + arr2(this->phirar1, this->nrar, ind) = 0.0; + for (j = 3; j <= this->nrar - 2; j++) { + arr2(this->phirar1, j, ind) = + ((arr2(this->phirar, j - 2, ind) - + arr2(this->phirar, j + 2, ind)) + + 8.0 * (arr2(this->phirar, j + 1, ind) - + arr2(this->phirar, j - 1, ind))) / 12.; } - for (j = 1; j <= meam_data.nrar - 1; j++) { - arr2(meam_data.phirar2, j, ind) = + for (j = 1; j <= this->nrar - 1; j++) { + arr2(this->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) - + (arr2(this->phirar, j + 1, ind) - arr2(this->phirar, j, ind)) - + 2.0 * arr2(this->phirar1, j, ind) - + arr2(this->phirar1, j + 1, ind); + arr2(this->phirar3, j, ind) = + arr2(this->phirar1, j, ind) + arr2(this->phirar1, j + 1, ind) - 2.0 * - (arr2(meam_data.phirar, j + 1, ind) - arr2(meam_data.phirar, j, ind)); + (arr2(this->phirar, j + 1, ind) - arr2(this->phirar, j, ind)); } - arr2(meam_data.phirar2, meam_data.nrar, ind) = 0.0; - arr2(meam_data.phirar3, meam_data.nrar, ind) = 0.0; + arr2(this->phirar2, this->nrar, ind) = 0.0; + arr2(this->phirar3, this->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; + for (j = 1; j <= this->nrar; j++) { + arr2(this->phirar4, j, ind) = arr2(this->phirar1, j, ind) / drar; + arr2(this->phirar5, j, ind) = + 2.0 * arr2(this->phirar2, j, ind) / drar; + arr2(this->phirar6, j, ind) = + 3.0 * arr2(this->phirar3, j, ind) / drar; } } @@ -1062,24 +1033,23 @@ interpolate_meam(int ind) // Compute Rose energy function, I.16 // double -compute_phi(double rij, int elti, int eltj) +MEAM::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; + ind = this->eltind[elti][eltj]; + pp = rij * this->rdrar + 1.0; kk = (int)pp; - kk = min(kk, meam_data.nrar - 1); + kk = min(kk, this->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)) * + double result = ((arr2(this->phirar3, kk, ind) * pp + + arr2(this->phirar2, kk, ind)) * pp + - arr2(meam_data.phirar1, kk, ind)) * + arr2(this->phirar1, kk, ind)) * pp + - arr2(meam_data.phirar, kk, ind); + arr2(this->phirar, kk, ind); return result; } -} diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index eac03277f2..dba4f40d18 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -1,4 +1,3 @@ -extern "C" { #include "meam.h" #include @@ -18,7 +17,7 @@ extern "C" { // void -meam_setup_global_(int* nelt, int* lat, double* z, int* ielement, double* atwt, +MEAM::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, @@ -28,72 +27,71 @@ meam_setup_global_(int* nelt, int* lat, double* z, int* ielement, double* atwt, int i; double tmplat[maxelt]; - meam_data.neltypes = *nelt; + this->neltypes = *nelt; for (i = 1; i <= *nelt; i++) { if (arr1v(lat, i) == 0) - meam_data.lattce_meam[i][i] = FCC; + this->lattce_meam[i][i] = FCC; else if (arr1v(lat, i) == 1) - meam_data.lattce_meam[i][i] = BCC; + this->lattce_meam[i][i] = BCC; else if (arr1v(lat, i) == 2) - meam_data.lattce_meam[i][i] = HCP; + this->lattce_meam[i][i] = HCP; else if (arr1v(lat, i) == 3) - meam_data.lattce_meam[i][i] = DIM; + this->lattce_meam[i][i] = DIM; else if (arr1v(lat, i) == 4) - meam_data.lattce_meam[i][i] = DIA; + this->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); + this->Z_meam[i] = arr1v(z, i); + this->ielt_meam[i] = arr1v(ielement, i); + this->alpha_meam[i][i] = arr1v(alpha, i); + this->beta0_meam[i] = arr1v(b0, i); + this->beta1_meam[i] = arr1v(b1, i); + this->beta2_meam[i] = arr1v(b2, i); + this->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); + this->Ec_meam[i][i] = arr1v(esub, i); + this->A_meam[i] = arr1v(asub, i); + this->t0_meam[i] = arr1v(t0, i); + this->t1_meam[i] = arr1v(t1, i); + this->t2_meam[i] = arr1v(t2, i); + this->t3_meam[i] = arr1v(t3, i); + this->rho0_meam[i] = arr1v(rozero, i); + this->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; + if (this->lattce_meam[i][i] == FCC) + this->re_meam[i][i] = tmplat[i] / sqrt(2.0); + else if (this->lattce_meam[i][i] == BCC) + this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0; + else if (this->lattce_meam[i][i] == HCP) + this->re_meam[i][i] = tmplat[i]; + else if (this->lattce_meam[i][i] == DIM) + this->re_meam[i][i] = tmplat[i]; + else if (this->lattce_meam[i][i] == DIA) + this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0; else { // 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; + this->rc_meam = 4.0; + this->delr_meam = 0.1; + setall2d(this->attrac_meam, 0.0); + setall2d(this->repuls_meam, 0.0); + setall3d(this->Cmax_meam, 2.8); + setall3d(this->Cmin_meam, 2.0); + setall2d(this->ebound_meam, pow(2.8, 2) / (4.0 * (2.8 - 1.0))); + setall2d(this->delta_meam, 0.0); + setall2d(this->nn2_meam, 0); + setall2d(this->zbl_meam, 1); + this->gsmooth_factor = 99.0; + this->augt1 = 1; + this->ialloy = 0; + this->mix_ref_t = 0; + this->emb_lin_neg = 0; + this->bkgd_dyn = 0; + this->erose_form = 0; } -} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index 4be38a0cb1..7b8b5e2cef 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -1,10 +1,9 @@ -extern "C" { #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) +MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr) { //: idx[0..2] *ierr = 0; @@ -58,7 +57,7 @@ meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr) // 20 = bkgd_dyn void -meam_setup_param_(int* which_p, double* value_p, int* nindex_p, +MEAM::meam_setup_param_(int* which_p, double* value_p, int* nindex_p, int* index /*index(3)*/, int* errorflag) { //: index[0..2] @@ -74,7 +73,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.Ec_meam[index[0]][index[1]] = value; + this->Ec_meam[index[0]][index[1]] = value; break; // 1 = alpha_meam @@ -82,7 +81,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.alpha_meam[index[0]][index[1]] = value; + this->alpha_meam[index[0]][index[1]] = value; break; // 2 = rho0_meam @@ -90,7 +89,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(1, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.rho0_meam[index[0]] = value; + this->rho0_meam[index[0]] = value; break; // 3 = delta_meam @@ -98,7 +97,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.delta_meam[index[0]][index[1]] = value; + this->delta_meam[index[0]][index[1]] = value; break; // 4 = lattce_meam @@ -109,23 +108,23 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, val = (int)value; if (val == 0) - meam_data.lattce_meam[index[0]][index[1]] = FCC; + this->lattce_meam[index[0]][index[1]] = FCC; else if (val == 1) - meam_data.lattce_meam[index[0]][index[1]] = BCC; + this->lattce_meam[index[0]][index[1]] = BCC; else if (val == 2) - meam_data.lattce_meam[index[0]][index[1]] = HCP; + this->lattce_meam[index[0]][index[1]] = HCP; else if (val == 3) - meam_data.lattce_meam[index[0]][index[1]] = DIM; + this->lattce_meam[index[0]][index[1]] = DIM; else if (val == 4) - meam_data.lattce_meam[index[0]][index[1]] = DIA; + this->lattce_meam[index[0]][index[1]] = DIA; else if (val == 5) - meam_data.lattce_meam[index[0]][index[1]] = B1; + this->lattce_meam[index[0]][index[1]] = B1; else if (val == 6) - meam_data.lattce_meam[index[0]][index[1]] = C11; + this->lattce_meam[index[0]][index[1]] = C11; else if (val == 7) - meam_data.lattce_meam[index[0]][index[1]] = L12; + this->lattce_meam[index[0]][index[1]] = L12; else if (val == 8) - meam_data.lattce_meam[index[0]][index[1]] = B2; + this->lattce_meam[index[0]][index[1]] = B2; break; // 5 = attrac_meam @@ -133,7 +132,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.attrac_meam[index[0]][index[1]] = value; + this->attrac_meam[index[0]][index[1]] = value; break; // 6 = repuls_meam @@ -141,7 +140,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.repuls_meam[index[0]][index[1]] = value; + this->repuls_meam[index[0]][index[1]] = value; break; // 7 = nn2_meam @@ -151,7 +150,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, return; i1 = min(index[0], index[1]); i2 = max(index[0], index[1]); - meam_data.nn2_meam[i1][i2] = (int)value; + this->nn2_meam[i1][i2] = (int)value; break; // 8 = Cmin_meam @@ -159,7 +158,7 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(3, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.Cmin_meam[index[0]][index[1]][index[2]] = value; + this->Cmin_meam[index[0]][index[1]][index[2]] = value; break; // 9 = Cmax_meam @@ -167,27 +166,27 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(3, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.Cmax_meam[index[0]][index[1]][index[2]] = value; + this->Cmax_meam[index[0]][index[1]][index[2]] = value; break; // 10 = rc_meam case 10: - meam_data.rc_meam = value; + this->rc_meam = value; break; // 11 = delr_meam case 11: - meam_data.delr_meam = value; + this->delr_meam = value; break; // 12 = augt1 case 12: - meam_data.augt1 = (int)value; + this->augt1 = (int)value; break; // 13 = gsmooth case 13: - meam_data.gsmooth_factor = value; + this->gsmooth_factor = value; break; // 14 = re_meam @@ -195,22 +194,22 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - meam_data.re_meam[index[0]][index[1]] = value; + this->re_meam[index[0]][index[1]] = value; break; // 15 = ialloy case 15: - meam_data.ialloy = (int)value; + this->ialloy = (int)value; break; // 16 = mixture_ref_t case 16: - meam_data.mix_ref_t = (int)value; + this->mix_ref_t = (int)value; break; // 17 = erose_form case 17: - meam_data.erose_form = (int)value; + this->erose_form = (int)value; break; // 18 = zbl_meam @@ -220,21 +219,20 @@ meam_setup_param_(int* which_p, double* value_p, int* nindex_p, return; i1 = min(index[0], index[1]); i2 = max(index[0], index[1]); - meam_data.zbl_meam[i1][i2] = (int)value; + this->zbl_meam[i1][i2] = (int)value; break; // 19 = emb_lin_neg case 19: - meam_data.emb_lin_neg = (int)value; + this->emb_lin_neg = (int)value; break; // 20 = bkgd_dyn case 20: - meam_data.bkgd_dyn = (int)value; + this->bkgd_dyn = (int)value; break; default: *errorflag = 1; } } -} \ No newline at end of file diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index f2c7da9239..a92c9c6b50 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -34,7 +34,6 @@ 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", @@ -78,7 +77,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp) PairMEAMC::~PairMEAMC() { - meam_cleanup_(); + meam_inst.meam_cleanup_(); memory->destroy(rho); memory->destroy(rho0); @@ -237,7 +236,7 @@ void PairMEAMC::compute(int eflag, int vflag) 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], + meam_inst.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], @@ -254,7 +253,7 @@ void PairMEAMC::compute(int eflag, int vflag) comm->reverse_comm_pair(this); - meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, + meam_inst.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, @@ -279,7 +278,7 @@ void PairMEAMC::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom, + meam_inst.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], @@ -368,7 +367,7 @@ void PairMEAMC::coeff(int narg, char **arg) // tell MEAM package that setup is done read_files(arg[2],arg[2+nelements+1]); - meam_setup_done_(&cutmax); + meam_inst.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 @@ -567,7 +566,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) 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 if (strcmp(words[1],"dia") == 0) lat[i] = DIA; else error->all(FLERR,"Unrecognized lattice type in MEAM file 1"); // store parameters @@ -600,7 +599,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) // pass element parameters to MEAM package - meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, + meam_inst.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 @@ -704,7 +703,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) 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],"dia") == 0) value = DIA; else if (strcmp(params[nparams-1],"b1") == 0) value = B1; else if (strcmp(params[nparams-1],"c11") == 0) value = C11; else if (strcmp(params[nparams-1],"l12") == 0) value = L12; @@ -716,7 +715,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) // pass single setting to MEAM package int errorflag = 0; - meam_setup_param_(&which,&value,&nindex,index,&errorflag); + meam_inst.meam_setup_param_(&which,&value,&nindex,index,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h index c4e03aef55..c9aee39364 100644 --- a/src/USER-MEAMC/pair_meamc.h +++ b/src/USER-MEAMC/pair_meamc.h @@ -20,41 +20,8 @@ PairStyle(meam/c,PairMEAMC) #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" +#include "meam.h" namespace LAMMPS_NS { @@ -76,6 +43,7 @@ class PairMEAMC : public Pair { double memory_usage(); private: + MEAM meam_inst; double cutmax; // max cutoff for all elements int nelements; // # of unique elements char **elements; // names of unique elements