傅碧娜，陳 俊，劉天輝，邵科杰，張東輝

中國科學(xué)院大連化學(xué)物理研究所，分子反應(yīng)動力學(xué)國家重點實驗室，理論與計算化學(xué)研究中心，遼寧 大連 116023

1 引言

化學(xué)反應(yīng)基本上是一個舊化學(xué)鍵斷裂，新化學(xué)鍵形成的過程。從最基本和微觀水平研究化學(xué)反應(yīng)已經(jīng)成為現(xiàn)代物理化學(xué)的重要目標(biāo)。分子反應(yīng)動力學(xué)是在原子、分子層次上研究化學(xué)反應(yīng)微觀動態(tài)和機(jī)理的一門學(xué)科，它所研究的基元化學(xué)反應(yīng)能夠解釋量子態(tài)分辨的反應(yīng)碰撞過程，并提供重要的反應(yīng)機(jī)理信息。隨著高靈敏性，高分辨率實驗技術(shù)和動力學(xué)理論研究的發(fā)展和應(yīng)用，反應(yīng)動力學(xué)領(lǐng)域取得了許多有重要意義的成果，并且已經(jīng)從一些研究很徹底的三原子反應(yīng)過渡到真正的多原子反應(yīng)1–30。

多原子反應(yīng)動力學(xué)理論研究在兩個方面受到了很大的挑戰(zhàn)，一是勢能面，二是動力學(xué)。在Born-Oppenheimer近似下，我們首先通過分離電子的運動和原子核的運動，采用從頭算得到不同原子核構(gòu)型下的勢能值，根據(jù)這些分立的勢能值構(gòu)造反應(yīng)的勢能面；然后在構(gòu)造好的勢能面上求解核運動Schr?dinger方程，就能得到反應(yīng)的動力學(xué)信息，如反應(yīng)幾率，微分和積分反應(yīng)截面，以及反應(yīng)速率常數(shù)等，這個過程也就是我們常說的動力學(xué)計算。由此可見，構(gòu)造一個精確可靠的勢能面對于反應(yīng)動力學(xué)的理論研究是頭等重要的，它是我們正確理解一個化學(xué)反應(yīng)機(jī)理的前提。

近幾十年來，隨著電子結(jié)構(gòu)理論和從頭算數(shù)據(jù)擬合技術(shù)的發(fā)展，使得精確的多維全域勢能面的構(gòu)造成為可能31–40。構(gòu)造一個精確的勢能面需要大量的高精度從頭算，這些從頭算能量點分布在一個很大的構(gòu)型空間，另外需要一個可信的函數(shù)形式來表征這些分立的能量點。由于計算機(jī)能力的不斷提升，高精度從頭算理論比如耦合簇(CC)41或者多參考組態(tài)相互作用(MRCI)42方法，結(jié)合較大的基組，基本上可以給出可靠的從頭算能量。然而，基于這些分立的能量點，如何構(gòu)造多維分子體系可靠的全域勢能面還是一個很大的挑戰(zhàn)。

樣條插值方法以及基于多體展開和函數(shù)形式的非線性擬合方法在三原子體系中運用是相當(dāng)成功的43,44。對于四原子以上的反應(yīng)體系，體系自由度的迅速增加使得如何能精確有效地構(gòu)造多原子反應(yīng)勢能面成為了動力學(xué)領(lǐng)域的關(guān)鍵問題和難點。如果按每個自由度分布固定數(shù)據(jù)點的方法來做樣條插值，就需要非常大的從頭算計算量，這對目前的計算機(jī)水平是無法實現(xiàn)的。移動插值最小二乘法(IMLS)31和改進(jìn)的Shepard插值方法32是典型的插值方法，用來構(gòu)造多原子分子體系勢能面。我們組曾經(jīng)用改進(jìn)的Shepard插值方法構(gòu)造了CH5和 OH3體系全維勢能面23,45。這些插值勢能面的缺點是計算量很大，很大程度上限制了動力學(xué)勢能值計算的速度。Bowman等人發(fā)展了置換不變多項式(PIP)擬合方法33,34，在構(gòu)造高維勢能面的工作上前進(jìn)了重要的一步。置換不變多項式擬合方法使用首要不變量(首要多項式)和次要不變量(次要多項式)作為基函數(shù)，生成所需的置換不變多項式，最后通過線性最小二乘法獲得各個多項式的系數(shù)。這個方法所用到的多項式基組考慮了體系的置換不變對稱性，是一種有效的復(fù)雜多原子反應(yīng)體系的全維勢能面構(gòu)造方法，并被成功應(yīng)用于一系列復(fù)雜體系的勢能面構(gòu)建46,47。

近幾年，神經(jīng)網(wǎng)絡(luò)(NN)方法48被廣泛應(yīng)用到多原子分子體系的高精度勢能面構(gòu)建中。神經(jīng)網(wǎng)絡(luò)的函數(shù)形式非常靈活，因此原理上可以利用NN擬合任意一種函數(shù)形式，而且擬合精度也能很高。在過去的二十年，神經(jīng)網(wǎng)絡(luò)方法在氣相分子與金屬表面反應(yīng)以及氣相體系勢能面的構(gòu)建中有不少研究49–55。然而，怎樣系統(tǒng)地選擇新的構(gòu)型用來從頭算，進(jìn)而開展有效的勢能面擬合以及動力學(xué)計算，最后得到可靠的動力學(xué)結(jié)果，一直以來是困難的工作。近幾年，我們組完善和發(fā)展了一套系統(tǒng)的方案用于有效地選擇構(gòu)型36,37,56，根據(jù)此方案，可以流程化地選擇足夠多且廣泛分布在勢能面重要區(qū)域的數(shù)據(jù)點集。最后能讓神經(jīng)網(wǎng)絡(luò)擬合得到的勢能面精度達(dá)到量子動力學(xué)的需求，并且得到可靠的量子動力學(xué)結(jié)果。接著，Guo等人提出了PIPNN方法38,39，他們不再使用分子的核間距作為神經(jīng)網(wǎng)絡(luò)的輸入層，而是采用 Bowman等提出的置換不變多項式，作為神經(jīng)網(wǎng)絡(luò)的輸入層。PIP-NN擬合方法將相同原子的置換對稱性嚴(yán)格考慮到神經(jīng)網(wǎng)絡(luò)方法擬合中。這個方法理論上要求所使用的多項式包含所有的首要不變量和次要不變量。但隨著體系的增大，輸入層多項式的個數(shù)隨次要不變量的最高次數(shù)呈非線性增長。

