Merge pull request #2 from OpenBioSim/fix_1

Fix issue #1

Author: lohedges

.github/workflows/pr.yaml

@@ -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:
@@ -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:

src/ghostly/_ghostly.py

@@ -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
@@ -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