蓝宝石（α-AL2O3）的磨削分子动力学仿真遇到模型错误以及程序运行不了的情况

ZJHHHHH

dimension  3                                      # 3维模拟
boundary   p p p                                  # 边界条件: x,y,z方向周期性
units metal                                       # 使用金属单位制
atom_style  full                                # 原子类型为

neighbor   2.5 bin                                # 邻居列表截断半径
neigh_modify  delay 5   check yes                          # 每5步更新一次邻居列表
timestep    0.001                                 # 时间步长(ps)

region     boundary1    block 0 14.274 0 90.596 0 116.91 units box      # 边界区域 - 底部固定层
region     boundary2    block 14.274 199.836 0 90.596 0 12.99 units box
region     boundary     union  2  boundary1  boundary2
region     temp_layer1  block 14.274 21.411 0 90.596 12.99 116.91 units box    # 恒温层区域
region     temp_layer2    block 21.411 199.836 0 90.596 12.99 25.98 units box
region     temp_layer     union  2  temp_layer1  temp_layer2
region     newton_layer block 21.411 199.836 0 90.596 25.98 116.91 units box  # 牛顿动力学层区域
region     piece union 3 boundary temp_layer newton_layer  # 整个工件区域

#球形金刚石磨粒区域 - 半径3.5nm，磨削深度1nm
region     tool      sphere 234.836 45.298 136.91  35 units box         # 球形磨粒区域，半径35Å(3.5nm)

#设置box区域
region     box block 0 300 0 91 0 300 units box        # 扩大BOX尺寸

#生成box，3种原子类型（Al, O, C）
create_box  3 box                                  # 创建包含3种原子类型的模拟盒子

#生成工件原子 Al2O3 (蓝宝石) - (0001)晶面取向
lattice custom 1.0 a1 4.758 0.0 0.0 a2 -2.379 4.118 0.0 a3 0.0 0.0 12.99 &
basis 0.0 0.0 0.352 &
basis 0.0 0.0 0.648 &
basis 0.0 0.0 0.852 &
basis 0.0 0.0 0.148 &
basis 0.333 0.666 0.018 &
basis 0.333 0.666 0.314 &
basis 0.333 0.666 0.518 &
basis 0.333 0.666 0.814 &
basis 0.666 0.333 0.685 &
basis 0.666 0.333 0.981 &
basis 0.666 0.333 0.185 &
basis 0.666 0.333 0.481 &
basis 0.306 0.0 0.25 &
basis 0.0 0.306 0.25 &
basis 0.694 0.694 0.25 &
basis 0.694 0.0 0.75 &
basis 0.0 0.694 0.75 &
basis 0.306 0.306 0.75 &
basis 0.639 0.666 0.916 &
basis 0.333 0.972 0.916 &
basis 0.027 0.360 0.916 &
basis 0.027 0.666 0.416 &
basis 0.333 0.360 0.416 &
basis 0.639 0.972 0.416 &
basis 0.972 0.333 0.583 &
basis 0.666 0.639 0.583 &
basis 0.360 0.027 0.583 &
basis 0.360 0.333 0.083 &
basis 0.666 0.027 0.083 &
basis 0.972 0.639 0.083 &

create_atoms  1 region piece &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 1 &
basis 5 1 &
basis 6 1 &
basis 7 1 &
basis 8 1 &
basis 9 1 &
basis 10 1 &
basis 11 1 &
basis 12 1 &
basis 13 2 &
basis 14 2 &
basis 15 2 &
basis 16 2 &
basis 17 2 &
basis 18 2 &
basis 19 2 &
basis 20 2 &
basis 21 2 &
basis 22 2 &
basis 23 2 &
basis 24 2 &
basis 25 2 &
basis 26 2 &
basis 27 2 &
basis 28 2 &
basis 29 2 &
basis 30 2 &

# 设置电荷
set type 1 charge 1.4175                           # Al电荷 +1.4175
set type 2 charge -0.945                          # O电荷 -0.945

#生成磨粒原子C(金刚石)
lattice       diamond 3.567
create_atoms  3 region tool                        # 在磨粒区域创建类型3原子(碳)

set type 3 charge 0.0                      # 新增：C 0

#设置原子质量
mass    1 26.981                                  # 类型1原子质量(铝)
mass    2 15.999                                  # 类型2原子质量(氧)
mass    3 12.01070                                # 类型3原子质量(碳-磨粒)

