How I used DynamO particle simulator

On my site, I host data on the dynamics of hard-sphere fluids obtained using the DynamO particle simulator. Below, I describe how I used DynamO to generate the data.

1. Simulation parameters
   1. Velocity autocorrelation timescale
   2. Free flight timescales
   3. Sound propagation timescale
   4. Optimal values for $t_\text{run}$ and $N$
2. DynamO simulation setup
   1. Initial configuration file
   2. Equilibrating the system
   3. Running the simulation
3. Outputting particle trajectories
4. Converting DynamO output to HDF5

Simulation parameters

At equilibrium, the monodisperse hard-sphere fluid is described by only two dimensionless parameters: the number of particles $N$ and the reduced number density $\rho^* = \sigma^3 n$. Here $\sigma$ is the diameter of a hard sphere, $n = N/V$ is the particle number density, $V$ is the system volume. We are interested in the thermodynamic limit ($N \to \infty$, finite $\rho^*$), so that $\rho^*$ becomes the only relevant thermodynamic parameter. To study the thermodynamic limit, we conduct simulations of the hard-sphere fluid using periodic boundary conditions and a relatively large number of particles $N$.

To run the simulation, in addition to $\rho^*$, we need to specify time and length scales: the simulation duration $t_\text{run}$ and the size of the domain (which corresponds to $N$). Both need to be large enough to capture temporal and spatial decorrelation of particle velocities. Further, to collect sufficient statistics on free flight times (or, equivalently, collision frequencies), the simulation has to be long enough for each particle to experience several collisions.

Below I discuss the characteristic scales of the problem and suggest values of $t_\text{run}$ and $N$ sufficient to resolve the statistics of interest. Because I am primarily interested in the gas-like rather than the liquid-like state of the hard-sphere fluid, all computations apply to a sufficiently dilute gas.

All the quantities are studied in reduced units, with $\sigma$ being the unit of length and $\sigma \sqrt{m / (k_\text{B} T)}$ being the unit of time. Here $m$ is the mass of a hard sphere, $T$ is the temperature and $k_\mathrm{B}$ is the Boltzmann constant. The unit of velocity is then $\sqrt{k_\text{B} T / m}$.

Velocity autocorrelation timescale

First, we discuss the velocity autocorrelation time $\tau_u$, defined as

The simulation duration $t_\text{run}$ shall satisfy $t_\text{run} \gg \tau_u$. Via a Green–Kubo relation $\tau_u$ is directly related to the self-diffusion coefficient $D$; in reduced units, $\tau_u^* = D^*$:

For a dilute hard-sphere gas:

In dense fluids, $D^*$ and thus $\tau_u^*$ are smaller than the dilute-gas estimation above.

Free flight timescales

Second, we discuss the free flight times of particles, based on the article Collisional statistics of the hard-sphere gas by P. Visco, F. van Wijland, and E. Trizac (Phys. Rev. E 2008). The mean free flight time $\langle t_\text{f} \rangle$ in a dilute hard-sphere gas is given by (needs to be checked)

or, in non-dmensional units,

For small free flight times, their probability density function (PDF) $f(t_\text{t})$ is close to the Poisson one,

However, for $t_\text{f} \gg \langle t_\text{f} \rangle$ the PDF deviates from the Poisson one, its tail is given by

To resolve statistics on the free flight times we want each molecule to collide several times during the simulation. For this we require the expected number $N_\text{f}$ of particles without collisions during the simulation to be much less than $1$:

Sound propagation timescale

Third, we discuss the time $t_c$ it takes for a sound wave to propagate through half the length of the simulation box, $L/2 = \frac{1}{2}\sqrt[3] V$. To avoid the influence of the periodic boundary conditions on the gathered time statistics, we should take $t_\text{run} < t_c$. The speed of sound $c$ in a dilute hard-sphere gas is

Given that $t_c = L/(2c)$, we explicitly calculate $t_c$ to be

Optimal values for $t_\text{run}$ and $N$

Let us calculate $N$ from the requirement $t_\text{run} = t_c$, so that

