大规模、量子精度的分子动力学模拟:以极端条件液态铁为例*

曾启昱 陈博 康冬冬† 戴佳钰‡

1) (国防科技大学理学院,长沙 410073)

2) (国防科技大学,极端条件物理及应用湖南省重点实验室,长沙 410073)

液态铁作为类地行星内核的主要组成成分,其在高温高压条件下的热力学、输运及动力学性质研究,对理解行星演化有着重要意义.极端条件物态物性在实验条件下产生困难且诊断手段有限,而理论模拟在动力学、输运性质计算方面面临着规模、精度的双重要求,极大限制了这方面的有效进展.本文结合深度学习技术,通过神经网络构造液态铁的高维相互作用势能面,在保证第一性原理计算精度的前提下,将计算规模从数百原子扩展到数十万原子体系.研究了从常压到核幔边界条件下液态铁的动力学及输运性质,并与X 射线衍射、非弹性X 射线散射实验对比,二者的一致性指出,深度学习技术与分子模拟的结合为我们高通量研究极端条件下真实体系的物态物性及动力学提供了有效手段.

1 引言

类地行星往往由致密的金属内核与硅酸盐、氧化物的地幔构成.以地球为例,地核分为液态外核与固态内核两个部分,其中铁占据了主要成分,而镍的占比只有5%–15%[1].地球的液态外核占据了地核总体积的96%,跨越内部2900–5150 km 深度的范围.随着深度增加,压强范围跨越135–330 GPa,温度范围在3500–6500 K[2,3].对地核的研究,有利于理解行星内部演化的诸多问题,包括地核余热的释放、地核的热对流与质量输运、行星的磁场产生等.因此,作为地核主要组成成分,铁及其掺杂体系在极端条件下的物态物性[4-8],一直以来是行星物理密切关心的重要主题之一.

在实验室条件下高温高压状态的产生,可以利用激光、离子束、超高速撞击等方式使物质状态沿着不同的热力学路径(等容线、冲击绝热线或准等熵线)演化[9].由于高能量密度物理的特性,极端条件物态往往是瞬态,其诊断严重依赖于时间分辨的实验技术.完备的热力学状态方程、静态结构、动态结构及输运性质测量,依然是实验技术前沿.因此从量子理论出发的数值模拟一直以来是不可或缺的重要研究手段,但是受限于效率-精度两难困境.

精度方面,是以密度泛函理论(density functional theory,DFT)为代表的第一性原理分子动力学(abinitiomolecular dynamics,AIMD),其核心思想是求解多电子的薛定谔方程.在获得了电子对离子的平均场之后,可以有效计算原子的受力、应力张量,从而获得离子系统的统计、关联信息.在Kohn-Sham DFT 的框架中,多电子系统的基态电荷密度等价于给定外场下单电子系统的基态电荷密度,通过引入交换关联泛函,将单电子近似引入的误差归到其中,这使得计算真实体系的电子结构成为可能.但是DFT 方法的主要计算开销在Kohn-Sham 方程的自洽场迭代求解上,其计算复杂度与自由度的立方成正比,即O(N3),这就使得模拟规模受限于算力,往往只能在周期性边界条件下计算包含数十、数百原子的小系统.

效率方面,经典分子动力学(classic molecular dynamics,CMD)通过先验知识构造离子间等效相互势能面的解析表达式,跳过了电子结构的求解,使得计算复杂度降低到了线性标度.结合高性能并行计算,数千万乃至上亿原子的大规模模拟成为可能,我们可以从原子尺度出发,去研究动理论甚至是流体力学尺度的科学问题.但是,由于势能面的准确性受限于模型复杂度,使得CMD 在描述真实体系的物态物性上精度有限,往往无法处理宽热力学区间、多相共存、多元素等复杂体系的物质结构及动力学.

