非絕熱分子動力學(xué)的量子路徑模擬

李曉克 馮偉

(天津大學(xué)物理系,天津 300350)

1 引 言

由于全量子動力學(xué)模擬計算量過大,混合經(jīng)典-量子分子動力學(xué)模擬方法已經(jīng)被廣泛應(yīng)用于氣體分子、勢能面系統(tǒng)、生物分子以及超分子系統(tǒng).在這類系統(tǒng)中,由于原子核質(zhì)量大小是電子質(zhì)量的103—104倍,電子速度比原子核速度大3—4個數(shù)量級,根據(jù)玻恩-奧本海默(BO)近似,可近似認為每一時刻電子運動在靜止的原子核勢場中,原子核幾乎不受電子位置的影響.因此在求解電子的波函數(shù)時可將核運動與電子運動分離.在該近似下,利用量子化學(xué)等電子結(jié)構(gòu)計算方法,可以得到系統(tǒng)的勢能面.若分子始終處于單一勢能面上,電子態(tài)之間的躍遷可以忽略,這就是所謂的絕熱過程.相反,如果勢能面相互靠近,在勢能面的(近)交叉區(qū)可能發(fā)生電子態(tài)的躍遷,就必須處理非絕熱過程.

混合經(jīng)典-量子分子動力學(xué)方法的基本思想是將原子核部分的運動通過牛頓運動方程來處理,而電子部分則用量子力學(xué)處理.混合經(jīng)典-量子分子動力學(xué)模擬方法主要有Ehrenfest方法[1?3]、勢能面躍遷法[4?9]、二者的混合方法[10?13]、CSDM方法(coherent switching with decay of mixing method)[14]以及最近提出的量子路徑方法[15]等.這些方法已被成功應(yīng)用于分子動力學(xué)[16?20]以及電荷傳輸研究[21].除以上混合模擬方法外,含時量子波包方法[22?24]也被廣泛應(yīng)用于三原子分子系統(tǒng)的分子動力學(xué)模擬[25?28].

本文對最近提出的量子路徑方案做進一步研究.量子路徑方案認為,由牛頓方程描述的核的經(jīng)典運動對電子的量子運動起到了量子測量的作用,應(yīng)當(dāng)用量子測量理論給出描述[29].除了在幾個典型的模型勢能面系統(tǒng)中獲得成功外,量子路徑方案還被用于模擬有機分子N2CO的光分解動力學(xué),獲得了令人滿意的結(jié)果[30].本文將進一步模擬檢驗量子路徑方案在更多的模型系統(tǒng)中的相關(guān)結(jié)果,包括單交叉模型、雙交叉模型、拓展耦合模型、啞鈴模型以及雙弓模型.由于難以在嚴格意義上得到退相干速率,我們特別模擬比較了從不同物理考慮得到的三個退相干速率公式,包括凍結(jié)高斯波包近似(frozen Gaussian approximation,FGA)[31]、能量分辨(energy discrimination,ED)[14,32]以及力分辨(force discrimination,FD)[15].發(fā)現(xiàn)對于結(jié)構(gòu)較簡單的勢能面模型,三種退相干速率都能得到較好的結(jié)果.然而對于較復(fù)雜的勢能面模型,由于復(fù)雜量子干涉的原因,與其他混合經(jīng)典-量子動力學(xué)方案類似,量子路徑方案仍然難以得到較準(zhǔn)確的結(jié)果.

2 量子路徑模擬方案

2.1 哈密頓量和量子路徑方程

在包含所有電子和原子核的總哈密頓量中除去原子核的動能部分,可以得到電子的哈密頓量:

式中,第一項描述了電子的動能,其中mj為電子質(zhì)量;第二項描述所有相互作用靜電勢能,其取決于量子部分電子坐標(biāo)r≡{rj;j=1,2,···}及經(jīng)典部分原子構(gòu)型R.選擇一組滿足正交歸一性和完備性的電子基函數(shù)?j(r,R)展開電子波函數(shù)ψ(r,R,t)=∑jcj(t)?j(r,R).通常此電子波函數(shù)的基函數(shù)選取電子哈密頓量Hel(r,R)的瞬時本征態(tài) (只考慮非簡并情況),即BO絕熱波函數(shù)Hel(r,R)?j(r,R)=εj(R)?j(r,R), 其中,εj(R)代表第j個絕熱電子態(tài)|?j(R)〉對應(yīng)的勢能面.

電子部分的動力學(xué)性質(zhì)由薛定諤方程i?˙ψ=Helψ確定. 我們定義電子哈密頓量矩陣元Vjk(R)=〈?j(r,R)|Hel(r,R)|?k(r,R)〉, 在絕熱表象下方程的形式為

