Problem with bond breakage and formation during ReaxFF simulation in LAMMPS

LAMMPS LAMMPS Beginners
,
1

Hi everyone,

I am trying to run a ReaxFF simulation for lignin and ozone in the water system. I am using ReaxFF potential from here, which I believe is good for me. I noticed that few bonds are stretched to unrealistic lengths without breaking, and a couple of atoms get very close without forming a bond. I am not able to understand why bond formation and breakage are not taking place. Moreover, during simulation, bonds from water molecules are running across the box, which I suppose is because of PBC. I tried using bond/break but it gave me an error that I need to define bond_style which I have not. I am attaching my input file below:

units real
boundary p p p

atom_style charge
read_data data.lignin+ozone+water.charge

pair_style reaxff lmp_control
pair_coeff * * ffield.reax.ozone_oxi C H O

compute reax all pair reaxff

variable eb equal c_reax[1]
variable ea equal c_reax[2]
variable elp equal c_reax[3]
variable emol equal c_reax[4]
variable ev equal c_reax[5]
variable epen equal c_reax[6]
variable ecoa equal c_reax[7]
variable ehb equal c_reax[8]
variable et equal c_reax[9]
variable eco equal c_reax[10]
variable ew equal c_reax[11]
variable ep equal c_reax[12]
variable efi equal c_reax[13]
variable eqeq equal c_reax[14]

neighbor 2.5 bin
neigh_modify every 10 delay 0 check yes

velocity all create 300.0 4928459 mom yes rot yes dist gaussian

fix 1 all nve
fix 2 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff
fix 4 all reaxff/bonds 5 bonds.reaxff.lignin+ozone+water
#fix 5 all bond/break 5 1 1.9

variable nqeq equal f_2

thermo 10
minimize 1.0e-5 1.0e-7 10000 100000
thermo_style custom step temp epair etotal press &
v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa &
v_ehb v_et v_eco v_ew v_ep v_efi v_eqeq v_nqeq

timestep 0.1

dump 1 all atom 10 dump.reaxff.lignin+ozone+water

run 1000

I would be grateful for your expertise.

Best Regards
Yuvam

2

This all sounds not like a problem of your simulation, but of your visualization. How are the bonds you “see” determined? Many visualization tools for molecular simulations determine them heuristically and won’t update them, because for typical simulations, there are no reactions and the bond topology is an input parameter of the simulation.

3

Hi,
Thank you for your reply.
I tried two different visualization software - VMD and Ovito and both visualized bonds in the same way.
For VMD I just used the dump file. Whereas for Ovito, I used a pipeline with the data file, the dump file, and the bond file (extracted from the simulation at every 5 steps). Can you suggest a better visualization software for reaxff simulation?
Thank you

4

No. But there may be ways to use the existing tools.
For VMD you could use the DynamicBonds representation. This won’t be overly accurate, but better than not updating bonds at all. Please check out the User’s Guide.
For Ovito you have to ask the Ovito developers. There is are Ovito categories on this site.