本文主要把之前学习石墨烯拉伸velocity拉伸做个总结。
一、模拟环境参数设置
units metal # 使用"metal"单位，这是 LAMMPS 中的一种长度单位
dimension 3 # 模拟的维度为三维空间
boundary s p p # 周期性边界条件：s 表示周期性边界条件，p 表示非周期性边界条件（periodic、permeable、periodic）
atom_style atomic # 原子类型为 atomic，表示每个原子作为一个粒子而不考虑其内部结构
atom_modify map array sort 100 2.0 # 原子修改：map 将原子映射到处理器上，array 创建数组以存储原子信息，sort 每隔 100 步排序一次原子，2.0 是排序截止半径
neighbor 3.0 bin # 设置邻居列表参数：3.0 是邻居搜索半径，bin 使用网格作为邻居列表构建方法
neigh_modify every 1 delay 0 check yes # 修改邻居列表行为：every 1 表示每个时间步更新一次，delay 0 表示立即执行更新，check yes 表示在每个时间步检查原子是否移出模拟区域
二、定义全局变量
variable temperature equal 300 #初始温度
variable relaxtemperature equal 300 #弛豫温度
variable tensiontemperature equal 300 #拉伸温度
variable tstep equal 0.001 #步长
variable thermalstep equal 100 #输出模拟结果信息
variable dumpstep equal 100 #参数输出步长
variable medoldumpstep equal 1000 #模型输出步长
variable relaxtime equal 10000 #弛豫时长
variable tensiontime equal 10000000#拉伸时长
variable pressure equal 0 #压强控制
variable strain equal 0.3 #拉伸应变
variable velocityrate equal 0.1 #velocity拉伸应变
variable fixed equal 10 #固定端三、创建/读取模型
# 边界层分组：将类型为 1 的原子分组为 graphene
group graphene type 1
# 使用变量定义边界层的最大和最小 x 坐标
variable xmax equal bound(graphene,xmax)-${fixed} # 边界层最大 x 坐标减去 ${fixed}
variable xmin equal bound(graphene,xmin)+${fixed} # 边界层最小 x 坐标加上 ${fixed}
# 确定边界层分界点坐标
region left block INF ${xmin} INF INF INF INF units box # 定义左边界层的区域
group left region left # 将左边界层的区域分组
region right block ${xmax} INF INF INF INF INF units box # 定义右边界层的区域
group right region right # 将右边界层的区域分组
# 边界层、变形层分组：将边界层定义为左右两个区域的并集
group boundary union left right # 将左右两个边界层的区域并集作为边界层分组
# 将边界层之外的所有原子作为变形层分组
group mobile subtract all boundary
# 设置边界层原子类型：将左边界层的原子类型设置为 2，将右边界层的原子类型设置为 3
set group left type 2
set group right type 3
# 设置所有原子的质量为 12.011150
mass * 12.011150
四、势函数
pair_style airebo 3.0
pair_coeff * * CH.airebo C C C五、能量最小化、温度初始化和步长
min_style cg #也可以选取sd
minimize 1.0e-10 1.0e-10 10000 10000
velocity all create ${temperature} 534567 dist gaussian units box
timestep ${tstep}六、固定两端
#固定边界层
velocity left set 0.0 0.0 0.0 units box
velocity right set 0.0 0.0 0.0 units box
fix left left setforce 0 0 0
fix right right setforce 0 0 0
七、npt弛豫
compute 2 all pe/atom #计算每个原子的势能
compute pe mobile reduce sum c_2 #对变形区域原子势能进行求和
variable PE equal "c_pe"
thermo ${thermalstep}
thermo_style custom step temp press lx ly lz pe ke etotal
fix 1 all npt temp ${relaxtemperature} ${relaxtemperature} 0.1 y ${pressure} ${pressure} 10 z ${pressure} ${pressure} 10 drag 1
fix pe all print ${dumpstep} "${PE} " file PE2.txt screen no
dump 1 all custom ${medoldumpstep} relax_x2.lammpstrj id type x y z
run ${relaxtime}
undump 1 #取消fix、dump设定，步数清零
unfix 1
reset_timestep 0八、输出参数计算
variable num equal count(mobile) #计算变形区域原子数目
variable toval equal (lx-2*${fixed})*ly*3.35 #计算系统总体积，单层石墨烯厚度为0.335nm
variable vol equal ${toval} #体系的初始体积
variable vatom equal v_vol/v_num #单个原子体积
variable l_x equal lx-2*${fixed}
variable lx0 equal ${l_x}
variable l_y equal ly
variable ly0 equal ${l_y}
variable fixed2 equal 2*${fixed}
variable xstrain equal "(lx - v_fixed2 - v_lx0) / v_lx0" #计算x方向应变
variable ystrain equal "(ly - v_ly0) / v_ly0" #计算y方向应变
compute stress mobile stress/atom NULL #计算体系中单原子应力
compute xalls mobile reduce sum c_stress[1] #计算体系中各方向应力
compute yalls mobile reduce sum c_stress[2]
compute zalls mobile reduce sum c_stress[3]
compute xyalls mobile reduce sum c_stress[4]
variable xstress equal "1.0e-4 * c_xalls/v_toval"
variable ystress equal "1.0e-4 * c_yalls/v_toval"
variable zstress equal "1.0e-4 * c_zalls/v_toval"
variable xystress equal "1.0e-4 * c_xyalls/v_toval"九、velocity拉伸
#拉伸结果输出
thermo ${thermalstep}
thermo_style custom step press temp lx ly lz vol v_xstrain v_ystrain v_xstress v_ystress pe ke etotal
#拉伸时系综设置
fix 2 all nvt temp ${tensiontemperature} ${tensiontemperature} 0.01
#设置两端边界层双向移动
velocity right set ${velocityrate} 0 0
velocity left set -${velocityrate} 0 0
#设置中间变形区域速度梯度
velocity mobile ramp vx -${velocityrate} ${velocityrate} x ${xmin} ${xmax} sum yes units box
#启动运算
fix 4 all ave/time 1 ${dumpstep} ${dumpstep} v_xstrain v_ystrain v_xstress v_ystress file x.txt
dump 2 all custom ${medoldumpstep} load_x.lammpstrj id type x y z c_stress[1] c_stress[2] c_stress[3] c_stress[4]
fix 5 all halt ${dumpstep} v_xstrain > ${strain} error continue #拉
run ${tensiontime}
十、完整代码及结果
1、代码
units metal #单位为lammps 中的metel 类型
dimension 3 #模拟的维度 三维
boundary s p p #周期边界条件
atom_style atomic #原子类型自动
atom_modify map array sort 100 2.0
neighbor 3.0 bin
neigh_modify every 1 delay 0 check yes
variable temperature equal 300 #初始温度
variable relaxtemperature equal 300 #弛豫温度
variable tensiontemperature equal 300 #拉伸温度
variable tstep equal 0.001 #步长
variable thermalstep equal 100 #输出模拟结果信息
variable dumpstep equal 100 #参数输出步长
variable medoldumpstep equal 1000 #模型输出步长
variable relaxtime equal 10000 #弛豫时长
variable tensiontime equal 10000000#拉伸时长
variable pressure equal 0 #压强控制
variable strain equal 0.3 #拉伸应变
variable velocityrate equal 0.1 #velocity拉伸应变
variable fixed equal 10 #固定端
read_data gp.lammpstrj extra/atom/types 2
#边界层分组
group graphene type 1
variable xmax equal bound(graphene,xmax)-${fixed}
variable xmin equal bound(graphene,xmin)+${fixed}
#确定边界层分界点坐标
region left block INF ${xmin} INF INF INF INF units box
group left region left
region right block ${xmax} INF INF INF INF INF units box
group right region right
#边界层、变形层分组
group boundary union left right
group mobile subtract all boundary
#设置边界层原子类型
set group left type 2
set group right type 3
mass * 12.011150
pair_style airebo 3.0
pair_coeff * * CH.airebo C C C
min_style cg #也可以选取sd
minimize 1.0e-10 1.0e-10 10000 10000
velocity all create ${temperature} 534567 dist gaussian units box
timestep ${tstep}
#固定边界层
velocity left set 0.0 0.0 0.0 units box
velocity right set 0.0 0.0 0.0 units box
fix left left setforce 0 0 0
fix right right setforce 0 0 0
compute 2 all pe/atom #计算每个原子的势能
compute pe mobile reduce sum c_2 #对变形区域原子势能进行求和
variable PE equal "c_pe"
thermo ${thermalstep}
thermo_style custom step temp press lx ly lz pe ke etotal
fix 1 all npt temp ${relaxtemperature} ${relaxtemperature} 0.1 y ${pressure} ${pressure} 10 z ${pressure} ${pressure} 10 drag 1
fix pe all print ${dumpstep} "${PE} " file PE2.txt screen no
dump 1 all custom ${medoldumpstep} relax_x2.lammpstrj id type x y z
run ${relaxtime}
undump 1 #取消fix、dump设定，步数清零
unfix 1
reset_timestep 0
variable num equal count(mobile) #计算变形区域原子数目
variable toval equal (lx-2*${fixed})*ly*3.35 #计算系统总体积，单层石墨烯厚度为0.335nm
variable vol equal ${toval} #体系的初始体积
variable vatom equal v_vol/v_num #单个原子体积
variable l_x equal lx-2*${fixed}
variable lx0 equal ${l_x}
variable l_y equal ly
variable ly0 equal ${l_y}
variable fixed2 equal 2*${fixed}
variable xstrain equal "(lx - v_fixed2 - v_lx0) / v_lx0" #计算x方向应变
variable ystrain equal "(ly - v_ly0) / v_ly0" #计算y方向应变
compute stress mobile stress/atom NULL #计算体系中单原子应力
compute xalls mobile reduce sum c_stress[1] #计算体系中各方向应力
compute yalls mobile reduce sum c_stress[2]
compute zalls mobile reduce sum c_stress[3]
compute xyalls mobile reduce sum c_stress[4]
variable xstress equal "1.0e-4 * c_xalls/v_toval"
variable ystress equal "1.0e-4 * c_yalls/v_toval"
variable zstress equal "1.0e-4 * c_zalls/v_toval"
variable xystress equal "1.0e-4 * c_xyalls/v_toval"
#拉伸结果输出
thermo ${thermalstep}
thermo_style custom step press temp lx ly lz vol v_xstrain v_ystrain v_xstress v_ystress pe ke etotal
#拉伸时系综设置
fix 2 all nvt temp ${tensiontemperature} ${tensiontemperature} 0.01
#设置两端边界层双向移动
velocity right set ${velocityrate} 0 0
velocity left set -${velocityrate} 0 0
#设置中间变形区域速度梯度
velocity mobile ramp vx -${velocityrate} ${velocityrate} x ${xmin} ${xmax} sum yes units box
#启动运算
fix 4 all ave/time 1 ${dumpstep} ${dumpstep} v_xstrain v_ystrain v_xstress v_ystress file x.txt
dump 2 all custom ${medoldumpstep} load_x.lammpstrj id type x y z c_stress[1] c_stress[2] c_stress[3] c_stress[4]
fix 5 all halt ${dumpstep} v_xstrain > ${strain} error continue #拉伸终止
run ${tensiontime}二、结果
若有问题欢迎讨论，会持续更新改正
![[嵌入式系统-52]：RT-Thread：RT-Thread Studio工程的文件结构](https://img-blog.csdnimg.cn/direct/8722ed3a7bf14b299065ec06936d930a.png)
![[华为OD] 给航天器一侧加装长方形或正方形的太阳能板 100](https://img-blog.csdnimg.cn/direct/51754d605de44c40aec2e6d0fecada95.png)
