求助：聚合物熔体NPT收缩出现空腔

dayu8278

各位老师、同学好：

我在使用 LAMMPS 模拟线性高分子熔体的组装过程时遇到了瓶颈。目前在 NPT 平衡阶段，体系内部出现了明显的空腔（物理意义上的真空泡），这与粗粒化模拟的经验不符。

具体情况如下：

• 体系构建：
  
  
  • 模型： 100 条线性高分子链，采用 OPLS-AA 全原子力场。高分子结构如下所示
  • 初始状态： 随机置于 200A˚ 的立方盒子内，使用 pair_style soft 缓慢推开重叠原子，避免能量发散。通过soft后，100分子链均匀分布在盒子中。
    
    
    

  

• 模拟流程：
  
  
  • NVT 预松弛： 在 800K 下运行 500,000 steps。
  • NPT 压缩平衡： 设置目标压力为 1atm（或常压），运行 6,000,000 steps。
    
    
• 存在问题：随着 NPT 运行，体系虽然在收缩，但内部出现了明显的空腔（如下图所示）。
  
  

• 疑惑点：此前在粗粒化（CG）模拟中从未遇到此类情况。全原子体系下出现空腔，是由于 **NPT 压力控制参数（P-damp）不当、初始密度过低导致的亚稳态，还是力场截断半径（cutoff）**或长程相互作用设置的问题？
  
  

恳请有经验的学长或老师指点迷津，非常感谢！

模拟脚本如下：
##################################################################
# ----------- VARIABLES -----------------
variable     T_START equal 500
variable     T_HIGH  equal 800
variable     P       equal 1.0

units        real
boundary     p p p

# ----------- FORCE FIELD SETUP ----------------
atom_style   full

bond_style   harmonic
angle_style  harmonic
dihedral_style opls
improper_style cvff

pair_style   lj/cut/coul/long 12.0
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333

kspace_style pppm 1.0e-4
pair_modify  mix arithmetic

# ----------- READ SYSTEM ----------------
read_data       ../../02_Soft/lammps_C1-OPLSL/pushed_off_100chains.data nocoeff
include         ../../Init_Conf/lammps_C1-OPLSL/forcefield.params

# ----------- NEIGHBOR -------------------
neighbor     2.0 bin
neigh_modify delay 0 every 1 check yes

# ----------- 1. 能量最小化 ----------------
thermo       1000
thermo_style custom step temp press density etotal ke pe evdwl ecoul elong emol

print        "--- Phase 1: Neutral Minimization ---"
min_style    cg
minimize     1.0e-6 1.0e-8 5000 50000

# ----------- 2. 速度初始化 & 热身 ----------------
print        "--- Phase 2: Dynamics Warmup ---"
reset_timestep 0

# 【安全建议】高温退火必须使用 0.5 fs 步长，防止键断裂
timestep     0.5

velocity     all create ${T_START} 4928459 dist gaussian mom yes rot yes

# NVT 预热 (让链动起来)
fix          nvt_heat all nvt temp ${T_START} ${T_START} 100.0 drag 0.2
run          500000  
unfix        nvt_heat

# ----------- 3. 退火流程 (Annealing Process) ----------------
print        "--- Phase 3: Annealing (Eliminating Voids) ---"
dump         traj all custom 20000 anneal_neutral.lammpstrj id mol type x y z
dump_modify  traj sort id

# [3.1] 升温阶段：500K -> 800K
# 给予分子链巨大的动能来冲破缠结
print        "--- 3.1 Heating: 500K to 800K ---"
fix          npt_heat all npt temp ${T_START} ${T_HIGH} 100 iso ${P} ${P} 1000 drag 1.0
run          1000000  # 500 ps 升温
unfix        npt_heat

# [3.2] 高温保持：在 800K 持续振荡
# 这是消除空洞、形成 Melt 的核心阶段
print        "--- 3.2 Holding: 800K Melt Relaxation ---"
fix          npt_hold all npt temp ${T_HIGH} ${T_HIGH} 100 iso ${P} ${P} 1000 drag 1.0
run          2000000  # 1 ns 保持
unfix        npt_hold

# [3.3] 降温阶段：800K -> 500K
# 缓慢收缩盒子，让密度平稳回到正常值
print        "--- 3.3 Cooling: 800K to 500K ---"
fix          npt_cool all npt temp ${T_HIGH} ${T_START} 100 iso ${P} ${P} 1000 drag 1.0
run          2000000  # 1 ns 降温
unfix        npt_cool

# ----------- 4. 最终常温弛豫 ----------------
print        "--- Phase 4: Final 500K Equilibration ---"
fix          npt_final all npt temp ${T_START} ${T_START} 100 iso ${P} ${P} 1000 drag 1.0
run          1000000  # 500 ps 稳定
unfix        npt_final

# ----------- FINAL OUTPUT ----------------
undump       traj
write_data   Annealed_Melt.data

print        "ANNEALING DONE"

yuzc

松弛问题，构象熵太大导致的，而粗粒化自由能面都磨平了完全没有这个问题。
退火可以不作单调处理，具体模拟方法建议参考这篇聚合物模拟领域的人都会看的文献：
1.Larsen, G. S., Lin, P., Hart, K. E. & Colina, C. M. Molecular simulations of PIM-1-like polymers of intrinsic microporosity. Macromolecules 44, 6944–6951 (2011)
2. Polymatic: a generalized simulated polymerization algorithm for amorphous polymers | Theoretical Chemistry Accounts | Springer Nature Link
可以说是一种protocol了。

dayu8278

yuzc 发表于 2026-2-27 14:27
松弛问题，构象熵太大导致的，而粗粒化自由能面都磨平了完全没有这个问题。
退火可以不作单调处理，具体模 ...

万分感谢，已经测试了21步松弛法，确实能有效消除熔体空腔现象。

弱弱的问一句，模拟聚合物全原子尺度动力学，采用gaff2加AM1-Bcc电荷可以吗？或者OPLS-AA加AM1-Bcc电荷，请问您觉得哪个比较好，我这聚合物包含比较多的氢键位点和芳香环

yuzc

dayu8278 发表于 2026-2-28 10:24
万分感谢，已经测试了21步松弛法，确实能有效消除熔体空腔现象。

弱弱的问一句，模拟聚合物全原子尺度 ...

有些御用搭配是不建议变的。
GAFF/GAFF2搭配RESP是amber开发者推荐的，而AM1-BCC是为了便宜地得到与RESP相近电荷数值的方案。所以AM1-BCC也是建议搭配GAFF的。
OPLS-AA系列的开发者的御用电荷是1.14*CM1A和1.2*CM5，这个力场的的各种开发方案也是围绕着CM系列电荷方案的。尽管有些文献也使用了RESP和AM1-BCC搭配OPLS-AA，但是我认为作通用型处理这是不妥的，case-by-case的话可以测试一下结构和动力学性质是否能够很好重现。

yuzc

dayu8278 发表于 2026-2-28 10:24
万分感谢，已经测试了21步松弛法，确实能有效消除熔体空腔现象。

弱弱的问一句，模拟聚合物全原子尺度 ...

如果是lammps用户可以尝试使用github上的一个项目：Radonpy 自动工作流产生线形聚合物体相性质的，这个软件对聚合物进行松弛时就内置了EQ21（21步循环退火）算法。
内置了gaff2+RESP自动计算方案。
同时由于其是通过if判据进行力场分配的，同样的思路可以利用chatgpt对这个项目进行二次开发，从moltemplate里把opls-AA的smarts规则转换为Radonpy风格的力场分配判据。 这样可以通过同一个架构实现OPLS-AA力场的分配和计算，电荷同理。