极端条件下材料化学和输运性质的第一原理研究

摘要

地球是一个动态系统。从常温常压到地核的360 GPa 和6600 K，地球的热力学状态随着深度的变化而发生很大变化。因此，地球成分（如硅酸盐和碳酸盐矿物）的物理和化学性质受其所处环境的显著影响。在过去的30 年中，在极端条件下的材料表征的实验技术和理论模拟方法都取得了巨大的进步。这些进步提升了我们对矿物特性的理解，这对于充分认识这个星球的形成及其生命起源至关重要。本文回顾了用于预测极端条件下材料性质的最新计算技术，本研究仅限于应用第一原理分子动力学（first principle molecular dynamics, FPMD）方法研究与地球科学相关的化学和热力学输运过程。

正文

《1. 背景介绍》

1. 背景介绍

具有讽刺意味的是，虽然我们对太阳系以及周围的恒星和星系都有了一定程度的了解，但对我们自己星球的成分及其相互作用的准确了解却相当有限。由于有史以来最深的钻孔仅穿透地壳12 km [1]，因此，原则上，我们没有直接证据证实深藏在地球内部物质的特性。我们关于地幔和地核组成的大部分信息都是间接获得的，如地震学、火山活动引起的出露矿物或钻石包裹体。因此，在极端条件下对材料的物理和化学性质的深入和精确认识是至关重要的，以便基于这些知识构建连贯、合理和可靠的认知。在早期，人们就利用各种岩石的弹性常数建立随深度变化的密度剖面。从那时起，理论计算和模拟成为地球科学中不可或缺的一部分[2]。基础物理和化学原理的融合导致了许多成功的发现[3]。其中一个例子是化学成键理论，它是在19世纪后期发展起来的，后来被许多科学家完善。通过这一理论，只需要几个关键结构参数和物理性质的信息就可以预测和理解常压下许多矿物的晶体结构和性质[4,5]。这一理论领域在过去20年中取得的最大进步可以说是现在许多物理量能够从第一原理中可靠地计算出来[6,7]，这种方法基于量子理论，只用输入组成元素来解薛定谔方程。最近，即使是晶体结构也能可靠地预测。计算机建模已经成为一种新的范式，它为研究材料在极端压力和温度下的行为开辟了一个新的维度，而这种压力和温度是不容易进行实验研究的。本文回顾了几种计算方法，重点介绍了它们在化学反应和输运性质研究中的应用。本次简短回顾的目的不是对这一领域的技术现状进行全面的调查，也不是收集所有相关的文章，相反，我们的目标是以作者最近在地球科学方面的工作为例，为广大读者提供这个领域的研究现状概要。

材料的性质完全由其原子的原子间力决定。一旦原子间力被定义，原则上，所有的性质都可以计算出来。例如，对应于元胞中原子最低能量排列的晶体结构；弹性是系统对外部应力的响应。利用经典分子动力学（molecular dynamics, MD）和统计力学理论[8]，还可以计算有限温度下的动力学和输运性质。例如，导热系数是能量通过材料的传递，而液体的剪切黏度是由应变速率相关的变形引起的内力所产生的弹性（黏性）应力。从数值模拟中计算材料性能并不新奇。事实上，自从MD方法在20世纪60年代被应用以来，几乎所有的理论形式都已经建立起来了[9−12]。因此，一旦通过拟合实验观测值（经验值）或者原子或分子团簇的高级量子化学计算值（非经验值）得到原子间势能函数并定义系统中原子间力，就可以预测该系统的各种性质。例如，Si−O、O−O和Si−Si原子间相互作用势已经构建得非常成功[13, 14]，并证明了它们能够在低压下重现晶体二氧化硅和非晶二氧化硅族的实验结构、弹性常数、振动光谱和结构相变[15]。多组分体系的势函数也为许多硅酸盐和铝硅酸盐固体和熔体的结构和动力学提供了有用的信息[16]。具有可靠原子间势的经典MD计算仍然是使用当今计算机技术的有力预测手段；在数以百万计的原子上，可以模拟10~100 ns级的持续时间。研究慢松弛液体和玻璃通常需要很长的模拟时间。经验原子间势的一个主要缺点是它只能在一定范围条件下应用。因此，大多数经验势是不可转移的，并且会影响预测的可靠性。在地球条件下，在较宽的温度和压力范围内，矿物的行为不能用经验原子势来描述。为了克服这一困难，有必要采用无可调参数的第一原理量子力学方法。这就是本文的主题。

《2. 理论背景》

2. 理论背景

根据量子力学理论[7,17]，多粒子系统（电子、原子核等）的电子状态由一个波函数，所唯一确定，而这个波函数取决于电子坐标和原子坐标。系统的总能量是哈密顿能量的期望值，。为哈密顿总能量，是动能和势能的简单相加。如果电子和原子核的运动（玻恩-奥本海默近似）能够分离，电子的哈密顿量就可以写成：

式中，m是电子的质量；势能U是电子与电子、电子与原子核、原子核与原子核库仑相互作用之和。用量子表达式替换动量算符，，哈密顿能量变为：

这就是系统的薛定谔方程。除了氢原子，电子波函数的确切形式是未知的，并且对于多粒子系统来说，波函数可能非常复杂。一种解决方法是假设系统中的每个电子都只是经历了其他电子的平均场作用，并且系统中的每个电子都有独特的波函数（独立粒子近似）。然后，多电子体系总的波函数，（由于泡利不相容原理而反对称）可以近似为单个和独立轨迹波函数乘积的行列式。将总波函数代入式（2），利用变分原理使电子轨迹总能量最小化，得到基态波函数和能量。这个过程产生了一组耦合微分方程，可以用迭代的方式求解。这种方法被称为Hartree-Fock方法，已经被证明对分子非常有效。然而，由于需要大量计算多中心双电子的库仑作用，，和交换积分，，这给体系的扩展带来了困难。