神經(jīng)網(wǎng)絡(luò)是一個較為理想的通用型擬合方法，若能在保證對稱性的前提下有效減少神經(jīng)網(wǎng)絡(luò)輸入層多項式的個數(shù)，那么神經(jīng)網(wǎng)絡(luò)就有可能用于更為復(fù)雜的體系的勢能面構(gòu)建上。為此，我們組最近提出了基本不變量神經(jīng)網(wǎng)絡(luò)(FI-NN)擬合方法40。FI-NN使用基本不變量作為神經(jīng)網(wǎng)絡(luò)的輸入，在保證勢能面的對稱性的同時能有效提高勢能面的計算速度。在PIP及PIP-NN方法中，由于多項式的數(shù)目增長迅速，很多情況下兩者需要使用按次數(shù)截斷的多項式。對于較小的分子體系，使用PIP和PIP-NN方法擬合勢能面時，通過選擇合適的截斷次數(shù)，可以得到較好的擬合結(jié)果，并同時確保所使用的多項式數(shù)目不至于過多。但當(dāng)分子體系增大時，只能選擇較低階的多項式，這會使勢能面的擬合結(jié)果達(dá)不到理想水平。同時，舍棄高階多項式意味著擬合時所使用的生成元不完備，這也將導(dǎo)致額外誤差的引入。為了減少神經(jīng)網(wǎng)絡(luò)輸入層多項式的數(shù)目，并同時保證生成元的完備性，我們可以使用基本不變量作為神經(jīng)網(wǎng)絡(luò)的輸入，這就是基本不變量神經(jīng)網(wǎng)絡(luò)。理論上，通過增加隱藏層神經(jīng)元的個數(shù)，神經(jīng)網(wǎng)絡(luò)可以用于擬合任意形式的函數(shù)。又由于基本不變量是完備的生成元，所以基本不變量神經(jīng)網(wǎng)絡(luò)可以用于擬合任意形式的具有交換不變性的勢能面?？傮w來說，F(xiàn)INN方法在數(shù)學(xué)形式上是精確而簡潔的，同時由于輸入層限制的減少，其神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)也更為合理。

本文主要介紹我們組在近幾年基于神經(jīng)網(wǎng)絡(luò)方法的勢能面構(gòu)造的進(jìn)展。首先簡單介紹基于神經(jīng)網(wǎng)絡(luò)的勢能面構(gòu)造方法，然后詳細(xì)討論這些方法在氣相多原子反應(yīng)，以及氣相分子在金屬表面反應(yīng)的應(yīng)用。這些反應(yīng)體系包括OH3，HOCO，CH5氣相反應(yīng)體系，以及 H2O在 Cu(111)表面的解離等。

2 神經(jīng)網(wǎng)絡(luò)(NN)

前饋型神經(jīng)網(wǎng)絡(luò)是一種較為簡單的神經(jīng)網(wǎng)絡(luò)函數(shù)，也是常用于化學(xué)反應(yīng)勢能面構(gòu)造的類型。我們組都是采用前饋型神經(jīng)網(wǎng)絡(luò)進(jìn)行勢能面擬合。圖 1(a)是具有兩個隱藏層的前饋型神經(jīng)網(wǎng)絡(luò)函數(shù)的結(jié)構(gòu)示意圖，圖中的每一個圓形節(jié)點稱為一個神經(jīng)元。如果輸入層(input layer)、第1隱藏層(1st hidden layer)、第2隱藏層(2nd hidden layer)、輸出層(output layer)的神經(jīng)元數(shù)目分別為I、J、K、L，那么這個神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)為I-J-K-L。圖1(b)是其中一個神經(jīng)元的運算規(guī)則。對于第m層的第i個神經(jīng)元，它有對應(yīng)的輸入值，輸出值，以及激發(fā)函數(shù)fm。如圖1所示，神經(jīng)網(wǎng)絡(luò)第1隱藏層的輸入與輸出分別為：

神經(jīng)網(wǎng)絡(luò)第2隱藏層的輸入與輸出分別為：

神經(jīng)網(wǎng)絡(luò)最終的輸出層為：

圖1 (a) 具有兩個隱藏層的前饋型神經(jīng)網(wǎng)絡(luò)函數(shù)結(jié)構(gòu)示意圖；(b) 隱藏層中神經(jīng)元的運算規(guī)則Fig.1 (a) Illustration of the functional structure of a feed forward neural networks;(b) the computing rule of a neuron in the hidden layer.

當(dāng)神經(jīng)網(wǎng)絡(luò)應(yīng)用到分子反應(yīng)勢能面的擬合時，神經(jīng)網(wǎng)絡(luò)的輸入層對應(yīng)分子構(gòu)型的坐標(biāo)，一般采用分子的內(nèi)坐標(biāo)或者原子核間距來表示，也可以使用精心構(gòu)造的對稱化坐標(biāo)。在擬合的過程中，我們會對隱藏層選取不同的神經(jīng)元個數(shù)來進(jìn)行嘗試擬合得到一個好的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，因為神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)對于擬合的好壞有很重要的影響。對于一個確定好的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，我們可以通過Levenberg-Marquardt運算法則57對公式(1)–(5)中的權(quán)重和偏差進(jìn)行調(diào)整，我們用均方根偏差(RMSE)來確定神經(jīng)網(wǎng)絡(luò)擬合的好壞，

