本文介绍了基于AttilaSzabo和NeilS.Ostlund著作的Hartree-Fock方法的Python实现流程，包括计算初始数据、构造Fock矩阵、自洽场迭代和能量求解，旨在帮助初学者理解量子化学计算基础。
一、前言
本文全部内容都基于 Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory （作者：Attila Szabo，Neil S. Ostlund）中关于 Hartree-Fock 的部分 （146 ~ 164页）。文中不涉及太多理论，仅仅介绍了一个自洽场计算的流程希望帮助读者理解，另外 Python 对于新手来说是非常简单易懂的语言，刚好可以用来做演示。事实上早在多年前网络上有关于这个章节的讨论，笔者在初学量子化学的时候受到了相应的启发。然而现在网络上已经找不到当时的文章，另外其他某些相关的文章居然开始收费……，以此为背景，笔者重新整理并写下这篇短文希望可以为初学量子化学的新人提供一些帮助，如果存在不严谨或错误的地方，欢迎大家批评指正。
二、流程简介
首先根据初始设定的条件计算初始数据包括核哈密顿、重叠积分等。
然后猜测初始电子密度矩阵构造 Fock 矩阵进行自洽场求解。
不断循环迭代直至达到收敛条件得出体系总能量。
三、程序实现
1、构建一个分子，就是定义这个分子的基本性质，如：原子核坐标 , 原子数量
, 电子数量
, 基函数
.(本文假定了一个
分子，那么就包含两个原子，两个电子，基函数由初始选择的基组决定）同时引入相关的 Python 库
import numpy as np2、初始基函数对不同的算符做积分可以得到初始用于计算的参数，在此不详细介绍积分过程，直接使用积分结果。
- 核哈密顿
，主要是单电子算符对应的哈密顿，包括 1、核的动能，2、每个电子与所有原子核之间的势能：
# Szabo. pp161(p3.233) H_core = np.array([[-1.1204, -0.9584], [-0.9584, -1.1204]])
- 轨道重叠积分
，当两个原子接近时，他们的电子云或轨道会相互重叠，如果两个轨道完全不重叠，S 矩阵的非对角元为0，对角元为1。在本文中：
# Szabo. pp160(p3.229) S_overlap = np.array([[1.0, 0.6593], [0.6593, 1.0]])
3、 对角化重叠矩阵 ，并且得到转换矩阵
, 这一步需要引入 numpy 库中的相关函数。自动完成对角化过程并返回相关的矩阵。
# Szabo. pp.143 (p3.167)
l, U = np.linalg.eigh(S_overlap)
s_1 = np.zeros((2, 2))
for i in range(2):
s_1[i,i] = 1.0 / l[i] ** 0.5
X_transfer = np.dot(np.dot(U, s_1), np.linalg.inv(U)) # 求转换矩阵 X
X_adjoint = np.matrix.getH(X_transfer) # 共轭转置4、随便猜一个密度矩阵 ，懒得猜全都设成 0 也可以，反正后面会进行迭代。
P_density = np.zeros((2, 2)) # 初始化密度矩阵以上都是固定过程，一旦生成后面就不需要再动了。接下来是自洽场部分，程序上就是不断的进行循环直到收敛。
5、根据密度矩阵计算双电子积分 ，其中对矩阵元的积分过程用数值代替了，感兴趣的读者可以自行探索。
def double_ele_integral(u,v,p,q): # u,v,p,q的取值范围为 0 或 1,算两电子积分
if u < v: # u一定要比v大或相等
u,v = v,u
if p < q: # p一定要大于等于q
p,q = q,p # 先 u, v 再 p, q 加一乘是怕有 0 值干扰
if (u+1)*(v+1) < (p+1)*(q+1):
u,v,p,q = p,q,u,v # Szabo. pp162(p3.235)
if (u,v,p,q) == (0,0,0,0) or (u,v,p,q) == (1,1,1,1):
return 0.7746
elif (u,v,p,q) == (1,1,0,0):
return 0.5697
elif (u,v,p,q) == (1,0,0,0) or (u,v,p,q) == (1,1,1,0):
return 0.4441
elif (u,v,p,q) == (1,0,1,0):
return 0.2970
else:
raise
def G_two_ele_integral(density): #Szabo., pp141(3.154)
G_matrix = np.zeros(density.shape)
dim = density.shape[0]
for u in range(dim):
for v in range(dim):
temp = 0.
for p in range(dim):
for q in range(dim):
doubleJ = double_ele_integral(u,v,p,q)
doubleK = 0.5 * double_ele_integral(u,q,p,v)
temp += density[p,q] * (doubleJ - doubleK)
G_matrix[u,v] = temp
return G_matrix
6、得到最终的 Fock 矩阵
Fock = H_core + G_matrix7、将变换矩阵 与
作用在
上得到
Fock_prime = X_adjoint.dot(Fock.dot(X_transfer))8、求 的本征矢量
结合变换矩阵
求
_, C_prime = np.linalg.eigh(Fock_prime)9、通过 和
得到新的密度矩阵
def calc_P_rhf(C_new):
density = np.zeros((2,2))
for i in range (2):
for j in range (2):
for x in range (1):
density[i,j] += 2*C_new[i,x]*C_new[j,x]
return density10、最后通过密度矩阵、核哈密顿、Fock矩阵求得体系的总能量
#Szabo. pp150(3.184)
def calc_E0_rhf(density, H_core, Fock):
H_F = H_core + Fock
E0 = 0.
for u in range(2):
for v in range(2):
E0 += density[v,u] * H_F[u,v]
E0 *= 0.5
return E0步骤 6 到步骤 10 为一个循环，每一次循环可以通过初始的密度矩阵得到新的密度矩阵，并得到当前状态下的电子能量。如果没有满足收敛条件，则可以根据新的密度矩阵作为输入再进行一次迭代得到下一个密度矩阵。
11、最后，体系的总能量还要加上原子核之间的排斥。
E1 = 1 / 1.4 # 两 H 原子间距离为 1.4 （原子单位）
Etot = E0 + E112、经过三次迭代之后，结果如下：单位为 a.u. 与真实情况一致。
**** iteration: 1 ****
Etot:, 0.7143
**** iteration: 2 ****
Etot:, -1.1168
**** iteration: 3 ****
Etot:, -1.1168四、总结
最终程序如下：
import numpy as np
# Szabo. pp161(p3.233)
H_core = np.array([[-1.1204, -0.9584],
[-0.9584, -1.1204]])
# Szabo. pp160(p3.229)
S_overlap = np.array([[1.0, 0.6593],
[0.6593, 1.0]])
# Szabo. pp.143 (p3.167)
l, U = np.linalg.eigh(S_overlap)
s_1 = np.zeros((2, 2))
for i in range(2):
s_1[i,i] = 1.0 / l[i] ** 0.5
X_transfer = np.dot(np.dot(U, s_1), np.linalg.inv(U)) # 求转换矩阵 X
X_adjoint = np.matrix.getH(X_transfer) # 共轭转置
P_density = np.zeros((2, 2)) # 初始化密度矩阵
def double_ele_integral(u,v,p,q): # u,v,p,q的取值范围为 0 或 1,算两电子积分
if u < v: # u一定要比v大或相等
u,v = v,u
if p < q: # p一定要大于等于q
p,q = q,p # 先 u, v 再 p, q 加一乘是怕有 0 值干扰
if (u+1)*(v+1) < (p+1)*(q+1):
u,v,p,q = p,q,u,v # Szabo. pp162(p3.235)
if (u,v,p,q) == (0,0,0,0) or (u,v,p,q) == (1,1,1,1):
return 0.7746
elif (u,v,p,q) == (1,1,0,0):
return 0.5697
elif (u,v,p,q) == (1,0,0,0) or (u,v,p,q) == (1,1,1,0):
return 0.4441
elif (u,v,p,q) == (1,0,1,0):
return 0.2970
else:
raise
#Szabo., pp141(3.154)
def G_two_ele_integral(density):
G_matrix = np.zeros(density.shape)
dim = density.shape[0]
for u in range(dim):
for v in range(dim):
temp = 0.
for p in range(dim):
for q in range(dim):
doubleJ = double_ele_integral(u,v,p,q)
doubleK = 0.5 * double_ele_integral(u,q,p,v)
temp += density[p,q] * (doubleJ - doubleK)
G_matrix[u,v] = temp
return G_matrix
#Szabo. pp150(3.184)
def calc_E0_rhf(density, H_core, Fock):
H_F = H_core + Fock
E0 = 0.
for u in range(2):
for v in range(2):
E0 += density[v,u] * H_F[u,v]
E0 *= 0.5
return E0
def calc_P_rhf(C_new):
density = np.zeros((2,2))
for i in range (2):
for j in range (2):
for x in range (1):
density[i,j] += 2*C_new[i,x]*C_new[j,x]
return density
for i in range(3):
G_matrix = G_two_ele_integral(P_density)
Fock = H_core + G_matrix # obtain Fock' by rotating Fock matrix # Szabo. pp.145 (3.177)
Fock_prime = X_adjoint.dot(Fock.dot(X_transfer)) # Solve Eigenvalue problem # Szabo. pp.145 (3.178)
_, C_prime = np.linalg.eigh(Fock_prime) # obtain C_new by rotating C_new_prime
C_new = X_transfer.dot(C_prime) # Calculate Density Matrix
P_density_new = calc_P_rhf(C_new)
E0 = calc_E0_rhf(P_density, H_core, Fock)
P_density = P_density_new
E1 = 1 / 1.4 # 两 H 原子间距离为 1.4 （原子单位）
Etot = E0 + E1
print("**** iteration: {} ****".format(i+1))
print(f"Etot:, {Etot:.4f}")
程序背后的理论背景，无需在文中详述。相关公式都给出了在书中的位置。至于书籍读者可以在网络上自行寻找，搭配程序，效果更佳。
代码方面有很大提升空间，很多运算可以用矩阵乘法进行简化，但会导致程序更加抽象难以理解，读者可以根据自身习惯进行修改。
为了简化程序使之更易于新手理解，文中省略了很多积分过程。事实上积分部分十分困难，感兴趣的读者可以挑战。
扫一扫
2190
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