Hohenberg和Kohn [18]在1964年对电子结构理论进行了另一种不同的表述，这种表述只基于电子密度，而没有显式地调用波函数。密度泛函理论（density functional theory, DFT）证明了电子密度是主要变量，一个相互作用的多电子系统的总能量是电子密度的函数。后来，Sham和Kohn [19]提出了一种计算电子系统总能量的实用方案。DFT理论消除了计算交换积分的负担，用密度函数代替了电子交换和相互关联作用，尽管这个函数的确切形式未知。系统总能量的Kohn-Sham方程为：

式中，n是位矢处的电子密度。函数F[n]是电子与电子、电子与原子核相互作用势能以及电子动能之和。同样，函数F[n]的精确形式依然未知。很容易将电子库仑相互作用的贡献从总电子相互作用中分离出来。F[n]可以进一步分解为：

式中，G[n]仍为含有密度n的非相互作用电子体系动能T[n]的泛函[19]，其余电子相互作用项归为E_xc [n]，称为交换相关联项（exchange correlation, XC）。

根据变分定理，在系统中粒子数N守恒的约束下，使电子密度对基态总能量最小，得到基态总能量。也就是说，。因此，取方程（3）的微分得到：

式中，μ是具有电子密度n的均匀气体的化学势。XC势定义为。

由于相互作用的电子系统的确切动能泛函是未知的，因此有必要将多体波函数近似为单电子波函数的乘积，如Hartree-Fock方法：。相互作用电子系统中的每一个电子都满足单电子Kohn-Sham (KS)方程：

其中，有效势定义为：

为了将KS DFT应用于实际问题，必须对XC泛函E_xc [n]作近似。第一个近似是局域密度近似（local-density approximation, LDA），其中非均匀真实系统中每一点的XC势[19]被设为相同电子密度[20]下量子蒙特卡罗计算的均匀电子气体的交换相关势。简单的LDA近似忽略了电子密度的不均匀性，因此往往高估了接近平衡的化学键的强度，低估了原子间的距离。通过加入电子密度梯度[广义梯度近似（generalized gradient approximation, GGA）] [21]来考虑电子密度的曲率，可以改进预测结构，有时也可以提高总能量，但不能两者兼顾。这是由于GGA不满足所有已知半局域函数的相关约束。为了解决这一问题，引入了动能密度来解释金属键电子密度的缓慢变化和共价键电子密度的紧密变化。这种类型的函数称为meta-GGA。最近，构造了一个“强应变近似赋范”（strongly strained and approximately normed, SCAN）meta-GGA函数，以满足半局域泛函[22]交换相关能所需的所有17个已知精确约束。由于轨迹φ_i 已经被计算出来，将非相互作用的动能项归入到半局域交换关联能中，几乎没有增加计算成本。SCAN泛函已证明，与所有其他可用的功能[23]相比，其具有更好的几何形状和总能量。它是迄今为止对主族化合物最有效和可靠的泛函。

DFT计算的一个关键因素是对作用于原子上的力的评估。这些信息对于几何优化和MD模拟是必不可少的。作用于原子上的力F_i 是能量对位移 λ_i 的负导数，可以通过赫尔曼-费曼定理[24]方便地求出：

式中，波函数导数的烦琐积分被简化为更易于处理的哈密顿能量导数的表达式。

一旦选择了泛函和基组[25]，就可以计算出孤立或扩展系统的结构和能量。自20世纪80年代以来，DFT在原子、分子和凝聚态电子结构计算中的应用取得了巨大的成功。

一旦作用于原子上的力和超胞模型被确定，原子的轨迹就可以根据玻耳兹曼统计进行抽样。在相空间任意温度压力下的原子轨迹作为时间的函数都可以通过求解经典的牛顿运动方程获得[26]，，式中，m_i 是原子质量，v_i (t)是瞬时速度。MD计算得到的有用结果是每个原子的瞬时位置、速度和力。从这些量可以计算出平衡态和热输运性质[8,26]。给定系综H以及温度T下的平衡热力学变量（如A），采用遍历假设的统计理论计算得到，即热力学平均可以近似为计算量的时间平均为玻耳兹曼常量[27]：

基于Green和Kubo[8−12]提出的时间相关形式，平衡MD模拟还可以计算热输运特性。从本质上讲，该理论假定输运过程的动力学是随机的（马尔可夫过程），热平衡态的扰动很小。随机过程（A）的耗散率的特征是由衰减时间相关函数（如果A是一个复动态量）所确定的。例如，扩散速率（D）是原子（扩散）的均方位移（位置）随时间变化的函数。可以看出，随着时间t→∞，速度自相关积分得到扩散系数[8]：

通过类比，可以建立黏度的时间关联函数。如上所述，黏度是液体对弹性变形的响应。变形模型系统中粒子的瞬时原子位置与参考量有关，由通过变形矩阵h的线性变换得到；h是参照系受无穷小应变ε，由h = I + ε变形而得到，其中，I为单位矩阵。在t = 0时刻，应变被切断，系统剪切速率将经历δ-函数峰，。剪切应力与峰的时间关系[8]为：

应力张量(σ)描述了基向量在应变ε_ij 的无穷小变形下总能量的变化。应用赫尔曼-费曼定理[式（9）]，可以估算应力分量：。稳态剪切产生的应力（如P_xy ）为随时间的积分，剪切黏度由式（13）[8]给出：

地球科学的另一个有用的输运性质是热导率κ。导热系数是累积的矢量热流 ——单位面积和单位时间内的能量动量传递速率，。导热系数的表达式为热通量自相关函数形式[8, 28]：

注意，计算扩散速率、黏度和导热系数所需的信息是瞬时位置、模型晶胞上的剪应力和粒子能量。所有这些量均可通过MD计算得到。时间关联函数是一种功能强大的技术手段，其中许多动态特性，如振动态密度（vibrational density of states, VDOS）、红外和拉曼光谱，都可以从原子轨迹的MD计算中获得[26]。在接下来的讨论中，我们将举例说明这些方法在极端条件下对地球科学中几个选定问题的应用。

《3. 举例说明》

3. 举例说明

