I wrote this minimal example to examine the Leapfrog integration algorithm.
However, I am not sure it is the correct algorithm, and the listing is giving the correct output.
Is this the Leapfrog algorithm?
NOTE: I am using the Lennard-Jones system.
#include "stdafx.h"
#include <cmath>
#include <iostream>
class Vec3 {
public:
Vec3() {}
Vec3(double x, double y, double z) : x(x), y(y), z(z) {}
Vec3(const Vec3& other) : x(other.x), y(other.y), z(other.z) {} // Copy constructor
Vec3& operator=(const Vec3& other) {
if (this != &other) {
x = other.x;
y = other.y;
z = other.z;
}
return *this;
} // Copy assignment operator
Vec3 operator+(const Vec3& b) const {
return Vec3(x + b.x, y + b.y, z + b.z);
}
Vec3 operator*(double scalar) const {
return Vec3(scalar * x, scalar * y, scalar * z);
}
// Friend function for scalar multiplication with scalar on the left-hand side
friend Vec3 operator*(double scalar, const Vec3& vector) {
return vector * scalar;
}
double getX() const { return x; }
double getY() const { return y; }
double getZ() const { return z; }
static Vec3 Zero() {
return Vec3(0.0, 0.0, 0.0);
}
private:
double x, y, z;
};
class Atom {
public:
Atom() {}
Atom(const Vec3& position, const Vec3& velocity)
: position(position), velocity(velocity) {}
Atom(const Atom& other)
: position(other.position), velocity(other.velocity) {}
// Copy constructor
Atom& operator=(const Atom& other) {
if (this != &other) {
position = other.position;
velocity = other.velocity;
}
return *this;
}
// Copy assignment operator
Vec3 getPosition() const { return position; }
Vec3 getVelocity() const { return velocity; }
void setPosition(const Vec3& position) { this->position = position; }
void setVelocity(const Vec3& velocity) { this->velocity = velocity; }
private:
Vec3 position;
Vec3 velocity;
};
class ArgonGasSimulation {
private:
static constexpr double Sigma = 3.4;
static constexpr double Epsilon = 1.0;
static Vec3 ComputeForce(const Vec3& position) {
double distance = sqrt(pow(position.getX(), 2) + pow(position.getY(), 2) + pow(position.getZ(), 2));
double r6 = pow(distance, 6);
double r12 = pow(r6, 2);
double magnitude = 48 * Epsilon * (pow((Sigma / distance), 6) - 0.5 * pow((Sigma / distance), 4)) / distance;
return position * magnitude;
}
static void Leapfrog(Atom& atom, double timeStep) {
Vec3 force = ComputeForce(atom.getPosition()); // Compute force based on the current position
Vec3 newPosition = atom.getPosition() + atom.getVelocity() * timeStep + 0.5 * force * pow(timeStep, 2);
Vec3 newForce = ComputeForce(newPosition);
Vec3 newVelocity = atom.getVelocity() + (force + newForce) * timeStep * 0.5;
// Update state
atom.setPosition(newPosition);
atom.setVelocity(newVelocity);
}
public:
static void Main() {
// Initial state
Vec3 position(1.0, 2.0, 3.0); // initial position
Vec3 velocity = Vec3::Zero(); // initial velocity
double timeStep = 0.01;
Atom atom(position, velocity);
// Leapfrog integration
Leapfrog(atom, timeStep);
// Output new position and velocity
std::cout << "New Position: (" << atom.getPosition().getX() << ", " << atom.getPosition().getY() << ", " << atom.getPosition().getZ() << ")" << std::endl;
std::cout << "New Velocity: (" << atom.getVelocity().getX() << ", " << atom.getVelocity().getY() << ", " << atom.getVelocity().getZ() << ")" << std::endl;
}
};
int main() {
ArgonGasSimulation::Main();
return 0;
}
Output:
New Position: (1.00014244382518, 2.00028488765035, 3.00042733147553)
New Velocity: (0.028470372373851, 0.0569407447477019, 0.0854111171215529)
Press any key to continue . . .
EDIT: The following is rewritten in C++ on the basis of the Lutz Lehman's answer.
#include "stdafx.h"
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
// Constants for Argon
constexpr double epsilon = 119.8; // Depth of the potential well (in K)
constexpr double sigma = 3.405; // Distance for zero potential (in Angstrom)
constexpr double mass = 39.948; // Mass of Argon (in amu)
struct Vector3D {
double x, y, z;
};
// Lennard-Jones potential function
double lj_potential(const Vector3D& 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
Vector3D lj_force(const Vector3D& r)
{
// Define a small distance for the derivative approximation
double dr = 1e-6;
Vector3D force;
std::vector<Vector3D> 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;
}
// Update the 'accel' function
void accel(std::vector<Vector3D>& a, const std::vector<Vector3D>& x)
{
int n = x.size(); // number of points
for (int i = 0; i < n; i++)
{
auto force = lj_force(x[i]);
// use Lennard-Jones force law
a[i].x = -force.x / mass;
a[i].y = -force.y / mass;
a[i].z = -force.z / mass;
}
}
void leapfrog_step(std::vector<Vector3D>& x, std::vector<Vector3D>& v, double dt) {
int n = x.size(); // number of points
std::vector<Vector3D> a(n);
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // advance vel by half-step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // advance vel by half-step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // advance vel by half-step
x[i].x = x[i].x + dt * v[i].x; // advance pos by full-step
x[i].y = x[i].y + dt * v[i].y; // advance pos by full-step
x[i].z = x[i].z + dt * v[i].z; // advance pos by full-step
}
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // and complete vel. step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // and complete vel. step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // and complete vel. step
}
}
// Initialize positions and velocities of particles
void initialize(std::vector<Vector3D>& x, std::vector<Vector3D>& v, int n_particles, double box_size, double max_vel) {
// Create a random number generator
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(-0.5, 0.5);
// Resize the vectors to hold the positions and velocities of all particles
x.resize(n_particles);
v.resize(n_particles);
// Initialize positions and velocities
for (int i = 0; i < n_particles; i++)
{
// Assign random initial positions within the box
x[i].x = box_size * distribution(generator);
x[i].y = box_size * distribution(generator);
x[i].z = box_size * distribution(generator);
// Assign random initial velocities up to max_vel
v[i].x = max_vel * distribution(generator);
v[i].y = max_vel * distribution(generator);
v[i].z = max_vel * distribution(generator);
}
}
int main() {
int n_particles = 100; // number of particles
double box_size = 10.0; // size of the simulation box
double max_vel = 0.1; // maximum initial velocity
double dt = 0.01; // time step
int n_steps = 10; // number of time steps
// Positions and velocities of the particles
std::vector<Vector3D> x, v;
// Initialize the particles
initialize(x, v, n_particles, box_size, max_vel);
// Run the simulation
for (int step = 0; step < n_steps; step++) {
leapfrog_step(x, v, dt);
}
return 0;
}
EDIT-2: I wrote the following listing after reading the comments under Lutz Lehman's post.
#include "stdafx.h"
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
// Constants for Argon
constexpr double epsilon = 119.8; // Depth of the potential well (in K)
constexpr double sigma = 3.405; // Distance for zero potential (in Angstrom)
constexpr double mass = 39.948; // Mass of Argon (in amu)
struct Vector3D {
double x, y, z;
};
// Derivative of the Lennard-Jones potential
Vector3D lj_force(const Vector3D& 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);
Vector3D force;
force.x = factor * r.x;
force.y = factor * r.y;
force.z = factor * r.z;
return force;
}
// Update the 'accel' function
void accel(std::vector<Vector3D>& a, const std::vector<Vector3D>& x)
{
int n = x.size(); // number of points
// Reset the acceleration to zero
for (int i = 0; i < n; i++)
{
a[i].x = 0.0;
a[i].y = 0.0;
a[i].z = 0.0;
}
// Compute the acceleration due to each pair of particles
for (int i = 0; i < n; i++)
{
for (int j = i + 1; j < n; j++)
{
Vector3D r;
r.x = x[j].x - x[i].x;
r.y = x[j].y - x[i].y;
r.z = x[j].z - x[i].z;
Vector3D force = lj_force(r);
// use Lennard-Jones force law
a[i].x += force.x / mass;
a[i].y += force.y / mass;
a[i].z += force.z / mass;
a[j].x -= force.x / mass;
a[j].y -= force.y / mass;
a[j].z -= force.z / mass;
}
}
}
void leapfrog_step(std::vector<Vector3D>& x, std::vector<Vector3D>& v, double dt) {
int n = x.size(); // number of points
std::vector<Vector3D> a(n);
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // advance vel by half-step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // advance vel by half-step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // advance vel by half-step
x[i].x = x[i].x + dt * v[i].x; // advance pos by full-step
x[i].y = x[i].y + dt * v[i].y; // advance pos by full-step
x[i].z = x[i].z + dt * v[i].z; // advance pos by full-step
}
accel(a, x);
for (int i = 0; i < n; i++)
{
v[i].x = v[i].x + 0.5 * dt * a[i].x; // and complete vel. step
v[i].y = v[i].y + 0.5 * dt * a[i].y; // and complete vel. step
v[i].z = v[i].z + 0.5 * dt * a[i].z; // and complete vel. step
}
}
// Initialize positions and velocities of particles
void initialize(std::vector<Vector3D>& x, std::vector<Vector3D>& v, int n_particles, double box_size, double max_vel) {
// Create a random number generator
std::default_random_engine generator;
std::uniform_real_distribution<double> distribution(-0.5, 0.5);
// Resize the vectors to hold the positions and velocities of all particles
x.resize(n_particles);
v.resize(n_particles);
// Initialize positions and velocities
for (int i = 0; i < n_particles; i++)
{
// Assign random initial positions within the box
x[i].x = box_size * distribution(generator);
x[i].y = box_size * distribution(generator);
x[i].z = box_size * distribution(generator);
// Assign random initial velocities up to max_vel
v[i].x = max_vel * distribution(generator);
v[i].y = max_vel * distribution(generator);
v[i].z = max_vel * distribution(generator);
}
}
int main() {
int n_particles = 100; // number of particles
double box_size = 10.0; // size of the simulation box
double max_vel = 0.1; // maximum initial velocity
double dt = 0.01; // time step
int n_steps = 10; // number of time steps
// Positions and velocities of the particles
std::vector<Vector3D> x, v;
// Initialize the particles
initialize(x, v, n_particles, box_size, max_vel);
// Run the simulation
for (int step = 0; step < n_steps; step++) {
leapfrog_step(x, v, dt);
}
return 0;
}
The standard Verlet method transforms the central second order difference quotient
(x(t+h)-2x(t)+x(t-h))/h^2 = x''(t) + O(h^2) = a(x(t)) + O(h^2)
into a time propagator
x[i+1] = 2*x[i]-x[i-1] + h^2*a[i]
The inherent double summation can be split into two chained summations using the helper variables
h*v[i+0.5] = x[i+1]-x[i]
These are at first just difference quotients or secant slopes for the point series of the above sequence. They approximate the derivative of the solution best where they are the central difference quotients, that is at the midpoints t[i+0.5]=t[i]+0.5*h. The reduced Verlet recursion is
h*a[i] = v[i+0.5] - v[i-0.5]
or as algorithm
a[i] = acc(x[i])
v[i+0.5] = v[i-0.5] + h*a[i]
x[i+1] = x[i] + h*v[i+0.5]
Methods that combine values from such staggered grids are given the moniker "leapfrog" or "checkerboard".
To get the full state of position and velocity at the original grid points t[i] the next-best central difference quotient is
2*h*v[i] = x[i+1]-x[i-1] = h*(v[i+0.5]+v[i-0.5])
There are now several ways to get the step equation for these velocities. They all end with
v[i+1]-v[i] = 0.5*h*(a[i+1]+a[i])
As algorithm in explicit steps this is realized as
v[i+0.5] = v[i] + 0.5*h*a[i]
x[i+1] = x[i] + h*v[i+0.5]
a[i+1] = acc(x[i+1])
v[i+1] = v[i+0.5] + 0.5*h*a[i+1]
This is now the fully formed symplectic method. As a short name it is called "velocity" Verlet method.
One gets a variant of this by applying a time shift of half a step, so that now the main state sequence is what previously were the states at the half-steps. This gives a different point sequence, but the difference is below the method error, so of no practical relevance
x[i+0.5] = x[i] + 0.5*h*v[i]
a[i+0.5] = acc(x[i+0.5])
v[i+1] = v[i] + h*a[i+0.5]
x[i+1] = x[i+0.5] + 0.5*h*v[i+1]
The acceleration computation and application can now be seen as one unit, the acceleration vector is only relevant for the current step, it does not need be saved for the next step.