我們在擬合過程中不斷地優(yōu)化參數(shù)，使誤差RMSE減小到足夠的精度。當(dāng)輸出的RMSE小于預(yù)先設(shè)定的閾值，或者RMSE的下降速率小于一定的標(biāo)準(zhǔn)，意味著擬合結(jié)束，我們就用神經(jīng)網(wǎng)絡(luò)得到了一個NN勢能面。值得指出的是，設(shè)定不同的隨機(jī)值神經(jīng)網(wǎng)絡(luò)會得到不同的擬合結(jié)果，這是因為神經(jīng)網(wǎng)絡(luò)擬合是通過尋找局域擬合誤差最優(yōu)的過程。在神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)確定的情況下，我們一般也會進(jìn)行多次擬合，最后取RMSE最小的一次或幾次擬合結(jié)果作為最終的勢能面。

Levenberg-Marquardt運算法則在神經(jīng)網(wǎng)絡(luò)擬合中是通過迭代來實現(xiàn)的，它的計算量和計算時間隨著要擬合的點的數(shù)目和我們所選取的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)中神經(jīng)元的數(shù)目的增加而快速上升。為了加快神經(jīng)網(wǎng)絡(luò)擬合的速度，我們一般采用對總的能量點分塊的方法(一般根據(jù)能量點構(gòu)型所處的位置如入口區(qū)，相互作用區(qū)，出口區(qū)等)，對這些部分分別進(jìn)行訓(xùn)練和測試來得到比較好的擬合，最后再把它們連在一起構(gòu)成一個全域的勢能面。這樣處理能夠減少神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)中神經(jīng)元的數(shù)目而減少計算量和擬合所需要的時間，而且能使我們方便地對各部分?jǐn)M合得到比較理想的RMSE，特別是對動力學(xué)反應(yīng)比較重要的區(qū)域如入口區(qū)和相互作用區(qū)，我們可以投入更多的時間和精力來把它們擬合得更好。

選擇合適的構(gòu)型用于從頭算是勢能面構(gòu)造的關(guān)鍵，完備的構(gòu)型集合需要覆蓋反應(yīng)涉及到的所有區(qū)域，并且在過渡態(tài)等重要區(qū)域包括更多的構(gòu)型。我們發(fā)展和完善了一套系統(tǒng)的構(gòu)型選擇方案，通過多次添加構(gòu)型和擬合，最終得到可靠的勢能面，從而計算出可靠的反應(yīng)動力學(xué)結(jié)果。我們一般通過低精度的從頭算動力學(xué)方法計算某些軌線，去搜索所有可能的反應(yīng)通道和分子能夠到達(dá)的構(gòu)型空間。并將這些軌線上的構(gòu)型經(jīng)過篩選(基于每兩個核構(gòu)型間的距離)，得到第一批構(gòu)型的集合。然后使用高精度從頭算方法計算這些構(gòu)型的能量，再使用NN方法進(jìn)行多次擬合，得到了最初版本的勢能面。得到最初版本的勢能面以后，就能很快開展大量的準(zhǔn)經(jīng)典軌線計算，在軌線經(jīng)歷的這些構(gòu)型中篩選新的構(gòu)型加入原來的從頭算數(shù)據(jù)集，進(jìn)行之后勢能面的改進(jìn)。我們需要從各個通道和不同的過渡態(tài)區(qū)域出發(fā)進(jìn)行經(jīng)典軌線計算，目的是讓勢能面能精確地描述每一個反應(yīng)通道。這樣經(jīng)過若干次重復(fù)之后，就可以進(jìn)行量子動力學(xué)計算，判斷勢能面是否收斂。

最終擬合得到的勢能面必須包括以下兩方面才是可靠的，一是擬合誤差RMSE足夠小，在多次擬合的勢能面上能計算得到一樣的動力學(xué)結(jié)果；二是用于擬合的從頭算構(gòu)型數(shù)目足夠多，并且覆蓋范圍足夠廣，再繼續(xù)增加或者刪除從頭算構(gòu)型重新擬合也能得到同樣的動力學(xué)結(jié)果。當(dāng)以上兩個要求都滿足時，表明勢能面的構(gòu)造過程已經(jīng)完全收斂，基于此勢能面能夠得到可靠的動力學(xué)結(jié)果。

2.1 OH3體系

H2+ OH ? H2O + H反應(yīng)是典型的四原子體系，它在大氣化學(xué)和燃燒化學(xué)中占有重要的位置58,59。由于OH3體系中有3個原子是H，從頭算或量子動力學(xué)計算量都相對比較小，所以是從頭算勢能面構(gòu)造和量子動力學(xué)計算的典型代表。在過去幾十年，人們用不同方法構(gòu)造了WDSE60，OC61,62，WSLFH63，YZCL264，XXZ23等該體系的勢能面。這些勢能面要么不夠精確，要么運算速度太慢。

因此，2013年，我們組使用神經(jīng)網(wǎng)絡(luò)方法，結(jié)合上述系統(tǒng)的構(gòu)型選擇方案，構(gòu)造了一個全新的從頭算勢能面(CXZ勢能面)36。CXZ勢能面的構(gòu)造將XXZ勢能面已有的8000多個構(gòu)型作為初始擬合數(shù)據(jù)集，采用UCCSD(T)-F12a/aug-ccpVTZ從頭算方法計算勢能，用多次神經(jīng)網(wǎng)絡(luò)擬合與準(zhǔn)經(jīng)典軌線計算相結(jié)合的方法來篩選構(gòu)型，最終選擇了 16814個構(gòu)型作為擬合數(shù)據(jù)集。在進(jìn)行擬合的時候，我們按照H原子與O原子的距離對3個H 原子進(jìn)行排序，使得 ROH1≤ ROH2≤ ROH3，以此來保證一個構(gòu)型在NN勢能面上能量的唯一性。

