嵌入原子势计算优化方法与流程

本发明设计大规模分子动力学模拟中EAM势及相应力的计算优化技术，特别涉及一种嵌入原子势计算优化方法。

背景技术：

在生物、化学、材料科学及计算机科学交叉领域，分子动力学模拟方法扮演着举足轻重的角色。而在分子动力学的模拟过程中，原子(或者分子)的势能及相互之间作用力的计算，占据了模拟的绝大部分时间。对此部分计算优化，能大幅度提高模拟效率。EAM势考虑粒子间的多体相互作用，包括两部分：一是粒子间的相互作用对势、二是粒子镶嵌在电子云中的嵌入能。

现有技术中，经典的计算方法需要两层循环，第一层循环遍历所有的粒子，以计算其电子云密度，第二层循环计算总能量、作用力及其他的信息。

技术实现要素：

本发明的目的在于，简化分子动力学模拟中的计算热点：粒子的势及粒子间相互作用力的计算，以此来提高模拟效率。以EAM势，尤其是其中的嵌入能计算部分作为优化对象。

为实现上述目的，本发明提供了一种嵌入原子势计算优化方法，嵌入原子势计算优化方法包括：通过初始化设置第一至少一个时间步及第一至少一个时间步对应的电子云密度；根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度；根据第一至少一个时间步及第一至少一个时间步对应的电子云密度和第二时间步及第二时间步对应的电子云密度，确定嵌入原子势。

优选地，根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的步骤之后还包括：记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数。

优选地，记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数的步骤之后还包括：根据预先设定的更新频率和记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数判断；当记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数为预先设定的更新频率的正整数倍时，更新邻居列表，设定第三至少一个时间步及第三至少一个时间步对应的电子云密度，并由第三至少一个时间步及第三至少一个时间步对应的电子云密度，确定第四时间步及第四时间步对应的电子云密度；当记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数不为预先设定的更新频率的正整数倍时，根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度。

优选地，根据预先设定的更新频率和记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数判断的步骤之后还包括：根据记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数与待模拟的时间步的个数进行判断；当记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数小于待模拟的时间步的个数时，根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度。

优选地，根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度：

ρ＝l_m，i(t)ρ_i+l_m，i+1(t)ρ_i+1+l_m，i+2(t)ρ_i+2+…+l_m，i+n(t)ρ_i+n

其中，t_i，t_i+1，...t_i+n为前第一时间步，ρ_i、ρ_i+1、ρ_i+2、ρ_i+n前第一时间步对应的电子云密度，m为插值多项式的项数。

优选地，根据第一至少一个时间步及第一至少一个时间步对应的电子云密度和第二时间步及第二时间步对应的电子云密度，确定嵌入原子势为：

其中，e_i为第i个粒子的对势，ρ_i为第i个粒子对应的电子云密度，F(ρi)为嵌入能。

优选地，ρ_i为：

ρ_i＝∑_i≠jf_ij(r_ij)

其中，r_ij项指原子i和原子j之间的距离，f_ij(r_ij)是第j个粒子核外电子在第i个粒子处贡献的电荷密度。

优选地，e_i为：

其中，Φ_ij(r_ij)为原子i和原子j之间的对势。

本发明通过简化动力学模拟中的计算热点：势及粒子间相互作用力的计算，以此来提高模拟效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种嵌入原子势计算优化方法的流程图；

图2为图1中的一个具体的实施例；

图3为BCC铁原子的电子云密度随时间变化图。

具体实施方式

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

图1为本发明实施例提供的一种嵌入原子势计算优化方法的流程图。如图1所示，嵌入原子势计算优化方法包括以下步骤：

步骤S101：通过初始化设置第一至少一个时间步及第一至少一个时间步对应的电子云密度；

步骤S102：根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二至少一个时间步及第二至少一个时间步对应的电子云密度；

根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二至少一个时间步及第二至少一个时间步对应的电子云密度为：

ρ＝l_m，i(t)ρ_i+l_m，i+1(t)ρ_i+1+l_m，i+2(t)ρ_i+2+…+l_m，i+n(t)ρ_i+n

其中，t_i，t_i+1，…t_i+n为时间步，ρ_i、ρ_i+1、ρ_i+2、ρ_i+n时间步对应的电子云密度。

需要说明的是，l_m，i(t)中的m为插值多项式的项数，可以采用二次插值多项式l_2，i(t)，也可以采用三次插值多项式l_3，i(t)，采用几次插值多项式也是与拟合分子有关，在此不赘述。

