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

Move rho/gamma arrays to fields of MEAM, remove arguments and arrdim macros

Browse Source
This commit is contained in:
Sebastian Hütter
2017-06-13 18:39:40 +02:00
parent 078f2a0a47
commit 29e60fa53a
7 changed files with 409 additions and 418 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -12,12 +12,8 @@ typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t;
class MEAM {
public:
MEAM(Memory *mem) :
memory(mem) {};
~MEAM() {
meam_cleanup();
}
MEAM(Memory *mem);
~MEAM();
private:
Memory *&memory;
@ -102,14 +98,29 @@ class MEAM {
int nr, nrar;
double dr, rdrar;
public:
int nmax;
double *rho,*rho0,*rho1,*rho2,*rho3,*frhop;
double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b;
double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave;
int maxneigh;
double *scrfcn,*dscrfcn,*fcpair;
protected:
void meam_checkindex(int, int, int, int*, int*);
void G_gam(double, int, double, double*, int*);
void dG_gam(double, int, double, double*, double*);
void getscreen(int, int, double*, double*, double*, double*, int, int*, int, int*, int, int*, int*);
void screen(int, int, int, double*, double, double*, int, int*, int, int*, int*);
void calc_rho1(int, int, int, int*, int*, double*, int, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*);
void dsij(int, int, int, int, int, int, double, double*, double*, int, int*, int*, double*, double*, double*);
void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair,
double** x, int numneigh, int* firstneigh, int numneigh_full,
int* firstneigh_full, int ntype, int* type, int* fmap);
void screen(int i, int j, double** x, double rijsq, double* sij,
int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap);
void calc_rho1(int i, int ntype, int* type, int* fmap, double** x,
int numneigh, int* firstneigh, double* scrfcn, double* fcpair);
void dsij(int i, int j, int k, int jn, int numneigh, double rij2,
double* dsij1, double* dsij2, int ntype, int* type, int* fmap, double** x,
double* scrfcn, double* fcpair);
void fcut(double, double*);
void dfcut(double, double*, double*);
void dCfunc(double, double, double, double*);
@ -133,9 +144,19 @@ class MEAM {
void meam_setup_global(int*, int*, double*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*);
void meam_setup_param(int*, double*, int*, int*, int*);
void meam_setup_done(double*);
void meam_dens_init(int*, int*, int*, int*, int*, double*, int*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*);
void meam_dens_final(int*, int*, int*, int*, int*, double*, double*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*);
void meam_force(int*, int*, int*, int*, int*, int*, double*, double*, int*, int*, int*, double*, int*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*);
void meam_dens_setup(int, int, int);
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, int* errorflag);
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);
void meam_force(int* iptr, 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, int* numneigh_full, int* firstneigh_full,
int fnoffset, double** f, double** vatom,
int* errorflag);
void meam_cleanup();
};
@ -165,25 +186,8 @@ Fortran Array Semantics in C.
(or be used with 0-based indexing)
*/
// we receive a pointer to the first element, and F dimensions is ptr(a,b,c)
// we know c data structure is ptr[c][b][a]
#define arrdim2v(ptr, a, b) \
const int DIM1__##ptr = a; \
const int DIM2__##ptr = b; \
(void)(DIM1__##ptr); \
(void)(DIM2__##ptr);
#define arrdim3v(ptr, a, b, c) \
const int DIM1__##ptr = a; \
const int DIM2__##ptr = b; \
const int DIM3__##ptr = c; \
(void)(DIM1__##ptr); (void)(DIM2__##ptr; (void)(DIM3__##ptr);
// access data with same index as used in fortran (1-based)
#define arr1v(ptr, i) ptr[i - 1]
#define arr2v(ptr, i, j) ptr[(DIM1__##ptr) * (j - 1) + (i - 1)]
#define arr3v(ptr, i, j, k) \
ptr[(i - 1) + (j - 1) * (DIM1__##ptr) + \
(k - 1) * (DIM1__##ptr) * (DIM2__##ptr)]
};
#define arr2v(ptr, i, j) ptr[j - 1][i - 1]
#endif

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

@ -21,23 +21,10 @@ using namespace LAMMPS_NS;
//
void
MEAM::meam_dens_final(int* nlocal, int* nmax, int* eflag_either, int* eflag_global,
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, double* Arho1, double* Arho2,
double* Arho2b, double* Arho3, double* Arho3b, double* t_ave,
double* tsq_ave, double* Gamma, double* dGamma1,
double* dGamma2, double* dGamma3, double* rho, double* rho0,
double* rho1, double* rho2, double* rho3, double* fp,
int* errorflag)
int* type, int* fmap, int* errorflag)
{
arrdim2v(Arho1, 3, *nmax);
arrdim2v(Arho2, 6, *nmax);
arrdim2v(Arho3, 10, *nmax);
arrdim2v(Arho3b, 3, *nmax);
arrdim2v(t_ave, 3, *nmax);
arrdim2v(tsq_ave, 3, *nmax);
int i, elti;
int m;
double rhob, G, dG, Gbar, dGbar, gam, shp[3 + 1], Z;
@ -49,23 +36,23 @@ MEAM::meam_dens_final(int* nlocal, int* nmax, int* eflag_either, int* eflag_glob
elti = arr1v(fmap, arr1v(type, i));
if (elti > 0) {
arr1v(rho1, i) = 0.0;
arr1v(rho2, i) = -1.0 / 3.0 * arr1v(Arho2b, i) * arr1v(Arho2b, i);
arr1v(rho2, i) = -1.0 / 3.0 * arr1v(arho2b, i) * arr1v(arho2b, i);
arr1v(rho3, i) = 0.0;
for (m = 1; m <= 3; m++) {
arr1v(rho1, i) =
arr1v(rho1, i) + arr2v(Arho1, m, i) * arr2v(Arho1, m, i);
arr1v(rho1, i) + arr2v(arho1, m, i) * arr2v(arho1, m, i);
arr1v(rho3, i) = arr1v(rho3, i) -
3.0 / 5.0 * arr2v(Arho3b, m, i) * arr2v(Arho3b, m, i);
3.0 / 5.0 * arr2v(arho3b, m, i) * arr2v(arho3b, m, i);
}
for (m = 1; m <= 6; m++) {
arr1v(rho2, i) =
arr1v(rho2, i) +
this->v2D[m] * arr2v(Arho2, m, i) * arr2v(Arho2, m, i);
this->v2D[m] * arr2v(arho2, m, i) * arr2v(arho2, m, i);
}
for (m = 1; m <= 10; m++) {
arr1v(rho3, i) =
arr1v(rho3, i) +
this->v3D[m] * arr2v(Arho3, m, i) * arr2v(Arho3, m, i);
this->v3D[m] * arr2v(arho3, m, i) * arr2v(arho3, m, i);
}
if (arr1v(rho0, i) > 0.0) {
@ -84,17 +71,17 @@ MEAM::meam_dens_final(int* nlocal, int* nmax, int* eflag_either, int* eflag_glob
}
}
arr1v(Gamma, i) = arr2v(t_ave, 1, i) * arr1v(rho1, i) +
arr1v(gamma, i) = arr2v(t_ave, 1, i) * arr1v(rho1, i) +
arr2v(t_ave, 2, i) * arr1v(rho2, i) +
arr2v(t_ave, 3, i) * arr1v(rho3, i);
if (arr1v(rho0, i) > 0.0) {
arr1v(Gamma, i) = arr1v(Gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i));
arr1v(gamma, i) = arr1v(gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i));
}
Z = this->Z_meam[elti];
G_gam(arr1v(Gamma, i), this->ibar_meam[elti],
G_gam(arr1v(gamma, i), this->ibar_meam[elti],
this->gsmooth_factor, &G, errorflag);
if (*errorflag != 0)
return;
@ -140,33 +127,33 @@ MEAM::meam_dens_final(int* nlocal, int* nmax, int* eflag_either, int* eflag_glob
rhob = arr1v(rho, i) / rho_bkgd;
denom = 1.0 / rho_bkgd;
dG_gam(arr1v(Gamma, i), this->ibar_meam[elti],
dG_gam(arr1v(gamma, i), this->ibar_meam[elti],
this->gsmooth_factor, &G, &dG);
arr1v(dGamma1, i) = (G - 2 * dG * arr1v(Gamma, i)) * denom;
arr1v(dgamma1, i) = (G - 2 * dG * arr1v(gamma, i)) * denom;
if (!iszero(arr1v(rho0, i))) {
arr1v(dGamma2, i) = (dG / arr1v(rho0, i)) * denom;
arr1v(dgamma2, i) = (dG / arr1v(rho0, i)) * denom;
} else {
arr1v(dGamma2, i) = 0.0;
arr1v(dgamma2, i) = 0.0;
}
// dGamma3 is nonzero only if we are using the "mixed" rule for
// dgamma3 is nonzero only if we are using the "mixed" rule for
// computing t in the reference system (which is not correct, but
// included for backward compatibility
if (this->mix_ref_t == 1) {
arr1v(dGamma3, i) = arr1v(rho0, i) * G * dGbar / (Gbar * Z * Z) * denom;
arr1v(dgamma3, i) = arr1v(rho0, i) * G * dGbar / (Gbar * Z * Z) * denom;
} else {
arr1v(dGamma3, i) = 0.0;
arr1v(dgamma3, i) = 0.0;
}
B = this->A_meam[elti] * this->Ec_meam[elti][elti];
if (!iszero(rhob)) {
if (this->emb_lin_neg == 1 && rhob <= 0) {
arr1v(fp, i) = -B;
arr1v(frhop, i) = -B;
} else {
arr1v(fp, i) = B * (log(rhob) + 1.0);
arr1v(frhop, i) = B * (log(rhob) + 1.0);
}
if (*eflag_either != 0) {
if (*eflag_global != 0) {
@ -186,9 +173,9 @@ MEAM::meam_dens_final(int* nlocal, int* nmax, int* eflag_either, int* eflag_glob
}
} else {
if (this->emb_lin_neg == 1) {
arr1v(fp, i) = -B;
arr1v(frhop, i) = -B;
} else {
arr1v(fp, i) = B;
arr1v(frhop, i) = B;
}
}
}
@ -198,38 +185,38 @@ MEAM::meam_dens_final(int* nlocal, int* nmax, int* eflag_either, int* eflag_glob
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* errorflag)
MEAM::G_gam(double gamma, int ibar, double gsmooth_factor, double* G, int* errorflag)
{
// Compute G(Gamma) based on selection flag ibar:
// 0 => G = sqrt(1+Gamma)
// 1 => G = exp(Gamma/2)
// Compute G(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+gamma)
// 1 => G = exp(gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+exp(-Gamma))
// 4 => G = sqrt(1+Gamma)
// -5 => G = +-sqrt(abs(1+Gamma))
// 3 => G = 2/(1+exp(-gamma))
// 4 => G = sqrt(1+gamma)
// -5 => G = +-sqrt(abs(1+gamma))
double gsmooth_switchpoint;
if (ibar == 0 || ibar == 4) {
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);
if (Gamma < gsmooth_switchpoint) {
if (gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/Gamma)**99
// G = 0.01*(-0.99/gamma)**99
*G = 1 / (gsmooth_factor + 1) *
pow((gsmooth_switchpoint / Gamma), gsmooth_factor);
pow((gsmooth_switchpoint / gamma), gsmooth_factor);
*G = sqrt(*G);
} else {
*G = sqrt(1.0 + Gamma);
*G = sqrt(1.0 + gamma);
}
} else if (ibar == 1) {
*G = MathSpecial::fm_exp(Gamma / 2.0);
*G = MathSpecial::fm_exp(gamma / 2.0);
} else if (ibar == 3) {
*G = 2.0 / (1.0 + exp(-Gamma));
*G = 2.0 / (1.0 + exp(-gamma));
} else if (ibar == -5) {
if ((1.0 + Gamma) >= 0) {
*G = sqrt(1.0 + Gamma);
if ((1.0 + gamma) >= 0) {
*G = sqrt(1.0 + gamma);
} else {
*G = -sqrt(-1.0 - Gamma);
*G = -sqrt(-1.0 - gamma);
}
} else {
*errorflag = 1;
@ -239,43 +226,43 @@ MEAM::G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* error
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* dG)
MEAM::dG_gam(double gamma, int ibar, double gsmooth_factor, double* G, double* dG)
{
// Compute G(Gamma) and dG(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+Gamma)
// 1 => G = MathSpecial::fm_exp(Gamma/2)
// Compute G(gamma) and dG(gamma) based on selection flag ibar:
// 0 => G = sqrt(1+gamma)
// 1 => G = MathSpecial::fm_exp(gamma/2)
// 2 => not implemented
// 3 => G = 2/(1+MathSpecial::fm_exp(-Gamma))
// 4 => G = sqrt(1+Gamma)
// -5 => G = +-sqrt(abs(1+Gamma))
// 3 => G = 2/(1+MathSpecial::fm_exp(-gamma))
// 4 => G = sqrt(1+gamma)
// -5 => G = +-sqrt(abs(1+gamma))
double gsmooth_switchpoint;
if (ibar == 0 || ibar == 4) {
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);
if (Gamma < gsmooth_switchpoint) {
if (gamma < gsmooth_switchpoint) {
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/Gamma)**99
// G = 0.01*(-0.99/gamma)**99
*G = 1 / (gsmooth_factor + 1) *
pow((gsmooth_switchpoint / Gamma), gsmooth_factor);
pow((gsmooth_switchpoint / gamma), gsmooth_factor);
*G = sqrt(*G);
*dG = -gsmooth_factor * *G / (2.0 * Gamma);
*dG = -gsmooth_factor * *G / (2.0 * gamma);
} else {
*G = sqrt(1.0 + Gamma);
*G = sqrt(1.0 + gamma);
*dG = 1.0 / (2.0 * *G);
}
} else if (ibar == 1) {
*G = MathSpecial::fm_exp(Gamma / 2.0);
*G = MathSpecial::fm_exp(gamma / 2.0);
*dG = *G / 2.0;
} else if (ibar == 3) {
*G = 2.0 / (1.0 + MathSpecial::fm_exp(-Gamma));
*G = 2.0 / (1.0 + MathSpecial::fm_exp(-gamma));
*dG = *G * (2.0 - *G) / 2;
} else if (ibar == -5) {
if ((1.0 + Gamma) >= 0) {
*G = sqrt(1.0 + Gamma);
if ((1.0 + gamma) >= 0) {
*G = sqrt(1.0 + gamma);
*dG = 1.0 / (2.0 * *G);
} else {
*G = -sqrt(-1.0 - Gamma);
*G = -sqrt(-1.0 - gamma);
*dG = -1.0 / (2.0 * *G);
}
}

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

@ -3,6 +3,85 @@
#include "math_special.h"
using namespace LAMMPS_NS;
void
MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
{
int i, j;
// grow local arrays if necessary
if (atom_nmax > nmax) {
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
nmax = atom_nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(rho0,nmax,"pair:rho0");
memory->create(rho1,nmax,"pair:rho1");
memory->create(rho2,nmax,"pair:rho2");
memory->create(rho3,nmax,"pair:rho3");
memory->create(frhop,nmax,"pair:frhop");
memory->create(gamma,nmax,"pair:gamma");
memory->create(dgamma1,nmax,"pair:dgamma1");
memory->create(dgamma2,nmax,"pair:dgamma2");
memory->create(dgamma3,nmax,"pair:dgamma3");
memory->create(arho2b,nmax,"pair:arho2b");
memory->create(arho1,nmax,3,"pair:arho1");
memory->create(arho2,nmax,6,"pair:arho2");
memory->create(arho3,nmax,10,"pair:arho3");
memory->create(arho3b,nmax,3,"pair:arho3b");
memory->create(t_ave,nmax,3,"pair:t_ave");
memory->create(tsq_ave,nmax,3,"pair:tsq_ave");
}
if (n_neigh > maxneigh) {
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
maxneigh = n_neigh;
memory->create(scrfcn,maxneigh,"pair:scrfcn");
memory->create(dscrfcn,maxneigh,"pair:dscrfcn");
memory->create(fcpair,maxneigh,"pair:fcpair");
}
// zero out local arrays
for (i = 0; i < nall; i++) {
rho0[i] = 0.0;
arho2b[i] = 0.0;
arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
for (j = 0; j < 6; j++) arho2[i][j] = 0.0;
for (j = 0; j < 10; j++) arho3[i][j] = 0.0;
arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
}
}
// Extern "C" declaration has the form:
//
// void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *,
@ -22,34 +101,27 @@ using namespace LAMMPS_NS;
//
void
MEAM::meam_dens_init(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x,
MEAM::meam_dens_init(int* i, int* ntype, int* type, int* fmap, double** x,
int* numneigh, int* firstneigh, int* numneigh_full,
int* firstneigh_full, double* scrfcn, double* dscrfcn,
double* fcpair, double* rho0, double* arho1, double* arho2,
double* arho2b, double* arho3, double* arho3b, double* t_ave,
double* tsq_ave, int* errorflag)
int* firstneigh_full, int fnoffset, int* errorflag)
{
*errorflag = 0;
// Compute screening function and derivatives
getscreen(*i, *nmax, scrfcn, dscrfcn, fcpair, x, *numneigh, firstneigh,
getscreen(*i, &scrfcn[fnoffset], &dscrfcn[fnoffset], &fcpair[fnoffset], x, *numneigh, firstneigh,
*numneigh_full, firstneigh_full, *ntype, type, fmap);
// Calculate intermediate density terms to be communicated
calc_rho1(*i, *nmax, *ntype, type, fmap, x, *numneigh, firstneigh, scrfcn,
fcpair, rho0, arho1, arho2, arho2b, arho3, arho3b, t_ave, tsq_ave);
calc_rho1(*i, *ntype, type, fmap, x, *numneigh, firstneigh, &scrfcn[fnoffset], &fcpair[fnoffset]);
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair,
double* x, int numneigh, int* firstneigh, int numneigh_full,
MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair,
double** x, int numneigh, int* firstneigh, int numneigh_full,
int* firstneigh_full, int ntype, int* type, int* fmap)
{
arrdim2v(x, 3, nmax);
int jn, j, kn, k;
int elti, eltj, eltk;
double xitmp, yitmp, zitmp, delxij, delyij, delzij, rij2, rij;
@ -92,8 +164,7 @@ MEAM::getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair
sij = 0.0;
} else {
rnorm = (this->rc_meam - rij) * drinv;
screen(i, j, nmax, x, rij2, &sij, numneigh_full, firstneigh_full,
ntype, type, fmap);
screen(i, j, x, rij2, &sij, numneigh_full, firstneigh_full, ntype, type, fmap);
dfcut(rnorm, &fc, &dfc);
fcij = fc;
dfcij = dfc * drinv;
@ -170,19 +241,9 @@ MEAM::getscreen(int i, int nmax, double* scrfcn, double* dscrfcn, double* fcpair
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x,
int numneigh, int* firstneigh, double* scrfcn, double* fcpair,
double* rho0, double* arho1, double* arho2, double* arho2b,
double* arho3, double* arho3b, double* t_ave, double* tsq_ave)
MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x,
int numneigh, int* firstneigh, double* scrfcn, double* fcpair)
{
arrdim2v(x, 3, nmax);
arrdim2v(arho1, 3, nmax);
arrdim2v(arho2, 6, nmax);
arrdim2v(arho3, 10, nmax);
arrdim2v(arho3b, 3, nmax);
arrdim2v(t_ave, 3, nmax);
arrdim2v(tsq_ave, 3, nmax);
int jn, j, m, n, p, elti, eltj;
int nv2, nv3;
double xtmp, ytmp, ztmp, delij[3 + 1], rij2, rij, sij;
@ -302,7 +363,7 @@ MEAM::calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x,
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::screen(int i, int j, int nmax, double* x, double rijsq, double* sij,
MEAM::screen(int i, int j, double** x, double rijsq, double* sij,
int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap)
// Screening function
// Inputs: i = atom 1 id (integer)
@ -310,11 +371,7 @@ MEAM::screen(int i, int j, int nmax, double* x, double rijsq, double* sij,
// rijsq = squared distance between i and j
// Outputs: sij = screening function
{
arrdim2v(x, 3, nmax)
int k,
nk /*,m*/;
int k, nk /*,m*/;
int elti, eltj, eltk;
double delxik, delyik, delzik;
double delxjk, delyjk, delzjk;
@ -375,8 +432,8 @@ MEAM::screen(int i, int j, int nmax, double* x, double rijsq, double* sij,
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2,
double* dsij1, double* dsij2, int ntype, int* type, int* fmap, double* x,
MEAM::dsij(int i, int j, int k, int jn, int numneigh, double rij2,
double* dsij1, double* dsij2, int ntype, int* type, int* fmap, double** x,
double* scrfcn, double* fcpair)
{
// Inputs: i,j,k = id's of 3 atom triplet
@ -385,7 +442,7 @@ MEAM::dsij(int i, int j, int k, int jn, int nmax, int numneigh, double rij2,
// Outputs: dsij1 = deriv. of sij w.r.t. rik
// dsij2 = deriv. of sij w.r.t. rjk
arrdim2v(x, 3, nmax) int elti, eltj, eltk;
int elti, eltj, eltk;
double rik2, rjk2;
double dxik, dyik, dzik;

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

@ -26,27 +26,13 @@ using namespace LAMMPS_NS;
//
void
MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
MEAM::meam_force(int* iptr, 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* ntype, int* type, int* fmap, double** x, int* numneigh,
int* firstneigh, int* numneigh_full, int* firstneigh_full,
double* scrfcn, double* dscrfcn, double* fcpair, double* dGamma1,
double* dGamma2, double* dGamma3, double* rho0, double* rho1,
double* rho2, double* rho3, double* fp, double* Arho1,
double* Arho2, double* Arho2b, double* Arho3, double* Arho3b,
double* t_ave, double* tsq_ave, double* f, double* vatom,
int fnoffset, double** f, double** vatom,
int* errorflag)
{
arrdim2v(x, 3, *nmax);
arrdim2v(Arho1, 3, *nmax);
arrdim2v(Arho2, 6, *nmax);
arrdim2v(Arho3, 10, *nmax);
arrdim2v(Arho3b, 3, *nmax);
arrdim2v(t_ave, 3, *nmax);
arrdim2v(tsq_ave, 3, *nmax);
arrdim2v(f, 3, *nmax);
arrdim2v(vatom, 6, *nmax);
int i, j, jn, k, kn, kk, m, n, p, q;
int nv2, nv3, elti, eltj, eltk, ind;
double xitmp, yitmp, zitmp, delij[3 + 1], rij2, rij, rij3;
@ -100,9 +86,9 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
j = arr1v(firstneigh, jn);
eltj = arr1v(fmap, arr1v(type, j));
if (!iszero(arr1v(scrfcn, jn)) && eltj > 0) {
if (!iszero(arr1v(scrfcn, fnoffset + jn)) && eltj > 0) {
sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn);
sij = arr1v(scrfcn, fnoffset + jn) * arr1v(fcpair, fnoffset + jn);
delij[1] = arr2v(x, 1, j) - xitmp;
delij[2] = arr2v(x, 2, j) - yitmp;
delij[3] = arr2v(x, 3, j) - zitmp;
@ -207,19 +193,19 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
for (p = n; p <= 3; p++) {
for (q = p; q <= 3; q++) {
arg = delij[n] * delij[p] * delij[q] * this->v3D[nv3];
arg1i3 = arg1i3 + arr2v(Arho3, nv3, i) * arg;
arg1j3 = arg1j3 - arr2v(Arho3, nv3, j) * arg;
arg1i3 = arg1i3 + arr2v(arho3, nv3, i) * arg;
arg1j3 = arg1j3 - arr2v(arho3, nv3, j) * arg;
nv3 = nv3 + 1;
}
arg = delij[n] * delij[p] * this->v2D[nv2];
arg1i2 = arg1i2 + arr2v(Arho2, nv2, i) * arg;
arg1j2 = arg1j2 + arr2v(Arho2, nv2, j) * arg;
arg1i2 = arg1i2 + arr2v(arho2, nv2, i) * arg;
arg1j2 = arg1j2 + arr2v(arho2, nv2, j) * arg;
nv2 = nv2 + 1;
}
arg1i1 = arg1i1 + arr2v(Arho1, n, i) * delij[n];
arg1j1 = arg1j1 - arr2v(Arho1, n, j) * delij[n];
arg3i3 = arg3i3 + arr2v(Arho3b, n, i) * delij[n];
arg3j3 = arg3j3 - arr2v(Arho3b, n, j) * delij[n];
arg1i1 = arg1i1 + arr2v(arho1, n, i) * delij[n];
arg1j1 = arg1j1 - arr2v(arho1, n, j) * delij[n];
arg3i3 = arg3i3 + arr2v(arho3b, n, i) * delij[n];
arg3j3 = arg3j3 - arr2v(arho3b, n, j) * delij[n];
}
// rho0 terms
@ -232,25 +218,25 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1;
a1 = 2.0 * sij / rij;
for (m = 1; m <= 3; m++) {
drho1drm1[m] = a1 * rhoa1j * arr2v(Arho1, m, i);
drho1drm2[m] = -a1 * rhoa1i * arr2v(Arho1, m, j);
drho1drm1[m] = a1 * rhoa1j * arr2v(arho1, m, i);
drho1drm2[m] = -a1 * rhoa1i * arr2v(arho1, m, j);
}
// rho2 terms
a2 = 2 * sij / rij2;
drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 -
2.0 / 3.0 * arr1v(Arho2b, i) * drhoa2j * sij;
2.0 / 3.0 * arr1v(arho2b, i) * drhoa2j * sij;
drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 -
2.0 / 3.0 * arr1v(Arho2b, j) * drhoa2i * sij;
2.0 / 3.0 * arr1v(arho2b, j) * drhoa2i * sij;
a2 = 4 * sij / rij2;
for (m = 1; m <= 3; m++) {
drho2drm1[m] = 0.0;
drho2drm2[m] = 0.0;
for (n = 1; n <= 3; n++) {
drho2drm1[m] = drho2drm1[m] +
arr2v(Arho2, this->vind2D[m][n], i) * delij[n];
arr2v(arho2, this->vind2D[m][n], i) * delij[n];
drho2drm2[m] = drho2drm2[m] -
arr2v(Arho2, this->vind2D[m][n], j) * delij[n];
arr2v(arho2, this->vind2D[m][n], j) * delij[n];
}
drho2drm1[m] = a2 * rhoa2j * drho2drm1[m];
drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m];
@ -274,16 +260,16 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
for (p = n; p <= 3; p++) {
arg = delij[n] * delij[p] * this->v2D[nv2];
drho3drm1[m] = drho3drm1[m] +
arr2v(Arho3, this->vind3D[m][n][p], i) * arg;
arr2v(arho3, this->vind3D[m][n][p], i) * arg;
drho3drm2[m] = drho3drm2[m] +
arr2v(Arho3, this->vind3D[m][n][p], j) * arg;
arr2v(arho3, this->vind3D[m][n][p], j) * arg;
nv2 = nv2 + 1;
}
}
drho3drm1[m] =
(a3 * drho3drm1[m] - a3a * arr2v(Arho3b, m, i)) * rhoa3j;
(a3 * drho3drm1[m] - a3a * arr2v(arho3b, m, i)) * rhoa3j;
drho3drm2[m] =
(-a3 * drho3drm2[m] + a3a * arr2v(Arho3b, m, j)) * rhoa3i;
(-a3 * drho3drm2[m] + a3a * arr2v(arho3b, m, j)) * rhoa3i;
}
// Compute derivatives of weighting functions t wrt rij
@ -358,32 +344,32 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
get_shpfcn(shpi, this->lattce_meam[elti][elti]);
get_shpfcn(shpj, this->lattce_meam[eltj][eltj]);
drhodr1 =
arr1v(dGamma1, i) * drho0dr1 +
arr1v(dGamma2, i) * (dt1dr1 * arr1v(rho1, i) + t1i * drho1dr1 +
arr1v(dgamma1, i) * drho0dr1 +
arr1v(dgamma2, i) * (dt1dr1 * arr1v(rho1, i) + t1i * drho1dr1 +
dt2dr1 * arr1v(rho2, i) + t2i * drho2dr1 +
dt3dr1 * arr1v(rho3, i) + t3i * drho3dr1) -
arr1v(dGamma3, i) *
arr1v(dgamma3, i) *
(shpi[1] * dt1dr1 + shpi[2] * dt2dr1 + shpi[3] * dt3dr1);
drhodr2 =
arr1v(dGamma1, j) * drho0dr2 +
arr1v(dGamma2, j) * (dt1dr2 * arr1v(rho1, j) + t1j * drho1dr2 +
arr1v(dgamma1, j) * drho0dr2 +
arr1v(dgamma2, j) * (dt1dr2 * arr1v(rho1, j) + t1j * drho1dr2 +
dt2dr2 * arr1v(rho2, j) + t2j * drho2dr2 +
dt3dr2 * arr1v(rho3, j) + t3j * drho3dr2) -
arr1v(dGamma3, j) *
arr1v(dgamma3, j) *
(shpj[1] * dt1dr2 + shpj[2] * dt2dr2 + shpj[3] * dt3dr2);
for (m = 1; m <= 3; m++) {
drhodrm1[m] = 0.0;
drhodrm2[m] = 0.0;
drhodrm1[m] =
arr1v(dGamma2, i) *
arr1v(dgamma2, i) *
(t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]);
drhodrm2[m] =
arr1v(dGamma2, j) *
arr1v(dgamma2, j) *
(t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]);
}
// Compute derivatives wrt sij, but only if necessary
if (!iszero(arr1v(dscrfcn, jn))) {
if (!iszero(arr1v(dscrfcn, fnoffset + jn))) {
drho0ds1 = rhoa0j;
drho0ds2 = rhoa0i;
a1 = 2.0 / rij;
@ -391,9 +377,9 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
drho1ds2 = a1 * rhoa1i * arg1j1;
a2 = 2.0 / rij2;
drho2ds1 =
a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arr1v(Arho2b, i) * rhoa2j;
a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arr1v(arho2b, i) * rhoa2j;
drho2ds2 =
a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arr1v(Arho2b, j) * rhoa2i;
a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arr1v(arho2b, j) * rhoa2i;
a3 = 2.0 / rij3;
a3a = 6.0 / (5.0 * rij);
drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3;
@ -460,35 +446,35 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
}
drhods1 =
arr1v(dGamma1, i) * drho0ds1 +
arr1v(dGamma2, i) * (dt1ds1 * arr1v(rho1, i) + t1i * drho1ds1 +
arr1v(dgamma1, i) * drho0ds1 +
arr1v(dgamma2, i) * (dt1ds1 * arr1v(rho1, i) + t1i * drho1ds1 +
dt2ds1 * arr1v(rho2, i) + t2i * drho2ds1 +
dt3ds1 * arr1v(rho3, i) + t3i * drho3ds1) -
arr1v(dGamma3, i) *
arr1v(dgamma3, i) *
(shpi[1] * dt1ds1 + shpi[2] * dt2ds1 + shpi[3] * dt3ds1);
drhods2 =
arr1v(dGamma1, j) * drho0ds2 +
arr1v(dGamma2, j) * (dt1ds2 * arr1v(rho1, j) + t1j * drho1ds2 +
arr1v(dgamma1, j) * drho0ds2 +
arr1v(dgamma2, j) * (dt1ds2 * arr1v(rho1, j) + t1j * drho1ds2 +
dt2ds2 * arr1v(rho2, j) + t2j * drho2ds2 +
dt3ds2 * arr1v(rho3, j) + t3j * drho3ds2) -
arr1v(dGamma3, j) *
arr1v(dgamma3, j) *
(shpj[1] * dt1ds2 + shpj[2] * dt2ds2 + shpj[3] * dt3ds2);
}
// Compute derivatives of energy wrt rij, sij and rij[3]
dUdrij = phip * sij + arr1v(fp, i) * drhodr1 + arr1v(fp, j) * drhodr2;
dUdrij = phip * sij + arr1v(frhop, i) * drhodr1 + arr1v(frhop, j) * drhodr2;
dUdsij = 0.0;
if (!iszero(arr1v(dscrfcn, jn))) {
dUdsij = phi + arr1v(fp, i) * drhods1 + arr1v(fp, j) * drhods2;
if (!iszero(arr1v(dscrfcn, fnoffset + jn))) {
dUdsij = phi + arr1v(frhop, i) * drhods1 + arr1v(frhop, j) * drhods2;
}
for (m = 1; m <= 3; m++) {
dUdrijm[m] =
arr1v(fp, i) * drhodrm1[m] + arr1v(fp, j) * drhodrm2[m];
arr1v(frhop, i) * drhodrm1[m] + arr1v(frhop, j) * drhodrm2[m];
}
// Add the part of the force due to dUdrij and dUdsij
force = dUdrij * recip + dUdsij * arr1v(dscrfcn, jn);
force = dUdrij * recip + dUdsij * arr1v(dscrfcn, fnoffset + jn);
for (m = 1; m <= 3; m++) {
forcem = delij[m] * force + dUdrijm[m];
arr2v(f, m, i) = arr2v(f, m, i) + forcem;
@ -530,8 +516,8 @@ MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global,
k = arr1v(firstneigh_full, kn);
eltk = arr1v(fmap, arr1v(type, k));
if (k != j && eltk > 0) {
dsij(i, j, k, jn, *nmax, *numneigh, rij2, &dsij1, &dsij2, *ntype,
type, fmap, x, scrfcn, fcpair);
dsij(i, j, k, jn, *numneigh, rij2, &dsij1, &dsij2, *ntype,
type, fmap, x, &scrfcn[fnoffset], &fcpair[fnoffset]);
if (!iszero(dsij1) || !iszero(dsij2)) {
force1 = dUdsij * dsij1;
force2 = dUdsij * dsij2;

64
src/USER-MEAMC/meam_impl.cpp Normal file
View File

@ -0,0 +1,64 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Sebastian Hütter (OvGU)
------------------------------------------------------------------------- */
#include "meam.h"
#include "memory.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MEAM::MEAM(Memory *mem) : memory(mem)
{
nmax = 0;
rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
maxneigh = 0;
scrfcn = dscrfcn = fcpair = NULL;
}
MEAM::~MEAM()
{
meam_cleanup();
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
}

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

@ -52,14 +52,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp)
one_coeff = 1;
manybody_flag = 1;
nmax = 0;
rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
maxneigh = 0;
allocated = 0;
scrfcn = dscrfcn = fcpair = NULL;
nelements = 0;
elements = NULL;
@ -81,29 +74,6 @@ PairMEAMC::~PairMEAMC()
{
delete meam_inst;
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] mass;
@ -128,48 +98,6 @@ void PairMEAMC::compute(int eflag, int vflag)
else evflag = vflag_fdotr = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
// grow local arrays if necessary
if (atom->nmax > nmax) {
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(rho0,nmax,"pair:rho0");
memory->create(rho1,nmax,"pair:rho1");
memory->create(rho2,nmax,"pair:rho2");
memory->create(rho3,nmax,"pair:rho3");
memory->create(frhop,nmax,"pair:frhop");
memory->create(gamma,nmax,"pair:gamma");
memory->create(dgamma1,nmax,"pair:dgamma1");
memory->create(dgamma2,nmax,"pair:dgamma2");
memory->create(dgamma3,nmax,"pair:dgamma3");
memory->create(arho2b,nmax,"pair:arho2b");
memory->create(arho1,nmax,3,"pair:arho1");
memory->create(arho2,nmax,6,"pair:arho2");
memory->create(arho3,nmax,10,"pair:arho3");
memory->create(arho3b,nmax,3,"pair:arho3b");
memory->create(t_ave,nmax,3,"pair:t_ave");
memory->create(tsq_ave,nmax,3,"pair:tsq_ave");
}
// neighbor list info
inum_half = listhalf->inum;
@ -194,29 +122,8 @@ void PairMEAMC::compute(int eflag, int vflag)
n = 0;
for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
if (n > maxneigh) {
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
maxneigh = n;
memory->create(scrfcn,maxneigh,"pair:scrfcn");
memory->create(dscrfcn,maxneigh,"pair:dscrfcn");
memory->create(fcpair,maxneigh,"pair:fcpair");
}
// zero out local arrays
for (i = 0; i < nall; i++) {
rho0[i] = 0.0;
arho2b[i] = 0.0;
arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
for (j = 0; j < 6; j++) arho2[i][j] = 0.0;
for (j = 0; j < 10; j++) arho3[i][j] = 0.0;
arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
}
meam_inst->meam_dens_setup(atom->nmax, nall, n);
double **x = atom->x;
double **f = atom->f;
@ -238,12 +145,10 @@ void PairMEAMC::compute(int eflag, int vflag)
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
ifort = i+1;
meam_inst->meam_dens_init(&ifort,&nmax,&ntype,type,fmap,&x[0][0],
meam_inst->meam_dens_init(&ifort,&ntype,type,fmap,x,
&numneigh_half[i],firstneigh_half[i],
&numneigh_full[i],firstneigh_full[i],
&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
rho0,&arho1[0][0],&arho2[0][0],arho2b,
&arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],
offset,
&errorflag);
if (errorflag) {
char str[128];
@ -255,11 +160,9 @@ void PairMEAMC::compute(int eflag, int vflag)
comm->reverse_comm_pair(this);
meam_inst->meam_dens_final(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
meam_inst->meam_dens_final(&nlocal,&eflag_either,&eflag_global,&eflag_atom,
&eng_vdwl,eatom,&ntype,type,fmap,
&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
@ -273,21 +176,19 @@ void PairMEAMC::compute(int eflag, int vflag)
// vptr is first value in vatom if it will be used by meam_force()
// else vatom may not exist, so pass dummy ptr
double *vptr;
if (vflag_atom) vptr = &vatom[0][0];
else vptr = &cutmax;
double **vptr;
if (vflag_atom) vptr = vatom;
else vptr = NULL;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
ifort = i+1;
meam_inst->meam_force(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom,
&vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
meam_inst->meam_force(&ifort,&eflag_either,&eflag_global,&eflag_atom,
&vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,x,
&numneigh_half[i],firstneigh_half[i],
&numneigh_full[i],firstneigh_full[i],
&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
&t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag);
offset,
f,vptr,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
@ -738,35 +639,35 @@ int PairMEAMC::pack_forward_comm(int n, int *list, double *buf,
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = rho0[j];
buf[m++] = rho1[j];
buf[m++] = rho2[j];
buf[m++] = rho3[j];
buf[m++] = frhop[j];
buf[m++] = gamma[j];
buf[m++] = dgamma1[j];
buf[m++] = dgamma2[j];
buf[m++] = dgamma3[j];
buf[m++] = arho2b[j];
buf[m++] = arho1[j][0];
buf[m++] = arho1[j][1];
buf[m++] = arho1[j][2];
buf[m++] = arho2[j][0];
buf[m++] = arho2[j][1];
buf[m++] = arho2[j][2];
buf[m++] = arho2[j][3];
buf[m++] = arho2[j][4];
buf[m++] = arho2[j][5];
for (k = 0; k < 10; k++) buf[m++] = arho3[j][k];
buf[m++] = arho3b[j][0];
buf[m++] = arho3b[j][1];
buf[m++] = arho3b[j][2];
buf[m++] = t_ave[j][0];
buf[m++] = t_ave[j][1];
buf[m++] = t_ave[j][2];
buf[m++] = tsq_ave[j][0];
buf[m++] = tsq_ave[j][1];
buf[m++] = tsq_ave[j][2];
buf[m++] = meam_inst->rho0[j];
buf[m++] = meam_inst->rho1[j];
buf[m++] = meam_inst->rho2[j];
buf[m++] = meam_inst->rho3[j];
buf[m++] = meam_inst->frhop[j];
buf[m++] = meam_inst->gamma[j];
buf[m++] = meam_inst->dgamma1[j];
buf[m++] = meam_inst->dgamma2[j];
buf[m++] = meam_inst->dgamma3[j];
buf[m++] = meam_inst->arho2b[j];
buf[m++] = meam_inst->arho1[j][0];
buf[m++] = meam_inst->arho1[j][1];
buf[m++] = meam_inst->arho1[j][2];
buf[m++] = meam_inst->arho2[j][0];
buf[m++] = meam_inst->arho2[j][1];
buf[m++] = meam_inst->arho2[j][2];
buf[m++] = meam_inst->arho2[j][3];
buf[m++] = meam_inst->arho2[j][4];
buf[m++] = meam_inst->arho2[j][5];
for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[j][k];
buf[m++] = meam_inst->arho3b[j][0];
buf[m++] = meam_inst->arho3b[j][1];
buf[m++] = meam_inst->arho3b[j][2];
buf[m++] = meam_inst->t_ave[j][0];
buf[m++] = meam_inst->t_ave[j][1];
buf[m++] = meam_inst->t_ave[j][2];
buf[m++] = meam_inst->tsq_ave[j][0];
buf[m++] = meam_inst->tsq_ave[j][1];
buf[m++] = meam_inst->tsq_ave[j][2];
}
return m;
@ -781,35 +682,35 @@ void PairMEAMC::unpack_forward_comm(int n, int first, double *buf)
m = 0;
last = first + n;
for (i = first; i < last; i++) {
rho0[i] = buf[m++];
rho1[i] = buf[m++];
rho2[i] = buf[m++];
rho3[i] = buf[m++];
frhop[i] = buf[m++];
gamma[i] = buf[m++];
dgamma1[i] = buf[m++];
dgamma2[i] = buf[m++];
dgamma3[i] = buf[m++];
arho2b[i] = buf[m++];
arho1[i][0] = buf[m++];
arho1[i][1] = buf[m++];
arho1[i][2] = buf[m++];
arho2[i][0] = buf[m++];
arho2[i][1] = buf[m++];
arho2[i][2] = buf[m++];
arho2[i][3] = buf[m++];
arho2[i][4] = buf[m++];
arho2[i][5] = buf[m++];
for (k = 0; k < 10; k++) arho3[i][k] = buf[m++];
arho3b[i][0] = buf[m++];
arho3b[i][1] = buf[m++];
arho3b[i][2] = buf[m++];
t_ave[i][0] = buf[m++];
t_ave[i][1] = buf[m++];
t_ave[i][2] = buf[m++];
tsq_ave[i][0] = buf[m++];
tsq_ave[i][1] = buf[m++];
tsq_ave[i][2] = buf[m++];
meam_inst->rho0[i] = buf[m++];
meam_inst->rho1[i] = buf[m++];
meam_inst->rho2[i] = buf[m++];
meam_inst->rho3[i] = buf[m++];
meam_inst->frhop[i] = buf[m++];
meam_inst->gamma[i] = buf[m++];
meam_inst->dgamma1[i] = buf[m++];
meam_inst->dgamma2[i] = buf[m++];
meam_inst->dgamma3[i] = buf[m++];
meam_inst->arho2b[i] = buf[m++];
meam_inst->arho1[i][0] = buf[m++];
meam_inst->arho1[i][1] = buf[m++];
meam_inst->arho1[i][2] = buf[m++];
meam_inst->arho2[i][0] = buf[m++];
meam_inst->arho2[i][1] = buf[m++];
meam_inst->arho2[i][2] = buf[m++];
meam_inst->arho2[i][3] = buf[m++];
meam_inst->arho2[i][4] = buf[m++];
meam_inst->arho2[i][5] = buf[m++];
for (k = 0; k < 10; k++) meam_inst->arho3[i][k] = buf[m++];
meam_inst->arho3b[i][0] = buf[m++];
meam_inst->arho3b[i][1] = buf[m++];
meam_inst->arho3b[i][2] = buf[m++];
meam_inst->t_ave[i][0] = buf[m++];
meam_inst->t_ave[i][1] = buf[m++];
meam_inst->t_ave[i][2] = buf[m++];
meam_inst->tsq_ave[i][0] = buf[m++];
meam_inst->tsq_ave[i][1] = buf[m++];
meam_inst->tsq_ave[i][2] = buf[m++];
}
}
@ -822,27 +723,27 @@ int PairMEAMC::pack_reverse_comm(int n, int first, double *buf)
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = rho0[i];
buf[m++] = arho2b[i];
buf[m++] = arho1[i][0];
buf[m++] = arho1[i][1];
buf[m++] = arho1[i][2];
buf[m++] = arho2[i][0];
buf[m++] = arho2[i][1];
buf[m++] = arho2[i][2];
buf[m++] = arho2[i][3];
buf[m++] = arho2[i][4];
buf[m++] = arho2[i][5];
for (k = 0; k < 10; k++) buf[m++] = arho3[i][k];
buf[m++] = arho3b[i][0];
buf[m++] = arho3b[i][1];
buf[m++] = arho3b[i][2];
buf[m++] = t_ave[i][0];
buf[m++] = t_ave[i][1];
buf[m++] = t_ave[i][2];
buf[m++] = tsq_ave[i][0];
buf[m++] = tsq_ave[i][1];
buf[m++] = tsq_ave[i][2];
buf[m++] = meam_inst->rho0[i];
buf[m++] = meam_inst->arho2b[i];
buf[m++] = meam_inst->arho1[i][0];
buf[m++] = meam_inst->arho1[i][1];
buf[m++] = meam_inst->arho1[i][2];
buf[m++] = meam_inst->arho2[i][0];
buf[m++] = meam_inst->arho2[i][1];
buf[m++] = meam_inst->arho2[i][2];
buf[m++] = meam_inst->arho2[i][3];
buf[m++] = meam_inst->arho2[i][4];
buf[m++] = meam_inst->arho2[i][5];
for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[i][k];
buf[m++] = meam_inst->arho3b[i][0];
buf[m++] = meam_inst->arho3b[i][1];
buf[m++] = meam_inst->arho3b[i][2];
buf[m++] = meam_inst->t_ave[i][0];
buf[m++] = meam_inst->t_ave[i][1];
buf[m++] = meam_inst->t_ave[i][2];
buf[m++] = meam_inst->tsq_ave[i][0];
buf[m++] = meam_inst->tsq_ave[i][1];
buf[m++] = meam_inst->tsq_ave[i][2];
}
return m;
@ -857,27 +758,27 @@ void PairMEAMC::unpack_reverse_comm(int n, int *list, double *buf)
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
rho0[j] += buf[m++];
arho2b[j] += buf[m++];
arho1[j][0] += buf[m++];
arho1[j][1] += buf[m++];
arho1[j][2] += buf[m++];
arho2[j][0] += buf[m++];
arho2[j][1] += buf[m++];
arho2[j][2] += buf[m++];
arho2[j][3] += buf[m++];
arho2[j][4] += buf[m++];
arho2[j][5] += buf[m++];
for (k = 0; k < 10; k++) arho3[j][k] += buf[m++];
arho3b[j][0] += buf[m++];
arho3b[j][1] += buf[m++];
arho3b[j][2] += buf[m++];
t_ave[j][0] += buf[m++];
t_ave[j][1] += buf[m++];
t_ave[j][2] += buf[m++];
tsq_ave[j][0] += buf[m++];
tsq_ave[j][1] += buf[m++];
tsq_ave[j][2] += buf[m++];
meam_inst->rho0[j] += buf[m++];
meam_inst->arho2b[j] += buf[m++];
meam_inst->arho1[j][0] += buf[m++];
meam_inst->arho1[j][1] += buf[m++];
meam_inst->arho1[j][2] += buf[m++];
meam_inst->arho2[j][0] += buf[m++];
meam_inst->arho2[j][1] += buf[m++];
meam_inst->arho2[j][2] += buf[m++];
meam_inst->arho2[j][3] += buf[m++];
meam_inst->arho2[j][4] += buf[m++];
meam_inst->arho2[j][5] += buf[m++];
for (k = 0; k < 10; k++) meam_inst->arho3[j][k] += buf[m++];
meam_inst->arho3b[j][0] += buf[m++];
meam_inst->arho3b[j][1] += buf[m++];
meam_inst->arho3b[j][2] += buf[m++];
meam_inst->t_ave[j][0] += buf[m++];
meam_inst->t_ave[j][1] += buf[m++];
meam_inst->t_ave[j][2] += buf[m++];
meam_inst->tsq_ave[j][0] += buf[m++];
meam_inst->tsq_ave[j][1] += buf[m++];
meam_inst->tsq_ave[j][2] += buf[m++];
}
}
@ -887,9 +788,9 @@ void PairMEAMC::unpack_reverse_comm(int n, int *list, double *buf)
double PairMEAMC::memory_usage()
{
double bytes = 11 * nmax * sizeof(double);
bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double);
bytes += 3 * maxneigh * sizeof(double);
double bytes = 11 * meam_inst->nmax * sizeof(double);
bytes += (3 + 6 + 10 + 3 + 3 + 3) * meam_inst->nmax * sizeof(double);
bytes += 3 * meam_inst->maxneigh * sizeof(double);
return bytes;
}

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

@ -52,14 +52,6 @@ class PairMEAMC : public Pair {
int *map; // mapping from atom types to elements
int *fmap; // Fortran version of map array for MEAM lib
int maxneigh;
double *scrfcn,*dscrfcn,*fcpair;
int nmax;
double *rho,*rho0,*rho1,*rho2,*rho3,*frhop;
double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b;
double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave;
void allocate();
void read_files(char *, char *);
void neigh_strip(int, int *, int *, int **);