diff --git a/CHANGELOG.md b/CHANGELOG.md index 54639e5..a53bdb2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,8 @@ Changelog * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR. * Stiffen the intraghost angle to 180 degrees for dual junctions with two ghost branches instead of removing it. With both ghosts confined to the normal of the physical plane, removal left the same-side arrangement degenerate at the end state, giving large free energy errors for e.g. methoxy and hydroxyl perturbations when the input geometry had both ghosts on the same side [#17](https://github.com/OpenBioSim/ghostly/pull/17). +* Fixed ``--stiffen-rotamers`` placing the minimum of every stiffened dihedral at 0 degrees (cis), which eclipsed the ghosts with the anchor and stacked multiple ghosts on the same bridge. The minimum is now the dihedral angle measured in the end state where the ghosts are physical [#19](https://github.com/OpenBioSim/ghostly/pull/19). +* Added ``--single-well-anchors`` and ``--k-single-well`` to replace surviving anchor dihedrals with a single-well cosine potential whose minimum is the geometry of the end state where the ghost is physical. Multi-well anchors, e.g. the n=2 term of a dummy hydrogen on a carbonyl oxygen, can trap ghosts in an alternative well at intermediate lambda, giving hysteresis between independently equilibrated windows. Anchors with zero force constants or an ill-defined reference dihedral are left unchanged. (Default: False, 5 kcal/mol) [#19](https://github.com/OpenBioSim/ghostly/pull/19). [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026 ---------------------------------------------------------------------------------------- diff --git a/src/ghostly/_cli.py b/src/ghostly/_cli.py index cd54251..6ecb93d 100644 --- a/src/ghostly/_cli.py +++ b/src/ghostly/_cli.py @@ -151,6 +151,31 @@ def run(): required=False, ) + parser.add_argument( + "--single-well-anchors", + action=argparse.BooleanOptionalAction, + help=""" + Whether to replace surviving anchor dihedrals (one ghost terminal + atom, three physical atoms) with a single-well cosine potential + whose minimum is the dihedral angle of the end state where the + ghost is physical. This prevents ghost atoms being trapped in an + alternative well of a multi-well dihedral at intermediate lambda. + """, + default=False, + required=False, + ) + + parser.add_argument( + "--k-single-well", + type=str, + help=""" + The force constant for the single-well anchor dihedrals. The + resulting barrier height is 2 * k_single_well. + """, + default="5 kcal/mol", + required=False, + ) + parser.add_argument( "--stiffen-rotamers", action=argparse.BooleanOptionalAction, @@ -356,6 +381,20 @@ def run(): logger.error("k-rotamer must have units of kcal/mol") sys.exit(1) + # Try to parse the k-single-well value. + try: + k_single_well = sr.u(args.k_single_well) + except Exception as e: + logger.error(f"An error occurred while parsing the k-single-well value: {e}") + sys.exit(1) + + if not k_single_well.has_same_units(u_energy): + logger.error("k-single-well must have units of kcal/mol") + sys.exit(1) + if k_single_well.value() <= 0: + logger.error("k-single-well must be positive") + sys.exit(1) + # Try to merge the reference and perturbed molecules. if args.system is None: try: @@ -380,6 +419,8 @@ def run(): optimise_angles=args.optimise_angles, num_optimise=args.num_optimise, soften_anchors=args.soften_anchors, + single_well_anchors=args.single_well_anchors, + k_single_well=k_single_well.value(), stiffen_rotamers=args.stiffen_rotamers, k_rotamer=k_rotamer.value(), stiffen_ring_bridges=args.stiffen_ring_bridges, diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index f7a0e6f..159e9fd 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -52,6 +52,8 @@ def modify( optimise_angles=False, num_optimise=10, soften_anchors=1.0, + single_well_anchors=False, + k_single_well=5, stiffen_rotamers=False, k_rotamer=50, stiffen_ring_bridges=False, @@ -113,13 +115,34 @@ def modify( constrained too tightly, particularly when multiple ghost groups share hub atoms. + single_well_anchors : bool, optional + Whether to replace each surviving anchor dihedral, i.e. one with a + single ghost terminal atom and three physical atoms, with a single-well + n=1 cosine potential whose minimum is the dihedral angle of the end + state where the ghost is physical. Multi-well anchor dihedrals can trap + ghost atoms in an alternative well at intermediate lambda, which + appears as hysteresis between independently equilibrated windows. + Applied to all junction types. Single-well terms replace, rather than + scale, softened anchors, and ``stiffen_rotamers`` takes precedence on + rotatable bonds. Anchors with zero force constants, or whose reference + dihedral is ill-defined, are left unchanged. Has no effect when + ``soften_anchors`` is 0.0, since all mixed dihedrals are removed. + Default is False. + + k_single_well : float, optional + The force constant for the single-well anchor dihedrals (in kcal/mol). + The resulting barrier height is ``2 * k_single_well``. Only used when + ``single_well_anchors`` is True. + stiffen_rotamers : bool, optional Whether to replace rotamer anchor dihedrals with a stiff single-well cosine potential. When a bridge--physical bond is rotatable (not in a ring, sp3 bridge), surviving anchor dihedrals can allow rotameric transitions of ghost atoms at intermediate lambda. Enabling this replaces those dihedrals with a single cosine well of depth - ``2 * k_rotamer``. Default is False (only log warnings). + ``2 * k_rotamer``, with its minimum at the dihedral angle of the end + state where the ghosts are physical. Default is False (only log + warnings). k_rotamer : float, optional The force constant for the replacement cosine well when stiffening @@ -198,6 +221,7 @@ def modify( "softened_angles": {}, "softened_dihedrals": [], "stiffened_dihedrals": [], + "single_well_dihedrals": [], "linearised_ring_break": [], "softened_bonds": {}, } @@ -208,6 +232,7 @@ def modify( "softened_angles": {}, "softened_dihedrals": [], "stiffened_dihedrals": [], + "single_well_dihedrals": [], "linearised_ring_break": [], "softened_bonds": {}, } @@ -429,24 +454,19 @@ def modify( mol, ghosts0, modifications, skip_ghosts=linearised0, is_lambda1=False ) - # Soften any surviving mixed ghost/physical dihedrals. - mol = _soften_mixed_dihedrals( - mol, - ghosts0, - modifications, - soften_anchors=soften_anchors, - is_lambda1=False, - ) - - # Check for potential rotamer anchor dihedrals. - mol = _check_rotamer_anchors( + # Soften, single-well or stiffen surviving mixed ghost/physical dihedrals. + mol = _modify_mixed_dihedrals( mol, bridges0, physical0, ghosts0, + ghosts1, modifications, is_lambda1=False, - stiffen=stiffen_rotamers, + soften_anchors=soften_anchors, + single_well_anchors=single_well_anchors, + k_single_well=k_single_well, + stiffen_rotamers=stiffen_rotamers, k_rotamer=k_rotamer, ) @@ -543,24 +563,19 @@ def modify( mol, ghosts1, modifications, skip_ghosts=linearised1, is_lambda1=True ) - # Soften any surviving mixed ghost/physical dihedrals. - mol = _soften_mixed_dihedrals( - mol, - ghosts1, - modifications, - soften_anchors=soften_anchors, - is_lambda1=True, - ) - - # Check for potential rotamer anchor dihedrals. - mol = _check_rotamer_anchors( + # Soften, single-well or stiffen surviving mixed ghost/physical dihedrals. + mol = _modify_mixed_dihedrals( mol, bridges1, physical1, ghosts1, + ghosts0, modifications, is_lambda1=True, - stiffen=stiffen_rotamers, + soften_anchors=soften_anchors, + single_well_anchors=single_well_anchors, + k_single_well=k_single_well, + stiffen_rotamers=stiffen_rotamers, k_rotamer=k_rotamer, ) @@ -2213,9 +2228,7 @@ def _linearise_ring_break( from sire.convert import to_rdkit as _to_rdkit from rdkit.Chem import ChiralType as _ChiralType - phys_mol = ( - _morph.link_to_reference(mol) if is_lambda1 else _morph.link_to_perturbed(mol) - ) + phys_mol = _physical_state(mol, is_lambda1) rdmol = _to_rdkit(phys_mol) linearised = set() @@ -2461,23 +2474,48 @@ def _remove_ghost_centre_angles( return mol -def _soften_mixed_dihedrals( - mol, ghosts, modifications, soften_anchors=1.0, is_lambda1=False +def _modify_mixed_dihedrals( + mol, + bridges, + physical, + ghosts, + ref_ghosts, + modifications, + is_lambda1=False, + soften_anchors=1.0, + single_well_anchors=False, + k_single_well=5, + stiffen_rotamers=False, + k_rotamer=50, ): r""" - Soften surviving mixed ghost/physical dihedral terms by scaling their - force constants. This is a post-processing step that runs after all - per-bridge junction handlers and residual removal passes. + Post-process surviving mixed ghost/physical dihedral terms. This runs + after all per-bridge junction handlers and residual removal passes. A "mixed" dihedral is one that involves at least one ghost atom and at - least one physical atom. These dihedrals couple the ghost and physical - regions and can cause dynamics crashes at small lambda when the ghost - atoms start gaining softcore nonbonded interactions but are constrained - too tightly by bonded terms. + least one physical atom. Each is modified at most once, with the first + applicable action taking precedence: + + 1. Removed, when ``soften_anchors`` is 0.0. + 2. Stiffened, when ``stiffen_rotamers`` is True and the dihedral spans a + rotatable sp3 bridge--physical bond with a ghost terminal atom. + 3. Replaced with a single well of ``k_single_well``, when + ``single_well_anchors`` is True and the dihedral is a non-zero anchor, + i.e. a single ghost terminal atom and three physical atoms. + 4. Scaled by ``soften_anchors``, when it is less than 1.0. + + Stiffened and single-well terms are n=1 cosine wells: - When ``soften_anchors`` is 1.0 (default), no modifications are made. - When 0.0, all mixed dihedrals are removed. Intermediate values scale - the force constants. + V(phi) = k [1 + cos(phi - phi_ref + pi)] + + with a single minimum at phi_ref, the dihedral angle measured in the end + state where the ghost is physical, and a barrier of 2k. When phi_ref is + ill-defined the term falls through to the next action. + + Rotatable sp3 bridge--physical bonds are always logged as warnings when + ``stiffen_rotamers`` is False, since rotameric transitions of ghost atoms + at intermediate lambda can degrade convergence (see Boresch et al., JCTC + 2021). Parameters ---------- @@ -2485,17 +2523,41 @@ def _soften_mixed_dihedrals( mol : sire.mol.Molecule The perturbable molecule. + bridges : dict + A dictionary mapping bridge atoms to their ghost neighbours. + + physical : dict + A dictionary mapping bridge atoms to their physical neighbours. + ghosts : List[sire.legacy.Mol.AtomIdx] The list of ghost atoms at the current end state. + ref_ghosts : List[sire.legacy.Mol.AtomIdx] + The list of ghost atoms at the opposite end state. + modifications : dict A dictionary to store details of the modifications made. + is_lambda1 : bool, optional + Whether to modify dihedrals at lambda = 1. + soften_anchors : float, optional Scale factor for mixed dihedral force constants (0.0 to 1.0). - is_lambda1 : bool, optional - Whether to modify dihedrals at lambda = 1. + single_well_anchors : bool, optional + Whether to replace anchor dihedrals with single-well terms. + + k_single_well : float, optional + The force constant for the single-well terms. (In kcal/mol) + + stiffen_rotamers : bool, optional + Whether to replace rotamer anchor dihedrals with stiff single-well + terms. + + k_rotamer : float, optional + The force constant for the stiffened rotamer terms. (In kcal/mol) + The default of 50 is a placeholder and has not been calibrated + against simulation data. It should be benchmarked before use. Returns ------- @@ -2504,8 +2566,26 @@ def _soften_mixed_dihedrals( The updated molecule. """ - # Nothing to do if there are no ghost atoms or no softening is requested. - if not ghosts or soften_anchors >= 1.0: + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + + ghosts = set(ghosts) + ref_ghosts = set(ref_ghosts) + + # The ghosts are physical at the opposite end state. + ref = _physical_state(mol, is_lambda1) + + rotatable_bonds = _find_rotatable_bonds( + ref, bridges, physical, is_lambda1=is_lambda1, warn=not stiffen_rotamers + ) + + # Nothing to do if no modifications are requested. + if ( + soften_anchors >= 1.0 + and not single_well_anchors + and not (stiffen_rotamers and rotatable_bonds) + ): return mol # Store the molecular info. @@ -2536,27 +2616,64 @@ def _soften_mixed_dihedrals( idx3 = info.atom_idx(p.atom3()) atoms = (idx0, idx1, idx2, idx3) - has_ghost = any(a in ghosts for a in atoms) - has_physical = any(a not in ghosts for a in atoms) - - if has_ghost and has_physical: - dih_idx_str = ",".join(str(a.value()) for a in atoms) - if soften_anchors > 0.0: - scaled = p.function() * soften_anchors - new_dihedrals.set(idx0, idx1, idx2, idx3, scaled) + is_ghost = [a in ghosts for a in atoms] + + # Leave dihedrals that are entirely ghost or entirely physical. + if all(is_ghost) or not any(is_ghost): + new_dihedrals.set(idx0, idx1, idx2, idx3, p.function()) + continue + + dih_idx_str = ",".join(str(a.value()) for a in atoms) + has_ghost_terminal = is_ghost[0] or is_ghost[3] + is_anchor = sum(is_ghost) == 1 and has_ghost_terminal + + if soften_anchors <= 0.0: + _logger.debug(f" Removing mixed dihedral: [{dih_idx_str}], {p.function()}") + modifications[mod_key]["removed_dihedrals"].append(dih_idx_str) + modified = True + continue + + function = None + if ( + stiffen_rotamers + and has_ghost_terminal + and ((idx1, idx2) in rotatable_bonds or (idx2, idx1) in rotatable_bonds) + ): + function = _single_well(ref, atoms, k_rotamer, ref_ghosts) + if function is not None: _logger.debug( - f" Softening mixed dihedral: [{dih_idx_str}], " - f"scale={soften_anchors}" + f" Stiffening rotamer anchor dihedral: [{dih_idx_str}], " + f"{p.function()} --> {function}" ) - modifications[mod_key]["softened_dihedrals"].append(dih_idx_str) - else: + modifications[mod_key]["stiffened_dihedrals"].append(dih_idx_str) + + if ( + function is None + and single_well_anchors + and is_anchor + and not p.function().is_zero() + ): + function = _single_well(ref, atoms, k_single_well, ref_ghosts) + if function is not None: _logger.debug( - f" Removing mixed dihedral: [{dih_idx_str}], {p.function()}" + f" Single-well anchor dihedral: [{dih_idx_str}], " + f"{p.function()} --> {function}" ) - modifications[mod_key]["removed_dihedrals"].append(dih_idx_str) - modified = True + modifications[mod_key]["single_well_dihedrals"].append(dih_idx_str) + + if function is None and soften_anchors < 1.0: + function = p.function() * soften_anchors + _logger.debug( + f" Softening mixed dihedral: [{dih_idx_str}], scale={soften_anchors}" + ) + modifications[mod_key]["softened_dihedrals"].append(dih_idx_str) + + if function is None: + function = p.function() else: - new_dihedrals.set(idx0, idx1, idx2, idx3, p.function()) + modified = True + + new_dihedrals.set(idx0, idx1, idx2, idx3, function) # Set the updated dihedrals. if modified: @@ -2566,39 +2683,16 @@ def _soften_mixed_dihedrals( .molecule() .commit() ) - n_softened = len(modifications[mod_key]["softened_dihedrals"]) - lam = f"{_lam_sym}={int(is_lambda1)}" - if soften_anchors > 0.0: - _logger.debug( - f"Softened {n_softened} mixed ghost/physical dihedrals at " - f"{lam} (scale={soften_anchors})" - ) - else: - _logger.debug(f"Removed all mixed ghost/physical dihedrals at {lam}") # Return the updated molecule. return mol -def _check_rotamer_anchors( - mol, - bridges, - physical, - ghosts, - modifications, - is_lambda1=False, - stiffen=False, - k_rotamer=50, -): +def _find_rotatable_bonds(ref, bridges, physical, is_lambda1=False, warn=True): r""" - Detect and optionally stiffen rotamer anchor dihedrals. - - After ghost modifications, each bridge atom retains anchor dihedral(s) - that constrain the ghost group's rotation. If the bridge--physical bond - is not in a ring and the bridge atom is sp3-hybridised, the anchor - dihedral is likely a rotamer. This can allow the ghost group to flip - between rotameric states at intermediate lambda, degrading convergence - (see Boresch et al., JCTC 2021). + Find the rotatable bridge--physical bonds, i.e. those that are not in a + ring and have an sp3 bridge atom, so that surviving anchor dihedrals + across them can allow rotameric transitions of ghost atoms. R1 DR1 \ / @@ -2608,19 +2702,12 @@ def _check_rotamer_anchors( / \ R2 DR3 - Rotamer anchors are always logged as warnings. When ``stiffen=True``, - affected dihedrals (those spanning the rotatable bond with at least one - ghost terminal) are replaced with a single n=1 cosine well: - - V(phi) = k [1 + cos(phi - pi)] - - which has a single minimum at phi = 0 (trans) and a barrier of 2k. - Parameters ---------- - mol : sire.mol.Molecule - The perturbable molecule. + ref : sire.mol.Molecule + The perturbable molecule linked to the end state where the ghosts + are physical. bridges : dict A dictionary mapping bridge atoms to their ghost neighbours. @@ -2628,57 +2715,34 @@ def _check_rotamer_anchors( physical : dict A dictionary mapping bridge atoms to their physical neighbours. - ghosts : List[sire.legacy.Mol.AtomIdx] - The list of ghost atoms at the current end state. - - modifications : dict - A dictionary to store details of the modifications made. - is_lambda1 : bool, optional - Whether to check/modify the lambda = 1 end state. - - stiffen : bool, optional - Whether to replace rotamer anchor dihedrals with a stiff cosine - well. If False (default), only log warnings. + Whether the bridges are at lambda = 1. - k_rotamer : float, optional - The force constant for the replacement cosine term. The resulting - barrier height is 2 * k_rotamer. (In kcal/mol) Only used when - ``stiffen=True``. The default of 50 is a placeholder and has not - been calibrated against simulation data. It should be benchmarked - before use. + warn : bool, optional + Whether to log a warning for each rotatable bond. Returns ------- - mol : sire.mol.Molecule - The updated molecule (unchanged if ``stiffen=False``). + rotatable_bonds : set + The (physical, bridge) atom index pairs of the rotatable bonds. """ + rotatable_bonds = set() + # Nothing to do if there are no bridges. if not bridges: - return mol + return rotatable_bonds from rdkit.Chem import HybridizationType from sire.convert import to_rdkit - lam = int(is_lambda1) - - # Check hybridisation and ring membership at the state where the ghost is - # physically present. For bridges0, the ghost is physically present at - # lambda=1 (it is ghost/virtual at lambda=0). For bridges1, the ghost is - # physically present at lambda=0. try: - if is_lambda1: - rdmol = to_rdkit(_morph.link_to_reference(mol)) - else: - rdmol = to_rdkit(_morph.link_to_perturbed(mol)) + rdmol = to_rdkit(ref) except Exception as e: _logger.warning(f"Failed to convert molecule to RDKit for rotamer check: {e}") - return mol + return rotatable_bonds - # Identify rotatable bridge--physical bond pairs. - rotatable_bonds = set() for bridge in bridges: for p in physical[bridge]: b_idx = bridge.value() @@ -2688,8 +2752,7 @@ def _check_rotamer_anchors( if rd_bond is None or rd_bond.IsInRing(): continue - rd_bridge = rdmol.GetAtomWithIdx(b_idx) - hybridisation = rd_bridge.GetHybridization() + hybridisation = rdmol.GetAtomWithIdx(b_idx).GetHybridization() if hybridisation in ( HybridizationType.SP3, @@ -2697,86 +2760,61 @@ def _check_rotamer_anchors( HybridizationType.SP3D2, ): rotatable_bonds.add((p, bridge)) - if not stiffen: + if warn: _logger.warning( - f"Potential rotamer anchor at {_lam_sym} = {lam}: " + f"Potential rotamer anchor at {_lam_sym} = {int(is_lambda1)}: " f"bond {p_idx}-{b_idx} is a rotatable sp3 bond. " f"Surviving anchor dihedrals may allow rotameric " f"transitions of ghost atoms." ) - # If not stiffening, or no rotatable bonds found, return early. - if not stiffen or not rotatable_bonds: - return mol - - # Stiffen the anchor dihedrals spanning the rotatable bonds. - - from math import pi + return rotatable_bonds - from sire.legacy.CAS import Symbol - # Get the end state property suffix. +def _physical_state(mol, is_lambda1=False): + """ + Return the molecule linked to the end state where the ghost atoms of the + given end state are physical. + """ if is_lambda1: - mod_key = "lambda_1" - suffix = "1" + return _morph.link_to_reference(mol) else: - mod_key = "lambda_0" - suffix = "0" - - # Store the molecular info. - info = mol.info() + return _morph.link_to_perturbed(mol) - # Get the end state dihedral functions. - dihedrals = mol.property("dihedral" + suffix) - # Create the stiff single-well replacement: k[1 + cos(phi - pi)]. - # Minimum at phi = 0 (trans), barrier height = 2k. - replacement = _SireMM.AmberDihedral( - _SireMM.AmberDihPart(k_rotamer, 1, pi) - ).to_expression(Symbol("phi")) - - # Initialise a container to store the updated dihedral functions. - new_dihedrals = _SireMM.FourAtomFunctions(mol.info()) +def _single_well(ref, idxs, k, ref_ghosts, max_angle=170): + """ + Return a single-well n=1 cosine dihedral term, k[1 + cos(phi - phi_ref + pi)], + with its minimum at the dihedral angle phi_ref of the atoms ``idxs`` in the + reference molecule ``ref``. Returns None when phi_ref is ill-defined, i.e. + an atom is a ghost in the reference state or a bond angle is near linear. + """ + from math import pi - modified = False + from sire.legacy.CAS import Symbol - 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()) + dih_idx_str = ",".join(str(a.value()) for a in idxs) - # Check if the central bond (idx1-idx2) is a rotatable bridge bond - # and at least one terminal atom is ghost. - bond_pair = (idx1, idx2) - bond_pair_rev = (idx2, idx1) - has_ghost_terminal = idx0 in ghosts or idx3 in ghosts + if any(a in ref_ghosts for a in idxs): + _logger.warning( + f"Keeping dihedral [{dih_idx_str}]: an atom is a ghost in the " + "reference state, so the reference dihedral is undefined." + ) + return None - if has_ghost_terminal and ( - bond_pair in rotatable_bonds or bond_pair_rev in rotatable_bonds - ): - _logger.debug( - f" Stiffening rotamer anchor dihedral: " - f"[{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], " - f"{p.function()} --> {replacement}" + for angle in (idxs[:3], idxs[1:]): + if ref.angle(*angle).size().to("degrees") > max_angle: + _logger.warning( + f"Keeping dihedral [{dih_idx_str}]: a bond angle is near linear, " + "so the reference dihedral is undefined." ) - new_dihedrals.set(idx0, idx1, idx2, idx3, replacement) - dih_idx = (idx0.value(), idx1.value(), idx2.value(), idx3.value()) - dih_idx = ",".join([str(i) for i in dih_idx]) - modifications[mod_key]["stiffened_dihedrals"].append(dih_idx) - modified = True - else: - new_dihedrals.set(idx0, idx1, idx2, idx3, p.function()) + return None - if modified: - mol = ( - mol.edit() - .set_property("dihedral" + suffix, new_dihedrals) - .molecule() - .commit() - ) + phi_ref = ref.dihedral(*idxs).size().to("radians") - return mol + return _SireMM.AmberDihedral( + _SireMM.AmberDihPart(k, 1, phi_ref - pi) + ).to_expression(Symbol("phi")) def _create_connectivity(mol): diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index 334d3be..0344012 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -717,6 +717,267 @@ def test_pfkfb3_48_to_47_skipped_neighbour(): assert str(p.function()) == original[idx] +def test_stiffen_rotamers_reference_geometry(): + """ + Test that stiffened rotamer anchor dihedrals have their minimum at the + geometry of the end state where the ghosts are physical. The three ghosts + on bridge 3 of hexane to propane are staggered at 60, 180 and 300 degrees + about the 2-3 bond. + """ + + import sire.morph as morph + + mols = sr.load_test_files("hex2prp.s3") + + new_mols, modifications = modify(mols, stiffen_rotamers=True, k_rotamer=50) + + stiffened = modifications["lambda_1"]["stiffened_dihedrals"] + assert set(stiffened) == {"1,2,3,4", "1,2,3,13", "1,2,3,14"} + + energies = dihedral_energies(new_mols[0], "dihedral1", morph.link_to_reference) + for dihedral in stiffened: + assert abs(energies[dihedral]) < 1e-6 + + +def test_single_well_anchors(): + """ + Test that anchor dihedrals are replaced with single-well terms whose + minimum is the geometry of the end state where the ghosts are physical. + The three ghosts on bridge 3 of hexane to propane are staggered at 60, + 180 and 300 degrees about the 2-3 bond. + """ + + import sire.morph as morph + + mols = sr.load_test_files("hex2prp.s3") + + default_mols, default = modify(mols) + new_mols, modifications = modify(mols, single_well_anchors=True) + + # Only the functional form changes, not which dihedrals survive. + assert ( + new_mols[0].property("dihedral1").num_functions() + == default_mols[0].property("dihedral1").num_functions() + ) + removed = modifications["lambda_1"]["removed_dihedrals"] + assert removed == default["lambda_1"]["removed_dihedrals"] + + single_well = modifications["lambda_1"]["single_well_dihedrals"] + assert set(single_well) == {"1,2,3,4", "1,2,3,13", "1,2,3,14"} + assert modifications["lambda_0"]["single_well_dihedrals"] == [] + + energies = dihedral_energies(new_mols[0], "dihedral1", morph.link_to_reference) + for dihedral in single_well: + assert abs(energies[dihedral]) < 1e-6 + + +def test_fac4c_to_edo_single_well(): + """ + Test the single-well anchor for the ring-breaking perturbation FAC 4c to + EDO. At lambda = 1, carbon 36 becomes a carbonyl oxygen with ghost 67 + attached, a terminal junction anchored through 10-12-36-67. The default + n=2 anchor has wells at both cis and trans, so the ghost can be trapped + trans at intermediate lambda. The single well keeps only cis. + """ + + import sire.morph as morph + + mols = sr.load_test_files("fac4c_edo.s3") + + new_mols, modifications = modify(mols, single_well_anchors=True, k_single_well=5) + + assert modifications["lambda_0"]["single_well_dihedrals"] == [] + assert modifications["lambda_1"]["single_well_dihedrals"] == ["10,12,36,67"] + assert modifications["lambda_1"]["removed_dihedrals"] == ["13,12,36,67"] + + energies = dihedral_energies(new_mols[0], "dihedral1", morph.link_to_reference) + assert abs(energies["10,12,36,67"]) < 1e-6 + + # Rotating the ghost to trans gives the full barrier of 2k. + mol = morph.link_to_reference(new_mols[0]) + cursor = mol.cursor() + cursor.dihedral("atomidx 10", "atomidx 12", "atomidx 36", "atomidx 67").change_size( + 180 * sr.units.degrees, move_all=False + ) + trans = dihedral_energies(cursor.commit(), "dihedral1", lambda m: m) + assert abs(trans["10,12,36,67"] - 10) < 1e-3 + + +def test_soften_anchors(): + """ + Test that surviving mixed dihedrals are scaled by soften_anchors, and + removed when it is zero. + """ + + mols = sr.load_test_files("hex2prp.s3") + + default_mols, _ = modify(mols) + soft_mols, soft = modify(mols, soften_anchors=0.5) + removed_mols, removed = modify(mols, soften_anchors=0.0) + + anchors = {"1,2,3,4", "1,2,3,13", "1,2,3,14"} + assert anchors <= set(soft["lambda_1"]["softened_dihedrals"]) + assert anchors <= set(removed["lambda_1"]["removed_dihedrals"]) + + default = dihedral_functions(default_mols[0], "dihedral1") + softened = dihedral_functions(soft_mols[0], "dihedral1") + remaining = dihedral_functions(removed_mols[0], "dihedral1") + + for dihedral in soft["lambda_1"]["softened_dihedrals"]: + for phi in (0.0, 1.0, 2.0, 3.0): + assert ( + abs( + evaluate(softened[dihedral], phi) + - 0.5 * evaluate(default[dihedral], phi) + ) + < 1e-6 + ) + + assert not anchors & set(remaining) + + +def test_single_well_overrides_soften_anchors(): + """ + Test that single-well and stiffened rotamer terms replace softened + anchors and that each dihedral is recorded under a single modification. + """ + + mols = sr.load_test_files("hex2prp.s3") + + anchors = {"1,2,3,4", "1,2,3,13", "1,2,3,14"} + + _, modifications = modify(mols, soften_anchors=0.5, single_well_anchors=True) + assert set(modifications["lambda_1"]["single_well_dihedrals"]) == anchors + assert not anchors & set(modifications["lambda_1"]["softened_dihedrals"]) + + _, modifications = modify( + mols, soften_anchors=0.5, single_well_anchors=True, stiffen_rotamers=True + ) + assert set(modifications["lambda_1"]["stiffened_dihedrals"]) == anchors + assert modifications["lambda_1"]["single_well_dihedrals"] == [] + assert not anchors & set(modifications["lambda_1"]["softened_dihedrals"]) + + +def test_single_well_skips_zero_anchor(): + """ + Test that an anchor dihedral with a zero force constant is not turned + into a single-well restraint. + """ + + from sire.legacy.CAS import Expression + from sire.legacy.Mol import AtomIdx + + mols = sr.load_test_files("hex2prp.s3") + + dihedrals = mols[0].property("dihedral1") + dihedrals.set(AtomIdx(1), AtomIdx(2), AtomIdx(3), AtomIdx(4), Expression(0)) + mol = mols[0].edit().set_property("dihedral1", dihedrals).molecule().commit() + mols.update(mol) + + _, modifications = modify(mols, single_well_anchors=True) + + assert set(modifications["lambda_1"]["single_well_dihedrals"]) == { + "1,2,3,13", + "1,2,3,14", + } + + +def test_single_well_skips_linear_reference(): + """ + Test that an anchor dihedral is left unchanged when a bond angle of its + reference geometry is linear, so that the reference dihedral is undefined. + """ + + from sire.legacy.Mol import AtomIdx + + mols = sr.load_test_files("hex2prp.s3") + + # Place ghost 4 on the extension of the 2-3 bond in the lambda = 0 state. + c2 = mols[0].atom(AtomIdx(2)).property("coordinates0") + c3 = mols[0].atom(AtomIdx(3)).property("coordinates0") + c4 = c3 + (c3 - c2).normalise() * 1.53 + mol = mols[0].edit().atom(AtomIdx(4)).set_property("coordinates0", c4) + mols.update(mol.molecule().commit()) + + new_mols, modifications = modify(mols, single_well_anchors=True) + + assert set(modifications["lambda_1"]["single_well_dihedrals"]) == { + "1,2,3,13", + "1,2,3,14", + } + + default = dihedral_functions(modify(mols)[0][0], "dihedral1") + single_well = dihedral_functions(new_mols[0], "dihedral1") + assert str(single_well["1,2,3,4"]) == str(default["1,2,3,4"]) + + +def test_single_well_skips_reference_ghost(): + """ + Test that no single-well term is built when one of the atoms is a ghost + in the reference state. + """ + + import sire.morph as morph + from sire.legacy.Mol import AtomIdx + + from ghostly._ghostly import _single_well + + mol = sr.load_test_files("hex2prp.s3")[0] + ref = morph.link_to_reference(mol) + idxs = (AtomIdx(1), AtomIdx(2), AtomIdx(3), AtomIdx(4)) + + assert _single_well(ref, idxs, 5, set()) is not None + assert _single_well(ref, idxs, 5, {AtomIdx(1)}) is None + + +def dihedral_functions(mol, prop): + """ + Return a dictionary mapping dihedral indices to their functions. + """ + + info = mol.info() + + return { + ",".join( + str(info.atom_idx(a).value()) + for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3()) + ): p.function() + for p in mol.property(prop).potentials() + } + + +def evaluate(function, phi): + """ + Evaluate a dihedral function at the given angle in radians. + """ + + from sire.legacy.CAS import Symbol, Values + + values = Values() + values.set(Symbol("phi"), phi) + + return function.evaluate(values) + + +def dihedral_energies(mol, prop, link): + """ + Return the energy of each dihedral in the given dihedral property, + evaluated at the coordinates of the linked end state. + """ + + dihedrals = mol.property(prop) + ref = link(mol).edit().set_property("dihedral", dihedrals).molecule().commit() + info = mol.info() + + energies = {} + for p in dihedrals.potentials(): + idx = [info.atom_idx(a) for a in (p.atom0(), p.atom1(), p.atom2(), p.atom3())] + key = ",".join(str(i.value()) for i in idx) + energies[key] = ref.dihedral(*idx).energy().value() + + return energies + + def check_angle(info, potentials, idx0, idx1, idx2): """ Check if an angle potential is in a list of potentials.