A harmonic bond¶

A harmonic bond, Hooke's law F = -k (r - r0) along the bond, positive apart: pulled together when stretched, pushed apart when compressed, nothing at rest.

Source: Allen & Tildesley, Computer Simulation of Liquids, 1.3. Checked against a reference in numpy, which has its own tests.

Source

Allen & Tildesley, Computer Simulation of Liquids, 1.3.

Proof

curve, limit, invariant.

bond_force.trait

namespace lusine.laws;

float<N> bond_force(float<m> r, float<N/m> stiffness, float<m> rest)
{
  return -stiffness * (r - rest);
}

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

The force of a spring between two bonded particles, Hooke's law F = -k (r - r0), positive apart.

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/bond_force.trait'))

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

// lusine.laws.bond_force: a harmonic bond, Hooke's law F = -k (r - r0) along the bond, positive apart: pulled
// together when stretched, pushed apart when compressed, nothing at rest.
// Source: Allen & Tildesley, Computer Simulation of Liquids, 1.3.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

float<N> bond_force(float<m> r, float<N/m> stiffness, float<m> rest)
{
  return -stiffness * (r - rest);
}

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/bond_force.trait', 'bond_force'))

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

namespace probe;

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

struct item {
  host float<m> r;
  host float<N/m> stiffness;
  host float<m> rest;
  float<N> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.laws.bond_force(items.r, items.stiffness, items.rest); }

Curve: a straight line through the rest length¶

The force against -k (r - r0), within binary32's rounding.

In [4]:

k, r0 = 50.0, 1.0
r = np.linspace(0.5, 1.5, 21)
F = trait_run.law('laws/bond_force.trait', 'bond_force', r=r, stiffness=k, rest=r0)
assert np.allclose(F, -k * (r - r0), rtol=4 * EPS, atol=4 * EPS * k)
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(r, F, 'o-')
ax.axhline(0, lw=0.5)
ax.set_xlabel('r (m)')
ax.set_ylabel('F (N)')
plt.show()

k, r0 = 50.0, 1.0 r = np.linspace(0.5, 1.5, 21) F = trait_run.law('laws/bond_force.trait', 'bond_force', r=r, stiffness=k, rest=r0) assert np.allclose(F, -k * (r - r0), rtol=4 * EPS, atol=4 * EPS * k) fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(r, F, 'o-') ax.axhline(0, lw=0.5) ax.set_xlabel('r (m)') ax.set_ylabel('F (N)') plt.show()

Limit: at rest, no force; stretched, pulled back¶

The rest length is the equilibrium, and the force always points back toward it.

In [5]:

F = trait_run.law('laws/bond_force.trait', 'bond_force', r=[1.0, 1.2, 0.8], stiffness=50.0, rest=1.0)
assert F[0] == 0.0 and F[1] < 0 < F[2], F
print('restoring:', F)

F = trait_run.law('laws/bond_force.trait', 'bond_force', r=[1.0, 1.2, 0.8], stiffness=50.0, rest=1.0) assert F[0] == 0.0 and F[1] < 0 < F[2], F print('restoring:', F)

restoring: [ -0.         -10.00000191   9.99999905]

Invariant: a bond oscillates at sqrt(k/m) and keeps its energy¶

A mass on the bond, stepped by the library's velocity Verlet (lusine.integrator.kick_drift, the bond's force, kick), oscillates with period 2 pi sqrt(m/k) and holds its energy. Newton's second law, the force over the mass, is the program's own kernel.

In [6]:

SPRING = '''namespace spring;

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<N/m> stiffness; float<m> rest; 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.bond_force(bodies.r, settings.stiffness, settings.rest);
  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.bond_force(bodies.r, settings.stiffness, settings.rest);
  bodies.acceleration = newton(bodies.force, settings.mass);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''
k, r0, m, dt, steps = 50.0, 1.0, 2.0, 1e-3, 6000
run = trait_run.Program(['laws/bond_force.trait', 'integrator/kick.trait'], SPRING, {'spring.bodies': [1]})
for name, value in (('stiffness', k), ('rest', r0), ('mass', m), ('dt', dt)):
    run.store('spring.settings', name, value)
run.store('spring.bodies', 'r', 1.1)
run.step('spring.start')
r, v = (a.ravel() for a in run.steps('spring.main', steps, [('spring.bodies', 'r'), ('spring.bodies', 'velocity')]))
E = 0.5 * m * v ** 2 + 0.5 * k * (r - r0) ** 2
E0 = 0.5 * k * 0.1 ** 2
assert np.max(np.abs(E - E0)) < 1e-4 * E0, np.max(np.abs(E - E0)) / E0
x = r - r0
down = np.nonzero((x[:-1] > 0) & (x[1:] <= 0))[0]
t = dt * (down + 1 + x[down] / (x[down] - x[down + 1]))       # each downward crossing, between its steps
T = (t[-1] - t[0]) / (len(t) - 1)
assert abs(T / (2 * np.pi * np.sqrt(m / k)) - 1) < 1e-4
print(f'period {T:.5f} s against {2 * np.pi * np.sqrt(m / k):.5f}; energy within '
      f'{np.max(np.abs(E - E0)) / E0:.1e} of its start')

SPRING = '''namespace spring; 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<N/m> stiffness; float<m> rest; 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.bond_force(bodies.r, settings.stiffness, settings.rest); 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.bond_force(bodies.r, settings.stiffness, settings.rest); bodies.acceleration = newton(bodies.force, settings.mass); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' k, r0, m, dt, steps = 50.0, 1.0, 2.0, 1e-3, 6000 run = trait_run.Program(['laws/bond_force.trait', 'integrator/kick.trait'], SPRING, {'spring.bodies': [1]}) for name, value in (('stiffness', k), ('rest', r0), ('mass', m), ('dt', dt)): run.store('spring.settings', name, value) run.store('spring.bodies', 'r', 1.1) run.step('spring.start') r, v = (a.ravel() for a in run.steps('spring.main', steps, [('spring.bodies', 'r'), ('spring.bodies', 'velocity')])) E = 0.5 * m * v ** 2 + 0.5 * k * (r - r0) ** 2 E0 = 0.5 * k * 0.1 ** 2 assert np.max(np.abs(E - E0)) < 1e-4 * E0, np.max(np.abs(E - E0)) / E0 x = r - r0 down = np.nonzero((x[:-1] > 0) & (x[1:] <= 0))[0] t = dt * (down + 1 + x[down] / (x[down] - x[down + 1])) # each downward crossing, between its steps T = (t[-1] - t[0]) / (len(t) - 1) assert abs(T / (2 * np.pi * np.sqrt(m / k)) - 1) < 1e-4 print(f'period {T:.5f} s against {2 * np.pi * np.sqrt(m / k):.5f}; energy within ' f'{np.max(np.abs(E - E0)) / E0:.1e} of its start')

period 1.25664 s against 1.25664; energy within 2.4e-05 of its start