# Where should we compute Energy and Trajectory in this listing?

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


I have three questions in this regard:

1. What is the best place to compute Trajectory and output it to a CSV file?
2. What is the best place to compute Energy and output it to a CSV file?
3. Since computation of LJ potential requires double looping, should I do #1 and #2 this inside accel() function?
• mattermodeling.stackexchange.com/questions/11628/… Commented Sep 15, 2023 at 2:51
• Due to disk space limitations you might want to reduce output in the future, i.e., only write output every X steps. I would keep the output methods separate and in the for loop of your main method. There you can conveniently inject some (if(n%10)==0)) Commented Sep 15, 2023 at 6:43
• Take care to keep track of the units. You calculate distances in au, so are your speeds in au/sec. You have to keep that in mind when calculating the kinetic energy etc. Commented Sep 15, 2023 at 6:47
• Asking the "best place" is poorly framed. If you have tried doing it (printing Trajectory, Energy) one way and got unsatisfactory results, explaining how that was so would make the Question more substantial. At this point you seem open to all sorts of suggestions about writing C++ code which is hardly the forte of this Community.
– hardmath
Commented Sep 15, 2023 at 13:32