Use OOP forcefields

[a-ws-m]

This implements several changes to make it easier for end-users to implement their own forcefields. Most of these changes revolve around the Forcefield class, which stores the path dictionaries that were previously separate entities. To use a non-default forcefield, make a new Forcefield instance and then pass it to Simulation using the ff_class argument. This argument is optional; passing a string value for forcefield will still work, it's simply converted behind the scenes. I also tidied up a few things, which may be worth discussing:

• The system.top and system_octwet.top files have been changed to be compatible with Python's built-in Template functionality. This makes the code a bit easier to understand.
• Water models are now attributes of the individual forcefields, and are read directly from the corresponding watermodels.dat file. As I understand it, GROMACS_WATER_MODELS is a master list of all water models, some of which are only available for OPLS-AA. If this is correct, it might be good to phase this out. I have changed the core functions that used this: get_solvent_model and get_solvent_identifier. The only other places in the code where GROMACS_WATER_MODELS is used is in the tests and in get_water_model, a function which is also unused in the code base except in tests.
• I have changed one test to check for equality rather than equivalence (is -> ==).

I'll write some example code that demonstrates how to use a Forcefield next week. In the meantime, I'm sharing this PR so that the changes can be discussed.

[a-ws-m]

@orbeckst here's the PR you asked me to ping you about!

[codecov]

Codecov Report

Merging #267 (3645776) into develop (f98424b) will increase coverage by 0.42%.
The diff coverage is 93.54%.

@@             Coverage Diff             @@
##           develop     #267      +/-   ##
===========================================
+ Coverage    81.26%   81.68%   +0.42%     
===========================================
  Files           15       15              
  Lines         1975     2026      +51     
  Branches       297      310      +13     
===========================================
+ Hits          1605     1655      +50     
+ Misses         278      276       -2     
- Partials        92       95       +3     

Files Changed	Coverage Δ
mdpow/equil.py	80.79% <88.23%> (+0.11%)	⬆️
mdpow/forcefields.py	92.56% <94.73%> (+4.75%)	⬆️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[orbeckst]

Thank you for the PR. In general I like your approach

• consolidate in Forcefield class
• make the templates real templates

I have two initial requests:

1. The original intention with the templates was that they should in principle be legal input form GROMACS. However, I like explicit templating better (also, with 10+ experience since I wrote the original code...) so could you please rename the files from *.top to *.top.template (or some other appropriate suffix... *.top.in?) to make clear that these are now template files?

2. You seem to be using an autoformatter which makes reading the diff difficult. Can you please either undo the formatting or submit a separate PR only with the formatting changes on the files that you touched. We can then merge that PR and then we can merge the main into this PR (or you can rebase) and we have a clean diff.

[a-ws-m]

Done! I always forget to disable black when I'm working on new code, sorry about that.

[orbeckst]

You gave me the motivation to try out black on a project so I opened #271. Let's do PR #269 first and then come back here.

[orbeckst]

@a-ws-m now that PR #269 is merged, could you please rebase against current develop so that we hopefully end up with a nice, clean diff? Thanks!

[orbeckst]

Thanks for the update.

There are failing tests. Please make changes so that the tests pass again.

You can also run the tests locally (which is faster) if you have GROMACS installed (there are conda packages).

[a-ws-m]

I've made a few important changes:

• The solvation errors arose because Gromacs expects a .top file suffix, but the copies of the templates were saved with identical names (i.e., including .template). I have added a few lines to ensure that the topology file does have the correct suffix.
• forcefields.get_solvent_identifier() is now consistent with get_solvent_model() and raises an error when no such identifier exists for the given force field.
• I removed the python.system_packages key from the read the docs configuration as it was just deprecated and was causing build errors.

[a-ws-m]

Regarding the error in the FEP tests, I'm unable to reproduce it on my machine and I think it's an issue with the MD workflow rather than any Python code that I've changed here.

[a-ws-m]

I've made those changes. The GitHub tests seem to be throwing some weird errors to do with cloning the repo.

[orbeckst]

I restarted CI and all jobs are running. It was a temporary problem.

[orbeckst]

This is all looking really good and I only have very few final requests:

1. I'd like to unify the API a bit by only using forcefield as a kwarg instead of ff_class and forcefield (see comment).

