一种适用于分子晶体从头算分子动力学轨迹的处理方法

摘要

一种可通过编程实现的适用于分子晶体从头算分子动力学轨迹的处理方法，针对体系内被周期性边界条件约束的原子坐标，通过第一性原理计算方法和原子坐标处理算法，还原受约束的原子坐标，消除轨迹中不合理的分子碎片，使得未发生化学反应的分子结构完整，协助科研人员更高效地分析化学反应路径。该方法具有效率高、可靠性好、指导性强等优点，减轻了对国外昂贵商业软件的依赖，提高了分子动力学研究水平。

说明书

技术领域

本发明属于分子动力学轨迹分析方法，具体涉及从头算分子晶体动力学轨迹的处理方法。

背景技术

分子动力学是一类综合数学、物理和化学的分子模拟方法，是一种模拟多体体系平衡性质、传递性质和微观动力学行为细节的技术。通常情况下，研究某些物质在高温和高压等极端条件下的物理化学变化就会采用分子动力学模拟的手段，以此获取极端条件下难以通过实验技术捕捉到的物质响应特性和化学反应过程。

分子动力学模拟依据其算法的不同可以分为基于经典分子力场的分子动力学和从头算分子动力学。前者只考虑多体系统中原子核的运动，而完全忽略了电子的运动和量子效应；后者重点关注电子的量子效应，尤其是Born-Oppenheimer分子动力学更是在每一个动力学步骤中运用迭代法求解Kohn-Sham方程得到模拟体系能量本征值，计算离子受力，通过拉格朗日运动方程得出下一个动力学步骤中原子核的运动速度和坐标。因此，Born-Oppenheimer分子动力学模拟化学反应过程更加准确可信，但也因采用密度泛函计算方法导致计算量十分庞大。近几年，随着计算机硬件的飞速发展，计算能力的快速提升，Born-Oppenheimer分子动力学渐渐受到国内外学者的重视，并广泛用于分子晶体的研究中。然而，在实际的研究过程中，对Born-Oppenheimer分子动力学轨迹的分析存在一些技术瓶颈：

1.原子成断键判断标准不统一。由于Born-Oppenheimer分子动力学计算量大，主流的VASP和CP2K等软件没有在动力学步中计算体系内键级的功能，因此无法采用键级快捷地判断原子间成断键情况，也无法高效地统计分解模拟过程中的分解产物和分子碎片；

2.分子动力学轨迹可视化效果受周期性边界条件约束。从头算分子动力学模拟软件常常采用分数坐标建模，模型内原子坐标受周期性边界条件的约束，超出晶胞的原子被认为处于下一个晶胞中，程序通过周期性坐标变换的方式将该部分原子放置于原始晶胞内，从而导致原始晶胞中出现大量分子碎片或孤立原子。然而，目前能够可视化分子动力学轨迹的开源程序不具有还原原始晶胞内被周期性边界条件约束的原子坐标的功能，从而导致科研人员难以区分分子碎片和孤立原子是由分解反应产生，还是由周期性坐标变换产生。这对科研人员高效地分析分子动力学轨迹产生了极大的阻碍。

3.主流的VASP、CP2K和Quantum Espresso等软件没有配套的分子动力学轨迹处理程序，而具有分子动力学模拟分析功能的Materials Studio和ADF等外国软件价格高昂且不兼容其它计算软件的轨迹文件。

受现有技术瓶颈限制，在采用从头算分子动力学模拟化学反应的研究中，对具体化学反应路径和分解机理的探讨鲜有报道。因此，需要发明一种可通过编程实现的适用于分子晶体从头算分子动力学轨迹的处理方法，能够还原分子动力学轨迹中被周期性边界条件约束的原子坐标，消除轨迹中不合理的分子碎片，使得未发生化学反应的分子结构完整，协助科研人员更高效地分析化学反应路径，更准确地捕捉化学反应机理，减轻对国外昂贵商业软件的依赖，提高分子动力学研究水平。

发明内容

针对现有技术存在的不足，本发明提供可通过编程实现的适用于分子晶体从头算分子动力学轨迹的处理方法，用来还原体系内原子坐标在模拟过程中受周期性边界条件约束而产生的分子碎片，协助科研人员分析动力学轨迹，包括但不限于观察化学反应路径和统计化学分解产物等。

