LAMMPS软件的主要算法

Mr. Material

已于 2023-03-05 11:06:23 修改

阅读量2.4k

点赞数 1

分类专栏： LAMMPS小技巧文章标签： 1024程序员节

于 2021-10-24 22:44:15 首次发布

本文链接：https://blog.csdn.net/qq_43689832/article/details/120943044

版权

LAMMPS小技巧专栏收录该内容

32 篇文章

订阅专栏

本文档介绍了LAMMPS中的空间分配算法，包括处理器间通信策略、数据交换优化，以及针对后处理的MATLAB集成和Python绘图技巧。涵盖了安装教程、实例分析和编程技巧等内容，适合对分子动力学感兴趣的读者。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

关注 $\color{Violet} \rm Mr.material\ ,$ $\color{red}{更}$ $\color{blue}{多}$ $\color{orange}{精}$ $\color{green}{彩}$ ！

主要专栏内容包括：
$\textbf{ \underline{\dag《LAMMPS小技巧》：}}$ 主要介绍采用分子动力学（ $L amm p s$ ）模拟相关安装教程、原理以及模拟小技巧（难度： $\bigstar$ ）
$\textbf{ \underline{\dag\dag《LAMMPS实例教程—In文件详解》：}}$ 主要介绍采用分子动力学（ $L amm p s$ ）模拟相关物理过程模拟。（包含：热导率计算、定压比热容计算，难度： $\bigstar$ $\bigstar$ $\bigstar$ ）
$\textbf{ \underline{\dag\dag\dag《Lammps编程技巧及后处理程序技巧》：}}$ 主要介绍针对分子模拟的动力学过程（轨迹文件）进行后相关的处理分析（需要一定编程能力。难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag《分子动力学后处理集成函数—Matlab》：}}$ 主要介绍针对后处理过程中指定函数，进行包装，方便使用者直接调用（需要一定编程能力，难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag\dag《SCI论文绘图—Python绘图常用模板及技巧》：}}$ 主要介绍针对处理后的数据可视化，并提供对应的绘图模板（需要一定编程能力，难度： $\bigstar$ $\bigstar$ $\bigstar$ $\bigstar$ ）。
$\textbf{ \underline{\dag\dag\dag\dag\dag\dag《分子模拟—Ovito渲染案例教程》：}}$ 主要采用 $\rm Ovito$ 软件，对 $\rm Lammps$ 生成的轨迹文件进行渲染（难度： $\bigstar$ $\bigstar$ ）。

$\color{red}{\textbf{ \underline{专栏说明（订阅后可浏览对应专栏全部博文）：}}}$
$\color{red} 注意：$ 如需只订阅某个单独博文，请联系博主邮箱咨询。 $\rm lammps\_materials@163.com$

$\spadesuit$ $\dag$ 开源后处理集成程序：请关注专栏《LAMMPS后处理——MATLAB子函数合集整理》
$\spadesuit$ $\dag$ $\dag$ 需要付费定制后处理程序请邮件联系： $\rm lammps\_materials@163.com$

LAMMPS软件的主要算法

在大多数情况下，LAMMPS数值积分原子，分子或者宏观颗粒等粒子的牛顿运动方程。这些粒子通过短程或者长程（一般指库仑力）相互作用。为了计算效率，LAMMPS采用邻居列表的方法来追踪某个粒子周围对其产生相互作用的粒子。这些列表根据粒子靠近时的相互排斥的特点进行优化，从而保证系统的某处的局部密度不会过大。在多核并行计算时，LAMMPS采用基于空间分配技术将模拟区域分割为小的三维子区域，然后将每一个子区域分配给一个处理器进行计算。处理器之间通过MPI进行消息通信并且将跨越子区域边界的原子记录为ghost原子以备调用。当原子跨越子区域边界后就会被分配给新的处理器。为了计算某个处理器负责的原子所受的力，该处理仅需要掌握其周围附近原子的位置即可。因此空间分配算法处理器之间的通信是局部的，这与基于原子分配和受力分配算法不同（这两种算法的通信是全局的）。
此处仅介绍三维长方体模拟区域的情形。此时，分配给每个处理的子区域的大小和形状取决于原子数N和处理器个数P以及模拟区域的长宽比。模拟区域三个方向上设定的处理器个数应当使得每个处理器负责子区域的形状越接近正方体越好。这样可以使不同处理器之间的通信量最小。在空间配算法中，通信消耗的时间与模拟子区域的表面积成正比。在空间分配算法中，每个处理器负责两个数据结构，一个存储该处理负责子区域内的N/P个原子，另一个存储周围子区域的原子。在第一个数据结构中，每个处理器存储原子的完整信息，如位置，速度，邻居列表等。这个数据结构以链表的形式存储，这样就允许原子在相邻子区域交换时，可以插入或删除链表中的原子信息数据。第二个数据结构只存储原子的位置。处理器之间在每一时间步的通信可以确保其负责的原子信息是最新的。图2.4.1为处理器与其周围处理器之间的通信示意图。第一步每个处理器与其东西相邻的处理进行通信。处理器2将其区域内落入处理器1中原子受力截断半径的原子位置存入缓存区。如果在东西方向上子区域的长度d小于截断半径r_s，则将处理器2的所有原子存入缓存区。现在每个处理器向西将信息传入相邻处理器（2传到1）并且结构来自向东方向的信息。每个处理器把接受到的信息存入第二个数据结构中。然后，这个过程反过来，每个处理器向东传出信息，接受来自西边的信息。如果d>r_s，那么每个处理器在东西方向上就完成了所需交换的信息（即原子位置）。如果d<r_s，那么东西方向上的信息传递会被重复直到所有落入相关截断半径的原子的信息传递给对应的处理器。然后类似的过程在南北方向上重复，如图2.4.1（b）所示。差别就在于发送给相邻的处理器的信息不仅包含此处理器负责子区域内的原子（在它的第一个数据结构中的原子），而且任何原子被相邻处理所需要的原子位置均需要发送。对于d=r_s，这就相当于把三个盒子的所有原子位置放在一条信息中传递出去了，如图2.4.1（b）所示。最后在信息传递过程在上下方向重复进行。现在原子位置一个平面9个盒子的信息都会被传递给相邻的处理器。