根據(jù)絕熱電子態(tài)的正交性可證絕熱耦合矢量具有反厄米性質(zhì):djk(R)=?d?kj(R).與djk(R)決定了各BO勢能面之間的非絕熱耦合強度.當(dāng)兩個絕熱電子態(tài)的能量差很大時,djk(R)可以忽略,則可以回到絕熱近似下的在單一勢能面上動力學(xué).反之,則需要求解電子與原子核互相耦合的運動.

由于電子和原子核存在著耦合,電子和原子核的運動也是互相影響的.討論電子態(tài)的動力學(xué)時,通常認為原子核對電子部分的影響在于電子態(tài)是原子核構(gòu)型R的函數(shù).如果從量子測量的觀點研究非絕熱分子動力學(xué)模擬問題,經(jīng)典的原子運動也會導(dǎo)致電子部分的量子態(tài)產(chǎn)生退相干[32?36].隨著量子態(tài)的不斷演化,原子部分也將受到與量子態(tài)對應(yīng)的經(jīng)典力的影響.如果將此經(jīng)典力視為一種可以獲得電子態(tài)信息的“測量儀器”,則該儀器對電子部分的量子態(tài)進行有效的連續(xù)測量.例如,原子在不同的勢能面上會感受到不同的經(jīng)典力,相當(dāng)于測量過程中獲得了某種“結(jié)果”,從而將引起量子跳躍或量子疊加態(tài)的逐漸塌縮.這種連續(xù)獲取信息的測量過程下的動力學(xué)演化,可以由量子路徑方程描述如下[29]:

式中ρ是電子態(tài)密度矩陣,其對角元描述各個BO態(tài)的占據(jù)概率,非對角元描述它們之間的相干性.

方程(4)右端的第一項描述了電子態(tài)的相干動力學(xué)演化,第二項和第三項描述了“信息獲取”(量子測量)的反作用.其中的Lindblad超算符定義為

式 中, 測 量 算 符Mjk(R)=|?j(R)〉〈?j(R)|?|?k(R)〉〈?k(R)|, 是與態(tài)?j(R)和?k(R)相關(guān)的測量算符.Γjk描述了電子由于原子力和其他因素導(dǎo)致的測量過程中的退相干速率,γF,jk是根據(jù)原子力測量對信息的獲取速率.在忽略其他環(huán)境影響下Γjk=γF,jk.描述 “信息獲取”導(dǎo)致的反作用的另外一個超算符,定義為

2.2 退相干速率

方程(4)中的退相干速率,難以在嚴格意義上得到,但可以通過定性的物理分析得到有明確物理意義的一些結(jié)果.例如,利用FGA[37],Schwartz等[33,38]通過計算不同勢能面上演化的凍結(jié)高斯波包之間的交疊,得到了如下的退相干速率公式:

式中,an是核波函數(shù)的有效寬度,其大小為

其中,ω是一個經(jīng)驗參數(shù),a0是玻爾半徑,Ekin是原子核的能量;Fi(j)(t)是?i(j)(R,r)態(tài)的Hellmann-Feynman力

進一步,Akimov等[31]對(7)式進行了簡化,得到了更加簡單和實用的退相干速率公式:

其中,p0是初始原子核動量大小,C是一個經(jīng)驗參數(shù).

作為另外一種不同的方案,Truhlar等[14,32]根據(jù)“退相干時間應(yīng)不小于最小的電子特征時間尺度”(基于能量-時間測不準(zhǔn)關(guān)系的考慮),唯像地提出了一個ED退相干速率公式:

其中,Ekin是原子核的能量,C=0.1 Hartree是一個經(jīng)驗參數(shù).值得注意的是,此退相干公式與原子核受力無關(guān),在勢能面平行區(qū)得到的退相干速率非零.

除了“ED”,原子在不同勢能面上的經(jīng)典運動,將感受到不同的“力”.基于對力的分辨,提出了力分辨(FD)退相干速率[15]:

式中,τc為核的特征運動時間,通常為電子運動特征時間的10?3倍;j(t)為原子處于勢能面Vj(R)感受到的粗?；骄?其具體形式將在下節(jié)中詳細討論.

2.3 量子路徑算法中原子核的受力形式

在混合經(jīng)典-量子動力學(xué)中,原子核由牛頓運動方程確定經(jīng)典運動規(guī)律.其中Erenfest平均場理論結(jié)合電子態(tài)演化的薛定諤方程,可以得到經(jīng)典力:

由于含時的|ψ(r;R)〉并不一定是電子哈密頓量的本征態(tài),式中第三個等號成立并不是由Hellmann-Feynman定理直接求得的結(jié)果[39],而是與Erenfest平均場理論的演化方程(2)有關(guān).此Erenfest力被許多勢能面跳躍理論采用作為原子力.但電子演化過程中含有隨機性質(zhì)的勢能跳躍,由于電子演化規(guī)律不遵循簡單的薛定諤方程,原子受力仍采用如上的形式就不再合適.

在量子路徑理論中,絕熱表象下決定原子核運動的經(jīng)典力為

