镎(IV)、镅(III)、鋦(III)的AMBER力場參數化及評估

劉子義，夏苗仁，柴之芳，王東琪

1中國科學院高能物理研究所多學科中心, 北京 100049

2中國科學院大學化學學院, 北京 100049

3蘇州大學放射醫學與輻射防護國家重點實驗室, 蘇州大學放射醫學與防護學院, 江蘇 蘇州 215123

1 引言

近年來，隨著計算能力的迅速提高，使用計算機模擬化學微觀體系并探究作用機制的方法逐漸成為化學研究的主要手段之一。在此方法中，合理地描述粒子間的相互作用是確保模擬結果可靠的關鍵。為了改進和完善對粒子間相互作用的描述，在過去的幾十年里，發展了多種理論模型，主要可分為兩類：(a)基于電子運動規律的量子力學描述；(b)基于分子力場的經典力學描述。其中，量子力學方法通過求解薛定諤方程獲得電子的運動規律，進而得出分子體系的性質，與實驗值有較好的吻合，但由于算法等限制，計算量大，以目前的計算能力難以處理復雜體系。分子力場的描述忽略了電子自由度，用經典力學的經驗方程取代迭代算法，將原子作為點模型處理，有效地提高了計算效率，可以處理數百萬原子的體系。在分子力場的描述中，粒子間的相互作用分解為鍵相互作用和非鍵相互作用。其中，鍵相互作用包括鍵伸縮、角彎曲、二面角扭轉等，非鍵相互作用包括庫倫作用和范德華作用。根據模擬體系和科學目標的側重點不同，相互作用的函數表達式不同，發展出不同版本的力場。目前常用的力場有AMBER1、GROMOS2及CHARMM3等。

研究錒系元素(An)在水溶液中的行為對于理解錒系元素的遷移、放射性廢物處置及核燃料的分離過程有著至關重要的作用。特別是針對不斷增加的核燃料廢料的科學挑戰，使用離子交換、溶劑萃取、電沉積等方法從混合物中分離出放射性核素(如U、Np、Pu、Am等)，都需要深入了解錒系元素在水溶液中的化學行為。通常，在水溶液及生物體系中，Th、U、Np、Pu等元素以四價單原子離子形式存在，其中U還能以鈾酰離子(UO22+)形式存在，而Am和Cm一般以三價單原子離子形式存在。由于錒系元素具有高放射性和化學毒性，在實驗上研究這些金屬離子的水溶液化學行為較為困難。而理論化學中的量子力學描述方法由于昂貴的計算成本，只能關注這些元素周圍第一配位層水分子的配位行為，同時，也很難提供動力學的信息。基于分子力場的模擬是一種可替代的方案4，而錒系金屬的力場參數是描述錒系溶液行為的關鍵。

由于錒系元素電子結構的復雜性，如高角動量、多組態、極化效應、相對論效應、旋軌耦合等，導致錒系元素具有復雜的化學行為，因此，相比于碳、氫、氧、氮等主族元素，錒系元素力場參數的擬合更具有挑戰性，需要考慮多方面的影響，如錒系元素的價態、種態和化學環境等5。從上個世紀90年代以來，多個課題組致力于錒系元素特定氧化態下力場參數的開發，包括有Wipff課題組報道的UO22+在水溶液中的AMBER力場參數6,7；Maginn課題組根據量子化學計算出的勢能面擬合出五價、六價錒酰離子(AnO2+/2+)的AMBER力場參數8；Merz等在AMBER力場框架下根據An―Ow的鍵長和An的溶劑化自由能擬合出了Th4+、U4+、Pu4+的12-6及12-6-4勢的力場參數9,10等。Am3+和Cm3+在極化力場中有相應的參數報道，如Gagliardi等11報道的Cm3+的NEMO力場參數、Spezia等12報道的Am3+、Cm3+的類Car-Parrinello力場參數等。這些參數能夠較好地重現實驗結果，但受限于較大的計算量、力場參數較少且各版本極化力場之間的不可移植性等因素，制約了極化力場的應用，有關錒系溶液動力學在力場水平上的研究仍以非極化力場為主。另外，目前Np4+、Am3+、Cm3+的非極化力場尚未報道。

本項工作旨在補充對于錒系關鍵離子Np4+、Am3+、Cm3+的AMBER力場參數，通過對比文獻報道的相關理論與實驗結果，評估參數的可靠性，進而結合文獻報道10的Th4+、U4+、Pu4+參數，系統研究這幾種重要的錒系離子在水溶液中的動力學。

2 研究方法

2.1 參數化

