add FFT-based Buccaneer-inspired C-alpha detection with ranked candidate output and validation on 6VXX and 1JYX

Author: rmeanapa

doc/implementation_notes/calpha_finder.md

@@ -0,0 +1,119 @@
+# C-alpha finder
+
+## Contract
+
+The first model-building slice adds `simple_exec prg=detect_calpha`. It locates
+oriented C-alpha seed candidates in a cryo-EM map using a built-in expected
+backbone density target.
+
+Inputs:
+
+- `vol1`: work map to search.
+- `smpd`: sampling distance in Angstroms per voxel.
+- `angstep`: coarse SO(3) angular sampling step.
+- `npeaks`: maximum number of non-overlapping candidates.
+- `thres`: minimum weighted normalized-correlation score.
+
+Outputs:
+
+- `pdbout`: candidate C-alpha coordinates in PDB format, with the score in
+  the B-factor field.
+- A same-stem CSV file containing coordinates, scores, and orientation
+  matrices. Each matrix maps the canonical target frame into the work
+  map frame.
+- `outvol`: maximum score over orientations at every work-map voxel.
+
+The initial coordinate contract follows the existing SIMPLE atom/map
+convention: voxel `(1,1,1)` corresponds to `(0,0,0)` Angstrom. MRC start and
+origin metadata are not interpreted by this slice. The FFT-searched work map
+must have even dimensions.
+
+## Numerical method
+
+1. Construct an ideal local `N-CA-C` geometry using 1.458 Angstrom `N-CA` and
+   1.525 Angstrom `CA-C` bonds and a 111.2 degree `N-CA-C` angle.
+2. Represent the three backbone atoms with Gaussian density. The C-alpha lobe
+   is slightly stronger to anchor the reported coordinate at the target center.
+   Gaussian sigma is `max(0.85 Angstrom, 0.75 * smpd)` to remain sampled on
+   coarser voxel grids.
+3. Sample this expected density inside a 4 Angstrom sphere and use a cosine
+   taper as the target reliability weight.
+4. Generate a coarse SO(3) grid from Fibonacci-sphere axes and uniformly
+   sampled roll angles.
+5. For every orientation, rotate the analytic target and evaluate weighted
+   normalized cross-correlation at every translation. Three 3-D correlations
+   provide the weighted data sum, squared-data sum, and target-data sum.
+6. Retain the best orientation per voxel and select score maxima with a
+   2 Angstrom exclusion radius.
+
+The normalized score fits a local scale and offset implicitly and is the first
+SIMPLE analogue of Buccaneer's cryo-EM correlation mode. Chain growth,
+fragment joining, sequence assignment, side-chain classification, and model
+refinement are deliberately outside this slice.
+
+## Review findings
+
+- `image%ccf` supplies the required FFT correlation primitive and preserves
+  SIMPLE's phase-origin convention. A new `image%ccf_into` variant writes into
+  an already allocated image so the search reuses FFT plans and buffers rather
+  than rebuilding them for every orientation.
+- An isolated spherical C-alpha Gaussian would make orientation search
+  redundant and provide poor discrimination. Including the expected adjacent
+  N and C density preserves a useful oriented Buccaneer-style target without
+  requiring `vol2` or a reference model.
+- The existing atom-centered validation code assumes the same zero-origin
+  coordinate convention adopted here.
+- The algorithm belongs in a new model-building domain. The commander owns
+  I/O and lifecycle; the domain object owns analytic target construction and
+  numerical search.
+
+## Validation criteria
+
+- Static Fortran checks report no overlong lines or invalid multiline
+  diagnostic macros.
+- The UI registry contains `detect_calpha` and the execution API routes it to
+  one commander.
+- `simple_test_exec test=detect_calpha` builds repeated analytic backbone
+  density in memory at an orientation in the search grid, searches it, and
+  verifies that a reported candidate is within 1.5 Angstrom of a known
+  C-alpha position.
+- User validation remains: compile the affected targets, run the synthetic
+  test, and evaluate precision/recall on representative SIMPLE maps.
+
+## Built-in molecule benchmark
+
+`simple_test_exec test=detect_calpha_molecules` is an opt-in accuracy benchmark
+using SIMPLE's hard-coded 6VXX spike-protein and 1JYX beta-galactosidase
+coordinates. The fixtures contain 2,916 and 4,044 protein C-alpha atoms,
+respectively. For each structure the benchmark:
+
+1. creates a tightly padded, even-dimension density map at the requested
+   sampling distance;
+2. runs the same analytic-target search used by `simple_exec`;
+3. performs greedy one-to-one matching of score-ordered predictions to the
+   ground truth within 2 Angstrom; and
+4. reports truth, prediction, unique-match, missed, false-positive, recall,
+   and precision counts while retaining the input map, centered model,
+   candidate PDB/CSV, and score volume.
+
+The default benchmark settings are `smpd=1.3`, `angstep=45`, and `thres=0.25`.
+They can be overridden on the test command line. No accuracy threshold is a
+test failure yet; the first observed run establishes the baseline. Candidate
+selection sorts threshold-eligible scores once before non-maximum suppression,
+avoiding the previous full-volume scan for every requested peak.
+
+### 2026-09-10 baseline
+
+The Release `simple_test_exec` benchmark was run with eight threads and the
+default settings. Matches are greedy, one-to-one, and within 2 Angstrom.
+
+| Structure | Truth | Top-N matches | Top-N recall/precision | Top-2N matches | Top-2N recall | Top-2N precision |
+| --- | ---: | ---: | ---: | ---: | ---: | ---: |
+| 6VXX | 2,916 | 1,753 | 0.601 | 2,755 | 0.945 | 0.472 |
+| 1JYX | 4,044 | 2,582 | 0.638 | 3,482 | 0.861 | 0.431 |
+
+Both structures reached the `2 * ntruth` output cap. At `thres=0.25`, 103,078
+6VXX voxels and 126,413 1JYX voxels were score-eligible, so this threshold does
+not determine the reported operating points. Follow-up calibration should use
+the score-ranked outputs to choose a precision/recall tradeoff before setting a
+production threshold.

