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

-el fix #559

merged 8 commits into from
Nov 13, 2023

Conversation

csbrasnett
Copy link
Collaborator

Currently, if elastic network rubber bands are applied and are 1) below the cutoff for -el and 2) not otherwise eliminated by the -ermd threshold for the force field, they are kept in. Ie. the -el cutoff is not actually currently applied. This PR applies the cutoff.

@fgrunewald
Copy link
Member

@csbrasnett do you dare to add a test for this behaviour? Editing 'tests/test_apply_rubber_bands' would be the place. Probably would make sense to add a new test for the compute_force_constants function. To get started you can ask chatGPT to write a test for example.

@csbrasnett
Copy link
Collaborator Author

csbrasnett commented Nov 1, 2023

@fgrunewald I was worried you'd ask... Have had a go now, that'll cover the testing of the different bounds. I haven't added anything to test the decay factors/powers because that gets in to testing floats and I didn't want to get in to that/ the point of this is to test that the cutoffs are actually applied, rather than the maths working correctly!

Copy link
Member

@fgrunewald fgrunewald left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For extendability purposes I suggest to use np.allclose and then you could also add a test for the decay factor but totally optional here.

minimum_force)

# Assert the result against the expected output
assert result.all() == expected_output.all()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you want to test floats you can use np.allclose(array1, array2) together with setting the absolute tolerance for the float. Or Peter seems to prefer the pytest build in using pytest.approx.

@pckroon
Copy link
Member

pckroon commented Nov 2, 2023

I'm not sure the behaviour you propose here is desired actually. The semantics between upper and lower bounds are annoyingly different in this context. The code you propose to change has been there since #54, the initial implementation for the EN.

If this change goes through we will produce a different elastic network than martinize1

@fgrunewald
Copy link
Member

@csbrasnett did you compare this to martinize1? Also @pckroon I don't really care if we diverge here. The implementation in martinize1 seems just incorrect or am I missing something?

@csbrasnett
Copy link
Collaborator Author

I haven't checked against martinize1, I can try and find time today.

But actually @pckroon is correct in that if the distance is < lower then it needs to just be set as the constant, rather than being deleted as this currently behaves. Which still wasn't happening before - eg. if you applied a decay and a lower cutoff, then the forces below the cutoff would be larger than the given constant. But that's easy enough to change. I can fix these and check against martinize1 if we agree that's the expected behaviour?

@pckroon
Copy link
Member

pckroon commented Nov 2, 2023

@csbrasnett Could you chase down the desired behaviour in the original EN paper? I'm not sure which one it is, either martini2 protein or elendyn I guess.

@csbrasnett
Copy link
Collaborator Author

Right, so this would change the behaviour against martinize1, which currently reproduces what the main branch does. However, I 1) cannot find where the concept of a lower cutoff has emerged from, as in both the papers and the citation trail that follows, only the upper cutoff is ever discussed.

That said, I still think that the expected behaviour is currently wrong based on the descriptions. This PR now actually only fixes the side case where decays are used. In the fig here, where ea = ep = 0, the force is constant below el for both old and new. But with decays (and with line 99 now reading constants[distances < lower_bound] = base_constant as it should) the force below the el is correctly set to the given force constant.

fig

pckroon
pckroon previously approved these changes Nov 2, 2023
Copy link
Member

@pckroon pckroon left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Excellent. Thanks for digging up the details!

@csbrasnett
Copy link
Collaborator Author

Thanks! nb don't merge yet, need to add the final condition, which isn't as trivial as I thought it'd be.

Copy link
Member

@fgrunewald fgrunewald left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@csbrasnett let us know once you're done.

@csbrasnett
Copy link
Collaborator Author

@fgrunewald now corrected. Based on the discussion and graph above, I realised the second condition wasn't actually necessary, because the deviation is constant > base_constant for distance < el, so it's simple enough to set anything larger as the base_constant.

fgrunewald
fgrunewald previously approved these changes Nov 3, 2023
@fgrunewald fgrunewald requested a review from pckroon November 3, 2023 14:47
pckroon
pckroon previously approved these changes Nov 6, 2023
Copy link
Member

@pckroon pckroon left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd be very happy if you could expand the docstring of compute_force_constants now that you know exactly what it does and is supposed to do.
LGTM otherwise, nice work :)

@csbrasnett csbrasnett dismissed stale reviews from pckroon and fgrunewald via a14cd23 November 6, 2023 13:26
@fgrunewald fgrunewald merged commit 899a8df into marrink-lab:master Nov 13, 2023
8 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants