Hartree-Fock 自洽场计算流程

本文档以 Attila Szabo and Neil S. Ostlund 所编写的 《Modern Quantum Chemistry - Introduction to Advanced Electronic Structure Theory》 一书中的 $He{H}^{+}$ 分子为例展示电子积分及 Hartree-Fock 自洽场的计算流程。或者说是自己基于 Python 对该书附录 B - Two-Electron Self-Consistent-Field Program 部分的实现。
主要内容包括：

• 构建STO-NG高斯基组（展示1G、2G、3G基组参数）

• 实现关键积分计算函数：

  • 重叠积分（公式A.9）
  • 动能积分（公式A.11）
  • 核吸引积分（公式A.33，使用Boys函数）

通过Gaussian类处理基函数操作，包括归一化和乘积计算矩阵构建方法生成重叠矩阵、动能矩阵和核吸引矩阵。程序可扩展用于更复杂分子的电子结构计算。
更多信息请关注原书。

##### HeH+ 分子案例 #####import numpy as np
from itertools import product
from scipy.special import erf# STO-NG 基组参数 (H原子)
sto_1g = {"d":[1.], # 收缩系数"alphas":[0.270950]} # 原始指数sto_2g = {"d":[0.678914, 0.430129],"alphas":[0.151623, 0.851819]}sto_3g = {"d":[0.444635, 0.535328, 0.154329],"alphas":[0.109818, 0.405771, 2.22766]}def get_scaled_basis(zeta: float, basis: dict):"""根据给定的zeta缩放基组指数 (公式 3.224)"""return {"d": basis["d"],"alphas": [alpha * zeta**2 for alpha in basis["alphas"]]}class Gaussian:"""高斯函数类"""def __init__(self, alpha: float = 1., center: np.ndarray = np.array([0., 0., 0.])):self.alpha = alpha  # 指数self.center = center # 中心坐标 [x, y, z]self.coeff = 1.  # 系数def normalize(self):"""将系数重设为归一化系数（归一化操作）"""new = Gaussian(self.alpha, self.center)new.coeff = (2*self.alpha/np.pi)**0.75return newdef __str__(self):return f"{self.coeff}*G[-{self.alpha}(r-{self.center})**2]"def gaussians_product(gaus0:Gaussian, gaus1:Gaussian):"""两个高斯函数的乘积为一个新的高斯函数 (公式 A.1 - A.5)"""alpha_sum = gaus0.alpha + gaus1.alphadist_square = np.sum((gaus0.center - gaus1.center)**2) # 两个高斯函数中心之间距离的平方new_center = (gaus0.alpha * gaus0.center + gaus1.alpha * gaus1.center) / alpha_sum # 新高斯函数的中心, 公式 A.4K = np.exp(-gaus0.alpha*gaus1.alpha / alpha_sum * dist_square) # 公式 A.3new_gaus = Gaussian(alpha_sum, new_center)new_gaus.coeff *= K * gaus0.normalize().coeff * gaus1.normalize().coeff # 乘上原来两个高斯函数的归一化系数return new_gausdef gaussian_integral(gaus):"""计算高斯函数的全空间积分"""return (np.pi/gaus.alpha)**1.5 * gaus.coeff # 根据gamma函数算 r 积分，再乘上角度积分 4*pidef overlap_integral(gaus0: Gaussian, gaus1: Gaussian):"""计算两个高斯函数的重叠积分 (公式 A.9)"""gaus0 = gaus0.normalize()gaus1 = gaus1.normalize()alpha_sum = gaus0.alpha + gaus1.alphadist_square = np.sum((gaus0.center - gaus1.center)**2) # 两个高斯函数中心之间距离的平方tem = (np.pi/alpha_sum)**1.5 * np.exp(-gaus0.alpha*gaus1.alpha/alpha_sum * dist_square) # 公式 A.9return gaus0.coeff * gaus1.coeff * tem # 要加上归一化系数def get_overlap_matrix(bases:list[dict], centers:list):"""计算分子的重叠积分矩阵"""n_basis_funs = len(bases)overlap_list = []for (basis0, center0), (basis1, center1) in product(zip(bases, centers), repeat=2):tem = 0.for i, j in product(range(len(basis0["alphas"])), repeat=2):gaus0 = Gaussian(basis0["alphas"][i], center0)gaus1 = Gaussian(basis1["alphas"][j], center1)tem += (basis0["d"][i] * basis1["d"][j]) * overlap_integral(gaus0, gaus1)overlap_list.append(tem)return np.array(overlap_list).reshape((n_basis_funs, n_basis_funs))def kinetic_integral(gaus0: Gaussian, gaus1: Gaussian):"""计算动能积分 (公式 A.11)"""alpha_prod = gaus0.alpha * gaus1.alphaalpha_sum = gaus0.alpha + gaus1.alphadist_square = np.sum((gaus0.center - gaus1.center)**2)S = overlap_integral(gaus0, gaus1) # 计算重叠积分部分, 其内部包含了归一化系数return (alpha_prod/alpha_sum) * (3 - 2*alpha_prod/alpha_sum * dist_square) * Sdef get_kinetic_matrix(bases:list[dict], centers:list):"""计算分子的动能积分矩阵"""n_basis_funs = len(bases)kinetic_list = []for (basis0, center0), (basis1, center1) in product(zip(bases, centers), repeat=2):tem = 0.for i, j in product(range(len(basis0["alphas"])), repeat=2):gaus0 = Gaussian(basis0["alphas"][i], center0)gaus1 = Gaussian(basis1["alphas"][j], center1)tem += (basis0["d"][i] * basis1["d"][j]) * kinetic_integral(gaus0, gaus1)kinetic_list.append(tem)return np.array(kinetic_list).reshape((n_basis_funs, n_basis_funs))def F_0(t):"""构造 Boys 函数, 公式 A.32"""if t < 1e-10:return 1.0return 0.5 * np.sqrt(np.pi/t) * erf(np.sqrt(t))def attraction_integral(gaus0: Gaussian, gaus1: Gaussian, center_c: np.ndarray, Zc: float):"""计算核吸引积分 (公式 A.33)"""# 计算乘积高斯函数alpha_prod = gaus0.alpha * gaus1.alphaalpha_sum = gaus0.alpha + gaus1.alphacenter_p = gaussians_product(gaus0, gaus1).centerdist_pc_square = np.sum((center_p - center_c)**2)dist_ab_square = np.sum((gaus0.center - gaus1.center)**2)tem = -2*np.pi / alpha_sum * Zc * np.exp(-alpha_prod/alpha_sum * dist_ab_square) * F_0(alpha_sum*dist_pc_square)return tem * gaus0.normalize().coeff * gaus1.normalize().coeffdef get_attraction_matrix(bases:list[dict], centers:list, center2:np.ndarray, Zc):"""计算分子的核吸引积分矩阵"""n_basis_funs = len(bases)attraction_list = []for (basis0, center0), (basis1, center1) in product(zip(bases, centers), repeat=2):tem = 0.for i, j in product(range(len(basis0["alphas"])), repeat=2):gaus0 = Gaussian(basis0["alphas"][i], center0)gaus1 = Gaussian(basis1["alphas"][j], center1)tem += (basis0["d"][i] * basis1["d"][j]) * attraction_integral(gaus0, gaus1, center2, Zc)attraction_list.append(tem)return np.array(attraction_list).reshape((n_basis_funs, n_basis_funs))def two_electron_integral(gaus0:Gaussian, gaus1:Gaussian, gaus2:Gaussian, gaus3:Gaussian):"""计算双电子积分 公式 A.41"""alpha_sum_p = gaus0.alpha + gaus1.alphaalpha_sum_q = gaus2.alpha + gaus3.alphaalpha_prod_p = gaus0.alpha * gaus1.alphaalpha_prod_q = gaus2.alpha * gaus3.alphacenter_p = gaussians_product(gaus0, gaus1).centercenter_q = gaussians_product(gaus2, gaus3).centerdist_ab_squrare =  np.sum((gaus0.center - gaus1.center)**2)dist_cd_square = np.sum((gaus2.center - gaus3.center)**2)dist_pq_square = np.sum((center_p - center_q)**2)term_0 = 2 * np.pi**2.5 / ((alpha_sum_p * alpha_sum_q * np.sqrt(alpha_sum_p + alpha_sum_q)))term_1 = np.exp(-alpha_prod_p/alpha_sum_p * dist_ab_squrare - alpha_prod_q/alpha_sum_q * dist_cd_square)term_2 = F_0(alpha_sum_p * alpha_sum_q / (alpha_sum_p + alpha_sum_q) * dist_pq_square)term_3 = gaus0.normalize().coeff * gaus1.normalize().coeff * gaus2.normalize().coeff * gaus3.normalize().coeff return term_0 * term_1 * term_2 *  term_3def get_two_electron_tensor(bases:list[dict], centers:list):"""计算双电子积分张量"""n_basis_funs = len(bases)tensor_list = []for (basis0, center0), (basis1, center1),(basis2, center2), (basis3, center3) in product(zip(bases,centers), repeat=4):tmp = 0.for m, n, p, q in product(range(len(basis0["alphas"])), repeat=4):gaus0 = Gaussian(basis0["alphas"][m], center0)gaus1 = Gaussian(basis1["alphas"][n], center1)gaus2 = Gaussian(basis2["alphas"][p], center2)gaus3 = Gaussian(basis3["alphas"][q], center3)tmp += (basis0["d"][m] * basis1["d"][n] * basis2["d"][p] * basis3["d"][q]) * two_electron_integral(gaus0, gaus1, gaus2, gaus3)tensor_list.append(tmp)return np.array(tensor_list).reshape([n_basis_funs]*4)def get_fock_two_electron(P:np.ndarray, tensor:np.ndarray):"""计算 Fock 矩阵的双电子部分"""J = np.einsum('yo,uvoy->uv', P, tensor) k = np.einsum('yo,uyov->uv', P, tensor) return J - 0.5 * kdef get_density_matrix(C:np.ndarray, n_occ:int):"""从HF系数矩阵计算闭壳层密度矩阵"""# 提取占据轨道的系数矩阵 (前n_occ列)C_occ = C[:, :n_occ]P = 2 * np.einsum('ua,va->uv', C_occ, C_occ)return Pdef canonical_orthogonalization(S:np.ndarray):"""正则正交化：计算变换矩阵 X 及其逆矩阵 X^{-1}"""assert np.allclose(S, S.T), "重叠矩阵 S 必须对称！"eigenvalues, U = np.linalg.eigh(S)  # 对角化 Seigenvalues = eigenvalues[::-1] # 本征值从大到小排U = U[:, ::-1]# 检查本征值是否为正（避免数值误差）if np.any(eigenvalues <= 1e-10):print("警告：重叠矩阵 S 接近奇异，可能需截断小本征值！")eigenvalues = np.maximum(eigenvalues, 1e-10)  # 截断小本征值# 构造 Λ^{-1/2} 和 Λ^{1/2}lambda_inv_sqrt = np.diag(1.0 / np.sqrt(eigenvalues))X = np.einsum('ij,jk->ik', U, lambda_inv_sqrt)return Xdef get_hf_energy(P:np.ndarray, h_core:np.ndarray, F:np.ndarray):"""计算 HF 能量 公式 3.184"""return 0.5 * np.einsum('vu,uv->', P, h_core + F)# 定义两个原子的位置 (沿z轴)
bond_length = 1.4632
center_he = np.array([0., 0., 0])
center_h = np.array([0., 0., bond_length])
centers = [center_he, center_h]zeta_he = 2.0925
zeta_h = 1.24
basis_he = get_scaled_basis(zeta_he, sto_3g)
basis_h = get_scaled_basis(zeta_h, sto_3g)
bases = [basis_he, basis_h]### 重叠积分矩阵
overlap_matrix = get_overlap_matrix(bases, centers)
print("重叠积分矩阵 公式 3.251")
print(overlap_matrix)

重叠积分矩阵 公式 3.251
[[1.00000143 0.45077041][0.45077041 1.00000143]]

#### 动能积分矩阵
kinetic_matrix = get_kinetic_matrix(bases, centers)
print("\n动能积分矩阵 公式 3.252")
print(kinetic_matrix)

动能积分矩阵 公式 3.252
[[2.16431256 0.16701287][0.16701287 0.76003294]]

### 核吸引积分矩阵
center2 = np.array([0., 0., 0])
center3 = np.array([0., 0., bond_length])
attraction_matrix_2 = get_attraction_matrix(bases, centers, center2, 2)
attraction_matrix_3 = get_attraction_matrix(bases, centers, center3, 1)
print("\n核吸引积分矩阵 公式 3.253")
print(attraction_matrix_2)
print("\n核吸引积分矩阵 公式 3.254")
print(attraction_matrix_3)

核吸引积分矩阵 公式 3.253
[[-4.1398272  -1.10291257][-1.10291257 -1.26524587]]核吸引积分矩阵 公式 3.254
[[-0.67723008 -0.41130547][-0.41130547 -1.22661547]]

### 芯矩阵
h_core = kinetic_matrix + attraction_matrix_2 + attraction_matrix_3
print("\n芯矩阵 公式 3.255")
print(h_core)

芯矩阵 公式 3.255
[[-2.65274472 -1.34720517][-1.34720517 -1.7318284 ]]

### 双电子积分张量
two_electron_tensor = get_two_electron_tensor(bases, centers)
print("\n双电子积分张量 P172 双电子积分")
print(two_electron_tensor[0,0,0,0])
print(two_electron_tensor[1,1,0,0])
print(two_electron_tensor[1,0,0,0])
print(two_electron_tensor[1,1,1,0])
print(two_electron_tensor[1,0,1,0])
print(two_electron_tensor[1,1,1,1])

双电子积分张量 P172 双电子积分
1.307151607555482
0.6057033663335994
0.4372793252654168
0.31179457036897396
0.1772671219506615
0.7746083600328786

### 变换矩阵 X
X = canonical_orthogonalization(overlap_matrix)
X[1,1] = -X[1,1] # 添加上这两行就和书上完全一样了
X[0,1] = -X[0,1]
print("\n正交化变换矩阵 X 公式 3.262")
print(X)

正交化变换矩阵 X 公式 3.262
[[ 0.58706399  0.95412983][ 0.58706399 -0.95412983]]

F = X.T @ h_core @ X
print("\n Fock 矩阵 公式 3.266")
print(F)

Fock 矩阵 公式 3.266
[[-2.43973011 -0.51583772][-0.51583772 -1.53866291]]

epsilon, C_prime = np.linalg.eigh(F)
print("\nC' 矩阵 公式 3.268")
print(C_prime)

C' 矩阵 公式 3.268
[[-0.91044566  0.4136287 ][-0.4136287  -0.91044566]]

print("\nepsilon 公式 3.269")
print(epsilon)

epsilon 公式 3.269
[-2.67408268 -1.30431033]

C = X @ C_primeprint("\n系数矩阵 C 公式 3.270")
print(C)

系数矩阵 C 公式 3.270
[[-0.92914535 -0.62585685][-0.13983438  1.11150988]]

P = get_density_matrix(C, 1)
print("\n密度矩阵 公式 3.271")
print(P)

密度矩阵 公式 3.271
[[1.72662215 0.25985293][0.25985293 0.03910731]]

print("\nFock 双电子部分 G 公式 3.273")
G = get_fock_two_electron(P, two_electron_tensor)
print(G)

Fock 双电子部分 G 公式 3.273
[[1.26232798 0.37400298][0.37400298 0.98895134]]

F = h_core + G
print("\n新的 Fock 矩阵 公式 3.273")
print(F)

新的 Fock 矩阵 公式 3.273
[[-1.39041674 -0.97320219][-0.97320219 -0.74287706]]

energy = get_hf_energy(P, h_core, F)
print(f"\nHF 能量为 {energy}")

HF 能量为 -4.141860256587706

### 整体迭代
print("================")
print("\nHF 迭代")
max_iter = 10
tol = 1e-6
n_occ = 1 # HeH+ 有两个电子，占据一个轨道
energy_old = 1.F = h_core + 0.  # G 初猜值为 0 矩阵
for iteration in range(max_iter):F_prime = X.T @ F @ Xepsilon, C_prime = np.linalg.eigh(F_prime)C = X @ C_primeP = get_density_matrix(C, n_occ)G = get_fock_two_electron(P, two_electron_tensor)F = h_core + Genergy = get_hf_energy(P, h_core, F)delta_e = abs(energy - energy_old)print(f"Iter {iteration+1}: Energy = {energy:.8f}, ΔE = {delta_e:.2e}")if delta_e < tol:print("\nSCF 收敛了!")break  energy_old = energy

================HF 迭代
Iter 1: Energy = -4.14186026, ΔE = 5.14e+00
Iter 2: Energy = -4.22648847, ΔE = 8.46e-02
Iter 3: Energy = -4.22751948, ΔE = 1.03e-03
Iter 4: Energy = -4.22752582, ΔE = 6.34e-06
Iter 5: Energy = -4.22752586, ΔE = 3.59e-08SCF 收敛了!