#原子分组
group      boundary region boundary                # 边界原子组
group      temp_layer region temp_layer            # 恒温层原子组
group      newton_layer region newton_layer        # 牛顿层原子组
group      tool region tool                        # 磨粒原子组
group      mobile union temp_layer newton_layer    # 可移动原子组

#保存模型文件
write_data    Sapphire_0001_small.data            # 输出初始结构数据文件
# 固定边界层 - z方向底部固定
fix      01 boundary setforce 0.0 0.0 0.0          # 固定边界层原子

#--------------------势函数设置------------------------------
pair_style hybrid buck/coul/long 12.0 lj/cut 12.0

kspace_style pppm 1.0e-5                           # 长程库仑处理

pair_coeff 1 1 buck/coul/long 284.80 0.068 14.07         # Al-Al

pair_coeff 1 2 buck/coul/long 6463.4 0.172 34.63         # Al-O

pair_coeff 2 2 buck/coul/long 31574.470 0.276 85.22       # O-O

pair_coeff 3 3 lj/cut 0.0 1.0                     # C-C (忽略交互，零势)

pair_coeff 1 3 lj/cut 3.0135 0.035078           # C-Al

pair_coeff 2 3 lj/cut 3.195 0.0032918      # C-O


#---------------------能量最小化--------------------------
dump                1 all atom 50 mini_sapphire_0001.xyz # 输出最小化过程轨迹
thermo              100                           # 每100步输出热力学信息
min_style           cg                            # 使用共轭梯度法最小化
minimize            1e-15 1e-15 5000 5000         # 能量最小化
undump              1                             # 停止输出轨迹
reset_timestep        0                           # 重置时间步计数器
write_restart restart0_sapphire_0001.equil          # 保存最小化后的重启文件

#--------------------初始化温度----------------------------
velocity            temp_layer create 300 8877423 # 初始化恒温层温度为300K

#--------------------恒温层驰豫-----------------------------------
fix                 1 temp_layer nvt temp 300 300 0.1  # NVT系综控制恒温层温度
thermo              100                           # 热力学输出频率
thermo_style     custom step temp etotal ke pe press pxx pyy pzz  # 输出内容
dump                1 all atom 1000 nvt1_sapphire_0001.xyz # 输出驰豫过程轨迹
run                 15000                         # 运行15000步
undump              1                             # 停止输出轨迹
unfix               1                             # 移除NVT约束
reset_timestep      0                             # 重置时间步计数器

# 输出驰豫后文件
write_restart restart1_sapphire_0001.equil          # 保存恒温层驰豫后的重启文件

#--------------------工件驰豫-----------------------------------
fix                 1 newton_layer nvt temp 0.1 0.1 0.1  # 低温驰豫牛顿层
thermo              100                           # 热力学输出频率
thermo_style     custom step temp etotal ke pe press pxx pyy pzz  # 输出内容
dump                1 all atom 1000 nvt2_sapphire_0001.xyz # 输出驰豫过程轨迹
run                 15000                         # 运行15000步
undump              1                             # 停止输出轨迹
unfix               1                             # 移除NVT约束
reset_timestep      0                             # 重置时间步计数器

# 输出驰豫后文件
write_restart restart2_sapphire_0001.equil          # 保存工件驰豫后的重启文件

#------------------磨削过程---------------------
# 设置磨粒运动：速度为10m/s (0.1 Å/ps)，沿x轴方向运动
velocity            tool set -0.1 0.0 0.0 sum yes units box  # 设置磨粒运动速度10m/s

compute        new_temp temp_layer temp/com        # 计算恒温层温度

thermo              200                           # 热力学输出频率
thermo_style     custom step temp etotal ke pe press pxx pyy pzz  # 输出内容

# 轨迹文件输出设置
dump                 1 all custom 1000 grinding_sapphire_0001.lammpstrj id type x y z vx vy vz fx fy fz

# 设置系综 - 根据新的边界条件调整
fix          10 tool nve/noforce                          # 磨粒使用NVE系综刚体
fix          20 newton_layer nve                  # 牛顿层使用NVE系综
# 对恒温层使用温度标定法进行控温
fix          30 temp_layer nve                    # 恒温层使用NVE系综
fix          31 temp_layer langevin 300.0 300.0 0.1 12345  # 温度重标定控温
fix_modify      31 temp new_temp                  # 使用计算的新温度
# 设置边界层为固定
fix          40 boundary   nve/noforce   # 边界层固定，不参与动力学

run          100000                               # 运行磨削过程100000步
write_restart  Sapphire_grinding_0001.restart       # 保存最终重启文件