diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp index 4d9007adb7..6e145efa10 100644 --- a/src/MOLECULE/angle_table.cpp +++ b/src/MOLECULE/angle_table.cpp @@ -609,18 +609,22 @@ double AngleTable::splint(double *xa, double *ya, double *y2a, int n, double x) void AngleTable::uf_lookup(int type, double x, double &u, double &f) { - int itable; - double fraction,a,b; + if (!ISFINITE(x)) { + error->one(FLERR,"Illegal angle in angle style table"); + } - Table *tb = &tables[tabindex[type]]; + double fraction,a,b; + const Table *tb = &tables[tabindex[type]]; + int itable = static_cast (x * tb->invdelta); + + if (itable < 0) itable = 0; + if (itable >= tablength) itable = tablength-1; if (tabstyle == LINEAR) { - itable = static_cast ( x * tb->invdelta); fraction = (x - tb->ang[itable]) * tb->invdelta; u = tb->e[itable] + fraction*tb->de[itable]; f = tb->f[itable] + fraction*tb->df[itable]; } else if (tabstyle == SPLINE) { - itable = static_cast ( x * tb->invdelta); fraction = (x - tb->ang[itable]) * tb->invdelta; b = (x - tb->ang[itable]) * tb->invdelta; @@ -640,17 +644,21 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f) void AngleTable::u_lookup(int type, double x, double &u) { - int itable; - double fraction,a,b; + if (!ISFINITE(x)) { + error->one(FLERR,"Illegal angle in angle style table"); + } - Table *tb = &tables[tabindex[type]]; + double fraction,a,b; + const Table *tb = &tables[tabindex[type]]; + int itable = static_cast ( x * tb->invdelta); + + if (itable < 0) itable = 0; + if (itable >= tablength) itable = tablength-1; if (tabstyle == LINEAR) { - itable = static_cast ( x * tb->invdelta); fraction = (x - tb->ang[itable]) * tb->invdelta; u = tb->e[itable] + fraction*tb->de[itable]; } else if (tabstyle == SPLINE) { - itable = static_cast ( x * tb->invdelta); fraction = (x - tb->ang[itable]) * tb->invdelta; b = (x - tb->ang[itable]) * tb->invdelta; diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp index 38cbe7e406..4f8db66757 100644 --- a/src/MOLECULE/bond_table.cpp +++ b/src/MOLECULE/bond_table.cpp @@ -590,29 +590,29 @@ double BondTable::splint(double *xa, double *ya, double *y2a, int n, double x) void BondTable::uf_lookup(int type, double x, double &u, double &f) { - int itable; + if (!ISFINITE(x)) { + error->one(FLERR,"Illegal bond in bond style table"); + } + double fraction,a,b; char estr[128]; - - Table *tb = &tables[tabindex[type]]; - if (x < tb->lo) { + const Table *tb = &tables[tabindex[type]]; + const int itable = static_cast ((x - tb->lo) * tb->invdelta); + if (itable < 0) { sprintf(estr,"Bond length < table inner cutoff: " "type %d length %g",type,x); error->one(FLERR,estr); - } - if (x > tb->hi) { + } else if (itable >= tablength) { sprintf(estr,"Bond length > table outer cutoff: " "type %d length %g",type,x); error->one(FLERR,estr); } if (tabstyle == LINEAR) { - itable = static_cast ((x - tb->lo) * tb->invdelta); fraction = (x - tb->r[itable]) * tb->invdelta; u = tb->e[itable] + fraction*tb->de[itable]; f = tb->f[itable] + fraction*tb->df[itable]; } else if (tabstyle == SPLINE) { - itable = static_cast ((x - tb->lo) * tb->invdelta); fraction = (x - tb->r[itable]) * tb->invdelta; b = (x - tb->r[itable]) * tb->invdelta; @@ -633,19 +633,28 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f) void BondTable::u_lookup(int type, double x, double &u) { - int itable; - double fraction,a,b; + if (!ISFINITE(x)) { + error->one(FLERR,"Illegal bond in bond style table"); + } - Table *tb = &tables[tabindex[type]]; - x = MAX(x,tb->lo); - x = MIN(x,tb->hi); + double fraction,a,b; + char estr[128]; + const Table *tb = &tables[tabindex[type]]; + const int itable = static_cast ((x - tb->lo) * tb->invdelta); + if (itable < 0) { + sprintf(estr,"Bond length < table inner cutoff: " + "type %d length %g",type,x); + error->one(FLERR,estr); + } else if (itable >= tablength) { + sprintf(estr,"Bond length > table outer cutoff: " + "type %d length %g",type,x); + error->one(FLERR,estr); + } if (tabstyle == LINEAR) { - itable = static_cast ((x - tb->lo) * tb->invdelta); fraction = (x - tb->r[itable]) * tb->invdelta; u = tb->e[itable] + fraction*tb->de[itable]; } else if (tabstyle == SPLINE) { - itable = static_cast ((x - tb->lo) * tb->invdelta); fraction = (x - tb->r[itable]) * tb->invdelta; b = (x - tb->r[itable]) * tb->invdelta; diff --git a/src/output.cpp b/src/output.cpp index ce593ec6ae..11c6fa073e 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -844,9 +844,9 @@ void Output::memory_usage() MPI_Reduce(&mbytes,&mbavg,1,MPI_DOUBLE,MPI_SUM,0,world); MPI_Reduce(&mbytes,&mbmin,1,MPI_DOUBLE,MPI_MIN,0,world); MPI_Reduce(&mbytes,&mbmax,1,MPI_DOUBLE,MPI_MAX,0,world); - mbavg /= comm->nprocs; if (comm->me == 0) { + mbavg /= comm->nprocs; if (screen) fprintf(screen,"Per MPI rank memory allocation (min/avg/max) = " "%.4g | %.4g | %.4g Mbytes\n",mbmin,mbavg,mbmax);