為了提高擬合效率，我們將整個勢能面分為三個區(qū)域，即(I) H2+ OH通道，(II) H + H2O通道，(III) OH3相互作用區(qū)，如圖 2中的紅色分割線所示。在分別對三個區(qū)域擬合時，為了保證連接處吻合程度，在連接處附近有一小部分構(gòu)型同時用于兩邊的區(qū)域擬合。整體勢能面則是用光滑的開關(guān)函數(shù)將這 3個單獨的擬合連接，是它們的權(quán)重之和，記為 NN1 (即 CXZ勢能面)，整體擬合誤差RMSE為1.61 meV。另外,我們也將所有16814個從頭算點同時用于一次擬合，采用6-50-80-1的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)得到了一個全域的勢能面，記為NN2。NN2勢能面的擬合誤差RMSE為1.50 meV。為了測試從頭算構(gòu)型是否足夠多，我們隨機(jī)選擇了16814個構(gòu)型中 80%的數(shù)據(jù)，即使用 13451個構(gòu)型，也按照整體擬合的方式構(gòu)建了NN3勢能面。NN1、NN2和NN3這三個勢能面對所有數(shù)據(jù)的擬合誤差隨勢能的分布如圖3所示。

我們利用量子動力學(xué)方法來測試勢能面構(gòu)造的收斂性質(zhì)?；谶@3個NN勢能面，我們使用全維量子波包方法計算了OH + H2→ H2O + H反應(yīng)幾率，以及交換反應(yīng)H + HOH → H’ + HOH的反應(yīng)幾率，其中反應(yīng)物都處于振轉(zhuǎn)基態(tài)，分別如圖4和5所示?？梢钥闯?，對于這兩個反應(yīng)通道,基于分區(qū)域擬合的 NN1勢能面和整體擬合的 NN2勢能面的量子動力學(xué)結(jié)果基本沒有差別，而且，采用80%的數(shù)據(jù)點擬合出來的NN3勢能面可以得到與NN1勢能面基本一致的動力學(xué)結(jié)果。這些結(jié)果表明我們得到的勢能面對擬合的數(shù)據(jù)點和擬合過程都能很好地收斂。基于這樣的勢能面，我們最終能得到可靠的動力學(xué)結(jié)果。

圖2 所有用于從頭算構(gòu)型在O-H2和O-H3坐標(biāo)空間上的分布，并以紅色分割線將整體區(qū)域分成OH + H2，H + H2O和OH3相互作用區(qū)這三個部分36Fig.2 The spacial distribution of ab initio data points as a function of O-H2 and O-H3,with the red dividing lines for the OH + H2,H + H2O,and OH3 interaction parts 36.

圖3 在NN1，NN2和NN3勢能面上所有構(gòu)型的擬合誤差隨著相對于OH + H2平衡構(gòu)型的從頭算能量的分布36Fig.3 The fitting errors for all the data points in NN1,NN2 and NN3 PESs,as a function of their corresponding ab initio energies with respect to OH + H2 36.

值得一提的是，CXZ(NN1)勢能面的計算速度比XXZ勢能面快挺多，并且要比XXZ勢能面光滑，因此它代表目前OH3反應(yīng)體系最精確的基態(tài)勢能面。進(jìn)一步的研究證明，基于CXZ勢能面的量子動力學(xué)計算結(jié)果能與交叉分子束實驗結(jié)果高度吻合。我們的勢能面能夠滿足高精度實驗和理論研究的需求。

圖4 用六維含時波包法在NN1、NN2和NN3勢能面上，在體系的總角動量Jtot = 0時計算得到的H2 + OH → H2O + H的反應(yīng)幾率36Fig.4 Reaction probabilities of H2 + OH → H2O + H from six-dimensional time dependent wave packet calculations on NN1,NN2 and NN3 PESs at the total angular momentum Jtot = 0 36.

圖5 用六維含時波包法在NN1、NN2和NN3勢能面上，在體系的總角動量Jtot = 0時計算得到的H2O + H →H2O + H′的反應(yīng)幾率36Fig.5 Reaction probabilities of H2O + H → H2O + H′from six-dimensional time dependent wave packet calculations on NN1,NN2 and NN3 PESs at the total angular momentum Jtot = 0 36.

2.2 HOCO體系

OH + CO → H + CO2反應(yīng)在大氣化學(xué)和燃燒化學(xué)中都扮演著重要的角色65,66。它與OH + H2直接反應(yīng)不同，是一個典型的在反應(yīng)路徑上存在長壽命中間體的復(fù)雜反應(yīng)。該反應(yīng)過程中存在較穩(wěn)定的HOCO順式和反式結(jié)構(gòu)，并且入口(OH + CO)與出口(H + CO2)通道都存在明顯的勢壘。由于HOCO體系包含3個非氫原子，并且反應(yīng)途徑上存在長壽命的中間體，使得這個體系的勢能面構(gòu)造和量子動力學(xué)計算都非常困難。

早期的一些勢能面都不能可靠完整地描述這個體系67–71。2012 年，Li與 Guo 等使用 UCCSD(T)-F12b/aug-cc-pVTZ方法計算了大約35000個從頭算能量點，通過置換不變多項式(PIP)擬合的方法構(gòu)造了HOCO體系的全域勢能面(CCSD-1/d和改進(jìn)的CCSD-2/d勢能面)72。這兩個PIP勢能面比以前的勢能面要更加可靠，并且被用來做一些動力學(xué)計算。2013年，我們基于大約 70000個UCCSD(T)-F12b/aug-cc-pVTZ能量點，用神經(jīng)網(wǎng)絡(luò)方法構(gòu)造了一個擬合精度更高的勢能面37。這個HOCO NN勢能面的構(gòu)建過程跟上述的OH3體系相似。

