feat(fold): tensor-network protein folding backend with TEBD solver

Add fold module to arvak-proj: MPS-based protein folding simulation
using the same adaptive bond dimension machinery (sin(C/2)) that
powers the quantum circuit projector.

Pipeline: PDB → ContactMap → ANM → Commensurability → Hamiltonian → TEBD

Modules:
- pdb: Minimal Cα parser (no external PDB library)
- contact: Native contact map, crossing profile, relative contact order
- anm: Anisotropic Network Model (Hessian eigendecomposition via nalgebra)
- commensurability: sin(C/2) scoring of contacts → adaptive χ allocation
- hamiltonian: Go-model with backbone, local (Ramachandran), and LR terms
- mpo: Finite State Automaton MPO construction with pruning
- dmrg: Two-site DMRG with adaptive bond dimension (baseline solver)
- tebd: Imaginary-time TEBD with SWAP networks (primary solver)
- trajectory: Folding pathway analysis, ⟨r⟩ intermediate detection

Performance (1FME BBA, 28 residues, d=3, χ∈[4,16], release):
  TEBD: 11.8s — 86× faster than DMRG baseline

20 unit tests + 17 integration tests on real PDB data.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: Daniel Hinderink

Cargo.lock

@@ -17,3 +17,4 @@ tracing = { workspace = true }
 serde = { workspace = true }
 rayon = "1.10"
 faer = "0.24"
+nalgebra = "0.33"

crates/arvak-proj/Cargo.toml

@@ -0,0 +1,131 @@
+use nalgebra::{DMatrix, SymmetricEigen};
+
+pub struct ANMResult {
+    pub eigenvalues: Vec<f64>,
+    pub frequencies: Vec<f64>,
+    pub eigenvectors: Vec<Vec<f64>>,
+    pub n_modes: usize,
+    pub n_residues: usize,
+}
+
+impl ANMResult {
+    pub fn compute(
+        chain: &super::pdb::ProteinChain,
+        cutoff: f64,
+        gamma: f64,
+        n_modes: Option<usize>,
+    ) -> Self {
+        let n = chain.len();
+        let dim = 3 * n;
+
+        // Build Hessian matrix (3N x 3N)
+        let mut hess = DMatrix::<f64>::zeros(dim, dim);
+
+        for i in 0..n {
+            for j in i + 1..n {
+                let d = chain.distance(i, j);
+                if d > cutoff || d < 1e-10 {
+                    continue;
+                }
+
+                let ci = &chain.residues[i].coords;
+                let cj = &chain.residues[j].coords;
+                let d2 = d * d;
+
+                // Off-diagonal 3x3 block: -gamma * (r_ij (x) r_ij) / |r_ij|^2
+                for a in 0..3 {
+                    for b in 0..3 {
+                        let val = -gamma * (ci[a] - cj[a]) * (ci[b] - cj[b]) / d2;
+                        hess[(3 * i + a, 3 * j + b)] = val;
+                        hess[(3 * j + b, 3 * i + a)] = val;
+                        // Diagonal blocks: subtract off-diagonal contributions
+                        hess[(3 * i + a, 3 * i + b)] -= val;
+                        hess[(3 * j + a, 3 * j + b)] -= val;
+                    }
+                }
+            }
+        }
+
+        // Eigendecomposition
+        let eigen = SymmetricEigen::new(hess);
+        let mut eig_pairs: Vec<(f64, Vec<f64>)> = eigen
+            .eigenvalues
+            .iter()
+            .enumerate()
+            .map(|(idx, &val)| {
+                let vec: Vec<f64> = eigen.eigenvectors.column(idx).iter().copied().collect();
+                (val, vec)
+            })
+            .collect();
+
+        // Sort by eigenvalue ascending
+        eig_pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
+
+        // Skip first 6 trivial modes (translation + rotation)
+        let non_trivial: Vec<_> = eig_pairs
+            .into_iter()
+            .filter(|(val, _)| *val > 1e-6)
+            .collect();
+
+        let keep = n_modes.unwrap_or(non_trivial.len()).min(non_trivial.len());
+
+        let eigenvalues: Vec<f64> = non_trivial[..keep].iter().map(|(v, _)| *v).collect();
+        let frequencies: Vec<f64> = eigenvalues.iter().map(|v| v.sqrt()).collect();
+        let eigenvectors: Vec<Vec<f64>> =
+            non_trivial[..keep].iter().map(|(_, v)| v.clone()).collect();
+
+        ANMResult {
+            eigenvalues,
+            frequencies,
+            eigenvectors,
+            n_modes: keep,
+            n_residues: n,
+        }
+    }
+
+    pub fn residue_participation(&self, residue: usize) -> Vec<(usize, f64)> {
+        let mut participations: Vec<(usize, f64)> = self
+            .eigenvectors
+            .iter()
+            .enumerate()
+            .map(|(mode_idx, evec)| {
+                let base = 3 * residue;
+                let p = evec[base].powi(2) + evec[base + 1].powi(2) + evec[base + 2].powi(2);
+                (mode_idx, p)
+            })
+            .collect();
+        participations.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
+        participations
+    }
+
+    pub fn local_frequencies(&self, residue: usize, n: usize) -> Vec<f64> {
+        self.residue_participation(residue)
+            .into_iter()
+            .take(n)
+            .map(|(idx, _)| self.frequencies[idx])
+            .collect()
+    }
+
+    pub fn level_spacing_ratio(&self) -> f64 {
+        if self.eigenvalues.len() < 3 {
+            return 0.0;
+        }
+        let spacings: Vec<f64> = self
+            .eigenvalues
+            .windows(2)
+            .map(|w| (w[1] - w[0]).abs())
+            .collect();
+        if spacings.len() < 2 {
+            return 0.0;
+        }
+        let ratios: Vec<f64> = spacings
+            .windows(2)
+            .map(|w| {
+                let (a, b) = (w[0], w[1]);
+                let (min, max) = if a < b { (a, b) } else { (b, a) };
+                if max < 1e-15 { 0.0 } else { min / max }
+            })
+            .collect();
+        ratios.iter().sum::<f64>() / ratios.len() as f64
+    }
+}

