熔融氯化钠晶体介电常数的模拟计算方法

陈 新 陈 卫 张现仁 任 瑛 李粮生

(1.中国科学院过程工程研究所 多相复杂系统国家重点实验室,北京 100190;2.北京化工大学 有机无机复合材料国家重点实验室,北京 100029;3.北京化工大学 北京市能源环境催化重点实验室,北京 100029;4.中国科学院大学,北京 100049;5.电磁散射重点实验室,北京 100854)

引 言

随着微波技术的不断发展,电介质材料的介电常数作为描述材料电磁特性的重要参数,反映了材料对外界电磁场的响应,目前已广泛应用于飞机隐身材料[1-2]、雷达探测[3-4]以及5G 通讯[5-6]等重要领域。介电弛豫光谱法[7-9]由于具有非侵入性、检测迅速等优点,常用于测定电介质材料的介电常数,但该方法往往受限于实际应用中的复杂环境以及缺乏较为精确的参数模型。

近年来,随着计算机技术的不断发展,分子动力学模拟方法为学者在原子或分子尺度上理解材料的介电特性提供了重要帮助。除了可以有效预测实际工作情况下电介质的介电特性,目前分子模拟方法已经将频率拓展到太赫兹(THz)光谱。Kirkwood[10]提出了介电常数的经典统计理论,这是利用计算机模拟介电常数的基石;Neumann 等[11-13]基于线性响应理论,利用分子动力学模拟将晶体的偶极矩涨落与周期性边界条件下的频率相关介电常数联系起来,得到频率相关介电常数的计算表达式;Caillol等[14-16]以离子-溶剂混合物体系为研究对象,根据电流和微观极化的相关函数,计算与频率相关的介电常数和电导率,得到相对复杂的复介电常数表达式,该表达式包含偶极子-偶极子、电流-偶极子、电流-电流相关函数的贡献;Thomas 等[17]提出了利用从头计算分子动力学的方法计算二丁醇的磁偶极矩;Chen 等[18]利用机器学习结合分子动力学的方法提高了介电常数的计算精度;Yang 等[19]研究了高分子溶液中浓度变化导致的聚合物结构转变对体系介电常数的影响。目前,对电介质材料介电常数的研究主要集中于固体材料和聚合物体系方面,但是针对体系中存在含净电荷的离子、分子或熔融盐晶体的计算机模拟介电常数的相关研究却很少。这是由于传统的偶极矩涨落方法无法适用于体系中存在含自由电荷的粒子的情况,其原因是:体系内存在导电介质,无法维持静电场平衡;周期性边界对离子偶极矩计算会产生数值跳跃,这是因为偶极矩定义为粒子所带电荷与位移矢量的乘积,当体系中施加周期性边界条件时,对于类似水的这种极性分子,其具有完整的分子信息,周期性边界条件并不会影响分子的偶极矩,然而对于含有自由电荷的体系,自由电荷可跨越周期性边界条件,导致体系的偶极矩值不能收敛[20]。针对上述问题,Sega 等[21-22]提出了利用电导率计算NaCl 溶液介电频谱的方法;Schröder 等[23-25]致力于计算离子液体介电常数的研究,建立了计算离子液体介电常数和电导率的理论基础,研究分析了极化的微观机制对离子液体广义介电常数或磁化率的贡献。

本文基于Neumann 介电理论公式[11-13],进一步推导了适用于离子溶液或熔融盐晶体体系的介电常数解析表达式;以熔融氯化钠(NaCl)晶体为研究算例,利用分子动力学方法测定了不同温度下NaCl晶体的介电常数和频率相关的介电谱,研究了NaCl晶体熔化过程中介电特性的变化;最后利用经典外加电场法计算了熔融NaCl 晶体的静态介电常数,从而对本文的方法进行了验证。

1 介电常数计算及分子模拟方法

1.1 介电常数计算

基于线性响应理论,Neumann 等[11-13]将分子动力学模拟得到的晶体偶极矩涨落与周期性边界条件下的频率相关介电常数联系起来,得到与频率ω相关的经典介电常数ε的计算表达式(式(1))。

式中:V为模拟体系的体积,kB为玻尔兹曼常数,T为模拟体系温度。

式中,φ(t)是总偶极矩Mtot的自相关函数。

考虑一个由N个带电粒子组成、体积为V、温度为T、粒子位移矢量为ri以及所带电荷量为Qi(i=1,2,…,N)的体系,根据Schröder 等[26]的研究,可以将体系的总偶极矩Mtot划分为平动偶极矩MJ和质心修正偶极矩MD,如公式(4)所示。

