Chore/compat arianna v0.2

Project.toml

@@ -18,7 +18,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.8"
-Arianna = "0.1"
+Arianna = "0.2"
 Comonicon = "1.0"
 ComponentArrays = "0.15"
 ConcreteStructs = "0.2"

examples/lj-mixture/run-validation.py

@@ -116,7 +116,12 @@ def create_params(parameters: dict, path_to_params: str) -> None:
 
         [[simulation.output]]
         algorithm = "StoreCallbacks"
-        callbacks = ["energy", "acceptance"]
+        callbacks = ["energy"]
+        scheduler_params = {{linear_interval = 100}}
+
+        [[simulation.output]]
+        algorithm = "StoreAcceptance"
+        dependencies = ["Metropolis"]
         scheduler_params = {{linear_interval = 100}}
 
         [[simulation.output]]
@@ -128,12 +133,12 @@ def create_params(parameters: dict, path_to_params: str) -> None:
 
 
 def run_simulations(output_path: str) -> None:
-    df = pd.read_csv("reference-data.csv")
+    df_ref = pd.read_csv("reference-data.csv")
 
     path_to_config = "config.exyz"
     path_to_params = "params.toml"
     data = []
-    for i, row in df.iterrows():
+    for i, row in df_ref.iterrows():
         workdir = f"./tmp/{i}"
         os.makedirs(workdir, exist_ok=True)
 
@@ -168,17 +173,31 @@ def run_simulations(output_path: str) -> None:
         )
 
         # Post-process the energies
-        energies = pd.read_csv(f"{workdir}/energy.dat", sep="\\s+", names=["i", "e"])[
+        energies = pd.read_csv(f"{workdir}/chains/1/energy.dat", sep="\\s+", names=["i", "e"])[
             "e"
         ]
         # Remove the first half as equilibration, just to be sure
         energies = energies[int(len(energies) / 2) :]
+
+        moves = {
+            1: "displacement",
+            2: "swap",
+        }
 
-        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])
+        dfs = []
+
+        for move_id, move_name in moves.items():
+            path = f"{workdir}/moves/{move_id}/acceptance.dat"
+            df = pd.read_csv(path, sep=r"\s+", names=["i", move_name])
+            dfs.append(df)
+
+        # Merge all on column "i"
+        df_acceptance_rates = dfs[0]
+        for df in dfs[1:]:
+            df_acceptance_rates = df_acceptance_rates.merge(df, on="i")
+
+        displacement_acceptance = float(df_acceptance_rates["displacement"].iloc[-1])
+        swap_acceptance = float(df_acceptance_rates["swap"].iloc[-1])
 
         # Compute long-range corrections from the cutoff
         # Formula from Gromacs https://manual.gromacs.org/current/reference-manual/functions/long-range-vdw.html
@@ -207,8 +226,7 @@ def run_simulations(output_path: str) -> None:
                 "acceptance_rate_swap": swap_acceptance,
             }
         )
-
-        df.merge(pd.DataFrame(data)).to_csv(output_path)
+        df_ref.merge(pd.DataFrame(data)).to_csv(output_path)
 
 
 def plot_results(path_to_energies: str) -> None:

examples/movie/params.toml

@@ -20,7 +20,12 @@ parameters = {sigma = 0.05}
 
 [[simulation.output]]
 algorithm = "StoreCallbacks"
-callbacks = ["energy", "acceptance"]
+callbacks = ["energy"]
+scheduler_params = {linear_interval = 500}
+
+[[simulation.output]]
+algorithm = "StoreAcceptance"
+dependencies = ["Metropolis"]
 scheduler_params = {linear_interval = 500}
 
 [[simulation.output]]

examples/ortho-terphenyl/2-equilibrate-at-different-temperatures/params-template.toml

@@ -69,7 +69,12 @@ policy = "DoubleUniform"
 
 [[simulation.output]]
 algorithm = "StoreCallbacks"
-callbacks = ["energy", "acceptance"]
+callbacks = ["energy"]
+scheduler_params = {linear_interval = 1000}
+
+[[simulation.output]]
+algorithm = "StoreAcceptance"
+dependencies = ["Metropolis"]
 scheduler_params = {linear_interval = 1000}
 
 [[simulation.output]]

examples/ortho-terphenyl/3-run-production/params-template.toml

@@ -69,7 +69,7 @@ policy = "DoubleUniform"
 
 [[simulation.output]]
 algorithm = "StoreCallbacks"
-callbacks = ["energy", "acceptance"]
+callbacks = ["energy"]
 scheduler_params = {linear_interval = 1000}
 
 [[simulation.output]]

examples/ortho-terphenyl/4-compute-correlation-functions/compute-observables.py

