diff --git a/doc/src/pair_meamc.rst b/doc/src/pair_meamc.rst index daa10906c7..ac858ce24b 100644 --- a/doc/src/pair_meamc.rst +++ b/doc/src/pair_meamc.rst @@ -233,15 +233,20 @@ where Cmin(I,J,K) = Cmin screening parameter when I-J pair is screened by K (I<=J); default = 2.0 lattce(I,J) = lattice structure of I-J reference structure: - dia = diamond (interlaced fcc for alloy) fcc = face centered cubic bcc = body centered cubic - dim = dimer - b1 = rock salt (NaCl structure) hcp = hexagonal close-packed + dim = dimer + dia = diamond (interlaced fcc for alloy) + dia3= diamond structure with primary 1NN and secondary 3NN interation + b1 = rock salt (NaCl structure) c11 = MoSi2 structure l12 = Cu3Au structure (lower case L, followed by 12) b2 = CsCl structure (interpenetrating simple cubic) + ch4 = methane-like structure, only for binary system + lin = linear structure (180 degree angle) + zig = zigzag structure with a uniform angle + tri = H2O-like structure that has an angle nn2(I,J) = turn on second-nearest neighbor MEAM formulation for I-J pair (see for example :ref:`(Lee) `). 0 = second-nearest neighbor formulation off @@ -254,6 +259,8 @@ where zbl(I,J) = blend the MEAM I-J pair potential with the ZBL potential for small atom separations :ref:`(ZBL) ` default = 1 + theta(I,J) = angle between three atoms in line, zigzag, and trimer reference structures in degrees + default = 180 gsmooth_factor = factor determining the length of the G-function smoothing region; only significant for ibar=0 or ibar=4. 99.0 = short smoothing region, sharp step diff --git a/doc/txt/pair_meamc.txt b/doc/txt/pair_meamc.txt index 7c42e9d2f2..2831600e08 100644 --- a/doc/txt/pair_meamc.txt +++ b/doc/txt/pair_meamc.txt @@ -139,6 +139,8 @@ potential parameters: lat = lattice structure of reference configuration z = number of nearest neighbors in the reference structure + This field is only read for compatibility, the correct + value is inferred from the lattice structure ielement = atomic number atwt = atomic weight alat = lattice constant of reference structure diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 4e5b298ffb..719e1af6f0 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -2,13 +2,15 @@ #define LMP_MEAM_H #include +#include +#include "math_const.h" #define maxelt 5 namespace LAMMPS_NS { class Memory; -typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t; +typedef enum { FCC, BCC, HCP, DIM, DIA, DIA3, B1, C11, L12, B2, CH4, LIN, ZIG, TRI } lattice_t; class MEAM { @@ -26,9 +28,7 @@ private: // Ec_meam = cohesive energy // re_meam = nearest-neighbor distance - // Omega_meam = atomic volume // B_meam = bulk modulus - // Z_meam = number of first neighbors for reference structure // ielt_meam = atomic number of element // A_meam = adjustable parameter // alpha_meam = sqrt(9*Omega*B/Ec) @@ -64,8 +64,11 @@ private: // nr,dr = pair function discretization parameters // nrar,rdrar = spline coeff array parameters + // theta = angle between three atoms in line, zigzag, and trimer reference structures + // stheta_meam = sin(theta/2) in radian used in line, zigzag, and trimer reference structures + // ctheta_meam = cos(theta/2) in radian used in line, zigzag, and trimer reference structures + double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt]; - double Omega_meam[maxelt], Z_meam[maxelt]; double A_meam[maxelt], alpha_meam[maxelt][maxelt], rho0_meam[maxelt]; double delta_meam[maxelt][maxelt]; double beta0_meam[maxelt], beta1_meam[maxelt]; @@ -108,6 +111,10 @@ public: int maxneigh; double *scrfcn, *dscrfcn, *fcpair; + //angle for trimer, zigzag, line reference structures + double stheta_meam[maxelt][maxelt]; + double ctheta_meam[maxelt][maxelt]; + protected: // meam_funcs.cpp @@ -191,9 +198,12 @@ protected: double embedding(const double A, const double Ec, const double rhobar, double& dF) const; static double erose(const double r, const double re, const double alpha, const double Ec, const double repuls, const double attrac, const int form); - static void get_shpfcn(const lattice_t latt, double (&s)[3]); - static int get_Zij(const lattice_t latt); - static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, double &a, double &S); + static void get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3]); + + static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, + const double sthe, double &a, double &S); + static int get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S); + protected: void meam_checkindex(int, int, int, int*, int*); void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh, @@ -204,6 +214,7 @@ protected: void alloyparams(); void compute_pair_meam(); double phi_meam(double, int, int); + const double phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat); void compute_reference_density(); void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double, double, double, double, int, int, lattice_t); @@ -212,7 +223,41 @@ protected: void interpolate_meam(int); public: - void meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha, + //----------------------------------------------------------------------------- + // convert lattice spec to lattice_t + // only use single-element lattices if single=true + // return false on failure + // return true and set lat on success + static bool str_to_lat(const char* str, bool single, lattice_t& lat) + { + if (strcmp(str,"fcc") == 0) lat = FCC; + else if (strcmp(str,"bcc") == 0) lat = BCC; + else if (strcmp(str,"hcp") == 0) lat = HCP; + else if (strcmp(str,"dim") == 0) lat = DIM; + else if (strcmp(str,"dia") == 0) lat = DIA; + else if (strcmp(str,"dia3") == 0) lat = DIA3; + else if (strcmp(str,"lin") == 0) lat = LIN; + else if (strcmp(str,"zig") == 0) lat = ZIG; + else if (strcmp(str,"tri") == 0) lat = TRI; + else { + if (single) + return false; + + if (strcmp(str,"b1") == 0) lat = B1; + else if (strcmp(str,"c11") == 0) lat = C11; + else if (strcmp(str,"l12") == 0) lat = L12; + else if (strcmp(str,"b2") == 0) lat = B2; + else if (strcmp(str,"ch4") == 0) lat = CH4; + else if (strcmp(str,"lin") == 0) lat =LIN; + else if (strcmp(str,"zig") == 0) lat = ZIG; + else if (strcmp(str,"tri") == 0) lat = TRI; + else return false; + } + return true; + } + + static int get_Zij(const lattice_t latt); + void meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* atwt, double* alpha, double* b0, double* b1, double* b2, double* b3, double* alat, double* esub, double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero, int* ibar); @@ -222,9 +267,9 @@ public: void meam_dens_init(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh, int numneigh_full, int* firstneigh_full, int fnoffset); void meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl, - double* eatom, int ntype, int* type, int* fmap, int& errorflag); + double* eatom, int ntype, int* type, int* fmap, double** scale, int& errorflag); void meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl, - double* eatom, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh, + double* eatom, int ntype, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh, int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom); }; @@ -234,6 +279,10 @@ static inline bool iszero(const double f) { return fabs(f) < 1e-20; } +static inline bool isone(const double f) { + return fabs(f-1.0) < 1e-20; +} + // Helper functions static inline double fdiv_zero(const double n, const double d) { diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index 867106df88..fe9e74ca7c 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -4,18 +4,20 @@ using namespace LAMMPS_NS; void MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl, - double* eatom, int /*ntype*/, int* type, int* fmap, int& errorflag) + double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, int& errorflag) { int i, elti; int m; double rhob, G, dG, Gbar, dGbar, gam, shp[3], Z; double denom, rho_bkgd, Fl; + double scaleii; // Complete the calculation of density for (i = 0; i < nlocal; i++) { elti = fmap[type[i]]; if (elti >= 0) { + scaleii = scale[type[i]][type[i]]; rho1[i] = 0.0; rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i]; rho3[i] = 0.0; @@ -52,12 +54,14 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ gamma[i] = gamma[i] / (rho0[i] * rho0[i]); } - Z = this->Z_meam[elti]; + Z = get_Zij(this->lattce_meam[elti][elti]); G = G_gam(gamma[i], this->ibar_meam[elti], errorflag); if (errorflag != 0) return; - get_shpfcn(this->lattce_meam[elti][elti], shp); + + get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shp); + if (this->ibar_meam[elti] <= 0) { Gbar = 1.0; dGbar = 0.0; @@ -113,6 +117,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]); if (eflag_either != 0) { + Fl *= scaleii; if (eflag_global != 0) { *eng_vdwl = *eng_vdwl + Fl; } diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index 155941422c..39e3ff9467 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -197,7 +197,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** // Now compute derivatives dscrfcn[jn] = 0.0; sfcij = sij * fcij; - if (!iszero(sfcij) && !iszero(sfcij - 1.0)) { + if (!iszero(sfcij) && !isone(sfcij)) { for (kn = 0; kn < numneigh_full; kn++) { k = firstneigh_full[kn]; if (k == j) continue; diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index 09aad90111..2b6832e155 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -8,7 +8,7 @@ using namespace LAMMPS_NS; void MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl, - double* eatom, int /*ntype*/, int* type, int* fmap, double** x, int numneigh, int* firstneigh, + double* eatom, int /*ntype*/, int* type, int* fmap, double** scale, double** x, int numneigh, int* firstneigh, int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom) { int j, jn, k, kn, kk, m, n, p, q; @@ -17,7 +17,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int double v[6], fi[3], fj[3]; double third, sixth; double pp, dUdrij, dUdsij, dUdrijm[3], force, forcem; - double r, recip, phi, phip; + double recip, phi, phip; double sij; double a1, a1i, a1j, a2, a2i, a2j; double a3i, a3j; @@ -42,6 +42,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3; double dsij1, dsij2, force1, force2; double t1i, t2i, t3i, t1j, t2j, t3j; + double scaleij; third = 1.0 / 3.0; sixth = 1.0 / 6.0; @@ -59,6 +60,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int for (jn = 0; jn < numneigh; jn++) { j = firstneigh[jn]; eltj = fmap[type[j]]; + scaleij = scale[type[i]][type[j]]; if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) { @@ -69,7 +71,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2]; if (rij2 < this->cutforcesq) { rij = sqrt(rij2); - r = rij; + recip = 1.0 / rij; // Compute phi and phip ind = this->eltind[elti][eltj]; @@ -78,17 +80,16 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int kk = std::min(kk, this->nrar - 2); pp = pp - kk; pp = std::min(pp, 1.0); - phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + - this->phirar[ind][kk]; + phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + this->phirar[ind][kk]; phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk]; - recip = 1.0 / r; if (eflag_either != 0) { + double phi_sc = phi * scaleij; if (eflag_global != 0) - *eng_vdwl = *eng_vdwl + phi * sij; + *eng_vdwl = *eng_vdwl + phi_sc * sij; if (eflag_atom != 0) { - eatom[i] = eatom[i] + 0.5 * phi * sij; - eatom[j] = eatom[j] + 0.5 * phi * sij; + eatom[i] = eatom[i] + 0.5 * phi_sc * sij; + eatom[j] = eatom[j] + 0.5 * phi_sc * sij; } } @@ -287,8 +288,9 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int } // Compute derivatives of total density wrt rij, sij and rij(3) - get_shpfcn(this->lattce_meam[elti][elti], shpi); - get_shpfcn(this->lattce_meam[eltj][eltj], shpj); + get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpi); + get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpj); + drhodr1 = dgamma1[i] * drho0dr1 + dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 + dt3dr1 * rho3[i] + t3i * drho3dr1) - @@ -379,6 +381,13 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int for (m = 0; m < 3; m++) { dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m]; } + if (!isone(scaleij)) { + dUdrij *= scaleij; + dUdsij *= scaleij; + dUdrijm[0] *= scaleij; + dUdrijm[1] *= scaleij; + dUdrijm[2] *= scaleij; + } // Add the part of the force due to dUdrij and dUdsij @@ -410,7 +419,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int // Now compute forces on other atoms k due to change in sij - if (iszero(sij) || iszero(sij - 1.0)) continue; //: cont jn loop + if (iszero(sij) || isone(sij)) continue; //: cont jn loop double dxik(0), dyik(0), dzik(0); double dxjk(0), dyjk(0), dzjk(0); @@ -429,7 +438,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int dsij1 = 0.0; dsij2 = 0.0; - if (!iszero(sij) && !iszero(sij - 1.0)) { + if (!iszero(sij) && !isone(sij)) { const double rbound = rij2 * this->ebound_meam[elti][eltj]; delc = Cmax - Cmin; dxjk = x[k][0] - x[j][0]; diff --git a/src/USER-MEAMC/meam_funcs.cpp b/src/USER-MEAMC/meam_funcs.cpp index 312caff686..da984cfeb7 100644 --- a/src/USER-MEAMC/meam_funcs.cpp +++ b/src/USER-MEAMC/meam_funcs.cpp @@ -198,7 +198,7 @@ MEAM::erose(const double r, const double re, const double alpha, const double Ec // Shape factors for various configurations // void -MEAM::get_shpfcn(const lattice_t latt, double (&s)[3]) +MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3]) { switch (latt) { case FCC: @@ -214,7 +214,9 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3]) s[1] = 0.0; s[2] = 1.0 / 3.0; break; + case CH4: // CH4 actually needs shape factor for diamond for C, dimer for H case DIA: + case DIA3: s[0] = 0.0; s[1] = 0.0; s[2] = 32.0 / 9.0; @@ -222,8 +224,20 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3]) case DIM: s[0] = 1.0; s[1] = 2.0 / 3.0; - // s(3) = 1.d0 - s[2] = 0.40; + // s(3) = 1.d0 // this should be 0.4 unless (1-legendre) is multiplied in the density calc. + s[2] = 0.40; // this is (1-legendre) where legendre = 0.6 in dynamo is accounted. + break; + case LIN: //linear, theta being 180 + s[0] = 0.0; + s[1] = 8.0 / 3.0; // 4*(co**4 + si**4 - 1.0/3.0) in zig become 4*(1-1/3) + s[2] = 0.0; + break; + case ZIG: //zig-zag + case TRI: //trimer e.g. H2O + s[0] = 4.0*pow(cthe,2); + s[1] = 4.0*(pow(cthe,4) + pow(sthe,4) - 1.0/3.0); + s[2] = 4.0*(pow(cthe,2) * (3*pow(sthe,4) + pow(cthe,4))); + s[2] = s[2] - 0.6*s[0]; //legend in dyn, 0.6 is default value. break; default: s[0] = 0.0; @@ -232,7 +246,7 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3]) } //----------------------------------------------------------------------------- -// Number of neighbors for the reference structure +// Number of first neighbors for reference structure // int MEAM::get_Zij(const lattice_t latt) @@ -244,18 +258,25 @@ MEAM::get_Zij(const lattice_t latt) return 8; case HCP: return 12; - case B1: - return 6; case DIA: + case DIA3: return 4; case DIM: return 1; + case B1: + return 6; case C11: return 10; case L12: return 12; case B2: return 8; + case CH4: // DYNAMO currenly implemented this way while it needs two Z values, 4 and 1 + return 4; + case LIN: + case ZIG: + case TRI: + return 2; // call error('Lattice not defined in get_Zij.') } return 0; @@ -263,11 +284,13 @@ MEAM::get_Zij(const lattice_t latt) //----------------------------------------------------------------------------- // Number of second neighbors for the reference structure -// a = distance ratio R1/R2 -// S = second neighbor screening function +// a = distance ratio R1/R2 (a2nn in dynamo) +// numscr = number of atoms that screen the 2NN bond +// S = second neighbor screening function (xfac, a part of b2nn in dynamo) // int -MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, double& a, double& S) +MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, + const double stheta, double& a, double& S) { double C, x, sijk; @@ -299,7 +322,7 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl numscr = 2; break; - case DIA: + case DIA: // 2NN Zij2 = 12; a = sqrt(8.0 / 3.0); numscr = 1; @@ -308,12 +331,35 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl } break; + case DIA3: // 3NN + Zij2 = 12; + a = sqrt(11.0 / 3.0); + numscr = 4; + if (cmin < 0.500001) { + // call error('can not do 2NN MEAM for dia') + } + break; + + case CH4: //does not have 2nn structure so it returns 0 + case LIN: //line case DIM: // this really shouldn't be allowed; make sure screening is zero a = 1.0; S = 0.0; return 0; + case TRI: //TRI + Zij2 = 1; + a = 2.0*stheta; + numscr=2; + break; + + case ZIG: //zig-zag + Zij2 = 2; + a = 2.0*stheta; + numscr=1; + break; + case L12: Zij2 = 6; a = sqrt(2.0); @@ -335,10 +381,38 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl } // Compute screening for each first neighbor - C = 4.0 / (a * a) - 1.0; + if (latt==DIA3){ // special case for 3NN diamond structure + C = 1.0; + } else { + C = 4.0 / (a * a) - 1.0; + } x = (C - cmin) / (cmax - cmin); sijk = fcut(x); // There are numscr first neighbors screening the second neighbors S = MathSpecial::powint(sijk, numscr); return Zij2; } + +int +MEAM::get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S) +{ + + double x, sijk, C; + int numscr = 0, Zij2 = 0; + switch (latt) { + case ZIG: //zig-zag for b11s and b22s term + case TRI: //trimer for b11s + Zij2 = 2; + numscr=1; + break; + default: + // unknown lattic flag in get Zij + // call error('Lattice not defined in get_Zij.') + break; + } + C = 1.0; + x = (C - cmin) / (cmax - cmin); + sijk = fcut(x); + S = MathSpecial::powint(sijk, numscr); + return Zij2; +} \ No newline at end of file diff --git a/src/USER-MEAMC/meam_impl.cpp b/src/USER-MEAMC/meam_impl.cpp index 74e8af1cde..0a0b95e14a 100644 --- a/src/USER-MEAMC/meam_impl.cpp +++ b/src/USER-MEAMC/meam_impl.cpp @@ -38,7 +38,7 @@ MEAM::MEAM(Memory* mem) neltypes = 0; for (int i = 0; i < maxelt; i++) { - Omega_meam[i] = Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] = + A_meam[i] = rho0_meam[i] = beta0_meam[i] = beta1_meam[i]= beta2_meam[i] = beta3_meam[i] = t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] = rho_ref_meam[i] = ibar_meam[i] = ielt_meam[i] = 0.0; diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp index 37bfce5873..99b330a83a 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -103,6 +103,9 @@ MEAM::alloyparams(void) this->alpha_meam[i][j] = this->alpha_meam[j][i]; this->lattce_meam[i][j] = this->lattce_meam[j][i]; this->nn2_meam[i][j] = this->nn2_meam[j][i]; + // theta for lin,tri,zig references + this->stheta_meam[i][j] = this->stheta_meam[j][i]; + this->ctheta_meam[i][j] = this->ctheta_meam[j][i]; // If i i) { @@ -161,10 +164,10 @@ void MEAM::compute_pair_meam(void) { - double r /*ununsed:, temp*/; + double r, b2nn, phi_val; int j, a, b, nv2; double astar, frac, phizbl; - int n, nmax, Z1, Z2; + int n, Z1, Z2; double arat, rarat, scrn, scrn2; double phiaa, phibb /*unused:,phitmp*/; double C, s111, s112, s221, S11, S22; @@ -215,7 +218,7 @@ MEAM::compute_pair_meam(void) if (this->nn2_meam[a][b] == 1) { Z1 = get_Zij(this->lattce_meam[a][b]); Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], - this->Cmax_meam[a][a][b], arat, scrn); + this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn); // The B1, B2, and L12 cases with NN2 have a trick to them; we need to // compute the contributions from second nearest neighbors, like a-a @@ -229,35 +232,26 @@ MEAM::compute_pair_meam(void) phiaa = phi_meam(rarat, a, a); Z1 = get_Zij(this->lattce_meam[a][a]); Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a], - this->Cmax_meam[a][a][a], arat, scrn); - 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); - } - } + this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn); + phiaa+= phi_meam_series(scrn, Z1, Z2, a, a, rarat, arat); // phi_bb phibb = phi_meam(rarat, b, b); Z1 = get_Zij(this->lattce_meam[b][b]); Z2 = get_Zij2(this->lattce_meam[b][b], this->Cmin_meam[b][b][b], - this->Cmax_meam[b][b][b], arat, scrn); - 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); - } - } + this->Cmax_meam[b][b][b], this->stheta_meam[b][b], arat, scrn); + phibb+= phi_meam_series(scrn, Z1, Z2, b, b, rarat, arat); if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 || this->lattce_meam[a][b] == DIA) { // Add contributions to the B1 or B2 potential Z1 = get_Zij(this->lattce_meam[a][b]); Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], - this->Cmax_meam[a][a][b], arat, scrn); + this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn); this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa; Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[b][b][a], - this->Cmax_meam[b][b][a], arat, scrn2); + this->Cmax_meam[b][b][a], this->stheta_meam[a][b], arat, scrn2); + this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb; } else if (this->lattce_meam[a][b] == L12) { @@ -279,11 +273,7 @@ MEAM::compute_pair_meam(void) } } else { - nmax = 10; - for (n = 1; n <= nmax; n++) { - this->phir[nv2][j] = - this->phir[nv2][j] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b); - } + this->phir[nv2][j]+= phi_meam_series(scrn, Z1, Z2, a, b, r, arat); } } @@ -328,11 +318,12 @@ MEAM::phi_meam(double r, int a, int b) double rho02, rho12, rho22, rho32; double scalfac, phiaa, phibb; double Eu; - double arat, scrn /*unused:,scrn2*/; + double arat, scrn, scrn2; int Z12, errorflag; - int n, nmax, Z1nn, Z2nn; + int n, Z1nn, Z2nn; lattice_t latta /*unused:,lattb*/; double rho_bkgd1, rho_bkgd2; + double b11s, b22s; double phi_m = 0.0; @@ -341,6 +332,8 @@ MEAM::phi_meam(double r, int a, int b) // get number of neighbors in the reference structure // Nref[i][j] = # of i's neighbors of type j + Z1 = get_Zij(this->lattce_meam[a][a]); + Z2 = get_Zij(this->lattce_meam[b][b]); Z12 = get_Zij(this->lattce_meam[a][b]); get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32); @@ -401,19 +394,17 @@ MEAM::phi_meam(double r, int a, int b) // If using mixing rule for t, apply to reference structure; else // use precomputed values if (this->mix_ref_t == 1) { - Z1 = this->Z_meam[a]; - Z2 = this->Z_meam[b]; if (this->ibar_meam[a] <= 0) G1 = 1.0; else { - get_shpfcn(this->lattce_meam[a][a], s1); + get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], s1); Gam1 = (s1[0] * t11av + s1[1] * t21av + s1[2] * t31av) / (Z1 * Z1); G1 = G_gam(Gam1, this->ibar_meam[a], errorflag); } if (this->ibar_meam[b] <= 0) G2 = 1.0; else { - get_shpfcn(this->lattce_meam[b][b], s2); + get_shpfcn(this->lattce_meam[b][b], this->stheta_meam[b][b], this->ctheta_meam[b][b], s2); Gam2 = (s2[0] * t12av + s2[1] * t22av + s2[2] * t32av) / (Z2 * Z2); G2 = G_gam(Gam2, this->ibar_meam[b], errorflag); } @@ -439,8 +430,8 @@ MEAM::phi_meam(double r, int a, int b) rho_bkgd2 = rho0_2; } else { if (this->bkgd_dyn == 1) { - rho_bkgd1 = this->rho0_meam[a] * this->Z_meam[a]; - rho_bkgd2 = this->rho0_meam[b] * this->Z_meam[b]; + rho_bkgd1 = this->rho0_meam[a] * Z1; + rho_bkgd2 = this->rho0_meam[b] * Z2; } else { rho_bkgd1 = this->rho_ref_meam[a]; rho_bkgd2 = this->rho_ref_meam[b]; @@ -475,20 +466,44 @@ MEAM::phi_meam(double r, int a, int b) // account for second neighbor a-a potential here... Z1nn = get_Zij(this->lattce_meam[a][a]); Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a], - this->Cmax_meam[a][a][a], arat, scrn); - 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); - } - } + this->Cmax_meam[a][a][a], this->stheta_meam[a][b], arat, scrn); + + + phiaa += phi_meam_series(scrn, Z1nn, Z2nn, a, a, r, arat); phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa; + } else if (this->lattce_meam[a][b] == CH4) { + phi_m = (5 * Eu - F1 - 4*F2)/4; + + } else if (this->lattce_meam[a][b] == ZIG){ + if (a==b){ + phi_m = (2 * Eu - F1 - F2) / Z12; + } else{ + Z1 = get_Zij(this->lattce_meam[a][b]); + Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn); + b11s = -Z2/Z1*scrn; + Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], scrn2); + b22s = -Z2/Z1*scrn2; + + phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a); + phibb = phi_meam(2.0*this->stheta_meam[a][b]*r, b, b); + phi_m = (2.0*Eu - F1 - F2 + phiaa*b11s + phibb*b22s) / Z12; + } + + } else if (this->lattce_meam[a][b] == TRI) { + if (a==b){ + phi_m = (3.0*Eu - 2.0*F1 - F2) / Z12; + } else { + Z1 = get_Zij(this->lattce_meam[a][b]); + Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn); + b11s = -Z2/Z1*scrn; + phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a); + phi_m = (3.0*Eu - 2.0*F1 - F2 + phiaa*b11s) / Z12; + } + } else { - // // potential is computed from Rose function and embedding energy phi_m = (2 * Eu - F1 - F2) / Z12; - // } // if r = 0, just return 0 @@ -499,6 +514,31 @@ MEAM::phi_meam(double r, int a, int b) return phi_m; } +//----------------------------------------------------------------------c +// Compute 2NN series terms for phi +// To avoid nan values of phir due to rapid decrease of b2nn^n or/and +// argument of phi_meam, i.e. r*arat^n, in some cases (3NN dia with low Cmin value) +// +const double +MEAM::phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat) +{ + double phi_sum = 0.0; + double b2nn, phi_val; + if (scrn > 0.0) { + b2nn = -Z2*scrn/Z1; + for (int n = 1; n <= 10; n++) { + phi_val = MathSpecial::powint(b2nn,n) * phi_meam(r * MathSpecial::powint(arat, n), a, b); + if (iszero(phi_val)) { + // once either term becomes zero at some point, all folliwng will also be zero + // necessary to avoid numerical error (nan or infty) due to exponential decay in phi_meam + break; + } + phi_sum += phi_val; + } + } + return phi_sum; +} + //----------------------------------------------------------------------c // Compute background density for reference structure of each element void @@ -510,11 +550,11 @@ MEAM::compute_reference_density(void) // loop over element types for (a = 0; a < this->neltypes; a++) { - Z = (int)this->Z_meam[a]; + Z = get_Zij(this->lattce_meam[a][a]); if (this->ibar_meam[a] <= 0) Gbar = 1.0; else { - get_shpfcn(this->lattce_meam[a][a], shp); + get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], shp); gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] + this->t3_meam[a] * shp[2]) / (Z * Z); Gbar = G_gam(gam, this->ibar_meam[a], errorflag); } @@ -529,7 +569,7 @@ MEAM::compute_reference_density(void) // screening) if (this->nn2_meam[a][a] == 1) { Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a], - this->Cmax_meam[a][a][a], arat, scrn); + this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn); rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1)); rho0 = rho0 + Z2 * rho0_2nn * scrn; } @@ -555,30 +595,43 @@ MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, dou *t12av = t12; *t22av = t22; *t32av = t32; - } else if (latt == FCC || latt == BCC || latt == DIA || latt == HCP || latt == B1 || latt == DIM || latt == B2) { - // all neighbors are of the opposite type - *t11av = t12; - *t21av = t22; - *t31av = t32; - *t12av = t11; - *t22av = t21; - *t32av = t31; - } else { - a1 = r / this->re_meam[a][a] - 1.0; - a2 = r / this->re_meam[b][b] - 1.0; - rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1); - rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2); - if (latt == L12) { - rho01 = 8 * rhoa01 + 4 * rhoa02; - *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01; + } else switch (latt) { + case FCC: + case BCC: + case DIA: + case DIA3: + case HCP: + case B1: + case DIM: + case B2: + case CH4: + case LIN: + case ZIG: + case TRI: + // all neighbors are of the opposite type + *t11av = t12; + *t21av = t22; + *t31av = t32; *t12av = t11; - *t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01; *t22av = t21; - *t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01; *t32av = t31; - } else { - // call error('Lattice not defined in get_tavref.') - } + break; + default: + a1 = r / this->re_meam[a][a] - 1.0; + a2 = r / this->re_meam[b][b] - 1.0; + rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1); + rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2); + if (latt == L12) { + rho01 = 8 * rhoa01 + 4 * rhoa02; + *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01; + *t12av = t11; + *t21av = (8 * t21 * rhoa01 + 4 * t22 * rhoa02) / rho01; + *t22av = t21; + *t31av = (8 * t31 * rhoa01 + 4 * t32 * rhoa02) / rho01; + *t32av = t31; + } else { + // call error('Lattice not defined in get_tavref.') + } } } @@ -600,7 +653,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* double a1, a2; double s[3]; lattice_t lat; - int Zij2nn; + int Zij,Zij2nn; double rhoa01nn, rhoa02nn; double rhoa01, rhoa11, rhoa21, rhoa31; double rhoa02, rhoa12, rhoa22, rhoa32; @@ -621,6 +674,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* lat = this->lattce_meam[a][b]; + Zij = get_Zij(lat); + *rho11 = 0.0; *rho21 = 0.0; *rho31 = 0.0; @@ -628,64 +683,133 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* *rho22 = 0.0; *rho32 = 0.0; - if (lat == FCC) { - *rho01 = 12.0 * rhoa02; - *rho02 = 12.0 * rhoa01; - } else if (lat == BCC) { - *rho01 = 8.0 * rhoa02; - *rho02 = 8.0 * rhoa01; - } else if (lat == B1) { - *rho01 = 6.0 * rhoa02; - *rho02 = 6.0 * rhoa01; - } else if (lat == DIA) { - *rho01 = 4.0 * rhoa02; - *rho02 = 4.0 * rhoa01; - *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32; - *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31; - } else if (lat == HCP) { - *rho01 = 12 * rhoa02; - *rho02 = 12 * rhoa01; - *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32; - *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31; - } else if (lat == DIM) { - get_shpfcn(DIM, s); - *rho01 = rhoa02; - *rho02 = rhoa01; - *rho11 = s[0] * rhoa12 * rhoa12; - *rho12 = s[0] * rhoa11 * rhoa11; - *rho21 = s[1] * rhoa22 * rhoa22; - *rho22 = s[1] * rhoa21 * rhoa21; - *rho31 = s[2] * rhoa32 * rhoa32; - *rho32 = s[2] * rhoa31 * rhoa31; - } else if (lat == C11) { - *rho01 = rhoa01; - *rho02 = rhoa02; - *rho11 = rhoa11; - *rho12 = rhoa12; - *rho21 = rhoa21; - *rho22 = rhoa22; - *rho31 = rhoa31; - *rho32 = rhoa32; - } else if (lat == L12) { - *rho01 = 8 * rhoa01 + 4 * rhoa02; - *rho02 = 12 * rhoa01; - if (this->ialloy == 1) { - *rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]); - denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]); - if (denom > 0.) - *rho21 = *rho21 / denom * *rho01; - } else - *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22); - } else if (lat == B2) { - *rho01 = 8.0 * rhoa02; - *rho02 = 8.0 * rhoa01; - } else { + switch (lat) { + case FCC: + *rho01 = 12.0 * rhoa02; + *rho02 = 12.0 * rhoa01; + break; + case BCC: + *rho01 = 8.0 * rhoa02; + *rho02 = 8.0 * rhoa01; + break; + case B1: + *rho01 = 6.0 * rhoa02; + *rho02 = 6.0 * rhoa01; + break; + case DIA: + case DIA3: + *rho01 = 4.0 * rhoa02; + *rho02 = 4.0 * rhoa01; + *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32; + *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31; + break; + case HCP: + *rho01 = 12 * rhoa02; + *rho02 = 12 * rhoa01; + *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32; + *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31; + break; + case DIM: + get_shpfcn(DIM, 0, 0, s); + *rho01 = rhoa02; + *rho02 = rhoa01; + *rho11 = s[0] * rhoa12 * rhoa12; + *rho12 = s[0] * rhoa11 * rhoa11; + *rho21 = s[1] * rhoa22 * rhoa22; + *rho22 = s[1] * rhoa21 * rhoa21; + *rho31 = s[2] * rhoa32 * rhoa32; + *rho32 = s[2] * rhoa31 * rhoa31; + break; + case C11: + *rho01 = rhoa01; + *rho02 = rhoa02; + *rho11 = rhoa11; + *rho12 = rhoa12; + *rho21 = rhoa21; + *rho22 = rhoa22; + *rho31 = rhoa31; + *rho32 = rhoa32; + break; + case L12: + *rho01 = 8 * rhoa01 + 4 * rhoa02; + *rho02 = 12 * rhoa01; + if (this->ialloy == 1) { + *rho21 = 8. / 3. * MathSpecial::square(rhoa21 * this->t2_meam[a] - rhoa22 * this->t2_meam[b]); + denom = 8 * rhoa01 * MathSpecial::square(this->t2_meam[a]) + 4 * rhoa02 * MathSpecial::square(this->t2_meam[b]); + if (denom > 0.) + *rho21 = *rho21 / denom * *rho01; + } else + *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22); + break; + case B2: + *rho01 = 8.0 * rhoa02; + *rho02 = 8.0 * rhoa01; + break; + case CH4: + *rho01 = 4.0 * rhoa02; //in assumption that 'a' represent carbon + *rho02 = rhoa01; //in assumption that 'b' represent hydrogen + + get_shpfcn(DIM, 0, 0, s); //H + *rho12 = s[0] * rhoa11 * rhoa11; + *rho22 = s[1] * rhoa21 * rhoa21; + *rho32 = s[2] * rhoa31 * rhoa31; + + get_shpfcn(CH4, 0, 0, s); //C + *rho11 = s[0] * rhoa12 * rhoa12; + *rho21 = s[1] * rhoa22 * rhoa22; + *rho31 = s[2] * rhoa32 * rhoa32; + break; + case LIN: + *rho01 = rhoa02*Zij; + *rho02 = rhoa01*Zij; + + get_shpfcn(LIN, this->stheta_meam[a][b], this->ctheta_meam[a][b], s); + *rho12 = s[0] * rhoa11 * rhoa11; + *rho22 = s[1] * rhoa21 * rhoa21; + *rho32 = s[2] * rhoa31 * rhoa31; + *rho11 = s[0] * rhoa12 * rhoa12; + *rho21 = s[1] * rhoa22 * rhoa22; + *rho31 = s[2] * rhoa32 * rhoa32; + break; + case ZIG: + *rho01 = rhoa02*Zij; + *rho02 = rhoa01*Zij; + + get_shpfcn(ZIG, this->stheta_meam[a][b], this->ctheta_meam[a][b], s); + *rho12 = s[0] * rhoa11 * rhoa11; + *rho22 = s[1] * rhoa21 * rhoa21; + *rho32 = s[2] * rhoa31 * rhoa31; + *rho11 = s[0] * rhoa12 * rhoa12; + *rho21 = s[1] * rhoa22 * rhoa22; + *rho31 = s[2] * rhoa32 * rhoa32; + break; + case TRI: + *rho01 = rhoa02; + *rho02 = rhoa01*Zij; + + get_shpfcn(TRI, this->stheta_meam[a][b], this->ctheta_meam[a][b], s); + *rho12 = s[0] * rhoa11 * rhoa11; + *rho22 = s[1] * rhoa21 * rhoa21; + *rho32 = s[2] * rhoa31 * rhoa31; + s[0] = 1.0; + s[1] = 2.0/3.0; + s[2] = 1.0 - 0.6*s[0]; + + *rho11 = s[0] * rhoa12 * rhoa12; + *rho21 = s[1] * rhoa22 * rhoa22; + *rho31 = s[2] * rhoa32 * rhoa32; + break; + + + // default: // call error('Lattice not defined in get_densref.') } if (this->nn2_meam[a][b] == 1) { - Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], arat, scrn); + + Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], + this->stheta_meam[a][b], arat, scrn); a1 = arat * r / this->re_meam[a][a] - 1.0; a2 = arat * r / this->re_meam[b][b] - 1.0; @@ -712,7 +836,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* *rho01 = *rho01 + Zij2nn * scrn * rhoa01nn; // Assume Zij2nn and arat don't depend on order, but scrn might - Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], arat, scrn); + Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], + this->stheta_meam[a][b], arat, scrn); *rho02 = *rho02 + Zij2nn * scrn * rhoa02nn; } } diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index b66aafa06c..ac5718914c 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -18,7 +18,7 @@ static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) { } void -MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* /*atwt*/, double* alpha, +MEAM::meam_setup_global(int nelt, lattice_t* lat, int* ielement, double* /*atwt*/, double* alpha, double* b0, double* b1, double* b2, double* b3, double* alat, double* esub, double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero, int* ibar) @@ -32,7 +32,6 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub for (i = 0; i < nelt; i++) { this->lattce_meam[i][i] = lat[i]; - this->Z_meam[i] = z[i]; this->ielt_meam[i] = ielement[i]; this->alpha_meam[i][i] = alpha[i]; this->beta0_meam[i] = b0[i]; @@ -49,17 +48,26 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub this->rho0_meam[i] = rozero[i]; this->ibar_meam[i] = ibar[i]; - if (this->lattce_meam[i][i] == FCC) - this->re_meam[i][i] = tmplat[i] / sqrt(2.0); - else if (this->lattce_meam[i][i] == BCC) - this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0; - else if (this->lattce_meam[i][i] == HCP) - this->re_meam[i][i] = tmplat[i]; - else if (this->lattce_meam[i][i] == DIM) - this->re_meam[i][i] = tmplat[i]; - else if (this->lattce_meam[i][i] == DIA) - this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0; - else { + switch(this->lattce_meam[i][i]) { + case FCC: + this->re_meam[i][i] = tmplat[i] / sqrt(2.0); + break; + case BCC: + this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 2.0; + break; + case HCP: + case DIM: + case CH4: + case LIN: + case ZIG: + case TRI: + this->re_meam[i][i] = tmplat[i]; + break; + case DIA: + case DIA3: + this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0; + break; + //default: // error } } @@ -82,4 +90,7 @@ MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, doub this->emb_lin_neg = 0; this->bkgd_dyn = 0; this->erose_form = 0; + // for trimer, zigzag, line refernece structure, sungkwang + setall2d(this->stheta_meam, 1.0); // stheta = sin(theta/2*pi/180) where theta is 180, so 1.0 + setall2d(this->ctheta_meam, 0.0); // stheta = cos(theta/2*pi/180) where theta is 180, so 0 } diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index 585b5b5712..70d8b63598 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -1,6 +1,7 @@ #include "meam.h" #include using namespace LAMMPS_NS; +using namespace MathConst; // // do a sanity check on index parameters @@ -46,6 +47,8 @@ MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr // 18 = zbl_meam // 19 = emb_lin_neg // 20 = bkgd_dyn +// 21 = theta + void MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag) @@ -203,6 +206,20 @@ MEAM::meam_setup_param(int which, double value, int nindex, int* index /*index(3 this->bkgd_dyn = (int)value; break; + // 21 = theta + // see alloyparams(void) in meam_setup_done.cpp + case 21: + // const double PI = 3.141592653589793238463; + meam_checkindex(2, neltypes, nindex, index, errorflag); + if (*errorflag != 0) + return; + i1 = std::min(index[0], index[1]); + i2 = std::max(index[0], index[1]); + // we don't use theta, instead stheta and ctheta + this->stheta_meam[i1][i2] = sin(value/2*MY_PI/180.0); + this->ctheta_meam[i1][i2] = cos(value/2*MY_PI/180.0); + break; + default: *errorflag = 1; } diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index a9034a1af3..9733ee2dad 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -33,13 +33,13 @@ using namespace LAMMPS_NS; #define MAXLINE 1024 -static const int nkeywords = 21; +static const int nkeywords = 22; static const char *keywords[] = { "Ec","alpha","rho0","delta","lattce", "attrac","repuls","nn2","Cmin","Cmax","rc","delr", "augt1","gsmooth_factor","re","ialloy", "mixture_ref_t","erose_form","zbl", - "emb_lin_neg","bkgd_dyn"}; + "emb_lin_neg","bkgd_dyn", "theta"}; /* ---------------------------------------------------------------------- */ @@ -56,6 +56,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp) elements = NULL; mass = NULL; meam_inst = new MEAM(memory); + scale = NULL; // set comm size needed by this Pair @@ -79,6 +80,7 @@ PairMEAMC::~PairMEAMC() if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); + memory->destroy(scale); delete [] map; } } @@ -143,7 +145,7 @@ void PairMEAMC::compute(int eflag, int vflag) comm->reverse_comm_pair(this); meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom, - &eng_vdwl,eatom,ntype,type,map,errorflag); + &eng_vdwl,eatom,ntype,type,map,scale,errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); @@ -164,7 +166,7 @@ void PairMEAMC::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; meam_inst->meam_force(i,eflag_either,eflag_global,eflag_atom, - vflag_atom,&eng_vdwl,eatom,ntype,type,map,x, + vflag_atom,&eng_vdwl,eatom,ntype,type,map,scale,x, numneigh_half[i],firstneigh_half[i], numneigh_full[i],firstneigh_full[i], offset,f,vptr); @@ -183,6 +185,7 @@ void PairMEAMC::allocate() memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(scale,n+1,n+1,"pair:scale"); map = new int[n+1]; } @@ -267,13 +270,16 @@ void PairMEAMC::coeff(int narg, char **arg) // set mass for i,i in atom class int count = 0; - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) + for (int i = 1; i <= n; i++) { + for (int j = i; j <= n; j++) { if (map[i] >= 0 && map[j] >= 0) { setflag[i][j] = 1; if (i == j) atom->set_mass(FLERR,i,mass[map[i]]); count++; } + scale[i][j] = 1.0; + } + } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -312,8 +318,10 @@ void PairMEAMC::init_list(int id, NeighList *ptr) init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairMEAMC::init_one(int /*i*/, int /*j*/) +double PairMEAMC::init_one(int i, int j) { + if (setflag[i][j] == 0) scale[i][j] = 1.0; + scale[j][i] = scale[i][j]; return cutmax; } @@ -430,13 +438,9 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) found[i] = true; // map lat string to an integer - - if (strcmp(words[1],"fcc") == 0) lat[i] = FCC; - else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC; - else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP; - else if (strcmp(words[1],"dim") == 0) lat[i] = DIM; - else if (strcmp(words[1],"dia") == 0) lat[i] = DIA; - else error->all(FLERR,"Unrecognized lattice type in MEAM file 1"); + + if (!MEAM::str_to_lat(words[1], true, lat[i])) + error->all(FLERR,"Unrecognized lattice type in MEAM file 1"); // store parameters @@ -458,8 +462,12 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) rozero[i] = atof(words[17]); ibar[i] = atoi(words[18]); - if (!iszero(t0[i]-1.0)) - error->all(FLERR,"Unsupported parameter in MEAM potential file"); + if (!isone(t0[i])) + error->all(FLERR,"Unsupported parameter in MEAM potential file: t0!=1"); + + // z given is ignored: if this is mismatched, we definitely won't do what the user said -> fatal error + if (z[i] != MEAM::get_Zij(lat[i])) + error->all(FLERR,"Mismatched parameter in MEAM potential file: z!=lat"); nset++; } @@ -471,7 +479,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) // pass element parameters to MEAM package - meam_inst->meam_setup_global(nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, + meam_inst->meam_setup_global(nelements,lat,ielement,atwt,alpha,b0,b1,b2,b3, alat,esub,asub,t0,t1,t2,t3,rozero,ibar); // set element masses @@ -523,6 +531,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) int which; double value; + lattice_t latt; int nindex,index[3]; int maxparams = 6; char **params = new char*[maxparams]; @@ -571,16 +580,9 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) // map lattce_meam value to an integer if (which == 4) { - if (strcmp(params[nparams-1],"fcc") == 0) value = FCC; - else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC; - else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP; - else if (strcmp(params[nparams-1],"dim") == 0) value = DIM; - else if (strcmp(params[nparams-1],"dia") == 0) value = DIA; - else if (strcmp(params[nparams-1],"b1") == 0) value = B1; - else if (strcmp(params[nparams-1],"c11") == 0) value = C11; - else if (strcmp(params[nparams-1],"l12") == 0) value = L12; - else if (strcmp(params[nparams-1],"b2") == 0) value = B2; - else error->all(FLERR,"Unrecognized lattice type in MEAM file 2"); + if (!MEAM::str_to_lat(params[nparams-1], false, latt)) + error->all(FLERR,"Unrecognized lattice type in MEAM file 2"); + value = latt; } else value = atof(params[nparams-1]); @@ -783,3 +785,12 @@ void PairMEAMC::neigh_strip(int inum, int *ilist, for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK; } } + +/* ---------------------------------------------------------------------- */ + +void *PairMEAMC::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str,"scale") == 0) return (void *) scale; + return NULL; +} diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h index 31dd8ba022..102985f0ca 100644 --- a/src/USER-MEAMC/pair_meamc.h +++ b/src/USER-MEAMC/pair_meamc.h @@ -35,6 +35,7 @@ class PairMEAMC : public Pair { void init_style(); void init_list(int, class NeighList *); double init_one(int, int); + virtual void *extract(const char *, int &); int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); @@ -50,6 +51,7 @@ class PairMEAMC : public Pair { double *mass; // mass of each element int *map; // mapping from atom types (1-indexed) to elements (1-indexed) + double **scale; // scaling factor for adapt void allocate(); void read_files(char *, char *); diff --git a/src/math_special.cpp b/src/math_special.cpp index bf11a1ad45..d4abc36f25 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -538,6 +538,8 @@ double MathSpecial::exp2_x86(double x) double MathSpecial::fm_exp(double x) { #if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__) + if (x < -1022.0/FM_DOUBLE_LOG2OFE) return 0; + if (x > 1023.0/FM_DOUBLE_LOG2OFE) return INFINITY; return exp2_x86(FM_DOUBLE_LOG2OFE * x); #else return ::exp(x);