CN110070918B - 基于分子间相互作用的粗粒化方法 - Google Patents

Abstract

本发明公开了一种基于分子间相互作用的粗粒化方法，包括以下步骤，对晶体材料进行量子力学(QM)计算、对三斜晶系进行校正、建立伪化学键模型、对Morse函数参数进行拟合、将参数、模型带入MD框架进行模拟、输出模拟结果。利用本发明提供的方法，能够计算材料颗粒之间的相互作用，能够显著的提高计算机模拟的效率，相对的提高计算准确度。

Description

基于分子间相互作用的粗粒化方法

技术领域

本发明涉及计算机模拟领域，尤其涉及计算机模拟中粗粒化模型的搭建和对应力场的计算。

背景技术

计算机模拟经过了几十年的发展，不论在理论方面还是在应用方面都取得了巨大成就。总的来说，计算机模拟的发展趋势是以更高的效率、模拟更大的体系、实现更长的演化时间、取得更精确的模拟结果为目的。为了实现这些目标，必须从计算机技术、模拟算法等方面进行广泛而深入的研究。

在经历了约半个世纪的指数式提高，计算机核心部件CPU在21世纪初期后，出现了停滞现象，失去了过去那种按摩尔定律快速提高的趋势。因此，要实现计算机模拟高效率、大体系、长时间、高精确的目的，一味的堆积硬件资源显然不可取。而改进模拟算法、简化分子模型可以达到任何计算设备的改进均无法实现的作用。

现目前主流的粗粒化方法有两种：

方法一：基于Martini力场的粗粒化方法，就是按照一定的规则方法将一个或几个原子团约化成一个粗粒化的粒子，在计算过程中，只考虑粒子的各种性质，而不考虑原先的原子团。基于Martini力场的粗粒化方法可以降低计算量，但是对于一些大体系、有着空间拓扑结构的物质，它的计算结果并不是很准确。

方法二：耗散粒子动力学(DPD)方法，其单个粒子代表整个分子或包含多个分子，或高分子的一个片段的流体区域，而不是单个原子，并且不考虑原子的行为细节，认为其与过程无关。粒子自身的自由度被整合，粒子间的受力由一对保守力、耗散力与随机力表示。耗散粒子动力学方法在流体模拟中有着较高的计算准确度，但是却不能准确的描述其他体系。

计算机模拟中消耗计算时间最多的部分就是分子间相互作用、静电相互作用、多体相互作用等。

发明内容

本发明的目的是，提供一种在分子动力学(MD)框架下的，基于分子间相互作用的粗粒化方法，利用本发明提供的方法，能够计算材料颗粒之间的相互作用，能够显著的提高计算机模拟的效率，相对的提高计算准确度。

为此本发明采用如下技术方案：基于分子间相互作用的粗粒化方法，包括下列步骤：

步骤1、对晶体材料进行量子力学(QM)计算：根据晶体材料的晶向对晶体材料的单体分子进行建模，并用Gaussian 09进行计算，计算完成后，改变两个单体分子之间的距离，重新计算；经过一系列的计算后，得到距离-能量矩阵。

步骤2、对三斜晶系进行校正：对于三斜晶系材料，对其晶体结构进行正交化处理，得到斜方晶系；正交化处理如下，a′＝a，b′＝b×cos(γ-90)，c′＝c×cos(α-90)，

构成斜方晶系。a、b、c、α、β、γ是指原三斜晶系(也就是TATB晶体)中的6个晶体参数，其中a、b、c是指晶胞的三组棱长，α、β、γ是指三组棱相互间的夹角(即晶体的轴角)，带′的是正交化为斜方晶系后的6个晶体参数。

步骤3、建立伪化学键模型：根据步骤2的斜方晶系建立伪化学键模型，建立一种专一的粗粒化模型；该模型x、y、z三个方向的距离依据步骤2中计算得到的斜方晶系晶格常数确定，然后将每一个材料单体粗粒化为一个珠子，珠子与珠子之间，使用伪化学键连接。伪化学键不是真正意义上的化学键，他的长度远远超出了化学键的长度，且不表示成对电子，而是一种分子间的相互作用，其包含了分子间相互作用、静电力等多种长程力，且不同方向的伪化学键的键长和类型也不同，如x方向的伪化学键长度为a′，类型为①，y方向的伪化学键长度为b′，类型为②，z方向的伪化学键长度为是c′，类型为③。

步骤4、对Morse函数参数进行拟合：对步骤1计算得到的距离-能量矩阵进行Morse函数参数拟合，得到D、α、r₀的参数矩阵；Morse势函数的公式为：

