【Lammps】拉伸模拟实验

一、代码

以下是使用lammps模拟对金属 Al 进行拉伸实验的 in文件。

# ------------------------ INITIALIZATION ----------------------------
units 		metal
dimension	3
boundary	p	p	p
atom_style	atomic
variable 	latparam equal 4.0320


# ----------------------- ATOM DEFINITION----------------------------
lattice		fcc ${latparam}
region whole block 0 10 0 10 0 10
create_box 1 whole
lattice fcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 region whole


# ------------------------ FORCE FIELDS------------------------------
pair_style eam/alloy
pair_coeff * * Al_EA.eam.alloy Al


# ------------------------EQUILIBRATION---------------------------------

variable tstp equal 0.001
timestep ${tstp}
variable tdamp equal "v_tstp * 100"

variable pdamp equal "v_tstp * 1000"
fix 1 all nve
fix 2 all langevin 300.0 300.0 ${pdamp} 904297
thermo 500
run 10000
unfix 1
unfix 2

fix 1 all npt temp 300 300 1 iso 0 0 1 drag 1
thermo 200

thermo_style custom step lx ly lz press pxx pyy pzz pe temp
run 10000
unfix 1




# ------------------------DEFORMATION---------------------------------
variable tmp equal "lx"
variable L0 equal ${tmp}

reset_timestep 0
fix 1 all npt temp 300 300 tdampy11
{pdamp} z 1 1 ${pdamp} drag 1
variable srate equal 1.0e10
variable srate1 equal "v_srate / 1.0e12"
fix 2 all deform 1 x erate ${srate1} units box remap x

# ------------------------OUTPUTS---------------------------------
# for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa]
# for units real, pressure is in [atm] = 101.325 [kPa] = 1.01325/10000 [GPa]

variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "-pxx/10000"
variable p2 equal "-pyy/10000"
variable p3 equal "-pzz/10000"


fix print all print 10 "strain
{p1} p2
{p3}" file stress-strain.txt screen no

dump     2 all atom 200 tensile_test_movie.lammpstrj
thermo 	100
thermo_style	custom step temp v_strain v_p1 v_p2 v_p3 ke pe press
run		20000

二、解析

（1）初始化

# ------------------------ INITIALIZATION ----------------------------
units 		metal
dimension	3
boundary	p	p	p
atom_style	atomic
variable 	latparam equal 4.0320

• units metal：设置单位系统为金属单位，其中压力单位是巴（1 bar = 100 kPa = 0.0001 GPa）。

• dimension 3：设置模拟为三维。

• boundary p p p：设置三个维度的边界条件为周期性边界。

• atom_style atomic：设置原子风格为原子，表示模拟中包含的是单个原子而不是分子。

• variable latparam equal 4.0320：定义一个变量latparam，其值为4.0320，这通常用于定义晶格参数。

注意： variable latparam equal 4.0320 表示定义一个变量latparam，他的值为4.032。

（2）建模

• lattice fcc ${latparam}：定义一个面心立方（FCC）晶格，晶格参数为之前定义的latparam。
• region whole block 0 10 0 10 0 10：定义一个长方体区域，尺寸为10x10x10。
• create_box 1 whole：在上述区域中创建一个模拟盒子。
• lattice fcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1：重新定义FCC晶格的方向，使其与坐标轴对齐。
• create_atoms 1 region whole：在上述区域中创建原子。

这一个板块，主要确定好初始晶胞，在创建一个盒子，在盒子里面建好模型。

其中，lattice fcc ${latparam}通过 lattice 关键字确定初始晶胞的晶格类型，如fcc（面心立方结构的晶胞） ，${latparam}是lammps语言中的一种变量引用，即我们需要调用之前定义的结果，需要使用 ${} 的格式。

（3）力场设置

pair_style eam/alloy：设置对势风格为EAM（嵌入原子法）合金势。
pair_coeff * * Al_EA.eam.alloy Al：为所有原子对指定EAM势的系数，这里使用的是铝合金的势能文件。

固定格式，确定设置模拟所需要的势场函数。

（4）热平衡和压力平衡

既然是进行模拟，一些特定的物理条件是十分必要的，比如时间步长、温度、压力、压强的设置。

设置时间单位

variable tstp equal 0.001：定义时间步长为0.001。
timestep ${tstp}：设置时间步长。
variable tdamp equal "v_tstp * 100"：定义热和压力弛豫时间。
variable pdamp equal "v_tstp * 1000"：定义压力弛豫时间。

1. variable tstp equal 0.001：定义了一个变量tstp（通常代表时间步长），并将其值设置为0.001。在模拟中，时间步长决定了模拟中时间的最小单位。

2. timestep ${tstp}：将模拟的时间步长设置为变量tstp的值，即0.001。这意味着在模拟中，每0.001单位时间将进行一次原子位置和速度的更新。

3. variable tdamp equal "v_tstp * 100"：定义了一个变量tdamp（通常代表热弛豫时间），并将其值设置为时间步长的100倍。热弛豫时间是温度达到平衡所需的时间尺度。

4. variable pdamp equal "v_tstp * 1000"：定义了一个变量pdamp（通常代表压力弛豫时间），并将其值设置为时间步长的1000倍。压力弛豫时间是压力达到平衡所需的时间尺度。

nve系综设置

fix 1 all nve：对所有原子应用NVE系综（恒定粒子数、体积、能量）。
fix 2 all langevin 300.0 300.0 ${pdamp} 904297：对所有原子应用Langevin恒温器和压力控制器。
thermo 500
run 10000
unfix 1
unfix 2

1. fix 1 all nve：创建了一个编号为1的fix组，对所有原子应用NVE系综。NVE系综是指恒定粒子数（N）、体积（V）和能量（E）的系综，也就是说，在这种系综下，系统的总能量是守恒的，但不控制温度或压力。

2. fix 2 all langevin 300.0 300.0 ${pdamp} 904297：创建了一个编号为2的fix组，对所有原子应用Langevin恒温器和压力控制器。Langevin恒温器是一种随机的、阻尼的恒温器，可以模拟原子在热浴中的动力学行为。这里设置了目标温度为300.0单位温度，压力控制器的压力弛豫时间为pdamp变量的值，最后一个参数904297是随机数种子，用于生成随机数。

3. thermo 500：设置每500步输出一次热力学信息，如温度、压力等。

4. run 10000：运行模拟10000步。

5. unfix 1和unfix 2：分别移除编号为1和2的fix组。

这段代码的目的是初始化一个模拟系统，设置时间步长和弛豫时间，然后对系统进行NVE系综模拟，并使用Langevin恒温器和压力控制器来控制温度和压力。在模拟过程中，每500步输出一次热力学信息，并在模拟结束后移除之前设置的fix组。

注意：需要重点理解 fix 关键字，它可以理解为进行时进行模拟条件的设置，因为模拟的条件可以是多种多样的，比如刚刚的例子中fix 1 和fix 2 分别表示对当前的系统施加了两种条件，在这两种条件的背景下，我们使用 run 10000，也就是模拟运行0.001*10000=10 ps（皮秒），执行完之后，我们设置的模拟条件为了防止与后续模拟条件的冲突，推荐使用 unfix 进行移除。

npt 系综设置

fix 1 all npt temp 300 300 1 iso 0 0 1 drag 1
thermo 200

thermo_style custom step lx ly lz press pxx pyy pzz pe temp
run 10000
unfix 1

这段代码是在运行NPT系综（恒定粒子数、压力、温度）的模拟。下面是对这段代码的详细解释，包括关键字的解释：

