diff --git a/molSimplify/Scripts/enhanced_structgen.py b/molSimplify/Scripts/enhanced_structgen.py
index 16c06e43..1a86d0f1 100644
--- a/molSimplify/Scripts/enhanced_structgen.py
+++ b/molSimplify/Scripts/enhanced_structgen.py
@@ -189,6 +189,7 @@ def generate_complex(
 
     # sterics report
     run_sterics: bool = True,
+    ann_bool: bool = True,
 ):
     """
     Build an complex from a list of ligands.
@@ -353,8 +354,7 @@ def generate_complex(
             # END NEW
 
             # clean bonds & optimize
-            core3D.convert2OBMol(force_clean=True)
-            replace_bonds(core3D.OBMol, core3D.bo_dict)
+            core3D = sync_obmol_from_bodict(core3D)
 
             optimized_coords = constrained_forcefield_optimization(
                 core3D,
@@ -376,8 +376,8 @@ def generate_complex(
                 if len(keep_piercings) != 0:
                     new_coords2, moved_atoms = correct_ring_piercings(core3D, keep_piercings)
                     core3D = set_new_coords(core3D, new_coords2)
-                    core3D.convert2OBMol(force_clean=True)
-                    replace_bonds(core3D.OBMol, core3D.bo_dict)
+                    core3D = sync_obmol_from_bodict(core3D)
+
 
                     optimized_coords = constrained_forcefield_optimization(
                         core3D,
@@ -386,8 +386,8 @@ def generate_complex(
                         ff_name=ff_name
                     )
                     core3D = set_new_coords(core3D, optimized_coords)
-                    core3D.convert2OBMol(force_clean=True)
-                    replace_bonds(core3D.OBMol, core3D.bo_dict)
+                    core3D = sync_obmol_from_bodict(core3D)
+
 
                     optimized_coords = constrained_forcefield_optimization(
                         core3D,
@@ -406,8 +406,8 @@ def generate_complex(
                         ff_name='GAFF'
                     )
                     core3D = set_new_coords(core3D, optimized_coords)
-                    core3D.convert2OBMol(force_clean=True)
-                    replace_bonds(core3D.OBMol, core3D.bo_dict)
+                    core3D = sync_obmol_from_bodict(core3D)
+
 
                     optimized_coords = constrained_forcefield_optimization(
                         core3D,
@@ -440,13 +440,6 @@ def generate_complex(
         )
         reapply_all_haptics_and_sync(core3D, bond_order=1, prefer_nearest_metal=True)
 
-    clashes = None
-    severity = None
-    # get sterics report
-    if run_sterics:
-        clashes, severity, fig = run_sterics_check(core3D, max_steps, ff_name)
-    else:
-        fig = None  # just to be explicit
 
     # NEW: produce the final ordered list of core indices per filled backbone site
     backbone_core_indices = [core for (site, core) in sorted(backbone_core_pairs, key=lambda x: x[0])]  # NEW
@@ -454,8 +447,8 @@ def generate_complex(
     # -------------------- FINAL: unconstrained FF relax --------------------
     try:
         # make sure OBMol matches our current coords/bonds before the FF pass
-        core3D.convert2OBMol(force_clean=True)
-        replace_bonds(core3D.OBMol, core3D.bo_dict)
+        core3D = sync_obmol_from_bodict(core3D)
+
 
         optimized_coords = constrained_forcefield_optimization(
             core3D,
@@ -465,23 +458,26 @@ def generate_complex(
         )
         core3D = set_new_coords(core3D, optimized_coords)
         # re-sync OBMol again for downstream writeout / sterics
-        core3D.convert2OBMol(force_clean=True)
-        replace_bonds(core3D.OBMol, core3D.bo_dict)
+        core3D = sync_obmol_from_bodict(core3D)
+
     except Exception as e:
         if verbose:
             print(f"[warn] Final unconstrained FF relax failed; keeping previous geometry. ({e})")
-    return core3D, clashes, severity, fig, batslist, backbone_core_indices
 
-def run_sterics_check(core3D, max_steps, ff_name):
-    optimized_coords, per_atom_ff_force = constrained_forcefield_optimization(
-        core3D,
-        max_steps=max_steps,
-        ff_name=ff_name,
-        return_per_atom_ff_force=True,
-        fd_delta=1e-3,
-        isolate_vdw=True,  # optional: emphasize sterics
-    )
 
+    core3D, per_atom_ff_force, optimized_coords = enforce_metal_ligand_distances_and_optimize(core3D, None, backbone_core_indices)
+    clashes = None
+    severity = None
+    # get sterics report
+    if run_sterics:
+        clashes, severity, fig = run_sterics_check(core3D, per_atom_ff_force, optimized_coords)
+    else:
+        fig = None  # just to be explicit
+
+
+    return core3D, clashes, severity, fig, batslist, backbone_core_indices
+
+def run_sterics_check(core3D, per_atom_ff_force, optimized_coords):
     elements = [at.sym for at in core3D.atoms]
     tree = KDTree(optimized_coords)
 
@@ -500,6 +496,9 @@ def run_sterics_check(core3D, max_steps, ff_name):
         clearance_HX=0.30,
         clearance_HH=0.35,
         exclude_hops=(1,2,3),
+        epsilon=0.10,
+        gamma=2.0,
+        pair_score_threshold=0.05,
     )
 
     # NOTE: pass clashes (list of pairs) to steric_pairs
@@ -899,8 +898,8 @@ def donor_metal_for(idx):
 
     # 6) Constrained FF optimization (freeze metals + donors)
     if constrain:
-        core3D.convert2OBMol(force_clean=True)
-        replace_bonds(core3D.OBMol, core3D.bo_dict)
+        core3D = sync_obmol_from_bodict(core3D)
+
         frozen = sorted(set(metal_indices + donor_idxs))
         optimized = constrained_forcefield_optimization(core3D, frozen, max_steps=max_steps, ff_name=ff_name)
         core3D = set_new_coords(core3D, optimized)
@@ -910,24 +909,26 @@ def donor_metal_for(idx):
         had_haptics = bool(getattr(core3D, "_haptic_groups_global", None) or [])
 
         # Always resync OBMol/bonds before FF
-        core3D.convert2OBMol(force_clean=True)
-        replace_bonds(core3D.OBMol, core3D.bo_dict)
+        core3D = sync_obmol_from_bodict(core3D)
+
 
         # Ask for fully-unconstrained relax.
         # If haptics exist, constrained_forcefield_optimization will automatically
         # freeze metal atoms only and disable rotor search (per our new implementation).
-        optimized2 = constrained_forcefield_optimization(
+        optimized2, per_atom_ff_force = constrained_forcefield_optimization(
             core3D,
-            fixed_atom_indices=None,
-            max_steps=final_relax_steps,
+            max_steps=max_steps,
             ff_name=ff_name,
+            return_per_atom_ff_force=True,
+            fd_delta=1e-3,
+            isolate_vdw=True,  # optional: emphasize sterics
         )
         core3D = set_new_coords(core3D, optimized2)
 
         # Final sync and reapply haptics for visual consistency
-        core3D.convert2OBMol(force_clean=True)
-        replace_bonds(core3D.OBMol, core3D.bo_dict)
+        core3D = sync_obmol_from_bodict(core3D)
+
         if had_haptics:
             reapply_all_haptics_and_sync(core3D, bond_order=1, prefer_nearest_metal=True)
 
-    return core3D
+    return core3D, per_atom_ff_force, optimized2
diff --git a/molSimplify/Scripts/enhanced_structgen_functionality.py b/molSimplify/Scripts/enhanced_structgen_functionality.py
index a0f81258..6ed713ee 100644
--- a/molSimplify/Scripts/enhanced_structgen_functionality.py
+++ b/molSimplify/Scripts/enhanced_structgen_functionality.py
@@ -32,6 +32,28 @@
 import matplotlib.colors as mcolors  # needed by visualize_molecule
 import os
 
+def apply_aromatic_flags_from_bodict(core3D):
+    obmol = core3D.OBMol
+    for (i, j), bo in core3D.bo_dict.items():
+        is_ar = (bo == "ar") or (isinstance(bo, float) and abs(bo - 1.5) < 1e-6)
+        if not is_ar:
+            continue
+        bond = obmol.GetBond(i+1, j+1)  # OB is 1-based
+        if bond is None:
+            continue
+        if hasattr(bond, "SetAromatic"):
+            bond.SetAromatic(True)
+        a = obmol.GetAtom(i+1); b = obmol.GetAtom(j+1)
+        if hasattr(a, "SetAromatic"): a.SetAromatic(True)
+        if hasattr(b, "SetAromatic"): b.SetAromatic(True)
+
+def sync_obmol_from_bodict(core3D):
+    core3D.convert2OBMol(force_clean=True)
+    replace_bonds(core3D.OBMol, core3D.bo_dict)
+    apply_aromatic_flags_from_bodict(core3D)
+    return core3D
+
+
 
 def add_ligand_to_complex(
     ligand_donor_coords,
@@ -2784,15 +2806,39 @@ def _shares_metal(u, v, donors_by_m):
 
 
 def _pca_frame(P):
-    """Return orthonormal principal axes for points P (N,3) and the centroid."""
+    """Return 3x3 orthonormal principal axes for points P (N,3) and the centroid."""
+    P = np.asarray(P, float)
     C = P.mean(axis=0)
     X = P - C
+
+    # Not enough information for a stable frame
     if X.shape[0] < 2 or np.linalg.matrix_rank(X) < 2:
         return np.eye(3), C
-    U, _, _ = np.linalg.svd(X, full_matrices=False)
-    if np.linalg.det(U) < 0:
-        U[:, -1] *= -1
-    return U, C
+
+    # SVD: X = U S Vt ; principal directions in coordinate space are columns of V = Vt.T
+    _, _, Vt = np.linalg.svd(X, full_matrices=False)
+    V = Vt.T  # shape (3,k) where k=min(N,3)
+
+    if V.shape[1] == 3:
+        R = V
+    elif V.shape[1] == 2:
+        # Complete to 3D with a right-handed third axis
+        a = V[:, 0]
+        b = V[:, 1]
+        c = np.cross(a, b)
+        cn = np.linalg.norm(c)
+        if cn < 1e-12:
+            return np.eye(3), C
+        c = c / cn
+        R = np.column_stack([a, b, c])
+    else:
+        return np.eye(3), C
+
+    # Ensure right-handed frame
+    if np.linalg.det(R) < 0:
+        R[:, -1] *= -1
+
+    return R, C
 
 
 def _pca_align(P, Q):
@@ -2993,32 +3039,25 @@ def check_sterics_with_ff_embedding(
     per_atom_nonbonded=None,
     per_atom_vdw=None,
     per_atom_ff_force=None,
-    energy_source='none',     # 'none' | 'nonbonded' | 'vdw' | 'ff_force'
-    energy_weighted=False,
-    # ---- NEW: pair-energy gate (option 4) ----
-    pair_energy_threshold: float = 0.5,  # require mean(e_i,e_j) >= this to consider a clash
-    #
+    energy_source='none',     # keep for backwards compat
+    energy_weighted=False,    # keep for backwards compat
+
+    # --- NEW (replace pair_energy_threshold gating) ---
+    epsilon: float = 0.10,
+    gamma: float = 2.0,
+    pair_score_threshold: float = 0.05,   # in Å-equivalent if E_geom in Å
+
+    # old args
+    pair_energy_threshold: float = 0.5,   # deprecated; ignore or keep for legacy
     scale=1.0,
     default_vdw=1.5,
     exclude_hops=(1, 2, 3),
-    # ---- BASE clearance for heavy-heavy ----
     clearance_heavy: float = 0.20,
-    # ---- NEW: H-aware per-pair clearances (option 3) ----
-    clearance_HH: float = 0.35,   # H···H pairs must penetrate > this to count
-    clearance_HX: float = 0.30,   # H···X (X != H) must penetrate > this to count
-    #
+    clearance_HH: float = 0.35,
+    clearance_HX: float = 0.30,
     infer_H_bonds: bool = True,
-    H_bond_max: float = 1.25      # Å: attach isolated H to nearest atom if bo_dict lacks H bonds
+    H_bond_max: float = 1.25
 ):
-    """
-    Steric clash detection with:
-      - robust bond-graph pruning (≤3 hops),
-      - H-aware per-pair clearance (HH/HX/heavy-heavy),
-      - optional pair-energy gating based on a per-atom FF embedding,
-      - optional gentle energy weighting of severity.
-
-    Returns: (clashes (list[(i,j)]), severity_scores dict[(i,j)]->float)
-    """
     import numpy as np
     import networkx as nx
 
@@ -3034,14 +3073,10 @@ def check_sterics_with_ff_embedding(
             G.add_edge(i, j)
             bond_pairs.add((min(i, j), max(i, j)))
 
-    # --- infer missing H–X bonds if Hs are isolated (common in heavy-only bo_dicts) ---
+    # --- infer missing H–X bonds (unchanged) ---
     if infer_H_bonds:
-        import numpy as np
         xyz = np.asarray(coords, float)
-        is_H: NDArray[np.bool_] = np.array(
-            [str(z).upper() == "H" for z in elements],
-            dtype=bool,
-        )
+        is_H = np.array([str(z).upper() == "H" for z in elements], dtype=bool)
         for h in range(N):
             if not is_H[h] or G.degree[h] > 0:
                 continue
@@ -3052,7 +3087,7 @@ def check_sterics_with_ff_embedding(
                 G.add_edge(h, j)
                 bond_pairs.add((min(h, j), max(h, j)))
 
-    # --- choose FF embedding for energy/force info ---
+    # --- choose per-atom embedding (we care mainly about ff_force here) ---
     energies = None
     if energy_source == 'ff_force':
         energies = per_atom_ff_force
@@ -3060,13 +3095,18 @@ def check_sterics_with_ff_embedding(
         energies = per_atom_nonbonded
     elif energy_source == 'vdw':
         energies = per_atom_vdw
+
     if energies is not None:
-        energies = np.clip(np.asarray(energies, float), 0.0, None)
+        w_raw = np.clip(np.asarray(energies, float), 0.0, None)
+        # robust normalization -> [0,1]
+        p10 = np.percentile(w_raw, 10)
+        p90 = np.percentile(w_raw, 90)
+        denom = (p90 - p10) + 1e-12
+        w = np.clip((w_raw - p10) / denom, 0.0, 1.0)
+    else:
+        w = None
 
-    is_H = np.array(
-        [str(z).upper() == "H" for z in elements],
-        dtype=bool,
-    )
+    is_H = np.array([str(z).upper() == "H" for z in elements], dtype=bool)
 
     clashes = set()
     severity_scores = {}
@@ -3075,24 +3115,17 @@ def check_sterics_with_ff_embedding(
 
     for i in range(N):
         ri = vdw_radii.get(elements[i], default_vdw)
-        e_i = (energies[i] if (energies is not None and i < len(energies)) else 0.0)
-
-        # Precompute hop distances up to cutoff (safe; no NodeNotFound)
         hop_len = nx.single_source_shortest_path_length(G, i, cutoff=max_hops) if max_hops > 0 else {}
 
         for j in tree.query_ball_point(coords[i], (ri + max_vdw) * scale):
             if i >= j:
                 continue
-
-            # bonded? skip
             if (min(i, j), max(i, j)) in bond_pairs:
                 continue
-
-            # prune 1–2/1–3/1–4 neighbors (if requested)
             if max_hops > 0 and (j in hop_len) and (hop_len[j] in exclude_hops):
                 continue
 
-            # dynamic, H-aware clearance
+            # dynamic clearance
             if is_H[i] and is_H[j]:
                 clearance = clearance_HH
             elif is_H[i] or is_H[j]:
@@ -3104,26 +3137,29 @@ def check_sterics_with_ff_embedding(
             cutoff = (ri + rj) * scale
             d = float(np.linalg.norm(coords[i] - coords[j]))
             penetration = cutoff - d
+
+            # must exceed clearance at all
             if penetration <= clearance:
                 continue
 
-            # pair-energy gate (mean of per-atom energies/forces)
-            if energies is not None and pair_energy_threshold is not None and pair_energy_threshold > 0.0:
-                e_j = energies[j] if j < len(energies) else 0.0
-                mean_e = 0.5 * (e_i + e_j)
-                if mean_e < pair_energy_threshold:
-                    continue  # too “quiet” energetically → not a real clash
+            # --- NEW: E_geom ---
+            E_geom = max(0.0, penetration - clearance)
 
-            # severity beyond the per-pair clearance
-            sev = penetration - clearance
-            if energy_weighted and (energies is not None):
-                # gentle up-weighting; avoids runaway from a single spiky atom
-                e_j = energies[j] if j < len(energies) else 0.0
-                sev *= np.log1p(0.5 * (e_i + e_j))
+            # --- NEW: pair weighting from per-atom ff_force (or other embedding) ---
+            if w is not None and i < len(w) and j < len(w):
+                w_ij = 0.5 * (w[i] + w[j])
+            else:
+                w_ij = 0.0
+
+            E_pair = E_geom * ((epsilon + w_ij) ** gamma)
+
+            # --- NEW: gate on pair score ---
+            if E_pair < pair_score_threshold:
+                continue
 
             key = (i, j)
-            severity_scores[key] = float(sev)
             clashes.add(key)
+            severity_scores[key] = E_pair
 
     return sorted(clashes), severity_scores
 
@@ -3261,8 +3297,7 @@ def _ensure_bond(i, j, order=1):
 
     # 5) write back, rebuild OBMol bonds
     core3D.bo_dict = bo
-    core3D.convert2OBMol(force_clean=True)
-    replace_bonds(core3D.OBMol, core3D.bo_dict)
+    core3D = sync_obmol_from_bodict(core3D)
 
 
 def map_ligand_local_to_core_indices_by_range(core_before_count, core_after_count, ligand_len):
@@ -3393,8 +3428,7 @@ def _ensure_bond(i, j, order=1):
 
     # write back and sync to OBMol
     core3D.bo_dict = bo
-    core3D.convert2OBMol(force_clean=True)
-    replace_bonds(core3D.OBMol, core3D.bo_dict)
+    core3D = sync_obmol_from_bodict(core3D)
 
 
 def reapply_all_haptics_and_sync(core3D, bond_order=1, prefer_nearest_metal=True):
diff --git a/molSimplify/utils/openbabel_helpers.py b/molSimplify/utils/openbabel_helpers.py
index cf80a128..e23da4e8 100644
--- a/molSimplify/utils/openbabel_helpers.py
+++ b/molSimplify/utils/openbabel_helpers.py
@@ -2,7 +2,15 @@
 from molSimplify.Classes.mol3D import mol3D
 import numpy as np
 
-def constrained_forcefield_optimization(
+def count_aromatic(obmol):
+    nb = 0
+    for bond in openbabel.OBMolBondIter(obmol):
+        if bond.IsAromatic():
+            nb += 1
+    return nb
+
+
+def constrained_forcefield_optimization_bu(
     mol,
     fixed_atom_indices=None,
     max_steps=250,
@@ -25,13 +33,26 @@ def constrained_forcefield_optimization(
     from openbabel import openbabel
     import numpy as np
 
+    print("aromatic BEFORE FF:", count_aromatic(mol.OBMol))
+
+    # --- hard assert: if bo_dict says aromatic but OBMol isn't, fix it right here ---
+    if hasattr(mol, "bo_dict"):
+        n_ar_tokens = sum(1 for v in mol.bo_dict.values() if str(v).lower() in ("ar", "am"))
+        if n_ar_tokens > 0 and count_aromatic(mol.OBMol) == 0:
+            replace_bonds(mol.OBMol, mol.bo_dict)
+            mol.OBMol.FindRingAtomsAndBonds()
+            # optional: print once
+            # print("reapplied aromaticity; aromatic bonds now:", count_aromatic(mol.OBMol))
+
+
+
     ff = openbabel.OBForceField.FindForceField(ff_name)
     if ff is None:
         raise RuntimeError(f"Forcefield '{ff_name}' not found.")
 
     obmol = mol.OBMol
     obmol.FindRingAtomsAndBonds()
-    openbabel.OBAtomTyper().AssignTypes(obmol)
+    #openbabel.OBAtomTyper().AssignTypes(obmol)
 
     if isolate_vdw:
         for a in openbabel.OBMolAtomIter(obmol):
@@ -150,11 +171,267 @@ def _apply_coords(xyz):
     if return_vdw_energy:
         results.append(ff.GetVDWEnergy())
 
+
+    print("aromatic AFTER FF:", count_aromatic(mol.OBMol))
+
     return tuple(results) if len(results) > 1 else results[0]
 
+def constrained_forcefield_optimization(
+    mol,
+    fixed_atom_indices=None,
+    max_steps=250,
+    ff_name="mmff94",
+    return_per_atom_nonbonded=False,   # deprecated path (delete-one-atom)
+    return_vdw_energy=False,
+    *,
+    return_per_atom_ff_force=False,    # preferred embedding
+    fd_delta=1e-3,                     # Å, finite-difference step
+    isolate_vdw=False,                 # True -> zero partial charges to emphasize vdW
+    freeze_aromatic_rings=True,        # NEW: freeze aromatic ring atoms (from bo_dict)
+    fallback_forcefields=("uff", "ghemical"),  # NEW: auto fallback if ff_name fails
+    verbose_ff=False                   # NEW: extra prints for debugging
+):
+    """
+    Run an OpenBabel forcefield optimization on mol.OBMol, optionally freezing atoms.
 
