ArrayStruct(N)

A structure that holds various arrays for storing particle properties and the simulation state. These arrays are typically updated at each time step of the simulation. The type parameter NB refers to the type of the neighbor list structure used (e.g., VoronoiNeighborList).

Fields

• positions::Vector{SVector{2, Float64}}: Current 2D positions of all particles.
• old_positions::Vector{SVector{2, Float64}}: Positions of all particles at the previous time step. Used for some integration schemes or state tracking.
• forces::Vector{SVector{2, Float64}}: Current 2D forces acting on all particles.
• old_forces::Vector{SVector{2, Float64}}: Forces acting on all particles at the previous time step. Used for some integration schemes.
• orientations::Vector{Float64}: Current orientations of all particles, typically represented as an angle.
• old_orientations::Vector{Float64}: Orientations of all particles at the previous time step.
• areas::Vector{Float64}: Current areas of the Voronoi cell corresponding to each particle.
• perimeters::Vector{Float64}: Current perimeters of the Voronoi cell corresponding to each particle.
• random_forces::Vector{Float64}: Stores random numbers or forces, often used for implementing stochastic elements like rotational diffusion.
• neighborlist::NB: A structure (e.g., VoronoiNeighborList) that holds neighbor information for each particle, essential for calculating interactions and Voronoi properties.

DumpInfo(; save=true, filename="dump_...", when_to_save_array=0:1000:1000000, save_r=true, save_F=false, save_Epot=false)

A mutable struct that holds parameters for controlling the dumping (saving) of simulation data to a file. The type parameter A corresponds to the type of when_to_save_array, which is typically an array or range of integers.

Fields

• save::Bool: A boolean flag to enable (true) or disable (false) the saving of simulation data. Defaults to true.
• filename::String: The name of the file where simulation data will be saved. Defaults to a randomly generated HDF5 filename (e.g., "dump_12345.h5").
• when_to_save_array::A: An array or range specifying the simulation steps at which data should be saved. A is the type of this field, typically AbstractArray. Defaults to saving every 1000 steps up to 1,000,000.
• save_r::Bool: A boolean flag indicating whether to save particle positions (true) or not (false). Defaults to true.
• save_F::Bool: A boolean flag indicating whether to save particle forces (true) or not (false). Defaults to false.
• save_Epot::Bool: A boolean flag indicating whether to save the total potential energy of the system (true) or not (false). Defaults to false.

Output()

A mutable struct that stores various simulation output quantities. These quantities are typically updated throughout the simulation and can be used for monitoring progress, logging, or analysis.

Fields

• potential_energy::Float64: The total potential energy of the system at the current state.
• steps_done::Int64: The number of simulation time steps that have been completed.
• N_voronoi_tesselations::Int64: The total number of Voronoi tesselations that have been performed during the simulation.

ParameterStruct(N =100, dt = 0.01, N_steps = 10000, kBT = 1.0, frictionconstant = 1.0,
                periodic_boundary_layer_depth = 1.0, verbose = false, box = SimulationBox(10.0, 10.0),
                particles = VoronoiCells(zeros(Float64, 100), zeros(Float64, 100), zeros(Float64, 100),
                                          zeros(Float64, 100), ones(Float64, 100), ones(Float64, 100)),
                dump_info = DumpInfo(), callback = x -> nothing,
                RNG = Random.MersenneTwister(1234))

A struct that holds all essential parameters for running a simulation. This includes physical parameters, numerical parameters, and settings for I/O and control.

Type parameters:

• P: The type of the particle properties structure (e.g., VoronoiCells), which must be a subtype of Particles.
• CB: The type of the callback function.

Fields

