地質流體中絡合物同位素平衡分餾系數的第一性原理計算*

孫明光，王佳新，袁順達

(1中國地質科學院礦產資源研究所自然資源部成礦作用與資源評價重點實驗室，北京 100037；2北京大學造山帶與地殼演化教育部重點實驗室，北京 100871)

流體的溫度、壓力、相態或絡合物的價態及成鍵環境等的變化可以導致具有多種同位素的元素在遷移、富集和沉淀成礦過程中發生顯著的同位素分餾，這些同位素分餾可以作為示蹤流體地球化學過程的重要手段（Kelly et al.,1979;Sheppard.,1981;Gallagher et al.,1992;Rye,1993;Bethke et al.,2005;Dreher et al.,2008;Wen et al.,2015;Ohmoto,2018;Gao et al.,2018），尤其在示蹤成礦物質來源、解釋成礦元素同位素分餾機制及進一步確定礦床成因和建立成礦過程中相關同位素組成變化的模型等方面具有重要意義：如錫同位素分餾可以記錄熱液流體的氧化還原反應，從而可以用來解釋含錫成礦流體的物理化學條件變化（Yao et al.,2018）；在熱液體系中銻同位素分餾可以用來指示含銻成礦流體的流動方向，從而指示銻礦的勘查等（Zhai et al.,2021）。同位素質譜分析技術的革命性進展，尤其是多接收器等離子體質譜儀(MC-ICP-MS)的出現，推動了非傳統穩定同位素尤其是金屬穩定同位素地球化學的快速發展，并在地質研究中得到了廣泛的應用。目前，鋰、鐵、鎂、鈣、鈦、釩、鎘等非傳統穩定同位素在地表風化、表生環境、熱液成礦等各類地質過程研究中發揮了越來越重要的作用，并取得了一系列重要進展（朱祥坤等，2013）。其中，對不同絡合物之間的穩定同位素平衡分餾系數的測定是開展穩定同位素地球化學研究的必要前提。然而，一些元素由于易揮發、易氧化、相對質量差值小以及溶液中絡合物難以分離測定等原因，導致實驗方法無法測得溶液中一些元素不同絡合物之間的平衡同位素分餾系數(Schauble,2004;Seo et al.,2007)，從而阻礙定量研究這些元素的同位素分餾機制。另外，在合理的時間尺度范圍內，一些礦物-溶液之間的同位素交換平衡分餾系數難以通過實驗手段獲得，在一定程度上限制了這類同位素地球化學的推廣應用（Schauble,2004）。

研究表明，基于密度泛函理論（DFT:density functional theory）的第一性原理計算是獲得各類元素在不同絡合物或礦物之間同位素分餾系數的一種重要手段（Cygan，2001）。目前，應用第一性原理計算溶液中絡合物的約化配分函數比（β:reduced isotope partition function ratio）的模型主要有2種：顯式溶劑模型（explicit solvation models）和隱式溶劑模型（implicit solvent model)。通過采用這2種模型，前人已經成功獲得了一系列元素在不同物種之間穩定同位素的平衡分餾系數，如O、H、Cl、Fe、Mg、Si、Cu、Ag、Cd、Ge、Se、B、Br等（Seo et al.,2007;Méheut et al.,2007;2009;2010;Domagal-Goldman et al.,2008;Pinilla et al.,2015;Yang et al.,2015;Gao et al.,2018;Wang et al.,2019;Fujii et al.,2010;2013;2014;2018;Gao et al.,2021），在很大程度上推動了同位素地球化學的發展。然而，采用上述2種模型對溶液中的絡合物進行模擬計算時，由于介質環境及絡合物構型等方面存在明顯不同，導致計算的結果往往會存在一定差異。本文對比研究了前人應用這2種模型對不同同位素體系的計算結果，進一步分析了這2種計算穩定同位素分餾系數方法的優、缺點，為今后開展第一性原理計算流體中絡合物平衡同位素分餾時，為合理選擇不同模型提供科學依據。

1 同位素分餾理論和計算方法

1.1 同位素分餾理論