为了实现上述目标，本发明采取如下的技术方案：一种适用于分子晶体从头算分子动力学轨迹的处理方法，针对体系内原子坐标在模拟过程中受周期性边界条件约束进行的坐标变换，其特征在于，所述方法包括以下步骤：

本发明与现有技术相比，具有以下技术优势：

1.一种适用于分子晶体从头算分子动力学轨迹的处理方法，针对体系内原子坐标在模拟过程中受周期性边界条件约束进行的坐标变换，所述方法包括以下步骤：

步骤1，获取晶体内不同种类原子能够两两成键的最大距离L

步骤2，获取分子晶体的晶格长度、分子晶体内的原子种类、每种原子数量、所有原子的分数坐标，将分子晶体内原子进行编号；

步骤3，将分子晶体内所有原子的分数坐标转换为笛卡尔坐标；

步骤4，根据分子晶体内原子间距离和L

步骤5，将能成键的原子归入一个原子组，将多个原子组生成一个原子组集合，进一步统计得到分子晶体内原子组的个数N1和各个原子组中原子的个数，将原子组集合中的原子组按照原子组中原子个数依次排序；

步骤6，将分子晶体内所有原子状态定义为移动或固定，将晶体内所有原子都标记为移动的状态；

步骤7，按照步骤5得到的原子组集合中的原子组的顺序选择一个原子组GA；

步骤8，从步骤7选择的原子组GA中选择一个原子A，检查该原子的状态；

步骤9，如果原子A的状态为固定则回到步骤7选择其他原子组；如果原子A的状态为移动，则需要根据原子A的初始分数坐标判断移动方向；

步骤10，分别根据原子A的移动方向，计算被移动后的原子A的新分数坐标，进一步计算出被移动后原子A的新笛卡尔坐标，根据移动后原子A与晶体内其它原子之间的距离判断是否有其它原子B与移动后的原子A成键；

步骤11，如果有其它原子B与移动后的原子A成键，按照原子A的移动方向，计算出原子A所在原子组GA内所有原子的新分数坐标，并替换掉原始的分数坐标，将GA内所有原子和原子B所在原子组GB内所有原子标记为固定，进入步骤12；

如果没有其它原子B与移动后的原子A成键，则返回步骤10选择原子A另一个移动方向，如果步骤10遍历了所有移动方向后没有其它原子与移动后的原子A成键，则返回步骤8；如果步骤8遍历了原子组GA中所有的原子后没有其它原子与选择并移动后的原子成键，则进入步骤12；

步骤12，循环步骤7～11直至遍历了所有原子组，得到晶体内所有原子新的分数坐标；

步骤13，将得到具有新的分数坐标的晶体内所有原子，按照步骤5进行归类，依次得到新的原子组集合，进一步统计得到新的原子组集合、分子晶体内原子组的个数N2和各个原子组中原子的个数，将新原子组集合中的原子组按照原子组中原子个数依次排序；

步骤14，如果N2＜N1，回到步骤S12，先将步骤S7读取的原子组集合替换成新的原子组集合且将晶体内所有原子都标记为移动的状态，接着按照步骤12的方式循环步骤7～11；如果N2≥N1，则将步骤12得到的晶体内所有原子新的分数坐标按照相应分子动力学轨迹文件的标准格式输出。

进一步地，所述的步骤3具体：晶格长度参数分别为a

优选地，步骤4：判断原子是否两两成键的条件为：分子晶体内两原子之间的距离小于Lmax则能成键，反之不能成键；

步骤5或步骤13，将原子组集合中的原子组按照原子组中原子个数的大小从小到大依次排序。

步骤9中，如果原子A的状态为移动，则需要根据原子A的初始分数坐标判断移动方向，如果原子A的分数坐标小于0.5，则通过分数坐标加1.0的方式获得移动原子的新分数坐标，反之则通过分数坐标减1.0的方式获得移动原子的新分数坐标。

步骤9中，原子A移动方向为将晶体的a轴、b轴和c轴三种方向进行移动方式的组合，即a轴动，b轴c轴不动、b轴动，a轴c轴不动、c轴动，a轴b轴不动、a轴b轴动，c轴不动、a轴c轴动，b轴不动、b轴c轴动，a轴不动、a轴b轴c轴同时动。

