深度學習模擬密度泛函計算初試筆記（一）

一、前言

量子力學（Quantum Mechanics）是當今計算化學之中常用的獲取體系性質的方法之一，其中以密度泛函理論（Density Functional Theory, DFT）最為流行。DFT可計算給定體系的各種性質，許多文獻通過各種對比充分肯定了DFT在預測化學體系性質中的有效性，其中又以單點能最為基本。然而，DFT的計算成本較高，對於較大體系可能會花費很久的時間。因此，使用成本較低的方法對DFT計算進行模擬便成為了計算化學的一個重要研究方向。

深度學習（Deep Learning）自從上世紀中葉被提出後，一直因為其強大的潛力而廣受關注。然而，顯而易見地，機器學習對於計算機的計算性能具有很高的要求，因此直到本世紀初，機器學習才隨著計算機計算性能的大幅提高而廣泛使用。機器學習的特點是，可以以較低的成本逼近複雜的非線性函數，且隨著模型的設置，可訓練得到極為複雜的模型。

DFT本質上也是一個函數，輸入原子坐標，輸出單點能。那麼，能否使用深度學習模擬DFT計算，從而以較低的成本得到近似精度的計算結果呢？筆者近期便做了一些嘗試。

二、基本思路

作為最初的嘗試，筆者採取非常簡單的單水分子（H₂O）作為體系，首先隨機獲取一些單水分子體系的坐標，通過DFT計算得到其單點能作為參考數據，隨後構建訓練集與測試集，進行深度學習計算。

三、參考數據獲取

筆者最初的嘗試是隨機獲取(-5, 5)內的浮點數，作為坐標中的一項（單位為埃，即10^-10米）。對於每個數據點，共計獲取9個浮點數，分為作為三個原子的xyz坐標。對於每個數據點，採取B3LYP 6-31G^*進行單點能計算，得到單點能數值作為參考值。

然而，實際操作時，筆者發現，由於隨機生成的體系很可能缺乏合理的化學含義（例如兩個原子距離過近，或者原子排布明顯不合理），從而使得進行單點能計算時，存在 SCF不收斂的問題。因此經過多次嘗試后，筆者選擇了固定氧原子坐標為(0.0, 0.0, 0.0)，兩個氫原子的xyz坐標範圍為(-2, 2)，且第一個氫原子的x坐標大於等於0，第二個氫原子的坐標小於等於0。實驗數據中水分子的H-O鍵鍵長約為1埃，H-O-H鍵角約為105°，因此這個份範圍依然可以涵蓋絕大多數合理情況。同時，為了消除初始對稱性的影響，以及促進SCF收斂，筆者加入了nosymm和SCF=conver=6關鍵字。示例輸入文件如下：

# b3lyp 6-31G* nosymm SCF=conver=6

template

0 1
O 0.0 0.0 0.0
H 1.2098128821231413 0.7498742314194025 1.2486296263802839
H -0.04312995252541918 1.161517161983439 1.4739122046307012

筆者共計獲取了1,500個體系的單點能作為參考數據。對於經過以上設置仍出現SCF不收斂的體系，將其刪掉重新生成，直至參考數據集為1500個數據為止。

隨後，筆者將兩個氫的xyz坐標（共計6個數）作為參考數據的輸入值X，將對應的單點能作為參考數據的輸出值y。參考數據集示意如下（可於文末下載）：

X[0] = [1.2098128821231413, 0.7498742314194025, 1.2486296263802839, -0.04312995252541918, 1.161517161983439, 1.4739122046307012]
y[0] = -76.1141702491

X[1] = [0.16316043179887973, 0.8455602443198043, 1.4291116311573102, -1.6856207179243436, -0.9464153792320902, -1.5885870521611554]
y[1] = -76.0434955957

…………

X[1499] = [0.40836400213078816, 1.2678932382271548, 1.1915913112555758, -1.7410125269177223, -1.0842450909696848, -1.6416433580996737]
y[1499] = -76.0138053880

四、訓練集與測試集的構建

由於深度學習對於輸入值的縮放較為敏感，因此筆者採用sklearn中的StandardScaler對輸入數據X進行處理：

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

隨後採取sklearn中的train_test_split對數據進行打亂分組：

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, random_state=24)

得到訓練集1125個數據點，測試集375個數據點。

