lammps的helloworld：in文件基础——分子动力学仿真及其在离子阱中运用学习与实践分享第2.1

基态清原子

已于 2024-03-15 00:42:10 修改

阅读量982

点赞数 19

分类专栏：分子动力学文章标签：学习

于 2024-03-15 00:29:32 首次发布

本文链接：https://blog.csdn.net/peace__mind/article/details/136518922

版权

分子动力学专栏收录该内容

5 篇文章 1 订阅

订阅专栏

基本调用方式

lammps软件只是一个求解器，所有的命令只能通过代码的方式输入到求解器进行求解计算。输入和输出都是文件形式，这就是为什么lammps官网上也给出了一些预处理/后处理程序。它们用于生成和处理这些文件。

lammps的运行形式是逐行读取并运行一个脚本文件，这个文件一般被称为in文件。

in文件简单语法

在合适的位置建立一个空文件in.helloworld，开始我们的第一次尝试吧。这个文件的基本结构分为四个部分，后面会进行逐段分析，首先学习几条基本的语法。

变量定义

lammps的变量定义有多种形式。其基本语法格式为：

variable name style args ...

其中的style包含“delete or atomfile or file or format or getenv or index or internal or loop or python or string or timer or uloop or universe or world or equal or vector or atom”多种，类似于其它语言变量定义中的变量类型但又略有区别，主要体现于这里的”variable“实在是非常广义的。其中最基本最常用的，是“equal“和”index“，其它类型或许会在后续进一步深入时遇到吧。

equal是将equal后面的表达式赋给前面定义的变量，equal是一种变量的类型，是由数学公式计算得出的结果，换言之就是数值。在网页的说明上，也给出了许多lammps当中支持的数学函数。例如可以写：

variable a equal ln(exp(1))

需要注意的是，equal定义的变量不会立即计算，而是在每次使用到该变量的时候再进行计算，也就是说这里定义的本质是一个函数或者计算法则。与其类似的还有atom、vector和python类型。
字符串也是常见的数据类型，需要使用string定义，定义格式是类似的，例如：

variable outstring string "hello world"

而index使用一个或者多个字符串定义，给出的是类似于python中list的一个列表，可以用于逐个调用。例如：

variable a index 1 2 3 45

index还有其他的妙用，后面还会提到它。

字符串与其输出

print函数是最基础的输出函数，其默认是向屏幕输出信息，也可以向文件中输出。其基本的结构为：

print string keyword value

其中keyworld和value可以没有。
在lammps当中字符串十分灵活，变量可以灵活的嵌入字符串中。变量的引用方法是$ 加上变量名即 $v 或者$ {value}。
根据上面的语法，我们就可以写出：

第一版hello world

完整代码如下：

# in.helloworld
variable outstring string "hello world"
print "${outstring}"

这里需要注意的细节是print的内容需要包含在双引号之内，因为直接引用变量时它的原理相当于将原先的字符串抄写在此处了，并没有双引号，print的调用会出问题。
在命令行终端中运行命令：

lmp -in in.helloworld

即可在终端中看到输出，其中包含代码中要求输出的hello world字符，其显示类似于：

LAMMPS (24 Dec 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94)
  using 1 OpenMP thread(s) per MPI task
hello world
Total wall time: 0:00:00

我们的lammps编程旅程正式开始！

条件结构

lammps当中的条件结构关键词为if、then、elif、else。其基本的语法与其它的高级语言十分相似：

if boolean then t1 t2 ... elif boolean f1 f2 ... elif boolean f1 f2 ... else e1 e2 ...

这里需要注意的点包括;

条件满足时的执行命令每一条都需要被包括在双引号中
这是一条单独的语句，因此换行时，需要在行末加上&来表示该语句没有结束

循环结构

lammps是采用其它语句中已经被“不推荐使用”的label-jump语法来实现循环结构的。
循环结构的本质就是一段代码重复跑好几遍，因此jump就是通过不包装语句，而是在运行中直接改变运行顺序，使程序从某一（前面的）位置开始运行的方式，实现循环的。label的作用是作为jump的目标，本质就是个传送锚点。
其一般的结构应当写作：

label loop
	variable  i loop $N
	~~~~~~
	if "BOOL" then "jump in.script break"
	next i
	jump 
label break
variable i delete

对于上述代码，做一些简单的注解：

