tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Re-Run clang-format

Browse Source
This commit is contained in:
Sebastian Hütter
2017-06-30 12:26:12 +02:00
parent 8fca667e4b
commit e0939ac795
7 changed files with 259 additions and 428 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,8 +1,8 @@
#ifndef LMP_MEAM_H
#define LMP_MEAM_H
#include <stdlib.h>
#include "memory.h"
#include <stdlib.h>
#define maxelt 5
@ -10,12 +10,14 @@ namespace LAMMPS_NS {
typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t;
class MEAM {
public:
MEAM(Memory *mem);
class MEAM
{
public:
MEAM(Memory* mem);
~MEAM();
private:
Memory *memory;
private:
Memory* memory;
// cutforce = force cutoff
// cutforcesq = force cutoff squared
@ -64,8 +66,7 @@ class MEAM {
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 A_meam[maxelt], alpha_meam[maxelt][maxelt], rho0_meam[maxelt];
double delta_meam[maxelt][maxelt];
double beta0_meam[maxelt], beta1_meam[maxelt];
double beta2_meam[maxelt], beta3_meam[maxelt];
@ -79,13 +80,11 @@ class MEAM {
int eltind[maxelt][maxelt];
int neltypes;
double **phir;
double** phir;
double **phirar, **phirar1, **phirar2, **phirar3, **phirar4, **phirar5,
**phirar6;
double **phirar, **phirar1, **phirar2, **phirar3, **phirar4, **phirar5, **phirar6;
double attrac_meam[maxelt][maxelt],
repuls_meam[maxelt][maxelt];
double attrac_meam[maxelt][maxelt], repuls_meam[maxelt][maxelt];
double Cmin_meam[maxelt][maxelt][maxelt];
double Cmax_meam[maxelt][maxelt][maxelt];
@ -101,26 +100,25 @@ class MEAM {
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;
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:
double *scrfcn, *dscrfcn, *fcpair;
protected:
void meam_checkindex(int, int, int, int*, int*);
void G_gam(double, int, double*, int*);
void dG_gam(double, int, 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 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*);
@ -131,7 +129,8 @@ public:
double phi_meam(double, int, int);
void compute_reference_density();
void get_shpfcn(double*, lattice_t);
void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double, double, double, double, int, int, lattice_t);
void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double,
double, double, double, int, int, lattice_t);
void get_Zij(int*, lattice_t);
void get_Zij2(int*, double*, double*, lattice_t, double, double);
void get_sijk(double, int, int, int, double*);
@ -140,46 +139,41 @@ public:
double erose(double, double, double, double, double, double, int);
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, 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);
public:
void meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha,
double* b0, double* b1, double* b2, double* b3, double* alat, double* esub,
double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero,
int* ibar);
void meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag);
void meam_setup_done(double* cutmax);
void meam_dens_setup(int atom_nmax, int nall, int n_neigh);
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 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, int numneigh_full, int* firstneigh_full,
int fnoffset, double** f, double** vatom,
int* errorflag);
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 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,
int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom,
int* errorflag);
};
// Functions we need for compat
#define iszero(f) (fabs(f) < 1e-20)
#define setall2d(arr, v) \
{ \
for (int __i = 0; __i < maxelt; __i++) \
for (int __j = 0; __j < maxelt; __j++) \
arr[__i][__j] = v; \
#define setall2d(arr, v) \
{ \
for (int __i = 0; __i < maxelt; __i++) \
for (int __j = 0; __j < maxelt; __j++) \
arr[__i][__j] = v; \
}
#define setall3d(arr, v) \
{ \
for (int __i = 0; __i < maxelt; __i++) \
for (int __j = 0; __j < maxelt; __j++) \
for (int __k = 0; __k < maxelt; __k++) \
arr[__i][__j][__k] = v; \
#define setall3d(arr, v) \
{ \
for (int __i = 0; __i < maxelt; __i++) \
for (int __j = 0; __j < maxelt; __j++) \
for (int __k = 0; __k < maxelt; __k++) \
arr[__i][__j][__k] = v; \
}
};
#endif

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

@ -1,6 +1,6 @@
#include "meam.h"
#include <math.h>
#include "math_special.h"
#include <math.h>
using namespace LAMMPS_NS;
// Extern "C" declaration has the form:
@ -21,9 +21,8 @@ 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)
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)
{
int i, elti;
int m;
@ -43,14 +42,10 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global,
rho3[i] = rho3[i] - 3.0 / 5.0 * arho3b[i][m] * arho3b[i][m];
}
for (m = 0; m < 6; m++) {
rho2[i] =
rho2[i] +
this->v2D[m] * arho2[i][m] * arho2[i][m];
rho2[i] = rho2[i] + this->v2D[m] * arho2[i][m] * arho2[i][m];
}
for (m = 0; m < 10; m++) {
rho3[i] =
rho3[i] +
this->v3D[m] * arho3[i][m] * arho3[i][m];
rho3[i] = rho3[i] + this->v3D[m] * arho3[i][m] * arho3[i][m];
}
if (rho0[i] > 0.0) {
@ -69,9 +64,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global,
}
}
gamma[i] = t_ave[i][0] * rho1[i] +
t_ave[i][1] * rho2[i] +
t_ave[i][2] * rho3[i];
gamma[i] = t_ave[i][0] * rho1[i] + t_ave[i][1] * rho2[i] + t_ave[i][2] * rho3[i];
if (rho0[i] > 0.0) {
gamma[i] = gamma[i] / (rho0[i] * rho0[i]);
@ -88,13 +81,9 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global,
dGbar = 0.0;
} else {
if (this->mix_ref_t == 1) {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] +
t_ave[i][2] * shp[2]) /
(Z * Z);
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
} else {
gam = (this->t1_meam[elti] * shp[0] +
this->t2_meam[elti] * shp[1] +
this->t3_meam[elti] * shp[2]) /
gam = (this->t1_meam[elti] * shp[0] + this->t2_meam[elti] * shp[1] + this->t3_meam[elti] * shp[2]) /
(Z * Z);
}
G_gam(gam, this->ibar_meam[elti], &Gbar, errorflag);
@ -106,9 +95,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global,
Gbar = 1.0;
dGbar = 0.0;
} else {
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] +
t_ave[i][2] * shp[2]) /
(Z * Z);
gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);
dG_gam(gam, this->ibar_meam[elti], &Gbar, &dGbar);
}
rho_bkgd = this->rho0_meam[elti] * Z * Gbar;
@ -196,8 +183,7 @@ MEAM::G_gam(double gamma, int ibar, double* G, int* errorflag)
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/gamma)**99
*G = 1 / (gsmooth_factor + 1) *
pow((gsmooth_switchpoint / gamma), gsmooth_factor);
*G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
*G = sqrt(*G);
} else {
*G = sqrt(1.0 + gamma);
@ -237,8 +223,7 @@ MEAM::dG_gam(double gamma, int ibar, double* G, double* dG)
// e.g. gsmooth_factor is 99, {:
// gsmooth_switchpoint = -0.99
// G = 0.01*(-0.99/gamma)**99
*G = 1 / (gsmooth_factor + 1) *
pow((gsmooth_switchpoint / gamma), gsmooth_factor);
*G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
*G = sqrt(*G);
*dG = -gsmooth_factor * *G / (2.0 * gamma);
} else {

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

@ -1,10 +1,9 @@
#include "meam.h"
#include <math.h>
#include "math_special.h"
#include <math.h>
using namespace LAMMPS_NS;
void
MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
{
@ -33,23 +32,23 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
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");
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) {
@ -57,9 +56,9 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
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");
memory->create(scrfcn, maxneigh, "pair:scrfcn");
memory->create(dscrfcn, maxneigh, "pair:dscrfcn");
memory->create(fcpair, maxneigh, "pair:fcpair");
}
// zero out local arrays
@ -68,20 +67,16 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
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;
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 *,
@ -101,9 +96,8 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
//
void
MEAM::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)
MEAM::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)
{
*errorflag = 0;
@ -118,9 +112,8 @@ MEAM::meam_dens_init(int i, int ntype, int* type, int* fmap, double** x,
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
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)
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)
{
int jn, j, kn, k;
int elti, eltj, eltk;
@ -241,8 +234,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair,
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x,
int numneigh, int* firstneigh, double* scrfcn, double* fcpair)
MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
double* scrfcn, double* fcpair)
{
int jn, j, m, n, p, elti, eltj;
int nv2, nv3;
@ -341,8 +334,8 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x,
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
MEAM::screen(int i, int j, double** x, double rijsq, double* sij,
int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap)
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)
// j = atom 2 id (integer)
@ -410,9 +403,8 @@ MEAM::screen(int i, int j, double** x, double rijsq, double* sij,
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
void
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)
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
// jn = id of i-j pair
@ -548,10 +540,8 @@ MEAM::dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, double* dCi
b = rik2 + rjk2;
denom = rij4 - a * a;
denom = denom * denom;
*dCikj1 = 4 * rij2 *
(rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom;
*dCikj2 = 4 * rij2 *
(rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom;
*dCikj1 = 4 * rij2 * (rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom;
*dCikj2 = 4 * rij2 * (rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom;
(void)(b);
}

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

@ -1,7 +1,7 @@
#include "meam.h"
#include <math.h>
#include <algorithm>
#include "math_special.h"
#include <algorithm>
#include <math.h>
using namespace LAMMPS_NS;
// Extern "C" declaration has the form:
@ -26,12 +26,10 @@ 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, int numneigh_full, int* firstneigh_full,
int fnoffset, double** f, double** vatom,
int* errorflag)
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,
int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom,
int* errorflag)
{
int j, jn, k, kn, kk, m, n, p, q;
int nv2, nv3, elti, eltj, eltk, ind;
@ -101,16 +99,9 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
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 +
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];
phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk];
recip = 1.0 / r;
if (eflag_either != 0) {
@ -221,10 +212,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
// rho2 terms
a2 = 2 * sij / rij2;
drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 -
2.0 / 3.0 * arho2b[i] * drhoa2j * sij;
drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 -
2.0 / 3.0 * arho2b[j] * drhoa2i * sij;
drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - 2.0 / 3.0 * arho2b[i] * drhoa2j * sij;
drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - 2.0 / 3.0 * arho2b[j] * drhoa2i * sij;
a2 = 4 * sij / rij2;
for (m = 0; m < 3; m++) {
drho2drm1[m] = 0.0;
@ -241,10 +230,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
rij3 = rij * rij2;
a3 = 2 * sij / rij3;
a3a = 6.0 / 5.0 * sij / rij;
drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 -
a3a * (drhoa3j - rhoa3j / rij) * arg3i3;
drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 -
a3a * (drhoa3i - rhoa3i / rij) * arg3j3;
drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 - a3a * (drhoa3j - rhoa3j / rij) * arg3i3;
drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 - a3a * (drhoa3i - rhoa3i / rij) * arg3j3;
a3 = 6 * sij / rij3;
a3a = 6 * sij / (5 * rij);
for (m = 0; m < 3; m++) {
@ -259,10 +246,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
nv2 = nv2 + 1;
}
}
drho3drm1[m] =
(a3 * drho3drm1[m] - a3a * arho3b[i][m]) * rhoa3j;
drho3drm2[m] =
(-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i;
drho3drm1[m] = (a3 * drho3drm1[m] - a3a * arho3b[i][m]) * rhoa3j;
drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i;
}
// Compute derivatives of weighting functions t wrt rij
@ -294,18 +279,12 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
if (!iszero(tsq_ave[j][2]))
a3j = drhoa0i * sij / tsq_ave[j][2];
dt1dr1 = a1i * (this->t1_meam[eltj] -
t1i * pow(this->t1_meam[eltj], 2));
dt1dr2 = a1j * (this->t1_meam[elti] -
t1j * pow(this->t1_meam[elti], 2));
dt2dr1 = a2i * (this->t2_meam[eltj] -
t2i * pow(this->t2_meam[eltj], 2));
dt2dr2 = a2j * (this->t2_meam[elti] -
t2j * pow(this->t2_meam[elti], 2));
dt3dr1 = a3i * (this->t3_meam[eltj] -
t3i * pow(this->t3_meam[eltj], 2));
dt3dr2 = a3j * (this->t3_meam[elti] -
t3j * pow(this->t3_meam[elti], 2));
dt1dr1 = a1i * (this->t1_meam[eltj] - t1i * pow(this->t1_meam[eltj], 2));
dt1dr2 = a1j * (this->t1_meam[elti] - t1j * pow(this->t1_meam[elti], 2));
dt2dr1 = a2i * (this->t2_meam[eltj] - t2i * pow(this->t2_meam[eltj], 2));
dt2dr2 = a2j * (this->t2_meam[elti] - t2j * pow(this->t2_meam[elti], 2));
dt3dr1 = a3i * (this->t3_meam[eltj] - t3i * pow(this->t3_meam[eltj], 2));
dt3dr2 = a3j * (this->t3_meam[elti] - t3j * pow(this->t3_meam[elti], 2));
} else if (this->ialloy == 2) {
@ -336,29 +315,19 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
// Compute derivatives of total density wrt rij, sij and rij(3)
get_shpfcn(shpi, this->lattce_meam[elti][elti]);
get_shpfcn(shpj, this->lattce_meam[eltj][eltj]);
drhodr1 =
dgamma1[i] * drho0dr1 +
dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 +
dt2dr1 * rho2[i] + t2i * drho2dr1 +
dt3dr1 * rho3[i] + t3i * drho3dr1) -
dgamma3[i] *
(shpi[0] * dt1dr1 + shpi[1] * dt2dr1 + shpi[2] * dt3dr1);
drhodr2 =
dgamma1[j] * drho0dr2 +
dgamma2[j] * (dt1dr2 * rho1[j] + t1j * drho1dr2 +
dt2dr2 * rho2[j] + t2j * drho2dr2 +
dt3dr2 * rho3[j] + t3j * drho3dr2) -
dgamma3[j] *
(shpj[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2);
drhodr1 = dgamma1[i] * drho0dr1 +
dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 +
dt3dr1 * rho3[i] + t3i * drho3dr1) -
dgamma3[i] * (shpi[0] * dt1dr1 + shpi[1] * dt2dr1 + shpi[2] * dt3dr1);
drhodr2 = dgamma1[j] * drho0dr2 +
dgamma2[j] * (dt1dr2 * rho1[j] + t1j * drho1dr2 + dt2dr2 * rho2[j] + t2j * drho2dr2 +
dt3dr2 * rho3[j] + t3j * drho3dr2) -
dgamma3[j] * (shpj[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2);
for (m = 0; m < 3; m++) {
drhodrm1[m] = 0.0;
drhodrm2[m] = 0.0;
drhodrm1[m] =
dgamma2[i] *
(t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]);
drhodrm2[m] =
dgamma2[j] *
(t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]);
drhodrm1[m] = dgamma2[i] * (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]);
drhodrm2[m] = dgamma2[j] * (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]);
}
// Compute derivatives wrt sij, but only if necessary
@ -369,10 +338,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
drho1ds1 = a1 * rhoa1j * arg1i1;
drho1ds2 = a1 * rhoa1i * arg1j1;
a2 = 2.0 / rij2;
drho2ds1 =
a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arho2b[i] * rhoa2j;
drho2ds2 =
a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arho2b[j] * rhoa2i;
drho2ds1 = a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arho2b[i] * rhoa2j;
drho2ds2 = a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arho2b[j] * rhoa2i;
a3 = 2.0 / rij3;
a3a = 6.0 / (5.0 * rij);
drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3;
@ -386,25 +353,25 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
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];
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 * (this->t1_meam[eltj] -
t1i * pow(this->t1_meam[eltj], 2));
dt1ds2 = a1j * (this->t1_meam[elti] -
t1j * pow(this->t1_meam[elti], 2));
dt2ds1 = a2i * (this->t2_meam[eltj] -
t2i * pow(this->t2_meam[eltj], 2));
dt2ds2 = a2j * (this->t2_meam[elti] -
t2j * pow(this->t2_meam[elti], 2));
dt3ds1 = a3i * (this->t3_meam[eltj] -
t3i * pow(this->t3_meam[eltj], 2));
dt3ds2 = a3j * (this->t3_meam[elti] -
t3j * pow(this->t3_meam[elti], 2));
dt1ds1 = a1i * (this->t1_meam[eltj] - t1i * pow(this->t1_meam[eltj], 2));
dt1ds2 = a1j * (this->t1_meam[elti] - t1j * pow(this->t1_meam[elti], 2));
dt2ds1 = a2i * (this->t2_meam[eltj] - t2i * pow(this->t2_meam[eltj], 2));
dt2ds2 = a2j * (this->t2_meam[elti] - t2j * pow(this->t2_meam[elti], 2));
dt3ds1 = a3i * (this->t3_meam[eltj] - t3i * pow(this->t3_meam[eltj], 2));
dt3ds2 = a3j * (this->t3_meam[elti] - t3j * pow(this->t3_meam[elti], 2));
} else if (this->ialloy == 2) {
@ -432,20 +399,14 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
dt3ds2 = aj * (this->t3_meam[elti] - t3j);
}
drhods1 =
dgamma1[i] * drho0ds1 +
dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 +
dt2ds1 * rho2[i] + t2i * drho2ds1 +
dt3ds1 * rho3[i] + t3i * drho3ds1) -
dgamma3[i] *
(shpi[0] * dt1ds1 + shpi[1] * dt2ds1 + shpi[2] * dt3ds1);
drhods2 =
dgamma1[j] * drho0ds2 +
dgamma2[j] * (dt1ds2 * rho1[j] + t1j * drho1ds2 +
dt2ds2 * rho2[j] + t2j * drho2ds2 +
dt3ds2 * rho3[j] + t3j * drho3ds2) -
dgamma3[j] *
(shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
drhods1 = dgamma1[i] * drho0ds1 +
dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + dt2ds1 * rho2[i] + t2i * drho2ds1 +
dt3ds1 * rho3[i] + t3i * drho3ds1) -
dgamma3[i] * (shpi[0] * dt1ds1 + shpi[1] * dt2ds1 + shpi[2] * dt3ds1);
drhods2 = dgamma1[j] * drho0ds2 +
dgamma2[j] * (dt1ds2 * rho1[j] + t1j * drho1ds2 + dt2ds2 * rho2[j] + t2j * drho2ds2 +
dt3ds2 * rho3[j] + t3j * drho3ds2) -
dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
}
// Compute derivatives of energy wrt rij, sij and rij[3]
@ -455,8 +416,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2;
}
for (m = 0; m < 3; m++) {
dUdrijm[m] =
frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];
dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];
}
// Add the part of the force due to dUdrij and dUdsij
@ -481,7 +441,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
v[4] = -0.25 * (delij[0] * fi[2] + delij[2] * fi[0]);
v[5] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]);
for (m = 0; m<6; m++) {
for (m = 0; m < 6; m++) {
vatom[i][m] = vatom[i][m] + v[m];
vatom[j][m] = vatom[j][m] + v[m];
}
@ -495,8 +455,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
k = firstneigh_full[kn];
eltk = fmap[type[k]];
if (k != j && eltk >= 0) {
dsij(i, j, k, jn, numneigh, rij2, &dsij1, &dsij2, ntype,
type, fmap, x, &scrfcn[fnoffset], &fcpair[fnoffset]);
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;
@ -522,14 +482,11 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
v[0] = -third * (delik[0] * fi[0] + deljk[0] * fj[0]);
v[1] = -third * (delik[1] * fi[1] + deljk[1] * fj[1]);
v[2] = -third * (delik[2] * fi[2] + deljk[2] * fj[2]);
v[3] = -sixth * (delik[0] * fi[1] + deljk[0] * fj[1] +
delik[1] * fi[0] + deljk[1] * fj[0]);
v[4] = -sixth * (delik[0] * fi[2] + deljk[0] * fj[2] +
delik[2] * fi[0] + deljk[2] * fj[0]);
v[5] = -sixth * (delik[1] * fi[2] + deljk[1] * fj[2] +
delik[2] * fi[1] + deljk[2] * fj[1]);
v[3] = -sixth * (delik[0] * fi[1] + deljk[0] * fj[1] + delik[1] * fi[0] + deljk[1] * fj[0]);
v[4] = -sixth * (delik[0] * fi[2] + deljk[0] * fj[2] + delik[2] * fi[0] + deljk[2] * fj[0]);
v[5] = -sixth * (delik[1] * fi[2] + deljk[1] * fj[2] + delik[2] * fi[1] + deljk[2] * fj[1]);
for (m = 0; m<6; m++) {
for (m = 0; m < 6; m++) {
vatom[i][m] = vatom[i][m] + v[m];
vatom[j][m] = vatom[j][m] + v[m];
vatom[k][m] = vatom[k][m] + v[m];

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

@ -22,7 +22,8 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MEAM::MEAM(Memory *mem) : memory(mem)
MEAM::MEAM(Memory* mem)
: memory(mem)
{
phir = phirar = phirar1 = phirar2 = phirar3 = phirar4 = phirar5 = phirar6 = NULL;
@ -69,5 +70,3 @@ MEAM::~MEAM()
memory->destroy(this->dscrfcn);
memory->destroy(this->fcpair);
}

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

@ -1,7 +1,7 @@
#include "meam.h"
#include <math.h>
#include <algorithm>
#include "math_special.h"
#include <algorithm>
#include <math.h>
using namespace LAMMPS_NS;
// Declaration in pair_meam.h:
@ -27,8 +27,7 @@ MEAM::meam_setup_done(double* cutmax)
// Augment t1 term
for (int i = 0; i < maxelt; i++)
this->t1_meam[i] =
this->t1_meam[i] + this->augt1 * 3.0 / 5.0 * this->t3_meam[i];
this->t1_meam[i] = this->t1_meam[i] + this->augt1 * 3.0 / 5.0 * this->t3_meam[i];
// Compute off-diagonal alloy parameters
alloyparams();
@ -116,28 +115,21 @@ MEAM::alloyparams(void)
if (iszero(this->Ec_meam[i][j])) {
if (this->lattce_meam[i][j] == L12)
this->Ec_meam[i][j] =
(3 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 4.0 -
this->delta_meam[i][j];
(3 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 4.0 - this->delta_meam[i][j];
else if (this->lattce_meam[i][j] == C11) {
if (this->lattce_meam[i][i] == DIA)
this->Ec_meam[i][j] =
(2 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 3.0 -
this->delta_meam[i][j];
(2 * this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
else
this->Ec_meam[i][j] =
(this->Ec_meam[i][i] + 2 * this->Ec_meam[j][j]) / 3.0 -
this->delta_meam[i][j];
(this->Ec_meam[i][i] + 2 * this->Ec_meam[j][j]) / 3.0 - this->delta_meam[i][j];
} else
this->Ec_meam[i][j] =
(this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 2.0 -
this->delta_meam[i][j];
this->Ec_meam[i][j] = (this->Ec_meam[i][i] + this->Ec_meam[j][j]) / 2.0 - this->delta_meam[i][j];
}
if (iszero(this->alpha_meam[i][j]))
this->alpha_meam[i][j] =
(this->alpha_meam[i][i] + this->alpha_meam[j][j]) / 2.0;
this->alpha_meam[i][j] = (this->alpha_meam[i][i] + this->alpha_meam[j][j]) / 2.0;
if (iszero(this->re_meam[i][j]))
this->re_meam[i][j] =
(this->re_meam[i][i] + this->re_meam[j][j]) / 2.0;
this->re_meam[i][j] = (this->re_meam[i][i] + this->re_meam[j][j]) / 2.0;
}
}
}
@ -160,8 +152,7 @@ MEAM::alloyparams(void)
for (i = 0; i < this->neltypes; i++) {
for (j = 0; j < this->neltypes; j++) {
for (k = 0; k < this->neltypes; k++) {
eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) /
(4.0 * (this->Cmax_meam[i][j][k] - 1.0));
eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) / (4.0 * (this->Cmax_meam[i][j][k] - 1.0));
this->ebound_meam[i][j] = std::max(this->ebound_meam[i][j], eb);
}
}
@ -203,44 +194,20 @@ MEAM::compute_pair_meam(void)
memory->destroy(this->phirar6);
// allocate memory for array that defines the potential
memory->create(this->phir,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phir");
memory->create(this->phir, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phir");
// allocate coeff memory
memory->create(this->phirar,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phirar");
memory->create(this->phirar1,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phirar1");
memory->create(this->phirar2,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phirar2");
memory->create(this->phirar3,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phirar3");
memory->create(this->phirar4,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phirar4");
memory->create(this->phirar5,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phirar5");
memory->create(this->phirar6,
(this->neltypes * (this->neltypes + 1)) / 2,
this->nr,
"pair:phirar6");
memory->create(this->phirar, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar");
memory->create(this->phirar1, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar1");
memory->create(this->phirar2, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar2");
memory->create(this->phirar3, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar3");
memory->create(this->phirar4, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar4");
memory->create(this->phirar5, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar5");
memory->create(this->phirar6, (this->neltypes * (this->neltypes + 1)) / 2, this->nr, "pair:phirar6");
// loop over pairs of element types
nv2 = 0;
nv2 = 0;
for (a = 0; a < this->neltypes; a++) {
for (b = a; b < this->neltypes; b++) {
// loop over r values and compute
@ -253,8 +220,8 @@ MEAM::compute_pair_meam(void)
// (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
if (this->nn2_meam[a][b] == 1) {
get_Zij(&Z1, this->lattce_meam[a][b]);
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][b],
this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b]);
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
this->Cmax_meam[a][a][b]);
// The B1, B2, and L12 cases with NN2 have a trick to them; we
// need to
@ -262,55 +229,43 @@ MEAM::compute_pair_meam(void)
// a-a
// pairs, but need to include NN2 contributions to those pairs as
// well.
if (this->lattce_meam[a][b] == B1 ||
this->lattce_meam[a][b] == B2 ||
if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||
this->lattce_meam[a][b] == L12) {
rarat = r * arat;
// phi_aa
phiaa = phi_meam(rarat, a, a);
get_Zij(&Z1, this->lattce_meam[a][a]);
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][a],
this->Cmin_meam[a][a][a],
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a]);
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);
phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a);
}
}
// phi_bb
phibb = phi_meam(rarat, b, b);
get_Zij(&Z1, this->lattce_meam[b][b]);
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[b][b],
this->Cmin_meam[b][b][b],
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[b][b], this->Cmin_meam[b][b][b],
this->Cmax_meam[b][b][b]);
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);
phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b);
}
}
if (this->lattce_meam[a][b] == B1 ||
this->lattce_meam[a][b] == B2) {
if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2) {
// Add contributions to the B1 or B2 potential
get_Zij(&Z1, this->lattce_meam[a][b]);
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][b],
this->Cmin_meam[a][a][b],
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][b], this->Cmin_meam[a][a][b],
this->Cmax_meam[a][a][b]);
this->phir[nv2][j] =
this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;
get_Zij2(&Z2, &arat, &scrn2, this->lattce_meam[a][b],
this->Cmin_meam[b][b][a],
this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;
get_Zij2(&Z2, &arat, &scrn2, this->lattce_meam[a][b], this->Cmin_meam[b][b][a],
this->Cmax_meam[b][b][a]);
this->phir[nv2][j] =
this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;
} else if (this->lattce_meam[a][b] == L12) {
// The L12 case has one last trick; we have to be careful to
@ -327,17 +282,14 @@ MEAM::compute_pair_meam(void)
get_sijk(C, b, b, a, &s221);
S11 = s111 * s111 * s112 * s112;
S22 = pow(s221, 4);
this->phir[nv2][j] = this->phir[nv2][j] -
0.75 * S11 * phiaa -
0.25 * S22 * phibb;
this->phir[nv2][j] = this->phir[nv2][j] - 0.75 * S11 * phiaa - 0.25 * S22 * phibb;
}
} 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] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b);
}
}
}
@ -349,16 +301,13 @@ MEAM::compute_pair_meam(void)
// potential is linear combination with zbl potential
// endif
if (this->zbl_meam[a][b] == 1) {
astar =
this->alpha_meam[a][b] * (r / this->re_meam[a][b] - 1.0);
astar = this->alpha_meam[a][b] * (r / this->re_meam[a][b] - 1.0);
if (astar <= -3.0)
this->phir[nv2][j] =
zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
this->phir[nv2][j] = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
else if (astar > -3.0 && astar < -1.0) {
fcut(1 - (astar + 1.0) / (-3.0 + 1.0), &frac);
phizbl = zbl(r, this->ielt_meam[a], this->ielt_meam[b]);
this->phir[nv2][j] =
frac * this->phir[nv2][j] + (1 - frac) * phizbl;
this->phir[nv2][j] = frac * this->phir[nv2][j] + (1 - frac) * phizbl;
}
}
}
@ -401,8 +350,7 @@ MEAM::phi_meam(double r, int a, int b)
// Nref[i][j] = # of i's neighbors of type j
get_Zij(&Z12, this->lattce_meam[a][b]);
get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22,
&rho32);
get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32);
// if densities are too small, numerical problems may result; just return zero
if (rho01 <= 1e-14 && rho02 <= 1e-14)
@ -419,40 +367,32 @@ MEAM::phi_meam(double r, int a, int b)
t32av = this->t3_meam[b];
} else {
scalfac = 1.0 / (rho01 + rho02);
t11av =
scalfac * (this->t1_meam[a] * rho01 + this->t1_meam[b] * rho02);
t11av = scalfac * (this->t1_meam[a] * rho01 + this->t1_meam[b] * rho02);
t12av = t11av;
t21av =
scalfac * (this->t2_meam[a] * rho01 + this->t2_meam[b] * rho02);
t21av = scalfac * (this->t2_meam[a] * rho01 + this->t2_meam[b] * rho02);
t22av = t21av;
t31av =
scalfac * (this->t3_meam[a] * rho01 + this->t3_meam[b] * rho02);
t31av = scalfac * (this->t3_meam[a] * rho01 + this->t3_meam[b] * rho02);
t32av = t31av;
}
} else {
// average weighting factors for the reference structure, eqn. I.8
get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av,
this->t1_meam[a], this->t2_meam[a], this->t3_meam[a],
this->t1_meam[b], this->t2_meam[b], this->t3_meam[b],
r, a, b, this->lattce_meam[a][b]);
get_tavref(&t11av, &t21av, &t31av, &t12av, &t22av, &t32av, this->t1_meam[a], this->t2_meam[a],
this->t3_meam[a], this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], r, a, b,
this->lattce_meam[a][b]);
}
// for c11b structure, calculate background electron densities
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 = 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 = sqrt(rhobar1);
rhobar2 = pow(Z12 * rho01, 2) + 2.0 / 3.0 * t21av * pow(rho21, 2);
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 = 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 = sqrt(rhobar2);
rhobar1 = pow(Z12 * rho02, 2) + 2.0 / 3.0 * t22av * pow(rho22, 2);
rhobar1 = sqrt(rhobar1);
@ -520,8 +460,7 @@ MEAM::phi_meam(double r, int a, int b)
if (this->emb_lin_neg == 1 && rhobar1 <= 0)
F1 = -this->A_meam[a] * this->Ec_meam[a][a] * rhobar1;
else
F1 =
this->A_meam[a] * this->Ec_meam[a][a] * rhobar1 * log(rhobar1);
F1 = this->A_meam[a] * this->Ec_meam[a][a] * rhobar1 * log(rhobar1);
}
if (iszero(rhobar2))
F2 = 0.0;
@ -529,13 +468,11 @@ MEAM::phi_meam(double r, int a, int b)
if (this->emb_lin_neg == 1 && rhobar2 <= 0)
F2 = -this->A_meam[b] * this->Ec_meam[b][b] * rhobar2;
else
F2 =
this->A_meam[b] * this->Ec_meam[b][b] * rhobar2 * log(rhobar2);
F2 = this->A_meam[b] * this->Ec_meam[b][b] * rhobar2 * log(rhobar2);
}
// compute Rose function, I.16
Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b],
this->Ec_meam[a][b], this->repuls_meam[a][b],
Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b], this->Ec_meam[a][b], this->repuls_meam[a][b],
this->attrac_meam[a][b], this->erose_form);
// calculate the pair energy
@ -552,14 +489,12 @@ MEAM::phi_meam(double r, int a, int b)
phiaa = phi_meam(r, a, a);
// account for second neighbor a-a potential here...
get_Zij(&Z1nn, this->lattce_meam[a][a]);
get_Zij2(&Z2nn, &arat, &scrn, this->lattce_meam[a][a],
this->Cmin_meam[a][a][a], this->Cmax_meam[a][a][a]);
get_Zij2(&Z2nn, &arat, &scrn, this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a]);
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);
phiaa = phiaa + pow((-Z2nn * scrn / Z1nn), n) * phi_meam(r * pow(arat, n), a, a);
}
}
phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;
@ -595,9 +530,7 @@ MEAM::compute_reference_density(void)
Gbar = 1.0;
else {
get_shpfcn(shp, this->lattce_meam[a][a]);
gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] +
this->t3_meam[a] * shp[2]) /
(Z * Z);
gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] + this->t3_meam[a] * shp[2]) / (Z * Z);
G_gam(gam, this->ibar_meam[a], &Gbar, &errorflag);
}
@ -610,10 +543,9 @@ MEAM::compute_reference_density(void)
// add on the contribution from those (accounting for partial
// screening)
if (this->nn2_meam[a][a] == 1) {
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][a],
this->Cmin_meam[a][a][a], this->Cmax_meam[a][a][a]);
rho0_2nn =
this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));
get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][a], this->Cmin_meam[a][a][a],
this->Cmax_meam[a][a][a]);
rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));
rho0 = rho0 + Z2 * rho0_2nn * scrn;
}
@ -652,10 +584,9 @@ MEAM::get_shpfcn(double* s /* s(3) */, lattice_t latt)
//------------------------------------------------------------------------------c
// Average weighting factors for the reference structure
void
MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av,
double* t22av, double* t32av, double t11, double t21, double t31,
double t12, double t22, double t32, double r, int a, int b,
lattice_t latt)
MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, double* t22av, double* t32av,
double t11, double t21, double t31, double t12, double t22, double t32, double r, int a,
int b, lattice_t latt)
{
double rhoa01, rhoa02, a1, a2, rho01 /*,rho02*/;
@ -668,8 +599,7 @@ MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av,
*t22av = t22;
*t32av = t32;
} else {
if (latt == FCC || latt == BCC || latt == DIA || latt == HCP ||
latt == B1 || latt == DIM || latt == B2) {
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;
@ -730,8 +660,7 @@ MEAM::get_Zij(int* Zij, lattice_t latt)
// neighbor screening function for lattice type "latt"
void
MEAM::get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin,
double cmax)
MEAM::get_Zij2(int* Zij2, double* a, double* S, lattice_t latt, double cmin, double cmax)
{
double /*rratio,*/ C, x, sijk;
@ -791,17 +720,15 @@ void
MEAM::get_sijk(double C, int i, int j, int k, double* sijk)
{
double x;
x = (C - this->Cmin_meam[i][j][k]) /
(this->Cmax_meam[i][j][k] - this->Cmin_meam[i][j][k]);
x = (C - this->Cmin_meam[i][j][k]) / (this->Cmax_meam[i][j][k] - this->Cmin_meam[i][j][k]);
fcut(x, sijk);
}
//------------------------------------------------------------------------------c
// Calculate density functions, assuming reference configuration
void
MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21,
double* rho31, double* rho02, double* rho12, double* rho22,
double* rho32)
MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* rho21, double* rho31,
double* rho02, double* rho12, double* rho22, double* rho32)
{
double a1, a2;
double s[3];
@ -878,11 +805,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. * 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);
if (denom > 0.)
*rho21 = *rho21 / denom * *rho01;
} else
@ -896,8 +820,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
if (this->nn2_meam[a][b] == 1) {
get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[a][a][b],
this->Cmax_meam[a][a][b]);
get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b]);
a1 = arat * r / this->re_meam[a][a] - 1.0;
a2 = arat * r / this->re_meam[b][b] - 1.0;
@ -924,8 +847,7 @@ 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
get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[b][b][a],
this->Cmax_meam[b][b][a]);
get_Zij2(&Zij2nn, &arat, &scrn, lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a]);
*rho02 = *rho02 + Zij2nn * scrn * rhoa02nn;
}
}
@ -959,8 +881,7 @@ MEAM::zbl(double r, int z1, int z2)
// Compute Rose energy function, I.16
//
double
MEAM::erose(double r, double re, double alpha, double Ec, double repuls,
double attrac, int form)
MEAM::erose(double r, double re, double alpha, double Ec, double repuls, double attrac, int form)
{
double astar, a3;
double result = 0.0;
@ -974,13 +895,11 @@ MEAM::erose(double r, double re, double alpha, double Ec, double repuls,
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) * pow(astar, 3)) * MathSpecial::fm_exp(-astar);
else if (form == 2)
result = -Ec * (1 + astar + a3 * pow(astar, 3)) * MathSpecial::fm_exp(-astar);
else
result =
-Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * MathSpecial::fm_exp(-astar);
result = -Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * MathSpecial::fm_exp(-astar);
}
return result;
}
@ -1005,37 +924,28 @@ MEAM::interpolate_meam(int ind)
}
this->phirar1[ind][0] = this->phirar[ind][1] - this->phirar[ind][0];
this->phirar1[ind][1] = 0.5 * (this->phirar[ind][2] - this->phirar[ind][0]);
this->phirar1[ind][this->nrar - 2] = 0.5 * (this->phirar[ind][this->nrar - 1] - this->phirar[ind][this->nrar - 3]);
this->phirar1[ind][this->nrar - 2] =
0.5 * (this->phirar[ind][this->nrar - 1] - this->phirar[ind][this->nrar - 3]);
this->phirar1[ind][this->nrar - 1] = 0.0;
for (j = 2; j < this->nrar - 2; j++) {
this->phirar1[ind][j] =
((this->phirar[ind][j - 2] -
this->phirar[ind][j + 2]) +
8.0 * (this->phirar[ind][j + 1] -
this->phirar[ind][j - 1])) /
12.;
this->phirar1[ind][j] = ((this->phirar[ind][j - 2] - this->phirar[ind][j + 2]) +
8.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j - 1])) /
12.;
}
for (j = 0; j < this->nrar - 1; j++) {
this->phirar2[ind][j] =
3.0 *
(this->phirar[ind][j + 1] - this->phirar[ind][j]) -
2.0 * this->phirar1[ind][j] -
this->phirar1[ind][j + 1];
this->phirar3[ind][j] =
this->phirar1[ind][j] + this->phirar1[ind][j + 1] -
2.0 *
(this->phirar[ind][j + 1] - this->phirar[ind][j]);
this->phirar2[ind][j] = 3.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]) -
2.0 * this->phirar1[ind][j] - this->phirar1[ind][j + 1];
this->phirar3[ind][j] = this->phirar1[ind][j] + this->phirar1[ind][j + 1] -
2.0 * (this->phirar[ind][j + 1] - this->phirar[ind][j]);
}
this->phirar2[ind][this->nrar - 1] = 0.0;
this->phirar3[ind][this->nrar - 1] = 0.0;
for (j = 0; j < this->nrar; j++) {
this->phirar4[ind][j] = this->phirar1[ind][j] / drar;
this->phirar5[ind][j] =
2.0 * this->phirar2[ind][j] / drar;
this->phirar6[ind][j] =
3.0 * this->phirar3[ind][j] / drar;
this->phirar5[ind][j] = 2.0 * this->phirar2[ind][j] / drar;
this->phirar6[ind][j] = 3.0 * this->phirar3[ind][j] / drar;
}
}
@ -1054,12 +964,9 @@ MEAM::compute_phi(double rij, int elti, int eltj)
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];
double result =
((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp +
this->phirar[ind][kk];
return result;
}

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

@ -18,11 +18,10 @@ using namespace LAMMPS_NS;
//
void
MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt,
double* alpha, double* b0, double* b1, double* b2,
double* b3, double* alat, double* esub, double* asub,
double* t0, double* t1, double* t2, double* t3,
double* rozero, int* ibar)
MEAM::meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, double* alpha,
double* b0, double* b1, double* b2, double* b3, double* alat, double* esub,
double* asub, double* t0, double* t1, double* t2, double* t3, double* rozero,
int* ibar)
{
int i;