diff --git a/validation/lj-mixture/README.md b/validation/lj-mixture/README.md
new file mode 100644
index 0000000..82dea55
--- /dev/null
+++ b/validation/lj-mixture/README.md
@@ -0,0 +1,34 @@
+# Binary mixture of LJ particles
+
+This directory contains scripts to run a simple simulation of a binary mixture of Lennard-Jones particles, in 3D.
+
+The average energy is recorded, for different conditions (temperature, density, ratio between the two species), and compared to published results.
+
+Reproducing the published results validates the implementation of the MC scheme, the LJ potential, and interaction counting.
+
+## How to run
+
+The `particlemc` exectuable should be in the PATH env variable, otherwise a specific path can
+be specified at the top of the run-validation.py script.
+
+The `run-validation.py` script is meant to be run with [uv](https://docs.astral.sh/uv/), with the following command:
+```bash
+uv run --script run-validation.py
+```
+which automatically install all dependencies.
+
+The script will run simulations at all desired parameters, save the results to csv, and make a correlation plot with published results.
+
+## Results
+
+ParticleMC reproduces published results very well:
+
+
+![correlation-plot](correlation-plot.jpeg)
+
+On this plot, `x` is the ratio between the two types of particles, and the energy is per particle.
+
+For some setups, the energies from ParticleMC are a bit higher than published ones. This actually happens at low density, because the simulation times are too short and we don't reach proper convergence (timeseries of energies not shown).
+
+## Reference data
+Comes from Monte Carlo Simulations of Binary Lennard–Jones Mixtures: A Test of the van der Waals One-Fluid Model, https://doi.org/10.1023/A:1022614200488
diff --git a/validation/lj-mixture/calculated-energies.csv b/validation/lj-mixture/calculated-energies.csv
new file mode 100644
index 0000000..05385a6
--- /dev/null
+++ b/validation/lj-mixture/calculated-energies.csv
@@ -0,0 +1,24 @@
+,t,x,density,u,energy,energy_err,lr_correction,acceptance_rate
+0,1.2183,0.5,0.8,-6.991912,-6.824332775509206,0.007332089000306406,-151.9222906615277,0.523858
+1,1.2183,0.5,0.7,-6.28399,-6.15696393793348,0.0115122047621387,-132.93200432883674,0.605872
+2,1.2183,0.5,0.65,-5.86756,-5.747027114561171,0.008533981399021708,-123.43686116249125,0.643768
+3,1.5096,0.25,0.8,-5.87326,-5.728584068606758,0.021466385837235483,-127.1877273361547,0.580572
+4,1.5096,0.25,0.7,-5.2830446,-5.1486481236643895,0.015199454304553882,-111.28926141913534,0.65148
+5,1.5096,0.25,0.6,-4.57915,-4.457067486522425,0.012901052999909882,-95.39079550211606,0.713526
+6,1.5096,0.25,0.5,-3.84601,-3.7550796383339944,0.016097072216878654,-79.4923295850967,0.767107
+7,1.5096,0.25,0.4,-3.145625,-3.031679704426388,0.01957993121178323,-63.59386366807736,0.812887
+8,1.5096,0.25,0.3,-2.446133,-2.29186093973726,0.019961432466160035,-47.69539775105803,0.853689
+9,1.5096,0.25,0.2,-1.71833,-1.5204237232517432,0.018705478102023994,-31.796931834038674,0.897772
+10,1.5096,0.25,0.1,-0.88341,-0.7731743940092469,0.012321114326961577,-15.898465917019347,0.94422
+11,1.5096,0.75,0.8,-7.528813,-7.387557821039766,0.019429067147871005,-178.9234450957788,0.515488
+12,1.5096,0.75,0.7,-6.89008,-6.705393398859641,0.024860213510585486,-156.55801445880638,0.598998
+13,1.5096,0.75,0.6,-6.0058145,-5.844619979911088,0.01211174105592481,-134.19258382183412,0.674251
+14,1.5096,0.75,0.55,-5.52616,-5.422084348147788,0.005698415488936215,-123.00986850334797,0.708598
+15,1.6438,0.5,0.8,-6.54401,-6.36678624267275,0.008473648040505592,-151.9222906615277,0.557943
+16,1.6438,0.5,0.7,-5.95657,-5.844425928860567,0.013831459037555685,-132.93200432883674,0.634995
+17,1.6438,0.5,0.6,-5.18836,-5.085670690270235,0.010455755569991681,-113.94171799614583,0.702152
+18,1.6438,0.5,0.5,-4.35859,-4.21187686179021,0.010561762660606716,-94.95143166345483,0.758437
+19,1.6438,0.5,0.4,-3.549019,-3.4182878823413763,0.011275107483718705,-75.96114533076386,0.807701
+20,1.6438,0.5,0.3,-2.768327,-2.575420052801317,0.010850764609473712,-56.97085899807291,0.853085
+21,1.6438,0.5,0.2,-1.936931,-1.7475874342217939,0.01920265179952064,-37.98057266538193,0.89591
+22,1.6438,0.5,0.1,-1.004111,-0.8560035870798322,0.012956204273870288,-18.99028633269097,0.942764
diff --git a/validation/lj-mixture/correlation-plot.jpeg b/validation/lj-mixture/correlation-plot.jpeg
new file mode 100644
index 0000000..abfe4c1
Binary files /dev/null and b/validation/lj-mixture/correlation-plot.jpeg differ
diff --git a/validation/lj-mixture/reference-data.csv b/validation/lj-mixture/reference-data.csv
new file mode 100644
index 0000000..039dec6
--- /dev/null
+++ b/validation/lj-mixture/reference-data.csv
@@ -0,0 +1,24 @@
+t,x,density,u
+1.2183,0.5,0.8,-6.991912
+1.2183,0.5,0.7,-6.28399
+1.2183,0.5,0.65,-5.86756
+1.5096,0.25,0.8,-5.87326
+1.5096,0.25,0.7,-5.2830446
+1.5096,0.25,0.6,-4.57915
+1.5096,0.25,0.5,-3.846010
+1.5096,0.25,0.4,-3.145625
+1.5096,0.25,0.3,-2.446133
+1.5096,0.25,0.2,-1.718330
+1.5096,0.25,0.1,-0.883410
+1.5096,0.75,0.8,-7.528813
+1.5096,0.75,0.7,-6.89008
+1.5096,0.75,0.6,-6.0058145
+1.5096,0.75,0.55,-5.52616
+1.6438,0.5,0.8,-6.544010
+1.6438,0.5,0.7,-5.95657
+1.6438,0.5,0.6,-5.18836
+1.6438,0.5,0.5,-4.35859
+1.6438,0.5,0.4,-3.549019
+1.6438,0.5,0.3,-2.768327
+1.6438,0.5,0.2,-1.936931
+1.6438,0.5,0.1,-1.004111
diff --git a/validation/lj-mixture/run-validation.py b/validation/lj-mixture/run-validation.py
new file mode 100644
index 0000000..992e09e
--- /dev/null
+++ b/validation/lj-mixture/run-validation.py
@@ -0,0 +1,224 @@
+# /// script
+# requires-python = ">=3.13"
+# dependencies = [
+#     "matplotlib>=3.10.8",
+#     "numpy>=2.4.1",
+#     "pandas>=2.3.3",
+#     "seaborn>=0.13.2",
+# ]
+# ///
+
+# Comes from Monte Carlo Simulations of Binary Lennard–Jones Mixtures: A Test of the van der Waals One-Fluid Model, https://doi.org/10.1023/A:1022614200488
+
+import itertools
+import json
+import math
+import os
+import subprocess
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import seaborn as sns
+
+# Different from publication, because we start from rectangular lattice instead of fcc
+# It shouldn't affect results, as we're not looking at crystallisation
+N = 1000
+path_to_exe = "particlesmc"
+
+epsilon_1 = 1
+epsilon_12 = 1.1523
+epsilon_2 = 1.3702
+sigma_1 = 1
+sigma_12 = 1.0339
+sigma_2 = 1.0640
+
+
+def create_config(n1: int, n2: int, box_length: float, filepath: str) -> None:
+    # We do a square lattice
+    # Species are put randomly amongst the sites
+
+    with open(filepath, "w") as f:
+        f.write(f"{N:d}\n")
+        f.write(
+            f'Lattice="{box_length:.4f} 0.0 0.0 0.0 {box_length:.4f} 0.0 0.0 0.0 {box_length:.4f}" Properties=:species:S:1:pos:R:3\n'
+        )
+
+        species_indices = [1] * n1 + [2] * n2
+        np.random.shuffle(species_indices)
+
+        number_of_particle_in_each_direction = round(N ** (1 / 3))
+        if number_of_particle_in_each_direction**3 != N:
+            raise RuntimeError("N is not x^3")
+        dxdydz = box_length / number_of_particle_in_each_direction
+        counter = 0
+        for i, j, k in itertools.product(
+            range(number_of_particle_in_each_direction),
+            range(number_of_particle_in_each_direction),
+            range(number_of_particle_in_each_direction),
+        ):
+            x = (i + 0.5) * dxdydz - box_length / 2
+            y = (j + 0.5) * dxdydz - box_length / 2
+            z = (k + 0.5) * dxdydz - box_length / 2
+            f.write(f"{species_indices[counter]} {x} {y} {z}\n")
+            counter += 1
+
+
+def create_params(parameters: dict, path_to_params: str) -> None:
+    with open(path_to_params, "w") as f:
+        f.write(f"""
+        [system]
+        config = "{parameters["config"]}"
+        temperature = {parameters["temperature"]}
+        density = {parameters["density"]}
+        list_type = "LinkedList"
+
+        [model]
+        [model."1-1"]
+        name = "LennardJones"
+        epsilon = {parameters["epsilon_1"]:.2f}
+        sigma = {parameters["sigma_1"]:.2f}
+        rcut = {parameters["rcut"]:.2f}
+        shift_potential = false
+
+        [model."1-2"]
+        name = "LennardJones"
+        epsilon = {parameters["epsilon_12"]}
+        sigma = {parameters["sigma_12"]}
+        rcut = {parameters["rcut"]:.2f}
+        shift_potential = false
+
+        [model."2-2"]
+        name = "LennardJones"
+        epsilon = {parameters["epsilon_2"]}
+        sigma = {parameters["sigma_2"]}
+        rcut = {parameters["rcut"]:.2f}
+        shift_potential = false
+
+        [simulation]
+        type = "Metropolis"
+        steps = {int(parameters["steps"]):d}
+        seed = 42
+        parallel = false
+        output_path = "./"
+
+        [[simulation.move]]
+        action = "Displacement"
+        probability = 1.0
+        policy = "SimpleGaussian"
+        parameters = {{sigma = 0.05}}
+
+        [[simulation.output]]
+        algorithm = "StoreCallbacks"
+        callbacks = ["energy", "acceptance"]
+        scheduler_params = {{linear_interval = 100}}
+
+        [[simulation.output]]
+        algorithm = "StoreLastFrames"
+        scheduler_params = {{linear_interval = 1000}}
+        fmt = "EXYZ"
+
+        """)
+
+
+def run_simulations(output_path: str) -> None:
+    df = pd.read_csv("reference-data.csv")
+
+    path_to_config = "config.exyz"
+    path_to_params = "params.toml"
+    data = []
+    for i, row in df.iterrows():
+        workdir = f"./tmp/{i}"
+        os.makedirs(workdir, exist_ok=True)
+
+        print(row["t"], row["x"], row["density"])
+
+        # density = N / box_length**3
+        box_length = (N / row["density"]) ** (1 / 3)
+
+        n2 = round(N * row["x"])
+        n1 = N - n2
+        create_config(n1, n2, box_length, f"{workdir}/{path_to_config}")
+
+        # Generate input, and run
+        rc = 4 * sigma_1
+        parameters = {
+            "config": path_to_config,
+            "temperature": row["t"],
+            "epsilon_1": epsilon_1,
+            "epsilon_12": epsilon_12,
+            "epsilon_2": epsilon_2,
+            "sigma_1": sigma_1,
+            "sigma_12": sigma_12,
+            "sigma_2": sigma_2,
+            "steps": 1e3,  # 1e4 equilibration in paper
+            "rcut": rc,
+            "density": N / box_length**3,
+        }
+        create_params(parameters, f"{workdir}/{path_to_params}")
+
+        subprocess.run(
+            [path_to_exe, path_to_params], cwd=workdir, stdout=subprocess.DEVNULL
+        )
+
+        # Post-process the energies
+        energies = pd.read_csv(f"{workdir}/energy.dat", sep="\\s+", names=["i", "e"])[
+            "e"
+        ]
+        # Remove the first half as equilibration, just to be sure
+        energies = energies[int(len(energies) / 2) :]
+
+        acceptance_rate = np.array(
+            pd.read_csv(f"{workdir}/acceptance.dat", sep="\\s+", names=["i", "a"])["a"]
+        )[-1]
+
+        # Compute long-range corrections from the cutoff
+        # Formula from Gromacs https://manual.gromacs.org/current/reference-manual/functions/long-range-vdw.html
+        lr_correction = 0
+
+        c6_11 = 4 * epsilon_1 * sigma_1**6
+        rho_11 = n1 / box_length**3
+        lr_correction += -2 / 3 * math.pi * n1 * rho_11 * c6_11 / rc**3
+
+        c6_22 = 4 * epsilon_2 * sigma_2**6
+        rho_22 = n2 / box_length**3
+        lr_correction += -2 / 3 * math.pi * n2 * rho_22 * c6_22 / rc**3
+
+        c6_12 = 4 * epsilon_12 * sigma_12**6
+        lr_correction += -2 / 3 * math.pi * (n1 * rho_22 + n2 * rho_11) * c6_12 / rc**3
+
+        data.append(
+            {
+                "t": row["t"],
+                "x": row["x"],
+                "density": row["density"],
+                "energy": np.mean(energies),
+                "energy_err": np.std(energies) / np.sqrt(len(energies)),
+                "lr_correction": lr_correction,
+                "acceptance_rate": json.loads(acceptance_rate)[0],
+            }
+        )
+
+        df.merge(pd.DataFrame(data)).to_csv(output_path)
+
+
+def plot_results(path_to_energies: str) -> None:
+    df = pd.read_csv(path_to_energies)
+
+    # u (from the paper) is actually total energy / epsilon_11 N
+    # epsilon_11 = 1, so it's the energy per particule
+    # Add long range corrections to the MC results, which is also in energy / particle
+    df["energy_lr"] = df["energy"] + df["lr_correction"] / N
+
+    sns.scatterplot(data=df, x="u", y="energy_lr", hue="density", style="x")
+    plt.axline((-5, -5), slope=1)
+    plt.xlabel("Published energy")
+    plt.ylabel("Re-computed from ParticleMC")
+    plt.savefig("correlation-plot.jpeg")
+    plt.show()
+
+
+if __name__ == "__main__":
+    path_to_energies = "calculated-energies.csv"
+    run_simulations(path_to_energies)
+    plot_results(path_to_energies)