[challenge]: Ab initio full-π electronic structure of fullerenes with DMRG and tree tensor networks

Released by

Guo CHEN, HKUST(GZ)

Contact email

guochen@hkust-gz.edu.cn

Method

MPS Based Algorithm

Challenge issue

The challenge

The electronic structure of the fullerenes is more subtle than these symmetric, single-element cages suggest. Density functional theory (DFT) does not reliably predict the energy ordering of the C₂₀ isomers (Grossman 1995). The singlet–triplet gap of C₃₆ is uncertain. The correlation strength of the C₆₀ π shell, the basic input to theories of fulleride superconductivity, has only been estimated (Gunnarsson 1997). The fullerenes are also a natural benchmark family for tensor networks: chemically uniform cages whose entanglement grows with size, playing the same role for closed surfaces that the polyacenes played for chains.

However, to our knowledge, no ab initio tensor-network calculation of a fullerene has been published. The many-body studies that exist use model Hamiltonians: a Hubbard model for C₂₀ (Lin 2007), a Heisenberg model on the C₆₀ geometry (Rausch 2021). This challenge is to do it with the real electronic Hamiltonian: one- and two-electron integrals over the π orbitals of the curved cage.

Consider a series of three fullerenes with increasing size: C₂₀, C₃₆, and C₆₀. C₂₀ is a CAS(20,20) problem — a complete active space (CAS) of 20 π electrons in 20 orbitals — small enough to solve near-exactly with selected CI, so you generate your own reference. C₃₆, CAS(36,36), is beyond the reach of FCI (exact diagonalization), and no DMRG result has been reported yet. C₆₀ is the stretch goal. The strained small cages have genuine diradical character (Grossman 1995), while neutral C₆₀ is only mildly multireference — the correlation strength is itself size-dependent, which is part of the story.

Cage	Active space	FCI dim	Goal
C20 (dodecahedron; Ih, Jahn–Teller distorted)	CAS(20,20)	~10¹⁰	generate the SHCI reference; validate MPS and TTNS against it
C36 (D6h)	CAS(36,36)	~10²⁰	first energy and S–T gap at this level
C60 (truncated icosahedron; Ih)	CAS(60,60)	~10³⁴	extrapolated energy (stretch)

The difficulty is that a cage cannot be unrolled into a chain. An MPS — the variational state that DMRG optimizes — arranges the orbitals along a line; a cut divides them into two blocks, and the bond dimension at the cut grows exponentially with the entanglement between the blocks. For localized π orbitals this entanglement is roughly proportional to the number of C–C bonds joining the two blocks. A linear molecule can be ordered so that every cut crosses only a few bonds, whereas a cage cannot; some cut must bisect it, crossing a ring of about √N bonds: 7 for C₂₀, 9 for C₃₆, 13 for C₆₀. A bond dimension of hundreds therefore suffices for C₂₀, C₃₆ needs thousands, and C₆₀ approaches the limit of an MPS. A tree tensor network state (TTNS) redistributes this entanglement, and the same DMRG sweeps optimize it (Task 3). The counting is only a guide, as the Coulomb interaction is long-ranged; the entanglement map of Task 1 tests it directly.

Observables: the neutral ground-state energy, the π-space correlation energy (CAS energy minus RHF), and the singlet–triplet gap.

Setup

Only the π electrons are active: one radial-p orbital per carbon. The 1s core and the σ framework are frozen, doubly occupied in every configuration. The active electrons see the Hamiltonian

$$
H = E_0 + \sum_{pq} \tilde{h}{pq} a^\dagger_p a_q + \frac{1}{2}\sum{pqrs} \langle pq|v|rs\rangle a^\dagger_p a^\dagger_q a_s a_r.
$$

The integrals, written in the Dirac notation, run over the active orbitals; $E_0$ contains the nuclear repulsion and the frozen-shell energy. The one-electron term is dressed by the frozen shells,

$$
\tilde{h}{pq} = h{pq} + \sum_{i \in \mathrm{frozen}} \left( 2 \langle pi|v|qi \rangle - \langle pi|v|iq \rangle \right).
$$

PySCF (Sun 2018) returns these objects from its CASCI machinery: get_h1eff() gives the dressed one-electron matrix together with $E_0$, and get_h2eff() gives the active-space two-electron integrals (in the Mulliken notation). Pass them to Block2 (Zhai 2023) in memory (pyblock2.driver.DMRGDriver, get_qc_mpo), and write one FCIDUMP file per system: it is the input of the Dice SHCI code (Holmes 2016; Sharma 2017) and the record of the Hamiltonian every solver shares.

