基本不變量神經(jīng)網(wǎng)絡(luò)解析梯度方法的研究

商辰堯，張東輝

（1.中國科學院大連化學物理研究所,分子反應動力學國家重點實驗室,大連116023；2.中國科學院大學,北京100049）

多原子體系的勢能面構(gòu)建是理論化學研究中的一個重要課題［1～4］.根據(jù)波恩-奧本海默近似，分子體系的薛定諤方程可以近似拆分為原子核薛定諤方程與電子薛定諤方程.由于電子質(zhì)量遠小于原子核質(zhì)量，在求解電子薛定諤方程時，可以近似認為原子核位置不變，由此可以求解出給定原子核坐標下的電子勢能.以原子核坐標為自變量的電子勢能函數(shù)即為勢能面.通過采用高精度的從頭算方法，如多參考組態(tài)相互作用（MRCI）方法或耦合簇（CC）方法，再結(jié)合函數(shù)插值或擬合等方法，就可以構(gòu)造出高精度的勢能面.構(gòu)造完成之后的勢能面可以用于進行一系列的化學動力學計算，如基于量子力學的量子動力學計算、以及基于牛頓力學的準經(jīng)典軌線計算.利用精確的動力學方法，可以求解分子振轉(zhuǎn)光譜、化學反應速率常數(shù)、以及態(tài)-態(tài)反應微分截面等信息.隨著計算機性能的提升，人們開始追求更為準確的動力學模擬結(jié)果.高精度勢能面不再僅被應用于小分子體系的計算中，而逐漸開始被應用于如分子動力學模擬這種研究體系宏觀熱力學性質(zhì)的領(lǐng)域［5～7］.

在實際的動力學模擬過程中，體系往往存在一些全同原子，交換兩個相同原子的坐標后，勢能面的輸出值不發(fā)生改變，這就是勢能面的置換不變性.全同原子置換不變性是勢能面構(gòu)造時必須要解決的一個問題.解決置換不變性最簡單的方法是將原子序號按照某種給定的規(guī)則進行排序［8］，這種排序規(guī)則通常以原子間距離為依據(jù).通過排序，分子坐標被映射到一個唯一的子空間，對這個子空間中的分子構(gòu)型進行擬合，即可以得到具有置換不變性的勢能面.然而，這種方案并沒有從根本上解決勢能面函數(shù)置換不變的問題，僅是通過排序避免了原子坐標發(fā)生置換的可能性.對于具有較高對稱性的體系，坐標空間中往往存在多個對稱面.由于擬合時僅僅考慮由對稱面包圍的唯一一個子空間，函數(shù)在對稱面處的一階導數(shù)往往是不連續(xù)的.對于基于牛頓力學的動力學方法——如準經(jīng)典軌線方法，這種不連續(xù)將導致原子在對稱面處的受力發(fā)生突變，從而導致錯誤的動力學結(jié)果.為了解決多原子體系中全同原子置換不變性問題，Bowman等［9～11］發(fā)展了置換不變多項式（PIP）方法.置換不變多項式是指在原子交換操作下數(shù)值不發(fā)生改變的一組多項式.置換不變多項式可以通過啟發(fā)式方法——如單項式對稱化方法（MSA）生成，也可以借助數(shù)學中的不變量理論，通過求解首要不變量和次要不變量系統(tǒng)地生成.在PIP方法中，勢能面被構(gòu)造為置換不變多項式的線性函數(shù)，在產(chǎn)生體系的置換不變多項式之后，通過最小二乘法擬合體系的勢能，進而得到具有置換不變性的勢能面.Guo等［12～15］在PIP的基礎(chǔ)上，更進一步將置換不變多項式作為神經(jīng)網(wǎng)絡(luò)的輸入，發(fā)展出了PIP-NN勢能面.在構(gòu)造PIP獲得置換不變性的基礎(chǔ)上，利用神經(jīng)網(wǎng)絡(luò)的強大函數(shù)擬合能力，可以精確地構(gòu)造更為復雜體系的勢能面.另一方面，Behler等［16～19］發(fā)展了高維神經(jīng)網(wǎng)絡(luò)（HDNN）方法.該方法基于以下假設(shè)：體系的電子勢能只與相距較近的原子間產(chǎn)生的相互作用有關(guān).在這條假設(shè)下，體系的總勢能可以表示為每個原子貢獻的局部勢能之和，即式中每個原子的勢能只與該原子附近的化學環(huán)境相關(guān).Behler等［17］設(shè)計了一組具有平移、轉(zhuǎn)動、置換不變性的對稱函數(shù)，用于描述給定原子周圍的化學環(huán)境.在計算每個原子的化學環(huán)境之后，將對稱函數(shù)作為神經(jīng)網(wǎng)絡(luò)的輸入，從而得到每個原子的局部勢能，最終通過對局部勢能求和，得到整個體系的總勢能.

