Skip to content

Fix issue #1 #2

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

Merged
lohedges merged 2 commits into from
Nov 28, 2025
Merged

Fix issue #1 #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d7fdc6f
Use bridge atom hybridisation to determine planarity of triple junction
lohedges Nov 28, 2025
bfdc43c
Remove Windows PR builds. [ci skip]
lohedges Nov 28, 2025
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
5 changes: 0 additions & 5 deletions .github/workflows/pr.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@ jobs:
matrix:
python-version: ["3.10", "3.11", "3.12"]
platform:
- { name: "windows", os: "windows-latest", shell: "pwsh" }
- { name: "linux", os: "ubuntu-latest", shell: "bash -l {0}" }
- { name: "macos", os: "macos-latest", shell: "bash -l {0}" }
exclude:
Expand All @@ -23,13 +22,9 @@ jobs:
- platform:
{ name: "macos", os: "macos-latest", shell: "bash -l {0}" }
python-version: "3.10"
- platform: { name: "windows", os: "windows-latest", shell: "pwsh" }
python-version: "3.10"
- platform:
{ name: "macos", os: "macos-latest", shell: "bash -l {0}" }
python-version: "3.12" # MacOS can't run 3.12 yet...
- platform: { name: "windows", os: "windows-latest", shell: "pwsh" }
python-version: "3.11"
environment:
name: ghostly-build
defaults:
Expand Down
273 changes: 81 additions & 192 deletions src/ghostly/_ghostly.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -680,11 +680,88 @@ def _triple(
except:
connectivity = mol.property("connectivity")

# Store the element of the bridge atom.
element = mol.atom(bridge).property("element" + suffix)
# Link the molecule to the desired end state.
if is_lambda1:
end_state_mol = _morph.link_to_perturbed(mol)
else:
end_state_mol = _morph.link_to_reference(mol)

from rdkit.Chem import HybridizationType
from sire.convert import to_rdkit

# Convert the molecule to RDKit to determine the hybridisation of the bridge atom.
try:
rdmol = to_rdkit(end_state_mol)
except Exception as e:
msg = "Failed to convert molecule to RDKit for hybridisation check: " + str(e)
_logger.error(msg)
raise RuntimeError(msg)

# Get the hybridisation of the bridge atom.
hybridisation = rdmol.GetAtomWithIdx(bridge.value()).GetHybridization()

# Non-planar junction.
if (
hybridisation == HybridizationType.SP3
or hybridisation == HybridizationType.SP3D
):
_logger.debug(" Non-planar junction.")

# First, modify the force constants of the angle terms between the ghost
# atoms and the physical system to be very low.

# Get the end state angle functions.
angles = mol.property("angle" + suffix)

# Planar junction.
if element == _SireMol.Element("C"):
# Initialise a container to store the updated angle functions.
new_angles = _SireMM.ThreeAtomFunctions(mol.info())

# Indices for the softened angle terms.
angle_idxs = []

for p in angles.potentials():
idx0 = info.atom_idx(p.atom0())
idx1 = info.atom_idx(p.atom1())
idx2 = info.atom_idx(p.atom2())

if (
idx0 in ghosts
and idx2 in physical
or idx2 in ghosts
and idx0 in physical
):
from sire.legacy.CAS import Symbol

theta = Symbol("theta")

# Cast the angle to an Amber angle to get the expression.
amber_angle = _SireMM.AmberAngle(p.function(), theta)

# Create a new Amber angle with a modified force constant.

# We'll optimise the equilibrium angle for the softened angle term.
if optimise_angles:
amber_angle = _SireMM.AmberAngle(0.05, amber_angle.theta0())
angle_idxs.append((idx0, idx1, idx2))
# Use the existing equilibrium angle.
else:
amber_angle = _SireMM.AmberAngle(k_soft, amber_angle.theta0())

# Generate the new angle expression.
expression = amber_angle.to_expression(theta)

# Set the force constant to a very low value.
new_angles.set(idx0, idx1, idx2, expression)

_logger.debug(
f" Softening angle: [{idx0.value()}-{idx1.value()}-{idx2.value()}], "
f"{p.function()} --> {expression}"
)

else:
new_angles.set(idx0, idx1, idx2, p.function())

else:
_logger.debug(" Planar junction.")

# First remove all bonded terms between one of the physical atoms
Expand Down Expand Up @@ -782,194 +859,6 @@ def _triple(
.commit()
)

# Non-planar junction.
elif element == _SireMol.Element("N"):
_logger.debug(" Non-planar junction.")

# First, modify the force constants of the angle terms between the ghost
# atoms and the physical system to be very low.

# Get the end state angle functions.
angles = mol.property("angle" + suffix)

# Initialise a container to store the updated angle functions.
new_angles = _SireMM.ThreeAtomFunctions(mol.info())

# Indices for the softened angle terms.
angle_idxs = []

for p in angles.potentials():
idx0 = info.atom_idx(p.atom0())
idx1 = info.atom_idx(p.atom1())
idx2 = info.atom_idx(p.atom2())

if (
idx0 in ghosts
and idx2 in physical
or idx2 in ghosts
and idx0 in physical
):
from sire.legacy.CAS import Symbol

theta = Symbol("theta")

# Cast the angle to an Amber angle to get the expression.
amber_angle = _SireMM.AmberAngle(p.function(), theta)

# Create a new Amber angle with a modified force constant.

# We'll optimise the equilibrium angle for the softened angle term.
if optimise_angles:
amber_angle = _SireMM.AmberAngle(0.05, amber_angle.theta0())
angle_idxs.append((idx0, idx1, idx2))
# Use the existing equilibrium angle.
else:
amber_angle = _SireMM.AmberAngle(k_soft, amber_angle.theta0())

# Generate the new angle expression.
expression = amber_angle.to_expression(theta)

# Set the force constant to a very low value.
new_angles.set(idx0, idx1, idx2, expression)

_logger.debug(
f" Softening angle: [{idx0.value()}-{idx1.value()}-{idx2.value()}], "
f"{p.function()} --> {expression}"
)

else:
new_angles.set(idx0, idx1, idx2, p.function())

# Next, remove all dihedral starting from the ghost atoms and ending in
# the physical system. Also, only preserve dihedrals terminating at one
# of the physical atoms.

# Get the end state dihedral functions.
dihedrals = mol.property("dihedral" + suffix)

# Initialise containers to store the updated dihedral functions.
new_dihedrals = _SireMM.FourAtomFunctions(mol.info())

for p in dihedrals.potentials():
idx0 = info.atom_idx(p.atom0())
idx1 = info.atom_idx(p.atom1())
idx2 = info.atom_idx(p.atom2())
idx3 = info.atom_idx(p.atom3())
idxs = [idx0, idx1, idx2, idx3]

# If there is one ghost atom, then this dihedral must begin or terminate
# at the ghost atom.
num_ghosts = len([x for x in idxs if x in ghosts])
if num_ghosts == 1:
_logger.debug(
f" Removing dihedral: [{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], {p.function()}"
)
# Remove the dihedral if includes a ghost and doesn't terminate at the first
# physical atom.
elif (_is_ghost(mol, [idx0], is_lambda1)[0] and idx3 in physical[1:]) or (
_is_ghost(mol, [idx3], is_lambda1)[0] and idx0 in physical[1:]
):
_logger.debug(
f" Removing dihedral: [{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], {p.function()}"
)
else:
new_dihedrals.set(idx0, idx1, idx2, idx3, p.function())

# Update the molecule.
mol = mol.edit().set_property("angle" + suffix, new_angles).molecule().commit()
mol = (
mol.edit()
.set_property("dihedral" + suffix, new_dihedrals)
.molecule()
.commit()
)

# Optimise the equilibrium value of theta0 for the softened angle terms.
if optimise_angles:
_logger.debug(" Optimising equilibrium values for softened angles.")

import sire.morph as _morph
from sire.units import radian as _radian

# Initialise the equilibrium angle values.
theta0s = {}
for idx in angle_idxs:
theta0s[idx] = []

# Perform multiple minimisations to get an average for the theta0 values.
for _ in range(num_optimise):
# Minimise the molecule.
min_mol = _morph.link_to_reference(mol)
minimiser = min_mol.minimisation(
lambda_value=1.0 if is_lambda1 else 0.0,
constraint="none",
platform="cpu",
)
minimiser.run()

# Commit the changes.
min_mol = minimiser.commit()

# Get the equilibrium angle values.
for idx in angle_idxs:
try:
theta0s[idx].append(min_mol.angles(*idx).sizes()[0].to(_radian))
except:
raise ValueError(f"Could not find optimised angle term: {idx}")

# Compute the mean and standard error.
import numpy as _np

theta0_means = {}
theta0_stds = {}
for idx in theta0s:
theta0_means[idx] = _np.mean(theta0s[idx])
theta0_stds[idx] = _np.std(theta0s[idx]) / _np.sqrt(len(theta0s[idx]))

# Get the existing angles.
angles = mol.property("angle" + suffix)

# Initialise a container to store the updated angle functions.
new_angles = _SireMM.ThreeAtomFunctions(mol.info())

# Update the angle potentials.
for p in angles.potentials():
idx0 = info.atom_idx(p.atom0())
idx1 = info.atom_idx(p.atom1())
idx2 = info.atom_idx(p.atom2())
idx = (idx0, idx1, idx2)

# This is the softened angle term.
if idx in angle_idxs:
# Get the optimised equilibrium angle.
theta0 = theta0_means[idx]
std = theta0_stds[idx]

# Create the new angle function.
amber_angle = _SireMM.AmberAngle(k_soft, theta0)

# Generate the new angle expression.
expression = amber_angle.to_expression(Symbol("theta"))

# Set the equilibrium angle to 90 degrees.
new_angles.set(idx0, idx1, idx2, expression)

_logger.debug(
f" Optimising angle: [{idx0.value()}-{idx1.value()}-{idx2.value()}], "
f"{p.function()} --> {expression} (std err: {std:.3f} radian)"
)

else:
new_angles.set(idx0, idx1, idx2, p.function())

# Update the molecule.
mol = (
mol.edit()
.set_property("angle" + suffix, new_angles)
.molecule()
.commit()
)

# Return the updated molecule.
return mol

Expand Down