Skip to content

Commit

Permalink
Merge pull request #2 from andrewkhardy/main
Browse files Browse the repository at this point in the history
updated the Sample codes and made them more uniform
  • Loading branch information
Anjishnubose authored Dec 19, 2023
2 parents ca7bb84 + b74c561 commit 7f173a1
Show file tree
Hide file tree
Showing 5 changed files with 287 additions and 222 deletions.
105 changes: 53 additions & 52 deletions Sample/KagomeHeis.jl
Original file line number Diff line number Diff line change
@@ -1,83 +1,84 @@
using MeanFieldToolkit, LinearAlgebra, TightBindingToolkit, Distributions, LaTeXStrings

using MeanFieldToolkit, TightBindingToolkit, FixedPointToolkit
using LinearAlgebra, Distributions, Distributions, LaTeXStrings, Base.Threads
##### YOU NEED TO CREATE THE FOLDER Sample/KagomeHeis_Data BEFORE RUNNING THIS SCRIPT
########## Defining Kagome lattice with the doubled unit cell
##### Primitives
const a1 = [4.0, 0.0]
const a2 = [1.0, sqrt(3)]
const a1 = [4.0, 0.0]
const a2 = [1.0, sqrt(3)]
##### 6 sublattices
const b1 = [0.0, 0.0]
const b2 = [1.0, 0.0]
const b3 = [0.5, sqrt(3)/2]
const b4 = [2.0 + 0.0, 0.0]
const b5 = [2.0 + 1.0, 0.0]
const b6 = [2.0 + 0.5, sqrt(3)/2]
const b1 = [0.0, 0.0]
const b2 = [1.0, 0.0]
const b3 = [0.5, sqrt(3) / 2]
const b4 = [2.0 + 0.0, 0.0]
const b5 = [2.0 + 1.0, 0.0]
const b6 = [2.0 + 0.5, sqrt(3) / 2]
##### On-site spin matrices
const InitialField = zeros(Float64, 4)
const OnSite = SpinMats(1//2)
const InitialField = zeros(Float64, 4)
const OnSite = SpinMats(1 // 2)
##### Unit cell with local dimensions of 2, and tracking rank-2 bonds
UC = UnitCell([a1, a2], 2, 2)
AddBasisSite!.(Ref(UC), [b1, b2, b3, b4, b5, b6], Ref(InitialField) , Ref(OnSite))
UC = UnitCell([a1, a2], 2, 2)
AddBasisSite!.(Ref(UC), [b1, b2, b3, b4, b5, b6], Ref(InitialField), Ref(OnSite))
##### Strength of the Heisenberg interaction
const J = 1.0
const J = 1.0
##### Heisenberg spin exchange matrix
Jmatrix = Matrix{Float64}(I, 3, 3)
Jmatrix = Matrix{Float64}(I, 3, 3)
##### Heisenberg interaction written in terms of 4-parton interactions --> rank-4 array.
U = SpinToPartonCoupling(Jmatrix, 1//2)
U = SpinToPartonCoupling(Jmatrix, 1 // 2)
###### Interaction parameter which is isotropic.
JParam = Param(J, 4)
JParam = Param(J, 4)
AddIsotropicBonds!(JParam, UC, 1.0, U, "Heisenberg Interaction")
Interactions= [JParam]
Interactions = [JParam]
##### Brillouin zone
const n = 10
const kSize = 6 * n + 3
bz = BZ([kSize, kSize])
const n = 10
const kSize = 6 * n + 3
bz = BZ([kSize, kSize])
FillBZ!(bz, UC)
##### Thermodynamic parameters
const T = 0.001
const filling = 0.5
const stat = -1
const T = 0.001
const filling = 0.5
const stat = -1
##### Mixing alpha used for self-consistency solver
const mixingAlpha = 0.5
const mixingAlpha = 0.5
##### Different flux configuration ansatzes to be tested.
phis = collect(LinRange(-pi/2, pi/2, 21))
phis = collect(LinRange(-pi / 2, pi / 2, 21))

for ϕ in phis

##### Hopping expectation params
##### Hopping with flux ϕ through the triangles, and π-2ϕ through the hexagons
t_flux = Param(1.0, 2)
AddAnisotropicBond!(t_flux, UC, 1, 2, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 2, 3, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 3, 1, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 4, 5, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 5, 6, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 6, 4, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 4, 2, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
t_flux = Param(1.0, 2)
AddAnisotropicBond!(t_flux, UC, 1, 2, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 2, 3, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 3, 1, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 4, 5, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 5, 6, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 6, 4, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 4, 2, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")

AddAnisotropicBond!(t_flux, UC, 2, 6, [ 0, -1], exp( im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 4, 6, [ 0, -1], exp(-im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 5, 3, [ 1, -1], exp( im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 5, 1, [ 1, 0], exp(-im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 3, 1, [ 0, 1], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 2, 6, [0, -1], exp(im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 4, 6, [0, -1], exp(-im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 5, 3, [1, -1], exp(im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 5, 1, [1, 0], exp(-im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")
AddAnisotropicBond!(t_flux, UC, 3, 1, [0, 1], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)")

######### Weiss fields for ordering
Ordering = Param(1.0, 2)
AddAnisotropicBond!(Ordering, UC, 1, 1, [ 0, 0], sum([ -sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 2, 2, [ 0, 0], sum([ sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 3, 3, [ 0, 0], sum([ 0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 4, 4, [ 0, 0], sum([ -sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 5, 5, [ 0, 0], sum([ sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 6, 6, [ 0, 0], sum([ 0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
Ordering = Param(1.0, 2)
AddAnisotropicBond!(Ordering, UC, 1, 1, [0, 0], sum([-sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 2, 2, [0, 0], sum([sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 3, 3, [0, 0], sum([0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 4, 4, [0, 0], sum([-sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 5, 5, [0, 0], sum([sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")
AddAnisotropicBond!(Ordering, UC, 6, 6, [0, 0], sum([0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering")

HoppingOrders = [t_flux, Ordering]
HoppingOrders = [t_flux, Ordering]
##### Build the model (trivial since we are looking at the pure spin Heisenberg model)
H = Hamiltonian(UC, bz)
H = Hamiltonian(UC, bz)
DiagonalizeHamiltonian!(H)
M = Model(UC, bz, H ; T=T, filling=filling, stat=stat)
M = Model(UC, bz, H; T=T, filling=filling, stat=stat)
##### Build the mean-field theory object, with chosen scaling such that on-site ordering is suppressed.
mft = TightBindingMFT(M, HoppingOrders, Interactions, Function[InterQuarticToHopping], Dict{String, Float64}("ij" => 1.0, "ii" => 0.0, "jj" => 0.0) )
mft = TightBindingMFT(M, HoppingOrders, Interactions, Function[InterQuarticToHopping], Dict{String,Float64}("ij" => 1.0, "ii" => 0.0, "jj" => 0.0))
##### File to save data to
fileName = "./KagomeHeis_Data/J=$(round(J, digits=3))_phi=$(round/pi, digits=3))Pi_New.jld2"
fileName = "./KagomeHeis_Data/J=$(round(J, digits=3))_phi=$(round/pi, digits=3))Pi_New.jld2"
##### Solve the mean-field theory and save the results in fileName
SolveMFT!(mft, fileName)

Expand Down
82 changes: 40 additions & 42 deletions Sample/RenormalizedSquaretJ.jl
Original file line number Diff line number Diff line change
@@ -1,90 +1,88 @@
include("../src/MeanFieldToolkit.jl")
using .MeanFieldToolkit

using LinearAlgebra, TightBindingToolkit, Distributions, Base.Threads

using MeanFieldToolkit, TightBindingToolkit, FixedPointToolkit
using LinearAlgebra, Distributions, Distributions, LaTeXStrings, Base.Threads
##### YOU NEED TO CREATE THE FOLDER Sample/SquaretJ_Data BEFORE RUNNING THIS SCRIPT
################### Square lattice with the doubles Unit Cell
##### Primitive vectors
const a1 = [1.0, 1.0]
const a2 = [1.0, -1.0]
##### Primitives
const a1 = [1.0, 1.0]
const a2 = [1.0, -1.0]

const b1 = [0.0, 0.0]
const b2 = [1.0, 0.0]
const b1 = [0.0, 0.0]
const b2 = [1.0, 0.0]

HoppingUC = UnitCell([a1, a2], 2, 2)
PairingUC = UnitCell([a1, a2], 2, 2)
HoppingUC = UnitCell([a1, a2], 2, 2)
PairingUC = UnitCell([a1, a2], 2, 2)

AddBasisSite!(HoppingUC, b1)
AddBasisSite!(HoppingUC, b2)

AddBasisSite!(PairingUC, b1)
AddBasisSite!(PairingUC, b2)

SpinVec = SpinMats(1//2)
SpinVec = SpinMats(1 // 2)

##### HoppingParams
const t = 1.0
t1Param = Param(t, 2)
const t = 1.0
t1Param = Param(t, 2)
AddIsotropicBonds!(t1Param, HoppingUC, 1.0, SpinVec[4], "t1")

const J = t/5
Jmatrix = Matrix{Float64}(I, 3, 3)
U = SpinToPartonCoupling(Jmatrix, 1//2)
JParam = Param(J, 4)
const J = t / 5
Jmatrix = Matrix{Float64}(I, 3, 3)
U = SpinToPartonCoupling(Jmatrix, 1 // 2)
JParam = Param(J, 4)
AddIsotropicBonds!(JParam, HoppingUC, 1.0, U, "Heisenberg Interaction")

const n = 10
const kSize = 6 * n + 3
bz = BZ([kSize, kSize])
FillBZ!(bz, HoppingUC)
const n = 10
const kSize = 6 * n + 3
bz = BZ([kSize, kSize])
FillBZ!(bz, HoppingUC)



##### Thermodynamic parameters
const T = 0.001
const stat = -1
const T = 0.001
const stat = -1

const mixingAlpha = 0.5
const mixingAlpha = 0.5

##### Renormalized MFT from https://arxiv.org/pdf/cond-mat/0311604.pdf

fillings = collect(LinRange(0.2, 0.49, 2))

for filling in fillings

delta = 1-2*filling
rnorm_hopping_factor = 2*delta/(1+delta)
rnorm_int_factor = 4/((1+delta)^2)
delta = 1 - 2 * filling
rnorm_hopping_factor = 2 * delta / (1 + delta)
rnorm_int_factor = 4 / ((1 + delta)^2)

push!(t1Param.value, -t * rnorm_hopping_factor)
push!(JParam.value, J * rnorm_int_factor)

ModifyUnitCell!(HoppingUC, [t1Param])
Interactions = [JParam]
Interactions = [JParam]

##### Hopping expectation params
t_s = Param(1.0, 2)
t_s = Param(1.0, 2)
AddIsotropicBonds!(t_s, HoppingUC, 1.0, SpinVec[4], "s Hopping")

HoppingOrders = [t_s]
HoppingOrders = [t_s]

##### Pairing expectation params
p_dx2y2 = Param(1.0, 2)
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [ 0, 0], (SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [-1 , -1], (SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [ 0, -1], -(SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [ -1, 0], -(SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")
p_dx2y2 = Param(1.0, 2)
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [0, 0], (SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [-1, -1], (SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [0, -1], -(SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")
AddAnisotropicBond!(p_dx2y2, PairingUC, 1, 2, [-1, 0], -(SpinVec[2]), 1.0, "d_x^2-y^2 Pairing")

PairingOrders = [p_dx2y2]
PairingOrders = [p_dx2y2]

bdgH = Hamiltonian(HoppingUC, PairingUC, bz)
bdgH = Hamiltonian(HoppingUC, PairingUC, bz)
DiagonalizeHamiltonian!(bdgH)

bdgModel = BdGModel(HoppingUC, PairingUC, bz, bdgH ; T=T, filling=filling, stat=stat)
bdgModel = BdGModel(HoppingUC, PairingUC, bz, bdgH; T=T, filling=filling, stat=stat)
SolveModel!(bdgModel)

bdgmft = BdGMFT(bdgModel, HoppingOrders, PairingOrders, Interactions, InterQuarticToHopping, InterQuarticToPairing)
fileName = "./Sample/SquaretJ_Data/filling=$(round(filling, digits=3))_t1=$(round(t1Param.value[end], digits=3))_J=$(round(J, digits=3))_wtWeiss.jld2"
bdgmft = BdGMFT(bdgModel, HoppingOrders, PairingOrders, Interactions, InterQuarticToHopping, InterQuarticToPairing)
fileName = "./SquaretJ_Data/filling=$(round(filling, digits=3))_t1=$(round(t1Param.value[end], digits=3))_J=$(round(J, digits=3))_wtWeiss.jld2"
SolveMFT!(bdgmft, fileName)

end
Loading

0 comments on commit 7f173a1

Please sign in to comment.