通常情況下，對于單原子離子，粒子間的相互作用只包含非鍵相互作用，即靜電相互作用和范德華相互作用13，其中，在非極化力場框架下，靜電相互作用采取固定點電荷模型的庫倫公式，范德華相互作用一般采取12-6 Lennard-Jones (LJ)勢。在AMBER力場中，其函數形式如下：

由上式可知，靜電作用取決于原子所帶電荷，對于單原子離子，原子電荷與其價態相等。LJ參數中的ε和R分別代表兩個原子之間LJ關系中的勢阱深度及對應的原子間距，是準確描述單原子金屬離子范德華相互作用的關鍵。

2.1.1R和ε的關系

1924年，Lennard-Jones 提出LJ勢14的初衷是準確描述惰性氣體中原子之間的相互作用，其中r-12項表示由于近距離分子軌道重疊而引起的泡利斥力，r-6項表示色散力產生的長程分子間的吸引力。之后，這種描述方法被證明也能夠很好地描述中性原子和分子之間的范德華相互作用。1964年，Stokes運用quantum mechanical scaling principle(QMSP)方法15計算了與惰性氣體原子具有相同電子結構的金屬離子的范德華半徑。2013年，Merz等9在一項二價金屬離子范德華參數優化工作中依照此范德華半徑數值，通過在參數空間內尋找最優參數，證實單原子金屬離子的R和ε的關系與惰性氣體一致。同時提出惰性氣體的R和ε的關系滿足：

其中，F(x)代表-log(ε)，x代表R/2。

基于此關系，單原子離子的LJ勢UvdW與R、ε的關系，可以轉換為UvdW僅與R的關系。

2.1.2 R和IOD的關系

由上可知，單原子金屬離子的Uvdw與R值是一一對應的。在水溶液中UvdW值對應于M―Ow原子對的距離(ion-oxygen distance，IOD)，即在水溶液動力學過程中，相同價態的單原子金屬離子的R值與IOD值有一一對應關系。根據Merz等10報道的常見M3+/4+的力場參數及特定價態的單原子金屬離子在水溶液動力學模擬中R和IOD的變化趨勢，可得M3+/4+的IOD和R/2存在以下關系(圖1)：

本文選取SPC/E作為水模型，根據各錒系離子的IOD值得到其R/2值。

2.2 IOD值的選取

依照圖1，對于Np4+、Am3+、Cm3+離子，只要確定合適的An―Ow的鍵長(IOD)，就可得出該離子的LJ參數。三種離子的IOD值確定方法如下：

(1) Np4+：有多項實驗研究報道了Np的水溶液化學性質，包括其熱力學16、電化學17、光譜學18及配合物結構19等。根據Allen等20對Cl-濃度極低的水溶液環境中Np4+的水合物的擴展X射線吸收精細結構(EXAFS)光譜的解析，Np4+與Ow的配位鍵長為0.240 nm，配位數(CN)為11.2。另一項工作21解析了Np4+在1 mol·L-1HClO4溶液中的EXAFS譜，得到相同的配位鍵長。

由于相關報道中22各項工作溶液離子環境的不同，Np4+―Ow鍵的鍵長測量在0.239-0.240 nm之間，稀溶液中Np4+―Ow的鍵長多為0.240 nm。由于稀溶液中IOD值受反離子的影響最小，本文選取Np4+在水溶液中的IOD值為0.240 nm。

(2) Am3+：根據Stumpf等23對水溶液中的Am3+的EXAFS研究，得到Am3+―Ow的鍵長在0.247-0.249 nm，配位數為8.3，這與Allen等24在0.25 mol·L-1HCl溶液中得到的Am3+―Ow鍵長(0.248 nm)一致。本文選取0.248 nm作為Am3+在水溶液中的IOD值。

圖1 SPC/E水環境中R/2和IOD之間關系的多項式方程10Fig. 1 IOD as a function of R/2 in water (SPC/E) 10.(a) M3+; (b) M4+.

(3) Cm3+：根據時間分辨激光熒光光譜研究25，Cm3+的高氯酸鹽在水溶液中存在8配位水合物和9配位水合物物種的平衡，且以9配位為主；另外一項高能X-射線與EXAFS聯用研究26表明1 mol·L-1HClO4溶液中，Cm3+的水分子配位數為8.8。兩項工作都表明Cm3+外層水分子配位數接近于9。Allen等24對Cm3+溶液的EXAFS研究報道Cm3+―Ow鍵長為0.245-0.246 nm，在較低的離子強度下測得IOD為0.245 nm，CN = 10.2；而在較高的離子強度下IOD為0.246 nm，CN = 9.3。為了盡可能準確模擬Cm3+的外層水分子配位數，本文選取0.246 nm作為Cm3+的IOD值。

