高性能计算服务
>
软件专区
>
最佳实践
>
高效使用LAMMPS进行合金材料分子动力学模拟
LAMMPS高效模拟中短程有序对中熵合金的强度和韧性的影响
如何高效使用LAMMPS进行合金材料分子动力学模拟?
LAMMPS是一款用于分子动力学模拟的开源软件，广泛应用于材料科学、物理、化学和生物领域。它能够对各种粒子体系进行大规模的并行计算，适用于研究材料的原子级结构、分子行为等。
目前，超算互联网已上线Lammps各个主要版本，支持开箱即用。
本次实操，我们将复现发表在力学顶刊《International Journal of Plasticity》上一项研究的模拟，详细介绍如何在超算互联网使用Lammps 7Feb2024版本，研究短程有序（SRO）对中熵合金的强度和韧性性能进行分子动力学模拟，揭示裂纹扩展过程中的纳米尺度变形机制。
论文链接：https://doi.org/10.1016/j.ijplas.2024.103980
一、仿真参数来源
文献2.1章节是本次实操获取对应仿真参数的重要来源。在后处理阶段使用北理工开发的mdapy模组进行后处理。
作者一共建立了9种不同的合金模型：三种随机固溶体（RSS）模型，RSS在lammps中可以通过set type/ratio命令轻松建立。与三个RSS模型对应的是6个具有SRO的模型。
并且作者还明确指出MC部分使用的是VCSGC系综。
二、计算准备
为了简化计算，分别从SRO和RSS模型中挑选出一个代表模型进行计算。这个两个模型分别为Ni33(CoCr)66和对应的SRO-650K模型。
在lammps中，为了使用VSCGC系综进行MC模拟，需要用户指定元素之间的化学势差，进行MC模拟的温度和MC步数：MD步数之间的比值。如果模型中有N种元素，则需要用户提供N-1个化学势差。此外，MC模拟需要在650K和1250K下进行，MD和MC模拟步数之间的频率为：每进行20次MD模拟，就进行一次MC模拟。整个模拟时长需要持续1000ps。
在MC系综弛豫完成之后，作者给出的拉伸参数为应变率为5 x 10^8/s,温度为1K，系综为NVT。且在拉伸之前需要在相同系综和温度下进行弛豫100ps。为了实现上述目标，撰写如下的in文件：
shell
processors 8 1 8
# 模型的基础设置
units metal
boundary p p p # 在建模阶段, 需要先设为PPP, 才能确保单晶的周期性结构不被破坏
timestep 0.001 # 单步仿真步长, 单位为ps
shell mkdir dump
# -------------建模参数-------------
variable lat equal 3.556633 # fcc单晶的晶格常数
variable latx equal v_lat*sqrt(6)/2 # 非标晶胞在x轴方向上的长度
variable laty equal v_lat*sqrt(2)/2 # 非标晶胞在y轴方向上的长度
variable latz equal v_lat*sqrt(3) # 非标晶胞在z轴方向上的长度
variable lengthx equal 650 # 文献中单晶在X轴上的长度, 以埃为单位
variable lengthy equal 25 # 文献中单晶在Y轴上的长度, 以埃为单位
variable lengthz equal 370 # 文献中单晶在Z轴上的长度, 以埃为单位
variable nx equal floor(${lengthx}/${latx}) # 用单晶在X轴的长度除以非标单晶晶胞的x长度, 并向下取整获得扩胞倍数
variable ny equal round(${lengthy}/${laty}) # 用单晶在Y轴的长度除以非标单晶晶胞的y长度, 并向下取整获得扩胞倍数
variable nz equal round(${lengthz}/${latz}) # 用单晶在Z轴的长度除以非标单晶晶胞的z长度, 并四舍五入获得扩胞倍数
variable crack_length equal 200 # 裂纹长度为200埃, 也就是20纳米
variable half_crack_width equal 5 # 裂纹宽度的一半为5埃, 也就是0.5纳米
variable fixed_layer_length equal 10 # 固定层的长度, 根据文献上确定为10埃米
# -------------MD仿真参数-------------
variable thermo_step equal 5000 # 多少次输出一次热力学信息
variable dump_step equal 2000 # 多少步输出一次轨迹文件
variable relax_time equal 100 # 弛豫多久时间, 单位为Ps
variable relax_step equal round(v_relax_time/dt) # 计算出来的需要弛豫多少步数
variable relax_temp equal 1 # 弛豫温度, 单位为K
variable strain_rate equal 5e8 # 应变率, 以s^-1为单位
variable final_strain equal 0.5 # 期望的最后模型达到的应变, 这里是20%
# -------------MC仿真参数-------------
variable mc_md_time equal 1000 # MC/MD混合模拟的持续时间, 以ps为单位
variable mc_md_steps equal round(${mc_md_time}/dt) # 总共MC/MD混合模拟的运行步数
variable mc_freq equal 20 # MC交换的频率, 即每多少步MD进行一次MC交换
variable mu_Ni_Co equal 0.021 # Ni和Co之间的化学势差, 单位为ev
variable mu_Ni_Cr equal -0.31 # Ni和Cr之间的化学势差, 单位为ev
variable kappa equal 1000 # vcsgc系综的方差, 越大模型中的浓度变化越小, 一般1000以上就可近似认为模型中的元素浓度不变
variable mc_temp equal 650 # 混合MC/MD温度, 单位为K
variable swap_fraction equal 0.25 # 随机抽取原子的比例
variable target_concentration equal "1.0/3.0" # 目标浓度, 这次做的是等原子比合金, 3种合金占比分别为1/3
# 建模设置, 单晶取向参考文献设置为X轴[11-2], Y轴[1-10], Z轴[111],
# 但是在lammps中要求三个晶向指数满足右手定则, 即X轴晶向指数叉乘Y轴
# 晶向指数结果需要等于Z轴晶向指数, 因此设置Z轴为[-1-1-1]
lattice fcc ${lat} orient x 1 1 -2 orient y 1 -1 0 orient z -1 -1 -1
region box block 0 ${latx} 0 ${laty} 0 ${latz} units box
# 模型中含有Ni, Co, Cr三种原子
create_box 3 box
# 向盒子中填充1号原子
create_atoms 1 box
# 由于lammps的非标晶向指数建模存在BUG, 因此只能通过建立单胞
# 之后, 再指定xyz方向上的扩胞次数来实现非标单晶的建模
replicate ${nx} ${ny} ${nz}
# 势函数设置
pair_style eam/alloy
pair_coeff * * NiCoCr.lammps.eam Ni Co Cr
# 随机替换原子
set type 1 type/ratio 2 0.666666 123
set type 2 type/ratio 3 0.5 123456
# --------------------------------------------------开始进行MC/MD模拟仿真--------------------------------------------------
# 创建初始温度
velocity all create ${mc_temp} 428459 dist gaussian
# 创建NPT和VCSGC对应的fix
fix integrate all npt temp ${mc_temp} ${mc_temp} $(100*dt) aniso 0.0 0.0 $(1000*dt)
fix mc all sgcmc ${mc_freq} ${swap_fraction} ${mc_temp} ${mu_Ni_Co} ${mu_Ni_Cr} &
randseed 324234 &
variance ${kappa} ${target_concentration} ${target_concentration}
# 设置对应的热力学输出
thermo ${thermo_step}
thermo_style custom step cpu temp time pe press lx ly lz
thermo_modify flush yes
# 运行指定的步数
run ${mc_md_steps}
# 取消所有的fix和dump
unfix integrate
unfix mc
reset_timestep 0
# 将模型从650K降温至1k, 在100ps以内
fix 1 all npt temp ${mc_temp} ${relax_temp} $(100*dt) aniso 0.0 0.0 $(1000*dt)
run 100000
unfix 1
reset_timestep 0
# 导出模型
write_data mc_650K_NiCoCr.lmp
# --------------------------------------------------开始进行拉伸--------------------------------------------------
# 构建边缘裂纹
region crack block $(xlo) $(xlo+v_crack_length) INF INF &
$((zlo+zhi)/2-v_half_crack_width) $((zlo+zhi)/2+v_half_crack_width) units box
delete_atoms region crack
# 将模型的边界调整到与文献中一样
change_box all boundary s p s
# 进行能量最小化
thermo 100
fix 1 all box/relax y 0.0 vmax 0.1
minimize 1e-8 1e-8 10000 10000
min_style cg
unfix 1
reset_timestep 0
# 划分固定组和变形组，固定组的原子将被施加恒定的速度以进行拉伸
region bottom block INF INF INF INF $(zlo) $(zlo+v_fixed_layer_length) units box
region top block INF INF INF INF $(zhi-v_fixed_layer_length) $(zhi) units box
group top region top
group bottom region bottom
group boundary union top bottom
group mobile subtract all boundary
# 导出模型进行查看固定组的分区是否正确
write_dump boundary atom boundary.xyz
# 为原子创建初始温度1K, 并确保整个模型线动量守恒
velocity all create ${relax_temp} 4928459 mom yes dist gaussian
# 定义热力学输出
thermo ${thermo_step}
thermo_style custom step cpu temp pxx pyy pzz press pe
# -------------开始弛豫-------------
# 在NVT系综下, 1K, 弛豫100ps,
fix 1 all nvt temp ${relax_temp} ${relax_temp} $(100*dt)
run ${relax_step}
unfix 1
reset_timestep 0
# 弛豫完成之后, 导出模型以进行后续的计算
write_data after_relaxed.lmp
# -------------开始拉伸-------------
# 计算拉伸过程会用到的中间变量
variable lzz equal lz # 由于模型在Z轴拉伸, 获取当前模型在z轴的长度
variable clzz equal ${lzz} # 将上面变量转为一个常量(constant)
variable strain_rate_ps equal v_strain_rate/1e12 # 计算在ps下的应变速率
variable velz equal v_strain_rate_ps*v_clzz # 计算应该给固定层在Z轴上的拉伸速率
variable tension_step equal round(v_final_strain/${strain_rate_ps}/dt) # 计算为了达到期望的拉伸应变需要运行多少步
variable stressz equal -pzz/10000 # 拉伸方向的应力, 单位为GPa
variable strainz equal v_strain_rate_ps*dt*step # 注意使用这个公式之前需要用reset_timestep为0
# 沿着Z轴拉伸, 所以忽略Z轴速度计算mobile组的温度
compute newtemp mobile temp/partial 1 1 0
# 固定原子, 将boundary组内的原子速度都归零
fix fixed boundary setforce 0.0 0.0 0.0
# 设置初始速度
velocity bottom set 0.0 0.0 0.0 units box
velocity top set 0.0 0.0 ${velz} units box
# 设置速度梯度
velocity mobile ramp vz 0 ${velz} z $(bound(bottom, zmax)) $(bound(top, zmin)) units box
# 定义热力学输出为修改之后的温度
thermo_modify temp newtemp
# 在NVT系综下保温拉伸
fix 1 all nvt temp ${relax_temp} ${relax_temp} $(100*dt)
# 将应力应变曲线打印出来
fix 3 all print 500 " ${strainz} ${stressz} " file unave_stress.txt screen no # 输出未平均的应力应变
fix 4 all ave/time 2000 5 10000 v_strainz v_stressz file ave_stress.txt off 1 # 不对变量1进行平均输出
# 修订控温对象为忽略Z轴速度计算模型温度
fix_modify 1 temp newtemp
# 导出轨迹文件
dump 1 all custom/gz ${dump_step} dump/tension_*.dump.gz id type mass x y z vx vy vz
run ${tension_step}接下来是对RSS模型进行拉伸，由于不需要进行MC/MD混合模拟，它的in文件会相对来说比较简洁一点：
shell
# 模型的基础设置
units metal
boundary p p p # 在建模阶段, 需要先设为PPP, 才能确保单晶的周期性结构不被破坏
timestep 0.001 # 单步仿真步长, 单位为ps
shell mkdir dump
# -------------建模参数-------------
variable lat equal 3.556633 # fcc单晶的晶格常数
variable latx equal v_lat*sqrt(6)/2 # 非标晶胞在x轴方向上的长度
variable laty equal v_lat*sqrt(2)/2 # 非标晶胞在y轴方向上的长度
variable latz equal v_lat*sqrt(3) # 非标晶胞在z轴方向上的长度
variable lengthx equal 650 # 文献中单晶在X轴上的长度, 以埃为单位
variable lengthy equal 25 # 文献中单晶在Y轴上的长度, 以埃为单位
variable lengthz equal 370 # 文献中单晶在Z轴上的长度, 以埃为单位
variable nx equal floor(${lengthx}/${latx}) # 用单晶在X轴的长度除以非标单晶晶胞的x长度, 并向下取整获得扩胞倍数
variable ny equal round(${lengthy}/${laty}) # 用单晶在Y轴的长度除以非标单晶晶胞的y长度, 并向下取整获得扩胞倍数
variable nz equal round(${lengthz}/${latz}) # 用单晶在Z轴的长度除以非标单晶晶胞的z长度, 并四舍五入获得扩胞倍数
variable crack_length equal 200 # 裂纹长度为200埃, 也就是20纳米
variable half_crack_width equal 5 # 裂纹宽度的一半为5埃, 也就是0.5纳米
variable fixed_layer_length equal 10 # 固定层的长度, 根据文献上确定为10埃米
# -------------仿真参数-------------
variable thermo_step equal 5000 # 多少次输出一次热力学信息
variable dump_step equal 2000 # 多少步输出一次轨迹文件
variable relax_time equal 100 # 弛豫多久时间, 单位为Ps
variable relax_step equal round(v_relax_time/dt) # 计算出来的需要弛豫多少步数
variable relax_temp equal 1 # 弛豫温度, 单位为K
variable strain_rate equal 5e8 # 应变率, 以s^-1为单位
variable final_strain equal 0.5 # 期望的最后模型达到的应变, 这里是50%
# 建模设置, 单晶取向参考文献设置为X轴[11-2], Y轴[1-10], Z轴[111],
# 但是在lammps中要求三个晶向指数满足右手定则, 即X轴晶向指数叉乘Y轴
# 晶向指数结果需要等于Z轴晶向指数, 因此设置Z轴为[-1-1-1]
lattice fcc ${lat} orient x 1 1 -2 orient y 1 -1 0 orient z -1 -1 -1
region box block 0 ${latx} 0 ${laty} 0 ${latz} units box
# 模型中含有Ni, Co, Cr三种原子
create_box 3 box
# 向盒子中填充1号原子
create_atoms 1 box
# 由于lammps的非标晶向指数建模存在BUG, 因此只能通过建立单胞
# 之后, 再指定xyz方向上的扩胞次数来实现非标单晶的建模
replicate ${nx} ${ny} ${nz}
# 势函数设置
pair_style eam/alloy
pair_coeff * * NiCoCr.lammps.eam Ni Co Cr
# 随机替换原子
set type 1 type/ratio 2 0.666666 123
set type 2 type/ratio 3 0.5 123456
write_data RSS.lmp
# 构建边缘裂纹
region crack block $(xlo) $(xlo+v_crack_length) INF INF &
$((zlo+zhi)/2-v_half_crack_width) $((zlo+zhi)/2+v_half_crack_width) units box
delete_atoms region crack
# 将模型的边界调整到与文献中一样
change_box all boundary s p s
# 进行能量最小化
thermo 100
fix 1 all box/relax y 0.0 vmax 0.1
minimize 1e-8 1e-8 10000 10000
min_style cg
unfix 1
reset_timestep 0
# 划分固定组和变形组，固定组的原子将被施加恒定的速度以进行拉伸
region bottom block INF INF INF INF $(zlo) $(zlo+v_fixed_layer_length) units box
region top block INF INF INF INF $(zhi-v_fixed_layer_length) $(zhi) units box
group top region top
group bottom region bottom
group boundary union top bottom
group mobile subtract all boundary
# 导出模型进行查看固定组的分区是否正确
write_dump boundary atom boundary.xyz
# 为原子创建初始温度1K, 并确保整个模型线动量守恒
velocity all create ${relax_temp} 4928459 mom yes dist gaussian
# 定义热力学输出
thermo ${thermo_step}
thermo_style custom step cpu temp pxx pyy pzz press pe
# -------------开始弛豫-------------
# 在NVT系综下, 1K, 弛豫100ps,
fix 1 all nvt temp ${relax_temp} ${relax_temp} $(100*dt)
run ${relax_step}
# run ${relax_step}
unfix 1
reset_timestep 0
# 弛豫完成之后, 导出模型以进行后续的计算
write_data after_relaxed.lmp
# -------------开始拉伸-------------
# 计算拉伸过程会用到的中间变量
variable lzz equal lz # 由于模型在Z轴拉伸, 获取当前模型在z轴的长度
variable clzz equal ${lzz} # 将上面变量转为一个常量(constant)
variable strain_rate_ps equal v_strain_rate/1e12 # 计算在ps下的应变速率
variable velz equal v_strain_rate_ps*v_clzz # 计算应该给固定层在Z轴上的拉伸速率
variable tension_step equal round(v_final_strain/${strain_rate_ps}/dt) # 计算为了达到期望的拉伸应变需要运行多少步
variable stressz equal -pzz/10000 # 拉伸方向的应力, 单位为GPa
variable strainz equal v_strain_rate_ps*dt*step # 注意使用这个公式之前需要用reset_timestep为0
# 沿着Z轴拉伸, 所以忽略Z轴速度计算mobile组的温度
compute newtemp mobile temp/partial 1 1 0
# 固定原子, 将boundary组内的原子速度都归零
fix fixed boundary setforce 0.0 0.0 0.0
# 设置初始速度
velocity bottom set 0.0 0.0 0.0 units box
velocity top set 0.0 0.0 ${velz} units box
# 设置速度梯度
velocity mobile ramp vz 0 ${velz} z $(bound(bottom, zmax)) $(bound(top, zmin)) units box
# 定义热力学输出为修改之后的温度
thermo_modify temp newtemp
# 在NVT系综下保温拉伸
fix 1 all nvt temp ${relax_temp} ${relax_temp} $(100*dt)
# 将应力应变曲线打印出来
fix 3 all print 500 " ${strainz} ${stressz} " file unave_stress.txt screen no # 输出未平均的应力应变
fix 4 all ave/time 2000 5 10000 v_strainz v_stressz file ave_stress.txt off 1 # 不对变量1进行平均输出
# 修订控温对象为忽略Z轴速度计算模型温度
fix_modify 1 temp newtemp
# 导出轨迹文件
dump 1 all custom/gz ${dump_step} dump/tension_*.dump.gz id type mass x y z vx vy vz
run ${tension_step}三、在线启动Lammps
不考虑MC/MD混合模拟，其余部分两个in文件是完全相同的。那么接下来就是提交任务了。
在超算互联网购买后，通过我的商品，在应用软件中找到或搜索该版本商品，点击“命令行”，即可快速使用Lammps软件。 
https://www.scnet.cn/ui/mall/detail/goods?type=software&common1=APP_SOFTWARE&id=1770722763948261377&resource=APP_SOFTWARE&keyword=LAMMPS2024
将写好的lammps输入文件上传到超算互联网平台并提交计算。
计算完成之后，从超算互联网平台下载轨迹文件。
在正确运行完成in文件之后，我们可以得到轨迹文件和包含应力应变数据的ave_stress.txt文件，注意在上述in文件中，输出的轨迹文件是压缩之后的，但ovito可以直接读取，可以节省大量的硬盘空间。
接下来就是运用我们得到的应力应变曲线，来获得峰值应力和平均名义应力。
首先新建一个excel文件并打开ave_stress.txt文件，将txt文件中的所有数据选中并复制，然后粘贴到刚刚新建立的excel文件中。对数据进行分列操作，注意此时的应力单位为GPa。之后选中应力那一列，在插入界面选择插入折线图，即可绘制出应力曲线。 峰值应力和平均名义应力的计算在excel可以很轻松的完成。对于峰值应力，可以调用excel中的max函数，指定数据范围为C1到C100之后，max函数就会返回指定范围内最大的数据，而这个数据就是峰值应力了。
对于平均名义应力计算来说，可以使用excel的average函数来轻松完成，指定数据范围为C26到C100。
四、表面能计算
接下来，作者为了量化不同样品的韧性值，分别计算了样品的表面能和不稳定堆垛参差能，并用它们的比值作为韧性参数。
利用lammps计算[111]晶面的表面能非常简单，in文件实现的具体的思路就是首先构建一块单晶，这个单晶有一个轴的晶向指数为[111]。首先在边界条件为PPP下计算出当前模型的势能P0，接着通过change_box命令增加盒子长度以引入真空层，使模型形成了两个[111]的自由表面。接着计算当前模型的势能P1。最后，代入公式之后就可以轻松求出表面能Es的值了。
GSFE曲线的计算思路和表面能大致相同。都可以概括为人工构建面缺陷，然后计算含有面缺陷的单晶和完美单晶之间的势能差，然后利用势能差除以面缺陷的面积即可。
只不过这一次需要的不是一个数值而是一整条曲线，只需要通过一个循环即可。
五、面缺陷识别
然后，作者使用了OVITO的CNA算法，以及基于OVITO的CNA算法开发的PDA算法展示了裂纹扩展过程中的原子微结构变化和面缺陷的变化，不过这里需要说明一下，PDA算法效率相当低下，因此这里推荐使用mdapy进行计算。
调用上述脚本之后，轨迹文件中会多出一列名为fault_types的原子属性，利用expression select选中，然后逐一assign color，即可达到论文中的可视化效果，而fault_types与面缺陷的对应关系如下。
接着，基于上述脚本加入一个for循环就可以获得面缺陷的变化曲线了。
关于位错密度的变化曲线，通过ovito左上角的导出文件，选择类型为table of value就可以导出每一帧位错线长度。这样利用excel我们就可以轻松计算出位错密度曲线。
紧接着作者统计并绘制了各个样品的裂纹表面积变化曲线，这一点利用ovito也可以轻松做到。
最后是模型的原子应变分布情况，通过联合调用atomic strain和color coding即可实现。此外，通过histgram可以统计观察原子应变的分布情况。
调用construct surface mesh，然后ovito就会统计出自由表面的总面积之和，然后与导出位错密度类似，将数据导出到txt中，利用excel或者origin之类的软件就可以绘制裂纹表面积的变化曲线了。
希望本期最佳实践为您提供一些LAMMPS的实战计算技巧。 注：本期最佳实践所提及的mp4视频均整理在飞书文档，可通过以下链接获取：https://bvjoh3z2qoz.feishu.cn/docx/X9f6dr7pKoSaKWxg2eecEwJFnXc