[WIP] Feature branch: Topology Biopolymer Refactor

[j-wags]

Feature branch plans: https://openforcefield.atlassian.net/wiki/spaces/IN/pages/1579122724/2021+Topology+Refactor

Do not edit this branch directly -- All changes should go through PR review!

Notes and to-dos inherited from #881

Adds cachetools as a requirement

• Addresses Adding lur_cache to toolkits #536
• (temporary, should be replaced by a more general solution) Adds a lightweight, index-sensitive hash method for molecules (FrozenMolecule.ordered_connection_table_hash)
• Adds caching for to_openeye and to_rdkit using cachelib
• Adds caching for Atom.molecule_atom_index, since the data structure makes self-index lookup an O(n_atoms_in_molecule) operation. Molecule._invalidate_cached_properties now recursively deletes this cached property on the molecule's Atoms.
• Removes @classmethod designation from OpenEyeToolkitWrapper.to_openeye and RDKitToolkitWrapper.to_rdkit (I don't see a need for them to be classmethods, this was probably from before when I knew what I was doing)
• Makes Molecule.to_openeye and to_rdkit accept a toolkit registry kwarg. It's unclear why they didn't before (a user could reasonably want to divert Molecule.to_rdkit through a custom ToolkitWrapper)
• Incorporate suggestions from Improving SMARTS matching performance with OpenEye backend #992
• Add tests
  • appropriate cached molecule hash invalidation (eg when an atom/bond is added)
  • appropriate cached atom index invalidation (necessary? Atoms can only be appended at present)
  • cache for to_openeye / to_rdkit actually detects when a aromaticity model is requested
  • test that operations on cached outputs don't mutate the object in the cache
• Update docstrings/documentation, if applicable
• Lint codebase
• Update changelog

Notes and to-dos inherited from #973

Remaining things that this doesn't do are:

• Standard caps
  • Partial support via [WIP] Perceive common caps #1072 . We need a way to programatically extract the information from CCD or something similar.
• uncapped C terminal
• Main-chain prolines
  • (possibly related) Figure out what's up with the proline substructure that gets an N- in it I'm pretty sure this is a problem with components.cif

Notes and to-dos inherited from #982

• Add pycifrw as dependency to some recipes
• Add tests for substructure dictionary creation

Notes and to-dos from #1072

• Better way to deal with strict chirality. Check code dev notes We won't do this - A partnering company has asked us to remove stereomarks from protein FFs.

Notes and to-dos from #1085

• Implement Molecule-level HierarchySchemes and HierarchyElements, providing ways to generate a static index where atoms are grouped according to their common metadata.
• In our original API design/plan and there is an option where we want to initialize a Molecule with some hierarchy flavor.
• Add an example that both details how the default behavior is implemented and also implements custom hierarchy scheme.
• Refactor to make Topology objects consist of explicit copies of Molecules, deprecating TopologyMolecule, TopologyAtom, TopologyVirtualSite, etc
• Implement access to HierarchySchemes from the Topology level.
• Implement to/from_dict methods for Topology
  protein_param_test.tar.gz

Notes and to-dos from #1097

• Refactor to make Topology objects consist of explicit copies of Molecules, removing TopologyMolecule, TopologyAtom, TopologyVirtualSite, etc
• Refactor test_topology.py to use basic pytest features
• Implement access to HierarchySchemes from the Topology level, using offtop.hierarchy_iterator (I'm happy to change the name of this in a future PR)
• Implement to/from_dict methods for Topology
• Add Molecule.particle_index(particle), atom(index), atom_index(atom), virtual_site(index), virtual_site_index(vsite), virtual_site_particle_start_index(vsite), and bond(index).
• Add Topology.atoms, atom_index(atom), particle_index(particle), virtual_site_particle_start_index(virtual_site), molecule_index(molecule), molecule_atom_start_index(molecule), molecule_virtual_particle_start_index(molecule)
• Fixes Topologies cannot be instantiated from other topologies #946 (but SHOULDN'T close it until the refactor is merged to master)
• Implement duplicate-molecule-detection to save time in ToolkitAM1BCCHandler (done in Improve topology refactor performance, offer API point to group molecules #1140)
• Implement duplicate-molecule-detection to save time during WBO assignment
• Test speed of atom/particle indexing operations on new Topology class and implement heuristics to recover original performance
• Have all deleted API points now either redirect to a equivalent new API point or raise a descriptive error
• Decide whether a topology should have an aromaticity model, or if it's redundant with the Molecule and ForceField aromaticity models
• Ensure that all API points for permutations of {Molecule,Topology}.{atom,particle,virtual_site}{,s,_index} and counters like n_atomsare exposed and tested.
• Wire up residue and chain information to play nice with Topology.{to,from}_openmm
• Update docstrings/documentation, if applicable
• Lint codebase

Notes and to-dos from #1105

• Prototype implementation of reading a molecule from a PDB consisting of entirely "standard" components and assigning bond orders and formal charges from aa_variants.cif (and some manually-added caps)
• Fix issue where our peptide bond SMARTS could mislabel PRO ring-closing carbon as CA (maybe go over metadata assignment substructures in some sort of order that puts peptide+disulfide bonds last, and make earlier matches take precedence)
• Code cleanup
• More extensive testing
• Ensure that perceive_residues still works with newly-generated substructure files
• Resolve whether the same file can be used for BOTH metadata assignment (perceive_residues) AND bond order/formal charge assignment (from_pdb)
  • It can't -- The atom-naming dict needs to remove imidazole+guanidinium bond orders and formal charges to catch resonance states. The chemistry-assigning dict needs to know that exact information.
• Make stereo warnings quiet (maybe don't do from_rdkit at all and just make offmol directly in OFFTK API?)
  • Loading process now ingests coordinates and perceives stereo from that.
• Remove direct rdkit dependence from molecule.py
• Add extensive tests for PDB loading (take PDBs from PLBenchmarks, COVID moonshot?, other sources?)
• Add logic to gracefully error for unrecognized residues in PDB (like TPO)
• Add API points to support adding additional residues to the substructure dict
• Handle loading PDBs with multiple components (even if they're a mix of polymer and small molecules) (do we really want to support this? It's easy to split PDBs) This is out-of-scope for the 0.11.0 release
• Offer API point for residue substructure perception when loading OpenMM topologies This is out of scope for the 0.11.0 release
• Fill residue_name, residue_number metadata
• Fill atom_name metadata
• Add tests for metadata loading
• Test roundtrips to/from PDB (Guarantee that we keep atom name, atom order, atom number(maybe not this one, if original indexing doesn't start at 0?), atom name, residue name, residue number, chain name, element, coordinates)
• Make sure perceive_residues assigns appropriate residue numbers at disulfide bridges (numbering should increment at each peptide bond, NOT jump disulfide bonds)
• Schedule meeting with RCSB to talk about ways to improve aa-variants file

Notes and to-dos from #1140

• Add functionality to efficiently group identical molecules
• Use new functionality to recover SMARTS matching performance
• Polish Topology.identical_molecule_groups docstring so actual humans can read it
• Use new functionality to deduplicate molecules in special charge assignment cases (like ToolkitAM1BCC and friends)
• Add tests
• Update docstrings/documentation, if applicable

Notes and to-dos from #1195

• Refactor PDB loading into toolkit wrappers
• Refactor (to/from) (openmm,rdkit,openeye) to handle hierarchy metadata
  • Document spec for hierarchy info transfer, edge cases, and roundtrip reliability
  • Document spec for other info transfer (basically get Preserve Molecule core properties through RDMol/OEMol conversions #135 out of issue tracker and into docs)
  • Update docstrings to give basic info about hierarchy info transfer spec
  • Make public API point for loading custom substructure dicts
  • (maybe?) Don't assign info to non-tagged atoms in substructure SMARTS

Notes and to-dos from #1285

• Cherry pick vsite refactor
• (blocking) Ensure that tests from test_parameters highlighted at this link get ported to Interchange
• (blocking) Re-enable vsite_showcase notebook and get it to a working state

[codecov]

Codecov Report

Merging #951 (088fc9f) into master (97462b8) will increase coverage by 2.78%.
The diff coverage is 64.70%.

❗ Current head 088fc9f differs from pull request most recent head 8acf6ad. Consider uploading reports for the commit 8acf6ad to get more accurate results

[lgtm-com]

This pull request introduces 2 alerts when merging 519f69c into b2b76bb - view on LGTM.com

new alerts:

• 1 for Unused local variable
• 1 for Syntax error

[lgtm-com]

This pull request introduces 2 alerts when merging 7fb46c7 into b2b76bb - view on LGTM.com

new alerts:

• 1 for Unused local variable
• 1 for Unused import

[mattwthompson]

Suggested change

	from cached_property import cached_property
	from functools import cached_property

?

[j-wags]

This may be related to cached_property only being available in 3.8

[mattwthompson]

Looks like you're usiing https://pypi.org/project/cached-property/ I'm not sure if https://pypi.org/project/backports.cached-property/ would be better, and the answer comes from how each works with 3.8+; a solution in the standard library should be used over a third-party backport. I have not looked into this but it should be a quick check.

[mattwthompson]

With the T4 lysozyme, I'm getting a single atom with missing metadata, causing issues with export:

In:

import re
from collections import Counter

from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.toolkit.utils import get_data_file_path
from openff.units import unit
from simtk import unit as simtk_unit
from simtk.openmm.app import PDBFile

from openff.interchange.components.interchange import Interchange
from openff.interchange.drivers import get_gromacs_energies, get_openmm_energies
from openff.interchange.utils import get_test_file_path


def _fix_openff_impropers(force_field: ForceField):
    impropers = force_field.get_parameter_handler("ImproperTorsions")
    for imp in impropers.parameters:
        imp.smirks = re.sub("(:[23])(.+)(:[23])", r"\3\2\1", imp.smirks)

    return force_field


protein = Molecule.from_file(get_data_file_path("proteins/T4-protein.sdf"))
protein.perceive_residues()

top = protein.to_topology()

print("Counter(Length of metadata dict : number of atoms)")
print(Counter([len(a.atom.metadata) for a in top.topology_atoms]))

print("The atom with no metadata:")
for top_atom in top.topology_atoms:
    if len(top_atom.atom.metadata) == 0:
        print(top_atom)

Out:

Counter(Length of metadata dict : number of atoms)
Counter({3: 2633, 0: 1})
The atom with no metadata:
TopologyAtom 2633 with reference atom <Atom name='' atomic number='8'> and parent TopologyMolecule <openff.toolkit.topology.topology.TopologyMolecule object at 0x11993baf0>

[lgtm-com]

This pull request introduces 2 alerts when merging 58c4284 into d78912b - view on LGTM.com

new alerts:

• 1 for Unused local variable
• 1 for Unused import

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[mattwthompson]

This can probably be closed now that all of these changes are in the default branch.

[mattwthompson]

I think we're forced to "close" this but it has been implicitly merged in - the default branch is now main, which was seeded from topology-biopolymer-refactor. (Updating the base branch to main would effectively move us back a month or two in time.) Changes that would have gone into this PR are now going straight into main, just as planned, and landing in beta/RC builds of version 0.11.0.

Of course, please re-open if I'm wrong! But I don't see any other ending here.

[j-wags]

You served us well, giant PR!