此原子力由含有隨機項的電子態(tài)演化方程(4)推導(dǎo).此時原子力形式較為復(fù)雜,且其表達式中含有隨機項.與我們將原子核受力作為對電子態(tài)的量子測量輸出的觀點十分吻合.表明在電子時間尺度觀察原子,其感受到的 “瞬時力”同電子演化一樣也含有隨機性.(13)式中第二項Vj?ρjj=Vj˙ρjj/˙R已經(jīng)看出“原子瞬時力”包含了隨機漲落成分.

由于原子運動相對電子運動比較緩慢,原子在t時刻受到特征時間τ內(nèi)“粗?；钠骄Α薄j(t)作用,且此力不再包含隨機性.此外,在量子理論中也可以將 “原子時間尺度 (τ)”內(nèi)的平均勢能變化率作為原子受力,其大小為

廣泛被勢能面跳躍方法采用的Erenfest力(12)式和平均勢能變化確定的平均力(14)式,其區(qū)別不僅在于是否包含隨機性.而且在轉(zhuǎn)向點處可能會出現(xiàn)兩種力的受力方向相反的矛盾(附錄A).因此,在勢能面跳躍方法中原子受力仍采用Erenfest受力是不適宜的.

2.4 關(guān)于能量守恒方案

在勢能面跳躍理論中,由于原子核在不同勢能面之間跳躍,原子核勢能的變化并不連續(xù).因此,在跳躍后調(diào)節(jié)原子核的運動速度以保證原子核的總能量守恒:

式中,dk采用有效平均場非絕熱耦合矢量[33],

速度調(diào)節(jié)參數(shù)γ通過跳躍前后的動能差等于負的勢能差確定.如果發(fā)生原子核跳躍至某勢能面,而且此勢能面能量高于原子核初始總能量的情況,則認為此次跳躍是“禁閉”的.此外,根據(jù)經(jīng)典力學(xué)理論,此點是原子核的“反向點”應(yīng)將原子核的動量轉(zhuǎn)向,同時保持電子部分的量子態(tài)不變[40].

與勢能面跳躍理論不同,量子路徑理論中原子核的勢能變化是連續(xù)的.在理論上,采用(13)式或在dt→0時本質(zhì)上與其一致的(14)式作為原子力,可以保證能量自動守恒.但由于實際計算時采用的時間步長dt/=0,導(dǎo)致演化中每一步的能量計算誤差均為正數(shù),即能量 “正誤差”(詳見附錄B).因此我們針對量子路徑理論設(shè)計以下方案處理能量守恒問題.

第一,在模擬過程中,始終確保原子勢能小于系統(tǒng)的初始時刻的總能量即V(R)=Tr[Hel(R)ρ(R)]＜E0,E0可以定義為演化初始時刻的能量,倘若原子處于某一個單勢能面Vj(R),同樣要求Vj(R)≤E0.

第二,對于滿足V(R1)=E0的位置R1,由于其動能為0無法繼續(xù)向前運動,因此需重新定置R2,但新平均勢能VM(R2)=E0仍成立,方法如下:首先確保R2與R1有相同的電子波函數(shù);再去掉最低勢能面的分量;最后重新歸一化波函數(shù)計算新的平均勢能VM(R2).

第三,從R2處開始新的演化模擬,且保持電子部分不會發(fā)生任何改變,但原子的動能會改變?yōu)镋0?V(R2),其動量的方向與原來動量方向相反(即反射).

第四,在非轉(zhuǎn)折點,設(shè)定能量誤差上限δ,當(dāng)某時刻原子核總能量與初始總能量之間的差值大于δ時,則保持電子態(tài)不變并重新調(diào)整動量使原子核總能量滿足能量守恒.

其中,在反向點處的處理方法與文獻[12]處理方式以及Subotnik[41]提出的能量守恒方法有類似之處.

2.5 量子路徑算法的實施步驟

量子路徑方案不僅在理論概念上有清晰的物理意義,同時也是一種容易實現(xiàn)的分子動力學(xué)模擬算法.與Tully等[9]的勢能面跳躍算法相比,該方案并不需要人為設(shè)計復(fù)雜的算法使系統(tǒng)“跳躍”到某個勢能面上.現(xiàn)將量子路徑方案的算法概述如下.

1)將每一條軌跡中原子核的位置、動量以及電子密度矩陣初始化.

2)計算t時刻核構(gòu)型R(t)所對應(yīng)的哈密頓量(1)式的本征值和本征矢,得到絕熱基態(tài).每隔一個 “原子時間尺度 (τ=1000dt)”,計算原子的受力:.在同一個原子時間尺度τ內(nèi),原子核的速度保持不變,在t+τ時刻,原子核的速度變?yōu)?t+τ)=˙(t)+(t)τ.

3)每經(jīng)過“電子時間尺度”dt,選取合適的測量速率形式,計算Wiener增量ΔWjk=ξjk(t)dt.并求解相應(yīng)的量子路徑方程(4)式.

