The Lennard-Jones pair force¶

The Lennard-Jones 12-6 force between two particles along the line between them, positive apart: V = 4 eps ((s/r)^12 - (s/r)^6), F = -dV/dr = 24 eps (2 (s/r)^12 - (s/r)^6) / r; zero at r = 2^(1/6) s, repulsive inside, attractive beyond.

Source: Allen & Tildesley, Computer Simulation of Liquids, 1.4. The law divides by r: the game keeps a pair apart (two bodies at one point have no direction to push along), as integrator.pendulum_constraint keeps its bob off the pivot. Checked against a reference in numpy, which has its own tests.

Source

Allen & Tildesley, Computer Simulation of Liquids, 1.4.

Proof

curve, limit, invariant.

pair_force.trait

namespace lusine.laws;

float<N> pair_force(float<m> r, float<J> epsilon, float<m> sigma)
{
  let s = sigma / r;
  let s2 = s * s;
  let s6 = s2 * s2 * s2;
  return 24.0 * epsilon * (2.0 * s6 * s6 - s6) / r;
}

Uses, without proving it here: integrator.kick_drift, integrator.kick.

The force between two particles, minus the derivative of the 12-6 potential, with argon's parameters.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
EPS = float(np.finfo(np.float32).eps)

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run EPS = float(np.finfo(np.float32).eps)

The kernel¶

The law as the library states it, and as this page runs it.

In [2]:

print(trait_run.source('laws/pair_force.trait'))

print(trait_run.source('laws/pair_force.trait'))

// lusine.laws.pair_force: the Lennard-Jones 12-6 force between two particles along the line between them, positive
// apart: V = 4 eps ((s/r)^12 - (s/r)^6), F = -dV/dr = 24 eps (2 (s/r)^12 - (s/r)^6) / r; zero at r = 2^(1/6) s,
// repulsive inside, attractive beyond.
// Source: Allen & Tildesley, Computer Simulation of Liquids, 1.4.
// The law divides by r: the game keeps a pair apart (two bodies at one point have no direction to push along), as
// integrator.pendulum_constraint keeps its bob off the pivot.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

float<N> pair_force(float<m> r, float<J> epsilon, float<m> sigma)
{
  let s = sigma / r;
  let s2 = s * s;
  let s6 = s2 * s2 * s2;
  return 24.0 * epsilon * (2.0 * s6 * s6 - s6) / r;
}

The program¶

Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.

In [3]:

print(trait_run.law_program('laws/pair_force.trait', 'pair_force'))

print(trait_run.law_program('laws/pair_force.trait', 'pair_force'))

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host float<m> r;
  host float<J> epsilon;
  host float<m> sigma;
  float<N> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.laws.pair_force(items.r, items.epsilon, items.sigma); }

Curve: the force is -dV/dr¶

The kernel against a central difference of V = 4 eps ((sigma/r)^12 - (sigma/r)^6), in double precision.

In [4]:

eps, sig = 1.65e-21, 3.4e-10
V = lambda r: 4 * eps * ((sig / r) ** 12 - (sig / r) ** 6)
r = np.linspace(0.9 * sig, 2.5 * sig, 60)
h = 1e-15
numeric = -(V(r + h) - V(r - h)) / (2 * h)
force = trait_run.law('laws/pair_force.trait', 'pair_force', r=r, epsilon=eps, sigma=sig)
assert np.max(np.abs(force - numeric)) < 1e-4 * np.max(np.abs(force))
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(r / sig, force, label='pair_force')
ax.plot(r / sig, numeric, '.', ms=3, label='-dV/dr')
ax.axhline(0, lw=0.5)
ax.set_xlabel('r / sigma')
ax.legend()
plt.show()

eps, sig = 1.65e-21, 3.4e-10 V = lambda r: 4 * eps * ((sig / r) ** 12 - (sig / r) ** 6) r = np.linspace(0.9 * sig, 2.5 * sig, 60) h = 1e-15 numeric = -(V(r + h) - V(r - h)) / (2 * h) force = trait_run.law('laws/pair_force.trait', 'pair_force', r=r, epsilon=eps, sigma=sig) assert np.max(np.abs(force - numeric)) < 1e-4 * np.max(np.abs(force)) fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(r / sig, force, label='pair_force') ax.plot(r / sig, numeric, '.', ms=3, label='-dV/dr') ax.axhline(0, lw=0.5) ax.set_xlabel('r / sigma') ax.legend() plt.show()

Limit: zero at 2^(1/6) sigma, and nothing far away¶

The force vanishes at the potential's minimum, to binary32's cancellation there, and falls toward nothing with distance.

In [5]:

eps, sig = 1.65e-21, 3.4e-10
r0 = 2 ** (1 / 6) * sig
F = trait_run.law('laws/pair_force.trait', 'pair_force', r=[r0, 20 * sig, 1.2 * sig, 0.95 * r0, 1.1 * r0],
                  epsilon=eps, sigma=sig)
scale = 24 * eps / sig
assert abs(F[0]) < 8 * EPS * scale and abs(F[1]) < 1e-7 * abs(F[2]) and F[3] > 0 > F[4], F
print(f'at 2^(1/6) sigma: {F[0]:.2e} N, of a scale {scale:.2e} N')

