Where should we compute `Trajectory` in this listing?

Question

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 Trajectory and output it to a CSV file?