Dancing particles

In this chapter we focus on developing the particle dynamics code. First we install required packages:

In [ ]:
import Pkg; Pkg.add(["Plots", "Printf", "Zygote"])
using LinearAlgebra

Plotting in Julia

In many day-to-day tasks there is not yet the one canonical way to do it in Julia. This also concerns plotting and visualisation of data. Right now quite a few different packages with their own APIs and feature sets exist. To give you some examples:

  • PyPlot (Plotting via python and matplotlib)
  • PGFPlots (Plotting via TeX, TikZ and pgf)
  • Gladify (Pure Julia plotting implementation, widespread)
  • PlotlyJS (Ploting via Javascript and Plotly)
  • GR (Plotting based on the GR framework, pretty fast)
  • Makie (Extremely feature-rich; can plot on the GPU directly)

What we will use in this course is Plots. I like it because it is a metapackage unifying some of the above under a common interface. Flipping the backend is just a single function call, so you can change between them as needed. Also it will probably not disappear any time soon ;). We load Plots first and explicitly switch to the default GR backend.

In [ ]:
using Plots
Plots.gr()

Plotting is then done with the plot function. Additional data series can be added using plot! and attributes can be changed via keywords arguments or standalone functions like xaxis!, title! and so on.

In [ ]:
Vho(r, a=0.5)  = (r - a)^2   # Shifted harmonic oscillator
dVho(r, a=0.5) = 2r -2a      # Derivative

r = collect(-5:0.1:5)        # Create series of x-values
p = plot(r, Vho.(r), label="Vho a=0.5")
plot!(p, r, dVho.(r), label="∂Vho")
xaxis!(p, "relative radial distance")

Of course scatter plots can be made as well ...

In [ ]:
scatter(r, Vho.(r), label="Vho", marker=:cross)
scatter!(r, dVho.(r), label="∂Vho", marker=:circle)
xaxis!("relative radial distance")

... and in fact there are plenty of things you can do to style your plot further, see attributes and layouts for some ideas if you are curious.

Taking derivatives without headaches

Recall that for our dynamics with a given potential, say $$ V\left( \begin{array}{c} x_1 \\ x_2 \end{array}\right) = x_1^2 + x_2^2 $$ we needed the acceleration $$ \vec{A}_V = - \nabla V. $$

For the harmonic oscillator this is kind of ok, but for more complicated potentials $V$, taking the derivative is error prone and most importantly quite boring. The solution is the great package Zygote, which is capable of adjoint-mode automatic differentiation. Without going into details, this means that the Julia compiler takes the derivatives for us and the result is (usually almost) as fast as if we did it ourselves. Let's try this in 1D on Vho as defined above:

In [ ]:
using Zygote

