# Am I doing leapfrog integration correctly?

I wrote this minimal example to examine the Leapfrog integration algorithm.

However, I am not sure it is the correct algorithm, and the listing is giving the correct output.

Is this the Leapfrog algorithm?

NOTE: I am using the Lennard-Jones system.

NOTE-2: I am using C#.

using System;

public class Vec3
{
public double X { get; set; }
public double Y { get; set; }
public double Z { get; set; }

public Vec3(double x, double y, double z)
{
X = x;
Y = y;
Z = z;
}

public static Vec3 operator +(Vec3 a, Vec3 b)
{
return new Vec3(a.X + b.X, a.Y + b.Y, a.Z + b.Z);
}

public static Vec3 operator *(Vec3 vector, double scalar)
{
return new Vec3(scalar * vector.X, scalar * vector.Y, scalar * vector.Z);
}

public static Vec3 operator *(double scalar, Vec3 vector)
{
return new Vec3(scalar * vector.X, scalar * vector.Y, scalar * vector.Z);
}

public static Vec3 Zero
{
get { return new Vec3(0.0, 0.0, 0.0); }
}
}

public class Atom
{
public Vec3 Position { get; set; }
public Vec3 Velocity { get; set; }

public Atom() { }

public Atom(Vec3 position, Vec3 velocity)
{
Position = position;
Velocity = velocity;
}
}

public class ArgonGasSimulation
{
private const double Sigma = 3.4;
private const double Epsilon = 1.0;

#region [unused]
private static double ComputePotential(double distance)
{
double r6 = Math.Pow(distance, 6);
double r12 = Math.Pow(r6, 2);
return 4 * Epsilon * ((Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance) - (Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance));
}
#endregion

private static Vec3 ComputeForce(Vec3 position)
{
double distance = Math.Sqrt(position.X * position.X + position.Y * position.Y + position.Z * position.Z);
double r6 = Math.Pow(distance, 6);
double r12 = Math.Pow(r6, 2);
double magnitude = 48 * Epsilon * ((Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance) - 0.5 * (Sigma / distance) * (Sigma / distance) * (Sigma / distance) * (Sigma / distance)) / distance;
return magnitude * position;
}

private static void Leapfrog(Atom atom, double timeStep)
{
Vec3 force = ComputeForce(atom.Position); // Compute force based on the current position
Vec3 newPosition = atom.Position + atom.Velocity * timeStep
+ 0.5 * force * timeStep * timeStep;

Vec3 newForce = ComputeForce(newPosition);

Vec3 newVelocity = atom.Velocity + (force + newForce) * timeStep * 0.5;

// Update state
atom.Position = newPosition;
atom.Velocity = newVelocity;
}

public static void Main()
{
// Initial state
Vec3 position = new Vec3(1.0, 2.0, 3.0); // initial position
Vec3 velocity = new Vec3(0.0, 0.0, 0.0); // initial velocity

double timeStep = 0.01;

Atom atom = new Atom(position, velocity);

// Leapfrog integration
Leapfrog(atom, timeStep);

// Output new position and velocity
Console.WriteLine("New Position: ({0}, {1}, {2})", atom.Position.X, atom.Position.Y, atom.Position.Z);
Console.WriteLine("New Velocity: ({0}, {1}, {2})", atom.Velocity.X, atom.Velocity.Y, atom.Velocity.Z);
}
}


Output:

New Position: (1.00014244382518, 2.00028488765035, 3.00042733147553)
New Velocity: (0.028470372373851, 0.0569407447477019, 0.0854111171215529)
Press any key to continue . . .

• At the end of main, you are printing the original position and velocity, rather than the new position and velocity that should be stored in the atom, e.g. atom.position.X. I didn't verify the logic otherwise, but just adding some print statements in Leapfrog shows that newPosition and newVelocity are different from the initial. Sep 10, 2023 at 1:39
• the main thing you are doing wrong is writing your own routine for integration rather than using a library Sep 10, 2023 at 5:19
• Really? For something as simple as Verlet Leapfrog, which should require no more than 2 or 3 lines of code? Out of curiosity which libraries offer an implementation of the Verlet algorithm or its variants? Sep 10, 2023 at 13:52
• @IanBush, Yeah? OOP works like this. This is not Fortran or Python. Sep 11, 2023 at 3:40
• Please mark or better avoid cross-postings, stackoverflow.com/questions/77074492/… Sep 11, 2023 at 8:38