-
Notifications
You must be signed in to change notification settings - Fork 22
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.
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
Fixes residue ID increment in GROMACS GRO files #1013
Conversation
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 |
There was a problem hiding this comment.
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
A few observations after having a poke around, based on the new unit test provided:
I tried to hack together a fix based on this, but encountered an Hopefully this has been useful in pinning down the issue! I'm happy to look some more otherwise |
Just a quick thought before I dig into this later - topology files need unique molecules and can't store information about copies, but coordinate files need all copies and don't care particularly about any de-duplication. I think the |
Okay, there's a lot going on here and I need to be extra careful to only fix what needs fixing. I've opened a new issue to track requirements for various use cases. I am asking for input at the linked issue, specifically
|
@@ -445,6 +445,7 @@ def _write_gro(self, gro, decimal: int): | |||
gro.write(f"{n_particles}\n") | |||
|
|||
count = 0 | |||
residue_id = 1 |
There was a problem hiding this comment.
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
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There was a problem hiding this comment.
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.
- https://github.com/openforcefield/protein-ligand-benchmark/blob/a0b67b93c4f12f63a8bb417a1f6e4df4dc10c6c5/data/tyk2/01_protein/crd/protein.pdb#L332-L337
- https://github.com/openforcefield/protein-ligand-benchmark/blob/a0b67b93c4f12f63a8bb417a1f6e4df4dc10c6c5/data/cdk2/01_protein/crd/protein.pdb#L623-L628
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.
There was a problem hiding this comment.
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:
- Additional loop logic is needed.
- The "writer" either only checks the final unique_molecule or is overwritten by each new identical molecule.
- 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.
There was a problem hiding this comment.
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
Description
Addresses failure to increment the residue ID of identical molecule types in .gro files. Similar to #1000.
Checklist