diff --git a/.github/workflows/pr.yaml b/.github/workflows/pr.yaml index 8cbb9d7..6b56abc 100644 --- a/.github/workflows/pr.yaml +++ b/.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: diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index c441c6c..a4b0dbe 100644 --- a/src/ghostly/_ghostly.py +++ b/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