diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 3ef6bf8a4..8cf51c773 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -1,219 +1,212 @@ #ifndef LMP_MEAM_H #define LMP_MEAM_H #include #define maxelt 5 extern "C" { double fm_exp(double); } +namespace LAMMPS_NS { + typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t; typedef struct { int dim1, dim2; double* ptr; } allocatable_double_2d; class MEAM { -private: + public: + MEAM() {}; + ~MEAM() { + meam_cleanup_(); + } + private: // cutforce = force cutoff // cutforcesq = force cutoff squared double cutforce, cutforcesq; // Ec_meam = cohesive energy // re_meam = nearest-neighbor distance // Omega_meam = atomic volume // B_meam = bulk modulus // Z_meam = number of first neighbors for reference structure // ielt_meam = atomic number of element // A_meam = adjustable parameter // alpha_meam = sqrt(9*Omega*B/Ec) // rho0_meam = density scaling parameter // delta_meam = heat of formation for alloys // beta[0-3]_meam = electron density constants // t[0-3]_meam = coefficients on densities in Gamma computation // rho_ref_meam = background density for reference structure // ibar_meam(i) = selection parameter for Gamma function for elt i, // lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j) // neltypes = maximum number of element type defined // eltind = index number of pair (similar to Voigt notation; ij = ji) // phir = pair potential function array // phirar[1-6] = spline coeffs // attrac_meam = attraction parameter in Rose energy // repuls_meam = repulsion parameter in Rose energy // nn2_meam = 1 if second nearest neighbors are to be computed, else 0 // zbl_meam = 1 if zbl potential for small r to be use, else 0 // emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0 // bkgd_dyn = 1 if reference densities follows Dynamo, else 0 // Cmin_meam, Cmax_meam = min and max values in screening cutoff // rc_meam = cutoff distance for meam // delr_meam = cutoff region for meam // ebound_meam = factor giving maximum boundary of sceen fcn ellipse // augt1 = flag for whether t1 coefficient should be augmented // ialloy = flag for newer alloy formulation (as in dynamo code) // mix_ref_t = flag to recover "old" way of computing t in reference config // erose_form = selection parameter for form of E_rose function // gsmooth_factor = factor determining length of G smoothing region // vind[23]D = Voight notation index maps for 2 and 3D // v2D,v3D = array of factors to apply for Voight notation // nr,dr = pair function discretization parameters // nrar,rdrar = spline coeff array parameters double Ec_meam[maxelt + 1][maxelt + 1], re_meam[maxelt + 1][maxelt + 1]; double Omega_meam[maxelt + 1], Z_meam[maxelt + 1]; double A_meam[maxelt + 1], alpha_meam[maxelt + 1][maxelt + 1], rho0_meam[maxelt + 1]; double delta_meam[maxelt + 1][maxelt + 1]; double beta0_meam[maxelt + 1], beta1_meam[maxelt + 1]; double beta2_meam[maxelt + 1], beta3_meam[maxelt + 1]; double t0_meam[maxelt + 1], t1_meam[maxelt + 1]; double t2_meam[maxelt + 1], t3_meam[maxelt + 1]; double rho_ref_meam[maxelt + 1]; int ibar_meam[maxelt + 1], ielt_meam[maxelt + 1]; lattice_t lattce_meam[maxelt + 1][maxelt + 1]; int nn2_meam[maxelt + 1][maxelt + 1]; int zbl_meam[maxelt + 1][maxelt + 1]; int eltind[maxelt + 1][maxelt + 1]; int neltypes; allocatable_double_2d phir; // [:,:] allocatable_double_2d phirar, phirar1, phirar2, phirar3, phirar4, phirar5, phirar6; // [:,:] double attrac_meam[maxelt + 1][maxelt + 1], repuls_meam[maxelt + 1][maxelt + 1]; double Cmin_meam[maxelt + 1][maxelt + 1][maxelt + 1]; double Cmax_meam[maxelt + 1][maxelt + 1][maxelt + 1]; double rc_meam, delr_meam, ebound_meam[maxelt + 1][maxelt + 1]; int augt1, ialloy, mix_ref_t, erose_form; int emb_lin_neg, bkgd_dyn; double gsmooth_factor; int vind2D[3 + 1][3 + 1], vind3D[3 + 1][3 + 1][3 + 1]; int v2D[6 + 1], v3D[10 + 1]; int nr, nrar; double dr, rdrar; -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*); + 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; \ - }) -#endif - -#ifndef min -#define min(a, b) \ - ({ \ - __typeof__(a) _a = (a); \ - __typeof__(b) _b = (b); \ - _a < _b ? _a : _b; \ - }) -#endif +#define MIN(A,B) ((A) < (B) ? (A) : (B)) +#define MAX(A,B) ((A) > (B) ? (A) : (B)) #define iszero(f) (fabs(f) < 1e-20) #define setall2d(arr, v) \ { \ for (int __i = 1; __i <= maxelt; __i++) \ for (int __j = 1; __j <= maxelt; __j++) \ arr[__i][__j] = v; \ } #define setall3d(arr, v) \ { \ for (int __i = 1; __i <= maxelt; __i++) \ for (int __j = 1; __j <= maxelt; __j++) \ for (int __k = 1; __k <= maxelt; __k++) \ arr[__i][__j][__k] = v; \ } /* Fortran Array Semantics in C. - Stack-Allocated and global arrays are 1-based, declared as foo[N+1] and simply ignoring the first element - Multi-Dimensional MUST be declared in reverse, so that the order when accessing is the same as in Fortran - arrays that are passed externally via the meam_* functions must use the arr*v() functions below (or be used with 0-based indexing) - allocatable arrays (only global phir*) are actually a struct, and must be accessed with arr2() */ // we receive a pointer to the first element, and F dimensions is ptr(a,b,c) // we know c data structure is ptr[c][b][a] #define arrdim2v(ptr, a, b) \ const int DIM1__##ptr = a; \ const int DIM2__##ptr = b; \ (void)(DIM1__##ptr); \ (void)(DIM2__##ptr); #define arrdim3v(ptr, a, b, c) \ const int DIM1__##ptr = a; \ const int DIM2__##ptr = b; \ const int DIM3__##ptr = c; \ (void)(DIM1__##ptr); (void)(DIM2__##ptr; (void)(DIM3__##ptr); // access data with same index as used in fortran (1-based) #define arr1v(ptr, i) ptr[i - 1] #define arr2v(ptr, i, j) ptr[(DIM1__##ptr) * (j - 1) + (i - 1)] #define arr3v(ptr, i, j, k) \ ptr[(i - 1) + (j - 1) * (DIM1__##ptr) + \ (k - 1) * (DIM1__##ptr) * (DIM2__##ptr)] // allocatable arrays via special type #define allocate_2d(arr, cols, rows) \ arr.dim1 = cols; \ arr.dim2 = rows; \ arr.ptr = (double*)calloc((size_t)(cols) * (size_t)(rows), sizeof(double)); #define allocated(a) (a.ptr != NULL) #define 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)] +}; -#endif \ No newline at end of file +#endif diff --git a/src/USER-MEAMC/meam_cleanup.cpp b/src/USER-MEAMC/meam_cleanup.cpp index b0ee013c4..1b1dc4595 100644 --- a/src/USER-MEAMC/meam_cleanup.cpp +++ b/src/USER-MEAMC/meam_cleanup.cpp @@ -1,15 +1,17 @@ #include "meam.h" +using namespace LAMMPS_NS; + void MEAM::meam_cleanup_(void) { 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_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index 8a63fbd82..778787a06 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -1,282 +1,283 @@ #include "meam.h" #include +using namespace LAMMPS_NS; // Extern "C" declaration has the form: // // void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, // int *, int *, int *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, int *); // // Call from pair_meam.cpp has the form: // // meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, // &eng_vdwl,eatom,ntype,type,fmap, // &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], // &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1, // dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); // void MEAM::meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, int* eflag_atom, double* eng_vdwl, double* eatom, int* ntype, int* type, int* fmap, double* Arho1, double* Arho2, double* Arho2b, double* Arho3, double* Arho3b, double* t_ave, double* tsq_ave, double* Gamma, double* dGamma1, double* dGamma2, double* dGamma3, double* rho, double* rho0, double* rho1, double* rho2, double* rho3, double* fp, int* errorflag) { arrdim2v(Arho1, 3, *nmax); arrdim2v(Arho2, 6, *nmax); arrdim2v(Arho3, 10, *nmax); arrdim2v(Arho3b, 3, *nmax); arrdim2v(t_ave, 3, *nmax); arrdim2v(tsq_ave, 3, *nmax); int i, elti; int m; double rhob, G, dG, Gbar, dGbar, gam, shp[3 + 1], Z; double B, denom, rho_bkgd; // Complete the calculation of density for (i = 1; i <= *nlocal; i++) { elti = arr1v(fmap, arr1v(type, i)); if (elti > 0) { arr1v(rho1, i) = 0.0; arr1v(rho2, i) = -1.0 / 3.0 * arr1v(Arho2b, i) * arr1v(Arho2b, i); arr1v(rho3, i) = 0.0; for (m = 1; m <= 3; m++) { arr1v(rho1, i) = arr1v(rho1, i) + arr2v(Arho1, m, i) * arr2v(Arho1, m, i); arr1v(rho3, i) = arr1v(rho3, i) - 3.0 / 5.0 * arr2v(Arho3b, m, i) * arr2v(Arho3b, m, i); } for (m = 1; m <= 6; m++) { arr1v(rho2, i) = arr1v(rho2, i) + this->v2D[m] * arr2v(Arho2, m, i) * arr2v(Arho2, m, i); } for (m = 1; m <= 10; m++) { arr1v(rho3, i) = arr1v(rho3, i) + this->v3D[m] * arr2v(Arho3, m, i) * arr2v(Arho3, m, i); } if (arr1v(rho0, i) > 0.0) { 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 (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); arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) / arr1v(rho0, i); } } arr1v(Gamma, i) = arr2v(t_ave, 1, i) * arr1v(rho1, i) + arr2v(t_ave, 2, i) * arr1v(rho2, i) + arr2v(t_ave, 3, i) * arr1v(rho3, i); if (arr1v(rho0, i) > 0.0) { arr1v(Gamma, i) = arr1v(Gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i)); } Z = this->Z_meam[elti]; G_gam(arr1v(Gamma, i), this->ibar_meam[elti], this->gsmooth_factor, &G, errorflag); if (*errorflag != 0) return; get_shpfcn(shp, this->lattce_meam[elti][elti]); if (this->ibar_meam[elti] <= 0) { Gbar = 1.0; dGbar = 0.0; } else { 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 = (this->t1_meam[elti] * shp[1] + this->t2_meam[elti] * shp[2] + this->t3_meam[elti] * shp[3]) / (Z * Z); } G_gam(gam, this->ibar_meam[elti], this->gsmooth_factor, &Gbar, errorflag); } arr1v(rho, i) = arr1v(rho0, i) * G; 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, this->ibar_meam[elti], this->gsmooth_factor, &Gbar, &dGbar); } rho_bkgd = this->rho0_meam[elti] * Z * Gbar; } else { if (this->bkgd_dyn == 1) { rho_bkgd = this->rho0_meam[elti] * Z; } else { rho_bkgd = this->rho_ref_meam[elti]; } } rhob = arr1v(rho, i) / rho_bkgd; denom = 1.0 / rho_bkgd; dG_gam(arr1v(Gamma, i), this->ibar_meam[elti], this->gsmooth_factor, &G, &dG); arr1v(dGamma1, i) = (G - 2 * dG * arr1v(Gamma, i)) * denom; if (!iszero(arr1v(rho0, i))) { arr1v(dGamma2, i) = (dG / arr1v(rho0, i)) * denom; } else { arr1v(dGamma2, i) = 0.0; } // dGamma3 is nonzero only if we are using the "mixed" rule for // computing t in the reference system (which is not correct, but // included for backward compatibility if (this->mix_ref_t == 1) { arr1v(dGamma3, i) = arr1v(rho0, i) * G * dGbar / (Gbar * Z * Z) * denom; } else { arr1v(dGamma3, i) = 0.0; } B = this->A_meam[elti] * this->Ec_meam[elti][elti]; if (!iszero(rhob)) { 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 (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 (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); } } } } else { if (this->emb_lin_neg == 1) { arr1v(fp, i) = -B; } else { arr1v(fp, i) = B; } } } } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void 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) // 1 => G = exp(Gamma/2) // 2 => not implemented // 3 => G = 2/(1+exp(-Gamma)) // 4 => G = sqrt(1+Gamma) // -5 => G = +-sqrt(abs(1+Gamma)) double gsmooth_switchpoint; if (ibar == 0 || ibar == 4) { gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1); if (Gamma < gsmooth_switchpoint) { // e.g. gsmooth_factor is 99, {: // gsmooth_switchpoint = -0.99 // G = 0.01*(-0.99/Gamma)**99 *G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / Gamma), gsmooth_factor); *G = sqrt(*G); } else { *G = sqrt(1.0 + Gamma); } } else if (ibar == 1) { *G = fm_exp(Gamma / 2.0); } else if (ibar == 3) { *G = 2.0 / (1.0 + exp(-Gamma)); } else if (ibar == -5) { if ((1.0 + Gamma) >= 0) { *G = sqrt(1.0 + Gamma); } else { *G = -sqrt(-1.0 - Gamma); } } else { *errorflag = 1; } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void 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) // 1 => G = fm_exp(Gamma/2) // 2 => not implemented // 3 => G = 2/(1+fm_exp(-Gamma)) // 4 => G = sqrt(1+Gamma) // -5 => G = +-sqrt(abs(1+Gamma)) double gsmooth_switchpoint; if (ibar == 0 || ibar == 4) { gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1); if (Gamma < gsmooth_switchpoint) { // e.g. gsmooth_factor is 99, {: // gsmooth_switchpoint = -0.99 // G = 0.01*(-0.99/Gamma)**99 *G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / Gamma), gsmooth_factor); *G = sqrt(*G); *dG = -gsmooth_factor * *G / (2.0 * Gamma); } else { *G = sqrt(1.0 + Gamma); *dG = 1.0 / (2.0 * *G); } } else if (ibar == 1) { *G = fm_exp(Gamma / 2.0); *dG = *G / 2.0; } else if (ibar == 3) { *G = 2.0 / (1.0 + fm_exp(-Gamma)); *dG = *G * (2.0 - *G) / 2; } else if (ibar == -5) { if ((1.0 + Gamma) >= 0) { *G = sqrt(1.0 + Gamma); *dG = 1.0 / (2.0 * *G); } else { *G = -sqrt(-1.0 - Gamma); *dG = -1.0 / (2.0 * *G); } } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index b94eec65c..f03bc5a12 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -1,522 +1,523 @@ #include "meam.h" #include +using namespace LAMMPS_NS; // Extern "C" declaration has the form: // // void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, // double *, // int *, int *, int *, int *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, int *); // // // Call from pair_meam.cpp has the form: // // meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0], // &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], // &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], // rho0,&arho1[0][0],&arho2[0][0],arho2b, // &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag); // void MEAM::meam_dens_init_(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x, int* numneigh, int* firstneigh, int* numneigh_full, int* firstneigh_full, double* scrfcn, double* dscrfcn, double* fcpair, double* rho0, double* arho1, double* arho2, double* arho2b, double* arho3, double* arho3b, double* t_ave, double* tsq_ave, int* errorflag) { *errorflag = 0; // Compute screening function and derivatives getscreen(*i, *nmax, scrfcn, dscrfcn, fcpair, x, *numneigh, firstneigh, *numneigh_full, firstneigh_full, *ntype, type, fmap); // Calculate intermediate density terms to be communicated calc_rho1(*i, *nmax, *ntype, type, fmap, x, *numneigh, firstneigh, scrfcn, fcpair, rho0, arho1, arho2, arho2b, arho3, arho3b, t_ave, tsq_ave); } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void 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) { arrdim2v(x, 3, nmax); int jn, j, kn, k; int elti, eltj, eltk; double xitmp, yitmp, zitmp, delxij, delyij, delzij, rij2, rij; double xjtmp, yjtmp, zjtmp, delxik, delyik, delzik, rik2 /*,rik*/; double xktmp, yktmp, zktmp, delxjk, delyjk, delzjk, rjk2 /*,rjk*/; double xik, xjk, sij, fcij, sfcij, dfcij, sikj, dfikj, cikj; double Cmin, Cmax, delc, /*ebound,*/ rbound, a, coef1, coef2; // double coef1a,coef1b,coef2a,coef2b; double dCikj; /*unused:double dC1a,dC1b,dC2a,dC2b;*/ double rnorm, fc, dfc, drinv; drinv = 1.0 / this->delr_meam; elti = arr1v(fmap, arr1v(type, i)); if (elti > 0) { xitmp = arr2v(x, 1, i); yitmp = arr2v(x, 2, i); zitmp = arr2v(x, 3, i); for (jn = 1; jn <= numneigh; jn++) { j = arr1v(firstneigh, jn); eltj = arr1v(fmap, arr1v(type, j)); if (eltj > 0) { // First compute screening function itself, sij xjtmp = arr2v(x, 1, j); yjtmp = arr2v(x, 2, j); zjtmp = arr2v(x, 3, j); delxij = xjtmp - xitmp; delyij = yjtmp - yitmp; delzij = zjtmp - zitmp; rij2 = delxij * delxij + delyij * delyij + delzij * delzij; rij = sqrt(rij2); if (rij > this->rc_meam) { fcij = 0.0; dfcij = 0.0; sij = 0.0; } else { rnorm = (this->rc_meam - rij) * drinv; screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full, ntype, type, fmap); dfcut(rnorm, &fc, &dfc); fcij = fc; dfcij = dfc * drinv; } // Now compute derivatives arr1v(dscrfcn, jn) = 0.0; sfcij = sij * fcij; if (iszero(sfcij) || iszero(sfcij - 1.0)) goto LABEL_100; rbound = this->ebound_meam[elti][eltj] * rij2; for (kn = 1; kn <= numneigh_full; kn++) { k = arr1v(firstneigh_full, kn); if (k == j) continue; eltk = arr1v(fmap, arr1v(type, k)); if (eltk == 0) continue; xktmp = arr2v(x, 1, k); yktmp = arr2v(x, 2, k); zktmp = arr2v(x, 3, k); delxjk = xktmp - xjtmp; delyjk = yktmp - yjtmp; delzjk = zktmp - zjtmp; rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk; if (rjk2 > rbound) continue; delxik = xktmp - xitmp; delyik = yktmp - yitmp; delzik = zktmp - zitmp; rik2 = delxik * delxik + delyik * delyik + delzik * delzik; if (rik2 > rbound) continue; xik = rik2 / rij2; xjk = rjk2 / rij2; a = 1 - (xik - xjk) * (xik - xjk); // if a < 0, then ellipse equation doesn't describe this case and // atom k can't possibly screen i-j if (a <= 0.0) continue; cikj = (2.0 * (xik + xjk) + a - 2.0) / a; Cmax = 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 // tolerance) if atoms are colinear, so don't reject that case // here // (other negative cikj cases were handled by the test on "a" // above) // Note that we never have 0cutforcesq) { eltj = arr1v(fmap, arr1v(type, j)); rij = sqrt(rij2); 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 (this->ialloy != 2) { arr2v(t_ave, 1, i) = arr2v(t_ave, 1, i) + this->t1_meam[eltj] * rhoa0j; arr2v(t_ave, 2, i) = arr2v(t_ave, 2, i) + this->t2_meam[eltj] * rhoa0j; arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) + this->t3_meam[eltj] * rhoa0j; arr2v(t_ave, 1, j) = arr2v(t_ave, 1, j) + this->t1_meam[elti] * rhoa0i; arr2v(t_ave, 2, j) = arr2v(t_ave, 2, j) + this->t2_meam[elti] * rhoa0i; arr2v(t_ave, 3, j) = arr2v(t_ave, 3, j) + this->t3_meam[elti] * rhoa0i; } if (this->ialloy == 1) { arr2v(tsq_ave, 1, i) = arr2v(tsq_ave, 1, i) + this->t1_meam[eltj] * this->t1_meam[eltj] * rhoa0j; arr2v(tsq_ave, 2, i) = arr2v(tsq_ave, 2, i) + this->t2_meam[eltj] * this->t2_meam[eltj] * rhoa0j; arr2v(tsq_ave, 3, i) = arr2v(tsq_ave, 3, i) + this->t3_meam[eltj] * this->t3_meam[eltj] * rhoa0j; arr2v(tsq_ave, 1, j) = arr2v(tsq_ave, 1, j) + this->t1_meam[elti] * this->t1_meam[elti] * rhoa0i; arr2v(tsq_ave, 2, j) = arr2v(tsq_ave, 2, j) + this->t2_meam[elti] * this->t2_meam[elti] * rhoa0i; arr2v(tsq_ave, 3, j) = arr2v(tsq_ave, 3, j) + this->t3_meam[elti] * this->t3_meam[elti] * rhoa0i; } arr1v(arho2b, i) = arr1v(arho2b, i) + rhoa2j; arr1v(arho2b, j) = arr1v(arho2b, j) + rhoa2i; A1j = rhoa1j / rij; A2j = rhoa2j / rij2; A3j = rhoa3j / (rij2 * rij); A1i = rhoa1i / rij; A2i = rhoa2i / rij2; A3i = rhoa3i / (rij2 * rij); nv2 = 1; nv3 = 1; for (m = 1; m <= 3; m++) { arr2v(arho1, m, i) = arr2v(arho1, m, i) + A1j * delij[m]; arr2v(arho1, m, j) = arr2v(arho1, m, j) - A1i * delij[m]; arr2v(arho3b, m, i) = arr2v(arho3b, m, i) + rhoa3j * delij[m] / rij; arr2v(arho3b, m, j) = arr2v(arho3b, m, j) - rhoa3i * delij[m] / rij; for (n = m; n <= 3; n++) { arr2v(arho2, nv2, i) = arr2v(arho2, nv2, i) + A2j * delij[m] * delij[n]; arr2v(arho2, nv2, j) = arr2v(arho2, nv2, j) + A2i * delij[m] * delij[n]; nv2 = nv2 + 1; for (p = n; p <= 3; p++) { arr2v(arho3, nv3, i) = arr2v(arho3, nv3, i) + A3j * delij[m] * delij[n] * delij[p]; arr2v(arho3, nv3, j) = arr2v(arho3, nv3, j) - A3i * delij[m] * delij[n] * delij[p]; nv3 = nv3 + 1; } } } } } } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void 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) // j = atom 2 id (integer) // rijsq = squared distance between i and j // Outputs: sij = screening function { arrdim2v(x, 3, nmax) int k, nk /*,m*/; int elti, eltj, eltk; double delxik, delyik, delzik; double delxjk, delyjk, delzjk; double riksq, rjksq, xik, xjk, cikj, a, delc, sikj /*,fcij,rij*/; double Cmax, Cmin, rbound; *sij = 1.0; eltj = arr1v(fmap, arr1v(type, j)); elti = arr1v(fmap, arr1v(type, j)); // if rjksq > ebound*rijsq, atom k is definitely outside the ellipse rbound = this->ebound_meam[elti][eltj] * rijsq; for (nk = 1; nk <= numneigh_full; nk++) { k = arr1v(firstneigh_full, nk); eltk = arr1v(fmap, arr1v(type, k)); if (k == j) continue; delxjk = arr2v(x, 1, k) - arr2v(x, 1, j); delyjk = arr2v(x, 2, k) - arr2v(x, 2, j); delzjk = arr2v(x, 3, k) - arr2v(x, 3, j); rjksq = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk; if (rjksq > rbound) continue; delxik = arr2v(x, 1, k) - arr2v(x, 1, i); delyik = arr2v(x, 2, k) - arr2v(x, 2, i); delzik = arr2v(x, 3, k) - arr2v(x, 3, i); riksq = delxik * delxik + delyik * delyik + delzik * delzik; if (riksq > rbound) continue; xik = riksq / rijsq; xjk = rjksq / rijsq; a = 1 - (xik - xjk) * (xik - xjk); // if a < 0, then ellipse equation doesn't describe this case and // atom k can't possibly screen i-j if (a <= 0.0) continue; cikj = (2.0 * (xik + xjk) + a - 2.0) / a; Cmax = 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 // tolerance) if atoms are colinear, so don't reject that case here // (other negative cikj cases were handled by the test on "a" above) else if (cikj <= Cmin) { *sij = 0.0; break; } else { delc = Cmax - Cmin; cikj = (cikj - Cmin) / delc; fcut(cikj, &sikj); } *sij = *sij * sikj; } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void 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) { // Inputs: i,j,k = id's of 3 atom triplet // jn = id of i-j pair // rij2 = squared distance between i and j // Outputs: dsij1 = deriv. of sij w.r.t. rik // dsij2 = deriv. of sij w.r.t. rjk arrdim2v(x, 3, nmax) int elti, eltj, eltk; double rik2, rjk2; double dxik, dyik, dzik; double dxjk, dyjk, dzjk; double rbound, delc, sij, xik, xjk, cikj, sikj, dfc, a; double Cmax, Cmin, dCikj1, dCikj2; sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn); elti = arr1v(fmap, arr1v(type, i)); eltj = arr1v(fmap, arr1v(type, j)); eltk = arr1v(fmap, arr1v(type, k)); Cmax = 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 * 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); dzjk = arr2v(x, 3, k) - arr2v(x, 3, j); rjk2 = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; if (rjk2 <= rbound) { dxik = arr2v(x, 1, k) - arr2v(x, 1, i); dyik = arr2v(x, 2, k) - arr2v(x, 2, i); dzik = arr2v(x, 3, k) - arr2v(x, 3, i); rik2 = dxik * dxik + dyik * dyik + dzik * dzik; if (rik2 <= rbound) { xik = rik2 / rij2; xjk = rjk2 / rij2; a = 1 - (xik - xjk) * (xik - xjk); if (!iszero(a)) { cikj = (2.0 * (xik + xjk) + a - 2.0) / a; if (cikj >= Cmin && cikj <= Cmax) { cikj = (cikj - Cmin) / delc; dfcut(cikj, &sikj, &dfc); dCfunc2(rij2, rik2, rjk2, &dCikj1, &dCikj2); a = sij / delc * dfc / sikj; *dsij1 = a * dCikj1; *dsij2 = a * dCikj2; } } } } } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void MEAM::fcut(double xi, double* fc) { // cutoff function double a; if (xi >= 1.0) *fc = 1.0; else if (xi <= 0.0) *fc = 0.0; else { a = 1.0 - xi; a = a * a; a = a * a; a = 1.0 - a; *fc = a * a; // fc = xi } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void MEAM::dfcut(double xi, double* fc, double* dfc) { // cutoff function and its derivative double a, a3, a4; if (xi >= 1.0) { *fc = 1.0; *dfc = 0.0; } else if (xi <= 0.0) { *fc = 0.0; *dfc = 0.0; } else { a = 1.0 - xi; a3 = a * a * a; a4 = a * a3; *fc = pow((1.0 - a4), 2); *dfc = 8 * (1.0 - a4) * a3; // fc = xi // dfc = 1.d0 } } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void MEAM::dCfunc(double rij2, double rik2, double rjk2, double* dCikj) { // Inputs: rij,rij2,rik2,rjk2 // Outputs: dCikj = derivative of Cikj w.r.t. rij double rij4, a, b, denom; rij4 = rij2 * rij2; a = rik2 - rjk2; b = rik2 + rjk2; denom = rij4 - a * a; denom = denom * denom; *dCikj = -4 * (-2 * rij2 * a * a + rij4 * b + a * a * b) / denom; } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc void 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 // dCikj2 = derivative of Cikj w.r.t. rjk double rij4, rik4, rjk4, a, b, denom; rij4 = rij2 * rij2; rik4 = rik2 * rik2; rjk4 = rjk2 * rjk2; a = rik2 - rjk2; b = rik2 + rjk2; denom = rij4 - a * a; denom = denom * denom; *dCikj1 = 4 * rij2 * (rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom; *dCikj2 = 4 * rij2 * (rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom; (void)(b); } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index 80e34c66d..f5d4e7579 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -1,595 +1,596 @@ #include "meam.h" #include +using namespace LAMMPS_NS; // Extern "C" declaration has the form: // // void meam_force_(int *, int *, int *, double *, int *, int *, int *, double // *, // int *, int *, int *, int *, double *, double *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, double *, int //*); // // Call from pair_meam.cpp has the form: // // meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom, // &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], // &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], // &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], // dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, // &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], // &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); // void MEAM::meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, int* eflag_atom, int* vflag_atom, double* eng_vdwl, double* eatom, int* ntype, int* type, int* fmap, double* x, int* numneigh, int* firstneigh, int* numneigh_full, int* firstneigh_full, double* scrfcn, double* dscrfcn, double* fcpair, double* dGamma1, double* dGamma2, double* dGamma3, double* rho0, double* rho1, double* rho2, double* rho3, double* fp, double* Arho1, double* Arho2, double* Arho2b, double* Arho3, double* Arho3b, double* t_ave, double* tsq_ave, double* f, double* vatom, int* errorflag) { arrdim2v(x, 3, *nmax); arrdim2v(Arho1, 3, *nmax); arrdim2v(Arho2, 6, *nmax); arrdim2v(Arho3, 10, *nmax); arrdim2v(Arho3b, 3, *nmax); arrdim2v(t_ave, 3, *nmax); arrdim2v(tsq_ave, 3, *nmax); arrdim2v(f, 3, *nmax); arrdim2v(vatom, 6, *nmax); int i, j, jn, k, kn, kk, m, n, p, q; int nv2, nv3, elti, eltj, eltk, ind; double xitmp, yitmp, zitmp, delij[3 + 1], rij2, rij, rij3; double delik[3 + 1], deljk[3 + 1], v[6 + 1], fi[3 + 1], fj[3 + 1]; double third, sixth; double pp, dUdrij, dUdsij, dUdrijm[3 + 1], force, forcem; double r, recip, phi, phip; double sij; double a1, a1i, a1j, a2, a2i, a2j; double a3i, a3j; double shpi[3 + 1], shpj[3 + 1]; double ai, aj, ro0i, ro0j, invrei, invrej; double rhoa0j, drhoa0j, rhoa0i, drhoa0i; double rhoa1j, drhoa1j, rhoa1i, drhoa1i; double rhoa2j, drhoa2j, rhoa2i, drhoa2i; double a3, a3a, rhoa3j, drhoa3j, rhoa3i, drhoa3i; double drho0dr1, drho0dr2, drho0ds1, drho0ds2; double drho1dr1, drho1dr2, drho1ds1, drho1ds2; double drho1drm1[3 + 1], drho1drm2[3 + 1]; double drho2dr1, drho2dr2, drho2ds1, drho2ds2; double drho2drm1[3 + 1], drho2drm2[3 + 1]; double drho3dr1, drho3dr2, drho3ds1, drho3ds2; double drho3drm1[3 + 1], drho3drm2[3 + 1]; double dt1dr1, dt1dr2, dt1ds1, dt1ds2; double dt2dr1, dt2dr2, dt2ds1, dt2ds2; double dt3dr1, dt3dr2, dt3ds1, dt3ds2; double drhodr1, drhodr2, drhods1, drhods2, drhodrm1[3 + 1], drhodrm2[3 + 1]; double arg; double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3; double dsij1, dsij2, force1, force2; double t1i, t2i, t3i, t1j, t2j, t3j; *errorflag = 0; third = 1.0 / 3.0; sixth = 1.0 / 6.0; //: aliased i = *iptr; // Compute forces atom i elti = arr1v(fmap, arr1v(type, i)); if (elti > 0) { xitmp = arr2v(x, 1, i); yitmp = arr2v(x, 2, i); zitmp = arr2v(x, 3, i); // Treat each pair for (jn = 1; jn <= *numneigh; jn++) { j = arr1v(firstneigh, jn); eltj = arr1v(fmap, arr1v(type, j)); if (!iszero(arr1v(scrfcn, jn)) && eltj > 0) { sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn); delij[1] = arr2v(x, 1, j) - xitmp; delij[2] = arr2v(x, 2, j) - yitmp; delij[3] = arr2v(x, 3, j) - zitmp; rij2 = delij[1] * delij[1] + delij[2] * delij[2] + delij[3] * delij[3]; if (rij2 < this->cutforcesq) { rij = sqrt(rij2); r = rij; // Compute phi and phip ind = this->eltind[elti][eltj]; pp = rij * this->rdrar + 1.0; kk = (int)pp; - kk = min(kk, this->nrar - 1); + kk = MIN(kk, this->nrar - 1); pp = pp - kk; - pp = min(pp, 1.0); + pp = MIN(pp, 1.0); phi = ((arr2(this->phirar3, kk, ind) * pp + arr2(this->phirar2, kk, ind)) * pp + arr2(this->phirar1, kk, ind)) * pp + arr2(this->phirar, kk, ind); phip = (arr2(this->phirar6, kk, ind) * pp + arr2(this->phirar5, kk, ind)) * pp + arr2(this->phirar4, kk, ind); recip = 1.0 / r; if (*eflag_either != 0) { if (*eflag_global != 0) *eng_vdwl = *eng_vdwl + phi * sij; if (*eflag_atom != 0) { arr1v(eatom, i) = arr1v(eatom, i) + 0.5 * phi * sij; arr1v(eatom, j) = arr1v(eatom, j) + 0.5 * phi * sij; } } // write(1,*) "force_meamf: phi: ",phi // write(1,*) "force_meamf: phip: ",phip // Compute pair densities and derivatives invrei = 1.0 / this->re_meam[elti][elti]; ai = rij * invrei - 1.0; 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 / this->re_meam[eltj][eltj]; aj = rij * invrej - 1.0; 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; rhoa1j = rhoa1i; drhoa1j = drhoa1i; rhoa2j = rhoa2i; drhoa2j = drhoa2i; rhoa3j = rhoa3i; drhoa3j = drhoa3i; } 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; nv3 = 1; arg1i1 = 0.0; arg1j1 = 0.0; arg1i2 = 0.0; arg1j2 = 0.0; arg1i3 = 0.0; arg1j3 = 0.0; arg3i3 = 0.0; arg3j3 = 0.0; for (n = 1; n <= 3; n++) { for (p = n; p <= 3; p++) { for (q = p; q <= 3; q++) { arg = delij[n] * delij[p] * delij[q] * 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] * this->v2D[nv2]; arg1i2 = arg1i2 + arr2v(Arho2, nv2, i) * arg; arg1j2 = arg1j2 + arr2v(Arho2, nv2, j) * arg; nv2 = nv2 + 1; } arg1i1 = arg1i1 + arr2v(Arho1, n, i) * delij[n]; arg1j1 = arg1j1 - arr2v(Arho1, n, j) * delij[n]; arg3i3 = arg3i3 + arr2v(Arho3b, n, i) * delij[n]; arg3j3 = arg3j3 - arr2v(Arho3b, n, j) * delij[n]; } // rho0 terms drho0dr1 = drhoa0j * sij; drho0dr2 = drhoa0i * sij; // rho1 terms a1 = 2 * sij / rij; drho1dr1 = a1 * (drhoa1j - rhoa1j / rij) * arg1i1; drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1; a1 = 2.0 * sij / rij; for (m = 1; m <= 3; m++) { drho1drm1[m] = a1 * rhoa1j * arr2v(Arho1, m, i); drho1drm2[m] = -a1 * rhoa1i * arr2v(Arho1, m, j); } // rho2 terms a2 = 2 * sij / rij2; drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - 2.0 / 3.0 * arr1v(Arho2b, i) * drhoa2j * sij; drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - 2.0 / 3.0 * arr1v(Arho2b, j) * drhoa2i * sij; a2 = 4 * sij / rij2; for (m = 1; m <= 3; m++) { drho2drm1[m] = 0.0; drho2drm2[m] = 0.0; for (n = 1; n <= 3; n++) { drho2drm1[m] = drho2drm1[m] + arr2v(Arho2, this->vind2D[m][n], i) * delij[n]; drho2drm2[m] = drho2drm2[m] - arr2v(Arho2, this->vind2D[m][n], j) * delij[n]; } drho2drm1[m] = a2 * rhoa2j * drho2drm1[m]; drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m]; } // rho3 terms rij3 = rij * rij2; a3 = 2 * sij / rij3; a3a = 6.0 / 5.0 * sij / rij; drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 - a3a * (drhoa3j - rhoa3j / rij) * arg3i3; drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 - a3a * (drhoa3i - rhoa3i / rij) * arg3j3; a3 = 6 * sij / rij3; a3a = 6 * sij / (5 * rij); for (m = 1; m <= 3; m++) { drho3drm1[m] = 0.0; drho3drm2[m] = 0.0; nv2 = 1; for (n = 1; n <= 3; n++) { for (p = n; p <= 3; p++) { arg = delij[n] * delij[p] * this->v2D[nv2]; drho3drm1[m] = drho3drm1[m] + arr2v(Arho3, this->vind3D[m][n][p], i) * arg; drho3drm2[m] = drho3drm2[m] + arr2v(Arho3, this->vind3D[m][n][p], j) * arg; nv2 = nv2 + 1; } } drho3drm1[m] = (a3 * drho3drm1[m] - a3a * arr2v(Arho3b, m, i)) * rhoa3j; drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * arr2v(Arho3b, m, j)) * rhoa3i; } // Compute derivatives of weighting functions t wrt rij t1i = arr2v(t_ave, 1, i); t2i = arr2v(t_ave, 2, i); t3i = arr2v(t_ave, 3, i); t1j = arr2v(t_ave, 1, j); t2j = arr2v(t_ave, 2, j); t3j = arr2v(t_ave, 3, j); if (this->ialloy == 1) { a1i = 0.0; a1j = 0.0; a2i = 0.0; a2j = 0.0; a3i = 0.0; a3j = 0.0; if (!iszero(arr2v(tsq_ave, 1, i))) a1i = drhoa0j * sij / arr2v(tsq_ave, 1, i); if (!iszero(arr2v(tsq_ave, 1, j))) a1j = drhoa0i * sij / arr2v(tsq_ave, 1, j); if (!iszero(arr2v(tsq_ave, 2, i))) a2i = drhoa0j * sij / arr2v(tsq_ave, 2, i); if (!iszero(arr2v(tsq_ave, 2, j))) a2j = drhoa0i * sij / arr2v(tsq_ave, 2, j); if (!iszero(arr2v(tsq_ave, 3, i))) a3i = drhoa0j * sij / arr2v(tsq_ave, 3, i); if (!iszero(arr2v(tsq_ave, 3, j))) a3j = drhoa0i * sij / arr2v(tsq_ave, 3, j); dt1dr1 = a1i * (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 (this->ialloy == 2) { dt1dr1 = 0.0; dt1dr2 = 0.0; dt2dr1 = 0.0; dt2dr2 = 0.0; dt3dr1 = 0.0; dt3dr2 = 0.0; } else { ai = 0.0; if (!iszero(arr1v(rho0, i))) ai = drhoa0j * sij / arr1v(rho0, i); aj = 0.0; if (!iszero(arr1v(rho0, j))) aj = drhoa0i * sij / arr1v(rho0, j); dt1dr1 = ai * (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, 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 + dt2dr1 * arr1v(rho2, i) + t2i * drho2dr1 + dt3dr1 * arr1v(rho3, i) + t3i * drho3dr1) - arr1v(dGamma3, i) * (shpi[1] * dt1dr1 + shpi[2] * dt2dr1 + shpi[3] * dt3dr1); drhodr2 = arr1v(dGamma1, j) * drho0dr2 + arr1v(dGamma2, j) * (dt1dr2 * arr1v(rho1, j) + t1j * drho1dr2 + dt2dr2 * arr1v(rho2, j) + t2j * drho2dr2 + dt3dr2 * arr1v(rho3, j) + t3j * drho3dr2) - arr1v(dGamma3, j) * (shpj[1] * dt1dr2 + shpj[2] * dt2dr2 + shpj[3] * dt3dr2); for (m = 1; m <= 3; m++) { drhodrm1[m] = 0.0; drhodrm2[m] = 0.0; drhodrm1[m] = arr1v(dGamma2, i) * (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]); drhodrm2[m] = arr1v(dGamma2, j) * (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]); } // Compute derivatives wrt sij, but only if necessary if (!iszero(arr1v(dscrfcn, jn))) { drho0ds1 = rhoa0j; drho0ds2 = rhoa0i; a1 = 2.0 / rij; drho1ds1 = a1 * rhoa1j * arg1i1; drho1ds2 = a1 * rhoa1i * arg1j1; a2 = 2.0 / rij2; drho2ds1 = a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arr1v(Arho2b, i) * rhoa2j; drho2ds2 = a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arr1v(Arho2b, j) * rhoa2i; a3 = 2.0 / rij3; a3a = 6.0 / (5.0 * rij); drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3; drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3; if (this->ialloy == 1) { a1i = 0.0; a1j = 0.0; a2i = 0.0; a2j = 0.0; a3i = 0.0; a3j = 0.0; if (!iszero(arr2v(tsq_ave, 1, i))) a1i = rhoa0j / arr2v(tsq_ave, 1, i); if (!iszero(arr2v(tsq_ave, 1, j))) a1j = rhoa0i / arr2v(tsq_ave, 1, j); if (!iszero(arr2v(tsq_ave, 2, i))) a2i = rhoa0j / arr2v(tsq_ave, 2, i); if (!iszero(arr2v(tsq_ave, 2, j))) a2j = rhoa0i / arr2v(tsq_ave, 2, j); if (!iszero(arr2v(tsq_ave, 3, i))) a3i = rhoa0j / arr2v(tsq_ave, 3, i); if (!iszero(arr2v(tsq_ave, 3, j))) a3j = rhoa0i / arr2v(tsq_ave, 3, j); dt1ds1 = a1i * (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 (this->ialloy == 2) { dt1ds1 = 0.0; dt1ds2 = 0.0; dt2ds1 = 0.0; dt2ds2 = 0.0; dt3ds1 = 0.0; dt3ds2 = 0.0; } else { ai = 0.0; if (!iszero(arr1v(rho0, i))) ai = rhoa0j / arr1v(rho0, i); aj = 0.0; if (!iszero(arr1v(rho0, j))) aj = rhoa0i / arr1v(rho0, j); dt1ds1 = ai * (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 = arr1v(dGamma1, i) * drho0ds1 + arr1v(dGamma2, i) * (dt1ds1 * arr1v(rho1, i) + t1i * drho1ds1 + dt2ds1 * arr1v(rho2, i) + t2i * drho2ds1 + dt3ds1 * arr1v(rho3, i) + t3i * drho3ds1) - arr1v(dGamma3, i) * (shpi[1] * dt1ds1 + shpi[2] * dt2ds1 + shpi[3] * dt3ds1); drhods2 = arr1v(dGamma1, j) * drho0ds2 + arr1v(dGamma2, j) * (dt1ds2 * arr1v(rho1, j) + t1j * drho1ds2 + dt2ds2 * arr1v(rho2, j) + t2j * drho2ds2 + dt3ds2 * arr1v(rho3, j) + t3j * drho3ds2) - arr1v(dGamma3, j) * (shpj[1] * dt1ds2 + shpj[2] * dt2ds2 + shpj[3] * dt3ds2); } // Compute derivatives of energy wrt rij, sij and rij[3] dUdrij = phip * sij + arr1v(fp, i) * drhodr1 + arr1v(fp, j) * drhodr2; dUdsij = 0.0; if (!iszero(arr1v(dscrfcn, jn))) { dUdsij = phi + arr1v(fp, i) * drhods1 + arr1v(fp, j) * drhods2; } for (m = 1; m <= 3; m++) { dUdrijm[m] = arr1v(fp, i) * drhodrm1[m] + arr1v(fp, j) * drhodrm2[m]; } // Add the part of the force due to dUdrij and dUdsij force = dUdrij * recip + dUdsij * arr1v(dscrfcn, jn); for (m = 1; m <= 3; m++) { forcem = delij[m] * force + dUdrijm[m]; arr2v(f, m, i) = arr2v(f, m, i) + forcem; arr2v(f, m, j) = arr2v(f, m, j) - forcem; } // Tabulate per-atom virial as symmetrized stress tensor if (*vflag_atom != 0) { fi[1] = delij[1] * force + dUdrijm[1]; fi[2] = delij[2] * force + dUdrijm[2]; fi[3] = delij[3] * force + dUdrijm[3]; v[1] = -0.5 * (delij[1] * fi[1]); v[2] = -0.5 * (delij[2] * fi[2]); v[3] = -0.5 * (delij[3] * fi[3]); v[4] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]); v[5] = -0.25 * (delij[1] * fi[3] + delij[3] * fi[1]); v[6] = -0.25 * (delij[2] * fi[3] + delij[3] * fi[2]); arr2v(vatom, 1, i) = arr2v(vatom, 1, i) + v[1]; arr2v(vatom, 2, i) = arr2v(vatom, 2, i) + v[2]; arr2v(vatom, 3, i) = arr2v(vatom, 3, i) + v[3]; arr2v(vatom, 4, i) = arr2v(vatom, 4, i) + v[4]; arr2v(vatom, 5, i) = arr2v(vatom, 5, i) + v[5]; arr2v(vatom, 6, i) = arr2v(vatom, 6, i) + v[6]; arr2v(vatom, 1, j) = arr2v(vatom, 1, j) + v[1]; arr2v(vatom, 2, j) = arr2v(vatom, 2, j) + v[2]; arr2v(vatom, 3, j) = arr2v(vatom, 3, j) + v[3]; arr2v(vatom, 4, j) = arr2v(vatom, 4, j) + v[4]; arr2v(vatom, 5, j) = arr2v(vatom, 5, j) + v[5]; arr2v(vatom, 6, j) = arr2v(vatom, 6, j) + v[6]; } // Now compute forces on other atoms k due to change in sij if (iszero(sij) || iszero(sij - 1.0)) continue; //: cont jn loop for (kn = 1; kn <= *numneigh_full; kn++) { k = arr1v(firstneigh_full, kn); eltk = arr1v(fmap, arr1v(type, k)); if (k != j && eltk > 0) { dsij(i, j, k, jn, *nmax, *numneigh, rij2, &dsij1, &dsij2, *ntype, type, fmap, x, scrfcn, fcpair); if (!iszero(dsij1) || !iszero(dsij2)) { force1 = dUdsij * dsij1; force2 = dUdsij * dsij2; for (m = 1; m <= 3; m++) { delik[m] = arr2v(x, m, k) - arr2v(x, m, i); deljk[m] = arr2v(x, m, k) - arr2v(x, m, j); } for (m = 1; m <= 3; m++) { arr2v(f, m, i) = arr2v(f, m, i) + force1 * delik[m]; arr2v(f, m, j) = arr2v(f, m, j) + force2 * deljk[m]; arr2v(f, m, k) = arr2v(f, m, k) - force1 * delik[m] - force2 * deljk[m]; } // Tabulate per-atom virial as symmetrized stress tensor if (*vflag_atom != 0) { fi[1] = force1 * delik[1]; fi[2] = force1 * delik[2]; fi[3] = force1 * delik[3]; fj[1] = force2 * deljk[1]; fj[2] = force2 * deljk[2]; fj[3] = force2 * deljk[3]; v[1] = -third * (delik[1] * fi[1] + deljk[1] * fj[1]); v[2] = -third * (delik[2] * fi[2] + deljk[2] * fj[2]); v[3] = -third * (delik[3] * fi[3] + deljk[3] * fj[3]); v[4] = -sixth * (delik[1] * fi[2] + deljk[1] * fj[2] + delik[2] * fi[1] + deljk[2] * fj[1]); v[5] = -sixth * (delik[1] * fi[3] + deljk[1] * fj[3] + delik[3] * fi[1] + deljk[3] * fj[1]); v[6] = -sixth * (delik[2] * fi[3] + deljk[2] * fj[3] + delik[3] * fi[2] + deljk[3] * fj[2]); arr2v(vatom, 1, i) = arr2v(vatom, 1, i) + v[1]; arr2v(vatom, 2, i) = arr2v(vatom, 2, i) + v[2]; arr2v(vatom, 3, i) = arr2v(vatom, 3, i) + v[3]; arr2v(vatom, 4, i) = arr2v(vatom, 4, i) + v[4]; arr2v(vatom, 5, i) = arr2v(vatom, 5, i) + v[5]; arr2v(vatom, 6, i) = arr2v(vatom, 6, i) + v[6]; arr2v(vatom, 1, j) = arr2v(vatom, 1, j) + v[1]; arr2v(vatom, 2, j) = arr2v(vatom, 2, j) + v[2]; arr2v(vatom, 3, j) = arr2v(vatom, 3, j) + v[3]; arr2v(vatom, 4, j) = arr2v(vatom, 4, j) + v[4]; arr2v(vatom, 5, j) = arr2v(vatom, 5, j) + v[5]; arr2v(vatom, 6, j) = arr2v(vatom, 6, j) + v[6]; arr2v(vatom, 1, k) = arr2v(vatom, 1, k) + v[1]; arr2v(vatom, 2, k) = arr2v(vatom, 2, k) + v[2]; arr2v(vatom, 3, k) = arr2v(vatom, 3, k) + v[3]; arr2v(vatom, 4, k) = arr2v(vatom, 4, k) + v[4]; arr2v(vatom, 5, k) = arr2v(vatom, 5, k) + v[5]; arr2v(vatom, 6, k) = arr2v(vatom, 6, k) + v[6]; } } } // end of k loop } } } // end of j loop } // else if elti=0, this is not a meam atom } } diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp index b05e3824b..bea6962de 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -1,1055 +1,1056 @@ #include "meam.h" #include +using namespace LAMMPS_NS; // Declaration in pair_meam.h: // // void meam_setup_done(double *) // // Call from pair_meam.cpp: // // meam_setup_done(&cutmax) // void MEAM::meam_setup_done_(double* cutmax) { int nv2, nv3, m, n, p; // Force cutoff this->cutforce = this->rc_meam; this->cutforcesq = this->cutforce * this->cutforce; // Pass cutoff back to calling program *cutmax = this->cutforce; // Augment t1 term for (int i = 1; i <= maxelt; 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(); // indices and factors for Voight notation nv2 = 1; nv3 = 1; for (m = 1; m <= 3; m++) { for (n = m; n <= 3; n++) { this->vind2D[m][n] = nv2; this->vind2D[n][m] = nv2; nv2 = nv2 + 1; for (p = n; p <= 3; p++) { 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; } } } this->v2D[1] = 1; this->v2D[2] = 2; this->v2D[3] = 2; this->v2D[4] = 1; this->v2D[5] = 2; this->v2D[6] = 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 <= this->neltypes; m++) { for (n = m; n <= this->neltypes; n++) { this->eltind[m][n] = nv2; this->eltind[n][m] = nv2; nv2 = nv2 + 1; } } // Compute background densities for reference structure compute_reference_density(); // Compute pair potentials and setup arrays for interpolation this->nr = 1000; this->dr = 1.1 * this->rc_meam / this->nr; compute_pair_meam(); } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc // Fill off-diagonal alloy parameters void MEAM::alloyparams(void) { int i, j, k; double eb; // Loop over pairs 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) { 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(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 this->Ec_meam[i][j] = (this->Ec_meam[i][i] + 2 * this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j]; } else this->Ec_meam[i][j] = (this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 2.0 - this->delta_meam[i][j]; } 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 <= 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]; } } } // ebound gives the squared distance such that, for rik2 or rjk2>ebound, // atom k definitely lies outside the screening function ellipse (so // there is no need to calculate its effects). Here, compute it for all // triplets [i][j][k] so that ebound[i][j] is the maximized over k for (i = 2; i <= 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); + this->ebound_meam[i][j] = MAX(this->ebound_meam[i][j], eb); } } } } //----------------------------------------------------------------------- // compute MEAM pair potential for each pair of element types // void MEAM::compute_pair_meam(void) { double r /*ununsed:, temp*/; int j, a, b, nv2; double astar, frac, phizbl; int n, nmax, Z1, Z2; double arat, rarat, scrn, scrn2; double phiaa, phibb /*unused:,phitmp*/; double C, s111, s112, s221, S11, S22; // check for previously allocated arrays and free them if (allocated(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(this->phir, this->nr, (this->neltypes * (this->neltypes + 1)) / 2); // allocate coeff memory 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 <= this->neltypes; a++) { for (b = a; b <= this->neltypes; b++) { nv2 = nv2 + 1; // loop over r values and compute for (j = 1; j <= this->nr; j++) { r = (j - 1) * this->dr; 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 (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 // compute the contributions from second nearest neighbors, like // a-a // pairs, but need to include NN2 contributions to those pairs as // well. 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, 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++) { phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a); } } // phi_bb phibb = phi_meam(rarat, b, b); get_Zij(&Z1, 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++) { phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b); } } if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2) { // Add contributions to the B1 or B2 potential 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 (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 // second-neighbor pairs are screened by 2 type 1 atoms and // two type // 2 atoms. 2-2 second-neighbor pairs are screened by 4 type // 1 // atoms. C = 1.0; get_sijk(C, a, a, a, &s111); get_sijk(C, a, a, b, &s112); get_sijk(C, b, b, a, &s221); S11 = s111 * s111 * s112 * s112; S22 = pow(s221, 4); arr2(this->phir, j, nv2) = arr2(this->phir, j, nv2) - 0.75 * S11 * phiaa - 0.25 * S22 * phibb; } } else { nmax = 10; for (n = 1; n <= nmax; n++) { arr2(this->phir, j, nv2) = arr2(this->phir, j, nv2) + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b); } } } // For Zbl potential: // if astar <= -3 // potential is zbl potential // else if -3 < astar < -1 // potential is linear combination with zbl potential // endif if (this->zbl_meam[a][b] == 1) { astar = this->alpha_meam[a][b] * (r / this->re_meam[a][b] - 1.0); if (astar <= -3.0) 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, this->ielt_meam[a], this->ielt_meam[b]); arr2(this->phir, j, nv2) = frac * arr2(this->phir, j, nv2) + (1 - frac) * phizbl; } } } // call interpolation interpolate_meam(nv2); } } } //----------------------------------------------------------------------c // Compute MEAM pair potential for distance r, element types a and b // double MEAM::phi_meam(double r, int a, int b) { /*unused:double a1,a2,a12;*/ double t11av, t21av, t31av, t12av, t22av, t32av; double G1, G2, s1[3 + 1], s2[3 + 1] /*,s12[3+1]*/, rho0_1, rho0_2; double Gam1, Gam2, Z1, Z2; double rhobar1, rhobar2, F1, F2; double rho01, rho11, rho21, rho31; double rho02, rho12, rho22, rho32; double scalfac, phiaa, phibb; double Eu; double arat, scrn /*unused:,scrn2*/; int Z12, errorflag; int n, nmax, Z1nn, Z2nn; lattice_t latta /*unused:,lattb*/; double rho_bkgd1, rho_bkgd2; double phi_m = 0.0; // Equation numbers below refer to: // I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615 // get number of neighbors in the reference structure // Nref[i][j] = # of i's neighbors of type j get_Zij(&Z12, this->lattce_meam[a][b]); get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32); // if densities are too small, numerical problems may result; just return zero if (rho01 <= 1e-14 && rho02 <= 1e-14) return 0.0; // calculate average weighting factors for the reference structure if (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 * (this->t1_meam[a] * rho01 + this->t1_meam[b] * rho02); t12av = t11av; t21av = scalfac * (this->t2_meam[a] * rho01 + this->t2_meam[b] * rho02); t22av = t21av; t31av = 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, 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 (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) + t21av / 6.0 * pow(rho22 + rho21, 2) + 121.0 / 40.0 * t31av * pow((rho32 - rho31), 2); rhobar1 = sqrt(rhobar1); rhobar2 = pow(Z12 * rho01, 2) + 2.0 / 3.0 * t21av * pow(rho21, 2); rhobar2 = sqrt(rhobar2); } else { rhobar2 = pow(((Z12 / 2) * (rho01 + rho02)), 2) + t12av * pow((rho11 - rho12), 2) + t22av / 6.0 * pow(rho21 + rho22, 2) + 121.0 / 40.0 * t32av * pow((rho31 - rho32), 2); rhobar2 = sqrt(rhobar2); rhobar1 = pow(Z12 * rho02, 2) + 2.0 / 3.0 * t22av * pow(rho22, 2); rhobar1 = sqrt(rhobar1); } } else { // for other structures, use formalism developed in Huang's paper // // composition-dependent scaling, equation I.7 // If using mixing rule for t, apply to reference structure; else // use precomputed values if (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, this->lattce_meam[a][a]); Gam1 = (s1[1] * t11av + s1[2] * t21av + s1[3] * t31av) / (Z1 * Z1); G_gam(Gam1, this->ibar_meam[a], this->gsmooth_factor, &G1, &errorflag); } if (this->ibar_meam[b] <= 0) G2 = 1.0; else { get_shpfcn(s2, this->lattce_meam[b][b]); Gam2 = (s2[1] * t12av + s2[2] * t22av + s2[3] * t32av) / (Z2 * Z2); G_gam(Gam2, this->ibar_meam[b], this->gsmooth_factor, &G2, &errorflag); } 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) Gam1 = 0.0; else Gam1 = Gam1 / (rho01 * rho01); Gam2 = (t12av * rho12 + t22av * rho22 + t32av * rho32); if (rho02 < 1.0e-14) Gam2 = 0.0; else Gam2 = Gam2 / (rho02 * rho02); G_gam(Gam1, this->ibar_meam[a], this->gsmooth_factor, &G1, &errorflag); G_gam(Gam2, this->ibar_meam[b], this->gsmooth_factor, &G2, &errorflag); if (this->mix_ref_t == 1) { rho_bkgd1 = rho0_1; rho_bkgd2 = rho0_2; } else { if (this->bkgd_dyn == 1) { rho_bkgd1 = this->rho0_meam[a] * this->Z_meam[a]; rho_bkgd2 = this->rho0_meam[b] * this->Z_meam[b]; } else { rho_bkgd1 = this->rho_ref_meam[a]; rho_bkgd2 = this->rho_ref_meam[b]; } } rhobar1 = rho01 / rho_bkgd1 * G1; rhobar2 = rho02 / rho_bkgd2 * G2; } // compute embedding functions, eqn I.5 if (iszero(rhobar1)) F1 = 0.0; else { if (this->emb_lin_neg == 1 && rhobar1 <= 0) F1 = -this->A_meam[a] * this->Ec_meam[a][a] * rhobar1; else F1 = this->A_meam[a] * this->Ec_meam[a][a] * rhobar1 * log(rhobar1); } if (iszero(rhobar2)) F2 = 0.0; else { if (this->emb_lin_neg == 1 && rhobar2 <= 0) F2 = -this->A_meam[b] * this->Ec_meam[b][b] * rhobar2; else F2 = this->A_meam[b] * this->Ec_meam[b][b] * rhobar2 * log(rhobar2); } // compute Rose function, I.16 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 (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; } else { phibb = phi_meam(r, b, b); phi_m = (3 * Eu - F1 - 2 * F2 - 5 * phibb) / Z12; } } 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, 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++) { phiaa = phiaa + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a); } } phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa; } else { // // potential is computed from Rose function and embedding energy phi_m = (2 * Eu - F1 - F2) / Z12; // } // if r = 0, just return 0 if (iszero(r)) { phi_m = 0.0; } return phi_m; } //----------------------------------------------------------------------c // Compute background density for reference structure of each element void 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 <= this->neltypes; a++) { Z = (int)this->Z_meam[a]; if (this->ibar_meam[a] <= 0) Gbar = 1.0; else { 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, 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 = 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 (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 = this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * (arat - 1)); rho0 = rho0 + Z2 * rho0_2nn * scrn; } this->rho_ref_meam[a] = rho0 * Gbar; } } //----------------------------------------------------------------------c // Shape factors for various configurations void MEAM::get_shpfcn(double* s /* s(3) */, lattice_t latt) { if (latt == FCC || latt == BCC || latt == B1 || latt == B2) { s[1] = 0.0; s[2] = 0.0; s[3] = 0.0; } else if (latt == HCP) { s[1] = 0.0; s[2] = 0.0; s[3] = 1.0 / 3.0; } else if (latt == DIA) { s[1] = 0.0; s[2] = 0.0; s[3] = 32.0 / 9.0; } else if (latt == DIM) { s[1] = 1.0; s[2] = 2.0 / 3.0; // s(3) = 1.d0 s[3] = 0.40; } else { s[1] = 0.0; // call error('Lattice not defined in get_shpfcn.') } } //------------------------------------------------------------------------------c // Average weighting factors for the reference structure void 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) { double rhoa01, rhoa02, a1, a2, rho01 /*,rho02*/; // For ialloy = 2, no averaging is done if (this->ialloy == 2) { *t11av = t11; *t21av = t21; *t31av = t31; *t12av = t12; *t22av = t22; *t32av = t32; } else { if (latt == FCC || latt == BCC || latt == DIA || latt == HCP || latt == B1 || latt == DIM || latt == B2) { // all neighbors are of the opposite type *t11av = t12; *t21av = t22; *t31av = t32; *t12av = t11; *t22av = t21; *t32av = t31; } else { a1 = r / 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; *t12av = t11; *t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01; *t22av = t21; *t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01; *t32av = t31; } else { // call error('Lattice not defined in get_tavref.') } } } } //------------------------------------------------------------------------------c // Number of neighbors for the reference structure void MEAM::get_Zij(int* Zij, lattice_t latt) { if (latt == FCC) *Zij = 12; else if (latt == BCC) *Zij = 8; else if (latt == HCP) *Zij = 12; else if (latt == B1) *Zij = 6; else if (latt == DIA) *Zij = 4; else if (latt == DIM) *Zij = 1; else if (latt == C11) *Zij = 10; else if (latt == L12) *Zij = 12; else if (latt == B2) *Zij = 8; else { // call error('Lattice not defined in get_Zij.') } } //------------------------------------------------------------------------------c // Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second // neighbor screening function for lattice type "latt" void MEAM::get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin, double cmax) { double /*rratio,*/ C, x, sijk; int numscr = 0; if (latt == BCC) { *Zij2 = 6; *a = 2.0 / sqrt(3.0); numscr = 4; } else if (latt == FCC) { *Zij2 = 6; *a = sqrt(2.0); numscr = 4; } else if (latt == DIA) { *Zij2 = 0; *a = sqrt(8.0 / 3.0); numscr = 4; if (cmin < 0.500001) { // call error('can not do 2NN MEAM for dia') } } else if (latt == HCP) { *Zij2 = 6; *a = sqrt(2.0); numscr = 4; } else if (latt == B1) { *Zij2 = 12; *a = sqrt(2.0); numscr = 2; } else if (latt == L12) { *Zij2 = 6; *a = sqrt(2.0); numscr = 4; } else if (latt == B2) { *Zij2 = 6; *a = 2.0 / sqrt(3.0); numscr = 4; } else if (latt == DIM) { // this really shouldn't be allowed; make sure screening is zero *Zij2 = 0; *a = 1; *S = 0; return; } else { // call error('Lattice not defined in get_Zij2.') } // Compute screening for each first neighbor C = 4.0 / (*a * *a) - 1.0; x = (C - cmin) / (cmax - cmin); fcut(x, &sijk); // There are numscr first neighbors screening the second neighbors *S = pow(sijk, numscr); } //------------------------------------------------------------------------------c void MEAM::get_sijk(double C, int i, int j, int k, double* sijk) { double x; 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 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) { double a1, a2; double s[3 + 1]; lattice_t lat; int Zij1nn, Zij2nn; double rhoa01nn, rhoa02nn; double rhoa01, rhoa11, rhoa21, rhoa31; double rhoa02, rhoa12, rhoa22, rhoa32; double arat, scrn, denom; double C, s111, s112, s221, S11, S22; a1 = r / 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); 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 = this->lattce_meam[a][b]; *rho11 = 0.0; *rho21 = 0.0; *rho31 = 0.0; *rho12 = 0.0; *rho22 = 0.0; *rho32 = 0.0; get_Zij(&Zij1nn, lat); if (lat == FCC) { *rho01 = 12.0 * rhoa02; *rho02 = 12.0 * rhoa01; } else if (lat == BCC) { *rho01 = 8.0 * rhoa02; *rho02 = 8.0 * rhoa01; } else if (lat == B1) { *rho01 = 6.0 * rhoa02; *rho02 = 6.0 * rhoa01; } else if (lat == DIA) { *rho01 = 4.0 * rhoa02; *rho02 = 4.0 * rhoa01; *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32; *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31; } else if (lat == HCP) { *rho01 = 12 * rhoa02; *rho02 = 12 * rhoa01; *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32; *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31; } else if (lat == DIM) { get_shpfcn(s, DIM); *rho01 = rhoa02; *rho02 = rhoa01; *rho11 = s[1] * rhoa12 * rhoa12; *rho12 = s[1] * rhoa11 * rhoa11; *rho21 = s[2] * rhoa22 * rhoa22; *rho22 = s[2] * rhoa21 * rhoa21; *rho31 = s[3] * rhoa32 * rhoa32; *rho32 = s[3] * rhoa31 * rhoa31; } else if (lat == C11) { *rho01 = rhoa01; *rho02 = rhoa02; *rho11 = rhoa11; *rho12 = rhoa12; *rho21 = rhoa21; *rho22 = rhoa22; *rho31 = rhoa31; *rho32 = rhoa32; } else if (lat == L12) { *rho01 = 8 * rhoa01 + 4 * rhoa02; *rho02 = 12 * rhoa01; if (this->ialloy == 1) { *rho21 = 8. / 3. * 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 *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22); } else if (lat == B2) { *rho01 = 8.0 * rhoa02; *rho02 = 8.0 * rhoa01; } else { // call error('Lattice not defined in get_densref.') } if (this->nn2_meam[a][b] == 1) { get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b]); a1 = arat * r / this->re_meam[a][a] - 1.0; a2 = arat * r / this->re_meam[b][b] - 1.0; 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 // the screening functions C = 1.0; get_sijk(C, a, a, a, &s111); get_sijk(C, a, a, b, &s112); get_sijk(C, b, b, a, &s221); S11 = s111 * s111 * s112 * s112; S22 = pow(s221, 4); *rho01 = *rho01 + 6 * S11 * rhoa01nn; *rho02 = *rho02 + 6 * S22 * rhoa02nn; } else { // For other cases, assume that second neighbor is of same type, // first neighbor may be of different type *rho01 = *rho01 + Zij2nn * scrn * rhoa01nn; // Assume Zij2nn and arat don't depend on order, but scrn might get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a]); *rho02 = *rho02 + Zij2nn * scrn * rhoa02nn; } } } //--------------------------------------------------------------------- // Compute ZBL potential // double MEAM::zbl(double r, int z1, int z2) { int i; const double c[] = { 0.028171, 0.28022, 0.50986, 0.18175 }; const double d[] = { 0.20162, 0.40290, 0.94229, 3.1998 }; const double azero = 0.4685; const double cc = 14.3997; double a, x; // azero = (9pi^2/128)^1/3 (0.529) Angstroms a = azero / (pow(z1, 0.23) + pow(z2, 0.23)); double result = 0.0; x = r / a; for (i = 0; i <= 3; i++) { result = result + c[i] * fm_exp(-d[i] * x); } if (r > 0.0) result = result * z1 * z2 / r * cc; return result; } //--------------------------------------------------------------------- // Compute Rose energy function, I.16 // double MEAM::erose(double r, double re, double alpha, double Ec, double repuls, double attrac, int form) { double astar, a3; double result = 0.0; if (r > 0.0) { astar = alpha * (r / re - 1.0); a3 = 0.0; if (astar >= 0) a3 = attrac; else if (astar < 0) a3 = repuls; if (form == 1) result = -Ec * (1 + astar + (-attrac + repuls / r) * pow(astar, 3)) * fm_exp(-astar); else if (form == 2) result = -Ec * (1 + astar + a3 * pow(astar, 3)) * fm_exp(-astar); else result = -Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * fm_exp(-astar); } return result; } // ----------------------------------------------------------------------- void MEAM::interpolate_meam(int ind) { int j; double drar; // map to coefficient space this->nrar = this->nr; drar = this->dr; this->rdrar = 1.0 / drar; // phir interp for (j = 1; j <= this->nrar; j++) { arr2(this->phirar, j, ind) = arr2(this->phir, j, 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 <= this->nrar - 1; j++) { arr2(this->phirar2, j, ind) = 3.0 * (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(this->phirar, j + 1, ind) - arr2(this->phirar, j, ind)); } arr2(this->phirar2, this->nrar, ind) = 0.0; arr2(this->phirar3, this->nrar, ind) = 0.0; 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; } } //--------------------------------------------------------------------- // Compute Rose energy function, I.16 // double MEAM::compute_phi(double rij, int elti, int eltj) { double pp; int ind, kk; ind = this->eltind[elti][eltj]; pp = rij * this->rdrar + 1.0; kk = (int)pp; - kk = min(kk, this->nrar - 1); + kk = MIN(kk, this->nrar - 1); pp = pp - kk; - pp = min(pp, 1.0); + pp = MIN(pp, 1.0); double result = ((arr2(this->phirar3, kk, ind) * pp + arr2(this->phirar2, kk, ind)) * pp + arr2(this->phirar1, kk, ind)) * pp + 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 dba4f40d1..6d390e2e5 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -1,97 +1,98 @@ #include "meam.h" #include +using namespace LAMMPS_NS; // // declaration in pair_meam.h: // // void meam_setup_global(int *, int *, double *, int *, double *, double *, // double *, double *, double *, double *, double *, // double *, double *, double *, double *, double *, // double *, double *, int *); // // call in pair_meam.cpp: // // meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, // alat,esub,asub,t0,t1,t2,t3,rozero,ibar); // // void MEAM::meam_setup_global_(int* nelt, int* lat, double* z, int* ielement, double* atwt, double* alpha, double* b0, double* b1, double* b2, double* b3, double* alat, double* esub, double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero, int* ibar) { int i; double tmplat[maxelt]; this->neltypes = *nelt; for (i = 1; i <= *nelt; i++) { if (arr1v(lat, i) == 0) this->lattce_meam[i][i] = FCC; else if (arr1v(lat, i) == 1) this->lattce_meam[i][i] = BCC; else if (arr1v(lat, i) == 2) this->lattce_meam[i][i] = HCP; else if (arr1v(lat, i) == 3) this->lattce_meam[i][i] = DIM; else if (arr1v(lat, i) == 4) this->lattce_meam[i][i] = DIA; else { // unknown } 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); 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 (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 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; } diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index 7b8b5e2ce..a04e5f818 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -1,238 +1,239 @@ #include "meam.h" +using namespace LAMMPS_NS; // // do a sanity check on index parameters void MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr) { //: idx[0..2] *ierr = 0; if (nidx < num) { *ierr = 2; return; } for (int i = 0; i < num; i++) { if ((idx[i] < 1) || (idx[i] > lim)) { *ierr = 3; return; } } } // // Declaration in pair_meam.h: // // void meam_setup_param(int *, double *, int *, int *, int *); // // in pair_meam.cpp // // meam_setup_param(&which,&value,&nindex,index,&errorflag); // // // // The "which" argument corresponds to the index of the "keyword" array // in pair_meam.cpp: // // 0 = Ec_meam // 1 = alpha_meam // 2 = rho0_meam // 3 = delta_meam // 4 = lattce_meam // 5 = attrac_meam // 6 = repuls_meam // 7 = nn2_meam // 8 = Cmin_meam // 9 = Cmax_meam // 10 = rc_meam // 11 = delr_meam // 12 = augt1 // 13 = gsmooth_factor // 14 = re_meam // 15 = ialloy // 16 = mixture_ref_t // 17 = erose_form // 18 = zbl_meam // 19 = emb_lin_neg // 20 = bkgd_dyn void MEAM::meam_setup_param_(int* which_p, double* value_p, int* nindex_p, int* index /*index(3)*/, int* errorflag) { //: index[0..2] int i1, i2, val; *errorflag = 0; int which = *which_p; double value = *value_p; int nindex = *nindex_p; switch (which) { // 0 = Ec_meam case 0: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->Ec_meam[index[0]][index[1]] = value; break; // 1 = alpha_meam case 1: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->alpha_meam[index[0]][index[1]] = value; break; // 2 = rho0_meam case 2: meam_checkindex(1, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->rho0_meam[index[0]] = value; break; // 3 = delta_meam case 3: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->delta_meam[index[0]][index[1]] = value; break; // 4 = lattce_meam case 4: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; val = (int)value; if (val == 0) this->lattce_meam[index[0]][index[1]] = FCC; else if (val == 1) this->lattce_meam[index[0]][index[1]] = BCC; else if (val == 2) this->lattce_meam[index[0]][index[1]] = HCP; else if (val == 3) this->lattce_meam[index[0]][index[1]] = DIM; else if (val == 4) this->lattce_meam[index[0]][index[1]] = DIA; else if (val == 5) this->lattce_meam[index[0]][index[1]] = B1; else if (val == 6) this->lattce_meam[index[0]][index[1]] = C11; else if (val == 7) this->lattce_meam[index[0]][index[1]] = L12; else if (val == 8) this->lattce_meam[index[0]][index[1]] = B2; break; // 5 = attrac_meam case 5: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->attrac_meam[index[0]][index[1]] = value; break; // 6 = repuls_meam case 6: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->repuls_meam[index[0]][index[1]] = value; break; // 7 = nn2_meam case 7: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - i1 = min(index[0], index[1]); - i2 = max(index[0], index[1]); + i1 = MIN(index[0], index[1]); + i2 = MAX(index[0], index[1]); this->nn2_meam[i1][i2] = (int)value; break; // 8 = Cmin_meam case 8: meam_checkindex(3, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->Cmin_meam[index[0]][index[1]][index[2]] = value; break; // 9 = Cmax_meam case 9: meam_checkindex(3, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->Cmax_meam[index[0]][index[1]][index[2]] = value; break; // 10 = rc_meam case 10: this->rc_meam = value; break; // 11 = delr_meam case 11: this->delr_meam = value; break; // 12 = augt1 case 12: this->augt1 = (int)value; break; // 13 = gsmooth case 13: this->gsmooth_factor = value; break; // 14 = re_meam case 14: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; this->re_meam[index[0]][index[1]] = value; break; // 15 = ialloy case 15: this->ialloy = (int)value; break; // 16 = mixture_ref_t case 16: this->mix_ref_t = (int)value; break; // 17 = erose_form case 17: this->erose_form = (int)value; break; // 18 = zbl_meam case 18: meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - i1 = min(index[0], index[1]); - i2 = max(index[0], index[1]); + i1 = MIN(index[0], index[1]); + i2 = MAX(index[0], index[1]); this->zbl_meam[i1][i2] = (int)value; break; // 19 = emb_lin_neg case 19: this->emb_lin_neg = (int)value; break; // 20 = bkgd_dyn case 20: this->bkgd_dyn = (int)value; break; default: *errorflag = 1; } } diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index a92c9c6b5..94b0df7f7 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -1,944 +1,946 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Greg Wagner (SNL) ------------------------------------------------------------------------- */ #include #include #include #include +#include "meam.h" #include "pair_meamc.h" #include "atom.h" #include "force.h" #include "comm.h" #include "memory.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define MAXLINE 1024 static const int nkeywords = 21; static const char *keywords[] = { "Ec","alpha","rho0","delta","lattce", "attrac","repuls","nn2","Cmin","Cmax","rc","delr", "augt1","gsmooth_factor","re","ialloy", "mixture_ref_t","erose_form","zbl", "emb_lin_neg","bkgd_dyn"}; /* ---------------------------------------------------------------------- */ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; restartinfo = 0; one_coeff = 1; manybody_flag = 1; nmax = 0; rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL; gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL; arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL; maxneigh = 0; allocated = 0; scrfcn = dscrfcn = fcpair = NULL; nelements = 0; elements = NULL; mass = NULL; + meam_inst = new MEAM; // set comm size needed by this Pair comm_forward = 38; comm_reverse = 30; } /* ---------------------------------------------------------------------- free all arrays check if allocated, since class can be destructed when incomplete ------------------------------------------------------------------------- */ PairMEAMC::~PairMEAMC() { - meam_inst.meam_cleanup_(); + delete meam_inst; memory->destroy(rho); memory->destroy(rho0); memory->destroy(rho1); memory->destroy(rho2); memory->destroy(rho3); memory->destroy(frhop); memory->destroy(gamma); memory->destroy(dgamma1); memory->destroy(dgamma2); memory->destroy(dgamma3); memory->destroy(arho2b); memory->destroy(arho1); memory->destroy(arho2); memory->destroy(arho3); memory->destroy(arho3b); memory->destroy(t_ave); memory->destroy(tsq_ave); memory->destroy(scrfcn); memory->destroy(dscrfcn); memory->destroy(fcpair); for (int i = 0; i < nelements; i++) delete [] elements[i]; delete [] elements; delete [] mass; if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); delete [] map; delete [] fmap; } } /* ---------------------------------------------------------------------- */ void PairMEAMC::compute(int eflag, int vflag) { int i,j,ii,n,inum_half,errorflag; int *ilist_half,*numneigh_half,**firstneigh_half; int *numneigh_full,**firstneigh_full; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = eflag_global = vflag_global = eflag_atom = vflag_atom = 0; // grow local arrays if necessary if (atom->nmax > nmax) { memory->destroy(rho); memory->destroy(rho0); memory->destroy(rho1); memory->destroy(rho2); memory->destroy(rho3); memory->destroy(frhop); memory->destroy(gamma); memory->destroy(dgamma1); memory->destroy(dgamma2); memory->destroy(dgamma3); memory->destroy(arho2b); memory->destroy(arho1); memory->destroy(arho2); memory->destroy(arho3); memory->destroy(arho3b); memory->destroy(t_ave); memory->destroy(tsq_ave); nmax = atom->nmax; memory->create(rho,nmax,"pair:rho"); memory->create(rho0,nmax,"pair:rho0"); memory->create(rho1,nmax,"pair:rho1"); memory->create(rho2,nmax,"pair:rho2"); memory->create(rho3,nmax,"pair:rho3"); memory->create(frhop,nmax,"pair:frhop"); memory->create(gamma,nmax,"pair:gamma"); memory->create(dgamma1,nmax,"pair:dgamma1"); memory->create(dgamma2,nmax,"pair:dgamma2"); memory->create(dgamma3,nmax,"pair:dgamma3"); memory->create(arho2b,nmax,"pair:arho2b"); memory->create(arho1,nmax,3,"pair:arho1"); memory->create(arho2,nmax,6,"pair:arho2"); memory->create(arho3,nmax,10,"pair:arho3"); memory->create(arho3b,nmax,3,"pair:arho3b"); memory->create(t_ave,nmax,3,"pair:t_ave"); memory->create(tsq_ave,nmax,3,"pair:tsq_ave"); } // neighbor list info inum_half = listhalf->inum; ilist_half = listhalf->ilist; numneigh_half = listhalf->numneigh; firstneigh_half = listhalf->firstneigh; numneigh_full = listfull->numneigh; firstneigh_full = listfull->firstneigh; // strip neighbor lists of any special bond flags before using with MEAM // necessary before doing neigh_f2c and neigh_c2f conversions each step if (neighbor->ago == 0) { neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full); } // check size of scrfcn based on half neighbor list int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; n = 0; for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]]; if (n > maxneigh) { memory->destroy(scrfcn); memory->destroy(dscrfcn); memory->destroy(fcpair); maxneigh = n; memory->create(scrfcn,maxneigh,"pair:scrfcn"); memory->create(dscrfcn,maxneigh,"pair:dscrfcn"); memory->create(fcpair,maxneigh,"pair:fcpair"); } // zero out local arrays for (i = 0; i < nall; i++) { rho0[i] = 0.0; arho2b[i] = 0.0; arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0; for (j = 0; j < 6; j++) arho2[i][j] = 0.0; for (j = 0; j < 10; j++) arho3[i][j] = 0.0; arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0; t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0; tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0; } double **x = atom->x; double **f = atom->f; int *type = atom->type; int ntype = atom->ntypes; // change neighbor list indices to Fortran indexing neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full); // 3 stages of MEAM calculation // loop over my atoms followed by communication int ifort; int offset = 0; errorflag = 0; for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_inst.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], rho0,&arho1[0][0],&arho2[0][0],arho2b, &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0], &errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); error->one(FLERR,str); } offset += numneigh_half[i]; } comm->reverse_comm_pair(this); - meam_inst.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, dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); error->one(FLERR,str); } comm->forward_comm_pair(this); offset = 0; // vptr is first value in vatom if it will be used by meam_force() // else vatom may not exist, so pass dummy ptr double *vptr; if (vflag_atom) vptr = &vatom[0][0]; else vptr = &cutmax; for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_inst.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], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], &t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); error->one(FLERR,str); } offset += numneigh_half[i]; } // change neighbor list indices back to C indexing neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- */ void PairMEAMC::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); map = new int[n+1]; fmap = new int[n]; } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairMEAMC::settings(int narg, char **arg) { if (narg != 0) error->all(FLERR,"Illegal pair_style command"); } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairMEAMC::coeff(int narg, char **arg) { int i,j,m,n; if (!allocated) allocate(); if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients"); // insure I,J args are * * if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->all(FLERR,"Incorrect args for pair coefficients"); // read MEAM element names between 2 filenames // nelements = # of MEAM elements // elements = list of unique element names if (nelements) { for (i = 0; i < nelements; i++) delete [] elements[i]; delete [] elements; delete [] mass; } nelements = narg - 4 - atom->ntypes; if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients"); elements = new char*[nelements]; mass = new double[nelements]; for (i = 0; i < nelements; i++) { n = strlen(arg[i+3]) + 1; elements[i] = new char[n]; strcpy(elements[i],arg[i+3]); } // read MEAM library and parameter files // pass all parameters to MEAM package // tell MEAM package that setup is done read_files(arg[2],arg[2+nelements+1]); - meam_inst.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 for (i = 4 + nelements; i < narg; i++) { m = i - (4+nelements) + 1; for (j = 0; j < nelements; j++) if (strcmp(arg[i],elements[j]) == 0) break; if (j < nelements) map[m] = j; else if (strcmp(arg[i],"NULL") == 0) map[m] = -1; else error->all(FLERR,"Incorrect args for pair coefficients"); } // clear setflag since coeff() called once with I,J = * * n = atom->ntypes; for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; // set setflag i,j for type pairs where both are mapped to elements // set mass for i,i in atom class int count = 0; for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) if (map[i] >= 0 && map[j] >= 0) { setflag[i][j] = 1; if (i == j) atom->set_mass(FLERR,i,mass[map[i]]); count++; } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairMEAMC::init_style() { if (force->newton_pair == 0) error->all(FLERR,"Pair style MEAM requires newton pair on"); // need full and half neighbor list int irequest_full = neighbor->request(this,instance_me); neighbor->requests[irequest_full]->id = 1; neighbor->requests[irequest_full]->half = 0; neighbor->requests[irequest_full]->full = 1; int irequest_half = neighbor->request(this,instance_me); neighbor->requests[irequest_half]->id = 2; // setup Fortran-style mapping array needed by MEAM package // fmap is indexed from 1:ntypes by Fortran and stores a Fortran index // if type I is not a MEAM atom, fmap stores a 0 for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1; } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use half or full ------------------------------------------------------------------------- */ void PairMEAMC::init_list(int id, NeighList *ptr) { if (id == 1) listfull = ptr; else if (id == 2) listhalf = ptr; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairMEAMC::init_one(int i, int j) { return cutmax; } /* ---------------------------------------------------------------------- */ void PairMEAMC::read_files(char *globalfile, char *userfile) { // open global meamf file on proc 0 FILE *fp; if (comm->me == 0) { fp = force->open_potential(globalfile); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open MEAM potential file %s",globalfile); error->one(FLERR,str); } } // allocate parameter arrays int params_per_line = 19; int *lat = new int[nelements]; int *ielement = new int[nelements]; int *ibar = new int[nelements]; double *z = new double[nelements]; double *atwt = new double[nelements]; double *alpha = new double[nelements]; double *b0 = new double[nelements]; double *b1 = new double[nelements]; double *b2 = new double[nelements]; double *b3 = new double[nelements]; double *alat = new double[nelements]; double *esub = new double[nelements]; double *asub = new double[nelements]; double *t0 = new double[nelements]; double *t1 = new double[nelements]; double *t2 = new double[nelements]; double *t3 = new double[nelements]; double *rozero = new double[nelements]; bool *found = new bool[nelements]; for (int i = 0; i < nelements; i++) found[i] = false; // read each set of params from global MEAM file // one set of params can span multiple lines // store params if element name is in element list // if element name appears multiple times, only store 1st entry int i,n,nwords; char **words = new char*[params_per_line+1]; char line[MAXLINE],*ptr; int eof = 0; int nset = 0; while (1) { if (comm->me == 0) { ptr = fgets(line,MAXLINE,fp); if (ptr == NULL) { eof = 1; fclose(fp); } else n = strlen(line) + 1; } MPI_Bcast(&eof,1,MPI_INT,0,world); if (eof) break; MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); // strip comment, skip line if blank if ((ptr = strchr(line,'#'))) *ptr = '\0'; nwords = atom->count_words(line); if (nwords == 0) continue; // concatenate additional lines until have params_per_line words while (nwords < params_per_line) { n = strlen(line); if (comm->me == 0) { ptr = fgets(&line[n],MAXLINE-n,fp); if (ptr == NULL) { eof = 1; fclose(fp); } else n = strlen(line) + 1; } MPI_Bcast(&eof,1,MPI_INT,0,world); if (eof) break; MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); if ((ptr = strchr(line,'#'))) *ptr = '\0'; nwords = atom->count_words(line); } if (nwords != params_per_line) error->all(FLERR,"Incorrect format in MEAM potential file"); // words = ptrs to all words in line // strip single and double quotes from words nwords = 0; words[nwords++] = strtok(line,"' \t\n\r\f"); while ((words[nwords++] = strtok(NULL,"' \t\n\r\f"))) continue; // skip if element name isn't in element list for (i = 0; i < nelements; i++) if (strcmp(words[0],elements[i]) == 0) break; if (i >= nelements) continue; // skip if element already appeared if (found[i] == true) continue; found[i] = true; // map lat string to an integer if (strcmp(words[1],"fcc") == 0) lat[i] = FCC; else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC; else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP; else if (strcmp(words[1],"dim") == 0) lat[i] = DIM; else if (strcmp(words[1],"dia") == 0) lat[i] = DIA; else error->all(FLERR,"Unrecognized lattice type in MEAM file 1"); // store parameters z[i] = atof(words[2]); ielement[i] = atoi(words[3]); atwt[i] = atof(words[4]); alpha[i] = atof(words[5]); b0[i] = atof(words[6]); b1[i] = atof(words[7]); b2[i] = atof(words[8]); b3[i] = atof(words[9]); alat[i] = atof(words[10]); esub[i] = atof(words[11]); asub[i] = atof(words[12]); t0[i] = atof(words[13]); t1[i] = atof(words[14]); t2[i] = atof(words[15]); t3[i] = atof(words[16]); rozero[i] = atof(words[17]); ibar[i] = atoi(words[18]); nset++; } // error if didn't find all elements in file if (nset != nelements) error->all(FLERR,"Did not find all elements in MEAM library file"); // pass element parameters to MEAM package - meam_inst.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 for (i = 0; i < nelements; i++) mass[i] = atwt[i]; // clean-up memory delete [] words; delete [] lat; delete [] ielement; delete [] ibar; delete [] z; delete [] atwt; delete [] alpha; delete [] b0; delete [] b1; delete [] b2; delete [] b3; delete [] alat; delete [] esub; delete [] asub; delete [] t0; delete [] t1; delete [] t2; delete [] t3; delete [] rozero; delete [] found; // done if user param file is NULL if (strcmp(userfile,"NULL") == 0) return; // open user param file on proc 0 if (comm->me == 0) { fp = force->open_potential(userfile); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open MEAM potential file %s",userfile); error->one(FLERR,str); } } // read settings // pass them one at a time to MEAM package // match strings to list of corresponding ints int which; double value; int nindex,index[3]; int maxparams = 6; char **params = new char*[maxparams]; int nparams; eof = 0; while (1) { if (comm->me == 0) { ptr = fgets(line,MAXLINE,fp); if (ptr == NULL) { eof = 1; fclose(fp); } else n = strlen(line) + 1; } MPI_Bcast(&eof,1,MPI_INT,0,world); if (eof) break; MPI_Bcast(&n,1,MPI_INT,0,world); MPI_Bcast(line,n,MPI_CHAR,0,world); // strip comment, skip line if blank if ((ptr = strchr(line,'#'))) *ptr = '\0'; nparams = atom->count_words(line); if (nparams == 0) continue; // words = ptrs to all words in line nparams = 0; params[nparams++] = strtok(line,"=(), '\t\n\r\f"); while (nparams < maxparams && (params[nparams++] = strtok(NULL,"=(), '\t\n\r\f"))) continue; nparams--; for (which = 0; which < nkeywords; which++) if (strcmp(params[0],keywords[which]) == 0) break; if (which == nkeywords) { char str[128]; sprintf(str,"Keyword %s in MEAM parameter file not recognized", params[0]); error->all(FLERR,str); } nindex = nparams - 2; for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]); // map lattce_meam value to an integer if (which == 4) { if (strcmp(params[nparams-1],"fcc") == 0) value = FCC; else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC; else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP; else if (strcmp(params[nparams-1],"dim") == 0) value = DIM; else if (strcmp(params[nparams-1],"dia") == 0) value = DIA; else if (strcmp(params[nparams-1],"b1") == 0) value = B1; else if (strcmp(params[nparams-1],"c11") == 0) value = C11; else if (strcmp(params[nparams-1],"l12") == 0) value = L12; else if (strcmp(params[nparams-1],"b2") == 0) value = B2; else error->all(FLERR,"Unrecognized lattice type in MEAM file 2"); } else value = atof(params[nparams-1]); // pass single setting to MEAM package int errorflag = 0; - meam_inst.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); error->all(FLERR,str); } } delete [] params; } /* ---------------------------------------------------------------------- */ int PairMEAMC::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { int i,j,k,m; m = 0; for (i = 0; i < n; i++) { j = list[i]; buf[m++] = rho0[j]; buf[m++] = rho1[j]; buf[m++] = rho2[j]; buf[m++] = rho3[j]; buf[m++] = frhop[j]; buf[m++] = gamma[j]; buf[m++] = dgamma1[j]; buf[m++] = dgamma2[j]; buf[m++] = dgamma3[j]; buf[m++] = arho2b[j]; buf[m++] = arho1[j][0]; buf[m++] = arho1[j][1]; buf[m++] = arho1[j][2]; buf[m++] = arho2[j][0]; buf[m++] = arho2[j][1]; buf[m++] = arho2[j][2]; buf[m++] = arho2[j][3]; buf[m++] = arho2[j][4]; buf[m++] = arho2[j][5]; for (k = 0; k < 10; k++) buf[m++] = arho3[j][k]; buf[m++] = arho3b[j][0]; buf[m++] = arho3b[j][1]; buf[m++] = arho3b[j][2]; buf[m++] = t_ave[j][0]; buf[m++] = t_ave[j][1]; buf[m++] = t_ave[j][2]; buf[m++] = tsq_ave[j][0]; buf[m++] = tsq_ave[j][1]; buf[m++] = tsq_ave[j][2]; } return m; } /* ---------------------------------------------------------------------- */ void PairMEAMC::unpack_forward_comm(int n, int first, double *buf) { int i,k,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { rho0[i] = buf[m++]; rho1[i] = buf[m++]; rho2[i] = buf[m++]; rho3[i] = buf[m++]; frhop[i] = buf[m++]; gamma[i] = buf[m++]; dgamma1[i] = buf[m++]; dgamma2[i] = buf[m++]; dgamma3[i] = buf[m++]; arho2b[i] = buf[m++]; arho1[i][0] = buf[m++]; arho1[i][1] = buf[m++]; arho1[i][2] = buf[m++]; arho2[i][0] = buf[m++]; arho2[i][1] = buf[m++]; arho2[i][2] = buf[m++]; arho2[i][3] = buf[m++]; arho2[i][4] = buf[m++]; arho2[i][5] = buf[m++]; for (k = 0; k < 10; k++) arho3[i][k] = buf[m++]; arho3b[i][0] = buf[m++]; arho3b[i][1] = buf[m++]; arho3b[i][2] = buf[m++]; t_ave[i][0] = buf[m++]; t_ave[i][1] = buf[m++]; t_ave[i][2] = buf[m++]; tsq_ave[i][0] = buf[m++]; tsq_ave[i][1] = buf[m++]; tsq_ave[i][2] = buf[m++]; } } /* ---------------------------------------------------------------------- */ int PairMEAMC::pack_reverse_comm(int n, int first, double *buf) { int i,k,m,last; m = 0; last = first + n; for (i = first; i < last; i++) { buf[m++] = rho0[i]; buf[m++] = arho2b[i]; buf[m++] = arho1[i][0]; buf[m++] = arho1[i][1]; buf[m++] = arho1[i][2]; buf[m++] = arho2[i][0]; buf[m++] = arho2[i][1]; buf[m++] = arho2[i][2]; buf[m++] = arho2[i][3]; buf[m++] = arho2[i][4]; buf[m++] = arho2[i][5]; for (k = 0; k < 10; k++) buf[m++] = arho3[i][k]; buf[m++] = arho3b[i][0]; buf[m++] = arho3b[i][1]; buf[m++] = arho3b[i][2]; buf[m++] = t_ave[i][0]; buf[m++] = t_ave[i][1]; buf[m++] = t_ave[i][2]; buf[m++] = tsq_ave[i][0]; buf[m++] = tsq_ave[i][1]; buf[m++] = tsq_ave[i][2]; } return m; } /* ---------------------------------------------------------------------- */ void PairMEAMC::unpack_reverse_comm(int n, int *list, double *buf) { int i,j,k,m; m = 0; for (i = 0; i < n; i++) { j = list[i]; rho0[j] += buf[m++]; arho2b[j] += buf[m++]; arho1[j][0] += buf[m++]; arho1[j][1] += buf[m++]; arho1[j][2] += buf[m++]; arho2[j][0] += buf[m++]; arho2[j][1] += buf[m++]; arho2[j][2] += buf[m++]; arho2[j][3] += buf[m++]; arho2[j][4] += buf[m++]; arho2[j][5] += buf[m++]; for (k = 0; k < 10; k++) arho3[j][k] += buf[m++]; arho3b[j][0] += buf[m++]; arho3b[j][1] += buf[m++]; arho3b[j][2] += buf[m++]; t_ave[j][0] += buf[m++]; t_ave[j][1] += buf[m++]; t_ave[j][2] += buf[m++]; tsq_ave[j][0] += buf[m++]; tsq_ave[j][1] += buf[m++]; tsq_ave[j][2] += buf[m++]; } } /* ---------------------------------------------------------------------- memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double PairMEAMC::memory_usage() { double bytes = 11 * nmax * sizeof(double); bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double); bytes += 3 * maxneigh * sizeof(double); return bytes; } /* ---------------------------------------------------------------------- strip special bond flags from neighbor list entries are not used with MEAM need to do here so Fortran lib doesn't see them done once per reneighbor so that neigh_f2c and neigh_c2f don't see them ------------------------------------------------------------------------- */ void PairMEAMC::neigh_strip(int inum, int *ilist, int *numneigh, int **firstneigh) { int i,j,ii,jnum; int *jlist; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK; } } /* ---------------------------------------------------------------------- toggle neighbor list indices between zero- and one-based values needed for access by MEAM Fortran library ------------------------------------------------------------------------- */ void PairMEAMC::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh) { int i,j,ii,jnum; int *jlist; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; for (j = 0; j < jnum; j++) jlist[j]--; } } void PairMEAMC::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh) { int i,j,ii,jnum; int *jlist; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jlist = firstneigh[i]; jnum = numneigh[i]; for (j = 0; j < jnum; j++) jlist[j]++; } } diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h index c9aee3936..e566f40a0 100644 --- a/src/USER-MEAMC/pair_meamc.h +++ b/src/USER-MEAMC/pair_meamc.h @@ -1,123 +1,123 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS PairStyle(meam/c,PairMEAMC) #else #ifndef LMP_PAIR_MEAMC_H #define LMP_PAIR_MEAMC_H #include "pair.h" -#include "meam.h" namespace LAMMPS_NS { +class MEAM; class PairMEAMC : public Pair { public: PairMEAMC(class LAMMPS *); ~PairMEAMC(); void compute(int, int); void settings(int, char **); void coeff(int, char **); void init_style(); void init_list(int, class NeighList *); double init_one(int, int); int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); int pack_reverse_comm(int, int, double *); void unpack_reverse_comm(int, int *, double *); double memory_usage(); private: - MEAM meam_inst; + class MEAM *meam_inst; double cutmax; // max cutoff for all elements int nelements; // # of unique elements char **elements; // names of unique elements double *mass; // mass of each element int *map; // mapping from atom types to elements int *fmap; // Fortran version of map array for MEAM lib int maxneigh; double *scrfcn,*dscrfcn,*fcpair; int nmax; double *rho,*rho0,*rho1,*rho2,*rho3,*frhop; double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b; double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave; void allocate(); void read_files(char *, char *); void neigh_strip(int, int *, int *, int **); void neigh_f2c(int, int *, int *, int **); void neigh_c2f(int, int *, int *, int **); }; } #endif #endif /* ERROR/WARNING messages: E: MEAM library error %d A call to the MEAM Fortran library returned an error. E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Incorrect args for pair coefficients Self-explanatory. Check the input script or data file. E: Pair style MEAM requires newton pair on See the newton command. This is a restriction to use the MEAM potential. E: Cannot open MEAM potential file %s The specified MEAM potential file cannot be opened. Check that the path and name are correct. E: Incorrect format in MEAM potential file Incorrect number of words per line in the potential file. E: Unrecognized lattice type in MEAM file 1 The lattice type in an entry of the MEAM library file is not valid. E: Did not find all elements in MEAM library file The requested elements were not found in the MEAM file. E: Keyword %s in MEAM parameter file not recognized Self-explanatory. E: Unrecognized lattice type in MEAM file 2 The lattice type in an entry of the MEAM parameter file is not valid. */