Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

-el fix #559

Merged
merged 8 commits into from
Nov 13, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions vermouth/processors/apply_rubber_band.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
47 changes: 47 additions & 0 deletions vermouth/tests/test_apply_rubber_band.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)