find_rotatable_bonds is biased toward input kekule form/bond orders

[j-wags]

Describe the bug

This molecule is symmetric, so it should have either 0 or 8 rotatable bonds.

The toolkit (both backends) reports that it has 4.

from openff.toolkit.topology import Molecule
mol = Molecule.from_smiles('[O-]C=CC=CC=CC=CC=O', allow_undefined_stereo=True)
print(len(mol.find_rotatable_bonds()))

4

Desired fixes

• Return either 0 or 8 rotatable bonds for this molecule
• Ensure that both OETKW and RDKTKW return the same result
• Document in the find_rotatable_bonds docstrings which behavior should be expected with conjugated bonds

Computing environment (please complete the following information):

• Mac 10.14.3
• rdkit 2021.03.1 py38hba24be9_0 conda-forge
• openff-toolkit 0.9.1

Additional context

The initial implementation, which works quite well for most cases, is here: https://github.com/openforcefield/openff-toolkit/pull/506

[j-wags]

On the plus side, both backends do seem to do the right thing for the resonance forms of amide bonds

from openff.toolkit.topology import Molecule
mol = Molecule.from_smiles('CNC(=O)C')
mol.visualize()
mol.chemical_environment_matches('[!$(*#*)&!D1:1]-&!@[!$(*#*)&!D1:2]')
mol2 = Molecule.from_smiles('C/[N+]=C(\[O-])C')
mol2.visualize()
mol.find_rotatable_bonds()