VQE — Technical Protocol & Benchmark Report

1. Executive summary

This VQE finds the lowest-energy 3D conformation of a short peptide on a tetrahedral lattice with a Miyazawa–Jernigan (MJ) contact potential. It is a faithful reimplementation of Robert et al. 2021: the exact Hamiltonian builder, a RealAmplitudes hardware-efficient ansatz, a CVaR cost, and a Differential-Evolution (DE) classical optimizer. Two facts drive everything below:

• Energy is reproduced exactly at the paper-narrow CVaR tail (α = 0.001): APRLRFY → −1.425, matching exact diagonalization to 0.000.
• Ground-state probability P(GS) — the concentration of measurement shots on the true ground state — does not reach the paper's >20% on harder peptides with a single fixed α. The fix is α-annealing (§6), which lifts P(GS) from ~6% to 44–61% while keeping energy exact.

2. Physical model & Hamiltonian

Built by build_qubit_op(sequence). Because H is diagonal, the diagonal is extracted in O(N_terms · 2ⁿ) without a dense matrix (compute_diagonal), giving the exact classical reference for free.

3. Algorithm pipeline

sequence ──► build_qubit_op() ──► qubit_op (diagonal SparsePauliOp)
                  │
                  ├─► compute_diagonal() ──► H_diag ──► classical_energy = min(H_diag)   [exact reference]
                  │
ansatz = RealAmplitudes(n_qubits, reps=2, entanglement="full")
sampled = ansatz.copy(); sampled.measure_all()
                  │
cost = cvar_cost(sampler, sampled, H_diag, alpha, shots=1024)   # CVaR over sampled energies
optimizer = DifferentialEvolution(pop=5·reps·n_qubits, gens=80, restarts=3, F=0.5, CR=0.7)
executor = VariationalExecutor(ansatz, cost, optimizer, seed=42)
best_params, best_energy = executor.run()
                  │
final_counts = sampler.run(sampled, best_params, shots=10240)
P(GS) = Σ counts where |H_diag[bit] - classical_energy| < 1e-6   /  total   (degeneracy-aware)

The quantum part is genuine gate-level simulation: a parameterized RealAmplitudes circuit is transpiled and executed on qiskit-aer (statevector method, GPU via cuStateVec or CPU). No numpy shortcut — the sampler returns real measurement counts that index H_diag.

4. Hyperparameters (verbatim from production)

5. CVaR objective & the three α-modes

Instead of the mean energy, the cost is the Conditional Value-at-Risk: sort the per-shot energies ascending, keep the lowest α-fraction, average them. Small α rewards the optimizer for putting anyamplitude on low-energy states (good for finding the GS); large α rewards putting most amplitude there (good for concentrating P(GS)).

CVaR_alpha(energies) = mean( sort(energies)[: ceil(alpha * shots)] )

eval α — for reporting, the final CVaR energy is always re-evaluated at the paper-narrow α = 0.001 regardless of the training α, so reported energies are comparable across modes.

6. α-annealing (the P(GS) fix)

Problem: finding the GS and concentrating probability on it are coupled in a single fixed-α run. A fresh DE at wide α cannot find the GS on a rugged landscape (energy diverges); a fresh DE at narrow α finds it but has no concentration pressure (P(GS) ≈ 6–12%).

Fix — decouple the two:

Phase 1  FIND        DE @ α=0.05, full paper budget          → lands on GS (energy exact)
Phase 2  CONCENTRATE  warm-start a DE population AT the Phase-1 optimum,
                      then optimize @ α=0.25, then @ α=0.40.
                      Population already sits on the GS, so wide α only SHARPENS
                      P(GS) — it does not have to re-find the state.

The warm-started DE (_warm_de) is the only non-paper logic: the project's stock DifferentialEvolution always inits from random, so Phase 2 seeds a Gaussian cloud (σ=0.4) around the Phase-1 parameters and keeps the exact warm point as member 0. Result on DRVYIHPF (13q): P(GS) 6.0% → 24.2% → 42.7% → 60.6% across α = 0.05/0.25/0.40/0.60, energy pinned at −1.925 throughout.

7. Code map

Sampler contract (the only shim): .run(circuits, parameter_values, shots) → list[dict[bitstring,int]]. The Aer adapter caches transpiled circuits, binds parameters per eval, and runs statevector on CPU or GPU (cuStateVec_enable=True on GPU).

8. Reproduction environment (fusion4)

9. Exact reproduction commands

# single run (paper mode), GPU
python run_vqe_exact.py --sequence APRLRFY --device GPU

# concentrate mode (wide α), CPU
python run_vqe_exact.py --sequence APRLRFY --device CPU --alpha 0.25 --eval-alpha 0.001

# two-phase α-anneal
python run_vqe_anneal.py --sequence DRVYIHPF --device CPU --alpha-find 0.05 --alpha-conc 0.25,0.4

# auto: ascend α, energy-gated, until P(GS) >= 20%
python run_vqe_exact.py --sequence DRVYIHPF --device GPU --auto-target 20

# THE BATCH (10 peptides × 3 modes), capped, crash-safe, nohup
OMP_NUM_THREADS=16 OPENBLAS_NUM_THREADS=16 MKL_NUM_THREADS=16 \
  nohup python run_vqe_batch.py --devices CPU --max-threads 16 \
  --json-out vqe_batch_CPU.json > vqe_batch_CPU.log 2>&1 &

The 10 length-7 peptides (all 9 qubits): APRLRFY, DRVYIHP, KLAKLAK, RGDSPAS, YGGFMRF, GHKGHKG, WLLWLLW, ACDEFGH, HEAGAWG, PPGFSPF. GPU half is deferred behind a watcher (watch_gpu_then_vqe.sh) that auto-launches it when the QAOA-GPU benchmark frees the GPU, so GPU timings are uncontended.

10. Benchmark results (fusion4)

10.1 Reference runs (completed, validated)

Reading these: energy matches exactly at narrow α on both peptides. Single-phase wide α diverges on the hard 13q peptide (E −0.71 vs −1.925). Only anneal clears >20% P(GS) on the hard one (60.6%). P(GS) at narrow α is intrinsically tiny and shot-noisy — it is not a failure, it is what α = 0.001 is supposed to do.

10.2 Generalization batch — 10 arbitrary length-7 peptides

Question being answered: does an arbitrary length-7 peptide still match exact energy, and what P(GS) does each α-mode reach? Every peptide is 9 qubits (length-determined). All match-flags use tol = max(0.05, 2%·|E|).

Conclusions. (1) Energy generalizes — every arbitrary length-7 peptide reproduces exact diagonalization at the paper α (10/10). (2) Anneal generalizes — it matches energy on all 10 and concentrates P(GS) into a tight 41–49% band regardless of sequence, confirming the two-phase route is not specific to APRLRFY/DRVYIHPF. (3) Single-phase concentrate is fragile — it diverged in energy on 2/10 (KLAKLAK −0.86 vs −0.97; PPGFSPF −1.34 vs −1.44), the same wide-α failure mode documented in §6, so anneal is the recommended concentration path. CPU wall-time ≈ 4 min per 9-qubit run; peak RSS ≈ 200 MB (see §11).

11. Resource utilization (time, CPU & memory)

Every run is instrumented by a background MemMonitor thread (psutil RSS sampler) plus time.process_time() (cumulative CPU-seconds) and nvidia-smi --query-gpu=memory.used. That separates three distinct things: wall-clock (how long you wait), CPU-time(how many core-seconds the work actually consumed), and peak RAM.

11.1 Measured profile (fusion4, capped at 16 threads)

Host RAM is tiny and flat: ~200 MB at 9 qubits, ~215 MB at 13 qubits — dominated by the fixed Python/Qiskit/Aer baseline, not the statevector. CPU-time ≫ wall (≈11× at the 16-thread cap) because Aer parallelizes the statevector; the box stayed far under the 65–70% target. The DRVYIHPF-anneal row shows ~66 effective cores because that historical run was launched uncapped — the current batch is capped at 16.

11.2 Why memory is not the bottleneck (the scaling law)

A statevector simulation stores one complex128 (16 bytes) per basis state, so the quantum-state RAM is:

statevector RAM  =  16 bytes x 2^n_qubits   (+ ~200 MB fixed Python/Qiskit/Aer baseline)

   n_qubits      states        statevector       total (≈ baseline + state)
   ─────────────────────────────────────────────────────────────────────────
      9 (len-7)      512           8 KB             ~200 MB   (measured 200.8)
     13 (len-8)    8,192         128 KB             ~215 MB   (measured 214.7)
     20           1.05 M          16 MB             ~220 MB
     26          67.1 M           1 GB              ~1.2 GB
     30           1.07 B          16 GB             ~16 GB    (statevector-bound)
     34          17.2 B          256 GB             needs MPS / tensor-network

For every peptide in this study (9–13 qubits) the statevector is KB-scale — RAM is irrelevant and runtime is set by the number of cost evaluations (~21,870 circuit runs at 9q), not by memory. Memory only becomes the binding constraint past ~30 qubits, where the run switches from statevector to a matrix-product-state backend.

11.3 A caveat on the GPU VRAM numbers

Practical takeaway: for these peptide sizes, CPU is both faster and simpler (no host↔device transfer overhead), and memory — host or device — is never the limiting resource.

12. Known deviations & next steps

• P(GS) vs paper: a single fixed α does not reproduce the paper's >20% on hard peptides. α-annealing does, and is now the recommended/default concentration path (wired into build_workload).
• GPU vs CPU for small circuits: for 9-qubit circuits, GPU is slower than CPU (host↔device overhead dominates; APRLRFY paper CPU 791s vs GPU 1316s). GPU only wins at higher qubit counts.
• Warm-DE is the only non-paper logic — everything else is the project's exact code. The warm-start is a strict superset (it can only help concentration).
• Next: finish the 10-peptide CPU+GPU batch; fill 10.2; then extend anneal sweep to α=0.6 for the length-7 set to confirm the 60% ceiling generalizes beyond DRVYIHPF.