Bond order of heteroatoms in aromatics

[alongd]

We should fix the treatment of heteroatoms (i.e., with lone pairs) that participate in aromatic rings.

Since RMG hard-codes a "benzene" bond as having a bond order of 1.5, we get the following error for sulfur, for example:
InvalidAdjacencyListError: Invalid valency for atom S (S4b) with 0 unpaired electrons, 1 pairs of electrons, 0 charge, and bonds [1.5,1.5].
One of the two lone pairs of sulfur participates in the aromaticity. However, with an aromatic (benzene) bond order of 1.5, RMG always gets the valency wrong:

• If we ascribe the S atom both lone pairs, its valence_electrons become 2 pairs * 2 electron/pair + 2 bonds * 1.5 bond order = 7 ≠ 6
• If we don't ascribe the lone pair that gets delocalozed to the S atom (as in the above example), its valence_electrons become 1 pairs * 2 electron/pair + 2 bonds * 1.5 bond order = 5 ≠ 6

We should consider the following cases:

• Carbon donates 3 electrons in total to the ring, hence its aromatic bond order is indeed 1.5
• Nitrogen has two cases:
  • It could either behave like carbon, what we currently call an N3b atomType, retaining its lone pair in an aromatic ring, donating 3 electrons in total to the ring, making its aromatic bond order 1.5, like carbon (e.g., pyridine)
  • But it could also share its lone pair with the ring (with a 5-member ring), what we currently call an N5b atomType, donating 4 electrons in total to the ring, making its aromatic bond order 2 (e.g., pyrrole)
• Oxygen in an aromatic ring donates one of its lone pairs, or 4 electrons in total, hence its aromatic bond order is 2.
• Sulfur behaves similarly to oxygen (aromatic bond order of 2). We currently call this an S4b atomType.

With these definitions implemented, the valency in the above example will be valid: Sulfur's valence_electrons in an aromatic ring become 1 pairs * 2 electron/pair + 2 bonds * 2 bond order = 6

I suggest we switch to calling these type of bonds aromatic, rather than benzene.
We should perhaps hard code the different aromatic bond orders according to the respective atomType.
We should also define an aromatic oxygen atomType.

[mliu49]

This is definitely a complex topic, and I think it needs to be considered from the point of view of the effect it has on downstream thermo/reaction/kinetics generation.

For benzenoid aromatics, there are defined Cb groups for thermo group additivity, such that benzene is accurately represented by the sum of 6 Cb-H groups without additional corrections. However, for heteroatom aromatics, we currently match non-aromatic groups, then apply a ring correction which captures the aromatic stabilization. If we switch to aromatic representations, then we would either change the ring correction, or maybe fit new groups for the particular heteroatom in the ring.

For reaction and kinetics generation, I'm not very familiar with reactivity of heteroatom aromatics. For benzene, all 6 bonds behave the same, in a way that's different from a normal single or double bond. Using a benzene bond representation is a straightforward way to characterize them differently. However, for heteroatom aromatics, the bonds don't seem like the would behave the same, although they would probably still be different from a typical single or double bond. Using a uniform aromatic bond type might not give the results that we want.

Another possibility is separating bond order and aromaticity into separate attributes, like RDKit does (and possibly OpenBabel too). Then you could have an aromatic single bond or an aromatic double bond. However, it might be complicated to integrate into the existing bond framework, and it's not clear to me exactly how we could take advantage of it.

[github-actions]

This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days.