一种基于力匹配的铝电解熔盐体系力场拟合方法与流程

本发明属于铝电解技术领域，具体涉及一种基于力匹配的铝电解熔盐体系力场拟合方法。

背景技术：

铝电解工业作为国民经济的支柱产业，一直以来受到人们的广泛关注，而铝电解中最为重要的就是铝电解质，作为电解铝的反应物之一，人们对其结构和性质的研究从未停止。但是由于实验条件的限制，熔盐电解质体系相关的性质很难通过实验获取，并且相关的作用机理也无法获悉。近年来，随着计算化学的逐渐发展，采用计算模拟的方法来研究物质的结构和性质受到越来越多人的认可。分子动力学模拟方法作为一种研究熔盐离子结构和输运性质较好的方法已经得到了广泛的应用，分子动力学模拟过程最为重要的就是找到合适的势函数以及势参数。但是适用于铝电解熔盐体系的势函数及相关的势参数匮乏，这成了研究的最大的阻碍。buckingham作为一种二体势已经应用于硅酸盐体系的分子模拟中，考虑到铝电解熔盐体系同样属于高温体系，buckingham势函数同样可以运用于铝电解熔盐体系的结构和性质的研究中。选择了合适的势函数以后，合适的势参数就成为计算能否成功的关键。

目前，急需开发一种简单有效的适合于铝电解熔盐体系分子动力学模拟的势参数构建方法，以适用于不同的铝电解熔盐体系，为铝电解生产提供有效指导。

技术实现要素：

本申请旨在至少解决现有技术中存在的技术问题之一。为此，本发明的目的之一在于提供一种基于力匹配的铝电解熔盐体系力场拟合方法。该方法通过力匹配结合非线性拟合方法得到适用于铝电解熔盐体系分子动力学模拟的势参数，进而通过计算得到铝电解熔盐体系的离子结构和输运性质。

为解决上述技术问题，本申请采用如下技术方案：

一种基于力匹配的铝电解熔盐体系力场拟合方法，包括如下步骤：

步骤1：根据铝电解熔盐体系的密度和熔点建立计算模型，对模型进行第一性原理分子动力学模拟，得到稳定状态下的熔盐中离子的运动轨迹；

步骤2：提取步骤1中稳定状态下每一帧运动轨迹中离子的受力矩阵b和位置矩阵c；

步骤3：设定周期性边界条件，将模拟体系在三维坐标系下进行n*n*n倍的扩大，提取扩大后整个体系中离子的位置坐标矩阵d；

步骤4：根据步骤3中扩大的整个体系中离子坐标计算离子间的距离矩阵e；

步骤5：根据步骤4中离子间的距离矩阵计算任意离子所受的库仑力；

步骤6：对步骤5中的库伦力进行分解得到非库仑力；

步骤7：将buckingham势函数对距离进行求导得到非库伦力的近似函数；

步骤8：结合步骤6和7进行非线性拟合得到适用于铝电解熔盐体系分子动力学模拟的势参数。

步骤9：根据所构建的势参数对熔盐体系进行分子动力学模拟，计算熔盐体系的离子结构和输运性质。

进一步的，步骤6中的库伦力分解过程按如下公式计算：

g＝b-f

其中，g表示体系中任意离子所受的非库伦力矩阵。

进一步的，步骤7中buckingham势函数以及相应的非库伦力的近似函数形式如下：

其中，v(rij)为体系的buckingham势能，aij，βij，cij是buckingham势参数，用来描述两离子间的相互作用。ff(rij)为任意两离子间的非库仑力，将离子间的非库伦力进行叠加得到所有离子所受的非库仑力矩阵记为h。

进一步的，步骤8中的非线性拟合方法根据如下公式计算：

min∑(h(x,y)-g(x,y))²

基于最小二乘原理，在matlab中对函数进行非线性拟合得到适用于熔盐体系分子动力学模拟的buckingham势参数。

进一步的，步骤1中第一性原理分子动力学模拟的体系大小一般设置为100-200个离子，模拟时间一般为10-20ps。稳定状态取最后的10ps的离子运动轨迹。

进一步优选，步骤2中离子的受力矩阵和位置矩阵来自于步骤1中稳定状态下每一帧运动轨迹文件，所有的矩阵均表示向量形式。其中，离子的受力矩阵设为b，位置矩阵设为c。

进一步优选，步骤3中n的值根据模拟盒子的边长来定，一般只考虑中心模拟盒子以内的范德华力和库仑力，n的取值按如下公式计算：

其中，表示向上取整符号，l表示模型盒子的边长。扩展后的离子位置坐标矩阵设为d。

进一步的，步骤4中离子距离r按如下公式计算：

将所有的离子对的距离按顺序存储在一个矩阵中e中。

进一步的，步骤5中的库仑力按如下公式计算：

其中，fc(rij)表示任意两离子间库仑力，k表示库伦常数，qi和qj分别表示离子所带电荷，rij表示离子间距离。将所有离子对的库仑力进行叠加后计算得到每个离子所受的库仑力，结果按顺序记录在矩阵f中。

