|
| 1 | +# /// script |
| 2 | +# requires-python = ">=3.13" |
| 3 | +# dependencies = [ |
| 4 | +# "matplotlib>=3.10.8", |
| 5 | +# "numpy>=2.4.1", |
| 6 | +# "pandas>=2.3.3", |
| 7 | +# "seaborn>=0.13.2", |
| 8 | +# ] |
| 9 | +# /// |
| 10 | + |
| 11 | +# 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 |
| 12 | + |
| 13 | +import itertools |
| 14 | +import json |
| 15 | +import math |
| 16 | +import os |
| 17 | +import subprocess |
| 18 | + |
| 19 | +import matplotlib.pyplot as plt |
| 20 | +import numpy as np |
| 21 | +import pandas as pd |
| 22 | +import seaborn as sns |
| 23 | + |
| 24 | +# Different from publication, because we start from rectangular lattice instead of fcc |
| 25 | +# It shouldn't affect results, as we're not looking at crystallisation |
| 26 | +N = 1000 |
| 27 | +path_to_exe = "particlesmc" |
| 28 | + |
| 29 | +epsilon_1 = 1 |
| 30 | +epsilon_12 = 1.1523 |
| 31 | +epsilon_2 = 1.3702 |
| 32 | +sigma_1 = 1 |
| 33 | +sigma_12 = 1.0339 |
| 34 | +sigma_2 = 1.0640 |
| 35 | + |
| 36 | + |
| 37 | +def create_config(n1: int, n2: int, box_length: float, filepath: str) -> None: |
| 38 | + # We do a square lattice |
| 39 | + # Species are put randomly amongst the sites |
| 40 | + |
| 41 | + with open(filepath, "w") as f: |
| 42 | + f.write(f"{N:d}\n") |
| 43 | + f.write( |
| 44 | + 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' |
| 45 | + ) |
| 46 | + |
| 47 | + species_indices = [1] * n1 + [2] * n2 |
| 48 | + np.random.shuffle(species_indices) |
| 49 | + |
| 50 | + number_of_particle_in_each_direction = round(N ** (1 / 3)) |
| 51 | + if number_of_particle_in_each_direction**3 != N: |
| 52 | + raise RuntimeError("N is not x^3") |
| 53 | + dxdydz = box_length / number_of_particle_in_each_direction |
| 54 | + counter = 0 |
| 55 | + for i, j, k in itertools.product( |
| 56 | + range(number_of_particle_in_each_direction), |
| 57 | + range(number_of_particle_in_each_direction), |
| 58 | + range(number_of_particle_in_each_direction), |
| 59 | + ): |
| 60 | + x = (i + 0.5) * dxdydz - box_length / 2 |
| 61 | + y = (j + 0.5) * dxdydz - box_length / 2 |
| 62 | + z = (k + 0.5) * dxdydz - box_length / 2 |
| 63 | + f.write(f"{species_indices[counter]} {x} {y} {z}\n") |
| 64 | + counter += 1 |
| 65 | + |
| 66 | + |
| 67 | +def create_params(parameters: dict, path_to_params: str) -> None: |
| 68 | + with open(path_to_params, "w") as f: |
| 69 | + f.write(f""" |
| 70 | + [system] |
| 71 | + config = "{parameters["config"]}" |
| 72 | + temperature = {parameters["temperature"]} |
| 73 | + density = {parameters["density"]} |
| 74 | + list_type = "LinkedList" |
| 75 | +
|
| 76 | + [model] |
| 77 | + [model."1-1"] |
| 78 | + name = "LennardJones" |
| 79 | + epsilon = {parameters["epsilon_1"]:.2f} |
| 80 | + sigma = {parameters["sigma_1"]:.2f} |
| 81 | + rcut = {parameters["rcut"]:.2f} |
| 82 | + shift_potential = false |
| 83 | +
|
| 84 | + [model."1-2"] |
| 85 | + name = "LennardJones" |
| 86 | + epsilon = {parameters["epsilon_12"]} |
| 87 | + sigma = {parameters["sigma_12"]} |
| 88 | + rcut = {parameters["rcut"]:.2f} |
| 89 | + shift_potential = false |
| 90 | +
|
| 91 | + [model."2-2"] |
| 92 | + name = "LennardJones" |
| 93 | + epsilon = {parameters["epsilon_2"]} |
| 94 | + sigma = {parameters["sigma_2"]} |
| 95 | + rcut = {parameters["rcut"]:.2f} |
| 96 | + shift_potential = false |
| 97 | +
|
| 98 | + [simulation] |
| 99 | + type = "Metropolis" |
| 100 | + steps = {int(parameters["steps"]):d} |
| 101 | + seed = 42 |
| 102 | + parallel = false |
| 103 | + output_path = "./" |
| 104 | +
|
| 105 | + [[simulation.move]] |
| 106 | + action = "Displacement" |
| 107 | + probability = 1.0 |
| 108 | + policy = "SimpleGaussian" |
| 109 | + parameters = {{sigma = 0.05}} |
| 110 | +
|
| 111 | + [[simulation.output]] |
| 112 | + algorithm = "StoreCallbacks" |
| 113 | + callbacks = ["energy", "acceptance"] |
| 114 | + scheduler_params = {{linear_interval = 100}} |
| 115 | +
|
| 116 | + [[simulation.output]] |
| 117 | + algorithm = "StoreLastFrames" |
| 118 | + scheduler_params = {{linear_interval = 1000}} |
| 119 | + fmt = "EXYZ" |
| 120 | +
|
| 121 | + """) |
| 122 | + |
| 123 | + |
| 124 | +def run_simulations(output_path: str) -> None: |
| 125 | + df = pd.read_csv("reference-data.csv") |
| 126 | + |
| 127 | + path_to_config = "config.exyz" |
| 128 | + path_to_params = "params.toml" |
| 129 | + data = [] |
| 130 | + for i, row in df.iterrows(): |
| 131 | + workdir = f"./tmp/{i}" |
| 132 | + os.makedirs(workdir, exist_ok=True) |
| 133 | + |
| 134 | + print(row["t"], row["x"], row["density"]) |
| 135 | + |
| 136 | + # density = N / box_length**3 |
| 137 | + box_length = (N / row["density"]) ** (1 / 3) |
| 138 | + |
| 139 | + n2 = round(N * row["x"]) |
| 140 | + n1 = N - n2 |
| 141 | + create_config(n1, n2, box_length, f"{workdir}/{path_to_config}") |
| 142 | + |
| 143 | + # Generate input, and run |
| 144 | + rc = 4 * sigma_1 |
| 145 | + parameters = { |
| 146 | + "config": path_to_config, |
| 147 | + "temperature": row["t"], |
| 148 | + "epsilon_1": epsilon_1, |
| 149 | + "epsilon_12": epsilon_12, |
| 150 | + "epsilon_2": epsilon_2, |
| 151 | + "sigma_1": sigma_1, |
| 152 | + "sigma_12": sigma_12, |
| 153 | + "sigma_2": sigma_2, |
| 154 | + "steps": 1e3, # 1e4 equilibration in paper |
| 155 | + "rcut": rc, |
| 156 | + "density": N / box_length**3, |
| 157 | + } |
| 158 | + create_params(parameters, f"{workdir}/{path_to_params}") |
| 159 | + |
| 160 | + subprocess.run( |
| 161 | + [path_to_exe, path_to_params], cwd=workdir, stdout=subprocess.DEVNULL |
| 162 | + ) |
| 163 | + |
| 164 | + # Post-process the energies |
| 165 | + energies = pd.read_csv(f"{workdir}/energy.dat", sep="\\s+", names=["i", "e"])[ |
| 166 | + "e" |
| 167 | + ] |
| 168 | + # Remove the first half as equilibration, just to be sure |
| 169 | + energies = energies[int(len(energies) / 2) :] |
| 170 | + |
| 171 | + acceptance_rate = np.array( |
| 172 | + pd.read_csv(f"{workdir}/acceptance.dat", sep="\\s+", names=["i", "a"])["a"] |
| 173 | + )[-1] |
| 174 | + |
| 175 | + # Compute long-range corrections from the cutoff |
| 176 | + # Formula from Gromacs https://manual.gromacs.org/current/reference-manual/functions/long-range-vdw.html |
| 177 | + lr_correction = 0 |
| 178 | + |
| 179 | + c6_11 = 4 * epsilon_1 * sigma_1**6 |
| 180 | + rho_11 = n1 / box_length**3 |
| 181 | + lr_correction += -2 / 3 * math.pi * n1 * rho_11 * c6_11 / rc**3 |
| 182 | + |
| 183 | + c6_22 = 4 * epsilon_2 * sigma_2**6 |
| 184 | + rho_22 = n2 / box_length**3 |
| 185 | + lr_correction += -2 / 3 * math.pi * n2 * rho_22 * c6_22 / rc**3 |
| 186 | + |
| 187 | + c6_12 = 4 * epsilon_12 * sigma_12**6 |
| 188 | + lr_correction += -2 / 3 * math.pi * (n1 * rho_22 + n2 * rho_11) * c6_12 / rc**3 |
| 189 | + |
| 190 | + data.append( |
| 191 | + { |
| 192 | + "t": row["t"], |
| 193 | + "x": row["x"], |
| 194 | + "density": row["density"], |
| 195 | + "energy": np.mean(energies), |
| 196 | + "energy_err": np.std(energies) / np.sqrt(len(energies)), |
| 197 | + "lr_correction": lr_correction, |
| 198 | + "acceptance_rate": json.loads(acceptance_rate)[0], |
| 199 | + } |
| 200 | + ) |
| 201 | + |
| 202 | + df.merge(pd.DataFrame(data)).to_csv(output_path) |
| 203 | + |
| 204 | + |
| 205 | +def plot_results(path_to_energies: str) -> None: |
| 206 | + df = pd.read_csv(path_to_energies) |
| 207 | + |
| 208 | + # u (from the paper) is actually total energy / epsilon_11 N |
| 209 | + # epsilon_11 = 1, so it's the energy per particule |
| 210 | + # Add long range corrections to the MC results, which is also in energy / particle |
| 211 | + df["energy_lr"] = df["energy"] + df["lr_correction"] / N |
| 212 | + |
| 213 | + sns.scatterplot(data=df, x="u", y="energy_lr", hue="density", style="x") |
| 214 | + plt.axline((-5, -5), slope=1) |
| 215 | + plt.xlabel("Published energy") |
| 216 | + plt.ylabel("Re-computed from ParticleMC") |
| 217 | + plt.savefig("correlation-plot.jpeg") |
| 218 | + plt.show() |
| 219 | + |
| 220 | + |
| 221 | +if __name__ == "__main__": |
| 222 | + path_to_energies = "calculated-energies.csv" |
| 223 | + run_simulations(path_to_energies) |
| 224 | + plot_results(path_to_energies) |
0 commit comments