Look into 2D/3D bond graph differences for angles and torsions

[mattwthompson]

Unfortunately I think the fudging we found in #174 is happening elsewhere as well, at least for the MM conformer:

In [43]: labels = ForceField("openff-2.1.0.offxml").label_molecules(qm_molecule.to_topology())[0]

In [44]: ic = get_internal_coordinates(qm_molecule, qm_molecule.conformers[0], mm_molecule.conformers[0])

In [45]: qm_molecule.n_bonds, len(labels['Bonds']), len(ic['Bond'])
Out[45]: (84, 84, 86)

In [46]: qm_molecule.n_angles, len(labels['Angles']), len(ic['Angle'])
Out[46]: (153, 153, 132)

In [47]: qm_molecule.n_propers, len(labels['ProperTorsions']), len(ic['Dihedral'])
Out[47]: (228, 228, 239)

In [48]: qm_molecule.n_impropers, len(labels['ImproperTorsions']), len(ic['Improper'])
Out[48]: (456, 24, 27)

It's easy enough to exclude terms which are added in, but less easy to add back in what's not there (the 19 angles in this case).

The code is hard to follow; there is a method which probably infers a 2D bond graph from 3D coordinates, but that's meant to not happen in the only place I think it's called

[chapincavender]

Thanks for the sleuthing.

I think the logic for inferring bonds from coordinates is on this line. So, even when force_bonds is False, geomeTRIC will still infer bonds from coordinates if geometric_molecule.top_settings["read_bonds"] is False.

Maybe a solution is to add geometric_molecule.top_settings["read_bonds"] = True in yammbs here.

[mattwthompson]

That's an improvement:

In [8]: qm_molecule.n_bonds, len(labels['Bonds']), len(ic['Bond'])
Out[8]: (84, 84, 84)

In [9]: qm_molecule.n_angles, len(labels['Angles']), len(ic['Angle'])
Out[9]: (153, 153, 128)

In [10]: qm_molecule.n_propers, len(labels['ProperTorsions']), len(ic['Dihedral'])
Out[10]: (228, 228, 228)

In [11]: qm_molecule.n_impropers, len(labels['ImproperTorsions']), len(ic['Improper'])
Out[11]: (456, 24, 25)

[mattwthompson]

But there's still confusing behavior. With commit fe96fa1, 25 angles are skipped by geomeTRIC:

Traceback (most recent call last):
  File "/Users/mattthompson/software/yammbs/script.py", line 10, in <module>
    ic = get_internal_coordinates(qm_molecule, qm_molecule.conformers[0], mm_molecule.conformers[0])
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/mattthompson/software/yammbs/yammbs/analysis.py", line 318, in get_internal_coordinates
    assert len(openff_angles) == 0, (
           ^^^^^^^^^^^^^^^^^^^^^^^
AssertionError: Some angles were not found (0-indexed, indices sorted ascending): [(28, 29, 30), (27, 28, 29), (16, 17, 18), (8, 9, 10), (1, 3, 40), (9, 10, 11), (10, 11, 12), (12, 13, 14), (22, 23, 24), (9, 13, 14), (41, 45, 46), (16, 24, 39), (24, 25, 26), (17, 21, 22), (42, 43, 44), (25, 29, 30), (6, 41, 42), (23, 24, 25), (11, 12, 13), (25, 26, 27), (41, 42, 43), (43, 44, 45), (44, 45, 46), (15, 16, 17), (26, 27, 28)]

This is the opposite behavior of what we originally saw; instead of inferring extra angles from a funky geometry, it's skipping some angles which existing in the topology and do not seem to have bad geometry. All of the angles I've cross-references with the molecule (remembering to add 1 to indices) are C-C-C bonds in aromatic rings.

[mattwthompson]

When I look more closely at angles, I'm seeing several angles skipped over by geomeTRIC:

>>> from yammbs.analysis import get_internal_coordinates
>>> from openff.toolkit import ForceField, Molecule
>>> from yammbs.analysis import get_internal_coordinates
>>> qm_molecule = Molecule("36966574-qm.sdf"); mm_molecule = Molecule("36966574-mm.sdf")
>>> _ = get_internal_coordinates(molecule=qm_molecule, reference=qm_molecule.conformers[0], target=mm_molecule.conformers[0])['Angle']
WARNING:yammbs.analysis:Some angles were not found (0-indexed, indices sorted ascending): [(25, 26, 27), (17, 21, 22), (16, 17, 18), (12, 13, 14), (26, 27, 28), (27, 28, 29), (44, 45, 46), (8, 9, 10), (6, 41, 42), (16, 24, 39), (22, 23, 24), (23, 24, 25), (28, 29, 30), (9, 13, 14), (15, 16, 17), (24, 25, 26), (42, 43, 44), (1, 3, 40), (43, 44, 45), (10, 11, 12), (41, 45, 46), (11, 12, 13), (41, 42, 43), (9, 10, 11), (25, 29, 30)]

I was seeing this earlier (128 < 153). I have no idea how to wrestle with this.

I think this has always been baked in; here is the code this derives from: https://github.com/openforcefield/openff-sage/blob/a6a061ba8ebc144e3fda46f5f7a9e7a352b9cd2f/inputs-and-results/benchmarks/qc-opt-geo/03-a-compute-metrics.py#L46C37-L65

[mattwthompson]

#184 (at 0a3c68c) fixes the skipping of linear angles (see CO2 and acetonitrile) but something still off, probably heavy atoms in rings (6 missed in benzene, many more missed in larger molecules)

from openff.toolkit import Molecule
from yammbs.analysis import get_internal_coordinates


adenosine = Molecule.from_smiles("n2c1c(ncnc1n(c2)[C@@H]3O[C@@H]([C@@H](O)[C@H]3O)CO)N")
remdesivir = Molecule.from_smiles(
    "CCC(COC(=O)[C@@H](NP(=O)(Oc1ccccc1)OC[C@H]1O[C@@]([C@@H]([C@@H]1O)O)(C#N)c1ccc2n1ncnc2N)C)CC",
    allow_undefined_stereo=True,
)
acetonitrile = Molecule.from_smiles("CC#N")
benzene = Molecule.from_smiles("c1ccccc1")
carbon_dioxide = Molecule.from_smiles("O=C=O")

for molecule in [
    carbon_dioxide,
    acetonitrile,
    benzene,
    adenosine,
    remdesivir,
]:
    # Source - https://stackoverflow.com/a
    # Posted by Aivar Paalberg, modified by community. See post 'Timeline' for change history
    # Retrieved 2025-11-10, License - CC BY-SA 4.0
    print("\n" + molecule.to_iupac())
    molecule.generate_conformers()

    print(
        f"{molecule.n_angles=}",
        f"Num angles from geomeTRIC: {len(get_internal_coordinates(
            molecule=molecule,
            reference=molecule.conformers[0],
            target=molecule.conformers[-1],
        )['Angle'])}",
    )

carbon dioxide
molecule.n_angles=1 Num angles from geomeTRIC: 1

acetonitrile
molecule.n_angles=7 Num angles from geomeTRIC: 7

benzene
molecule.n_angles=18 Num angles from geomeTRIC: 12

(2~{R},3~{R},4~{S},5~{R})-2-(6-aminopurin-9-yl)-5-(hydroxymethyl)tetrahydrofuran-3,4-diol
molecule.n_angles=58 Num angles from geomeTRIC: 51

2-ethylbutyl (2~{S})-2-[[[(2~{R},3~{S},4~{R},5~{R})-5-(4-aminopyrrolo[2,1-f][1,2,4]triazin-7-yl)-5-cyano-3,4-dihydroxy-tetrahydrofuran-2-yl]methoxy-phenoxy-phosphoryl]amino]propanoate
Warning: : Failed due to unspecified stereochemistry
molecule.n_angles=141 Num angles from geomeTRIC: 126