Choosing the π orbitals. On a curved cage the σ–π separation is not exact. Curvature pyramidalizes each carbon — about 11.6° for C₆₀, more for C₂₀ (POAV analysis; Haddon 1987, 1993) — so the π orbitals acquire s-character and tilt, and no global axis defines them. Because the σ framework is frozen, this choice affects the energy, not just the representation. The construction has two steps:

1. Assign each atom a radial axis from the geometry: the vector from the cage center, or its POAV refinement from the three neighbors. This is a direction, not an orbital.
2. Get the orbitals by projection or localization. AVAS (Sayfutyarova 2017): build one reference orbital per atom, the p function along its radial axis; project the MO space onto them; keep the eigenvectors with large overlap eigenvalues. PySCF's avas module accepts only AO labels, so supply the directed reference orbitals yourself; the projection is a few lines. Alternatively, localize the occupied and virtual valence spaces separately (IBO, Knizia 2013; Pipek & Mezey 1989) and keep the most radial orbital per carbon; the axis enters only in this final selection, not in the localization itself.

Checks: the AVAS overlap spectrum should show a clean gap between kept and discarded orbitals (a smeared spectrum means the σ–π separation is failing; report it); one active orbital per atom; degeneracies consistent with the point group; an RHF energy unchanged by the freeze. Get one recipe working, then test the other and report how the correlation energy moves.

Fix one geometry per cage and commit the XYZ file. DFT geometries (B3LYP-D3/6-31G*) for every isomer up to C₆₀ are in the supplementary information of Liu 2024. A small basis (def2-SVP or cc-pVDZ) is fine. The symmetric I_h C₂₀ cage is Jahn–Teller unstable, not a minimum; use a distorted stationary point — the Liu structures already are distorted, or re-optimize once at a stated level. Every solver and both spin states use the same file.

Tasks

First, validate the pipeline. Match DMRG to SHCI on C₂₀ CAS(20,20) to sub-mEh. SHCI (Dice) gives the near-FCI reference: variational selection plus Epstein–Nesbet PT2, extrapolated to zero threshold. Agreement of two independent solvers on one FCIDUMP validates the pipeline; nothing downstream is trustworthy without it.

Mandatory

1. C₂₀ CAS(20,20): a converged DMRG energy matching SHCI, the π-space correlation energy, and the entanglement map (single-orbital entropies, orbital-pair mutual information). Then an ordering study: energy error versus bond dimension for three orbital orderings, including Fiedler (from the mutual-information graph; the default, built into Block2), latitude (atoms sorted along a cage axis, so each cut is a small ring), and naive (atom order of the XYZ file; the control). If the entanglement follows the cage graph, Fiedler and latitude should nearly coincide; check whether they do.
2. C₃₆ CAS(36,36) with DMRG: the first C₃₆ total energy with the full π space correlated, reported as a discarded-weight extrapolation with an error bar, not a single-bond-dimension number. The discarded weight is the density-matrix weight dropped at each truncation; at small truncation the energy error is linear in it, and extrapolating to zero is standard practice (Olivares-Amaya 2015). Use spin-adapted DMRG and ramp the bond dimension to 2000–4000 multiplets (one HPC cluster node, overnight), fit the largest-D points, and check that the fit is linear. The result is complete even if the bar is a few mEh. Add the S–T gap if time allows; it is a second run of the same cost.

Optional

3. Tree tensor networks (Nakatani & Chan 2013). The question is whether a tree requires a smaller bond dimension than an MPS to reach the same accuracy on a cage, and how the answer changes from C₂₀ to C₃₆. Begin with the Heisenberg model on the cage graph, the setting of Rausch 2021, in ITensorNetworks.jl; these calculations are cheap and cannot fail. Compare an MPS, a tree that follows the structure of the cage (one subtree per symmetry-equivalent cluster of atoms: the twelve pentagons of C₆₀, the D_6h rings of C₃₆), and a tree whose connectivity is optimized by local reconnection (Hikihara 2023). Then repeat the comparison for the ab initio Hamiltonians of CAS(20,20) and CAS(36,36). The main obstacle in the ab initio case is the tree tensor network operator (TTNO), the tree analog of the matrix product operator (MPO): opsum_to_ttn in ITensorNetworks.jl constructs one automatically but has not been tested at this scale, while T3NS (Gunst 2018, 2019) performs ab initio tree calculations directly from an FCIDUMP. The advantage should not be expected to be large: some edge of any tree still bisects the cage, so a tree does not evade the bisection bound; it can only shorten the effective interaction range and place the dominant cuts more favorably. Tree tensor networks for quantum chemistry are well established (Nakatani & Chan 2013), as is their symmetry adaptation (Gunst 2019); what is new is the fullerene application and the size dependence of the advantage. Measure it and report it, whether or not it is significant.

Stretch

4. C₆₀ CAS(60,60) with an MPS or a TTNS: an extrapolated energy with an error bar.

Compute reality check

The challenge is posted before the school: the environment, the Dice build, and the full C₂₀ pipeline can be done in advance on a laptop. The large runs are overnight cluster jobs. Estimated costs, calibrating on C₂₀ before trusting them:

Run	Hardware	Walltime
C20 SHCI, threshold-extrapolated	laptop	under an hour (several threshold points)
C20 DMRG, D up to 1000	laptop	under an hour
C36 DMRG, D 500–1000	laptop	1–3 h
C36 DMRG, D 2000–4000	one cluster node, 32–64 cores	overnight
C60 DMRG, D 1000–2000	one cluster node, 32–64 cores	overnight to a day

All runs fit on one cluster node; multi-node parallelism is unnecessary. Give Block2 generous memory and scratch disk, tens of GB at the largest bond dimensions. Convergence will be slower than for quasi-1D π systems such as polyacenes. That slowdown is data: it is what the tree comparison probes.

Hints & pitfalls

• DMRG practice: SU(2) spin adaptation, Fiedler or genetic orbital ordering, bond-dimension ramp with decreasing noise, discarded-weight extrapolation. The energy is variational and must fall monotonically with bond dimension.
• S–T gap: same orbitals and same ordering for both spin sectors, or the error cancellation is lost. Target each spin sector explicitly; the frontier is near-degenerate and Jahn–Teller active.
• A cheap CASCI on a small subspace is an exact check DMRG must reproduce.
• If you impose point-group symmetry on the tensors, verify the target quantum numbers are conserved: a silent symmetry bug shows up as a too-high energy.
• Report the active-space construction and its effect on the energy; it is a real uncertainty.
• Tooling: build Dice first; compilation is the main friction (pyscf.hci is a fallback for the variational stage). Block2 installs from PyPI and includes the pyblock2 API; the dmrgscf plugin wraps it as a PySCF fcisolver, but prefer explicit integrals, so that SHCI and DMRG solve the same recorded Hamiltonian. Block2 has no TTNS; tree calculations use ITensorNetworks.jl or T3NS. Pin package versions.

Tools

• PySCF — HF, integrals / FCIDUMP, CASSCF, and an SHCI interface.
• Dice — SHCI / semistochastic heat-bath CI (the near-FCI reference).
• Block2 — ab initio DMRG, DMRG-SCF, orbital reordering, NEVPT2.
• ITensorNetworks.jl (Julia) — tree tensor networks; ITensorMPS.jl — DMRG.
• T3NS — ab initio three-legged tree tensor networks; reads FCIDUMP.

Deliverables

Fork and work under the MPS/DMRG track's solutions/<your-team>/ directory (e.g. tracks/mps/solutions/<your-team>/; confirm the track path at the help desk). Open a PR containing:

1. Harness improvement: a reusable routine from a fullerene geometry to (active-space FCIDUMP, DMRG/SHCI energy, single-orbital entropies, point-group orbit data, candidate tree) for any C_n. Add a .knowledge model card ("ab initio π active space of fullerenes") and a benchmark card recording the C₂₀ SHCI energy.
2. Solution scripts + README stating what was computed and why it matters.
3. One reproduction prompt (PR comment) via the challenge-report skill, reproducing the C₂₀ DMRG↔SHCI agreement.

Generated data goes to a gitignored results/.

References

