Add files via upload

BoseHubbard/Calculations.jl

@@ -0,0 +1,134 @@
+using LinearAlgebra
+
+function time_slice_spectrum()
+    H_matrix = kinetic_matrix_imgtime * auxiliary_field_matrix() * kinetic_matrix_imgtime
+end
+
+function sub_partition_weight(L)
+    matrix_product = Matrix{Complex{Float64}}(I, num_sites, num_sites)
+    for i = 1 : L
+        matrix_product *= time_slice_spectrum()
+    end
+    #return matrix_product
+    #if flag_background_subtraction == 1
+        #matrix_product *= exp(sum(β * (U/2) * background_dist.^2))
+        #return matrix_product
+    #end
+    effective_spectrum, u = eigen(matrix_product) # spectrum element: e^(-βϵᵢ)
+    u_inv = inv(u)
+    U_array = complex(zeros(num_sites, num_sites, num_sites))
+    for i = 1 : num_sites
+        for j = 1 : num_sites
+            for k = 1 : num_sites
+                U_array[i, j, k] = u_inv[k, i] * u[j, k]
+            end
+        end
+    end
+    return effective_spectrum, U_array
+end
+
+function recursion_partition_function(effective_spectrum)
+    Z = complex(zeros(num_bosons))
+    zₖ = complex(zeros(num_bosons))
+    for k = 1 : num_bosons
+        zₖ[k] = sum(effective_spectrum .^ k)
+    end
+    for N = 1 : num_bosons
+        for M = 1 : N
+            if M == N
+                Z[N] += 1 * zₖ[M]
+            else
+                Z[N] += Z[N-M] * zₖ[M]
+            end
+        end
+        Z[N] = Z[N] / N
+    end
+    return Z
+end
+
+function recursion_OccNum(Z, effective_spectrum) # Calculating <nᵢ>
+    n_i = complex(zeros(num_bosons, num_sites))
+    n_i[1, :] = effective_spectrum
+    for N = 2 : num_bosons
+        n_i[N, :] = effective_spectrum .* (Z[N-1] .+ n_i[N-1, :])
+    end
+    return n_i
+end
+
+function recursion_double_OccNum_eq(Z, effective_spectrum, n_i) # Calculating <nᵢ^2>
+    n²_i = complex(zeros(num_bosons, num_sites)) # <nᵢ^2> array
+    n²_i[1, :] = effective_spectrum
+    for N = 2 : num_bosons
+        n²_i[N, :] = effective_spectrum .* (Z[N-1] .+ 2 * n_i[N-1, :] .+ n²_i[N-1, :])
+    end
+    return n²_i[num_bosons, :]
+end
+
+function recursion_double_OccNum_neq(Z, effective_spectrum, n_i) # Calculating <nᵢnⱼ>
+    n²_i = complex(zeros(num_bosons, num_sites, num_sites)) # <nᵢnⱼ> array
+    for i = 1 : num_sites
+        for j = i + 1 : num_sites
+            for N = 2 : num_bosons
+                n²_i[N, i, j] = effective_spectrum[j] / 2 * (n_i[N-1, i] + n²_i[N-1, i, j]) +
+                            + effective_spectrum[i] / 2 * (n_i[N-1, j] + n²_i[N-1, i, j])
+            end
+            n²_i[num_bosons, j, i] = n²_i[num_bosons, i, j]
+        end
+    end
+    return n²_i[num_bosons, :, :]
+end
+
+#function measure_condensate_fraciton()
+#end
+
+function measure_TotEnergy(num_iterations)
+    partition_func_array = Vector{Complex{Float64}}()
+    Ek_array = Vector{Complex{Float64}}()
+    Ep_array = Vector{Complex{Float64}}()
+    E_tot_array = Vector{Complex{Float64}}()
+    for N = 1 : num_iterations
+        effective_spectrum, U_array = sub_partition_weight(L)
+        Z = recursion_partition_function(effective_spectrum)
+        n_i = recursion_OccNum(Z, effective_spectrum)
+        n²_i_diag = recursion_double_OccNum_eq(Z, effective_spectrum, n_i)
+        n²_i_offdiag = recursion_double_OccNum_neq(Z, effective_spectrum, n_i)
+        E_onebody, E_twobody = 0, 0
+        # one-body energy calculation
+        for i = 1 : num_sites
+            for j = 1 : num_sites
+                if abs(i - j) == 1 || abs(i - j) == num_sites - 1
+                    D_ij = sum(U_array[i, j, :] .* n_i[num_bosons, :])
+                    E_onebody += -t * D_ij
+                end
+            end
+        end
+        # two-body energy calculation
+        for i = 1 : num_sites
+            U_matrix = complex(zeros(num_sites, num_sites))
+            D_iiii = 0
+            for p = 1 :num_sites
+                for r = 1 : num_sites
+                    if p != r
+                        U_matrix[p, r] = U_array[i, i, p] * U_array[i, i, r]
+                        D_iiii += U_matrix[p, r] * n_i[num_bosons, p]
+                    end
+                end
+                D_iiii += U_array[i, i, p]^2 * n²_i_diag[p]
+            end
+            D_iiii += sum( 2 * U_matrix .* n²_i_offdiag)
+            E_twobody += U/2 * D_iiii
+        end
+        push!(partition_func_array, Z[num_bosons])
+        push!(Ek_array, E_onebody)
+        push!(Ep_array, E_twobody)
+        push!(E_tot_array, E_onebody + E_twobody)
+    end
+    Z_mean = mean(partition_func_array)
+    Ek_mean = mean(Ek_array) / Z_mean
+    Ek_error = std(Ek_array) / (sqrt(num_iterations) * real(Z_mean))
+    Ep_mean = mean(Ep_array) / Z_mean
+    Ep_error = std(Ep_array) / (sqrt(num_iterations) * real(Z_mean))
+    E_tot_mean = mean(E_tot_array) / Z_mean
+    E_tot_error_bar = std(E_tot_array) / (sqrt(num_iterations) * real(Z_mean))
+    return Ek_mean, Ek_error, Ep_mean, Ep_error, E_tot_mean, E_tot_error_bar
+end

