From efef3db8cdf449dd0f3eb4e9cd6654ab8b0de589 Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Mon, 4 Dec 2023 14:27:10 +0000 Subject: [PATCH 1/9] Attempted fix for #320 AdvancedMH or one of its dependencies may have updated a keyword from "initial_config" to "initial_params". This caused NQCD to not deliver initial positions to the sampling chain correctly. --- src/InitialConditions/ThermalMonteCarlo.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/InitialConditions/ThermalMonteCarlo.jl b/src/InitialConditions/ThermalMonteCarlo.jl index 00f837344..55b6a2810 100644 --- a/src/InitialConditions/ThermalMonteCarlo.jl +++ b/src/InitialConditions/ThermalMonteCarlo.jl @@ -27,7 +27,7 @@ using NQCDynamics: masses """ - run_advancedhmc_sampling(sim, r, steps, σ; move_ratio=0.0, internal_ratio=0.0) + run_advancedmh_sampling(sim, r, steps, σ; move_ratio=0.0, internal_ratio=0.0) Sample the configuration space for the simulation `sim` starting from `r`. @@ -47,7 +47,8 @@ function run_advancedmh_sampling( steps::Real, σ::Dict{Symbol,<:Real}; move_ratio=0.0, - internal_ratio=0.0 + internal_ratio=0.0, + kwargs... ) density = get_density_function(sim) @@ -60,7 +61,7 @@ function run_advancedmh_sampling( initial_config = reshape_input(sim, copy(r)) chain = AdvancedMH.sample(density_model, sampler, convert(Int, steps); - init_params=initial_config) + initial_params=initial_config, kwargs...) return reshape_output(sim, chain) end From f22e003409f8ca213559d6079a9e6d13e500cb6d Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Mon, 4 Dec 2023 15:32:15 +0000 Subject: [PATCH 2/9] Initial test for #321 --- src/InitialConditions/ThermalMonteCarlo.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/InitialConditions/ThermalMonteCarlo.jl b/src/InitialConditions/ThermalMonteCarlo.jl index 55b6a2810..d360c73dd 100644 --- a/src/InitialConditions/ThermalMonteCarlo.jl +++ b/src/InitialConditions/ThermalMonteCarlo.jl @@ -165,7 +165,7 @@ end function Base.rand(rng::Random.AbstractRNG, p::ClassicalProposal) result = map(x -> rand(rng, x), p.proposal) - result[Random.randsubseq(eachindex(result), p.move_ratio)] .= 0 + result[Random.randsubseq(findall(x -> x!=0.0, result), p.move_ratio)] .=0 return result end @@ -174,12 +174,12 @@ function Base.rand(rng::Random.AbstractRNG, p::RingPolymerProposal) reshaped_result = reshape(result, size(p.sim)) # Zero some of the centroid moves - sequence = Random.randsubseq(CartesianIndices(size(p.sim)[1:2]), p.move_ratio) + sequence = Random.randsubseq(findall(x -> x != 0, CartesianIndices(size(p.sim)[1:2])), p.move_ratio) reshaped_result[sequence, 1] .= 0 # Zero some of the internal mode moves for i = 2:nbeads(p.sim) - sequence = Random.randsubseq(CartesianIndices((ndofs(p.sim), natoms(p.sim))), p.internal_ratio) + sequence = Random.randsubseq(findall(x -> x != 0, CartesianIndices((ndofs(p.sim), natoms(p.sim)))), p.internal_ratio) reshaped_result[sequence, i] .= 0 end From f16d792872578e98269ae03295823c56314eb078 Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Tue, 5 Dec 2023 11:49:16 +0000 Subject: [PATCH 3/9] Modify AdvancedMH.jl keyword handling #321 - Update ThermalMonteCarlo.jl with new keywords movement_ratio and stop_ratio, where movement_ratio=1-stop_ratio - Add warning for deprecated kwargs. - Update documentation with information about new keywords --- .../initialconditions/metropolishastings.md | 15 ++++- src/InitialConditions/ThermalMonteCarlo.jl | 63 +++++++++++++------ 2 files changed, 57 insertions(+), 21 deletions(-) diff --git a/docs/src/initialconditions/metropolishastings.md b/docs/src/initialconditions/metropolishastings.md index e188b4295..aabfe599c 100644 --- a/docs/src/initialconditions/metropolishastings.md +++ b/docs/src/initialconditions/metropolishastings.md @@ -33,14 +33,25 @@ steps = 1e4 step_size = Dict(:H=>1) ``` -Now we can run the sampling. The extra keyword argument `move_ratio` is used to specify +Now we can run the sampling. The extra keyword argument `movement_ratio` is used to specify the fraction of the system moved during each Monte Carlo step. If we attempt to move the entire system at once, we can expect a very low acceptance ratio, whereas is we move only a single atom, the sampling will take much longer. You will likely have to experiment with this parameter to achieve optimal sampling. + +!!! note + + Keyword arguments relating to how much of the system to sample were recently changed. The existing `move_ratio` and `internal_ratio` arguments are no longer used. + Instead, there are now two options to specify how much of the system to move: + `movement_ratio`: Defines which fraction of the system to move. 1 moves the entire system. + `stop_ratio`: Defines which fraction of the system *not* to move. 1 stops the entire system. + `movement_ratio_internal`: Defines which proportion of ring polymer normal modes to perturb. 1 moves the entire system. + `stop_ratio_internal`: Defines which proportion of ring polymer normal modes not to perturb. 1 stops the entire system. + + ```@example mh using NQCDynamics.InitialConditions: ThermalMonteCarlo -chain = ThermalMonteCarlo.run_advancedmh_sampling(sim, r0, steps, step_size; move_ratio=0.5) +chain = ThermalMonteCarlo.run_advancedmh_sampling(sim, r0, steps, step_size; movement_ratio=0.5) ``` Now that our sampling is complete we can evaluate the potential energy expectation value. diff --git a/src/InitialConditions/ThermalMonteCarlo.jl b/src/InitialConditions/ThermalMonteCarlo.jl index d360c73dd..71e5a4e5f 100644 --- a/src/InitialConditions/ThermalMonteCarlo.jl +++ b/src/InitialConditions/ThermalMonteCarlo.jl @@ -27,34 +27,59 @@ using NQCDynamics: masses """ - run_advancedmh_sampling(sim, r, steps, σ; move_ratio=0.0, internal_ratio=0.0) + run_advancedmh_sampling(sim, r, steps, σ; movement_ratio=nothing, movement_ratio_internal=nothing) Sample the configuration space for the simulation `sim` starting from `r`. Total number of steps is given by `steps` and `σ` is the dictionary of step sizes for each species. -`move_ratio` defaults to `0.0` and denotes the fraction of system moved each step. -If `move_ratio = 0`, every degree of freedom is moved at each step. -If `move_ratio = 1`, then nothing will happen. Experiment with this parameter to achieve -optimal sampling. +`movement_ratio` denotes the fraction of system moved each step. +`internal_ratio` works as for `movement_ratio` but for the internal modes of the ring polymer. +For `movement_ratio = 0`, every degree of freedom is moved at each step, if `movement_ratio = 1`, then nothing will happen. -`internal_ratio` works as for `move_ratio` but for the internal modes of the ring polymer. +If neither arguments are defined, default behaviour is to move one atom (and one ring polymer normal mode) per step on average. """ function run_advancedmh_sampling( sim::AbstractSimulation, r, steps::Real, σ::Dict{Symbol,<:Real}; - move_ratio=0.0, - internal_ratio=0.0, + movement_ratio=nothing, + stop_ratio=nothing, + movement_ratio_internal=nothing, + stop_ratio_internal=nothing, + move_ratio=nothing, # deprecated + internal_ratio=nothing, # deprecated kwargs... ) + # Give a warning if move_ratio or internal_ratio are used + if move_ratio !== nothing || internal_ratio !== nothing + @warn "move_ratio and internal_ratio kwargs are deprecated and may be removed in future. More information: https://nqcd.github.io/NQCDynamics.jl/stable/initialconditions/metropolishastings/" + stop_ratio=move_ratio === nothing ? nothing : move_ratio + stop_ratio_internal=internal_ratio === nothing ? nothing : internal_ratio + end + + # If no kwargs for system fraction to move are given, perturb one atom and one normal mode at a time + if movement_ratio===nothing && stop_ratio===nothing + @debug "No movement restriction for atoms provided, automatically setting to move one atom per step on average." + movement_ratio=1/size(sim)[2] + end + + if movement_ratio_internal===nothing && stop_ratio_internal===nothing + @debug "No movement restriction for ring polymer normal modes provided, automatically setting to move one mode per step on average." + movement_ratio_internal=length(size(sim))==3 ? 1/size(sim)[3] : 0.0 + end + + # Set atom movement ratio by using whichever keyword is defined + stop_ratio = stop_ratio===nothing ? 1-movement_ratio : stop_ratio + stop_ratio_internal = stop_ratio_internal===nothing ? 1-movement_ratio_internal : stop_ratio_internal + density = get_density_function(sim) density_model = AdvancedMH.DensityModel(density) - proposal = get_proposal(sim, σ, move_ratio, internal_ratio) + proposal = get_proposal(sim, σ, stop_ratio, stop_ratio_internal) sampler = AdvancedMH.MetropolisHastings(proposal) @@ -117,29 +142,29 @@ end struct ClassicalProposal{P,T,S<:Simulation} <: AdvancedMH.Proposal{P} proposal::P - move_ratio::T + stop_ratio::T sim::S end struct RingPolymerProposal{P,T,S<:RingPolymerSimulation} <: AdvancedMH.Proposal{P} proposal::P - move_ratio::T - internal_ratio::T + stop_ratio::T + stop_ratio_internal::T sim::S end const MolecularProposal{P,T,S} = Union{RingPolymerProposal{P,T,S},ClassicalProposal{P,T,S}} -function ClassicalProposal(sim::Simulation, σ, move_ratio) +function ClassicalProposal(sim::Simulation, σ, stop_ratio) proposals = Matrix{UnivariateDistribution}(undef, size(sim)) for (i, symbol) in enumerate(sim.atoms.types) # Position proposals distribution = σ[symbol] == 0 ? Dirac(0) : Normal(0, σ[symbol]) proposals[:, i] .= distribution end - ClassicalProposal(proposals[:], move_ratio, sim) + ClassicalProposal(proposals[:], stop_ratio, sim) end -function RingPolymerProposal(sim::RingPolymerSimulation, σ, move_ratio, internal_ratio) +function RingPolymerProposal(sim::RingPolymerSimulation, σ, stop_ratio, stop_ratio_internal) ωₖ = RingPolymers.get_matsubara_frequencies(nbeads(sim), sim.beads.ω_n) proposals = Array{UnivariateDistribution}(undef, size(sim)) for i = 1:nbeads(sim) @@ -160,12 +185,12 @@ function RingPolymerProposal(sim::RingPolymerSimulation, σ, move_ratio, interna end end end - RingPolymerProposal(proposals[:], move_ratio, internal_ratio, sim) + RingPolymerProposal(proposals[:], stop_ratio, stop_ratio_internal, sim) end function Base.rand(rng::Random.AbstractRNG, p::ClassicalProposal) result = map(x -> rand(rng, x), p.proposal) - result[Random.randsubseq(findall(x -> x!=0.0, result), p.move_ratio)] .=0 + result[Random.randsubseq(findall(x -> x!=0.0, result), p.stop_ratio)] .=0 return result end @@ -174,12 +199,12 @@ function Base.rand(rng::Random.AbstractRNG, p::RingPolymerProposal) reshaped_result = reshape(result, size(p.sim)) # Zero some of the centroid moves - sequence = Random.randsubseq(findall(x -> x != 0, CartesianIndices(size(p.sim)[1:2])), p.move_ratio) + sequence = Random.randsubseq(findall(x -> x != 0, CartesianIndices(size(p.sim)[1:2])), p.stop_ratio) reshaped_result[sequence, 1] .= 0 # Zero some of the internal mode moves for i = 2:nbeads(p.sim) - sequence = Random.randsubseq(findall(x -> x != 0, CartesianIndices((ndofs(p.sim), natoms(p.sim)))), p.internal_ratio) + sequence = Random.randsubseq(findall(x -> x != 0, CartesianIndices((ndofs(p.sim), natoms(p.sim)))), p.stop_ratio_internal) reshaped_result[sequence, i] .= 0 end From 043da11906daba1297d1ca8659526a484165f17b Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Tue, 5 Dec 2023 11:51:13 +0000 Subject: [PATCH 4/9] Drop compat for AdvancedMH.jl<0.8 for #320,#321 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index dca911953..ebb4ed394 100644 --- a/Project.toml +++ b/Project.toml @@ -46,7 +46,7 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" [compat] AdvancedHMC = "0.5, 0.6" -AdvancedMH = "0.6, 0.7, 0.8" +AdvancedMH = "0.8" ComponentArrays = "0.11, 0.12, 0.13, 0.14, 0.15" DEDataArrays = "0.2" Dictionaries = "0.3" From 06a49ba5e957aad4de4ae06e96c7afe46598a040 Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Mon, 8 Jan 2024 14:10:11 +0000 Subject: [PATCH 5/9] Improved docs --- docs/src/initialconditions/metropolishastings.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/src/initialconditions/metropolishastings.md b/docs/src/initialconditions/metropolishastings.md index aabfe599c..b7cb0faf7 100644 --- a/docs/src/initialconditions/metropolishastings.md +++ b/docs/src/initialconditions/metropolishastings.md @@ -42,10 +42,15 @@ You will likely have to experiment with this parameter to achieve optimal sampli !!! note Keyword arguments relating to how much of the system to sample were recently changed. The existing `move_ratio` and `internal_ratio` arguments are no longer used. + Instead, there are now two options to specify how much of the system to move: + `movement_ratio`: Defines which fraction of the system to move. 1 moves the entire system. + `stop_ratio`: Defines which fraction of the system *not* to move. 1 stops the entire system. + `movement_ratio_internal`: Defines which proportion of ring polymer normal modes to perturb. 1 moves the entire system. + `stop_ratio_internal`: Defines which proportion of ring polymer normal modes not to perturb. 1 stops the entire system. From 4bca8e11f5eca3b5a803671c2fd763bab3691791 Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Tue, 9 Jan 2024 09:43:52 +0000 Subject: [PATCH 6/9] Document kwarg passing to AdvancedMH.sample() --- .../initialconditions/metropolishastings.md | 58 +++++++++---------- src/InitialConditions/ThermalMonteCarlo.jl | 4 +- 2 files changed, 32 insertions(+), 30 deletions(-) diff --git a/docs/src/initialconditions/metropolishastings.md b/docs/src/initialconditions/metropolishastings.md index b7cb0faf7..04f3b8d40 100644 --- a/docs/src/initialconditions/metropolishastings.md +++ b/docs/src/initialconditions/metropolishastings.md @@ -14,7 +14,17 @@ random walk starting from an initial configuration. These are accepted or rejected based upon the Metropolis-Hastings criteria. The result is a Markov chain that samples the canonical distribution. -## Example +!!! Legacy version + + Prior to the use of [`AdvancedMH.jl`](https://github.com/TuringLang/AdvancedMH.jl), + an alternative version of the algorithm was implemented that works for both classical + and ring polymer systems: `MetropolisHastings.run_monte_carlo_sampling(sim, R0, Δ, passes)` + + This is currently still included in the code but should be regarded as deprecated and + will likely be removed/combined with the [`AdvancedMH.jl`](https://github.com/TuringLang/AdvancedMH.jl) + version. + +## Example 1 We can perform the sampling by setting up a classical simulation in the usual way and providing an appropriate initial configuration. @@ -69,17 +79,10 @@ Estimators.@estimate potential_energy(sim, chain) sim.temperature / 2 * 5 ``` -## Legacy version - -Prior to the use of [`AdvancedMH.jl`](https://github.com/TuringLang/AdvancedMH.jl), -an alternative version of the algorithm was implemented that works for both classical -and ring polymer systems. -This is currently still included in the code but should be regarded as deprecated and -will likely be removed/combined with the [`AdvancedMH.jl`](https://github.com/TuringLang/AdvancedMH.jl) -version. -Here, we use the legacy version to obtain a thermal distribution in a simple -model system. +## Example 2 +Here, we obtain a thermal distribution in a simple model system with some additional tweaks to +try and sample a larger configuration space. ```@setup monte using NQCDynamics @@ -99,40 +102,37 @@ nothing # hide ``` Then we have to specify the parameters for the Monte Carlo simulation and perform the sampling. -`Δ` contains the step sizes for each of the species, `R0` the initial geometry and `passes` the -number of monte carlo passes we perform (`passes*n_atoms` steps total). +`Δ` contains the step sizes for each of the species, `R0` the initial geometry and `samples` the +number of configurations we want to obtain. + +*AdvancedMH.jl* provides some additional options to control the sampling process. To (hopefully) include + a larger configuration space in the final results, we set a `thinning` of 10, meaning that we only keep + every 10th proposed configuration. In addition, we discard the first 100 samples, since our initial configuration might not +lie in the equilibrium. This is set with `discard_initial`. +Further explanations of the keyword arguments can be found in the *AbstractMCMC.jl* [documentation](https://turinglang.org/AbstractMCMC.jl/dev/api/#Common-keyword-arguments). + ```@example monte -Δ = Dict([(:N, 0.1), (:O, 0.1)]) +Δ = Dict(:N => 0.1, :O => 0.1) R0 = [1.0 0.0; 0.0 0.0; 0.0 0.0] -passes = 1000 -output = InitialConditions.MetropolisHastings.run_monte_carlo_sampling(sim, R0, Δ, passes) +samples = 1000 +output = ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, samples, Δ; movement_ratio=0.5, thinning=thinning, discard_initial=discard_initial) nothing # hide ``` -Output has three fields: the acceptance rates for each species and the energies and geometries -obtained during sampling. -```@repl monte -output.acceptance -``` -```@example monte -plot(output.energy) -xlabel!("Step") # hide -ylabel!("Energy") # hide -``` - We can calculate the distance between each atom and plot the bond length throughout the sampling. ```@example monte using LinearAlgebra -plot([norm(R[:,1] .- R[:,2]) for R in output.R]) +plot([norm(R[:,1] .- R[:,2]) for R in output]) xlabel!("Step") # hide ylabel!("Bond length") # hide ``` The result of this simulation seamlessly interfaces with the `DynamicalDistribution` -presented in the previous section and `output.R` can be readily passed to provide +presented in the previous section and `output` can be readily passed to provide the position distribution. The Monte Carlo sampling does not include velocities but these can be readily obtained from the Maxwell-Boltzmann distribution. + ```@setup logging runtime = round(time() - start_time; digits=2) @info "...done after $runtime s." diff --git a/src/InitialConditions/ThermalMonteCarlo.jl b/src/InitialConditions/ThermalMonteCarlo.jl index 71e5a4e5f..8101e8e24 100644 --- a/src/InitialConditions/ThermalMonteCarlo.jl +++ b/src/InitialConditions/ThermalMonteCarlo.jl @@ -27,7 +27,7 @@ using NQCDynamics: masses """ - run_advancedmh_sampling(sim, r, steps, σ; movement_ratio=nothing, movement_ratio_internal=nothing) + run_advancedmh_sampling(sim, r, steps, σ; movement_ratio=nothing, movement_ratio_internal=nothing, kwargs...) Sample the configuration space for the simulation `sim` starting from `r`. @@ -39,6 +39,8 @@ step sizes for each species. For `movement_ratio = 0`, every degree of freedom is moved at each step, if `movement_ratio = 1`, then nothing will happen. If neither arguments are defined, default behaviour is to move one atom (and one ring polymer normal mode) per step on average. + +Further kwargs are passed to `AdvancedMH.sample` to allow for [extra functionality](https://turinglang.org/AbstractMCMC.jl/dev/api/#Common-keyword-arguments). """ function run_advancedmh_sampling( sim::AbstractSimulation, From 017c61767cd7e4f1f37e0e8b1eba887067ba1f57 Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Tue, 9 Jan 2024 13:00:54 +0000 Subject: [PATCH 7/9] Docs example had missing variables --- docs/src/initialconditions/metropolishastings.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/initialconditions/metropolishastings.md b/docs/src/initialconditions/metropolishastings.md index 04f3b8d40..512a59160 100644 --- a/docs/src/initialconditions/metropolishastings.md +++ b/docs/src/initialconditions/metropolishastings.md @@ -115,7 +115,7 @@ Further explanations of the keyword arguments can be found in the *AbstractMCMC. Δ = Dict(:N => 0.1, :O => 0.1) R0 = [1.0 0.0; 0.0 0.0; 0.0 0.0] samples = 1000 -output = ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, samples, Δ; movement_ratio=0.5, thinning=thinning, discard_initial=discard_initial) +output = ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, samples, Δ; movement_ratio=0.5, thinning=10, discard_initial=100) nothing # hide ``` From 98259ab62404e68db296b7c8621a56d54529b5de Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Mon, 15 Jan 2024 10:53:46 +0000 Subject: [PATCH 8/9] Fixed mistake in docs --- docs/src/initialconditions/metropolishastings.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/initialconditions/metropolishastings.md b/docs/src/initialconditions/metropolishastings.md index 512a59160..3e848fee7 100644 --- a/docs/src/initialconditions/metropolishastings.md +++ b/docs/src/initialconditions/metropolishastings.md @@ -115,7 +115,7 @@ Further explanations of the keyword arguments can be found in the *AbstractMCMC. Δ = Dict(:N => 0.1, :O => 0.1) R0 = [1.0 0.0; 0.0 0.0; 0.0 0.0] samples = 1000 -output = ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, samples, Δ; movement_ratio=0.5, thinning=10, discard_initial=100) +output = InitialConditions.ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, samples, Δ; movement_ratio=0.5, thinning=10, discard_initial=100) nothing # hide ``` From 54ab1e6bbd79d0d75e409bc16c91c41e75c2deb1 Mon Sep 17 00:00:00 2001 From: Alexander Spears <39826690+Alexsp32@users.noreply.github.com> Date: Mon, 15 Jan 2024 17:01:26 +0000 Subject: [PATCH 9/9] Version bump for #325 and #329 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ebb4ed394..d6b6b70c4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NQCDynamics" uuid = "36248dfb-79eb-4f4d-ab9c-e29ea5f33e14" authors = ["James "] -version = "0.13.4" +version = "0.13.5" [deps] AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d"