在本节中，将阐述上述方法在几个与地质学有关的问题上的应用。首先， “静态”DFT计算用于研究硅（Si）和氧化锗玻璃[29]的结构变化。结果表明，传统的原子尺度刚性球概念在高压下已不再适用。在玻璃结构中，在任何给定的压力下，都没有独特的Si−O或Ge−O配位。因此，任何试图将配位数（coordination number, CN）与玻璃属性联系起来的尝试都是无效的。本文作者的研究团队先前描述的DFT计算方法，结合先进结构搜索技术，被用来预测高压[30]下铁（Fe）卤化物的稳定多态性。结果表明，地核中可能存在富铁卤化物，这表明重卤素可能被分异在地核中，为陆地挥发性卤素的耗尽提供了一个合理的解释。另外，采用第一原理分子动力学（first-principles molecular dynamics, FPMD）计算方法，研究了二氧化碳（CO₂ ）与二氧化硅（SiO₂ ）[31]之间可能发生的反应。研究发现，在相对较低的压力和温度下，分子筛孔隙中的CO₂ 与SiO₂ 发生化学反应。但是在过渡区的压力和温度条件下，液态CO₂ 与石英晶体和辉长石界面接触时不发生反应。在另一项研究中，我们发现H₂在较低的压力和温和的温度下很容易与石英反应生成水[32]。剩下的讨论将致力于计算传输属性。我们考察了纯碳酸钙（CaCO₃ ）熔体的黏度和碳酸钙（CaCO₃ ）熔体在少量水[33]存在下的黏度，发现含水熔体的黏度低于纯熔体。我们还讨论了一种新的计算有序和无序固体在高压和高温[34]下导热系数的有效近似方法。

《3.1. SiO2 和 GeO2 的玻璃结构》

3.1. SiO₂ 和 GeO₂ 的玻璃结构

由于高温高压下矿物熔体的性质难以测量，淬火玻璃的结构常被作为理解其高温高压下熔体结构的近似。如上所述，原子大小和原子价的概念最初是由Pauling [4]及O’keeffee和Navrotsky [5]提出的，是阐明各种矿物在常压下晶体结构的一个非常有力的工具。将这一强大的概念推广到无序玻璃并非易事。然而，根据一定的观测参数对高压液体和玻璃的性质进行分类是值得的。最近，人们试图将各种氧化物玻璃的体积特性与氧填充率联系起来[35]。这一想法似乎不无道理，因为在高压下，氧化玻璃中体积较大的氧（O）原子（氧化阴离子）相互压得很近，从而改变了局域环境，并且会影响某些体积特性，如压缩性和黏度。难点在于确定一个唯一的O离子半径，这是计算高压玻璃中氧填充率所需要的。为此，人们认为离子半径与阳离子-氧CN和键长有关，并可由一组经验关系式导出[35]。由于局域结构（如CN）可能会随着压力的变化而变化，所以在不同的压力状态下必须使用不同的公式。要应用该方法，必须先确定氧化玻璃中的局域结构（即4、5或6配位）。这就提出了一个问题：如何指定氧化玻璃在给定的压力下单个CN的有效性。此外，由于单个原子的“大小”受周围环境（如填充率）的影响，应不同于环境压力下的原子“大小”，因此没有理由认为给定的化学键类型在较宽压力范围内的原子半径比是恒定的。为了解决这个问题，利用FPMD计算得到的理论结构[26]，分析了不同压力下SiO₂ [36]和锗（GeO₂ ）[29]玻璃的结构及相关配位参数。玻璃中CN的定义有些武断，因为它取决于截止分离半径的选择。在研究中，选取Si–O和Ge–O径向分布函数中第一个最小值的距离作为截止分离半径。图1和图2简洁地总结了GeO₂ 的结果。同样的结果也出现在SiO₂ 玻璃上。显然，除了在最低压力下，没有任何单一的CN可以代表玻璃结构。结构总是由不同CN的混合物组成。图1（a）为不同CN在GeO₂ 玻璃中的比例与压力的关系。从0到15 GPa，随着5和6配位的增加，4配位稳步下降。在25~70 GPa，玻璃已经混合了5、6配位和7配位。大于7的配合量随着压力的进一步增加而增加。这些理论结果表明，对于任何压力下的玻璃，没有任何依据假设其CN是唯一的。下一个问题是：原子的“大小”是如何随着压强而变化的？不论压力大小，Ge(Si)和O原子的原子大小的变化是均匀的吗？当Ge(Si)–O距离随压力变化时，Ge(Si)/O原子大小的比值是否仍然是常数？即使在晶体中，由于不同的键长，在扩展固体中原子大小的分配也是模糊不清的。为了消除这种随机性，Gibbs等[37]基于Bader[38]首创的原子-分子（atom-in-molecule, AIM）量子理论，从电子密度分布的数学分析出发，提出了计算原子大小的方法。根据该理论，分子或固体中的原子区域是由电子拓扑结构所唯一定义的。原子的Bader体积是围绕原子的零电子通量表面所包围的空间。这个严格的定义基于量子力学理论，消除了在极端压力下，为玻璃状固体的严重扭曲的原子密度分配球形半径的任意性和模糊性。例如，这种划分原子体积的方法已被证明成功地解释了结构变化对金属玻璃电子性质和结构转变的影响[39]。图1（b）比较了锗（Ge）和O在GeO 2 玻璃中的Bader体积随压力的变化。正如所料，原子“体积”随压力的增大而减小。平均O Bader体积在0 GPa时是12.5 Å³ ，在120 GPa时下降到7.5 Å³ ，下降了40%。相比之下，Ge Bader体积的变化要小得多：从0 GPa时的6.2 Å³ 到160 GPa时的5.2 Å³ ，下降了16%。最重要的是，Bader体积比随压力的变化而变化，并不是一个常数。该理论研究结果驳斥了氧填充率假说。

《图1》

图1.（a）研究了GeO₂ 玻璃中不同Ge−O配位含量随压力的变化；（b）GeO₂ 玻璃中Ge和O Bader体积随压力的变化。

