2
\$\begingroup\$

Please review this C++ listing of an implementation of Leapfrog integration.

Was the algorithm implemented correctly?

#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 leapstep(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 = 10000;    // 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++) {
        leapstep(x, v, dt);
    }

    return 0;
}
\$\endgroup\$
2
  • \$\begingroup\$ Correct me if I am wrong, but I don't see how the particles interact. Which acceleration does accel compute? \$\endgroup\$
    – vnp
    Sep 12, 2023 at 21:25
  • \$\begingroup\$ Please don't add your updated code to the question. That can make it harder to match reviews to the code that they are reviewing. See What should I do when someone answers my question? - it's better to ask a new question (if you want further review), post an answer (if you have additional insights) or comment with a link to updated code. I've rolled back to the previous version and recommend you take one of those actions instead. \$\endgroup\$ Sep 15, 2023 at 7:52

1 Answer 1

3
\$\begingroup\$

The leapfrog part itself seems correct.

Use a vector math library

Instead of creating your own Vector3D, use a library that implements a similar type for you, along with overloads for all the mathematical operations you want to do on them. That avoids you from having to spell out all the operations done on the individual components of those vectors. I personally favor GLM (although it is oriented towards graphics it works fine for 3D vectors), but there are many more.

Use exact calculations where possible

We use numerical integration because it is the only practical thing to do when you have a system of many particles. However, you should always favor exact solutions where possible. For example, instead of using a numerical derivative of the Lennard-Jones potential, you should be able to write an exact version. This will then avoid any issues like the choice of dr, which may or may not be right depending on r itself.

What kind of system is this simulating?

vnp already mentioned in the comments that he doesn't see how the particles interact. Indeed, the acceleration of each particle only depends in its distance to the origin of the system, not on any other particle.

What are the units? I see some comments after the declarations of sigma, epsilon and mass, which is great, and from that I can assume that particle positions are also in ångström? But then the box seems very small compared to the size of the well in the Lennard-Jones potential. Even more important, what about time?

Are particles meant to be kept inside the box? Are there periodic boundary conditions? There are lots of questions here. You should document all this in comments in your source code, and/or refer to a paper or other document describing exactly what you are trying to simulate.

Note that these questions are important for deciding whether your use of the leapfrog algorithm is correct: it is just an approximation, and there will are sources of errors. Whether those errors will dominate your results depend on the velocities, forces and the size of your timesteps.

Performance optimizations

While the leapfrom algorithm is implemented correctly, you are calculating the accelarations twice as often as necessary. Consider that the result of the second call to accel() will be the same as the first call to accel() in the next iteration of the for-loop in main().

Related to that: avoid creating a new std::vector holding accelerations for each timestep. Just declare a in main() at the same place you are declaring x and v, and pass it on to leapfrog() by reference.

Create structs to organize your data

There are several ways you can group related data together in structs. For example, you could create a struct System that holds all the data related to the system you are simulating, in particular the particles' positions, velocities and accelerations, a struct Parameters to group all the parameters like sigma, epsilon, n_particles, box_size and so on. You can then easily pass these to initialize() and leapfrog().

You can also create a struct Particle to hold all the data for one particle, so that you only need one std::vector<Particle> particles instead of multiple separate vectors. Usually this make for more intuitive code, but this may or may not be faster.

\$\endgroup\$
5
  • \$\begingroup\$ In the Performance optimizations section "Consider that the result of the second call to accel() will be the same as the first call to accel() in the next iteration of the for-loop in main()" is not clear to me. \$\endgroup\$
    – user366312
    Sep 13, 2023 at 12:37
  • \$\begingroup\$ accel() only depends on the positions of the particles. The positions don't change between the second accel() of the first iteration and the first accel() of the next iteration. \$\endgroup\$
    – G. Sliepen
    Sep 13, 2023 at 12:39
  • \$\begingroup\$ Does that mean I can remove the 2nd call to the accel() function? \$\endgroup\$
    – user366312
    Sep 15, 2023 at 0:30
  • \$\begingroup\$ You can't just delete the line, you have to restructure your code so that you remember the accelerations of the particles between calls to leapstep(). \$\endgroup\$
    – G. Sliepen
    Sep 15, 2023 at 8:38
  • 1
    \$\begingroup\$ codereview.stackexchange.com/questions/287045/… \$\endgroup\$
    – user366312
    Sep 15, 2023 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.