Merge pull request #1636 from glotzerlab/test/sdf-changes

Fix SDF computation with patchy potentials.
glotzerlab · Oct 23, 2023 · 9953645 · 9953645
2 parents 5183c32 + 7f98259
commit 9953645
Show file tree

Hide file tree

Showing 2 changed files with 56 additions and 10 deletions.
diff --git a/hoomd/hpmc/ComputeSDF.h b/hoomd/hpmc/ComputeSDF.h
@@ -616,8 +616,15 @@ template<class Shape> void ComputeSDF<Shape>::countHistogramLinearSearch(uint64_
                                     if (u_ij_new != u_ij_0)
                                         {
                                         min_bin_compression = bin_to_sample;
-                                        hist_weight_ptl_i_compression
-                                            = 1.0 - fast::exp(-(u_ij_new - u_ij_0));
+                                        if (u_ij_new < u_ij_0)
+                                            {
+                                            hist_weight_ptl_i_compression = 0;
+                                            }
+                                        else
+                                            {
+                                            hist_weight_ptl_i_compression
+                                                = 1.0 - fast::exp(-(u_ij_new - u_ij_0));
+                                            }
                                         }
                                     } // end if (!hard_overlap && m_mc->getPatchEnergy())
                                 }     // end loop over bins for compression
@@ -668,8 +675,15 @@ template<class Shape> void ComputeSDF<Shape>::countHistogramLinearSearch(uint64_
                                     if (u_ij_new != u_ij_0)
                                         {
                                         min_bin_expansion = bin_to_sample;
-                                        hist_weight_ptl_i_expansion
-                                            = 1.0 - fast::exp(-(u_ij_new - u_ij_0));
+                                        if (u_ij_new < u_ij_0)
+                                            {
+                                            hist_weight_ptl_i_expansion = 0;
+                                            }
+                                        else
+                                            {
+                                            hist_weight_ptl_i_expansion
+                                                = 1.0 - fast::exp(-(u_ij_new - u_ij_0));
+                                            }
                                         }
                                     } // end if (!hard_overlap && m_mc->getPatchEnergy())
                                 }     // end loop over histogram bins for expansions

diff --git a/hoomd/hpmc/pytest/test_compute_sdf.py b/hoomd/hpmc/pytest/test_compute_sdf.py
@@ -277,26 +277,58 @@ def test_linear_search_path(simulation_factory, two_particle_snapshot_factory):
         assert (numpy.count_nonzero(sdf_result) == 1)
 
     # add pair potential
-    square_well = rf'''float rsq = dot(r_ij, r_ij);
-                    float rcut = {r_patch};
+    square_well = r'''float rsq = dot(r_ij, r_ij);
+                    float rcut = param_array[1];
                     if (rsq > 2*rcut)
                         return 0.0f;
                     else
-                        return -{epsilon};'''
+                        return -param_array[0];'''
     patch = hoomd.hpmc.pair.user.CPPPotential(r_cut=2.5 * r_patch,
                                               code=square_well,
-                                              param_array=[])
+                                              param_array=[epsilon, r_patch])
     mc.pair_potential = patch
 
-    # assert we have an overlap in the first bin with the expected weight and
-    # that the SDF is zero everywhere else
+    # for a soft overlap with a negative change in energy, the weight of the
+    # overlap is zero, so we still expect all zeros in the sdf array
+    sim.run(1)
+    neg_mayerF = 1 - numpy.exp(epsilon)
+    sdf_result = sdf.sdf_compression
+    if sim.device.communicator.rank == 0:
+        assert (numpy.count_nonzero(sdf_result) == 0)
+
+    # change sign of epsilon, so that now there is a positive energy soft
+    # overlap in bin 0
+    epsilon *= -1
+    patch.param_array[0] = epsilon
     sim.run(1)
     neg_mayerF = 1 - numpy.exp(epsilon)
     sdf_result = sdf.sdf_compression
     if sim.device.communicator.rank == 0:
         assert (numpy.count_nonzero(sdf_result) == 1)
         assert (sdf_result[0] == neg_mayerF * norm_factor)
 
+    # extend patches so that there is a soft overlap with positive energy
+    # in the configuration and there will be a change in energy on expansions
+    # the change in energy is < 0 so the weight should be 0 and
+    # sdf_expansion should be all zeros
+    patch.param_array[1] += 2 * dx
+    sim.run(1)
+    sdf_result = sdf.sdf_expansion
+    if sim.device.communicator.rank == 0:
+        assert (numpy.count_nonzero(sdf_result) == 0)
+
+    # change sign of epsilon so now the change in energy on expansion is
+    # positive and the weight is nonzero and one of the sdf_expansion counts
+    # is nonzero
+    epsilon *= -1
+    patch.param_array[0] = epsilon
+    sim.run(1)
+    neg_mayerF = 1 - numpy.exp(-epsilon)
+    sdf_result = sdf.sdf_expansion
+    if sim.device.communicator.rank == 0:
+        assert (numpy.count_nonzero(sdf_result) == 1)
+        assert (sdf_result[-1] == neg_mayerF * norm_factor)
+
 
 @pytest.mark.cpu
 @pytest.mark.serial