2023-07-19 FF Fitting Meeting

I had a couple of short updates, one about redundant parameters in Sage and one about conformera with high energy sterics.

Lexie

• Slides:

• TG – 6X2 shouldn’t really exist

• CBy – argues for having generic parameters

  • Cby – there’s a strong history of pulling parameters out of the air. We could have a generic parameter in there as a stop gap measure, that we’ve trained to fit to all the training data

  • CBy – could run FF through huge database like eMolecules (stuff like enamine is too realistic) to evaluate coverage. If nothing matches, we could have a case for getting rid of t123 because we have nothing that matches it

  • PB – we do have some generic parameters

  • CBy – we don’t want to have a Jenga stack of a FF

  • PB – there is a molecule in PubChem that’s quite similar (although t123 wouldn’t actually apply due to double bond)

    • https://pubchem.ncbi.nlm.nih.gov/compound/102015867

  • CBy – could we have routine test on some large dataset to test parameter coverage

  • PB – yes, did suggest something like ChEMBL30

  • CBy – could look at all parameters that have low population

  • PB – in 2.1.0 I added in parameters to cover rare chemistries in Chembl and pubchem

• TG – we can go through parameters and sanitise things, e.g. torsion with N connected to two things.

  • PB – also, could start with k value of 0.1 and fitting the force field from there, so we don’t get ridiculous force constants

    • PB – could be equivalent of MSM for torsion parameters. Currently starting from previous fits for initial values

    • Would need something like simulated annealing for proper sampling

  • TG – did you make t123a PB?

  • PB – yes

  • TG/PB – changed periodicity from 1 to 3 in changing from t123 to t123a

• (conformer sterics)

  • CC – are you including the 1-4 LJ interactions?

    • General – this is probably included by default

    • PB – @Chapin Cavender Force groups in openmm
      'HarmonicAngleForce'
      'HarmonicBondForce'
      'NonbondedForce'
      'PeriodicTorsionForce'

      and Lexie is zeroing out the charges to make electrostatic contribution zero

  • CBy – atropisomerism is common in medicinal chemistry

    • e.g. https://en.wikipedia.org/wiki/Indometacin

    • The interactions that hinder rotation aren’t 1-4s, they are 1-5s, 1-6s

    • where the high steric interactions are a problem is where we parameterize torsions. The high sterics need to be screened out for the torsions, but it’s critical that the FF can account for these energies.

    • CBy likes the dexp form for softening the hard repulsions that lead to the high LJ (e.g. 12 kcal/mol). To get the minimum correct, we end up with the excessively high energies we see here

    • CBy – we need a way to encounter high steric interactions with pharmaceutical relevance, e.g. torsional transitions that we should get right

  • BS – we saw a lot of C=O interacting with S. It looked a lot like they were electrostatically attractive when the FF would make them repulsive. We found a lot of FF energies were messed up with these interactions

    • CBy – virtual sites would solve this – aromatic sulfur sigma holes are important here

  • CBy – you probably don’t want LJ interactions beyond 1-4, 1-5s for training torsions. The exception is if the SMIRKS specifies more – e.g. if there’s a methyl coming off an aromatic N, you know there could be key LJ 1-6+ interactions.

NEB update

Bill

• LPW – @Bill, one comment I have on your contour plots, (I really like them!), is that the TD path might not be perpendicular to the contours. I think maybe the contours should be perpendicular to the driving angle in the TD pathway, because that is the constrained DoF. Does that make sense?

  • I think maybe the contours should be vertical or parallel to the y-axis at the TD path. (Sorry again, I can't speak without waking people up.)

• CBy – use WBO for TD scan?

  • Potentially for every geometry with a WBO that changes significantly, we drive in the other direction instead. This also lets us avoid doing a full 2D scan

  • Would expect all conjugated single bonds to exhibit this behaviour

• LPW – I wonder if it is possible to compute a 2D PES by scanning both the driving angle and the pucker (which could be done with geometric or with TD)? This would depend on there not being other DoFs in the system that could flip.

  • CBy – one issue could be that the improper term could become dependent on the proper, and we may need cross-terms for this

  • LPW – I agree with Christopher that we would want to do it as a case study here, and potentially as something to guide improving the functional form, but not for all the training molecules because the cost would be too high.

  • CBy – BS/TG, could you do regular driving from the other direction?

    • TG – yes

    • CBy – do you compute WBOs?

    • TG – yes, they’re free

• TG – this initially started as a comparison between TD and NEB

  • I did AIMD on this molecule and I do see this bending at 300 K, so I know this is legitimate behaviour on part of TD, but have not been able to reproduce the bending with NEB, even with 500(?) images

  • In general I think TD is doing the right thing here, or at least sampling what we need

  • BS – what angles are you not sampling?

  • TG – if I only fit an FF on the NEB data, I wouldn’t capture any pucker

  • TG – current hypothesis is the geometric is the optimizer finding these images, not Psi4. Luckily geometric has NEB implemented so will see if it can find these puckers