From d0feb7fa3c6dc6ec931b7bb217c4079a189631b3 Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Tue, 31 Oct 2023 16:38:34 +0100 Subject: [PATCH 1/8] ensured el flag is used --- vermouth/processors/apply_rubber_band.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index 9342631be..7da4b13ce 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -339,7 +339,7 @@ def apply_rubber_band(molecule, selector, to_key = idx_to_node[selection[to_idx]] force_constant = constants[from_idx, to_idx] length = distance_matrix[from_idx, to_idx] - if force_constant > minimum_force: + if (force_constant > minimum_force) and (force_constant < base_constant): molecule.add_interaction( type_='bonds', atoms=(from_key, to_key), From 8ef0e4244b4fc43f329d41cb0bae10043b83fc13 Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Tue, 31 Oct 2023 17:17:08 +0100 Subject: [PATCH 2/8] fixed el threshold application --- vermouth/processors/apply_rubber_band.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index 7da4b13ce..4cabd0d1c 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -93,7 +93,9 @@ def compute_force_constants(distance_matrix, lower_bound, upper_bound, 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] = 0 constants[distance_matrix > upper_bound] = 0 return constants @@ -339,7 +341,7 @@ def apply_rubber_band(molecule, selector, to_key = idx_to_node[selection[to_idx]] force_constant = constants[from_idx, to_idx] length = distance_matrix[from_idx, to_idx] - if (force_constant > minimum_force) and (force_constant < base_constant): + if (force_constant > minimum_force): molecule.add_interaction( type_='bonds', atoms=(from_key, to_key), From edbfb6f9781098be8b5814373a4206920d21ef12 Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Tue, 31 Oct 2023 17:22:00 +0100 Subject: [PATCH 3/8] removed brackets --- vermouth/processors/apply_rubber_band.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index 4cabd0d1c..5ee6b672e 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -341,7 +341,7 @@ def apply_rubber_band(molecule, selector, to_key = idx_to_node[selection[to_idx]] force_constant = constants[from_idx, to_idx] length = distance_matrix[from_idx, to_idx] - if (force_constant > minimum_force): + if force_constant > minimum_force: molecule.add_interaction( type_='bonds', atoms=(from_key, to_key), From bcf46994656162473fcb7b53716360353d1b098a Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Tue, 31 Oct 2023 17:34:53 +0100 Subject: [PATCH 4/8] didn't account for constant forces --- vermouth/processors/apply_rubber_band.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index 5ee6b672e..4f0f23ef4 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -93,10 +93,10 @@ def compute_force_constants(distance_matrix, lower_bound, upper_bound, 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] = 0 constants[distance_matrix > upper_bound] = 0 + constants[distance_matrix < lower_bound] = 0 return constants From 99c04520b79e55060698621521a5c747015e7ae7 Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Wed, 1 Nov 2023 14:15:54 +0100 Subject: [PATCH 5/8] added tests for rubber band application --- vermouth/tests/test_apply_rubber_band.py | 35 ++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/vermouth/tests/test_apply_rubber_band.py b/vermouth/tests/test_apply_rubber_band.py index b7fa5cad1..17ab609e1 100644 --- a/vermouth/tests/test_apply_rubber_band.py +++ b/vermouth/tests/test_apply_rubber_band.py @@ -537,3 +537,38 @@ 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,500]])], # 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 + ) + ) +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.all() == expected_output.all() From 9104ceca4f73e6a817b867de830f199709d56017 Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Wed, 1 Nov 2023 15:12:33 +0100 Subject: [PATCH 6/8] added decay tests --- vermouth/tests/test_apply_rubber_band.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/vermouth/tests/test_apply_rubber_band.py b/vermouth/tests/test_apply_rubber_band.py index 17ab609e1..131525157 100644 --- a/vermouth/tests/test_apply_rubber_band.py +++ b/vermouth/tests/test_apply_rubber_band.py @@ -543,19 +543,31 @@ def test_bail_out_on_nan(caplog, test_molecule): ([1, 2, 0, 0, 500, 400, np.array([[ 0,500,500, 0], [500, 0,500,500], [500,500, 0,500], - [ 0,500,500,500]])], # no decays, return the base constant within the bounds + [ 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 + [ 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): @@ -571,4 +583,4 @@ def test_compute_force_constants(lower_bound, upper_bound, decay_factor, decay_p minimum_force) # Assert the result against the expected output - assert result.all() == expected_output.all() + assert result == pytest.approx(expected_output) \ No newline at end of file From 1617e842e98db5cb5de32a0ab57f4654f3265308 Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Fri, 3 Nov 2023 11:40:44 +0100 Subject: [PATCH 7/8] corrected behaviour --- vermouth/processors/apply_rubber_band.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index 4f0f23ef4..f897aae86 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -94,9 +94,8 @@ def compute_force_constants(distance_matrix, lower_bound, upper_bound, np.fill_diagonal(constants, 0) constants *= base_constant constants[constants < minimum_force] = 0 - constants[constants > base_constant] = 0 + constants[constants > base_constant] = base_constant constants[distance_matrix > upper_bound] = 0 - constants[distance_matrix < lower_bound] = 0 return constants From a14cd235faef71f45295f9fae6e2769758d819ea Mon Sep 17 00:00:00 2001 From: csbrasnett Date: Mon, 6 Nov 2023 14:25:52 +0100 Subject: [PATCH 8/8] expanded docstring of compute_force_constants --- vermouth/processors/apply_rubber_band.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index f897aae86..100ffc134 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -89,6 +89,14 @@ 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)