Skip to content

Conversation

akalpokas
Copy link
Contributor

Introduction

This PR expands the restraints functionality of sire in order to include Morse potential restraints, which can be used in normal and alchemical simulations within sire. This restraints expansion in theory enables alternative route to scaffold-hopping & ring-breaking FEP already implemented in sire as part of #272.

Proposed Changes

Specifically this PR does the following:

  • Adds support for Morse restraints in the C++ layer (as part of SireMM corelib).
  • Updates sire to openMM system conversion code to create a CustomBondForce based on these restraints.
  • Updates Lambda lever code (lambdalever.cpp) to allow scaling of the restraints force with rho.
  • Builds python wrappers for all of the modifications above.
  • Exposes these restraints to the sire python API layer.
  • Expands the restraints tutorial showcasing the creation of these restraints for perturbation involving transformation of cyclopentane to cyclohexane (needs to be uploaded), and adds basic unit testing for the python code.

In the end the user can create Morse restraints, for example:

>>> mols = sr.load_test_files("cyclopentane_cyclohexane.bss")
>>> morse_restraints = sr.restraints.morse_potential(mols, atoms0=mols["molecule property is_perturbable and atomidx 0"], atoms1=mols["molecule property is_perturbable and atomidx 4"], k="100 kcal mol-1 A-2", r0="1.5 A", de="50 kcal mol-1")
>>> morse_restraint = morse_restraints[0]
>>> print(morse_restraint)
MorsePotentialRestraint( 0 <=> 4, k=100 kcal mol-1 Å-2 : r0=1.5 Å : de=50 kcal mol-1 )

creates a Morse potential restraint between atoms 0 and 4.
If the molecule contains a bond that is being alchemically annihilated at λ=1, automatic parametrisation can be used to detect the atoms and match Morse restraint parameters to the Harmonic bond as closely as possible.

>>> morse_restraints = sr.restraints.morse_potential(mols, auto_parametrise=True, de="50 kcal mol-1")
>>> morse_restraints = morse_restraints[0]
>>> print(morse_restraints)
MorsePotentialRestraint( 0 <=> 4, k=228.89 kcal mol-1 Å-2 : r0=1.5354 Å : de=50 kcal mol-1 )
 )

Because most of the restraints code is adapted from already existing restraints classes, these restraints support the same functionality as others (custom force values, usage of containers for passing atoms, multiple restraints, etc.)

Potential further additions

This PR only introduces basic Morse potential restraint functionality, and can be potentially expanded (or submitted in future as a PR) with further modifications which are currently highly experimental:

  • REST2 with angle tempering
  • Bonded term annihilation helper function (needed if merge does not introduce bond annihilation)
  • Standard ring-breaking λ-schedule
  • Direct Morse replacement functionality.

Checklists

  • 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]

Suggested reviewers:

@lohedges

Copy link
Contributor

@lohedges lohedges left a comment

Choose a reason for hiding this comment

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

Thanks, this looks great and all the tests are passing. Let me know when to take another look.

@akalpokas
Copy link
Contributor Author

akalpokas commented Oct 3, 2025

Thanks for the review, I implemented all of your feedback now. For the build tests it appears that build (3.10, linux) and build (3.12, windows) never finish for some reason.

@akalpokas akalpokas requested a review from lohedges October 3, 2025 14:19
@lohedges
Copy link
Contributor

lohedges commented Oct 3, 2025

Thanks, I'll be able to look again on Monday. Don't worry about the CI: Windows has been broken for a long time and the Linux Python 3.10 build has just become incredibly slow in the past few runs.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants