Commit fa2e5ac2 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

handle lookup exceptions consistently across energy and energy+force lookup in...

handle lookup exceptions consistently across energy and energy+force lookup in bond/angle style table
parent b7c74926
Loading
Loading
Loading
Loading
+15 −5
Original line number Diff line number Diff line
@@ -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<int> (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<int> ( 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];
+16 −9
Original line number Diff line number Diff line
@@ -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<int> ((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<int> ((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;