基于上述的各金屬離子的IOD值，優化得到相應的LJ參數(表1)。采用該套參數計算得到的金屬離子與SPC/E水分子模型相互作用的Lennard-Jones勢函數曲線如圖2所示。

2.3 計算方法

本文的分子動力學工作采用AmberTools1627程序完成，An3+/4+的范德華參數采用表1中的數值。水分子采用SPC/E模型28處理，Cl-的參數取自文獻29，CO32-、NO3-的力場參數使用antechamber程序30擬合，其電荷采用由半經驗量化程序sqm工具生成的AM1-bcc電荷31,32。動力學過程中的非鍵相互作用截斷距離為1 nm，長程靜電相互作用采用smooth PME的計算方法33。H原子與重原子之間的化學鍵采用LINCS限制算法34加以約束以保證1 fs時間步長的合理性。所有的模擬體系均包含約1500個水分子，分布在大小約2.9 nm × 2.9 nm × 2.9 nm的正方體盒中。各搭建好的模型體系經過以下四個階段獲取可用于分析的平衡動力學數據：(1)初始構型能量最小化；(2)升溫到目標溫度298 K：從250 K分兩階段至298 K，總時長100 ps；(3)預平衡過程：NPT系綜，T= 298 K，P= 1.013 ×105Pa，總時長500 ps；(4)采樣：在298 K、1.013 ×105Pa的NPT系綜下，采樣20 ns用于分析。溫控采用Berendsen弱耦合技術35，其中τT= 0.5 ps；壓力控制采用Berendsen壓浴方法，其中τP= 2 ps。

表1 幾種三價(Am3+、Cm3+)和四價(Th4+、U4+、Np4+、Pu4+) a錒系離子的IOD及擬合的LJ勢參數R/2和εTable 1 The IOD and LJ parameters (R/2 and ε) of several tri- (Am3+, Cm3+) and tetra-valent (Th4+, U4+, Np4+,Pu4+) a actinide cations.

圖2 根據表1中的LJ參數得出的An3+/4+-水的Lennard-Jones勢函數曲線Fig. 2 The Lennard-Jones potentials of An3+/4+-Ow using the parameters in Table 1.

金屬離子溶劑化自由能的計算采用自由能微擾(Free Energy Perturbation，FEP)的方法，將金屬離子從氣相微擾到水溶液相，每個微擾狀態由λvdW和λcoul兩個變量控制，分別表示金屬與水環境的范德華、庫倫相互作用微擾系數，每個微擾變量從0到1的步長為0.05，對應每個金屬離子的微擾過程共有41個窗口，每個窗口在NPT系綜下采樣1 ns，最后基于每個窗口間的交叉能量關系，采用Bennett acceptance ratio方法36計算出兩個微擾窗口間的自由能差異：

其中f是氣體常數R和溫度T的函數：

金屬第一溶劑化層水分子的平均駐留時間的計算首先根據下式計算相關函數R(t)37：

其中N0為時刻離子第一溶劑化層的總水分子數，Pi(t,t*)為階梯函數，當t時刻第i個水分子與0時刻第i個水分子相關，則為1。當且僅當第i個水分子同時在t和0時刻在第一溶劑化層，同時不能連續脫離第一溶劑化層超過t*時間長度，Pi(t,t*)才為1，否則為0。計算中使用的t*為10 ps，第一溶劑化層的截斷距離為0.282 nm。然后使用如下的公式計算駐留時間τ：

當R= e-1≈ 0.3679時，τ=t，即該值為駐留時間。

配位層水分子與體相水的交換速率與其反應活化能的關系用Eyring方程38表示，其形式如下：

其中kB為波爾茲曼常數(1.381 × 10-23J·K-1)，T為絕對溫度，h為普朗克常數，k為反應的速率，ΔG為反應的活化自由能。

3 結果與討論

3.1 An3+/4+水溶液動力學

3.1.1 An3+/4+溶劑化結構

