2022-11-23 BespokeFit Meeting notes

General updates

OpenMMForceFields issue

• DC + JH – DEXP fits were OK to opts, but bad for torsion scans. It turns out that torsion scans were being parameterized through OpenMMForceFields, which silently ignore custom nonbonded forces. I made some local changes to make custom offxmls loadable and fix the customnonbondedforce export and the torsiondrives+fits look good now. So OpenMMForceFields could be updated for this, or we could have QCEngine directly call the OpenFF toolkit.

• MT – OMMFFs is very powerful but should probably be temporary. Eventually we’d like for this to be handled by OpenFF’s infrastructure. But it will be a ton of work to replace OMMFFs, so it’s not likely to be replaced in the short term.

• MT – I’ve opened a PR to raise an exception if a customnonbondedforce comes out from openmmforcefields.

• MT – Hard to decide what to do with OMMFFs - We could commit to maintaining it but then we’d shoulder a lot. Chodera lab doesn’t really pay attention to it until things get pretty bad.

• DC – One question is how badly we need this given the situation where this came up.

• JH – So this didn’t come up until we tried to do fitting to torsiondrives. But for straightforward parameterization the toolkit should be fine, but for protein-ligand stuff we’ll probably need it. I could try putting a fix in QCEngine.

• DM – The QCEngine route may be the best thing here.

• JW – Nervous that we’d end up also maintaining QCEngine then too.

• JH – That’s pretty much already the case for the QCE OMM interface.

• MT – … Also changes will be needed for the new interchange plugin interface, so all of this is around 0.10.X OFFTK line. So that’s another dimension to keep track of.

• JH – True, and OMMFFs gives us access to GAFF, which is really useful for experimentation. So unless we get GAFF in SMIRNOFF format we’ll also want that.

• DM – Porting GAFF will be a ton of work.

• MT – Interchange may be able to help here - Would eventually let us import files directly from AMBER parameterization pipeline. Then we could make GAFF systems in AMBER and combine them.

• DM – I wonder how much we want to invest in a future where we continue comparing to GAFF - They don’t even give us adequate instruction on using GAFF correctly or version it clearly. Pharma partners aren’t asking us to do GAFF comparison, so maybe we don’t need this.

• DC – So QCEngine fix may be the way to go? JH, could you open a PR to QCE?

  • JH – Sure, will do.

• PB (chat) – I think we may need openmmforcefields to use/test espaloma as well

  • DC – I’m remembering the torsion scan problems with espaloma - Might these be related to this torsion scan problem?

  • JH – I can look into this.

Macrocycle fragmenting

• DC – I’m not convinced that BespokeFit was intended to treat macrocycles. It seems like a science project to decide what to do with them

• JW – Seems like there’s three options:

  • What we’re doing now

  • Skip making bespoke torsions if the fragment would involve macrocycle

  • Something clever where we cut up the macrocycle

• V – We noticed two things WRT macrocycles:

  • Fragmentation – We’re ending up with cases where fragments are larger than they should be. One part of this is keeping rotatable substitutents in the fragment. Another is some unexpected behavior from WBO fragmentation.

    • DC – Details?

    • V – Heuristic is “path length”. So the “WBO” setting grows the fragment until the WBO of the target bond stops changing. “Path length” counts the number of bonds away from the target bond.

    • V – Also, we use xTB to get WBOs in our code. Seems to give very similar results to AmberTools and runs much faster

    • DM – I think we’d seen this before when using xTB to get WBOs.

    • PB – Re: path length - I think Fragmenter fragments based on WBO, and the heuristic you’re talking about is how bespokefit generates the new torsion parameters?

      • V – The heuristic here may be “on what basis am I adding the substituent to grow the fragment”

      • JH – I think that’s right, it calculates the WBO for the proposed one, then chooses which path length to grow to.

    • PB – It looks like the problem in the GH issue may not have arised in the BespokeFit stack.

      • V – Right, so we may need to open issues on fragmenter and/or forceblaance

    •  

  • FB fitting

    • V – It seems like there’s an issue when there are a lot of flexible substituents hanging off, where the parameter fitting takes a long time to run. JW sent a reproducing example from his end and it ran in ~700 seconds for both of us. Then I tweaked my implementation/settings and got about 700 seconds from my workflow as well.

    • DC – This agrees with my intuition - FB will take a long time to work on big potential energy surfaces. So smaller fragments/fewer dofs makes a big difference.

    • V – Agree, it seems like we get better performance+parameters from reducing dofs/flexible substituents.

  • V – We have another non-macrocycle moleucle that we can’t share that takes a long time in ForceBalance, and so I need to check whether these changes will improve runtime for that. JW suggested that there may be representation issues with 5-membered aromatic rings but that didn’t seem to be a case.

• V – Another thing we hit was proton transfer. During some optimizations the connectivity of the input molecule would change.

  • JW – We’ve added changes to raise an exception when this happends https://github.com/openforcefield/openff-bespokefit/pull/193

    • Intend to make a release with current main branch, which includes https://github.com/openforcefield/openff-bespokefit/compare/0.1.2...main at commit https://github.com/openforcefield/openff-bespokefit/commit/3abdcca136cf47d33420081ea7c4e4fb5c3f052e

  • V – Right now I use GeomeTRIC to constrain NH bonds. This avoids the proton transfer issues, but I’m not sure how much it affects the results when there isn’t a transfer

  • JH + DC – It may be cleaner to not use points where proton transfers occurred.

  • V – This requires avoiding cases where we lose too many points, and that’s a science question. It also depends on whether the missing points are something like a minimum, in which case it’s really important to keep this around.

  • JW – I’d be in favor of “raise an exception if there’s a proton transfer” or “don’t make a new parameter”, followed by “constrain NH bonds”

  • V – One thing is of you put constraints to “auto”, then one opt has an NH bond go from 1 to 1.2, then the next goes from 1.2 to 1.4, etc, then you can eventually get a proton transfer. So that’s one thing to watch out for

  • DC – Yeah, that would also involve re-tuning the settings, which would be a lot of work.

  • JW – For now, we’ll push out the “raise an exception” change, and if we have extra engineer-time next year we can look at a better solution. Could you share code snippet for constaining NHs in torsiondrive+xtb?

  • V – I only have it for xtb, it’s a simple “iterate over atoms in molecule and find NH bond indices”

  • JH – Chapin’s done bond length constraints in torsiondrive, so he or I could provide code snippets/example inputs.