I hope the following question will not be perceived as to vague. I am trying to ask the question directly, without going into too much information regarding my code. I am currently developing a molecular dynamics code for a custom application, as part of my PhD thesis. I have checked over and over again the validity of the force subroutines (both non bonded and bonded) by comparing numerical values directly from lammps. The integration routines were the same as I used previously to reproduce LJ equation of state data, so I am operating on the assumption they are a correct form of the verlet algorithm. However, I still have a sneaking suspicion I am experiencing more energy drift than I should be. For instance, a recent test I just ran is 1 glycine molecule in an 92.1*92.1*92.1 angstrom box. I am using OPLS force field with wolf summation for the electrostatics. I defined the drift as
drift = (total potential + kinetic - total initial energy)/total initial energy
After 50,000 steps (1.0 fs time step, no shake, just brute force MD), the drift value is upwards of 0.07 (7%). To test whether or not my assumption about the integrator being correct was true, I ran a simulation for 100k iterations with a tilmestep of 0.5 fs and one with 1.0 fs. I have provided the energy drift plots. For these simulations, the optimization flag was set to -O0. The forces are converted internally from kcal/mol/A to amu*A/fs^2 by multiplying dthalf by 0.0004184.
Given that I am running in full double precision,
1) How is it possible to determine possible sources of energy drift in the code base? For instance, if I were comparing the same simulation ran with lammps and my own, saw that my drift was greater, how can I go about mechanistically to determine what possible sources of this error are in my lines of code? What tests should I be running?
2) What is the rule of thumb for acceptable energy drift? For instance, how can I quantify whether a simulation is following "good" energy drift or "poor" energy drift. Is this metric system size dependent?