考虑到模拟系统中存在Nmol个分子,分子i包含ni个电荷,每个自由离子也被视为ni=1 的分子,因此分别表示分子i上的第j个电荷的电荷量和位置向量,则总偶极矩可以表示为

公式(5)中,等号右边第一项表示质心矫正的分子偶极矩μi,可以表示为

公式(5)中,等号右边第二项是平动偶极矩MJ,

将公式(4)和(8)带入公式(1)的傅里叶-拉普拉斯变换中,可以得到频率相关介电常数的表达式(式(9))。

由于公式(9)考虑了体系中分子偶极贡献与离子电流贡献,故可广泛作为体系的介电常数解析表达式。对于本文考虑的熔融NaCl 晶体体系,随着温度升高,NaCl 晶体结构被破坏,体系中含有自由移动的离子电流,离子i的坐标位置等于其质心坐标位置,那么体系中分子偶极矩μi=0,即体系中不存在分子偶极矩贡献,故仅考虑体系中离子电流贡献,其频率相关介电常数的解析表达式为

当频率小于103Hz 时,体系的介电常数称为静态介电常数,利用电流法计算静态介电常数εstatic,其表达式为

1.2 分子模拟方法

利用分子建模脚本Moltemplate[27]搭建了包含有1 000 个原子的NaCl 晶体的初始构型,立方体模拟盒在3 个方向上的大小均为2.8 nm。对模拟盒的3 个方向都应用周期性边界条件,选用经典原子间相互作用势函数Born-Mayer -Huggins(BMH)势[28-30]描述熔融NaCl 晶体离子之间的相互作用势。

2月8日，徐云天让弟弟给父亲打电话，谎称得了重病。徐河来了，见小儿子根本没病，他转身就要走。徐云天鼓足勇气劝道：“爸，弟弟今年就要中考了，就算你跟我妈的缘份到头了，也要替弟弟着想啊！”

式中:uij为NaCl 晶体离子之间的相互作用势,Aij为BMH 势排斥参数,Cij和Dij为范德华参数,rij为粒子i和j之间的距离,ρij为离子对ij的长度参数,qi和qj分别为i和j的电荷量,ε0为真空介电常数。

选用经典分子模拟软件Lammps 进行分子动力学模拟[31],选取随机数种子为825 577,使体系中粒子的速度在给定的温度下满足玻尔兹曼分布;之后运行100 ps 的能量最小化过程,防止不合理的初始结构导致模拟过程崩溃。采用NPT 系综,分别使用Nosé-Hoover 热浴法和Berendsen 耦合方法使体系维持给定的温度和1 个标准大气压。使用Ewlad 求和法计算长程库仑相互作用,截断半径为12 Å,采用Velocity-Verlet 算法求解运动方程的积分,时间步长设置为1 fs;每次进行MD 模拟过程时首先运行2 ns 以达到平衡状态,之后进行1 ns 的数据采样过程,每隔0.005 ps 进行一次体系平均,从而计算静态介电常数、频率相关介电谱以及进行后处理。使用可视化软件VMD[32]对轨迹文件进行观察。本研究的模拟温度区间为300 ～1 500 K,调节不同温度大小是为了揭示NaCl 晶体的相变演化过程,研究相变过程中晶体材料介电特性的变化。NaCl 晶体的初始构型由VMD 绘制[32],如图1 所示。

图1 NaCl 晶体的初始结构模型Fig.1 Initial structure model of NaCl crystals

2 结果与讨论

2.1 NaCl 晶体相变演化

在1 个标准大气压下,分别在300、500、700、900、1 100、1 180、1 200、1 220、1 240、1 300、1 500 K的温度下计算体系的体积V,结果如图2 所示。可以看出,体系体积随着温度的升高而增加。当温度升高至1 240K 时,体积开始快速膨胀,这表明固态NaCl 开始熔化为液态,这与实验测得的NaCl 的熔点(1 074 K)[33]基本吻合。当温度达到1 300 K 时,体系已经完全成为液态。

图2 NaCl 的体积-温度相图Fig.2 Volume-temperature phase diagram of NaCl

