使用reaxff力场模拟出来丙烯酸酯分子密度比文献值偏高

aofchas

大家好，本人想用reaxff模拟丙烯酸脂交联的性质。但是在对丙烯酸酯单体弛豫之后获得的密度一直比论文原文偏高。我担心是我的弛豫过程有问题 各位老师能帮忙看一下吗
我用的reaxff是从以下文献来的

Development of a Reactive Force Field for Simulating Photoinitiated Acrylate Polymerization

文中对MD的过程是这么描述的

All MD simulations were performed in LAMMPS, and all initial configurations with different monomers were independently generated using Packmol [size=17.008px][url=](24)[/url] where molecules are randomly distributed, and atom overlaps were minimized. The initial box sizes were determined assuming the densities of the liquid monomers were the same as experimental values. All initial configurations were subjected to equilibration. During equilibration, a series of short trajectories with short time steps were performed and the temperature was carefully ramped up from 100 to 300 K in the NVT (canonical) ensemble followed by 500 ps in the NPT (isothermal–isobaric) ensemble at 300 K and 1 atm to eliminate any high-energy configurations and relax the systems across the periodic boundaries. The liquid monomer properties were then taken from the equilibrated systems at ambient conditions. For radical polymerization reactions, CH4 was first introduced to the monomers and once the systems were equilibrated, one of the hydrogens on CH4 was removed to generate the CH3• radical and initiate the reaction. Reactions were run in NPT for 20 ns. For all MD simulations, the temperature was maintained using a Nosé–Hoover thermostat with a damping parameter of 25 fs, and pressure was maintained using a Nosé–Hoover barostat with a damping parameter of 250 fs. Except during the temperature ramping process, time steps of 0.25 fs were used for all runs.


文中弛豫的部分（蓝色文字）我用下面这个input文件进行了复现
我也使用了文中使用的Methyl Acrylate 120个MA用packmol装进了27A的盒子里面 lammps的版本是 2Aug2023 我通过kokkos使用了GPU加速
原文中模拟出来的密度是0.99 ± 0.02 可是我模拟出来的密度却有1.02左右 各位老师能帮我看看我写的input文件有什么问题吗 模拟出来密度偏高的原因有可能是哪些


# Define Variables
variable  Initial_T equal 100.0
variable  Final_T equal 300.0
variable  Reaction_T equal 700.0
variable  T_damp equal 25.0
variable  P equal 1.0
variable  P_damp equal 250.0
variable  Reaction_P equal 3000.0
variable  Base_Name string "ma_equ"

# Initialize simulation
package kokkos neigh half comm device
newton on
units real
atom_style full
boundary p p p

# Read data file
read_data ma_box.data

# create velocity
velocity all create ${Initial_T} 4928459 dist gaussian

# Define force field
pair_style reaxff/kk NULL checkqeq yes
pair_coeff * * SI_ffield.reax.C_H_O_N_.txt C H O

# Settings for neighbor lists
neighbor 2.5 bin
neigh_modify delay 0 every 10 check no

# Calculate properties
compute msd all msd

##############################################################################################################

# NVT ensemble: Ramp up temperature from 100 K to 300 K
reset_timestep 0
timestep 0.25
fix 1 all nvt/kk temp ${Initial_T} ${Final_T} ${T_damp}
fix 2 all qeq/reaxff/kk 1 0.0 10.0 1e-6 reaxff

# Print out thermo data
thermo_style custom step cpu time temp press density etotal c_msd[4] cpuremain
thermo 1000

# Save Output and Log
variable Ramp_up_file string "${Base_Name}_Ramp_up.lammpstrj"
variable Ramp_up_log string "${Base_Name}_Ramp_up.log"
dump 1 all custom 1000 ${Ramp_up_file} id type x y z ix iy iz
dump_modify 1 sort id
log ${Ramp_up_log}

# Run
run 50000
variable Ramp_up_restart string "${Base_Name}_Ramp_up.restart"
write_restart ${Ramp_up_restart}

##############################################################################################################

# NPT ensemble: Equilibrate at 300 K and 1 atm (500 ps)
unfix 1
undump 1
fix 3 all npt/kk temp ${Final_T} ${Final_T} ${T_damp} iso ${P} ${P} ${P_damp}

# Print out thermo data
thermo_style custom step cpu time temp press density etotal c_msd[4] cpuremain
thermo 1000

# Save Output and Log
variable Equilibration_file string "${Base_Name}_Equilibration.lammpstrj"
variable Equilibration_log string "${Base_Name}_Equilibration.log"
dump 2 all custom 5000 ${Equilibration_file} id type x y z ix iy iz
dump_modify 2 sort id
log ${Equilibration_log}

# Run
run 2000000
variable Equilibration_restart string "${Base_Name}_Equilibration.restart"
write_restart ${Equilibration_restart}

wildon

0.99和1.02差距不能接受嘛

aofchas

wildon 发表于 2024-7-25 11:24
0.99和1.02差距不能接受嘛

谢谢回复
1. 因为作者没有给出input file 我对lammps也没有那么熟悉 只是会用而已 我担心有我没有考虑到的地方
2. 作者在文章中说实际密度是0.96 模拟出来是0.99 我复现出来是1.02 我复现出来的error直接从作者的结果3%到6%了 我担心这个会影响后面聚合物的性质