上述空间分配和消息传递算法具有以下几个特点来降低计算过程中所需要的通信。首先当d\geq r_s时，本来一个处理器要与周围26个处理器进行通信，现在仅需要六次数据交换。而且，如果并行计算机是一个超立方配置的，那么模拟区域的分配可以使得所有的六个这些处理器可以直接与中心的处理器联系。这样通信就会快速而且无争用的。其次，当d<r_s时，所需的原子信息需要来自更远处的子区域。这仅会增加少量的数据交换且这些交换也是与相邻的六个处理器进行的。这个特性使得在较小的体系下采用较多的处理器进行计算也会有很好的效率。
第三个特点是数据交换量是最少的。每个处理器所需要的只是落入这个处理器负责子区域的截断半径内的原子。第四处理器所有接受到的原子信息都可以连续的存储在该处理器的第二个数据结构中。除非需要创建拟发送的缓存信息，否则并不需要花时间重排收到的信息。最后，拟发送的缓存信息可以被很快的创建。对两个数据结构的扫描只需要隔几个时间步才进行（因为此时邻居列表会被重建）以确定哪些原子需要被发送。扫描程序会创建需要被通信的原子列表。在后续的几个时间步中，这个列表仍可以采用，从而直接索引需要的原子，快速缓存通信信息。
Fig. 2.4.2为空间分配算法S1的大纲。设子区域z分配给了处理器P_z，其中z从0到P-1。处理器P_z在x_z中存储了该处理器负责的N/P个原子的位置，在f_z中存储了这些原子的受力。步骤（1a）到（1c）是邻居列表的建立。邻居列表建立每个几个时间步才会建立。在其他算法中用于通信的原子列表在每一时间步都需要建立。在步骤1（a）中有原子从子区域z中移动出去，此时需要删除这些原子在子区域z对应的x_z（第一个数据结构）中的信息，并把这些信息存储为缓存信息。这些原子与6个相邻的处理器通过Fig.2.4.1所示的通信方式进行交换。同时，处理器P_z接受到新的原子并把它们添加至x_z中。接着在步骤（1b）中，所有子区域z截断半径内的原子与相邻处理器进行通信，并将通信原子建立列表存储。

