import numpy as np
[docs]
def ghf(S, Hcore, Eri, nalpha, nbeta, energy_nuc, I_thld=500, E_thld=1e-10, D_thld=1e-8):
"""Generalized Hartree-Fock (GHF) algorithm for non-collinear spin systems.
Parameters
----------
S: np.ndarray
Overlap matrix.
Hcore: np.ndarray
Core Hamiltonian matrix.
Eri: np.ndarray
Electron repulsion integrals (in chemists' notation, aosym='s8').
nalpha: int
Number of spin alpha electrons.
nbeta: int
Number of spin beta electrons.
energy_nuc: float
Nuclear repulsion energy.
I_thld: int
Maximum number of iterations.
E_thld: float
Threshold for energy convergence.
D_thld: float
Threshold for density matrix convergence.
Returns
-------
results : dict
Hartree Fock calculation results.
"""
nbf = S.shape[0]
n2 = 2 * nbf
nelec = nalpha + nbeta
# Build orthogonalization matrix for GHF (block-diagonal)
eigval, eigvec = np.linalg.eigh(S)
X = eigvec @ np.diag(eigval**(-0.5)) @ eigvec.T
X_ghf = np.kron(np.eye(2), X)
# Extended Hückel Theory (EHT) initial guess for GHF
# Wolfsberg-Helmholtz approximation: F_ij = 0.5 * K * S_ij * (H_ii + H_jj)
Hdiag = np.diag(Hcore)
K = 1.75 # Wolfsberg-Helmholtz constant
F_eht = 0.5 * K * S * (Hdiag[:, None] + Hdiag[None, :])
np.fill_diagonal(F_eht, Hdiag) # Diagonal elements are H_ii
# Build GHF Fock matrix with spin mixing for true GHF
# Introduce off-diagonal spin blocks to break UHF symmetry
F_eht_ghf = np.zeros((n2, n2))
F_eht_ghf[:nbf, :nbf] = F_eht
F_eht_ghf[nbf:, nbf:] = F_eht
# Add small spin-orbit mixing term to break symmetry
# This is crucial for getting true GHF solutions (not UHF)
mixing_factor = 1e-3
F_eht_ghf[:nbf, nbf:] = mixing_factor * F_eht
F_eht_ghf[nbf:, :nbf] = mixing_factor * F_eht
F_ortho = X_ghf.T @ F_eht_ghf @ X_ghf
eps, C_ortho = np.linalg.eigh(F_ortho)
C = X_ghf @ C_ortho
# Initial density matrix (occupy lowest nelec orbitals)
C_occ = C[:, :nelec]
D = C_occ @ C_occ.T
D_old = D.copy()
# Extract spin blocks
Da = D[:nbf, :nbf]
Db = D[nbf:, nbf:]
Dab = D[:nbf, nbf:]
Dba = D[nbf:, :nbf]
# Initial Fock matrices
D_total = Da + Db
J = np.einsum('pqrs,rs->pq', Eri, D_total)
Kaa = np.einsum('prqs,rs->pq', Eri, Da)
Kbb = np.einsum('prqs,rs->pq', Eri, Db)
Kab = np.einsum('prqs,rs->pq', Eri, Dab)
Kba = np.einsum('prqs,rs->pq', Eri, Dba)
Fa = Hcore + J - Kaa
Fb = Hcore + J - Kbb
Fab = -Kab
Fba = -Kba
# Initial energy
E_elec = 0.5 * np.einsum('pq,pq->', Da + Db, Hcore) + \
0.5 * (np.einsum('pq,pq->', Fa, Da) + np.einsum('pq,pq->', Fb, Db) +
np.einsum('pq,pq->', Fab, Dab) + np.einsum('pq,pq->', Fba, Dba))
E_old = E_elec + energy_nuc
# SCF Iterations
for Iter in range(1, I_thld + 1):
# Recompute Fock matrices
D_total = Da + Db
J = np.einsum('pqrs,rs->pq', Eri, D_total)
Kaa = np.einsum('prqs,rs->pq', Eri, Da)
Kbb = np.einsum('prqs,rs->pq', Eri, Db)
Kab = np.einsum('prqs,rs->pq', Eri, Dab)
Kba = np.einsum('prqs,rs->pq', Eri, Dba)
Fa = Hcore + J - Kaa
Fb = Hcore + J - Kbb
Fab = -Kab
Fba = -Kba
# Build GHF Fock matrix (2nbf x 2nbf)
F_ghf = np.zeros((n2, n2))
F_ghf[:nbf, :nbf] = Fa
F_ghf[nbf:, nbf:] = Fb
F_ghf[:nbf, nbf:] = Fab
F_ghf[nbf:, :nbf] = Fba
# Diagonalize in orthogonal basis
F_ortho = X_ghf.T @ F_ghf @ X_ghf
eps, C_ortho = np.linalg.eigh(F_ortho)
C = X_ghf @ C_ortho
# New density matrix (occupy lowest nelec orbitals)
C_occ = C[:, :nelec]
D_new = C_occ @ C_occ.T
D_old = D.copy()
D = D_new
# Extract new spin blocks
Da_new = D[:nbf, :nbf]
Db_new = D[nbf:, nbf:]
Dab_new = D[:nbf, nbf:]
Dba_new = D[nbf:, :nbf]
# Recompute Fock with new density for energy
D_total_new = Da_new + Db_new
J_new = np.einsum('pqrs,rs->pq', Eri, D_total_new)
Kaa_new = np.einsum('prqs,rs->pq', Eri, Da_new)
Kbb_new = np.einsum('prqs,rs->pq', Eri, Db_new)
Kab_new = np.einsum('prqs,rs->pq', Eri, Dab_new)
Kba_new = np.einsum('prqs,rs->pq', Eri, Dba_new)
Fa_new = Hcore + J_new - Kaa_new
Fb_new = Hcore + J_new - Kbb_new
Fab_new = -Kab_new
Fba_new = -Kba_new
# Electronic energy (GHF expression)
E_elec = 0.5 * np.einsum('pq,pq->', Da_new + Db_new, Hcore) + \
0.5 * (np.einsum('pq,pq->', Fa_new, Da_new) + np.einsum('pq,pq->', Fb_new, Db_new) +
np.einsum('pq,pq->', Fab_new, Dab_new) + np.einsum('pq,pq->', Fba_new, Dba_new))
E_total = E_elec + energy_nuc
# Convergence criteria
dE = abs(E_total - E_old)
dRMS = np.mean((D - D_old) ** 2) ** 0.5
# Update
Da, Db, Dab, Dba = Da_new, Db_new, Dab_new, Dba_new
E_old = E_total
# Check convergence
if (dE < E_thld) and (dRMS < D_thld):
break
if Iter == I_thld:
raise Exception(f"Maximum number of SCF cycles ({I_thld}) exceeded.")
# Prepare results
results = {
"energy_elec": float(E_elec),
"energy_nuc": float(energy_nuc),
"energy_total": float(E_total),
"mo_energy": eps,
"mo_coeff": C,
"rdm1": [Da, Db, Dab, Dba], # Full GHF density matrix blocks
}
return results