-
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
Update topology/system particle bookkeeping in OpenMM #1051
base: develop
Are you sure you want to change the base?
Changes from all commits
5553c0b
e9b52ca
35d4d51
c8f44b0
aa9a724
694449f
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,30 @@ | ||
import pytest | ||
from openff.toolkit import Topology | ||
|
||
|
||
@pytest.mark.parametrize("collate", [True, False]) | ||
def test_collate_or_not_virtual_site_ordering( | ||
tip4p, | ||
water, | ||
collate, | ||
): | ||
pytest.importorskip("openmm") | ||
|
||
openmm_topology = tip4p.create_interchange( | ||
Topology.from_molecules([water, water]), | ||
).to_openmm_topology( | ||
collate=collate, | ||
) | ||
|
||
if collate: | ||
# particles should be O H H VS O H H VS | ||
# 0 1 2 3 4 5 6 7 | ||
virtual_site_indices = (3, 7) | ||
|
||
else: | ||
# particles should be O H H O H H VS VS | ||
# 0 1 2 3 4 5 6 7 | ||
virtual_site_indices = (6, 7) | ||
|
||
for particle_index, particle in enumerate(openmm_topology.atoms()): | ||
assert (particle.element is None) == (particle_index in virtual_site_indices) |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -15,9 +15,22 @@ | |
|
||
def to_openmm_topology( | ||
interchange: "Interchange", | ||
collate: bool = False, | ||
ensure_unique_atom_names: str | bool = "residues", | ||
) -> "openmm.app.Topology": | ||
"""Create an OpenMM Topology containing some virtual site information (if appropriate).""" | ||
""" | ||
Create an OpenMM Topology containing some virtual site information (if appropriate). | ||
|
||
Parameters | ||
---------- | ||
interchange | ||
The Interchange object to convert to an OpenMM Topology. | ||
collate | ||
If False, the default, all virtual sites will be added to a single residue at the end of the topology. | ||
If True, virtual sites will be collated with their associated molecule and added to the residue of the last | ||
atom in the molecule they belong to. | ||
|
||
""" | ||
# Heavily cribbed from the toolkit | ||
# https://github.com/openforcefield/openff-toolkit/blob/0.11.0rc2/openff/toolkit/topology/topology.py | ||
|
||
|
@@ -126,7 +139,7 @@ def to_openmm_topology( | |
last_chain = chain | ||
last_residue = residue | ||
|
||
if has_virtual_sites: | ||
if has_virtual_sites and collate: | ||
virtual_sites_in_this_molecule: list[VirtualSiteKey] = molecule_virtual_site_map[molecule_index] | ||
for this_virtual_site in virtual_sites_in_this_molecule: | ||
virtual_site_name = this_virtual_site.name | ||
|
@@ -173,6 +186,21 @@ def to_openmm_topology( | |
order=bond_order, | ||
) | ||
|
||
if has_virtual_sites and not collate: | ||
virtual_site_chain = openmm_topology.addChain(atom_chain_id) | ||
virtual_site_residue = openmm_topology.addResidue("VS", virtual_site_chain) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. From a practical standpoint, a reference to the original residue would possibly be useful - i.e. when I'm doing "solvent within N of Y" I'm looking for a residue name that matches SOL | WAT | HOH not VS. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think this would be solved by your suggestion later - that each "virtual site residue" be named according to the molecule/residue it's from? |
||
|
||
for molecule_index, molecule in enumerate(topology.molecules): | ||
virtual_sites_in_this_molecule: list[VirtualSiteKey] = molecule_virtual_site_map[molecule_index] | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Only briefly looking at things, but it looks to me like that we're really needing is this map. If we could map things to/from the collated form then yes we have to plan for this map being optionally around, but at the same time it's not crazy extra amounts of work to work around it. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I'd like to expose a map like this to the user - I have to re-create it a number of times while exporting, might as well cache it on the |
||
for this_virtual_site in virtual_sites_in_this_molecule: | ||
virtual_site_name = this_virtual_site.name | ||
|
||
openmm_topology.addAtom( | ||
virtual_site_name, | ||
virtual_site_element, | ||
virtual_site_residue, | ||
) | ||
|
||
if interchange.box is not None: | ||
from openff.units.openmm import to_openmm | ||
|
||
|
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'm sure I'm missing something important here - is there a practical advantage to doing a single residue rather than a residue per molecule's worth of virtual sites?
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.
A little easier to group everything together but nothing significant. I suppose I could have each molecule's virtual site(s) be given their own residue?
Is there an upper bound on residue numbers for some formats? A large-but-not-impractical number of waters, say 20,000, would create an extra 20,000 residues just for the virtual sites of each
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 could see AMBER and GROMACS being issues here :/ (at least iirc AMBER is still technically fixed width in parm7).