From: Gilles Gouaillardet Date: Wed, 2 Sep 2020 09:01:44 +0000 (+0000) Subject: Fix improved SIMD angles() X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=132f5c9f1ad409fa3f54d18f082a3271c82d1c31;p=alexxy%2Fgromacs.git 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 93a8bcd2cd..39876051f8 100644 --- 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 caaa6e09ad..3ab4b4fd66 100644 --- a/src/gromacs/listed_forces/tests/bonded.cpp +++ b/src/gromacs/listed_forces/tests/bonded.cpp @@ -675,6 +675,16 @@ std::vector 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 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 c_InputAnglesZeroAngle = { + { iListInput(2e-3, 1e-8).setHarmonic(F_ANGLES, 0.0, 50.0) }, +}; + //! Coordinates for testing std::vector> 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> 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> 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> 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 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 index 0000000000..bab7c13cfc --- /dev/null +++ b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_0.xml @@ -0,0 +1,38 @@ + + + + + + 84.797965428018102 + 0 + + -2.1316282072803006e-14 + -1.7763568394002505e-15 + 0 + + + 4 + + 0 + 0 + 0 + + + 861.88775797312974 + 318.05244565695551 + -43.094387898656521 + + + -896.43676003712528 + -333.82837010100394 + 209.27290807871069 + + + 34.549002063995523 + 15.775924444048426 + -166.17852018005416 + + + + + 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 index 0000000000..bab7c13cfc --- /dev/null +++ b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_1.xml @@ -0,0 +1,38 @@ + + + + + + 84.797965428018102 + 0 + + -2.1316282072803006e-14 + -1.7763568394002505e-15 + 0 + + + 4 + + 0 + 0 + 0 + + + 861.88775797312974 + 318.05244565695551 + -43.094387898656521 + + + -896.43676003712528 + -333.82837010100394 + 209.27290807871069 + + + 34.549002063995523 + 15.775924444048426 + -166.17852018005416 + + + + + 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 index 0000000000..bab7c13cfc --- /dev/null +++ b/src/gromacs/listed_forces/tests/refdata/AngleZero_ListedForcesTest_Ifunc_2.xml @@ -0,0 +1,38 @@ + + + + + + 84.797965428018102 + 0 + + -2.1316282072803006e-14 + -1.7763568394002505e-15 + 0 + + + 4 + + 0 + 0 + 0 + + + 861.88775797312974 + 318.05244565695551 + -43.094387898656521 + + + -896.43676003712528 + -333.82837010100394 + 209.27290807871069 + + + 34.549002063995523 + 15.775924444048426 + -166.17852018005416 + + + + + 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 index 0000000000..160a61dd25 --- /dev/null +++ b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_0.xml @@ -0,0 +1,38 @@ + + + + + + 2.5062500000000005 + 0 + + 0 + 0 + 0 + + + 4 + + 0 + 0 + 0 + + + 2.5 + 0 + 50 + + + -2.5 + 0 + -50 + + + 0 + 0 + 0 + + + + + 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 index 0000000000..160a61dd25 --- /dev/null +++ b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_1.xml @@ -0,0 +1,38 @@ + + + + + + 2.5062500000000005 + 0 + + 0 + 0 + 0 + + + 4 + + 0 + 0 + 0 + + + 2.5 + 0 + 50 + + + -2.5 + 0 + -50 + + + 0 + 0 + 0 + + + + + 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 index 0000000000..160a61dd25 --- /dev/null +++ b/src/gromacs/listed_forces/tests/refdata/BondZeroLength_ListedForcesTest_Ifunc_2.xml @@ -0,0 +1,38 @@ + + + + + + 2.5062500000000005 + 0 + + 0 + 0 + 0 + + + 4 + + 0 + 0 + 0 + + + 2.5 + 0 + 50 + + + -2.5 + 0 + -50 + + + 0 + 0 + 0 + + + + +