近年来,机器学习技术的发展为建模高维、非线性的势能面提供了有效工具.2007 年,Behler 和Parrinello[10]首先将人工神经网络用于势能面建模问题,提出了单原子能量分解、原子环境描述符以满足势能面建模对广延性、对称性的要求.后续陆续涌现出了高斯近似势(Gaussian approximation potential,GAP)[11]、谱近邻分析势(spectral neighbor analysis potential,SNAP)[12]、深度势能(deep potential,DP)[13]等一系列机器学习势能面构造的应用,不同技术方案的差异在局域原子环境的描述、拟合描述子与能量之间的技术上.在Behler-Parrinello 网络中,局域原子环境描述符往往需要针对不同的材料体系,人为地构造和选取相适应的原子对称函数,限制了该方法的适用性和泛化性.DP 方法则提出了端到端的描述子方案[14],构造嵌入网络以描述复杂的局域原子环境,使得神经网络势能面可以有效应用到金属键、共价键、离子键等不同键合的晶体、团簇、无序体系.随着不断发展,基于机器学习的势能面方法通过学习第一性原理计算标签好的数据集(能量、受力、应力张量),逐步成为建模多相、多元素系统在不同温度、密度条件下相互作用的有效方案[15,16].在保持量子精度的同时,神经网络势能面把计算开销从求解单电子薛定谔方程转变为神经网络模型的推断,有效解决了建模领域的效率精度两难问题,使得研究大尺度、长时间的性质及动力学成为可能[17-21].

本文基于深度势能分子动力学(deep potential molecular dynamics,DPMD),探究了常压到核幔边界条件(p=136 GPa,T=4100 K)下液态铁的动力学与输运性质,并与X 射线非弹性散射实验[22-24]进行对比.二者的一致性使得从实验到理论之间建立了有效的桥梁,使得大规模研究极端条件下真实复杂体系的理化性质成为可能.

2 方 法

在分子动力学模拟中,体系中的原子处理为经典的质点,原子间相互作用由势能面V主导,可以写出相应的牛顿运动方程,即

其中M为原子质量,R表示原子的坐标.该原子间相互作用势能面,由核-核之间的库仑相互作用与电子提供的平均场构成.在传统的KS-DFT 框架中,后者需要通过求解Kohn-Sham 方程来获得.在经典分子动力学的框架中,则是构造解析形式的力场模型去建模势能面.相比之下,在DPMD 中,我们采用了神经网络形式的势能面,它具有更强的描述能力与泛化性能.

深度势能为势能面构造提供了端到端的对称性保持方案,通过嵌入网络提取原子局域环境的特征,以构造保持平移、旋转、置换不变性的描述子.进一步通过拟合网络,建立描述子到单原子能量的高维关系,其模型的数学表述为[14]

其中,Dαi表示第i个原子(元素类型为α)的局域环境描述子,Nαi表示通过神经网络来建立原子i局域环境的描述子到单原子能量的高维映射.有了能量与原子坐标的关系,就可以良好定义力矩阵F(R)≡F∈RN×3和 3×3 的维里张量(Virial tensor)Ξ(R)≡Ξ∈R3×3:

在求解多离子的牛顿运动方程中,每步根据原子在其截断半径内所有近邻原子的相对位置,计算受力、维里张量,并更新下一步每个原子的位置、动量.在获得每个原子的位置与动量({ri},{pi})后,就可以获得系统的任意力学量A(r,p) .通过统计理论可以获得系统的关联与热力学信息.

本文的研究对象是液态铁,初始结构为包含54 个原子的无序结构.首先通过VASP 软件包[25]在一定的温度、密度条件下产生了长度为2 ps 的第一性原理分子动力学轨迹,时间步长为1 fs.每50 步选取一个构型,得到初始训练集.如图1 所示,在初始训练集的基础上,结合同步学习方案,同时训练4 个初始随机数种子不同的DP 模型,沿着铁的熔化曲线探索,从p=0 GPa,T=1873 K 到p=136 GPa,T=4250 K,总共采样了1280 个液相构型,图中黑色实线为Morard 等[26]确定的相边界.

图1 液态铁的训练集采样示意图,紫色空心圆圈表示在该(p,T)条件下采样的液态结构,黑色实线为Morard 等[26]确定的相边界Fig.1.Atomic configurations sampling in P-T space,where the black solid line denotes the experimental measured phase boundaries by Morard et al[26].

