# Understanding leapfrog integration algorithm

The leapfrog.cpp is an implementation of leapfrog integration algorithm where f() function is being integrated:

## leapfrog.cpp

#include <cmath>
#include <vector>

struct Particle
{
double x, y, z; // positions
double vx, vy, vz; // velocities
};

// Function to calculate f(x, y, z)
void f(const Particle& particle, double& fx, double& fy, double& fz) {
// Compute forces based on the particle's position
// Modify this function to calculate the forces based on your specific problem
fx = -particle.x;
fy = -particle.y;
fz = -particle.z;
}

void leapfrogIntegration(std::vector<Particle>& particles, int k, int N)
{
double h = 2 * M_PI / N; // step size

// Initial conditions
for (auto& particle : particles)
{
particle.vx = 1.0;
particle.vy = 0.0;
particle.vz = 0.0;
}

// Leapfrog method
for (int i = 1; i <= k * N; i++)
{
for (auto& particle : particles)
{
double fx, fy, fz;

f(particle, fx, fy, fz);

particle.x += particle.vx * h + 0.5 * fx * h * h;
particle.y += particle.vy * h + 0.5 * fy * h * h;
particle.z += particle.vz * h + 0.5 * fz * h * h;

particle.vx += 0.5 * fx * h;
particle.vy += 0.5 * fy * h;
particle.vz += 0.5 * fz * h;

f(particle, fx, fy, fz);

particle.vx += 0.5 * fx * h;
particle.vy += 0.5 * fy * h;
particle.vz += 0.5 * fz * h;
}
}
}

The routine clearly demonstrates how the integrator should work, and the function being integrated is also obvious.

.

Now please take a look at the MD simulation listing MD.cpp.

There you will find a routine integrate() which integrated Newton's equation of motion.

void integrate(const bool isStepOne, const double timeStep, std::vector<Atom>& atoms)
{
const double timeStepHalf = timeStep * 0.5;

for (Atom& atom : atoms)
{
const double ax = atom.fx / atom.mass;
const double ay = atom.fy / atom.mass;
const double az = atom.fz / atom.mass;
if (isStepOne)
{
atom.vx += ax * timeStepHalf;
atom.vy += ay * timeStepHalf;
atom.vz += az * timeStepHalf;
atom.x += atom.vx * timeStep;
atom.y += atom.vy * timeStep;
atom.z += atom.vz * timeStep;
}
else
{
atom.vx += ax * timeStep;
atom.vy += ay * timeStep;
atom.vz += az * timeStep;
}
}
}

I have a few confusions in this regard:

1 Is this a correct implementation of leapfrog integration algorithm?
2. Why does this routine look so radically different than the first one?
3. Where and how is the Newton's equation of motion being integrated (coz, I don't see any such function being called inside the routine)?

## MD.cpp

#include <cmath>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <sstream>
#include <string>
#include <vector>

const int Ns = 100;
const double K_B = 8.617343e-5;
const double TIME_UNIT_CONVERSION = 1.018051e+1;

struct Atom {
int number;
double box[6];
double pe;
double mass, x, y, z, vx, vy, vz, fx, fy, fz;
};

struct Configuration {
int diagnosticLevel;
int numStepsThermalization;
int numStepsHeating;
int numStepsCooling;
int energyOutputFrequency;
int coordinateOutputFrequency;
};

double computeKineticEnergy(const std::vector<Atom>& atoms)
{
double kineticEnergy = 0.0;
for (const Atom& atom : atoms) {
double v2 = atom.vx * atom.vx + atom.vy * atom.vy + atom.vz * atom.vz;
kineticEnergy += atom.mass * v2;
}
return kineticEnergy * 0.5;
}

void scaleVelocity(const double T0, std::vector<Atom>& atoms)
{
const double temperature =
computeKineticEnergy(atoms) * 2.0 / (3.0 * K_B * atoms.size());
double scaleFactor = sqrt(T0 / temperature);
for (Atom& atom : atoms) {
atom.vx *= scaleFactor;
atom.vy *= scaleFactor;
atom.vz *= scaleFactor;
}
}

