Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixes residue ID increment in GROMACS GRO files #1013

Closed
wants to merge 3 commits into from
Closed

Fixes residue ID increment in GROMACS GRO files #1013

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
357aa31
Update _export.py
PEFrankel Jul 17, 2024
ce25fc4
Update test_export.py
PEFrankel Jul 17, 2024
c39effa
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jul 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 23 additions & 0 deletions openff/interchange/_tests/unit_tests/interop/gromacs/export/test_export.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,6 +193,29 @@ def test_atom_names_pdb(self):

assert openmm_atom_names == pdb_atom_names

def test_residue_id_increment(self):
"""Test that residue IDs increment properly for multiple molecules."""

mol1 = MoleculeWithConformer.from_smiles("CCO")
mol2 = MoleculeWithConformer.from_smiles("CCO")

for atom in mol1.atoms:
atom.metadata["residue_name"] = "MOL1"
for atom in mol2.atoms:
atom.metadata["residue_name"] = "MOL2"

interchange = Interchange.from_smirnoff(
ForceField("openff-2.0.0.offxml"), [mol1, mol2]
)
interchange.to_gro("resid.gro")

with open("resid.gro") as gro_file:
lines = gro_file.readlines()
residue_ids = {int(line[:5].strip()) for line in lines[2:-2]}

# expecting residue IDs 1 and 2 for MOL1 and MOL2 respectively
assert residue_ids == {1, 2}


class TestGROMACS(_NeedsGROMACS):
@pytest.mark.slow
Expand Down
5 changes: 4 additions & 1 deletion openff/interchange/interop/gromacs/export/_export.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -445,6 +445,7 @@ def _write_gro(self, gro, decimal: int):
gro.write(f"{n_particles}\n")

count = 0
residue_id = 1
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We cant' always start residue ID at 1, this is often a different value in user inputs and ought to be respected, i.e. this protein which starts a 3 for whatever reason https://github.com/openforcefield/protein-ligand-benchmark/blob/a0b67b93c4f12f63a8bb417a1f6e4df4dc10c6c5/data/shp2/01_protein/crd/protein.pdb#L591-L621

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Believed you had moved it to #1016 or #1024 / planned on encountering it there. I'll review this again then.

It doesn't matter, but I would want to know "whatever reason" this example starts at 3 if this were the fix.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know why it's the case, but it's common enough in biophysics that we need to support it, i.e.

PLB is a relatively high-quality dataset that's valued, if not endorsed in some capacity, by consortium alumni and industry partners. It confuses me, like plenty of PDB files, but it's at least a concrete starting point to work on.

Note that this applies to residue IDs defined in biopolymers, not necessarily applicable to non-biological polymers and certainly not applicable to cases in which Interchange applies "residue IDs" on the fly just to make GROMACS happy. One of the details here that's been a thorn in my side for a decade or so is that these engines inconsistently handle residue and molecule indices, almost none handling them both at the same time.

Copy link
Author

@PEFrankel PEFrankel Aug 13, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mattwthompson
An additional loop appears to be necessary within the existing one which covers unique_molecules in order to loop through the topology_mol_idx in the lists:
example: 0: [. (0, {0: 0, 1: 1, 2: 2}), <--the existing loop only accesses the unique_molecule but we need the list int.
(1, {0: 1, 1: 0, 2: 2})]...
This can be achieved via:

        list_of_mol_tuples = unique_molecule_map[unique_molecule_index]
        for (topology_mol_idx, atom_map) in list_of_mol_tuples:
            unique_molecule_index = topology_mol_idx    # lazy diction assignment

...then in GROMACSAtom:
residue_index=unique_molecule_index + 1

This change prints the maximum molecule int in the system (originally all the minimum) which could indicate a few things:

  1. Additional loop logic is needed.
  2. The "writer" either only checks the final unique_molecule or is overwritten by each new identical molecule.
  3. GROMACSMolecule contains one GROMACSAtom set/molecule instance that overwrites itself before sending to the system.

I checked the last possibility by reviewing the contents of GROMACSMolecule after each append below GROMACSAtom (I checked this a number of ways, but this seems to be the most telling. Printing the GROMACSMolecule ID is a little confusing):

            print(f"After append: unique_molecule.name = {unique_molecule.name}")
            for a in molecule.atoms:
                print(a)

This result, interestingly, does increment the residue_index correctly (testing with the original loop will only throw one instance of course) which suggests it may be the writer or a process nearby upstream. Regardless, I think the loop here does need to be changed, possibly in addition to a writer fix.

Also, both loops in _gromacs.py may need to be changed (line 113), but the errors I got from not doing so were inconsistent.

I think I'm quite close and I could try to track down the writer in comp/interchange which my debugger accesses as the next path, but you may be able to solve this far faster than me. Let me know your thoughts. I spent a while on this and to little avail since I accidently tracked deepcopy and add_mol_keep_cache for molecules before returning to this file thinking it would help me understand the background.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be frank, there's too much going on here for me to follow, and I don't fully understand the behavior you're trying to fix. If you can open an issue with a minimally reproducing example demonstrating a difference between expectation and behavior, that would greatly help me understand the objective here.

I think #1024 there are changes you might be after? In particular I think this test looks at residue IDs in .gro files when multiple copies of a unique molecule are present

https://github.com/openforcefield/openff-interchange/pull/1024/files#diff-796291360ccaaab03d27aab793674a9c3b447a64c3f6fd2402b600611f1e36f5R200-R214

for molecule_name, molecule in self.system.molecule_types.items():
n_copies = self.system.molecules[molecule_name]

Expand All @@ -456,7 +457,8 @@ def _write_gro(self, gro, decimal: int):
f"%{decimal + 5}.{decimal}f"
f"%{decimal + 5}.{decimal}f\n"
% (
atom.residue_index, # This needs to be looked up from a different data structure
residue_id = 1,
# (hidden: 1013) atom.residue_index, # This needs to be looked up from a different data structure
Comment on lines -459 to +461
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you help me understand why the residue number should be hardcoded instead of looked up from the topology? Sometimes residue information is meaningless but it's often useful, and modifying it at export time is something I'd like to avoid

atom.residue_name[:5],
atom.name[:5],
(count + 1) % 100000,
Expand All @@ -467,6 +469,7 @@ def _write_gro(self, gro, decimal: int):
)

count += 1
residue_id += 1

if self.system.box is None:
warnings.warn(
Expand Down
Loading