专栏首页Bottom-up Material DesignLAMMPS的C++源代码的内部结构

LAMMPS的C++源代码的内部结构

↑分子模拟交流,请点击蓝字关注

LAMMPS作为分子动力学模拟软件,最大的优点在于其开源的C++代码。C++作为典型的面向对象的程序语言,拥有很好的封装性和继承性,便于全世界的研究者按照自己的需求进行代码扩展和补充。本文将对LAMMPS的源代码内进行介绍,并介绍LAMMPS的每一个时间步内是如何进行的。主要内容如下:

1、LAMMPS源代码中的类层级

2、每一个时间步怎样进行?

3、LAMMPS代码扩展


LAMMPS源代码中的类层级

LAMMPS虽然有很多的源文件,但其类层级是很简单的,如下图所示。每一个方框中代表C++中的类封装,每一类对应两个代码文件分别是.cpp源文件和.h头文件。

LAMMPS中的类层级总图

lammps类是整个代码的最高层级,其实例化的对象可供LAMMPS的其他代码多次使用。比如LAMMPS的主程序main.cpp中就简单调用lammps的实例化对象来输入一个in文件。

LAMMPS源程序main.cpp主要代码

蓝色框中的类也为LAMMPS的“顶级类”,但不是真正意义上的顶级类。之所以称之为顶级类,是因为在其他类中也可以随意访问这些类中的函数或变量。在LAMMPS编程时,这里运用了一些小技巧,即建立了一个叫做Pointers的类,蓝色框中的顶级类都为Pointers的继承类。

Pointers类中的指针定义

红色框中标注的类是一些包括虚函数的母类。虚函数是C++中的一种重要函数,简而言之,虚函数中的变量在程序编译时并不能确定具体值,而是在运行中才能确定,需要边运行边编译。LAMMPS中需要让用户选择类型的命令基本上都用的是虚函数的概念。比如fix操作命令的母类下有100多个子类,分别是fix所进行的具体的操作,例如:fix nve, fix shake, fix ave/time等。

其中一个例外是command类,其子类用于配置in文件中各个命令,可以在run命令之前,之后,之间调用。例如CreateBox, Minimize, Run, Velocity都是command母类的子类。

一个时间步内LAMMPS如何工作

深入学习LAMMPS最重要的一步就是理解一个时间步的内部结构。时间步的推进是依靠LAMMPS中的Update类中的Integrate类完成的。Intergrate仍然是一个母类,对用in文件中不同的run_style,它仍有不同的子类。以下就以其中的Verlet子类为例进行描述,另一个子类Respa的配置方式类似。

Verlet类在源代码中的verlet.cpp和verlet.h文件中,该类中的变量和函数实现了速度-verlet积分算法,该算法是进行牛顿运动力学积分兼顾精度和效率的算法。在Verlet类中,承担主要计算任务的是Verlet::run()函数。

Verlet::run()函数中计算力的部分代码(红线为关键代码)

在具体描述Verlet:run()函数之前,先介绍该类中另外几个重要的函数:init(), setup(), force_clear().

init()函数在每一次run之前都会被调用,该函数设置类一些基本的标志变量,这些标志变量是根据用户的一些设置而确定的。

setup()函数在每一次run之前也会被调用。速度积分velocity-Verlet算法在计算每一步的位置和速度的时候需要当前的力,这些例行的计算会在该函数中计算。在该函数中,某些fix操作也会被调用。各种各样的计数器也会在函数中被初始化。

force_clear() 函数初始化原子受力的向量并初始化值为0.

下面具体讨论Verlet::run函数

Verlet::run()函数(非C++代码,简化以方便理解)

ev_set()函数(位于类Integrate中的函数)设置了两个标志位变量eflagvflag,分别代表能量和位力函数的计算标志。这两个标志变量标志了在该时间步内是否进行了能量和位力函数的计算。对于某些时间步,能量和位力函数不需要计算,设置标志位可以增加计算效率。

一个Timestep中LAMMPS的执行结构:

1、首先是fix-->initial_integrate(),这个函数其实是属于Modify类的。Modify类保存了所有Fix的对象和在一个时间步内被调用的顺序。例如 fix nve, fix, nvt, fix npt这些操作都是在每一步的初始完成的。该函数将会完成半步的速度verlet算法,即得到半步的是速度和整步的位置,而整步的速度需要后面一般速度verlet算法来计算

2、函数post_integrate()会执行一些后处理操作,比如粒子被模拟盒子壁反弹等等。

3、属于Neighbor类的decide()函数确定当前时间步中是否需要确定邻域列表。如果不需要,Comm类中的forward_comm()函数会把镜像原子的坐标划分到每个处理器中。如果邻域列表需要重新建立,则会进行一系列建立邻域列表的命令,在上图的if条件语句中可以看到。

4、在建立新的邻域列表前,pre_exchange()函数首先会被调用,该函数对应in文件中一些添加或删除原子的操作fix。

5、Domain类中的pbc()函数此时会进行周期性边界条件的确认,将移出周期性边界外的原子重新移动到模拟盒子内部。注意该函数不是每一步都调用的,只有在邻域列表被建立时,该函数才会被重新调用。这也就是为什么有时导出的原子坐标的文件中会有原子的坐标超出模拟盒子的范围的原因。

