"""Slater-Condon rules for matrix elements between determinants.
Implements the standard Slater-Condon rules for evaluating
<D|H|D'> where H = S_{pq} h_{pq} a+_p a_q + 1/4 S_{pqrs} g_{pqrs} a+_p a+_q a_s a_r
The rules are:
* **Diagonal** (D = D'):
S_i h_{ii} + 1/2 S_{ij} (g_{ijij} - g_{ijji})
* **Single excitation** (D' = a+_a a_i D):
sign x [h_{ia} + S_j (g_{ijaj} - g_{ijja})]
* **Double excitation** (D' = a+_a a_i a+_b a_j D, i<j, a<b):
sign x (g_{ijab} - g_{ijba})
All other matrix elements are zero for a two-body Hamiltonian.
References
----------
* Szabo & Ostlund, *Modern Quantum Chemistry*, Sec.2.3.
* Helgaker, Jorgensen & Olsen, *Molecular Electronic-Structure Theory*, Sec.4.2.
"""
from __future__ import annotations
from typing import Optional
import numpy as np
from ._determinant import Det, excitation_rank
# -- Phase factor ----------------------------------------------------------
def _phase_factor(occ: Det, excitation_path: list[tuple[int, int]]) -> int:
"""Compute the fermionic phase (-1)^P for a sequence of creations/annihilations.
Parameters
----------
occ : Det
Sorted tuple of occupied orbital indices.
excitation_path : list[tuple[int, int]]
List of (hole, particle) pairs, ordered as applied.
Returns
-------
sign : int (+1 or -1)
"""
# Build bitmask and compute via permutation count
occ_set = set(occ)
sign = 1
for hole, particle in excitation_path:
# Count occupied orbitals strictly between hole and particle
lo = min(hole, particle)
hi = max(hole, particle)
for idx in occ_set:
if lo < idx < hi:
sign *= -1
# Update occupation
occ_set.discard(hole)
occ_set.add(particle)
return sign
# -- Diagonal --------------------------------------------------------------
def diagonal_matrix_element(
occ: Det,
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> float:
"""<D|H|D> for a closed-shell (spin-restricted) determinant.
Each orbital in ``occ`` is doubly occupied (a + b).
Formula:
E_D = 2 S_i h_{ii} + S_{ij} (2 g_{ijij} - g_{ijji})
"""
occ_arr = np.asarray(occ, dtype=np.intp)
E = 2.0 * float(np.sum(h1e[occ_arr, occ_arr]))
ii, jj = np.meshgrid(occ_arr, occ_arr, indexing="ij")
E += float(np.sum(2.0 * h2e_phys[ii, jj, ii, jj] - h2e_phys[ii, jj, jj, ii]))
return E
[docs]
def diagonal_matrix_element_unrestricted(
alpha_occ: Det,
beta_occ: Det,
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> float:
"""<D|H|D> for an unrestricted determinant.
E_D = S_{iina} h_{ii} + S_{iinb} h_{ii}
+ 1/2 S_{iina} S_{jina} (g_{ijij} - g_{ijji})
+ 1/2 S_{iinb} S_{jinb} (g_{ijij} - g_{ijji})
+ S_{iina} S_{jinb} g_{ijij}
"""
a_arr = np.asarray(alpha_occ, dtype=np.intp)
b_arr = np.asarray(beta_occ, dtype=np.intp)
E = float(np.sum(h1e[a_arr, a_arr]) + np.sum(h1e[b_arr, b_arr]))
# aa
ii, jj = np.meshgrid(a_arr, a_arr, indexing="ij")
E += 0.5 * float(np.sum(h2e_phys[ii, jj, ii, jj] - h2e_phys[ii, jj, jj, ii]))
# bb
ii, jj = np.meshgrid(b_arr, b_arr, indexing="ij")
E += 0.5 * float(np.sum(h2e_phys[ii, jj, ii, jj] - h2e_phys[ii, jj, jj, ii]))
# ab
E += float(
np.sum(h2e_phys[a_arr[:, None], b_arr[None, :], a_arr[:, None], b_arr[None, :]])
)
return E
def batch_diagonal_elements(
det_list: list[Det],
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> np.ndarray:
"""Vectorized diagonal elements for a list of closed-shell determinants.
Parameters
----------
det_list : list[Det]
List of determinants.
h1e : (norb, norb) ndarray
h2e_phys : (norb, norb, norb, norb) ndarray
Returns
-------
diag : (ndet,) ndarray
"""
diag = np.zeros(len(det_list))
for idx, occ in enumerate(det_list):
occ_arr = np.asarray(occ, dtype=np.intp)
d = 2.0 * np.sum(h1e[occ_arr, occ_arr])
ii, jj = np.meshgrid(occ_arr, occ_arr, indexing="ij")
d += np.sum(2.0 * h2e_phys[ii, jj, ii, jj] - h2e_phys[ii, jj, jj, ii])
diag[idx] = d
return diag
def batch_diagonal_elements_unrestricted(
det_list: list[SpinDet],
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> np.ndarray:
"""Vectorized diagonal elements for a list of unrestricted SpinDet.
Parameters
----------
det_list : list[SpinDet]
List of (alpha_occ, beta_occ) tuples.
h1e : (norb, norb) ndarray
h2e_phys : (norb, norb, norb, norb) ndarray
Returns
-------
diag : (ndet,) ndarray
"""
diag = np.zeros(len(det_list))
for idx, (a_occ, b_occ) in enumerate(det_list):
diag[idx] = diagonal_matrix_element_unrestricted(a_occ, b_occ, h1e, h2e_phys)
return diag
# -- Single excitation -----------------------------------------------------
def single_excitation_matrix_element(
occ: Det,
i: int, # hole
a: int, # particle
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> float:
"""<D^{a}_{i}|H|D> for a spin-restricted (closed-shell) determinant.
D^{a}_{i} has the a electron from orbital i promoted to a.
The phase is (-1)^{#occupied between i and a}.
Formula:
h_{ia} + S_j (2 g_{ijaj} - g_{ijja})
"""
sign = _phase_factor(occ, [(i, a)])
occ_arr = np.asarray(occ, dtype=np.intp)
val = h1e[i, a]
val += float(
np.sum(
2.0 * h2e_phys[i, occ_arr, a, occ_arr] - h2e_phys[i, occ_arr, occ_arr, a]
)
)
return sign * val
# -- Double excitation -----------------------------------------------------
def double_excitation_matrix_element(
occ: Det,
i: int, # hole 1
j: int, # hole 2 (i < j)
a: int, # particle 1
b: int, # particle 2 (a < b)
h2e_phys: np.ndarray,
) -> float:
"""<D^{ab}_{ij}|H|D> for a spin-restricted determinant.
Both electrons are a-electrons promoted from i,j to a,b.
Phase: (-1)^{#occupied between i,a + #occupied between j,b after first excitation}.
Formula:
g_{ijab} - g_{ijba}
"""
# Phase for two successive single excitations
sign = _phase_factor(occ, [(i, a), (j, b)])
return sign * (h2e_phys[i, j, a, b] - h2e_phys[i, j, b, a])
# -- Hamiltonian matrix element (general) ----------------------------------
def hamiltonian_matrix_element(
bra_occ: Det,
ket_occ: Det,
h1e: np.ndarray,
h2e_phys: np.ndarray,
*,
nelec: Optional[int] = None,
) -> float:
"""Compute <D_bra|H|D_ket> for spin-restricted determinants.
Dispatches to diagonal / single / double rules based on excitation rank.
Returns 0 for rank > 2.
"""
rank = excitation_rank(bra_occ, ket_occ)
if rank == 0:
return diagonal_matrix_element(bra_occ, h1e, h2e_phys)
bra_set = set(bra_occ)
ket_set = set(ket_occ)
holes = sorted(ket_set - bra_set) # orbitals in ket but not bra
particles = sorted(bra_set - ket_set) # orbitals in bra but not ket
if rank == 1:
i, a = holes[0], particles[0]
return single_excitation_matrix_element(ket_occ, i, a, h1e, h2e_phys)
if rank == 2:
i, j = holes
a, b = particles
return double_excitation_matrix_element(ket_occ, i, j, a, b, h2e_phys)
return 0.0
# -- Hamiltonian-vector product (s-vector) ---------------------------------
def hamiltonian_dot(
coeffs: np.ndarray,
det_list: list[Det],
h1e: np.ndarray,
h2e_phys: np.ndarray,
*,
diagonal_precomputed: Optional[np.ndarray] = None,
) -> np.ndarray:
"""Compute H.c (s-vector) for a CI wavefunction.
.. deprecated:: 0.15
This Python O(ndet^2) implementation is superseded by the C++
``casci_direct_sigma`` kernel (OpenMP-parallel, O(ndet) pair-driven
Slater-Condon). The C++ path is the production default for all CI
s-vector builds; this function is retained for validation parity only.
"""
import warnings
warnings.warn(
"hamiltonian_dot is deprecated; use the C++ casci_direct_sigma kernel",
DeprecationWarning,
stacklevel=2,
)
"""Parameters
----------
coeffs : (ndet,) ndarray
CI coefficient vector.
det_list : list[Det]
List of determinants in the CI expansion.
h1e : (norb, norb) ndarray
h2e_phys : (norb, norb, norb, norb) ndarray
diagonal_precomputed : (ndet,) ndarray, optional
Precomputed diagonal elements for speed.
Returns
-------
sigma : (ndet,) ndarray
H.c.
"""
import warnings
warnings.warn(
"hamiltonian_dot is deprecated; use the C++ casci_direct_sigma kernel",
DeprecationWarning,
stacklevel=2,
)
ndet = len(det_list)
sigma = np.zeros(ndet)
# Precompute diagonals if not provided
if diagonal_precomputed is None:
diag = batch_diagonal_elements(det_list, h1e, h2e_phys)
else:
diag = diagonal_precomputed
sigma += diag * coeffs
# Off-diagonal: loop over all pairs of determinants
# For efficiency, we precompute the singles/doubles connections
norb = h1e.shape[0]
for I in range(ndet):
occ_I = det_list[I]
occ_set = set(occ_I)
vir_set = set(range(norb)) - occ_set
vir_list = sorted(vir_set)
# Singles from I
for i in occ_set:
for a in vir_list:
# Build the single-excited determinant
new_occ = tuple(sorted((occ_set - {i}) | {a}))
try:
J = det_list.index(new_occ)
val = single_excitation_matrix_element(occ_I, i, a, h1e, h2e_phys)
sigma[I] += val * coeffs[J]
except ValueError:
pass # not in the CI space
# Doubles from I
occ_list = sorted(occ_set)
for idx_i in range(len(occ_list)):
i = occ_list[idx_i]
for idx_j in range(idx_i + 1, len(occ_list)):
j = occ_list[idx_j]
for idx_a in range(len(vir_list)):
a = vir_list[idx_a]
for idx_b in range(idx_a + 1, len(vir_list)):
b = vir_list[idx_b]
new_occ = tuple(sorted((occ_set - {i, j}) | {a, b}))
try:
J = det_list.index(new_occ)
val = double_excitation_matrix_element(
occ_I, i, j, a, b, h2e_phys
)
sigma[I] += val * coeffs[J]
except ValueError:
pass
return sigma
# -- Build full Hamiltonian matrix -----------------------------------------
def build_hamiltonian_matrix(
det_list: list[Det],
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> np.ndarray:
"""Build the full CI Hamiltonian matrix H_{IJ} = <D_I|H|D_J>.
Parameters
----------
det_list : list[Det]
List of determinants.
h1e : (norb, norb) ndarray
h2e_phys : (norb, norb, norb, norb) ndarray
Returns
-------
H : (ndet, ndet) ndarray
Symmetric Hamiltonian matrix.
"""
ndet = len(det_list)
H = np.zeros((ndet, ndet))
norb = h1e.shape[0]
# Precompute diagonal
diag = batch_diagonal_elements(det_list, h1e, h2e_phys)
H[np.diag_indices(ndet)] = diag
# Off-diagonal via connected determinants
# Build index map for O(1) lookups instead of O(ndet) list.index()
det_to_idx = {d: i for i, d in enumerate(det_list)}
for I in range(ndet):
occ_I = det_list[I]
occ_set = set(occ_I)
vir_set = set(range(norb)) - occ_set
vir_list = sorted(vir_set)
# Singles
for i in occ_set:
for a in vir_list:
new_occ = tuple(sorted((occ_set - {i}) | {a}))
J = det_to_idx.get(new_occ)
if J is not None and J > I:
val = single_excitation_matrix_element(occ_I, i, a, h1e, h2e_phys)
H[I, J] = val
H[J, I] = val
# Doubles
occ_list = sorted(occ_set)
for idx_i in range(len(occ_list)):
i = occ_list[idx_i]
for idx_j in range(idx_i + 1, len(occ_list)):
j = occ_list[idx_j]
for idx_a in range(len(vir_list)):
a = vir_list[idx_a]
for idx_b in range(idx_a + 1, len(vir_list)):
b = vir_list[idx_b]
new_occ = tuple(sorted((occ_set - {i, j}) | {a, b}))
J = det_to_idx.get(new_occ)
if J is not None and J > I:
val = double_excitation_matrix_element(
occ_I, i, j, a, b, h2e_phys
)
H[I, J] = val
H[J, I] = val
return H
# -- Unrestricted Slater-Condon rules -------------------------------------
from ._determinant import SpinDet
def _unrestricted_single_excitation(
occ_ref: Det,
occ_other: Det,
i: int,
a: int,
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> float:
"""<...a+_a a_i...|H|D> for an excitation in one spin sector.
The other spin sector's occupation is ``occ_other`` and is
unchanged by the excitation. The phase is computed from
``occ_ref`` (the sector being excited).
Formula:
sign x [h_{ia} + S_{kinocc_ref} (g_{ikak} - g_{ikka})
+ S_{kinocc_other} g_{ikak}]
"""
sign = _phase_factor(occ_ref, [(i, a)])
ref_arr = np.asarray(occ_ref, dtype=np.intp)
other_arr = np.asarray(occ_other, dtype=np.intp)
val = h1e[i, a]
val += float(
np.sum(h2e_phys[i, ref_arr, a, ref_arr] - h2e_phys[i, ref_arr, ref_arr, a])
)
val += float(np.sum(h2e_phys[i, other_arr, a, other_arr]))
return sign * val
def hamiltonian_matrix_element_unrestricted(
alpha_bra: Det,
beta_bra: Det,
alpha_ket: Det,
beta_ket: Det,
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> float:
"""<a_bra b_bra|H|a_ket b_ket> for unrestricted determinants.
Dispatches based on excitation ranks in each spin sector.
Total rank = rank_a + rank_b <= 2 for non-zero matrix elements.
"""
rank_a = excitation_rank(alpha_bra, alpha_ket)
rank_b = excitation_rank(beta_bra, beta_ket)
total = rank_a + rank_b
if total == 0:
return diagonal_matrix_element_unrestricted(alpha_bra, beta_bra, h1e, h2e_phys)
if total > 2:
return 0.0
holes_a = sorted(set(alpha_ket) - set(alpha_bra))
parts_a = sorted(set(alpha_bra) - set(alpha_ket))
holes_b = sorted(set(beta_ket) - set(beta_bra))
parts_b = sorted(set(beta_bra) - set(beta_ket))
# Case 1: Single excitation in a only
if rank_a == 1 and rank_b == 0:
return _unrestricted_single_excitation(
alpha_ket, beta_ket, holes_a[0], parts_a[0], h1e, h2e_phys
)
# Case 2: Single excitation in b only
if rank_a == 0 and rank_b == 1:
return _unrestricted_single_excitation(
beta_ket, alpha_ket, holes_b[0], parts_b[0], h1e, h2e_phys
)
# Case 3: Double excitation, both in a
if rank_a == 2 and rank_b == 0:
return double_excitation_matrix_element(
alpha_ket, holes_a[0], holes_a[1], parts_a[0], parts_a[1], h2e_phys
)
# Case 4: Double excitation, both in b
if rank_a == 0 and rank_b == 2:
return double_excitation_matrix_element(
beta_ket, holes_b[0], holes_b[1], parts_b[0], parts_b[1], h2e_phys
)
# Case 5: One excitation in each spin sector
if rank_a == 1 and rank_b == 1:
sign = _phase_factor(alpha_ket, [(holes_a[0], parts_a[0])])
sign *= _phase_factor(beta_ket, [(holes_b[0], parts_b[0])])
return sign * h2e_phys[holes_a[0], holes_b[0], parts_a[0], parts_b[0]]
return 0.0
[docs]
def build_hamiltonian_matrix_unrestricted(
det_list: list[SpinDet],
h1e: np.ndarray,
h2e_phys: np.ndarray,
) -> np.ndarray:
"""Build the full CI Hamiltonian matrix in the unrestricted
determinant basis.
Uses a connection map to avoid O(ndet^2) all-pairs iteration --
only computes matrix elements for determinants connected by
<= 2 excitations.
Parameters
----------
det_list : list[SpinDet]
List of (alpha_occ, beta_occ) tuples.
h1e : (norb, norb) ndarray
h2e_phys : (norb, norb, norb, norb) ndarray
Returns
-------
H : (ndet, ndet) ndarray
Symmetric Hamiltonian matrix.
"""
ndet = len(det_list)
H = np.zeros((ndet, ndet))
norb = h1e.shape[0]
# Precompute diagonal
diag = batch_diagonal_elements_unrestricted(det_list, h1e, h2e_phys)
H[np.diag_indices(ndet)] = diag
# Build index map for O(1) lookup
det_to_idx = {d: i for i, d in enumerate(det_list)}
# Off-diagonal: only connected pairs via explicit excitation generation
for I in range(ndet):
aI, bI = det_list[I]
occ_a = set(aI)
occ_b = set(bI)
vir_a = sorted(set(range(norb)) - occ_a)
vir_b = sorted(set(range(norb)) - occ_b)
# a singles
for i in occ_a:
for a in vir_a:
new_a = tuple(sorted((occ_a - {i}) | {a}))
J = det_to_idx.get((new_a, bI))
if J is not None and J > I:
val = hamiltonian_matrix_element_unrestricted(
new_a, bI, aI, bI, h1e, h2e_phys
)
if abs(val) > 1e-16:
H[I, J] = val
H[J, I] = val
# b singles
for i in occ_b:
for a in vir_b:
new_b = tuple(sorted((occ_b - {i}) | {a}))
J = det_to_idx.get((aI, new_b))
if J is not None and J > I:
val = hamiltonian_matrix_element_unrestricted(
aI, new_b, aI, bI, h1e, h2e_phys
)
if abs(val) > 1e-16:
H[I, J] = val
H[J, I] = val
# aa doubles
a_list = sorted(occ_a)
for idx_i in range(len(a_list)):
for idx_j in range(idx_i + 1, len(a_list)):
i, j = a_list[idx_i], a_list[idx_j]
for idx_p in range(len(vir_a)):
for idx_q in range(idx_p + 1, len(vir_a)):
a, b = vir_a[idx_p], vir_a[idx_q]
new_a = tuple(sorted((occ_a - {i, j}) | {a, b}))
J = det_to_idx.get((new_a, bI))
if J is not None and J > I:
val = hamiltonian_matrix_element_unrestricted(
new_a, bI, aI, bI, h1e, h2e_phys
)
if abs(val) > 1e-16:
H[I, J] = val
H[J, I] = val
# bb doubles
b_list = sorted(occ_b)
for idx_i in range(len(b_list)):
for idx_j in range(idx_i + 1, len(b_list)):
i, j = b_list[idx_i], b_list[idx_j]
for idx_p in range(len(vir_b)):
for idx_q in range(idx_p + 1, len(vir_b)):
a, b = vir_b[idx_p], vir_b[idx_q]
new_b = tuple(sorted((occ_b - {i, j}) | {a, b}))
J = det_to_idx.get((aI, new_b))
if J is not None and J > I:
val = hamiltonian_matrix_element_unrestricted(
aI, new_b, aI, bI, h1e, h2e_phys
)
if abs(val) > 1e-16:
H[I, J] = val
H[J, I] = val
# ab doubles
for i in occ_a:
for j in occ_b:
for a in vir_a:
for b in vir_b:
new_a = tuple(sorted((occ_a - {i}) | {a}))
new_b = tuple(sorted((occ_b - {j}) | {b}))
J = det_to_idx.get((new_a, new_b))
if J is not None and J > I:
val = hamiltonian_matrix_element_unrestricted(
new_a, new_b, aI, bI, h1e, h2e_phys
)
if abs(val) > 1e-16:
H[I, J] = val
H[J, I] = val
return H