深物质能量与极端条件下材料化学的第一性原理方法

工程5（2019）421研究深物质能量极端条件下材料化学和输运性质研究的第一性原理方法John S.谢a，ba加拿大萨斯卡通萨斯喀彻温大学物理和工程物理系，SK S7N 5E2b高压科学技术前沿研究中心，上海201203阿提奇莱因福奥文章历史记录：2018年6月30日收到2018年8月23日修订2019年1月16日接受在线预订2019年保留字：高压高温地幔分子动力学密度泛函方法输运性质化学反应性A B S T R A C T地球是一个动态系统。地球的热力学条件因深度而异，从表面的环境温度和压力到核心的360 GPa和6600 K因此，地球成分的物理和化学性质（例如，硅酸盐和碳酸盐矿物）受其直接环境的强烈影响。在过去的30年里，在极端条件下材料表征的实验技术和理论建模方法方面都取得了巨大的进展。这些进展提高了我们对矿物性质的认识，这对于充分理解地球的形成和生命起源是必不可少的。本文综述了最近用于预测极端条件下材料行为的计算技术。本文主要介绍第一性原理分子动力学（FPMD）方法在地球科学中的应用。©2019 The Corner.Elsevier LTD代表中国工程院出版，高等教育出版社有限公司。这是一篇CC BY-NC-ND许可下的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0/）。1. 介绍具有讽刺意味的是，虽然我们对太阳系以及周围一定距离的恒星和星系了解很多，但对我们自己星球的组成及其相互作用的确切知识却相当有限。有史以来最深的钻孔只穿透地壳12公里[1];因此，原则上，我们没有直接证据表明更深处埋藏的物质及其行为。我们关于地幔和地核组成的大多数信息都是间接地从地震学、火山活动产生的挤压矿物或钻石中的矿物中获得的。 在这方面，在极端条件下材料的物理和化学的精确知识是至关重要的，以便根据这些发现构建连贯，逻辑和可靠的解释。在早期，各种岩石类型的弹性常数被用来解释密度剖面作为深度的函数。从那时起，理论和建模在地质研究中就不可或缺了[2]。基本物理学和化学原理的融合导致了许多成功和发现[3]。一个例子是化学键理论，它是在19世纪后期发展起来的，此后，许多后来的科学家对该理论进行了改进;该理论产生了预测原理和规则，可用于解释环境压力下许多矿物的晶体结构，仅需要几个关键结构参数和物理性质的信息[4，5]。可以说，在过去20年中，这个理论领域取得的最大进步是许多量现在可以从第一原理可靠地计算出来[6，7];该方法基于量子理论，只需要输入组成元素的知识就可以求解薛定谔方程。最近，甚至可以可靠地预测晶体结构计算机建模已经成为一种新的范式，它为研究材料在极端压力和温度下的行为开辟了一个新的维度，而这些压力和温度不容易通过实验进行修正。本文综述了几种计算方法，重点介绍了它们在化学反应和输运性质研究中的应用。这篇简短的评论并不打算对这一领域的最新技术进行全面的调查，也不打算收集所有相关的文章。相反，本书的目的是向一般读者提供该领域现状的概要电子邮件地址：john. usask.cahttps://doi.org/10.1016/j.eng.2019.01.0082095-8099/©2019 THE CONDITOR.由爱思唯尔有限公司代表中国工程院和高等教育出版社有限公司出版。这是一个在CC BY-NC-ND许可证下的开放获取文章（http://creativecommons.org/licenses/by-nc-nd/4.0/）。可在ScienceDirect上获得目录列表工程杂志主页：www.elsevier.com/locate/eng422J.S. 谢/工程 5 （2019）421. Σ..1.阿克萨斯交流e，K 1/2/1/2/2/3/1/2/3/1/3积分。i jijij. .- 是的- 是的- 是的半]里伊半]小时 j ji..半]半]半]2m12312312312132123.Σ材料的性质完全取决于材料之间的相互作用，原子之间的原子力一旦原子间的作用力2. ！re;！re;！re;：U. ！re;！re;！re;：2原则上，所有的属性都可以计算。例如，诸如晶体结构的静态性质对应于晶胞中原子的最低能量排列;刚性是系统对外部应力的响应动态这就是系统的薛定谔方程。除了氢原子，电子波的确切形式功能W. ！re;！re;！re;：是未知的，并且对于在有限温度下的输运性质也可以被计算。使用经典分子动力学（MD）结合统计力学理论[8]。例如，热导率是能量通过材料的传输，而液体的剪切粘度是由于与应变率相关的变形而引起的内力引起的弹性（粘性）应力。通过数值模拟计算材料特性并不是什么新鲜事。事实上，自从MD方法在20世纪60年代被构思和实施以来，几乎所有的理论形式主义都已经建立[9 因此，一旦原子间的作用力通过将原子间势能函数拟合到实验观测值（经验的）或原子或分子簇的高级量子化学计算（非经验的）来定义系统，可以预测各种性质。例如，已经构建了几个非常成功的Si-O、O-O和Si-Si相互作用的原子间势多组分体系的势函数也是可用的，并提供了许多硅酸盐和铝硅酸盐固体和熔体的结构和动力学的有用信息具有可靠的原子间势的经典MD计算仍然是使用当今计算机技术的强大预测工具;可以对数百万个原子执行10至100纳秒量级的持续时间的模拟。研究缓慢弛豫的液体和玻璃通常需要很长的模拟时间。预定义的原子间势的一个主要缺点是它只能应用于参数化范围内的环境。因此，大多数潜在参数是不可传递的。当然，这会影响预测的可靠性。因此，矿物在地球条件下在广泛的温度和压力范围内的行为为了克服这一困难，必须求助于无可调参数的第一性原理量子力学方法。这就是这次小型审查的主题。2. 理论背景根据量子力学的基本假设[7，17]，多粒子系统（电子、原子核等）是唯一地定义通过一波功能W. ！re;！re;多粒子系统一种方法是假设每个电子在系统中只经历其他电子的平均场，并且每个电子都有自己独特的波函数f！Ri（独立粒子近似）然后，多电子全波函数W！re;！re;！re;：/！r1/！r2/！r3：，（由于泡利不相容而反对称原理）可以通过单个和独立轨道波函数的乘积的行列式来近似将总波函数代入Eq.（2）利用变分原理，使总能量相对于电子轨道最小这个过程导致一组耦合的微分方程，可以用迭代的方式求解这种方法被称为库仑上的二选r，Jij¼. /i1/j2。 1. /i1/j2，以及里伊1964年，霍恩伯格和科恩提出了一种完全基于电子密度而不明确调用波函数的电子结构理论的不同表述[18]。密度泛函理论（DFT）证明了电子密度是主要变量，并且相互作用的多电子系统的总能量是电子密度的函数。后来，Sham和Kohn[19]提出了一个计算电子系统总能量密度泛函理论消除了计算交换积分的负担，并用一个密度函数的通用泛函代替了电子交换和相关相互作用，尽管这个泛函的确切形式是未知的。系统总能量的Kohn-ShamEehn。！t. ！rn。！rd！rF½n]3其中n是电子密度。R. 函数alFn是电子-电子和电子-核相互作用的势能与电子动能的和。再一次，F n的确切形式未知。很容易将电子库仑相互作用的贡献从总电子相互作用中分离出来。可以进一步分解为以下：！re;：;！re;！rN;！rN;！rN;：;！rN，它显式地依赖于. ！- 是的！0Σ3n1 2 3MZZnRn r电子和原子核的坐标。的总能量1..！啊！0ijF½n]¼2.！啊！0。dr dr[2019-04 - 24]系统是能量哈密顿量的期望值EWH^W。总能量哈密顿量，简单地说就是动能和势能。如果电子和核的运动（2R-R其中G n仍然是一个泛函，包含a的动能T n，系统的非相互作用的电子密度n[19]，和剩余电子相互作用术语分为Excn，这被称为好的。！re;！re;！re;：12mG½n]<$T½n]Exc½n]5J.S. 谢/工程 5 （2019）421423我R其中m是电子的质量，势能U是总和electron–electron, electron–nuclei and库仑相互作用。 将动量算符替换为根据变分原理，基态总能量是通过最小化总能量来获得的，电子密度受到粒子数，量子表示，变成：p^！呃，的能源哈密尔顿·尼安系统中的N是守恒的，即，424J.S. 谢/工程 5 （2019）421Eq的微分 （3）产量V，n. ！rd！r1/4N。因此，J.S. 谢/工程 5 （2019）421425. - 是的Σ1/1我¼ ðÞR.半]！12019年12月28日ðÞðÞi我Pr1我...ð Þð Þð ÞðÞð@@@@eijVkBTXYXY.CUP. . - 是的2iiZn. ！r≤8V.！rZn.！rdExchn.！鲁伊dT½n]l9=d！r10的一旦选择了泛函和基集[25]，dext：1995年。！r-！r0. þdn. ！r被遗弃的; ！r-;¼隔离或扩展系统的结构和能量可以..ð6Þ这里，l是电子密度为n的均匀气体的化学势。 XC电位定义为V XC！R1/4dEx c½n]=dn！R.由于相互作用电子系统的精确动能泛函是未知的，因此必须将多体波函数的动能近似为单电子的乘积。波函数，作为在的Hartree–Fock计算。DFT在电子结构计算中的应用自20世纪80年代以来，在原子、分子和凝聚态物质方面取得了巨大的成功一旦确定了作用在原子和模型超原胞上的力，就可以根据相空间中作为时间函数的玻尔兹曼统计对原子的轨迹进行采样在任何选定的温度和压力下，牛顿运动方程[26]，Fitmi½dvit=dt]，其中mi是原子质量，vi是时间速度。有用n！R¼N u！R2. 相互作用电子ic系统中的每个电子-1r2 ！你好。！我是你。！r7其中有效电势Veff. ！r的定义如下：从MD计算获得的结果是时间位置，每个原子的运动速度和作用力根据这些量，可以计算平衡和热输运性质[8，26]。平衡热力学变量（例如，A）是根据统计理论应用遍历假设计算的，即热力学平均值可以近似为计算量的时间平均值（b1=kBT，kB是玻尔兹曼常数）。【二十七】：.！Σ.！Zþn.！r.！ΣhAi¼RdCAtexp-bHð10 ÞVeff R1/4V输出V.！R！r0.Vxcrð8ÞdCexpression-bHexpression热输运性质也可以从等式计算为了将KS DFT应用于实际问题，必须对XC泛函EXCn进行近似。第一近似是局域密度近似（LDA），其中XC势[19]在非均匀真实系统中的每一点处，被设置为在相同电子密度下通过量子蒙特卡罗计算的均匀电子气的交换相关势[20]。简单的LDA近似忽略了电子密度的不均匀性，因此，倾向于高估接近平衡的化学键的强度，低估原子间的距离。通过包括电子密度梯度（广义梯度近似，GGA）[21]来考虑电子密度的曲率，基于Green和Kubo开创的时间相关形式主义的library MD模拟[8本质上，该理论假设输运过程的动力学是随机的（马尔可夫随机过程（A）的耗散率的一腐烂时间相关性函数C t Aωt A0（如果A是复动力学量）。例如，扩散速率（D）是作为时间的函数的原子（扩散）的均方（位置）位移可以表明，随着时间t，速度自相关的积分给出了扩散系数[8]：预测的结构或有时总能量，但不能两者兼而有之。这个问题与GGA不满足半局部泛函的所有已知约束的事实相关联。所有人-.h！R 不！r0i26Dt2tZtD！v0t！v0Edt11为了解决这个问题，引入动能密度S1=2OCCJUIJ2，以解释金属键的电子密度的缓慢变化和共价键的更紧密的密度这种类型的功能被称为元GGA。最近，构造了一个由于轨道ui已经计算，将非相互作用动能项包括到半局域交换相关能中与所有其他可用的泛函相比，SCAN泛函已被证明具有优越的几何形状和总能量[23]。它是迄今为止最有效和可靠的主族化合物的功能。DFT计算的一个关键组成部分是评估作用在原子上的力。这些信息对于几何优化和MD模拟是必不可少的。的受力情况一个原子，Fi，是能量对通过类比，可以为粘度构造时间相关函数。如上所述，粘度是液体对弹性变形的响应。 时间原子的位置，！变形模型系统中的！首先，通过变形矩阵H的线性变换，！别这样！r<$t <$;h从参考系变形了无穷小应变e，h<$I<$e，I是单位矩阵。在t1/40处，应变被切断，并且系统将经历剪切速率的d函数尖峰，@v x=@y1/4-edtt=：以下[8]：.Px yt-eV kBT.Pxy0Pxyt12应力张量（r）描述了在以下条件下总能量的变化一个无穷小失真的的基向量到位移kj，并且可以通过Hellman-Feynman定理方便地进行评估h0k1/4Pj. dijeijh下的应变eij。 应用赫尔曼-E HFFeynman theorem（费曼定理）（9）），可以对应力分量进行估计。日期：Pij¼-hWj@HhkjWi. 稳态应力（E。例如，在一个实施例中，Pxy）i¼ -@k¼ -@khWjHjWi¼-hWj@kjWi由稳定剪切产生的只是时间上的积分，j j j和剪切粘度g. ^rxy=ex由以下公式给出[8]：其中繁琐的对波的导数的积分函数被简化为更易于管理的表达式，能量哈密顿量的导数。g¼ 1Z1d t.PtP0130¼¼Þ-426J.S. 谢/工程 5 （2019）421PDtN第1页！JJ3VkBT20另一个对地球科学有用的输运性质是热导率j。热导率是累积的矢量热通量！单位面积和单位时间内的能量动量转移率！我不知道他会不会杀了我。的表达式压力将这个强大的概念扩展到无序玻璃并不简单。然而，根据某些可观察的参数对高压液体和玻璃的性质进行分类是很有吸引力的。已尝试最近，为了关联各种氧化物的体性质，热流自相关函数变为[8，28]：j¼1Z1dtD！我的天！日本语简体中文繁体中文注意，计算扩散率、粘度和热导率所需的信息是时间位置、模型单元上的剪切应力和粒子能量。 所有这些量都可以从MD计算中获得。时间相关函数是一种强大的技术，其中许多动力学性质，如振动态密度（VDOS，红外（IR）和拉曼光谱），可以从MD计算的原子轨迹计算[26]。在下面的讨论中，我们将举例说明这些方法在极端条件下对地球科学中几个选定问题的应用。3. 说明性实例在这一节中，说明了上述方法在几个地质相关问题中的应用。首先，它将表明，传统的刚性球的概念，确定的原子大小不再适用于高压。在玻璃结构中，在任何给定的压力下，不存在唯一的Si-O或Ge-O配位。因此，任何试图将配位数（CN）与玻璃性质联系起来的尝试都是无效的。作者结果表明，富铁卤化物可以存在于地球的核心;这表明，重卤素可能被隔离到核心，提供了一个合理的解释地球挥发性卤素的耗尽。第一性原理分子动力学（FPMD）计算用于研究二氧化碳（CO2）和结晶二氧化硅（SiO2）之间可能的反应[31]。结果表明，在较低的压力和温度下，CO2会与SiO2发生化学反应。然而，如果液态CO2与界面接触，在过渡区的压力和温度条件下，预期不会发生反应。 在另一项研究中，我们发现H2可以很容易地与石英反应在低压和温和的温度下形成水[32]。其余的讨论是专门计算的运输性能。我们研究了纯碳酸钙（CaCO3）熔体的粘度和少量水存在下CaCO3熔体的粘度[33]，发现水合熔体的粘度低于纯液体。我们还讨论了一种新的有效近似计算有序和无序固体在高压和高温下的热导率[34]。3.1. SiO2和GeO2玻璃的结构由于液体矿物的性质在高温高压下很难测量，淬火玻璃常被用作解释液体结构的第一近似值如前所述，原子大小和原子价的概念最初由Pauling[4]和O'keeffee和Navrotsky [5]引入，已成为阐明环境温度下各种矿物晶体结构的非常有力的工具。[35]第三十五章：我的女人这个想法似乎并不不合理，因为在高压下，氧化物玻璃中体积较大的氧（O）原子（氧化物阴离子）彼此压得很近，这改变了局部环境，因此应该会影响某些体积性质，如压缩性和粘度。困难在于分配唯一的O离子半径，这是计算高压玻璃中的氧填充分数所需要的。为此，有人认为离子半径与阳离子-氧CN和键长有关，并且可以从一组经验关系式中导出[35]。由于局部几何形状（即，CN）随压力变化，因此在不同的压力状态下必须使用不同的公式为了应用该方法，有必要先验地选择局部几何形状（即，4-、5-或6-配位）。这就提出了在给定压力下在氧化物玻璃中分配单个CN的有效性的问题此外，对于给定的键类型，在宽的压力范围内假设恒定的原子半径比是没有道理的，因为单个原子的填料），并应不同于在环境压力下。为了解决这个问题，使用FPMD计算[26]获得的理论结构分析了不同压力下SiO2[36]和二氧化锗（GeO2）[29]玻璃的结构和相关几何参数。玻璃中CN的定义有些随意，因为它取决于截止分离的选择。在研究中，选择Si-O和Ge-O径向分布函数中第一极小值的距离GeO2的结果简要地总结在图2和图3中。 1和2. 对于SiO2玻璃也发现了类似的显然，除了在最低压力下，没有单一的唯一CN可以归因于玻璃结构。该结构总是由氯化萘的混合物组成 图图1（a）显示了不同氯化萘在GeO 2玻璃中的比例与压力的关系。从0到15GPa，随着5和6倍配位的增加，四倍配位稳定下降。在25-大于7-配位的CN的量随着压力的进一步增加而增加这些理论结果表明，假设玻璃在任何压力下都有唯一的CN是没有根据的下一个Ge（Si）和O原子的原子尺寸的变化是否与压力无关当Ge（Si）-O距离随压力变化时，Ge（Si）/O原子尺寸的比率是否仍然是常数即使在结晶固体中，由于多个不同的键长，延伸固体中原子大小的分配也是不明确的。为了消除任意性，吉布斯等人。[37]提出根据Bader[38]提出的分子中原子（AIM）量子理论，通过对电子密度分布的数学分析计算原子大小。根据AIM理论，分子或固体中的原子区域原子的巴德体积是由包围原子的零电子通量表面 这个基于量子力学理论的严格定义消除了在极端压力下将球半径分配给玻璃态固体的严重扭曲的原子密度的任意性和模糊性。例如，这种划分原子体积的方法已被证明成功地解释了金属玻璃中电子性质和结构转变的结构变化[39]。图2显示了锗（Ge）和O在GeO 2玻璃中的巴德体积随压力变化的比较。 1（b）. 正如预期的那样，“原子”的体积随着压力的增加而0 GPa时OBader体积的平均值为J.S. 谢/工程 5 （2019）421427Fig. 1. (a)GeO 2玻璃中不同Ge-O配位的百分比贡献的压力依赖性图2. (a)在选定的压力下具有不同化学计量的二元铁碘化物多晶型物的稳定性的凸壳图;（b）从结构搜索中揭示的最稳定的高压相的晶体结构。12.5在120 GPa时，下降到7.5 × 10 -3，下降了40%。相比之下，GeBader体积的变化要小得多：从0 GPa时的6.2 μ3到160 GPa时的5.2μ3，仅减少了16%最重要的是，巴德体积的比率随压力而变化，并且不是常数。理论结果反驳了氧填充分数假说.3.2. 失踪的卤素悖论挥发性元 素的贫化与星云 的挥发性有关 。碳质陨石（英 语：Carbonaceous meteorite）是一组原始陨石，其元素分布与星云的组成非常相似。地球中元素的浓度与C1-C10的丰度的比较为地球早期的吸积提供了有价值的信息。令人惊讶的是，地球上重卤素（溴（Br）和碘（I））的丰度被发现比它们在CI陨石中的丰度低一个数量级。这种缺乏与地球的形成和成分演变有关，并可能影响当今海洋的盐度。对于这种“丢失卤素悖论”的假设为了限制地核中卤素的浓度，重要的是要知道在地幔和地核的压力下是否能形成稳定的为此，本文研究了铁卤化物已经使用基于智能的结构搜索技术进行了探索。这种方法采用第一原理DFT计算结合结构搜索技术[40，41]，并且是一种非常强大的策略，可以以无偏的方式定位最稳定的晶体结构，只需在给定的热力学条件下了解组成原子。这种计算已被证明提供了新的洞察力的晶体类型和化学键的材料在压力条件下无法验证实验。搜索过程如下：随机生成具有不同Fe/卤素化学计量的初始结构集。试验结构首先使用DFT计算进行优化，然后使用粒子群优化算法进行改进[42]。选择具有更好能量的那些结构作为下一个优化周期的试验结构。重复该过程，直到没有发现新的能量上有利的结构 在铁卤化物的情况下，这种方法成功地再现了已知的环境结构的FeCl 3（菱铁矿），FeBr 3，FeCl 2（水铝榴石），FeBr 2，和FeI 2。为了研究高压多晶型体的结构和稳定性，在过渡区（25 GPa）、幔核边界（70 GPa）和深核（135和360 GPa）的压力下，对Fe：卤化物的化学计量比为4：1 ~1：4的Fe卤化物进行了计算。在结构搜索中发现的多晶型物的相对稳定性描绘在凸壳图中。凸包线连接了不同化学计量的预测结构，这些结构对分解成其组分是稳定的。428J.S. 谢/工程 5 （2019）421≤≤元素或化合物的混合物。为了说明结果，Fe-I系统的凸包图通过计算发现了一系列稳定的不同化学计量比的铁碘化合物。在低压下，富碘二元化合物，如FeI4和FeI3，被发现是最稳定的。1：1的FeI在125GPa时最稳定。在高压下，富Fe化合物是优选的，并且Fe2I在320GPa下最稳定。高压下稳定铁卤化物的化学计量变化规律与二元Fe-Cl和Fe-Br化合物的化学演化规律在低压下（Fe-Cl和Fe-Br系统为135 GPa，Fe-I系统为70 GPa），富卤素化合物更稳定。稳定的铁卤化物仅在高压下富铁条件下形成此外，铁卤化物和组成元素，固体铁和卤素之间的焓差，被发现在较高的压力下较大这些结果表明，二元铁卤化物，特别是碘化物，可以稳定在地球值得注意的是，理论结果证实，稳定的铁卤化物也可以在压力小于100 GPa的条件下以富卤素相的形式存在，这与早期地球的条件接近。因此，描述挥发性卤化物耗尽的一个合理的假设是，铁卤化物在吸积的早期阶段形成，然后随着地球的生长和压力的增加而被隔离到在早期阶段在低压下形成的富卤素Fe-卤化物可能经历了化学演化，释放出过量的卤素以满足高压下化学稳定性富铁卤化物的残留物存在于现今地球3.3. 地幔条件下CO2和SiO2的化学行为CO2和SiO2都是IV族氧化物，在地球上含量丰富 在正常条件下，SiO2不会与CO2发生反应.众所周知，化学键的性质受压力的强烈影响例如，Si在室温下是半导体，但在11 GPa下转变为金属超导体！一个更戏剧性的例子是钠，一种自由电子金属的原型模型，在极端压力下成为一个大的带隙绝缘体更相关的是，当压缩到40 GPa并加热到2000 K时，固体CO2聚合成完全4配位的扩展固体，然后在20 ℃下转化为金属液体。60 GPa和4000 K。因此，在高压和高温下，可以改变CO2和SiO2的化学性质，以促进化学反应。此外，在这些条件下，CO2是一种超临界流体，可以渗透到多孔材料中.如果反应确实发生，这一知识可能会对我们对深层碳循环的理解Santoro等人[43]第一次观察到碳酸硅的形成从CO2负载将微孔SiO2（沸石）样品压缩到18-26 GPa并加热到600和980 K。光谱学证据明确地表明发生了反应。然而，没有提供详细的结构信息或机理分析[43]。考虑到这些反应在地球科学中的潜在重要性和相关性，采用等温等压系综（恒压恒温，NPT）分子动力学方法，对1：1CO2填充的SSZ-56沸石模型体系进行了理论研究.使用MD的动机在于它能够通过在不同温度下向模型系统施加外部压力来模拟实验，并在最终达到热力学平衡之前监测CO2和SiO2之间的相互作用（即，所有的反应都停止了）。对平衡原子轨迹的分析将提供所形成物质的振动光谱信息，这些信息可以直接与实验IR和拉曼测量进行比较[43]。在0GPa和300K的初始热平衡期，未观察到SiO2骨架与被捕获的CO2当压缩到26 GPa及以上时，SiO2骨架开始收缩并变得扭曲，如晶胞轴的缩短所示，特别是在结晶学b方向上。 当压缩结构在26 GPa下加热到1000 K时，CO2和通道表面暴露的SiO2之间开始发生化学反应（图1和图2）。3（a）和（b））。因此，可以从产物中鉴定出单齿、双齿和桥接碳酸盐（图3）。预测的物种是一致的建议来自实验振动光谱分析。在8ps后，没有观察到进一步的反应，尽管计算继续到25ps。暴露的SiO2表面完全被与CO2反应的产物所覆盖，腔中剩余的自由CO2不再与之相互作用，模型系统达到了热平衡，可以用附加的平衡原子轨道计算系统的动力学性质.对模型最终结构的考察表明，反应产物主要由位于沸石大孔道中的3-配位碳酸酯（CO 3）的聚合物（-O令人惊讶的是，沸石骨架在反应完成后仅轻微变形。这是因为化学反应仅发生在空通道的表面上，即使在加压和热处理后，整体框架也不会受到严重影响[43]。原子动态模拟的一个明显优点是，反应机制可以从原子轨迹的检查中解开。在这种情况下（图4），分析表明化学反应分三步进行。步骤i：高温下，CO2分子受到热激发，进行大幅度的第二步：减少图三.在1000 K和26 GPa下NPT模拟后的SiO2-CO2模型的结构J.S. 谢/工程 5 （2019）421429图4. 三个主要反应步骤的方案。图五.前线轨道相互作用图描绘了孤对电子从CO 2插入空的Si-O r * 轨道，引发CO 2和SiO 2之间的化学反应。经学术出版社许可，复制自参考文献[44]，©1971。模型的压力体积使CO2和SiO2更接近。弯曲的CO2削弱了C@ O键，并引发了与形成沸石通道表面的SiO4单元的Si在前线轨道描述[44]中，CO2的O原子上的孤对电子转移到SiO2的空r *轨道（图5），并伴随着削弱随后，碳（C）与桥接Si-O-Si键合第三步：与空腔中的其他游离CO2分子发生连续反应，生成聚合物结构。上述理论结果说明了MD方法的通用性分子动力学除了研究平衡性质外，还可用于研究非平衡化学过程，并从原子轨道分析中揭示反应的原子机制3.4. 地球外地幔中的水形成SiO2以SiO2和硅酸盐的形式存在，是地球和其他类地行星上最丰富的物质之一石英是最常见的晶体形式，在环境条件下是化学惰性的然而，Shinozaki等人[45]最近报道了在相对温和的压力和温度条件下，SiO2在气态氢（H2）地球的上地幔。此外，原位拉曼光谱和红外光谱显示了硅烷（SiH4）和水（H2O）形成的证据.在这种热力学条件下对水的观测是令人兴奋的，因为它可能对地球上水的起源有重要的意义。此外，微量水存在的可能性可能会影响触发和孕育大陆地幔岩石圈深震的地幔动力学。为了理解这一新现象，在100 K下，对由含有1296SiO2的a-石英板和填充有3294H2分子的空区域目前在这样一个大的系统上进行第一性原理DFT MD是不可行的。相反，使用了反作用力场ReaxFF[46，47]，一种能够模拟化学变化的多体势。ReaxFF根据以下参数通过电荷平衡（QEq）算法确定原子电荷原子间的相互作用在反应过程中根据化学和结构的变化进行动态调整。模型系统在1700 K和2GPa下缓慢地进入实验条件。最初，发现氢分子扩散到石英板的次表面区域（图6）。由于石英中的通道相对宽敞，H2分子可以渗透到板的深处;最终，图第六章在MD模拟的不同阶段的SiO2平板模型的快照：（a）初始结构：（b）初始模型平衡到100 K和1GPa;(c)在580 K和1.3 GPa下的结构;（d）在1220 K和1.7 GPa下的结构;（e）在1700 K和2 GPa下的最终结构。所有结构都朝向（100）方向。Si原子显示为绿色，O为红色，H为浅蓝色。430J.S. 谢/工程 5 （2019）421·ð Þ42捕 获 的 氢 导 致 半 均 匀 的 浓 度 分 布 （ 图 1 和 2 ） 。 （ 6 ） （ b ） -（d））。在此阶段，氢保持其分子结构，但H-H键从0.76埃的分子键长延长至1.0-1.1埃（图1和2）。6（b）-（d））。随后，被削弱的氢分子解离（图1）。 6（e））转化为原子氢，并与石英的O原子发生化学相互作用。该反应削弱了水分子的分解[48]。对水结构的目视检查和对水域中存在的化学物质的分析证实了这一建议。发现的分子种类在图中以图形方式总结。 8（b）. 在水区域可以清楚地看到OH-、H2 O、H3 O +和H4 O2+的存在这些理论结果是出乎意料的，可能会有一些重大的后果。由于水的离解而产生大量的离子种类，如H3O+、H2O+板。这个过程一直持续到一个水囊被限制在板内（图6（e））。反应后的石英转变为非晶态SiOxHy。总之，SiO2/H2溶解反应分三个不同的阶段进行（图7）。首先，H2渗透到石英的多孔通道中.接下来，H2与其中一个角共享O原子和一个Si原子相互作用。反应导致H2解离，形成然后，第二个H2以类似的方式与相邻的SiO4分子中的Si反应，导致H2O的形成和SiO4四面体键的断裂。MD计算成功地完全再现了实验中观察到的所有特征。 这项研究支持了这样一种观点，即地幔中的SiO2（也许还有硅酸盐）和H2之间的反应可能是地球上产生水的一个可能的本土来源对承压水结构的分析揭示了一个令人惊讶的现象。水区域中的特别是，没有明确的第一配位壳层--也就是说，在距离r = 1.25 μ m处的O-H径向分布中的第一个最小值见图7。氢与石英反应生成水。从左到右：一个H2分子首先吸附到晶格中的在异常高的围压下不仅可以改变其化学性质，例如增强的反应性，而且可以改变其物理性质（例如，使其比类似条件下的纯对应物具有更高的热扩散率和电导率）;反过来，这将对从地幔交代作用到部分熔融和岩浆作用的各种地质过程产生深远的影响3.5. 碳酸钙熔体的粘度CaCO3是岩石矿物中常见的化学成分.在自然界中，它以方解石和文石的纯净形式存在。在石灰石、白垩和大理石中也发现了大量的碳酸盐因此，碳酸盐矿物是潜在的碳库。了解地幔矿物熔体的粘度是定量研究火山作用和岩浆流动的一个重要参数然而，在高压和高温下测量粘度具有挑战性。因此，在过渡带或更深的条件下，这是一个深震活跃区，很少有实验粘度已知。最近，利用同步辐射X射线超快照相技术测量了CaCO 3熔体在2000-2500 K的粘度，粘度为6.2GPa。结果发现，碳酸盐熔体的粘度与环境条件下的液态水相当，并且比硅酸钙熔体的粘度低约一个数量级[49]。为了研究这一令人惊讶的低粘度的起源，并获得与地球过渡区相关的信息粘度由格林-（13））从MD计算在恒定相关函数法的一个缺点是相关函数在长时间内会出现持续的小振荡一个准确的见图8。(a)SiO2板中水层的O-O、O-H和H-H径向分布函数;（b）承压水层中确定的ID：已识别。J.S. 谢/工程 5 （2019）421431··hðÞðÞðÞihi粘度的估计需要大于50 ps的平衡MD模拟。当熔体接近液态玻璃边界时，需要更长的在32.6GPa和3000K下计算的CaCO3熔体的剪切应力时间自相关函数，如图9（a）所示，清楚地显示了长的振荡尾。该体系的粘度为（14 ± 1）mPa·s。在2000和3000 K的选定压力下进行类似的粘度计算。结果总结在图9（b）中。在在 2063K 和 6.2GPa 下 测 得 的 粘 度 值 约 为 6mPa·s ， 理 论 粘 度 为5.8mPa·s，与实验值符合较在2000K下，当压力超过8GPa时，CaCO3熔体变为粘性液体，自扩散系数很低由于需要很长的模拟时间，因此计算粘度不再实用。另一方面，CaCO3熔体在3000K下是流体状的，并且粘度可以计算到32.6GPa。在0.6GPa和3000K下，计算的粘度为2mPa·s。这个值是一个数量级低于大多数硅酸盐熔体在类似的热力学条件下，顺便说一句，是接近液态水在环境条件下的粘度。然而，这种比较是偶然的。熔体的结构分析表明，碳酸根阴离子之间没有聚合的迹象。在低压下，CaCO3熔体表现得像理想的分子液体实际上，直接计算的粘度与使用计算的扩散系数的Stoke-Einstein方程的估计值合理一致与实验趋势类似[49]，在3000 K下计算的熔体粘度仅在小于10 GPa时缓慢增加一个潜在的重要预测是粘度在大于10 GPa的压力下突然增加，表明液体中短程有序的变化这一预测得到了原子时间构型的直观检验的证实，结果表明，当压力大于10GPa时，碳酸根阴离子与Ca2+离子之间存在短暂的相互作用这种相互作用可以从钙（Ca）与碳酸盐的C和O原子之间的速度互相关（VCC）的计算（VCa0VCt和VCa0VOt）来量化。 计算的Ca-O和Ca-C VCC在3000 K在几个压力下进行了比较图。 10（a）和（b）。在高达7.1 GPa的压力下，相关函数在短时间内的快速和平滑的衰减，以及VCC曲线在长时间内不存在弱振荡，表明了正常的扩散流体行为。然而，在压力大于11.2 GPa时，VCC中开始出现100和350 fs之间的振荡特征这些特点是由于相关案例. O和Ca…. C动议。 这些振荡特征图第九章（a）在3000 K和32.6 GPa下CaCO3的剪切应力自相关函数（黑色实线）、累积积分（蓝色实线）和粘度（红色虚线）(b)在2000和3000 K下熔融CaCO3的粘度和在30 GPa下15%水-CaCO3溶液的粘度的压力依赖性。实验：实验值图10. CaCO 3熔体中（a）Ca-C和（b）Ca-O的压力相关VCC函数（c）Ca在3000 K时的投影速度态密度432J.S. 谢/工程 5 （2019）421eN“X我X1型 我N3× ×2吨！12tVkTB¼-··Nlationfunction，tx····！你好！v100米宽的t. 事实上，一个新的乐队在1ixt0vieidtuijrij千分之四LLLð2我2预计将在功率谱中形成一个独特的振动带。 单粒子VDOS可以从傅里叶变换（功率谱）的速度时间自相关-代入Eq。（16），热通量可以写为：研发与工程1“XN1天.编号400 - 500 cm-1，由于Ca.. . CO 3振动出现在Ca投影的VDOS中（图1）。 10（c）），压力大于10 GPa。的振动能量与VCC中的振荡很好地匹配，具有约71 fs或470 cm-1的周期。可以得出结论我我1NNveV2IXjRij。Fijvij#ð18 Þ毫无疑问，在高压下（2000 K> 4GPa和3000 K> 10 GPa），Ca2+和CO2-相互压缩，它们的运动是相关的。结果，原子扩散速率降低，熔体粘度降低突然增加此外，粘度的增加影响热传递（导热性）。这些效应可能对地球上地幔熔体的流动性产生研究了水杂质对CaCO3熔体粘度的影响.在30GPa和3000K下，对CaCO3熔体中水的摩尔分数为15%的模型体系进行了NVT计算的粘度为10mPa·s，在Eq中的第一项（17）是由于对流而导致的导热率第二项是由于传导引起的热导率。在固体中，对流项为零，可以忽略不计。固体的能量动量R可以定义如下：R¼Xeiri 19我注意，热流J是R的时间导数，即，J¼dR= dt。粒子k的能量动量为：比相同温度下纯CaCO3熔体低10%动态条件RXrZlFvdt20L03.6. 方镁石（MgO）导热系数是材料的重要热物理性质。矿物在高温高压下的热导率与我们对岩心热损失速率和热演化的认识有关。地球在极端条件下直接测量导热系数是一项艰巨的任务。近年来，已经做出了许多努力来开发有效的计算方案，以根据第一原理电子计算来计算热导率[50，51]。对于可以通过多体成对加性势充分描述的系统，热导率最方便地从热通量时间相关函数使用平衡MD采用Green-Kubo形式主义（Eq. （14）[28]。对于成对加性势，热通量的动能和势能项可以被划分为每个原子的贡献。然而，这种方法不能扩展到基于第一原理电子结构的分子动力学，因为没有有效的方法来定义量子力学系统中单个原子的势能。最近[34]，研究表明，如果忽略热传导率的对流项，例如在固体（晶体，无序和无定形）中，热导率可以从单个原子的能量动量之和估计（见等式10）。（14）[28]）。热传导的傅立叶定律指出，材料的导热系数k是传递的热通量q与热梯度的比率（oT=ox）：不是使用第1页lim1D½Rt-