《3.2. 卤素缺失悖论》

3.2. 卤素缺失悖论

挥发性元素的耗竭与星云挥发性有关。碳质（carbonaceous, CI）球粒陨石是一组原始陨石，其元素分布与星云的组成极为相似。比较地球元素浓度与CI年代学资料丰度，为研究地球早期的元素聚集提供了有价值的资料。令人惊讶的是，地球上的重卤素溴（Br）和碘（I）的丰度比它们在CI年代仪中的丰度低一个数量级。这一缺陷与地球的形成和成分演化有关，并可能影响当今海洋的盐度。对于这个“卤素缺失悖论”的假设是，卤素可能溶解在早期地球的液态铁中，然后在结晶后被分异到地核中。为了限制卤素在地核中的浓度，了解稳定的卤化铁化合物能否在地幔和地核的压力下形成是很重要的。为此，利用基于智能的结构搜索技术，对卤化铁的相图进行了研究。这种方法结合了第一原理DFT计算和结构搜索技术[40,41]，是一种非常强大的手段，可以在给定热力学条件下，仅考虑组成原子的信息，以无偏向的方式定位最稳定的晶体结构。研究表明，这种计算方法为研究无法通过实验修正的压力条件下材料的晶体类型和化学键形成提供了新的思路。搜索过程如下：随机生成一组不同Fe/卤素化学计量的初始结构。首先利用DFT算法对试验结构进行优化，然后利用粒子群优化算法对试验结构进行改进[42]。选择能量较好的结构作为下一个优化周期的试验结构。重复这个过程，直到没有发现新的有利的能量结构。以卤化铁为例，该方法成功地再现了已知的FeCl₃ 、FeBr₃ 、FeCl₂ 、FeBr₂ 和FeI₂ 的环境结构。为了研究高压多晶型的结构和稳定性，在过渡区（25 GPa）、核幔边界（70 GPa）和深部地核（135 GPa和360 GPa）的压力下，采用Fe:卤化物化学比法，从4:1到1:4对Fe卤化物进行了计算。在结构搜索中发现的多态性的相对稳定性被描述在一个凹凸的凸包图中。凸包线连接了不同化学计量比的预测结构，这些结构在给定压力下是稳定的，不会分解为其组成元素或混合物。为了说明结果，Fe-I系统的凸包图如图2所示。通过计算得到了一系列稳定的不同化学计量的Fe-I化合物。在低压下，富含碘的二元化合物如FeI₄ 和FeI₃ 被发现是最稳定的。1:1的FeI在125 GPa时变得最稳定。在高压下，富铁化合物是首选，而Fe₂I在320 GPa时是最稳定的。在二元Fe-Cl和Fe-Br化合物的化学演化过程中，高压下稳定卤化铁的化学计量比变化顺序相同。在低压下（Fe-Cl和Fe-Br体系≤135 GPa，Fe-I体系≤70 GPa），富卤素化合物更稳定。稳定的卤化铁只能在高铁以及高压条件下形成。此外，在较高的压力下，卤化铁与组成元素（固态铁和卤素）的活化焓差异较大。这些结果表明，二元卤化铁特别是碘化物，可以稳定在地球的内核中。值得注意的是，理论结果证实，在低于100 GPa的压力下，稳定的卤化铁也可以以富卤素相的形式存在，这一条件与早期地球的条件相近。因此，描述挥发性卤化物耗竭的一种可能的情况是，卤化铁形成于积聚的早期阶段，然后随着地球的增长和压力的增加而分异到地核中。早期在低压下形成的富卤素卤化铁，可能由于过量卤素的释放而发生化学演化，以满足高压下化学稳定性的计量比要求。富铁卤化物的残骸可能仍然存在于现今地球的内核中。

《图2》

图2.（a）选定压力下具有不同化学计量比的二元铁-碘多晶型的凸包图；（b）从结构搜索研究中发现的最稳定高压相的晶体结构。

《3.3. 地幔条件下 CO2 和 SiO2 的化学行为》

3.3. 地幔条件下 CO₂ 和 SiO₂ 的化学行为

CO₂ 和SiO₂ 都是第四主族氧化物，在地球上含量丰富。在正常情况下，SiO₂ 不会与CO₂ 发生反应。众所周知，压力对化学键的性质有很大的影响。例如，Si是一种在室温下的半导体，但在11 GPa时转变为金属超导体。一个更引人注目的例子是钠，它是自由电子金属的一个典型模型，在极端压力下它会变成一个宽带隙绝缘体。更与主题相关的是，当压缩到40 GPa并加热到2000 K时，固体CO₂ 聚合成一个完全4配位的扩展固体，然后在60 GPa和4000 K时转变成金属液体。因此，在高压和高温下，可以通过改变CO₂ 和SiO₂ 的化学性质来促进化学反应。此外，在这些条件下，CO₂ 是一种超临界流体，可以渗透到多孔材料中。如果确实发生了反应，这一知识可能会对我们理解深层碳循环产生深远的影响。Santoro等[43]首先观察到，在压缩至18~26 GPa以及加热至600 K和980 K的条件下，光谱证据明确表明加载CO₂ 的微孔SiO₂ （沸石）样品会生成碳酸硅。然而，尚未有详细的结构信息或力学分析[43]。考虑到这些反应对地球科学的潜在重要性和相关性，利用恒压-恒温（constant pressure–constant temperature, NPT）MD计算方法，对1:1 CO₂ 填充沸石SSZ-56体系进行了理论研究。使用MD的动机在于它能够通过对模型系统施加不同温度下的外部压力来模拟实验，并在最终达到热力学平衡（即所有的反应都停止了）。平衡原子轨迹的分析将提供所形成物种振动光谱的有关信息，可以直接与实验红外（IR）和拉曼测量进行比较[43]。在0 GPa和300 K的初始热平衡阶段，SiO₂ 骨架与捕获的CO₂ 之间没有发生反应。当压缩到26 GPa及以上时，可以从晶体细胞轴的缩短来看出，尤其是在晶体b方向，SiO₂ 骨架开始收缩并变形。当压缩结构在26 GPa下加热至1000 K时，CO₂ 与通道表面暴露的SiO₂ 发生化学反应[图3（a）和（b）]。因此，从产品中可以鉴别出单齿、双齿和桥接碳酸盐[图3（c）]。预测结果与实验结果吻合较好。虽然继续计算到25 ps，但8 ps后没有观察到进一步的反应，暴露在外的SiO₂ 表面现在完全被与CO₂ 反应的产物所覆盖，空腔内剩余的游离的CO₂ 不再与之相互作用。该模型系统已经达到热平衡，可以利用附加的平衡原子轨迹来计算系统的动态特性。通过对模型最终结构的考察，发现反应产物主要由位于沸石大通道的3配位碳酸盐（CO₃）的聚合物[–O–(C=O)–O–]链和未反应的CO₂ 分子组成。令人惊讶的是，反应完成后沸石框架只发生了轻微扭曲。这是因为化学反应只发生在空通道表面，即使经过加压和热处理[43]，整体框架也没有受到严重影响。