首先計算了Th4+、U4+、Np4+、Pu4+、Am3+、Cm3+等錒系離子與水中氧原子Ow之間的徑向分布函數(RDF)及其積分值以分析其溶劑化結構，并與文獻值10,12,20,24,39-50進行對比，詳見圖3和表2。可以看出六種金屬都有明顯的第一、二水配位層，其第一配位層水的RDF峰值分別在0.243 nm (Th4+)、0.239 nm (U4+)、0.239 nm (Np4+)、0.237 nm (Pu4+)、0.248 nm (Am3+)、0.246 nm (Cm3+) (表2，圖3)，除了Th4+的配位數為10外其余都為9，這與參數擬合所用的IOD值20,24吻合。同時，我們的模擬結果也與已報道的相關理論工作吻合，如Yang和Bursten39采用量子化學的方法計算得出Cm3+與第一配位層水的配位鍵長為0.248 nm；Gagliardi等11采用NEMO力場在300 K下模擬得到的Cm3+―Ow的鍵長為0.255 nm，高于Yang等39、EXAFS測量值24及本套力場模擬值；Spezia等12采用極化力場得到Am3+―Ow和Cm3+―Ow的鍵長分別為0.248和0.245 nm，配位數均為9，與本文的非極化力場模擬結果一致；Rode等40采用量子力學電荷場(QMCF)動力學報道第一配位層U4+―Ow的鍵長為0.245 nm，略高于實驗值41(0.242 nm)及本文計算值；Jong等42使用CPMD方法計算得到的Pu4+―Ow的鍵長為0.247 nm，高于實驗值43(0.239 nm)，配位數(水分子)為9；Macedo等44采用BP86泛函方法計算出在COSMO溶劑模型下，Np4+、Th4+的水合物中金屬離子的配位數為9時，Th4+―Ow和Np4+―Ow的距離分別為0.247和0.239 nm，與本項工作吻合。

圖3 20 ns模擬過程的An4+/3+周圍水分子氧原子的徑向分布函數(黑色實線)及其積分(紅色虛線)Fig.3 Radial distribution function (solid black line) and its integral (dashed red line) of water molecules around An4+/3+ during 20 ns MD simulations(a) Th4+, (b) U4+, (c) Np4+, (d) Pu4+ (e) Am3+, (f) Cm3+. Color online.

實驗報道Th4+和Np4+的配位數分別是9-10和8-1020，且不同配位數的配位結構間能量差很小45。Yang等報道46的采用B3LYP泛函的[Th·9H2O]4+·3H2O與[Th·10H2O]4+·2H2O水溶液結合能僅差1.256 kJ·mol-1, 表明了9與10配位結構可以在溶液中共存，本項工作在力場水平上模擬出了一致的結果。

第二配位層的水分子受金屬的影響相對較弱。由于有關Am3+的第二水分子配位層的配位距離報道較少，而U4+、Np4+、Pu4+則缺乏相關的實驗數據，我們主要討論Cm3+的第二配位層結構。根據對錒系離子周圍水分子徑向分布函數(RDF)的分析(表2，圖3)，Cm3+與第二水分子層水的Cm3+―Ow鍵長為0.456 nm，這較接近于Spezia等12極化力場下的動力學結果(0.462 nm)。第二水分子層配位數為16，與Atta-Fynn等47的AIMD (15.8)、3-body MD (17.4)結果接近，小于Spezia等12(20.8)、Yang和Bursten等39(21)的模擬結果，說明第二水分子層的溶劑結構對于模擬體系、模擬方法的依賴性較高。

3.1.2 金屬離子配位層的水分子的運動性質

由上述金屬離子的溶劑化結構可知，An4+比An3+具有更緊湊的溶劑結構，表明An4+與溶劑的相互作用更強(圖4)。這種較強的相互作用不僅表現為更短的An―Ow距離，而且對金屬周圍水分子的空間排序、平動及轉動性質都帶來較大的影響。

An3+/4+第一配位層的水分子的空間取向高度有序，其偶極矩的朝向與An―Ow鍵軸方向總是保持一致。圖5表明，在1 nm的非鍵相互作用截斷距離內，An3+/An4+對外層水分子有序性都有影響，并隨著距離的增大水分子的有序性減弱。An4+對第一、二層水分子的有序性影響要強于An3+，表現在更小的偶極矩與An―Ow方向的夾角(第一配位層水分子：0°-33°vs0°-45°，第二配位層水分子：0°-70°vs0°-90°)。

金屬外圍水分子的高度有序給水分子的平動、轉動性質帶來了影響。由表3可得，在同樣九配位的An4+/An3+體系中，An4+前兩水分子層的擴散速度要慢于An3+體系，同時該部分水分子的轉動弛豫時間也要慢于An3+體系。這歸因于An4+與外圍水分子具有更強的靜電相互作用(圖4)。

3.1.3 第一配位層水分子的駐留時間與交換活化能

在20 ns的分子動力學模擬中，Am3+、Cm3+的第一配位層水分子較少發生與第二層的水分子進行交換，說明Am3+、Cm3+對第一配位層的水有較強的束縛能力，通過計算相關函數得出Am3+的第一層水駐留時間為18.4 ns，而對于Cm3+，該值無法在20 ns的時長采樣下準確測得，但根據前20 ns的軌跡分析，估算Cm3+的第一層水駐留時間約為34 ns。這與Dognon等51、Spezia等12的An3+在數納秒的溶液動力學內都未觀測到第一層水與第二層水的交換現象吻合。

