forked from lijiext/lammps
Merge pull request #842 from martok/meamc-bugfix
MEAM/C Improvements/fixes
This commit is contained in:
Steve Plimpton
GitHub
commit
db5a6b272b
|
|
@ -208,7 +208,6 @@ protected:
|
|||
void get_sijk(double, int, int, int, double*);
|
||||
void get_densref(double, int, int, double*, double*, double*, double*, double*, double*, double*, double*);
|
||||
void interpolate_meam(int);
|
||||
double compute_phi(double, int, int);
|
||||
|
||||
public:
|
||||
void meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha,
|
||||
|
|
@ -248,5 +247,13 @@ static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) {
|
|||
arr[i][j][k] = v;
|
||||
}
|
||||
|
||||
// Helper functions
|
||||
|
||||
static inline double fdiv_zero(const double n, const double d) {
|
||||
if (iszero(d))
|
||||
return 0.0;
|
||||
return n / d;
|
||||
}
|
||||
|
||||
};
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -33,9 +33,9 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
|
|||
|
||||
if (rho0[i] > 0.0) {
|
||||
if (this->ialloy == 1) {
|
||||
t_ave[i][0] = (tsq_ave[i][0] != 0.0 ) ? t_ave[i][0] / tsq_ave[i][0] : 0.0;
|
||||
t_ave[i][1] = (tsq_ave[i][1] != 0.0 ) ? t_ave[i][1] / tsq_ave[i][1] : 0.0;
|
||||
t_ave[i][0] = (tsq_ave[i][2] != 0.0 ) ? t_ave[i][2] / tsq_ave[i][2] : 0.0;
|
||||
t_ave[i][0] = fdiv_zero(t_ave[i][0], tsq_ave[i][0]);
|
||||
t_ave[i][1] = fdiv_zero(t_ave[i][1], tsq_ave[i][1]);
|
||||
t_ave[i][2] = fdiv_zero(t_ave[i][2], tsq_ave[i][2]);
|
||||
} else if (this->ialloy == 2) {
|
||||
t_ave[i][0] = this->t1_meam[elti];
|
||||
t_ave[i][1] = this->t2_meam[elti];
|
||||
|
|
|
|||
|
|
@ -245,31 +245,19 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
|
|||
|
||||
if (this->ialloy == 1) {
|
||||
|
||||
a1i = 0.0;
|
||||
a1j = 0.0;
|
||||
a2i = 0.0;
|
||||
a2j = 0.0;
|
||||
a3i = 0.0;
|
||||
a3j = 0.0;
|
||||
if (!iszero(tsq_ave[i][0]))
|
||||
a1i = drhoa0j * sij / tsq_ave[i][0];
|
||||
if (!iszero(tsq_ave[j][0]))
|
||||
a1j = drhoa0i * sij / tsq_ave[j][0];
|
||||
if (!iszero(tsq_ave[i][1]))
|
||||
a2i = drhoa0j * sij / tsq_ave[i][1];
|
||||
if (!iszero(tsq_ave[j][1]))
|
||||
a2j = drhoa0i * sij / tsq_ave[j][1];
|
||||
if (!iszero(tsq_ave[i][2]))
|
||||
a3i = drhoa0j * sij / tsq_ave[i][2];
|
||||
if (!iszero(tsq_ave[j][2]))
|
||||
a3j = drhoa0i * sij / tsq_ave[j][2];
|
||||
a1i = fdiv_zero(drhoa0j * sij, tsq_ave[i][0]);
|
||||
a1j = fdiv_zero(drhoa0i * sij, tsq_ave[j][0]);
|
||||
a2i = fdiv_zero(drhoa0j * sij, tsq_ave[i][1]);
|
||||
a2j = fdiv_zero(drhoa0i * sij, tsq_ave[j][1]);
|
||||
a3i = fdiv_zero(drhoa0j * sij, tsq_ave[i][2]);
|
||||
a3j = fdiv_zero(drhoa0i * sij, tsq_ave[j][2]);
|
||||
|
||||
dt1dr1 = a1i * (t1mj - t1i * t1mj*t1mj);
|
||||
dt1dr2 = a1j * (t1mi - t1j * t1mi*t1mi);
|
||||
dt2dr1 = a2i * (t2mj - t2i * t2mj*t2mj);
|
||||
dt2dr2 = a2j * (t2mi - t2j * t2mi*t2mi);
|
||||
dt3dr1 = a3i * (t3mj - t3i * t3mj*t3mj);
|
||||
dt3dr2 = a3j * (t3mi - t3j * t3mi*t3mi);
|
||||
dt1dr1 = a1i * (t1mj - t1i * MathSpecial::square(t1mj));
|
||||
dt1dr2 = a1j * (t1mi - t1j * MathSpecial::square(t1mi));
|
||||
dt2dr1 = a2i * (t2mj - t2i * MathSpecial::square(t2mj));
|
||||
dt2dr2 = a2j * (t2mi - t2j * MathSpecial::square(t2mi));
|
||||
dt3dr1 = a3i * (t3mj - t3i * MathSpecial::square(t3mj));
|
||||
dt3dr2 = a3j * (t3mi - t3j * MathSpecial::square(t3mi));
|
||||
|
||||
} else if (this->ialloy == 2) {
|
||||
|
||||
|
|
@ -331,32 +319,19 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
|
|||
drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3;
|
||||
|
||||
if (this->ialloy == 1) {
|
||||
a1i = fdiv_zero(rhoa0j, tsq_ave[i][0]);
|
||||
a1j = fdiv_zero(rhoa0i, tsq_ave[j][0]);
|
||||
a2i = fdiv_zero(rhoa0j, tsq_ave[i][1]);
|
||||
a2j = fdiv_zero(rhoa0i, tsq_ave[j][1]);
|
||||
a3i = fdiv_zero(rhoa0j, tsq_ave[i][2]);
|
||||
a3j = fdiv_zero(rhoa0i, tsq_ave[j][2]);
|
||||
|
||||
a1i = 0.0;
|
||||
a1j = 0.0;
|
||||
a2i = 0.0;
|
||||
a2j = 0.0;
|
||||
a3i = 0.0;
|
||||
a3j = 0.0;
|
||||
if (!iszero(tsq_ave[i][0]))
|
||||
a1i = rhoa0j / tsq_ave[i][0];
|
||||
if (!iszero(tsq_ave[j][0]))
|
||||
a1j = rhoa0i / tsq_ave[j][0];
|
||||
if (!iszero(tsq_ave[i][1]))
|
||||
a2i = rhoa0j / tsq_ave[i][1];
|
||||
if (!iszero(tsq_ave[j][1]))
|
||||
a2j = rhoa0i / tsq_ave[j][1];
|
||||
if (!iszero(tsq_ave[i][2]))
|
||||
a3i = rhoa0j / tsq_ave[i][2];
|
||||
if (!iszero(tsq_ave[j][2]))
|
||||
a3j = rhoa0i / tsq_ave[j][2];
|
||||
|
||||
dt1ds1 = a1i * (t1mj - t1i * pow(t1mj, 2));
|
||||
dt1ds2 = a1j * (t1mi - t1j * pow(t1mi, 2));
|
||||
dt2ds1 = a2i * (t2mj - t2i * pow(t2mj, 2));
|
||||
dt2ds2 = a2j * (t2mi - t2j * pow(t2mi, 2));
|
||||
dt3ds1 = a3i * (t3mj - t3i * pow(t3mj, 2));
|
||||
dt3ds2 = a3j * (t3mi - t3j * pow(t3mi, 2));
|
||||
dt1ds1 = a1i * (t1mj - t1i * MathSpecial::square(t1mj));
|
||||
dt1ds2 = a1j * (t1mi - t1j * MathSpecial::square(t1mi));
|
||||
dt2ds1 = a2i * (t2mj - t2i * MathSpecial::square(t2mj));
|
||||
dt2ds2 = a2j * (t2mi - t2j * MathSpecial::square(t2mi));
|
||||
dt3ds1 = a3i * (t3mj - t3i * MathSpecial::square(t3mj));
|
||||
dt3ds2 = a3j * (t3mi - t3j * MathSpecial::square(t3mi));
|
||||
|
||||
} else if (this->ialloy == 2) {
|
||||
|
||||
|
|
|
|||
|
|
@ -162,11 +162,11 @@ MEAM::erose(const double r, const double re, const double alpha, const double Ec
|
|||
a3 = repuls;
|
||||
|
||||
if (form == 1)
|
||||
result = -Ec * (1 + astar + (-attrac + repuls / r) * pow(astar, 3)) * MathSpecial::fm_exp(-astar);
|
||||
result = -Ec * (1 + astar + (-attrac + repuls / r) * MathSpecial::cube(astar)) * MathSpecial::fm_exp(-astar);
|
||||
else if (form == 2)
|
||||
result = -Ec * (1 + astar + a3 * pow(astar, 3)) * MathSpecial::fm_exp(-astar);
|
||||
result = -Ec * (1 + astar + a3 * MathSpecial::cube(astar)) * MathSpecial::fm_exp(-astar);
|
||||
else
|
||||
result = -Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * MathSpecial::fm_exp(-astar);
|
||||
result = -Ec * (1 + astar + a3 * MathSpecial::cube(astar) / (r / re)) * MathSpecial::fm_exp(-astar);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -375,16 +375,20 @@ MEAM::phi_meam(double r, int a, int b)
|
|||
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 = MathSpecial::square((Z12 / 2) * (rho02 + rho01))
|
||||
+ t11av * MathSpecial::square(rho12 - rho11)
|
||||
+ t21av / 6.0 * MathSpecial::square(rho22 + rho21)
|
||||
+ 121.0 / 40.0 * t31av * MathSpecial::square(rho32 - rho31);
|
||||
rhobar1 = sqrt(rhobar1);
|
||||
rhobar2 = pow(Z12 * rho01, 2) + 2.0 / 3.0 * t21av * pow(rho21, 2);
|
||||
rhobar2 = MathSpecial::square(Z12 * rho01) + 2.0 / 3.0 * t21av * MathSpecial::square(rho21);
|
||||
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 = MathSpecial::square((Z12 / 2) * (rho01 + rho02))
|
||||
+ t12av * MathSpecial::square(rho11 - rho12)
|
||||
+ t22av / 6.0 * MathSpecial::square(rho21 + rho22)
|
||||
+ 121.0 / 40.0 * t32av * MathSpecial::square(rho31 - rho32);
|
||||
rhobar2 = sqrt(rhobar2);
|
||||
rhobar1 = pow(Z12 * rho02, 2) + 2.0 / 3.0 * t22av * pow(rho22, 2);
|
||||
rhobar1 = MathSpecial::square(Z12 * rho02) + 2.0 / 3.0 * t22av * MathSpecial::square(rho22);
|
||||
rhobar1 = sqrt(rhobar1);
|
||||
}
|
||||
} else {
|
||||
|
|
@ -560,31 +564,29 @@ 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 {
|
||||
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;
|
||||
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 {
|
||||
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.')
|
||||
}
|
||||
// call error('Lattice not defined in get_tavref.')
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -677,8 +679,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
|
|||
*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);
|
||||
*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
|
||||
|
|
@ -768,25 +770,3 @@ MEAM::interpolate_meam(int ind)
|
|||
this->phirar6[ind][j] = 3.0 * this->phirar3[ind][j] / 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;
|
||||
kk = (int)pp;
|
||||
kk = std::min(kk, this->nrar - 2);
|
||||
pp = pp - kk;
|
||||
pp = std::min(pp, 1.0);
|
||||
double result =
|
||||
((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp +
|
||||
this->phirar[ind][kk];
|
||||
|
||||
return result;
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue