Chemically equivalent angles in arginine are assigned different parameters

[epretti]

Describe the bug

Some chemically equivalent C-N-H angles in arginine are assigned different harmonic angle parameters.

To Reproduce

Save arginine.pdb.txt as arginine.pdb and run

from openff.toolkit import ForceField, Topology
from openmm import HarmonicAngleForce, unit

force_field_path = "openff-2.3.0.offxml"

topology = Topology.from_pdb("arginine.pdb")
molecule = topology.molecule(0)
system = ForceField(force_field_path).create_openmm_system(topology)
angles, = (force for force in system.getForces() if isinstance(force, HarmonicAngleForce))

for angle_index in range(angles.getNumAngles()):
    *indices, theta, k = angles.getAngleParameters(angle_index)
    if indices in ([9, 10, 20], [9, 10, 21], [9, 11, 22], [9, 11, 23]):
        print(
            *(molecule.atom(index).name for index in indices),
            theta.in_units_of(unit.degree),
            k.in_units_of(unit.kilocalorie_per_mole / unit.radian ** 2),
            sep="\t"
        )

Output

CZ      NH1     HH11    119.104777675 deg       97.22432164919 kcal/(mol rad**2)
CZ      NH1     HH12    119.104777675 deg       97.22432164919 kcal/(mol rad**2)
CZ      NH2     HH21    117.15139543910001 deg  128.6229511823 kcal/(mol rad**2)
CZ      NH2     HH22    117.15139543910001 deg  128.6229511823 kcal/(mol rad**2)

Both NH* nitrogen atoms attached to the CZ carbon atom should be chemically equivalent, as well as the HH** hydrogen atoms attached to those nitrogen atoms. However, the C-N-H angles get assigned quite different parameters.

It seems like these two different SMIRKS are matching the two pairs of angles:

openff-forcefields/openforcefields/offxml/openff-2.3.0.offxml

Lines 123 to 124 in 1ddd930

	<Angle smirks="[*:1]~[#7X3$(*~[#6X3,#6X2,#7X2+0]):2]~[*:3]" angle="119.104777675 * degree ** 1" k="97.22432164919 * kilocalorie_per_mole ** 1 * radian ** -2" id="a20"></Angle>
	<Angle smirks="[#1:1]-[#7X3$(*~[#6X3,#6X2,#7X2+0]):2]-[*:3]" angle="117.1513954391 * degree ** 1" k="128.6229511823 * kilocalorie_per_mole ** 1 * radian ** -2" id="a21"></Angle>

When loaded, the resonance form of the molecule is such that one C-N bond is a single bond and the other is a double bond:

But the different resonance forms of the guanidinium group should make all four of these angles equivalent.

Incidentally, switching out the force field for openff_no_water-3.0.0-alpha0.offxml gives

CZ      NH1     HH11    119.75608417580001 deg  156.29492916679996 kcal/(mol rad**2)
CZ      NH1     HH12    119.75608417580001 deg  156.29492916679996 kcal/(mol rad**2)
CZ      NH2     HH21    118.04593704200002 deg  70.86640775433 kcal/(mol rad**2)
CZ      NH2     HH22    118.04593704200002 deg  70.86640775433 kcal/(mol rad**2)

i.e., the issue is still present, and in fact, the difference between the force constants more significant.

Computing environment (please complete the following information):
macOS 15.7.4, conda 25.11.1, environment is conda_list.txt

Additional context
This seems related to #56, although that issue discusses torsions, and this one demonstrates a problem with angles.

[chapincavender]

Nice find! When we discovered #56, I checked the NE-CZ-NH angles but not the CZ-NH-HH angles. We've run into persistent issues with how SMILES represents delocalized charges.

Note that you get the correct behavior if you write a SMILES for arginine with the NE=CZ resonance form (first molecule below) rather than one of the CZ=NH forms (second molecule below).

from openff.toolkit import ForceField, Molecule
from openmm import HarmonicAngleForce, unit

force_field_path = "openff-2.3.0.offxml"

print("NE=CZ resonance")

molecule = Molecule.from_smiles("[NH3+][C@@H](CCC[NH+]=C(N)N)C(=O)[O-]")
topology = molecule.to_topology()
system = ForceField(force_field_path).create_openmm_system(topology)
angles, = (force for force in system.getForces() if isinstance(force, HarmonicAngleForce))

for angle_index in range(angles.getNumAngles()):
    *indices, theta, k = angles.getAngleParameters(angle_index)
    if indices in ([6, 7, 23], [6, 7, 24], [6, 8, 25], [6, 8, 26]):
        print(
            *(molecule.atom(index).name for index in indices),
            theta.in_units_of(unit.degree),
            k.in_units_of(unit.kilocalorie_per_mole / unit.radian ** 2),
            sep="\t"
        )

print("CZ=NH resonance")

molecule = Molecule.from_smiles("[NH3+][C@@H](CCCNC(=[NH2+])N)C(=O)[O-]")
topology = molecule.to_topology()
system = ForceField(force_field_path).create_openmm_system(topology)
angles, = (force for force in system.getForces() if isinstance(force, HarmonicAngleForce))

for angle_index in range(angles.getNumAngles()):
    *indices, theta, k = angles.getAngleParameters(angle_index)
    if indices in ([6, 7, 23], [6, 7, 24], [6, 8, 25], [6, 8, 26]):
        print(
            *(molecule.atom(index).name for index in indices),
            theta.in_units_of(unit.degree),
            k.in_units_of(unit.kilocalorie_per_mole / unit.radian ** 2),
            sep="\t"
        )

Produces:

NE=CZ resonance
			117.15139543910001 deg	128.6229511823 kcal/(mol rad**2)
			117.15139543910001 deg	128.6229511823 kcal/(mol rad**2)
			117.15139543910001 deg	128.6229511823 kcal/(mol rad**2)
			117.15139543910001 deg	128.6229511823 kcal/(mol rad**2)
CZ=NH resonance
			119.104777675 deg	97.22432164919 kcal/(mol rad**2)
			119.104777675 deg	97.22432164919 kcal/(mol rad**2)
			117.15139543910001 deg	128.6229511823 kcal/(mol rad**2)
			117.15139543910001 deg	128.6229511823 kcal/(mol rad**2)

So a workaround for arginine specifically could be to patch Topology.from_pdb (and probably pablo.topology_from_pdb) to make sure we always have the NE=CZ resonance form. But this is not a general solution.

Long term, we probably we want all of these angles to get the stiffer force constant, so we need to make a new angle parameter a21a with SMIRKS [#1:1]~!@[#7X3:2]~!@[#6X3:3](~!@[#7X3])~!@[#7X3]

[lilyminium]

Thanks for the report Evan! I agree with Chapin that while there are some workarounds for now, ultimately the answer would be to fix or add the SMIRKS of the correct angle type, and ideally re-train. We'll add this to our to-do list for our next force field!