Replicate a LAMMPS simulation


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
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.


Please note that your question is not really about LAMMPS but about your own MD code and thus it is off-topic for the LAMMPS category.

Studying the LAMMPS code to learn basic MD timestepping is in my personal opinion not a good idea. Because LAMMPS has so many features it is quite an effort to put all the bits and pieces together (e.g. you need not only look at verlet.cpp but also fix_nve.cpp, not to mention the neighbor list and domain decomposition and parallel communication code.
There are loads of MD program tutorials on the web. I have written a minimal (single atom type, cubic box) MD code for teaching purposes myself (to teach code optimization, parallelization, and collaborative software development at the HPC master program at GitHub - akohlmey/ljmd-c: Lennard-Jones MD code in C for teaching purposes (an equivalent Fortran version exists as well).

In your description you didn’t make any statement about how you compute the forces, and most importantly there is no mention of how you implement handling of (periodic?) boundaries.

1 Like

Thanks Axel, and my apology for being off-topic here. I am looking at your MD code to see what’s wrong with mine.