BoseHubbard/InputBoseHubbard

@@ -0,0 +1,14 @@
+beta    0.2
+number_sites	3
+numeber_bosons	3
+delta_tau	0.005
+t	1
+U	2
+num_iteration   1e3
+flag_background_subtraction 1
+flag_density_matrix 1
+# flag for obsevables #
+KineticEnergy 1
+PotentialEnergy 0
+TotalEnergy 0
+CondensationFraction 0

BoseHubbard/Main.jl

@@ -0,0 +1,67 @@
+using Statistics
+using DelimitedFiles
+
+include("MatrixGenerator.jl")
+include("Calculations.jl")
+include("MeanFieldSolver.jl")
+include("Test.jl")
+
+# get inputs, initialization
+lines = readlines("InputBoseHubbard")
+# inverse temperature
+temp_name, β_str = split(lines[1])
+const β = parse(Float64, β_str)
+# number of sites
+temp_name, num_sites_str = split(lines[2])
+const num_sites = parse(Int64, num_sites_str)
+# number of bosons
+temp_name, num_bosons_str = split(lines[3])
+const num_bosons = parse(Int64, num_bosons_str)
+# length of the imaginary time interval
+temp_name, Δτ_str = split(lines[4])
+const Δτ = parse(Float64, Δτ_str)
+# hopping constant
+temp_name, t_str = split(lines[5])
+const t = parse(Float64, t_str)
+# on-site repulsion constant
+temp_name, U_str = split(lines[6])
+const U = parse(Float64, U_str)
+# number of iterations
+temp_name, num_iterations_str = split(lines[7])
+const num_iterations = convert(Int64, parse(Float64, num_iterations_str))
+# flag for the background subtraction
+temp_name, flag_background_subtraction_str = split(lines[8])
+const flag_background_subtraction = parse(Int16, flag_background_subtraction_str)
+# flag for calculating the density matrix
+temp_name, flag_density_matrix_str = split(lines[9])
+const flag_density_matrix = convert(Int64, parse(Int16, flag_density_matrix_str))
+
+const background_dist = ones(num_sites) / num_sites #  Equally distributed among sites
+const residue_term = sum(β * (U/2) * background_dist.^2)
+
+# Constants for calculations
+# Kinetic matrix in each imaginary time slice is fixed throughout the calculation.
+const kinetic_matrix_imgtime = complex(exp(-kinetic_matrix(num_sites, t) * Δτ / 2))
+# number of imaginary time steps
+const L = convert(Int64, round(β / Δτ))
+# coefficient of the H-S transformed potential term
+const AF_coefficient = sqrt(complex(-Δτ * U))
+
+#mean_value, error_bar = Z_iteration(num_iteration, L)
+
+#write_to_file(mean_value, error_bar, Δτ)
+
+#function write_to_file(mean_value, error_bar, Δτ)
+    #output_file = open("result","a+")
+    #write(output_file, string(β), " ", string(Δτ), " ", string(mean_value), " ", string(error_bar), "\n")
+    #close(output_file)
+#end
+
+Ek_mean, Ek_error, Ep_mean, Ep_error, E_tot_mean, E_tot_error_bar = measure_TotEnergy(num_iterations)
+println("U = ", U)
+println("Beta = ", β)
+println("Ns = ", num_sites)
+println("NBosons = ", num_bosons)
+println("Ek = ", Ek_mean, ", Error = ", Ek_error)
+println("Ep = ", Ep_mean, ", Error = ", Ep_error)
+println("Etot = ", E_tot_mean, ", Error = ", E_tot_error_bar)