《图3》

图3. 在1000 K和26 GPa条件下，NPT模拟SiO₂ -CO₂ 模型的结构，显示了单齿、双齿和桥接碳酸盐结构的存在。

原子动力学模拟的一个显著优点是，反应机理可以从对原子运动轨迹的研究中得到阐明。在这种情况下（图4），分析表明化学反应分为三个步骤。① 在高温下，CO₂ 分子受到热激发，并进行大幅度的O-C-O弯曲运动。② 通过降低模型体积来增加压力，使CO₂和SiO₂ 更接近。弯曲的CO₂ 使C=O键减弱，并与沸石通道表面的SiO₄ 单元的Si发生反应。在前沿轨道描述中[44]，CO₂ 的O原子上的孤对电子被转移到SiO₂ 的空σ ^* 轨道上（图5），同时削弱了Si–O键。随后，碳（C）与桥接Si–O–Si成键，形成CO₃ 。③ 随后腔内其他游离CO₂ 分子的连续反应产生聚合物结构。

《图4》

图4. 三种主要反应步骤示意图。

《图5》

上述理论结果说明了MD方法的通用性。除了平衡性质的研究，MD还可以用来研究非平衡化学过程，并从原子轨迹的分析中揭示反应的原子机理。

《3.4. 地球外地幔中水的形成》

3.4. 地球外地幔中水的形成

硅以SiO₂ 和硅酸盐的形式存在，是地球和其他类地行星上最丰富的物质之一。石英是最常见的晶体，在环境条件下具有化学惰性。然而，Shinozaki等[45]最近报道了SiO₂ 在与地球上地幔对应的压力和温度条件下在气态氢（H₂ ）中发生溶解。此外，原位拉曼光谱和红外光谱显示了硅烷（SiH₄ ）和水（H₂O）形成的证据。在这种热力学条件下观察到水是令人兴奋的，因为它可能对地球上水的起源有重要的意义。此外，微量水存在的可能性会影响地幔动力学的过程，引发大陆地幔岩石圈深部地震的成核与发生。为了理解这一新的现象，我们在含有1296个SiO₂ 的α石英板片上建立了一个模型，该模型的真空区域充满了3294个H₂ 分子，使气体温度达到100 K，同时压力达到1 GPa。在如此大的系统上执行第一原理DFT MD目前是不可行的。取而代之，使用了反应力场ReaxFF [46,47]，这是一种能够模拟化学变化的多体势场。ReaxFF基于电荷平衡（charge-equilibration, QEq）算法确定原子电荷，继而动态计算化学键顺序和原子价态。在这个过程中，原子间的相互作用会随着化学的和结构的变化而动态调整。模型体系被缓慢地置于1700 K和2 GPa的实验模拟条件。最初发现氢分子扩散到石英板片的亚表面区域（图6），由于石英内部通道相对较宽，H₂ 分子可以深入到石英板片的内部，最终，被捕获的H形成了半均匀的浓度分布结构[图6（b）~（d）]。在这个阶段，H保持了分子结构，但是氢键从0.76 Å的分子键长延长到了1.0~1.1 Å [图6（b）~（d）]。随后，被削弱的氢分子解离[图6（e）]成氢原子，与石英的O原子发生化学作用。这个反应削弱了硅氧键，导致了板片的内部水分子的形成。这一过程一直持续到大量水被限制在板片内[图6（e）]。反应后的石英转变为非晶态的二氧化硅。综上所述，SiO₂ /H₂ 的溶解反应经历了三个不同的阶段（图7）。首先，H₂ 侵入石英多孔通道。接下来H₂ 和共享角O原子之一以及Si原子相互反应。然后，第二个H₂ 以类似的方式与邻近SiO₄ 团簇的Si反应，导致H₂O的形成和SiO₄ 四面体键的断裂。MD计算成功地再现了实验中观察到的所有特征。这项研究支持了地幔中SiO₂ （可能还有硅酸盐）和H₂ 之间的反应可能是地球上产生水的一个本地来源的观点。

《图6》

图6. SiO₂ 模型在MD模拟不同阶段的结构。（a）初始结构；（b）初始模型平衡至100 K和1 GPa；（c）580 K和1.3 GPa的结构；（d）1220 K和1.7 GPa的结构；（e）最终结构为1700 K及2 GPa。所有结构都朝向（100）方向。Si原子用绿色表示，O原子用红色表示，H原子用浅蓝表示。

《图7》

图7. 氢与石英反应生成水。从左到右： H₂ 分子首先吸附在晶格中的Si−O对上，H原子解离，然后与Si原子上的其余H原子形成水分子。

