From a15e829c9e554d0eb67494c3793277c4bf700b79 Mon Sep 17 00:00:00 2001 From: Joachim Brand Date: Sat, 25 Nov 2023 23:14:59 +1300 Subject: [PATCH] Rework BHM-example.jl, new function `default_starting_vector` --- docs/make.jl | 8 +- scripts/BHM-example.jl | 160 +++++++++--------- src/Rimu.jl | 1 + src/lomc.jl | 50 ++++-- src/strategies_and_params/poststepstrategy.jl | 6 +- 5 files changed, 134 insertions(+), 91 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 6613cbe79..621c87de6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -44,7 +44,11 @@ end makedocs(; modules=[Rimu,Rimu.RimuIO], - format=Documenter.HTML(prettyurls = false), + format=Documenter.HTML( + prettyurls = false, + size_threshold=500_000, # 500 kB + size_threshold_warn=200_000, # 200 kB + ), pages=[ "Guide" => "index.md", "Examples" => EXAMPLES_PAIRS[sortperm(EXAMPLES_NUMS)], @@ -69,7 +73,7 @@ makedocs(; sitename="Rimu.jl", authors="Joachim Brand ", checkdocs=:exports, - doctest=false # Doctests are done while testing. + doctest=false, # Doctests are done while testing. ) deploydocs( diff --git a/scripts/BHM-example.jl b/scripts/BHM-example.jl index a99b13dff..b5858e22d 100644 --- a/scripts/BHM-example.jl +++ b/scripts/BHM-example.jl @@ -11,103 +11,109 @@ # `Rimu` for FCIQMC calculation; using Rimu -using Random +using Plots # for plotting + +# ## Setting up the model # Now we define the physical problem: -# Setting the number of lattice sites `m = 6`; -# and the number of particles `n = 6`: -m = n = 6 -# Generating a configuration that particles are evenly distributed: -aIni = near_uniform(BoseFS{n,m}) -# where `BoseFS` is used to create a bosonic system. -# The Hamiltonian is defined based on the configuration `aIni`, -# with additional onsite interaction strength `u = 6.0` -# and the hopping strength `t = 1.0`: +# Generating a configuration where 6 particles are evenly distributed in 6 lattice sites: +aIni = near_uniform(BoseFS{6,6}) +# where `BoseFS` is used to create a bosonic Fock state. +# The Hamiltonian is defined by specifying the model and the parameters. +# Here we use the Bose Hubbard model in one dimension and real space: Ĥ = HubbardReal1D(aIni; u = 6.0, t = 1.0) +# ## Parameters of the calculation + +# Now let's setup the Monte Carlo calculation. +# We need to decide the number of walkers to use in this Monte Carlo run, which is +# equivalent to the average one-norm of the coefficient vector: +targetwalkers = 1_000; + +# It is good practice to equilibrate the time series before taking statistics. +steps_equilibrate = 1_000; +steps_measure = 2_000; + +# The appropriate size of the time step is problem dependent. +dτ = 0.001; + +# ## Defining an observable + +# Now let's set up an observable to measure. Here we will measure the projected energy. +# In additon to the `shift`, the projected energy is a second estimator for the energy. + +# We first need to define a projector. Here we use the function `default_starting_vector` +# to generate a vector with only a single occupied configuration, the same vector that we +# will use as starting vector for the FCIQMC calculation. +svec = default_starting_vector(aIni; style=IsStochasticInteger()) +# The choice of the `style` argument already determines the FCIQMC algorithm to use +# (while it is irrelevant for the projected energy). -# Now let's setup the Monte Carlo settings. -# The number of walkers to use in this Monte Carlo run: -targetwalkers = 1_000 -# The number of time steps before doing statistics, -# i.e. letting the walkers to sample Hilbert and to equilibrate: -steps_equilibrate = 1_000 -# And the number of time steps used for getting statistics, -# e.g. time-average of shift, projected energy, walker numbers, etc.: -steps_measure = 2_000 - - -# Set the size of a time step -dτ = 0.001 -# and we report QMC data every k-th step, -# set the interval to record QMC data: -reporting_interval = 1 - -# Now we prepare initial state and allocate memory. -# The initial address is defined above as `aIni = near_uniform(Ĥ)`. -# Putting one of walkers into the initial address `aIni` -svec = DVec(aIni => 1) -# Let's plant a seed for the random number generator to get consistent result: -Random.seed!(17) - -# Now let's setup all the FCIQMC strategies. - -# Passing dτ and total number of time steps into params: -params = RunTillLastStep(dτ = dτ, laststep = steps_equilibrate + steps_measure) -# Strategy for updating the shift: -s_strat = DoubleLogUpdate(targetwalkers = targetwalkers, ζ = 0.08) -# Strategy for reporting info: -r_strat = ReportDFAndInfo(reporting_interval = reporting_interval, info_interval = 100) -# Strategy for updating dτ: -τ_strat = ConstantTimeStep() -# set up the calculation and reporting of the projected energy -# in this case we are projecting onto the starting vector, -# which contains a single configuration -post_step = ProjectedEnergy(Ĥ, copy(svec)) - -# Print out info about what we are doing: -println("Finding ground state for:") -println(Ĥ) -println("Strategies for run:") -println(params, s_strat) -println(τ_strat) +# Observables are passed into the `lomc!` function with the `post_step` keyword argument. +post_step = ProjectedEnergy(Ĥ, svec) + +# ## Running the calculation + +# Seeding the random number generator is sometimes useful in order to get reproducible +# results +using Random +Random.seed!(17); # Finally, we can start the main FCIQMC loop: -df, state = lomc!(Ĥ,svec; - params, - laststep = steps_equilibrate + steps_measure, - s_strat, - r_strat, - τ_strat, +df, state = lomc!(Ĥ, svec; + laststep = steps_equilibrate + steps_measure, # total number of steps + dτ, + targetwalkers, post_step, ); -# Here is how to save the output data stored in `df` into a `.arrow` file, -# which can be read in later: -println("Writing data to disk...") -save_df("fciqmcdata.arrow", df) +# `df` is a `DataFrame` containing the time series data. + +# ## Analysing the results -# Now let's look at the calculated energy from the shift. -# Loading the equilibrated data -qmcdata = last(df,steps_measure); +# We can plot the norm of the coefficient vector as a function of the number of steps: +hline([targetwalkers], label="targetwalkers", color=:red, linestyle=:dash) +plot!(df.steps, df.norm, label="norm", ylabel="norm", xlabel="steps") +# After some equilibriation steps, the norm fluctuates around the target number of walkers. -# compute the average shift and its standard error -se = shift_estimator(qmcdata) +# Now let's look at estimating the energy from the shift. +# The mean of the shift is a useful estimator of the shift. Calculating the error bars +# is a bit more involved as correlations have to be removed from the time series. +# The following code does that: +se = shift_estimator(df; skip=steps_equilibrate) # For the projected energy, it a bit more complicated as it's a ratio of two means: -pe = projected_energy(qmcdata) +pe = projected_energy(df; skip=steps_equilibrate) # The result is a ratio distribution. Let's get its median and lower and upper error bars # for a 95% confidence interval v = val_and_errs(pe; p=0.95) +# Let's visualise these estimators together with the time series of the shift +plot(df.steps, df.shift, ylabel="energy", xlabel="steps", label="shift") + +plot!(x->se.mean, df.steps[steps_equilibrate+1:end], ribbon=se.err, label="shift mean") +plot!( + x -> v.val, df.steps[steps_equilibrate+1:end], ribbon=(v.val_l,v.val_u), + label="projected_energy", +) +# In this case the projected energy and the shift are close to each other an the error bars +# are hard to see on this scale. + +# The problem was just a toy example, as the dimension of the Hamiltonian is rather small: +dimension(Ĥ) + +# In this case is easy (and more efficient) to calculate the exact ground state energy +# using standard linear algebra: +using LinearAlgebra +exact_energy = eigvals(Matrix(Ĥ))[1] + +# Comparing our results for the energy: println("Energy from $steps_measure steps with $targetwalkers walkers: - Shift: $(se.mean) ± $(se.err); - Projected Energy: $(v.val) ± ($(v.val_l), $(v.val_u))") + Shift: $(se.mean) ± $(se.err) + Projected Energy: $(v.val) ± ($(v.val_l), $(v.val_u)) + Exact Energy: $exact_energy") -# Finished! using Test #hide -@test isfile("fciqmcdata.arrow") #hide -@test se.mean ≈ -4.0215 rtol=0.1 #hide -rm("fciqmcdata.arrow", force=true) #hide +@test se.mean ≈ -4.0215 rtol=0.1; #hide diff --git a/src/Rimu.jl b/src/Rimu.jl index 77019bd77..856756e34 100644 --- a/src/Rimu.jl +++ b/src/Rimu.jl @@ -49,6 +49,7 @@ include("StatsTools/StatsTools.jl") @reexport using .StatsTools export lomc! +export default_starting_vector export FciqmcRunStrategy, RunTillLastStep export ShiftStrategy, LogUpdate, LogUpdateAfterTargetWalkers export DontUpdate, DoubleLogUpdate, DoubleLogUpdateAfterTargetWalkers diff --git a/src/lomc.jl b/src/lomc.jl index f33f590c8..922a37c9a 100644 --- a/src/lomc.jl +++ b/src/lomc.jl @@ -279,6 +279,43 @@ function report_default_metadata!(report::Report, state::QMCState) return report end +""" + default_starting_vector(hamiltonian::AbstractHamiltonian; kwargs...) + default_starting_vector( + address=starting_address(hamiltonian); + style=IsStochasticInteger(), + threading=nothing + ) +Return a default starting vector for [`lomc!`](@ref). The default choice for the starting +vector is +```julia +v = PDVec(starting_address => 10; style) +``` +if threading is available or + +```julia +v = DVec(starting_address => 10; style) +``` +otherwise. See [`PDVec`](@ref), [`DVec`](@ref) and [`StochasticStyle`](@ref). +""" +function default_starting_vector(hamiltonian::AbstractHamiltonian; kwargs...) + return default_starting_vector(starting_address(hamiltonian); kwargs...) +end +function default_starting_vector(address; + style=IsStochasticInteger(), + threading=nothing, +) + if isnothing(threading) + threading = Threads.nthreads() > 1 + end + if threading + v = PDVec(address => 10; style) + else + v = DVec(address => 10; style) + end + return v +end + """ lomc!(ham::AbstractHamiltonian, [v]; kwargs...) -> df, state lomc!(state::QMCState, [df]; kwargs...) -> df, state @@ -363,16 +400,8 @@ julia> metadata(df2, "hamiltonian") # some metadata is automatically added be passed for merging with the report `df`. The default choice for the starting vector is -```julia -v = PDVec(starting_address => 10; style) -``` -if threading is available or - -```julia -v = DVec(starting_address => 10; style) -``` -otherwise. It triggers the integer walker FCIQMC algorithm. See [`PDVec`](@ref), -[`DVec`](@ref) and [`StochasticStyle`](@ref). +`v = default_starting_vector(; address, style, threading)`. +See [`default_starting_vector`](@ref). """ function lomc!(ham, v; df=DataFrame(), name="lomc!", metadata=nothing, kwargs...) state = QMCState(ham, v; kwargs...) @@ -385,6 +414,7 @@ function lomc!( address=starting_address(ham), kwargs... ) + v = default_starting_vector(ham; style, threading, address) if Threads.nthreads() > 1 && (threading ≠ false) || threading == true v = PDVec(address => 10; style) else diff --git a/src/strategies_and_params/poststepstrategy.jl b/src/strategies_and_params/poststepstrategy.jl index 4f293ddef..2dadc8d67 100644 --- a/src/strategies_and_params/poststepstrategy.jl +++ b/src/strategies_and_params/poststepstrategy.jl @@ -24,8 +24,10 @@ abstract type PostStepStrategy end """ post_step(::PostStepStrategy, ::ReplicaState) -> kvpairs -Compute statistics after FCIQMC step. Should return a tuple of `:key => value` pairs. -This function is only called every [`reporting_interval`](@ref) steps, as defined by the `ReportingStrategy`. +Compute statistics after FCIQMC step. Should return a tuple of `:key => value` pairs. +This function is only called every [`reporting_interval`](@ref) steps, as defined by the +`ReportingStrategy`. + See also [`PostStepStrategy`](@ref), [`ReportingStrategy`](@ref). """ post_step