mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7196 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
51aef530de
commit
cbbb66194c
|
|
@ -13,6 +13,8 @@
|
|||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Ase Henry (MIT)
|
||||
Bugfixes and optimizations:
|
||||
Marcel Fallet & Steve Stuart (Clemson), Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "math.h"
|
||||
|
|
@ -45,7 +47,6 @@ using namespace MathConst;
|
|||
PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
single_enable = 0;
|
||||
restartinfo = 0;
|
||||
one_coeff = 1;
|
||||
ghostneigh = 1;
|
||||
|
||||
|
|
@ -221,7 +222,7 @@ void PairAIREBO::init_style()
|
|||
|
||||
pgsize = neighbor->pgsize;
|
||||
oneatom = neighbor->oneatom;
|
||||
if (maxpage == 0) add_pages(0);
|
||||
if (maxpage == 0) add_pages();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -304,7 +305,7 @@ void PairAIREBO::REBO_neigh()
|
|||
int nlocal = atom->nlocal;
|
||||
int nall = nlocal + atom->nghost;
|
||||
|
||||
if (nall > maxlocal) {
|
||||
if (atom->nmax > maxlocal) {
|
||||
maxlocal = atom->nmax;
|
||||
memory->destroy(REBO_numneigh);
|
||||
memory->sfree(REBO_firstneigh);
|
||||
|
|
@ -325,7 +326,7 @@ void PairAIREBO::REBO_neigh()
|
|||
// store all REBO neighs of owned and ghost atoms
|
||||
// scan full neighbor list of I
|
||||
|
||||
npage = 0;
|
||||
int npage = 0;
|
||||
int npnt = 0;
|
||||
|
||||
for (ii = 0; ii < allnum; ii++) {
|
||||
|
|
@ -334,7 +335,7 @@ void PairAIREBO::REBO_neigh()
|
|||
if (pgsize - npnt < oneatom) {
|
||||
npnt = 0;
|
||||
npage++;
|
||||
if (npage == maxpage) add_pages(npage);
|
||||
if (npage == maxpage) add_pages();
|
||||
}
|
||||
neighptr = &pages[npage][npnt];
|
||||
n = 0;
|
||||
|
|
@ -1196,60 +1197,6 @@ void PairAIREBO::TORSION(int eflag, int vflag)
|
|||
}
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------
|
||||
// S'(t) and S(t) cutoff functions
|
||||
// ----------------------------------------------------------------------
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
cutoff function Sprime
|
||||
return cutoff and dX = derivative
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairAIREBO::Sp(double Xij, double Xmin, double Xmax, double &dX)
|
||||
{
|
||||
double cutoff;
|
||||
|
||||
double t = (Xij-Xmin) / (Xmax-Xmin);
|
||||
if (t <= 0.0) {
|
||||
cutoff = 1.0;
|
||||
dX = 0.0;
|
||||
}
|
||||
else if (t >= 1.0) {
|
||||
cutoff = 0.0;
|
||||
dX = 0.0;
|
||||
}
|
||||
else {
|
||||
cutoff = 0.5 * (1.0+cos(MY_PI*t));
|
||||
dX = (-MY_PI2*sin(MY_PI*t)) / (Xmax-Xmin);
|
||||
}
|
||||
return cutoff;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
LJ cutoff function Sp2
|
||||
return cutoff and dX = derivative
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairAIREBO::Sp2(double Xij, double Xmin, double Xmax, double &dX)
|
||||
{
|
||||
double cutoff;
|
||||
|
||||
double t = (Xij-Xmin) / (Xmax-Xmin);
|
||||
if (t <= 0.0) {
|
||||
cutoff = 1.0;
|
||||
dX = 0.0;
|
||||
}
|
||||
if (t >= 1.0) {
|
||||
cutoff = 0.0;
|
||||
dX = 0.0;
|
||||
}
|
||||
if (t>0.0 && t<1.0) {
|
||||
cutoff = (1.0-(t*t*(3.0-2.0*t)));
|
||||
dX = 6.0*(t*t-t) / (Xmax-Xmin);
|
||||
}
|
||||
return cutoff;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Bij function
|
||||
------------------------------------------------------------------------- */
|
||||
|
|
@ -2123,11 +2070,6 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag,
|
|||
NjiC = nC[atomj]-(wij*kronecker(itype,0));
|
||||
NjiH = nH[atomj]-(wij*kronecker(itype,1));
|
||||
|
||||
rij[0] = rij0[0];
|
||||
rij[1] = rij0[1];
|
||||
rij[2] = rij0[2];
|
||||
rijmag = rij0mag;
|
||||
|
||||
bij = 0.0;
|
||||
tmp = 0.0;
|
||||
tmp2 = 0.0;
|
||||
|
|
@ -3133,7 +3075,7 @@ double PairAIREBO::PijSpline(double NijC, double NijH, int typei, int typej,
|
|||
if (typei == 0 && typej == 1){
|
||||
if (NijC < pCHdom[0][0]) NijC=pCHdom[0][0];
|
||||
if (NijC > pCHdom[0][1]) NijC=pCHdom[0][1];
|
||||
if (NijH < pCHdom[0][0]) NijH=pCHdom[1][0];
|
||||
if (NijH < pCHdom[1][0]) NijH=pCHdom[1][0];
|
||||
if (NijH > pCHdom[1][1]) NijH=pCHdom[1][1];
|
||||
|
||||
if (fabs(NijC-floor(NijC)) < TOL && fabs(NijH-floor(NijH)) < TOL) {
|
||||
|
|
@ -3338,27 +3280,17 @@ double PairAIREBO::TijSpline(double Nij, double Nji,
|
|||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Kronecker delta function
|
||||
add pages to REBO neighbor list
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairAIREBO::kronecker(int a, int b)
|
||||
{
|
||||
double kd;
|
||||
if (a == b) kd = 1.0;
|
||||
else kd = 0.0;
|
||||
return kd;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
add pages to REBO neighbor list, starting at npage
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairAIREBO::add_pages(int npage)
|
||||
void PairAIREBO::add_pages()
|
||||
{
|
||||
int toppage = maxpage;
|
||||
maxpage += PGDELTA;
|
||||
|
||||
pages = (int **)
|
||||
memory->srealloc(pages,maxpage*sizeof(int *),"AIREBO:pages");
|
||||
for (int i = npage; i < maxpage; i++)
|
||||
for (int i = toppage; i < maxpage; i++)
|
||||
memory->create(pages[i],pgsize,"AIREBO:pages[i]");
|
||||
}
|
||||
|
||||
|
|
@ -3954,17 +3886,23 @@ void PairAIREBO::read_file(char *filename)
|
|||
|
||||
double PairAIREBO::Sp5th(double x, double coeffs[6], double *df)
|
||||
{
|
||||
double f;
|
||||
int i;
|
||||
i = 0;
|
||||
f = 0.0;
|
||||
*df = 0.0;
|
||||
double f, d;
|
||||
const double x2 = x*x;
|
||||
const double x3 = x2*x;
|
||||
|
||||
for (i = 0; i<6; i++) {
|
||||
f += coeffs[i]*pow(x,((double) i));
|
||||
if (i > 0) *df += coeffs[i]*((double) i)*pow(x,((double) i-1.0));
|
||||
}
|
||||
f = coeffs[0];
|
||||
f += coeffs[1]*x;
|
||||
d = coeffs[1];
|
||||
f += coeffs[2]*x2;
|
||||
d += 2.0*coeffs[2]*x;
|
||||
f += coeffs[3]*x3;
|
||||
d += 3.0*coeffs[3]*x2;
|
||||
f += coeffs[4]*x2*x2;
|
||||
d += 4.0*coeffs[4]*x3;
|
||||
f += coeffs[5]*x2*x3;
|
||||
d += 5.0*coeffs[5]*x2*x2;
|
||||
|
||||
*df = d;
|
||||
return f;
|
||||
}
|
||||
|
||||
|
|
@ -3975,25 +3913,28 @@ double PairAIREBO::Sp5th(double x, double coeffs[6], double *df)
|
|||
double PairAIREBO::Spbicubic(double x, double y,
|
||||
double coeffs[16], double df[2])
|
||||
{
|
||||
double f;
|
||||
int i,j,cn;
|
||||
double f,xn,yn,xn1,yn1,c;
|
||||
int i,j;
|
||||
|
||||
f = 0.0;
|
||||
df[0] = 0.0;
|
||||
df[1] = 0.0;
|
||||
cn = 0;
|
||||
|
||||
xn = 1.0;
|
||||
for (i = 0; i < 4; i++) {
|
||||
yn = 1.0;
|
||||
for (j = 0; j < 4; j++) {
|
||||
f += coeffs[cn]*pow(x,((double) i))*pow(y,((double) j));
|
||||
if (i > 0) df[0] +=
|
||||
(coeffs[cn]*((double) i)*pow(x,((double) i-1.0)) *
|
||||
pow(y,((double) j)));
|
||||
if (j > 0) df[1] +=
|
||||
(coeffs[cn]*((double) j)*pow(x,((double) i)) *
|
||||
pow(y,((double) j-1.0)));
|
||||
cn++;
|
||||
c = coeffs[i*4+j];
|
||||
|
||||
f += c*xn*yn;
|
||||
if (i > 0) df[0] += c * ((double) i) * xn1 * yn;
|
||||
if (j > 0) df[1] += c * ((double) j) * xn * yn1;
|
||||
|
||||
yn1 = yn;
|
||||
yn *= y;
|
||||
}
|
||||
xn1 = xn;
|
||||
xn *= x;
|
||||
}
|
||||
|
||||
return f;
|
||||
|
|
@ -4006,31 +3947,36 @@ double PairAIREBO::Spbicubic(double x, double y,
|
|||
double PairAIREBO::Sptricubic(double x, double y, double z,
|
||||
double coeffs[64], double df[3])
|
||||
{
|
||||
double f,ir,jr,kr;
|
||||
int i,j,k,cn;
|
||||
double f,ir,jr,kr,xn,yn,zn,xn1,yn1,zn1,c;
|
||||
int i,j,k;
|
||||
|
||||
f = 0.0;
|
||||
df[0] = 0.0;
|
||||
df[1] = 0.0;
|
||||
df[2] = 0.0;
|
||||
cn = 0;
|
||||
|
||||
xn = 1.0;
|
||||
for (i = 0; i < 4; i++) {
|
||||
ir = (double) i;
|
||||
yn = 1.0;
|
||||
for (j = 0; j < 4; j++) {
|
||||
jr = (double) j;
|
||||
zn = 1.0;
|
||||
for (k = 0; k < 4; k++) {
|
||||
kr = (double) k;
|
||||
f += (coeffs[cn]*pow(x,ir)*pow(y,jr)*pow(z,kr));
|
||||
if (i > 0) df[0] +=
|
||||
(coeffs[cn]*ir*pow(x,ir-1.0)*pow(y,jr)*pow(z,kr));
|
||||
if (j > 0) df[1] +=
|
||||
(coeffs[cn]*jr*pow(x,ir)*pow(y,jr-1.0)*pow(z,kr));
|
||||
if (k > 0) df[2] +=
|
||||
(coeffs[cn]*kr*pow(x,ir)*pow(y,jr)*pow(z,kr-1.0));
|
||||
cn++;
|
||||
c = coeffs[16*i+4*j+k];
|
||||
f += c*xn*yn*zn;
|
||||
if (i > 0) df[0] += c * ir * xn1 * yn * zn;
|
||||
if (j > 0) df[1] += c * jr * xn * yn1 * zn;
|
||||
if (k > 0) df[2] += c * kr * xn * yn * zn1;
|
||||
zn1 = zn;
|
||||
zn *= z;
|
||||
}
|
||||
yn1 = yn;
|
||||
yn *= y;
|
||||
}
|
||||
xn1 = xn;
|
||||
xn *= x;
|
||||
}
|
||||
|
||||
return f;
|
||||
|
|
|
|||
|
|
@ -21,6 +21,8 @@ PairStyle(airebo,PairAIREBO)
|
|||
#define LMP_PAIR_AIREBO_H
|
||||
|
||||
#include "pair.h"
|
||||
#include "math.h"
|
||||
#include "math_const.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
|
|
@ -28,7 +30,7 @@ class PairAIREBO : public Pair {
|
|||
public:
|
||||
PairAIREBO(class LAMMPS *);
|
||||
virtual ~PairAIREBO();
|
||||
void compute(int, int);
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
|
|
@ -36,15 +38,15 @@ class PairAIREBO : public Pair {
|
|||
double memory_usage();
|
||||
|
||||
protected:
|
||||
int **pages; // neighbor list pages
|
||||
int *map; // 0 (C), 1 (H), or -1 (NULL) for each type
|
||||
|
||||
int me;
|
||||
int ljflag,torflag; // 0/1 if LJ,torsion terms included
|
||||
int maxlocal; // size of numneigh, firstneigh arrays
|
||||
int **pages; // neighbor list pages
|
||||
int maxpage; // # of pages currently allocated
|
||||
int pgsize; // size of neighbor page
|
||||
int oneatom; // max # of neighbors for one atom
|
||||
int npage; // current page in page list
|
||||
int *map; // 0 (C), 1 (H), or -1 (NULL) for each type
|
||||
|
||||
double cutlj; // user-specified LJ cutoff
|
||||
double cutljrebosq; // cut for when to compute
|
||||
|
|
@ -77,12 +79,12 @@ class PairAIREBO : public Pair {
|
|||
|
||||
double PCCf[5][5],PCCdfdx[5][5],PCCdfdy[5][5],PCHf[5][5];
|
||||
double PCHdfdx[5][5],PCHdfdy[5][5];
|
||||
double piCCf[5][5][10],piCCdfdx[5][5][10];
|
||||
double piCCdfdy[5][5][10],piCCdfdz[5][5][10];
|
||||
double piCHf[5][5][10],piCHdfdx[5][5][10];
|
||||
double piCHdfdy[5][5][10],piCHdfdz[5][5][10];
|
||||
double piHHf[5][5][10],piHHdfdx[5][5][10];
|
||||
double piHHdfdy[5][5][10],piHHdfdz[5][5][10];
|
||||
double piCCf[5][5][11],piCCdfdx[5][5][11];
|
||||
double piCCdfdy[5][5][11],piCCdfdz[5][5][11];
|
||||
double piCHf[5][5][11],piCHdfdx[5][5][11];
|
||||
double piCHdfdy[5][5][11],piCHdfdz[5][5][11];
|
||||
double piHHf[5][5][11],piHHdfdx[5][5][11];
|
||||
double piHHdfdy[5][5][11],piHHdfdz[5][5][11];
|
||||
double Tf[5][5][10],Tdfdx[5][5][10],Tdfdy[5][5][10],Tdfdz[5][5][10];
|
||||
|
||||
void REBO_neigh();
|
||||
|
|
@ -94,17 +96,12 @@ class PairAIREBO : public Pair {
|
|||
double bondorderLJ(int, int, double *, double, double,
|
||||
double *, double, double **, int);
|
||||
|
||||
double Sp(double, double, double, double &);
|
||||
double Sp2(double, double, double, double &);
|
||||
|
||||
double gSpline(double, double, int, double *, double *);
|
||||
double PijSpline(double, double, int, int, double *);
|
||||
double piRCSpline(double, double, double, int, int, double *);
|
||||
double TijSpline(double, double, double, double *);
|
||||
|
||||
double kronecker(int, int);
|
||||
|
||||
void add_pages(int);
|
||||
void add_pages();
|
||||
void read_file(char *);
|
||||
|
||||
double Sp5th(double, double *, double *);
|
||||
|
|
@ -113,6 +110,63 @@ class PairAIREBO : public Pair {
|
|||
void spline_init();
|
||||
|
||||
void allocate();
|
||||
|
||||
// ----------------------------------------------------------------------
|
||||
// S'(t) and S(t) cutoff functions
|
||||
// added to header for inlining
|
||||
// ----------------------------------------------------------------------
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
cutoff function Sprime
|
||||
return cutoff and dX = derivative
|
||||
no side effects
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
inline double Sp(double Xij, double Xmin, double Xmax, double &dX) const {
|
||||
double cutoff;
|
||||
|
||||
double t = (Xij-Xmin) / (Xmax-Xmin);
|
||||
if (t <= 0.0) {
|
||||
cutoff = 1.0;
|
||||
dX = 0.0;
|
||||
} else if (t >= 1.0) {
|
||||
cutoff = 0.0;
|
||||
dX = 0.0;
|
||||
} else {
|
||||
cutoff = 0.5 * (1.0+cos(t*MathConst::MY_PI));
|
||||
dX = (-0.5*MathConst::MY_PI*sin(t*MathConst::MY_PI)) / (Xmax-Xmin);
|
||||
}
|
||||
return cutoff;
|
||||
};
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
LJ cutoff function Sp2
|
||||
return cutoff and dX = derivative
|
||||
no side effects
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
inline double Sp2(double Xij, double Xmin, double Xmax, double &dX) const {
|
||||
double cutoff;
|
||||
|
||||
double t = (Xij-Xmin) / (Xmax-Xmin);
|
||||
if (t <= 0.0) {
|
||||
cutoff = 1.0;
|
||||
dX = 0.0;
|
||||
} else if (t >= 1.0) {
|
||||
cutoff = 0.0;
|
||||
dX = 0.0;
|
||||
} else {
|
||||
cutoff = (1.0-(t*t*(3.0-2.0*t)));
|
||||
dX = 6.0*(t*t-t) / (Xmax-Xmin);
|
||||
}
|
||||
return cutoff;
|
||||
};
|
||||
|
||||
/* kronecker delta function returning a double */
|
||||
inline double kronecker(const int a, const int b) const {
|
||||
return (a == b) ? 1.0 : 0.0;
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue