Fix improved SIMD angles()

avoid division by zero when the angle is (too close to) zero.

Fix a regression introduced in fb7a6480db5914a50a97618b3e6316663041aafb

diff --git a/src/gromacs/listed_forces/bonded.cpp b/src/gromacs/listed_forces/bonded.cpp
index 93a8bcd2cd325b9a42475120eb242a12264b0417..39876051f87c8ecf2560a3a2be053c8e39ff5720 100644 (file)
--- a/src/gromacs/listed_forces/bonded.cpp
+++ b/src/gromacs/listed_forces/bonded.cpp
@@ -1089,7 +1089,7 @@ angles(int             nbonds,
     SimdReal                                 rijx_S, rijy_S, rijz_S;
     SimdReal                                 rkjx_S, rkjy_S, rkjz_S;
     SimdReal                                 one_S(1.0);
-    SimdReal min_one_plus_eps_S(-1.0 + 2.0 * GMX_REAL_EPS); // Smallest number > -1
+    SimdReal one_min_eps_S(1.0_real - GMX_REAL_EPS); // Largest number < 1
 
     SimdReal                         rij_rkj_S;
     SimdReal                         nrij2_S, nrij_1_S;
@@ -1174,14 +1174,14 @@ angles(int             nbonds,
          */
         cos2_S = rij_rkj_S * rij_rkj_S / (nrij2_S * nrkj2_S);
 
-        /* To allow for 180 degrees, we take the max of cos and -1 + 1bit,
-         * so we can safely get the 1/sin from 1/sqrt(1 - cos^2).
-         * This also ensures that rounding errors would cause the argument
-         * of simdAcos to be < -1.
-         * Note that we do not take precautions for cos(0)=1, so the outer
-         * atoms in an angle should not be on top of each other.
+        /* To allow for 0 and 180 degrees, we need to avoid issues when
+         * the cosine is close to -1 and 1. For acos() the argument needs
+         * to be in that range. For cos^2 we take the min of cos^2 and 1 - 1bit,
+         * so we can safely compute 1/sin as 1/sqrt(1 - cos^2).
          */
-        cos_S = max(cos_S, min_one_plus_eps_S);
+        cos_S  = max(cos_S, -one_S);
+        cos_S  = min(cos_S, one_S);
+        cos2_S = min(cos2_S, one_min_eps_S);
 
         theta_S = acos(cos_S);
 

diff --git a/src/gromacs/listed_forces/tests/bonded.cpp b/src/gromacs/listed_forces/tests/bonded.cpp
index caaa6e09ad0be4359c6b52d7596550aa61a581a5..3ab4b4fd6659a66138d1ab9b6b9f9f81fd0ee843 100644 (file)
--- a/src/gromacs/listed_forces/tests/bonded.cpp
+++ b/src/gromacs/listed_forces/tests/bonded.cpp
@@ -675,6 +675,16 @@ std::vector<iListInput> c_InputRestraints = {
     { iListInput(2e-3, 1e-8).setHarmonic(F_RESTRANGLES, 100.0, 50.0, 110.0, 45.0) }
 };
 
+//! Function types for testing bond with zero length, has zero reference length to make physical sense.
+std::vector<iListInput> c_InputBondsZeroLength = {
+    { iListInput().setHarmonic(F_BONDS, 0.0, 500.0) },
+};
+
+//! Function types for testing angles with zero angle, has zero reference angle to make physical sense.
+std::vector<iListInput> c_InputAnglesZeroAngle = {
+    { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 0.0, 50.0) },
+};
+
 //! Coordinates for testing
 std::vector<PaddedVector<RVec>> c_coordinatesForTests = {
     { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.2 }, { 0.005, 0.0, 0.1 }, { -0.001, 0.1, 0.0 } },
@@ -682,6 +692,16 @@ std::vector<PaddedVector<RVec>> c_coordinatesForTests = {
     { { -0.1143, -0.0282, 0.0 }, { 0.0, 0.0434, 0.0 }, { 0.1185, -0.0138, 0.0 }, { -0.0195, 0.1498, 0.0 } }
 };
 
+//! Coordinates for testing bonds with zero length
+std::vector<PaddedVector<RVec>> c_coordinatesForTestsZeroBondLength = {
+    { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.005, 0.0, 0.1 }, { 0.005, 0.0, 0.1 } }
+};
+
+//! Coordinates for testing bonds with zero length
+std::vector<PaddedVector<RVec>> c_coordinatesForTestsZeroAngle = {
+    { { 0.005, 0.0, 0.1 }, { 0.0, 0.0, 0.0 }, { 0.005, 0.0, 0.1 }, { 0.5, 0.18, 0.22 } }
+};
+
 //! PBC values for testing
 std::vector<PbcType> c_pbcForTests = { PbcType::No, PbcType::XY, PbcType::Xyz };
 
@@ -716,6 +736,18 @@ INSTANTIATE_TEST_CASE_P(Restraints,
                         ::testing::Combine(::testing::ValuesIn(c_InputRestraints),
                                            ::testing::ValuesIn(c_coordinatesForTests),
                                            ::testing::ValuesIn(c_pbcForTests)));
+
+INSTANTIATE_TEST_CASE_P(BondZeroLength,
+                        ListedForcesTest,
+                        ::testing::Combine(::testing::ValuesIn(c_InputBondsZeroLength),
+                                           ::testing::ValuesIn(c_coordinatesForTestsZeroBondLength),
+                                           ::testing::ValuesIn(c_pbcForTests)));
+
+INSTANTIATE_TEST_CASE_P(AngleZero,
+                        ListedForcesTest,
+                        ::testing::Combine(::testing::ValuesIn(c_InputAnglesZeroAngle),
+                                           ::testing::ValuesIn(c_coordinatesForTestsZeroAngle),
+                                           ::testing::ValuesIn(c_pbcForTests)));
 #endif
 } // namespace
 

diff --git a/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_0.xml b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_0.xml
new file mode 100644 (file)
index 0000000..bab7c13
--- /dev/null
+++ b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_0.xml
@@ -0,0 +1,38 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <FunctionType Name="ANGLES">
+    <FEP Name="No">
+      <Real Name="Epot ">84.797965428018102</Real>
+      <Real Name="dVdlambda ">0</Real>
+      <Vector Name="Central shift forces">
+        <Real Name="X">-2.1316282072803006e-14</Real>
+        <Real Name="Y">-1.7763568394002505e-15</Real>
+        <Real Name="Z">0</Real>
+      </Vector>
+      <Sequence Name="Forces">
+        <Int Name="Length">4</Int>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">861.88775797312974</Real>
+          <Real Name="Y">318.05244565695551</Real>
+          <Real Name="Z">-43.094387898656521</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">-896.43676003712528</Real>
+          <Real Name="Y">-333.82837010100394</Real>
+          <Real Name="Z">209.27290807871069</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">34.549002063995523</Real>
+          <Real Name="Y">15.775924444048426</Real>
+          <Real Name="Z">-166.17852018005416</Real>
+        </Vector>
+      </Sequence>
+    </FEP>
+  </FunctionType>
+</ReferenceData>

diff --git a/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_1.xml b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_1.xml
new file mode 100644 (file)
index 0000000..bab7c13
--- /dev/null
+++ b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_1.xml
@@ -0,0 +1,38 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <FunctionType Name="ANGLES">
+    <FEP Name="No">
+      <Real Name="Epot ">84.797965428018102</Real>
+      <Real Name="dVdlambda ">0</Real>
+      <Vector Name="Central shift forces">
+        <Real Name="X">-2.1316282072803006e-14</Real>
+        <Real Name="Y">-1.7763568394002505e-15</Real>
+        <Real Name="Z">0</Real>
+      </Vector>
+      <Sequence Name="Forces">
+        <Int Name="Length">4</Int>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">861.88775797312974</Real>
+          <Real Name="Y">318.05244565695551</Real>
+          <Real Name="Z">-43.094387898656521</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">-896.43676003712528</Real>
+          <Real Name="Y">-333.82837010100394</Real>
+          <Real Name="Z">209.27290807871069</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">34.549002063995523</Real>
+          <Real Name="Y">15.775924444048426</Real>
+          <Real Name="Z">-166.17852018005416</Real>
+        </Vector>
+      </Sequence>
+    </FEP>
+  </FunctionType>
+</ReferenceData>

diff --git a/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_2.xml b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_2.xml
new file mode 100644 (file)
index 0000000..bab7c13
--- /dev/null
+++ b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_2.xml
@@ -0,0 +1,38 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <FunctionType Name="ANGLES">
+    <FEP Name="No">
+      <Real Name="Epot ">84.797965428018102</Real>
+      <Real Name="dVdlambda ">0</Real>
+      <Vector Name="Central shift forces">
+        <Real Name="X">-2.1316282072803006e-14</Real>
+        <Real Name="Y">-1.7763568394002505e-15</Real>
+        <Real Name="Z">0</Real>
+      </Vector>
+      <Sequence Name="Forces">
+        <Int Name="Length">4</Int>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">861.88775797312974</Real>
+          <Real Name="Y">318.05244565695551</Real>
+          <Real Name="Z">-43.094387898656521</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">-896.43676003712528</Real>
+          <Real Name="Y">-333.82837010100394</Real>
+          <Real Name="Z">209.27290807871069</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">34.549002063995523</Real>
+          <Real Name="Y">15.775924444048426</Real>
+          <Real Name="Z">-166.17852018005416</Real>
+        </Vector>
+      </Sequence>
+    </FEP>
+  </FunctionType>
+</ReferenceData>

diff --git a/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_0.xml b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_0.xml
new file mode 100644 (file)
index 0000000..160a61d
--- /dev/null
+++ b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_0.xml
@@ -0,0 +1,38 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <FunctionType Name="BONDS">
+    <FEP Name="No">
+      <Real Name="Epot ">2.5062500000000005</Real>
+      <Real Name="dVdlambda ">0</Real>
+      <Vector Name="Central shift forces">
+        <Real Name="X">0</Real>
+        <Real Name="Y">0</Real>
+        <Real Name="Z">0</Real>
+      </Vector>
+      <Sequence Name="Forces">
+        <Int Name="Length">4</Int>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">2.5</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">50</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">-2.5</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">-50</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+      </Sequence>
+    </FEP>
+  </FunctionType>
+</ReferenceData>

diff --git a/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_1.xml b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_1.xml
new file mode 100644 (file)
index 0000000..160a61d
--- /dev/null
+++ b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_1.xml
@@ -0,0 +1,38 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <FunctionType Name="BONDS">
+    <FEP Name="No">
+      <Real Name="Epot ">2.5062500000000005</Real>
+      <Real Name="dVdlambda ">0</Real>
+      <Vector Name="Central shift forces">
+        <Real Name="X">0</Real>
+        <Real Name="Y">0</Real>
+        <Real Name="Z">0</Real>
+      </Vector>
+      <Sequence Name="Forces">
+        <Int Name="Length">4</Int>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">2.5</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">50</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">-2.5</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">-50</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+      </Sequence>
+    </FEP>
+  </FunctionType>
+</ReferenceData>

diff --git a/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_2.xml b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_2.xml
new file mode 100644 (file)
index 0000000..160a61d
--- /dev/null
+++ b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_2.xml
@@ -0,0 +1,38 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <FunctionType Name="BONDS">
+    <FEP Name="No">
+      <Real Name="Epot ">2.5062500000000005</Real>
+      <Real Name="dVdlambda ">0</Real>
+      <Vector Name="Central shift forces">
+        <Real Name="X">0</Real>
+        <Real Name="Y">0</Real>
+        <Real Name="Z">0</Real>
+      </Vector>
+      <Sequence Name="Forces">
+        <Int Name="Length">4</Int>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">2.5</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">50</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">-2.5</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">-50</Real>
+        </Vector>
+        <Vector>
+          <Real Name="X">0</Real>
+          <Real Name="Y">0</Real>
+          <Real Name="Z">0</Real>
+        </Vector>
+      </Sequence>
+    </FEP>
+  </FunctionType>
+</ReferenceData>