商辰堯，張東輝

（1.中國(guó)科學(xué)院大連化學(xué)物理研究所,分子反應(yīng)動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,大連116023；2.中國(guó)科學(xué)院大學(xué),北京100049）

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

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

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

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

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

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

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

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

1 理論方法

1.1 FI-NN算法

首先，簡(jiǎn)要介紹一下FI-NN勢(shì)能面的計(jì)算流程，然后再對(duì)勢(shì)能面的解析梯度算法進(jìn)行推導(dǎo).

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

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

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

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

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

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

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

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

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

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

1.2 FI-NN解析梯度算法

現(xiàn)在開始求解勢(shì)能面在笛卡爾坐標(biāo)下的解析梯度.要得到的目標(biāo)變量是能量對(duì)笛卡爾坐標(biāo)的導(dǎo)數(shù)dE/dx.根據(jù)鏈?zhǔn)角髮?dǎo)法則可以得出：

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

式中：向量a與向量b是任意兩個(gè)向量.

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

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

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

通過自動(dòng)微分程序，可以將上述多項(xiàng)式的微分輸出為

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

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

利用求導(dǎo)法則不難得出：

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

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

2 結(jié)果與討論

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

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

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

圖1示出了在111772個(gè)兩體水分子構(gòu)型上分別使用數(shù)值梯度和解析梯度計(jì)算得到的O原子受力（Fanal.-Fnum.）的差值.可以看出，絕大部分的差值在10-4eV/nm以下，考慮到在兩體水分子相互作用能中，O原子受力的絕對(duì)值通常在±100 eV/nm范圍內(nèi)，解析梯度與數(shù)值梯度之間的誤差小于十萬分之一，可以認(rèn)為上述解析梯度算法已經(jīng)得到了正確的結(jié)果.從圖中可以看到，O原子之間的距離（ROO）越近，解析梯度與數(shù)值梯度的差值就越大.這是因?yàn)楫?dāng)兩個(gè)分子逐漸接近時(shí)，分子間的相互作用越強(qiáng)，勢(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è)水分子的平衡構(gòu)型上，數(shù)值梯度與解析梯度在O原子上的差值隨數(shù)值梯度中Δx的變化趨勢(shì).可見，10-5nm附近是Δx的最佳取值范圍，當(dāng)Δx值大于10-4nm后，數(shù)值梯度引入的誤差開始急劇增大，當(dāng)Δx值小于10-6nm后，由計(jì)算機(jī)浮點(diǎn)數(shù)引入的不穩(wěn)定性開始逐漸明顯起來.總體而言，在水的平衡構(gòu)型處，當(dāng)選擇了較為合適的Δx時(shí)，數(shù)值梯度與解析梯度的誤差處于10-5eV/nm的數(shù)量級(jí)內(nèi)，屬于可以接受的范圍.但同時(shí)需要指出的是，相比于最佳的Δx取值，隨意地選擇一個(gè)不合適的Δx可能會(huì)導(dǎo)致數(shù)十倍的誤差.因此，在單純使用數(shù)值梯度勢(shì)能面的場(chǎng)合下，有必要事先對(duì)Δx進(jìn)行掃描，從而尋找一個(gè)合適的取值.

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

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

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

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

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

3 結(jié) 論

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