@@ -20,7 +20,7 @@ def compute_fskt() -> pd.DataFrame:
     df_list = []
     for T in temperatures:
         print(f"T = {T}")
-        traj = Trajectory(f"../3-run-production/{T}/trajectories/1/trajectory.xyz")
+        traj = Trajectory(f"../3-run-production/{T}/chains/1/trajectory.xyz")
 
         cf = pp.SelfIntermediateScatteringFast(
             traj,

src/ParticlesMC.jl

@@ -112,7 +112,7 @@ function compute_energy_particle(system::Particles, ids::AbstractVector)
 end
 
 
-export callback_energy
+export energy
 #export nearest_image_distance
 export Model, GeneralKG, JBB, BHHP, SoftSpheres, KobAndersen, Trimer
 export NeighbourList, LinkedList, CellList, EmptyList, VerletList
@@ -249,6 +249,7 @@ ParticlesMC implemented in Comonicon.
     for output in sim["output"]
         alg = output["algorithm"]
         scheduler_params = output["scheduler_params"]
+        dependencies = get(output, "dependencies", nothing)
         callbacks = get(output, "callbacks", [])
         fmt = get(output, "fmt", "XYZ")
         interval = scheduler_params["linear_interval"]
@@ -259,12 +260,19 @@ ParticlesMC implemented in Comonicon.
             sched = build_schedule(steps, burn, interval)
         end
         if alg == "StoreCallbacks"
-            callbacks = map(c -> eval(Meta.parse("callback_$c")), callbacks)
+            callbacks = map(c -> eval(Meta.parse("$c")), callbacks)
             algorithm = (
                 algorithm = eval(Meta.parse(alg)),
                 callbacks = callbacks,
                 scheduler = sched,
             )
+        elseif alg == "StoreAcceptance"
+            dependencies = map(d -> eval(Meta.parse("$d")), dependencies)
+            algorithm = (
+                algorithm = eval(Meta.parse(alg)),
+                dependencies = dependencies,
+                scheduler = sched,
+            )
         elseif alg == "StoreTrajectories" || alg == "StoreLastFrames"
              algorithm = (
                 algorithm = eval(Meta.parse(alg)),

src/molecules.jl

@@ -240,3 +240,7 @@ function compute_chain_correlation(system::Molecules)
     end
     return sum(correlation_array.^2)
 end
+
+Arianna.@callback function chain_correlation(system::Molecules)
+    return compute_chain_correlation(system)
+end

src/utils.jl

@@ -48,12 +48,8 @@ function SpeciesList(species)
     return SpeciesList(ids, heads)
 end
 
-function callback_energy(simulation)
-    return mean(system.energy[1] / length(system) for system in simulation.chains)
-end
-
-function callback_chain_correlation(simulation)
-    return mean(compute_chain_correlation(system) for system in simulation.chains)
+Arianna.@callback function energy(system::Particles)
+    return system.energy[1]
 end
 
 function volume_sphere(r, d::Int)

test/gerhard_energy_distribution.jl

@@ -29,7 +29,6 @@ steps = 1000
 burn = 0
 block = [0, 10]
 sampletimes = build_schedule(steps, burn, block)
-callbacks = (callback_energy, callback_acceptance)
 
 # NO SWAPS
 pswap = 0.0
@@ -41,7 +40,8 @@ pool = (
 
 algorithm_list = (
     (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=system.N),
-    (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes),
+    (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+    (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes),
     (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=EXYZ()),
     (algorithm=StoreLastFrames, scheduler=[steps], fmt=EXYZ()),
     (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10)),
@@ -72,7 +72,8 @@ pool = (
 )
 algorithm_list = (
     (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=system.N),
-    (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes),
+    (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+    (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes),
     (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=LAMMPS()),
     (algorithm=StoreLastFrames, scheduler=[steps], fmt=EXYZ()),
     (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10)),
@@ -105,7 +106,8 @@ pool = (
 )
 algorithm_list = (
     (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=system.N),
-    (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes),
+    (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+    (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes),
     (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=LAMMPS()),
     (algorithm=StoreLastFrames, scheduler=[steps], fmt=LAMMPS()),
     (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10)),

test/molecules_test.jl

@@ -43,8 +43,8 @@ callbacks = (callback_energy, callback_acceptance)
 path = "data/test/particles/Molecules/T$temperature/N$N/M$M/seed$seed"
 algorithm_list = (
     (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=N),
-    (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes),
-    (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=XYZ()),
+    (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+    (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes),    (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=XYZ()),
     (algorithm=StoreLastFrames, scheduler=[steps], fmt=XYZ()),
     (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10), fmt=XYZ())
     )

test/params.toml

@@ -20,9 +20,14 @@ parameters = {sigma = 0.05}
 
 [[simulation.output]]
 algorithm = "StoreCallbacks"
-callbacks = ["energy", "acceptance"]
+callbacks = ["energy"]
 scheduler_params = {linear_interval = 100}
 
+[[simulation.output]]
+algorithm = "StoreAcceptance"
+dependencies = ["Metropolis"]
+scheduler_params = {linear_interval = 500}
+
 [[simulation.output]]
 algorithm = "StoreTrajectories"
 scheduler_params = {linear_interval = 100}

test/runtests.jl

@@ -46,7 +46,6 @@ end
     burn = 0
     block = [0, 1, 2, 4, 8]
     sampletimes = build_schedule(steps, burn, block)
-    callbacks = (callback_energy, callback_acceptance)
 
     # NO SWAPS
     pswap = 0.0
@@ -57,7 +56,8 @@ end
     )
     algorithm_list = (
     (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=system_el.N),
-    (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes),
+    (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+    (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes), 
     (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=EXYZ()),
     (algorithm=StoreLastFrames, scheduler=[steps], fmt=LAMMPS()),
     (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10)),
@@ -81,9 +81,9 @@ end
     run!(simulation)
 
     ## Read energy data and compare
-    path_energy_el = joinpath(path_el, "energy.dat")
-    path_energy_ll = joinpath(path_ll, "energy.dat")
-    path_energy_vl = joinpath(path_vl, "energy.dat")
+    path_energy_el = joinpath(path_el, "chains/1/energy.dat")
+    path_energy_ll = joinpath(path_ll, "chains/1/energy.dat")
+    path_energy_vl = joinpath(path_vl, "chains/1/energy.dat")
     energy_el= readdlm(path_energy_el)[:, 2]
     energy_ll = readdlm(path_energy_ll)[:, 2]
     energy_vl = readdlm(path_energy_vl)[:, 2]
@@ -103,7 +103,8 @@ end
     )
     algorithm_list = (
         (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=system_el.N),
-        (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes, fmt=XYZ()),
+        (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+        (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes),
         (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=XYZ()),
         (algorithm=StoreLastFrames, scheduler=[steps], fmt=XYZ()),
         (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10), fmt=XYZ()),
@@ -121,8 +122,8 @@ end
     run!(simulation)
 
     ## Read energy data and compare
-    path_energy_el = joinpath(path_el, "energy.dat")
-    path_energy_ll = joinpath(path_ll, "energy.dat")
+    path_energy_el = joinpath(path_el, "chains/1/energy.dat")
+    path_energy_ll = joinpath(path_ll, "chains/1/energy.dat")
     energy_el= readdlm(path_energy_el)[:, 2]
     energy_ll = readdlm(path_energy_ll)[:, 2]
     @test isapprox(energy_el, energy_ll, atol=1e-6)
@@ -154,7 +155,6 @@ end
     burn = 0
     block = [0, 1, 2, 4, 8]
     sampletimes = build_schedule(steps, burn, block)
-    callbacks = (callback_energy, callback_acceptance)
 
     # NO SWAPS
     pswap = 0.0
@@ -165,7 +165,8 @@ end
     )
     algorithm_list = (
     (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=system_el.N),
-    (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes),
+    (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+    (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes),        
     (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=EXYZ()),
     (algorithm=StoreLastFrames, scheduler=[steps], fmt=EXYZ()),
     (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10)),
@@ -183,8 +184,8 @@ end
     run!(simulation)
 
     ## Read energy data and compare
-    path_energy_el = joinpath(path_el, "energy.dat")
-    path_energy_ll = joinpath(path_ll, "energy.dat")
+    path_energy_el = joinpath(path_el, "chains/1/energy.dat")
+    path_energy_ll = joinpath(path_ll, "chains/1/energy.dat")
     energy_el= readdlm(path_energy_el)[:, 2]
     energy_ll = readdlm(path_energy_ll)[:, 2]
     @test isapprox(energy_el, energy_ll, atol=1e-6)

test/simple_test.jl

@@ -49,8 +49,8 @@ callbacks = (callback_energy, callback_acceptance)
 
 algorithm_list = (
     (algorithm=Metropolis, pool=pool, seed=seed, parallel=false, sweepstep=length(system_ll)),
-    (algorithm=StoreCallbacks, callbacks=(callback_energy, callback_acceptance), scheduler=sampletimes),
-    (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=XYZ()),
+    (algorithm=StoreCallbacks, callbacks=(energy,), scheduler=sampletimes),
+    (algorithm=StoreAcceptance, dependencies=(Metropolis,), scheduler=sampletimes),    (algorithm=StoreTrajectories, scheduler=sampletimes, fmt=XYZ()),
     (algorithm=StoreLastFrames, scheduler=[steps], fmt=XYZ()),
     (algorithm=PrintTimeSteps, scheduler=build_schedule(steps, burn, steps ÷ 10)),
 )