[lammps-users] questions about reax force field and its use

Hi,
I have used the scripts as showed below to take advantage of Reax force field to run molecular dynamic simulations, but the run end up with an error like this "ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after 709517 of 1634650 pairs completed.STOP Too many bonds on atom. Increase MBONDDEF ". I have no idea about how this comes up and how to fix it. Can anybody give me some advice? Shoud it be neccessary for me to combine commands like delete_bonds or some other commands related to my script? Or could there be something else I should pay attention to?Any sugesstion will be appreciated! Thank you in advance!
The scipt I have used is as follows:

read_restart sio.2975.min

pair_style reax 10.0 1.0e-6
pair_coeff * * ffield.reax 2 1 1 3 2 1 4 3

neighbor 2.0 multi
neigh_modify every 5 delay 0 check yes

restart 10 sio.*.min

thermo 1
thermo_style custom step etotal ke temp pe ebond eangle edihed eimp evdwl ecoul elong press enthalpy etail vol
thermo_modify flush yes

timestep 0.25

dump 1 all xtc 10 min.xtc
dump 2 all atom 10 min.reax..atom
dump 3 all xyz 10 min.reax.
.xyz
dump 4 all dcd 10 min.reax.dcd
dump_modify 2 scale no
dump_modify 4 unwrap yes
dump_modify 1 precision 10000

velocity all create 298 4928459 dist gaussian
fix 1 all reax/bonds 10 minbonds.tatb
fix 2 all nvt 298 298 10 chain no

run 1000
run 1000

Best Wishes!

Shuping Huang

I'll let Aidan answer this. It's been discussed on the mail list
before. It's an internal ReaxFF library error that has to
do with the memory it allocates.

Steve

I'll let Aidan answer this. It's been discussed on the mail list
before. It's an internal ReaxFF library error that has to
do with the memory it allocates.

Steve

This error occurred because at that point in the simulation, one of the
atoms had 21 bonds, whereas memory was only allocated 20 bonds. You can
increase the memory allocation by increasing MBONDDEF (as the message says).
See lib/reax/README for more info.

Thanks, Steve and Adian. But how could this happen? Why didn’t any bonds on that specific atom break? Although this is a stage for the system to reach equilibrium, it can never be physics that an atom could be connected with such a number of bonds. What is the point of such an operation? Could you please give me some spark to get through this?
Thanks again.

–Shuping

Sometimes a large value of mbond is due to an unphysical geometry. Other
times, it is just a spontaneous occurrence of a very low probability
geometry that will quickly dissappear. Finally, it could be how ReaxFF
handles certain chemical environments, for example, some metals (where
second and third shell neighbors are treated as bonded). How do you tell the
difference?

"ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after

709517

of 1634650 pairs completed.STOP Too many bonds on atom. Increase
MBONDDEF ".

The anser is hinted at by the fact that you processed more than half the
bonds before encountering a problem. This suggest it was just a very low
probability geometry. Increase MBONDDEF to 30 and you should be okay.

Aidan

Thanks, Aidan. I’ll have a try again, Thank you very much!

Cheers,
Shuping