編寫一個量化計算程序,求解簡單限制性HartreeFock方程(python編寫)

上次說完HF方程需要自洽場迭代求解，這次我們來寫一個程序實現自洽場過程。下面是上次解讀的傳送門，建議先看理論再來看這次如何實現自洽場迭代程序

淺談量子化學（2）Hartree Fock

解Hartree-Fock-Roothaan方程:

SCF的過程為：

1.明確一個分子（原子的坐標集合{RA},原子數目{ZA},電子的數量{N},和一個基組）

2.計算所有需求的分子積分

3.重疊積分S對角化，並獲得一個轉換的矩陣X（有兩種方式：對稱正交化和標準正交化）

4.獲得一個初始的密度矩陣P

5.通過密度矩陣P和兩電子積分，帶入方程算出fock矩陣的兩電子積分項G

6.把前面得到的Fock矩陣的兩電子積分項與單電子的哈密頓算法積分相加得到Fock矩陣

7.計算轉換後的Fock矩陣，使用公式

8.對角化F'獲得新的係數矩陣C'和能量本徵值

9.將新的係數矩陣與老的係數矩陣通過X聯繫起來

10.把C帶入方程

獲得其密度矩陣，其表徵電子在分子各個原子中或原子之間的分布情況。

11.判斷是否收斂，可用能量或密度矩陣來判斷---新的密度矩陣（能量）是否與先前的密度矩陣（能量）達到了我們設定的誤差值之內

12.如果程序收斂了，就可以用得到的解去算其他的量子性質了

編寫代碼部分

從上面我們可以得到，第4步往前，都是我們在進行自洽場過程之前需要有的結果，明確你要計算的分子，求出它的，及初始密度矩陣P

為了尋求簡單，我們只算H2分子。為什麼算H2分子呢，因為以上要先求出的值呢，書上都已經給我們求好了。

import mathimport numpy as np

def calc_Hcore_H2():    dim = 2    Hcore = np.zeros( (dim,dim) )        Hcore[0,0] = Hcore[1,1] = -1.1204    Hcore[0,1] = Hcore[1,0] = -0.9584    return Hcore
def calc_S_H2():    dim = 2    S = np.zeros( (dim, dim) )    S[0, 0] = S[1,1] = 1.    S[0, 1] = S[1,0] = 0.6593    return S
將S對角化得X：
def symmetric_orthogonalization(S):            l, U = np.linalg.eigh(S)    l_rt = np.zeros( S.shape )    for i in range(len(l)):        l_rt[i,i] = 1./(math.sqrt(l[i]))    X =  U.dot(l_rt.dot(np.conjugate(U)))     return X
def eri_table(u,v,p,q):        if u < v:        u,v = v,u        if p < q:        p,q = q,p        if (u+1)*(v+1) < (p+1)*(q+1):        u,v,p,q = p,q,u,v            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:                print(u,v,p,q)        raise
初始的密度矩陣呢，就寫成零矩陣，反正後期迭代也會產生新的密度矩陣。
def guess_D():    dim = 2    D = np.zeros( (dim,dim) )    return D
前四步我們已經做完了，接下來做第五步計算G，公式見上面
def calc_G_H2(D):
    G = np.zeros(D.shape)    dim = D.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 = eri_table(u,v,p,q)                    K = 0.5 * eri_table(u,q,p,v)                    temp += D[p,q] * (doubleJ - K)            G[u,v] = temp    return G
第六步，相加得到Fock矩陣。這一步由於都是變量，所以需要定義在函數當中迭代（待會一起寫）
第七步為計算轉換後的Fock矩陣，並且第八步對角化F'獲得新的係數矩陣C'和能量本徵值，由於是迭代過程也得寫在迭代的函數當中。
第九步將新矩陣與原矩陣聯繫，第十步計算新的密度矩陣，判斷密度矩陣與原矩陣是否收斂
首先要寫一個判斷是否收斂的函數用rmsd均方根誤差來判斷吧。
def calc_rmsd(m1,m2):    d = m1 - m2    ret =  np.sum(d**2)/d.size    return ret
得到新的密度矩陣D，計算均方誤差，如果小於設定值就收斂，將這個值賦給D，如果沒有小於就進行下一次循環（也得寫在迭代函數中，待會一併給出）。
我們再寫一個計算能量的函數
def calc_E0_rhf(D,Hcore,F):        H_F = Hcore + F    row, col = H_F.shape    E0 = 0.    for u in range(row):        for v in range(col):            E0 += D[v,u] * H_F[u,v]    E0 *= 0.5;    return E0
def calc_E1():                                                        E1 = 1/1.4    return E1
最後整合一下以上需要迭代的量寫一個函數：
def rhf(nelec):    n_occ_orbitals = 1 # int(nelec/2)    max_iteration = 10    convergence_threshold = 1.e-4  #最小收斂值    dim = 2 # The number of Basis functions
    S = calc_S_H2()    X = symmetric_orthogonalization(S)    X_adj = np.matrix.getH(X)  # np.matrix.getH() returns the adjoint matrix 返回伴隨矩陣
    assert np.allclose( X_adj.dot(S.dot(X)), np.identity(dim) ), \            "X does not satisfy (X_adjoint * S * X)"
    Hcore = calc_Hcore_H2()    D = guess_D()
    for i in range(max_iteration):        print("**** iteration: {} ****".format(i+1))        G = calc_G_H2(D)        F = Hcore + G
        F_prime = X_adj.dot(F.dot(X))        e, C_new_prime = np.linalg.eigh(F_prime)  #返回厄密矩陣或對稱矩陣的特徵值和特徵向量。
        C_new = X.dot(C_new_prime)        D_new = calc_D_rhf(C_new, n_occ_orbitals)                E0 = calc_E0_rhf(D, Hcore, F)        rmsd = calc_rmsd(D_new, D)        D = D_new        E1 = calc_E1()        Etot = E0 + E1        print("Etot:  ", Etot)        print("RMSD:", rmsd)        if rmsd < convergence_threshold:            break    return E0
我們只需要運行這個函數就可以得到H2分子的能量了。
結果為：
我們使用高斯計算一下H2，做比較
可以說結果相當的接近了！
 
 
參考網站：量子化學計算コードを自作してみた記録　その１https://qiita.com/swakamoto/items/5694fdb099cb96b050e8#%E9%87%8D%E3%81%AA%E3%82%8A%E7%A9%8D%E5%88%86s%E3%81%AE%E5%AF%BE%E8%A7%92%E5%8C%96
參考書籍：Szabo A , Ostlund N S . Modern quantum chemistry : introduction to advanced electronic structure theory[J]. 1989.
完
如有概念錯誤，還請指出，我將甚為感謝，如需轉載，請在公眾號聊天框說明。
免責聲明：部分圖片及資料來源於網絡，轉載的目的在於傳遞更多信息及分享，如涉及侵權，請聯繫我及時修改或刪除。