LAMMPS GPU包加速下CO₂–壁面相互作用出现能量突变与bond lost错误

c2jeremy

各位老师好，我在进行 CO₂ 驱替烷烃体系的分子动力学模拟 时遇到了稳定性问题。
模拟体系包括石英孔壁、烷烃以及CO₂，使用的力场与软件环境如下：

• CO₂： EPM2
• 石英： CLAYFF
• 烷烃： GAFF2
• 软件版本： LAMMPS (22 Jul 2025)
• 加速方式： GPU包
  
  

模拟设置概述：

边界条件： p p p

pair_style： lj/cut/coul/long 10.0 14.0

kspace_style： pppm 5.0e-5

积分步长： timestep 1.0

恒温器： NVT (T = 330 K)

驱替方式： 通过对上下游活塞（wall groups）施加不同方向的恒定力实现压力差驱替（200 atm vs. 50 atm）

具体问题：

在GPU加速的计算中，模拟运行到一定时间后，体系出现能量突增（energy spike），随后发生“bond lost”错误，导致CO₂分子解体、体系“爆炸”。

通过对崩溃前的 restart 文件重算，发现：

• 崩溃前出现 CO₂ 分子速度极大化（velocity spike）；
• 增大 CO₂–壁面之间 LJ 平衡距离（通过调整 ε 和 σ 参数）未能解决问题；
• 当 关闭 GPU 加速（纯 CPU 计算）后，体系可以安全跨过该时刻；
• 同时观察到：库伦截断距离 (cutoff) 越大，体系越容易崩溃，但 GPU 利用率随 cutoff 增大而提升。
  
  

相关in文件：

1. units real
   
2. dimension 3
   
3. boundary p p p
   
4. atom_style full
   
5. bond_style harmonic
   
6. pair_style lj/cut/coul/long 10.0 14.0
   
7. angle_style harmonic
   
8. dihedral_style fourier
   
9. special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333
   
10. kspace_style pppm 5.0e-5
    
11. pair_modify mix arithmetic
    
12. 
    
13. read_data combined_system.data
    
14. include combined_system.in.settings
    
15. change_box all z final -150 700.0
    
16. 
    
17. # groups
    
18. group upper_wall id 213292:230217
    
19. group lower_wall id 230218:247143
    
20. group mobile subtract all upper_wall lower_wall freeze_group
    
21. group oil type 1 2
    
22. group CO2 type 8 9
    
23. group quartz_surface subtract mobile oil CO2
    
24. group liquid union oil CO2
    
25. 
    
26. neighbor 2.5 bin
    
27. neigh_modify delay 0 every 1 check yes one 4000
    
28. 
    
29. velocity mobile create 330.0 12345 dist gaussian
    
30. fix thermo_liquid liquid nvt temp 330.0 330.0 100.0
    
31. fix thermo_surface quartz_surface nvt temp 330.0 330.0 100.0
    
32. fix freeze_atoms freeze_group setforce 0.0 0.0 0.0
    
33. 
    
34. # piston forces
    
35. variable area equal "111.32 * 114.01"
    
36. variable P_hi equal 200.0
    
37. variable P_lo equal 50.0
    
38. variable P2F equal 1.4583e-5
    
39. variable F_hi_final equal "-v_P_hi * v_area * v_P2F"
    
40. variable F_lo_final equal "v_P_lo * v_area * v_P2F"
    
41. 
    
42. variable N_hi equal count(upper_wall)
    
43. variable N_lo equal count(lower_wall)
    
44. variable f_hi_per_atom equal v_F_hi_final/v_N_hi
    
45. variable f_lo_per_atom equal v_F_lo_final/v_N_lo
    
46. 
    
47. fix f_force_hi_ramp upper_wall addforce 0.0 0.0 v_f_hi_per_atom
    
48. fix f_force_lo_ramp lower_wall addforce 0.0 0.0 v_f_lo_per_atom
    
49. fix wall all_walls nve
    
50. 
    
51. run 2000000
    

复制代码

已尝试的解决方法

• 调整 LJ 参数（增加壁面与 CO₂ 的平衡距离）→ 仍出现爆炸；
• 减小库伦 cutoff → 稳定性提升，但 GPU 利用率下降；
• 关闭 GPU 包 → 体系稳定运行。
  
  

TROBRIZE

个人感觉跟用了GPU加速关系不大，关键还是模拟不够稳定导致的，但从描述来看也没有出现构型很不合理的（也可能是图太小了or截图不完整）。直接用刚性的CO2和CH4力场就可以了吧，不仅不会出现lost bond还不会被审稿人挑刺（我之前用非刚性的力场被审稿人质疑timestep太大），非要用bond力场的话timestep得取<1.0 fs，用run_style respa 2 2计算。

c2jeremy

TROBRIZE 发表于 2025-10-28 14:25
个人感觉跟用了GPU加速关系不大，关键还是模拟不够稳定导致的，但从描述来看也没有出现构型很不合理的（也 ...

感谢回复，这里主要烷烃不是纯甲烷，更长的烷烃将C-H键固定好像不太合适吧

TROBRIZE

c2jeremy 发表于 2025-10-28 15:46
感谢回复，这里主要烷烃不是纯甲烷，更长的烷烃将C-H键固定好像不太合适吧

抱歉没看清楚以为是甲烷……烷烃的话就只能先试试用respa跑到稳定构象了。不过我感觉CO2密度是不是太大了，如果不是做临界态的话密度改小点试试？或者活塞压力改小一点试试。