4)在同一個原子時間尺度τ內(nèi),原子核的速度保持不變,原子核的位置為R(t+dt)=R(t)+(t)dt.如果在t時刻,將進入下一個τ時段,則從t到t+τ時段內(nèi)速度為˙(t+τ)=˙(t)+(t)τ.

5)查看系統(tǒng)總能量,如果不需要進行動量調(diào)整則返回第二步繼續(xù)計算.如需要則按照上節(jié)的討論進行調(diào)節(jié).

3 對若干模型的數(shù)值模擬

本節(jié)分別使用FGA,ED和FD得到的退相干速率公式,利用量子路徑方案模擬五個模型,包括三個 Tully模型[9]和Subotnikhe和Shenvi[42]提出的另外兩個模型,并與嚴格的量子動力學(xué)結(jié)果做比較.對每一個模型,我們都將平均2000條隨機量子路徑.對于每條路徑,粒子都從勢能面的左端出發(fā),向右演化.原子質(zhì)量設(shè)為M=2000 a.u..

3.1 單交叉模型

在非絕熱表象中,該模型勢能面為[9]

其中,A=0.01,B=1.6,C=0.005,D=1.0.該模型的絕熱勢能面和耦合強度如圖1(a)所示,模擬結(jié)果見圖1(b)—(d).可以看出,有3個特殊的動量區(qū)域值得討論.第一,動量k＜4.5時,經(jīng)典的粒子幾乎全部返回到低勢能面上,這主要是因為經(jīng)典粒子的動能較小,無法提供其克服勢壘繼續(xù)向前演化的能量.第二,動量5＜k＜10時,處于低勢能面的粒子會隧穿勢壘繼續(xù)沿原方向演化,造成隧穿低勢能概率在k≈5時有較大幅度的梯度增加(圖1(b)所示);動量5＜k＜8時粒子動能較小無法躍遷到高勢能面上,所以反射概率為0;動量k≈8時,反射到低勢能面的概率略有增加(圖1(c)所示),因為在8＜k＜10時,粒子有足夠的動能向高勢能躍遷,但無足夠的動能使其隧穿高勢能面的勢阱,所以會以一定概率返回到低勢能面上.第三,動量k＞10時,在耦合區(qū)域低勢能面的粒子以一定的概率躍遷到高勢能面并沿原方向繼續(xù)演化.

圖1 (網(wǎng)刊彩色)單交叉模型 (a)絕熱勢能(實線)和非絕熱耦合強度(虛線)隨原子位置坐標(biāo)x的變化;(b),(c),(d)分別為隧穿至低勢能面、反射至低勢能面和隧穿至高勢能面的概率隨初始動量k的變化;利用FD退相干速率(藍色方塊)、ED退相干速率(綠色三角形)以及FGA退相干速率(紅色空心圓)通過量子路徑理論進行模擬,并與全量子動力學(xué)(Exact)(黑色實心點)模擬做了對比Fig.1. (color online)Single avoided crossing model:(a)Adiabatic potential energy level(solid line)and nonadiabatic coupling strength(dashed line)as a function of position x;panels(b),(c)and(d)are,respectively,initial momentum k effects on probabilities of transmission to the lower energy surface,ref l ection to the lower energy surface,and transmission to the upper energy surface.We present the quantum trajectory results by using the FD(blue square),ED(green triangle)and FGA(open red circles)decoherence rates,and compare them with the exact result from full quantum dynamics(Exact)simulation(black solid circle).

通過對三種退相干速率模擬的結(jié)果的比較,發(fā)現(xiàn)它們與全量子動力學(xué)的模擬結(jié)果幾乎相同.在幾個特殊動量區(qū)域也都符合,尤其是使用FD退相干速率和FGA退相干速率,與全量子動力學(xué)的模擬結(jié)果完全吻合.

3.2 雙交叉模型

雙交叉模型將導(dǎo)致量子干涉,并將造成隧穿概率的Stueckelberg振蕩,所以該模型更具挑戰(zhàn)性.該模型勢能面為[9]

其中,A=0.10,B=0.28,C=0.015,D=0.06和E=0.05.從圖2(a)可以看出,此模型由一個直線和兩個反轉(zhuǎn)的高斯曲線組成,并且有兩個連續(xù)的強耦合區(qū).

