# Molecular Dynamics (MD) Simulation: energy fluctuations in NVE ensemble

I'm writing my first MD simulation (ever) for liquid Argon. The code is up and running. I am supposed to do the calculations in the NVE ensemble. Having implemented a 4th order symplectic integrator (forest-ruth) the total energy of my system (approx 1000 atoms) oscillates (as expected).

Now I am "measuring" the heat capacity of my system, which I do by calculating the mean total energy and it's fluctuations for several distinct simulations (using the same parameters).

This is exactly the point that I don't get (conceptually). By definition, the total energy of a system is fixed in the microcanonical ensemble. Am I in the microcanonical or in the canonical ensemble when I do an MD simulation? The literature (I have) is somewhat sloppy when it comes to this.

• ok, so after fixing the code, the total energy does not oscillate (as mentioned by Nick below). I was confused about this one before, because my reference was a harmonic oscillator, duhh... – seb Mar 23 '13 at 7:05
• Way too late, but Symplectic integrators will have an oscillating energy because they conserve an approximate energy. Not the exact hamiltonian, but one with a few perturbation terms. – aidan.plenert.macdonald Sep 7 '17 at 3:12

A month (and lots of lines of code) later, I have finally gained a better understanding of what I have been asking here previously. The book by Frenkel and Smit shed some light into this: when you are running an MD simulation, you are doing this in in something that is very close to the NVE ensemble. It is possible to run an MD simulation in something that is very close to the NVT ensemble by carefully implementing a thermostat.

Contrarily to what Nick has answered, I did measure the heat capacity through fluctuations, which produced very good results and is also recommended in the literature [Frenkel,Smit: Understanding MolecularSimulation ; Allen,Tildesley: Computer Simulation of Liquids].

I compared those measurements to measuring the heat capacity by interpolating energy measurements at different temperatures. The results were clearly inferior to computing the heat capacity from fluctuations.

If you are interested, here is my full report on the project and here is the source code in Python.

What I would do is set up a comparison to a well-known MD code so you can be sure that your code matches up with something that has been time-tested (I use LAMMPS personally).

You say the total energy of your system oscillates. This isn't supposed to happen in the NVE ensemble. I am guessing your timestep may be too big.

If you were to run an NVT (canonical) ensemble, then you would need to implement a thermostat of some sort, like the Nose-Hoover or the Berendsen thermostat. These constantly adjust the particles' velocities so to stay near a certain temperature.

Also, for a system of 1000 atoms, measuring the heat capacity using fluctuations isn't the best way because your error bars will probably be on the same order as your results.

• thanks for the advice. I'm checking out LAMMPS right now. the only other way of measuring the heat capacity that comes to my mind right now is running the simulation for different temperatures and determine the heat capacity from that. would you say that this is a better choice? do you have any other suggestions? – seb Feb 15 '13 at 17:37
• @Nick I wish this answer was around back when I was doing MD half a decade ago! :P – dearN Mar 5 '13 at 16:26