我們首先在早期的 LTSH勢能面上開展了少量的準(zhǔn)經(jīng)典軌線計算，通過判斷距離的方式選擇了第一批6000個數(shù)據(jù)點，進(jìn)行高精度的從頭算。從這些構(gòu)型出發(fā)，通過多次 NN擬合和準(zhǔn)經(jīng)典軌線計算的方式來產(chǎn)生新的構(gòu)型，構(gòu)型篩選過程使用了與OH3體系相同的方案。基于NN勢能面進(jìn)行了大量的量子動力學(xué)計算來判斷勢能面的擬合過程是否收斂。最終的 NN勢能面一共使用了～74400個從頭算構(gòu)型。

為了減少NN擬合計算量，提高擬合精度，我們又采用了 OH3體系中使用的分區(qū)域擬合的方案，將全域勢能面分為HO + CO入口區(qū)域，HOCO相互作用區(qū)域，以及H + CO2出口區(qū)域，并對這三個區(qū)域分別擬合，分別測試收斂性。同時，為了提高低能量區(qū)域的精確性，在使用 Levenberg-Marquardt算法進(jìn)行NN擬合的過程中，我們將高能部分適當(dāng)減少權(quán)重，以保證低能部分更加精確。由于相互作用區(qū)的勢能面形狀非常復(fù)雜，我們使用多個勢能面作平均的方式進(jìn)一步降低擬合誤差。最終的勢能面在整體的能量范圍內(nèi)，擬合誤差RMSE為9.2 meV。

我們用量子動力學(xué)計算了基于兩個平均勢能面以及去掉10%從頭算點后擬合的勢能面，HO +CO → H + CO2反應(yīng)的總反應(yīng)幾率，如圖6所示。從圖中比較發(fā)現(xiàn)，反應(yīng)幾率基本相似，從而說明勢能面相對擬合數(shù)據(jù)點和擬合方法都收斂得比較好。

相對于以前的勢能面，NN勢能面不僅使用了更大量的從頭算數(shù)據(jù)，而且覆蓋了所有的反應(yīng)通道，在每個反應(yīng)通道上的擬合精度都相當(dāng)高，是HOCO體系最為可靠的勢能面。

圖6 (a) OH + CO → H + CO2 反應(yīng)在2套勢能面上的反應(yīng)幾率比較；(b) 在其中一套勢能面上的反應(yīng)幾率和去掉10%從頭算點重新擬合的勢能面上計算得到的反應(yīng)幾率37Fig.6 (a) Comparison of total reaction probabilities for OH + CO → H + CO2 on two sets of PESs averaged over three fittings and (b) comparison of total reaction probabilities on one of the sets shown in (a) and a PES averaged over three fitting based on 10% less data points 37.

2.3 CH5體系

H + CH4→ H2+ CH3以及其逆反應(yīng)在甲烷燃燒過程中非常重要57。作為 6原子體系反應(yīng)的代表，人們對這個體系的研究已有很長的歷史，這些研究大大提高了我們對多原子反應(yīng)的認(rèn)識。因為這個反應(yīng)的6個原子中有5個H原子，所以成為很好的測試勢能面和各種動力學(xué)方法的典型。

科學(xué)家們?yōu)榱藰?gòu)建這個體系的精確勢能面做出了巨大的努力。然而早期的半經(jīng)典勢能面都不能定量描述該反應(yīng)73,74。2004年，Manthe等基于CCSD(T)從頭算能量，采用改進(jìn)的Shepard插值方法構(gòu)造了第一個比較可靠的勢能面75，并且用該勢能面精確地計算了反應(yīng)的速率常數(shù)。但是這個勢能面是個局域的勢能面，主要限制在過渡態(tài)區(qū)域，不能被應(yīng)用于詳細(xì)的全域動力學(xué)計算。Bowman等利用置換不變多項式方法，基于RCCSD(T)/augcc-pVTZ水平的高精度從頭算數(shù)據(jù)，構(gòu)建了ZBB1、ZBB2和ZBB376,77等一系列全域勢能面。這些勢能面包括了抽取和交換反應(yīng)通道，比之前的勢能面都更加精確。我們利用改進(jìn)的Shepard插值方法，基于UCCSD(T)/aug-cc-pVTZ水平的從頭算數(shù)據(jù)，構(gòu)造了全維的 ZFWCZ 勢能面45。測試表明，ZFWCZ 勢能面精確度 ZBB3 勢能面差不多，并且基于這兩個勢能面的量子動力學(xué)結(jié)果也基本一致。ZFWCZ勢能面作為高維插值型的勢能面，其計算速度比較慢，很難應(yīng)用于進(jìn)一步的動力學(xué)研究。

2014年，我們使用神經(jīng)網(wǎng)絡(luò)擬合方法，基于UCCSD(T)-F12a/aug-cc-pVTZ的從頭算數(shù)據(jù)，構(gòu)造了一個精確的勢能面。這個勢能面稱為XCZ勢能面54。XCZ勢能面的構(gòu)造采用ZFWCZ勢能面和ZBB3勢能面用到的50000多個構(gòu)型，進(jìn)而采用構(gòu)型之間的歐氏距離進(jìn)行篩選，保留了約8000個構(gòu)型用于初始數(shù)據(jù)集。坐標(biāo)采用15個鍵長作為神經(jīng)網(wǎng)絡(luò)的輸入層，并且將5個H原子排序使得C―H鍵長呈遞增關(guān)系。為了提高擬合和構(gòu)型選擇的效率,將整個勢能面分為如圖7所示的H + CH4漸近區(qū)(asymptotic region)、H2+ CH3漸近區(qū)(asymptotic region)以及 CH5相互作用區(qū)(interaction region)。與OH3體系類似，我們循環(huán)地增加新構(gòu)型，并進(jìn)行量子動力學(xué)結(jié)果測試。最終在全部的空間內(nèi)一共選擇了47783個構(gòu)型，它們的空間分布如圖7所示，這些構(gòu)型同時覆蓋了抽取和交換反應(yīng)通道。

