[BUG] somd-freenrg is expanding the water box
nigel-palmer-cresset opened this issue · 9 comments
Describe the bug
I have been seeing the following error from analyse_freenrg.py mbar.
Traceback (most recent call last):
File "/data/tmp/nigel.palmer/ceb3/conan/p/sirebb1417ebb1a62/p/lib/python3.11/site-packages/sire/legacy/Tools/FreeEnergyAnalysis.py", line 136, in run_mbar
results = MBAR_obj.compute_free_energy_differences()
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/tmp/nigel.palmer/ceb3/conan/p/sirebb1417ebb1a62/p/lib/python3.11/site-packages/pymbar/mbar.py", line 698, in compute_free_energy_differences
Theta_ij = self._computeAsymptoticCovarianceMatrix(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/tmp/nigel.palmer/ceb3/conan/p/sirebb1417ebb1a62/p/lib/python3.11/site-packages/pymbar/mbar.py", line 1793, in _computeAsymptoticCovarianceMatrix
check_w_normalized(W, N_k)
File "/data/tmp/nigel.palmer/ceb3/conan/p/sirebb1417ebb1a62/p/lib/python3.11/site-packages/pymbar/utils.py", line 371, in check_w_normalized
raise ParameterError(
pymbar.utils.ParameterError: Warning: Should have \sum_n W_nk = 1. Actual column sum for state 0 was 10.971530. 11 other columns have similar problems.
This generally indicates the free energies are not converged.
While investigating it I have found that the size of the water box is significantly larger between the beginning and end of the somd-freenrg calculation.
In the image below the input waterbox is about 30A afterwards its about 350A.
Do you know what could be causing this?
I have included the input, output and config files in
somd-freenrg_30~ABFE_Free_Discharge_0.000_2630_1.zip
(please complete the following information):
- OS: AlmaLinux release 8.9 (Midnight Oncilla)
- Version of Python: 3.11
- Version of sire: 2024.1.0
- I confirm that I have checked this bug still exists in the latest released version of sire: no
The command was "bin/sire_python" --ignore-ipython "share/Sire/scripts/somd-freenrg.py" -C Config.cfg
Additional context
Add any other context about the problem here.
Just to confirm that I see the same thing locally with the latest development version of Sire. I also see the same thing for both PME and RF cutoffs. I have also run the same perturbation using somd2 and see no such explosion.
I've seen similar things so thought I'd jump in.
Box Expansion
I think the apparent expansion is due to the coordinates being unwrapped. You can wrap your coordinates with cpptraj - here's a wrap_coordinates.in script:
parm Input.prm7
trajin traj000000001.dcd
autoimage
trajout wrapped_trajectory.dcd
run
which can be run with cpptraj -i wrap_coordinates.in. When I rerun your input and use this script, wrapped_coordinates.dcd shows no apparent expansion.
PyMBAR Convergence
I assume you're using pymbar 4? I have a memory of seeing this error after upgrading from pymbar 3 to 4 on a perturbation which previously converged with pymbar 3, but didn't have time to dig into exactly why. Could you try downgrading to pymbar 3?
Thanks, @fjclark. I can confirm that things look fine once the trajectory is wrapped.
I've checked with sire 2023.1.3 and the wrapped trajectory is saved by default with the above input files.
Thanks. I think Sire has a wrap boolean option that can be passed via the map kwarg. I'll see if I can get it to image correctly using that. I'm assuming the default behaviour for the DCDWriter must have changed for 2024.1.0, i.e. it writes the raw, unimaged coordinates. (This makes sense if you want to analyse things like particle displacements, etc.)
The change in imaging is the result of this fix to OpenMM. Previously the MonteCarloBarostat would automatically image the molecules, which isn't desired if you wish to calculate displacements along a trajectory.
Downgrading pymbar to version 3 did fix the convergence error.
Many thanks for all your help.
No problem, thanks for confirming.