[lammps-users] Implementation of ReaxFF on water

Dear Lammps users,

As a new user, I had attempted to implement the reaxFF potential on a simple system of water molecules before proceeding to a complex system. I have reviewed the reax example files provided with Lammps. For the simulation, a datafile of 60 water molecules, packed in a cube, was generated using packmol and the VMD command ‘writelammpsdata’. However, on implementing any ‘fix’ commands associated with reaxFF, the code-dumped error is returned. Any help and pointers on correcting this are greatly appreciated. The simplified input file, data and charge equilibrium parameter files defined are attached below.




  1. Input File:

---------- Initialize Simulation -----------

units real
dimension 3
boundary p p p
atom_style charge

---------- System ----------

read_data waterbox60_charge.dat

---------- Interatomic_Potential ----------

pair_style reax/c lmp_control
pair_coeff * * ffield.reax H O

compute reax all pair reax/c

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 bin
neigh_modify every 10 delay 0 check no

reset_timestep 0
timestep 0.2

fix 1 all qeq/reax 1 0.0 10.0 1e-6 param.qeq
velocity all create 300.0 4712398

thermo 1
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

fix 2 all nve
run 10
unfix 2

---------- End of Simulation -----------

  1. Param.qeq:

1 5.8678 14.0000 0.9000
2 8.5000 17.9978 1.0500

  1. Data File:

LAMMPS data file. CGCMM style. atom_style charge generated by VMD/TopoTools v1.7 on Tue May 11 15:46:45 IST 2021
180 atoms
0 bonds
0 angles
0 dihedrals
0 impropers
2 atom types
0 bond types
0 angle types
0 dihedral types
0 improper types
-0.500000 0.500000 xlo xhi
-0.500000 0.500000 ylo yhi
-0.500000 0.500000 zlo zhi

Pair Coeffs

there are several issues.

  • your data file has invalid box dimensions. the telltale sign is the warning “WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (src/domain.cpp:964)”
    your box is 1x1x1 angstrom, clearly not correct for 60 water molecules.
  • usually you won’t need a control file and qeq parameters but go with the defaults and the settings taken from the force field file
  • you have to be very careful in picking the suitable reaxff parameterization. e.g. using a combustion force field on bulk water is a mistake and will give bogus results.