在保證勢能面精度的情況下，勢能面的計算速度會直接影響動力學模擬的效率.盡量減小勢能面模型的大小，提高計算速度，是研究人員一直以來的努力方向.盡管PIP方法能夠解決置換不變性問題，但是隨著體系的增大，多項式的項數(shù)會急劇增長，將全部的多項式用于擬合勢能面會導致勢能面的計算成本過于龐大.為了控制勢能面的計算規(guī)模，通常的辦法是按照多項式的最高次數(shù)對多項式進行截斷，從而控制多項式的總數(shù)目.然而，這種對PIP項數(shù)的限制方法具有一定的隨意性.在進行截斷前，無法預先知道哪部分多項式對勢能面的擬合更有幫助，而哪部分的多項式實際上存在冗余，應當被除去.為此，Shao等［20］和Chen等［21］從基本不變量理論出發(fā)，發(fā)展了基本不變量神經(jīng)網(wǎng)絡(luò)（FI-NN）方法，以數(shù)學理論為依據(jù)減少了置換不變多項式的項數(shù).如前所述，置換不變多項式可以通過計算首要不變量與次要不變量來生成.事實上，次要不變量可以進一步分解為可約次要不變量與不可約次要不變量，其中不可約次要不變量可以作為可約次要不變量的生成元.這就意味著，置換不變多項式環(huán)的生成元存在一個最小集合，這個集合是首要不變量和不可約次要不變量的子集，也就是所謂的基本不變量（FI）.依據(jù)計算不變量理論以及King等［22］提出的算法，基本不變量可以通過程序直接生成.將生成后的基本不變量作為神經(jīng)網(wǎng)絡(luò)的輸入，就可以構(gòu)造出FI-NN勢能面.相對于啟發(fā)式地構(gòu)建置換不變多項式，基于基本不變量的多項式天然具有項數(shù)更少的優(yōu)點，并且在數(shù)學上保證了這種項數(shù)的減少不會影響多項式對構(gòu)型空間的表達能力.

除了減小勢能面模型的規(guī)模之外，直接在勢能面調(diào)用步驟上進行優(yōu)化，也能夠提高勢能面的調(diào)用速度.在以經(jīng)典力學為基礎(chǔ)的動力學模擬——如準經(jīng)典軌線方法和分子動力學模擬方法中，計算每個原子受到的相互作用力往往是整個模擬的決速步驟.在保守力系統(tǒng)中，原子間的相互作用力等于原子間勢能對坐標偏導的負數(shù)，可以通過求解勢能面對笛卡爾坐標的偏導得到.由于勢能面的函數(shù)形式往往十分復雜，對于簡單體系，勢能面對坐標的偏導數(shù)可以簡單采用數(shù)值差分的方法近似得到.考慮到勢能在笛卡爾坐標向量x中第i個分量的偏導數(shù)，根據(jù)多元函數(shù)微分原理，可以近似地將該偏導數(shù)用數(shù)值方法求解為

式中：E為勢能面的輸出能量；x代表體系的笛卡爾坐標向量，其由體系內(nèi)每個原子的笛卡爾坐標拼接而成，其維度為3×原子數(shù)；x i代表笛卡爾坐標向量的第i個分量；Δx是數(shù)值微分時的距離間隔，通常是一個任意設(shè)置的、足夠小的正實數(shù).

