validate: charged ligand → warning (with biologist-actionable advice), not error

Previously a charged ligand (formal_charge != 0) hard-FAILED
validation, which blocked any bench run on ionic carboxylates,
ammoniums, etc. — a huge swath of real medicinal chemistry. The
'Phase-2 item' message offered no path forward.

The actual physics: PME on a net-charged box adds a uniform
neutralising background, which is a finite-size artefact that
corrects out via the Rocklin 2013 5-term correction (EMP + NET
+ RIP + DSC + ANA). Without that correction, ABSOLUTE ΔG_bind
on a charged ligand carries a 5-15 kcal/mol systematic error
that depends on box size and protein charge. BUT — and this is
the biologist-actionable point — it CANCELS in ΔΔG within a
same-charge congeneric series. So Kendall-τ rank-correlation
gates remain valid for ionic ligand series; only absolute MAE
is meaningful when the correction is applied.

New validator behaviour:
  ? acetate: formal charge -1 — absolute ΔG_bind carries a
    5-15 kcal/mol PBC self-interaction error (Rocklin 2013
    J Chem Phys 139:184103). For ΔΔG within a same-charge
    congeneric series the error cancels; for absolute ΔG
    apply Rocklin post-hoc.

Warning, exit 0 — biologist proceeds informed.

Also added a stub comment block where the analytical
_rocklin_pbc_correction would live, documenting the 5 terms
and why the full implementation is Phase-3 (needs protein
charge + ion-pairing distribution as auxiliary state).

Test test_charged_ligand_flagged updated to assert the new
warning shape including the 'ΔΔG cancels' caveat — that
sentence is the load-bearing biologist-actionability claim.

19/19 validator smoke tests still PASS.

Author: Rockman6

src/fep/binding.py

@@ -396,6 +396,28 @@ def _harmonic_restraint_free_energy_kcalmol(
     return dG_kJ / 4.184
 
 
+# Charged-ligand finite-size PBC correction (Rocklin 2013 J Chem
+# Phys 139:184103) is NOT yet implemented. The full correction is
+# 5 terms (EMP + NET + RIP + DSC + ANA) and requires knowing the
+# protein's net charge + ion pairing distribution — too much
+# auxiliary state to derive from the YAML alone in a way I trust
+# to ship publication-grade numbers.
+#
+# Practical guidance for biologists hitting the validator's
+# 'formal_charge ≠ 0' message:
+#   - For a CONGENERIC SERIES of equally-charged ligands (e.g.
+#     all-anionic carboxylates), the PBC self-interaction error
+#     cancels in ΔΔG. Rank-correlation gates (Kendall τ) are still
+#     valid; absolute MAE is offset by a constant.
+#   - For absolute ΔG_bind on a charged ligand vs neutral protein,
+#     extend the Rocklin correction post-hoc per the cited paper.
+#     A typical -1 ligand in a 5 nm cubic TIP3P box adds ~5-15
+#     kcal/mol depending on protein charge.
+#   - When the Phase-3 implementation lands, this comment block
+#     gets replaced with an analytical _rocklin_pbc_correction
+#     function + autowiring into compute_absolute_binding_dg.
+
+
 def _build_complex_alchemical_system(
     smiles: str,
     receptor_pdb: str | Path,
@@ -1795,15 +1817,23 @@ def _similar(a, b):
                         "est_sampled_wall_cpu_s": sampled_cpu_s,
                         "est_sampled_wall_gpu_s": sampled_gpu_s,
                     })
-                    # Biologist gotchas: formal charge ≠ 0 means the
-                    # alchemical FEP needs a counter-ion correction
-                    # we haven't wired; flag it.
+                    # Charged ligand: PBC self-interaction error in
+                    # absolute ΔG. Cancels in ΔΔG within a same-charge
+                    # series (so still valid for rank-correlation
+                    # gates) but makes absolute MAE off by 5-15
+                    # kcal/mol depending on protein charge. Warning,
+                    # not error — biologists studying ionic ligand
+                    # series shouldn't be blocked, just informed.
                     if formal_charge != 0:
-                        issues.append(
+                        warnings.append(
                             f"{name}: formal charge "
-                            f"{formal_charge:+d} — alchemical FEP "
-                            "needs a Rocklin / Warren correction "
-                            "(not wired in Phase-1); Phase-2 item")
+                            f"{formal_charge:+d} — absolute ΔG_bind "
+                            "carries a 5-15 kcal/mol PBC self-"
+                            "interaction error (Rocklin 2013 J Chem "
+                            "Phys 139:184103). For ΔΔG within a "
+                            "same-charge congeneric series the error "
+                            "cancels; for absolute ΔG apply Rocklin "
+                            "post-hoc.")
                     if n_heavy > 60:
                         issues.append(
                             f"{name}: {n_heavy} heavy atoms > 60 — "

tests/fep/test_validate_smoke.py

@@ -203,8 +203,11 @@ def test_pinned_stereo_passes():
 
 
 def test_charged_ligand_flagged():
-    """Formal-charge != 0 needs a Rocklin/Warren correction we
-    haven't wired. Validator must flag so the biologist knows."""
+    """Formal-charge != 0 carries a PBC self-interaction error in
+    absolute ΔG (Rocklin 2013). Validator emits a WARNING (not
+    a hard error) so biologists studying ionic series aren't
+    blocked — the error cancels in same-charge ΔΔG.
+    """
     yaml = """
     receptor:
       pdb_path: benchmarks/dock/1stp.pdb
@@ -215,8 +218,14 @@ def test_charged_ligand_flagged():
         dG_bind_kcalmol: -5.0
     """
     rc, out = _run(yaml)
-    assert rc != 0, "charged ligand should be flagged as an issue"
-    assert "formal charge" in out.lower() or "rocklin" in out.lower()
+    # Warning, not hard fail — proceed.
+    assert rc == 0, (
+        f"charged ligand should warn but pass; got rc={rc}\n{out}")
+    assert "formal charge" in out.lower()
+    assert "rocklin" in out.lower()
+    assert "PBC self-interaction error" in out
+    assert "ΔΔG within a same-charge" in out, (
+        "biologist actionability: must mention ΔΔG-cancels caveat")
 
 
 def test_duplicate_names_flagged():