7
$\begingroup$

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. enter image description here enter image description here

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?

Improve this question
$\endgroup$
5
  • 8
    $\begingroup$ To check if your energy drift stems from your time integrator (you assume it does not, but might want to test anyway), I would try to find out whether the energy drift depends on the time step size. If smaller steps lead to less energy drift, your integrator may be the culprit. $\endgroup$
    – Daniel
    Commented Jul 22, 2015 at 7:28
  • 1
    $\begingroup$ Thanks, I turned off all optimization flags, and tried two tests as you suggested. I edited my question, and provided the plots. It doesn't appear to be the integrator, unless 0.5 fs isn't a small enough to notice the difference from 1.0 fs. I am so perplexed. $\endgroup$
    – user3225087
    Commented Jul 22, 2015 at 20:12
  • $\begingroup$ Yeah, this looks somewhat inconclusive. However, given these plots, are you sure you really have an energy drift? What you plot is more the energy error over time, not really the drift, correct? In the sense that drift would be something like the change of energy error with time? Just asking because a symplectic integrator like Verlet only gives you an energy error that is bounded over time - which your plots do not really rule out, do they? Sorry, while I know a bit about time integrators, I know very little about the specific problem you are solving, so I am not sure if this is helpful. $\endgroup$
    – Daniel
    Commented Jul 23, 2015 at 6:39
  • 1
    $\begingroup$ To elaborate a little more on my comment above: The quantity you define should probably better be called something like relative energy error instead of energy drift. Energy drift would then be the rate of change of the energy error over time. $\endgroup$
    – Daniel
    Commented Jul 23, 2015 at 7:00
  • $\begingroup$ Drift in MD is indeed commonly define as a rate of change of total energy. One should only make such measurements excluding equilibration time, unless you have reason to believe your starting configuration is a sample from the equilibrium ensemble. $\endgroup$
    – mabraham
    Commented Aug 26, 2015 at 20:54

1 Answer 1

Reset to default
2
$\begingroup$

You've measured the relative deviation from the initial energy, whereas "energy drift" in MD is only usefully measured as a rate of absolute deviation, and only after the system has equilibrated. (Neither of your simulations looks long enough to measure that equilibration has occurred, e.g. by measuring statistical properties of windows of the energy distribution.)

To answer your questions

1) You can systematically vary parameters in the program, observe the effect on drift and infer things about the numerical properties of the implementation. (This is what you have already done with the size of the time step.) Other parameters of interest are the atom ordering (e.g. in a water box, if you accumulate all the O-O forces on the O atoms before any of the O-H forces on the O atoms then you will observe different numerical properties than if you accumulate molecule by molecule), simulation cell size (the numerical precision at which you can compute the distance between two atoms around (10,10,10) is rather lower than that you can compute for two atoms around (1,1,1)), and atom masses.

2) As above, such metrics are size dependent, and drifts are often reported per atom or degree of freedom. Whether drift matters depends on what you are observing. Usually you are observing the statistical properties of some observable, and if those are indistinguishable from correct, then your job is done. A useful approach can be seen here https://scholar.google.com/scholar?oi=bibs&cluster=1069342445614725659&btnI=1&hl=en. Most people talk about drift because it's easy to measure and has a well-defined expected property, but if it walks and talks like a duck, it doesn't matter whether it has fleas unless you're going to use the down.

Improve this answer
$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.