void initializeVelocity(const double T0, std::vector<Atom>& atoms)
{
#ifndef DEBUG
srand(time(NULL));
#endif
double totalMass = 0.0;
double vSquareSum = 0.0;

for (Atom& atom : atoms) {
atom.vx = -6.0 + (rand() * 12.0) / RAND_MAX;
atom.vy = -6.0 + (rand() * 12.0) / RAND_MAX;
atom.vz = -6.0 + (rand() * 12.0) / RAND_MAX;
vSquareSum += atom.vx * atom.vx + atom.vy * atom.vy + atom.vz * atom.vz;
totalMass += atom.mass;
}

const double scaleFactor = sqrt((2.0 * T0 * K_B * totalMass) / vSquareSum);

for (Atom& atom : atoms) {
atom.vx *= scaleFactor;
atom.vy *= scaleFactor;
atom.vz *= scaleFactor;
}

double centerOfMassVelocity[3] = {0.0, 0.0, 0.0};
for (const Atom& atom : atoms) {
centerOfMassVelocity[0] += atom.mass * atom.vx;
centerOfMassVelocity[1] += atom.mass * atom.vy;
centerOfMassVelocity[2] += atom.mass * atom.vz;
}
centerOfMassVelocity[0] /= totalMass;
centerOfMassVelocity[1] /= totalMass;
centerOfMassVelocity[2] /= totalMass;
for (Atom& atom : atoms) {
atom.vx -= centerOfMassVelocity[0];
atom.vy -= centerOfMassVelocity[1];
atom.vz -= centerOfMassVelocity[2];
}
}

void applyMicOne(const double length, const double halfLength, double& x12)
{
if (x12 < -halfLength)
x12 += length;
else if (x12 > +halfLength)
x12 -= length;
}

void applyMic(const double box[6], double& x12, double& y12, double& z12)
{
applyMicOne(box[0], box[3], x12);
applyMicOne(box[1], box[4], y12);
applyMicOne(box[2], box[5], z12);
}

void computeForce(std::vector<Atom>& atoms)
{
const double epsilon = 1.032e-2;
const double sigma = 3.405;
const double cutoff = 9.0;
const double cutoffSquare = cutoff * cutoff;
const double sigma3 = sigma * sigma * sigma;
const double sigma6 = sigma3 * sigma3;
const double sigma12 = sigma6 * sigma6;
const double e24s6 = 24.0 * epsilon * sigma6;
const double e48s12 = 48.0 * epsilon * sigma12;
const double e4s6 = 4.0 * epsilon * sigma6* pow(rinv, 6);
const double e4s12 = 4.0 * epsilon * sigma12;

for (Atom& atom : atoms)
{
atom.pe = 0.0;
atom.fx = atom.fy = atom.fz = 0.0;
}

for (int i = 0; i < atoms.size() - 1; ++i) {
for (int j = i + 1; j < atoms.size(); ++j) {
double xij = atoms[j].x - atoms[i].x;
double yij = atoms[j].y - atoms[i].y;
double zij = atoms[j].z - atoms[i].z;
applyMic(atoms[0].box, xij, yij, zij);
const double r2 = xij * xij + yij * yij + zij * zij;
if (r2 > cutoffSquare)
continue;

const double r = sqrt(r2);
const double peTerm = (r > balloonRadius) ? (1.0 / B) * (r - balloonRadius) * (r - balloonRadius) : 0.0;
const double fijTerm = (r > balloonRadius) ? -B * (r - balloonRadius) / r : 0.0;

const double rinv = 1.0 / r;
const double f_ij = e24s6 * pow(rinv, 8) - e48s12 * pow(rinv, 14) + fijTerm;
atoms[i].pe += e4s12 * pow(rinv, 12) - e4s6 * pow(rinv, 6) + peTerm;
atoms[j].pe += e4s12 * pow(rinv, 12) - e4s6 * pow(rinv, 6) + peTerm;
atoms[i].fx += f_ij * xij;
atoms[j].fx -= f_ij * xij;
atoms[i].fy += f_ij * yij;
atoms[j].fy -= f_ij * yij;
atoms[i].fz += f_ij * zij;
atoms[j].fz -= f_ij * zij;
}
}
}

