carstenbauer / MonteCarlo.jl

Classical and quantum Monte Carlo simulations in Julia
https://carstenbauer.github.io/MonteCarlo.jl/dev/
Other
185 stars 18 forks source link

different interaction for multi-orbital model? #174

Open kfkq opened 1 year ago

kfkq commented 1 year ago

I am not sure whether we can implement a multi-orbital model with its specified interaction in current version of package. For example, given orbital A at a particular lattice position, it has its own μ_A, t_AA for same-orbital hopping, t_AB for different-orbital hopping, and U_A for Hubbard interaction. It is easy to implement the hopping matrix, but I am unsure about the interaction if I have more than one orbital {A, B, C}.

Due to a lack of documentation, it is difficult to figure out whether this feature is already part of the package or not. However, from what I can tell by reading the source code, it seems that this feature is not yet implemented and only supports Hubbard interaction for a single orbital. It does not require a lot of modification though. I only need to define my own interaction_matrix_exp and propose_local for different orbitals.

So, I need to construct the U matrix, specifically, the α.

α[i] = acosh(exp(-0.5 * param.delta_tau * ComplexF64(U[i])))

where i is positions of sites.

then finally, the interaction term need to be written to be like this

@inline function interaction_matrix_exp!(f::MagneticHirschField, result::Diagonal, slice, power)
    N = size(f.conf, 1)
    @inbounds for i in axes(f.conf, 1)
        result.diag[i]   = exp( power * α[i] * f.conf[i, slice])
    end
    @inbounds for i in axes(f.conf, 1)
        result.diag[i+N] = exp(-power * α[i] * f.conf[i, slice])
    end
    nothing
end

function propose_local(mc, f::MagneticHirschField, i, slice)
    ΔE_Boson = -2.0 * α[i] * f.conf[i, slice]
    mc.stack.field_cache.Δ[1] = exp(+ΔE_Boson) - 1.0
    mc.stack.field_cache.Δ[2] = exp(-ΔE_Boson) - 1.0
    detratio = calculate_detratio!(mc.stack.field_cache, mc.stack.greens, i)

    # There is no bosonic part (exp(-ΔE_Boson)) to the partition function.
    # Therefore pass 0.0
    return detratio, 0.0, nothing
end

let me know if this is correct approach.