Definition of force constants inconsistent between Boresch restraints correction and applied restraints
Opened this issue · 2 comments
Describe the bug
This bug was my mistake and I'll fix it once we agree on the best approach.
The issue is an inconsistency between the definitions of the force constants in the Boresch restraints Python class/ standard state correction, where I've assumed that the definition is consistent with E = 0.5 * k * x**2, and the implementation, where the ks apply as E = k * x**2. The definition as E = k * x**2 is consistent with other restraints implemented in Sire e.g. here.
One fix would be simply to fix the standard state correction I implemented and keep the current definition consistent with E = k * x**2. However, in this case the ks are no longer force constants (because F = 2kx), so I should be careful to emphasise this in the documentation. The downside here is that most Boresch restraints are specified in the literature according to the E = 0.5 * k * x**2 definition, so we'd have to be careful to make this clear. The plus is that the force constants would be consistent with SOMD1.
Alternatively, the OpenMM restraints could be created using the E = 0.5 * k * x**2 definition. Then the ks would correspond to force constants and the Boresch restraints could be supplied in a way which is consistent with most of the literature (but not SOMD1). However, the implementations of other restraints should be changed to be consistent with this, which could be a pain for people currently using the other restraints.
Please let me know what would be best, and I'll make the changes.
Thanks.
I don't have a strong opinion either way. The logic of using k rather than 0.5 k is because it is easier to be internally consistent across all restraints, i.e. we can say for all restraints that k is the scaling constant, and it is up to the user to map the higher level equation (e.g. 0.5 k) to the value of k used at the lower level. This way, they don't need to read the code to find out what factor has been used to scale k - it is always unscaled.
I believe that is also how we handle force constants for e.g. bond and angle terms too. It makes things easier not having to remember to multiply or divide by 0.5 within the code.
Great, thanks for the comment - I'll just update the restraints correction.