Why the lennard-jones potential calculated by lammps seems incorrect?


I’m trying to incorporate my own potential function in lammps. In order to test the accuracy of my code, I did an “experiment”. I put two methane molecules in a box and run a simulation for just 1 timestep. That way I have the potential energy in lammps output. The potential energy is -0.4201 kcal/mol for initial configuration and -0.4193 kcal/mol for timestep 1. Meanwhile, I calculated the potential energy by my pocket calculator and I used the same potential function, epsilon, sigma ,etc. But, now the potential energy is -0.34079kcal/mol for initial configuration and -0.34017kcal/mol for timestep 1. This suggests that my code has a bug and does not calculate the potential correctly.
BUT, I did the same procedure again, and this time I used the regular Lennard-Jones potential instead of my own potential. However, still the results does not match.

Can anybody tell me what is wrong with this experiment that I did? Following, you can see the data file and input file that I used.

Thank you

Data file:

My guess is you didn’t consider periodic boundary conditions in your manual calculation.

Can you seriously imagine a program like lammps being around for such a long time and nobody noticing such a fundamental issue? The problem has to be on your side.