r = collect(-5:0.1:5)
p = plot(r, Vho.(r), label="Vho a=0.5")
xaxis!(p, "relative radial distance")
plot!(p, r, Vho'.(r), label="∂Vho")   # Notice the dash

With ease we add the second derivative to the plot:

In [ ]:
plot!(p, r, Vho''.(r), label="∂∂Vho")

This is extremely helpful for more complicated potentials, for example the Morse potential, which is a common model for a chemical bond: $$ V_\text{Morse} = D (1 - exp(-\alpha (x-x_0)))^2 - D$$

In [ ]:
Vmorse(r; re=0.7, α=1.3, D=10) = D*(1 - exp(-α * (r - re)))^2 - D
In [ ]:
r = collect(0:0.1:7)
p = plot(r, Vmorse.(r), label="Vmorse", ylim=(-25, 25))
xaxis!(p, "Radial distance")
plot!(p, r, Vmorse'.(r), label="Vmorse'")
plot!(p, r, Vmorse''.(r), label="Vmorse''")

By the way: This works for higher dimensions and more complicated expressions, too, we will use this in a sec.

Inspecting Euler dynamics in 1D

Now we return to our forward-euler implementation, which works for both 1D and higher dimensions:

In [ ]:
euler(A, Δt, xn, vn) = xn + vn * Δt, vn + A(xn) * Δt

where $A_V$ (or in the code $A$) was the negative gradient of a potential $V$. With plotting at hand, let us actually see the dynamics. We define a plot function to plot the potential and animate the particle over time:

In [ ]:
using Printf

# The arguments after the ; are so-called keyword arguments, they can be omitted
# and then the default value after the = is used
function plot_dynamics_line(V, Δt, n_steps; x0=0, v0=randn(), integrator=euler,
                            xlim=(-5, 5), ylim=(0, 10))
    t, x, v = 0, x0, v0  # Initialise state
    
    # Compute potential values (for plotting)
    xrange = xlim[1]:0.1:xlim[2]
    Vrange = V.(xrange)
    
    # @gif is a macro to "record" an animation of the dynamics,
    # each loop iteration gives a frame
    @gif for i in 1:n_steps
        # Propagate dynamics (notice the derivative)
        x, v = integrator(x -> -V'(x), Δt, x, v)
        t += Δt

        # Plot potential
        p = plot(xrange, Vrange, label="Potential", xlim=xlim, ylim=ylim)
        
        # Plot the particle and its velocity (as an arrow)
        timestr = @sprintf "%.2f" t   # Format time as a nice string
        scatter!(p, [x], [V(x)], label="Particle at t=$timestr")
        plot!(p, [(x, V(x)), (x + 0.5v, V(x))], arrow=1.0,
              label="particle velocity / 2")
    end
end
In [ ]:
# Now let's actually see it ....
plot_dynamics_line(Vho, 0.1, 200)

Wow ... that's strange. Our plot points at the well-known problem that energy is not conserved for a simple forward Euler scheme. We can also diagnose this using a phase-space plot, where the time evolution of particle position $x$ and particle velocity $v$ is shown. This should be a closed curve if energy is conserved, so let's see ...

In [ ]:
function plot_phase_space(A, Δt, n_steps; x0=0, v0=randn(), integrator=euler,
                          xlim=(-7, 7), ylim=(-7, 7))
    x, v = x0, v0
    p = plot([x], [v], xlim=xlim, ylim=ylim, label="")  # Plot first position
    xaxis!(p, "position")
    yaxis!(p, "velocity")
    @gif for N in 1:n_steps
        x, v = integrator(A, Δt, x, v)
        push!(p, x, v)  # Add new positions to the plot ...
    end
end
In [ ]:
plot_phase_space(x -> -Vho'(x), 0.1, 200)

Exercise:

A much more stable integrator than the euler we used so far is the verlocity Verlet:

$$\left\{\begin{array}{l} x^{(n+1)} = x^{(n)} + v^{(n)} \Delta t + \frac{A_V(x^{(n)})}{2} \Delta t^2\\ v^{(n+1)} = v^{(n)} + \frac{A_V(x^{(n))} + A_V(x^{(n+1)})}{2} \Delta t\\ \end{array}\right. $$
  • Program this algorithm, taking care that it supports multi-dimensional positions and velocities as well. In practice we would like to avoid recomputing $A_V(x)$ as much as possible, since this is usually the expensive step of the dynamics. For our purposes there is no need to keep an eye on that.
  • How does the previous dynamics look like in this example. Does this algorithm conserve energy (phase-space plot)?
  • Also look at the Morse potential
In [ ]:
# solution
In [ ]:
# Some example plots and parameters:
# plot_dynamics_line(Vho, 0.1, 200, integrator=verlet, ylim=(0, 2.5), xlim=(-1, 3))
# plot_phase_space(x -> -Vho'(x), 0.1, 200, integrator=verlet, xlim=(-2, 2), ylim=(-2, 2))
# plot_dynamics_line(Vmorse, 0.03, 200, integrator=verlet, xlim=(-1, 4), ylim=(-10, 5), x0=2.0)

Multiple particles

Now we want to simulate multiple identical partices in 2D. For a system of $N$ particles in 2D, we define the particle positions as the matrix $$ \textbf{x} = \left( \begin{array}{cccc} x_{1x} & x_{2x} & \cdots & x_{Nx} \\ x_{1y} & x_{2y} & \cdots & x_{Ny} \end{array} \right) = \left( \vec{x}_1 \vec{x}_2 \cdots \vec{x}_N \right). $$ with the individual particle vectors as columns. We assume our particles interact via the idential pair potential $V_\text{pair}(r)$ depending only on particle distance $r$. The total potential is therefore $$ V_\text{tot}(\textbf{x}) = \sum_{i=1}^N \sum_{j=i+1}^N V_\text{pair}(\| \vec{x}_i - \vec{x}_j \|). $$

Exercise

Program the total potential function for a matrix $\textbf{x}$. A useful function is norm from the LinearAlgebra package.

In [ ]:
using LinearAlgebra
# You're code here

Now we want to plot the dynamics in a plane. In the following function the acceleration is computed using automatically generated derivatives.

In [ ]:
function plot_dynamics_plane(Vpair, Δt, n_steps; x0=randn(2, 2), v0=zero(x0), integrator=verlet,
                             xlim=(-3, 3), ylim=(-3, 3))
    # Total acceleration via automatic differentiation
    V(x) = Vtot(Vpair, x)
    Atot(x) = -V'(x)
    
    t, x, v = 0, x0, v0  # Initialise state
    @gif for i in 1:n_steps
        # Propagate dynamics
        x, v = integrator(Atot, Δt, x, v)
        t += Δt
        timestr = @sprintf "%.2f" t   # Format time
        
        # Plot the particles and their velocities
        p = scatter(x[1, :], x[2, :], xlim=xlim, ylim=ylim, label="Particles at t=$timestr")
        label = "Velocity / 4"
        for (xi, vi) in zip(eachcol(x), eachcol(v))
            plot!(p, [Tuple(xi), Tuple(xi + 0.25vi)], arrow=1.0, colour=:red, label=label)
            label = ""
        end
    end
end
In [ ]:
Δt = 0.07
n_steps = 300
x0 = [[0.; 0.] [1.; 0.] [-1.; 0.0]]
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0)

A few nice examples

In [ ]:
Δt = 0.07
n_steps = 300
x0 = [[0.; 0.] [1.; 0.] [-1.; 0.15]]
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0)
In [ ]:
Δt = 0.05
n_steps = 300
x0 = [[0.; 1.] [1.; 0.] [-1.; 0] [0; -1.2]]
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0, ylim=(-10, 10))
In [ ]:
Δt = 0.05
n_steps = 300
x0 = 4randn(2, 10)
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0, xlim=(-10, 10), ylim=(-10, 10))