# Where should we compute Energy in this listing?

Cross-posted here.

I have the following listing for a barebone MD simulation:

#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)
{
//TODO
}

void writeTrajectoryToFile(const std::string& filename, const Trajectory& trajectory)
{
//TODO
}

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


What is the best place to compute Energy and output it to a CSV file?