在KS-DFT 的计算中,采用了Perdew-Burke-Ernzerhof (PBE) 交换关联泛函[27].赝势则采用了投影缀加平面波 (projector augmented-wave,PAW)形式[28],包含了16 个价电子( 3s23p63d64s2).采用了900 eV 的动能截断.Brillouin 区的k点采样网格为 2×2×2 .此外,考虑到铁的价电子在d 轨道是不满壳层占据,因此计算时考虑了自旋极化.只有考虑自旋极化时,液态铁的静态结构才具有相对可靠的描述[29,30].

在势函数训练步骤,采用开源的深度学习势包DeePMD-kit[31,32].训练过程中,嵌入网格(embedding net)包含三层隐藏全连接层,每层的神经元分别为25,50 和100 个.拟合网络(fitting net)也包含三层隐藏层,每层用了240 个神经元节点,同样为全连接网络.将能量、力和维里项的训练权重初始时刻设置为0.02,1000 和0.02,在训练后期逐渐变化为1.0,1.0 和1.0.训练权重随训练步数变化的主要目的是使神经网络在训练前期优先拟合受力,即势能面的梯度,加快训练的收敛.这是因为梯度内禀包含了能量的相对值以及维里的信息.在训练后期再逐步增加对能量、维里项的权重,使得三者具有相同的拟合权重,保证神经网络对于能量、压强的拟合精度.深度势的截断半径取6.0 Å,训练步数设置为40 万步.最终得到的液态铁多体相互作用势在数学形式上是一个权值冻结的深度神经网络,嵌入网络、拟合网络中各神经元的权重、偏置,在训练中通过反向传播算法得到充分优化.在分子动力学模拟中,神经网络只保留前向传播过程,以进行能量、受力、维里的推断.如图2 所示,以测试集的均方根偏差(root mean squared error,RMSE)来评估模型的准确性,单原子的能量预测偏差为 7.36 meV/atom ,受力偏差为 0.36 eV/Å,压强偏差为 0.41 GPa,结果表现出良好的一致性.

图2 液态铁的能量、受力与压强在训练集上的预测偏差 (a) σE=7.36 meV/atom ;(b) σf=0.36 eV/Å ;(c)σp=0.41 GPaFig.2.DP-predicted energy per atom,force,and pressures versus the true KS-DFT calculations in the testing dataset: (a)σE=7.36 meV/atom ;(b) σf=0.36 eV/Å ;(c) σp=0.41 GPa .

3 结果

3.1 静态结构

液体作为动态无序系统,其原子在中短程的分布特征通常通过径向分布函数g(r),或者是静态结构 因子(static structure factor,SSF)S(q) 来描述.二者互为傅里叶变换,等价描述了粒子数密度的空间关联信息.通常,由于后者可以通过X 射线衍射(X-ray diffraction,XRD)测量直接获得.近年来,Kuwayama 等[23]测量了从常压到核幔边界条件下液态铁的静态结构因子,为验证理论计算的准确性提供了基准.

图3 所示为不同热力学条件下液态铁的静态结构因子,表现出典型的无序系统特征,失去了长程有序性,并具有明显的壳层结构.可以看到,在2.6—2.8 Å-1的波数范围内,SSF展示出一个对称的主峰,其峰值位置随着压强的增加而增加,描述了第一近邻的分布强度.结合径向分布函数可以发现,液态铁作为典型的液体金属系统,配位数在12.6 左右,是典型的密排结构.而非对称的次近邻峰则说明,液态铁具有非常明确的局域二十面体短程结构.与Kuwayama 等[23]的漫散射XRD 测量信号对比,DPMD 的结果表现出良好的一致性,说明PBE 交换关联泛函与自旋极化的考虑准确描述了液态铁的短程结构,而神经网络模型准确还原了第一性原理计算的精度.作为对比,选取Sun等[33]针对地球内核热力学条件优化的嵌入原子势(embedded-atom-method,EAM),该EAM 势在核幔边界条件及高温高压区间(p≥53 GPa,T≥3300 K)同样表现出较好的精度,但是随着压强过渡到常压,对于液态铁的结构描述偏差变大.考虑到经验势在计算效率和规模上的优势,我们可以认为,经验势在针对某些特定条件下的物理问题依然有优势.但是由于经验势的模型复杂度有限,通常无法兼顾从常态到高温高压状态的统一描述.