由上式可見，勢能面梯度的數(shù)值解法具有求解簡單的優(yōu)點.只要能夠得到勢能面的能量值，就可以輕易地得到勢能面梯度的近似解.然而同樣可以注意到，該方法在求解時，需要對體系內(nèi)每一個原子的每一個維度調(diào)用2次勢能面，對于一個有N個原子的系統(tǒng)來說，求解一次數(shù)值梯度需要進行6N次勢能面的調(diào)用.隨著體系規(guī)模N的增大，這種調(diào)用帶來的成本也會成比例的增大.更為嚴重的是，隨著體系原子數(shù)目的增加，勢能面所描述的構(gòu)型空間顯著增大.為了維持較高的勢能面精度，研究者往往需要使用更大規(guī)模的神經(jīng)網(wǎng)絡(luò)，這就意味著體系規(guī)模的增大使單次調(diào)用勢能面的成本也相應增加.顯然，在上述兩種因素的疊加影響下，簡單采用數(shù)值方法來求解體系的能量梯度已經(jīng)無法滿足更大體系的動力學研究.因此，需要尋求更為高效的勢能面梯度求解方法.

事實上，由于在某一特定體系下，勢能面往往具有靜態(tài)的數(shù)學形式，可以預先求解出勢能面對笛卡爾坐標偏導的解析解.將勢能面的求解算法和勢能面梯度的解析公式一起作為編譯期確定的函數(shù)，在調(diào)用前預先編譯成機器指令，可以極大地提高勢能面梯度的求解效率.通過解析方法求解勢能面只需調(diào)用一次解析梯度版本的勢能面函數(shù)，函數(shù)的調(diào)用次數(shù)不再隨原子數(shù)N線性增加.盡管解析梯度勢能面函數(shù)的單次調(diào)用時間會有所增加，但是在總體上依然能提高勢能梯度的求解效率.

本文發(fā)展了基于FI-NN的解析梯度勢能面方法.通過求解FI-NN勢能面函數(shù)的解析導數(shù)，精確求解了勢能面的梯度，同時極大地提升了勢能面的調(diào)用速度.

1 理論方法

1.1 FI-NN算法

首先，簡要介紹一下FI-NN勢能面的計算流程，然后再對勢能面的解析梯度算法進行推導.

在進行動力學模擬時，通常采用笛卡爾坐標系來描述體系中原子的坐標，由于笛卡爾坐標不具備平移、轉(zhuǎn)動、置換不變性，首先將笛卡爾坐標轉(zhuǎn)換成內(nèi)坐標：

式中：r ij表示原子i與原子j之間的距離；x i和x j分別表示原子i與原子j的笛卡爾坐標.

通過計算每兩個原子間的距離，可以構(gòu)建內(nèi)坐標向量r.對于有N個原子的體系，內(nèi)坐標向量的維度為N×(N-1)/2.構(gòu)建好內(nèi)坐標之后，體系描述只剩下置換不變性需要解決，可以用FI多項式來解決這個問題，但是在此之前，一般會先將內(nèi)坐標變換成倒數(shù)或e指數(shù)形式.研究發(fā)現(xiàn)，如果將原子間核間距變換為類似于Morse力場形式，在勢能面擬合中往往會取得更好的擬合結(jié)果［9］，這可能是由于Morse形式本身已經(jīng)是分子間作用力的一種較為粗略的近似.類似地，將內(nèi)坐標變換為倒數(shù)形式也能起到與Morse型變換類似的擬合提升效果.將變換后的向量表示為y：

式中：r是內(nèi)坐標向量；λ是一個與體系有關(guān)的常數(shù)，用于控制衰減系數(shù).隨后，將向量y變換為具有置換不變性的FI多項式p：

式中：生成FI多項式的過程由算符F^表示，具體的實現(xiàn)步驟可以參考Shao等［20］和Chen等［21］的研究結(jié)果.至此，已經(jīng)得到了具有平移、轉(zhuǎn)動和置換不變性的分子構(gòu)型描述方式.

這里需要注意的是，在訓練神經(jīng)網(wǎng)絡(luò)時，通常會將由訓練集計算出的FI多項式和從頭算的能量進行歸一化處理，使得神經(jīng)網(wǎng)絡(luò)的輸入數(shù)據(jù)和目標輸出均介于[-1，1]之間.這樣確保了神經(jīng)網(wǎng)絡(luò)的輸入輸出的數(shù)值不會過大，使神經(jīng)網(wǎng)絡(luò)能有更好的擬合效果.

