2021-10-14 Force Field Release meeting notes

Smirnoff plugin update

@Joshua Horton

• JH will post slides here

• Update on fitting force field parameters for alternative non-bonded functional forms.

• https://docs.google.com/presentation/d/161rEOSXGJY9xV0UHqS-W33M-zpVLfJENZLCvvlfN0xA/edit?usp=sharing

• VIDEO LINK: https://drive.google.com/file/d/1koAIGeeHqRAVNY1uYBg7KPVFWPM4V02B/view?usp=sharing

• JH – We can’t do free energies with soft-core potentials, can only use YANK. We’ll need infrastructure for this.

• JH – 6-8 potential coded as a customnonbondedforce is much slower than standard nonbondedforce. Double exponential is just about as fast, though.

• JH – Compared our B68 implementation to Rowley’s previous work - We found some sizable discrepancies so we’ll need to talk to him.

• DM: Double exponential results?

• JH: I had to refit tip4p for this, and got results close to TIP4p-DE-FB. Although the parameters alpha, beta are different the final objective is accurate enough.
  It’s probably because of a noisy PES.

• JH: Next steps, keeping the alpha, beta constant we will refit sage epsilon, etc., with physical property data.

• DM – Looking at final slide, it’s not clear that changing the functional form should lead to dramatically better/worse performance, and the plots show about the same

• JH –

• JW: You mentioned some infra needs to use other than yank.

• JH: Yeah, for free energies

• DM: Openmm doesn’t offer free energy calcs directly, so Yank is a wrapper. Could talk to Perses devs - It’s focused on nonequilibrium FE calcs / lambda scaling - I think Perses' lambda scaling doesn’t care about using an LJ potential specifically, so a double exponential could be easily supported.

• SB – I could see YANK being flexible enough to do something other than LJ. But Perses is possibly more likely to have the flexibility that we need in the future. Alternatively, we could “brute force” it, and just do the free energy calc ourselves, scaling the nonbonded energies in a few steps.

• DM – Agree, this could be done manually.

• DC – This problem is on my roadmap as well, and I’m willing to contribute resources/person-time to fixing it.

  • DM – Should discuss how to interact with Perses/YANK devs. But in the short term I’d recommend doing local/manual changes to get the needed numbers

• DM – Does GROMACS support double-exponential?

  • JH – Not sure, but I do recall that they implemented a new sort of long-range correction that may have been driven by the need to use double-exponential potentials.

  • DM: GROMACS only has LJ 6-12 or Buckingham exp-6. (There is support for others via tabulated nonbonded potentials, but that’s more complex/doesn’t necessarily work with free energy infra.

• DC – In final plots, are valence terms refit?

  • JH – No

  • SB – CBayly would be asking whether the new potential form would give better accuracy for the steep short-range “wall”, and thereby improve the torsions.

  • SB – Generally when fitting to physical parameters, we’d not want to let forcebalance change torsions. So we’d need to set dtorsion_parameter/dphysical_observable=0 , but still let QM data affect the torsion parameters

  • JH – Right, but at this point we haven’t done any torsion fitting with the double-exponential.

  • DC – Will we be able to fit alpha/beta?

  • SB – Yes, can implement this when needed. It’s not currently there but won’t be hard to add.

  • DC – I’m hoping that this can help during mixed solvent fits.

  • DM – The multiple-optima-for-water-models point would be great to publicly state at some point - A lot of the field assumes that there’s a single optimum so this could spark interesting discussion

  • DC – Agree

  • SB – It would be timely to revisit the question of “what makes a good water model?” - We currently show the “big 6” series of plots on bulk water, but I wonder if there are new targets that we should consider that include water-nonwater interactions, and treat the “big 6” as a prior.

  • DC – I wonder about things like the dielectric constant - We’re not modeling electrons explicitly so why would we expect to get that right?

  • DM – Agree, ??? has written that we should never be able to reproduce the dipole of cyclohexane using MM, since the assumption of fixed partial charges will never capture the induced-dipole induced-dipole interactions that would be the major contributor.

  • JH – (Something about protein properties vs pure water properties)

• DC – ETA for first fit?

  • SB – Turnaround for these sorts phys prop fits is 4-7 days on lilac. Turnaround for valence props is about another 7 days on lilac. We’re losing some time on the first iteration since this is the first time we’re doing vsite fits but we should only have to fix these once.

• DM – This is very exciting, and it’s great that we can move forward quickly to do cool things and not stop to publish. But let’s make sure to write about the multiple-minima thing in some public form.

• JW – Timeline for double-exponential getting into smirnoff spec?

  • SB – I’d like to see the results of the feasibility study before we put this in the spec/toolkit

  • DM – Agree

• DM – TGokey’s work on alkanes has revealed two new parameters that may significantly improve FF accuracy. So that may come up in a future meeting. One is a C-C-H angle, new angle after a2, "[#6X4:1]-[#6X4:2]-[#1:3]"; new bond "[#6X4;r4,r5:1]-[#6X4:2]"

• SB – General issue: To make rapid improvement in FF iterations we need the benchmark refactor. Until then I’m not sure how to decide whether to accept these proposed updates.

• JW: Biopolymer refactor is the first thing I handle, depending on availability I can see how benchmark refactor fits in.