• Grossman 1995 — J. C. Grossman, L. Mitas, K. Raghavachari, Structure and stability of molecular carbon: importance of electron correlation, Phys. Rev. Lett. 75, 3870 (1995). doi:10.1103/PhysRevLett.75.3870 — why C₂₀ is a canonical strong-correlation case (cage/bowl/ring near-degeneracy breaks DFT).
• Lin 2007 — F. Lin, E. S. Sørensen, C. Kallin, A. J. Berlinsky, Strong correlation effects in the fullerene C₂₀, Phys. Rev. B 76, 033414 (2007). doi:10.1103/PhysRevB.76.033414 — prior C₂₀ many-body study, but a one-band Hubbard model, not ab initio.
• Rausch 2021 — R. Rausch, C. Plorin, M. Peschke, The antiferromagnetic S=1/2 Heisenberg model on the C₆₀ fullerene geometry, SciPost Phys. 10, 087 (2021). arXiv:2011.12083 — prior fullerene tensor-network study, but a Heisenberg spin model, not ab initio.
• Gunnarsson 1997 — O. Gunnarsson, Superconductivity in fullerides, Rev. Mod. Phys. 69, 575 (1997). doi:10.1103/RevModPhys.69.575 — correlation in the C₆₀ π shell as the input to fulleride superconductivity.
• Holmes 2016 — A. A. Holmes, N. M. Tubman, C. J. Umrigar, Heat-bath configuration interaction, J. Chem. Theory Comput. 12, 3674 (2016). doi:10.1021/acs.jctc.6b00407.
• Sharma 2017 — S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, C. J. Umrigar, Semistochastic heat-bath configuration interaction (SHCI), J. Chem. Theory Comput. 13, 1595 (2017). doi:10.1021/acs.jctc.6b01028.
• Sun 2018 — Q. Sun et al., PySCF: the Python-based simulations of chemistry framework, WIREs Comput. Mol. Sci. 8, e1340 (2018). doi:10.1002/wcms.1340.
• Zhai 2023 — H. Zhai et al., Block2: a comprehensive open source framework to develop and apply state-of-the-art DMRG algorithms in electronic structure and beyond, J. Chem. Phys. 159, 234801 (2023). doi:10.1063/5.0180424.
• Knizia 2013 — G. Knizia, Intrinsic atomic orbitals, J. Chem. Theory Comput. 9, 4834 (2013). doi:10.1021/ct400687b — active-space (radial-p) construction.
• Pipek & Mezey 1989 — J. Pipek, P. G. Mezey, A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions, J. Chem. Phys. 90, 4916 (1989). doi:10.1063/1.456588.
• Sayfutyarova 2017 — E. R. Sayfutyarova, Q. Sun, G. K.-L. Chan, G. Knizia, Automated construction of molecular active spaces from atomic valence orbitals, J. Chem. Theory Comput. 13, 4063 (2017). doi:10.1021/acs.jctc.7b00128 — AVAS.
• Liu 2024 — B. Liu, J. Jin, M. Liu, Mapping structure-property relationships in fullerene systems: a computational study from C₂₀ to C₆₀, npj Comput. Mater. 10, 227 (2024). doi:10.1038/s41524-024-01410-7 — DFT (B3LYP-D3/6-31G*) geometries of all C₂₀–C₆₀ isomers, XYZ files in the SI.
• Olivares-Amaya 2015 — R. Olivares-Amaya, W. Hu, N. Nakatani, S. Sharma, J. Yang, G. K.-L. Chan, The ab-initio density matrix renormalization group in practice, J. Chem. Phys. 142, 034102 (2015). doi:10.1063/1.4905329 — DMRG practice, including the discarded-weight extrapolation.
• Nakatani & Chan 2013 — N. Nakatani, G. K.-L. Chan, Efficient tree tensor network states for quantum chemistry, J. Chem. Phys. 138, 134113 (2013). arXiv:1302.2298.
• Gunst 2018, 2019 — K. Gunst, F. Verstraete, S. Wouters, Ö. Legeza, D. Van Neck, T3NS: three-legged tree tensor network states, J. Chem. Theory Comput. 14, 2026 (2018), arXiv:1801.09998; …with SU(2) and molecular point group symmetry, JCTC 15, 2996 (2019), arXiv:1901.08926 — prior art on symmetry in trees.
• Hikihara 2023 — T. Hikihara, H. Ueda, K. Okunishi, K. Harada, T. Nishino, Automatic structural optimization of tree tensor networks, Phys. Rev. Research 5, 013031 (2023). doi:10.1103/PhysRevResearch.5.013031 — the entanglement-optimized-tree baseline.
• Haddon 1987 — R. C. Haddon, Rehybridization and π-orbital overlap in nonplanar conjugated organic molecules: π-orbital axis vector (POAV) analysis and 3D-HMO theory, J. Am. Chem. Soc. 109, 1676 (1987). doi:10.1021/ja00240a013 — the POAV method.
• Haddon 1993 — R. C. Haddon, Chemistry of the fullerenes: the manifestation of strain in a class of continuous aromatic molecules, Science 261, 1545 (1993). doi:10.1126/science.261.5128.1545 — POAV rehybridization/strain on fullerenes.

Start with the C20 pipeline described in the issue: PySCF CASCI get_h1eff()/get_h2eff(), Block2 DMRGDriver/get_qc_mpo, and Dice using a shared FCIDUMP. For the optional tree work, begin with the Heisenberg cage model in ITensorNetworks.jl and inspect opsum_to_ttn. Done means sub-mEh C20 validation against SHCI and a discarded-weight-extrapolated C36 DMRG energy with an error bar.