ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

MEAM/C: remove unused vars, refactoring for extensibility

Browse Source
This commit is contained in:
Sebastian Hütter
2019-08-21 18:38:45 +02:00
parent cd6d2c55d1
commit 204529bcaf
5 changed files with 151 additions and 116 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
src/USER-MEAMC/meam.h
View File

@ -2,6 +2,7 @@
#define LMP_MEAM_H
#include <cmath>
#include <cstring>
#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,

2
src/USER-MEAMC/meam_impl.cpp
View File

@ -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;

184
src/USER-MEAMC/meam_setup_done.cpp
View File

@ -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.')
}

28
src/USER-MEAMC/meam_setup_global.cpp
View File

@ -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
}
}

24
src/USER-MEAMC/pair_meamc.cpp
View File

@ -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]);