• N::Int: The total number of particles in the simulation.
• dt::Float64: The time step size for the integration algorithm.
• N_steps::Int: The total number of simulation steps to perform.
• kBT::Float64: The thermal energy, given by the product of Boltzmann's constant (kB) and temperature (T).
• frictionconstant::Float64: The friction coefficient, determining the damping of particle motion.
• periodic_boundary_layer_depth::Float64: The depth of the layer around the simulation box used for periodic boundary condition checks and neighbor finding.
• verbose::Bool: A boolean flag to enable (true) or disable (false) verbose output during the simulation.
• box::SimulationBox: A SimulationBox struct defining the dimensions and properties of the simulation domain.
• particles::P: A struct (of type P) holding the specific properties of the particles (e.g., VoronoiCells which includes target areas, perimeters, etc.).
• dump_info::DumpInfo: A DumpInfo struct containing parameters for controlling how simulation data is saved to files.
• callback::CB: A callback function (of type CB) that can be executed at specified intervals during the simulation (e.g., for custom analysis or logging).
• RNG::Random.MersenneTwister: An instance of a random number generator (specifically MersenneTwister) used for stochastic processes in the simulation.

SimulationBox(Lx, Ly)

Defines the simulation domain, typically a rectangular box with periodic boundary conditions.

Fields

• box_sizes::SVector{2, Float64}: A 2D static vector representing the dimensions (Lx, Ly) of the simulation box. Lx is the width and Ly is the height.

VoronoiCells(target_perimeters, target_areas, K_P, K_A, active_force_strengths, rotational_diffusion_constants)

Represents the properties of particles modeled as Voronoi cells. This struct stores parameters related to the mechanical properties (target area and perimeter, and associated spring constants) and active behavior (active force strength and rotational diffusion) of each particle/cell.

Fields

• target_perimeters::Vector{Float64}: Target perimeters for each Voronoi cell. Deviations from these values incur an energy penalty.
• target_areas::Vector{Float64}: Target areas for each Voronoi cell. Deviations from these values incur an energy penalty.
• K_P::Vector{Float64}: Spring constants for perimeter deviations. Determines how strongly a cell resists changes to its perimeter.
• K_A::Vector{Float64}: Spring constants for area deviations. Determines how strongly a cell resists changes to its area.
• active_force_strengths::Vector{Float64}: Magnitudes of the active force for each particle. This force propels the particle in the direction of its orientation.
• rotational_diffusion_constants::Vector{Float64}: Rotational diffusion rates for each particle. This determines how quickly the particle's orientation changes randomly.

NeighborList

A structure to hold the neighbor list and Voronoi edges for each particle. The neighbors field is a vector of vectors, where each inner vector contains the indices of the neighboring particles for each particle. The voronoi_vertices field is a vector of vectors, where each inner vector contains the index of the Voronoi vertices for each particle. The coordinates are accessed by delauney_circumcenters and are stored in counterclockwise order The delauney_circumcenters field is a vector of vectors, where each inner vector contains the coordinates of the Delauney circumcenters for each particle. The Delauney circumcenters are represented as SVector{2, Float64}.

applyperiodicboundaryconditions(position, boxsizes)

Applies periodic boundary conditions to a given position, ensuring it wraps around the simulation box. For a position coordinate x and a box dimension L, the new coordinate x' is x - floor(x/L) * L. This maps x to the interval [0, L). The same logic applies to all dimensions.

Arguments

• position: The original position, typically an SVector or any AbstractVector representing coordinates (e.g., [x, y]).
• box_sizes: The dimensions of the simulation box, typically an SVector or AbstractVector (e.g., [Lx, Ly]).

Returns

• new_position: The position after applying periodic boundary conditions, of the same type as position.

circumcenter(a::SVector{2, T}, b::SVector{2, T}, c::SVector{2, T}) where T

Calculates the circumcenter of a triangle defined by three 2D vertices a, b, and c. The circumcenter is the point where the perpendicular bisectors of the triangle's sides intersect, and it is equidistant from all three vertices. This point is also the center of the triangle's circumcircle.

The calculation is based on a standard formula for the circumcenter's coordinates (ux, uy) derived from the coordinates of the vertices (ax, ay), (bx, by), (cx, cy): D = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) ux = (1/D) * ((ax^2 + ay^2) * (by - cy) + (bx^2 + by^2) * (cy - ay) + (cx^2 + cy^2) * (ay - by)) uy = (1/D) * ((ax^2 + ay^2) * (cx - bx) + (bx^2 + by^2) * (ax - cx) + (cx^2 + cy^2) * (bx - ax)) This function uses norm2 which computes ax^2 + ay^2 (squared length from origin).

