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

第 50 卷 第 3 期2023 年北京化工大学学报(自然科学版)Journal of Beijing University of Chemical Technology(Natural Science)Vol.50,No.32023引用格式:陈新,陈卫,张现仁,等.熔融氯化钠晶体介电常数的模拟计算方法J.北京化工大学学报(自然科学版),2023,50(3):59-68.CHEN Xin,CHEN Wei,ZHANG XianRen,et al.A simulation method for calculating the dielectric constant of molten so鄄dium chloride crystalsJ.Journal of Beijing University of Chemical Technology(Natural Science),2023,50(3):59-68.熔融氯化钠晶体介电常数的模拟计算方法陈摇 新1,2,3摇 陈摇 卫1,4*摇 张现仁2,3摇 任摇 瑛1,4摇 李粮生5(1.中国科学院过程工程研究所 多相复杂系统国家重点实验室,北京摇 100190;2.北京化工大学 有机无机复合材料国家重点实验室,北京摇 100029;3.北京化工大学 北京市能源环境催化重点实验室,北京摇 100029;4.中国科学院大学,北京摇 100049;5.电磁散射重点实验室,北京摇 100854)摘摇 要:由于分子模拟过程中存在周期性边界的影响,离子偶极矩定义模糊,导致传统的偶极矩涨落方法在计算熔融盐晶体或含有自由电荷体系的介电常数时存在一定阻碍。根据经典介电理论,通过电流关联函数计算熔融氯化钠(NaCl)晶体体系的介电常数,研究了不同温度下氯化钠的频率相关介电谱,并分析了晶体相变对介电特性的影响。结果表明:当 NaCl 为晶体状态时,介电谱显示出明显的共振吸收峰;在 NaCl 熔化以后,介电谱的共振吸收峰的带宽变宽,波形发生了明显变化。分别利用电流法和外加电场法计算了不同温度下 NaCl 的静态介电常数,结果显示:两种方法的计算结果的一致性良好,从而验证了本文方法的准确性;NaCl 晶体在发生固液相变时,静态介电常数会在相变点发生跳变。关键词:电流关联函数;介电常数;熔融氯化钠晶体;相变中图分类号:TQ015郾 9;TP391郾 9;O6鄄39摇 摇 DOI:10.13543/j.bhxbzr.2023.03.007收稿日期:2023-01-05基金项目:多相复杂系统国家重点实验室开放研究基金(MPCS鄄2021鄄D鄄10);国 家 自 然 科 学 基 金(21973097/21978007);国 家 自 然 科 学 基 金 创 新 研 究 群 体 项 目(21821005)第一作者:男,1998 年生,硕士生*通信联系人E鄄mail:chenwei 引摇 言随着微波技术的不断发展,电介质材料的介电常数作为描述材料电磁特性的重要参数,反映了材料对外界电磁场的响应,目前已广泛应用于飞机隐身材料1-2、雷达探测3-4以及 5G 通讯5-6等重要领域。介电弛豫光谱法7-9由于具有非侵入性、检测迅速等优点,常用于测定电介质材料的介电常数,但该方法往往受限于实际应用中的复杂环境以及缺乏较为精确的参数模型。近年来,随着计算机技术的不断发展,分子动力学模拟方法为学者在原子或分子尺度上理解材料的介电特性提供了重要帮助。除了可以有效预测实际工作情况下电介质的介电特性,目前分子模拟方法已经将频率拓展到太赫兹(THz)光谱。Kirkwood10提出了介电常数的经典统计理论,这是利用计算机模拟介电常数的基石;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)。摇 摇 着(棕)=1+4仔3VkBTF-d渍(t)dt(1)式中:V 为模拟体系的体积,kB为玻尔兹曼常数,T为模拟体系温度。摇 摇F-d渍(t)dt=乙+肄0e-J2仔棕t-d渍(t)dtdt=渍(0)-J2仔棕乙+肄0e-J2仔棕t渍(t)dt,J2=-1(2)式中,渍(t)是总偶极矩 Mtot的自相关函数。摇 摇 渍(t)=骉Mtot(0)Mtot(t)骍(3)考虑一个由 N 个带电粒子组成、体积为 V、温度为 T、粒子位移矢量为 ri以及所带电荷量为 Qi(i=1,2,N)的体系,根据 Schr觟der 等26的研究,可以将体系的总偶极矩 Mtot划分为平动偶极矩 MJ和质心修正偶极矩 MD,如公式(4)所示。摇 摇 Mtot(t)=MD+MJ(4)考虑到模拟系统中存在 Nmol个分子,分子 i 包含 ni个电荷,每个自由离子也被视为 ni=1 的分子,因此移Nmoli=1ni=N。Qij和 rij分别表示分子 i 上的第 j 个电荷的电荷量和位置向量,则总偶极矩可以表示为摇 摇 Mtot=移Ni=1Qirij=移Nmoli=1移nij=1Qijrij=移Nmoli=1移nij=1Qij(rij-ricm)+Qijricm=移Nmoli=1滋i+移Nmoli=1Qiricm(5)公式(5)中,等号右边第一项表示质心矫正的分子偶极矩 滋i,可以表示为摇 摇 滋i=移nij=1Qij(rij-ricm)(6)式中,ricm表示分子 i 的质心位置。故考虑分子上每个电荷的坐标都表示为相对于质心的相对位置。由于分子偶极矩 MD不受限于周期性边界条件,因此可以对其进行有效定义和计算。对于带电荷的自由离子,由于离子位置与质心位置重合,即 rij=ricm,因此对分子偶极矩不提供贡献。公式(5)中,等号右边第二项是平动偶极矩MJ,摇 摇 MJ=移Nmoli=1Qiricm(7)式中,Qi=移nij=1Qij,是分子 i 的净电荷。由于周期性边界条件的存在,体系中所带电荷离子穿过周期性边界条件导致偶极矩值发生非物理抖动20。因此本文利用电流法20来计算离子对介电常数的贡献,由于电流 J 的定义是粒子速度 vicm和电荷量 Qi的乘积(式(8),摇 摇 J=移Nmoli=1Qivicm=移Nmoli=1Qidricmdt=dMJdt(8)将公式(4)和(8)带入公式(1)的傅里叶-拉普拉斯06北京化工大学学报(自然科学版)摇 摇 摇 摇 摇 摇摇摇 摇 摇 摇 摇 摇 摇 摇 2023 年变换中,可以得到频率相关介电常数的表达式(式(9)。摇 摇 着(棕)=1+4仔3VkBT掖M2D业+J棕 F 掖MD(0)MD(t)业-掖MD(0)J(t)业+J棕4仔3VkBTJ掖J(0)J(t)业+J棕 F 掖J(0)MD(t)业(9)由于公式(9)考虑了体系中分子偶极贡献与离子电流贡献,故可广泛作为体系的介电常数解析表达式。对于本文考虑的熔融 NaCl 晶体体系,随着温度升高,NaCl 晶体结构被破坏,体系中含有自由移动的离子电流,离子 i 的坐标位置 rij等于其质心坐标位置 ricm,那么体系中分子偶极矩 滋i=0,即体系中不存在分子偶极矩贡献,故仅考虑体系中离子电流贡献,其频率相关介电常数的解析表达式为着(棕)=1+J棕4仔3VkBTJ掖J(0)J(t)业(10)当频率小于 103Hz 时,体系的介电常数称为静态介电常数,利用电流法计算静态介电常数 着static,其表达式为着static=1-4仔3VkBT乙+肄0t掖J(0)J(t)业dt(11)1郾 2摇 分子模拟方法利用分子建模脚本 Moltemplate27搭建了包含有 1 000 个原子的 NaCl 晶体的初始构型,立方体模拟盒在 3 个方向上的大小均为 2郾 8 nm。对模拟盒的 3 个方向都应用周期性边界条件,选用经典原子间相互作用势函数 Born-Mayer-Huggins(BMH)势28-30描述熔融 NaCl 晶体离子之间的相互作用势。摇 摇 uij=Aije-rij籽ij-Cijr6ij-Dijr8ij+qiqj4仔着0rij(12)式中: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 进行一次体系平均,从而计算静态介电常数、频率相关介电谱以及进行后处理。使用可视化软件 VMD32对轨迹文件进行观察。本研究的模拟温度区间为 300 1 500 K,调节不同温度大小是为了揭示 NaCl 晶体的相变演化过程,研究相变过程中晶体材料介电特性的变化。NaCl 晶体的初始构型由 VMD 绘制32,如图 1 所示。绿色代表氯离子,黄色代表钠离子。图 1摇 NaCl 晶体的初始结构模型Fig.1摇 Initial structure model of NaCl crystals摇2摇 结果与讨论2郾 1摇 NaCl