Merge branch 'master' into documentation

marrink-lab · Nov 20, 2023 · 571283a · 571283a
2 parents a67a1d3 + 899a8df
commit 571283a
Show file tree

Hide file tree

Showing 2 changed files with 56 additions and 0 deletions.
diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py
@@ -89,11 +89,20 @@ def compute_force_constants(distance_matrix, lower_bound, upper_bound,
 
     The force constant can be modified with a decay function, and it can be
     bounded with a minimum threshold, or a distance upper and lower bonds.
+    
+    If decay_factor = decay_power = 0 all forces applied are = base_constant
+    
+    Forces applied to distances above upper_bound are removed.
+    Forces below minimum_force are removed.
+    
+    If decay_factor or decay_power != 0, forces below lower_bound are greater 
+    than base_constant, in which case they are set back to = base_constant    
     """
     constants = compute_decay(distance_matrix, lower_bound, decay_factor, decay_power)
     np.fill_diagonal(constants, 0)
     constants *= base_constant
     constants[constants < minimum_force] = 0
+    constants[constants > base_constant] = base_constant
     constants[distance_matrix > upper_bound] = 0
     return constants
 

diff --git a/vermouth/tests/test_apply_rubber_band.py b/vermouth/tests/test_apply_rubber_band.py
@@ -537,3 +537,50 @@ def test_bail_out_on_nan(caplog, test_molecule):
     assert record.getMessage() == required_warning
     assert len(caplog.records) == 1
     assert test_molecule.interactions['bonds'] == []
+
+
+@pytest.mark.parametrize('lower_bound, upper_bound, decay_factor, decay_power, base_constant, minimum_force, expected_output', 
+                          ([1, 2, 0, 0, 500, 400, np.array([[  0,500,500,  0],
+                                                            [500,  0,500,500],
+                                                            [500,500,  0,500],
+                                                            [  0,500,500,  0]])], # no decays, return the base constant within the bounds
+                           [1, 2, 0, 0, 500, 600, np.array([[  0,  0,  0,  0],
+                                                            [  0,  0,  0,  0],
+                                                            [  0,  0,  0,  0],
+                                                            [  0,  0,  0,  0]])], # all forces less than the minumum force
+                           [1, 0.5, 0, 0, 500, 600, np.array([[  0,  0,  0,  0],
+                                                              [  0,  0,  0,  0],
+                                                              [  0,  0,  0,  0],
+                                                              [  0,  0,  0,  0]])], # all distances larger than the upper bound
+                           [4, 2, 0, 0, 500, 600, np.array([[  0,  0,  0,  0],
+                                                            [  0,  0,  0,  0],
+                                                            [  0,  0,  0,  0],
+                                                            [  0,  0,  0,  0]])], # all distances less than the lower bound
+                           [1, 3, 0.1, 0, 500, 400, np.array([[  0.        , 452.41870902, 452.41870902, 452.41870902],
+                                                              [452.41870902,   0.        , 452.41870902, 452.41870902],
+                                                              [452.41870902, 452.41870902,   0.        , 452.41870902],
+                                                              [452.41870902, 452.41870902, 452.41870902,   0.        ]])], # with a decay factor
+                           [1, 3, 0, 2, 500, 400, np.array([[  0,500,500,500],
+                                                            [500,  0,500,500],
+                                                            [500,500,  0,500],
+                                                            [500,500,500,  0]])], # with a decay factor
+                           [1, 3, 0.1, 2, 500, 400, np.array([[  0.        , 500.        , 452.41870902,   0.        ],
+                                                              [500.        ,   0.        , 500.        , 452.41870902],
+                                                              [452.41870902, 500.        ,   0.        , 500.        ],
+                                                              [  0.        , 452.41870902, 500.        ,   0.        ]])], # with a complex decay
+                          )
+                        )
+def test_compute_force_constants(lower_bound, upper_bound, decay_factor, decay_power, base_constant, minimum_force, expected_output):
+    # Define the constant distance_matrix for the test cases
+    distance_matrix = np.array([[0, 1, 2, 3],
+                                [1, 0, 1, 2],
+                                [2, 1, 0, 1],
+                                [3, 2, 1, 0]])
+
+    # Call the compute_force_constants function with the constant distance_matrix and varied parameters
+    result = vermouth.processors.apply_rubber_band.compute_force_constants(distance_matrix, lower_bound, upper_bound, 
+                                                                           decay_factor, decay_power, base_constant, 
+                                                                           minimum_force)
+
+    # Assert the result against the expected output
+    assert result == pytest.approx(expected_output)