For this value of $N$, we compute $N_\text{f}$ and plot both of them as functions of $t_\text{run}^*$ for various values of $\rho^*$:

We observe that the number of particles $N$ increases exponentially as we enforce higher accuracy, ensuring each particle experiences a collision ($N_\text{f} \ll 1$) while minimizing the effect of periodic boundaries ($t_\text{run} = t_c$). Moreover, this effect is exacerbated in dilute systems, as $N$ increases exponentially with decreasing $\rho^*$.

We set our goal to accurately simulate systems with $\rho^*$ as low as $0.01$. We require $N_\text{f} < 10^{-2}$ and $t_\text{run} = t_c$. At $\rho^* = 0.01$,following the previous calculations, we get $N > 7.1 \times 10^6$ and $t_\text{run} > 16.3 \tau_u$.

Thus for all systems with $\rho^* \geqslant 0.01$, we use the same fixed value of $N$

and the value of $t_\text{run}$ tied to $\tau_u$ as

DynamO simulation setup

The DynamO simulation proceeds in three main steps:

• Create an initial configuration file.

• Equilibrate the system by running it for some time, to obtain an equilibrated state.

• Simulate the system from the equilibrated state for $t_\text{run}$. During this step DynamO also computes various statistical quantities.

I describe the three steps using a system at $\rho^* = 0.1$. For subsequent processing and conversion to HDF5, it is assumed that the commands below are executed in a folder named 1e-1, ensuring that the generated data is written to this directory.

Initial configuration file

First, I create the initial configuration file config.start.xml, in which particles are packed into a grid of face-centered cubic (FCC) unit cells (4 particles per cell), with the unit cells arranged in a cubic grid of side length $128$. This results a total of $N = 4 \times 128^3$ particles.

The initial configuration file for $\rho^* = 0.1$ and $N = 4 \times 128^3$ is created with the command:

dynamod --pack-mode 0 \
        --NCells 128 \
        --density 0.1 \
        --i1 0 \
        --out-config-file config.start.xml

Equilibrating the system

Once the initial configuration file was created, we equilibrate the system. For the chosen density, $\rho^* = 0.1$, we have $t_\text{run}^* = 36.0$. We equilibrate for $4t_\text{run} = 144$. Thus, the DynamO command used for equilibration is:

dynarun --config-file config.start.xml \
        --sim-end-time 144 \
        --out-data-file output.xml \
        --out-config-file config.equilibrated.xml

Running the simulation

Now we use the equilibrated state to launch the main simulation for $t_\text{run}^* = 36.0$. To save particle trajectories we launch dynarun with the Trajectory output plugin:

dynarun --config-file config.equilibrated.xml \
        --sim-end-time 36 \
        --load-plugin Trajectory \
        --load-plugin MSD \
        --unwrapped \
        --out-data-file output.xml \
        --out-config-file config.end.xml

Once the simulation is finished, you can find the detailed particle trajectories in the trajectory.out file. However, the Trajectory plugin, at the time of writing, is more of a diagnostic rather than investigatory tool. For it to output the desired data, I had to modify it, as I describe below.

Outputting particle trajectories

I output trajectories data as a CSV (technically, space-separated values) file with the following columns

To achieve this output from DynamO, I had to modify the Trajectory output plugin located in DynamO/src/dynamo/dynamo/outputplugins/trajectory.cpp. Specifically, I modified two methods, OPTrajectory::initialise and OPTrajectory::eventUpdate as follows:

void OPTrajectory::initialise() {
  if (logfile.is_open())
    logfile.close();

  logfile.open("trajectory.out", std::ios::out | std::ios::trunc);
  logfile.precision(std::numeric_limits<double>::max_digits10);
  logfile.setf(std::ios::scientific);

  logfile << "event "
          << "t "
          << "dt "
          << "p1 "
          << "p2 "
          << "x1 "
          << "y1 "
          << "z1 "
          << "x2 "
          << "y2 "
          << "z2 "
          << "u1 "
          << "v1 "
          << "w1 "
          << "u2 "
          << "v2 "
          << "w2\n";
}