图3 液态铁的静态结构因子 (a) p=0 GPa,T=1873 K;(b) p=21 GPa,T=2600 K;(c) p=40 GPa,T=3000 K;(d) p=53 GPa,T=3300 K;(e) p=74 GPa,T=3600 K;(f) p=106 GPa,T=4250 K;彩色点为DPMD 计算结果,灰色点为Inui 等[24](常压)和Kuwayama 等[23] (高温高压)的实验测量结果,灰色虚线为基于EAM 势的CMD 计算结果Fig.3.Static structure factor of liquid iron: (a) p=0 GPa,T=1873 K;(b) p=21 GPa,T=2600 K;(c) p=40 GPa,T=3000 K;(d) p=53 GPa,T=3300 K;(e) p=74 GPa,T=3600 K;(f) p=106 GPa,T=4250 K.Colored circles indicate the results from DPMD simulation,the gray square denotes the experimental measurements by Inui et al.[24] and Kuwayaka et al.[23],the gray dashed line denotes the results from CMD simulation with EAM potential.

3.2 动态结构

动态无序系统的动力学性质,包含了其输运性质与集体激发模式.离子-离子动态结构因子(ionion dynamic structure factor,DSF),通过密度涨落的时空关联,描述了一个系统内在不同空间尺度上的弛豫过程与传播过程,是研究离子集体动力学的关键物理量.对于均匀的流体系统,通过对密度-密度时空关联函数做两次傅里叶变换,从时域t到频域ω,从实空间r到倒空间k,可以写出DSF 的表达式为

其中,ρ(k,t) 表示对粒子数倒空间密度的含时演化.DSF 可以通过非弹性X 射线散射(inelastic Xray scattering,IXS)或者是非弹性中子散射(inelastic neutron scattering,INS)实验测量得到,因此它同样为标定方法进而提供了有效基准.对此,我们实现了包含32000 个铁原子的大尺度模拟,与常压条件(p=0 GPa,T=1873 K)下Hosokawa等[22]通过IXS 实验测量的液态铁数据进行对比.

如图4 所示,可以看到DSF 由零频的弹性峰与有限频率的非弹性峰构成,后者描述了入射X 射线光子或者中子与离子集体激发模式发生能量交换的非弹性散射过程.随着q值从小到大,描述离子集体模式的非弹性峰频率发生蓝移.从中可以提取出相应的色散关系ω=ω(q),其斜率为绝热声速.在常压下,DPMD 通过提取纵向集体模式的色散关系得到声速为 3891 m/s,与实验测量的值( (3820±150)m/s)非常接近.DPMD 计算的DSF与IXS 实验信号整体符合也非常好,尤其是小波矢下的弹性峰与非弹性峰的高度、宽度.当波矢足够小的情况下,倒空间无穷小对应实空间无穷大,此时的DSF 过渡到流体力学区间,描述了宏观过程的弛豫与传播过程.在流体力学极限下,DSF 表现出明显的Lorentz 线形,其高度、宽度具有明确的物理意义,可以从中提取出热扩散系数、声衰减系数和比热容比等信息[17].在3.1 节中,EAM 势在极端条件下对液态铁的静态结构同样具有良好的描述,因此我们进一步对比了DP 势与EAM 势在核幔边界条件下的动态结构.有趣的一点是,在小波矢情况下(q=0.38Å-1),EAM势准确描述了液态铁在长波下的集体行为.但随着波矢增加,在动理论区间,可以明显看到EAM 势相比DP 势所描述的液态铁具有相对较高的集体模式能量,相应的纵波声速分别为 7934 和 7600 m/s .这说明EAM 势虽然很好地描述了体系的静态结构,但是在离子动力学上依然有差异.

图4 (a)常压和(b)核幔边界条件下液态铁的动态结构因子,彩色实线为DPMD 计算结果,黑色圆圈为Hosokawa 等[22]的非弹性X 射线散射实验测量结果,黑色虚线为基于Sun 等[33]开发的EAM 势计算的结果Fig.4.Dynamic structure factor of liquid iron under (a) ambient pressure condition and (b) core-mantle boundary condition.Colored lines indicate the results from DPMD simulation,the black circles denote the experimental measurements by Hosokawa et al.[22],the black dashed lines denote the CMD simulation with EAM potential developed by Sun et al[33].

