From cbbb66194cebc013fc4d313a77d2176689ffda7e Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 25 Oct 2011 16:17:27 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7196 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MANYBODY/pair_airebo.cpp | 168 ++++++++++++----------------------- src/MANYBODY/pair_airebo.h | 86 ++++++++++++++---- 2 files changed, 127 insertions(+), 127 deletions(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 6f8ef42511..c2625072cc 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -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; diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index 166e1ed4ed..0c4795ab0d 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -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; + }; + }; }