void OPTrajectory::eventUpdate(const Event &eevent, const NEventData &SDat) {
  if (SDat.L2partChanges.empty()) { 
    return; // skips single-particle events (VIRTUAL for our case)
  }

  logfile << std::setw(12) << std::setfill('0') << Sim->eventCount // event
          << " " << Sim->systemTime / Sim->units.unitTime()        // t
          << " " << eevent._dt / Sim->units.unitTime();            // dt


  for (const PairEventData &pData : SDat.L2partChanges) {
    const size_t id1 = std::min(pData.particle1_.getParticleID(),
                                pData.particle2_.getParticleID());
    const size_t id2 = std::max(pData.particle1_.getParticleID(),
                                pData.particle2_.getParticleID());

    Vector r1 = Sim->particles[id1].getPosition() / Sim->units.unitLength();
    Vector r2 = Sim->particles[id2].getPosition() / Sim->units.unitLength();

    Vector v1 = Sim->particles[id1].getVelocity() / Sim->units.unitVelocity();
    Vector v2 = Sim->particles[id2].getVelocity() / Sim->units.unitVelocity();

    logfile << " " << id1    // p1
            << " " << id2    // p2
            << " " << r1[0]  // x1
            << " " << r1[1]  // y1
            << " " << r1[2]  // z1
            << " " << r2[0]  // x2
            << " " << r2[1]  // y2
            << " " << r2[2]  // z2
            << " " << v1[0]  // u1
            << " " << v1[1]  // v1
            << " " << v1[2]  // w1
            << " " << v2[0]  // u2
            << " " << v2[1]  // v2
            << " " << v2[2]; // w2
  }
  logfile << "\n";
}

Converting DynamO output to HDF5

Once the simulation is finished I convert the trajectory data to the HDF5 format structured as described preivously.

Assuming that the data was written into the folder 1e-1, the following code creates hsf_1e-1.h5 with processed data:

using HDF5
using CSV
using DataFrames
using EzXML
using Statistics
using LinearAlgebra
using StaticArrays
using ProgressMeter


function get_config_x_u(config_file, N)
  x = Vector{SVector{3, Float64}}(undef, N)
  u = Vector{SVector{3, Float64}}(undef, N)

  config_xml = EzXML.readxml(config_file)
  particles = findall("//Pt", config_xml)
  for pt in particles
    i = parse(Int64, pt["ID"]) + 1
    x_node = findfirst("P", pt)
    u_node = findfirst("V", pt)
    x[i] = SA[parse(Float64, x_node["x"]),
              parse(Float64, x_node["y"]),
              parse(Float64, x_node["z"])]
    u[i] = SA[parse(Float64, u_node["x"]),
              parse(Float64, u_node["y"]),
              parse(Float64, u_node["z"])]
  end

  return x, u
end


function get_parameters(output_file)
  output_xml = EzXML.readxml(output_file)
  ρ_node = findfirst("//Density", output_xml)
  N_node = findfirst("//ParticleCount", output_xml)
  duration_node = findfirst("//Duration", output_xml)
  size_node = findfirst("//PrimaryImageSimulationSize", output_xml)
  
  ρ = parse(Float64, ρ_node["val"])
  N = parse(Int64, N_node["val"])
  L = parse(Float64, size_node["x"])
  t_end = parse(Float64, duration_node["Time"])
  
  return ρ, N, L, t_end
end


function check_consistency(t_trajectory, β_trajectory, x_trajectory, u_trajectory)
  ################################################################################
  println("\nChecking position-velocity consistency...")
  
  x_from_u = [vcat( [x[1]], [x[1]] .+ cumsum(u[1:end-1] .* diff(t)) )
    for (t, x, u) in zip(t_trajectory, x_trajectory, u_trajectory)]
  
  Δx_max = 0.0
  for (x₁, x₂) in collect(zip(x_trajectory, x_from_u))
    Δx = maximum(norm.(x₁ .- x₂))
    Δx_max = Δx > Δx_max ? Δx : Δx_max
  end
  println("maximum Δx = ", Δx_max, " (should be ≈ +0)")

  ################################################################################
  println("\nChecking velocity-position consistency...")
  
  u_from_x = [diff(x) ./ diff(t) for (t, x) in zip(t_trajectory, x_trajectory)]
  for u in u_from_x push!(u, u[end]) end
  
  Δu_max = 0.0
  for (u₁, u₂) in zip(u_trajectory, u_from_x)
    Δu = maximum(norm.(u₁ .- u₂))
    Δu_max = Δu > Δu_max ? Δu : Δu_max
  end
  println("maximum Δu = ", Δu_max, " (should be ≈ +0)\n")
  
