Minimalistic implementation of Leapfrog integration algorithm (2)

Question

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

This C++ listing is rewritten according to this review.

#include "stdafx.h"
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
#include <fstream>

typedef double real;

// Constants for Argon
constexpr real epsilon = 119.8;   // Depth of the potential well (in K)
constexpr real sigma = 3.405;     // Distance for zero potential (in Angstrom)
constexpr real mass = 39.948;     // Mass of Argon (in amu)

struct Trajectory
{
    int step_no;
    real position_x;
    real position_y;
    real position_z;
    real velocity_x;
    real velocity_y;
    real velocity_z;
    real temperature;
};

struct Energy
{
    int step_no;
    real Total_energy;
    real Poten_energy;
    real Pot_engy_repulsive;
    real Pot_engy_attractive;
    real Pot_engy_balloon;
    real Kinetic_engy;
};

struct Vec3
{
    real x;
    real y;
    real z;
};

struct DataSet
{
    std::vector<real> VecX;
    std::vector<real> VecY;
    std::vector<real> VecZ;
};

void writeEnergyToFile(const std::string& filename, const Energy& energy)
{
    std::ofstream outputFile(filename, std::ios_base::app); // Open the file in append mode

    if (outputFile.is_open())
    {
        // Write the energy values to the file
        outputFile << "Step: " << energy.step_no << std::endl;
        outputFile << "Total Energy: " << energy.Total_energy << std::endl;
        outputFile << "Potential Energy (Repulsive): " << energy.Pot_engy_repulsive << std::endl;
        outputFile << "Potential Energy (Attractive): " << energy.Pot_engy_attractive << std::endl;
        outputFile << "Potential Energy: " << energy.Poten_energy << std::endl;
        outputFile << "Kinetic Energy: " << energy.Kinetic_engy << std::endl;

        outputFile.close(); // Close the file
    }
    else
    {
        std::cout << "Unable to open the file: " << filename.c_str() << std::endl;
    }
}

// Initialize positions and velocities of particles
void initialize(int n_particles, DataSet& positionData, DataSet& velocityData, real boxSize, real maxVelocity)
{
    // Create a random number generator
    std::default_random_engine generator;
    std::uniform_real_distribution<real> distribution(-0.5, 0.5);

    // Initialize positions and velocities
    for (int i = 0; i < n_particles; i++)
    {
        // Assign random initial positions within the box
        positionData.VecX[i] = boxSize * distribution(generator);
        positionData.VecY[i] = boxSize * distribution(generator);
        positionData.VecZ[i] = boxSize * distribution(generator);

        // Assign random initial velocities up to max_vel
        velocityData.VecX[i] = maxVelocity * distribution(generator);
        velocityData.VecY[i] = maxVelocity * distribution(generator);
        velocityData.VecZ[i] = maxVelocity * distribution(generator);
    }
}

// Derivative of the Lennard-Jones potential
Vec3 lj_force(Vec3 posVec)
{
    real r_mag = std::sqrt(posVec.x * posVec.x + posVec.y * posVec.y + posVec.z * posVec.z);
    real s_over_r = sigma / r_mag;
    real s_over_r6 = s_over_r * s_over_r * s_over_r * s_over_r * s_over_r * s_over_r;
    real s_over_r12 = s_over_r6 * s_over_r6;
    real factor = 24.0 * epsilon * (2.0 * s_over_r12 - s_over_r6) / (r_mag * r_mag * r_mag);

    Vec3 force;
    force.x = factor * posVec.x;
    force.y = factor * posVec.y;
    force.z = factor * posVec.z;
    return force;
}

// Update the 'accel' function
void accel(int n_particles, DataSet& accelData, DataSet& posData)
{
    // Reset the acceleration to zero
    for (int i = 0; i < n_particles; i++)
    {
        accelData.VecX[i] = 0.0;
        accelData.VecY[i] = 0.0;
        accelData.VecZ[i] = 0.0;
    }

    // Compute the acceleration due to each pair
    for (int i = 0; i < n_particles; i++)
    {
        for (int j = i + 1; j < n_particles; j++)
        {
            Vec3 posVec;
            posVec.x = posData.VecX[j] - posData.VecX[i];
            posVec.y = posData.VecY[j] - posData.VecY[i];
            posVec.z = posData.VecZ[j] - posData.VecZ[i];
            Vec3 force = lj_force(posVec);
            // use Lennard-Jones force law
            accelData.VecX[i] += force.x / mass;
            accelData.VecY[i] += force.y / mass;
            accelData.VecZ[i] += force.z / mass;

            accelData.VecX[j] -= force.x / mass;
            accelData.VecY[j] -= force.y / mass;
            accelData.VecZ[j] -= force.z / mass;
        }
    }
}