+    Haptics-aware behavior:
+      - If mol has _haptic_groups_global and caller requested an unconstrained relax
+        (fixed_atom_indices is None/empty), we automatically freeze metal atoms only.
+      - Rotor search is skipped when metals or haptics are present (prevents "peeling").
 
-def bond_order_from_str(bo_str):
+    NEW robustness:
+      - Re-assert aromaticity from mol.bo_dict if OBMol has 0 aromatic bonds.
+      - Freeze aromatic ring atoms (optional) to keep rings planar/rigid for steric screening.
+      - Auto fallback to UFF/Ghemical if requested FF cannot Setup() (MMFF94 often fails on TMs).
+    """
+    from openbabel import openbabel
+    import numpy as np
+
+    # -------------------- small local helpers --------------------
+    def _count_aromatic(obmol):
+        nb = 0
+        for b in openbabel.OBMolBondIter(obmol):
+            try:
+                if b.IsAromatic():
+                    nb += 1
+            except Exception:
+                pass
+        return nb
+
+    def _aromatic_atom_indices_from_bodict(bo_dict):
+        ar_atoms = set()
+        for (i, j), bo in bo_dict.items():
+            s = str(bo).strip().lower()
+            if s in ("ar", "am", "5", "1.5"):
+                ar_atoms.add(int(i))
+                ar_atoms.add(int(j))
+        return sorted(ar_atoms)
+
+    def _find_forcefield(name_list):
+        for nm in name_list:
+            ff = openbabel.OBForceField.FindForceField(nm)
+            if ff is not None:
+                return nm, ff
+        return None, None
+
+    # -------------------- sanity: aromaticity present in OBMol --------------------
+    obmol = mol.OBMol
+
+    # If bo_dict has aromatic tokens but OBMol has none, reapply bonds+aromatic flags
+    if hasattr(mol, "bo_dict"):
+        bo = mol.bo_dict
+        n_ar_tokens = sum(1 for v in bo.values() if str(v).strip().lower() in ("ar", "am", "5", "1.5"))
+        if n_ar_tokens > 0 and _count_aromatic(obmol) == 0:
+            if verbose_ff:
+                print(f"[FF] Reapplying aromaticity from bo_dict (ar_edges={n_ar_tokens})")
+            replace_bonds(obmol, bo)
+            if hasattr(obmol, "FindRingAtomsAndBonds"):
+                obmol.FindRingAtomsAndBonds()
+            if verbose_ff:
+                print("[FF] aromatic after reapply:", _count_aromatic(obmol))
+
+    # Some OB builds benefit from ring finding before typing/setup
+    if hasattr(obmol, "FindRingAtomsAndBonds"):
+        obmol.FindRingAtomsAndBonds()
+
+    # -------------------- choose a forcefield (with fallback) --------------------
+    ff_try_order = (ff_name,) + tuple(fallback_forcefields or ())
+    chosen_name, ff = _find_forcefield(ff_try_order)
+    if ff is None:
+        raise RuntimeError(f"No usable forcefield found. Tried: {ff_try_order}")
+
+    if isolate_vdw:
+        for a in openbabel.OBMolAtomIter(obmol):
+            a.SetPartialCharge(0.0)
+
+    if not ff.Setup(obmol):
+        # Try fallbacks if the chosen one failed setup
+        if verbose_ff:
+            print(f"[FF] ff.Setup('{chosen_name}') failed; trying fallbacks…")
+        ok = False
+        for nm in tuple(fallback_forcefields or ()):
+            ff2 = openbabel.OBForceField.FindForceField(nm)
+            if ff2 is not None and ff2.Setup(obmol):
+                ff = ff2
+                chosen_name = nm
+                ok = True
+                break
+        if not ok:
+            raise RuntimeError(f"Failed to set up forcefield on molecule. Tried: {ff_try_order}")
+
+    if verbose_ff:
+        print(f"[FF] Using forcefield: {chosen_name}")
+
+    # -------------------- haptics-aware defaults --------------------
+    groups = getattr(mol, "_haptic_groups_global", None) or []
+    has_haptics = bool(groups)
+
+    # Detect metals (prefer mol3D metal finder; fall back to OB atomic numbers)
+    try:
+        metals = mol.findMetal(transition_metals_only=False) or []
+    except Exception:
+        tm_atomic_nums = {
+            21,22,23,24,25,26,27,28,29,30,
+            39,40,41,42,43,44,45,46,47,48,
+            57,72,73,74,75,76,77,78,79,80
+        }
+        metals = []
+        for i, a in enumerate(openbabel.OBMolAtomIter(obmol)):
+            if a.GetAtomicNum() in tm_atomic_nums:
+                metals.append(i)
+    has_metal = bool(metals)
+
+    # Normalize fixed list
+    if fixed_atom_indices:
+        effective_fixed = [int(i) for i in fixed_atom_indices]
+    else:
+        effective_fixed = None
+
+    # If caller asked for "fully free" but haptics exist, freeze metal(s) only.
+    if (effective_fixed is None) and has_haptics and has_metal:
+        effective_fixed = sorted(set(int(i) for i in metals))
+
+    # NEW: optionally freeze aromatic ring atoms to keep rings planar/rigid
+    if freeze_aromatic_rings and hasattr(mol, "bo_dict"):
+        ar_atoms = _aromatic_atom_indices_from_bodict(mol.bo_dict)
+        if ar_atoms:
+            if effective_fixed is None:
+                effective_fixed = []
+            effective_fixed = sorted(set(effective_fixed) | set(ar_atoms))
+            if verbose_ff:
+                print(f"[FF] Freezing aromatic atoms: {len(ar_atoms)}")
+
+    # -------------------- constraints (always set; avoids sticky constraints) --------------------
+    constraints = openbabel.OBFFConstraints()
+    if effective_fixed:
+        for idx in effective_fixed:
+            constraints.AddAtomConstraint(int(idx) + 1)  # OB is 1-based
+    ff.SetConstraints(constraints)
+
+    # -------------------- optimization schedule --------------------
+    sd_steps = min(500, max_steps // 4 if max_steps >= 4 else max_steps)
+    cg_steps = max(1, max_steps - sd_steps)
+
+    ff.SteepestDescent(sd_steps)
+
+    # Rotor search can unstick torsions, but it's risky for TM complexes / haptics.
+    if (not has_haptics) and (not has_metal):
+        try:
+            ff.WeightedRotorSearch(50, 25)
+        except Exception:
+            pass
+
+    ff.ConjugateGradients(cg_steps)
+    ff.GetCoordinates(obmol)
+
+    optimized_coords = np.array(
+        [[a.GetX(), a.GetY(), a.GetZ()] for a in openbabel.OBMolAtomIter(obmol)],
+        dtype=float
+    )
+
+    # Early exit
+    if not (return_per_atom_nonbonded or return_vdw_energy or return_per_atom_ff_force):
+        return optimized_coords
+
+    results = [optimized_coords]
+
+    if return_per_atom_nonbonded:
+        # backwards-compat (delete-one-atom)
+        base_energy = ff.Energy()
+        per_atom_nonbonded = []
+        for i in range(obmol.NumAtoms()):
+            tempmol = openbabel.OBMol(obmol)
+            tempmol.DeleteAtom(tempmol.GetAtom(i + 1))
+            fftemp = openbabel.OBForceField.FindForceField(chosen_name)
+            if not fftemp or not fftemp.Setup(tempmol):
+                # try fallbacks here too
+                ok = False
+                for nm in tuple(fallback_forcefields or ()):
+                    ff2 = openbabel.OBForceField.FindForceField(nm)
+                    if ff2 is not None and ff2.Setup(tempmol):
+                        fftemp = ff2
+                        ok = True
+                        break
+                if not ok:
+                    raise RuntimeError("Failed to set up temporary forcefield.")
+            e = fftemp.Energy()
+            per_atom_nonbonded.append(base_energy - e)
+        results.append(per_atom_nonbonded)
+
+    if return_per_atom_ff_force:
+        # finite-difference gradient magnitude per atom -> "steric pressure"
+        N = optimized_coords.shape[0]
+        force_mag = np.zeros(N, dtype=float)
+
+        def _apply_coords(xyz):
+            for k, atom in enumerate(openbabel.OBMolAtomIter(obmol)):
+                atom.SetVector(float(xyz[k, 0]), float(xyz[k, 1]), float(xyz[k, 2]))
+            ff.SetCoordinates(obmol)
+
+        _apply_coords(optimized_coords)
+
+        fixed_set = set(int(i) for i in (effective_fixed or []))
+
+        for i in range(N):
+            # If atom is constrained/frozen, report 0 force (prevents misleading "pressure")
+            if i in fixed_set:
+                force_mag[i] = 0.0
+                continue
+
+            gi = np.zeros(3, dtype=float)
+            for d in range(3):
+                x_f = optimized_coords.copy()
+                x_f[i, d] += fd_delta
+                _apply_coords(x_f)
+                Ef = ff.Energy()
+
+                x_b = optimized_coords.copy()
+                x_b[i, d] -= fd_delta
+                _apply_coords(x_b)
+                Eb = ff.Energy()
+
+                gi[d] = (Ef - Eb) / (2.0 * fd_delta)
+
+            force_mag[i] = float(np.linalg.norm(gi))
+
+        # restore minimized coords
+        _apply_coords(optimized_coords)
+        results.append(force_mag)
+
+    if return_vdw_energy:
+        # Some FFs may not implement VDWEnergy; guard it.
+        vdw = None
+        try:
+            vdw = ff.GetVDWEnergy()
+        except Exception:
+            pass
+        results.append(vdw)
+
+    return tuple(results) if len(results) > 1 else results[0]
+
+
+
+def bond_order_from_str_bu(bo_str):
     """
     Convert bond order representation to OpenBabel-compatible values.
     - '1', '2', '3' ? integer bond orders