crates/arvak-proj/src/fold/anm.rs

@@ -0,0 +1,82 @@
+pub struct CommensurabilityResult {
+    pub contact_scores: Vec<f64>,
+    pub bond_budget: Vec<f64>,
+    pub adaptive_chi: Vec<usize>,
+}
+
+impl CommensurabilityResult {
+    pub fn compute(
+        anm: &super::anm::ANMResult,
+        contacts: &super::contact::ContactMap,
+        n_local_modes: usize,
+    ) -> Self {
+        let mut contact_scores = Vec::with_capacity(contacts.contacts.len());
+
+        for contact in &contacts.contacts {
+            let freqs_i = anm.local_frequencies(contact.i, n_local_modes);
+            let freqs_j = anm.local_frequencies(contact.j, n_local_modes);
+            let parts_i = anm.residue_participation(contact.i);
+            let parts_j = anm.residue_participation(contact.j);
+
+            let mut score = 0.0;
+            let mut norm = 0.0;
+
+            for (a, &fi) in freqs_i.iter().enumerate() {
+                let pi = parts_i.get(a).map_or(0.0, |(_, p)| *p);
+                for (b, &fj) in freqs_j.iter().enumerate() {
+                    let pj = parts_j.get(b).map_or(0.0, |(_, p)| *p);
+                    let w = pi * pj;
+                    if fi > 1e-10 && fj > 1e-10 {
+                        let ratio = fi / fj;
+                        let frac = ratio - ratio.floor();
+                        let sin_val = (std::f64::consts::PI * frac).sin().abs();
+                        score += w * (1.0 - sin_val);
+                    }
+                    norm += w;
+                }
+            }
+
+            let c = if norm > 1e-15 { score / norm } else { 0.0 };
+            contact_scores.push(c.clamp(0.0, 1.0));
+        }
+
+        // Bond budget: sum of commensurability scores crossing each bond
+        let n = contacts.n_residues;
+        let mut bond_budget = vec![0.0; n.saturating_sub(1)];
+        for (idx, contact) in contacts.contacts.iter().enumerate() {
+            for k in contact.i..contact.j.min(bond_budget.len()) {
+                bond_budget[k] += contact_scores[idx];
+            }
+        }
+
+        CommensurabilityResult {
+            contact_scores,
+            bond_budget,
+            adaptive_chi: Vec::new(), // filled by to_adaptive_chi
+        }
+    }
+
+    pub fn to_adaptive_chi(&mut self, chi_min: usize, chi_max: usize) -> Vec<usize> {
+        if self.bond_budget.is_empty() {
+            return Vec::new();
+        }
+
+        let max_budget = self.bond_budget.iter().copied().fold(0.0_f64, f64::max);
+        if max_budget < 1e-15 {
+            self.adaptive_chi = vec![chi_min; self.bond_budget.len()];
+            return self.adaptive_chi.clone();
+        }
+
+        self.adaptive_chi = self
+            .bond_budget
+            .iter()
+            .map(|&b| {
+                let frac = (b / max_budget).sqrt();
+                let chi = chi_min + ((chi_max - chi_min) as f64 * frac) as usize;
+                chi.clamp(chi_min, chi_max)
+            })
+            .collect();
+
+        self.adaptive_chi.clone()
+    }
+}

crates/arvak-proj/src/fold/commensurability.rs

@@ -0,0 +1,60 @@
+pub struct Contact {
+    pub i: usize,
+    pub j: usize,
+    pub distance: f64,
+    pub seq_sep: usize,
+}
+
+pub struct ContactMap {
+    pub contacts: Vec<Contact>,
+    pub n_residues: usize,
+}
+
+impl ContactMap {
+    pub fn from_chain(chain: &super::pdb::ProteinChain, cutoff: f64, min_sep: usize) -> Self {
+        let n = chain.len();
+        let mut contacts = Vec::new();
+        for i in 0..n {
+            for j in i + min_sep..n {
+                let d = chain.distance(i, j);
+                if d <= cutoff {
+                    contacts.push(Contact {
+                        i,
+                        j,
+                        distance: d,
+                        seq_sep: j - i,
+                    });
+                }
+            }
+        }
+        ContactMap {
+            contacts,
+            n_residues: n,
+        }
+    }
+
+    pub fn contacts_crossing_bond(&self, k: usize) -> Vec<&Contact> {
+        self.contacts
+            .iter()
+            .filter(|c| c.i <= k && c.j > k)
+            .collect()
+    }
+
+    pub fn crossing_profile(&self) -> Vec<usize> {
+        let mut profile = vec![0usize; self.n_residues.saturating_sub(1)];
+        for c in &self.contacts {
+            for k in c.i..c.j.min(profile.len()) {
+                profile[k] += 1;
+            }
+        }
+        profile
+    }
+
+    pub fn relative_contact_order(&self) -> f64 {
+        if self.contacts.is_empty() || self.n_residues == 0 {
+            return 0.0;
+        }
+        let total_sep: f64 = self.contacts.iter().map(|c| c.seq_sep as f64).sum();
+        total_sep / (self.contacts.len() as f64 * self.n_residues as f64)
+    }
+}