From b96028eaf2490554cfffa5021455dc737c81fe72 Mon Sep 17 00:00:00 2001 From: Christoph Scherer Date: Sat, 22 Oct 2022 16:38:37 +0200 Subject: [PATCH 1/3] Update pair_threebody_table.cpp Correcting for hard coded ntheta = 79 in the extreme case that theta is exactly equal to 180.0 degrees. --- src/MANYBODY/pair_threebody_table.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MANYBODY/pair_threebody_table.cpp b/src/MANYBODY/pair_threebody_table.cpp index 044f69a8da..74367defba 100644 --- a/src/MANYBODY/pair_threebody_table.cpp +++ b/src/MANYBODY/pair_threebody_table.cpp @@ -708,7 +708,7 @@ void PairThreebodyTable::uf_lookup(Param *pm, double r12, double r13, double the if (r13 == (pm->mltable->rmin - 0.5 * dr)) { nr13 = 0; } nr13 -= nr12; ntheta = (theta - 0.00000001) / dtheta; - if (theta == 180.0) { ntheta = 79; } + if (theta == 180.0) { ntheta = (pm->mltable->ninput * 2)-1; } itable = 0; for (i = 0; i < nr12; i++) { itable += (pm->mltable->ninput - i); } itable += nr13; @@ -721,7 +721,7 @@ void PairThreebodyTable::uf_lookup(Param *pm, double r12, double r13, double the nr13 = (r13 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; if (r13 == (pm->mltable->rmin - 0.5 * dr)) { nr13 = 0; } ntheta = (theta - 0.00000001) / dtheta; - if (theta == 180.0) { ntheta = 79; } + if (theta == 180.0) { ntheta = (pm->mltable->ninput * 2)-1; } itable = nr12 * (pm->mltable->ninput); itable += nr13; itable *= (pm->mltable->ninput * 2); From 94627b3ef72ba3881bfd9baae46c0a18a81c3157 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Oct 2022 14:54:53 -0400 Subject: [PATCH 2/3] remove redundant curly braces --- src/MANYBODY/pair_threebody_table.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/MANYBODY/pair_threebody_table.cpp b/src/MANYBODY/pair_threebody_table.cpp index 74367defba..17254c83a4 100644 --- a/src/MANYBODY/pair_threebody_table.cpp +++ b/src/MANYBODY/pair_threebody_table.cpp @@ -445,7 +445,7 @@ void PairThreebodyTable::read_table(Table *tb, char *file, char *keyword, bool s char *line = reader.find_section_start(keyword); - if (!line) { error->one(FLERR, "Did not find keyword in table file"); } + if (!line) error->one(FLERR, "Did not find keyword in table file"); // read args on 2nd line of section // allocate table arrays for file values @@ -703,25 +703,25 @@ void PairThreebodyTable::uf_lookup(Param *pm, double r12, double r13, double the // if it is a symmetric threebody interaction, less table entries are required if (pm->symmetric) { nr12 = (r12 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; - if (r12 == (pm->mltable->rmin - 0.5 * dr)) { nr12 = 0; } + if (r12 == (pm->mltable->rmin - 0.5 * dr)) nr12 = 0; nr13 = (r13 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; - if (r13 == (pm->mltable->rmin - 0.5 * dr)) { nr13 = 0; } + if (r13 == (pm->mltable->rmin - 0.5 * dr)) nr13 = 0; nr13 -= nr12; ntheta = (theta - 0.00000001) / dtheta; - if (theta == 180.0) { ntheta = (pm->mltable->ninput * 2)-1; } + if (theta == 180.0) ntheta = (pm->mltable->ninput * 2) - 1; itable = 0; - for (i = 0; i < nr12; i++) { itable += (pm->mltable->ninput - i); } + for (i = 0; i < nr12; i++) itable += (pm->mltable->ninput - i); itable += nr13; itable *= (pm->mltable->ninput * 2); itable += ntheta; } else { // else, more (full) table entries are required nr12 = (r12 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; - if (r12 == (pm->mltable->rmin - 0.5 * dr)) { nr12 = 0; } + if (r12 == (pm->mltable->rmin - 0.5 * dr)) nr12 = 0; nr13 = (r13 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; - if (r13 == (pm->mltable->rmin - 0.5 * dr)) { nr13 = 0; } + if (r13 == (pm->mltable->rmin - 0.5 * dr)) nr13 = 0; ntheta = (theta - 0.00000001) / dtheta; - if (theta == 180.0) { ntheta = (pm->mltable->ninput * 2)-1; } + if (theta == 180.0) ntheta = (pm->mltable->ninput * 2) - 1; itable = nr12 * (pm->mltable->ninput); itable += nr13; itable *= (pm->mltable->ninput * 2); From cd621e74f5c34475ded5ceb5bf1326e2bba704f8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 22 Oct 2022 15:03:37 -0400 Subject: [PATCH 3/3] it is safer to do a >= comparison instead of == for floating point numbers --- src/MANYBODY/pair_threebody_table.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MANYBODY/pair_threebody_table.cpp b/src/MANYBODY/pair_threebody_table.cpp index 17254c83a4..e2b58975b7 100644 --- a/src/MANYBODY/pair_threebody_table.cpp +++ b/src/MANYBODY/pair_threebody_table.cpp @@ -708,7 +708,7 @@ void PairThreebodyTable::uf_lookup(Param *pm, double r12, double r13, double the if (r13 == (pm->mltable->rmin - 0.5 * dr)) nr13 = 0; nr13 -= nr12; ntheta = (theta - 0.00000001) / dtheta; - if (theta == 180.0) ntheta = (pm->mltable->ninput * 2) - 1; + if (theta >= 180.0) ntheta = (pm->mltable->ninput * 2) - 1; itable = 0; for (i = 0; i < nr12; i++) itable += (pm->mltable->ninput - i); itable += nr13; @@ -721,7 +721,7 @@ void PairThreebodyTable::uf_lookup(Param *pm, double r12, double r13, double the nr13 = (r13 - pm->mltable->rmin + 0.5 * dr - 0.00000001) / dr; if (r13 == (pm->mltable->rmin - 0.5 * dr)) nr13 = 0; ntheta = (theta - 0.00000001) / dtheta; - if (theta == 180.0) ntheta = (pm->mltable->ninput * 2) - 1; + if (theta >= 180.0) ntheta = (pm->mltable->ninput * 2) - 1; itable = nr12 * (pm->mltable->ninput); itable += nr13; itable *= (pm->mltable->ninput * 2);