OL

src/Ansatz.jl

@@ -5,15 +5,91 @@
     g::Float64
     G::Float64
     Og::Float64
+    OL::Float64
 end
 
+"""
+Fast Gutzwiller Factor update technique from Becca and Sorella 2017
+"""
+function fast_G_update(newconf::BitVector, oldconf::BitVector, g::Float64, n_mean::Float64)
+    L = length(newconf) ÷ 2
+    occup = newconf[1:L] + newconf[L+1:end]
+    diff = newconf - oldconf
+    @assert sum(abs.(diff)) == 2 "Only one electron move"
+    ratio = exp(-g / 2 * sum(@. diff * ((1 - n_mean)^2) - n_mean^2))
+    return ratio
+end
+
+"""
+Fast computing technique from Becca and Sorella 2017
+"""
+function fast_update(
+    U::AbstractMatrix,
+    Uinvs::AbstractMatrix,
+    newconf::BitVector,
+    oldconf::BitVector,
+)
+    diff = newconf - oldconf
+    Rl = findfirst(==(-1), diff) # the old position of the l-th electron
+    l = sum(oldconf[1:Rl]) # l-th electron 
+    k = findfirst(==(1), diff) # the new position of the l-th electron K = R_l'
+    ratio = sum(U[k, :] .* Uinvs[:, l])
+    return ratio
+end
+
+"""
+The observable OL = <x|H|Ψ_G>/<x|Ψ_G> 
+"""
+function getOL(orb::AHmodel, conf_up::BitVector, conf_down::BitVector, g)
+    conf = vcat(conf_up, conf_down)
+    L = length(conf) ÷ 2
+    n_mean = (orb.N_up + orb.N_down) / orb.lattice.ns
+    OL = 0.0
+    U_upinvs = inv(orb.U_up[conf_up, :]) # might be optimized by column slicing
+    U_downinvs = inv(orb.U_down[conf_down, :])
+    xprime = getxprime(orb, LongBitStr(conf))
+    for confstr in keys(xprime)
+        new_conf = BitVector([confstr...])
+        coff = xprime[confstr]
+        if new_conf == conf
+            OL += coff
+        elseif new_conf[1:L] == conf_up
+            OL +=
+                coff *
+                fast_update(orb.U_down, U_downinvs, new_conf[L+1:end], conf_down) *
+                fast_G_update(new_conf[L+1:end], conf_down, g, n_mean)
+        elseif new_conf[L+1:end] == conf_down
+            OL +=
+                coff *
+                fast_update(orb.U_up, U_upinvs, new_conf[1:L], conf_up) *
+                fast_G_update(new_conf[1:L], conf_up, g, n_mean)
+        else
+            OL +=
+                coff *
+                fast_update(orb.U_up, U_upinvs, new_conf[1:L], conf_up) *
+                fast_update(orb.U_down, U_downinvs, new_conf[L+1:end], conf_down) *
+                fast_G_update(new_conf[1:L], conf_up, g, n_mean) *
+                fast_G_update(new_conf[L+1:end], conf_down, g, n_mean)
+        end
+
+    end
+    return OL
+end
+
+
 """
 add Gutzwiller Ansatz where G  = exp(-g/2 ∑_i (n_i - n_mean)^2), Ψ_G = G Ψ_0
 """
-function Gutzwiller(orbitals::AHmodel, conf_up::BitVector, conf_down::BitVector, g::Float64)
+function Gutzwiller(
+    orbitals::AHmodel{B},
+    conf_up::BitVector,
+    conf_down::BitVector,
+    g::Float64,
+) where {B}
     occupation = conf_up + conf_down
     n_mean = (orbitals.N_up + orbitals.N_down) / orbitals.lattice.ns
     Og = -1 / 2 * sum(@. (occupation - n_mean)^2)
     G = exp(g * Og)
-    return Gutzwiller(conf_up, conf_down, g, G, Og)
-end
+    OL = getOL(orbitals, conf_up, conf_down, g)
+    return Gutzwiller(conf_up, conf_down, g, G, Og, OL)
+end

src/Lattices.jl

@@ -26,7 +26,7 @@
     left = circshift(rectl, (-1, 0))
 
     # Store neighbor information in a nested vector
-    for i in 1:ns
+    @inbounds for i in 1:ns
         neigh[i] = [up[i], right[i], down[i], left[i]]
     end
     return LatticeRectangular{B}(nx, ny, ns, neigh)
@@ -36,7 +36,7 @@
     ns = nx * ny
     neigh = Vector{Vector{Int}}(undef, ns)
     # Loop through each site
-    for i in 1:ns
+    @inbounds for i in 1:ns
         neighbors = Int[]
 
         # Identify valid neighbors for each direction
@@ -59,6 +59,3 @@
     return LatticeRectangular{B}(nx, ny, ns, neigh)
 
 end
-
-
-

src/Orbitals.jl

@@ -30,14 +30,20 @@ struct AHmodel{B} <: AbstractOrbitals
     N_up::Int
     N_down::Int
     omega::Vector{Float64}
-    U_up::Matrix{ComplexF64}
-    U_down::Matrix{ComplexF64}
+    U_up::Matrix{Float64}
+    U_down::Matrix{Float64}
 end
 
 """
 Get the non-interacting Anderson model Hamiltonian Matrix to construct HF states
 """
-function getHmat(lattice::LatticeRectangular{B}, t::Float64, omega::Vector{Float64}, N_up::Int, N_down::Int) where {B}
+function getHmat(
+    lattice::LatticeRectangular{B},
+    t::Float64,
+    omega::Vector{Float64},
+    N_up::Int,
+    N_down::Int,
+) where {B}
     ns = lattice.ns
     N = N_up + N_down
     @assert ns / N == 1 "Should be hall-filling"
