jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1798 from martok/meamc-dec19 

Collected MEAM/C additions
This commit is contained in:
Axel Kohlmeyer 2020-01-06 15:07:08 -05:00 committed by GitHub
parent 402f5585ff cb20cb9f25
commit 1be8109618
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 516 additions and 202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
doc/src/pair_meamc.rst
View File

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

2
doc/txt/pair_meamc.txt
View File

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

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

@ -2,13 +2,15 @@
#define LMP_MEAM_H
#include <cmath>
#include <cstring>
#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) {

11
src/USER-MEAMC/meam_dens_final.cpp
View File

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

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

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

35
src/USER-MEAMC/meam_force.cpp
View File

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

96
src/USER-MEAMC/meam_funcs.cpp
View File

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

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] =
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;

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

@ -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<j and term is unset, use default values (e.g. mean of i-i and
// j-j)
} else if (j > 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;
}
}

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

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

17
src/USER-MEAMC/meam_setup_param.cpp
View File

@ -1,6 +1,7 @@
#include "meam.h"
#include <algorithm>
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;
}

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

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

2
src/USER-MEAMC/pair_meamc.h
View File

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

2
src/math_special.cpp
View File

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