File size: 5,304 Bytes
2a87ede | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 | """Restricted Hartree-Fock self-consistent field loop (Roothaan-Hall).
Standard textbook algorithm (e.g. Szabo & Ostlund, "Modern Quantum
Chemistry", ch. 3): orthogonalize the AO basis via S^(-1/2), build the
Fock matrix F = H_core + 2J - K from the current density, diagonalize
in the orthogonal basis, form a new density, repeat to convergence.
This part is genuinely simple compared to the integral evaluation and
doesn't need vectorizing -- a closed-shell molecule's SCF loop is a few
dozen matrix multiplies on an N x N matrix where N is a few tens at
most for STO-3G, nowhere near where PennyLane's implementation loses
its time (which is entirely in building H_core/repulsion tensor, done
once in assembly.py, not in this loop).
BUG FOUND (Si2, minimal 4-electron/4-orbital active space, R=2.184 A):
plain (undamped) density substitution never converged for this system
-- 100/100 iterations, still oscillating -- because two pairs of
orbitals near the active-space boundary are numerically degenerate
(HOMO-1/HOMO and LUMO/LUMO+1 each split by <1e-9 Ha), so each iteration
flips which member of a near-tied pair gets occupied, and the density
never settles. Confirmed this is a real oscillation, not just slow
convergence: three separate machines/runs of the undamped loop each
hit the iteration cap at a DIFFERENT total energy (-571.63, -570.69,
-571.02 Ha), all physically meaningless artifacts of whatever step the
loop happened to be on. Fixed with standard linear density damping
(P_next = alpha*P_new + (1-alpha)*P_old) -- textbook remedy for exactly
this oscillation failure mode (Szabo & Ostlund ch. 3.4.9). Verified
the damped loop converges to the SAME energy (-570.874032094871 Ha,
agreeing to 10 significant figures) across alpha in {0.1, 0.2, 0.3,
0.5, 0.7} -- the fix's correctness doesn't hinge on the specific alpha
chosen, only that damping is applied at all.
"""
import dataclasses
import numpy as np
@dataclasses.dataclass
class HFResult:
converged: bool
n_iterations: int
electronic_energy: float
nuclear_repulsion_energy: float
total_energy: float
orbital_energies: np.ndarray
orbital_coefficients: np.ndarray # C, shape (n_basis, n_basis)
density_matrix: np.ndarray
def nuclear_repulsion_energy(nuclear_charges: list[float], nuclear_positions: np.ndarray) -> float:
energy = 0.0
n = len(nuclear_charges)
for i in range(n):
for j in range(i + 1, n):
r = np.linalg.norm(nuclear_positions[i] - nuclear_positions[j])
energy += nuclear_charges[i] * nuclear_charges[j] / r
return energy
def _orthogonalizer(S: np.ndarray) -> np.ndarray:
w, v = np.linalg.eigh(S)
return v @ np.diag(1.0 / np.sqrt(w)) @ v.T
def _density_from_coefficients(C: np.ndarray, n_occupied_pairs: int) -> np.ndarray:
C_occ = C[:, :n_occupied_pairs]
return C_occ @ C_occ.T
def run_scf(
S: np.ndarray,
H_core: np.ndarray,
repulsion: np.ndarray,
n_electrons: int,
nuclear_charges: list[float],
nuclear_positions: np.ndarray,
max_iterations: int = 200,
convergence_tol: float = 1e-10,
damping: float = 0.5,
) -> HFResult:
if n_electrons % 2 != 0:
raise ValueError("Only closed-shell (even electron count) systems are supported.")
n_occupied_pairs = n_electrons // 2
X = _orthogonalizer(S)
orbital_energies, C_ortho = np.linalg.eigh(X.T @ H_core @ X)
C = X @ C_ortho
P = _density_from_coefficients(C, n_occupied_pairs)
converged = False
for iteration in range(1, max_iterations + 1):
J = np.einsum("pqrs,rs->pq", repulsion, P)
K = np.einsum("prqs,rs->pq", repulsion, P)
F = H_core + 2.0 * J - K
orbital_energies, C_ortho = np.linalg.eigh(X.T @ F @ X)
C = X @ C_ortho
P_new = _density_from_coefficients(C, n_occupied_pairs)
if np.linalg.norm(P_new - P) < convergence_tol:
P = P_new
converged = True
break
# Linear density damping (Szabo & Ostlund ch. 3.4.9): plain
# substitution (damping=1.0) oscillates indefinitely for systems
# with near-degenerate orbitals at the occupied/virtual boundary
# (see module docstring, Si2) instead of converging -- mixing in
# a fraction of the old density every step damps that
# oscillation. Verified to reproduce the same converged energy
# to 10 significant figures across damping in [0.1, 0.7] for the
# system that motivated this, and to leave already-convergent
# systems' converged energies unchanged (just costs a handful of
# extra iterations, well inside max_iterations).
P = damping * P_new + (1.0 - damping) * P
electronic_energy = float(np.sum(P * (H_core + F)))
e_nuc = nuclear_repulsion_energy(nuclear_charges, nuclear_positions)
return HFResult(
converged=converged,
n_iterations=iteration,
electronic_energy=electronic_energy,
nuclear_repulsion_energy=e_nuc,
total_energy=electronic_energy + e_nuc,
orbital_energies=orbital_energies,
orbital_coefficients=C,
density_matrix=P,
)
|