fix(orbit): use spatial orbital optimization to preserve spin symmetry - #223
Conversation
The previous implementation diagonalized the full spin-orbital RDM, which mixes alpha and beta spin-orbitals and breaks spin symmetry in restricted calculations. This caused energy to increase after orbital optimization. Changes: - NaturalOrbitCalculator.calculate_rdm() now returns spatial orbital RDM by averaging alpha-alpha and beta-beta blocks - Integrals are transformed at the spatial orbital level first, then expanded to spin-orbital basis with identical transformations for alpha and beta spins - Store block-diagonal spin-orbital unitary for tracking
cb9721b to
4741de9
Compare
There was a problem hiding this comment.
Pull request overview
Updates orbital optimization to operate on spatial (spin-restricted) orbitals so that restricted calculations preserve spin symmetry, avoiding the previous α/β mixing from diagonalizing the full spin-orbital 1-RDM.
Changes:
- Compute a spatial 1-RDM by combining the αα and ββ spin blocks, then diagonalize it to obtain a spatial-orbital rotation.
- Transform spatial one- and two-electron integrals first, then expand to a spin-orbital Hamiltonian using the same rotation for α and β blocks.
- Store a block-diagonal spin-orbital unitary (
U_spin) for tracking.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
| qmp/algorithms/orbit.py | Switches natural orbital construction and integral transformation to spatial-orbital form; stores block-diagonal spin unitary. |
| qmp/algorithms/haar_with_orbit.py | Aligns iterative orbital-optimization step with spatial-orbital (spin-restricted) transformation and storage. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| Calculate the spatial 1-RDM (One-Particle Reduced Density Matrix) for restricted calculations. | ||
|
|
||
| For restricted (closed-shell) systems, we compute the spatial orbital RDM by summing | ||
| contributions from both alpha and beta spin-orbitals, ensuring spin symmetry is preserved. | ||
|
|
There was a problem hiding this comment.
The updated calculate_rdm doc/comments say the spatial RDM is formed by “summing” / D_aa + D_bb, but the implementation actually averages the blocks: (rdm_aa + rdm_bb) / 2. Please make the docstring and inline comments consistent with the chosen convention (sum vs average), and clarify what the resulting eigenvalues represent (per-spin occupations vs spatial occupations).
| # 3. Load original integrals | ||
| # 3. Load original integrals (these are spatial orbital integrals from FCIDUMP) | ||
| (norb, nelec, nspin), e0, h1, h2 = _read_fcidump_tensors(self.src_fcidump) | ||
|
|
There was a problem hiding this comment.
This code now assumes a spin-restricted setting (same spatial-orbital rotation applied to both α/β, and RDM built by combining αα/ββ blocks). However, nspin (FCIDUMP MS2) / electron parity isn’t validated before proceeding, so running this on open-shell or otherwise non-restricted inputs will silently apply a potentially incorrect transformation. Consider explicitly checking for the restricted condition (e.g., MS2==0 and even electron count) and either raising a clear error or falling back to a spin-orbital/unrestricted path.
| # Validate that the FCIDUMP corresponds to a closed-shell, spin-restricted reference. | |
| # For a restricted calculation we require MS2 == 0 and an even number of electrons. | |
| if nspin != 0 or (nelec % 2) != 0: | |
| raise ValueError( | |
| "Orbit optimization currently assumes a spin-restricted (closed-shell) reference, " | |
| f"but FCIDUMP header reports nelec={nelec}, MS2={nspin}. " | |
| "Support for open-shell / unrestricted inputs is not yet implemented." | |
| ) |
| # 4. Load and Transform SPATIAL Integrals | ||
| (norb, nelec, nspin), e0, h1, h2 = _read_fcidump_tensors(self.src_fcidump) | ||
|
|
||
| device = context.device | ||
| U_spatial = U_spatial.to(device=device, dtype=torch.complex128) | ||
| h1 = h1.to(device=device, dtype=torch.complex128) | ||
| h2 = h2.to(device=device, dtype=torch.complex128) |
There was a problem hiding this comment.
Orbital optimization here is now implemented as spin-restricted (shared U_spatial for α/β and spatial-RDM-based rotation), but the FCIDUMP nspin (MS2) / electron parity isn’t checked before applying it. To avoid silently producing incorrect results on open-shell/non-restricted cases, add an explicit validation (or a fallback to the previous spin-orbital/unrestricted transformation) when MS2 != 0 / when α and β occupations differ.
The previous implementation diagonalized the full spin-orbital RDM, which mixes alpha and beta spin-orbitals and breaks spin symmetry in restricted calculations. This caused energy to increase after orbital optimization.
Changes: