lammps_jiayou的博客

fix recenter命令可以把设置的原子组pka调整到后面三个坐标点为止，同时，整个模型也会自动地进行偏移。最近几天有几个粉丝咨询如何把随机选择的辐照PKA原子设置到模型的中心位置，本文介绍一下这个方法。后面的参数0 0 0为模型的中心坐标，这个也需要根据自己的模型进行修改。686为data文件或者模型中的最大原子数，可根据自己的模型修改。这种建模，有两个关键点：随机选择和原子居中。蓝色的pka原子已经位于模型的中心位置。大家好，我是小马老师。

lammps reaxff反应势可以模拟分子间的化学反应，fix reaxff/species命令可以输出反应过程中的产物信息。在ovito中，查看总原子数，也可以看出，随着模拟过程的进行，总原子数越来越少，说明此命令确实能够动态删除产物分子。在2022之后的lammps版本中，这个命令新增了delete关键词，使用这个命令，可以动态删除反应过程中的分子。filedel：删除的分子信息保存到该文件中。大家好，我是小马老师。

在新版的lammps命令中，set命令有了较多的扩展，比如可以对group和region进行原子类型的替换。晶界偏析这篇推文是2021年写的，当时的set替换命令功能有限，当时只能对type类型进行替换。在之前的案例中，介绍过set命令随机替换原子进行高熵合金建模，也介绍过晶界偏析的建模方法。把中间区域mid中的type 1按比例替换为type 2。把top组的type 1原子全部替换为type 2。把type 1中的100个原子替换为type 2。因此，本文汇总最新的set替换命令，以供参考。

（2）弹出旋转对话框，第一行设置旋转轴，分别对应x、y和z轴，设置z轴为1，设置Angle角度为180，第三行设置旋转中心，保存不变。在ovito后处理中，如果需要对模型进行旋转，可以使用ovito的“Affine transformation"。（4）选择菜单，导出data文件，会发现新的data文件结构已经实现了沿z轴的180度旋转。以摩擦模型为例，球在x轴左侧，下面使用这个命令把结构沿z轴旋转180度。使用这个方法，也可以实现流体模型上下壁面的对称建模，有兴趣的可以试试。

分别统计红色小球周围蓝色原子和红色原子的数量，需要用到count(type,region)命令，改命令可统计区域region内type类型原子的数量。以小球的随机碰撞为例，原子类型为1的红色小球（id也为1）周围分布原子类型为2、3的蓝色和黄色小球。

先设置中间区域mid，调用fcm()命令计算，最后用thermo命令输出结果。在前面的推文中，介绍过lammps动态统计区域内原子数量、区域压力的方法。这个函数可以计算区域region内group原子组在dim方向的受力和。dim：力的方向，可设置x、y、z。大家好，我是小马老师。group：原子组id。region：区域id。

两个方法都可以实现区域压力的计算，第一个方法可以计算区域压力的分布，如果需要同时计算多个区域，可以采用第一个方法，如果只需要计算一个特定区域，第二个方法更灵活一些。在lammps模拟中，计算某一个固定区域内气体或者液态的压力，可以先计算该区域内所有单个原子的应力，然后把区域内原子的应力值求和再除以体积。如下面代码，在x方向划分了坐标为[20,40]的一个区域，计算完成后，会之间输出该区域的压力值。只有处于该区域内的原子参与压力的计算，当原子移动出该区域后，就不再计算该原子的应力。大家好，我是小马老师。

在上一次案例代码的基础上，稍加修改，添加一个for循环遍历所有的帧即可实现多帧孔洞体积的计算。在前面的专栏中，已经介绍了单帧孔洞体积的计算方法，有不少粉丝朋友咨询多帧孔洞体积的计算方法。本文介绍lammps后处理技巧：多帧孔洞体积和孔隙率的计算方法。计算的结果保存到Vol.txt文本中，便于绘图。

这种方法只能统计一个静态的区域内原子的数量变化，如果要统计某一个运动中的原子周围原子数量的变化，就需要根据该原子的运动状态，以小球的随机碰撞为例，原子类型为1的红色小球（id也为1）周围分布原子类型为2、3的蓝色和黄色小球。在之前的专栏中，曾介绍过动态统计某一个固定区域内原子数量的方法。根据这个原理，也可以统计该原子周围不同原子类型的数量。

ovito中有一个Construct Surface mesh命令，可以识别原子表面，使用这个功能，就可以直观的显示空洞的大小和形态。ovito免费版可以直观的显示孔洞的表面积，但是不能计算体积，计算体积需要ovito pro版本。这个命令还可以输出其他的量，如下图所示，具体可到ovito官网查看手册说明。在拉伸、辐照或者一些力学的lammps模拟中，材料内部会出现孔洞等缺陷。但是，使用python调用ovito模块，则可以很方便的计算孔洞的体积。

在lammps中，和MC相关的模拟命令有三个：fix gcmc、fix atom/swap、fix sgcmc。在一个lammps模拟中，mc和md可同时存在，间隔运行，所以设置fix sgcmc命令的同时，可以同时设置系综（nve、nvt或者npt）。lammps官网自带案例中也有一个fix sgcmc的案例，主要模拟FeCu合金的MC建模。sgcmc命令可应用到合金的模拟中，如短程有序模拟。MC运行间隔为100步，最终的浓度比例为0.02，也就是2%的Fe原子交换为Cu原子。大家好，我是小马老师。

上式中的compass_published对应lammps自带的力场文件compass_published.frc，windows版本的lammps中，这个文件在lammps安装目录的frc_files文件夹，ubuntu版本的lammps中，这个文件在安装目录的tools/msi2lmp/frc_files,如果以上命令不能用，先查找一下这个文件是否存在。compass力场在ms中用的比较多，但是在lammps中用的较少，因为很多参数不是公开的，限制了这个力场的应用范围。大家好，我是小马老师。

主要原理和之前的锐角刀具建模原理一样，用region plane设置6个面，用这个6个面组成一个封闭区域，在封闭区域内生成原子。lammps自带的形状有长方体、圆柱体、圆锥体和球体，其他特殊的情况大多需要用组合建模方法。如下图截面为六边形的纳米线，很难用一条命令建模。大家好，我是小马老师。

在lammps模拟时，因为键被拉的过长，会造成“bond atoms missing”错误，这种错误在聚合物模拟中经常会出现，这个错误的解决也是聚合物模拟的一个难点。主要原因是，上述结构单独弛豫过程中，当原子移动到边界处，因周期性边界原因，会被算法折叠到盒子的另一侧，造成了被键链接的两个原子分别位于盒子的两侧。在lammps模拟中，聚合物模拟应该算是比较复杂的一种模拟，不仅建模复杂、势参数较多，而且在模拟过程中也会经常出现各种意想不到的错误。本文主要解决聚合物建模过程中常遇到的一个问题：多成分的组合。

在lammps模拟中，所有的键默认是永久存在的，非正常情况下，不能断开，否则会产生"bond atoms missing”错误。如果要模拟拉伸过程中键的断裂，就需要加一个额外的命令：fix bond/break。这个命令可以设定每N步判断一次键长，如果键长超过设定的极限长度，则认为键被拉断。聚合物的拉伸模拟过程中，聚合物链并没有被拉断，而只是把不同的聚合物链分离。通过对比初始pe.data，可以发现，在拉伸过程中，一共拉断了6条键。f_4[1]、f_4[2]表示当前步数断裂键的数量和总的断裂数量。

本文分享lammps gcmc吸附的一个案例代码，上下壁面的通道材料为石墨烯，吸附的分子为tip4p水分子。大家好，我是小马老师。

以上命令可以单独的建出一个金刚石刀具模型，可以使用read_data命令把刀具模型合并到切削基体模型，也可以直接复制上面的代码，在基体模型的上方直接建立刀具模型。主要原理是通过三个斜面或者三个斜面加一个平面，合并为一个封闭区域，在区域内生成原子。lammps切削模拟的刀具形状有很多，如球形、锐角、钝角、三棱锥等刀具。球形、锐角、钝角等刀具建模已经在公众号发过，本文介绍三棱锥的建模。斜面和平面都是用region plane命令。本文介绍lammps三棱锥刀具建模方法。大家好，我是小马老师。

所有的reaxff模拟，如：氧化、腐蚀、热裂解、燃烧等模拟，均可采用此模拟框架，只需要更换data文件、势文件，代码中的control文件可以不写，电荷平衡文件也可以用reaxff替换。该案例主要模拟Fe(OH)3在高温下的反应过程，主要代码来自lammps自带的案例。lammps自带案例没有产物输出，故在此基础上稍加修改，增加了产物输出命令。如果需要系统学习反应势的模拟，可考虑报名培训课程。本文分享一个reaxff反应势的案例。in文件代码已经注释，以供参考。大家好，我是小马老师。

因为内容太多，查找起来也不是很方便，建议使用“搜索”功能，当需要查找某个方面的内容时，“Ctrl+F”打开搜索框，输入查找内容即可。5.24 “atoms are time integrated more than once”警告，常见三种原因及解决办法。1.31 compute temp/partial ,compute temp/com 温度修正。教程大约半年左右更新一次，本次为第3版，总页数600页左右，分上中下三册，每册200页。1.41 msi2lmp不能用，mpi不能连接主机，解决办法看这里。

和abacus一样，cp2k也是一款开源的第一性原理模拟软件，模拟的数据也可作为机器学习势的训练数据集。cp2k安装方法也很多，常见的有docker安装、代码编译安装。但是这两种方法对于新手或者小白来说，有点复杂。如果仅仅是用于机器学习势的训练，只需要cp2k计算能量和力的计算，不需要复杂的功能。所以，这里介绍一种非常简单的安装方法，只需要一个命令即可完成安装。该方法仅限于ubuntu系统，如果是其他linux系统，还是建议采用编译安装。

在一对一咨询过程中，发现有不少同学对于其它软件建好并通过read_data读入的模型，就不会分组了，也有同学知道分组方法，但是不会确定分组的region坐标范围，从而导致模拟结果出错。通过第三方软件建的模型，原子坐标不一定就是从0开始的，所以，在设置region命令的时候不能想当然的认为原子坐标就是从0开始，需要准确地确定坐标范围。各组大体比例如图所示，当然，用上面的代码分组后并不能改变原子的类型，也不能改变在ovito中显示的颜色，如需改变颜色，还需要用set命令更改原子类型。大家好，我是小马老师。

解决方法相对简单，即在弛豫过程中先取消固定代码，使用"fix all npt aniso"方式对整个体系进行充分的弛豫。然而，如果在弛豫之前就固定了底部的边界，导致底部原子无法调整原子间距，这将导致上层未固定原子无法完全调整到平衡状态，从而出现晃动现象。造成这种现象的根本原因是模型未能弛豫到平衡状态，导致原子间的能量过大，结构无法稳定。具体表现为在模拟过程中，模型无法保持稳定，可能会上下晃动或左右晃动。最近，在咨询过程中，多次遇到模型晃动的问题。一般来说，晶体模型的初始状态应该是平衡状态。

最近几个学员在做球形化合物相关的模拟，为方便大家学习，本文以Al2O3为例介绍一下球形化合物的建模方法。假如球的半径为20埃米，扩胞时要保证Al2O3在各个方向的尺寸超过40埃米，这样才能切出完整的一个球体。只要能找到初始的cif文件，都能够用这个方法建出球体、圆柱、圆锥、圆台等结构。先在长方体结构中选择一个球形的区域，然后删掉球区域外侧的原子。下载的Al2O3模型是一个斜交结构，先进行正交化。首先到晶体库网站下载相应的cif文件。来源：公号 lammps加油站。大家好，我是小马老师。

xcm（group-id，dim）函数可以计算group的质心的三个坐标，其中，dim表示质心坐标的方向（x、y、z）以上命令计算原子组Cu的质心，质心可通过c_1[*]或者c_1[1]、c_1[2]、c_1[3]输出。在润湿、蒸发、烧结或者冲击模拟中，可能需要计算某个原子组的重心，以判断其移动速度或者相对运动的情况。下面的案例分别使用以上两种方法计算Cu球的质心x坐标，并对两个结果进行对比。compute com命令同样可以计算原子组的质心，计算结果是一个矢量坐标。来源：公众号<lammps加油站>

写完in文件代码后，点击“Run”-“Run LAMMPS"，可以直接启动运算，不需要输入lmp -in in.melt.lmp这些命令，非常的方便。运行效果如下图所示，这个版本带有独立的编辑界面，可以使用菜单新建或者打开in文件，也可以使用菜单运行in文件，甚至可以直接查看运行结果。lammps一直没有编辑界面，对新手来说特别的不友好，不过，今年8月4号lammps推出了一款包含界面的版本。打开之后可以发现，lammps编辑界面非常清爽，不同的命令使用不同的颜色进行区分，类似vs code的效果。

npt iso方式下，三个方向同比例调整压力，如果三个方向尺寸不一致，也可能无法消除拉伸方向的应力。这个问题比较常见，在平时一对一指导过程中遇见多次，使用以上方法基本上都能解决初始应力不为零的问题。通常情况下，初始应力为零，以获得准确的应力-应变曲线。有的同学可能会说，明明已经进行了长时间的弛豫，可还是不能解决问题。然而，在某些情况下，模拟得到的初始应力大于零或者小于零。虽然进行了弛豫，但弛豫代码不合理的话，也会出现这种情况。如沿x方向进行拉伸，可单独对x方向进行弛豫。学lammps有问题找小马老师。

fix ave/histo命令是LAMMPS中一个非常有用的命令，可以帮助我们对系统的物理性质进行直方图统计分析。fix ave/histo命令则是LAMMPS中一个非常有用的命令，它可以帮助我们对系统进行直方图统计分析。通过该命令，我们可以获得物理量的分布情况，例如原子速度的分布、能量的分布等，进而深入了解系统的性质。通过该命令，我们可以深入了解系统中物理量的分布情况，为后续的分析和研究提供有效的数据支持。它可以将系统中某个物理量的取值范围划分成一系列的区间，并统计每个区间内的取值个数或概率。

通过调整模型的参数，可以研究不同条件下水的行为，如温度、压力和溶质浓度等，揭示其在各种环境条件下的性质变化。在模拟过程中，需要使用刚体命令来保证TIP5P水分子结构的稳定性，因此，在TIP5P分子模板中，不需要设置键角参数，仅需给出5个原子的坐标即可。TIP5P水分子是在TIP3P模型基础上添加了两个额外的虚原子L，通常这两个虚原子质量忽略不计，其中氧原子的电荷被放置在这两个虚原子L上。在lammps模拟中，常见的水分子结构为SPC/E、TIP3P和TIP4P，本文介绍另外一种水分子结构：TIP5P。

其中，data文件是lammps中最常用的输入文件之一，用于描述模拟系统的原子类型、坐标、力场参数等信息。atomsk建的模型默认的为atomic类型，如果atomic类型的data文件和full类型的data文件进行合并，需要把atomic类型转换为full类型。无论是将data文件转换为其他格式，还是使用转换后的文件进行进一步的模拟或分析，都需要谨慎处理数据，确保结果的准确性和可靠性。如果data文件不大，可以直接把data文件中的原子坐标部分复制到excel中，添加分子id和电荷信息两列即可。

minimize命令也会造成原子的移动，所以，固定的命令一般要放到minimize之前，不然，在minimize阶段也可能会丢失原子。以上是f边界条件下原子丢失的几种解决办法，没有一个方法是万能的，用哪个方法，还需要具体分析，不能盲目的套用。以切削模拟为例，边界条件设为ppf，在模拟的过程中，箭头处的原子有可能会移动到边界外面造成原子丢失。既然f边界处的原子在弛豫的过程中会移动，可增加一个真空层，给边界原子一个移动的空间。这就存在一个问题，过早的固定边界可能会影响整个体系的弛豫。

在指定的区域内，每隔一定步数将特定数量的原子或分子插入box中，直到有N个原子或分子被插入。通常用来模拟颗粒在重力作用下倒入容器的过程。下面案例是lammps自带的一个fix pour案例，主要模拟在体系内插入3000个原子，原子在重力作用下沉积到底部的过程。在绝大多数的分子动力学模拟过程中，原子数量保持恒定。为突出重点，删掉了若干不需要的语句。也可以使用fix pour命令。公众号：lammps加油站。

默认情况下，ovito会给出所有原子对的rdf，如果仅需要个别原子的rdf，可以先使用“select type”命令选择特定原子，然后再选择“coordination analysis”如下图仅计算了1-1的rdf。ovito计算rdf比较方便，但也有个缺点，只能输出某一帧的rdf值，无法对一定步数内rdf值进行平均（pro版本可以），如需平均值，可选择之前介绍的计算方法。Number of histogram bins：计算rdf时分的块数为100上图中的参数和下面的lammps命令是等效的。

