diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 9be4628022..698f620694 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -2,6 +2,7 @@ #define LMP_MEAM_H #include +#include #define maxelt 5 @@ -26,7 +27,6 @@ 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 @@ -65,7 +65,7 @@ private: // nrar,rdrar = spline coeff array parameters double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt]; - double Omega_meam[maxelt], Z_meam[maxelt]; + double 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]; @@ -212,6 +212,31 @@ protected: void interpolate_meam(int); public: + //----------------------------------------------------------------------------- + // 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 (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 return false; + } + return true; + } + void meam_setup_global(int nelt, lattice_t* 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, diff --git a/src/USER-MEAMC/meam_impl.cpp b/src/USER-MEAMC/meam_impl.cpp index 74e8af1cde..96ef364fa3 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] = + Z_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..1c0fc3a53c 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -161,10 +161,10 @@ void MEAM::compute_pair_meam(void) { - double r /*ununsed:, temp*/; + double r; 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; @@ -230,9 +230,8 @@ MEAM::compute_pair_meam(void) 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++) { + for (n = 1; n <= 10; n++) { phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a); } } @@ -242,9 +241,8 @@ MEAM::compute_pair_meam(void) 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++) { + for (n = 1; n <= 10; n++) { phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b); } } @@ -279,8 +277,7 @@ MEAM::compute_pair_meam(void) } } else { - nmax = 10; - for (n = 1; n <= nmax; n++) { + for (n = 1; n <= 10; n++) { this->phir[nv2][j] = this->phir[nv2][j] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b); } @@ -330,7 +327,7 @@ MEAM::phi_meam(double r, int a, int b) double Eu; double arat, scrn /*unused:,scrn2*/; int Z12, errorflag; - int n, nmax, Z1nn, Z2nn; + int n, Z1nn, Z2nn; lattice_t latta /*unused:,lattb*/; double rho_bkgd1, rho_bkgd2; @@ -476,9 +473,8 @@ MEAM::phi_meam(double r, int a, int b) 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++) { + for (n = 1; n <= 10; n++) { phiaa = phiaa + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a); } } @@ -555,30 +551,38 @@ 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 HCP: + case B1: + case DIM: + case B2: + // 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.') + } } } @@ -627,59 +631,69 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* *rho12 = 0.0; *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: + *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, 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; + // default: // call error('Lattice not defined in get_densref.') } diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index b66aafa06c..2306ac53c6 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -49,17 +49,23 @@ 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: + this->re_meam[i][i] = tmplat[i]; + break; + case DIM: + this->re_meam[i][i] = tmplat[i]; + break; + case DIA: + this->re_meam[i][i] = tmplat[i] * sqrt(3.0) / 4.0; + break; + //default: // error } } diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index dc2e194dea..bc366854b6 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -430,13 +430,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 @@ -523,6 +519,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 +568,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]);