Add files via upload

FermiHubbard/Calculations.jl

@@ -0,0 +1,177 @@
+using LinearAlgebra
+
+
+function recursion_partition_function(num_elec, eff_spec_SpinUp, eff_spec_SpinDn)
+    Z_up, Z_dn = zeros(num_elec), zeros(num_elec)
+    zₖ_up, zₖ_dn = zeros(num_elec), zeros(num_elec)
+    for k = 1 : num_elec
+        zₖ_up[k] = sum(eff_spec_SpinUp .^ k)
+        zₖ_dn[k] = sum(eff_spec_SpinDn .^ k)
+    end
+    for N = 1 : num_elec
+        for k = 1 : N
+            sgn_k = 2 * isodd(k) - 1
+            if k == N
+                Z_up[N] += sgn_k * zₖ_up[k]
+                Z_dn[N] += sgn_k * zₖ_dn[k]
+            else
+                Z_up[N] += sgn_k * Z_up[N-k] * zₖ_up[k]
+                Z_dn[N] += sgn_k * Z_dn[N-k] * zₖ_dn[k]
+            end
+        end
+        Z_up[N] = Z_up[N] / N
+        Z_dn[N] = Z_dn[N] / N
+    end
+    return Z_up, Z_dn
+end
+
+function recursion_OccNum(Z_up, eff_spec_SpinUp, Z_dn, eff_spec_SpinDn, num_elec, num_sites) # Calculating <nᵢ> * Z_N
+    n_up_mean = zeros(num_elec, num_sites)
+    n_dn_mean = zeros(num_elec, num_sites)
+    n_up_mean[1, :], n_dn_mean[1, :] = eff_spec_SpinUp, eff_spec_SpinDn
+    for N = 2 : num_elec
+        n_up_mean[N, :] = eff_spec_SpinUp .* (Z_up[N-1] .- n_up_mean[N-1, :])
+        n_dn_mean[N, :] = eff_spec_SpinDn .* (Z_dn[N-1] .- n_dn_mean[N-1, :])
+    end
+    return n_up_mean ./ Z_up, n_dn_mean ./ Z_dn
+end
+
+function measure_tot_energy(num_sites, t, U, U_array_SpinUp, U_array_SpinDn, n_up_mean, n_dn_mean, num_SpinUp, num_SpinDn, Z_up, Z_dn)
+    E_onebody, E_twobody = 0, 0
+    for i = 1 : num_sites
+        for j = 1 : num_sites
+            if abs(i - j) == 1 || abs(i - j) == num_sites - 1
+                D_ij_up = sum(U_array_SpinUp[i, j, :] .* n_up_mean[num_SpinUp, :])
+                D_ij_dn = sum(U_array_SpinDn[i, j, :] .* n_dn_mean[num_SpinDn, :])
+                E_onebody += -t * ( D_ij_up + D_ij_dn )
+            end
+        end
+    end
+    for i = 1 : num_sites
+        U_matrix = U_array_SpinUp[i, i, :] * U_array_SpinDn[i, i, :]'
+        n_matrix = n_up_mean[num_SpinUp, :] * n_dn_mean[num_SpinDn, :]'
+        D_iiii = sum(U_matrix .* n_matrix)
+        E_twobody += U * D_iiii
+    end
+    return E_onebody, E_twobody
+end
+
+
+function recursion(self::Lattice, expϵup::Array{Float64,1}, expϵdn::Array{Float64,1})
+    nup_mean = zeros(Float64, self.nsites, self.nsites)
+    ndn_mean = zeros(Float64, self.nsites, self.nsites)
+    Znup_ratio = zeros(Float64, self.ntot)
+    Zndn_ratio = zeros(Float64, self.ntot)
+    expϵup = expϵup ./ expϵup[self.nsites]
+    expϵdn = expϵdn ./ expϵdn[self.nsites]
+    Znup_ratio[1], Zndn_ratio[1] = sum(expϵup), sum(expϵdn)
+    nup_mean[:, 1], ndn_mean[:, 1] = expϵup / Znup_ratio[1], expϵdn / Zndn_ratio[1]
+    for N = 2 : max(self.nup, self.ndn)
+        Znup_ratio[N] = sum(expϵup .* (1 .- nup_mean[:, N-1])) / N
+        Zndn_ratio[N] = sum(expϵdn .* (1 .- ndn_mean[:, N-1])) / N
+        nup_mean[:, N] = expϵup .* (1 .- nup_mean[:, N - 1]) / Znup_ratio[N]
+        ndn_mean[:, N] = expϵdn .* (1 .- ndn_mean[:, N - 1]) / Zndn_ratio[N]
+    end
+    return nup_mean, ndn_mean, Znup_ratio, Zndn_ratio
+end
+
+function measure_energy(
+    self::Lattice, Uup_array::Array{Float64,3}, Udn_array::Array{Float64,3},
+    nup_mean::Array{Float64,2}, ndn_mean::Array{Float64,2}
+)
+    Eonebody, Etwobody = 0, 0
+    for i = 1 : self.nsites
+        for j = 1 : self.nsites
+            if abs(i - j) == 1 || abs(i - j) == self.nsites - 1
+                Dij_up = sum(Uup_array[:, i, j] .* nup_mean[:, self.nup])
+                Dij_dn = sum(Udn_array[:, i, j] .* ndn_mean[:, self.ndn])
+                Eonebody += -self.t * ( Dij_up + Dij_dn )
+            end
+        end
+    end
+    for i = 1 : self.nsites
+        Umatrix = Uup_array[:, i, i] * Udn_array[:, i, i]'
+        nmatrix = nup_mean[:, self.nup] * ndn_mean[:, self.ndn]'
+        Diiii = sum(Umatrix .* nmatrix)
+        Etwobody += self.U * Diiii
+    end
+    return Eonebody, Etwobody
+end
+
+# rank-one perturbation calculator
+function rankone_perturb!(
+    γ::Float64, σil::Int64, Ns::Int64,
+    MatProdUp_new::Array{Float64,2}, MatProdDn_new::Array{Float64,2}
+)
+    Δup, Δdn = Matrix{Float64}(I, size(MatProdUp_new)), Matrix{Float64}(I, size(MatProdDn_new))
+    Δup[Ns, Ns] = exp(-4 * γ * σil)
+    Δdn[Ns, Ns] = exp(4 * γ * σil)
+    MatProdUp_new = Δup * MatProdUp_new
+    MatProdDn_new = Δdn * MatProdDn_new
+    return MatProdUp_new, MatProdDn_new
+end
+
+# determine the sign of the accept_ratio
+@inline function Rsign(R::Float64)
+    Rsgn = 2 * (R > 0) - 1
+end
+
+function MCIntegration(lattice::Lattice, qmc::QMC, opt::Opt)
+    # initialize vectors that store MC data
+    Ek_array = Vector{Float64}() # kinetic energy
+    Ep_array = Vector{Float64}() # potential energy
+    E_array = Vector{Float64}() # total energy
+    sgn_array = Vector{Int64}() # sign
+    NiMax_array = Vector{Float64}()
+    # initialize a random configuration
+    σlist = 2 * (rand(Float64, (lattice.nsites, qmc.L)) .< 0.5) .- 1
+    MatProdUp, MatProdDn, EigensUp, EigensDn = matrix_product_QR(σlist, lattice, qmc, opt.stab_interval)
+    nup_mean, ndn_mean, Znup_ratio, Zndn_ratio = recursion(lattice, real(EigensUp.values), real(EigensDn.values))
+    sgn = Rsign(prod(Znup_ratio[1 : lattice.nup]) * prod(Zndn_ratio[1 : lattice.nup]))
+    for n = 1 : qmc.samples
+        σlist_new = copy(σlist)
+        # random flip a spin
+        #flipL, flipNs = rand(1 : qmc.L), rand(1 : lattice.nsites)
+        flip = rand(1 : qmc.L * lattice.nsites)
+        σlist_new[flip] *= -1
+        #rankone_perturb!(qmc.γ, σlist[flipNs, flipL], flipNs, MatProdUp_new, MatProdDn_new)
+        MatProdUp_new, MatProdDn_new, EigensUp_new, EigensDn_new = matrix_product_QR(σlist_new, lattice, qmc, opt.stab_interval)
+        nup_mean_new, ndn_mean_new, Znup_ratio_new, Zndn_ratio_new = recursion(lattice, real(EigensUp_new.values), real(EigensDn_new.values))
+        Rup = prod(Znup_ratio_new[1 : lattice.nup] ./ Znup_ratio[1 : lattice.nup])
+        Rdn = prod(Zndn_ratio_new[1 : lattice.ndn] ./ Zndn_ratio[1 : lattice.ndn])
+        sgn_new = Rsign(Rup * Rdn) * sgn
+        R = abs(Rup * Rdn)
+        accept_ratio = R / (1 + R)
+        if rand() <= accept_ratio   # accept
+            σlist = σlist_new
+            Znup_ratio, Zndn_ratio = Znup_ratio_new, Zndn_ratio_new
+            EigensUp, EigensDn = EigensUp_new, EigensDn_new
+            nup_mean, ndn_mean = nup_mean_new, ndn_mean_new
+            sgn = sgn_new
+            Uup_array, Udn_array = overlap_coefficient(lattice.nsites, EigensUp.vectors, EigensDn.vectors)
+            Eonebody, Etwobody = measure_energy(lattice, Uup_array, Udn_array, nup_mean, ndn_mean)
+            push!(sgn_array, sgn)
+            push!(Ek_array, Eonebody * sgn)
+            push!(Ep_array, Etwobody * sgn)
+            push!(E_array, Eonebody * sgn + Etwobody * sgn)
+            push!(NiMax_array, maximum(nup_mean[:, lattice.nup]))
+        else    # reject
+            Uup_array, Udn_array = overlap_coefficient(lattice.nsites, EigensUp.vectors, EigensDn.vectors)
+            Eonebody, Etwobody = measure_energy(lattice, Uup_array, Udn_array, nup_mean, ndn_mean)
+            push!(sgn_array, sgn)
+            push!(Ek_array, Eonebody * sgn)
+            push!(Ep_array, Etwobody * sgn)
+            push!(E_array, Eonebody * sgn + Etwobody * sgn)
+            push!(NiMax_array, maximum(nup_mean[:, lattice.nup]))
+        end
+    end
+    E_mean = mean(E_array[opt.nburn_in : qmc.samples])
+    E_error = std(E_array[opt.nburn_in : qmc.samples]) / sqrt(qmc.samples - opt.nburn_in)
+    Ek_mean = mean(Ek_array[opt.nburn_in : qmc.samples])
+    Ek_error = std(Ek_array[opt.nburn_in : qmc.samples]) / sqrt(qmc.samples - opt.nburn_in)
+    Ep_mean = mean(Ep_array[opt.nburn_in : qmc.samples])
+    Ep_error = std(Ep_array[opt.nburn_in : qmc.samples]) / sqrt(qmc.samples - opt.nburn_in)
+    sgn_avg = mean(sgn_array[opt.nburn_in : qmc.samples])
+    NiMax_avg = mean(NiMax_array[opt.nburn_in : qmc.samples])
+    return E_mean, E_error, Ek_mean, Ek_error, Ep_mean, Ep_error, sgn_avg, NiMax_avg
+end