src/main/commanders/simple/simple_commanders_atoms.f90

@@ -4,6 +4,7 @@ module simple_commanders_atoms
 use simple_nanoparticle_utils
 use simple_nanoparticle, only: nanoparticle
 use simple_atoms,        only: atoms
+use simple_calpha_finder, only: calpha_finder
 implicit none
 #include "simple_local_flags.inc"
 
@@ -42,6 +43,11 @@ module simple_commanders_atoms
     procedure :: execute      => exec_detect_atoms
 end type commander_detect_atoms
 
+type, extends(commander_base) :: commander_detect_calpha
+  contains
+    procedure :: execute      => exec_detect_calpha
+end type commander_detect_calpha
+
 type, extends(commander_base) :: commander_map2model_fsc
   contains
     procedure :: execute      => exec_map2model_fsc
@@ -416,6 +422,36 @@ subroutine exec_detect_atoms( self, cline )
         call simple_end('**** SIMPLE_DETECT_ATOMS NORMAL STOP ****')
     end subroutine exec_detect_atoms
 
+    subroutine exec_detect_calpha(self, cline)
+        class(commander_detect_calpha), intent(inout) :: self
+        class(cmdline),                 intent(inout) :: cline
+        type(parameters)    :: params
+        type(image)         :: workvol
+        type(calpha_finder) :: finder
+        integer :: work_ldim(3), nsections
+
+        if(.not.cline%defined('angstep')) call cline%set('angstep', 45)
+        if(.not.cline%defined('npeaks'))  call cline%set('npeaks', 100)
+        if(.not.cline%defined('thres'))   call cline%set('thres', 0.25)
+        if(.not.cline%defined('pdbout'))  call cline%set('pdbout', 'calpha_candidates.pdb')
+        if(.not.cline%defined('outvol'))  call cline%set('outvol', 'calpha_scores.mrc')
+        call params%new(cline)
+        if(params%smpd <= 0.) THROW_HARD('smpd must be positive for detect_calpha')
+
+        call find_ldim_nptcls(params%vols(1), work_ldim, nsections)
+        if(work_ldim(3) <= 1) THROW_HARD('vol1 must contain a 3D search volume')
+
+        call workvol%new(work_ldim, params%smpd)
+        call workvol%read(params%vols(1))
+        call finder%new(params%smpd, 4.0)
+        call finder%search(workvol, real(params%angstep), params%npeaks, params%thres, &
+            params%pdbout, params%outvol)
+
+        call finder%kill()
+        call workvol%kill()
+        call simple_end('**** SIMPLE_DETECT_CALPHA NORMAL STOP ****')
+    end subroutine exec_detect_calpha
+
     subroutine exec_map2model_fsc( self, cline )
         use simple_commanders_resolest, only: commander_fsc
         class(commander_map2model_fsc), intent(inout) :: self