3.3 输运性质

由于DPMD 计算的效率,可以进一步研究对模拟体系大小非常敏感的输运性质.自扩散系数描述了动态无序系统中原子的迁移能力,它可以通过两种方法获得,一个是通过均方根位移(meansquare displacement,MSD):

其中,Ri(t) 表示原子i在t时刻的位移.另一种方式是通过对速度自关联函数(velocity autocorrelation function,VACF)做积分:

其中,v表示原子的速度.本文采用这两种方法,计算了从108 原子体系到16384 原子不同体系大小下的扩散系数.其中时间步长为1 fs,总时间长度为10 ps.每个模拟独立运行了50 次,每次体系有着独立的初始速度分布,目的是做足够多的采样以获得充分系综平均的结果.

如图5 所示,扩散系数在常压和极端条件下都表现出明显的尺寸效应.注意到,基于MSD 和VACF 两种不同的方法在小体系下表现出较大的差异,但是随着体系的增加而逐渐收敛到相同结果.在108 原子体系时,基于VACF 的计算结果,与之前González 等[7]同样采用100 原子体系结合VACF 方法的计算结果非常接近.这说明,DPMD准确还原了AIMD 的精度,同时DPMD 的计算效率使得可以扩大计算规模,消除尺寸效应,获得更为合理的结果.

图5 液态铁的扩散系数 (a) p=0 GPa,T=1873 K;(b) p=96 GPa,T=3800 K;彩色点为DPMD 计算结果,实线为线性关系的拟合结果,黑色三角为González 等[7]的AIMD 计算结果Fig.5.Self-dissufion coefficient of liquid iron: (a) p=0 GPa,T=1873 K;(b) p=96 GPa,T=3800 K.Colored circles denote the results from DPMD simulation,colored solid lines denote the fitting curve for removal of size effect,and the balck triangle denotes the previous AIMD calcualtion by González et al[7].

由于扩散系数表现出强烈的尺度依赖关系,这里引入线性的尺度修正公式Dcorr=Dcalc+kBTξ/(6πηL)来获得无穷大体系下的扩散系数数值.其中,Dcorr为修正值,即无穷大体系的扩散系数,Dcalc为有限尺寸L下计算得到的值,η为剪切黏度,kB为玻尔兹曼常数,T为离子温度,ξ为未知系数,用于拟合.经过拟合,得到在常态(p=0 GPa,T=1873 K)条件下,液态铁的扩散系数为 0.52 Å2/ps,在接近核幔边界条件下,液态铁的扩散系数为 0.48 Å2/ps .该结果指出,DPMD 在研究输运性质时,相比AIMD,表现出计算规模上的优越性.

4 结论

本文构造了液态铁势能面的神经网络模型,热力学条件从常压到核幔边界.在此基础上,结合有限的实验测量数据,验证了DP 模型对液态铁的静态结构、动态结构、热力学性质及输运性质的描述能力.

经过数值实验发现,DP 可以准确描述极端条件下液态铁静态结构,使其与漫散射XRD信号一致.在常压下,液态铁的集体动力学对相互作用描述要求精度更高,要求准确描述密度涨落的时空关联,而DP 模型计算的动态结构因子与IXS 实验基本符合,并给出准确的绝热声速.采用针对核幔条件下液态铁优化的EAM 势,同样计算了常压到核幔边界液态铁的静态结构与动态结构.结果发现,即使是对于液态铁这种结构沿着熔化曲线并不敏感的液态金属体系,传统的EAM也难以兼顾常压到高温高压的统一描述.

进一步讨论了扩散系数的尺寸效应,DPMD使得我们可以把计算尺度从传统的数百原子体系扩展到数万原子体系.结果发现扩散系数强烈依赖于模拟尺寸,并与模拟尺寸的倒数表现出线性依赖关系.这意味着传统AIMD 的计算能力还远远没有达到收敛的条件,也意味着在输运性质的研究上,深度学习与分子模拟的结合会发挥出更大的优势.