2

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

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

3 Answers3

3

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

Juan M. Bello-Rivas
  • 3,994
  • 16
  • 24
2

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

Community
  • 1
mabraham
  • 297
  • 2
  • 7
1

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

JohnS
  • 129
  • 4