1. fix 1 all npt temp 300 300 1 iso 0 0 1 drag 1：

   • fix 1 all npt：创建了一个编号为1的fix组，对所有原子应用NPT系综。
   • temp 300 300 1：设置目标温度为300K，温度弛豫时间为1时间单位（由之前定义的tstp变量决定）。
   • iso 0 0 1：设置压力控制的方式为各向同性（iso），即在x、y、z三个方向上的压力相等，但这里指定只在z方向上控制压力（因为x和y方向的控制被设置为0）。
   • drag 1：设置系统的压力阻尼系数为1，这个参数影响压力弛豫的速率。

2. thermo 200：设置每200步输出一次热力学信息，如温度、压力等。

3. thermo_style custom step lx ly lz press pxx pyy pzz pe temp：

   • thermo_style custom：设置自定义的热力学输出信息。
   • step：输出当前的步数。
   • lx ly lz：输出当前模拟盒子在x、y、z方向上的长度。
   • press：输出当前的总压力。
   • pxx pyy pzz：输出当前在x、y、z方向上的压力分量。
   • pe：输出当前的势能。
   • temp：输出当前的温度。

4. run 10000：运行模拟10000步。

5. unfix 1：移除编号为1的fix组。

这段代码的目的是初始化一个模拟系统，设置为NPT系综，并在模拟过程中控制温度和压力。在模拟中，每200步输出一次自定义的热力学信息，包括模拟步数、盒子尺寸、总压力及其分量、势能和温度。模拟运行10000步后，移除之前设置的fix组。

注意：在这个板块内，需要根据自身体系进行系综的确定！！上述的参考代码中，采用的是 先nve系综下弛豫，在在npt系综下进行弛豫，这种处理方式并不都适用，需要结合自己的体系查阅参考文献确定。 **************

系综的选择推荐阅读下述文章：
https://blog.csdn.net/weixin_43848614/article/details/135510950

（5）拉伸模拟操作

variable tmp equal "lx"
variable L0 equal ${tmp}

reset_timestep 0
fix 1 all npt temp 300 300 tdampy11
{pdamp} z 1 1 ${pdamp} drag 1
variable srate equal 1.0e10
variable srate1 equal "v_srate / 1.0e12"
fix 2 all deform 1 x erate ${srate1} units box remap x

variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "-pxx/10000"
variable p2 equal "-pyy/10000"
variable p3 equal "-pzz/10000"

fix print all print 10 "strain
{p1} p2
{p3}" file stress-strain.txt screen no

dump     2 all atom 200 tensile_test_movie.lammpstrj
thermo 	100
thermo_style	custom step temp v_strain v_p1 v_p2 v_p3 ke pe press
run		20000

这段代码用于设置模拟盒子的尺寸、重置时间步长、应用NPT系综以及进行变形模拟的。下面是对这段代码的详细解释：

1. variable tmp equal "lx"：

   • 定义了一个变量tmp，并将其值设置为当前模拟盒子在x方向上的长度（lx）。lx是一个LAMMPS内置的热力学输出变量，它代表模拟盒子在x方向的尺寸。

2. variable L0 equal ${tmp}：

   • 定义了另一个变量L0，并将其值设置为变量tmp的值，即模拟盒子在x方向的初始长度。这个值可能用于后续的计算，比如应力-应变分析。

3. reset_timestep 0：

   • 重置模拟的时间步长为0。这通常在开始新的模拟运行之前进行，以确保从零开始计数。

4. fix 1 all npt temp 300 300 tdampy11：

   • 创建了一个编号为1的fix组，对所有原子应用NPT系综，目标温度设置为300K，温度弛豫时间设置为tdampy11变量的值，这个变量需要在脚本的其他地方定义。

5. {pdamp} z 1 1 ${pdamp} drag 1：

   • 这部分代码应该是fix命令的一部分，用于设置压力控制的方式。这里设置了压力弛豫时间为pdamp变量的值，只在z方向上控制压力，并且设置了阻尼系数为1。

6. variable srate equal 1.0e10：

   • 定义了一个变量srate，并将其值设置为1.0e10。这个变量可能用于定义应变率，即模拟中盒子尺寸变化的速率。

