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?


  • $\begingroup$ Please don't cross-post on two Stack Exchange sites. Cross-posting is frowned upon. $\endgroup$ – Geoff Oxberry Aug 7 '13 at 22:41
  • 1
    $\begingroup$ 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.. $\endgroup$ – RNs_Ghost Aug 7 '13 at 22:49
  • 1
    $\begingroup$ 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. $\endgroup$ – Geoff Oxberry Aug 7 '13 at 22:58
  • $\begingroup$ Don't worry; I won't cross post again :) $\endgroup$ – RNs_Ghost 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;

Graph validating simulation

| cite | improve this answer | |

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).

See also discussion at Reproducibility in molecular dynamics trajectories

| cite | improve this answer | |

This is not my field but perhaps you could numerically check conservation laws, such as ensemble averaged energy.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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