Complete remaining tasks: LOCCDT hook, NATURaL CF, docking.py::run()

LIB/read_input.cpp (Task 1.3 read_input hook):
- Add #include "CavityDetect/CavityDetect.h"
- Add LOCCDT RNGOPT branch: "RNGOPT LOCCDT [cleft_id] [min_r] [max_r]"
  Runs CavityDetector::load_from_fa() + detect() on the loaded receptor,
  writes loccdt_cleft.pdb for inspection, calls to_flexaid_spheres() to
  get sphere* linked list, then feeds it to generate_grid() + calc_cleftic()
  — identical pipeline to LOCCEN/LOCCLF but cavity auto-detected.
  Exits with error if no clefts found or requested cleft_id not present.

LIB/NATURaL/NATURaLDualAssembly.cpp (Phase 3 CF placeholder removal):
- Add #include "../gaboom.h" for vcfunction() declaration
- compute_partial_cf(): replace -0.1*count placeholder with a real
  vcfunction() call on the current Cartesian atom coords.  Temporarily
  trims FA_->num_optres to n_grown_residues so only grown residues are
  scored; restores the original count after the call.  Returns 0.0 on
  vcfunction error rather than crashing the growth simulation.

python/flexaidds/docking.py (Phase 2 docking.py::run()):
- Add subprocess, shutil, re imports
- BindingMode.free_energy: fall back to min(pose.energy) when _cpp_mode
  is None (Python-only path used by PDB parser)
- BindingPopulation.__init__: accept optional modes list + temperature
  so run() can construct it directly
- Docking.run(binary, timeout): full implementation —
  1. _find_binary(): search PATH → build/ → ../build/
  2. subprocess.run(FlexAID config.inp) with capture + timeout
  3. Discover *_*.pdb output files by mtime order
  4. _parse_remark_pdb(): extract mode index, CF, RMSD, frequency from
     REMARK lines written by output_BindingMode(); compute Boltzmann weight
  5. Aggregate poses into BindingMode objects, sort by free_energy,
     return BindingPopulation

https://claude.ai/code/session_01PQFoVRutnkDNUQTFX1JQAw

Author: claude

LIB/NATURaL/NATURaLDualAssembly.cpp

@@ -27,6 +27,7 @@
 #include "NATURaLDualAssembly.h"
 #include "RibosomeElongation.h"
 #include "TransloconInsertion.h"
+#include "../gaboom.h"   // vcfunction / cfstr / get_apparent_cf_evalue
 
 #include <cstring>
 #include <cmath>