式中：a0為真正輸入到神經(jīng)網(wǎng)絡(luò)的輸入層數(shù)據(jù)；t為神經(jīng)網(wǎng)絡(luò)擬合時的待學習數(shù)據(jù)；p是由數(shù)據(jù)集中的某個構(gòu)型計算得到的FI多項式；pmin和pmax是與p相同維度的向量，分別代表p向量在數(shù)據(jù)集中每個維度能取到的最小值和最大值；E是由從頭算方法得到的體系勢能；Emin和Emax分別是數(shù)據(jù)集中所有從頭算能量的最小值和最大值.

與深度學習領(lǐng)域常選用深層神經(jīng)網(wǎng)絡(luò)不同的是，在勢能面構(gòu)建時，通常只選用簡單的雙隱層前饋型神經(jīng)網(wǎng)絡(luò).一方面，雙隱層的神經(jīng)網(wǎng)絡(luò)已經(jīng)足以用于表達任何復雜函數(shù)形式，因此這種神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)已經(jīng)足夠用于勢能面的構(gòu)建.另一方面，動力學計算的速度嚴重依賴于神經(jīng)網(wǎng)絡(luò)的運算速度，在保證足夠精度的情況下，會盡可能選用計算效率最高的機器學習模型.神經(jīng)網(wǎng)絡(luò)的計算公式如下：

式中：下角標代表變量所在的神經(jīng)網(wǎng)絡(luò)層數(shù)：下角標0代表輸入層，下角標1，2代表第1，2層隱藏層，下角標3代表輸出層；向量a代表神經(jīng)網(wǎng)絡(luò)的神經(jīng)元數(shù)據(jù)；向量z代表神經(jīng)網(wǎng)絡(luò)前向傳播過程中，線性變換部分的計算結(jié)果；矩陣W i是連接第（i-1）層和第i層的權(quán)重矩陣，設(shè)第（i-1）層有m個神經(jīng)元，第i層有n個神經(jīng)元，那么矩陣W i的維度為(n×m)；向量b i為偏置向量，維度為n.在這里，輸出層的數(shù)據(jù)較為特殊；由于神經(jīng)網(wǎng)絡(luò)的輸出層只經(jīng)過線性變換，因此z3就是神經(jīng)網(wǎng)絡(luò)的輸出數(shù)據(jù)，又由于輸出數(shù)據(jù)z3為標量，因此，矩陣W3實際上是一個行向量，偏置值b3是一個標量.函數(shù)f(z)為神經(jīng)網(wǎng)絡(luò)的激活函數(shù)，這個函數(shù)用于對神經(jīng)元間傳遞的數(shù)據(jù)做一個非線性變換.常用的激活函數(shù)為sigmoid函數(shù)和tanh函數(shù).在這里為了提升勢能面的調(diào)用速度，采用作為激活函數(shù)，這個函數(shù)由于不需要進行指數(shù)計算，在具有較高擬合精度的同時，調(diào)用速度稍快于上述2個函數(shù).由于在訓練神經(jīng)網(wǎng)絡(luò)時，對神經(jīng)網(wǎng)絡(luò)的待學習數(shù)據(jù)進行了歸一化處理.因此，在需要實際調(diào)用神經(jīng)網(wǎng)絡(luò)勢能面時，需要對神經(jīng)網(wǎng)絡(luò)的輸出進行逆變換，以得到勢能面的輸出能量：

以上就是從笛卡爾坐標系下的分子坐標得到分子構(gòu)型能量的完整步驟.

1.2 FI-NN解析梯度算法

現(xiàn)在開始求解勢能面在笛卡爾坐標下的解析梯度.要得到的目標變量是能量對笛卡爾坐標的導數(shù)dE/dx.根據(jù)鏈式求導法則可以得出：

式中：E為勢能面的輸出能量；x為笛卡爾坐標向量；其余變量的含義均與上文一致.任意兩個向量間的導數(shù)以雅克比矩陣方式排列：

式中：向量a與向量b是任意兩個向量.

根據(jù)該定義［式（13）］，dE/dx是一個行向量，維度為3×原子數(shù).公式中第一項dE/dp可以使用神經(jīng)網(wǎng)絡(luò)中的反向傳播算法計算：