其中E是能量，r是两个分子之间的距离，D、α、r₀是需要拟合的参数。参数拟合使用Python中Scipy包的curve_fit函数，从而得到D、α、r₀的参数矩阵。D、α、r₀这三个参数是Morse函数中的三个参数，其中D是势井深度，r₀是作用势等于0时两个粒子间的距离，α是调和因子。

步骤5、将步骤4获得的参数矩阵和步骤3建立的模型带入MD框架进行模拟；MD模拟框架使用LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)程序。将步骤3建立好的模型和步骤4得到的参数矩阵转化后作为LAMMPS的输入文件输入，然后根据具不同的情况使用不同的系综进行模拟。

MD框架模拟中将键角设置为90°或180°。且MD模拟中，在一个晶体颗粒内，使用伪化学键描述分子之间的相互作用，而在不同的晶体颗粒之间使用Morse对势来描述相互作用。

步骤6、输出模拟结果：输出需要的计算结果或者可视化轨迹。

本发明的优点及有益效果如下：

1、计算速度快。为了检验该方法的计算效率，建立了两个测试模型，分别包含有108个、1728个TATB分子，在MD框架下一共测试3万步(0.03ns)、30万步(0.3ns)、300万步(3ns)；CPU为Intel Xeon E7-4820，使用8个核心进行并行计算。结果如图9、图10，可以看到，在粗粒化模型下运行效率是很高，平均是全原子力场效率的100倍左右。如果体系更大，那么效率提升的就越多。

2、在高计算效率的基础上，有着较高的准确度。如图7、图8，使用本方法得到的数据与实验数据曲线几乎重合，得到的结果与实验结果相差无几，因此有着较高的准确度。

附图说明

图1为QM平面内计算模型示意图；

图2为QM平面间计算模型示意图；

图3为伪化学键模型示意图；

图4为TATB分子粗粒化示意图；

图5为TATB伪化学键粗粒化计算模型；

图6为由机器学习拟合的Morse曲线；

图7为TATB晶体体积随温度变化关系图；

图8为TATB的P-V曲线；

图9、图10为本方法与全原子模型的效率对比。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明做进一步详细说明。

以下说明以含能材料TATB(三氨基三硝基苯)为例。

101、首先对TATB晶体单体进行QM计算，计算模型如图1、2所示，根据不同的晶向，建立不同的模型。计算程序使用Gaussian 09，方法和基组选择cam-B3LYP/6-31+G(d)。

102、一次计算完成后，改变两个分子单体之间的距离，重新计算。多次计算后，得到距离-能量矩阵，如表下表。

距离/Angstrom	能量/Kcal·mol<sup>-1</sup>	距离/Angstrom	能量/Kcal·mol<sup>-1</sup>
				4.63416	6.19	8.154247	56.81
4.711678	4.01	8.280947	32.73
				4.823614	1.31	8.439322	14.52
4.918824	0.06	8.561264	6.38
				4.963248	-0.12	8.629373	3.37
5.001009	-0.33	8.633444	3.19
				5.065078	-0.72	8.768269	-0.56
5.096421	-0.98	8.917987	-2.66
				5.097134	-0.88	8.928059	-2.74
5.100695	-0.81	8.942418	-2.85
				5.134998	-1.12	8.956777	-2.95
5.173588	-1.23	8.957657	-2.95
				5.17909	-1.13	8.985495	-3.12
5.225485	-1.3	9.00073	-3.2
				5.250805	-1.39	9.027723	-3.32
……	……	……	……

201、为了进一步减少计算量和模型复杂度，将所有的晶系正交化为斜方晶系。如TATB晶体为三斜晶系，a＝9.010，b＝9.028，c＝6.812，α＝108.59，β＝91.82，γ＝119.97，根据公式a′＝a，b′＝b×cos(γ-90)，c′＝c×cos(α-90)，α′＝β′＝γ′＝90得到斜方晶系的6个晶格常数：a′＝9.010，b′＝9.028×cos(119.97-90)＝7.82,c′＝6.812×cos(108.59-90)＝6.457,a′＝β′＝γ′＝90。

301、根据在201中计算的6个斜方晶系的晶格常数建立伪化学键模型。伪化学键不是真正意义上的化学键，他的长度远远超出了化学键的长度，且不表示成对电子，而是一种分子间的相互作用，其包含了分子间相互作用、静电力等多种长程力，且不同方向的伪化学键的键长和类型也不同，如图3，x方向的伪化学键长度为a′，类型为①，y方向的伪化学键长度为b′，类型为②，z方向的伪化学键长度为是c′，类型为③。①②③是LAMMPS软件中的1类型键、2类型键、3类型键，也分别指粗粒化模型中xyz三个方向的三种不同的键。如图4是TATB晶体的示意图；使用伪化学键将TATB分子连接起来，表示分子之间的各种相互作用力；然后将一个TATB单体分子进行粗粒化，用一个刚性珠子代替；