Arguments

• a::SVector{2, T}: The 2D coordinates of the first vertex of the triangle. T is its element type.
• b::SVector{2, T}: The 2D coordinates of the second vertex of the triangle. T is its element type.
• c::SVector{2, T}: The 2D coordinates of the third vertex of the triangle. T is its element type.

Returns

• SVector{2, Float64}: The 2D coordinates (ux, uy) of the circumcenter. The coordinates are always returned as Float64, regardless of the input type T.

Notes

• The function relies on norm2(v) to calculate the squared magnitude of vectors a, b, and c (from the origin), which are used in the circumcenter formula.
• If the three points a, b, c are collinear, the denominator D will be zero, leading to Inf or NaN coordinates. This case is not explicitly handled here beyond the standard floating-point behavior.

compute_energy(parameters, arrays, output)

Computes the total potential energy of the system. This energy is typically derived from the deviations of individual Voronoi cell areas and perimeters from their target values, penalized by corresponding spring constants.

This function first calls update_areas! and update_perimeters! to ensure that the current areas and perimeters in arrays are up-to-date before calculating the energy.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct.
• arrays::ArrayStruct: The struct holding simulation arrays. It's passed to update_areas! and update_perimeters!, which use the current areas and perimeters.
• output::Output: The simulation output struct. Passed to update_areas! and update_perimeters!.

Returns

• potential_energy::Float64: The total calculated potential energy of the system. This is the sum of energy contributions from each cell, where each cell's energy is K_A * (area - target_area)^2 + K_P * (perimeter - target_perimeter)^2.

compute_energy_i(i, parameters, arrays, output)

Computes the potential energy contribution of a single Voronoi cell indexed by i. This function calculates the energy based on the deviation of the cell's area and perimeter from their target values, penalized by the corresponding spring constants.

Arguments

• i::Int: The index of the Voronoi cell for which to compute the energy.
• parameters::ParameterStruct: The main simulation parameter struct. Used to access target areas, target perimeters, and spring constants.
• arrays::ArrayStruct: The struct holding simulation arrays. Used to access the current areas and perimeters of the Voronoi cells.
• output::Output: The simulation output struct. Not directly used in this function but included for consistency in function signatures across the module.

Returns

• E::Float64: The computed energy for the Voronoi cell indexed by i. This is calculated as:

KA[i] * (areas[i] - targetareas[i])^2 + KP[i] * (perimeters[i] - targetperimeters[i])^2

compute_forces_GCM!(parameters, arrays, output)

Computes the forces on each particle using a Gaussian Core Model (GCM). The GCM is a soft, repulsive pair potential where the force between two particles i and j is given by F_ij = -2 * r_ij * exp(-r_ij^2), where r_ij is the distance vector between them. This function calculates the net force on each particle by summing these pairwise forces. Periodic boundary conditions are handled using compute_pair_distance_vector.

This function might serve as an alternative or simpler test case for force calculations compared to the more complex SPV model.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct. It provides:
  • parameters.N: The number of particles.
  • parameters.box.box_sizes: Dimensions of the simulation box, used by compute_pair_distance_vector.
• arrays::ArrayStruct: The struct holding simulation arrays.
  • arrays.positions: A vector of SVector{2, Float64} representing the current positions of all particles.
  • arrays.forces: A vector of SVector{2, Float64} where the computed force for each particle will be stored. This array is modified in-place.
• output::Output: The simulation output struct. Not directly used in this function but included for consistency in function signatures.

Notes

• The function modifies arrays.forces in-place with the newly computed forces.
• The interaction is considered only if the distance r is less than 2^(1/6) (a common cutoff for this potential, related to its inflection point).

compute_forces_SPV!(parameters, arrays, output)

Computes the forces on each particle based on the Self-Propelled Voronoi (SPV) model. The forces arise from deviations of cell areas and perimeters from their target values, effectively modeling area and perimeter elasticity.

The function first ensures the Voronoi tessellation is current by calling verify_tesselation. If the tessellation is not valid, it's updated by voronoi_tesselation! or update_delauney_vertices!. Subsequently, it calls update_perimeters! and update_areas! to ensure these geometric properties are current before force calculation.