这里的"loop"和"break"并没有特殊性，可以将它们换成合适的名字。
注意then后面的jump语句需要包裹在双引号里面
变量的定义类型是loop。实际上这里也可以使用前文提到过的index，另有一种uloop。这些变量定义类型不会重复定义（在删除前）。因此在这个模板中，将参数变量定义在loop标签前面还是后面在一般情况下都是可以的。同时我们注意到由于index等类型不能重复定义，与equal类型对比，它完全可以被当作常量类型使用。
变量有$X和v_X两种引用模式，详情可见这篇博客。在下面的程序当中涉及这一问题。

第二版hello world

完整代码如下：

# in.helloworld
variable outstring string "hello world"
variable a index 1 2 3 45 67
variable sum equal 1
label loop
	print "${outstring} ${a} times"
	variable sum equal ${sum}+$a
	if "${sum}>10" then "jump in.helloworld break"
	next a
	jump in.helloworld loop
label break
variable a delete

可以得到循环运行的结果，参考输出为：

LAMMPS (24 Dec 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94)
  using 1 OpenMP thread(s) per MPI task
hello world 1 times
hello world 2 times
hello world 3 times
hello world 45 times
Total wall time: 0:00:00

最后一个输出是45，之后程序退出。

Hello world a trillion times！

in文件基本结构

接下来，正式进入仿真实用的in文件内容。一般的，in文件可以分为模型初始化、系统初始化、仿真设置、仿真运行四个部分。

模型初始化

单位制

第一行通常是单位制的定义，根据官网相关说明，lammps给出了一些常用的单位制，这些单位制影响了输入输出的单位以及很多内部参数的值，在计算机计算的舍入误差累积效应下，这种偏差导致的影响常常是不可忽略的，需要尽可能选择合适的单位制并保持连续一致。
当然，我们现在是helloworld，当然选择最通用而明晰的单位制：国际单位制，也就是第一行代码：

unis si

维度

lammps可以运行2维或3维的模型。我们采用最常见的3维，这是默认值，实际上不需要写出。

dimension 3

边界类型

boundary命令用于定义仿真盒子的边界类型。在仿真过程中，无论是什么系综，粒子的运动范围不会是无限的，这个命令给出了在6个方向上的边界上会受到的约束。
其中p即周期性边界条件，意味着边界是周期性的，在这一方向穿越边界，会从对方向返回，其它三种fsm都是非周期性的。如果xyz某个方向两侧的边界不同，则采用连续的两个字母定义。
我们是helloworld嘛，采用最简单的全周期条件，即：

boundary p p p

原子属性及设定

atom_style是一个用于指定用于选择 LAMMPS 模拟中哪些原子属性与原子相关联的命令。
不同的模型需要访问特定的原子属性，例如要计算库仑相互作用，原子必须具有 “电荷”（又称 “q”）属性，这些属性需要与这些原子一起存储和通信。为了节约计算资源，选择合适的原子属性十分重要。有许多不同的原子样式可以将属性组合在一起。有些原子样式是其他原子样式的超集。之后也可以通过混合样式（提供子样式属性的组合）或fix property/atom命令为原子添加更多属性。
原子属性集有很多种类型，适用于不同的问题，它们之间的区别可见于官网介绍。在第二列，介绍了它们之间的包含关系，而最后一列写出了它们对应用于求解的问题。第三列给出的包要求是对于lammps的安装build过程的，如果遇到问题，可以检查一下。
如果用于离子体系，那么显然就是"atomic systems with charges"了，选择charge，即：

atom_style charge

而atom_modify命令用于修改一些 LAMMPS 中定义和存储的原子的某些属性。这些属性大多与具体的计算过程有关，在初学时大概还不需要修改。

其它属性

还包括pair_style, bond_style, angle_style, dihedral_style, improper_style等，在此先不予考虑。

newton

该命令用于打开或关闭成对和成键相互作用的牛顿第三定律。对大多数问题而言，将牛顿第三定律设置为开启意味着计算量的适度节省，但通信量要增加两倍。这是否更快取决于问题的大小、力截止长度、机器的计算/通信比率以及使用的处理器数量。将 "成对牛顿 "标志设置为 "关闭 "意味着，如果两个相互作用的原子位于不同的处理器上，则两个处理器都会计算它们的相互作用，而不会通信所产生的力信息。同样，对于键相互作用，关闭牛顿标志意味着如果一个键、角、二面体或不恰当的相互作用包含 2 个或更多处理器上的原子，则相互作用由每个处理器计算。
除了舍入问题外，LAMMPS 在任何牛顿标志设置下都会产生相同的答案。在运行方式为 respa 且仅在最内层时间步中计算成键相互作用（键、角等）的情况下，关闭成键相互作用的 newton 标志可能会更快，以避免最内层循环中的额外通信。
由于我们目前的问题是小问题，因此根据其说明我想当然的认为关闭会更好一点，并且似乎区别不大：

