# How can I easily verify my molecular dynamics simulation?

I have written a molecular dynamics simulator with a single potential function (lennard jones). I initialise the temperature of the system using a Boltzmann velocity distribution.

Is there a standard way of verifying that it works? E.g. some pressure/temperature calculation for a known element?

Thanks

• Please don't cross-post on two Stack Exchange sites. Cross-posting is frowned upon. Aug 7 '13 at 22:41
• I have deleted the post on physics.SE. The only reason I x-posted is due to 1)the low traffic volume on this SE 2)In the past Questions have gone unanswered on this SE necessitating a post on physics.SE.. Aug 7 '13 at 22:49
• If your post is going unanswered for some period of time (say, a week or two), I think it's reasonable to cross-post to another relevant SE site (along with a comment below the post explaining why you're cross-posting). A better practice may be to ask a mod to migrate your question if it goes unanswered after a while. Cross-posting immediately after posting your question is going to cause moderators to point out that you shouldn't cross-post. Aug 7 '13 at 22:58
• Don't worry; I won't cross post again :) Aug 7 '13 at 23:06

A common test consists of computing some ensemble averages that you can check analytically and repeating the experiment a number of times while decreasing the time step length. You should observe that, when you plot the time step length against the error in log-log scale, the data points can be easily fitted to a line with slope equal to the order of the numerical scheme's global error (that would be two in the case of Verlet).

The following code computes an integral curve of a Lennard-Jones oscillator using the Verlet algorithm for different time step lengths and collecting the average total energy (which should be constant). In the plot below a simple linear fitting reveals that the slope is almost two, as expected from theory.

#include <cmath>
#include <iostream>

using namespace std;

const double total_time = 1e3;

const double sigma = 1.0;
const double sigma6 = pow(sigma, 6.0);
const double epsilon = 1.0;
const double four_epsilon = 4.0 * epsilon;

inline double force(double q, double& potential) {
const double r2 = q * q;
const double r6 = r2 * r2 * r2;
const double factor6  = sigma6 / r6;
const double factor12 = factor6 * factor6;

potential = four_epsilon * (factor12 - factor6);
return -four_epsilon * (6.0 * factor6 - 12.0 * factor12) / r2 * q;
}

int main() {
const double q0 = 1.5, p0 = 0.1;
double potential;
const double f0 = force(q0, potential);
const double total_energy_exact = p0 * p0 / 2.0 + potential;

for (double dt = 1e-3; dt <= 5e-2; dt *= 2.0) {
const long steps = long(total_time / dt);

double q = q0, p = p0, f = f0;
double total_energy_average = total_energy_exact;

for (long step = 1; step <= steps; ++step) {
p += dt / 2.0 * f;
q += dt * p;
f = force(q, potential);
p += dt / 2.0 * f;

total_energy_average += p * p / 2.0 + potential;
}

total_energy_average /= double(steps);

const double err = fabs(total_energy_exact - total_energy_average);
cout << log10(dt) << "\t"
<< log10(err) << endl;
}

return 0;
} One simple approach is to recognize that the distribution of (say) potential energies over the trajectory is known a priori, and varies in a known way with (say) temperature. A statistically significant observation that the same (and the expected) temperature dependence of the distribution is observed with the two integrators is an excellent start to demonstrating their equivalence. See http://dx.doi.org/10.1021/ct300688p (article also available on arxiv.org; code on https://simtk.org/home/checkensemble).