Merge pull request #51 from HITS-MBM/threshold_averaging

Threshold value of force and averaging
HITS-MBM · May 17, 2023 · 3dfde04 · 3dfde04
2 parents c359ccf + c946b1d
commit 3dfde04
Show file tree

Hide file tree

Showing 26 changed files with 6,238 additions and 4,679 deletions.
diff --git a/.vscode/launch.json b/.vscode/launch.json
@@ -9,7 +9,7 @@
             "type": "cppdbg",
             "request": "launch",
             "program": "${workspaceFolder}/build/bin/fda-integration-tests",
-            "args": [],
+            "args": ["--gtest_filter=AllFDATests/FDATest.Basic/4"],
             "stopAtEntry": false,
             "cwd": "${workspaceFolder}",
             "environment": [],

diff --git a/src/gromacs/fda/DistributedForces.cpp b/src/gromacs/fda/DistributedForces.cpp
@@ -153,6 +153,7 @@ void DistributedForces::write_summed_scalar(std::ostream& os, gmx::PaddedHostVec
             uint num_j = summed[i].size();
             os.write(reinterpret_cast<char*>(&num_j), sizeof(uint));
             for (uint p = 0; p != num_j; ++p) {
+                if (std::abs(scalar[i][p].force) < fda_settings.threshold) continue;
                 uint j = indices[i][p];
                 os.write(reinterpret_cast<char*>(&j), sizeof(uint));
                 real scalar = vector2signedscalar(summed[i][p].force.get_pointer(), x[i], x[j], box, fda_settings.v2s);
@@ -167,9 +168,12 @@ void DistributedForces::write_summed_scalar(std::ostream& os, gmx::PaddedHostVec
             for (size_t p = 0; p != summed_i.size(); ++p) {
                 size_t j = indices_i[p];
                 auto const& summed_j = summed_i[p];
-                os << i << " " << j << " "
-                   << vector2signedscalar(summed_j.force.get_pointer(), x[i], x[j], box, fda_settings.v2s) << " "
-                   << summed_j.type << std::endl;
+                auto&& summed_j_scalar_force = vector2signedscalar(summed_j.force.get_pointer(), x[i], x[j], box, fda_settings.v2s);
+                if (std::abs(summed_j_scalar_force) >= fda_settings.threshold) {
+                    os << i << " " << j << " "
+                       << summed_j_scalar_force << " "
+                       << summed_j.type << std::endl;
+                }
             }
         }
     }
@@ -186,6 +190,7 @@ void DistributedForces::write_scalar(std::ostream& os) const
             uint num_j = scalar[i].size();
             os.write(reinterpret_cast<char*>(&num_j), sizeof(uint));
             for (uint p = 0; p != num_j; ++p) {
+                if (std::abs(scalar[i][p].force) < fda_settings.threshold) continue;
                 uint j = indices[i][p];
                 os.write(reinterpret_cast<char*>(&j), sizeof(uint));
                 os.write(reinterpret_cast<const char*>(&scalar[i][p].force), sizeof(real));
@@ -199,9 +204,11 @@ void DistributedForces::write_scalar(std::ostream& os) const
             for (size_t p = 0; p != scalar_i.size(); ++p) {
                 size_t j = scalar_indices_i[p];
                 auto const& scalar_j = scalar_i[p];
-                os << i << " " << j << " "
-                   << scalar_j.force << " "
-                   << scalar_j.type << std::endl;
+                if (std::abs(scalar_j.force) >= fda_settings.threshold) {
+                    os << i << " " << j << " "
+                       << scalar_j.force << " "
+                       << scalar_j.type << std::endl;
+                }
             }
         }
     }

diff --git a/src/gromacs/nbnxm/kernels_reference/kernel_ref_inner.h b/src/gromacs/nbnxm/kernels_reference/kernel_ref_inner.h
@@ -353,13 +353,7 @@
 
 #ifdef CALC_ENERGIES
                 /* pairwise forces */
-				if (fabs(fcoul) > fda_threshold && fabs(fvdw) > fda_threshold) {
-					fda->add_nonbonded(cellInv[ai], cellInv[aj], fcoul, fvdw, dx, dy, dz);
-				} else if (fabs(fcoul) > fda_threshold) {
-					fda->add_nonbonded_single(cellInv[ai], cellInv[aj], fda::InteractionType_COULOMB, fcoul, dx, dy, dz);
-				} else if (fabs(fvdw) > fda_threshold) {
-					fda->add_nonbonded_single(cellInv[ai], cellInv[aj], fda::InteractionType_LJ, fscal, dx, dy, dz);
-				}
+                fda->add_nonbonded(cellInv[ai], cellInv[aj], fcoul, fvdw, dx, dy, dz);
 #endif
             }
 #    ifdef HALF_LJ
@@ -369,9 +363,7 @@
 
 #ifdef CALC_ENERGIES
 				/* pairwise forces */
-				if (fabs(fcoul) > fda_threshold) {
-					fda->add_nonbonded_single(cellInv[ai], cellInv[aj], fda::InteractionType_COULOMB, fcoul, dx, dy, dz);
-				}
+                fda->add_nonbonded_single(cellInv[ai], cellInv[aj], fda::InteractionType_COULOMB, fcoul, dx, dy, dz);
 #endif
             }
 #    endif
@@ -380,9 +372,7 @@
 
 #ifdef CALC_ENERGIES
             /* pairwise forces */
-            if (fabs(fscal) > fda_threshold) {
-            	fda->add_nonbonded_single(cellInv[ai], cellInv[aj], fda::InteractionType_LJ, fscal, dx, dy, dz);
-            }
+            fda->add_nonbonded_single(cellInv[ai], cellInv[aj], fda::InteractionType_LJ, fscal, dx, dy, dz);
 #endif
 #endif
             fx = fscal * dx;

diff --git a/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h b/src/gromacs/nbnxm/kernels_reference/kernel_ref_outer.h
@@ -130,7 +130,6 @@ void
     real qi[UNROLLI];
 
 #ifdef CALC_ENERGIES
-    real fda_threshold = fda->get_settings().threshold;
 #    ifndef ENERGY_GROUPS
 
     real Vvdw_ci, Vc_ci;

diff --git a/src/programs/mdrun/fda_tests/FDATest.cpp b/src/programs/mdrun/fda_tests/FDATest.cpp
@@ -150,8 +150,10 @@ std::vector<TestDataStructure> get_tests()
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_summed_renumbered", "pfa", "pfr"});
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_summed_atom_based", "pfa", ""});
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_summed_no_residue_based", "pfa", ""});
+    tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_threshold", "pfa", "pfr"});
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_average_0", "pfa", "pfr"});
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_average_3", "pfa", "pfr"});
+    tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_average_0_threshold", "pfa", "pfr"});
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_detailed", "pfa", "pfr"});
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_vector", "pfa", "pfr", "traj.trr", true});
     tests.push_back({"alagly_verlet", "alagly_pairwise_forces_scalar_detailed_nonbonded", "pfa", "pfr"});