当步骤（1a）和（1b）完成后，处理器的两个数据结构都只最新的，此时就可以建立邻居列表了。如果原子i和原子j都在子区域z内，则（ij）原子对只需在列表中存储一次。如果原子i和原子j在两个相邻的子区域中，则两个处理器分别存储原子i和原子j的列表。如果不这样做，那么处理器计算的力必须互相通信。下面会阐述一个修正的算法，避免这样两次计算原子受力。当d<2r_s时，直接对子区域z内的每个原子针对两个数据结构扫描找到其邻居相互作用。当d>2r_s时，对于围绕着子区域z的薄壳，在每个方向上大概有四个或者更多的薄层。此时，装入薄层法构建邻居列表会有更好的效率。在两个数据结构中的原子均被装入尺寸为r_s的区域内。在子区域z内的每个原子周围的区域都会被检测，以查找可能的邻居。
处理器P_z现在就可以基于建立的邻居列表在步骤（2）中计算其负责原子的受力。当发生相互作用的两个原子都在子区域z内时，受力会被存储两次，一次是原子i，一次是原子j。当发生作用的两个原子处在不用的子区域内时，各个处理器会存储其各自负责原子的受力。在受力f_z计算完成后就可以在步骤（4）中更新原子的位置。最后这些更新后的原子位置必须和相邻的处理器进行通信为下一个时间步的操作准备。这发生在步骤（5）采用图2.4.1所示的通信方式并利用之间建立的列表。在这个操作中的数据交换量是受力计算的截断半径和子区域长度相对值的函数。同时，需要指出的是在那些邻居列表已经建立的时间步中，步骤（5）不是一定要执行的。
S1算法中的通信操作发生在步骤1（a），1（b）和（5）中。其中步骤1（b）和（5）的通信是相同的。这些步骤的通信消耗的时间与数据交换量有关。对于步骤（5），如果我们假设原子的密度时均匀的，那么数据量正比于数据所占的体积即，(d+2r_s)^3-d^3。在体积d^3的区域内大于有N/P个原子。有三种情况需要考虑。第一种，当d<r_s时，来自多个周围的子区域的数据均要被考虑且这一操作量级为8r_s^3；第二种，当d\approx r_s时，来自周围26个子区域的数据都要被考虑，此时操作的数据量量级为27N/P；最后当d\gg r_s时，只有在子区域z六个面附近的原子需要被交换。此时交换的数据量量级比例于子区域的表面积，也就是{6r_sN/P}^{2/3}。这三种情况都列在了步骤（5）中。在图2.4.2中，使用了∆来代表在给定N，P和r_s时，三个中的哪一个会被应用。在步骤（1a）中涉及到的数据交换比较少因为并非所有都会跨越子区域的边界。但是这个操作的数据量依旧比例于子区域的表面积，因此也标示其比例于∆。
在步骤1（c）、（2）和（4）中是计算部分。所有这些操作都比例于N/P，除了在1（c）和（2）中需要额外处理存储在第二个数据结构中的来自周围子区域的原子。这些原子的个数比例于∆，因此也它也包含在了这些步骤的比例量级中。在步骤1（c）和（2）中的主要对计算的消耗贡献的量级为N/2P，因为在盒子内部每一对原子的受力只计算和存储一次。当系统的规模也即整个模拟区域比增大时，也即d相对于r_s增大，∆对计算的贡献减小，整体计算的量级趋向于N/2P。此时，每个处理器主要工作用于计算其内部原子的信息，并且与周围的处理器只交换少量的信息用以更新自己边界处的原子信息。
算法S1的一个很重要的特点就是两个数据结构只有隔几个时间步才在重新建立邻居列表时更新。特别的，即使原子移动超出了子区域的边界，它也不会被分配给新的处理器直到步骤（1a）被执行。处理器P_z依旧会正确的计算原子受力，只要满足以下两个条件。第一，在两次邻居列表被构建时，任一个原子都不会移动超过子区域边界d。第二，所有在r_s范围内的临近原子在每一时间步都被更新。另一种方法是讲更新后跨越边界的原子立刻交换给相邻的处理器。这样做的好处是当邻居列表没有构建时，仅需要交换r_c范围的原子。这可以减少通信量，因为r_s>r_c。但是，这样的话被交换原子的邻居列表的信息也需要发送。邻居列表中的信息是原子对存储其相邻原子当地内存位置的索引。如果原子被频繁的发送到新的处理器，那么这些索引就没有意义了。为了克服这一点，我们的操作是给每个原子分配一个全局指针，该指针随着原子在处理器之间交换。全局指针和局部内存之间的对应必须被每个处理器存储在长度为N的向量中。或者，当全局指针在邻居列表中被引用是，这个全局指针必须被排序和检查以找到它正确对应的原子。前一个方案限制模拟体系的规模；后一个方案需要额外的排序和排查。在算法S1中应用Tamayo和Giles的方法可以是代码不那么复杂，并且降低超量的计算和通信。这不会影响大规模的体系并且会提高算法在中等规模下的表现。
利用牛顿第三定律而对S1算法进行改进，改进算法称为S2。如果处理器P_z只从西面，南面和上面方向接受原子而向东面，北面和下面送出原子，这样相互作用只需要被计算一次，即使两个原子在不同的子区域。这要求将计算的力从相反方向送还给负责该原子的处理器，也即在算法中的第三步。这样操作并不会减小通信的数据量，因为原来一半的信息会被通信两次，但是这确实避免了对跨越两个盒子的相互作用的重复计算。关于这个需要指出两点。第一，从整体来看算法S2比S1节省的时间并不多，特别是对较大的体系。因此只有在（1c）和（2）步节省了∆项的时间。第二，对于大规模的体系，空间分配算法可以通过优化单个处理器计算力的效率来提升算法的效率。在矢量机器上，对于计算力和构建邻居列表的过程中需要特别注意数据结构和循环的顺序来获得高效的浮点率。实施S2算法需要额外注意子区域边界和角落处的原子来保证这些原子的受力只被计算一次，而这会妨碍单处理器的优化。
最后一个问题是在空间分配算法中负载平衡很重要。当模拟区域内的原子数密度大致均匀分布时，算法S1的负载才会平衡。如果模拟区域的形状不是长方体，那么很难讲模拟区域均匀的分成多份。此时，更复杂的算法被开发出来，来应对非长方体的模拟区域和原子密度非均匀分布的体系，但是这会使得子区域不是长方体或者以非长方体的方式连接相邻的子区域。这些情况都会使得向子区域分配原子和相邻子区域的通信变得低效和复杂。如果在分子动力学模拟过程中，体系的物理密度会一直发生改变，那么负载平衡问题是复杂的。任何动态的负载平衡都要求额外的计算资源和数据迁移。
空间分配算法是最高效的，但是其效率依赖于模拟的几何形状并且在具体实施时程序也是较复杂的。