newton off

processors

与处理器相关，暂时不必考虑

系统初始化

这一部分用于建立用于模拟的初始原子系统。模拟必须有一个起点，用下面的方式可以告知系统所涉及的粒子和它们之间（如果有）的相互作用、化学键等的初始状态。从此状态开始，系统进行运算。

内置初始化操作

region函数（建立区域）

region函数会给出一个集合区域定义。其具有如下基础语法：

region ID style args keyword arg ...

ID一项是一个自定义的字符串，不一定需要是数字。第二个style是区域形状类型，可选择的选项包括：“block or cone or cylinder or ellipsoid or plane or prism or sphere or delete or union or intersect”。
上述最后3个属于功能性操作，删除和取交集并集。其它则是可以选择的区域形状。其具体参数需要参考官方参考文档。
keyword包括side or units or move or rotate or open，是对于定义的进一步修改，可选。

最简单的是block，即立方体形状。它只需要6个参数，也就是x，y，z坐标的下界和上界。例如：

region theregion block -5 5 -5 5 0 2

create_box函数（建立模拟盒子）

create_box函数的基本语法为：

create_box N region-ID keyword value ...

其前两个参数是选取的原子类型编号和region创建的集合区域ID，而后面是对于其它内容（如化学键等）进行补充定义。
我们helloworld程序不希望过于复杂，因此采用最为基本的，即例：

create_box 1 theregion

lattice函数（建立晶格）

lattice函数只负责生成晶格，并不填充原子。

其语句基本结构为：
lattice style scale keyword values …
style是点阵的类型，包括了none 、sc 、 bcc、 fcc 、hcp、 diamond、 sq 、sq2、 hex和custom
scale是网格与模拟框之间的比例因子
keyword是后面追加的关键字，这些关键字包括 origin、 orient、 spacing、 a1、 a2、 a3、 basis，它们用于定制晶格参数。

为了后面的helloworld简单，我们采用bcc晶格，并进行最简单的初始化，追加关键字暂不予考虑，则采用语句：

lattice bcc 2

create_atom函数（创建原子）

在前面的基础上，我们就可以使用create_atom命令创建原子了。
其基本语法为：

create_atoms type style args keyword values ...

这里type是原子的type编号，从1开始的数字。style是创建原子方式的类型，包括“box or region or single or mesh or random”，也就是使用不同的方式产生原子。而keyword相关是进一步的细致设置。具体的可见于官网教程，这里不展开了。

这里先采用最基础的方式，即填充整个区域的方式：

create_atoms 1 region theregion ratio 0.5 23333

ratio 0.5代表随机填充一半的格子，后面跟着的是随机种子。

原子质量

顾名思义，mass，给定原子质量。我们采用我们碳基生物喜欢的碳原子做粒子，因此：

mass 1 12

文件初始化操作

LAMMPS使用data文件来定义系统的初始状态和模拟参数。
data文件是一个文本文件，包含了描述模拟系统的各种信息，包括原子类型、原子坐标、分子拓扑、力场参数等。第一行通常是描述系统的概要信息，比如模拟步数、原子数、边界条件等。接下来的行包含了系统中每个原子的信息。每一行用空格或制表符分隔，包含了原子的唯一标识号、原子类型、原子坐标（x、y、z轴坐标）以及可能的额外属性。之后的行定义了系统中的分子拓扑，如键连接、键角等。这些信息可以通过定义每个原子的连接列表来指定原子之间的相互作用。最后几行是力场参数的定义，包括原子类型的质量、电荷等。这些参数将影响模拟中的力场计算。

上面提到的系统定义和初始化方式过于基础，因此不能用于更大的系统。对于更大的系统，一定需要使用文件操作。我们首先可以看一些成熟的data文件作为格式的参考，最简单的是输出上述定义对应的data文件：

data的输出

到目前为止，我们在程序中已经采用的语句包括：

# part1 Initialization
unis si
dimension 3
boundary p p p
atom_style charge
newton off

# part2 System definition
region theregion block -5 5 -5 5 0 2
create_box 1 theregion
lattice bcc 2
create_atoms 1 region theregion ratio 0.5 23333
mass 1 12

将当前的仿真状态输出为data文件的函数为write_data，例如：

write_data data.helloworld

输出的data文件如下：

LAMMPS data file via write_data, version 24 Dec 2020, timestep = 0

25 atoms
1 atom types

-5 5 xlo xhi
-5 5 ylo yhi
0 2 zlo zhi

Masses