式中:diag[a]的作用是將向量a轉(zhuǎn)化成對角矩陣,其中a是任意一個向量:

式（13）中dp/dy即為對FI多項式的求導.不同的勢能面所采用的FI多項式的形式可能各不相同，因此無法像推導神經(jīng)網(wǎng)絡(luò)反向傳播公式那樣得到一個統(tǒng)一的解析解.但是，F(xiàn)I的函數(shù)形式至少在編譯期是確定的，因此可以編寫一個自動計算FI解析梯度的程序.該程序以FI的多項式公式作為輸入，然后將FI的解析梯度輸出成一個編譯期確定的函數(shù).以一個AB3體系為例，體系中共有4個原子，因此總的內(nèi)坐標個數(shù)為6個.通過基本不變量生成方法，可以將內(nèi)坐標轉(zhuǎn)化為共有9項的基本不變量多項式.本文只截取前5項為例：

通過自動微分程序，可以將上述多項式的微分輸出為

最后來求解dy/dx，在這里僅考慮y=1/r的情況，對于y=exp(-r/λ)的情況用類似的方法不難得到.考慮向量y的第k個分量y k：

式中：r k是內(nèi)坐標向量r的第k個分量；x i1代表第i個原子的笛卡爾坐標在x軸方向上的分量，其余變量的含義與之類似.

利用求導法則不難得出：

在對一個N原子體系的內(nèi)坐標求解時，i和j的遍歷滿足1≤i＜j≤N，在這種情況下，i，j和k存在如下關(guān)系：

根據(jù)式（31）和（32）就可以求解dy/dx矩陣中的每一個元素.需要補充的是，在算法的具體實現(xiàn)中存在一些優(yōu)化的辦法.盡管可以根據(jù)上述公式編寫一個二重循環(huán)逐次地求解dy/dx矩陣中的每一個元素，然而可以注意到，上述求解公式在編譯期同樣是確定的.因此，依然可以編寫一個程序，自動生成矩陣的求解代碼.這樣，i，j和k的數(shù)值在編譯期就能確定下來，從而避免了每一步必須動態(tài)求解的時間成本.

2 結(jié)果與討論

選用水分子的多體相互作用勢能面來對解析梯度方法進行測試.對于一個由N個水分子構(gòu)成的團簇，體系總勢能可以分解為若干多體相互作用能的疊加［23-26］：

式中：E N表示N個分子的總勢能；VNB表示N體相互作用能；x i代表第i個水分子的笛卡爾坐標.根據(jù)上述方法，構(gòu)造了水分子的兩體、三體、四體勢能面.其中兩體勢能面采用UCCSD（T）/CBS方法進行從頭算［27］，三體、四體由于體系較大，采用雙雜化泛函XYGJ-OS/AVTZ進行從頭算［28］.研究表明，大于四體的水分子相互作用影響十分有限，因此并沒有構(gòu)造更高階的相互作用勢能面［29］.由式（34）可以看出，N體相互作用勢能面需要使用所有N個水分子的笛卡爾坐標，勢能面的規(guī)模勢必會隨著N值相應增大.另一方面，在計算N體相互作用能的貢獻時，需要遍歷所有由N個水分子構(gòu)成的組合，當體系的分子數(shù)較多時，這樣的組合數(shù)就會十分龐大.由于上述2個原因，降低多體相互作用勢能面的計算成本變得十分重要，這就是為什么用這一組勢能面來驗證解析梯度方法的效果.

下面來驗證上述解析梯度方法的正確性.由于數(shù)值梯度的求解十分簡單，能夠用數(shù)值方法幾乎無誤地得到一個勢能面梯度的近似值.如果解析梯度的結(jié)果與數(shù)值結(jié)果十分接近，就可以認為采用的解析梯度方法是正確的.另一方面，當確認上述推導的解析梯度算法有效之后，可以用解析梯度的結(jié)果反過來評判數(shù)值梯度的結(jié)果.根據(jù)前面所述，數(shù)值梯度是通過給笛卡爾坐標變動一個較小的位移Δx來近似求得能量梯度，然而這個Δx的最優(yōu)數(shù)值無法預先知道.理論上，Δx越小所得到的結(jié)果越接近于真實值，但由于計算機在處理浮點數(shù)時存在精度問題，過小的Δx會使最終的結(jié)果不穩(wěn)定.同樣的，如果選用的Δx過大，在能量變化十分劇烈的位置求得的數(shù)值梯度可能會與該位置的真實梯度產(chǎn)生較大偏差.通過對比數(shù)值梯度與解析梯度的結(jié)果就可以看出所選用的Δx是否合適.