end


function process_DynamO_data(folder_name)
  _, N, _, t_end = get_parameters(folder_name * "/output.xml")

  println("Reading initial and final configurations...")
  x_start, u_start = get_config_x_u(folder_name * "/config.equilibrated.xml", N)
  x_end,   u_end   = get_config_x_u(folder_name * "/config.end.xml", N) 

  println("Reading trajectory data...")
  trajectory_df  = CSV.read(folder_name * "/trajectory.out", DataFrame, delim = ' ')

  println("Processing trajectory data...")
  
  t_trajectory = [[0.0] for _ in 1:N]
  x_trajectory = [[x] for x in x_start]
  u_trajectory = [[u] for u in u_start]
  β_trajectory = [Array{Int64, 1}([]) for _ in 1:N]

  @showprogress for r in eachrow(trajectory_df)
    push!(t_trajectory[r["p1"] + 1], r["t"])
    push!(t_trajectory[r["p2"] + 1], r["t"])

    push!(β_trajectory[r["p1"] + 1], r["p2"] + 1)
    push!(β_trajectory[r["p2"] + 1], r["p1"] + 1)
    
    push!(x_trajectory[r["p1"] + 1], SA[r["x1"], r["y1"], r["z1"]])
    push!(x_trajectory[r["p2"] + 1], SA[r["x2"], r["y2"], r["z2"]])

    push!(u_trajectory[r["p1"] + 1], SA[r["u1"], r["v1"], r["w1"]])
    push!(u_trajectory[r["p2"] + 1], SA[r["u2"], r["v2"], r["w2"]])
  end

  for i in 1:N
    push!(t_trajectory[i], t_end)
    push!(x_trajectory[i], x_end[i])
    push!(u_trajectory[i], u_end[i])
  end

  check_consistency(t_trajectory, β_trajectory, x_trajectory, u_trajectory)
  
  return t_trajectory, β_trajectory, x_trajectory, u_trajectory
end


function process_and_write_DynamO_data(folder_name)
  ρ, N, L, t_end = get_parameters(folder_name * "/output.xml")
  
  t_trajectory, β_trajectory, x_trajectory, _ = process_DynamO_data(folder_name)

  file_name = folder_name * "/hsf_" * folder_name * ".h5"
  println("Writing trajectory data into " * file_name * "...")
  h5open(file_name, "w") do f
    f["t_start"] = 0.0
    f["t_end"]   = t_end
    f["ρ"] = ρ
    f["N"] = N
    f["L"] = L
    f["README"] = "This file contains particle trajectories from the simulation \
    of a hard-sphere fluid at thermal equilibrium. The simulation was performed \
    at reduced density ρ with N particles in a periodic cubic domain of side L, \
    using the DynamO particle simulator; for details, see https://sarnitsky.se/."

    i_to_str(i) = lpad(i, 9, '0')
    
    println("Writing t data...")
    t_data = create_group(f, "t")
    @showprogress for (i, t) in enumerate(t_trajectory)
      t_data[i_to_str(i)] = t
    end

    println("Writing x data...")
    x_data = create_group(f, "x")
    @showprogress for (i, x) in enumerate(x_trajectory)
      x_data[i_to_str(i)] = collect(reinterpret(reshape, Float64, x))
    end

    println("Writing collision data...")
    β_data = create_group(f, "collision_partner")
    @showprogress for (i, β) in enumerate(β_trajectory)
      β_data[i_to_str(i)] = β
    end
    
  end
end


process_and_write_DynamO_data("1e-1")