同一元素不同同位素間質量差可以導致在各類地質過程中發生同位素分餾，因而在同位素交換反應達到平衡時，由于生成物和反應物的能量差別而發生同位素的平衡分餾。這種能量差極小，主要是由振動能引起的，因而可以依賴振動配分函數來計算2種絡合物之間的同位素平衡分餾系數(Urey et al.,1933;Bigeleisen et al.,1947;Urey,1947)。2種物質A和B中的同位素交換反應：

A和A*分別代表的是含有同一元素的輕、重同位素的物質。平衡常數KA-B可以通過2種物質的約化配分函數比（β）的比值獲得

反應式（1）中涉及到的是1個原子(n=1)的交換反應，因此分餾系數α可以替代K（α=K1/n），α=K。

β可以通過公式（3）獲得：

其中

h和k分別是普朗克常數和玻爾茲曼常數，T是溫度（Kelvin），υi表示原子的振動頻率。通過公式（3）和公式（4），可得到只與振動頻率υi有關的RPFR(β)。因此如果能夠獲得體系的振動頻率，便可獲得β及同位素平衡分餾系數。計算公式的詳細推導及解釋見Chacko等(2001)。

1.2 β值的計算方法

1.2.1 隱式溶劑模型

隱式溶劑模型是通過極化連續介質模型(PCM:Polarizable Continuum Models)的方法模擬溶液環境（Cances et al.,1997），該方法是將溶質置于一個被溶劑包圍的反應場空腔中，且空腔外部的溶劑被定義為具有均勻介電常數的連續介質（widerek et al.,2019;Wilson et al.,2015;Li et al.,2009），圖1為隱式溶劑模型，[Sn4+Cl4]（圖1a）和[Sn2+Cl3]-（圖1b）被置于空腔中，這種空腔簡化了水溶液的動態特性，因而計算耗時降低。例如，Sun等（2022）應用

圖1 優化后隱式溶劑模型的絡合物分子結構a.[Sn4+Cl4];b.[Sn2+Cl3]-（據Sun et al.,2022）Fig.1 Complexes molecular structure of optimized implicit solvent model a.[Sn4+Cl4];b.[Sn2+Cl3]-（after Sun et al.,2022）