BoseHubbard/MatrixGenerator.jl

@@ -0,0 +1,33 @@
+using LinearAlgebra
+using Distributions
+
+function kinetic_matrix(num_sites, t)   # one-body matrix, periodic boundary condition
+    K_matrix = zeros(num_sites, num_sites)
+    for i =  1 : num_sites
+        for j = 1 : num_sites
+            if abs(i - j) == 1 || abs(i - j) == num_sites - 1
+                K_matrix[i, j] = -t
+            end
+        end
+    end
+    if flag_background_subtraction == 1
+        K_matrix += Diagonal(background_dist) * U
+    end
+    return K_matrix
+end
+
+function auxiliary_field_matrix()
+    ϕ = rand(Normal(), num_sites)
+    coefficent_vector = ϕ * AF_coefficient
+    AF_matrix = complex(zeros(num_sites, num_sites))
+    for i = 1 : num_sites
+        AF_matrix[i, i] = coefficent_vector[i]
+    end
+    if flag_background_subtraction == 1
+        background_exponent = sum(-coefficent_vector .* background_dist)
+        background_weight = exp(background_exponent)
+        return background_weight * exp(AF_matrix)
+    elseif flag_background_subtraction == 0
+        return exp(AF_matrix)
+    end
+end

BoseHubbard/MeanFieldSolver.jl

@@ -0,0 +1,24 @@
+using LinearAlgebra
+
+include("MatrixGenerator.jl")
+
+function mean_field_solver(β, t, U, num_sites)
+    convergence_tolerance = 1e-6
+    kinetic_val = kinetic_matrix(num_sites, t)
+    n_set = ones(num_sites) / num_sites
+    # each element in the array represents the density of particles on each site
+    trial_matrix = kinetic_val + Diagonal(n_set * U)
+    n_set_old = zeros(num_sites)
+    while sum(abs.(n_set_old - n_set)) > convergence_tolerance
+        energy_spectrum, basis_vector = eigen(trial_matrix)
+        density_matrix = zeros(num_sites, num_sites)
+        n_set_old = n_set
+        Z = sum(exp.(- β * energy_spectrum))
+        for i = 1:num_sites
+            density_matrix += (exp.(- β * energy_spectrum[i])  / Z) * (basis_vector[:, i] * basis_vector[:, i]')
+        end
+        n_set = diag(density_matrix)
+        trial_matrix = kinetic_val + Diagonal(n_set * U)
+    end
+    return n_set
+end

BoseHubbard/Test.jl

@@ -0,0 +1,47 @@
+function test_density_matrix()
+    effective_spectrum, U_array = sub_partition_weight(L)
+    Z = recursion_partition_function(effective_spectrum)
+    D = complex(zeros(num_sites, num_sites))
+    for i = 1 : num_sites
+        for j = 1 : num_sites
+            n_i = recursion_OccNum(Z, effective_spectrum)
+            D[i, j] = sum(U_array[i, j, :] .* n_i[num_bosons, :])
+        end
+    end
+    return Z[num_bosons], D
+end
+
+function test_density_matrix_sampling(num_iterations) # tests for one-body density matrices
+    partition_func_array = Vector{Complex{Float64}}()
+    density_mat_array = Matrix{Complex{Float64}}[]
+    for N = 1 : num_iterations
+        partition_func_sample, density_mat_sample = test_density_matrix()
+        push!(partition_func_array, partition_func_sample)
+        push!(density_mat_array, density_mat_sample)
+    end
+    Z_mean = mean(partition_func_array)
+    Z_error_bar = stdm(partition_func_array, Z_mean) / sqrt(num_iterations)
+    D_mean = mean(density_mat_array) / Z_mean
+    return Z_mean, Z_error_bar, D_mean
+end
+
+function test_matrix_product(num_iterations)
+    mat_product_array = Vector{Array{Complex{Float64},1}}()
+    for N = 1 : num_iterations
+        mat_product_sample = sub_partition_weight(L)
+        push!(mat_product_array, mat_product_sample)
+    end
+    mat_product_mean = mean(mat_product_array)
+end
+
+function test_partition_function(num_iterations)
+    partition_func_array = Vector{Complex{Float64}}()
+    for i = 1 : num_iterations
+        effective_spectrum, U_array = sub_partition_weight(L)
+        Z = recursion_partition_function(effective_spectrum)
+        push!(partition_func_array, Z[num_bosons])
+    end
+    Z_mean = mean(partition_func_array)
+    Z_error_bar = stdm(partition_func_array, Z_mean) / num_iterations
+    return Z_mean, Z_error_bar
+end