2
$\begingroup$

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;
  double balloonRadius;
  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;
  const double balloonRadius = 10.0; // Adjust this value as needed
  const double B = 1.0 / balloonRadius / balloonRadius;

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

  return tokens;
}

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 readConfig()
{
  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") {
        config.balloonRadius = getDouble(tokens[1]);
        std::cout << "Balloon radius = " << config.balloonRadius << std::endl;
      } 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();
}

void readXyz(std::vector<Atom>& atoms)
{
  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)
{
  Configuration config = readConfig();
  int numSteps = config.numStepsThermalization + config.numStepsHeating + config.numStepsCooling;
  double temperature;
  double timeStep;
  readRun(numSteps, timeStep, temperature);

  timeStep /= TIME_UNIT_CONVERSION;

  std::vector<Atom> atoms;
  readXyz(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;
}
```
$\endgroup$
1

2 Answers 2

6
$\begingroup$

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

$\endgroup$
5
$\begingroup$

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

$\endgroup$

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.