圖7 H + CH4勢能面的從頭算點空間分布，其中H +CH4，H2 + CH3和CH5相互作用區(qū)分別由分割線隔開56Fig.7 The spacial distribution of ab initio data points as a funcion of RC-H4 and RC-H5,with the dividing lines for the H + CH4,H2 + CH3,and CH5 interaction parts 56.

為了進(jìn)一步降低擬合誤差，我們使用Agrafiotis曾提出的神經(jīng)網(wǎng)絡(luò)集合的方式78，將若干次擬合的勢能面取平均從而減小隨機(jī)誤差。我們?nèi)MSE最小的5次擬合的平均，作為此反應(yīng)體系的最終版本神經(jīng)網(wǎng)絡(luò)勢能面XCZ。最后，我們?yōu)榱藴y試從頭算構(gòu)型的數(shù)目是否足夠,只選擇了所有構(gòu)型中 90%的數(shù)據(jù)，擬合得到了另一個勢能面。這個勢能面與XCZ勢能面的量子動力學(xué)結(jié)果如圖 8所示，兩者基本沒有差異，表明勢能面已經(jīng)收斂得很好。

圖8 基于所有從頭算點擬合的NN勢能面和基于90%從頭算點擬合的勢能面，H + CH4 → H2 + CH3的總反應(yīng)幾率比較56Fig.8 The reaction probabilities for the H + CH4 →H2 + CH3 reaction on the NN PES and another NN PES fitted with 90% of all data points 56.

2.4 H2O + Cu(111)體系

氣相分子在過渡態(tài)金屬表面的解離吸附在異相催化反應(yīng)中起著非常重要的作用，是很多化學(xué)反應(yīng)的決速步驟79。但是由于在構(gòu)建勢能面和發(fā)展新的理論方法存在很多困難，只能在非常少的體系上面開展量子動力學(xué)計算80,81。水在過渡態(tài)金屬表面的解離吸附是蒸汽重整和水煤氣轉(zhuǎn)化過程中的重要一步82。基于剛性表面近似，H2O在金屬表面的解離總共需要考慮 9個自由度。這對構(gòu)建全維勢能面和開展全維量子動力學(xué)計算是很大的挑戰(zhàn)。

我們基于大量的DFT能量點，采用神經(jīng)網(wǎng)絡(luò)方法首次構(gòu)建了H2O + Cu(111)體系擬合精度很高的全維(9維)勢能面83,84。我們總共計算了 81102個DFT能量點，采用基于贗勢和平面波的VASP軟件包(Vienna ab initio simulation package)。DFT計算中采用廣義梯度近似(GGA)下的 PW91泛函形式來描述交換相關(guān)能。

由于Cu(111)表面具有C3v對稱性，我們只需要計算由top、fcc、hcp位構(gòu)成的最小不可約三角形內(nèi)的構(gòu)型的能量值，三角形外面的構(gòu)型的擬合的更好能量可以通過三角形內(nèi)部的構(gòu)型對稱得到。我們只擬合三角形內(nèi)部的數(shù)據(jù)點，考慮到特殊位置如top、bridge、fcc、hcp位置擬合的精確性以及邊界的連續(xù)性，我們特意在這四個位置以及三角形的邊界線上多選取了一些能量點，同時為了把邊界擬合得更好，我們把三角形邊界外部非?？拷吔绲囊徊糠帜芰奎c也包含進(jìn)來，這部分能量點的能量不需要重新計算，只需要按對稱性從三角形內(nèi)部對稱出來就行。這樣擬合總共包含了93908個數(shù)據(jù)點，采用H2O分子的9個笛卡爾坐標(biāo)作為神經(jīng)網(wǎng)絡(luò)擬合的輸入層。

為了提高神經(jīng)網(wǎng)絡(luò)擬合的速度，我們把全部的數(shù)據(jù)點分成四個部分：漸進(jìn)區(qū)(I)、相互作用區(qū)(II + III)和產(chǎn)物區(qū)(IV)，每個部分都有重疊。這樣分塊以后，每部分的能量點相比總的能量點以及需要的神經(jīng)元數(shù)目大為減少，在不損失擬合精度的條件下加快了擬合的進(jìn)度。對于每一部分，我們一般都要擬合幾十次。為了測試這四個部分?jǐn)M合的收斂性，我們每一部分選取三個RMSE比較小的擬合，它們一般是神經(jīng)網(wǎng)絡(luò)擬合的結(jié)構(gòu)不同或者結(jié)構(gòu)相同的時候，權(quán)重和偏差不一樣。這十二個擬合可以組成81個勢能面，我們把a(bǔ)1 + b1 + c1 +d1標(biāo)示成NN1勢能面，它的I，II，III，IV四部分的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)分別是9-80-75-1，9-80-75-1，9-80-75-1和9-85-70-1，RMSE分別是3.2 meV，10.1 meV,9.8 meV 和 13.1 meV。NN1 勢能面總的RMSE 為9.0 meV，對于對動力學(xué)研究比較重要的部分(能量小于2.0 eV)RMSE 只有6.0 meV。我們又從剩下的80個勢能面中隨機(jī)選取了三個，用7維含時波包法分別計算了top位和bridge位基于這三個勢能面和 NN1勢能面的基態(tài)反應(yīng)幾率，如圖9所示，可以看到，隨著碰撞能的增加，反應(yīng)幾率是上升的，而且這三個勢能面的反應(yīng)幾率和NN1勢能面的反應(yīng)幾率基本上是重合的，也就是說NN1勢能面對于神經(jīng)網(wǎng)絡(luò)擬合的過程是收斂的。

