Example 3: Calculating observables

This is an example calculation of the two-body correlation function $G_2$.

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

First, we load the reqired packages. Rimu for FCIQMC calculation, and DataFrames for maniplating the output.

using Rimu
using Random
using DataFrames

We use the same Hamiltonian as the first example, a Bose-Hubbard model with 6 particles in 6 sites, with strong interactions (we expect a Mott insulating state).

m = n = 6
initial_address = near_uniform(BoseFS{n,m})
H = HubbardReal1D(initial_address; u = 6.0, t = 1.0)

HubbardReal1D(BoseFS{6,6}(1, 1, 1, 1, 1, 1); u=6.0, t=1.0)

Now, we define the operators for the observables we wish to calculate.

dvals = 0:m-1
G2list = ((G2RealCorrelator(d) for d in dvals)...,)

(G2RealCorrelator(0), G2RealCorrelator(1), G2RealCorrelator(2), G2RealCorrelator(3), G2RealCorrelator(4), G2RealCorrelator(5))

This is a tuple of G2RealCorrelators, subtypes of AbstractHamiltonian. It calculates the density-density correlation function on a lattice

\[ \hat{G}^{(2)}(d) = \frac{1}{M} \sum_i^M \hat{n}_i (\hat{n}_{i+d} - \delta_{0d}).\]

with normalisation

\[ \sum_{d=0}^{M-1} \langle \hat{G}^{(2)}(d) \rangle = \frac{N (N-1)}{M}.\]

Observables that are defined by expectation values are calculated using the "replica trick". Thereby several independent copies or "replicas" of the state vector are propagated simultaneously. The reason is to have two (or more) stochastically independent copies of the state vector available such that we can calculate bias-free overlaps. We enable this by defining a ReplicaStrategy. Each replica has its own state and FCIQMC is effectively performed independently on each one. For calculating observables, we use AllOverlaps for the ReplicaStrategy. At each timestep, after the FCIQMC step is performed on, this strategy calculates the overlaps of every operator with the wavefunctions from each pair of replicas.

To obtain an unbiased result, at least two replicas should be used. One can also use more than two to improve the statistics. This is particularly helpful when the walker number is low.

num_replicas = 3
replica = AllOverlaps(num_replicas; operator = G2list)

AllOverlaps{3, 6, Tuple{G2RealCorrelator{0}, G2RealCorrelator{1}, G2RealCorrelator{2}, G2RealCorrelator{3}, G2RealCorrelator{4}, G2RealCorrelator{5}}, true}((G2RealCorrelator(0), G2RealCorrelator(1), G2RealCorrelator(2), G2RealCorrelator(3), G2RealCorrelator(4), G2RealCorrelator(5)))

Other FCIQMC parameters and strategies can be set in the same way as before.

steps_equilibrate = 1_000
steps_measure = 5_000
targetwalkers = 100;
dτ = 0.001

Now, we run FCIQMC. Note that passing an initial vector is optional - if we only pass the style, a vector with the appropriate style is created automatically.

df, state = lomc!(
    H; style=IsDynamicSemistochastic(),
    dτ,
    laststep = steps_equilibrate + steps_measure,
    targetwalkers,
    replica,
);

The output DataFrame has FCIQMC statistics for each replica (e.g. shift, norm),

println(filter(startswith("shift_"), names(df)))

["shift_1", "shift_2", "shift_3"]

as well as vector-vector overlaps (e.g. c1_dot_c2),

println(filter(contains("dot"), names(df)))

["c1_dot_c2", "c1_dot_c3", "c2_dot_c3"]

and operator overlaps (e.g. c1_Op1_c2) between the replicas.

println(filter(contains("Op"), names(df)))

["c1_Op1_c2", "c1_Op2_c2", "c1_Op3_c2", "c1_Op4_c2", "c1_Op5_c2", "c1_Op6_c2", "c1_Op1_c3", "c1_Op2_c3", "c1_Op3_c3", "c1_Op4_c3", "c1_Op5_c3", "c1_Op6_c3", "c2_Op1_c3", "c2_Op2_c3", "c2_Op3_c3", "c2_Op4_c3", "c2_Op5_c3", "c2_Op6_c3"]

The vector-vector and operator overlaps go into calculating the Rayleigh quotient for an observable

\[ \langle \hat{G}^{(2)}(d) \rangle = \frac{\sum_{a<b} \mathbf{c}_a^\dagger \cdot \hat{G}^{(2)}(d) \cdot \mathbf{c}_b}{\sum_{a<b} \mathbf{c}_a^\dagger \cdot \mathbf{c}_b }\]

The sum over all replica pairs (a,b), especially in the denominator, helps to avoid errors from poor sampling if the number of walkers is too low.

We use the function rayleigh_replica_estimator to calculate the Rayleigh quotient using all replicas in df, returning a RatioBlockingResult. Using the keyword skip will ignore the initial equilibration steps.

Now, we can calculate the correlation function for each value of $d$.

println("Two-body correlator from $num_replicas replicas:")
for d in dvals
    r = rayleigh_replica_estimator(df; op_name = "Op$(d+1)", skip=steps_equilibrate)
    println("   G2($d) = $(r.f) ± $(r.σ_f)")
end

Two-body correlator from 3 replicas:
   G2(0) = 0.210397630724869 ± 0.0029454115601256865
   G2(1) = 0.9179375998196901 ± 0.0014487882545570095
   G2(2) = 0.9824618365485462 ± 0.0008078855891766442
   G2(3) = 0.9888034965386578 ± 0.0012000192916385307
   G2(4) = 0.9824618365485462 ± 0.0008078855891766442
   G2(5) = 0.9179375998196901 ± 0.0014487882545570095

As expected, the onsite correlation at $d=0$ is low since this is a Mott insulating state with unit filling fraction, and is close to $1.0$ for all other values of the displacement $d$.

Since we ran multiple independent replicas, we also have multiple estimates of the shift energy.

println("Shift energy from $num_replicas replicas:")
for i in 1:num_replicas
    se = shift_estimator(df; shift="shift_$i", skip=steps_equilibrate)
    println("   Replica $i: $(se.mean) ± $(se.err)")
end

Shift energy from 3 replicas:
   Replica 1: -3.982824328882403 ± 0.1374885675286172
   Replica 2: -3.974753826250174 ± 0.12466355058904437
   Replica 3: -4.011054302602587 ± 0.12438161612595744

This page was generated using Literate.jl.