Discrepancy in Pressure, Coulomb Energy and Charge Equilibration Energy using ReaxFF with restart files


I am using LAMMPS (LAMMPS (28 Jun 2014)) to simulate RDX using the ReaxFF potential described in Lianchi Liu, J. Phys. Chem. A 2011, 115, 11016–11022. In order to run longer simulations, I was checking the write restart functionality and the 2 input files used for the same are attached along with the initial RDX atom definition file used for 1st simulation.

The first input file runs for 1000 timesteps and writes a restart file. The second one reads the restart file and runs for another 1000 timesteps. I find that the thermodynamic output for the last timestep from the first run and first timestep from the second run have discrepancies for all components of pressure, Coulomb energy and Charge equilibration energy computed from the pair style reax/c.

From the read_restart manual page, I find that there should not be immediate discrepancy in the restarted simulation. Below are the values of pressure components and energy components with discrepancy.

Can someone explain if the discrepancy is expected and fine or
is there some issue with my input files or
is this some issue with read_restart and ReaxFF potential ?

1st run, 2nd run,
1000th timestep 1st timestep Difference
Pxx 28351.65 28351.616 0.034
Pyy -47134.505 -47134.522 0.017
Pzz -58322.263 -58322.31 0.047
Pxy -683.09681 -683.09186 0.00495
Pyz 1352.3359 1352.3431 0.0072
Pxz 425.85928 425.85878 0.0005
Press -25701.706 -25701.739 0.033
ep -24422.635 -24422.661 0.026
eqeq 11763.065 11763.091 0.026

The units are real units. ep is the Coulomb energy and eqeq is the Charge equilibration energy.

RDX_atom_definition.txt (5.68 KB)

in1.lmp (1.66 KB)

in2.lmp (1.63 KB)

Ray can likely answer this.



This is related to QEq (fix qeq/reax) whose arrays are instantiated to zeros at the beginning of each run. Therefore you get slight differences in solved charges. If you turn of QEq in your second run, you will find exactly same results.