1
\$\begingroup\$

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


\$\endgroup\$

2 Answers 2

3
\$\begingroup\$

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\$.

\$\endgroup\$
6
  • \$\begingroup\$ 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. --- this is a too small project to use an external library. \$\endgroup\$
    – user366312
    Commented Sep 15, 2023 at 9:36
  • 1
    \$\begingroup\$ I don't see why the size of your project should matter for that. \$\endgroup\$
    – G. Sliepen
    Commented Sep 15, 2023 at 9:47
  • \$\begingroup\$ Deployment becomes complicated. \$\endgroup\$
    – user366312
    Commented Sep 15, 2023 at 9:50
  • 1
    \$\begingroup\$ That's a valid point, although nowadays it shouldn't be that hard to add some external dependencies to your code. But if it really is problematic, consider using a header-only library, like GLM for example. \$\endgroup\$
    – G. Sliepen
    Commented Sep 15, 2023 at 9:53
  • 3
    \$\begingroup\$ Aha, indeed if it needs to compile for CUDA as well, there are some restrictions. It would have been helpful if you mentioned that in your question. But I also recommend then that you wait until you have something that also compiles for CUDA, and only then post it for review, since it's hard to review unfinished code. \$\endgroup\$
    – G. Sliepen
    Commented Sep 15, 2023 at 9:56
3
\$\begingroup\$

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;
\$\endgroup\$
7
  • \$\begingroup\$ Say something about the leap_step() routine. I failed to modify it as per the previous review. \$\endgroup\$
    – user366312
    Commented Sep 15, 2023 at 10:16
  • \$\begingroup\$ Vect3 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: --- the aim was to keep the code as simple as possible. \$\endgroup\$
    – user366312
    Commented Sep 15, 2023 at 10:17
  • 2
    \$\begingroup\$ The question said nothing about CUDA. \$\endgroup\$ Commented Sep 15, 2023 at 10:30
  • 1
    \$\begingroup\$ Adding operators to Vec3 will make the code simpler - that's why I recommended it. \$\endgroup\$ Commented Sep 15, 2023 at 10:31
  • 1
    \$\begingroup\$ Note that CUDA supports OOP as well (it supports most of C++). Also, CUDA has vector types that work on both the GPU and CPU, so if you are going to target CUDA anyway, just use those. Note that an array of 3D vectors is perfectly fine in CUDA. Arrays of more complicated objects should also work, it's just that it might be less efficient depending how memory is going to be accessed. \$\endgroup\$
    – G. Sliepen
    Commented Sep 15, 2023 at 12:31

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.