302、根据此方法，建立一个含有13824个TATB分子的伪化学键模型，如图5。

401、对102计算得到的距离-能量矩阵进行Morse函数参数拟合，使用Python中Scipy包的curve_fit，代码如下：

拟合结果如图6。这样就得到了D、α、r₀的参数矩阵。

501、MD框架模拟使用LAMMPS程序包。一般的，LAMMPS程序需要两个输入文件：模型文件(DATA_FILE)和命令脚本(Commands Script)。将建立好的伪化学键模型进行转化，生成模型文件，使其能被LAMMPS程序识别。

502、在MD模拟中，由于伪化学键表示的是分子间的相互作用，所以将对势和伪化学键类型都设置为Morse函数，在一个晶体颗粒内，不使用Morse对势而使用伪化学键描述分子之间的相互作用，而在不同的晶体颗粒之间使用Morse对势来描述相互作用。将之前拟合出的D、α、r₀参数矩阵导入命令脚本：

设置伪化学键键角为90°和180°：

angle_coeff 1 20 90

angle_coeff 2 20 180

503、使用NPT系综进行控温控压模拟：fix 1 all npt temp 298 298 100 iso 11 1000。共模拟500000步。

601、使用dump命令输出可视化结果，使用thermo命令输出数据结果。处理数据后即可得到相应的结果。如图7，根据本方法计算得到的TATB的V-T图像，温度为200K-600K，压强为一个标准大气压。根据公式

计算得到热膨胀系数为7.25×10^-5K^-1，与实验值6.6×10^-5K^-1几乎一致。如图8，根据本方法计算得到TATB的高压P-V图，温度设置为298K，压强从0Gpa到10Gpa，每次增加0.2Gpa，横坐标是当前体积与原始体积的比值，纵坐标是压强，单位是GPa，可以看到，使用本方法得到的结果与stevens的实验结果几乎一致，与一些第一性原理算法的结果相差较大，原因可能是第一性原理计算不适合这种大体系的计算。

Claims (5)

1.基于分子间相互作用的粗粒化方法，其特征在于：包括以下步骤：

步骤1、对晶体材料进行量子力学计算：根据晶体材料的晶向对晶体材料的单体分子进行建模，并用QM方法进行计算，计算完成后，改变两个单体分子之间的距离，重新计算；经过一系列的计算后，得到距离-能量矩阵；

步骤2、对三斜晶系进行校正：对于三斜晶系材料，对其晶体结构进行正交化处理，得到斜方晶系；

步骤3、建立伪化学键模型：根据步骤2的斜方晶系建立伪化学键模型，具体为，模型x、y、z三个方向的距离依据校正后的斜方晶系晶格常数确定，然后将每一个晶体材料单体粗粒化为一个珠子，珠子与珠子之间，使用伪化学键连接，伪化学键不是真正意义上的化学键，他的长度远远超出了化学键的长度，表示成一种分子间的相互作用；

步骤4、对Morse函数参数进行拟合：对步骤1计算得到的距离-能量矩阵进行Morse函数参数拟合，得到D、α、r₀的参数矩阵；

步骤5、将步骤4获得的参数矩阵和步骤3建立的模型带入MD框架进行模拟；

步骤6、输出模拟结果：输出需要的计算结果或者可视化轨迹。

2.根据权利要求1所述基于分子间相互作用的粗粒化方法，其特征在于：步骤2所述对于三斜晶系材料的正交化处理为，a′＝a,b′＝b×cos(γ-90),c′＝c×cos(α-90),

a、b、c、α、β、γ是指原三斜晶系中的6个晶体参数，其中a、b、c是指晶胞的三组棱长，α、β、γ是指三组棱相互间的夹角，带′的是正交化为斜方晶系后的6个晶体参数。

3.根据权利要求1所述基于分子间相互作用的粗粒化方法，其特征在于：所述伪化学键的长度远远超出化学键的长度，且不表示成对电子，而表示为一种分子间的相互作用，且不同方向的伪化学键的键长和类型不同。

4.根据权利要求1-3任一项所述基于分子间相互作用的粗粒化方法，其特征在于：步骤5MD框架模拟中将键角设置为90°或180°。

5.根据权利要求4所述基于分子间相互作用的粗粒化方法，其特征在于：在MD模拟中，在一个晶体颗粒内，使用伪化学键描述分子之间的相互作用，而在不同的晶体颗粒之间使用Morse对势来描述相互作用。

Images