Skip to content

Conversation

@lohedges
Copy link
Contributor

@lohedges lohedges commented Nov 3, 2025

This PR adds support for InverseDistanceRestraint(s) and InverseBondRestraint(s), which allows two atoms to be kept apart. This is useful for charge change perturbations where it is desirable to keep alchemical ions away from the perturbing ligand.

The PR also closes #379 by allowing the user to specify whether periodic boundary conditions should be used when creating restraints via the new Sire Python API. By default, positional and distance restraints use PBC, and bonded restraints don't. This is needed to avoid issues with bonded restraints being incorrectly satisfied across the periodic boundary, e.g. when a broken bound is being restrained.

I decided to only fix this via the new Python API for simplicity, i.e. by adding a _use_pbc attribute to the restraints objects returned by the functions in sire.restraints, then passing this through to the Sire-to-OpenMM conversion layer via the map. This avoided the need to update all of the corresponding restraint classes in Sire::MM, along with adding new versions for all streaming operators. This could still be done in future, since the Python API wouldn't need to change.

I've added a unit test that confirms that the PBC flag is set correctly for different restraints and that this propagates through to the OpenMM force. In particular, I've tested that the BondRestraintForce, which can be used for distance and bond restraints, is set an appropriate default depending on whether it is used as a distance or bond restraint, i.e. distance does use periodicity, bond doesn't. If we strictly want these to be synonyms, then I can standardise the behaviour. My interpretation was that bond refers to physically bonded atoms, whereas distance can refer to any two atoms, although I guess they could still be part of the same molecule, so issues could occur with PBC. (Maybe we could check for this if needed, then set a default based on whether they are in the same molecule?)

This will be the last PR before the 2025.3.0 release. Following this, I'll create a new BioSimSpace release, then initial conda packages for our other tools, i.e. ghostly, loch, and somd2. This will allow everyone to work with a consistent set of packages for the purposes of running initial benchmarks.

  • I confirm that I have merged the latest version of devel into this branch before issuing this pull request (e.g. by running git pull origin devel): [y]
  • I confirm that I have added a test for any new functionality in this pull request: [y]
  • I confirm that I have added documentation (e.g. a new tutorial page or detailed guide) for any new functionality in this pull request: [y]
  • I confirm that I have added a changelog entry to the changelog (we will add a link to this PR as part of the review): [y]
  • I confirm that I have permission to release this code under the GPL3 license: [y]

@akalpokas: Tagging you for reference.

@lohedges
Copy link
Contributor Author

lohedges commented Nov 4, 2025

Not sure what happened, but it appears no test input files could be read. Everything worked locally before pushing. Guess GiHub might have had a connection issue.

@lohedges
Copy link
Contributor Author

lohedges commented Nov 4, 2025

The PR also adds support for specifying the OpenCL platform index when creating an OpenMM context via the opencl_platform_index map option, which will be exposed via somd2.

@mark-mackey-cresset: Tagging you for reference. Let me know if there are any other platform related options that need exposing, other than those that are already available. (See here for a list.)

@lohedges
Copy link
Contributor Author

lohedges commented Nov 4, 2025

The ._use_pbc approach doesn't work since it's no longer possible to pickle the restraint objects. I'll need to do this the hard way, i.e. update all of the restraint classes in Sire::MM and rebuild the wrappers. Bah.

@lohedges lohedges force-pushed the feature_inverse_bond_restraint branch from 5c3d45e to 3189d66 Compare November 4, 2025 16:01
@lohedges
Copy link
Contributor Author

lohedges commented Nov 4, 2025

I've now done this the hard way. In addition, I realised that the automatic coalchemical restraints generated by SOMD2 need to be passed through in the map as a separate property, i.e. coalchemical_restraints. This is because a use might be performing a charge change perturbation for a system where a restraint is already present, so it can't be handled via the regular restraint dynamics kwarg. This new option doesn't need to be exposed to the user, since the restraint will be automatic. The user can tune the distance of restraint, though.

@lohedges
Copy link
Contributor Author

lohedges commented Nov 5, 2025

This is now working and I've tested the InverseBondRestraint via somd2, where it can be automatically applied for each alchemical ion that is inserted. Remaining questions are:

  1. Are the defaults okay? We now use use_pbc=False for bonded restraints and use_pbc=True for position and distance restraints.
  2. Do we want to add any logic to autodetect whether PBC should be used, or to raise an exception/warning if the user-specified value is likely to cause issues. For example, for a bond restraint we could check whether the atoms are part of the same molecule, e.g.:
for a0, a1 in zip(atoms0, atoms1):
    if use_pbc and a0.molecule() == a1.molecule():
        warnings.warn(
            "You are adding a bond restraint between two atoms in the same molecule. "
            "Using periodic boundary conditions is not recommended."
        )

@mark-mackey-cresset
Copy link

@mark-mackey-cresset: Tagging you for reference. Let me know if there are any other platform related options that need exposing, other than those that are already available. (See here for a list.)

No, that should be fine. We've never had an issue with OpenCL context options, it's just there were some comments in our code about somd1 behaviour and I wanted to check that they were still valid for somd2.

@lohedges
Copy link
Contributor Author

lohedges commented Nov 6, 2025

I'll merge so I can get on with the release. We can evaluate defaults over the next release cycle. Everything can be update via the Python API, or updated defaults in the class headers, so no rebuilding of wrappers etc.

@lohedges lohedges merged commit 98246ff into devel Nov 6, 2025
3 of 5 checks passed
@lohedges lohedges deleted the feature_inverse_bond_restraint branch November 6, 2025 09:42
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working cresset related to work with cresset enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

[BUG] Periodic boundary conditions for OpenMM restraints

3 participants