對比之下，四價錒系離子的第一配位層的水分子在動力學過程中表現地比較活躍，并且隨原子序數的增大其平均駐留時間也增大。計算得到的An4+的駐留時間分別為0.19 (Th4+), 0.27 (U4+),0.5 (Np4+)和0.83 (Pu4+) ns。根據相關實驗報道52數據換算得到的Th4+與U4+的駐留時間分別是＜ 20 ns和180 ns。因此相對于實驗值，分子模擬過程的水分子交換過快，但是其變化趨勢是一致的。對于Np4+和Pu4+體系，文獻中鮮有相關的實驗和理論報道。從上述分子模擬的結果看，對于相同價態的An3+/4+，隨著原子序數的增大，對水分子的束縛能力增強，與An―Ow的鍵長值順序一致。

圖4 An4+/An3+與軌跡中穩定的第一水分子層間的相互作用能分解Fig. 4 The decomposition of the interaction energy between An4+/An3+ and the water molecules in the first coordination shell during MD simulations.(a) Th4+, (b) U4+, (c) Np4+, (d) Pu4+, (e) Am3+, (f) Cm3+.

圖5 在20 ns的動力學中An-Ow矢量與水分子偶極矩的夾角與An―Ow距離的關系，及構型間的Gibbs自由能差異(kJ·mol-1)Fig. 5 The population of the angle between An-Ow vector and water dipole moment as a function of the An―Ow distance,the populations were normalized and shown in the form of relative Gibbs free energy (kJ·mol-1).

由于An4+外圍水分子高度有序，可以認為熵效應并不是影響An4+第一水分子層較短駐留時間的主要原因。通過分析溶劑間的相互作用，發現An4+緊湊的溶劑化結構會帶來更多溶劑間的范德華互斥和靜電互斥作用(圖6)，比An3+溶劑間的能量高約84-155 kJ·mol-1，這促進了第一水分子層向第二水分子層的轉移。

表3 金屬離子擴散系數(DAn，10-5 cm2·s-1)，第一、二配位層水分子的擴散系數(Dwat，10-5 cm2·s-1)和轉動弛豫時間(τ，ps)，第一配位層水分子的平均駐留時間(τr，ns)和根據 Eyring方程 38計算的活化能(ΔG≠，kJ·mol-1)Table 3 The diffusion coefficients of An4+/An3+ (DAn, 0-5 cm2·s-1), the diffusion coefficients (Dwat, 10-5 cm2·s-1) and rotational correlation time (τ, ps) of water in the first and second hydration shells, the mean residence time (τr, ns) of water in the first hydration shell and corresponding activation energies (ΔG≠, kJ·mol-1) calculated by Eyring equation 38.

圖6 動力學過程中穩定水合結構的連續軌跡中第一配位層水分子內部相互作用的能量分解為靜電相互作用(黑色)和范德華相互作用(紅色)Fig. 6 The energy decomposition of interactions between intra-cluster water molecules in the first hydration shell: electrostatic (black) and vdW (red) interactions. Only continuous trajectory of stable hydrated structures without water exchange were analyzed.Average values of these energy components (kJ·mol-1) are shown. Color online.

3.1.4 配合結構立體構型

根據對軌跡的分析，在298 K，除Th4+外的五種錒系金屬離子的配位數都保持為9。從立體化學的角度，金屬及其外層的9個配位原子可形成包括單帽四方反棱柱體(capped square antiprism，CSAP)和三帽三角棱柱體(tricapped trigonal prism，TCTP)的立體結構53。

根據對動力學模擬軌跡的分析(圖7)，[ThIV(H2O)10]4+的第一配位層立體化學結構為雙帽四方反棱柱體(bicapped square antiprism，BCSAP)；而[UIV(H2O)9]4+，[NpIV(H2O)9]4+，[PuIV(H2O)9]4+采取的立體化學構型為CSAP，與Yang等46報道的采用量子化學的方法證明九配位的[ThIV(H2O)9]4+的CSAP構型比TCTP構型更穩定的結論相似。Am3+的相關報道較少，但一般認為An3+與Ln3+水合物的立體結構相似，均采取TCTP構型。一項An3+/Ln3+的配位化學研究54提及[AmIII(H2O)9][CF3SO3]3中為TCTP構型，但我們的純水溶液動力學中觀察到Am3+保持CSAP構型，這兩者的差異可能是因為單晶X-ray的測量方法中加入的[CF3SO3]-反離子分布在第二配位層，對第一配位層水分子的排列產生的影響導致的55。[CmIII(H2O)9]3+為TCTP構型，與Yang和Bursten等39運用量子化學方法得出相同的構型，表明該套參數的動力學模擬可以較好描述金屬第一配位原子層的立體化學結構。