圖10顯示的是H2O的質(zhì)心處于TS、Top、Bridge、Hcp位時的NN1勢能面等高線圖?？梢钥闯?，在這三個位置得到的L形狀的反應(yīng)路徑是非常光滑的，并且從勢壘的位置說明 H2O在Cu(111)解離是一個晚期勢壘體系，意味著振動激發(fā)比平動能的增加對于反應(yīng)會有更大的促進(jìn)作用。我們在這個高精度擬合的全維全域勢能面上成功在全維(9維)水平計算了H2O在Cu(111)表面的解離吸附幾率，從而首次實現(xiàn)了一個三原子分子在固體表面反應(yīng)的全維量子動力學(xué)研究，證明以前的減維模型會給計算帶來較大的誤差，代表著分子-表面量子散射動力學(xué)研究的一個重要突破82。

圖9 H2O處于振轉(zhuǎn)基態(tài)時，在top位和bridge位的7維解離幾率，其中四個勢能面從總共的81個勢能面隨機(jī)挑選產(chǎn)生83Fig.9 The seven-dimensional dissociation probabilities for H2O initially in the ground rovibrational state at fixed top and bridge sites calculated on four PESs we selected randomly from the total of 81 PESs (a,b,c,and d denote part I,II,III,and IV,respectively) 83.

圖10 在固定位置上，H2O在Cu(111)表面解離吸附的勢能面等高線圖83Fig.10 Fixed-site contour plots of the PES as a function of the vertical distance of H2O (Zcom) and the distance between the dissociating H atom and the center mass of OH (r2),with other coordinates fixed at the corresponding saddle-point geometries.The saddle point geometries are inserted in the right upper corner 83.

3 置換不變多項式神經(jīng)網(wǎng)絡(luò)(PIP-NN)

對于含有相同原子的化學(xué)反應(yīng)，勢能值應(yīng)該對相同原子的置換保持不變。因此，置換對稱性是勢能面構(gòu)建中需要考慮的重要因素。由于 NN函數(shù)關(guān)于相同原子的置換是不對稱的，在兩個或者兩個以上相等核間距的情況下，勢能值會出現(xiàn)不光滑的情況。這可能在準(zhǔn)經(jīng)典軌線計算中引入一些誤差，但對量子反應(yīng)動力學(xué)結(jié)果不會造成任何影響。我們在上述 NN的擬合過程中使用了對 H原子排序的方式來保證每個構(gòu)型能量的唯一性。

為了嚴(yán)格處理這種相同原子的交換對稱性問題，Jiang與Guo提出了PIP-NN方法38,39。他們從分子的鍵長出發(fā)，產(chǎn)生Bowman提出的一套置換不變多項式33,34，然后將這些多項式作為NN的輸入層用于擬合。NN的輸入層由低階的置換不變多項式組成，也就是對稱化的單項式，即其中pij= exp(-arij)，rij是核間距，lij是pij的階數(shù)，是對稱算符。是每個單項式的總階數(shù)。值得提出的是，輸入層中置換不變多項式的數(shù)量要足夠多，以保證體系的置換對稱性。基于PIP-NN方法，我們重新構(gòu)造了HOCO體系和CH5體系的PIP-NN勢能面85,86。這兩個PIP-NN勢能面更加適合于準(zhǔn)經(jīng)典軌線計算。盡管在構(gòu)建PIP-NN勢能面時又增加了一些從頭算構(gòu)型，但是量子動力學(xué)結(jié)果與之前NN勢能面上的結(jié)果都吻合得非常好。

并且，我們還和復(fù)旦大學(xué)徐昕教授等合作，利用新一代泛函XYG387,88計算所有結(jié)構(gòu)的能量，并用 PIP-NN方法擬合了一個新的 XYG3勢能面89。我們發(fā)現(xiàn)XYG3能基本達(dá)到CCSD(T)的精度，計算得到的反應(yīng)勢壘和CCSD(T)結(jié)果很接近?；?XYZG3勢能面的量子動力學(xué)反應(yīng)幾率和速率常數(shù)與 CCSD(T)勢能面的結(jié)果也很吻合。由于XYG3泛函遠(yuǎn)比 CCSD(T)計算速度快，可以將XYG3應(yīng)用到更大的體系，尤其是H的一系列抽取反應(yīng)。

4 基本不變量神經(jīng)網(wǎng)絡(luò)(FI-NN)

在PIP-NN方法中，神經(jīng)網(wǎng)絡(luò)使用置換不變多項式作為網(wǎng)絡(luò)的輸入層，理論上要求所使用的多項式包含所有的首要不變量和次要不變量。但隨著體系的增大，輸入層多項式的個數(shù)隨次要不變量的最高次數(shù)呈非線性增長。以簡單的A3B2體系為例，次要不變量的最高次數(shù)為 11，此時共有14984個多項式的次數(shù)小于等于 11。將這 14984個多項式全部作為神經(jīng)網(wǎng)絡(luò)的輸入將導(dǎo)致網(wǎng)絡(luò)參數(shù)過多，所以并不切合實際。所以在很多時候，PIPNN方法使用按一定次數(shù)截斷的多項式作為神經(jīng)網(wǎng)絡(luò)的輸入。對于A3B2體系，可以將次數(shù)限制在6次以下，這樣便可以將輸入層的多項式減少到525個。在A3B2體系中，由于截斷次數(shù)小于次要不變量的最高次數(shù)，擬合時有一半左右的次要不變量將不能被考慮在內(nèi)。

