[WIP] Custom GBSA

[j-wags]

Description

Adding GBSAHandler(s) that allow varying the alpha, beta, gamma, and possibly kappa globals in OBC models. Possibly also adding support for the GBn1+2 models.

Todos

• Get OBC1 reference energy tests matching
  • Document whether the radius particle parameter SHOULD or SHOULD NOT have the offset (0.009 nm) applied in each method
• Get OBC2 reference energy tests matching
  • In the OBC2 OFF vs. AMBER comparison tests, we're often getting an absolute GBSA energy discrepancy of 1e-5 (rel error ~1e-6 to 1e-7). Is this acceptable? Could be due to sig figs in conversion factors? Is it significant?
• Get some test with nonzero kappa matching known references (?)
• IF we allow changing the surface_area_penalty, add tests matching to a known model with a different SA term.
• Add support for GBn and GBn2 models

Questions

• Should we allow people to change the offset (0.009 nm)?
  • @jeff231li: yes, the offset is exposed in the offxml file

Status

• Ready to go

[codecov-commenter]

Codecov Report

Base: 79.27% // Head: 62.88% // Decreases project coverage by -16.39% ⚠️

Coverage data is based on head (208a741) compared to base (6118a88).
Patch coverage: 32.07% of modified lines in pull request are covered.

Additional details and impacted files

@@             Coverage Diff             @@
##             main      #11       +/-   ##
===========================================
- Coverage   79.27%   62.88%   -16.40%     
===========================================
  Files           5        6        +1     
  Lines         275      458      +183     
===========================================
+ Hits          218      288       +70     
- Misses         57      170      +113     

Flag	Coverage Δ
unittests	62.88% <32.07%> (-16.40%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
smirnoff_plugins/handlers/customgbsa.py	31.41% <31.41%> (ø)
smirnoff_plugins/handlers/nonbonded.py	93.92% <66.66%> (-0.99%)	⬇️
smirnoff_plugins/utilities/openmm.py	46.93% <0.00%> (-1.49%)	⬇️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[lgtm-com]

This pull request introduces 5 alerts when merging 966ec41 into 7131a2d - view on LGTM.com

new alerts:

• 5 for Unused import

[lgtm-com]

This pull request introduces 5 alerts when merging d3b2ea0 into 7131a2d - view on LGTM.com

new alerts:

• 5 for Unused import

[jeff231li]

@j-wags I ended up implementing a handler class that can assign parameters for all Amber-based implicit solvent models in OpenMM (HCT, OBC1-2, GBn1-2). All GB models use CustomGBForce even for OBC2 (the Toolkit uses GBSAOBCForce). However, for GBn2, I didn't implement the radius correction for nucleic acid. I tried to replace all of the hard-coded parameters but I still don't know what unit for 138.9354825 is in the energy expression. The parameters below are now exposed as attributes in the handler.

• alpha, beta, gamma - effective radius parameters
• neck_scale and neck_cutoff - for GBn and GBn2 models
• solvent_radius and surface_area_energy - for surface area calculation
• kappa - Debye-Huckel screening parameter

The regression tests work for all GB models + FreeSolv molecules but for OBC2 there were 5 molecules that failed. They pass when I increase the tolerance to 2.0e-4.

Question:

• Do you think the energy discrepancy is a significant issue to warrant further investigation?
• Should we update the current GBSAHandler class in the Toolkit and increase the GBSA spec to version 0.4 or leave it as a plugin for now?