src/main/commanders/test/simple_commanders_test_single.f90

@@ -13,6 +13,16 @@ module simple_commanders_test_single
     procedure :: execute      => exec_test_detect_atoms
 end type commander_test_detect_atoms
 
+type, extends(commander_base) :: commander_test_detect_calpha
+  contains
+    procedure :: execute      => exec_test_detect_calpha
+end type commander_test_detect_calpha
+
+type, extends(commander_base) :: commander_test_detect_calpha_molecules
+  contains
+    procedure :: execute      => exec_test_detect_calpha_molecules
+end type commander_test_detect_calpha_molecules
+
 type, extends(commander_base) :: commander_test_simulate_nanoparticle
   contains
     procedure :: execute      => exec_test_simulate_nanoparticle
@@ -93,6 +103,202 @@ subroutine exec_test_detect_atoms( self, cline )
     call simple_end('**** SIMPLE_TEST_DETECT_ATOMS NORMAL STOP ****')
 end subroutine exec_test_detect_atoms
 
+subroutine exec_test_detect_calpha( self, cline )
+    use simple_atoms,         only: atoms
+    use simple_calpha_finder, only: calpha_finder
+    class(commander_test_detect_calpha), intent(inout) :: self
+    class(cmdline),                      intent(inout) :: cline
+    character(len=*), parameter :: PDB_FILE = 'test_calpha_candidates.pdb'
+    character(len=*), parameter :: CSV_FILE = 'test_calpha_candidates.csv'
+    character(len=*), parameter :: MRC_FILE = 'test_calpha_scores.mrc'
+    integer,          parameter :: TEST_BOX = 32, NRES = 3
+    type(image)         :: workvol
+    type(atoms)         :: candidates
+    type(calpha_finder) :: finder
+    real(kind=c_float), pointer :: density(:,:,:)
+    real    :: centers(3,NRES), atom_sites(3,3), amplitudes(3), rotation(3,3)
+    real    :: xyz(3), site(3), delta(3), distance, closest
+    integer :: ldim(3), ires, iatom, ix, iy, iz
+
+    write(logfhandle,'(A)') '>>> TEST_DETECT_CALPHA:'
+    ldim            = [TEST_BOX,TEST_BOX,TEST_BOX]
+    centers(:,1)    = [8.,8.,8.]
+    centers(:,2)    = [16.,16.,16.]
+    centers(:,3)    = [24.,24.,24.]
+    atom_sites(:,1) = [0.,0.,0.]
+    atom_sites(:,2) = [1.458*cos(111.2*PI/180.), 1.458*sin(111.2*PI/180.), 0.]
+    atom_sites(:,3) = [1.525,0.,0.]
+    amplitudes      = [1.25,1.0,1.0]
+    rotation(:,1)   = [0.,1.,0.]
+    rotation(:,2)   = [-0.5,0.,sqrt(0.75)]
+    rotation(:,3)   = [sqrt(0.75),0.,0.5]
+    call workvol%new(ldim, 1.0)
+    call workvol%get_rmat_ptr(density)
+    density         = 0.
+    do ires         = 1, NRES
+        do iatom = 1, size(atom_sites,2)
+            site = centers(:,ires) + matmul(rotation, atom_sites(:,iatom))
+            do iz = 1, TEST_BOX
+                do iy = 1, TEST_BOX
+                    do ix = 1, TEST_BOX
+                        xyz = real([ix,iy,iz] - 1)
+                        delta = xyz - site
+                        density(ix,iy,iz) = density(ix,iy,iz) + amplitudes(iatom) * &
+                            exp(-0.5 * sum(delta * delta) / 0.85**2)
+                    enddo
+                enddo
+            enddo
+        enddo
+    enddo
+
+    call finder%new(1.0, 4.0)
+    call finder%search(workvol, 180.0, 10, 0.5, string(PDB_FILE), string(MRC_FILE))
+    if(nlines(string(PDB_FILE)) == 0) THROW_HARD('TEST_DETECT_CALPHA FAILED: no candidates')
+    call candidates%new(string(PDB_FILE))
+    closest = huge(1.)
+    do iatom = 1, candidates%get_n()
+        do ires = 1, NRES
+            distance = sqrt(sum((candidates%get_coord(iatom) - centers(:,ires))**2))
+            closest  = min(closest, distance)
+        enddo
+    enddo
+    if( closest > 1.5 ) THROW_HARD('TEST_DETECT_CALPHA FAILED: peak is displaced')
+
+    call candidates%kill()
+    call finder%kill()
+    call workvol%kill()
+    if(file_exists(PDB_FILE)) call del_file(PDB_FILE)
+    if(file_exists(CSV_FILE)) call del_file(CSV_FILE)
+    if(file_exists(MRC_FILE)) call del_file(MRC_FILE)
+    write(logfhandle,'(A,F7.3,A)') '>>> TEST_DETECT_CALPHA: PASS (closest peak ', closest, ' A)'
+    call simple_end('**** SIMPLE_TEST_DETECT_CALPHA NORMAL STOP ****')
+
+end subroutine exec_test_detect_calpha
+
+subroutine exec_test_detect_calpha_molecules( self, cline )
+    use simple_atoms,         only: atoms
+    use simple_calpha_finder, only: calpha_finder
+    use simple_molecule_data, only: molecule_data, betagal_1jyx, sars_cov2_spkgp_6vxx
+    class(commander_test_detect_calpha_molecules), intent(inout) :: self
+    class(cmdline),                                intent(inout) :: cline
+    type(parameters)    :: params
+    type(molecule_data) :: mol
+
+    if( .not.cline%defined('smpd') )    call cline%set('smpd', 1.3)
+    if( .not.cline%defined('angstep') ) call cline%set('angstep', 45)
+    if( .not.cline%defined('thres') )   call cline%set('thres', 0.25)
+    call params%new(cline)
+
+    write(logfhandle,'(A)') '>>> C-ALPHA MOLECULE BENCHMARK:'
+    write(logfhandle,'(A,F6.2,A,I0,A,F6.3)') '    smpd=', params%smpd, &
+        ' A, angstep=', params%angstep, ' degrees, threshold=', params%thres
+    mol = sars_cov2_spkgp_6vxx()
+    call evaluate_molecule('6VXX', mol, 2916, params%smpd, params%angstep, params%thres)
+    mol = betagal_1jyx()
+    call evaluate_molecule('1JYX', mol, 4044, params%smpd, params%angstep, params%thres)
+    call simple_end('**** SIMPLE_TEST_DETECT_CALPHA_MOLECULES NORMAL STOP ****')
+
+contains
+
+    subroutine evaluate_molecule( label, molecule_data_in, expected_truth, smpd, angstep, threshold )
+        character(len=*),    intent(in) :: label
+        type(molecule_data), intent(in) :: molecule_data_in
+        integer,             intent(in) :: expected_truth, angstep
+        real,                intent(in) :: smpd, threshold
+        real, parameter      :: MAP_PADDING = 12.0, MATCH_RADIUS = 2.0
+        type(atoms)          :: molecule, candidates
+        type(calpha_finder)  :: finder
+        type(image)          :: workvol
+        type(string)         :: source_file, truth_file, vol_file, candidate_file, score_file
+        real, allocatable    :: truth_xyz(:,:)
+        logical, allocatable :: truth_matched(:)
+        real    :: span(3), delta(3), best_distance_sq, recall, precision
+        real    :: recall_top_n, precision_top_n
+        integer :: ldim(3), iatom, itruth, ipred, ntruth, npred, nmatched, best_truth
+        integer :: top_n_count, nmatched_top_n
+
+        source_file    = trim(label)//'.pdb'
+        truth_file     = trim(label)//'_calpha_truth.pdb'
+        vol_file       = trim(label)//'_calpha_input.mrc'
+        candidate_file = trim(label)//'_calpha_candidates.pdb'
+        score_file     = trim(label)//'_calpha_scores.mrc'
+        span           = maxval(molecule_data_in%xyz, dim=1) - minval(molecule_data_in%xyz, dim=1)
+        ldim           = max(round2even((span + 2. * MAP_PADDING) / smpd), 16)
+        call molecule%pdb2mrc(pdbfile=source_file, volfile=vol_file, smpd=smpd, &
+            center_pdb=.true., pdb_out=truth_file, vol_dim=ldim, mol=molecule_data_in)
+
+        ntruth = 0
+        do iatom = 1, molecule%get_n()
+            if(molecule%get_name(iatom) == ' CA ' .and. molecule%get_element(iatom) == 'C ') &
+                ntruth = ntruth + 1
+        enddo
+        if(ntruth /= expected_truth) THROW_HARD('Unexpected built-in C-alpha count')
+        allocate(truth_xyz(3,ntruth), source=0.)
+        allocate(truth_matched(ntruth), source=.false.)
+        itruth = 0
+        do iatom = 1, molecule%get_n()
+            if(molecule%get_name(iatom) /= ' CA ' .or. molecule%get_element(iatom) /= 'C ') cycle
+            itruth = itruth + 1
+            truth_xyz(:,itruth) = molecule%get_coord(iatom)
+        enddo
+
+        call workvol%new(ldim, smpd)
+        call workvol%read(vol_file)
+        call finder%new(smpd, 4.0)
+        call finder%search(workvol, real(angstep), 2 * ntruth, threshold, candidate_file, score_file)
+
+        npred = 0
+        if(nlines(candidate_file) > 0)then
+            call candidates%new(candidate_file)
+            npred = candidates%get_n()
+        endif
+        nmatched       = 0
+        nmatched_top_n = 0
+        top_n_count = min(ntruth, npred)
+        do ipred = 1, npred
+            best_truth       = 0
+            best_distance_sq = huge(1.)
+            do itruth = 1, ntruth
+                if(truth_matched(itruth)) cycle
+                delta = candidates%get_coord(ipred) - truth_xyz(:,itruth)
+                if(sum(delta * delta) < best_distance_sq)then
+                    best_distance_sq = sum(delta * delta)
+                    best_truth       = itruth
+                endif
+            enddo
+            if(best_truth > 0 .and. best_distance_sq <= MATCH_RADIUS**2)then
+                truth_matched(best_truth) = .true.
+                nmatched                  = nmatched + 1
+            endif
+            if(ipred == top_n_count) nmatched_top_n = nmatched
+        enddo
+        recall_top_n    = real(nmatched_top_n) / real(ntruth)
+        precision_top_n = 0.
+        if(top_n_count > 0) precision_top_n = real(nmatched_top_n) / real(top_n_count)
+        recall    = real(nmatched) / real(ntruth)
+        precision = 0.
+        if(npred > 0) precision = real(nmatched) / real(npred)
+
+        write(logfhandle,'(A,A)') '>>> ', trim(label)
+        write(logfhandle,'(A,I0,A,I0)') '    truth=', ntruth, ', candidate cap=', 2 * ntruth
+        write(logfhandle,'(A,I0,A,I0,A,F7.3,A,F7.3)') '    top-N: predicted=', top_n_count, &
+            ', matched=', nmatched_top_n, ', recall=', recall_top_n, ', precision=', precision_top_n
+        write(logfhandle,'(A,I0,A,I0,A,I0,A,F7.3,A,F7.3)') '    top-2N: predicted=', npred, &
+            ', matched=', nmatched, ', missed=', ntruth - nmatched, ', recall=', recall, &
+            ', precision=', precision
+        write(logfhandle,'(A,3(I0,1X))') '    map dimensions=', ldim
+        write(logfhandle,'(A,A)') '    candidates: ', candidate_file%to_char()
+        write(logfhandle,'(A,A)') '    score volume: ', score_file%to_char()
+
+        if(npred > 0) call candidates%kill()
+        call finder%kill()
+        call workvol%kill()
+        call molecule%kill()
+        deallocate(truth_xyz, truth_matched)
+    end subroutine evaluate_molecule
+
+end subroutine exec_test_detect_calpha_molecules
+
 subroutine exec_test_simulate_nanoparticle( self, cline )
     use simple_commanders_sim, only: commander_simulate_nanoparticle
     class(commander_test_simulate_nanoparticle), intent(inout) :: self