粒子從左端較低的勢能面開始,向右演化.當(dāng)進入第一強耦合區(qū)可能發(fā)生電子態(tài)的跳躍,使粒子離開第一強耦合區(qū)域時會以一定的概率繼續(xù)沿低勢能面向前演化,或發(fā)生量子跳躍以一定概率沿高勢能面向前演化,亦或以一定概率返回高低勢能面.若從波包演化的角度分析,在此演化過程中可能發(fā)生波包的分離與相遇.據(jù)圖2(b)—(d)全量子動力學(xué)的模擬結(jié)果顯示:當(dāng)粒子初始動能較低(lnE＜?3)時,絕大多數(shù)粒子會繼續(xù)沿低勢能面向前演化,少數(shù)部分粒子會返回到低勢能面.造成此現(xiàn)象的主要原因是在第一個耦合區(qū)部分粒子會向高勢能面躍遷.但是,處于高勢能面的粒子勢能增加,由于動能較低,無法使其隧穿繼續(xù)向前演化.當(dāng)粒子初始動能較大(lnE＞?3),粒子會隧穿到低勢能面或高勢能面上繼續(xù)向前演化.高低勢能面波包在第二耦合區(qū)會發(fā)生相干疊加,從而造成振蕩現(xiàn)象.隨著初始動能的增加,高低勢能面波包演化速率差也會較小,促使振蕩振幅變大.

圖2 (網(wǎng)刊彩色)雙交叉模型 (a)絕熱勢能(實線)和非絕熱耦合強度(虛線)隨原子位置坐標(biāo)x的變化;(b),(c),(d)分別為隧穿至低勢能面、反射至低勢能面和隧穿至高勢能面的概率隨初始能量E的變化;利用FD退相干速率(藍色方塊)、ED退相干速率(綠色三角形)以及FGA退相干速率(紅色空心圓)通過量子路徑理論進行模擬,并與全量子動力學(xué)(Exact)(黑色實心點)模擬做了對比Fig.2.(color online)Dual avoided crossing model:(a)Adiabatic potential energy(solid lines)and nonadiabatic coupling strength(dashed line)as a function of position x;panels(b),(c)and(d)are,respectively,initial energy E effects on probabilities of transmission to the lower energy surface,ref l ection to the lower energy surface,and transmission to the upper energy surface.We present the quantum trajectory results by using the FD(blue square),ED(green triangle)and FGA(open red circles)decoherence rates,and compare them with the exact result from full quantum(Exact)dynamics simulation(black solid circle).

對三種退相干速率模擬結(jié)果比較發(fā)現(xiàn),利用FD和FGA退相干速率模擬的數(shù)值結(jié)果與全量子動力學(xué)模擬結(jié)果符合良好.然而利用ED退相干速率進行的模擬在高動能區(qū)域無振蕩現(xiàn)象產(chǎn)生;在低初始動能區(qū)域也過高估測了返回的概率(大約高一個量級),與全量子動力學(xué)結(jié)果相差較大.

3.3 擴展耦合模型

擴展耦合模型是一個非常重要的模型,其勢能面形式為[9]

其中的系數(shù)選取為A=6×10?4,B=0.1和C=0.9.相應(yīng)的絕熱勢能面如圖3(a)所示,可以發(fā)現(xiàn)在此模型中包含了一個很寬的耦合區(qū).

當(dāng)粒子以低動量從低勢能面開始演化時,在強耦合區(qū)變?yōu)榛旌席B加態(tài),但由于能量守恒,處在高勢能面的波包由于動能無法提供其繼續(xù)沿原方向演化的能量,可能造成高勢能面的波包返回,但低勢能面的波包可以繼續(xù)向前演化,造成波包演化的退相干.

從圖3(b)—(d)可以看出:當(dāng)初始動量k＜28時,處在高勢能面的波包返回,而處在低勢能面的波包繼續(xù)隧穿,利用FD退相干速率與FGA退相干速率模擬的演化結(jié)果與全量子動力學(xué)模擬結(jié)果幾乎一致,然而利用ED退相干速率模擬的結(jié)果與全量子模擬結(jié)果差異較大,原因是過高地估計了返回到高低勢能面的概率;對于初始動量k＞28,處在高勢能面的粒子有足夠能量使其向前演化,所以在此動量區(qū),利用三種速率公式模擬的結(jié)果幾乎與全量子結(jié)果相似.

圖3 (網(wǎng)刊彩色)擴展耦合模型 (a)絕熱勢能(實線)和非絕熱耦合強度(虛線)隨原子位置坐標(biāo)x的變化;(b),(c)和(d)分別為隧穿至低勢能面、反射至低勢能面和反射至高勢能面的概率隨初始動量k的變化;FD退相干速率(藍色方塊)、ED退相干速率(綠色三角形)以及FGA退相干速率(紅色空心圓)通過量子路徑理論進行模擬,并與全量子動力學(xué)(Exact)(黑色實心點)模擬做了對比Fig.3.(color online)Dual avoided crossing model:(a)Adiabatic potential energy(solid lines)and nonadiabatic coupling strength(dashed line)as a function of position x;panels(b),(c)and(d)are,respectively,initial momentum k effects on probabilities of transmission to the lower energy surface,ref l ection to the lower energy surface,and ref l ection to the upper energy surface.We present the quantum trajectory results by using the FD(blue square),ED(green triangle)and FGA(open red circles)decoherence rates,and compare them with the exact result from full quantum(Exact)dynamics simulation(black solid circle).