@@ -175,7 +452,23 @@ def bond_order_from_str(bo_str):
     else:
         raise ValueError(f"Unknown bond order string: {bo_str}")
 
-def replace_bonds(obmol, bo_dict):
+def bond_order_from_str(bo):
+    # numeric tokens can encode aromatic (common in pipelines)
+    if isinstance(bo, (int, float)):
+        v = float(bo)
+        if abs(v - 5.0) < 1e-6 or abs(v - 1.5) < 1e-6:
+            return 1, True      # aromatic -> order 1 + aromatic flag
+        return int(round(v)), False
+
+    s = str(bo).strip().lower()
+    if s in {"ar", "am", "5", "1.5"}:
+        return 1, True
+    if s in {"1", "2", "3"}:
+        return int(s), False
+    raise ValueError(f"Unknown bond token: {bo!r}")
+
+
+def replace_bonds_bu(obmol, bo_dict):
     """
     Replaces all bonds in the OBMol object based on a bond order dictionary.
 
@@ -203,8 +496,63 @@ def replace_bonds(obmol, bo_dict):
             bond.SetAromatic(True)
 
     # Recalculate connectivity and bonding
-    obmol.ConnectTheDots()
-    obmol.PerceiveBondOrders()
+    #obmol.ConnectTheDots()
+    #obmol.PerceiveBondOrders()
+
+
+from openbabel import openbabel
+
+_AR_TOKENS = {"ar", "am", "aro", "aromatic"}
+
+def _is_aromatic_token(bo):
+    # normalize strings
+    if isinstance(bo, str):
+        s = bo.strip().lower()
+        return s in _AR_TOKENS or s == "5" or s == "5.0" or s == "1.5"
+    # numeric conventions sometimes used for aromatic
+    if isinstance(bo, (int, float)):
+        return (bo == 5) or (abs(float(bo) - 1.5) < 1e-6)
+    return False
+
+def _numeric_bond_order(bo):
+    # for OBMol we keep aromatic as single + aromatic flag
+    if _is_aromatic_token(bo):
+        return 1
+    if isinstance(bo, str):
+        s = bo.strip()
+        if s.isdigit():
+            return int(s)
+        try:
+            return int(float(s))
+        except Exception:
+            pass
+    if isinstance(bo, (int, float)):
+        return int(round(float(bo)))
+    raise ValueError(f"Unrecognized bond order token: {bo!r}")
+
+def replace_bonds(obmol, bo_dict):
+    for bond in list(openbabel.OBMolBondIter(obmol)):
+        obmol.DeleteBond(bond)
+
+    aromatic_atoms = set()
+
+    for (a1, a2), bo in bo_dict.items():
+        i, j = a1 + 1, a2 + 1
+        order, is_ar = bond_order_from_str(bo)
+        obmol.AddBond(i, j, order)
+        if is_ar:
+            b = obmol.GetBond(i, j)
+            if b:
+                b.SetAromatic(True)
+            aromatic_atoms.add(i); aromatic_atoms.add(j)
+
+    for idx in aromatic_atoms:
+        a = obmol.GetAtom(idx)
+        if a:
+            a.SetAromatic(True)
+
+    obmol.FindRingAtomsAndBonds()
+
 
 def get_all_bonds(obmol):
     """