High force constant for bond [#6:1]=[#8X1+0,#8X2+1:2]

[pbuslaev]

I noticed that the force constant used in OpenFF 2.2.1 for [#6:1]=[#8X1+0,#8X2+1:2] bond is quite high. Compared to protein force field (Amber14 SB), where such bond is present in Asparagine, the force constant used for small molecule is 1.5 times higher. In some of my runs this leads to problems, when running with HMR and time step 4 ps. I am just wondering, why would the force constant for this double bond would be so much higher than for a similar double bond in protein force field? Is this something expected? What are the potential drawbacks of reducing this constant a bit, to gain simulation stability?

[lilyminium]

Hi Pavel, thanks for the report. Would you be able to provide a test HMR case for us to look at? Our bond parameters are currently only trained to gas-phase QM data, which mostly supports a force constant of this magnitude -- more below. While we have manually changed parameters in the past, we mostly aim to avoid that and let the data guide the parameters.

QM data

It looks like the force constant for this bond mostly comes from our small molecule dataset -- the values we get using the modified Seminario method range from ~1300 to ~1900 kcal/ (mol A2) across 50k molecules.

It's been suggested that some of the starting force constants we get from MSM are a bit high for some of our double/triple bonds, but from a preliminary look it doesn't seem that this is one of them -- in our most recent trial fits for the upcoming 2.3.0 force field, the force constant actually increases with training from ~1500 to ~1650.

A different force field fit to a very different dataset (SPICE2, which are single-points computed at a different level of theory and I believe includes amino acids) also sets this force constant quite high (https://github.com/openforcefield/yammbs-dataset-submission/blob/e498c5b6e2abeca7fcc817075a55cf570f519580/submissions/2025-08-13-smee-spice2-systematic-torsion-generation/smee-spice2-systematic-torsion-generation.offxml#L31).

[pbuslaev]

Hi @lilyminium. Unfortunately, I can't provide the molecule. But as I mentioned in my initial message, my main concern is that this constant is significantly higher than similar constant in AMBER FF (~1100 kcal/mol), which is out of the range 1300-1900, that you are observing. I was somehow under impression, that for proteins force constants, especially for bonded terms, should be highly optimized, and I was a bit puzzled with such a big difference.

The issue with HMR, is in a way quite clear: for numerical integrators to have reasonable accuracy, at least 5 numerical integration steps should occur during the bond oscillation period. Given the force constants and atomic masses (especially when the carbon mass can get reduced during HMR), it is clear that with force constant used and time-step of 4 fs, we are getting very close to the cutoff of 5 integration steps (essentially 4.9 with carbon mass reduced due to HMR, and 5.1 with normal carbon mass). This can lead to some instabilities, which I observe as increased simulation failure rates with both GROMACS and OpenMM. Also, in case of FEP for the changes that include such bonded term, I normally observe higher errors. I will need to run with smaller time-step to double check that the problem is due to integration inaccuracies, but this is my current working hypothesis.

[lilyminium]

Hi Pavel, at present Sage 2.2.1 is only a small molecule force field and fit to small molecule data -- are you using it on a protein or peptide system? We are in the process of fitting a protein force field at the moment to additional protein data, so the force constants there may change -- either way, our philosophy is to let the data guide the parameters.

Please do keep us updated on your findings, we'd be interested to see if it is indeed due to integration issues!

[pbuslaev]

Hi, maybe I was a bit unclear in my message. What I do, I use Sage 2.2.1 for ligands and Amber14SB for protein. I observed some higher failure rates for some of the ligands, and the common basis for all those ligand was [#6:1]=[#8X1+0,#8X2+1:2] bond. I then checked the predicted oscillation period for this bond given the force field and noticed that it was at the edge of what numerical integrator support. Then, I checked what is the force constant for the same group in established protein force field (Amber14SB in my case) and it was significantly lower, so I raised the question.

Hope that my story is clearer now. And essentially, I don't say that there is a problem or something. I am just confused by the fact that there is such a difference between what we have in protein for this bond, and what we have in Sage. I believe, the overall approach for predicting the force constant is more or less the same for fitting Amber14SB parameters and Sage parameters, but still the 1.5 factor difference leaves me puzzled.