FermiHubbard/InputFermiHubbard

@@ -0,0 +1,13 @@
+Beta    4
+NumOfSites	18
+NumOfSpinUp	9
+NumOfSpinDn 9
+delta_tau	0.05
+t	1
+U	4
+NumBlocks   1
+NumOfSamples  1e4
+
+## Optimization Parameters ##
+NumOfStablizationInterval   80
+NumofBurn-ins   1

FermiHubbard/Main.jl

@@ -0,0 +1,116 @@
+using Statistics
+using DelimitedFiles
+using Random
+using Printf
+
+struct Lattice
+    nsites::Int64
+    nup::Int64
+    ndn::Int64
+    ntot::Int64
+    t::Float64
+    U::Float64
+end
+
+struct QMC
+    nblocks::Int64
+    samples::Int64
+    Δτ::Float64
+    L::Int64
+    γ::Float64
+    Kmatrix::Array{Float64,2}
+end
+
+struct Opt  # Optimization
+    stab_interval::Int64
+    nburn_in::Int64
+end
+
+const FloatType = Union{Float64, Complex{Float64}}
+
+include("./MatrixGenerator.jl")
+include("./Calculations.jl")
+
+function getinputs()
+    lines = readlines("./InputFermiHubbard")
+
+    # basic parameters
+    β = parse(Float64, split(lines[1])[2])
+    num_sites = parse(Int64, split(lines[2])[2])
+    num_spinup = parse(Int64, split(lines[3])[2])
+    num_spindn = parse(Int64, split(lines[4])[2])
+    num_elec = num_spinup + num_spinup
+
+    Δτ = parse(Float64, split(lines[5])[2])
+    t = parse(Float64, split(lines[6])[2])
+    U = parse(Float64, split(lines[7])[2])
+    num_blocks = convert(Int64, parse(Float64, split(lines[8])[2]))
+    num_samples = convert(Int64, parse(Float64, split(lines[9])[2]))
+
+    # optimization parameters
+    stab_interval = parse(Int64, split(lines[12])[2])
+    num_burnins = convert(Int64, parse(Float64, split(lines[13])[2]))
+
+    # constants for calculations
+    L = convert(Int64, round(β / Δτ))
+    γ = atanh(sqrt(tanh(Δτ * U / 4)))
+
+    lattice = Lattice(num_sites, num_spinup, num_spindn, num_elec, t, U)
+
+    K_matrix = kinetic_matrix(lattice)
+    qmc = QMC(num_blocks, num_samples, Δτ, L, γ, exp(-K_matrix * Δτ / 2))
+
+    opt = Opt(stab_interval, num_burnins)
+    return lattice, qmc, opt
+end
+
+function mainbody()
+    lattice, qmc, opt = getinputs()
+    SubEk_array = Vector{Float64}()
+    SubEp_array = Vector{Float64}()
+    SubE_array = Vector{Float64}()
+    SubEk_error_array = Vector{Float64}()
+    SubEp_error_array = Vector{Float64}()
+    SubE_error_array = Vector{Float64}()
+    filename = string("output_Ns", lattice.nsites, "_U", lattice.U, "_beta", qmc.Δτ * qmc.L)
+    outputfile = open(filename, "a+")
+    @printf(outputfile, "       Canonical Ensemble AFQMC for Hubbard Model     \n\n")
+    @printf(outputfile, "Number of Blocks: %d\n", qmc.nblocks)
+    @printf(outputfile, "Length of Img Timestep: %.4f\n", qmc.Δτ)
+    @printf(outputfile, "Stablization Interval: %d\n\n", opt.stab_interval)
+    for i = 1 : qmc.nblocks
+        SubE, SubE_error, SubEk, SubEk_error, SubEp, SubEp_error, sgn_mean, NiMax = MCIntegration(lattice, qmc, opt)
+        SubE = SubE / sgn_mean
+        SubEk = SubEk / sgn_mean
+        SubEp = SubEp / sgn_mean
+
+        @printf(outputfile, "Block # %d out of %d\n", i, qmc.nblocks)
+        @printf(outputfile, "--------------------\n")
+        @printf(outputfile, "Average Sign: %.2f\n", sgn_mean)
+        @printf(outputfile, "Maximum Occupancy: %.2f\n", NiMax)
+        @printf(outputfile, "Ek            Ep            Etot\n")
+        @printf(outputfile, "%.5f     %.5f       %.5f\n", SubEk, SubEp, SubE)
+        @printf(outputfile, "Error:\n")
+        @printf(outputfile, "%.5f     %.5f       %.5f\n\n", SubEk_error, SubEp_error, SubE_error)
+
+        push!(SubEk_array, SubEk)
+        push!(SubEk_error_array, SubEk_error)
+        push!(SubEp_array, SubEp)
+        push!(SubEp_error_array, SubEp_error)
+        push!(SubE_array, SubE)
+        push!(SubE_error_array, SubE_error)
+    end
+    Ek_mean, Ek_error = mean(SubEk_array), mean(SubEk_error_array)
+    Ep_mean, Ep_error = mean(SubEp_array), mean(SubEp_error_array)
+    E_mean, E_error = mean(SubE_array), mean(SubE_error_array)
+
+    @printf(outputfile, "--------------------\n")
+    @printf(outputfile, "Results for this run:\n")
+    @printf(outputfile, "Ek            Ep            Etot\n")
+    @printf(outputfile, "%.5f     %.5f       %.5f\n", Ek_mean, Ep_mean, E_mean)
+    @printf(outputfile, "Error:\n")
+    @printf(outputfile, "%.5f     %.5f       %.5f\n\n", Ek_error, Ep_error, E_error)
+    close(outputfile)
+end
+
+mainbody()

