ABFE restraints torsional force constants are too high

[fjclark]

Describe the bug

This is a fairly minor, but I've found that my restraints fitting code is selecting torsional force constants ~ 3 times stronger than is consistent with normal fluctuations in the unrestrained state. This does not affect the correctness of any ABFE calculations (e.g. the restraint corrections are still correct and lower torsional force constants will not affect stability).

I came across this when I realised that I could test the restraint search code by feeding it a trajectory of a completely decoupled state (with some Boresch restraints applied), and the restraint code should recover exactly the restraints used to simulate the decoupled state. I found that this worked, other than the force constants being ~ 3 times too large. The issue is here:

                        for val in boresch_dof_data[pair][dof]["values"]:
                            dtheta = abs(val - circmean)
                            corrected_values.append(min(dtheta, 2 * _np.pi - dtheta))
                        corrected_values = _np.array(corrected_values)
                        boresch_dof_data[pair][dof]["var"] = corrected_values.var()

Here, I'm shifting the deviations between the measured torsions and the mean value so that they are in the range [-pi, pi], consistent with their periodicity. However, I've forgotten that because I've taken the absolute value, the deviations are all on the same side of the mean, and so corrected_values.var() on the last line is incorrect. Instead, it should be np.mean(corrected_values**2), which avoids erroneously re-centring through subtracting the mean of the deviations.

Assuming that the deviations in the unrestrained state used to fit the restraints are normal (which has been a pretty good approximation in my experience, at least for well-behaved systems), then this bug produces an overestimate of the variance by a factor of ~3:

>>> data = np.random.normal(size=1_000_000, scale=1)
>>> 1/np.var(abs(data))
2.75

(whereas 1/np.mean(abs(data)**2) = 1)

Boresch restraints are subject to instabilities if the angular force constants are too low, but lowering the torsional force constants should not affect stability.

I'll fix this along with the other discussed tweaks to the restraint selection code.