# energy drift in molecular dynamics

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?

• 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. – Daniel Jul 22 '15 at 7:28
• 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. – user3225087 Jul 22 '15 at 20:12
• 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. – Daniel Jul 23 '15 at 6:39
• 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. – Daniel Jul 23 '15 at 7:00
• 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. – mabraham Aug 26 '15 at 20:54
{"ft_lang_label":"__label__en","ft_lang_prob":0.9350725,"math_prob":0.82889193,"size":2056,"snap":"2019-51-2020-05","text_gpt3_token_len":472,"char_repetition_ratio":0.10818713,"word_repetition_ratio":0.0,"special_character_ratio":0.22762646,"punctuation_ratio":0.103365384,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.95061874,"pos_list":[0,1,2,3,4],"im_url_duplicate_count":[null,3,null,3,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2019-12-13T11:44:38Z\",\"WARC-Record-ID\":\"<urn:uuid:1e960275-68d1-4e41-8591-02f24dfcfb51>\",\"Content-Length\":\"139827\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:fc899335-cb01-49d3-810a-88d290651321>\",\"WARC-Concurrent-To\":\"<urn:uuid:18717727-0504-441b-88e6-9e87f46187fa>\",\"WARC-IP-Address\":\"151.101.1.69\",\"WARC-Target-URI\":\"https://scicomp.stackexchange.com/questions/20252/energy-drift-in-molecular-dynamics\",\"WARC-Payload-Digest\":\"sha1:TGOXJKJWZB5O7OUUOFJY4U5SFLA6DLWB\",\"WARC-Block-Digest\":\"sha1:DIYQHER27ZKYLCRI2C735LPHMIXJV7JG\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2019/CC-MAIN-2019-51/CC-MAIN-2019-51_segments_1575540553486.23_warc_CC-MAIN-20191213094833-20191213122833-00455.warc.gz\"}"}