void leapfrog_step(int n_particles, DataSet& posData, DataSet& velocData, real dt)
{
    DataSet a;

    accel(n_particles, a, posData); //compute acceleration
    for (int i = 0; i < n_particles; i++)
    {
        velocData.VecX[i] = velocData.VecX[i] + 0.5 * dt * a.VecX[i];    // advance vel by half-step
        velocData.VecY[i] = velocData.VecY[i] + 0.5 * dt * a.VecY[i];    // advance vel by half-step
        velocData.VecZ[i] = velocData.VecZ[i] + 0.5 * dt * a.VecZ[i];    // advance vel by half-step

        posData.VecX[i] = posData.VecX[i] + dt * velocData.VecX[i];      // advance pos by full-step
        posData.VecY[i] = posData.VecY[i] + dt * velocData.VecY[i];      // advance pos by full-step
        posData.VecZ[i] = posData.VecZ[i] + dt * velocData.VecZ[i];      // advance pos by full-step
    }

    /////////accel(n_particles, a, posData); //compute acceleration
    for (int i = 0; i < n_particles; i++)
    {
        velocData.VecX[i] = velocData.VecX[i] + 0.5 * dt * a.VecX[i];    // and complete vel. step
        velocData.VecY[i] = velocData.VecY[i] + 0.5 * dt * a.VecY[i];    // and complete vel. step
        velocData.VecZ[i] = velocData.VecZ[i] + 0.5 * dt * a.VecZ[i];    // and complete vel. step
    }
}

int main()
{
    int n_particles = 10;  // number of particles
    real box_size = 10.0; // size of the simulation box
    real max_vel = 0.1;   // maximum initial velocity
    real dt = 0.01;   // time step
    int n_steps = 10000;   // number of time steps

    DataSet posData; // Positions of the particles
    DataSet velData; // Velocities of the particles

                     // Initialize the particles
    initialize(n_particles, posData, velData, box_size, max_vel);

    // Run the simulation
    for (int step = 0; step < n_steps; step++)
    {
        leapfrog_step(n_particles, posData, velData, dt);
    }

    return 0;
}

Answer 1

Use a vector math library

As already mentioned in the previous review, you should really go and use a library that defines vector types for you. This will make your life much easier and will simplify the code a lot. It will probably have performance benefits as well.

For example, using the Eigen library, you would write:

using Vec3 = Eigen::Vector3D;

Vec3 lj_force(Vec3 posVec)
{
    real r_mag = posVec.norm();
    …
    real factor = …;
    return factor * posVec;
}
…
void accel(…)
{
    …
    auto posVec = posData[j] - posData[i];
    auto force = lj_force(posVec);
    accelData[i] += force / mass;
    accelData[j] -= force / mass;
    …
}

Note how you no longer have to access .x, .y and .z individually anymore.

Organizing your data

You have struct Vec3, but you are hardly using it. You should have used this (or better, a vector type from a library like mentioned above) everywhere that you have a 3D vector. So for example:

struct Trajectory 
{
    …
    Vec3D position;
    Vec3D velocity;
    …
};

using DataSet = std::vector<Vec3D>;

The struct Trajectory does not contain a trajectory, it just contains the current state of a single particle. So Particle would have been a better name for this. I would have added a Vec3D acceleration member variable, and removed step_no and temperature: the latter two can be derived from the context and from the particle's velocity.

What is the unit of time used in the simulation?

I see that you have some specific numbers for the properties of Argon atoms, along with their units. However, I don't see anywhere mentioned what the unit of time is, and the values for the maximum velocity and timestep seem to be chosen rather arbitrarily. Make sure you define that, and ensure that the velocities and timesteps are chosen such that the errors in your simulation are not too large.

The leapstep algorithm is now incorrect

In my previous review I mentioned that you only need to update the acceleration of the particles in leapfrog_step() once. However, you are now doing it at the wrong place. The acceleration should match that of where the particles are. So, you should update the acceleration right after you have updated all the positions.

It might be good to revisit the leapfrog integration algorithm, and pay attention to the subscripts \$i\$, \$i+1/2\$ and \$i+1\$ used for the position \$x\$, velocity \$v\$ and acceleration \$a\$.

Answer 2

I don't know what "stdafx.h" is - I can't find a package providing it. However, it doesn't seem to be necessary - the code compiles here if I remove that #include.

---

GCC gives an Effective C++ warning for uninitialised members in DataSet. It's a bit over-enthusiastic, because std::vector default-initialises, but it doesn't hurt to add = {} to those declarations to silence the warning.

That said, it looks like we have parallel arrays of components there - usually better represented as a single array of objects, probably std::vector<Vec3>.

---

There's no need to flush outputFile (using std::endl) so frequently. Just write an ordinary newline, and let close() do the flushing.

---

After outputFile.close(), we should test its state to take appropriate action if unsuccessful. BTW // Close the file is what we call a pointless comment.

---

std::sqrt(posVec.x * posVec.x + posVec.y * posVec.y + posVec.z * posVec.z) can be rewritten std::hypot(posVec.x, posVec.y, posVec.z), and that function is more robust against overflow and underflow.

---

Vec3 could do with some operators, so that we can write return posVec * factor instead of this C-like code, and simplify similar member-wise code elsewhere:

Vec3 force;
force.x = factor * posVec.x;
force.y = factor * posVec.y;
force.z = factor * posVec.z;
return force;