此外,径向分布函数(radial distribution function,RDF)[34]也可以作为分析体系结构变化的重要参数用来研究物质的有序性,其定义为一个系统在距离一个标记粒子的r处找到其他粒子的概率。图3 为不同温度下体系中Na-Cl 的RDF 图,其中g(r)Na-Cl表示钠离子周围距离半径r到r+dr处发现氯离子的概率。当温度由1 180 K 升高至1 240 K时,RDF 在不同位置呈现出4 个离散峰,并且随着温度升高,g(r)Na-Cl几乎没有增加,最大峰值大约在3.4 左右,表示此时熔融态NaCl 晶体中仍有较强的约束作用,结构只存在轻微扰动。随着温度升高到1 300 K 时,原来有序的熔融NaCl 晶体结构被破坏,相互作用力减弱,分子结构变得杂乱无序,RDF 峰值发生了显著变化,并且仅保留两个较为明显的离散峰,表明体系已经转变为液态。

图3 不同温度下Na-Cl 的径向分布函数Fig.3 Radial distribution function of Na-Cl at different temperatures

2.2 介电频谱

根据上文推导的介电常数的计算方法,为了获得不同温度下熔融NaCl 晶体体系的频率相关性,首先需要获得体系中离子电流-电流随时间的变化规律。由于时间相关函数本身就是信号与自身之间相关性程度的度量,因此深入分析相关函数的行为有助于理解频率相关介电谱的特性[21]。图4 给出了1 180 ～1 300 K 的温度下熔融NaCl 晶体的电流相关函数。可以看出,NaCl 电流自相关函数在约10-1ps处呈现极大的振荡行为,随后电流自相关函数逐渐衰减为零,表明该模拟体系在给定的时间尺度内达到收敛。值得注意的是,当温度升高至1 300 K 时,不仅电流相关函数迅速降低,而且峰值出现的频率位置也发生了显著变化,这可能表明随着体系相态结构的改变,共振频率发生了变化。由于时间相关函数十分依赖于较长的模拟时间尺度,并且在分子模拟过程中,由于原子轨迹上的每个点只提供一个值,难以对每个粒子残留物进行平均,这导致即使在长时间模拟下,无穷远处的自相关函数仍然在零附近不停振荡。此外,使用傅里叶-拉普拉斯变换的方法计算解析介电频谱时,离散傅里叶变换算法将从分子动力学模拟中获得的有限数据点进行周期性扩展,这并不是公式中简单的连续函数积分。为了减少数学算法差异,往往需要更多的数据点,这增加了模拟的负担,而离散傅里叶变换算法中使用的矩形积分方法将会引入显著的误差。因此,本文利用非线性最小二乘法将模拟数据拟合成特定的连续函数[35],并直接解析计算这些函数的傅里叶-拉普拉斯变换,这大大平滑了介电谱中的统计粗糙度。

图4 不同温度下NaCl 的电流相关函数Fig.4 Current correlation function of NaCl at different temperatures

由于电流自相关函数的行为类似于阻尼振荡,通常在几皮秒内迅速振荡,随后逐渐衰减到零,因此可以使用多指数余弦相关函数来拟合相应的模拟结果。

表1 电流相关函数的拟合参数Table 1 Fitting parameters of the current correlation function

图5 不同温度下NaCl 电流相关函数的拟合结果Fig.5 Fitting results of the current correlation function of NaCl at different temperatures

将电流拟合函数代入公式(10),利用Python脚本可以得到熔融NaCl 晶体的频率相关介电谱。图6 为对数坐标系中熔融NaCl 晶体在不同温度下的频率相关介电谱。在1 180 ～1 240 K 的温度范围内,随着频率的增加,大约在3 THz 处外场与NaCl 晶体的偶极子振荡频率发生共振,介电常数的实部εreal开始加速增大,达到峰值后快速减小,演化为接近于无外场作用下的极化状态。同样地,在达到共振频率时,虚部介电损耗εimag从0 逐渐增大到峰值,随后又减小到0。当温度进一步升高到1 300 K 时,NaCl 开始熔化发生相变,与NaCl为晶体状态时相比,熔融NaCl 晶体的介电谱发生了显著变化。在NaCl 熔化以后,由于体系中存在不同组分,各组分的介电响应相互重叠,介电谱的共振吸收峰的带宽变宽,波形发生了明显变化,导致体系的介电频谱十分复杂。可以发现,此时介电损耗峰值降低到约1.5,而实部处的共振峰值几乎消失。这是因为NaCl 晶体结构随体系温度改变而发生了变化,体系中离子组分随频率增加发生离子极化,分子的固有偶极矩的转向极化不再响应。