圖7 MD捕捉到的An4+/3+的水合物的優勢立體幾何結構Fig. 7 The dominant conformations of the hydrated An4+/3+ during MD sampling.(a) Th4+, (b) U4+, (c) Np4+, (d) Pu4+, (e) Am3+, (f) Cm3+.

3.1.5 溶劑化自由能(HFE)

通過運用自由能微擾方法，計算得出六種錒系金屬離子在該套參數描述下的溶劑化自由能ΔG(FEP)。但由于非極化力場函數式中未引入極化項，忽略了誘導偶極作用等相互作用，使用FEP算出的自由能往往會低估溶劑化自由能13。Merz等29對比一系列M3+與M4+的非極化力場的IOD值參數計算出的自由能與實驗值比較，總結出在非極化力場下運用FEP的方法，對于三價單原子金屬離子，往往會低估約9%的自由能；對于四價單原子金屬離子，會低估約16%的自由能13。我們依據此結論，對ΔG(FEP)進行校正，得出校正后的自由能ΔG(CORR)，見表4，并與文獻值進行對比56-59。

根據計算結果，Np4+和Pu4+的溶劑化自由能最大，Cm3+比Am3+的溶劑化自由能略高，這與文獻報道的結果一致，表明本項工作擬合的An3+/An4+參數能夠定性描述An3+/An4+與水的相互作用強弱。

從結果也可看出，隨著原子序數的升高，相同價態的An3+/4+的溶劑化效應增強。這主要是由于隨著原子序數的增大，An―Ow之間的斥力減小(圖2)，有助于An與水的相互作用。相比于An3+，四價離子較高的溶劑化自由能得益于較大的靜電相互作用的貢獻(圖8)。

3.2 An3+/4+與Cl-、NO3-、CO32-的溶液動力學

在鑭錒分離及其分析過程，限于實驗條件或測量技術，溶液當中會伴隨著反離子的存在，如Cl-、NO3-、CO32-等，這些無機陰離子會影響液-液萃取平衡及分析測量的結果，因此有必要研究這些無機陰離子與錒系金屬離子共存體系的溶液動力學。在本項工作的模擬體系中，為了探究金屬離子與該類無機陰離子的特征配位化學行為，金屬離子與無機陰離子均放置1個，以簡化環境因素的影響，分析結果與文獻值進行對比20,21,24,60-69。

Cl-在錒系單原子離子測量方面是一種常見的反離子，在之前Np4+、Am3+、Cm3+的相關EXAFS工作20,24中給出了不同Cl-濃度下An3+/4+―Cl的鍵長，我們的模擬結果與實驗值有較好吻合(表5)。An3+/4+與Cl-的RDF結果(圖S1，見Supporting Information)表明Cl-與Cm3+、Th4+表現出較弱的相互作用，在動力學過程中Cl-在前4配位層之間不斷交換，尤其是Cm3+，在2-3配位層的時間長于停留在第一配位層的時間，表明了Cm3+―Cl-弱于Cm3+―H2O的相互作用。其他四種金屬離子均可以將Cl-束縛在第一配位層，且An4+對Cl-有較短的配位鍵長，其作用力要強于An3+。Seminario等60運用DFT的方法計算了10配位的[ThCl3(H2O)7]+配合物在溶劑環境下Th4+―Cl的鍵長為0.272 nm，在同樣條件下8配位的PuCl4(H2O)4中Pu4+―Cl的配位鍵長為0.262 nm，與本項工作的動力學模擬結果誤差在±0.002 nm。而U4+與Cl-在水溶液配位結構的表征工作較少，根據一項1200 K高溫下UCl4-NaCl熔鹽體系的第一性原理動力學模擬(AIMD)的工作61，U4+―Cl的距離為0.263 nm。整體來看，相同價態下，An3+/4+與Cl-的相互作用隨著原子序數的升高而增強。

表4 An3+/4+的溶劑化自由能(kJ·mol-1)Table 4 The calculated hydration free energy (HFE) of An3+/4+ (kJ·mol-1).

圖8 軌跡最后5 ns的金屬與水相互作用的能量分解：靜電相互作用(黑色)和范德華相互作用(紅色)Fig. 8 The energy decomposition of non-bonded interaction into electrostatic (black) and vdW (red) interactions during the last 5 ns of production run.Average energies (kJ·mol-1) are showed. Color online.

表5 根據基于20 ns軌跡的RDF計算的An3+/4+與Cl-、NO-3、CO23-配位原子的配位鍵長Table 5 The coordination bond distance (nm) between actinide ion and coordinating atoms of Cl-, NO-3, and CO23- derived from calculated RDF over the 20 ns production run.