进一步的，步骤9中离子结构和输运性质的计算过程如下：

(1)根据熔盐体系的密度和熔点建立计算模型；

(2)对计算模型进行分子动力学模拟，收集稳定状态下的粒子轨迹数据；

(3)对得到的粒子轨迹数据进行统计分析，得到熔盐体系的离子结构和输运性质(键长，键角，配位数，粘度和电导率等)。

与现有技术相比，本发明的优点有：

本发明通过力匹配结合非线性拟合方法得到适用于铝电解熔盐体系分子动力学模拟的势参数，除了基本的密度和熔点数据外，大部分数据均来源于第一性原理计算得到，从而无需进行高温熔盐实验，节约了实验成本，避免了实验过程中对设备的损害。

2、本发明构建的势参数可以运用到铝电解熔盐体系的分子动力学模拟研究中，计算熔盐体系的离子结构和输运性质，通过与实际工况相比较，发现离子结构和输运性质与实际试验吻合度高，从而可以为铝电解生产上选择合适的电解质体系提供有益参考。

3、本发明在非线性拟合过程中的数据充足，足以根据第一性原理计算得到的离子所受非库仑力拟合得到相应的势参数，计算过程严谨，避免了数据量过少带来的拟合误差。

附图说明

图1为本发明的流程图；

图2为naf的计算模型图；

图3为基于力匹配的非线性拟合结果图；

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，一种基于力匹配的铝电解熔盐体系力场拟合方法，包括如下步骤：

步骤1：对铝电解熔盐体系进行第一性原理分子动力学模拟，得到稳定状态下的熔盐中离子的运动轨迹；以naf体系为例，naf的密度为1.125g/cm³，熔点为993℃，建立含100个离子的naf计算模型，模型的尺寸为如图2所示。第一性原理分子动力学模拟时间为20ps，取10-20ps的模拟轨迹进行后续的计算。

步骤2：提取步骤1中稳定状态下每一帧运动轨迹中离子的受力矩阵和位置矩阵；离子的受力矩阵和位置矩阵来自于步骤1中稳定状态下的运动轨迹文件，所有的矩阵均表示向量形式。其中，离子的受力矩阵设为b，位置矩阵设为c。

步骤3：设定周期性边界条件，n的计算过程如下：

因此，将模拟体系在三维坐标系下进行2*2*2倍的扩大，提取扩大后整个体系中离子的位置坐标矩阵设为d。

步骤4：根据步骤3中扩大的整个体系中离子坐标计算离子间的距离矩阵，离子间距离r按如下公式计算：

将所有的离子对的距离按顺序存储在一个矩阵中e中。

步骤5：根据步骤4中离子间的距离矩阵计算任意离子所受的库仑力，库仑力按如下公式计算：

其中，fc(rij)表示任意两离子间库仑力，k表示库伦常数，qi和qj分别表示离子所带电荷，rij表示i离子和j离子的距离。将所有离子对的库仑力进行叠加后计算得到每个离子所受的库仑力，结果按顺序记录在矩阵f中。

步骤6：对步骤5中的库伦力进行分解得到非库仑力，库伦力分解过程按如下公式计算：

g＝b-f

其中，g表示体系中任意离子所受的非库伦力矩阵。

步骤7：将buckingham势函数对距离进行求导得到非库伦力的近似函数，buckingham势函数以及相应的非库伦力的近似函数形式如下：

其中，v(rij)为体系的buckingham势能，aij，βij，cij是buckingham势参数，用来描述两离子间的相互作用。ff(rij)为任意两离子间的非库仑力，将离子间的非库伦力进行叠加得到所有离子所受的非库仑力矩阵记为h。

步骤8：结合步骤6和7进行非线性拟合得到适用于铝电解熔盐体系分子动力学模拟的势参数，非线性拟合方法根据如下公式计算：

min∑(h(x,y)-g(x,y))²

其中，x表示矩阵的行，y表示矩阵的列，基于最小二乘原理，在matlab中对函数进行非线性拟合得到适用于熔盐体系分子动力学模拟的buckingham势参数。拟合结果如图3所示，得到相应的势参数如表1所示。

表1naf体系中不同离子对的buckingham势参数

步骤9：根据所构建的势参数对naf熔盐体系进行分子动力学模拟，计算熔盐体系的离子结构和输运性质。模型的构建过程如步骤1所示，所不同的是采用本文所构建的势参数可以计算更长的时间，这意味着统计平均的数据更多，得到的结果也更为可靠。对图2的模型进行分子动力学模拟，总的模拟时间为300ps，取200-300ps时间段的粒子轨迹数据进行后续的离子结构和输运性质计算。这里采用matlab软件对粒子轨迹数据进行统计分析，所得的部分结果如表2所示。通过将表2的数据与试验测量数据进行对比，发现通过本方法获得的模拟数据与实际试验数据吻合度高，从而可以为铝电解生产上选择合适的电解质体系提供有益参考。

表2naf体系的离子结构和输运性质部分计算结果

需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。