# Computing LJ force from LJ potential, or not?

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?

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.