Fix redmine #1109: not-a-number values in punctual stress

src/gromacs/fda/FDA.cpp

@@ -235,6 +235,12 @@ void FDA::add_angle(int ai, int aj, int ak, rvec f_i, rvec f_j, rvec f_k)
 {
     rvec uf_i, uf_j, uf_k, f_j_i, f_j_k, f_i_k;
     real nf_j_i, nf_j_k;
+
+    // below computation can sometimes return before finishing to avoid division with very small numbers;
+    // this situation can occur f.e. when all f_i, f_j, f_k and f_l are (almost) zero;
+    // in this case there is no call to fda_add_bonded, no pairwise forces are recorded (which is different from recording zero forces!)
+    if (norm(f_i) + norm(f_j) + norm(f_k) == 0.0) return;
+
     unitv(f_i, uf_i);
     unitv(f_j, uf_j);
     unitv(f_k, uf_k);

src/gromacs/fda/tests/CMakeLists.txt

@@ -3,6 +3,7 @@ set(exename "fda-tests")
 
 gmx_add_gtest_executable(
     ${exename}
+    FDATest.cpp
     PairwiseForcesTest.cpp
 )
 

src/gromacs/fda/tests/FDATest.cpp

@@ -0,0 +1,49 @@
+/*
+ * FDA.cpp
+ *
+ *  Created on: Sep 15, 2017
+ *      Author: Bernd Doser, HITS gGmbH <[email protected]>
+ */
+
+#include <iostream>
+#include <gtest/gtest.h>
+#include "gromacs/fda/FDA.h"
+#include "testutils/integrationtests.h"
+
+namespace fda {
+
+//! Test fixture for PairwiseForces
+class FDATest : public ::gmx::test::IntegrationTestFixture
+{};
+
+TEST_F(FDATest, add_angle)
+{
+    FDASettings fda_settings;
+    fda_settings.atom_based_result_type = ResultType::NO;
+    fda_settings.residue_based_result_type = ResultType::PUNCTUAL_STRESS;
+    fda_settings.one_pair = OnePair::SUMMED;
+    fda_settings.v2s = Vector2Scalar::NORM;
+    fda_settings.residues_renumber = ResiduesRenumber::AUTO;
+    fda_settings.no_end_zeros = false;
+    fda_settings.syslen_atoms = 0;
+    fda_settings.syslen_residues = 0;
+    fda_settings.time_averaging_period = 1;
+    fda_settings.type = InteractionType_NONE;
+    fda_settings.nonbonded_exclusion_on = true;
+    fda_settings.bonded_exclusion_on = true;
+    fda_settings.index_group1 = -1;
+    fda_settings.index_group2 = -1;
+    fda_settings.groups = nullptr;
+    fda_settings.groupnames = nullptr;
+
+    FDA fda(fda_settings);
+    rvec f_i{0.0, 0.0, 0.0};
+    rvec f_j{0.0, 0.0, 0.0};
+    rvec f_k{0.0, 0.0, 0.0};
+    fda.add_angle(0, 1, 2, f_i, f_j, f_k);
+
+}
+
+
+
+} // namespace fda