@@ -56,7 +62,14 @@ end
 """
 generate Anderson-Hubbard model and get the sampling ensemble
 """
-function AHmodel(lattice::LatticeRectangular{B}, t::Float64, W::Float64, U::Float64, N_up::Int, N_down::Int) where {B}
+function AHmodel(
+    lattice::LatticeRectangular{B},
+    t::Float64,
+    W::Float64,
+    U::Float64,
+    N_up::Int,
+    N_down::Int,
+) where {B}
     omega = randn(Float64, lattice.ns) * W / 2
     H_mat = getHmat(lattice, t, omega, N_up, N_down)
     # get sampling ensemble U_up and U_down
@@ -76,7 +89,7 @@ function getxprime(orb::AHmodel{B}, x::BitStr{N,T}) where {B,N,T}
     L = length(x) ÷ 2  # Int division
     xprime = Dict{typeof(x),Float64}()
     # consider the spin up case
-    @inbounds for i in 1:L
+    @inbounds for i = 1:L
         if readbit(x, i) == 1
             xprime[x] = get!(xprime, x, 0.0) + orb.omega[i] # On-site energy
             if readbit(x, i) == 1 && readbit(x, i + L) == 1 #occp[i] == 2
@@ -92,11 +105,11 @@ function getxprime(orb::AHmodel{B}, x::BitStr{N,T}) where {B,N,T}
             end
         end
     end
-    @inbounds for i in L+1:length(x)
+    @inbounds for i = L+1:length(x)
         if readbit(x, i) == 1
             xprime[x] = get!(xprime, x, 0.0) + orb.omega[i-L] # On-site energy
             for neigh in orb.lattice.neigh[i-L]
-                if readbit(x, neigh) == 0
+                if readbit(x, neigh + L) == 0
                     _x = x
                     _x &= ~indicator(T, i)
                     _x |= indicator(T, neigh + L)
@@ -106,4 +119,4 @@ function getxprime(orb::AHmodel{B}, x::BitStr{N,T}) where {B,N,T}
         end
     end
     return xprime
-end
+end

test/Ansatz.jl

@@ -1,4 +1,5 @@
 using Test, FastFermionSampling
+using LinearAlgebra
 
 @testset "Gutzwiller" begin
     g = 1.0
@@ -14,4 +15,16 @@ using Test, FastFermionSampling
     @test ansatz.conf_down == conf_down
     ansatz2 = Gutzwiller(orb, conf_up, conf_down, 0.0)
     @test ansatz2.G == 1.0
+end
+
+@testset "fast_update" begin
+    conf_up = BitVector([1, 0, 1, 0])
+    conf_down = BitVector([0, 1, 0, 1])
+    orb = AHmodel(LatticeRectangular(2, 2, Periodic()), 1.0, 1.0, 1.0, 2, 2)
+    U_upinvs = inv(orb.U_up[conf_up, :])
+    U_downinvs = inv(orb.U_down[conf_down, :])
+    new_conf_up = BitVector([0, 1, 1, 0])
+    new_conf_down = BitVector([1, 0, 0, 1])
+    @test FastFermionSampling.fast_update(orb.U_up, U_upinvs, new_conf_up, conf_up) ≈ det(orb.U_up[new_conf_up, :]) / det(orb.U_up[conf_up, :])
+    @test FastFermionSampling.fast_update(orb.U_down, U_downinvs, new_conf_down, conf_down) ≈ det(orb.U_down[new_conf_down, :]) / det(orb.U_down[conf_down, :])
 end

test/Orbitals.jl

@@ -16,19 +16,24 @@ end
     x = LongBitStr([true, true, false, false, false, true, true, false])
     xprime = getxprime(orb, x)
     occp = x[1:4] + x[5:end]
-    @test length(keys(xprime)) == 5
+    @test length(keys(xprime)) == 7 
     @test xprime[x] ≈ sum(occp .* orb.omega) + 1.0
-    @test xprime[LongBitStr([false, true, true, false, false, true, true, false])]  == -2
-    @test xprime[LongBitStr([true, false, false, true, false, true, true, false])]  == -2
-    @test xprime[LongBitStr([true, true, false, false, false, true, false, true])]  == -2
-    @test xprime[LongBitStr([true, true, false, false, false, false, true, true])]  == -2
+    @test xprime[LongBitStr([false, true, true, false, false, true, true, false])] == -2
+    @test xprime[LongBitStr([true, false, false, true, false, true, true, false])] == -2
+    @test xprime[LongBitStr([true, true, false, false, false, true, false, true])] == -2
+    @test xprime[LongBitStr([true, true, false, false, false, false, true, true])] == -2
+    @test xprime[LongBitStr([true, true, false, false, true, true, false, false])] == -2
+    @test xprime[LongBitStr([true, true, false, false, true, false, true, false])] == -2
+    for conf in keys(xprime)
+        @test sum([conf...]) == 4
+    end
 end
 
 using BenchmarkTools
 function bm(n::Int)
     @assert n % 2 == 0
     lat = LatticeRectangular(n, n, Periodic())
-    orb = AHmodel(lat, 1.0, 1.0, 1.0, n^2÷2, n^2÷2)
+    orb = AHmodel(lat, 1.0, 1.0, 1.0, n^2 ÷ 2, n^2 ÷ 2)
     #x = rand(Bool, 2*n^2)
     x = LongBitStr(rand(0:1, 2n^2))
     @btime getxprime($orb, $x)

test/runtests.jl

@@ -1,4 +1,3 @@
-using FastFermionSampling
 using Test
 
 @testset "FastFermionSampling.jl" begin