硝酸根離子在鑭錒分離的萃取過程當中經常出現。在水溶液中，NO3-的配位能力較弱，我們的模擬結果表明六種金屬離子在20 ns的采樣過程中始終保持著第一配位層的單齒配位。Guibaud等62利用傘形采樣的方法計算了NO3-與鑭系金屬Nd3+和Dy3+的單配位與雙配位的能量，發現硝酸根在水溶液中從單配位到雙配位需要跨越較大能壘，并且雙配位模式不穩定，這在我們的模擬中也有所體現。從模擬的平均鍵長，硝酸根的配位氧原子(ON)與 An4+的 配 位 鍵 (An4+―ON)的 鍵 長 較Am3+/Cm3+―ON的鍵長短約0.010-0.018 nm。根據密度泛函理論(BP/SVP)研究63，在An(NO3)3·5(H2O)的三價錒系配合物中，2個NO3-是雙齒配位，一個NO3-是單齒配位時，總配位數為10，Am3+與其單齒配位的O的鍵長為0.255 nm，Cm的相應鍵長為0.251 nm。而對Np4+的HNO3水溶液的EXAFS光譜研究21測得在5.0 mol·L-1HNO3水溶液中，與Np4+配位的氧原子中有約4.5個來自硝酸根，另有7個水分子，配位數為11.5，Np4+和硝酸根配位O原子的配位距離為0.252 nm，這比我們的模擬結果略高，可能是由于研究體系的硝酸根數量、配位模式及金屬總的配位數等差異導致。通常情況下，金屬的高配位數會弱化配位鍵而致使配位鍵增長。另外兩項工作64,65分別表征了[Th(NO3)6]2-和[U(NO3)6]2-的晶體結構，NO3-均采用κ2的配位模式，由于金屬的高配位數12，其An4+―ON的距離大于我們的水溶液模擬結果。整體來看，同樣價態下，隨著原子序數的增高，An―ON的配位距離變短。

圖9 計算 An3+/4+與 Cl-、NO-3、結合自由能所采用的熱力學循環示意圖Fig. 9 The schematic of thermodynamic cycle used for calculating binding energy between metal ions and anion.

表6 熱力學循環計算過程中Cl-、、溶劑化能ΔG2及其與An3+/4+的復合化能ΔG1、結合自由能ΔGbind (單位：kJ·mol-1)Table 6 The anion’s complexation free energy ΔG1 and binding free energy ΔGbind (in kJ·mol-1) with An3+/4+.

表6 熱力學循環計算過程中Cl-、、溶劑化能ΔG2及其與An3+/4+的復合化能ΔG1、結合自由能ΔGbind (單位：kJ·mol-1)Table 6 The anion’s complexation free energy ΔG1 and binding free energy ΔGbind (in kJ·mol-1) with An3+/4+.

a The anion’s hydration free energy ΔG2 were calculated to be -371.8 ± 0.24 kJ·mol-1 for Cl-, -370.1 ± 0.21 kJ ·mol-1 for NO-3,and -1416.8 ± 0.52 kJ·mol-1 for, respectively.

Cl- NO3- CO3 2-ΔG1 ΔGbind ΔG1 ΔGbind ΔG1 ΔGbind Am3+ -380.2 ± 0.29 -0.84 -394.8 ± 0.27 -24.7 -1552.9 ± 0.59 -136.1 Cm3+ -371.0 ± 0.29 0.4 -399.4 ± 0.31 -28.9 -1544.9 ± 0.61 -128.1 Th4+ -370.1 ± 0.32 1.3 -411.6 ± 0.26 -41.4 -1671.4 ± 0.52 -254.6 U4+ -417.0 ± 0.28 -45.2 -428.3 ± 0.28 -58.2 -1671.4 ± 0.48 -254.6 Np4+ -417.8 ± 0.23 -46.5 -427.9 ± 0.33 -57.8 -1676.4 ± 0.44 -259.6 Pu4+ -410.3 ± 0.26 -38.9 -425.4 ± 0.25 -55.3 -1692.3 ± 0.52 -275.5

