Hi Ouyang,
Just wanted to quickly write up what I've found (some of this is from
memory, some I did last week):
The code does not show large discrepancies when comparing various
implementations on the LAMMPS side with your system --
however the forces deviate from what you have in your code (possibly
due to https://github.com/lammps/lammps/issues/726).
I've been experimenting with various possible fixes in areas where I
had a suspicion that a problem might exist.
I did some relatively long runs (4ns) for each of the possible
implementations, with timesteps varying from 1fs to 1/16 fs.
None of these "fixes" seemed to matter. See the table below for raw data:
Timestep | Min. energy spread | Max. energy spread (across all the
implementation variants I tested)
timestep (fs) | min | max
1 | crash | crash
0.5 | 1.2685 | 2.1206
0.25 | 0.0863 | 0.1744
0.125 | 0.015 | 0.0442
0.0625 | 0.0053 | 0.0088
Given the relative difference when halving is suspiciously close to
four (which is the expected value for verlet integration),
I'm inclined to say that all you're seeing is integration error,
amplified in the 1fs (and 0.5fs) case by the fact that it doesn't
resolve all oscillations.
This doesn't mean that the REBO code in LAMMPS is perfect (on the
contrary, see the issue referenced above), but that it, as far as I
can tell, is decently self-consistent.
However, admittedly, I am not an expert in running simulations, so I'm
open if there's a suggestion for a better test.
Cheers, Markus