I am trying to replicate LAMMPS’ result of MD simulation using Lennard-jones potential and verlet intergration. I did refer to source code of verlet.cpp in LAMMPS github and follow it to write my verlet as 5 steps:
(1): update velocities using half timestep: v += a * (0.5dt)
(2): update positions using full timestep: x += v * dt
(3): compute force, set it to accelerations as I am using lj unit (m=1)
(4): update velocities uisng half timestep and computed a: v+= a(0.5*dt)
To eliminate the fact that errors may belong to my neighbor-list implementation, I am not using neighbor list and just using naive implementation for force-pair calculations, i.e. traverse all pairs of atoms using nested for loop (i=0->natoms , j=i+1->natoms).
To get the same initial configs, I dumped a file containing postions, velocites, and accelerations from LAMMPS (step 0) and used them to initialize my data structures, respectively. Box length is also matched.
I am seeing some numerical miss match between my simulation and LAMMPS results, say atoms accelerations after just 1 step.
Am I missing something related to verlet integrator here? Any idea is very appreciated. Thank you.