本发明与现有技术相比，具有以下技术效果：

1.本发明采用的原子分组法和多重循环算法大幅度提高了分子动力学轨迹的处理效率，有效地消除了周期性边界条件对轨迹可视化的负面影响；

2.本发明形成的处理方法能够协助科研人员高效地分析从头算分子动力学轨迹，包括但不限于观察化学反应路径和统计化学分解产物等，减轻对国外昂贵商业软件的依赖，提高分子动力学研究水平。

附图说明

图1是实施例1中高氯酸三碳酰肼合锌的分子动力学轨迹处理前后的可视化效果对比，其中图1a为实施例1处理前可视化轨迹图，图1b为实施例1处理后可视化轨迹图。

图2是实施例2中高氯酸三碳酰肼合锌的分解产物轨迹处理前后的可视化效果对比，图2a为实施例2处理前可视化轨迹图，图2b为实施例2处理后可视化轨迹图。

图3是实施例3中碳酰肼的分子动力学轨迹处理前后的可视化效果对比，图3a为实施例3处理前可视化轨迹图，图3b为实施例3处理后可视化轨迹图。

具体实施方式

下面结合实施例对本发明做进一步说明，但本发明不受下列实施例的限制。

实施例1

本实施例选用的是过渡金属配合物高氯酸三碳酰肼合锌(Zn(CHZ)

步骤1，获取GTX分子晶体内不同种类原子能够两两成键的最大距离L

步骤2，获取GTX分子晶体的晶格长度、分子晶体内的原子种类、每种原子数量、所有原子的分数坐标，将分子晶体内原子按顺序编号；GTX分子晶体的晶格长度分别为a

步骤3，将GTX分子晶体内所有原子的分数坐标转换为笛卡尔坐标；以上一步选择的Zn原子为例，将上一步的晶格长度和Zn原子的分数坐标(A＝0.39485585，B＝0.18463756，C＝0.62613754)带入公式X＝a

步骤4，根据GTX分子晶体内原子间距离和L

步骤5，将能成键的原子归入一个原子组，将多个原子组生成一个原子组集合，进一步统计得到分子晶体内原子组的个数N1和各个原子组中原子的个数，将原子组集合中的原子组按照原子组中原子个数依次排序；本实施例中在处理前原子组集合个数N1＝67；

步骤6，将GTX分子晶体内所有原子状态定义为移动或固定，将GTX分子晶体内所有原子都标记为移动的状态；

步骤7，按照步骤5得到的原子组集合中的原子组的顺序选择一个原子组GA；本实施例以只有一个N原子的原子组GA为例；

步骤8，从步骤7选择的原子组GA中选择一个原子A，检查该原子的状态；本实施例的原子组GA中只有一个N原子，即为原子A，其分数坐标为(0.04973235，0.01250552，0.25003133)；

步骤9，如果原子A的状态为固定则回到步骤7选择其他原子组；如果原子A的状态为移动，则需要根据原子A的初始分数坐标判断移动方向；步骤8中选择的N原子状态为移动；

步骤10，分别根据原子A的移动方向，计算被移动后的原子A的新分数坐标，进一步计算出被移动后原子A的新笛卡尔坐标，根据移动后原子A与晶体内其它原子之间的距离判断是否有其它原子B与移动后的原子A成键；步骤8中选择的N原子沿a轴方向的分数坐标0.04973235，则该N原子可沿a轴方向将移动至分数坐标为1.04973235处；类似地，该N原子可沿b轴方向将移动至分数坐标1.01250552处，可沿c轴方向将移动至分数坐标1.25003133处；在上述移动条件的约束下，移动方向为：a轴动，b轴c轴不动、b轴动，a轴c轴不动、c轴动，a轴b轴不动、a轴b轴动，c轴不动、a轴c轴动，b轴不动、b轴c轴动，a轴不动、a轴b轴c轴同时动。

步骤11，如果有其它原子B与移动后的原子A成键，按照原子A的移动方向，计算出原子A所在原子组GA内所有原子的新分数坐标，并替换掉原始的分数坐标，将GA内所有原子和原子B所在原子组GB内所有原子标记为固定，进入步骤12；

如果没有其它原子B与移动后的原子A成键，则返回步骤10选择原子A另一个移动方向，如果步骤10遍历了所有移动方向后没有其它原子与移动后的原子A成键，则返回步骤8，如果步骤8遍历了原子组GA中所有的原子后没有其它原子与选择并移动后的原子成键，则进入步骤12。本实施例中，步骤8中选择的N原子在移动到分数坐标为(0.04973235，1.01250552，0.25003133)的位置时与分数坐标为(0.11244351，0.97758895，0.20575611)的C原子成键，此时将N原子和C原子所在的原子组GB内所有原子标记为固定；

步骤12，循环步骤7～11直至遍历了所有原子组，得到晶体内所有原子新的分数坐标；

步骤13，将得到具有新的分数坐标的晶体内所有原子，按照步骤5进行归类，依次得到新的原子组集合，进一步统计得到新的原子组集合、分子晶体内原子组的个数N2和各个原子组中原子的个数，将新原子组集合中的原子组按照原子组中原子个数依次排序；本实施例中在第一次处理至步骤13时，GTX分子晶体内原子组的个数N2＝26，原子组集合个数相比于处理前大幅度减少，说明大量残缺的GTX分子逐渐被补全；

步骤14，如果N2＜N1，回到步骤S12，先将步骤S7读取的原子组集合替换成新的原子组集合且将晶体内所有原子都标记为移动的状态，接着按照步骤12的方式循环步骤7～11；如果N2≥N1，则将步骤12得到的晶体内所有原子新的分数坐标按照相应分子动力学轨迹文件的标准格式输出。本实施例中，第一次处理至步骤13时，N2＝26

用ASE(Atomic Simulation Environment)程序将轨迹文件可视化(图1)，从前后的可视化效果对比结果来看，处理前没有完整的GTX分子，处理后该分子结构完整且没有不合理的分子碎片。通过实施例1可以说明本发明能够有效地解决周期性边界条件对轨迹可视化的负面影响，协助科研人员高效地分析动力学轨迹。

实施例2

本实施例选用的是过渡金属配合物高氯酸三碳酰肼合锌(Zn(CHZ)

步骤1，获取GTX分子晶体内不同种类原子能够两两成键的最大距离L

步骤2，获取GTX分子晶体的晶格长度、分子晶体内的原子种类、每种原子数量、所有原子的分数坐标，将分子晶体内原子按顺序编号；GTX分子晶体的晶格长度分别为a

步骤3，将GTX分子晶体内所有原子的分数坐标转换为笛卡尔坐标；以上一步选择的Zn原子为例，将上一步的晶格长度和Zn原子的分数坐标(A＝0.61562229，B＝0.92577598，C＝0.66034684)带入公式X＝a

步骤4，根据GTX分子晶体内原子间距离和L

步骤5，将能成键的原子归入一个原子组，将多个原子组生成一个原子组集合，进一步统计得到分子晶体内原子组的个数N1和各个原子组中原子的个数，将原子组集合中的原子组按照原子组中原子个数依次排序；本实施例中在处理前原子组集合个数N1＝143；

步骤6，将GTX分子晶体内所有原子状态定义为移动或固定，将GTX分子晶体内所有原子都标记为移动的状态；

步骤7，按照步骤5得到的原子组集合中的原子组的顺序选择一个原子组GA；本实施例以只有一个N原子的原子组GA为例；

步骤8，从步骤7选择的原子组GA中选择一个原子A，检查该原子的状态；本实施例的原子组GA中只有一个N原子，即为原子A，其分数坐标为(0.06651380，0.83617123，0.99387827)；

步骤9，如果原子A的状态为固定则回到步骤7选择其他原子组；如果原子A的状态为移动，则需要根据原子A的初始分数坐标判断移动方向；步骤8中选择的N原子状态为移动；

步骤10，分别根据原子A的移动方向，计算被移动后的原子A的新分数坐标，进一步计算出被移动后原子A的新笛卡尔坐标，根据移动后原子A与晶体内其它原子之间的距离判断是否有其它原子B与移动后的原子A成键；步骤8中选择的N原子沿a轴方向的分数坐标为0.06651380，则该N原子可沿a轴方向将移动至分数坐标为1.06651380处；相反地，该N原子沿b轴方向的分数坐标为0.83617123，则该N原子可沿b轴方向将移动至分数坐标为-0.16382877处；类似地，该N原子可沿c轴方向将移动至分数坐标-0.00612173处；在上述移动条件的约束下，移动方向为：a轴动，b轴c轴不动、b轴动，a轴c轴不动、c轴动，a轴b轴不动、a轴b轴动，c轴不动、a轴c轴动，b轴不动、b轴c轴动，a轴不动、a轴b轴c轴同时动。

步骤11，如果有其它原子B与移动后的原子A成键，按照原子A的移动方向，计算出原子A所在原子组GA内所有原子的新分数坐标，并替换掉原始的分数坐标，将GA内所有原子和原子B所在原子组GB内所有原子标记为固定，进入步骤12；

如果没有其它原子B与移动后的原子A成键，则返回步骤10选择原子A另一个移动方向，如果步骤10遍历了所有移动方向后没有其它原子与移动后的原子A成键，则返回步骤8，如果步骤8遍历了原子组GA中所有的原子后没有其它原子与选择并移动后的原子成键，则进入步骤12。本实施例中，步骤8中选择的N原子(原子A)在移动到分数坐标为(0.06651380，0.83617123，-0.00612173)的位置时与分数坐标为(0.09804023，0.84921393，0.04690389)的N原子(原子B)成键，此时将N原子(原子A)和N原子(原子B)所在的原子组GB内所有原子标记为固定，两个N原子实则为一个氮气分子(N

步骤12，循环步骤7～11直至遍历了所有原子组，得到晶体内所有原子新的分数坐标；

步骤13，将得到具有新的分数坐标的晶体内所有原子，按照步骤5进行归类，依次得到新的原子组集合，进一步统计得到新的原子组集合、分子晶体内原子组的个数N2和各个原子组中原子的个数，将新原子组集合中的原子组按照原子组中原子个数依次排序；本实施例中在第一次处理至步骤13时，GTX分子晶体内原子组的个数N2＝110，原子组集合个数相比于处理前略有减少，这说明模型中分解产物多为零散的小分子化合物；

步骤14，如果N2＜N1，回到步骤S12，先将步骤S7读取的原子组集合替换成新的原子组集合且将晶体内所有原子都标记为移动的状态，接着按照步骤12的方式循环步骤7～11；如果N2≥N1，则将步骤12得到的晶体内所有原子新的分数坐标按照相应分子动力学轨迹文件的标准格式输出。本实施例中，第一次处理至步骤13时，N2＝110

用ASE(Atomic Simulation Environment)程序将轨迹文件可视化(图2)，从前后的可视化效果对比结果来看，处理前分解产物仍然是分子碎片，这导致统计分解产物的种类和数量会遇到极大的困难；采用本发明的方法处理后，每个分解产物的结构均完整，所有的非常见结构的分子碎片均是分解产物或是分解产物的聚合物。通过实施例2可以说明本发明能够协助科研人员高效地分析从头算分子动力学轨迹，包括但不限于观察化学反应路径和统计化学分解产物等。

实施例3

本实施例选用的是碳酰肼分子(H

步骤1，获取CHZ分子晶体内不同种类原子能够两两成键的最大距离L

步骤2，获取CHZ分子晶体的晶格长度、分子晶体内的原子种类、每种原子数量、所有原子的分数坐标，将分子晶体内原子按顺序编号；CHZ分子晶体的晶格长度分别为a

步骤3，将CHZ分子晶体内所有原子的分数坐标转换为笛卡尔坐标；以上一步选择的O原子为例，将上一步的晶格长度和O原子的分数坐标(A＝0.07775692，B＝0.22337864，C＝0.89482989)带入公式X＝a

步骤4，根据CHZ分子晶体内原子间距离和L

步骤5，将能成键的原子归入一个原子组，将多个原子组生成一个原子组集合，进一步统计得到分子晶体内原子组的个数N1和各个原子组中原子的个数，将原子组集合中的原子组按照原子组中原子个数依次排序；本实施例中在处理前原子组集合个数N1＝40；

步骤6，将CHZ分子晶体内所有原子状态定义为移动或固定，将CHZ分子晶体内所有原子都标记为移动的状态；

步骤7，按照步骤5得到的原子组集合中的原子组的顺序选择一个原子组GA；本实施例以有2个N原子和2个H原子的原子组GA为例；

步骤8，从步骤7选择的原子组GA中选择一个原子A，检查该原子的状态；本实施例的原子组GA中的一个N原子，其分数坐标为(0.09578715，0.03920234，0.19603486)；

步骤9，如果原子A的状态为固定则回到步骤7选择其他原子组；如果原子A的状态为移动，则需要根据原子A的初始分数坐标判断移动方向；步骤8中选择的N原子状态为移动；

步骤10，分别根据原子A的移动方向，计算被移动后的原子A的新分数坐标，进一步计算出被移动后原子A的新笛卡尔坐标，根据移动后原子A与晶体内其它原子之间的距离判断是否有其它原子B与移动后的原子A成键；步骤8中选择的N原子沿a轴方向的分数坐标0.09578715，则该N原子可沿a轴方向将移动至分数坐标为1.09578715处；类似地，该N原子可沿b轴和c轴方向将移动至分数坐标1.03920234和1.19603486处；在上述移动条件的约束下，移动方向为：a轴动，b轴c轴不动、b轴动，a轴c轴不动、c轴动，a轴b轴不动、a轴b轴动，c轴不动、a轴c轴动，b轴不动、b轴c轴动，a轴不动、a轴b轴c轴同时动。

步骤11，如果有其它原子B与移动后的原子A成键，按照原子A的移动方向，计算出原子A所在原子组GA内所有原子的新分数坐标，并替换掉原始的分数坐标，将GA内所有原子和原子B所在原子组GB内所有原子标记为固定，进入步骤12；

如果没有其它原子B与移动后的原子A成键，则返回步骤10选择原子A另一个移动方向，如果步骤10遍历了所有移动方向后没有其它原子与移动后的原子A成键，则返回步骤8；如果步骤8遍历了原子组GA中所有的原子后没有其它原子与选择并移动后的原子成键，则进入步骤12。

本实施例中，步骤8中选择的N原子在移动在所有可能移动的位置上都未找到与之成键的原子，则返回步骤8选择一个新的N原子，新选择的N原子坐标为(0.06565819，0.11710711，0.08599817)。新选择的N原子移动到分数坐标为(0.06565819，0.11710711，1.08599817)的位置时与分数坐标为(0.10228022，0.14103637，0.99350166)的C原子(原子B)成键。此时，按照新选择的N原子的移动方向，计算出GA内所有原子的新分数坐标，并替换掉原始的分数坐标，将GA内所有原子和原子B所在原子组GB内所有原子标记为固定。

步骤12，循环步骤7～11直至遍历了所有原子组，得到晶体内所有原子新的分数坐标；

步骤13，将得到具有新的分数坐标的晶体内所有原子，按照步骤5进行归类，依次得到新的原子组集合，进一步统计得到新的原子组集合、分子晶体内原子组的个数N2和各个原子组中原子的个数，将新原子组集合中的原子组按照原子组中原子个数依次排序；本实施例中在第一次处理至步骤13时，CHZ分子晶体内原子组的个数N2＝26，原子组集合个数减少，残缺的CHZ分子逐渐被补全；

步骤14，如果N2＜N1，回到步骤S12，先将步骤S7读取的原子组集合替换成新的原子组集合且将晶体内所有原子都标记为移动的状态，接着按照步骤12的方式循环步骤7～11；如果N2≥N1，则将步骤12得到的晶体内所有原子新的分数坐标按照相应分子动力学轨迹文件的标准格式输出。本实施例中，第一次处理至步骤13时，N2＝26

用ASE(Atomic Simulation Environment)程序将轨迹文件可视化(图1)，从前后的可视化效果对比结果来看，处理前没有完整的CHZ分子，处理后该分子结构完整，说明本发明能够有效地消除周期性边界条件对轨迹可视化的负面影响。

尽管这里参照本发明的多个解释性实施例对本发明进行了描述，但是，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。