在開放體系中，空氣中的CO2會溶于水相并隨著體系pH的變化表現出不同的物種分布，CO32-是最為常見的物種之一。在動力學模擬中，溶液中的CO32-均以κ2模式與六種金屬配位，且在整個動力學模擬時間段內保持穩定，另有7個水分子在第一配位層分別與除Th4+外的五種金屬離子配位形成9配位結構，對于Th4+則有8個水分子在第一配位層形成10配位結構。根據Spezia等50基于極化力場對[Cm(CO3)3]3-的水溶液動力學的模擬研究，三個碳酸根均以雙齒配位與Cm3+結合，另有3個水分子與Cm3+共形成9配位結構，配位模式與配位數與本項工作一致。作者還根據不同計算方法得出Cm3+與CO32-上的配位O原子的配位鍵長分別為0.234-0.241 nm (POL)、 0.241 nm (DFT)、 0.241 nm(CPMD)，略高于本項模擬工作。一項實驗工作66合成并表征了[C(NH2)3]6[NpIV(CO3)5]·4H2O配合物，發現共有5個CO32-以雙齒配位的模式與Np4+配位形成10配位結構，[C(NH2)3]+與H2O分布在第二配位層，并統計出10個配位氧原子與Np4+配位的平均距離為0.244 nm，比使用本文報道的參數的模擬結果大0.023 nm，但由于配合物結構相差較大，CO32-對金屬配位具有較強的影響，兩組數據無法進行直接對比。

為了分析An3+/4+與Cl-、NO3-、CO32-在水溶液中的結合強弱，我們采用圖9所示的熱力學循環來計算陰陽離子的結合自由能。由于有限尺寸給自由能計算帶來的誤差與微擾粒子的電荷成正相關70，為了盡可能減小該部分誤差，我們選擇微擾具有低價態的陰離子，而不是高價態的An3+/4+。在圖9中，ΔG1和ΔG2表示目標陰離子從dummy態分別微擾到復合物水溶液和純水溶液的自由能變化，ΔG3表示該陰離子在復合物水溶液體系及純水體系中的空穴能差異，由于該能量差很小往往忽略不計。陰陽離子的結合自由能可以計算為：

自由能計算結果表明(表6)，CO32-對An3+/4+的結合能力要遠強于NO3-及Cl-，自由能微擾過程中的能量積分表明這主要來源于高電荷帶來的強靜電相互作用貢獻(圖S2)。NO3-的結合能力略強于Cl-，這主要由于Cl-與An3+/4+較高的范德華斥力導致。整體看來，靜電相互作用是An3+/4+與三種陰離子的結合的主要貢獻，范德華作用均表現為斥力，因此，對于同種陰離子，An4+的結合能力要高于An3+。對于An4+，Th4+與Cl-、NO3-的結合自由能要明顯弱于其他三種An4+，這與Th4+較高的配位數有關(表2)。

4 結論

本文基于惰性氣體的LJ勢中R和ε的函數關系，以及文獻報道的Np4+/Am3+/Cm3+―Ow鍵長的實驗值，擬合出該三種錒系金屬離子的AMBER力場參數并評估了其可靠性。結合已發表的Th4+、U4+、Pu4+參數，系統研究了六種錒系離子的溶液配位化學，獲得了與文獻一致的六種金屬離子在水溶液中的水合結構。從立體化學的角度，10配位的Th4+的水合結構采取BCASP構型，9配位的U4+、Np4+、Pu4+、Am3+的第一配位層原子采用CASP的構型，而9配位的Cm3+相應原子采用TCTP構型。研究進而采用FEP方法計算了六種錒系離子的溶劑化自由能，得到了與之前工作吻合的結果，并從能量分解的角度解釋了An4+的溶劑化效應要遠強于Am3+及Cm3+的原因，提出這種差異主要來源于An4+與水之間較強的靜電吸引作用。

采用該套力場參數，本文也分析了金屬對外圍水分子結構、平動與轉動性質的影響，發現可以提供金屬離子近程水分子的分布和運動的合理圖像，并可得出與文獻一致的Cm3+和Am3+的第一配位層水分子的駐留時間，但對于四價錒系離子，計算得到的駐留時間顯著小于文獻提供的推測值。雖然文獻報道的數值自身存在較大的誤差，我們謹慎建議本文報道的參數適用配位結構和結合自由能的研究；在用于動力學性質的研究時，僅可用于相同價態的錒系離子溶液體系的定性比較。另外，由于該套參數基于金屬與水之間的相互作用擬合得到，宜用于金屬與軟硬酸堿理論中較“硬”配體的相互作用的模擬；當用于金屬與較“軟”配體的配合物的模擬時，可能會帶來較大誤差，需審慎對待。

此外，我們運用此套參數模擬了六種金屬離子與重要的無機陰離子Cl-、NO3-和CO32-的溶液動力學的模擬研究，發現NO3-均以κ1模式配位，CO32-均以κ2模式配位。并通過熱力學定性分析陰陽離子的結合強弱。通過對比文獻報道的實驗及理論研究結果，發現采用該套參數可以獲得合理的配位模式、鍵長及結合強度等配合物結構信息，豐富我們對錒系溶液化學的認識。

Supporting Information:available free of chargeviathe internet at http://www.whxb.pku.edu.cn.