diff --git a/examples/ortho-terphenyl/1-create-config-right-density/create-initial-config.py b/examples/ortho-terphenyl/1-create-config-right-density/create-initial-config.py
new file mode 100644
index 0000000..27229ce
--- /dev/null
+++ b/examples/ortho-terphenyl/1-create-config-right-density/create-initial-config.py
@@ -0,0 +1,69 @@
+# /// script
+# requires-python = ">=3.13"
+# dependencies = [
+#     "numpy>=2.4.1",
+# ]
+# ///
+import itertools
+
+import numpy as np
+
+M = 1000
+N = M * 3
+number_density = 0.2  # This is per particle, not per molecule
+# rho = N/V
+V = N / number_density
+L = V ** (1 / 3)
+
+sigmas = [0.9, 1, 1.1]
+
+
+def main(output_file_name: str) -> None:
+    """Create a regular lattice of molecules. These planar molecules are in the x-y plane."""
+
+    number_in_each_direction = round(M ** (1 / 3))
+    dxdydz = L / number_in_each_direction
+    assert number_in_each_direction**3 == M, "M is not an int to power 3"
+
+    r_ab = (sigmas[0] + sigmas[1]) / 2
+    r_ac = (sigmas[0] + sigmas[2]) / 2
+    cos_alpha = np.cos(60 / 180)
+    sin_alpha = np.sin(60 / 180)
+
+    f = open(output_file_name, "w")
+
+    f.write(f"{N}\n")
+    f.write(f"columns:molecule,species,position cell:{L},{L},{L}\n")
+
+    counter = 1
+    for i, j, k in itertools.product(
+        range(number_in_each_direction),
+        range(number_in_each_direction),
+        range(number_in_each_direction),
+    ):
+        r_a = (i * dxdydz, j * dxdydz, k * dxdydz)
+        r_b = (i * dxdydz, j * dxdydz + r_ab, k * dxdydz)
+        r_c = (i * dxdydz + r_ac * cos_alpha, j * dxdydz + r_ac * sin_alpha, k * dxdydz)
+
+        f.write(f"{counter} 1 {r_a[0]} {r_a[1]} {r_a[2]}\n")
+        f.write(f"{counter} 2 {r_b[0]} {r_b[1]} {r_b[2]}\n")
+        f.write(f"{counter} 3 {r_c[0]} {r_c[1]} {r_c[2]}\n")
+
+        counter += 1
+
+    number_of_bonds = M * 3
+    f.write(f"{number_of_bonds}\n")
+    f.write("columns:bond\n")
+    for i in range(M):
+        index_of_a = 1 + i * 3
+        # AB
+        f.write(f"{index_of_a} {index_of_a + 1}\n")
+        # BC
+        f.write(f"{index_of_a + 1} {index_of_a + 2}\n")
+        # AC
+        f.write(f"{index_of_a} {index_of_a + 2}\n")
+
+
+if __name__ == "__main__":
+    output_file_name = "inputframe.xyz"
+    main(output_file_name)
diff --git a/examples/ortho-terphenyl/1-create-config-right-density/params-template.toml b/examples/ortho-terphenyl/1-create-config-right-density/params-template.toml
new file mode 100644
index 0000000..8c417ee
--- /dev/null
+++ b/examples/ortho-terphenyl/1-create-config-right-density/params-template.toml
@@ -0,0 +1,82 @@
+[system]
+config = "./inputframe.xyz"
+temperature = 2.0
+density = DENSITY
+list_type = "LinkedList"
+
+[model]
+
+[model."1-1"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 0.9
+k = 37.03703703703703
+r0 = 1.35
+
+[model."1-2"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 0.95
+k = 33.24099722991689
+r0 = 1.425
+
+[model."1-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.0
+k = 30.0
+r0 = 1.5
+
+[model."2-2"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.0
+k = 30.0
+r0 = 1.5
+
+[model."2-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.05
+k = 27.2108843537415
+r0 = 1.575
+
+[model."3-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.1
+k = 24.79338842975207
+r0 = 1.65
+
+
+
+[simulation]
+type = "Metropolis"
+steps = 10000
+seed = 10
+parallel = false
+
+[[simulation.move]]
+action = "Displacement"
+probability = 0.8
+policy = "SimpleGaussian"
+parameters = {sigma = 0.05}
+
+[[simulation.move]]
+action = "MoleculeFlip"
+probability = 0.2
+policy = "DoubleUniform"
+
+[[simulation.output]]
+algorithm = "StoreCallbacks"
+callbacks = ["energy", "acceptance"]
+scheduler_params = {linear_interval = 1}
+
+[[simulation.output]]
+algorithm = "StoreLastFrames"
+scheduler_params = {linear_interval = 10000}
+fmt = "XYZ"
+
+[[simulation.output]]
+algorithm = "PrintTimeSteps"
+scheduler_params = {linear_interval = 1}
diff --git a/examples/ortho-terphenyl/1-create-config-right-density/run.sh b/examples/ortho-terphenyl/1-create-config-right-density/run.sh
new file mode 100644
index 0000000..b54339f
--- /dev/null
+++ b/examples/ortho-terphenyl/1-create-config-right-density/run.sh
@@ -0,0 +1,16 @@
+#!/bin/bash
+
+uv run --script create-initial-config.py
+
+density=0.2
+
+while (( $(echo "$density < 1.2" | bc -l ) ))
+do
+    echo "Density" $density
+
+    sed "s/DENSITY/$density/" params-template.toml > params.toml
+    particlesmc params.toml
+    cp trajectories/1/lastframe.xyz inputframe.xyz
+
+    density=$(echo "$density" | awk '{printf "%f", $1 * 1.1}')
+done
diff --git a/examples/ortho-terphenyl/2-equilibrate-at-different-temperatures/create-folders.sh b/examples/ortho-terphenyl/2-equilibrate-at-different-temperatures/create-folders.sh
new file mode 100644
index 0000000..6c28da3
--- /dev/null
+++ b/examples/ortho-terphenyl/2-equilibrate-at-different-temperatures/create-folders.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+for temperature in 3.0 2.0 1.6 1.4 1.25 1.2 1.15 1.1 1.05 1.0
+do
+    mkdir -p $temperature
+    sed "s/TEMPERATURE/$temperature/" params-template.toml > $temperature/params.toml
+done
diff --git a/examples/ortho-terphenyl/2-equilibrate-at-different-temperatures/params-template.toml b/examples/ortho-terphenyl/2-equilibrate-at-different-temperatures/params-template.toml
new file mode 100644
index 0000000..72a2d86
--- /dev/null
+++ b/examples/ortho-terphenyl/2-equilibrate-at-different-temperatures/params-template.toml
@@ -0,0 +1,87 @@
+[system]
+config = "../../1-create-config-right-density/inputframe.xyz"
+temperature = TEMPERATURE
+density = 1.2
+list_type = "LinkedList"
+
+[model]
+
+[model."1-1"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 0.9
+k = 37.03703703703703
+r0 = 1.35
+
+[model."1-2"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 0.95
+k = 33.24099722991689
+r0 = 1.425
+
+[model."1-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.0
+k = 30.0
+r0 = 1.5
+
+[model."2-2"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.0
+k = 30.0
+r0 = 1.5
+
+[model."2-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.05
+k = 27.2108843537415
+r0 = 1.575
+
+[model."3-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.1
+k = 24.79338842975207
+r0 = 1.65
+
+
+
+[simulation]
+type = "Metropolis"
+steps = 5000000
+seed = 10
+parallel = false
+
+[[simulation.move]]
+action = "Displacement"
+probability = 0.8
+policy = "SimpleGaussian"
+parameters = {sigma = 0.05}
+
+[[simulation.move]]
+action = "MoleculeFlip"
+probability = 0.2
+policy = "DoubleUniform"
+
+[[simulation.output]]
+algorithm = "StoreCallbacks"
+callbacks = ["energy", "acceptance"]
+scheduler_params = {linear_interval = 1000}
+
+[[simulation.output]]
+algorithm = "StoreLastFrames"
+scheduler_params = {linear_interval = 10000}
+fmt = "XYZ"
+
+[[simulation.output]]
+algorithm = "StoreTrajectories"
+scheduler_params = {linear_interval = 1000}
+fmt = "XYZ"
+
+[[simulation.output]]
+algorithm = "PrintTimeSteps"
+scheduler_params = {linear_interval = 1000}
diff --git a/examples/ortho-terphenyl/3-run-production/create-folders.sh b/examples/ortho-terphenyl/3-run-production/create-folders.sh
new file mode 100644
index 0000000..6c28da3
--- /dev/null
+++ b/examples/ortho-terphenyl/3-run-production/create-folders.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+for temperature in 3.0 2.0 1.6 1.4 1.25 1.2 1.15 1.1 1.05 1.0
+do
+    mkdir -p $temperature
+    sed "s/TEMPERATURE/$temperature/" params-template.toml > $temperature/params.toml
+done
diff --git a/examples/ortho-terphenyl/3-run-production/params-template.toml b/examples/ortho-terphenyl/3-run-production/params-template.toml
new file mode 100644
index 0000000..8dd1bb8
--- /dev/null
+++ b/examples/ortho-terphenyl/3-run-production/params-template.toml
@@ -0,0 +1,87 @@
+[system]
+config = "../../2-equilibrate-at-different-temperatures/TEMPERATURE/trajectories/1/lastframe.xyz"
+temperature = TEMPERATURE
+density = 1.2
+list_type = "LinkedList"
+
+[model]
+
+[model."1-1"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 0.9
+k = 37.03703703703703
+r0 = 1.35
+
+[model."1-2"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 0.95
+k = 33.24099722991689
+r0 = 1.425
+
+[model."1-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.0
+k = 30.0
+r0 = 1.5
+
+[model."2-2"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.0
+k = 30.0
+r0 = 1.5
+
+[model."2-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.05
+k = 27.2108843537415
+r0 = 1.575
+
+[model."3-3"]
+name = "GeneralKG"
+epsilon = 1.0
+sigma = 1.1
+k = 24.79338842975207
+r0 = 1.65
+
+
+
+[simulation]
+type = "Metropolis"
+steps = 5000000
+seed = 10
+parallel = false
+
+[[simulation.move]]
+action = "Displacement"
+probability = 0.8
+policy = "SimpleGaussian"
+parameters = {sigma = 0.05}
+
+[[simulation.move]]
+action = "MoleculeFlip"
+probability = 0.2
+policy = "DoubleUniform"
+
+[[simulation.output]]
+algorithm = "StoreCallbacks"
+callbacks = ["energy", "acceptance"]
+scheduler_params = {linear_interval = 1000}
+
+[[simulation.output]]
+algorithm = "StoreLastFrames"
+scheduler_params = {linear_interval = 10000}
+fmt = "XYZ"
+
+[[simulation.output]]
+algorithm = "StoreTrajectories"
+scheduler_params = {linear_interval = 4096, log_base=2.0}
+fmt = "XYZ"
+
+[[simulation.output]]
+algorithm = "PrintTimeSteps"
+scheduler_params = {linear_interval = 1000}
diff --git a/examples/ortho-terphenyl/4-compute-correlation-functions/compute-observables.py b/examples/ortho-terphenyl/4-compute-correlation-functions/compute-observables.py
new file mode 100644
index 0000000..35d4dd7
--- /dev/null
+++ b/examples/ortho-terphenyl/4-compute-correlation-functions/compute-observables.py
@@ -0,0 +1,56 @@
+# /// script
+# requires-python = ">=3.14"
+# dependencies = [
+#     "atooms-pp>=4.2.1",
+#     "fastparquet>=2025.12.0",
+#     "pandas>=3.0.0",
+# ]
+# ///
+
+import atooms.postprocessing as pp
+import numpy as np
+import pandas as pd
+from atooms.trajectory import Trajectory
+
+temperatures = [3.0, 2.0, 1.6, 1.4, 1.25, 1.2, 1.15, 1.1, 1.05, 1.0]
+
+
+def compute_fskt() -> pd.DataFrame:
+    print("F_s(k,t)")
+    df_list = []
+    for T in temperatures:
+        print(f"T = {T}")
+        traj = Trajectory(f"../3-run-production/{T}/trajectories/1/trajectory.xyz")
+
+        cf = pp.SelfIntermediateScatteringFast(
+            traj,
+            ksamples=1,
+            kmin=7.4,
+            kmax=7.4,
+            tgrid=sorted({round(x) for x in np.logspace(0, 6, num=55)}),
+        )
+        cf.compute()
+
+        print(f"k = {cf.grid[0][0]}")
+
+        N = len(cf.grid[1])
+        df_list.append(
+            pd.DataFrame(
+                {
+                    "t": cf.grid[1],
+                    cf.qualified_name: np.array(cf.values)[0],
+                    "T": [T] * N,
+                }
+            )
+        )
+
+    return pd.concat(df_list).reset_index(drop=True)
+
+
+def main() -> None:
+    df = compute_fskt()
+    df.to_parquet("fskt.parquet")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/ortho-terphenyl/4-compute-correlation-functions/fskt.parquet b/examples/ortho-terphenyl/4-compute-correlation-functions/fskt.parquet
new file mode 100644
index 0000000..5341d9f
Binary files /dev/null and b/examples/ortho-terphenyl/4-compute-correlation-functions/fskt.parquet differ
diff --git a/examples/ortho-terphenyl/4-compute-correlation-functions/fskt.png b/examples/ortho-terphenyl/4-compute-correlation-functions/fskt.png
new file mode 100644
index 0000000..8c48b01
Binary files /dev/null and b/examples/ortho-terphenyl/4-compute-correlation-functions/fskt.png differ
diff --git a/examples/ortho-terphenyl/4-compute-correlation-functions/plot-comparison.py b/examples/ortho-terphenyl/4-compute-correlation-functions/plot-comparison.py
new file mode 100644
index 0000000..cce41f2
--- /dev/null
+++ b/examples/ortho-terphenyl/4-compute-correlation-functions/plot-comparison.py
@@ -0,0 +1,47 @@
+# /// script
+# requires-python = ">=3.13"
+# dependencies = [
+#     "fastparquet>=2025.12.0",
+#     "seaborn>=0.13.2",
+# ]
+# ///
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import seaborn as sns
+
+
+def main() -> None:
+    df_here = pd.read_parquet("fskt.parquet")
+    df_here["method"] = "ParticlesMC"
+
+    temperatures = [3, 2, 1.6, 1.4, 1.25, 1.2, 1.15, 1.1, 1.05, 1]
+    assert len(temperatures) == len(df_here["T"].unique())
+    df_list = []
+    for T in temperatures:
+        df = pd.read_csv(
+            f"4/fsqt_flip/temp_{T}.dat", sep=" ", skiprows=1, names=["t", "F_s(k,t)"]
+        )
+        df["T"] = T
+        df_list.append(df)
+    df_paper = pd.concat(df_list)
+    df_paper["method"] = "Paper"
+
+    df_tot = pd.concat([df_here, df_paper]).reset_index(drop=True)
+    df_tot["Temperature"] = df_tot["T"].apply(lambda x: str(x))
+
+    ax = sns.lineplot(
+        data=df_tot,
+        x="t",
+        y="F_s(k,t)",
+        hue="Temperature",
+        style="method",
+        markers={"ParticlesMC": "o", "Paper": ","},
+    )
+    ax.set_xscale("log")
+    plt.savefig("fskt.png")
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/ortho-terphenyl/README.md b/examples/ortho-terphenyl/README.md
new file mode 100644
index 0000000..bbefaf6
--- /dev/null
+++ b/examples/ortho-terphenyl/README.md
@@ -0,0 +1,41 @@
+# Simulation of ortho-terphenyl (MC swap between atoms of the same molecule)
+
+In this example, we simulate 1000 molecules of ortho-terphenyl (modeled as molecules containing three particles) at a number density of 1.2, for multiple temperatures.
+This mirrors what was done in the publication [Molecular motion at the experimental glass transition](https://doi.org/10.1103/4twk-33j7).
+
+In every case, we will run MC simulations that have swap moves between atoms of the same molecule.
+
+This example is split into multiple steps:
+
+## 1. Create a configuration at the correct density
+
+We want to run our simulations at a relatively high density.
+
+The easiest way to create a decent starting point is to first generate a structure at a low density (we can put molecules on a simple lattice, all in the same orientation), and slowly increase the density until we reach the desired one.
+This density increase is done by running multiple consecutive simulations, downscaling the simulation box between each simulation.
+
+The scripts to do this are in the [1-create-config-right-density](./1-create-config-right-density) subfolder. The process is to first run the python script, preferably using [uv](https://doi.org/10.1103/4twk-33j7). then the run.sh script, which will run the MC simulations.
+
+## 2 - Create equilibrated configurations at various temperatures
+
+The configuration generated at the previous step needs to be further equilibrated if we want representative starting points for production runs at different temperatures.
+
+The next step is therefore to run long simulations, one for each target temperature, to produce inputs for the production runs.
+
+This is done by running the bash script in the [2-equilibrate-at-different-temperatures](./2-equilibrate-at-different-temperatures) subfolder.
+
+## 3 - Run production simulations
+
+For each of the temperature, we run a long simulation to accumulate data, sung the bash script in the [3-run-production](./3-run-production) subfolder.
+
+## 4 - Analyze results
+
+Many analyses could be done. Here we look at a single time correlation function as an example. It corresponds to Fig. 4a of the paper listed above.
+
+The correlation function is computed using an external python package, [atooms-pp](https://framagit.org/atooms/postprocessing).
+The results are compared to the published results, to obtain the following plot
+
+
+![fsqt](./4-compute-correlation-functions/fskt.png)
+
+There is a very good agreement between two sets of results. As these are MC simulations, the time has no proper physical meaning, and changing MC parameters (like amplitude of displacements, ratio between types of moves) will affect the time scales. So actual quantitative comparison can't be done.
diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl
index 8f8dfd2..4199345 100644
--- a/src/ParticlesMC.jl
+++ b/src/ParticlesMC.jl
@@ -271,7 +271,7 @@ ParticlesMC implemented in Comonicon.
         elseif alg == "PrintTimeSteps"
             algorithm = (
                 algorithm = eval(Meta.parse(alg)),
-                scheduler = build_schedule(steps, burn, steps ÷ 10),
+                scheduler = sched,
             )
         else
             error("Unsupported output algorithm: $alg")