步骤S102中通过根据第一至少一个时间步及第一至少一个时间步对应的电子云密度，使用抛物插值多项式得到第二至少一个时间步对应的电子云密度。此处的抛物插值多项式只是最优的方法，能实现同等效果的算法都可以替代抛物插值多项式，在此不赘述。此处的第一至少一个时间步和第二时间步，只表示第二时间步是由第一至少一个时间步外推得到后面的时间步。

步骤S103：根据第一至少一个时间步及第一至少一个时间步对应的电子云密度和第二时间步及第二时间步对应的电子云密度，确定嵌入原子势。

具体地，根据第一至少一个时间步及第一至少一个时间步对应的电子云密度和第二时间步及第二时间步对应的电子云密度，确定嵌入原子势为：

其中，e是对势(pair potential)，F是嵌入能(embedding energy.)，电子云密度(electron density)ρ是周围原子作用的和。

ρ_i＝∑_i≠jf_ij(r_ij) (2-3)

其中，e_i为第i个粒子的对势，ρ_i为第i个粒子对应的电子云密度，F(ρi)为嵌入能；r_ij项指原子i和原子j之间的距离，f_ij(r_ij)是第j个粒子核外电子在第i个粒子处贡献的电荷密度；Φ_ij(r_ij)为原子i和原子j之间的对势。

可选的，由于拟合的过程中，原子的轨迹会在多次执行步骤S102发生变化，没有按照原先的规律。因此，根据执行几次步骤S102的规律，更新邻居列表，重新设定第三至少一个时间步。

另外，为了能够判断是否满足更新邻居列表的要求，可以考虑设置一个计数器，需模拟的时间步的个数为计数器的周期，在每一次执行步骤S102，计数器加一，直到计数器的值为预先设定更新频率的正整数倍时，更新邻居列表，也可以通过程序来设置。具体的实施方式比较灵活，本发明实施例对此不做限定。该计数器用于记录执行步骤S102的次数。

可选地，在步骤S102之后，嵌入原子势计算优化方法还包括：

步骤S104：记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数。记录执行步骤S102的次数，在初始化的时候可以初始化为零。

通过计数器记录执行步骤S102的次数和预先设定更新频率判断需不需要对邻居列表进行更新。

具体地，嵌入原子势计算优化方法还包括：

步骤S105：根据记录执行步骤S102的次数和预先设定的更新频率判断；

当记录执行步骤S102的次数为预先设定的更新频率的正整数倍时，更新邻居列表。

更新完邻居列表，由重新设定的第三至少一个时间步及第三至少一个时间步对应的电子云密度外推得到后面的时间步及时间步对应的电子云密度；也就是说每更新一次邻居列表，需要重新设定时间步及时间步对应的电子云密度。

当记录执行步骤S102的次数不为预先设定的更新频率的正整数倍时，继续执行步骤S102，由第一至少一个时间步及第一至少一个时间步对应的电子云密度外推得到后面时间步对应的电子云密度。

在具体实施模拟的时候，会预先设定一个需拟合的最大时间步的个数，通过记录执行步骤S102的次数和需模拟的最大时间步的次数进行比较，只在记录执行步骤S102的次数小于需模拟的最大时间步的次数时，继续执行由第一至少一个时间步及第一至少一个时间步对应的电子云密度外推得到后面时间步及时间步对应的电子云密度；当满足更新邻居列表的频率为记录执行步骤S102的次数的正整数倍时，更新邻居列表。

在邻居列表更新完后，重新设定至少一个时间步及至少一个时间步对应的电子云密度，并由重新设定的至少一个时间步及重新设定的至少一个时间步对应的电子云密度外推得到后面的时间步对应的电子云密度。

具体地，嵌入原子势计算优化方法还包括：

步骤S106：根据记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数与需模拟时间步的个数判断；

当记录执行步骤S102的次数小于待模拟时间步的个数，继续执行由第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度。

在记录执行步骤S102的次数小于待模拟时间步的个数时，继续执行步骤S102，由第一至少一个时间步及第一至少一个时间步对应的电子云密度外推得到后面时间步及时间步对应的电子云密度；当记录执行步骤S102的次数满足为预先设定的更新频率的正整数倍时，更新邻居列表。

在更新完邻居列表，需要设定第三至少一个时间步及第三至少一个时间步对应的电子云密度。