void integrate(const bool isStepOne, const double timeStep, std::vector<Atom>& atoms)
{
const double timeStepHalf = timeStep * 0.5;

for (Atom& atom : atoms)
{
const double ax = atom.fx / atom.mass;
const double ay = atom.fy / atom.mass;
const double az = atom.fz / atom.mass;
if (isStepOne)
{
atom.vx += ax * timeStepHalf;
atom.vy += ay * timeStepHalf;
atom.vz += az * timeStepHalf;
atom.x += atom.vx * timeStep;
atom.y += atom.vy * timeStep;
atom.z += atom.vz * timeStep;
}
else
{
atom.vx += ax * timeStep;
atom.vy += ay * timeStep;
atom.vz += az * timeStep;
}
}
}

std::vector<std::string> getTokens(std::ifstream& input)
{
std::string line;
std::getline(input, line);
std::istringstream iss(line);
std::vector<std::string> tokens{
std::istream_iterator<std::string>{iss},
std::istream_iterator<std::string>{}};

}

int getInt(std::string& token)
{
int value = 0;
try {
value = std::stoi(token);
} catch (const std::exception& e) {
std::cout << "Standard exception:" << e.what() << std::endl;
exit(1);
}
return value;
}

double getDouble(std::string& token)
{
float value = 0;
try {
value = std::stod(token);
} catch (const std::exception& e) {
std::cout << "Standard exception:" << e.what() << std::endl;
exit(1);
}
return value;
}

{
Configuration config;
std::ifstream input("config.ini");
if (!input.is_open()) {
std::cout << "Failed to open config.ini." << std::endl;
exit(1);
}

while (input.peek() != EOF) {
std::vector<std::string> tokens = getTokens(input);
if (tokens.size() > 0) {
if (tokens[0] == "diagnostic_level") {
config.diagnosticLevel = getInt
config.diagnosticLevel = getInt(tokens[1]);
std::cout << "Diagnostic level = " << config.diagnosticLevel << std::endl;
} else if (tokens[0] == "num_steps_thermalization") {
config.numStepsThermalization = getInt(tokens[1]);
std::cout << "Number of steps in thermalization phase = " << config.numStepsThermalization << std::endl;
} else if (tokens[0] == "num_steps_heating") {
config.numStepsHeating = getInt(tokens[1]);
std::cout << "Number of steps in heating phase = " << config.numStepsHeating << std::endl;
} else if (tokens[0] == "num_steps_cooling") {
config.numStepsCooling = getInt(tokens[1]);
std::cout << "Number of steps in cooling phase = " << config.numStepsCooling << std::endl;
} else if (tokens[0] == "balloon_radius") {
} else if (tokens[0] == "energy_output_frequency") {
config.energyOutputFrequency = getInt(tokens[1]);
std::cout << "Energy output frequency = " << config.energyOutputFrequency << std::endl;
} else if (tokens[0] == "coordinate_output_frequency") {
config.coordinateOutputFrequency = getInt(tokens[1]);
std::cout << "Coordinate output frequency = " << config.coordinateOutputFrequency << std::endl;
} else if (tokens[0][0] != '#') {
std::cout << tokens[0] << " is not a valid keyword." << std::endl;
exit(1);
}
}
}

input.close();
return config;
}