圖1示出了在111772個兩體水分子構(gòu)型上分別使用數(shù)值梯度和解析梯度計算得到的O原子受力（Fanal.-Fnum.）的差值.可以看出，絕大部分的差值在10-4eV/nm以下，考慮到在兩體水分子相互作用能中，O原子受力的絕對值通常在±100 eV/nm范圍內(nèi)，解析梯度與數(shù)值梯度之間的誤差小于十萬分之一，可以認為上述解析梯度算法已經(jīng)得到了正確的結(jié)果.從圖中可以看到，O原子之間的距離（ROO）越近，解析梯度與數(shù)值梯度的差值就越大.這是因為當兩個分子逐漸接近時，分子間的相互作用越強，勢能面的變化越劇烈，在Δx不變的情況下，數(shù)值梯度引入的誤差就越明顯.

Fig.1 Errors between numerical and analytical gradient on two-body water structures

Fig.2 Errors between numerical and analytical gradient on different choices ofΔx

圖2示出了在兩個水分子的平衡構(gòu)型上，數(shù)值梯度與解析梯度在O原子上的差值隨數(shù)值梯度中Δx的變化趨勢.可見，10-5nm附近是Δx的最佳取值范圍，當Δx值大于10-4nm后，數(shù)值梯度引入的誤差開始急劇增大，當Δx值小于10-6nm后，由計算機浮點數(shù)引入的不穩(wěn)定性開始逐漸明顯起來.總體而言，在水的平衡構(gòu)型處，當選擇了較為合適的Δx時，數(shù)值梯度與解析梯度的誤差處于10-5eV/nm的數(shù)量級內(nèi)，屬于可以接受的范圍.但同時需要指出的是，相比于最佳的Δx取值，隨意地選擇一個不合適的Δx可能會導致數(shù)十倍的誤差.因此，在單純使用數(shù)值梯度勢能面的場合下，有必要事先對Δx進行掃描，從而尋找一個合適的取值.

解析梯度勢能面方法最重要的意義在于提高勢能面的調(diào)用速度.為此，對勢能面的運行時間進行了測量.由于測量的目的僅在于對比解析梯度與數(shù)值梯度的速度差異，測試程序的絕對運行時間可以根據(jù)需要進行調(diào)節(jié).在測量時，如果程序運行時間過短，測得的數(shù)據(jù)就不夠準確.反之，如果運行時間過長，則沒有必要.為了能夠方便地調(diào)節(jié)測試程序的運行時間，預先將待測試的分子構(gòu)型加載到一個數(shù)組中，再選定一個合適的迭代次數(shù)，對數(shù)組中的分子構(gòu)型進行勢能面調(diào)用.這樣，程序的運行時間與自行設(shè)定的迭代次數(shù)相關(guān)，可以任意調(diào)節(jié).為了避免計算機緩存帶來的運行速度干擾，并沒有按數(shù)組順序依次進行遍歷，而是每次生成一個隨機數(shù)，隨機地挑選數(shù)據(jù)集中的某個分子構(gòu)型進行勢能面讀取.由于勢能面的運行速度與參與計算的數(shù)據(jù)值無關(guān)，這種隨機挑選構(gòu)型的方法是合理的.

表1示出了在多個水分子多體相互作用勢能面上分別進行勢能計算和勢能梯度計算所用的時間.對于兩體勢能面，選用了具有66621個兩體水分子構(gòu)型的數(shù)據(jù)集，每次測試進行了107次迭代.除此之外，分別測試了多個具有不同結(jié)構(gòu)的三體、四體勢能面.其中，所有三體勢能面的測試均采用了139929個三體水分子的數(shù)據(jù)集，每次測試進行了106次迭代；所有四體勢能面的測試均采用了152490個四體水分子的數(shù)據(jù)集，每次測試進行了105次迭代.在這種安排下，由于三體、四體勢能面沒有選擇相同的迭代次數(shù)，二者的絕對運行時間沒有可比性，但是同屬于三體或四體的不同勢能面之間具有比較價值.由此可以看出不同勢能面結(jié)構(gòu)對勢能面運行速度的影響.

