1
$\begingroup$

Lennard-Jones force can be computed in two different ways:

Version #1:

Here, the force computation is a self-contained routine.

Vec3 lj_force(const Vec3& r)
{
    double r_mag = std::sqrt(r.x * r.x + r.y * r.y + r.z * r.z);
    double s_over_r = sigma / r_mag;
    double s_over_r6 = s_over_r * s_over_r * s_over_r * s_over_r * s_over_r * s_over_r;
    double s_over_r12 = s_over_r6 * s_over_r6;
    double factor = 24.0 * epsilon * (2.0 * s_over_r12 - s_over_r6) / (r_mag * r_mag * r_mag);

    Vec3 force;
    force.x = factor * r.x;
    force.y = factor * r.y;
    force.z = factor * r.z;
    return force;
}

Version #2:

Here, the force is computed by reusing the potential routine.

// Lennard-Jones potential function
double lj_potential(const Vec3& r)
{
    double r_mag = std::sqrt(r.x * r.x + r.y * r.y + r.z * r.z);
    double s_over_r = sigma / r_mag;
    double s_over_r6 = pow(s_over_r, 6);
    return 4.0 * epsilon * (s_over_r6 * s_over_r6 - s_over_r6);
}

// Derivative of the Lennard-Jones potential
Vec3 lj_force(const Vec3& r)
{
    // Define a small distance for the derivative approximation
    double dr = 1e-6;
    Vec3 force;
    std::vector<Vec3> r_plus_dr = { r, r, r };
    r_plus_dr[0].x += dr;
    r_plus_dr[1].y += dr;
    r_plus_dr[2].z += dr;
    // The force is the negative derivative of the potential energy
    force.x = -(lj_potential(r_plus_dr[0]) - lj_potential(r)) / dr;
    force.y = -(lj_potential(r_plus_dr[1]) - lj_potential(r)) / dr;
    force.z = -(lj_potential(r_plus_dr[2]) - lj_potential(r)) / dr;
    return force;
}

Questions:

  1. Which one is computationally more accutate and why?
  2. Which ine is preferable in Molecular Dynamics simulation and why?
$\endgroup$
1

1 Answer 1

3
$\begingroup$

Which one is computationally more accutate [sic] and why?

The first one is more accurate. It uses the exact derivative of the LJ potential.

By contrast, the second one not only uses a numerical approximation for the derivative, but uses a particular fixed value for the finite difference denominator dr, which may not always be suitable, depending on the value of r and the quirks of floating-point arithmetic.

Which ine [sic] is preferable in Molecular Dynamics simulation and why?

Also the first one, because it is at least as performant as the second version and doesn't require an approximation.

$\endgroup$
1
  • $\begingroup$ I agree completely. He may use the second one to double-check the analytical implementation! $\endgroup$
    – MPIchael
    Commented Sep 14, 2023 at 12:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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