Example 1: 1D Bose-Hubbard Model

This is an example calculation finding the ground state of a 1D Bose-Hubbard chain with 6 particles in 6 lattice sites.

A runnable script for this example is located here. Run it with julia BHM-example.jl.

Firstly, we load all needed modules. Rimu for FCIQMC calculation;

using Rimu
using Plots # for plotting

Setting up the model

Now we define the physical problem: Generating a configuration where 6 particles are evenly distributed in 6 lattice sites:

aIni = near_uniform(BoseFS{6,6})

BoseFS{6,6}(1, 1, 1, 1, 1, 1)

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)

HubbardReal1D(BoseFS{6,6}(1, 1, 1, 1, 1, 1); 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())

DVec{BoseFS{6, 6, BitString{11, 1, UInt16}},Int64} with 1 entry, style = IsStochasticInteger{Int64}()
  fs"|1 1 1 1 1 1⟩" => 10

The choice of the style argument already determines the FCIQMC algorithm to use (while it is irrelevant for the projected energy).

Observables are passed into the lomc! function with the post_step keyword argument.

post_step = ProjectedEnergy(Ĥ, svec)

ProjectedEnergy{HubbardReal1D{Float64, BoseFS{6, 6, BitString{11, 1, UInt16}}, 6.0, 1.0}, Rimu.DictVectors.FrozenDVec{BoseFS{6, 6, BitString{11, 1, UInt16}}, Int64}, Rimu.DictVectors.FrozenDVec{BoseFS{6, 6, BitString{11, 1, UInt16}}, Float64}}(:vproj, :hproj, HubbardReal1D(BoseFS{6,6}(1, 1, 1, 1, 1, 1); u=6.0, t=1.0), Rimu.DictVectors.FrozenDVec{BoseFS{6, 6, BitString{11, 1, UInt16}}, Int64}(Pair{BoseFS{6, 6, BitString{11, 1, UInt16}}, Int64}[fs"|1 1 1 1 1 1⟩" => 10]), Rimu.DictVectors.FrozenDVec{BoseFS{6, 6, BitString{11, 1, UInt16}}, Float64}(Pair{BoseFS{6, 6, BitString{11, 1, UInt16}}, Float64}[fs"|1 1 1 1 2 0⟩" => -14.142135623730951, fs"|0 2 1 1 1 1⟩" => -14.142135623730951, fs"|1 1 1 1 0 2⟩" => -14.142135623730951, fs"|1 2 0 1 1 1⟩" => -14.142135623730951, fs"|2 0 1 1 1 1⟩" => -14.142135623730951, fs"|1 1 1 2 0 1⟩" => -14.142135623730951, fs"|1 1 2 0 1 1⟩" => -14.142135623730951, fs"|1 1 0 2 1 1⟩" => -14.142135623730951, fs"|1 1 1 0 2 1⟩" => -14.142135623730951, fs"|1 0 2 1 1 1⟩" => -14.142135623730951, fs"|2 1 1 1 1 0⟩" => -14.142135623730951, fs"|0 1 1 1 1 2⟩" => -14.142135623730951]))

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;
            laststep = steps_equilibrate + steps_measure, # total number of steps
            dτ,
            targetwalkers,
            post_step,
);

df is a DataFrame containing the time series data.

Analysing the results

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.

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 with blocking transformations:

se = shift_estimator(df; skip=steps_equilibrate)

BlockingResult{Float64}
  mean = -3.815 ± 0.099
  with uncertainty of ± 0.008970614342568576
  from 62 blocks after 5 transformations (k = 6).

For the projected energy, it a bit more complicated as it's a ratio of fluctuationg quantities:

pe = projected_energy(df; skip=steps_equilibrate)

RatioBlockingResult{Float64,MonteCarloMeasurements.Particles{Float64, 2000}}
  ratio = -4.1774 ± (0.0474639, 0.051944) (MC)
  95% confidence interval: [-4.28066, -4.07921]) (MC)
  linear error propagation: -4.1726 ± 0.049528
  |δ_y| = |0.00975764| (≤ 0.1 for normal approx)
  Blocking successful with 31 blocks after 6 transformations (k = 7).

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)

(val = -4.177400207765482, val_l = 0.1032568959464415, val_u = 0.09819112324625667)

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(Ĥ)

462

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]

-4.021502406906392

Read more about Rimu.jls capabilities for exact diagonalisation in the example "Exact diagonalisation".

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))
         Exact Energy: $exact_energy")

Energy from 2000 steps with 1000 walkers:
         Shift: -3.814556830153461 ± 0.09908367395962789
         Projected Energy: -4.177400207765482 ± (0.1032568959464415, 0.09819112324625667)
         Exact Energy: -4.021502406906392

This page was generated using Literate.jl.