Table 1 Comparison of calculation time of forward and backward propagation on FI-NN PESs

表1的第一列表示勢能面描述的體系，分別是兩體、三體和四體相互作用能.第三列示出了勢能面采用的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu).對于三體和四體勢能面，為了得到較高的擬合精度，分別選用了不同的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)進行擬合.顯然，神經(jīng)網(wǎng)絡(luò)規(guī)模越大，勢能面計算成本也會越大.神經(jīng)網(wǎng)絡(luò)第一層的大小既表示神經(jīng)網(wǎng)絡(luò)的輸入層維度，也代表基本不變量的項數(shù).表中第四、五列示出了勢能面前向傳播與反向傳播的計算時間.前向傳播是指從分子笛卡爾坐標計算得到分子勢能所用的時間，反向傳播指求解勢能梯度所用時間.由圖2可見，反向傳播過程的耗時往往數(shù)倍于前向傳播，這是因為在反向傳播中需要計算多個向量間的雅克比矩陣.根據(jù)上述兩個時間，可以在理論上估算解析梯度方法帶來的加速效果.考慮一個N原子體系的勢能面，如果需要得到該體系的總勢能和勢能梯度，設(shè)勢能面的前向傳播耗時為t1，反向傳播耗時t2，采用解析方法計算梯度的加速比可以按如下公式計算：

表1中第六列示出了式（35）的計算結(jié)果，第七列示出了實際勢能面調(diào)用時的加速比.注意到二者數(shù)值并不完全相同.這是因為，為了測試勢能面的前向傳播與反向傳播，勢能面程序需要進行一些結(jié)構(gòu)上的調(diào)整，并在計算過程中插入計時函數(shù)，這些操作一定程度上都會影響勢能面的調(diào)用速度.但是可以看到，兩列數(shù)據(jù)具有大致相同的變化趨勢.由表1可見，盡管反向傳播的計算成本較高，相比于數(shù)值方法，采用解析方法計算勢能面梯度依然能帶來10倍以上的性能提升.并且，相同體系的勢能面采用的FI項數(shù)越多，解析梯度方法帶來的性能提升也就越大.在處理較大體系時，為了保證勢能面擬合的精度，往往需要增加FI項數(shù).因此，隨著研究體系的增大，采用解析梯度方法能夠帶來更加明顯的性能收益.

3 結(jié) 論

本文在基本不變量神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上，提出了基本不變量解析梯度勢能面的計算方法.計算解析梯度勢能面的基本思想是利用鏈式求導法則進行反向傳播.計算大致分為神經(jīng)網(wǎng)絡(luò)部分、基本不變量部分和內(nèi)坐標部分3個部分.神經(jīng)網(wǎng)絡(luò)部分的反向傳播可以通過矩陣運算得到，基本不變量部分和內(nèi)坐標部分的計算沒有統(tǒng)一的求解公式，而是由研究體系決定.本文提出的思路是利用程序?qū)υ摬糠止竭M行符號微分計算，再將計算后的結(jié)果直接輸出為可編譯的代碼.這種方法提供了一種方便、快捷的勢能面構(gòu)建思路，同時最大程度上保證了勢能面的計算速度.在誤差測試中發(fā)現(xiàn)，在全部構(gòu)型中，解析梯度與數(shù)值梯度之間的誤差不會大于十萬分之一.數(shù)值梯度結(jié)果的準確性取決于勢能面變化的劇烈程度和Δx的取值，其中不同Δx的取值會對誤差產(chǎn)生數(shù)十倍的影響.在速度測試中，發(fā)現(xiàn)對于大部分勢能面，解析梯度方法能夠提供10倍以上的加速效果，且這種速度提升會隨勢能面規(guī)模的增加而愈發(fā)明顯，使得在將來有能力研究更大規(guī)模的分子體系.