Single-well anchor dihedrals

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

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,