I wrote this minimal example to examine the Leapfrog integration algorithm.

我编写了这个最小的示例来检查LeapFrog积分算法。

However, I am not sure it is the correct algorithm, and the listing is giving the correct output.

但是，我不确定这是正确的算法，并且列表给出了正确的输出。

Is this the Leapfrog algorithm?

这是LeapFrog算法吗？

NOTE: I am using the Lennard-Jones system.

注：我使用的是Lennard-Jones系统。

#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.

编辑：以下内容是根据Lutz Lehman的答案用C++重写的。

#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.

编辑2：在阅读了Lutz Lehman的帖子下的评论后，我写了以下清单。

#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;
}

On the question of why the code does not implement leapfrog Verlet and what that should be

The standard Verlet method transforms the central second order difference quotient

标准的Verlet方法将中心二阶差商

(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

使用helper变量可以将固有的双重求和拆分成两个链式求和

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

首先，这些只是上述序列的点序列的差商或割斜率。它们最好地逼近解的导数，其中它们是中心差商，也就是在中点t[i+0.5]=t[i]+0.5*h。

h*a[i] = v[i+0.5] - v[i-0.5]

or as algorithm

或AS算法

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".

将这种交错网格中的值组合在一起的方法被命名为“LeapFrog”或“Checkerboard”。

To get the full state of position and velocity at the original grid points t[i] the next-best central difference quotient is

要获得原始网格点t[i]上的位置和速度的完整状态，次优中心差商为

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.

这就是现在完全形成的辛方法。简称为“速度”Verlet法。

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.

现在可以将加速度计算和应用看作一个单元，加速度向量只与当前步骤相关，不需要保存以供下一步使用。