Add basic orbital optimization for a fixed 1-/2-RDM (restricted case)

[SemyonAndreyev]

Summary

Adds optimize_orbitals() (hamiltonians/orbital_optimization.py) for Issue #711: given a fixed 1-/2-particle RDM (e.g. from a truncated active-space CI/CASCI calculation, or recovered from a quantum device via sample-based diagonalization), finds a restricted (closed-shell) orbital rotation that minimizes the energy those RDM coefficients would give under a different one-particle basis, holding the RDM fixed.

Reuses the existing restricted-orbital parametrization already in hartree_fock.py (rhf_params_to_matrix — an antihermitian generator restricted to occupied/virtual blocks, exponentiated into a unitary) and general_basis_change (already used by HartreeFockFunctional.__init__) to rotate the Hamiltonian integrals into the trial basis at each step. scipy.optimize.minimize drives the search; the gradient is scipy's numerical one, not analytic — deliberately "basic," per the issue's own wording.

Correctness property, verified rather than assumed: for a full (untruncated) active space, the identity rotation is provably the global minimum of this objective — any rotation keeps the trial state inside the same complete N-electron Fock space that full CI already minimizes over exactly, so no rotation can score below the FCI energy. Tested this directly against real FCI diagonalization data (not asserted): an earlier version of this PR held the Hamiltonian fixed and rotated the RDM instead, expecting the energy landscape to be flat/invariant for a full space — that assumption failed against real H2 FCI data, which is what surfaced the correct formulation (rotate the Hamiltonian, keep the RDM fixed) and the correct test criterion (global minimum, not invariance).

Test plan

• test_optimize_orbitals_at_identity_reproduces_full_ci_energy — energy at the identity rotation exactly matches the FCI energy (sanity check that the objective is wired correctly)
• test_full_space_fci_rdm_is_never_beaten_by_any_rotation — for a full-space RDM, explicit nonzero rotations never score below the FCI floor, and the optimizer started away from the identity converges back down to it (not below)
• test_optimize_orbitals_improves_a_truncated_active_space — the realistic use case: a CASCI-style truncated-active-space RDM is improved (or at worst unchanged) by orbital rotation between active and excluded space
• test_optimize_orbitals_rejects_degenerate_orbital_split — raises ValueError when there's no occupied/virtual split to rotate over
• Full hamiltonians/ test suite: 161/161 passing, no regressions
• check/format-incremental, check/pylint-changed-files, check/mypy (pinned tool versions) all clean

Scope note

Test fixtures use H2 (sto-3g and 6-31g, both already shipped in the package's test data — no new dependency). Did not generate a LiH fixture: doing so needs PySCF, which requires a C/C++ toolchain not available in my dev environment; scoped the tests to H2 at multiple bond lengths/basis sets instead, which already exercises both the full-space and truncated-active-space code paths.

Fixes #711.

[gemini-code-assist]

Code Review

This pull request introduces a basic restricted (closed-shell) orbital optimization feature (optimize_orbitals) to find a rotation of one-particle orbitals that minimizes the total energy for a fixed 1- and 2-particle reduced density matrix (RDM). It includes comprehensive unit tests verifying physical correctness and active-space truncation improvements. The review feedback suggests adding robust input validation checks at the entry point of optimize_orbitals, specifically verifying the shapes of the integral and RDM arrays, ensuring n_electrons is even, and validating the size of the initial_guess parameter vector.

[gemini-code-assist]

The function optimize_orbitals lacks validation for the shapes of the input arrays (one_body_integrals, two_body_integrals, one_rdm, two_rdm) and the parity of n_electrons. Since this is a restricted closed-shell orbital optimization, n_electrons must be even, and the dimensions of the RDMs must match the number of orbitals. Adding defensive validation checks at the beginning of the function prevents runtime errors (such as shape mismatches or index errors) and provides clear, actionable error messages to the user.

if one_body_integrals.ndim != 2 or one_body_integrals.shape[0] != one_body_integrals.shape[1]:
        raise ValueError("one_body_integrals must be a square 2D array.")
    n_orbitals = one_body_integrals.shape[0]

    if two_body_integrals.shape != (n_orbitals, n_orbitals, n_orbitals, n_orbitals):
        raise ValueError(
            f"two_body_integrals must have shape {(n_orbitals,) * 4}, "
            f"but got {two_body_integrals.shape}."
        )
    if one_rdm.shape != (2 * n_orbitals, 2 * n_orbitals):
        raise ValueError(
            f"one_rdm must have shape {(2 * n_orbitals, 2 * n_orbitals)}, "
            f"but got {one_rdm.shape}."
        )
    if two_rdm.shape != (2 * n_orbitals,) * 4:
        raise ValueError(
            f"two_rdm must have shape {(2 * n_orbitals,) * 4}, "
            f"but got {two_rdm.shape}."
        )

    if n_electrons % 2 != 0:
        raise ValueError(
            f"Restricted closed-shell orbital optimization requires an even number of electrons, "
            f"but got n_electrons={n_electrons}."
        )

    nocc = n_electrons // 2
    nvirt = n_orbitals - nocc
    if nocc <= 0 or nvirt <= 0:

References

1. Validate and cast input parameters early in the class constructor or entry points, rather than performing complex type validation and casting in downstream utility functions.

[gemini-code-assist]

If initial_guess is provided, its size is not validated against the expected number of parameters (nocc * nvirt). If a user passes an initial_guess with an incorrect size, it will lead to an IndexError or silent bugs inside rhf_params_to_matrix. Adding a validation check ensures that the size of initial_guess matches the expected parameter count.

Suggested change

	if initial_guess is None:
	init_params = np.zeros(nocc * nvirt)
	else:
	init_params = np.asarray(initial_guess).flatten()
	if initial_guess is None:
	init_params = np.zeros(nocc * nvirt)
	else:
	init_params = np.asarray(initial_guess).flatten()
	if init_params.size != nocc * nvirt:
	raise ValueError(
	f"initial_guess has size {init_params.size}, but expected {nocc * nvirt} "
	f"(nocc={nocc}, nvirt={nvirt})."
	)

References

1. Validate and cast input parameters early in the class constructor or entry points, rather than performing complex type validation and casting in downstream utility functions.

[mhucka]

@SemyonAndreyev Thank you for this contribution. The comments by Gemini Code Assist actually seem reasonable. When you get a chance, can you address them?

[SemyonAndreyev]

Addressed both high-priority comments from Gemini Code Assist:

• optimize_orbitals now validates one_body_integrals/two_body_integrals/one_rdm/two_rdm shapes against each other and against n_electrons (must be even — this is the restricted/closed-shell case) before anything reaches energy(), where a mismatch would otherwise surface as an opaque broadcast/index error deep inside general_basis_change/rhf_params_to_matrix.
• initial_guess size is now validated against the expected nocc * nvirt parameter count for the same reason.

4 new regression tests cover the new checks (test_optimize_orbitals_rejects_mismatched_shapes, test_optimize_orbitals_rejects_odd_electron_count, test_optimize_orbitals_rejects_mismatched_initial_guess). All 7 tests in orbital_optimization_test.py pass, black/isort/pylint clean. Pushed as 6847abb.