B3LYP(Becke-style 3-parameter density functional theory with the Lee-Yang-Parr correlation functional)(Lee et al.,1988;Becke,1993a;1993b)計算含錫絡合物的分子結構和振動頻率，使用混合基組對[Sn4+Cl4]和[Sn2+Cl3]-進行結構優化和頻率計算,其中對Cl原子使用6-311+G(d,p)基組，Sn原子使用SDD(Stuttgart-Dresden relativistic effective core potential)(Schwerdtfeger et al.,1989），上述計算均基于Gaussian 16軟件（Frisch et al.,2016）。通過Gaussian View可視化軟件，分別搭建[Sn4+Cl4]和[Sn2+Cl3]-絡合物合理的分子結構，然后依次選擇上述的計算基組提交計算任務（依次計算124Sn、116Sn），計算完成后從輸出文件中檢查振動頻率是否有虛頻，若無虛頻則可通過公式（3），求得β值；若有虛頻表明此時的分子結構不穩定，體系能量不是最低，因而需要進一步調整分子結構重新提交計算任務，直至無虛頻為止。

1.2.2 顯式溶劑模型

為精確處理溶液中水分子的相互作用、水的動態特征和氫鍵對β值的影響，前人通過在溶質周圍添加水分子的方法對溶液環境進行模擬，故稱為顯式溶劑模型，此方法能夠很好地描述水溶液的動態特征（Liu et al.,2005;Li et al.,2009;2011;Wu et al.,2015;Gao et al.,2018;Wang et al.,2019;Gao et al.,2021）。圖2a、b為基于顯式溶劑模型對[Sn4+Cl4]和[Sn2+Cl3]-的結構優化圖，圖中可以顯示[Sn4+Cl4]和[Sn2+Cl3]-被周圍的水分子包裹。Sun等（2022）應用第一性原理分子動力學方法（first-principles molecular dynamics(FPMD)）對[Sn4+Cl4]和[Sn2+Cl3]-絡合物進行模擬計算。首先，對2個含錫絡合物在溶液中的動態特征進行模擬，計算采用NVT熱動力學系綜（NVT thermodynamic ensemble）(即在模擬計算過程中原子數(n),體積(V)和溫度(T)固定不變，壓力(P)和 能量(E)允許改變)。面波能量截取設置為600 eV，布里淵區k點 為1×1×1。已有研究表明，在模擬計算過程中50個水分子就能夠真實反映溶質所處的溶液環境，且獲得精確的β值(Wang et al.,2019;Gao et al.,2021)，因而筆者的計算也在[Sn4+Cl4]和[Sn2+Cl3]-溶質周圍添加了50個水分子。為保持整體電荷平衡，在[Sn2+Cl3]-溶液中添加了1個Na離子。參數設置確保了靜態壓力接近0，密度設置為水的密度1 g/cm3。模擬2個絡合物在溶液中的動力學特征的時間步長為1 fs，模擬的總時間是50 ps（50 000 fs），以獲得足夠的構型去計算每個絡合物的β值。應用VASP軟件能夠較簡便的建立團簇模型，可以通過設置時間步長(1 fs)和模擬計算的總時間（ps）來獲得不同的構型，比Gaussian軟件的處理過程相對簡單，但計算耗時較長，比如，通過改變[Sn4+Cl4]和[Sn2+Cl3]-絡合物的輸入文件，分別獲得[Sn4+Cl4]和[Sn2+Cl3]-的50 000個構型，體系的溫度和能量變化（圖3a～d）趨勢表明，在30 ps至50 ps之間體系能量變得穩定，因此在30 000至50 000構型之間每隔2 ps(2000構型)取1個構型，對獲得的10個不同構型結構優化，使用優化后的構型進行頻率計算，最終通過公式(3)得出每個絡合物的β平均值。

圖2 優化后顯式溶劑模型的絡合物分子結構a.[Sn4+Cl4](H2O)50(擴胞：2×2×2);b.[Sn2+Cl3]-(H2O)50(擴胞：2×2×2)綠色表示Cl原子；紅色表示O原子；白色表示H原子；紫色表示Sn原子；黃色表示Na原子；虛線為分子動力學模擬完成后H2O分子之間形成的氫鍵（據Sun et al.,2022）Fig.2 Complexes molecular structure of the optimized explicit solvent model a.[Sn4+Cl4](H2O)50(super cell：2×2×2);b.[Sn2+Cl3]-(H2O)50(super cell：2×2×2)Colorsfor atoms:chlorine,green;hydrogen,white;tin,purple;sodium,yellow;thedotted lineshowsthehydrogen bond formed between H2O molecules after molecular dynamics simulation(after Sun et al.,2022)

圖3 [Sn4+Cl4]和[Sn2+Cl3]-分子動力學(FPMD)模擬過程中溫度和能量隨時間的變化趨勢Fig.3 Variation trend of temperature and energy with time of[Sn4+Cl4]and[Sn2+Cl3]-complexes during molecular dynamics(FPMD)simulation

2 隱式溶劑模型與顯式溶劑模型對比

為了對比2種模型對絡合物平衡同位素分餾系數計算結果差異，本研究系統收集了前人通過隱式或顯式溶劑模型計算的B、Se和Sn同位素計算結果（表1）。從0～350℃，以隱式溶劑模型計算的[Sn2+Cl3]-的β值范圍為1.002 084～1.000 406，[Sn4+Cl4]的β值范圍為1.005 445～1.001 082，以顯式溶劑模型計算的[Sn2+Cl3]-和[Sn4+Cl4]的β值范圍分別為1.002 385～1.000 473和1.006 877～1.001 394（Sun et al.,2022），且在對Sn2+Cl3-絡合物執行分子動力學過程中，有1個水分子與Sn結合，此時的Sn2+Cl3-絡合物的配位數為4，配位元素為Cl和O，而隱式溶劑模型對此無法模擬。根據計算結果表明，在恒定溫度條件下同一種模型計算的β值的變化趨勢為[Sn4+Cl4]＞[Sn2+Cl3]-，這與前人得出高價態同位素富集重同位素的結論相一致(She et al.,2020;Wang et al.,2021),且隨著溫度的增加[Sn4+Cl4]和[Sn2+Cl3]-的β值都逐漸降低（圖4a、b）。雖然這種變化趨勢是一樣的，但是2種模型計算的β值有明顯的差異，如在25℃條件下，[Sn2+Cl3]-在隱式和顯式溶劑模型下的β值分別為1.001 755和1.002 014，[Sn4+Cl4]的β值分別為1.0046和1.005 829(表1)。此外，前人采用隱式溶解模型和顯式溶劑模型對同一絡合物β值的計算結果也存在一定差異：如在25℃，對流體中B（硼）絡合物B(OH)3，以隱式溶劑計算的β值為1.2313，以顯式溶劑模型（B(OH)3-(H2O)34）計算的β值為1.2267（表1）（Liu et al.,2005）；同樣，在25℃時，流體中Se絡合物以隱式溶劑模型計算的β值為1.0350，而顯式溶劑模型計算的β值則為1.0243（Li et al.,2011）。下面從以下3個方面對比研究這2種模型的主要差異。

圖4 兩種模型計算的Sn2+Cl3-和Sn4+Cl4絡合物的1000lnβ值（a）和1000lnα值（b）Fig.4 1000lnβ(a)and 1000lnα(b)values of Sn2+Cl3-and Sn4+Cl4 complexes calculated by two models

表1 2種模型計算的不同元素的β值(0～350℃)Table 1 βvalues for different elements calculated by the two models(0～350℃)

2.1 兩種模型計算β值耗時對比

隱式溶劑模型計算絡合物的β值并不增加溶質周圍的原子數目（圖1a、b）。已有研究表明，通常情況下，無機絡合物的單體溶質原子的數目一般較少（1～20），如HSe-（2）、Fe(H2O)62+（19）、Fe(H2O)63+（19）、Ge(OH)4（9）和V4+O(H2O)52+（17）等（Jarzecki et al.,2004;Li et al.,2009;2011;Wu et al.,2015），對隱式溶劑模型下的絡合物分子結構優化和頻率計算的過程中，若每個絡合物分子分配32個核數進行計算，大約耗時數個小時即可計算完成，且計算完成后輸出的結果文件大小僅為數兆（M）。而顯式溶劑模型需要在絡合物周圍添加不同數量的水分子,計算的原子數目增加，如B(表1)，Sn(圖2)，Ge和Se（表2），每個絡合物計算的原子數目范圍為20～154（表2）。而且由于絡合物構型對β值有影響，使得每個絡合物必須計算不同的構型（2～10個），因而在使用與隱式模型計算相同的基組水平對顯式溶劑模型下的絡合物進行計算時（Gaussian軟件），耗時比隱式溶劑模型所需更久，根據絡合物數量可能需要數月，計算完成后輸出的結果文件為數十兆至數百兆，如VASP計算的Sn2+Cl3--(H2O)50和Sn4+Cl4-(H2O)50的結果文件可以達到2 G～3 G。

2.2 兩種模型計算時水分子數和構型對β值的影響

隱式溶劑模型的方法計算的絡合物只有一種構型，因而不能評估不同構型(configurations)對絡合物β值的影響，如前人對B、Se和Sn絡合物β值的計算即只有一個構型（表1）。相對而言，采用顯式溶劑模型的方法對絡合物β值計算時，可以建立多種初始構型（如表2中A、B、C、D為4種不同的構型），且不同構型其β值也有一定差異（表2），表明構型是影響β值的因素之一，因而在計算流體中絡合物的平衡分餾系數需考慮構型的影響。如前人計算B、Se、Sn和Ge絡合物的β值，分別通過建立多種構型來計算β值，并最終采用不同構型獲得β值的平均值作為絡合物的β值，從而降低構型對計算結果的影響（表2）（Liu et al.,2005;Li et al.,2009;2011;Wu et al.,2015）。

表2 不同水分子數和不同構型的Ge、Se和B絡合物的β值（A、B、C和D表示不同的構型）Table 2 βvalues of Ge,Se and B complexes with different numbersof water molecules and different configurations(A,B,C and D represent different configurations)

隱式溶劑模型對溶質周圍水環境的模擬方法較為簡單，因而該方法不能系統評估水分子數的變化對計算結果的影響。顯式溶劑模型則可通過在溶質周圍添加一系列的水分子建立不同的團簇模型的方法（如6、12、18、24和30）來評估水分子數量差異對β值的影響（表2）：例如，Liu等(2005)從計算B(硼)最小的構型開始，逐漸增加溶質周圍的水分子數量（6、12、8、24和30），計算的結果表明，隨著水分子數的增加，各構型之間的β值差異逐漸降低，如30個水分子模型B(OH)3-(H2O)30_C和34個水分子模型B(OH)3-(H2O)34之間的β值相差為0.0004（表2），這種差異對于目前自然體系B同位素分餾水平可忽略不計，因此，30或34個水分子足以獲得B絡合物穩定的β值，無需進一步計算更多水分子的硼絡合物的β值。同樣，Li等（2009）建立了一系列的模型研究Ge絡合物的β值，在每次優化“水滴團簇模型”時確保了最小的模型（6個H2O分子）的能量為最低，結構最穩定，在上一個模型的基礎上增加6個水分子直至30，隨著水分子數的增加，計算獲得的β值變化逐漸減小。Li等(2011)計算Se絡合物之間的β值，也通過建立了一系列的模型（6、12、18和24）計算Se絡合物的平均β值（表2）?；谕瑯拥姆椒?，前人也成功預測了V（釩）絡合物的平衡同位素分餾系數（Wu et al.,2015）。因此，在計算流體中絡合物同位素分餾時絡合物的構型和水分子數量對其β值的計算產生一定影響，相比于隱式溶劑模型，顯式溶解模型在計算絡合物β值時更能充分的評估絡合物構型及水分子數對計算結果的影響。

2.3 兩種模型模擬的水溶液環境對β值的影響

根據前人應用隱式溶劑模型和顯式溶劑模型對B，Se和Sn絡合物的β值計算結果表明，不同模型計算絡合物的β值有明顯的差異（表1、表2），如B(OH)3絡合物的β值分別為1.2313（B(OH)3）和1.2267（B(OH)3-(H2O)34），SnCl3-的β值分別為1.001 755（SnCl3-）和1.002 014（SnCl3--(H2O)50），SnCl4的β值分別為1.0046（SnCl4）和1.005 829（SnCl4-(H2O)50），的β值分別為和1.038 29造成這種差異的原因是兩種模型對水環境的模擬方式不同。

隱式溶劑模型雖然可通過應用PCM方法對溶質周圍的水環境進行模擬，但該模型把溶質周圍的水定義成具有均勻介電常數的連續介質，處理水的動態特性及水分子之間形成的氫鍵精確度稍低（Liu et al.,2005;Li et al.,2011;2009）。盡管可通過Gaussian軟件中的多種PCM(UAHF,UA0,CPCM和SCIPCM)模型，可以提高計算過程中對溶液水環境中氫鍵模擬的精確度，已有學者通過PCM UA0模型計算獲的B(OH)3的振動頻率值（885 cm-1）與實驗測定值(880 cm-1)具有較好的一致性（Liu et al.,2005），但計算的β值與實驗測定值仍存在明顯差異（Klochko et al.,2006;Byrne et al.,2006;Rustad et al.,2007）。為了更好的模擬溶液的動態特征，并合理表達溶質與周圍水分子之間的相互作用，Liu等(2005)采用顯式溶劑模型（水滴團簇模型：water-droplet），在硼絡合物周圍添加了不同數量的H2O分子來模擬溶質周圍的水環境（表1）。雖然這種模型對計算的要求高，但是只要團簇模型足夠大，計算的結果精度就越高（Liu et al.,2005）?；谠擄@式溶劑模型，能夠使溶質周圍形成至少2～3層水分子，且每個水分子在自身的位置形成兩個或更多的氫鍵，最終形成由4、5和6個水分子組成的多個環狀團簇，從而在溶質周圍的水分子之間形成了更真實的氫鍵環境。圖2展示出在計算完成后，水分子由于相互作用形成的氫鍵（圖中虛線所示）及團簇?；陲@式溶劑模型，Liu等(2005)成功預測了B(OH)3和B(OH)4-之間的同位素分餾系數，并且被后續的實驗和計算結果所證實(Klochko et al.,2006;Byrne et al.,2006;Rustad et al.,2007)。

除此之外，Li等(2009)通過PCM方法和顯式溶劑模型方法計算同一Ge絡合物的1000lnβ值，計算結果表明雖然2種模型計算出的Ge-O鍵長幾乎一致，但是1000lnβ值仍有明顯的差異，尤其是對部分高負價態的絡合物，差異更加明顯，而且隱式溶劑模型（PCM方法）計算時假設絡合物中Ge-O鍵長相等的結構是高度對稱的，即使在輸入文件中寫入關鍵詞“nosymm”(not assume any symmertry),但是從輸出的結果可以看出分子結構仍然是對稱的。然而，分子結構的高度對稱是不真實的，因為溶質和水分子外層形成的氫鍵在不同的方向結構不可能是相同的。因此，雖然隱式溶劑模型對絡合物分子結構的優化能夠得到與顯式溶劑模型一致的鍵長，但是由于隱式溶劑模型對水環境的模擬沒有顯式溶劑模型精確，仍然導致2種模型計算的結果不一致(Li et al.,2009)。

綜上所述，應用隱式溶劑模型計算溶液中絡合物的β值有一定的局限性，雖然計算量相對小，能明顯節省計算時間，但是不能精確處理溶液中水分子之間的相互作用、強氫鍵，也不能有效評估構型對β的影響。而顯式溶劑模型通過在溶質周圍添加水分子的方法，可以更好的模擬真實溶液環境，而且隨著水分子數量的增加，計算的結果越精確，隨著水分子數增加，能夠逐漸真實反映溶質所處的溶液環境，并獲得精確的β值（Wang et al.,2019;Gao et al.,2021）。由于溶液的動態特征導致絡合物存在各種構型，為了降低構型對計算結果的影響，前人在保持溶質周圍水分子數量相同的條件下，對每個絡合物建立不同的構型，最終計算獲得β值的平均值。因此應用顯式溶劑模型計算水溶液中絡合物的同位素分餾系數時需要考慮溶質周圍添加的H2O分子數量和不同構型對計算結果的影響，從而盡可能降低計算誤差。如圖4a所示，雖然顯式溶劑模型和隱式溶劑模型計算的Sn4+Cl4絡合物的1000lnβ值相差不大，但是Sn2+Cl3-的1000lnβ有較大的差異，而且Sn4+Cl4和Sn2+Cl3-之間的平衡錫同位素分餾系數1000lnα（圖4b）也有較大差異，這主要是由于隱式模型忽略了水的構型和水分子與溶質之間相互作用對β值的影響，而顯式模型則充分模擬了構型和水分子之間的相互作用。

3 結論

顯式溶劑模型和隱式溶劑模型是計算溶液中各絡合物之間的平衡同位素分餾系數的2種主要模型，但是由于這2種模型對溶液模擬的方式不同，導致計算耗時及β值有差異，具體表現為：①隱式溶劑模型計算的原子數量較少、計算的速度快、耗時少。顯式溶劑模型涉及的水分子數量和構型多，因此需要計算的原子數量多且需要反復計算不同的構型，導致耗時久；②溶質周圍水分子數和構型對β值的影響：隱式溶劑模型對水溶液環境模擬的較為簡單，導致絡合物只有一個構型；顯式溶劑模型通過建立一系列的水分子數和構型評估它們對計算結果的影響，最終通過求各構型的平均值降低這些因素引起的誤差；③隱式溶劑模型把溶質周圍的水定義成具有均勻介電常數的連續介質，這不能很好的處理水的動態特性及水分子之間形成的氫鍵，導致計算的結果與實驗值偏差較大；顯式溶劑模型能更好的模擬溶液的動態特征，且表達溶質與周圍水分子之間的相互作用，能獲得更加精確的β值。因此，目前已有研究表明：顯式溶劑模型是計算溶液中絡合物之間平衡同位素分餾的更為有效的手段。

致謝感謝審稿人百忙之中抽出時間審閱本文，為本文提出的建設性意見，從而使得文章語言表達通暢，邏輯更加清楚，文章有了質的提升。