图6 不同温度下NaCl 的频率相关介电谱Fig.6 Frequency dependent dielectric spectra of NaCl at different temperatures

2.3 外加电场分子动力学模拟

为了进一步验证上述推导的介电常数计算的准确性,本文使用外加电场法计算静态介电常数。所谓静态介电常数是指当频率很低或接近零时,介电常数是与电场频率无关的量。利用外加电场法计算体系的静态介电常数εstatic,其表达式为

式中:E为外电场强度,ME为施加电场后介质的偶极矩,V为模拟体系的体积。因此,可以利用分子动力学模拟的方法在平衡态的NaCl 分子构型的Z方向上施加静电场EZ,电场强度以0.01 V/Å 为间隔,从0.01 V/Å 逐渐增大至0.1 V/Å。然而当T=1 100 K 时,如果继续向体系中施加0.1 V/Å 的电场强度,此时NaCl 晶体的偶极矩大小受电场强度的影响会发生剧烈波动,这会无法保证电场强度与诱导偶极矩之间高度一致的线性关系,从而影响静态介电常数的估算。因此本文选择从0.001 V/Å 逐渐增大至0.01 V/Å。根据电场诱导偶极矩和电场强度之间的线性关系,利用最小二乘法拟合可得到不同温度下NaCl 晶体的静态介电常数。

图7 为诱导偶极矩和电场强度之间的线性拟合结果。可以发现,随着体系施加的电场强度的增加,诱导偶极矩增大。在300 ～1 100 K 的温度范围内,电场强度和诱发偶极矩之间保持着良好的线性关系。随着温度进一步升高至1 300 K,NaCl 晶体开始熔化,采样点的波动明显增大,这可能与NaCl 晶体结构被破坏有关。

图7 不同温度下电场诱导偶极矩MZ和电场强度EZ之间的线性关系Fig.7 Linear relationship between the electric field induced dipole moment and the electric field intensity at different temperatures

分别采用电流法和外加电场法计算NaCl 的静态介电常数与温度的相关性,结果如图8 所示。可以看出:当NaCl 为固态晶体时,静态介电常数随着温度的升高而逐渐增加;当温度升高到固液临界点(T=1 100 K)时,静态介电常数达到最大;随着温度进一步升高,NaCl 发生固-液相变过程,体系的静态介电常数迅速下降。电流法和外加电场法的计算结果有着良好的一致性,电流法能够有效反映NaCl 体系相变过程中静态介电常数的变化规律,广泛适用于熔融盐晶体以及溶液中含有自由电荷的体系。根据文献[36]的实验测定结果,室温下NaCl 晶体的静态介电常数在5.45 ～5.9 之间。本文用电流法计算得到的静态介电常数为4.5(T=300 K),小于实验观测值,这可能与描述NaCl 晶体的力场参数有关。在今后的研究中,本文将采用机器学习结合量子力学第一性原理的计算方式来提高力场的精度。此外,由公式(11)可知,对于电流法,如果要准确计算静态介电常数,需要从0 积分到无穷大,而实际上本文首先用最小二乘法将电流关联函数拟合成一个解析函数,然后对解析函数进行积分,这样会引入误差,导致计算结果偏低。

图8 采用电流法和外加电场法计算得到的NaCl静态介电常数与温度的相关性Fig.8 Correlation between the static dielectric constant of NaCl and the temperature calculated by the current method and the external electric field method

3 结论

(1)基于经典介电理论公式,进一步推导了适用于熔融盐晶体以及溶液中含有自由电荷的体系的介电常数表达式,该公式可以用于计算熔融NaCl 晶体的介电性质。

(2)利用电流拟合相关函数分析了不同温度下NaCl 的频率相关介电谱,结果表明:随着温度升高(300 ～1 100 K),在共振吸收峰处,介电常数的实部先增大,达到峰值后快速减小;当温度升高至1 300 K 时,NaCl 晶体发生相变,介电损耗显著下降,在共振频率处介电常数的实部峰值几乎消失。

(3)分别利用电流法和外加电场法计算了不同温度下NaCl 的静态介电常数,结果显示:随着体系温度升高(300 ～1 100 K),NaCl 晶体的介电常数呈上升趋势;当温度高于相变点时,NaCl 的静态介电常数快速下降。电流法和外加电场法的计算结果有着良好的一致性,表明电流法可以有效描述NaCl 的相变过程。