斜交转正交后，晶格的划分方法变了，所以晶格的尺寸和晶格内原子个数均发生变化，在扩胞的时候，需要按正交后的晶格尺寸计算最终的盒子尺寸。虽然lammps支持斜交模型的模拟，但是在模拟过程中有一定的限制，如斜交模型的变形量不能太大，超过限定值后会提示错误。斜交改正交后，只是改变了整体盒子的形状，并没有改变原子间的相对位置，材料的性质是完全一样的。设定新的晶格，对于120度夹角的晶格，直接按下图参数更改。所以，多数情况下，把斜交的模型变成正交模型会更好一些。在ms中，斜交改正交，可使用重定义晶格的方法。

如果使用ms建模，聚合物的链会伸出到盒子外面，在导出data文件后，盒子外面的原子被周期性算法折叠进盒子对面，这就造成了边界附近的bond原子可能被分布到了盒子的两侧。并且，这个data文件不能使用npt系综进行压缩，npt系综必须使用周期性边界，在压缩时，bond原子同样会移动到对面的边界。所以，最好的解决办法是使用一个反射墙墙进行压缩，设置一个反射墙墙既可以实现压缩效果，也可以组织原子跨越box边界移动到对面。解决这个问题的原则只有一个：导出data时确保所有的bond原子没有跨越边界。

如果对编程比较熟悉，可以使用python或者MATLAB编个脚本提取，技术含量有点高，对于不熟悉编程的同学来说，可能有点难。计算表面粗糙度，需要用到表面原子的坐标，计算的一个关键技术是如果把表面的原子和其它原子区分开，或者单独把表面原子提取出来。（3）选择右上角“Add modification”-“Inver selection”，使用反向选择功能，选择非表面原子。Construct surface mesh主要功能是自动识别模型的表面，并且能标注表面上的原子。这样，表面上的原子就全部被选中。

lammps中文教程

在MS建模后，需要使用msi2lmp把car类型的文件转为lammps data文件，有时候模型在ms中显示没有问题，但是转换为data文件后结构会错位。在ovito打开转换后的data文件，如下图所示，红框内的结构应该在最下面，把鼠标放到红框内左下角的原子，ovito窗口的左下角会显示该原子的坐标。运行in文件后，会得到新的结构new.data，在ovito中打开new.data，检查结构是否符合预期，如果移动距离偏大或偏小，再进行微调。大家好，我是小马老师。（1）确定移动的距离。

如果dump文件是用dump custom 加通配符星号(*)输出的，文件名除了步数不同，其余部分是相同的，如：frame.0.dump, frame.1000.dump, frame.2000.dump。如果文件名没有规律，或者位于不同的文件，也可以新建一个文件列表file_list，把所有的文件全部输入到列表中，import_file()命令读入列表即可实现加载所有文件。在前面的文章中，介绍过python ovito可以把一个dump文件按帧拆分为多个dump文件。这种情况，可使用通配符读入所有文件。

如果模拟的体系没有键角，用这种方式是可以的，但如果模拟的是聚合物体系，用这种方法导出的data文件只有原子坐标，没有键角等信息。python ovito模块不仅可以一次性的读入dump文件，也可以分别读入包含键角信息的拓扑文件和包含模拟结果的轨迹文件。平常用的最多的方法是用ovito打开dump文件，点击“file”-“export file”的方式导出data文件。在本例中，包含键角的拓扑文件可以用初始结构文件pp.data，轨迹文件就是模拟产生的dump.xyz文件。公众号：lammps加油站。

方法很简单，使用python ovito模块，先读入dump文件，然后设置一个循环，依次读入每一帧的数据并导出到单独的文件。代码原理比较简单，import_file读入dump数据，在for循环中，export_file输出每一帧数据。在lammps模拟中，如果体系比较大，或者模拟时间比较长，得到的dump文件有可能会非常大。如果在模拟计算时，没有使用这个方法，所有的数据都输出到一个dump文件，还有一个补救的办法。使用这个方法，也可以输出某一特定步数的dump文件，比如只输出第1000步的dump文件。