eps, sig = 1.65e-21, 3.4e-10 r0 = 2 ** (1 / 6) * sig F = trait_run.law('laws/pair_force.trait', 'pair_force', r=[r0, 20 * sig, 1.2 * sig, 0.95 * r0, 1.1 * r0], epsilon=eps, sigma=sig) scale = 24 * eps / sig assert abs(F[0]) < 8 * EPS * scale and abs(F[1]) < 1e-7 * abs(F[2]) and F[3] > 0 > F[4], F print(f'at 2^(1/6) sigma: {F[0]:.2e} N, of a scale {scale:.2e} N')

at 2^(1/6) sigma: 2.47e-17 N, of a scale 1.16e-10 N

Invariant: two argon atoms keep their energy¶

Two atoms, as their separation r and its rate with the reduced mass mu = m / 2, stepped by the library's velocity Verlet under the pair force, exchange kinetic and potential energy and keep their sum. Newton's second law is the program's own kernel.

In [6]:

PAIR = '''namespace pair;

lusine.version version = { { number = 1 } };

struct body { float<m> r; float<m/s> velocity; float<m/s2> acceleration; float<N> force; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct { float<J> epsilon; float<m> sigma; float<kg> mass; float<s> dt; } settings;

float<m/s2> newton(float<N> force, float<kg> mass) { return force / mass; }

void start()
{
  bodies.force = lusine.laws.pair_force(bodies.r, settings.epsilon, settings.sigma);
  bodies.acceleration = newton(bodies.force, settings.mass);
}

void main()
{
  (bodies.r, bodies.velocity) = lusine.integrator.kick_drift(position = bodies.r, velocity = bodies.velocity,
                                                             acceleration = bodies.acceleration, dt = settings.dt);
  bodies.force = lusine.laws.pair_force(bodies.r, settings.epsilon, settings.sigma);
  bodies.acceleration = newton(bodies.force, settings.mass);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''
eps, sig, mu, dt = 1.65e-21, 3.4e-10, 6.63e-26 / 2, 2e-15
run = trait_run.Program(['laws/pair_force.trait', 'integrator/kick.trait'], PAIR, {'pair.bodies': [1]})
for name, value in (('epsilon', eps), ('sigma', sig), ('mass', mu), ('dt', dt)):
    run.store('pair.settings', name, value)
run.store('pair.bodies', 'r', 1.5 * sig)
run.store('pair.bodies', 'velocity', -200.0)
run.step('pair.start')
r, v = (a.ravel() for a in run.steps('pair.main', 4000, [('pair.bodies', 'r'), ('pair.bodies', 'velocity')]))
V = lambda r: 4 * eps * ((sig / r) ** 12 - (sig / r) ** 6)
E = 0.5 * mu * v ** 2 + V(r)
E0 = 0.5 * mu * float(np.float32(-200.0)) ** 2 + V(float(np.float32(1.5 * sig)))
worst = np.max(np.abs(E - E0)) / abs(E0)
assert worst < 1e-3 and r.min() < 1.2 * sig, (worst, r.min() / sig)
print(f'closest {r.min() / sig:.3f} sigma; energy within {worst:.1e} of its start')

PAIR = '''namespace pair; lusine.version version = { { number = 1 } }; struct body { float<m> r; float<m/s> velocity; float<m/s2> acceleration; float<N> force; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { float<J> epsilon; float<m> sigma; float<kg> mass; float<s> dt; } settings; float<m/s2> newton(float<N> force, float<kg> mass) { return force / mass; } void start() { bodies.force = lusine.laws.pair_force(bodies.r, settings.epsilon, settings.sigma); bodies.acceleration = newton(bodies.force, settings.mass); } void main() { (bodies.r, bodies.velocity) = lusine.integrator.kick_drift(position = bodies.r, velocity = bodies.velocity, acceleration = bodies.acceleration, dt = settings.dt); bodies.force = lusine.laws.pair_force(bodies.r, settings.epsilon, settings.sigma); bodies.acceleration = newton(bodies.force, settings.mass); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' eps, sig, mu, dt = 1.65e-21, 3.4e-10, 6.63e-26 / 2, 2e-15 run = trait_run.Program(['laws/pair_force.trait', 'integrator/kick.trait'], PAIR, {'pair.bodies': [1]}) for name, value in (('epsilon', eps), ('sigma', sig), ('mass', mu), ('dt', dt)): run.store('pair.settings', name, value) run.store('pair.bodies', 'r', 1.5 * sig) run.store('pair.bodies', 'velocity', -200.0) run.step('pair.start') r, v = (a.ravel() for a in run.steps('pair.main', 4000, [('pair.bodies', 'r'), ('pair.bodies', 'velocity')])) V = lambda r: 4 * eps * ((sig / r) ** 12 - (sig / r) ** 6) E = 0.5 * mu * v ** 2 + V(r) E0 = 0.5 * mu * float(np.float32(-200.0)) ** 2 + V(float(np.float32(1.5 * sig))) worst = np.max(np.abs(E - E0)) / abs(E0) assert worst < 1e-3 and r.min() < 1.2 * sig, (worst, r.min() / sig) print(f'closest {r.min() / sig:.3f} sigma; energy within {worst:.1e} of its start')

closest 0.997 sigma; energy within 6.9e-04 of its start