@@ -326,17 +327,28 @@ std::vector<DualAssemblyEngine::GrowthStep> DualAssemblyEngine::run() {
 
 // ─── compute_partial_cf ───────────────────────────────────────────────────────
 double DualAssemblyEngine::compute_partial_cf(int n_grown_residues) const {
-    // Simplified CF estimate: each grown residue contributes an attractive
-    // contact energy. Production code calls the full cffunction()/vcfunction()
-    // pipeline (ic2cf.cpp) with the current partially-grown complex.
-    if (!FA_ || !atoms_ || !residues_) return 0.0;
-
-    // Scale CF by accessible surface area proxy:
-    // more residues → more surface buried → stronger (negative) CF.
-    // Using -0.1 kcal/mol/residue as a placeholder; production uses the
-    // actual FlexAID Contact Function grid evaluation.
-    double count = std::min(n_grown_residues, n_residues_);
-    return -0.1 * count; // kcal/mol (attractive)
+    if (!FA_ || !VC_ || !atoms_ || !residues_) return 0.0;
+
+    // Temporarily mark grown residues as scorable and the rest as non-scorable
+    // by adjusting FA_->num_optres to cover only the grown portion.
+    // vcfunction() scores the current Cartesian atom coordinates directly —
+    // no IC → CC rebuild needed since DualAssembly keeps atoms_ in Cartesian.
+    const int saved_num_optres = FA_->num_optres;
+    const int n_score = std::min(n_grown_residues, n_residues_);
+
+    // Limit scoring to the first n_score residues (those grown so far).
+    // optres[] is ordered by residue index; trim the active count.
+    FA_->num_optres = std::min(n_score, saved_num_optres);
+
+    std::vector<std::pair<int,int>> intraclashes;
+    bool error = false;
+    double cf_val = vcfunction(FA_, VC_, atoms_, residues_, intraclashes, &error);
+
+    // Restore original optres count
+    FA_->num_optres = saved_num_optres;
+
+    if (error) return 0.0;
+    return cf_val; // kcal/mol (attractive values are negative)
 }
 
 // ─── compute_growth_entropy ───────────────────────────────────────────────────

LIB/read_input.cpp

@@ -1,5 +1,6 @@
 #include "flexaid.h"
 #include "fileio.h"
+#include "CavityDetect/CavityDetect.h"
 
 /***************************************************************************** 
  * SUBROUTINE read_input reads input file.
@@ -433,14 +434,62 @@ void read_input(FA_Global* FA,atom** atoms, resid** residue,rot** rotamer,gridpo
 		calc_cleftic(FA,*cleftgrid);
 
 	}else if(!strcmp(rngopt,"LOCCLF")){
-        
+
 		//RNGOPT LOCCLF filename.pdb
 		strcpy(FA->rngopt,"locclf");
 		strcpy(clf_file,&rngoptline[14]);
-        
+
 		printf("read binding-site grid <%s>\n",clf_file);
 		spheres = read_spheres(clf_file);
-        
+
+		(*cleftgrid) = generate_grid(FA,spheres,(*atoms),(*residue));
+		calc_cleftic(FA,*cleftgrid);
+
+	}else if(!strcmp(rngopt,"LOCCDT")){
+
+		// RNGOPT LOCCDT [cleft_id] [min_radius] [max_radius]
+		// Use native CavityDetector (SURFNET + AVX-512/OpenMP/Metal) to
+		// locate the largest binding cleft automatically.
+		// Optional args: cleft_id (default 1), probe radii (default 1.4–4.0 Å).
+		strcpy(FA->rngopt,"loccdt");
+
+		int   cdt_cleft_id  = 1;
+		float cdt_min_r     = 1.4f;
+		float cdt_max_r     = 4.0f;
+		sscanf(rngoptline,"%*s %d %f %f", &cdt_cleft_id, &cdt_min_r, &cdt_max_r);
+		if(cdt_cleft_id < 1) cdt_cleft_id = 1;
+
+		printf("LOCCDT: running native CavityDetector (cleft %d, probe %.2f–%.2f Å)\n",
+		       cdt_cleft_id, cdt_min_r, cdt_max_r);
+
+		{
+			cavity_detect::CavityDetector detector;
+			detector.load_from_fa(*atoms, *residue, FA->res_cnt);
+			detector.detect(cdt_min_r, cdt_max_r);
+
+			if(detector.clefts().empty()){
+				fprintf(stderr,"ERROR: LOCCDT found no clefts — "
+				        "check probe radii or fall back to LOCCLF/LOCCEN\n");
+				Terminate(2);
+			}
+
+			// Write the detected cleft as a CLF sphere PDB for inspection
+			detector.write_sphere_pdb("loccdt_cleft.pdb", cdt_cleft_id);
+			const auto& cleft_list = detector.clefts();
+			size_t cdt_idx = static_cast<size_t>(cdt_cleft_id - 1);
+			printf("LOCCDT: cleft %d detected (%zu spheres) — "
+			       "written to loccdt_cleft.pdb\n",
+			       cdt_cleft_id,
+			       cdt_idx < cleft_list.size() ? cleft_list[cdt_idx].spheres.size() : 0u);
+
+			spheres = detector.to_flexaid_spheres(cdt_cleft_id);
+		}
+
+		if(!spheres){
+			fprintf(stderr,"ERROR: LOCCDT cleft %d not found\n", cdt_cleft_id);
+			Terminate(2);
+		}
+
 		(*cleftgrid) = generate_grid(FA,spheres,(*atoms),(*residue));
 		calc_cleftic(FA,*cleftgrid);
 	}

python/flexaidds/docking.py

@@ -3,6 +3,9 @@
 Provides Pythonic API for molecular docking workflows.
 """
 
+import subprocess
+import shutil
+import re
 import numpy as np
 from pathlib import Path
 from typing import List, Optional, Dict, Any
@@ -85,6 +88,9 @@ def free_energy(self) -> float:
         """Helmholtz free energy F = H - TS (kcal/mol)."""
         if self._cpp_mode:
             return self._cpp_mode.get_free_energy()
+        # Python-only path: use best pose energy as proxy
+        if self._poses:
+            return min(p.energy for p in self._poses)
         return float('inf')
 
     @property
@@ -123,9 +129,10 @@ class BindingPopulation:
     Provides ensemble-level analysis and ranking of binding modes.
     """
 
-    def __init__(self):
-        self._modes: List[BindingMode] = []
-        self._temperature: float = 300.0
+    def __init__(self, modes: Optional[List[BindingMode]] = None,
+                 temperature: float = 300.0):
+        self._modes: List[BindingMode] = list(modes) if modes else []
+        self._temperature: float = temperature
 
     def add_mode(self, mode: BindingMode) -> None:
         """Add a binding mode to the population."""
@@ -286,21 +293,166 @@ def optimization_method(self) -> Optional[str]:
         """Optimization method, e.g. 'GA' (METOPT keyword)."""
         return self._config.get("METOPT")
 
-    def run(self, **kwargs) -> BindingPopulation:
-        """Execute docking simulation.
-        
+    def run(self, binary: Optional[str] = None,
+            timeout: int = 3600, **kwargs) -> BindingPopulation:
+        """Execute docking via the FlexAID C++ binary and parse results.
+
+        Locates the ``FlexAID`` binary (in PATH, project build/, or explicit
+        *binary* argument), invokes it with this config file, waits for
+        completion, then parses all ``*_N_M.pdb`` output files written by
+        ``output_Population()`` to reconstruct a ``BindingPopulation``.
+
+        Args:
+            binary:  Path to FlexAID executable.  If *None*, searches PATH and
+                     common build locations (``build/FlexAID``,
+                     ``../build/FlexAID``).
+            timeout: Wall-clock timeout in seconds (default 3600).
+            **kwargs: Ignored; reserved for future keyword overrides.
+
         Returns:
-            BindingPopulation with ranked binding modes
-        
-        Note:
-            Full implementation requires integration with C++ FlexAID GA engine.
-            This is a stub for Phase 2 development.
+            BindingPopulation populated from the PDB REMARK lines written by
+            ``output_BindingMode()`` / ``output_Population()``.
+
+        Raises:
+            FileNotFoundError: binary not found.
+            RuntimeError:      FlexAID exited non-zero or produced no output.
         """
-        raise NotImplementedError(
-            "Full docking pipeline integration in progress. "
-            "For now, use C++ FlexAID binary directly and parse output with "
-            "BindingMode/BindingPopulation wrappers."
+        # ── 1. Locate binary ─────────────────────────────────────────────────
+        exe = self._find_binary(binary)
+
+        # ── 2. Invoke FlexAID ────────────────────────────────────────────────
+        cmd = [str(exe), str(self.config_file)]
+        try:
+            result = subprocess.run(
+                cmd,
+                cwd=self.config_file.parent,
+                capture_output=True,
+                text=True,
+                timeout=timeout,
+            )
+        except subprocess.TimeoutExpired as exc:
+            raise RuntimeError(
+                f"FlexAID timed out after {timeout}s"
+            ) from exc
+
+        if result.returncode != 0:
+            raise RuntimeError(
+                f"FlexAID exited with code {result.returncode}.\n"
+                f"stdout: {result.stdout[-2000:]}\n"
+                f"stderr: {result.stderr[-2000:]}"
+            )
+
+        # ── 3. Discover output PDBs ──────────────────────────────────────────
+        # output_BindingMode writes files named <prefix>_<minPoints>_<mode>.pdb
+        # Collect all candidate PDB files in the working directory.
+        work_dir = self.config_file.parent
+        pdb_files = sorted(work_dir.glob("*_*.pdb"),
+                           key=lambda p: p.stat().st_mtime)
+
+        if not pdb_files:
+            raise RuntimeError(
+                "FlexAID completed but no PDB output files were found in "
+                f"{work_dir}. Check the config file NRGOUT / output settings."
+            )
+
+        # ── 4. Parse PDB REMARK lines into BindingModes ──────────────────────
+        temperature = self._config.get("TEMPER", 300) or 300
+        modes: List[BindingMode] = []
+        seen_modes: Dict[int, BindingMode] = {}
+
+        for pdb_path in pdb_files:
+            mode_info = self._parse_remark_pdb(pdb_path, temperature)
+            if mode_info is None:
+                continue
+            mode_idx, pose = mode_info
+            if mode_idx not in seen_modes:
+                seen_modes[mode_idx] = BindingMode()
+            seen_modes[mode_idx]._poses.append(pose)
+
+        # Sort modes by free energy (ascending → most favourable first)
+        modes = sorted(seen_modes.values(),
+                       key=lambda m: m.free_energy)
+
+        return BindingPopulation(modes, temperature=float(temperature))
+
+    # ── helpers ───────────────────────────────────────────────────────────────
+
+    def _find_binary(self, binary: Optional[str]) -> Path:
+        """Locate the FlexAID executable."""
+        if binary is not None:
+            p = Path(binary)
+            if not p.is_file():
+                raise FileNotFoundError(f"Specified FlexAID binary not found: {binary}")
+            return p
+
+        # Search order: PATH → project-relative build dirs
+        in_path = shutil.which("FlexAID")
+        if in_path:
+            return Path(in_path)
+
+        candidates = [
+            self.config_file.parent / "FlexAID",
+            self.config_file.parent / "build" / "FlexAID",
+            Path(__file__).parents[3] / "build" / "FlexAID",
+        ]
+        for c in candidates:
+            if c.is_file():
+                return c
+
+        raise FileNotFoundError(
+            "FlexAID binary not found in PATH or build/. "
+            "Build with 'cmake --build build' or pass binary= argument."
+        )
+
+    @staticmethod
+    def _parse_remark_pdb(
+            pdb_path: Path, temperature: float) -> Optional[tuple]:
+        """Parse a single output PDB written by output_BindingMode().
+
+        Extracts mode index, CF, RMSD, and per-pose energy from REMARK lines.
+        Returns (mode_index, Pose) or None if the file lacks FlexAID remarks.
+        """
+        mode_idx   = None
+        cf_val     = None
+        rmsd_val   = None
+        freq       = 1
+
+        try:
+            text = pdb_path.read_text(errors="replace")
+        except OSError:
+            return None
+
+        for line in text.splitlines():
+            if not line.startswith("REMARK"):
+                continue
+            # "REMARK Binding Mode:N Best CF in Binding Mode:X …"
+            m = re.search(
+                r"Binding Mode:(\d+).*?Best CF in Binding Mode:\s*([-\d.]+)"
+                r".*?Binding Mode Frequency:(\d+)",
+                line)
+            if m:
+                mode_idx = int(m.group(1))
+                cf_val   = float(m.group(2))
+                freq     = int(m.group(3))
+            # "REMARK 0.12345 RMSD to ref. structure …"
+            m2 = re.search(r"REMARK\s+([\d.]+)\s+RMSD to ref\.", line)
+            if m2 and rmsd_val is None:
+                rmsd_val = float(m2.group(1))
+
+        if mode_idx is None or cf_val is None:
+            return None
+
+        import math
+        beta = 1.0 / (0.001987206 * float(temperature))
+        bw   = math.exp(-beta * cf_val)
+
+        pose = Pose(
+            index=mode_idx,
+            energy=cf_val,
+            rmsd=rmsd_val,
+            boltzmann_weight=bw,
         )
+        return mode_idx, pose
 
     def __repr__(self) -> str:
         return f"<Docking config={self.config_file.name}>"