為此，我們提出了基本不變量神經(jīng)網(wǎng)絡(luò)(FINN)擬合方法40。這個方法能在保證對稱性的前提下有效減少神經(jīng)網(wǎng)絡(luò)輸入層多項式的個數(shù)，在保證勢能面的對稱性的同時能有效提高勢能面的計算速度。

對于AiBj···Xp分子體系，我們使用原子核間距(r1,r2,···,rn)作為勢能面的自變量，其中 n 是分子體系的鍵長個數(shù)。群 G = Si× Sj× ··· × Sp，是多個對稱群的直積，其中Sn是n階對稱群。定義?g為置換操作，為群G的某個元素。

假設(shè)R是一個集合，包含所有關(guān)于r = (r1,r2,···,rn)的多項式。當(dāng)R在加法和乘法運算下為阿貝爾群，且乘法對加法滿足左右分配律時，R被稱作多項式環(huán)。設(shè)RG是R的子環(huán)，其元素在?g的作用下保持不變，則RG是一個置換不變多項式環(huán)。RG是有限生成的并且有一個最小生成元集合，其最小生成元集合也可以稱作基本不變量(FI)。有了基本不變量，我們就可以生成任意形式的置換不變多項式。下面我們將首先給出基本不變量的計算方法。

對于特定的分子體系，其基本不變量可以使用 Singular等軟件計算90，計算時所用的算法是由King在2013 年提出的91。這里我們以A2B和A3B 體系為例，使用Singular計算其基本不變量，更多體系的基本不變量可在https://github.com/kjshao/FI中獲得。

(1) 對于 A2B分子體系，基本不變量為(xi代表分子核間距)：

(2) 對于 A3B分子體系，基本不變量為(xi代表分子核間距)：

使用PIP-NN方法擬合勢能面時，通過選擇合適的截斷次數(shù)，可以得到較好的擬合結(jié)果，并同時確保所使用的多項式數(shù)目不至于過多。但當(dāng)分子體系增大時，只能選擇較低階的多項式，這會使是勢能面的擬合結(jié)果達(dá)不到理想水平。同時，舍棄高階多項式意味著擬合時所使用的生成元不完備，這也將導(dǎo)致額外誤差的引入。相較于 PIP-NN 方法，F(xiàn)I-NN中輸入層的多項式數(shù)目明顯減少，將有效提高勢能面的解析速度。

用PIP-NN方法構(gòu)建的OH3體系的勢能面，輸入層總共使用了50個多項式，最高次數(shù)為4。但實際上，A3B體系的基本不變量只有上述的 9個，最高次數(shù)為3。新的OH3FI-NN勢能面的擬合總共用了16814個ROHF-UCCSD(T)-F12a/augcc-pVTZ個從頭算數(shù)據(jù)點，與之前用于擬合 OH3CXZ NN勢能面所用的數(shù)據(jù)點相同。最終 FI-NN勢能面的整體 RMSE、最大偏差為 1.18和 28.18 meV。并且，基于FI-NN和NN勢能面的H + H2O以及其逆反應(yīng)的反應(yīng)幾率都符合得很好，說明了FI-NN方法的可靠性以及擬合所得的勢能面的精確性。

我們還比較了FI-NN和PIP-NN勢能面的計算速度。對于較小的體系，多項式的計算基本不占時間，而神經(jīng)網(wǎng)絡(luò)的計算速度基本相同，所以整體而言，兩種方法的速度相近。FI-NN方法的優(yōu)勢主要體現(xiàn)在對較大體系的計算上，對于更大的體系，F(xiàn)I-NN 的速度提升將會更為顯著。

5 結(jié)論

精確的勢能面是反應(yīng)動力學(xué)理論研究的基礎(chǔ)，只有可靠的勢能面才能得出可靠的動力學(xué)計算結(jié)果。隨著反應(yīng)體系的增大和自由度的增加，精確的勢能面構(gòu)建給人們帶來了很大的挑戰(zhàn)92。本文主要介紹了近幾年來，用神經(jīng)網(wǎng)絡(luò)方法構(gòu)建多維體系的精確勢能面的進(jìn)展。這些勢能面包括氣相分子體系OH3,HOCO,CH5，以及H2O在Cu(111)表面的解離等。我們發(fā)展和完善了一套系統(tǒng)選擇構(gòu)型的方案，通過多次添加構(gòu)型以及擬合，最終的基于神經(jīng)網(wǎng)絡(luò)構(gòu)造的勢能面擬合精度相當(dāng)高，能夠得到收斂的量子反應(yīng)動力學(xué)結(jié)果。

我們基于置換不變多項式神經(jīng)網(wǎng)絡(luò)(PIP-NN)方法，重新構(gòu)造了HOCO體系和CH5體系的勢能面。這兩個PIP-NN勢能面更加適合于準(zhǔn)經(jīng)典軌線計算，但是量子動力學(xué)結(jié)果與之前NN勢能面上的結(jié)果都吻合得非常好。

并且，我們通過使用基本不變量作為神經(jīng)網(wǎng)絡(luò)的輸入，提出了一種新的置換不變勢能面的擬合方法，稱為基本不變量神經(jīng)網(wǎng)絡(luò)方法(FI-NN)。相較于其他方法，基本不變量的使用極大地減少了神經(jīng)網(wǎng)絡(luò)輸入層多項式的個數(shù)，有效提高了勢能面的計算速度。與其他方法相比，使用FI-NN方法可以得到很好的收斂結(jié)果，擬合誤差可以達(dá)到相同水平或更小，同時還能有效加快勢能面的計算速度。

總之，這些勢能面構(gòu)建方法的發(fā)展有望將反應(yīng)動力學(xué)理論研究推廣到更大更復(fù)雜的動力學(xué)體系。