3.4 啞鈴模型

啞鈴模型是一個對稱拓展耦合模型,其形狀類似于啞鈴.其勢能面和耦合強度,在x軸大于零部分(反應(yīng)結(jié)束部分)與擴展耦合模型的形狀一致,在x軸小于零部分(反應(yīng)開始階段)相當(dāng)于反轉(zhuǎn)的擴展耦合模型.其勢能面為[42]

A,B和C系數(shù)選取與模型(19)一致,Z=10.從圖4(a)可以看出,絕熱勢能面上有兩個強耦合區(qū).

在圖4(b)—(d)中,從全量子動力學(xué)演化結(jié)果可以看出:首先,在初始動量k＜28時,原子無法隧穿,幾乎全部反射到低勢能面上;其次,當(dāng)32＜k＜40時,部分原子隧穿至低勢面,部分反射到低勢能面,從波包演化的角度分析,當(dāng)初始動量在此間隔時,原子可以隧穿障礙到達高勢能面和低勢能面,但隧穿到高勢能面的原子在演化過程中勢能會不斷地增加,根據(jù)能量守恒,其無法提供繼續(xù)向前演化的動能,隧穿到高勢能面的原子會反射到低勢能面上;最后,當(dāng)k＞40時,原子完全隧穿至低勢能面或高勢能面上.

比較發(fā)現(xiàn),利用FGA退相干速率和FD退相干速率,通過量子路徑方法可以準(zhǔn)確地模擬系統(tǒng)的演化,尤其利用FGA退相干速率 (常數(shù)C=1.0)可以得到與全量子動力學(xué)模擬幾乎完全吻合的結(jié)果.但在初始動量32＜k＜40時,利用ED退相干速率模擬得到的結(jié)果與全量子動力學(xué)結(jié)果相差略大.

圖4 (網(wǎng)刊彩色)啞鈴模型 (a)絕熱勢能(實線)和非絕熱耦合強度(虛線)隨原子位置坐標(biāo)x的變化;(b),(c),(d)為隧穿、反射的概率隨初始動量k的變化(與圖3類似);利用ED退相干速率 (綠色虛線)、FD退相干速率(藍色實線)以及FGA退相干速率(紅色點劃線)分別通過量子路徑理論進行模擬,并與精確的量子動力學(xué)(Exact)(黑色虛線)模擬做了對比Fig.4.(color online)Dumbbell potential model:(a)Adiabatic potential energy(solid lines)and nonadiabatic coupling strength(dash line)as a function of position x;panels(b),(c)and(d)stand for initial momentum k effects on transmission and ref l ection probabilities as described in Fig.3.The quantum trajectory results by using FD(blue solid lines),ED(green dashed lines)and FGA(red chain line)decoherence rates are compared with the exact one from the full quantum(Exact)dynamics simulation(black dashed lines).

圖5 (網(wǎng)刊彩色)雙弓模型 (a)絕熱勢能(實線)和非絕熱耦合強度(虛線)隨原子位置坐標(biāo)x的變化;(b),(c),(d)為隧穿、反射的概率隨初始動量k的變化(與圖3類似);利用ED退相干速率(綠色虛線)、FD退相干速率(藍色實線)以及FGA退相干速率(紅色點劃線)分別通過量子路徑理論進行模擬,并與精確的量子動力學(xué)(Exact)(黑色虛線)模擬做了對比Fig.5.(color online)Double arch geometry:(a)Adiabatic potential energy(solid lines)and nonadiabatic coupling strength(dash line)as a function of position x;panels(b),(c)and(d)stand for initial momentum k effect on transmission and ref l ection probabilities as described in Fig.3.The quantum trajectory results by using FD(blue solid lines),ED(green dashed lines)and FGA(red chain line)decoherence rates are compared with the exact one from the full quantum(Exact)dynamics simulation(black dashed lines).

3.5 雙弓模型

雙弓模型在開始和反應(yīng)結(jié)束時,兩個勢能面基本重合,形式為[42]

A,B,C系數(shù)選取與模型(19)一致,Z=4.從圖5(a)可看出,該模型的非絕熱勢能面形狀像兩個弓.基于全量子動力學(xué)結(jié)果,此模型中有三個特殊的動能區(qū)域.第一,在低能量區(qū)域(k＜28),如圖5所示,僅僅只有處于低勢能面的波包進行了隧穿,而處于高勢能面的波包由于動能無法提供其繼續(xù)向前演化能量而被返回.造成其退相干率為100%.第二,在中能量區(qū)域(k剛大于28),波包可以隧穿到高勢能面或者低勢能面.然而低勢能面的波包由于能量守恒具有較大的速度,造成高低勢能面的波包在第二個非絕熱耦合區(qū)之前發(fā)生分離,所以勢能面的占據(jù)概率出現(xiàn)部分振蕩的現(xiàn)象.第三,在高能量區(qū)域(k?28),由于弓形勢壘所造成的速度差相對較小,上能級和下能級的波包發(fā)生了干涉現(xiàn)象,在此動量區(qū)域發(fā)現(xiàn)有較強的振蕩.