1 12

Atoms # charge

1 1 0 -1 -5 1 0 0 0
2 1 0 3 -5 1 0 0 0
3 1 0 -3 -3 1 0 0 0
4 1 0 -2 -4 0 0 0 0
5 1 0 -1 -3 1 0 0 0
6 1 0 3 -3 1 0 0 0
7 1 0 4 -4 0 0 0 0
8 1 0 -4 -2 0 0 0 0
9 1 0 -2 -2 0 0 0 0
10 1 0 0 -2 0 0 0 0
11 1 0 1 -1 1 0 0 0
12 1 0 3 -1 1 0 0 0
13 1 0 -1 1 1 0 0 0
14 1 0 0 0 0 0 0 0
15 1 0 2 0 0 0 0 0
16 1 0 3 1 1 0 0 0
17 1 0 4 0 0 0 0 0
18 1 0 -4 2 0 0 0 0
19 1 0 -3 3 1 0 0 0
20 1 0 -1 3 1 0 0 0
21 1 0 2 2 0 0 0 0
22 1 0 4 2 0 0 0 0
23 1 0 -4 4 0 0 0 0
24 1 0 2 4 0 0 0 0
25 1 0 4 4 0 0 0 0


Velocities

1 0 0 0
2 0 0 0
3 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0
7 0 0 0
8 0 0 0
9 0 0 0
10 0 0 0
11 0 0 0
12 0 0 0
13 0 0 0
14 0 0 0
15 0 0 0
16 0 0 0
17 0 0 0
18 0 0 0
19 0 0 0
20 0 0 0
21 0 0 0
22 0 0 0
23 0 0 0
24 0 0 0
25 0 0 0

可以看到包括了仿真区域的信息，原子的基本信息和所有原子的位置和速度信息。
data文件的格式内容丰富，将在下一篇文章里专门学习，可以通过官方教程的相关内容学习了解。

读入data文件

我们可以通过其它方式生成data文件，并读入data文件以建立原子系统的模型。由于实际的系统往往很大，不使用其它工具而采用手动方式建立模型是困难的。

读入data文件采用的函数是read_data，其基本的语法是：

read_data file keyword args ...

data文件的进一步理解将写在下一篇，这里简单略过。

重新启动操作

可以通过read_restart读入重启文件来重启和继续仿真进程。这些文件可以通过restart命令从运行过程中保存。

仿真设置

力场设置

模型建好之后，需要设置力场参数，也就是常说的势函数设置。最基本的是两两之间非键结相互作用，对应的是pair_style。

从物理模型角度，力场无疑是分子动力学仿真中最重要的部分，包含了对于整个系统物理性质的基本假设。因此选择合适、正确的力场模型至关重要。lammps当中给出了非常多的预设模型可供选用，可见于此链接。
在其中还包括一个python类型，这意味着可以通过python代码的方式自定义力场的数学形式。、

对于一般的中性原子分子，常常采用lj类型的力场，而我们这里使用的是离子，库伦场占据主导地位，因此采用库伦场，对应的是coul。而coul/cut，coul/debye等选项是对原先的库伦公式进行了各种各样的修改，以符合现实的多种干扰因素。我们选用最简单的coul/cut，其基本形式为：

pair_style coul/cut cutoff

每一种力场都需要相关的参数进一步给定，这个过程用到的命令是pair_coeff，至少会存在的内容是这种力场作用与哪些原子之间，如果作用于所有原子之间，可以用*代替。因此这里可以采用的代码：

pair_style coul/cut 3.0
pair_coeff * *

在此之上，还有很多种相互作用，包括化学键的作用，三个或者四个原子之间的作用类型，其定义方式与两个原子间的力场定义是相似的。这是可选的，因此此处先不选。

约束条件设置

在lammps当中，约束条件使用fix命令给出。fix的类型非常多样，默认给出的选项包括这个列表链接。如果不再需要某个fix，可以使用unfix命令撤销。

其中较为基础的一类是系综条件，其物理意义可以见本系列序号1.2的文章。在helloworld阶段，我们就采用微正则系综即可，可以写作：

fix 1 all nve

其中1是我们为这一条约束所起的名字，all的位置是原子组的编号，我们这里没有进行分组，因此写all。

仿真计算设置

compute命令可以进行单独的计算，其基本语法为：

compute ID group-ID style args

例如：

compute 1 all temp

仿真选项设置

timestep

这里给出仿真时每一步对应的时间。这个参量一般需要设定，但也存在默认值，默认值与所选单位制直接相关。
国际单位制的默认值为1e-8s即10ns。