The force on particle i is calculated as the negative gradient of the total potential energy with respect to its position r_i. This involves complex derivatives of cell areas and perimeters with respect to vertex positions, and then of vertex positions with respect to particle positions.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct. This provides access to:
  • parameters.N: Number of particles.
  • parameters.particles: Contains target_areas, target_perimeters, K_A (area spring constant), K_P (perimeter spring constant).
• arrays::ArrayStruct: The struct holding simulation arrays.
  • arrays.neighborlist: Used to access Voronoi neighbor information, vertex indices, vertex positions per particle, and positions with periodic boundary conditions.
  • arrays.areas: Current areas of Voronoi cells (read after being updated by update_areas!).
  • arrays.perimeters: Current perimeters of Voronoi cells (read after being updated by update_perimeters!).
  • arrays.forces: A vector of SVector{2, Float64} where the computed force for each particle will be stored. This array is modified in-place.
• output::Output: The simulation output struct. Passed to helper functions like voronoi_tesselation!, update_perimeters!, etc.

Notes

• The function modifies arrays.forces in-place with the newly computed forces.
• The calculation involves contributions from the particle's own cell energy (dEi/dri) and from the energy of its neighboring cells (dEj/dri).
• The mathematical details of the force calculation (derivatives of geometry) can be complex and are implemented within the loops.

compute_pair_distance_vector(p1, p2, box_sizes)

Computes the shortest vector (displacement) from particle p1 to particle p2 in a system with periodic boundary conditions (PBC). This is often referred to as the minimum image convention. For a coordinate difference dx and box dimension L, the PBC-aware difference dx' is dx - round(dx/L) * L.

Arguments

• p1: Position of the first particle (e.g., an SVector or AbstractVector like [x1, y1]).
• p2: Position of the second particle (e.g., an SVector or AbstractVector like [x2, y2]).
• box_sizes: Dimensions of the simulation box (e.g., an SVector or AbstractVector like [Lx, Ly]).

Returns

• delta: The displacement vector from p1 to p2 after applying the minimum image convention due to PBC. This will be of the same type as p1 and p2.

do_time_step(parameters, arrays, output)

Advances the simulation by a single time step by calling an appropriate integration scheme.

This function acts as a dispatcher for different integration methods. Currently, it selects:

• do_time_step_Euler_Murayama for the very first step of the simulation (output.steps_done == 0).
• do_time_step_Euler_Heun! for all subsequent steps.

This strategy might be used, for example, to ensure proper initialization of "old" states required by more complex integrators like Heun's method.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct, containing settings like dt, particle properties, etc. Passed through to the chosen integrator function.
• arrays::ArrayStruct: The struct holding simulation arrays (positions, forces, orientations, etc.). These are modified in-place by the chosen integrator function.
• output::Output: The simulation output struct. output.steps_done is used to determine which integration scheme to use. It's also passed through to the chosen integrator function.

Returns

• The function itself doesn't return a value but calls one of the integrator functions (do_time_step_Euler_Murayama or do_time_step_Euler_Heun!), which in turn modify arrays (and potentially output) in-place.

do_time_step_Euler_Heun!(parameters, arrays, output)

Performs a single time step using the Euler-Heun predictor-corrector integration scheme. This method is suitable for stochastic differential equations (SDEs) and offers better stability and accuracy than the simpler Euler-Maruyama method for some systems. It involves predicting a future state and then correcting it using an average of forces/drifts.

The operational steps are:

1. Compute forces at current positions: compute_forces_SPV! is called to get forces F(x(t_n)) based on current arrays.positions (state at t_n).
2. Update orientations: Particle orientations θ(t_n) are updated to θ(t_{n+1}) by adding a stochastic rotational diffusion term. arrays.orientations now stores θ(t_{n+1}). The orientations from the start of the step, θ(t_n), are available in arrays.old_orientations (assuming they were stored from a previous step or initialization).
3. Predictor step:
   • A predicted position x_tilde(t_{n+1}) is calculated using forces F(x(t_n)) and active forces based on orientations θ(t_n) (from arrays.old_orientations).
   • arrays.positions are updated to these predicted positions.
