Possible issue with angle parameter involving five-membered rings

[hjuinj]

Hi @bannanc, a colleague of mine has noticed a particular pattern seemed to be causing strains in some systems he was looking at.

The SMIRKS in questions is <Angle smirks="[*:1]~;!@[*;r5:2]~;@[*;r5:3]" angle="125." id="a14" k="140.0"/> in the latest smirnoff99Frosst.offxml. He shared with me calculations from cccbdb.nist.gov, for example in the image I show dft with cc-pVTZ calculation of cyclopentane which has the relevant angles at equilibrium values of ~109 degrees rather than 125 (similar angles when using other method and basis sets). Is it possible that there might be an issue with this parameter?

[bannanc]

Hi @hjuinj

I wanted to post a quick update on what I have found looking into this parameter. I would like to remind you that smirnoff99Frosst is a translation of the parm99/parm@frosst force field into the SMIRNOFF format. That means this is intended to represent an angle parameter taken from that set of parameters.

That being said, I this it is possible we had a transcription error because this parameter is assigned to pentane type rings (those with tetrahedral carbons) and pentene/aromtic-ish type rings (those with trivalent carbons). It doesn't make sense to me that these groups would receive the same parameters, so you are probably correct that there is a problem. This parameter would make more sense in the latter group where the extra-cyclic angles should be greater than 120 due to the ring strain.

Here are some example angles highlighted by this parameter from the MiniDrugBank set with the corresponding parm99/parm@frosst atom types:

Tetrahedral carbon

CT~CT~H1

Planar carbons:

CW~CC~H4

NA~CR~H5

[bannanc]

In this comment I have my excavation into the original parameter files:

Atom types in the examples above:

Here are the specific parameters for the examples I shared above from parm99.dat.

H5-CR-NA    50.0      120.00    AA his
CW-CC-NA    70.0      120.00    AA his
CT-CT-H1    50.0      109.50    changed based on NMA nmodes

Similar parameters in parm99

When those failed I looked for a carbon type specific for 5-membered rings. There are a lot, mostly focused on carbon in aromatic heterocycles. So, my next step is to find the angle="125." k="140.0" parameter in parm99. Keeping in mind we intentionally rounded some parameters so I considered those as well. I tried to separate these by the chemistry in the group. Keep in mind that SMIRNOFF ks are 2* amber ks.

Angles between the oxygen and other neighbors in carbonyl groups:

CB-C -O     80.0      128.80
CM-C -O     80.0      125.30
O -C -OS    80.0      125.00    Junmei et al, 1999
O2-C -O2    80.0      126.00    AA GLU            (SCH JPC 79,2379)

Pyrimidine carbons (these are 6-membered rings so they don't apply here):

CA-CM-HA    50.0      123.30
CA-CM-H4    50.0      123.30
CM-CM-OS    80.0      125.00    Junmei et al, 1999

Bridging carbons between 5&6 membered rings (like tryptophan)
We have a separate parameter for these a13

N*-CB-NC    70.0      126.20
CB-CB-NC    70.0      127.70
CA-CN-CB    63.0      122.70    changed from 85.0  bsd on C6H6 nmodes; AA trp

The 5-membered ring in purine type groups.

H5-CK-N*    50.0      123.05
H5-CK-NB    50.0      123.05

Perhaps this is the parameter we're looking for, this is "aromatic" 5-membered rings with a substituent. aromatic is in quotes since these molecules won't be aromatic in the MDL model.

CT-C*-CW    70.0      125.00    AA trp
CB-C*-CT    70.0      128.60    AA trp

Similar parameters in parm@Frosst:

There are a significant number of parameters added with the parm@frosst extension which include the angle="125." k="140.0" parameter. I assume of these types is what we were trying to capture with a14. Here are those categories:

Addition of more sulfur parameters, I assume these have nothing to do with the ring size:

O -C -S      70.000 124.000 thioester thet0 from b3pw91/6-311Gdp 11/00
CB-C*-SO     70.000 125.000 guess from C*-C*-CT
NB-C*-S      70.000 125.000 RHF/6-31G(d) Mar 7 2007
H4-CA-S      70.000 122.300 calc RHF/6-31G(d) (estimated force constant)'
CB-CC-SO     70.000 129.000 RHF/6-31G(d) opt June 29 2007
CT-CR-S      70.000 122.000 calc RHF/6-31G(d) (estimated force constant)
CC-CW-S      70.000 125.000
CC-CW-SO     70.000 125.000

Specific biphenyl CP parameters (these are for bridged 6-membered aromatic rings so not relevant here):

CA-CP-CA     63.000 123.000 ff94 CA-CA-CA ; for within the ring
CC-CP-CP     70.000 125.000 ff94 CT-C*-CW; AA trp

More 5&6 ring junctions:

CA-CB-CB     63.000 122.700 ff94 CA-CN-CB; AA trp
CA-CB-CP     63.000 122.700 ff94 CA-CN-CB; AA trp

Extension of the CM atom type for pyrimidines:

CL-CM-CM     70.000 123.000 fit to AM1 opt july 10 2003
CM-CM-F      70.000 123.000 fit to AM1 opt July 10 2003
N -CM-NB     70.000 124.200 from RHF/6-31G(d) opt estimated force constant april 11 2000

"aromatic" 5-membered rings such as histidine:

CA-CC-O2     70.000 124.300 calc RHF/6-31G(d) (estimated force constant)
CT-CC-CW     70.000 125.000 ff94 CT-C*-CW; AA trp
NB-CC-O2     70.000 127.800 calc RHF/6-31G(d) (estimated force constant)

and

BR-CW-CC     70.000 125.000
CC-CW-CR     70.000 125.000
CC-CW-CW     70.000 125.000
CC-CW-N*     70.000 125.000
CC-CW-NB     70.000 125.000
CC-CW-OS     70.000 125.000
CT-CW-N*     70.000 125.000 ff94 CT-C*-CW

Conclusions:

It seems clear this parameter is being assigned too generally. Perhaps there was a mis-ordering that lead to it being placed too far down on the hierarchy. I assume this parameter should only be applied to the last category, or at least only trivalent carbons so make the change from smirks="[*:1]~;!@[*;r5:2]~;@[*;r5:3]" to smirks="[*:1]~;!@[#6X3;r5:2]~;@[*;r5:3]". These aren't applied to non-carbon atoms since they are above all the hetro-atom parameters.

All that being said, I'd like to get @cbayly13's eye's on this logic to see if I'm understanding the reasoning behind this parameter.

[hjuinj]

Hi Caitlin, thanks for the swift and thorough investigation. I too would very much like to hear Christopher's take on this.

I will just briefly add that I did look at the ordering of this SMIRKS. It has moved from a15 in the pre-offxml era version (which includes the omini-matching [*:1]~[*:2]~[*:3]) up to a12 in the latest version. For the extra-cyclic angle in cyclopentane, the only prior angle matching is to [*:1]~[#6X4:2]-[*:3], which currently is a1.

[cbayly13]

I agree with the analysis of @davidlmobley and @bannanc above. The equilibrium angle of ~125 is for external bond angles around an sp2 atom in a 5-membered ring; as @hjuinj noted this angle is wrong for tetrahedral centers in a 5-membered ring. Perhaps a minimal change to the offending parameter could be to introduce a "trivalent" requirement along the lines of <Angle smirks="[*:1]~;!@[*X3;r5:2]~;@[*;r5:3]" angle="125." id="a14" k="140.0"/> (didn't check syntax); this would address angles around carbon but it could still leave an issue with aliphatic nitrogen in a 5-membered ring (as in pyrrolidine). How shall we distinguish a methyl pyrrole angle from a methyl pyrrolidine? This is a classic "planarity of the nitrogen" question. When we can examine the bond order of the bond inside the ring it will be helpful.

[bannanc]

@cbayly13 I'm happy with using smirks="[*:1]~;!@[*X3;r5:2]~;@[*;r5:3]" but if we want it to apply to nitrogens as well as carbons then we need to move it down in the hierarchy, currently this parameter comes before ALL of the nitrogen centered angle parameters. Do you think we should be reordering?

I can put in a pull request in a little bit with the X3 added to this parameter.

[bannanc]

Also, I will double check the syntax matches the atoms we expect, but that SMIRKS looks good to me.

[cbayly13]

@bannanc if that change fixes only carbons it will be sufficient to address most of the concerns raised by @hjuinj, @davidlmobley and you above. For nitrogens I am still hoping we can address the planarity through the OOP as we discussed a while back (summer?), ie valence angle always wants tetrahedral but OOP forces planarity more or less depending on chemistry, so they "fight" in the case of planar N's. Those hopes don't help us right now with smirnoff99Frosst or whatever we are calling our current development version. Currently perhaps the external N angle needs to look on both sides to see if there is an sp2 center, as in [*:1]-;!@[#7X3;r5:2](-;@[#6X3,#7X2;r5])~;@[*;r5:3] ; what do youthink? We should also make sure it is ~125 if we have an external double bond.

[davidlmobley]

Agree we should eventually deal with nitrogens separately.

@bannanc can you address?

[bannanc]

Yes, I got a bit distracted away from this yesterday, PRs coming momentarily.

[bannanc]

I'm closing this since we fixed it in PR #85