FermiHubbard/MatrixGenerator.jl

@@ -0,0 +1,89 @@
+using LinearAlgebra
+
+# one-body matrix, nearest neighbour, PBC
+function kinetic_matrix(self::Lattice)
+    K_matrix = zeros(self.nsites, self.nsites)
+    for i =  1 : self.nsites
+        for j = 1 : self.nsites
+            if abs(i - j) == 1 || abs(i - j) == self.nsites - 1
+                K_matrix[i, j] = -self.t
+            end
+        end
+    end
+    return K_matrix
+end
+
+# σ-dependent generator
+function auxiliary_field_matrix_stepwise(σ::Array{Int64,1}, γ::Float64, Δτ::Float64, U::Float64)
+    AF_matrix_SpinUp = diagm(0 => exp.(2 * γ * σ .- Δτ * U / 2))
+    AF_matrix_SpinDn = diagm(0 => exp.(-2 * γ * σ .- Δτ * U / 2))
+    return AF_matrix_SpinUp, AF_matrix_SpinDn
+end
+
+function matrix_product_stepwise(σlist::Array{Int64,2}, lattice::Lattice, qmc::QMC)
+    MatProdUp = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    MatProdDn = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    for i = 1 : qmc.L
+        σ = σlist[:, i]
+        AF_matrix_SpinUp, AF_matrix_SpinDn = auxiliary_field_matrix_stepwise(σ, qmc.γ, qmc.Δτ, lattice.U)
+        MatProdUp = Symmetric(qmc.Kmatrix * AF_matrix_SpinUp * qmc.Kmatrix) * MatProdUp
+        MatProdDn = Symmetric(qmc.Kmatrix * AF_matrix_SpinDn * qmc.Kmatrix) * MatProdDn
+    end
+    EigensUp = eigen(MatProdUp)
+    EigensDn = eigen(MatProdDn)
+    return MatProdUp, MatProdDn, EigensUp, EigensDn
+end
+
+# Use QR decomposition to stablize the propogation
+function matrix_product_QR(σlist::Array{Int64,2}, lattice::Lattice, qmc::QMC, n_stab_interval::Int64)
+    # Initialize QR matrices
+    QRUp_Q = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    QRUp_R = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    QRDn_Q = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    QRDn_R = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    # Final Product Matrix
+    MatProdUp = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    MatProdDn = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+    for i = 1 : div(qmc.L, n_stab_interval)
+        matprodUp = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+        matprodDn = Matrix{Float64}(I, lattice.nsites, lattice.nsites)
+        for j = 1 : n_stab_interval
+            σ = σlist[:, (i - 1) * n_stab_interval + j]
+            AF_matrix_SpinUp, AF_matrix_SpinDn = auxiliary_field_matrix_stepwise(σ, qmc.γ, qmc.Δτ, lattice.U)
+            matprodUp = Symmetric(qmc.Kmatrix * AF_matrix_SpinUp * qmc.Kmatrix) * matprodUp
+            matprodDn = Symmetric(qmc.Kmatrix * AF_matrix_SpinDn * qmc.Kmatrix) * matprodDn
+        end
+        # QR decomposition of the matrix product
+        # leave here to see if can be replaced with rmul!
+        QRUp_Q = matprodUp * QRUp_Q
+        DecompUp = qr(QRUp_Q)
+        QRDn_Q = matprodDn * QRDn_Q
+        DecompDn = qr(QRDn_Q)
+        # Remaining R_i goes to multiply R_1
+        QRUp_Q = DecompUp.Q
+        QRUp_R = DecompUp.R * QRUp_R
+        QRDn_Q = DecompDn.Q
+        QRDn_R = DecompDn.R * QRDn_R
+    end
+    mul!(MatProdUp, QRUp_Q, QRUp_R)
+    mul!(MatProdDn, QRDn_Q, QRDn_R)
+    EigensUp = eigen(MatProdUp)
+    EigensDn = eigen(MatProdDn)
+    return MatProdUp, MatProdDn, EigensUp, EigensDn
+end
+
+# generate overlap "tensor" / column major
+function overlap_coefficient(nsites::Int64, Uup::Array{T1,2}, Udn::Array{T2,2}) where {T1<:FloatType, T2<:FloatType}
+    Uup_inv, Udn_inv = inv(Uup), inv(Udn)
+    Uup_array = complex(zeros(nsites, nsites, nsites))
+    Udn_array = complex(zeros(nsites, nsites, nsites))
+    for i = 1 : nsites
+        for j = 1 : nsites
+            for k = 1 : nsites
+                Uup_array[k ,i, j] = Uup_inv[k, i] * Uup[j, k]
+                Udn_array[k, i, j] = Udn_inv[k, i] * Udn[j, k]
+            end
+        end
+    end
+    real(Uup_array), real(Udn_array)
+end