我們的模擬結(jié)果顯示:在初始具有較大動量時,利用FD退相干速率與ED退相干速率進行的模擬,最終兩個勢能面的占據(jù)概率幾乎是一致的,與全量子模擬的振蕩結(jié)果差別略大.利用FGA退相干速率的模擬結(jié)果,雖然在隧穿部分也有振蕩現(xiàn)象,但其振幅隨著初始動量的增加變化緩慢,與全量子動力學(xué)模擬也不完全符合.在初始動量較低時,利用三種速率公式模擬也都無法得到滿意的結(jié)果.原因是在“混合量子-經(jīng)典”方法中,核的運動由牛頓方程描述,難以對有復(fù)雜量子干涉現(xiàn)象的系統(tǒng)給出較理想的結(jié)果.

4 總 結(jié)

基于近期發(fā)展的經(jīng)典-量子混合模擬非絕熱分子動力學(xué)的量子路徑方法,本文對5個典型勢能面模型進行了模擬,包括單交叉模型、雙交叉模型、拓展耦合模型、啞鈴模型以及雙弓模型.在模擬過程中,我們對比檢驗了3個不同的退相干速率公式,分別為FGA,ED以及FD退相干速率.比較發(fā)現(xiàn),在單交叉模型中,三種退相干速率的模擬效果幾乎一致,均與全量子結(jié)果符合良好.對于雙交叉模型、拓展模型以及啞鈴模型,利用FGA與FD退相干速率,模擬的結(jié)果更好 (優(yōu)于ED退相干速率).在雙弓模型中,發(fā)現(xiàn)使用三種退相干速率,模擬結(jié)果與精確值差別較大,都不能很好地模擬精確的量子動力學(xué)過程.原因是“混合量子-經(jīng)典”方法中,核的運動由牛頓方程描述,難以納入復(fù)雜的量子干涉因素.如何發(fā)展更加有效的混合經(jīng)典-量子模擬方案,是未來研究的重要課題.

附錄A 瞬時力與Erenfest力“反向”問題

以一維勢能V1=?V2的Tully模型[9]一為例.在非耦合區(qū),定義V1為上能級勢能面,若粒子從位置R1點向R2點演化,通過平均勢能面梯度得到力的大小:

在任一R點處,其 Erenfest力大小則為(第二個等號是考慮非耦合區(qū),dji(R)近似為零的特例)

若從R1點演化到R2,電子態(tài)的密度矩陣不發(fā)生變化,二力可以近似相同.但考慮電子態(tài)由于測量投影可以發(fā)生如下變化:?V1始終是正值,在位置R1,R2處,粒子所處上能級概率分別為ρ1(R1),ρ1(R2),并且(2ρ1(R2)?1)V1(R2)?(2ρ1(R1)?1)V1(R1)＜0,2ρ1(R2)?1＞0.則

原子感受到的平均勢能面確實下降,所以根據(jù) “平均勢能的變化率”確定的力是正方向.但是,由于基本可以認為原子處于上勢能面,Erenfest力應(yīng)該是負方向.其原因在于量子路徑方法模擬過程中的∑jVj?Rρjj/= ∑j,kρkj(Vk?Vj)dkj,即Erenfest力不能描述電子態(tài)占據(jù)概率變化產(chǎn)生的力對系統(tǒng)勢能的改變.并且Vj?Rρjj=Vj˙ρjj/˙R,在速度較低的情況,勢能梯度產(chǎn)生的力?∑j(ρjj?RVj)小于電子態(tài)變化產(chǎn)生的力?∑j(Vj?Rρjj),在此情況下的Erenfest力忽略了起主導(dǎo)作用的電子態(tài)變化產(chǎn)生的力,導(dǎo)致了上述問題.

附錄B 借瞬時力(14)式實現(xiàn)能量守恒的誤差分析

采用新的原子受力的表達式,理論上可以自動實現(xiàn)能量守恒.在一維模型中,原子從t時刻,經(jīng)歷時間間隔Δt的過程中,速度保持不變,新的位置

t+Δt時刻的速度

定義時間間隔Δt內(nèi)的“平均加速度”

則容易推出t+Δt時刻的能量

在小量可以忽略的情況下,能量自動守恒

如果原子的初始動量較小,原子在跳躍至高勢能面后要轉(zhuǎn)向,發(fā)生“反射”現(xiàn)象.此時,原子速度要減小至接近0,然后轉(zhuǎn)向.在速度接近于0的過程中,(B6)式中的分母較小,因此能量中的增量仍可能具有一定的數(shù)值.并且1)即使這個增量較小,但由于是累加,因此總能量只會上升,無法減小;2)離開了耦合區(qū)后,只要投影沒有完畢,此增量就一直在累加,因此,總能量在實際計算時會有一定的能量增量誤差.