五、深度學習訓練

作為初試，筆者採取多層感知機（Multilayer Perceptron, MLP）進行訓練。MLP是一種簡單且常用的深度學習模型，包括輸入層、輸出層和隱藏層，層與層之間以一定的權重結合特定的激活函數進行連接。

筆者採取tanh作為激活函數，隱藏層設置為[25, 25]，正則化參數alpha設置為0.01：

from sklearn.neural_network import MLPRegressor

mlp = MLPRegressor(activation="tanh", hidden_layer_sizes=[25, 25], alpha=0.01, solver='lbfgs', max_iter=10000, random_state=45)
mlp.fit(X_train, y_train)

在筆者的MacBook Pro M1上訓練4.6秒即得到結果，首先看一下模型在訓練集和測試集上的表現：

mlp.score(X_train, y_train)
0.9372755336259925

mlp.score(X_test, y_test)
0.7523257261302045

可以看到模型在訓練集上表現很好，在測試集上表現不如訓練集那麼好（說明可能存在一定的過擬合），但0.75的擬合度也是可以接受的。筆者接下來測試一下該模型對水分子結構優化的性能。

六、利用深度學習對水分子進行結構優化

結構優化是DFT計算中很關鍵的一類任務，傳統的結構優化是通過DFT計算採用梯度下降法進行。這裡筆者採取另一種優化算法——協方差自適應進化策略（Covariance Matrix Adaptation – Evolution Strategy, CMA-ES）進行優化。

首先定義需要優化的能量函數：

def energy(xyz):
    return mlp.predict(scaler.transform([xyz]))[0]

隨後採取CMA-ES算法進行優化：

import cma

init_xyz = [0.2, 0.4, 0.7, -0.4, -0.4, -0.5]

es = cma.CMAEvolutionStrategy(init_xyz, 0.5, {"seed":74})
es.optimize(energy)

可以看到，經過0.3秒即優化完成。接下來採用兩個O-H鍵長與H-O-H鍵角作為指標檢驗一下優化結果的合理性：

import math

r1 = (es.result.xbest[0]**2+es.result.xbest[1]**2+es.result.xbest[2]**2)**0.5
r2 = (es.result.xbest[3]**2+es.result.xbest[4]**2+es.result.xbest[5]**2)**0.5
r3 = ((es.result.xbest[0]-es.result.xbest[3])**2+(es.result.xbest[1]-es.result.xbest[4])**2+(es.result.xbest[2]-es.result.xbest[5])**2)**0.5
a1 = math.acos((r1**2+r2**2-r3**2)/(2*r1*r2))*180/math.pi

其結果為：

r1 = 0.8778885328377475
r2 = 1.0248547218080677
a1 = 78.25493239118879

作為參考，筆者採用與參考數據集相同的DFT級別（即B3LYP/6-31G*）對水分子進行優化，固定O坐標為(0.0, 0.0, 0.0)，得到兩個氫原子坐標分別為(0.944156, 0.185549, -0.112264)和(-0.179397, -0.682453, -0.663734)，其鍵長與鍵角為：

r1 = 0.9687425816144349
r2 = 0.968745685706006
a1 = 103.65158225677882

對比如下：

	r1	r2	a1
深度學習	0.8779	1.0249	78.25
DFT	0.9687	0.9687	103.65
相對偏差	-9.4%	+5.8%	-24.5%

可見O-H鍵長擬合度較好，而H-O-H角度擬合度較差。此模型仍有進一步優化的空間，但考慮到水分子中O-H作用較強，H…H作用較弱，因此鍵長對化學性質的影響比起鍵角更大。考慮到僅僅是這樣一個簡單的模型依然可以給出定性正確的模型，可見採取深度學習擬合DFT計算依舊是具有可行性的。

七、下一步的嘗試

其實，從化學角度分析，影響單個水分子體系能量的因素只有兩個O-H鍵長與H-O-H鍵角三個變量，採取純坐標作為輸入是一種並不經濟的做法（例如分子整體旋轉會導致坐標發生大量變化，但分子能量並無變化），因此，可否使用這三個變量作為輸入數據呢（本質上就是採取內坐標代替笛卡爾坐標輸入）？這其實就屬於數據預處理的內容了，筆者將在下一篇筆記中進行嘗試。

八、參考數據集下載

Weida

2022年8月20日