内置建模功能¶
当我们迈出进行分子动力学模拟的第一步之前，必须首先对系统进行建模。建模的范围涵盖了分子结构的创建，以及力场参数和系统参数的配置。在这个教程中，我们将专注于体系结构的创建，也就是如何获得一个初始构象。这一步骤对于成功的模拟至关重要，因为它不仅影响了弛豫时间，有时还可能阻碍系统的正常运行。
体系建模的思路可以概括为以下三个关键步骤：规则、模板和优化。规则指的是确定分子如何被排列的方法。模板用于生成多个分子的蓝图。优化阶段依据分子间距等指标来度量体系结构，从而调整以使其更加均匀和自然。例如，生物体系中常用的初始化构象生成工具packmol，它接受PDB格式的分子模板，将模板复制并摆放到指定的区域，然后基于分子间距离进行优化，从而得到合理的初始构象。
举例来说，如果想生成晶体结构，首先需要确定晶格类型和晶格参数，这些规则将成为构建晶体的基础。然后，通过在晶格点上复制和平移摆放分子或原子，可以使用模板生成晶体的初始结构。
在 LAMMPS 中内置了一系列通用的规则和分子模板等功能，可以用来方便地构建常见的初始构象。本教程将按照这个顺序逐一介绍这些关键步骤，理解和应用 LAMMPS 的内置建模工具。
规则¶
LAMMPS中的规则分为几类，如一系列的晶体结构，随即摆放以及基于物理的摆放等等。
让我们从最简单的规则入手，即如何生成一个晶体结构。在这个规则中，我们需要指定一个区域，在这个区域中绘制一个晶格结构，将我们提供的模板复制并摆放到每一个格点上。lammps中，选中一个区域是通过region命令实现的：
region ID style args keywords args ...
ID = user-assigned name for the region
style = delete or block or cone or cylinder or ellipsoid or plane or prism or sphere or union or intersect
keyword = side or units or move or rotate or open
例如，如果我们想选中一个球形区域，命令可以写作region NAME sphere 0 0 0 4.1，即以原点为中心，选中一个半径为4.1的球形区域。同理，可以选择其他的几何形状，例如圆锥圆柱，并且对这些区域进行旋转。甚至甚至可以通过“交并补”的集合操作更加精细地修剪区域形状。
现在，我们选中了一个区域，其名称为NAME。接下来，我们需要在区域内添加格点。lammps内置了多种晶格类型，我们可以通过lattice命令进行操作：
lattice style scale keyword values ...
style = none or sc or bcc or fcc or hcp or diamond or sq or sq2 or hex or custom
scale = scale factor between lattice and simulation box
keyword = origin or orient or spacing or a1 or a2 or a3 or basis
lattice fcc 3.52建立了一个fcc晶格。注意，如果此时体系选择的单位为LJ单位，则约化密度为3.52；如果选择其他单位，3.52则是这种单位制下的晶格常数。
基本的建模规则我们已经定义完成。我们首先选中了一个区域，指定以fcc的晶型向其中添加原子，最后需要做的就是将原子摆放到晶格上。由于我们这里仅仅是向晶格上摆放一个原子，因此可以使用create_atoms 1 region NAME，其中1是指我们摆放的原子类型总数为1。
以上是使用lammps对体系建模的基本流程，我们可以对其中的某些步骤进行替换。如果我们对晶体不感兴趣，我们可以用其他的规则，如：
- create_atoms提供的一系列方法：
- single：在某个坐标点新建一个原子
- mesh：通过(STL)[https://www.adobe.com/creativecloud/file-types/image/vector/stl-file.html#:~:text=The%20name%20STL%20is%20an,a%203D%20model%20or%20object.]文件提供网格
- random: 将分子在空间中随机摆放
- fix pour 模拟倾倒的方式填充体系
- fix deposit 模拟沉积以固定间隔将分子插入到体系中
- fix gcmc 通过蒙特卡洛的方式插入分子
create_atoms提供了将分子随机放入体系的方法。例如，create_atoms 1 random 100 12345 NULL将在体系中随机放入100个类型为1的原子，其中12345是随机数种子，NULL代指在整个体系中随机放置，也可以用region-ID指定特定的区域。如果需要将一个具有多个原子和拓扑连接的分子作为模板放入体系，需要用到molecule命令和分子模板功能。我们将在下一节中详细介绍。
如果已经完成了对一个模拟盒子的建模，可以使用replicate命令对整个盒子进行整体的复制。例如，replicate 2 2 2将整个盒子复制8份，即体系中的原子数与体积将增加8倍。需要注意的是，如果你的体系是周期性的，且分子跨过了边界，此时复制的体系会出现断键的问题。一种行之有效的方式是在随机向盒子中摆放分子时，将分子的初始位置限制在盒子的中心区域，在盒子的四周留出足够的空白区域，从而避免分子跨边界。
create_atoms甚至支持采用STL)格式的文件作为输入。在CAD软件中完成对网格的编辑，生成STL文件，然后通过create_atoms命令将网格导入到lammps中。这种方法在建模复杂的几何形状时非常有用，例如在模拟微流控芯片时，可以直接将芯片的几何形状导入到lammps中，而不需要手动输入每一个原子的坐标。
在data文件中，我们会在文件头的部分指定原子、键、角的数量以及类型数量，之后lammps会在根据这些数量分配相应的内存。但是，当使用create_atoms新建原子时，lammps并不知道我们需要多少内存，更不知道此时模拟盒子的尺寸。因此需要我们通过create_box命令实现：
create_box N region-ID keyword value ...
* N = # of atom types to use in this simulation
* region-ID = ID of region to use as simulation domain
* zero or more keyword/value pairs may be appended
* keyword = bond/types or angle/types or dihedral/types or improper/types or extra/bond/per/atom or extra/angle/per/atom or extra/dihedral/per/atom or extra/improper/per/atom or extra/special/per/atom
除了需要手动指定原子类型的数量N之外，我们还需要指定键、角、二面角等的数量。当我们还需要添加分子时，我们需要额外指定extra/bond/per/atom等参数，lammps将根据这些参数为键接的原子分配额外的内存。例如，使用create_box 1 box可以满足晶体硅的模拟，但是如果我们需要继续添加水分子，则需要使用create_box 2 box extra/bond/per/atom 2 extra/angle/per/atom 1，其中extra/bond/per/atom 2表示每个原子需要额外分配2个键的内存，extra/angle/per/atom 1表示每个原子需要额外分配1个角的内存。当然，我们这里并不需要精确地设置这些数量，我们完全可以给一个足够大的数值。相对于模拟过程中的内存，这里的数量可以忽略不记。
模板¶
上一节中，我们知道可以直接使用create_atoms命令将原子放入体系，那么如何将分子放入体系呢？这就需要用到分子模板的功能。所谓molecule template实际上是与data文件很相近的文本文件，其中记录了一个分子的坐标、原子类型、电荷、拓扑连接等数据。这个文件将通过molecule命令读入，并赋予一个molecule-ID，之后可以在其他命令中引用。例如，create_atoms 1 random 100 12345 mol h2o 123将在体系中随机放入ID为h2o的分子，其中123同样也是随机数种子。
例如，一个长得像水分子模板文件的格式如下：
# water molecule
3 atoms
2 bonds
1 angles
Coords
1 0.000000 0.000000 0.000000
2 0.957200 0.000000 0.000000
3 0.000000 0.957200 0.000000
Types
1 1
2 2
3 2
Charges
1 -1.0
2 0.5
3 0.5
Bonds
1 1 1 2
2 1 1 3
Angles
1 1 2 3 1 2
如果读取多个分子模板，我们需要手动处理type-id的问题。例如，如果我们读取了两个分子，第一个分子的type-id为1到2，第二个分子的type-id为1到4。这样的type-id是有重叠的，因此我们需要手动将type-id进行平移，使得type-id连续。这个过程可以通过offset Toff Boff Aoff Doff Ioff参数实现。例如，molecule mol2 offset 2 2 2 2 2将所有第二个分子的所有type-id都增加2。
优化¶
通过上述非物理的方法将分子塞入体系中所得到的结果是不可信的。使用这种结构直接进行分子动力学模拟通常会出现一系列的问题，例如原子重叠或者距离过远等。我们需要对体系进行优化。最常见的优化方法是minimize。minimize模块是用于执行分子动力学模拟中的结构能量最小化的工具，是通过调整原子坐标来寻找系统的能量最低点或稳定结构，这对于研究分子、晶体、蛋白质和材料等各种体系的性质和稳定性非常重要。minimize模块通常通过在LAMMPS的输入文件中添加以下命令来使用：
minimize [style] [tolerance] [max_iter] [max_eval] [etol] [ftol]
- `[style]`：选择用于能量最小化的优化算法。LAMMPS支持不同的最小化算法，如Conjugate Gradient（共轭梯度法）、L-BFGS（Limited-memory Broyden-Fletcher-Goldfarb-Shanno法）等。
- `[tolerance]`：用于控制最小化过程的收敛的容忍度。当能量变化小于该值时，最小化过程将停止。
- `[max_iter]`：最大迭代次数。如果达到最大迭代次数而未达到收敛条件，则最小化过程将停止。
- `[max_eval]`：最大能量评估次数。如果达到最大评估次数而未达到收敛条件，则最小化过程将停止。
- `[etol]`：能量容忍度。当系统的能量变化小于该值时，最小化过程将停止。
- `[ftol]`：力容忍度。当系统的最大原子力小于该值时，最小化过程将停止。
另一种基于动力学的优化方法是使用 fix nve/limit。这个方法可以与诸如 pair_style soft 等势函数一起使用，其主要作用是限制每个原子在每一步模拟中的移动距离，以防止原子之间的距离缩短得过快，导致产生过大的相互作用力，并防止原子移动距离过大。这样的控制可以确保模拟系统的动力学行为更加稳定。