当记录执行步骤S102的次数大于或等于待模拟时间步的个数时，结束第一至少一个时间步及第一至少一个时间步对应的电子云密度外推得到后面时间步对应的电子云密度的操作，由外推得到的所有的电子云密度，得到嵌入原子势。

需要说明的是，对于步骤S106还需要设置在步骤S104之后，和设置在步骤S105之后，且步骤S106设置在步骤S104之后，与设置在步骤S105后相互独立，都是根据执行步骤S102的次数和需模拟的时间步的个数进行判断；当执行步骤S102的次数大于或等于需模拟的时间步的个数时，由外推得到的所有的时间步及所有时间步对应的电子云密度，得到嵌入原子势。

需要说明的是，第一次执行步骤S102之前，也需要根据记录执行根据设置的第一至少一个时间步及第一至少一个时间步对应的电子云密度，确定第二时间步及第二时间步对应的电子云密度的次数与需模拟时间步的个数判断。

步骤S106和步骤S103处计算嵌入原子势的公式相同，只是设定的外推起点对应的三个时间步及前三个时间步对应的电子云密度不相同，在此处不赘述。

图2为图1中的一个具体的实施例。如图2所示，图2为BCC铁的级联碰撞分子动力学示意图。

图2中的是一个粒子的电子云密度随时间变化为例，采用前三个时间步的电子云密度外推得到后面时间步的电子云密度，每执行5次外推得到后面时间步的电子云密度的操作，更新邻居列表；更新邻居列表后，需重新设定外推启动点，得到后面时间步对应的电子云密度，而不是沿用没更新之前的数据。

具体的，嵌入原子势计算优化方法的步骤包括：

步骤S201：通过初始化设置前三个时间步及前三个时间步对应的电子云密度；

如，前三个时间步分别为t₀、t₁、t₂，前三个时间步对应的电子云密度分别为ρ₀、ρ₁、ρ₂。

步骤S202：根据设置的前三个时间步及前三个时间步对应的电子云密度，确定前三个时间步后的时间步及前三个时间步后的时间步对应的电子云密度；

步骤S202使用抛物插值多项式得到当前时间步对应的电子云密度。此处的抛物插值多项式只是最优的方法，能实现同等效果的算法都可以替代抛物插值多项式，在此不赘述。

基于前三个时间步t₀、t₁、t₂，及前三个时间步分别对应的电子云密度为ρ₀、ρ₁、ρ₂。根据抛物插值算法构建关于前三个时间步t₀、t₁、t₂的二次二次插值基函数l_2，0(t)、l_2，1(t)、l_2，2(t)：

公式(1)、(2)和(3)需要满足公式(4)的要求。

然后，通过点(t₀，ρ₀)、(t₁，ρ₁)、(t₂，ρ₂)的二次插值多项式得到前三个时间步后的时间步对应的电子云密度ρ₃为：

ρ₃＝l_2，0(t)ρ₀+l_2，1(t)ρ₁+l_2，2(t)ρ₂公式(5)

t₀、t₁、t₂时间步后的对应电子云密度由公式(5)外推得到。

步骤S203：根据前三个时间步及前三个时间步对应的电子云密度和前三个时间步后的时间步及前三个时间步后的时间步对应的电子云密度，确定嵌入原子势。

根据前三个时间步后的时间步及前三个时间步后的时间步对应的电子云密度，确定当前时间步对应电子的嵌入原子势为：

其中，e是对势(pair potential)，F是嵌入能(embedding energy.)，电子云密度(electron density)ρ是周围原子作用的和。

ρ_i＝∑_i≠jf_ij(r_ij) (2-3)

其中，e_i为第i个粒子的对势，ρ_i为第i个粒子对应的电子云密度，F(ρi)为嵌入能；r_ij项指原子i和原子j之间的距离，f_ij(r_ij)是第j个粒子核外电子在第i个粒子处贡献的电荷密度；Φ_ij(r_ij)为原子i和原子j之间的对势。

本发明实施例中的嵌入原子势为EAM(Embedded-Atom Method)势，其他的嵌入原子势若计算规则如EAM势的计算规则相同，则认定为相同，在此不赘述。

可选的，由于拟合的过程中，原子的轨迹会在多次执行步骤S202发生变化，没有按照原先的规律。因此，根据执行几次步骤S202的规律，更新邻居列表，重新设定前第二至少一个时间步。

