Download this example as a Julia script.
Honeycomb Hubbard model
In this script we simulate the Hubbard model on a Honeycomb lattice, with a Hamiltonian given by
\[\begin{align*} \hat{H} = & -t \sum_{\sigma,\langle i, j \rangle} (\hat{c}^{\dagger}_{\sigma,i}, \hat{c}^{\phantom \dagger}_{\sigma,j} + {\rm h.c.}) -\mu \sum_{\sigma,i}\hat{n}_{\sigma,i}\\ & + U \sum_{i} (\hat{n}_{\uparrow,i}-\tfrac{1}{2})(\hat{n}_{\downarrow,i}-\tfrac{1}{2}), \end{align*}\]
where $\hat{c}^\dagger_{\sigma,i} \ (\hat{c}^{\phantom \dagger}_{\sigma,i})$ creates (annihilates) a spin $\sigma$ electron on site $i$ in the lattice, and $\hat{n}_{\sigma,i} = \hat{c}^\dagger_{\sigma,i} \hat{c}^{\phantom \dagger}_{\sigma,i}$ is the spin-$\sigma$ electron number operator for site $i$. The nearest-neighbor hopping amplitude is $t$ and $\mu$ is the chemical potential. The strength of the repulsive Hubbard interaction is controlled by $U>0$.
A short test simulation using the script associated with this example can be run as
> julia hubbard_honeycomb.jl 0 6.0 0.0 4.0 3 2000 10000 50This simulates a half-filled $(\mu = 0.0)$ Hubabrd model on a $3 \times 3$ unit cell honeycomb lattice, with a Hubbard interaction of $U = 6.0$, at inverse temperature $\beta = 4.0$.
Below you will find the source code from the julia script linked at the top of this page, but with additional comments giving more detailed explanations for what certain parts of the code are doing.
using LinearAlgebra
using Random
using Printf
using SmoQyDQMC
import SmoQyDQMC.LatticeUtilities as lu
import SmoQyDQMC.JDQMCFramework as dqmcf
import SmoQyDQMC.JDQMCMeasurements as dqmcm
# Define top-level function for running DQMC simulation.
function run_hubbard_chain_simulation(sID, U, μ, β, L, N_burnin, N_updates, N_bins; filepath = ".")
# Construct the foldername the data will be written to.
datafolder_prefix = @sprintf "hubbard_honeycomb_U%.2f_mu%.2f_L%d_b%.2f" U μ L β
# Initialize an instance of the SimulationInfo type.
simulation_info = SimulationInfo(
filepath = filepath,
datafolder_prefix = datafolder_prefix,
sID = sID
)
# Initialize the directory the data will be written to.
initialize_datafolder(simulation_info)
# Initialize a random number generator that will be used throughout the simulation.
seed = abs(rand(Int))
rng = Xoshiro(seed)
# Set the discretization in imaginary time for the DQMC simulation.
Δτ = 0.10
# Calculate the length of the imaginary time axis, Lτ = β/Δτ.
Lτ = dqmcf.eval_length_imaginary_axis(β, Δτ)
# This flag indicates whether or not to use the checkboard approximation to
# represent the exponentiated hopping matrix exp(-Δτ⋅K)
checkerboard = false
# Whether the propagator matrices should be represented using the
# symmetric form B = exp(-Δτ⋅K/2)⋅exp(-Δτ⋅V)⋅exp(-Δτ⋅K/2)
# or the asymetric form B = exp(-Δτ⋅V)⋅exp(-Δτ⋅K)
symmetric = false
# Set the initial period in imaginary time slices with which the Green's function matrices
# will be recomputed using a numerically stable procedure.
n_stab = 10
# Specify the maximum allowed error in any element of the Green's function matrix that is
# corrected by performing numerical stabiliziation.
δG_max = 1e-6
# Calculate the bins size.
bin_size = div(N_updates, N_bins)
# Initialize a dictionary to store additional information about the simulation.
additional_info = Dict(
"dG_max" => δG_max,
"N_burnin" => N_burnin,
"N_updates" => N_updates,
"N_bins" => N_bins,
"bin_size" => bin_size,
"local_acceptance_rate" => 0.0,
"n_stab_init" => n_stab,
"symmetric" => symmetric,
"checkerboard" => checkerboard,
"seed" => seed,
)
#######################
### DEFINE THE MODEL ##
#######################
# Initialize an instance of the type UnitCell.
unit_cell = lu.UnitCell(
lattice_vecs = [[3/2,√3/2],
[3/2,-√3/2]],
basis_vecs = [[0.,0.],
[1.,0.]]
)
# Initialize an instance of the type Lattice.
lattice = lu.Lattice(
L = [L, L],
periodic = [true, true]
)
# Initialize an instance of the ModelGeometry type.
model_geometry = ModelGeometry(unit_cell, lattice)
# Define the first nearest-neighbor bond in a honeycomb lattice.
bond_1 = lu.Bond(orbitals = (1,2), displacement = [0,0])
# Add the first nearest-neighbor bond in a honeycomb lattice to the model.
bond_1_id = add_bond!(model_geometry, bond_1)
# Define the second nearest-neighbor bond in a honeycomb lattice.
bond_2 = lu.Bond(orbitals = (1,2), displacement = [-1,0])
# Add the second nearest-neighbor bond in a honeycomb lattice to the model.
bond_2_id = add_bond!(model_geometry, bond_2)
# Define the third nearest-neighbor bond in a honeycomb lattice.
bond_3 = lu.Bond(orbitals = (1,2), displacement = [0,-1])
# Add the third nearest-neighbor bond in a honeycomb lattice to the model.
bond_3_id = add_bond!(model_geometry, bond_3)
# Define nearest-neighbor hopping amplitude, setting the energy scale for the system.
t = 1.0
# Define the tight-binding model
tight_binding_model = TightBindingModel(
model_geometry = model_geometry,
t_bonds = [bond_1, bond_2, bond_3], # defines hopping
t_mean = [t, t, t], # defines corresponding hopping amplitude
μ = μ, # set chemical potential
ϵ_mean = [0.0, 0.0] # set the (mean) on-site energy
)
# Initialize the Hubbard interaction in the model.
hubbard_model = HubbardModel(
shifted = false, # If true, Hubbard interaction instead parameterized as U⋅nup⋅ndn
U_orbital = [1, 2],
U_mean = [U, U],
)
# Write the model summary to file.
model_summary(
simulation_info = simulation_info,
β = β, Δτ = Δτ,
model_geometry = model_geometry,
tight_binding_model = tight_binding_model,
interactions = (hubbard_model,)
)
#########################################
### INITIALIZE FINITE MODEL PARAMETERS ##
#########################################
# Initialize tight-binding parameters.
tight_binding_parameters = TightBindingParameters(
tight_binding_model = tight_binding_model,
model_geometry = model_geometry,
rng = rng
)
# Initialize Hubbard interaction parameters.
hubbard_parameters = HubbardParameters(
model_geometry = model_geometry,
hubbard_model = hubbard_model,
rng = rng
)
# Apply Ising Hubbard-Stranonvich (HS) transformation, and initialize
# corresponding HS fields that will be sampled in DQMC simulation.
hubbard_ising_parameters = HubbardIsingHSParameters(
β = β, Δτ = Δτ,
hubbard_parameters = hubbard_parameters,
rng = rng
)
##############################
### INITIALIZE MEASUREMENTS ##
##############################
# Initialize the container that measurements will be accumulated into.
measurement_container = initialize_measurement_container(model_geometry, β, Δτ)
# Initialize the tight-binding model related measurements, like the hopping energy.
initialize_measurements!(measurement_container, tight_binding_model)
# Initialize the Hubbard interaction related measurements.
initialize_measurements!(measurement_container, hubbard_model)
# Initialize the single-particle electron Green's function measurement.
initialize_correlation_measurements!(
measurement_container = measurement_container,
model_geometry = model_geometry,
correlation = "greens",
time_displaced = true,
pairs = [(1, 1), (2, 2), (1, 2)]
)
# Initialize density correlation function measurement.
initialize_correlation_measurements!(
measurement_container = measurement_container,
model_geometry = model_geometry,
correlation = "density",
time_displaced = false,
integrated = true,
pairs = [(1, 1), (2, 2), (1, 2)]
)
# Initialize the pair correlation function measurement.
initialize_correlation_measurements!(
measurement_container = measurement_container,
model_geometry = model_geometry,
correlation = "pair",
time_displaced = false,
integrated = true,
pairs = [(1, 1), (2, 2), (1, 2)]
)
# Initialize the spin-z correlation function measurement.
initialize_correlation_measurements!(
measurement_container = measurement_container,
model_geometry = model_geometry,
correlation = "spin_z",
time_displaced = false,
integrated = true,
pairs = [(1, 1), (2, 2), (1, 2)]
)
# Initialize the sub-directories to which the various measurements will be written.
initialize_measurement_directories(
simulation_info = simulation_info,
measurement_container = measurement_container
)
#############################
### SET-UP DQMC SIMULATION ##
#############################
# Allocate FermionPathIntegral type for both the spin-up and spin-down electrons.
fermion_path_integral_up = FermionPathIntegral(tight_binding_parameters = tight_binding_parameters, β = β, Δτ = Δτ)
fermion_path_integral_dn = FermionPathIntegral(tight_binding_parameters = tight_binding_parameters, β = β, Δτ = Δτ)
# Initialize the FermionPathIntegral type for both the spin-up and spin-down electrons.
initialize!(fermion_path_integral_up, fermion_path_integral_dn, hubbard_parameters)
initialize!(fermion_path_integral_up, fermion_path_integral_dn, hubbard_ising_parameters)
# Initialize the imaginary-time propagators for each imaginary-time slice for both the
# spin-up and spin-down electrons.
Bup = initialize_propagators(fermion_path_integral_up, symmetric=symmetric, checkerboard=checkerboard)
Bdn = initialize_propagators(fermion_path_integral_dn, symmetric=symmetric, checkerboard=checkerboard)
# Initialize FermionGreensCalculator for the spin-up and spin-down electrons.
fermion_greens_calculator_up = dqmcf.FermionGreensCalculator(Bup, β, Δτ, n_stab)
fermion_greens_calculator_dn = dqmcf.FermionGreensCalculator(Bdn, β, Δτ, n_stab)
# Allcoate matrices for spin-up and spin-down electron Green's function matrices.
Gup = zeros(eltype(Bup[1]), size(Bup[1]))
Gdn = zeros(eltype(Bdn[1]), size(Bdn[1]))
# Initialize the spin-up and spin-down electron Green's function matrices, also
# calculating their respective determinants as the same time.
logdetGup, sgndetGup = dqmcf.calculate_equaltime_greens!(Gup, fermion_greens_calculator_up)
logdetGdn, sgndetGdn = dqmcf.calculate_equaltime_greens!(Gdn, fermion_greens_calculator_dn)
# Allocate matrices for various time-displaced Green's function matrices.
Gup_ττ = similar(Gup) # G↑(τ,τ)
Gup_τ0 = similar(Gup) # G↑(τ,0)
Gup_0τ = similar(Gup) # G↑(0,τ)
Gdn_ττ = similar(Gdn) # G↓(τ,τ)
Gdn_τ0 = similar(Gdn) # G↓(τ,0)
Gdn_0τ = similar(Gdn) # G↓(0,τ)
# Initialize variables to keep track of the largest numerical error in the
# Green's function matrices corrected by numerical stabalization.
δG = zero(typeof(logdetGup))
δθ = zero(typeof(sgndetGup))
####################################
### BURNIN/THERMALIZATION UPDATES ##
####################################
# Iterate over burnin/thermalization updates.
for n in 1:N_burnin
# Perform a sweep through the lattice, attemping an update to each Ising HS field.
(acceptance_rate, logdetGup, sgndetGup, logdetGdn, sgndetGdn, δG, δθ) = local_updates!(
Gup, logdetGup, sgndetGup, Gdn, logdetGdn, sgndetGdn,
hubbard_ising_parameters,
fermion_path_integral_up = fermion_path_integral_up,
fermion_path_integral_dn = fermion_path_integral_dn,
fermion_greens_calculator_up = fermion_greens_calculator_up,
fermion_greens_calculator_dn = fermion_greens_calculator_dn,
Bup = Bup, Bdn = Bdn, δG_max = δG_max, δG = δG, δθ = δθ, rng = rng
)
# Record the acceptance rate for the attempted local updates to the HS fields.
additional_info["local_acceptance_rate"] += acceptance_rate
end
################################
### START MAKING MEAUSREMENTS ##
################################
# Re-initialize variables to keep track of the largest numerical error in the
# Green's function matrices corrected by numerical stabalization.
δG = zero(typeof(logdetGup))
δθ = zero(typeof(sgndetGup))
# Iterate over the number of bin, i.e. the number of time measurements will be dumped to file.
for bin in 1:N_bins
# Iterate over the number of updates and measurements performed in the current bin.
for n in 1:bin_size
# Perform a sweep through the lattice, attemping an update to each Ising HS field.
(acceptance_rate, logdetGup, sgndetGup, logdetGdn, sgndetGdn, δG, δθ) = local_updates!(
Gup, logdetGup, sgndetGup, Gdn, logdetGdn, sgndetGdn,
hubbard_ising_parameters,
fermion_path_integral_up = fermion_path_integral_up,
fermion_path_integral_dn = fermion_path_integral_dn,
fermion_greens_calculator_up = fermion_greens_calculator_up,
fermion_greens_calculator_dn = fermion_greens_calculator_dn,
Bup = Bup, Bdn = Bdn, δG_max = δG_max, δG = δG, δθ = δθ, rng = rng
)
# Record the acceptance rate for the attempted local updates to the HS fields.
additional_info["local_acceptance_rate"] += acceptance_rate
# Make measurements, with the results being added to the measurement container.
(logdetGup, sgndetGup, logdetGdn, sgndetGdn, δG, δθ) = make_measurements!(
measurement_container,
logdetGup, sgndetGup, Gup, Gup_ττ, Gup_τ0, Gup_0τ,
logdetGdn, sgndetGdn, Gdn, Gdn_ττ, Gdn_τ0, Gdn_0τ,
fermion_path_integral_up = fermion_path_integral_up,
fermion_path_integral_dn = fermion_path_integral_dn,
fermion_greens_calculator_up = fermion_greens_calculator_up,
fermion_greens_calculator_dn = fermion_greens_calculator_dn,
Bup = Bup, Bdn = Bdn, δG_max = δG_max, δG = δG, δθ = δθ,
model_geometry = model_geometry, tight_binding_parameters = tight_binding_parameters,
coupling_parameters = (hubbard_parameters, hubbard_ising_parameters)
)
end
# Write the average measurements for the current bin to file.
write_measurements!(
measurement_container = measurement_container,
simulation_info = simulation_info,
model_geometry = model_geometry,
bin = bin,
bin_size = bin_size,
Δτ = Δτ
)
end
# Calculate acceptance rate for local updates.
additional_info["local_acceptance_rate"] /= (N_updates + N_burnin)
# Record the final numerical stabilization period that the simulation settled on.
additional_info["n_stab_final"] = fermion_greens_calculator_up.n_stab
# Record the maximum numerical error corrected by numerical stablization.
additional_info["dG"] = δG
# Write simulation summary TOML file.
save_simulation_info(simulation_info, additional_info)
#################################
### PROCESS SIMULATION RESULTS ##
#################################
# Process the simulation results, calculating final error bars for all measurements,
# writing final statisitics to CSV files.
process_measurements(simulation_info.datafolder, N_bins)
return nothing
end
# Only excute if the script is run directly from the command line.
if abspath(PROGRAM_FILE) == @__FILE__
# Read in the command line arguments.
sID = parse(Int, ARGS[1]) # simulation ID
U = parse(Float64, ARGS[2])
μ = parse(Float64, ARGS[3])
β = parse(Float64, ARGS[4])
L = parse(Int, ARGS[5])
N_burnin = parse(Int, ARGS[6])
N_updates = parse(Int, ARGS[7])
N_bins = parse(Int, ARGS[8])
# Run the simulation.
run_hubbard_chain_simulation(sID, U, μ, β, L, N_burnin, N_updates, N_bins)
end