6、Domain类中的reser_box()函数会将模拟盒子重新调整,例如采用最小收缩的固定边界条件时,该函数就会被调用。模拟盒子的变化需要镜像原子交流的内部信息(Comm class)和邻域列表网格(Neighbor class)重新更新。Neighbor类中的setup()函数和setup_bins()函数就是更新这些信息的。

7、在模拟盒子变化之后,需要将原子重新划分处理器。Comm类中的exchange()函数会执行该操作,该类中的另一个函数borders()会确认每个处理器中的原子的镜像原子和镜像原子与实际原子之间的对应关系。

8、现在每个处理器都拥有其实际原子和镜像原子,接下来LAMMPS会调Neighbor类中的build()函数建立邻域列表。该操作是通过给所有的实际原子或镜像原子划分网格(bin)实现的,网格的大小一般与力和邻域列表的截断半径有关。网格和邻域列表的建立是LAMMPS能够快速计算原子位置、速度的重要原因。

9、LAMMPS开始计算受力。首先force_clear()函数会把所有原子的力清零。如果在in文件中newton变量的标志开启,也就是原子间的受力满足牛顿第三定律,实际原子和镜像原子的受力都会被计算并累加,然后镜像原子的受力会对应回其实际原子(Comm类中的reverse_comm()命令实现了这一操作),这样可以节省一半的计算量。对势首先被计算。然后是分子内的相互作用,最后是长程库伦力。

10、此时每个原子的受力已经计算,下面就是用户外加的受力。外加的力是由函数post_force()实现的。然后另外半步速度verlet算法会完成计算,更新原子的速度。

11、在一个Timestep的最后,由end_of_step()函数定义的操作将开始进行,主要是一些诊断类型的计算,例如ave/time,ave/spatial等。最后Output类中的write()函数会把得到热力学信息、各个时间步的原子信息和restart文件写出。

LAMMPS代码的扩展

扩展LAMMPS代码最好的方法是在现有的代码中寻找类似功能的源文件和头文件,然后弄清楚其中每一个函数实现了什么功能,例如用户需要写一个新的势函数,就可以参照最简单的势函数源文件Pair_lj_cut.cpp,首先需要进行预编译,使自己写的势函数能够被LAMMPS识别,以下代码即可实现:

#ifdef PAIR_CLASS

PairStyle(lj/cut,PairLJCut)

/*lj/cut即为势函数类型, PairLJCut是用户自己写的类*/

#else

接下来需要一系列函数实现计算:

以上只是对LAMMPS扩展的一个例子,LAMMPS具有很好的可扩展性,具备一定的C++语言基础后,基本上LAMMPS所有的功能都可以被扩展,包括粒子类型、键合类型,compute类型,fix类型,输入in文件命令,输出dump类型,区域region类型,势函数pair_style类型等等。

本文分享自微信公众号 - Bottom2top(gh_a3885137171a)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2018-11-04

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 为什么C语言仍然占据统治地位?

    导读:C语言五十年来一直是软件开发的一种主力语言。本文介绍它在如今的2019年与C++,Java,C#,Go,Rust和Python抗衡的方式。

    华章科技
  • PCL 可视化

    可视化(visualization)是利用计算机图形学和图像处理技术,将数据转换图像在屏幕上显示出来,并进行交互处理的的理论,方法和技术,

    点云PCL博主
  • Sharding-JDBC:单库分表的实现

    通过上面的优化,已经能满足大部分的需求了。只有一种情况需要我们再次进行优化,那就是单表的数量急剧上升,超过了1千万以上,这个时候就要对表进行水平拆分了。

    猿天地
  • Android通过OpenCV和TesserartOCR实时进行识别

    最近一系列的文章都是用Android利用OpenCV NDK的方法通过摄像头实时获取图像进行图像处理,在上一篇《Android使用Tesseract-ocr进行...

    Vaccae
  • 连接两个点云中的字段或数据形成新点云以及Opennni Grabber初识

    (1)学习如何连接两个不同点云为一个点云,进行操作前要确保两个数据集中字段的类型相同和维度相等,同时了解如何连接两个不同点云的字段(例如颜色 法线)这种操作的强...

    点云PCL博主
  • 浅谈.Net反射 7

    emberInfo是一个基类,它包含的是类型的各种成员共有的一组信息。对于字段、属性、方法、事件等具体成员类型来说,它们包含的信息显然是不一样的,因此,在.NE...

    小蜜蜂
  • Spring Cloud Gateway 扩展支持动态限流

    之前分享过 一篇 《Spring Cloud Gateway 原生的接口限流该怎么玩》, 核心是依赖Spring Cloud Gateway 默认提供的限流过...

    冷冷
  • 入门任意一种编程语言所必须的几道习题

      每当学习一门计算机语言,我们也要做一些练习以便逐步熟悉。随着我们对这种编程语言本身支持的抽象手段理解的过程,以下这些问题,基本可以在几乎每门编程语言学习的过...

    窗户
  • 基础扩展 | 18. 静态链表

    通常,我们使用指针来构建链表。然而,对于没有指针的的编程语言来说,还可以使用数组来构建链表。可以让每个数组元素由两个元素项组成:data和cur,data用来存...

    fanjy
  • Android NDK OpenCV级联方式实时进行人脸检测

    前面的文章《Android通过OpenCV和TesserartOCR实时进行识别》我们已经搭好一个利用NDK方式实时处理摄像头数据的程序了,今天我们就在看看Op...

    Vaccae

扫码关注云+社区

领取腾讯云代金券