Problem with bond breakage and formation during ReaxFF simulation in LAMMPS

OVITO
, , ,
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 use Ovito for visualization and followed the same procedure mentioned in the website i.e., 1) Load data file; 2) Load trajectory (dump file); 3) Again load trajectory (bond file). The bond file was obtained from the command:
fix 4 all reaxff/bonds 5 bonds.reaxff.lignin+ozone+water

Can anyone help me fix this issue and improve the visualization?

Thanks
Yuvam

1 Like
2

Hi Yuvam,

Could you please describe some more how you dump the atomic trajectories in LAMMPS? Do you include the atom IDs in the trajectory dump file(s) that you write?

And how do you import the atomic trajectories in OVITO? Are you using a second Load Trajectory modifier for this? If so, did you place it before or after the Load Trajectory modifier for loading the bond information in the pipeline?

-Alex

3

Hi Alex,
Thank you for your reply.
So I use the following commands to obtain the dump file and bond file:
dump 1 all atom 10 dump.reaxff.lignin+ozone+water
fix 4 all reaxff/bonds 5 bonds.reaxff.lignin+ozone+water

And the process by which I upload it in Ovito is:

  1. Uploading data file first
  2. Using Load Trajectory to upload the dump file.
  3. Using one more Load Trajectory to upload the bond file.

The process mentioned above is mentioned on the Ovito website.

This work for a simple system like methane in an empty box but in a system with the extended movement of atoms, the problem of bond stretching occurs.

Regards
Yuvam

4

Hi Yuvam,
It sounds correct what you are doing. Hard to tell for me in this case what is going wrong. I should take a first-hand look at it. Do you think you can share the simulation with me (LAMMPS input script + all necessary input data files)? Then I will try to recreate the situation you describe and can investigate what is happening.

5

Hi Alex,
Thank you for your help.
I am attaching the data file, input file, and FF file.

Regards
Yuvam
ffield.reax.ozone_oxi2 (9.4 KB)
in.lignin+ozone+water (1.6 KB)
data.lignin+ozone+water.charge (602.3 KB)

6

Thanks. But I am having trouble downloading the file in.lignin+ozone+water. For some reason the MatSci forum says that it is private or does not exist. Could you please upload it again - perhaps under a different filename?

7

Hi,
Sorry for the late reply.

I am uploading it again and also writing all the commands 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_oxi2 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 hope it will work.
Thank you again for your help

Regards
Yuvam
in.lignin+ozone+water (1.6 KB)

8

Yuvam,

Your LAMMPS script has the following structure (leaving out the irrelevant parts):

fix 4 all reaxff/bonds 5 bonds.reaxff.lignin+ozone+water
minimize 1.0e-5 1.0e-7 10000 100000
dump 1 all atom 10 dump.reaxff.lignin+ozone+water
run 1000

Notice how fix reaxff/bonds starts writing snapshots at intervals of 5 steps already during the initial minimization, whereas the dump command starts writing after the minimization at intervals of 10 timesteps. That means the sequence of bond snapshots and the sequence of atomic trajectory snapshots will be out of sync from the beginning and will further diverge during the course of the MD simulation.

OVITO doesn’t use the LAMMPS timestep information as a common time basis. It rather assumes that all input files loaded via the Load Trajectory modifier contain simulation snapshots taken at the exact same times. This should explain why the bonds you see in OVITO do not really fit to the atomic positions.

Please check yourself if you obtain correct results once you align the output timesteps of the fix reaxff/bonds and dump commands.

9

Hi Alex,
Thank you very much for pointing it out.
I ran the simulation again and used the same interval for both bonds and dump this time i.e., 10 for both. But the issue is still there. Bonds are getting stretched without breaking and long bonds are getting reflected too across the box.
Is there something else that is causing it?

Regards
Yuvam

10

So far I couldn’t reproduce the problem you describe. In order for me to do that, please also provide:

  • The file lmp_control referenced in your LAMMPS script,
  • The visualization setup you have created in OVITO, saved as a .ovito session state file,
  • The program version number of OVITO you are using.

Thanks!