feat: enable discrete swaps

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
 

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
+end

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
 
-###############################################################################
+###############################################################################

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:

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

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,
             }
         )