void readRun(int& numSteps, double& timeStep, double& temperature)
{
std::ifstream input("run.in");
if (!input.is_open()) {
std::cout << "Failed to open run.in." << std::endl;
exit(1);
}

while (input.peek() != EOF) {
std::vector<std::string> tokens = getTokens(input);
if (tokens.size() > 0) {
if (tokens[0] == "time_step") {
timeStep = getDouble(tokens[1]);
if (timeStep < 0) {
std::cout << "timeStep should >= 0." << std::endl;
exit(1);
}
std::cout << "timeStep = " << timeStep << " fs." << std::endl;
} else if (tokens[0] == "run") {
numSteps = getInt(tokens[1]);
if (numSteps < 1) {
std::cout << "numSteps should >= 1." << std::endl;
exit(1);
}
std::cout << "numSteps = " << numSteps << std::endl;
} else if (tokens[0] == "velocity") {
temperature = getDouble(tokens[1]);
if (temperature < 0) {
std::cout << "temperature >= 0." << std::endl;
exit(1);
}
std::cout << "temperature = " << temperature << " K." << std::endl;
} else if (tokens[0][0] != '#') {
std::cout << tokens[0] << " is not a valid keyword." << std::endl;
exit(1);
}
}
}

input.close();
}

{
std::ifstream input("xyz.in");
if (!input.is_open()) {
std::cout << "Failed to open xyz.in." << std::endl;
exit(1);
}

std::vector<std::string> tokens = getTokens(input);

if (tokens.size() != 1) {
std::cout << "The first line of xyz.in should have one item." << std::endl;
exit(1);
}
int numCells = getInt(tokens[0]);
int totalAtoms = numCells * 4; // Each elementary cell has 4 atoms
std::cout << "Number of atoms = " << totalAtoms << std::endl;

atoms.resize(totalAtoms);

tokens = getTokens(input);
if (tokens.size() != 3) {
std::cout << "The second line of xyz.in should have 3 items." << std::endl;
exit(1);
}
std::cout << "box length = ";
for (int d = 0; d < 3; ++d) {
atoms[0].box[d] = getDouble(tokens[d]);
atoms[0].box[d + 3] = atoms[0].box[d] * 0.5;
std::cout << atoms[0].box[d] << " ";
}
std::cout << std::endl;

// Generate the atom positions in the grid of elementary cells
int atomIndex = 0;
for (int i = 0; i < numCells; ++i) {
for (int j = 0; j < numCells; ++j) {
for (int k = 0; k < numCells; ++k) {
double x = atoms[0].box[0] * (i + 0.5);
double y = atoms[0].box[1] * (j + 0.5);
double z = atoms[0].box[2] * (k + 0.5);
for (int l = 0; l < 4; ++l) {
atoms[atomIndex].x = x;
atoms[atomIndex].y = y;
atoms[atomIndex].z = z;
atoms[atomIndex].mass = 39.948; // Mass of Argon atom
atomIndex++;
// Update coordinates based on the arrangement of atoms in the elementary cell
if (l == 0)
y += atoms[0].box[1];
else if (l == 1)
z += atoms[0].box[2];
else if (l == 2)
x += atoms[0].box[0];
else if (l == 3) {
y -= atoms[0].box[1];
z += atoms[0].box[2];
}
}
}
}
}
}

