From cee94da85e1675bde9928c49d673f5d8331be179 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 13 Oct 2017 00:15:15 -0400 Subject: [PATCH 1/4] bugfix. avoiding operating on uninitialize data. closes #690 --- src/output.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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); From b7c74926088515c20e45350bf352697d44d57bd2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 13 Oct 2017 09:50:56 -0400 Subject: [PATCH 2/4] handle invalid lookup for bond/angle tabulation --- src/MOLECULE/angle_table.cpp | 16 +++++++--------- src/MOLECULE/bond_table.cpp | 24 +++++++++++++----------- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp index 4d9007adb7..b9f467f463 100644 --- a/src/MOLECULE/angle_table.cpp +++ b/src/MOLECULE/angle_table.cpp @@ -609,18 +609,18 @@ 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; - Table *tb = &tables[tabindex[type]]; + const Table *tb = &tables[tabindex[type]]; + const int itable = static_cast (x * tb->invdelta); - if (tabstyle == LINEAR) { - itable = static_cast ( x * tb->invdelta); + if ((itable < 0) || (itable >= tablength) || (!ISFINITE(itable))) { + error->one(FLERR,"Illegal angle in angle style table"); + } else if (tabstyle == LINEAR) { 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 +640,15 @@ 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; - Table *tb = &tables[tabindex[type]]; + const Table *tb = &tables[tabindex[type]]; + const int itable = static_cast ( x * tb->invdelta); 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..a641f41a6f 100644 --- a/src/MOLECULE/bond_table.cpp +++ b/src/MOLECULE/bond_table.cpp @@ -590,29 +590,28 @@ 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; 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); + } else if (!ISFINITE(itable)) { + error->one(FLERR,"Illegal bond length in bond style table"); } 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 +632,22 @@ 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; - Table *tb = &tables[tabindex[type]]; + if (!ISFINITE(x)) { + u = 0.0; + return; + } + + const Table *tb = &tables[tabindex[type]]; x = MAX(x,tb->lo); x = MIN(x,tb->hi); + const int itable = static_cast ((x - tb->lo) * tb->invdelta); 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; From fa2e5ac2d90234ae6b6d76b0b5d39616c0a5f76c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 13 Oct 2017 10:13:34 -0400 Subject: [PATCH 3/4] handle lookup exceptions consistently across energy and energy+force lookup in bond/angle style table --- src/MOLECULE/angle_table.cpp | 20 +++++++++++++++----- src/MOLECULE/bond_table.cpp | 25 ++++++++++++++++--------- 2 files changed, 31 insertions(+), 14 deletions(-) diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp index b9f467f463..a3201153da 100644 --- a/src/MOLECULE/angle_table.cpp +++ b/src/MOLECULE/angle_table.cpp @@ -609,14 +609,18 @@ 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) { - double fraction,a,b; + if (!ISFINITE(x)) { + error->one(FLERR,"Illegal angle in angle style table"); + } + double fraction,a,b; const Table *tb = &tables[tabindex[type]]; const int itable = static_cast (x * tb->invdelta); - if ((itable < 0) || (itable >= tablength) || (!ISFINITE(itable))) { - error->one(FLERR,"Illegal angle in angle style table"); - } else if (tabstyle == LINEAR) { + if (itable < 0) itable = 0; + if (itable >= tablength) itable = tablength-1; + + if (tabstyle == LINEAR) { fraction = (x - tb->ang[itable]) * tb->invdelta; u = tb->e[itable] + fraction*tb->de[itable]; f = tb->f[itable] + fraction*tb->df[itable]; @@ -640,11 +644,17 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f) void AngleTable::u_lookup(int type, double x, double &u) { - double fraction,a,b; + if (!ISFINITE(x)) { + error->one(FLERR,"Illegal angle in angle style table"); + } + double fraction,a,b; const Table *tb = &tables[tabindex[type]]; const int itable = static_cast ( x * tb->invdelta); + if (itable < 0) itable = 0; + if (itable >= tablength) itable = tablength-1; + if (tabstyle == LINEAR) { fraction = (x - tb->ang[itable]) * tb->invdelta; u = tb->e[itable] + fraction*tb->de[itable]; diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp index a641f41a6f..4f8db66757 100644 --- a/src/MOLECULE/bond_table.cpp +++ b/src/MOLECULE/bond_table.cpp @@ -590,9 +590,12 @@ 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) { + if (!ISFINITE(x)) { + error->one(FLERR,"Illegal bond in bond style table"); + } + 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) { @@ -603,8 +606,6 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f) sprintf(estr,"Bond length > table outer cutoff: " "type %d length %g",type,x); error->one(FLERR,estr); - } else if (!ISFINITE(itable)) { - error->one(FLERR,"Illegal bond length in bond style table"); } if (tabstyle == LINEAR) { @@ -632,17 +633,23 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f) void BondTable::u_lookup(int type, double x, double &u) { - double fraction,a,b; - if (!ISFINITE(x)) { - u = 0.0; - return; + error->one(FLERR,"Illegal bond in bond style table"); } + double fraction,a,b; + char estr[128]; const Table *tb = &tables[tabindex[type]]; - x = MAX(x,tb->lo); - x = MIN(x,tb->hi); 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) { fraction = (x - tb->r[itable]) * tb->invdelta; From 3df9caf435339585459b7bc00a54ab8f2029f3de Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 13 Oct 2017 10:29:49 -0400 Subject: [PATCH 4/4] drop const qualifier to allow bracketing of lookup index --- src/MOLECULE/angle_table.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp index a3201153da..6e145efa10 100644 --- a/src/MOLECULE/angle_table.cpp +++ b/src/MOLECULE/angle_table.cpp @@ -615,7 +615,7 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f) double fraction,a,b; const Table *tb = &tables[tabindex[type]]; - const int itable = static_cast (x * tb->invdelta); + int itable = static_cast (x * tb->invdelta); if (itable < 0) itable = 0; if (itable >= tablength) itable = tablength-1; @@ -650,7 +650,7 @@ void AngleTable::u_lookup(int type, double x, double &u) double fraction,a,b; const Table *tb = &tables[tabindex[type]]; - const int itable = static_cast ( x * tb->invdelta); + int itable = static_cast ( x * tb->invdelta); if (itable < 0) itable = 0; if (itable >= tablength) itable = tablength-1;