另外，为了能够判断是否满足更新邻居列表的要求，设置一个计数器，需模拟的时间步的个数为计数器的周期，在每一次执行步骤S202，计数器加一，直到更新邻居列表的频率(如5)为计数器的值的倍数(如，5,10,15等)，更新邻居列表。也可以通过编程控制更新邻居列表，具体的实施方式比较灵活，本发明实施例对此不做限定。

具体地，执行完步骤S202后，嵌入原子势计算优化方法还包括：

步骤S204：记录执行根据设置的前三个时间步及前三个时间步对应的电子云密度，确定前三个时间步后的时间步及前三个时间步后的时间步对应的电子云密度的次数；

通过计数器记录执行步骤S202的次数和更新邻居列表的频率(预先设定的更新频率)判断需不需要对邻居列表进行更新。

具体地，嵌入原子势计算优化方法还包括：

步骤S205：根据记录执行步骤S102的次数和预先设定的更新频率判断；

当记录执行步骤S102的次数为更新邻居列表的频率(如5)的正整数倍(如，5、10、15等)时，更新邻居列表。

在更新完邻居列表后，重新设定三个时间步及三个时间步对应的电子云密度。并由重新设置的三个时间步及三个时间步对应的电子云密度外推得到后面时间步对应的电子云密度。并在满足邻居列表更新频率是对之前设定的三个时间步及三个时间步对应的电子云密度重新设定。

当记录执行步骤S202的次数不为更新邻居列表的频率(如5)的正整数倍(如，5、10、15等)时，按照之前设定的三个时间步及三个时间步对应的电子云密度外推得到后面的时间步对应的电子云密度。

需要说明的是，每执行一次由三个时间步及三个时间步对应的电子云密度外推得到后面时间步及时间步对应的电子云密度，计数器都需要加1。

在具体实施模拟的时候，会预先设定一个需拟合的最大时间步的个数，通过计数器的数值和需模拟的最大时间步的次数进行比较，只在计数器的数值小于需模拟的最大时间步的次数时，继续执行由三个时间步及三个时间步对应的电子云密度外推得到后面时间步及时间步对应的电子云密度；当满足更新邻居列表的频率(如5)为计数器的数值的倍数时，更新邻居列表。

在邻居列表更新完后，重新设定三个时间步及三个时间步对应的电子云密度，并由重新设定的三个时间步及重新设定的三个时间步对应的电子云密度外推得到后面的时间步对应的电子云密度。

具体地，嵌入原子势计算优化方法还包括：

步骤S206：根据记录执行步骤S202的次数与待模拟时间步的个数判断；

当记录执行步骤S202的次数小于待模拟时间步的个数，继续执行由三个时间步及三个时间步对应的电子云密度外推得到后面时间步及时间步对应的电子云密度；当计数器的数值满足为更新邻居列表的整数倍时，更新邻居列表。

在邻居列表更新完后，重新设定新的三个时间步及新的三个时间步对应的电子云密度。

当记录执行步骤S202的次数大于或等于待模拟时间步的个数时，结束由三个时间步及三个时间步对应的电子云密度外推得到后面时间步对应的电子云密度的操作，由外推得到的所有的电子云密度，得到嵌入原子势。

需要说明的是，在每一次执行完步骤S202、步骤S204之后，需要判断将记录执行步骤S202的次数和待模拟的时间步的个数，因此，步骤S206也需要设置在步骤S204之后，与设置在步骤S205后相互独立。

需要说明的是，第一次执行步骤S202之前，也需要根据记录执行步骤S202的次数与待模拟时间步的个数判断。

步骤S206和步骤S203处计算嵌入原子势的公式相同，只是设定的外推起点对应的三个时间步及前三个时间步对应的电子云密度不相同，在此处不赘述。

图3为BCC铁原子的电子云密度随时间变化图。如图3所示，BCC铁的级联碰撞分子动力学模拟中，以某一原子的电子云密度随时间步变化来说明。图3中每隔5步曲线出现拐点，在出现拐点之前用前三个时间步对应的电子云密度外推得到前三个时间步后的时间步对应的电子云密度。出现拐点的原因是由于，更新了邻居列表。因此，需要在邻居列表更新后，需要重新设置新的前三个时间步及行的前三个时间步对应的电子云密度，使得根据新的前三个时间步及行的前三个时间步对应的电子云密度进行外推，而不是沿用之前的数据。具体更新邻居列表的频率，根据材料相关。

本发明通过简化动力学模拟中的计算热点：势及粒子间相互作用力的计算，以此来提高模拟效率。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。