Skip angle stiffening for bridges.

Author: lohedges

src/ghostly/_cli.py

@@ -99,19 +99,6 @@ def run():
         required=False,
     )
 
-    parser.add_argument(
-        "--k-hard-ring",
-        type=str,
-        help="""
-             The force constant to use when setting angle terms involving ghost
-             atoms to 90 degrees for bridge atoms that are in a ring. A lower
-             value reduces strain where the 90 degree target fights the natural
-             ring geometry.
-             """,
-        default="75 kcal/mol/rad**2",
-        required=False,
-    )
-
     parser.add_argument(
         "--k-soft",
         type=str,
@@ -292,12 +279,6 @@ def run():
         logger.error(f"An error occurred while parsing the k-hard value: {e}")
         sys.exit(1)
 
-    try:
-        k_hard_ring = sr.u(args.k_hard_ring)
-    except Exception as e:
-        logger.error(f"An error occurred while parsing the k-hard-ring value: {e}")
-        sys.exit(1)
-
     try:
         k_soft = sr.u(args.k_soft)
     except Exception as e:
@@ -308,9 +289,6 @@ def run():
     if not k_hard.has_same_units(u):
         logger.error("k-hard must have units of kcal/mol/rad**2")
         sys.exit(1)
-    if not k_hard_ring.has_same_units(u):
-        logger.error("k-hard-ring must have units of kcal/mol/rad**2")
-        sys.exit(1)
     if not k_soft.has_same_units(u):
         logger.error("k-soft must have units of kcal/mol/rad**2")
         sys.exit(1)
@@ -352,7 +330,6 @@ def run():
         system, modifications = modify(
             system,
             k_hard.value(),
-            k_hard_ring=k_hard_ring.value(),
             k_soft=k_soft.value(),
             optimise_angles=args.optimise_angles,
             num_optimise=args.num_optimise,

src/ghostly/_ghostly.py

@@ -47,7 +47,6 @@
 def modify(
     system,
     k_hard=100,
-    k_hard_ring=75,
     k_soft=5,
     optimise_angles=True,
     num_optimise=10,
@@ -69,13 +68,6 @@ def modify(
         The force constant to use to when setting angle terms involving ghost
         atoms to 90 degrees to avoid flapping. (In kcal/mol/rad^2)
 
-    k_hard_ring : float, optional
-        The force constant to use when setting angle terms involving ghost
-        atoms to 90 degrees for bridge atoms that are in a ring. Ring bridges
-        suffer strain from the 90 degree target fighting the natural ring
-        geometry (~110 degrees). A lower force constant reduces this strain
-        while still preventing flapping. (In kcal/mol/rad^2)
-
     k_soft : float, optional
         The force constant to use when setting angle terms involving ghost atoms
         for non-planar triple junctions. (In kcal/mol/rad^2)
@@ -322,7 +314,6 @@ def modify(
                     connectivity0,
                     modifications,
                     k_hard=k_hard,
-                    k_hard_ring=k_hard_ring,
                     bridge_indices=bridge_indices0,
                 )
 
@@ -336,7 +327,6 @@ def modify(
                     connectivity0,
                     modifications,
                     k_hard=k_hard,
-                    k_hard_ring=k_hard_ring,
                     k_soft=k_soft,
                     optimise_angles=optimise_angles,
                     num_optimise=num_optimise,
@@ -417,7 +407,6 @@ def modify(
                     connectivity1,
                     modifications,
                     k_hard=k_hard,
-                    k_hard_ring=k_hard_ring,
                     is_lambda1=True,
                     bridge_indices=bridge_indices1,
                 )
@@ -431,7 +420,6 @@ def modify(
                     connectivity1,
                     modifications,
                     k_hard=k_hard,
-                    k_hard_ring=k_hard_ring,
                     k_soft=k_soft,
                     optimise_angles=optimise_angles,
                     num_optimise=num_optimise,
@@ -741,7 +729,6 @@ def _dual(
     connectivity,
     modifications,
     k_hard=100,
-    k_hard_ring=75,
     is_lambda1=False,
     bridge_indices=None,
 ):
@@ -784,10 +771,6 @@ def _dual(
         The force constant to use when setting angle terms involving ghost
         atoms to 90 degrees to avoid flapping. (In kcal/mol/rad^2)
 
-    k_hard_ring : float, optional
-        The force constant to use when the bridge atom is in a ring, where
-        the 90 degree target fights the natural ring geometry. (In kcal/mol/rad^2)
-
     is_lambda1 : bool, optional
         Whether the junction is at lambda = 1.
 
@@ -859,9 +842,6 @@ def _is_heavy(atom_idx):
         best_phys_score = 0
         heavy_phys = set(physical)
 
-    # Choose the force constant for angle stiffening.
-    k = k_hard_ring if bridge_in_ring else k_hard
-
     # Single branch.
     if len(ghosts) == 1:
         _logger.debug("  Single branch:")
@@ -926,6 +906,19 @@ def _is_heavy(atom_idx):
                 # Identify the physical atom in this angle.
                 phys_atom = idx2 if idx0 in ghosts else idx0
 
+                # Ring bridges: the ring geometry already constrains the ghost
+                # position and prevents flapping (Boresch et al. note this is
+                # "acceptable if the physical molecule is a rigid ring system").
+                # Stiffening to 90° would introduce strain without benefit.
+                if bridge_in_ring:
+                    new_angles.set(idx0, idx1, idx2, p.function())
+                    _logger.debug(
+                        f"  Skipping stiffening for angle "
+                        f"[{idx0.value()}-{idx1.value()}-{idx2.value()}]: "
+                        f"bridge atom {bridge.value()} is in a ring."
+                    )
+                    continue
+
                 # Skip stiffening through poorly-scoring atoms if better ones
                 # exist AND at least one heavy atom would remain as anchor.
                 # Without a heavy anchor, the ghost is under-constrained.
@@ -950,7 +943,7 @@ def _is_heavy(atom_idx):
                 theta0 = pi / 2.0
 
                 # Create the new angle function.
-                amber_angle = _SireMM.AmberAngle(k, theta0)
+                amber_angle = _SireMM.AmberAngle(k_hard, theta0)
 
                 # Generate the new angle expression.
                 expression = amber_angle.to_expression(Symbol("theta"))
@@ -966,7 +959,7 @@ def _is_heavy(atom_idx):
                 ang_idx = ",".join(
                     [str(i) for i in (idx0.value(), idx1.value(), idx2.value())]
                 )
-                modifications[mod_key]["stiffened_angles"][ang_idx] = {"k": k}
+                modifications[mod_key]["stiffened_angles"][ang_idx] = {"k": k_hard}
 
             else:
                 new_angles.set(idx0, idx1, idx2, p.function())
@@ -1051,7 +1044,6 @@ def _is_heavy(atom_idx):
                 connectivity,
                 modifications,
                 k_hard=k_hard,
-                k_hard_ring=k_hard_ring,
                 is_lambda1=is_lambda1,
                 bridge_indices=bridge_indices,
             )
@@ -1068,7 +1060,6 @@ def _triple(
     connectivity,
     modifications,
     k_hard=100,
-    k_hard_ring=75,
     k_soft=5,
     optimise_angles=True,
     num_optimise=10,
@@ -1113,10 +1104,6 @@ def _triple(
         The force constant to use when setting angle terms involving ghost
         atoms to 90 degrees to avoid flapping. (In kcal/mol/rad^2)
 
-    k_hard_ring : float, optional
-        The force constant to use when the bridge atom is in a ring, where
-        the 90 degree target fights the natural ring geometry. (In kcal/mol/rad^2)
-
     k_soft : float, optional
         The force constant to use when setting angle terms involving ghost
         atoms for non-planar triple junctions. (In kcal/mol/rad^2)
@@ -1535,7 +1522,6 @@ def _is_heavy(atom_idx):
             connectivity,
             modifications,
             k_hard=k_hard,
-            k_hard_ring=k_hard_ring,
             is_lambda1=is_lambda1,
             bridge_indices=bridge_indices,
         )

tests/test_ghostly.py

@@ -74,11 +74,8 @@ def test_toluene_to_pyridine():
     angles = mols[0].property("angle1")
     dihedrals = mols[0].property("dihedral1")
 
-    # Force constant for stiffening angles where the bridge atom is in a ring.
-    k_kard_ring = 75
-
     # Apply the ghost atom modifications.
-    new_mols, _ = modify(mols, k_hard_ring=k_kard_ring)
+    new_mols, _ = modify(mols)
 
     # Get the new angles and dihedrals.
     new_angles = new_mols[0].property("angle1")
@@ -119,33 +116,22 @@ def test_toluene_to_pyridine():
         for dihedral in missing_dihedrals
     )
 
-    # Create a list of angle IDs for the modified angles.
-    modified_angles = [
+    # Bridge atom 1 is in an aromatic ring. Boresch notes that for rigid ring
+    # systems the coupling is acceptable, and the ring geometry already
+    # constrains the ghost position, so no stiffening is applied. The angles
+    # involving the ring-bridge ghost should be unchanged.
+    ring_bridge_angles = [
         (AtomIdx(0), AtomIdx(1), AtomIdx(2)),
         (AtomIdx(0), AtomIdx(1), AtomIdx(6)),
     ]
 
-    # Functional form of the modified angles.
-    # Bridge atom 1 is in a ring, so k_hard_ring (75) is used.
-    expression = f"{k_kard_ring} [theta - 1.5708]^2"
-
-    # Check that the original angles don't have the modified functional form.
-    for p in angles.potentials():
-        idx0 = info.atom_idx(p.atom0())
-        idx1 = info.atom_idx(p.atom1())
-        idx2 = info.atom_idx(p.atom2())
+    for p_orig, p_new in zip(angles.potentials(), new_angles.potentials()):
+        idx0 = info.atom_idx(p_orig.atom0())
+        idx1 = info.atom_idx(p_orig.atom1())
+        idx2 = info.atom_idx(p_orig.atom2())
 
-        if (idx0, idx1, idx2) in modified_angles:
-            assert str(p.function()) != expression
-
-    # Check that the modified angles have the correct functional form.
-    for p in new_angles.potentials():
-        idx0 = info.atom_idx(p.atom0())
-        idx1 = info.atom_idx(p.atom1())
-        idx2 = info.atom_idx(p.atom2())
-
-        if (idx0, idx1, idx2) in modified_angles:
-            assert str(p.function()) == expression
+        if (idx0, idx1, idx2) in ring_bridge_angles:
+            assert str(p_orig.function()) == str(p_new.function())
 
 
 def test_acetone_to_propenol():