对承压水结构的分析令人吃惊。含水区域的O−O径向分布函数g_OO (r) [图8（a）]与之前FPMD计算[48]中报道的密度为1.95 g·cm^–3 的高压水非常相似。特别是，没有一个清晰的第一配位壳层，即距离r = 1.25 Å处的O−H径向分布的第一个最小值为非零，这表明水分子发生了解离[48]。对水体结构的目测和对水体中化学物质的分析证实了这一观点。图8（b）对发现的分子种类进行了图形化总结。含水区域中OH^– 、H₂O、H³O⁺和H₄O²⁺ 的存在较为明显。这些理论结果出乎意料，可能会产生一些重大的后果。存在大量的离子种类，如H₃O + 、H₄O₂⁺ 和OH，从反常高压下解离的水不仅可以改变其化学性质，如提高反应活性，而且可提高其物理属性（例如，相比在类似条件下，它的纯配对物给它一个更高的热扩散率和导电性）；反过来，这将对从地幔交代作用到部分熔融和岩浆作用的各种地质过程产生深远的影响。

《图8》

图8.（a）SiO₂ 结构中含水层的O–O、O–H、H–H径向分布函数；（b）在封闭含水层中鉴定出的化学物质（OH^– 、H₂O、H₃O⁺ 和H₄O²⁺ 化合物）。ID：确认。

《3.5. CaCO3 熔体的黏度》

3.5. CaCO₃ 熔体的黏度

碳酸钙是岩石矿物中常见的化学成分。在自然界中，它以方解石和文石的形式存在。石灰岩、白垩土和大理石中也含有大量的碳酸盐。因此，碳酸盐岩矿物是潜在的碳库。了解地幔中富碳酸盐岩熔体的迁移规律对认识深部碳循环及相关地球化学、地球物理现象具有重要意义。矿物熔体黏度是火山活动和岩浆流动的一个重要定量参数。然而，在高压和高温下测量黏度是一项具有挑战性的工作。因此，在地震活动活跃的过渡带或更深层条件下，很少有实验测量的黏度数据。近年来，同步辐射X射线超快成像技术被用于测量CaCO₃ 熔体的黏度，温度和压力达到2000~2500 K和6.2 GPa。结果表明，碳酸盐类熔体的黏度与环境条件下的液态水的黏度相当，比硅酸钙熔体的黏度低一个数量级[49]。为了研究这种令人惊讶的低黏度的来源，并获得与地球过渡区有关的信息，FPMD计算了CaCO₃ 熔体在2000 K 和3000 K、35 GPa时的性质。黏度由Green-Kubo时间自相关函数计算得到，即恒定体积-恒定温度（constant volume-constant temperature, NVT）正则系综中MD计算得到的瞬时应力张量的非对角分量[式（13）]。相关函数法的缺点之一是在相关函数中存在长时间的持续小振荡。因此，准确估计黏度需要平衡MD模拟大于50 ps，当熔体接近液体-玻璃边界时，需要更长的模拟时间。从图9（a）中可以看出，在32.6 GPa和3000 K条件下计算出的CaCO₃ 熔体剪切应力时间自相关函数清晰地显示出长振荡尾巴。对于该体系，黏度估计为（14± 1）mPa·s。在2000 K和3000 K的选定压力下进行了类似的黏度计算，结果如图9（b）所示。在4.1 GPa和2000K时，5.8 mPa·s的理论黏度与2063 K和6.2 GPa时测得的约6 mPa·s的实验值基本一致。在2000 K大于8 GPa时，CaCO₃ 熔体变为黏性液体，自扩散系数非常低。由于需要很长时间的模拟计算，黏度的计算已不实用。另一方面，在3000 K时，CaCO₃ 熔体呈流体状，黏度可计算到32.6 GPa。在0.6 GPa和3000 K时，计算得到的黏度为2 mPa·s。这个数值比大多数硅酸盐熔体在相似的热力学条件下的数值低一个数量级，并且接近于液态水在环境条件下的黏度。然而，这种比较只是巧合，对熔体结构的分析表明，碳酸盐阴离子之间没有聚合的迹象。在低压下，CaCO₃ 熔体表现出类似理想的分子液体的性质。事实上，直接计算的黏度与斯托克斯-爱因斯坦方程计算的扩散系数基本一致。与实验趋势相似[49]，计算得到的熔体黏度在3000 K时仅在小于10 GPa时缓慢增加。一个潜在的重要预测是，当压力大于10 GPa时，黏度会突然增加，表明液体的短期顺序发生了变化。这一预测是通过对时间原子构型的可视化检查得到证实的，这表明在压力大于10 GPa时，碳酸盐岩阴离子和Ca²⁺ 之间存在短暂的相互作用。这种相互作用可以通过计算钙（Ca）与碳酸盐的C、O原子之间的速度互相关（velocity cross-correlation, VCC）来量化。在图10（a）和（b）中，我们比较了在3000 K 和几个压力下计算得到的Ca−O和Ca−C的VCC。在压力达到7.1 GPa时，短时间内快速平滑衰减的相关函数和较长时间内没有弱震荡的VCC结构都表征了通常的流体扩散行为。然而，当压力大于11.2 GPa时，VCC中开始出现100~350 fs的振荡特征。这些特征是由于相关的Ca…O和Ca…C运动。这些振荡特性有望在功率谱中形成一个独特的振动带。可由速度-时间自相关函数的傅里叶变换（功率谱）计算单粒子VDOS。事实上，在压力超过10 GPa时，一个新的波段在400~500 cm^–1由于Ca…CO₃ 振动出现在了Ca投影VDOS中[图10（c）]。振动能与VCC中的振动很好地匹配，周期约为71 fs或470 cm^–1 。可以毫不含糊地得出结论，在高压下（在2000 K大于4 GPa和在3000 K大于10 GPa），Ca²⁺ 和相互压缩，它们的运动是相关的。结果，原子扩散速率降低，熔体黏度突然增大。此外，黏度的增加影响热传输（导热系数）。这些影响可能对地球上地幔熔体的流动性产生重大影响。我们也研究了水杂质对碳酸钙熔体黏度的影响。在碳酸钙熔体中加入15%（摩尔分数）水，在30 GPa和3000 K条件下，对模型系统进行了NVT MD仿真。计算得到的黏度为10 mPa·s，比相同热力学条件下纯CaCO₃ 熔体的黏度低10%。

《图9》

