曙光智算应用发展部用户7072
最佳实践| 如何使用LAMMPS模拟中熵合金的强度和韧性
输入“/”快速插入内容
最佳实践| 如何使用LAMMPS模拟中熵合金的强度和韧性
最佳实践| 如何使用LAMMPS模拟中熵合金的强度和韧性
用户7072用户7072
用户4525用户45252025年11月25日修改
如何高效使用LAMMPS进行合金材料分子动力学模拟?
Lammps是一款用于分子动力学模拟的开源软件，广泛应用于材料科学、物理、化学和生物领域。它能够对各种粒子体系进行大规模的并行计算，适用于研究材料的原子级结构、分子行为等。
目前，超算互联网已上线Lammps各个主要版本，支持开箱即用。
本次实操，我们将复现发表在力学顶刊《International Journal of Plasticity》上的一项模拟，详细介绍如何在超算互联网使用Lammps 7Feb2024版本，研究短程有序（SRO）对中熵合金的强度和韧性性能的影响，揭示裂纹扩展过程中的纳米尺度变形机制。
一、仿真参数来源
文献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文件：
代码块
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_psdtstep # 注意使用这个公式之前需要用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}