圖B1 (網(wǎng)刊彩色)通過量子路徑方法,分別模擬了使用能量守恒(CE)方案(綠色實線)以及未使用能量守恒(NCE)方案(藍色虛線)兩種情況總能量隨時間的變化Fig.B1.(color online)Energy variations along the time.The results obtained by using(blue dashed lines)and not-using(green solid lines)the rule of energy conservation.為了更清楚地反映此現(xiàn)象,我們以拓展耦合模型為例,用量子路徑方法模擬了單條路徑總能量隨演化時間的變化(圖B1).設(shè)置初始動量k0=8.1578 a.u.(有反射現(xiàn)象產(chǎn)生).從圖B1我們發(fā)現(xiàn)其模擬結(jié)果與以上理論分析完全吻合.

[1]Gerber R B,Buch V,Ratner M A 1982J.Chem.Phys.77 3022

[2]Micha D A 1983J.Chem.Phys.78 7138

[3]Li X S,Tully J C,Schlegel H B,Frisch M J 2005J.Chem.Phys.123 084106

[4]Tully J C,Preston P K 1971J.Chem.Phys.55 562

[5]Miller W H,George T F 1972J.Chem.Phys.56 5637

[6]Kuntz P J,Kendrick J,Whitton W N 1979Chem.Phys.38 147

[7]Blais N C,Truhlar D G 1983J.Chem.Phys.79 1334

[8]Ali D P,Miller W H 1983J.Chem.Phys.78 6640

[9]Tully J C 1990J.Chem.Phys.93 1061

[10]Kuntz P J 1991J.Chem.Phys.95 141

[11]Webster F,Wang E T,Rossky P J,Friesner R A 1994J.Chem.Phys.100 4835

[12]Prezhdo O V,Rossky P J 1997J.Chem.Phys.107 825

[13]Zhu C Y,Jasper A W,Truhlar D G 2004J.Chem.Phys.120 5543

[14]Zhu C Y,Nangia S,Jasper A W,Truhlar D G 2004J.Chem.Phys.121 7658

[15]Feng W,Xu L T,Li X Q,Fang W H,Yan Y J 2014AIP Adv.4 077131

[16]Li B,Han K L 2009J.Phys.Chem.A113 10189

[17]Li B,Chu T S,Han K L 2010J.Comput.Chem.31 362

[18]Yang M H,Huo C Y,Li A Y,Lei Y B,Yu L,Zhu C Y 2017Phys.Chem.Chem.Phys.19 12185

[19]Lu J F,Zhou Z N 2016J.Chem.Phys.145 124109

[20]Schubert A,Falvo C,Meier C 2016J.Chem.Phys.145 054108

[21]Wang L J,Prezhdo O V,Beljonne D 2015Phys.Chem.Chem.Phys.17 12395

[22]KosloffR 1988J.Phys.Chem.92 2087

[23]Schatz G C 1996J.Phys.Chem.100 12839

[24]Zhang J Z H,Dai J,Zhu W 1997J.Phys.Chem.A101 2746

[25]Guo H,Yarkony D R 2016Phys.Chem.Chem.Phys.18 26335

[26]Chu T S,Zhang Y,Han K L 2006Int.Rev.Phys.Chem.25 201

[27]Chu T S,Han K L 2008Phys.Chem.Chem.Phys.10 2431

[28]Zhang S B,Wu Y,Wang J G 2016J.Chem.Phys.145 224306

[29]Jacobs K,Steck D A 2006Contemp.Phys.47 279

[30]Xie B B,Liu L H,Cui G L,Fang W H,Cao J,Feng W,Li X Q 2015J.Chem.Phys.143 194107

[31]Akimov A V,Long R,Prezhdo O V 2014J.Chem.Phys.140 194107

[32]Zhu C Y,Jasper A W,Truhlar D G 2005J.Chem.Theory Comput.1 527

[33]Bedard-Hearn M J,Larsen R E,Schwartz B J 2005J.Chem.Phys.123 234106

[34]Prezhdo O V 1999J.Chem.Phys.111 8366

[35]Granucci G,Persico M 2007J.Chem.Phys.126 134114

[36]Thachuk M,Ivanov M Y,Wardlaw D M 1998J.Chem.Phys.109 5747

[37]Heller E J 1981J.Chem.Phys.75 2923

[38]Schwartz B J,Bittner E R,Prezhdo O V,Rossky P J 1996J.Chem.Phys.104 5942

[39]Lan Z G,Shao J S 2012Prog.Chem.24 1105(in Chinese)[蘭崢崗,邵久書2012化學(xué)進展 24 1105]

[40]Hammes-Schiffer S,Tully J C 1994J.Chem.Phys.101 4657

[41]Subotnik J E 2010J.Chem.Phys.132 134112

[42]Subotnik J E,Shenvi N 2011J.Chem.Phys.134 024105