图9.（a）3000 K、32.6 GPa时CaCO₃ 剪切应力自相关函数（黑色实线）、累积积分（蓝色实线）、黏度（红色虚线）；（b）2000 K和3000 K时熔融CaCO₃ 的黏度和35 GPa时含水15%的CaCO₃ 溶液的黏度随压力的变化。Exp：实验值。

《图10》

图10. Ca-C（a）和Ca-O（b）在CaCO₃ 熔体中的压力依赖性VCC函数；（c）Ca在3000 K时预测态的速度密度。

《3.6. 方镁石（MgO）的热导率》

3.6. 方镁石（MgO）的热导率

导热系数是材料的重要热物理性质。对于矿物在高压和高温下的导热性的认识关乎我们对地核热损失率和地球演化的理解。在极端条件下直接测量导热系数是一项艰巨的任务。近年来，人们做了很多努力来开发有效的计算方案来从第一原理电子计算来计算导热系数[50,51]。对于可以由多体双叠加势来充分描述的系统，利用平衡MD和Green-Kubo形式[式（14）][28]从热流密度时间相关函数中计算得到导热系数最为方便。对于成对的叠加势，热流的动能项和势能项可以分解为单个原子的贡献。然而，这种方法不能推广到基于第一原理的电子结构量子力学中，因为没有有效的方法来确定单个原子在量子力学系统中的势能。最近的研究表明[34]，如果忽略导热系数的对流项，如在固体（晶体、无序和非晶态）中，导热系数可以由单个原子的能量动量之和来估计（见式（14）[28]）。推导过程如下：导热傅里叶定律表明，材料的导热系数是热流密度与热梯度的比值

热流J是单位时间内每个粒子ε_i通过横截面A（垂直于热流）表面的能量率之和：

每个粒子对应的能量是动能和势能(u)之和：

代入式（16），热通量变为：

式（17）中的第一项为对流导热系数，第二项为传导导热系数。在固体中，对流项为零，可以忽略。固体的能量动量R可以定义为：

注意，热流J是R的时间导数，。粒子k的能量动量为：

不用Green-Kubo公式[8]对热流J的时间自相关积分来计算导热系数k，这个量可以用能量动量扩散的等效爱因斯坦关系来计算，即能量动量[52]的均方位移：

从原理上讲，瞬时能量动量[式（20）]可以由对每个粒子的时间位置、速度和合力的知识来计算。这些量在MD计算中是可用的。

通过对各种有序和无序固体[34]的研究，验证了该简化方法[式（20）和式（21）]的有效性。矿物MgO是地球下地幔中常见的组分。鉴于其重要性、简单的立方结构，以及可获得的准确的高温下单晶[53]和高压下质量较差的粉末样品上的导热系数数据[54,55]，通常选择方镁石作为测试用例来验证新的计算方案。利用NVT MD计算了MgO在0 GPa、100 GPa、300 GPa、600 GPa、900 GPa和1200 K时的导热系数。图11（a）中采用4×4×4超胞的数值计算结果，与实验数据以及之前的非谐晶格动力学理论计算结果进行了比较。采用新方法计算了MgO在300 K时的导热系数为（70.3± 8.9）W·K ^–1 ·m^–1 ，与较为分散的实测数据相吻合，其范围为36~70 W·K^–1 ·m^–1 。这一数值也符合其他理论估计的65~111 W·K^–1 ·m^–1 。应当指出，在计算和（或）测量样品中镁（Mg）和O同位素的自然分布中使用的密度泛函的缺陷被忽略了。在一定的压力下，由于Umklapp 声子-声子散射过程[56]的主导作用，导热系数应随温度的升高而降低。计算结果再现了这一趋势，导热系数从100 K时的153 W·K^–1 ·m^–1 下降到1200 K时的13.7 W·K^–1 ·m^–1 。在较宽的温度范围内，计算的导热系数绝对值与实验值吻合较好。MD方法的另一个优点是，计算导热系数相关的统计误差可以估计为不同时间起点但时间长度相同的MD轨迹的时间段与计算平均值的标准差。如图11（a）所示，温度越高，统计误差越小。这表明在高温下对构型空间的采样变得更有效。此外，对于一个足够大的模型，由于非谐大幅度原子运动明确增强了Umklapp散射过程[56]，如此MD结果在高温下将会更加准确。这与非调和晶格动力学方法是相反的，在非调和晶格动力学方法中，温度的影响只通过玻色-爱因斯坦统计量来改变声子总体，而不修正非谐力常数和谐波声子频率，两者都是在0 K时计算出来的。实际上，与实验结果和MD结果相比，非谐晶格动力学计算得到的导热系数在高温下下降得更快[52,53]。在熔点附近，高阶力常数将变得重要，谐波频率的使用不再有效。理论上，所有阶次的非谐效应都包含在MD方法中。为了测试长波声子和长程相互作用的影响，即预测的导热系数与系统尺寸的收敛性，在300 K时使用一个更大的5×5×5超胞模型进行了另一个计算。较大模型的计算导热系数为（73.3±9.2）W·K^–1 ·m^–1 ，属于较小模型（70.3±8.9）W·K^–1 ·m^–1 （见下文）的4×4×4模型的误差估计范围。因此，至少MgO测试用例中的系统大小依赖性不是太严重。

《图11》

图11. 高温环境压力（a）和高压300 K（b）下实验导热系数和理论导热系数的比较 [(a) Reproduced from Ref. [53] with permission of National Academy of Sciences, 2010; (b) reproduced from Ref. [54] with permission of Springer Nature Publishing AG, 2013 and Ref. [55] with permission of American Geophysical Union, 2014]。

