diff --git a/README.md b/README.md index ff5c0e8..125cdbd 100644 --- a/README.md +++ b/README.md @@ -93,7 +93,7 @@ rcut = 2.5 [simulation] type = "Metropolis" -steps = 500 +steps = 500 seed = 10 parallel = false output_path = "./" @@ -120,7 +120,21 @@ This example defines a minimal Monte Carlo simulation setup: - The `[[simulation.move]]` section describes the Monte Carlo move to use: a displacement move with probability 1.0, guided by a simple Gaussian policy with a standard deviation (`sigma`) of 0.05. - The `[[simulation.output]]` section configures the output: trajectories will be stored every 50 steps in the XYZ format. -By executing `particlesmc params.toml` you will run a basic Metropolis Monte Carlo simulation of particles interacting via the Lennard-Jones potential, using displacement moves, and periodically saving the system's trajectory. +By executing `particlesmc params.toml** you will run a basic Metropolis Monte Carlo simulation of particles interacting via the Lennard-Jones potential, using displacement moves, and periodically saving the system's trajectory. + +**Extending beyond this simple example:** + +More models and moves are available. + +For instance, in a simulation that contains multiple species, __DiscreteSwap__ moves can be added to swap species between to random particles: +```toml +[[simulation.move]] +action = "DiscreteSwap" +probability = 0.1 +policy = "DoubleUniform" +parameters = {species = [1, 2]} +``` +is a move that will be attempted 10% of the time (probability = 0.1). It will select a random particle _i_ of species _1_, a random particle _j_ of species _2_, and will attempt to transform _i_ into species _2_, and _j_ into species _1_. ## Contributing diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index 8120b1e..8f8dfd2 100644 --- a/src/ParticlesMC.jl +++ b/src/ParticlesMC.jl @@ -153,8 +153,8 @@ ParticlesMC implemented in Comonicon. error("Configuration file '$config' does not exist in the current path.") end list_type = get(system, "list_type", "LinkedList") # optional field - bonds = get(system, "bonds", nothing) - + bonds = get(system, "bonds", nothing) + # Extract simulation parameters sim = params["simulation"] steps = sim["steps"] @@ -175,7 +175,7 @@ ParticlesMC implemented in Comonicon. "list_type" => list_type, "bonds" => bonds, )) - else + else chains = load_chains(config, args=Dict( "temperature" => [temperature], "density" => [density], @@ -214,8 +214,27 @@ ParticlesMC implemented in Comonicon. else error("Unsupported policy: $policy for action: $action") end + elseif action == "DiscreteSwap" + if "species" in keys(parameters) + species = parameters["species"] + if length(species) != 2 || eltype(species) != Int + error("'species' for action: $action must be an array of two ints") + end + else + error("Missing parameter 'species' for action: $action") + end + + # Use a system to initialize (chains[1]) + # This is because the action needs the number of particles for each species + action_obj = DiscreteSwap(species, chains[1]) + param_obj = Vector{Float64}() + if policy == "DoubleUniform" + policy_obj = DoubleUniform() + else + error("Unsupported policy: $policy for action: $action") + end else - error("Unsupported action: $action") + error("Unsupported action: $action") end # Build move move_obj = Move(action_obj, policy_obj, param_obj, prob) @@ -257,15 +276,15 @@ ParticlesMC implemented in Comonicon. else error("Unsupported output algorithm: $alg") end - push!(algorithm_list, algorithm) + push!(algorithm_list, algorithm) end M=1 path = joinpath(output_path) simulation = Simulation(chains, algorithm_list, steps; path=path, verbose=true) - + # Run the simulation run!(simulation) -end +end -end \ No newline at end of file +end diff --git a/src/moves.jl b/src/moves.jl index b357eda..f3e25ca 100644 --- a/src/moves.jl +++ b/src/moves.jl @@ -146,6 +146,12 @@ mutable struct DiscreteSwap <: Action δe::Float64 end +function DiscreteSwap(species::Vector{Int}, system::Atoms) + particles_per_species_1 = count(i -> (i == species[1]), system.species) + particles_per_species_2 = count(i -> (i == species[2]), system.species) + return DiscreteSwap(1, 1, (species[1], species[2]), (particles_per_species_1, particles_per_species_2), 0.0) +end + """ Swap the species identifiers of particles `i` and `j` and return the total pre- and post-swap energies for those two particles. @@ -349,4 +355,4 @@ function Arianna.sample_action!(action::MoleculeFlip, ::DoubleUniform, parameter action.i, action.j = i, j end -############################################################################### \ No newline at end of file +############################################################################### diff --git a/validation/lj-mixture/README.md b/validation/lj-mixture/README.md index 82dea55..16e52e8 100644 --- a/validation/lj-mixture/README.md +++ b/validation/lj-mixture/README.md @@ -1,14 +1,14 @@ # Binary mixture of LJ particles -This directory contains scripts to run a simple simulation of a binary mixture of Lennard-Jones particles, in 3D. +This directory contains scripts to run a simple simulation of a binary mixture of Lennard-Jones particles, in 3D. Swap Monte Carlo moves are attempted on top of normal displacements. 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. +Reproducing the published results validates the implementation of the MC scheme, the LJ potential, the discrete swap algorithm, and energy evaluation. ## How to run -The `particlemc` exectuable should be in the PATH env variable, otherwise a specific path can +The `particlemc` executable 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: diff --git a/validation/lj-mixture/calculated-energies.csv b/validation/lj-mixture/calculated-energies.csv index 05385a6..1c4d964 100644 --- a/validation/lj-mixture/calculated-energies.csv +++ b/validation/lj-mixture/calculated-energies.csv @@ -1,24 +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 +,t,x,density,u,energy,energy_err,lr_correction,acceptance_rate_displacement,acceptance_rate_swap +0,1.2183,0.5,0.8,-6.991912,-6.840534104086828,0.017637447754306167,-151.9222906615277,0.524164659893668,0.39982694488095 +1,1.2183,0.5,0.7,-6.28399,-6.165202896818967,0.016206321890262023,-132.93200432883674,0.606704734205417,0.496359875082052 +2,1.2183,0.5,0.65,-5.86756,-5.755000382698143,0.013358234226850594,-123.43686116249125,0.645592770725351,0.539404849521612 +3,1.5096,0.25,0.8,-5.87326,-5.770095515865939,0.01683749104035944,-127.1877273361547,0.581054728757668,0.45646768643208085 +4,1.5096,0.25,0.7,-5.2830446,-5.18929850304146,0.013574981514565761,-111.28926141913534,0.6520411271727,0.5359636385336065 +5,1.5096,0.25,0.6,-4.57915,-4.496678165866352,0.004256448674720936,-95.39079550211606,0.71337500305741,0.6050265550096473 +6,1.5096,0.25,0.5,-3.84601,-3.744194994789512,0.009028764596754972,-79.4923295850967,0.768620741027334,0.6672667236886599 +7,1.5096,0.25,0.4,-3.145625,-3.0625659840628465,0.009621436608541251,-63.59386366807736,0.812516259864317,0.7167863465478488 +8,1.5096,0.25,0.3,-2.446133,-2.315936149881792,0.007309758603592841,-47.69539775105803,0.85386912504697,0.766972331072345 +9,1.5096,0.25,0.2,-1.71833,-1.5634718658838844,0.013177941604522733,-31.796931834038674,0.897846916017939,0.8230561136196368 +10,1.5096,0.25,0.1,-0.88341,-0.7540551691012838,0.016305680822504467,-15.898465917019347,0.946107305098426,0.8910946233564736 +11,1.5096,0.75,0.8,-7.528813,-7.337683683940237,0.010453572023241963,-178.9234450957788,0.513666068526016,0.4023630974877171 +12,1.5096,0.75,0.7,-6.89008,-6.74127765861963,0.013307090789325909,-156.55801445880638,0.59955150569123,0.4993833668171782 +13,1.5096,0.75,0.6,-6.0058145,-5.886840679436459,0.014665817674222618,-134.19258382183412,0.675006170410048,0.5875420205677003 +14,1.5096,0.75,0.55,-5.52616,-5.410653702890696,0.014852028782687345,-123.00986850334797,0.70719013979592,0.6243311519105683 +15,1.6438,0.5,0.8,-6.54401,-6.394357885128574,0.019019860433670554,-151.9222906615277,0.557749479128449,0.4403258210172458 +16,1.6438,0.5,0.7,-5.95657,-5.836950816506536,0.013006657681474675,-132.93200432883674,0.633923469126825,0.5291707278260697 +17,1.6438,0.5,0.6,-5.18836,-5.103247816622887,0.009429572650158194,-113.94171799614583,0.702751891703188,0.6058918306048973 +18,1.6438,0.5,0.5,-4.35859,-4.2967131610134635,0.012212677851235598,-94.95143166345483,0.757340564386839,0.6687187953772402 +19,1.6438,0.5,0.4,-3.549019,-3.441614177830049,0.008689906174574758,-75.96114533076386,0.80890073311142,0.7258170389672389 +20,1.6438,0.5,0.3,-2.768327,-2.6187924148641315,0.01460481900316066,-56.97085899807291,0.85185679312116,0.7742326895152467 +21,1.6438,0.5,0.2,-1.936931,-1.7582716394839677,0.016702400788528176,-37.98057266538193,0.896708447569303,0.8260597139617687 +22,1.6438,0.5,0.1,-1.004111,-0.8719412344268423,0.018406946799003905,-18.99028633269097,0.943625799651788,0.8937103415352177 diff --git a/validation/lj-mixture/correlation-plot.jpeg b/validation/lj-mixture/correlation-plot.jpeg index abfe4c1..9bed755 100644 Binary files a/validation/lj-mixture/correlation-plot.jpeg and b/validation/lj-mixture/correlation-plot.jpeg differ diff --git a/validation/lj-mixture/run-validation.py b/validation/lj-mixture/run-validation.py index 992e09e..c147c64 100644 --- a/validation/lj-mixture/run-validation.py +++ b/validation/lj-mixture/run-validation.py @@ -104,10 +104,16 @@ def create_params(parameters: dict, path_to_params: str) -> None: [[simulation.move]] action = "Displacement" - probability = 1.0 + probability = 0.9 policy = "SimpleGaussian" parameters = {{sigma = 0.05}} + [[simulation.move]] + action = "DiscreteSwap" + probability = 0.1 + policy = "DoubleUniform" + parameters = {{species = [1, 2]}} + [[simulation.output]] algorithm = "StoreCallbacks" callbacks = ["energy", "acceptance"] @@ -168,9 +174,11 @@ def run_simulations(output_path: str) -> None: # 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] + df_acceptance_rates = pd.read_csv( + f"{workdir}/acceptance.dat", sep="\\s+", names=["i", "move", "swap"] + ) + displacement_acceptance = float(df_acceptance_rates["move"].iloc[-1][1:-2]) + swap_acceptance = float(df_acceptance_rates["swap"].iloc[-1][:-1]) # Compute long-range corrections from the cutoff # Formula from Gromacs https://manual.gromacs.org/current/reference-manual/functions/long-range-vdw.html @@ -195,7 +203,8 @@ def run_simulations(output_path: str) -> None: "energy": np.mean(energies), "energy_err": np.std(energies) / np.sqrt(len(energies)), "lr_correction": lr_correction, - "acceptance_rate": json.loads(acceptance_rate)[0], + "acceptance_rate_displacement": displacement_acceptance, + "acceptance_rate_swap": swap_acceptance, } )