4. Store current forces: The forces F(x(t_n)) (used in the predictor step) are stored in arrays.old_forces.
5. Compute forces at predicted positions: compute_forces_SPV! is called again to calculate forces F(x_tilde(t_{n+1})) based on the predicted arrays.positions. arrays.forces now stores these new forces.
6. Corrector step:
   • The final positions x(t_{n+1}) are calculated using an average of the "old" forces/orientations and the "new" (predicted) forces/orientations. Specifically, it uses 0.5 * [ (F(x(t_n)) + v_active(θ(t_n))) + (F(x_tilde(t_{n+1})) + v_active(θ(t_{n+1}))) ].
   • arrays.positions are updated to these corrected positions.
7. Update "old" states: arrays.old_positions are updated with the corrected arrays.positions, and arrays.old_orientations are updated with arrays.orientations (θ(t_{n+1})) to prepare for the next time step.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct, providing:
  • dt::Float64: The time step size.
  • frictionconstant::Float64: Used to calculate mobility (1 / frictionconstant).
  • particles.active_force_strengths::Vector{Float64}: Magnitudes of the active force.
  • particles.rotational_diffusion_constants::Vector{Float64}: Rotational diffusion rates.
  • N::Int: The total number of particles.
  • box.box_sizes::SVector{2, Float64}: Dimensions for periodic boundary conditions.
• arrays::ArrayStruct: The struct holding simulation arrays. The following are used and modified in-place:
  • positions::Vector{SVector{2, Float64}}: Initially x(t_n), then x_tilde(t_{n+1}), then corrected x(t_{n+1}).
  • forces::Vector{SVector{2, Float64}}: Initially stores F(x(t_n)), then F(x_tilde(t_{n+1})).
  • orientations::Vector{Float64}: Updated to θ(t_{n+1}) after rotational diffusion.
  • old_positions::Vector{SVector{2, Float64}}: Stores x(t_n) at the start, updated to corrected x(t_{n+1}) at the end.
  • old_forces::Vector{SVector{2, Float64}}: Updated to store F(x(t_n)).
  • old_orientations::Vector{Float64}: Stores θ(t_n) at the start, updated to θ(t_{n+1}) at the end.
• output::Output: The simulation output struct. Passed to compute_forces_SPV!.

Notes

• All relevant fields in arrays are modified in-place throughout the function.
• This scheme aims for better accuracy by averaging drifts over the interval.

do_time_step_Euler_Murayama(parameters, arrays, output)

Performs a single time step using the Euler-Maruyama integration scheme, suitable for stochastic differential equations (SDEs). This method updates particle positions and orientations considering deterministic forces, active self-propulsion forces, and stochastic rotational diffusion.

The function first calls compute_forces_SPV! to calculate the current deterministic forces on all particles based on their positions and Voronoi cell properties. Then, it updates orientations by adding a random component scaled by the rotational diffusion constant. Finally, it updates positions based on the sum of deterministic forces and active forces (directed by the orientation at the beginning of the step), scaled by mobility and time step dt. Periodic boundary conditions are applied to the new positions.

At the end of the step, the old_positions, old_forces, and old_orientations arrays are updated to store the state for potential use in subsequent, more complex integration schemes (though this specific integrator is often used as a simpler baseline or for the first step).

Arguments

• parameters::ParameterStruct: The main simulation parameter struct, providing:
  • dt::Float64: The time step size.
  • frictionconstant::Float64: Used to calculate mobility (1 / frictionconstant).
  • particles.active_force_strengths::Vector{Float64}: Magnitudes of the active force for each particle.
  • particles.rotational_diffusion_constants::Vector{Float64}: Rotational diffusion rates for each particle.
  • N::Int: The total number of particles.
  • box.box_sizes::SVector{2, Float64}: Dimensions of the simulation box for periodic boundary conditions.