直接原位实验测定材料在高温高压下的绝对导热系数是不可行的。在这些极端条件下，大部分的信息都是通过间接测量粉末样品的热扩散系数得到的，然后考虑恒压热容将其转化为热导率。热容通常是由状态方程的外推与其他已知的弹性性质相结合来估计的。近年来，人们尝试用瞬态热反射光谱法测定金刚石压腔的绝对导热系数[54]。这种方法的一个潜在缺点是必须使用（半）经验校正来校正来自样品围压的贡献。采用该方法，在室温[54]下，单晶MgO的导热系数可达60 GPa。采用相同的4×4×4超胞模型，计算了3300 K和选定压力下的导热系数。理论计算结果与图11（b）中的实验值进行了比较。该理论正确地预测了热导率随着压力增大而增大的观测趋势。在137 GPa的最高压力下，计算的导热系数为（357±41）W·K^–1 ·m^–1 ，几乎是环境压力（70.3±8.9）W·K^–1 ·m^–1 的5倍。一般情况下，预测的导热系数略高于观测值。例如，47 GPa和65 GPa下（178±25）W·K^–1 ·m^–1 和（211±30）W·K^–1 ·m^–1 的计算值分别大于140 W·K^–1 ·m^–1 和162 W·K^–1 ·m^–1 的实测值。考虑到高压下的较大的实验误差（47 GPa时误差为±20 W·K^–1 ·m^–1 ，60 GPa时误差为±30 W·K^–1 ·m^–1 ），采用新MD方法计算的理论导热系数并非不合理。系统尺寸依赖性也较小，使用较大的5×5×5超胞在65 GPa下，计算的导热系数（218±31） W·K ^–1 ·m^–1 在较小的4×4×4模型计算的误差范围（211±39） W·K^–1 ·m^–1 内。与高温情况相比，压力越大，统计误差越大。这是因为当原子相互压缩时，构型空间的热采样受到限制。为了获得更高的数值精度，需要更长的MD运行时间。

虽然计算正确地再现了导热系数的压力趋势，但其绝对值始终被高估，并且随着压力[55]的增大，与实验结果的偏差也增大。然而，MD结果与早期的晶格非谐动力学计算[53]非常相似。造成这种差异的原因不太可能是理论上的缺陷。由于实验测量的物理量为热扩散系数[55]，因此需要等压下的热容信息C_p ，将该物理量转化为导热系数。虽然在高压下C_p 不是一个可测量的量，但利用在低压下的弹性性质的信息[57]，可以通过估计已知的实验状态方程获得。把这个量外推到高压下可能是不正确的。通过将状态方程[54]得到的MgO在不同压力下的经验热容C_p 与准简谐近似（quasi-harmonic approximation, QHA）晶格动力学计算结果进行比较，证实了这一推测（图12）。结果表明，在环境压力下的估计值C_p 与计算值C_p 是一致的，但在较高压力下，状态方程（“经验”）的估计值低于准非谐理论的计算值。压力越大，这种差异就越大。因此，“实验”导热系数始终低于理论值。如果用适当的（QHA/经验）C_p 比来衡量“实验”导热系数（图12），则其与理论结果的绝对一致性将大大提高[图11（b）]。

《图12》

图12. 用不同方法得到的MgO热容量（C_p ）比较。红线：从准谐近似理论（QHA）计算得到；黑线：由实验得到的状态方程估算；蓝线：由参考文献[54]的补充材料报道的二阶多项式拟合得到。

该计算方法忽略了对流的贡献，不适用于液体。在一个完整的量子力学理论中，不可能把能量动量分成单个的原子分量。预测和观测到的大类别固体[34]的导热系数之间的非常合理的一致性可以通过以下直观的论证加以合理化：从平衡势能面出发，系统总能量的时间变化可以近似为单个原子上的扰动之和。与每个原子相关的能量变化只是位移乘以作用在该原子上的有效力，该力可通过Hellmann-Feynman定理计算得出。对于固体，式（21）是一个合理的近似。尽管存在不精确的地方，但经验证据支持了MD方法的准确性，表明该方法可以应用于高温环境，在高温环境下，高阶非谐效应是不可忽视的。

《4. 展望》

4. 展望

理论建模已经并将继续成为高压/高温研究不可缺少的工具。在本文中，只说明了小部分可能的应用。利用有限的算例，证明了用FPMD可以计算出可靠的热物理输运性质，如黏度和导热系数。该方法的优点是将温度效应和相关的非谐原子振动间接地纳入算法中。在理论上，DFT晶格动力学计算对应于非热条件（如0 K）。为了将理论结果与实验测量结果进行比较，需要用准谐理论等微扰处理对热效应进行校正[58]。只有当温度不远超固体德拜温度时，这些修正才有效。不幸的是，矿物中有许多有趣的异常现象发生在熔点附近，熔点的温度比德拜温度高得多。实验在极端温度和压力条件下存在重大挑战，无法直接测量黏度和导热系数。因此，FPMD可能是首选的计算方法。例如，该方法的一个潜在应用是研究六方体致密铁（hcp-Fe）在内核条件下接近熔化时的剪切弱化[59]。虽然这里没有提到，但弹性模量可以通过等压-等焓（P、H、N）系综中应变的波动动态计算获得[60–62]。本文还表明FPMD可以用来探索新的化学过程，并揭示目前在极端条件下的实验无法修正的原子机制。这里演示的FPMD计算需要大量的计算资源。

为了使用时间相关函数方法计算获得可靠的统计数据，需要进行长时间的模拟。就像计算导热系数一样，可以用相应的等效爱因斯坦扩散方程来减少热力学平衡后数据分析所需的时间步长。在这里描述的例子中，使用当前的计算机技术对一个热力学数据点上的几百个原子进行FPMD模拟，将需要使用100~200个核的大规模并行计算机同时计算1~2个月的时间。鉴于大规模电子结构计算在计算能力和能力方面的迅速发展，以及在更可靠和更有效的算法方面的积极发展，笔者乐观地认为，MD方法将在不久的将来得到广泛应用。

《Acknowledgements》

Acknowledgements

The author wishes to acknowledge his collaborators, Drs. N.J. English, Y. Pan, and T. Iitaka, and previous postdocs and students, Drs. X. Du, M. Wu, and X. Yong, for their contributions to the studies reported here. He also thanks WestGrid (Canada), RIKEN (Japan), and HPSTAR (China) for the allocation of computational resources.

展示更多