int main(int argc, char** argv)
{
int numSteps = config.numStepsThermalization + config.numStepsHeating + config.numStepsCooling;
double temperature;
double timeStep;

timeStep /= TIME_UNIT_CONVERSION;

std::vector<Atom> atoms;
initializeVelocity(temperature, atoms);

const int heatingSteps = config.numStepsHeating;
const double initialTemperature = 20.0;
const double finalTemperature = 120.0;
const double deltaT = finalTemperature - initialTemperature;
const double energyChangePerAtom = 3.0 * K_B * deltaT / atoms.size();

const clock_t tStart = clock();
std::ofstream ofile("thermo.out");
ofile << std::fixed << std::setprecision(16);

std::ofstream trajFile("argon_trajektoria.xyz");
trajFile << std::fixed << std::setprecision(6);

double totalEnergy = 0.0;
double prevTotalEnergy = 0.0;
int stabilityCount = 0;

for (int step = 0; step < numSteps; ++step) {
integrate(true, timeStep, atoms);
computeForce(atoms);
integrate(false, timeStep, atoms);
if (step % Ns == 0) {
const double kineticEnergy = computeKineticEnergy(atoms);
const double T = kineticEnergy / (1.5 * K_B * atoms.size());
ofile << T << " " << kineticEnergy << " " << atoms[0].pe << std::endl;
}
if (step == config.numStepsThermalization) {
scaleVelocity(initialTemperature, atoms);
}
if (step >= config.numStepsThermalization && step < config.numStepsThermalization + heatingSteps) {
scaleVelocity(temperature + energyChangePerAtom, atoms);
}
if (step % config.energyOutputFrequency == 0) {
std::ofstream energyFile("argon_energia.csv", std::ios::app); // Open the file in append mode
if (energyFile.is_open()) {
const double kineticEnergy = computeKineticEnergy(atoms);
const double potentialEnergy = atoms[0].pe;
const double repulsiveTerm = e4s12 * pow(rinv, 12);
const double attractiveTerm = -e4s6 * pow(rinv, 6);
const double balloonTerm = (1.0 / B) * (r - config.balloonRadius) * (r - config.balloonRadius);
const double totalEnergy = kineticEnergy + potentialEnergy;

energyFile << std::scientific << std::setprecision(config.diagnosticLevel)
<< step << "," << totalEnergy << "," << potentialEnergy << ","
<< repulsiveTerm << "," << attractiveTerm << ","
<< balloonTerm << "," << kineticEnergy << std::endl;
energyFile.close();
} else {
std::cout << "Failed to open argon_energia.csv for writing." << std::endl;
}
}
if (step % config.coordinateOutputFrequency == 0) {
trajFile << atoms.size() << std::endl;
trajFile << "Step: " << step << std::endl;
for (const Atom& atom : atoms) {
trajFile << "ATOM " << atom.number << " " << std::scientific << std::setprecision(config.diagnosticLevel)
<< atom.x << " " << atom.y << " " << atom.z << std::endl;
}
}

totalEnergy = computeKineticEnergy(atoms) + atoms[0].pe;
if (step >= config.numStepsThermalization) {
if (std::abs(totalEnergy - prevTotalEnergy) < 1e-8)
stabilityCount++;
else
stabilityCount = 0;
prevTotalEnergy = totalEnergy;
}
if (stabilityCount >= 10000) {
std::cout << "Stability condition met. Total energy is stable for 10000 steps." << std::endl;
break;
}
}
ofile.close();
trajFile.close();
const clock_t tStop = clock();
const float tElapsed = float(tStop - tStart) / CLOCKS_PER_SEC;
std::cout << "Time used = " << tElapsed << " s" << std::endl;

return 0;
}
$$`$$

In the second code the full time stepping is given by three lines in main()

integrate(true, timeStep, atoms);
computeForce(atoms);
integrate(false, timeStep, atoms);

This first call of integrate results in $$v_{i+1/2} = v_i + a_i \Delta t/2 \\ x_{i+1} = x_i + v_{i+1/2} \Delta t$$ The call of computeForce updates $$a_{i}$$ to $$a_{i+1}$$ and the second call of integrate updates the velocities $$v_{i+1} = v_{i+1/2} + a_{i+1} \Delta t$$

This looks to be the "kick drift" form of the algorithm, however, the second update to the velocity should use a half time step, $$\Delta t/2$$ instead of $$\Delta t$$.

Thus in the else branch of integrate timeStepHalf should be used instead of timeStep

The first code is wrong if you have multiple particles that interact. Due to the structure of the loop, the interaction forces with particles at the start of the list are computed with the updated positions, that is, you mix positions at time $$t$$ and at time $$t+dt$$. This usually results in an unphysical drift of the whole system, see https://stackoverflow.com/questions/23586195/n-body-gravity-simulation-in-javascript to see this effect in the original jsfiddle.

Thus the necessity to split the process into updates before the force calculation, using the old forces, and updates after the force calculations.

By the way, this is variant of the Verlet principle not leapfrog Verlet, but velocity Verlet. In the leapfrog variant you only keep track of the $$p_i$$ and $$v_{i+1/2}$$, disregarding the $$v_i$$.