• arrays::ArrayStruct: The struct holding simulation arrays. The following fields are used and modified:
  • positions::Vector{SVector{2, Float64}}: Updated in-place with new particle positions.
  • forces::Vector{SVector{2, Float64}}: Updated in-place by compute_forces_SPV!.
  • orientations::Vector{Float64}: Updated in-place with new particle orientations.
  • old_positions::Vector{SVector{2, Float64}}: Updated in-place to store the new positions.
  • old_forces::Vector{SVector{2, Float64}}: Updated in-place to store the new forces.
  • old_orientations::Vector{Float64}: Updated in-place to store the new orientations. (Note: the position update uses orientations from before this step's rotational diffusion for the active force direction).
• output::Output: The simulation output struct. Passed to compute_forces_SPV!.

Notes

• All relevant fields in arrays (positions, forces, orientations, oldpositions, oldforces, old_orientations) are modified in-place.

norm2(v::SVector{2, T}) where T

Computes the squared Euclidean norm (also known as squared length or squared magnitude) of a 2D vector v.

For a vector v = [v1, v2], this is calculated as v1^2 + v2^2. This is often used in computations where the actual distance (square root) is not needed, as it avoids a potentially costly square root operation.

Arguments

• v::SVector{2, T}: The 2D input vector. T is its element type (e.g., Float64).

Returns

• ::T: The squared Euclidean norm of the vector v, equal to v[1]^2 + v[2]^2. The return type is the same as the element type of the input vector.

replace_or_push!(array, value, index)

Modifies an array in-place by either replacing an element at a specified index or pushing a new value onto the end of the array.

The behavior is as follows:

• If index is within the current bounds of the array (i.e., 1 <= index <= length(array)), the element at array[index] is replaced with value.
• If index is one greater than the current length of the array (i.e., index == length(array) + 1), the value is pushed onto the end of the array using push!.
• If index is out of these bounds (e.g., less than 1, or greater than length(array) + 1), an ArgumentError is thrown.

Arguments

• array: The array to be modified. This array is changed in-place.
• value: The value to be inserted into the array or pushed onto its end.
• index::Integer: The index at which to replace the existing element, or the position (equal to length(array) + 1) where the new element should be pushed.

Returns

• The function does not explicitly return a value (nothing is implicitly returned if an ArgumentError is not thrown). The primary effect is the modification of the input array.

Throws

• ArgumentError: If index is not within the valid range for replacement or pushing (i.e., if index < 1 or index > length(array) + 1).

run_simulation!(parameters, arrays, output, N_steps)

Runs the main simulation loop for a specified number of time steps.

Before starting the loop, it performs an initial Voronoi tessellation by calling voronoi_tesselation!(parameters, arrays, output) to ensure the system's geometric properties are up-to-date.

The main loop then proceeds as follows for each step:

1. Verbose Output: If parameters.verbose is true and the current step is a multiple of 100, the step number is printed to the console.
2. Callback Invocation: The user-provided parameters.callback function is called, passing parameters, arrays, and output. This allows for custom actions or analysis during the simulation.
3. Time Step Advancement: do_time_step(parameters, arrays, output) is called to advance the simulation state by one time step, updating particle positions, orientations, and forces.
4. Data Dumping: If parameters.dump_info.save is true, it checks if the current step is in parameters.dump_info.when_to_save_array. If so, save_simulation_state!(parameters, arrays, output) is called to save the current simulation state.
5. Step Counting: The local step counter is incremented, and output.steps_done is updated to reflect the completed step.
6. Termination: The loop breaks if step exceeds N_steps.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct. Key fields used:
  • verbose::Bool: Flag to enable/disable verbose output.
  • dump_info::DumpInfo: Contains parameters for saving simulation data, including save (flag), when_to_save_array (steps for saving).
  • callback: A user-defined function called at each step.
• arrays::ArrayStruct: The struct holding simulation arrays (positions, forces, orientations, etc.). This is modified in-place by the functions called within the loop (e.g., do_time_step, voronoi_tesselation!).
• output::Output: The simulation output struct.
  • output.steps_done: Used to initialize the local step counter and is updated at each iteration. It is also passed to and modified by functions called within the loop.
• N_steps::Int: The total number of time steps to run the simulation.

Notes

• The arrays and output structs are modified in-place throughout the simulation by the various functions called.
• The simulation starts from the step number indicated by output.steps_done at the beginning of the function call.

save_simulation_state!(parameters::ParameterStruct, arrays::ArrayStruct, output::Output)

Saves the current state of the simulation to an HDF5 file.

The function uses the filename specified in parameters.dump_info.filename. It only supports files with .h5 or .hdf5 extensions.

File Structure:

If the HDF5 file does not exist, it will be created, and initial simulation parameters will be saved under the /parameters group. These include:

• /parameters/N: Number of particles
• /parameters/dt: Timestep size
• /parameters/kBT: Boltzmann temperature
• /parameters/frictionconstant: Friction constant
• /parameters/box_sizes: Dimensions of the simulation box (SVector)
• /parameters/particles/target_perimeters: Target perimeters for Voronoi cells
• /parameters/particles/target_areas: Target areas for Voronoi cells
• /parameters/particles/K_P: Perimeter stiffness constants
• /parameters/particles/K_A: Area stiffness constants
• /parameters/particles/active_force_strengths: Active force strengths
• /parameters/particles/rotational_diffusion_constants: Rotational diffusion constants

For both new and existing files, the simulation data for the current step (obtained from output.steps_done) is saved in a new group named after the step number (e.g., /0, /1000, etc.).

Each step group (/<step_number>/) contains the following datasets, depending on the flags in parameters.dump_info:

• positions: Vector of SVector{2, Float64} (saved if dump_info.save_r is true)
• orientations: Vector{Float64} (saved if dump_info.save_r is true)
• forces: Vector of SVector{2, Float64} (saved if dump_info.save_F is true)
• potential_energy: Float64 (saved if dump_info.save_Epot is true)
• areas: Vector{Float64} (always saved)
• perimeters: Vector{Float64} (always saved)

Usage Example:

# Assuming `params`, `arrs`, `outs` are populated ParameterStruct, ArrayStruct, and Output instances
# and params.dump_info is configured.
save_simulation_state!(params, arrs, outs)

The function returns nothing.

update_areas!(parameters, arrays, output)

Calculates and updates the areas of all Voronoi cells. The area of each cell is computed using the shoelace formula (also known as Gauss's area formula), based on the coordinates of its Voronoi vertices. The results are stored in arrays.areas.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct. Used here to get N, the number of particles.
• arrays::ArrayStruct: The struct holding simulation arrays.
• output::Output: The simulation output struct. Not directly used in this function but included for consistency in function signatures across the module.

Notes

• This function assumes that the Voronoi vertex positions for each cell are already calculated and available in arrays.neighborlist.voronoi_vertices.
• The arrays.areas vector is updated in-place with the new area values. The shoelace formula calculates signed area, so the division by 2.0 yields the geometric area assuming a consistent vertex ordering.

update_perimeters!(parameters, arrays, output)

Calculates and updates the perimeters of all Voronoi cells. The perimeter of each cell is computed by summing the lengths of the segments connecting its Voronoi vertices in sequence. The results are stored in arrays.perimeters.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct. Used here to get N, the number of particles.
• arrays::ArrayStruct: The struct holding simulation arrays.
• output::Output: The simulation output struct. Not directly used in this function but included for consistency in function signatures across the module.

Notes

• This function assumes that the Voronoi vertex positions for each cell are already calculated and available in arrays.neighborlist.voronoi_vertices.
• The arrays.perimeters vector is updated in-place with the new perimeter values.

update_positions_with_pbcs!(parameters, arrays, output)

Creates an extended list of particle positions that includes the original particles and their periodic images located within a defined boundary layer around the simulation box. This extended list is essential for Voronoi tessellation algorithms to correctly construct Voronoi cells for particles near the boundaries of the periodic domain.

The function first adds all original particle positions. Then, for each original particle, it checks if it's close to any of the box boundaries (left, right, bottom, top) or corners. If a particle is within pbc_layer_depth of a boundary, its corresponding periodic image(s) across that boundary (or corner) are added to the list.

The generated list of positions (arrays.neighborlist.positions_with_pbc) and a corresponding list of original particle indices (arrays.neighborlist.position_indices) are stored in the arrays.neighborlist structure.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct, providing:
• arrays::ArrayStruct: The struct holding simulation arrays.
• output::Output: The simulation output struct. Not directly used in this function but included for consistency in function signatures.

Notes

• The fields arrays.neighborlist.positions_with_pbc and arrays.neighborlist.position_indices are modified in-place.

verify_tesselation(parameters, arrays, output)

Verifies the validity of the Voronoi tessellation by checking if all particles are correctly positioned with respect to the Delaunay facets of the tessellation. This function checks that for each Delaunay facet (triangle) defined by a triplet of particles, the circumcircle of the triangle does not contain any other particle that is not one of the triangle's vertices.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct, providing the number of particles N.
• arrays::ArrayStruct: The struct holding simulation arrays, including positions and neighbor lists.
• output::Output: The simulation output struct. Not directly used in this function but included for consistency in function signatures.

Returns

• Bool: Returns true if the tessellation is valid (no violations found), or false if any particle is found within the circumcircle of a Delaunay facet that is not one of its vertices.

Notes

• The function assumes that the Voronoi tessellation has already been computed and that the necessary fields in arrays.neighborlist are populated, including delaunay_facet_triplets, positions_with_pbc, and voronoi_neighbors.
• The field arrays.neighborlist.positionswithpbc is updated to include periodic boundary conditions.
• If the flag arrays.neighborlist.check_tesselation is set to false, the function immediately returns false, indicating that tessellation verification is not required or has been disabled.

voronoi_tesselation!(parameters, arrays, output)

Performs a full Voronoi tessellation of the system based on current particle positions. This function orchestrates several steps to compute the Voronoi diagram, which is dual to the Delaunay triangulation. The results, including lists of Voronoi neighbors, vertex positions, and topological information, are stored in arrays.neighborlist.

The process involves:

1. Update Positions with PBC: Calls update_positions_with_pbcs! to generate an extended list of particle positions, including periodic images necessary for correct tessellation across boundaries. This list is stored in arrays.neighborlist.positions_with_pbc.
2. Delaunay Triangulation: Computes the Delaunay triangulation of the extended particle positions using Quickhull.delaunay.
3. Process Delaunay Facets: Iterates through each facet (triangle) of the Delaunay triangulation:
   • The three particles forming the triangle are identified as Voronoi neighbors of each other. These are added to arrays.neighborlist.voronoi_neighbors.
   • The circumcenter of the Delaunay triangle is calculated using the circumcenter function. This circumcenter is a Voronoi vertex. It's added to a global list of voronoi_vertices and also to arrays.neighborlist.voronoi_vertex_positions_per_particle for each of the three particles forming the triangle.
   • The index of this new Voronoi vertex is added to arrays.neighborlist.voronoi_vertex_indices for each of the three particles.
   • A tuple of the three particle indices forming the triangle (cell centers sharing this Voronoi vertex) is stored in arrays.neighborlist.cell_centers_that_share_a_vertex.
4. Sort Voronoi Vertices: For each of the original N particles, the collected Voronoi vertices associated with it are sorted in counter-clockwise order around the particle's actual position using sort_indices_counter_clockwise. This ensures a consistent representation of Voronoi cells. The sorted vertex indices and positions replace the unsorted ones in arrays.neighborlist.voronoi_vertex_indices and arrays.neighborlist.voronoi_vertex_positions_per_particle for the original particles.
5. Store Results: All computed lists (global Voronoi vertices, neighbors per particle, vertex indices per particle, vertex positions per particle, and cell centers sharing a vertex) are stored in the corresponding fields of arrays.neighborlist.

Arguments

• parameters::ParameterStruct: The main simulation parameter struct.
• arrays::ArrayStruct: The struct holding simulation arrays.
• output::Output: The simulation output struct. Passed to update_positions_with_pbcs!.

Notes

• This function heavily modifies the arrays.neighborlist structure in-place.
• It relies on Quickhull.jl for the Delaunay triangulation and helper functions like update_positions_with_pbcs!, circumcenter, and sort_indices_counter_clockwise.

compute_msd(displacement_array)

returns [1,2,3,4,5,6,7,8,9,10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, etc] Computes the mean squared displacement (MSD) of particles based on their displacement vectors. The MSD is a measure of the average squared distance that particles have moved from their initial positions.