7. variable srate1 equal "v_srate / 1.0e12"：

   • 定义了一个变量srate1，并将其值设置为变量srate的值除以1.0e12。这可能是为了将应变率转换为特定的单位。

8. fix 2 all deform 1 x erate ${srate1} units box remap x：

   • 创建了一个编号为2的fix组，对所有原子应用变形。这里指定了在x方向上以应变率srate1进行变形，units box表示应变率的单位与模拟盒子的单位一致，remap x表示在每次变形后重新计算x方向的盒子尺寸。

（6）结果输出与存储

variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "-pxx/10000"
variable p2 equal "-pyy/10000"
variable p3 equal "-pzz/10000"


fix print all print 10 "strain {p1} p2 {p3}" file stress-strain.txt screen no

dump     2 all atom 200 tensile_test_movie.lammpstrj
thermo 	100
thermo_style	custom step temp v_strain v_p1 v_p2 v_p3 ke pe press
run		20000

这段代码是用于监控和记录模拟过程中的应力-应变曲线以及其他相关物理量的。下面是对这段代码的详细解释：

这段代码的目的是在一个拉伸测试模拟中，记录和输出模拟的物理量，如应变、压力、温度、动能、势能等。通过这些信息，研究人员可以分析材料在拉伸过程中的力学响应，绘制应力-应变曲线，并观察材料的变形和破坏过程。

1. variable strain equal "(lx - v_L0)/v_L0"：

   • 定义了一个变量strain（应变），其值为当前模拟盒子在x方向的长度lx与初始长度L0之差，除以初始长度L0。这是一个常用的工程应变定义。

2. variable p1 equal "-pxx/10000"：

   • 定义了一个变量p1，并将其值设置为当前模拟系统中x方向的压力pxx的负值除以10000。这个变量可能用于将压力从LAMMPS的内部单位转换为更常用的单位，如GPa。

3. variable p2 equal "-pyy/10000"：

   • 类似地，定义了一个变量p2，其值为当前模拟系统中y方向的压力pyy的负值除以10000。

4. variable p3 equal "-pzz/10000"：

   • 定义了一个变量p3，其值为当前模拟系统中z方向的压力pzz的负值除以10000。

5. fix print all print 10 "strain ${p1} p2 ${p3}" file stress-strain.txt screen no：

   • 创建了一个编号为print的fix组，用于每10步将应变和压力信息打印到文件stress-strain.txt中，并且不显示在屏幕上。打印的信息包括应变strain和三个方向上的压力p1、p2、p3。

6. dump 2 all atom 200 tensile_test_movie.lammpstrj：

   • 设置了一个dump命令，每200步将原子的信息（位置、速度等）保存到文件tensile_test_movie.lammpstrj中。这个文件可以用于后续的可视化，观察模拟过程中原子的移动。

7. thermo 100：

   • 设置每100步输出一次热力学信息。

8. thermo_style custom step temp v_strain v_p1 v_p2 v_p3 ke pe press：

   • 自定义热力学输出信息，包括模拟步数、温度、应变、三个方向上的压力、动能、势能和总压力。

9. run 20000：

   • 运行模拟20000步。

三、总结

针对lammps的拉伸模拟，只需要 更换步骤（2）建模，（3）力场设置，其它模拟环境的参数设置可以根据上述解释自己替换。总的来说，建模部分最难，如果需要通过读取 cif 格式的晶体文件。可以首先对cif文件的结构进行扩胞，在转换为lammps可以读取的data文件。

常见势函数的获取方法可以参考下述文章：
https://blog.csdn.net/weixin_43848614/article/details/135686128

将 cif 转换为 data格式，可以使用Python的ase库，可以参考下述文章：

https://blog.csdn.net/weixin_43848614/article/details/134489979

具体的lammps建模操作，可以参考下述文章，后面将有空时更新建模的总结版。
https://blog.csdn.net/weixin_43848614/article/details/133755233

参考博客：

https://blog.sciencenet.cn/blog-3437453-1300250.html