0
$\begingroup$

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 . . .
$\endgroup$
9
  • $\begingroup$ 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. $\endgroup$
    – Tyberius
    Sep 10, 2023 at 1:39
  • $\begingroup$ the main thing you are doing wrong is writing your own routine for integration rather than using a library $\endgroup$ Sep 10, 2023 at 5:19
  • 2
    $\begingroup$ 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? $\endgroup$
    – Ian Bush
    Sep 10, 2023 at 13:52
  • $\begingroup$ @IanBush, Yeah? OOP works like this. This is not Fortran or Python. $\endgroup$
    – user366312
    Sep 11, 2023 at 3:40
  • 1
    $\begingroup$ Please mark or better avoid cross-postings, stackoverflow.com/questions/77074492/… $\endgroup$ Sep 11, 2023 at 8:38

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.