From 526f22fcf56f84552f04f789ebbe02aec75e6e09 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 22 Oct 2026 12:26:25 +0100
Subject: [PATCH 1/3] Use geometric reference phase for stiffened rotamer
 dihedrals.

---
 CHANGELOG.md            |  1 +
 src/ghostly/_ghostly.py | 49 +++++++++++++++++++++++++++++++----------
 tests/test_ghostly.py   | 41 ++++++++++++++++++++++++++++++++++
 3 files changed, 79 insertions(+), 12 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 54639e5..8711aa8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@ 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).
 
 [2026.2.0](https://github.com/openbiosim/ghostly/compare/2026.1.0...2026.2.0) - Sep 2026
 ----------------------------------------------------------------------------------------
diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py
index f7a0e6f..9ed9251 100644
--- a/src/ghostly/_ghostly.py
+++ b/src/ghostly/_ghostly.py
@@ -119,7 +119,9 @@ def modify(
         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
@@ -2612,9 +2614,10 @@ def _check_rotamer_anchors(
     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)]
+        V(phi) = k [1 + cos(phi - phi_ref + pi)]
 
-    which has a single minimum at phi = 0 (trans) and a barrier of 2k.
+    which has a single minimum at phi_ref, the dihedral angle measured in the
+    end state where the ghost is physical, and a barrier of 2k.
 
     Parameters
     ----------
@@ -2711,10 +2714,6 @@ def _check_rotamer_anchors(
 
     # Stiffen the anchor dihedrals spanning the rotatable bonds.
 
-    from math import pi
-
-    from sire.legacy.CAS import Symbol
-
     # Get the end state property suffix.
     if is_lambda1:
         mod_key = "lambda_1"
@@ -2729,11 +2728,8 @@ def _check_rotamer_anchors(
     # 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"))
+    # The ghosts are physical at the opposite end state.
+    ref = _physical_state(mol, is_lambda1)
 
     # Initialise a container to store the updated dihedral functions.
     new_dihedrals = _SireMM.FourAtomFunctions(mol.info())
@@ -2755,6 +2751,7 @@ def _check_rotamer_anchors(
         if has_ghost_terminal and (
             bond_pair in rotatable_bonds or bond_pair_rev in rotatable_bonds
         ):
+            replacement = _single_well(ref, (idx0, idx1, idx2, idx3), k_rotamer)
             _logger.debug(
                 f"  Stiffening rotamer anchor dihedral: "
                 f"[{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], "
@@ -2779,6 +2776,34 @@ def _check_rotamer_anchors(
     return mol
 
 
+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:
+        return _morph.link_to_reference(mol)
+    else:
+        return _morph.link_to_perturbed(mol)
+
+
+def _single_well(ref, idxs, k):
+    """
+    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``.
+    """
+    from math import pi
+
+    from sire.legacy.CAS import Symbol
+
+    phi_ref = ref.dihedral(*idxs).size().to("radians")
+
+    return _SireMM.AmberDihedral(
+        _SireMM.AmberDihPart(k, 1, phi_ref - pi)
+    ).to_expression(Symbol("phi"))
+
+
 def _create_connectivity(mol):
     """
     Create a connectivity object for an end state molecule.
diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py
index 334d3be..b3ea39b 100644
--- a/tests/test_ghostly.py
+++ b/tests/test_ghostly.py
@@ -717,6 +717,47 @@ 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 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.

From 74342162ae7c79c266488a6d6de6e469fae692fa Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 22 Oct 2026 12:27:45 +0100
Subject: [PATCH 2/3] Add option to replace anchor dihedrals with a single-well
 potential.

---
 CHANGELOG.md            |   1 +
 src/ghostly/_cli.py     |  38 +++++++++++
 src/ghostly/_ghostly.py | 142 ++++++++++++++++++++++++++++++++++++++++
 tests/test_ghostly.py   |  64 ++++++++++++++++++
 4 files changed, 245 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8711aa8..93558ac 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,7 @@ 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. (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..e025f54 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,17 @@ 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)
+
     # Try to merge the reference and perturbed molecules.
     if args.system is None:
         try:
@@ -380,6 +416,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 9ed9251..38d9309 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,6 +115,22 @@ 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, after ``soften_anchors`` and before
+        ``stiffen_rotamers``, which overrides it on rotatable bonds. 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
@@ -200,6 +218,7 @@ def modify(
         "softened_angles": {},
         "softened_dihedrals": [],
         "stiffened_dihedrals": [],
+        "single_well_dihedrals": [],
         "linearised_ring_break": [],
         "softened_bonds": {},
     }
@@ -210,6 +229,7 @@ def modify(
         "softened_angles": {},
         "softened_dihedrals": [],
         "stiffened_dihedrals": [],
+        "single_well_dihedrals": [],
         "linearised_ring_break": [],
         "softened_bonds": {},
     }
@@ -440,6 +460,16 @@ def modify(
             is_lambda1=False,
         )
 
+        # Optionally replace anchor dihedrals with single-well potentials.
+        if single_well_anchors:
+            mol = _single_well_anchors(
+                mol,
+                ghosts0,
+                modifications,
+                k_single_well=k_single_well,
+                is_lambda1=False,
+            )
+
         # Check for potential rotamer anchor dihedrals.
         mol = _check_rotamer_anchors(
             mol,
@@ -554,6 +584,16 @@ def modify(
             is_lambda1=True,
         )
 
+        # Optionally replace anchor dihedrals with single-well potentials.
+        if single_well_anchors:
+            mol = _single_well_anchors(
+                mol,
+                ghosts1,
+                modifications,
+                k_single_well=k_single_well,
+                is_lambda1=True,
+            )
+
         # Check for potential rotamer anchor dihedrals.
         mol = _check_rotamer_anchors(
             mol,
@@ -2582,6 +2622,108 @@ def _soften_mixed_dihedrals(
     return mol
 
 
+def _single_well_anchors(mol, ghosts, modifications, k_single_well=5, is_lambda1=False):
+    r"""
+    Replace surviving anchor dihedrals with single-well n=1 cosine terms.
+
+    An anchor dihedral has a single ghost terminal atom and three physical
+    atoms, e.g. P2-P1-X-DR1, where X is the bridge atom. Each is replaced with:
+
+        V(phi) = k [1 + cos(phi - phi_ref + pi)]
+
+    where phi_ref is the dihedral angle measured in the end state where the
+    ghost is physical. The functional form of the anchor does not affect the
+    factorisation of the ghost partition function, so the free energy
+    contribution still cancels between legs.
+
+    Parameters
+    ----------
+
+    mol : sire.mol.Molecule
+        The perturbable molecule.
+
+    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.
+
+    k_single_well : float, optional
+        The force constant for the single-well terms. (In kcal/mol)
+
+    is_lambda1 : bool, optional
+        Whether to modify dihedrals at lambda = 1.
+
+    Returns
+    -------
+
+    mol : sire.mol.Molecule
+        The updated molecule.
+    """
+
+    # Nothing to do if there are no ghost atoms.
+    if not ghosts:
+        return mol
+
+    # Store the molecular info.
+    info = mol.info()
+
+    # Get the end state property.
+    if is_lambda1:
+        mod_key = "lambda_1"
+        suffix = "1"
+    else:
+        mod_key = "lambda_0"
+        suffix = "0"
+
+    # Get the end state dihedral functions.
+    dihedrals = mol.property("dihedral" + suffix)
+
+    # The ghosts are physical at the opposite end state.
+    ref = _physical_state(mol, is_lambda1)
+
+    # Initialise a container to store the updated dihedral functions.
+    new_dihedrals = _SireMM.FourAtomFunctions(mol.info())
+
+    # Track whether any modifications were made.
+    modified = False
+
+    # Loop over the dihedral potentials.
+    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())
+
+        atoms = (idx0, idx1, idx2, idx3)
+        is_ghost = [a in ghosts for a in atoms]
+
+        if sum(is_ghost) == 1 and (is_ghost[0] or is_ghost[3]):
+            replacement = _single_well(ref, atoms, k_single_well)
+            dih_idx_str = ",".join(str(a.value()) for a in atoms)
+            _logger.debug(
+                f"  Single-well anchor dihedral: [{dih_idx_str}], "
+                f"{p.function()} --> {replacement}"
+            )
+            new_dihedrals.set(idx0, idx1, idx2, idx3, replacement)
+            modifications[mod_key]["single_well_dihedrals"].append(dih_idx_str)
+            modified = True
+        else:
+            new_dihedrals.set(idx0, idx1, idx2, idx3, p.function())
+
+    # Set the updated dihedrals.
+    if modified:
+        mol = (
+            mol.edit()
+            .set_property("dihedral" + suffix, new_dihedrals)
+            .molecule()
+            .commit()
+        )
+
+    # Return the updated molecule.
+    return mol
+
+
 def _check_rotamer_anchors(
     mol,
     bridges,
diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py
index b3ea39b..9099e62 100644
--- a/tests/test_ghostly.py
+++ b/tests/test_ghostly.py
@@ -739,6 +739,70 @@ def test_stiffen_rotamers_reference_geometry():
         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 dihedral_energies(mol, prop, link):
     """
     Return the energy of each dihedral in the given dihedral property,

From bf572e6cd519856ee3f5af3ddec3c6556a537703 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 22 Oct 2026 12:43:32 +0100
Subject: [PATCH 3/3] Unify mixed dihedral post-processing and guard
 single-well references.

---
 CHANGELOG.md            |   2 +-
 src/ghostly/_cli.py     |   3 +
 src/ghostly/_ghostly.py | 511 +++++++++++++++-------------------------
 tests/test_ghostly.py   | 156 ++++++++++++
 4 files changed, 351 insertions(+), 321 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 93558ac..a53bdb2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,7 +7,7 @@ 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. (Default: False, 5 kcal/mol) [#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 e025f54..6ecb93d 100644
--- a/src/ghostly/_cli.py
+++ b/src/ghostly/_cli.py
@@ -391,6 +391,9 @@ def run():
     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:
diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py
index 38d9309..159e9fd 100644
--- a/src/ghostly/_ghostly.py
+++ b/src/ghostly/_ghostly.py
@@ -122,9 +122,12 @@ def modify(
         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, after ``soften_anchors`` and before
-        ``stiffen_rotamers``, which overrides it on rotatable bonds. Default
-        is False.
+        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).
@@ -451,34 +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,
-        )
-
-        # Optionally replace anchor dihedrals with single-well potentials.
-        if single_well_anchors:
-            mol = _single_well_anchors(
-                mol,
-                ghosts0,
-                modifications,
-                k_single_well=k_single_well,
-                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,
         )
 
@@ -575,34 +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,
-        )
-
-        # Optionally replace anchor dihedrals with single-well potentials.
-        if single_well_anchors:
-            mol = _single_well_anchors(
-                mol,
-                ghosts1,
-                modifications,
-                k_single_well=k_single_well,
-                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,
         )
 
@@ -2255,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()
@@ -2503,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:
 
-    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.
+    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:
+
+        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
     ----------
@@ -2527,132 +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.
 
-    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.
 
-    Returns
-    -------
-
-    mol : sire.mol.Molecule
-        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:
-        return mol
-
-    # Store the molecular info.
-    info = mol.info()
-
-    # Get the end state property.
-    if is_lambda1:
-        mod_key = "lambda_1"
-        suffix = "1"
-    else:
-        mod_key = "lambda_0"
-        suffix = "0"
-
-    # Get the end state dihedral functions.
-    dihedrals = mol.property("dihedral" + suffix)
-
-    # Initialise a container to store the updated dihedral functions.
-    new_dihedrals = _SireMM.FourAtomFunctions(mol.info())
-
-    # Track whether any modifications were made.
-    modified = False
-
-    # Loop over the dihedral potentials.
-    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())
-
-        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)
-                _logger.debug(
-                    f"  Softening mixed dihedral: [{dih_idx_str}], "
-                    f"scale={soften_anchors}"
-                )
-                modifications[mod_key]["softened_dihedrals"].append(dih_idx_str)
-            else:
-                _logger.debug(
-                    f"  Removing mixed dihedral: [{dih_idx_str}], {p.function()}"
-                )
-                modifications[mod_key]["removed_dihedrals"].append(dih_idx_str)
-            modified = True
-        else:
-            new_dihedrals.set(idx0, idx1, idx2, idx3, p.function())
-
-    # Set the updated dihedrals.
-    if modified:
-        mol = (
-            mol.edit()
-            .set_property("dihedral" + suffix, new_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 _single_well_anchors(mol, ghosts, modifications, k_single_well=5, is_lambda1=False):
-    r"""
-    Replace surviving anchor dihedrals with single-well n=1 cosine terms.
-
-    An anchor dihedral has a single ghost terminal atom and three physical
-    atoms, e.g. P2-P1-X-DR1, where X is the bridge atom. Each is replaced with:
-
-        V(phi) = k [1 + cos(phi - phi_ref + pi)]
-
-    where phi_ref is the dihedral angle measured in the end state where the
-    ghost is physical. The functional form of the anchor does not affect the
-    factorisation of the ghost partition function, so the free energy
-    contribution still cancels between legs.
-
-    Parameters
-    ----------
-
-    mol : sire.mol.Molecule
-        The perturbable molecule.
-
-    ghosts : List[sire.legacy.Mol.AtomIdx]
-        The list of ghost atoms at the current end state.
+    soften_anchors : float, optional
+        Scale factor for mixed dihedral force constants (0.0 to 1.0).
 
-    modifications : dict
-        A dictionary to store details of the modifications made.
+    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)
 
-    is_lambda1 : bool, optional
-        Whether to modify dihedrals at lambda = 1.
+    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
     -------
@@ -2665,6 +2570,24 @@ def _single_well_anchors(mol, ghosts, modifications, k_single_well=5, is_lambda1
     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.
     info = mol.info()
 
@@ -2679,9 +2602,6 @@ def _single_well_anchors(mol, ghosts, modifications, k_single_well=5, is_lambda1
     # Get the end state dihedral functions.
     dihedrals = mol.property("dihedral" + suffix)
 
-    # The ghosts are physical at the opposite end state.
-    ref = _physical_state(mol, is_lambda1)
-
     # Initialise a container to store the updated dihedral functions.
     new_dihedrals = _SireMM.FourAtomFunctions(mol.info())
 
@@ -2698,18 +2618,62 @@ def _single_well_anchors(mol, ghosts, modifications, k_single_well=5, is_lambda1
         atoms = (idx0, idx1, idx2, idx3)
         is_ghost = [a in ghosts for a in atoms]
 
-        if sum(is_ghost) == 1 and (is_ghost[0] or is_ghost[3]):
-            replacement = _single_well(ref, atoms, k_single_well)
-            dih_idx_str = ",".join(str(a.value()) 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"  Stiffening rotamer anchor dihedral: [{dih_idx_str}], "
+                    f"{p.function()} --> {function}"
+                )
+                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"  Single-well anchor dihedral: [{dih_idx_str}], "
+                    f"{p.function()} --> {function}"
+                )
+                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"  Single-well anchor dihedral: [{dih_idx_str}], "
-                f"{p.function()} --> {replacement}"
+                f"  Softening mixed dihedral: [{dih_idx_str}], scale={soften_anchors}"
             )
-            new_dihedrals.set(idx0, idx1, idx2, idx3, replacement)
-            modifications[mod_key]["single_well_dihedrals"].append(dih_idx_str)
-            modified = True
+            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:
@@ -2724,25 +2688,11 @@ def _single_well_anchors(mol, ghosts, modifications, k_single_well=5, is_lambda1
     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
           \   /
@@ -2752,20 +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 - phi_ref + pi)]
-
-    which has a single minimum at phi_ref, the dihedral angle measured in the
-    end state where the ghost is physical, 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.
@@ -2773,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.
+        Whether the bridges are at lambda = 1.
 
-    stiffen : bool, optional
-        Whether to replace rotamer anchor dihedrals with a stiff cosine
-        well. If False (default), only log warnings.
-
-    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()
@@ -2833,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,
@@ -2842,80 +2760,15 @@ 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.
-
-    # Get the end state property suffix.
-    if is_lambda1:
-        mod_key = "lambda_1"
-        suffix = "1"
-    else:
-        mod_key = "lambda_0"
-        suffix = "0"
-
-    # Store the molecular info.
-    info = mol.info()
-
-    # Get the end state dihedral functions.
-    dihedrals = mol.property("dihedral" + suffix)
-
-    # The ghosts are physical at the opposite end state.
-    ref = _physical_state(mol, is_lambda1)
-
-    # Initialise a container to store the updated dihedral functions.
-    new_dihedrals = _SireMM.FourAtomFunctions(mol.info())
-
-    modified = False
-
-    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())
-
-        # 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 has_ghost_terminal and (
-            bond_pair in rotatable_bonds or bond_pair_rev in rotatable_bonds
-        ):
-            replacement = _single_well(ref, (idx0, idx1, idx2, idx3), k_rotamer)
-            _logger.debug(
-                f"  Stiffening rotamer anchor dihedral: "
-                f"[{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], "
-                f"{p.function()} --> {replacement}"
-            )
-            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())
-
-    if modified:
-        mol = (
-            mol.edit()
-            .set_property("dihedral" + suffix, new_dihedrals)
-            .molecule()
-            .commit()
-        )
-
-    return mol
+    return rotatable_bonds
 
 
 def _physical_state(mol, is_lambda1=False):
@@ -2929,16 +2782,34 @@ def _physical_state(mol, is_lambda1=False):
         return _morph.link_to_perturbed(mol)
 
 
-def _single_well(ref, idxs, k):
+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``.
+    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
 
     from sire.legacy.CAS import Symbol
 
+    dih_idx_str = ",".join(str(a.value()) for a in idxs)
+
+    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
+
+    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."
+            )
+            return None
+
     phi_ref = ref.dihedral(*idxs).size().to("radians")
 
     return _SireMM.AmberDihedral(
diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py
index 9099e62..0344012 100644
--- a/tests/test_ghostly.py
+++ b/tests/test_ghostly.py
@@ -803,6 +803,162 @@ def test_fac4c_to_edo_single_well():
     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,