Hi,

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.5*dt)
(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.

Regards,

Q