neibor

分子动力学仿真中，为了计算量的减少，通常不会计算所有原子之间的相互作用，而只会计算“相邻”的原子。这个相邻所对应的范围就需要由neibor给出。其基本语法为：

neighbor skin style

skin对应着能够被算作“相邻”的距离对应的范围，而style是建立邻居表的方式。如果在计算上没有特殊的要求，一般选择bin即可。这里同样有默认值，skin = 0.001即1mm，style默认为bin。

内容输出设置

可以默认地输出热力学量，语法为：

thermo N

即每隔N个时间步输出一次热力学量。这种输出默认是输出到命令行的，其输出格式等可以通过thermo_style, thermo_modify进行修改。这里就不详细展开了。

dump也会输出信息，但通常用于向文件夹中输出更加详细的信息，包括绘制图像等。它的设置可以非常复杂，但如果只是想看到原子某一时刻的位置，并不复杂，只需要：

dump d0 all image 100 dump.*.jpg type type

仿真运行

run！
最后一部分，我们就只需要run了，run后面跟的参数是run的步数。

第三版hello world

完整版的代码如下，此处做了一些修改（晶格改成了全填充，方便图像的辨识）：

# part1 Initialization
unis si
dimension 3
boundary p p p
atom_style charge
newton off

# part2 System definition
region theregion block -5 5 -5 5 0 2
create_box 1 theregion
lattice bcc 2
create_atoms 1 region theregion ratio 1 23333
mass 1 12

# part3 Simulation settings
pair_style coul/cut 3.0
pair_coeff * *
fix 1 all nve
thermo 10
dump d0 all image 50 dump.*.jpg type type

# part4 Run a simulation
run 100

很显然，我们这里选用的是国际单位制，因此回头看前面的定义，我们实际上定义了一个10米乘10米乘2米的空间，我们在这个或许可以在北京卖一个亿的空间中生成了50个原子，这真是非常理想的真空，原子之间的距离也远大于cutoff，它们根本就是不动的。主要还是提示自己和大家——注意单位！

你是故意的还是不小心的
我是，，故意不小心的

虽然不动，但至少这里我们看到了原子，我们在运行时会输出三张图，它们是：
在这里插入图片描述

根据图片也可以验证：它们确实一点没动。因此它们的温度也将始终为0

在过程中，还输出了内容，其内容列于下方，供参考，可以看到温度和能量确实是0：

LAMMPS (24 Dec 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:94)
  using 1 OpenMP thread(s) per MPI task
Created orthogonal box = (-5.0000000 -5.0000000 0.0000000) to (5.0000000 5.0000000 2.0000000)
  1 by 1 by 1 MPI processor grid
Lattice spacing in x,y,z = 2.0000000 2.0000000 2.0000000
Created 50 atoms
  create_atoms CPU = 0.000 seconds
System init for write_data ...
WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (src/comm_brick.cpp:167)
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 3.001
  ghost atom cutoff = 3.001
  binsize = 1.5005, bins = 7 7 2
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair coul/cut, perpetual
      attributes: half, newton off
      pair build: half/bin/newtoff
      stencil: half/bin/3d/newtoff
      bin: standard
Setting up Verlet run ...
  Unit style    : si
  Current step  : 0
  Time step     : 1e-08
Per MPI rank memory allocation (min/avg/max) = 3.449 | 3.449 | 3.449 Mbytes
Step Temp E_pair E_mol TotEng Press 
       0            0            0            0            0            0 
      10            0            0            0            0            0 
      20            0            0            0            0            0 
      30            0            0            0            0            0 
      40            0            0            0            0            0 
      50            0            0            0            0            0 
      60            0            0            0            0            0 
      70            0            0            0            0            0 
      80            0            0            0            0            0 
      90            0            0            0            0            0 
     100            0            0            0            0            0 
Loop time of 0.00509082 on 1 procs for 100 steps with 50 atoms

93.3% CPU use with 1 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.00047795 | 0.00047795 | 0.00047795 |   0.0 |  9.39
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 4.9218e-05 | 4.9218e-05 | 4.9218e-05 |   0.0 |  0.97
Output  | 0.0044977  | 0.0044977  | 0.0044977  |   0.0 | 88.35
Modify  | 3.1709e-05 | 3.1709e-05 | 3.1709e-05 |   0.0 |  0.62
Other   |            | 3.42e-05   |            |       |  0.67

Nlocal:        50.0000 ave          50 max          50 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:        594.000 ave         594 max         594 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:        1075.00 ave        1075 max        1075 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 1075
Ave neighs/atom = 21.500000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00