2. There are some other minor things in the comments.

3. Please update the docs at the top of forcefields.py where it says

   In the principle it is possible to switch to a different force field by supplying alternative template files.

   The idea with this PR is that it now becomes easier to add custom force fields. Can you show an example? You could add a section "Adding your own forcefield". That would really help (1) for maintenance because it tells developers how it works and (2) for users so that we don't need to answer questions but can just point to docs.

Sorry it's taking a while. It's just that this is a sizable (and very good!!) change and I want to make sure that everything continues to work well together.

[a-ws-m]

No problem! I'll work on those. I have actually been working on an example, it's in the martini-example branch of my fork. I have some questions before it's ready that I'll raise on the Q&A.

[orbeckst]

awesome!

[a-ws-m]

I've merged the example using Martini, you can find it in doc/examples/martini-example. I imagine there's scope for improving the run parameters to get a more accurate value, and perhaps fixing the issue whereby the nan error for the coulomb contribution propagates and gives a nan value for the overall error. I'm not sure how to improve these, though -- I'd be happy to hear suggestions!

I also need to go back and clear some of the old, erroneous output from that notebook. The plots, for example, are completely wrong and don't reflect the latest results.

[orbeckst]

Your MARTINI/MDPOW example notebook is fantastic!

(I also noted that some of the GROMACS runs actually didn't run but just restarted from the last frame. However, it's nice to know that it's not an issue re-running the workflow as it won't destroy data.)

For nan value in the Coulomb calculation I am not sure where it comes from. I assume that there is a valid Coulomb leg in the FEP calculation because MARTINI benzene has partial charges, doesn't it?

One thing that's a bit odd is that alchemlyb and mdpow report different values for statistical inefficiency, for example alchemlyb states 15.84 and mdpow states 7.9840.

mdpow.fep   : INFO     Performing statistical inefficiency analysis for window coulomb 0000
2023-09-04 22:41:16.809 | DEBUG    | alchemlyb.preprocessing.subsampling:statistical_inefficiency:520 - Running statistical inefficiency analysis.
2023-09-04 22:41:16.812 | DEBUG    | alchemlyb.preprocessing.subsampling:statistical_inefficiency:522 - Statistical inefficiency: 15.84.
2023-09-04 22:41:16.813 | DEBUG    | alchemlyb.preprocessing.subsampling:statistical_inefficiency:528 - Number of uncorrelated samples: 312.
mdpow.fep   : INFO     The statistical inefficiency value is 7.9840.
mdpow.fep   : INFO     The data are subsampled every 8 frames.

I think that's because we have to infer the actual value from what alchemlyb subsample in
https://github.com/Becksteinlab/MDPOW/blob/6449879a9bc5fbc80210290eea84fd40ad164dff/mdpow/fep.py#L1192C20-L1192C20
and there's also a suspicious looking factor of 2 involved. (The whole question about getting the raw SI from alchemlyb was the topic of alchemistry/alchemlyb#295 ). @VOD555 can you please check that the MDPOW output is consistent with the one from alchemlyb? I assume that other users that read the log output will have the same question as I

[VOD555]

@orbeckst I've checked the code. The logger in the mdpow outputs the correlation time (g/2) but not the statistical inefficiency (g). And as the one calculated in mdpow is based on the length output of the subsampled dataframe, there is small difference between the mdpow one and alchemlyb one.

As currently alchemlyb has a logger output for the SI, I think we can remove the mdpow logger information.

[VOD555]

Open an issue #275

[a-ws-m]

Thank you! The Gromacs output is for a restart because it was the quickest way to re-generate most of the useful information after clearing the cell outputs. I'll work on those other points today.

[orbeckst]

I am very happy with your contribution. It's a really neat re-organization that's also backwards compatible.

Your example notebook is also really good. Nice work!!

[orbeckst]

@a-ws-m I let you look at my one comment about "sedstate" in the MDP file but if you tell me that this is what it ought to be, I'll squash-merge. (If you want rewrite the history on the PR into less than 5 commits then you're also welcome to do so, e.g. code changes+tests, example, but if not then that's fine too.)

[orbeckst]

Thank you very much @a-ws-m , great contribution!!