|
1 |
| -using MeanFieldToolkit, LinearAlgebra, TightBindingToolkit, Distributions, LaTeXStrings |
2 |
| - |
| 1 | +using MeanFieldToolkit, TightBindingToolkit, FixedPointToolkit |
| 2 | +using LinearAlgebra, Distributions, Distributions, LaTeXStrings, Base.Threads |
| 3 | +##### YOU NEED TO CREATE THE FOLDER Sample/KagomeHeis_Data BEFORE RUNNING THIS SCRIPT |
3 | 4 | ########## Defining Kagome lattice with the doubled unit cell
|
4 | 5 | ##### Primitives
|
5 |
| -const a1 = [4.0, 0.0] |
6 |
| -const a2 = [1.0, sqrt(3)] |
| 6 | +const a1 = [4.0, 0.0] |
| 7 | +const a2 = [1.0, sqrt(3)] |
7 | 8 | ##### 6 sublattices
|
8 |
| -const b1 = [0.0, 0.0] |
9 |
| -const b2 = [1.0, 0.0] |
10 |
| -const b3 = [0.5, sqrt(3)/2] |
11 |
| -const b4 = [2.0 + 0.0, 0.0] |
12 |
| -const b5 = [2.0 + 1.0, 0.0] |
13 |
| -const b6 = [2.0 + 0.5, sqrt(3)/2] |
| 9 | +const b1 = [0.0, 0.0] |
| 10 | +const b2 = [1.0, 0.0] |
| 11 | +const b3 = [0.5, sqrt(3) / 2] |
| 12 | +const b4 = [2.0 + 0.0, 0.0] |
| 13 | +const b5 = [2.0 + 1.0, 0.0] |
| 14 | +const b6 = [2.0 + 0.5, sqrt(3) / 2] |
14 | 15 | ##### On-site spin matrices
|
15 |
| -const InitialField = zeros(Float64, 4) |
16 |
| -const OnSite = SpinMats(1//2) |
| 16 | +const InitialField = zeros(Float64, 4) |
| 17 | +const OnSite = SpinMats(1 // 2) |
17 | 18 | ##### Unit cell with local dimensions of 2, and tracking rank-2 bonds
|
18 |
| -UC = UnitCell([a1, a2], 2, 2) |
19 |
| -AddBasisSite!.(Ref(UC), [b1, b2, b3, b4, b5, b6], Ref(InitialField) , Ref(OnSite)) |
| 19 | +UC = UnitCell([a1, a2], 2, 2) |
| 20 | +AddBasisSite!.(Ref(UC), [b1, b2, b3, b4, b5, b6], Ref(InitialField), Ref(OnSite)) |
20 | 21 | ##### Strength of the Heisenberg interaction
|
21 |
| -const J = 1.0 |
| 22 | +const J = 1.0 |
22 | 23 | ##### Heisenberg spin exchange matrix
|
23 |
| -Jmatrix = Matrix{Float64}(I, 3, 3) |
| 24 | +Jmatrix = Matrix{Float64}(I, 3, 3) |
24 | 25 | ##### Heisenberg interaction written in terms of 4-parton interactions --> rank-4 array.
|
25 |
| -U = SpinToPartonCoupling(Jmatrix, 1//2) |
| 26 | +U = SpinToPartonCoupling(Jmatrix, 1 // 2) |
26 | 27 | ###### Interaction parameter which is isotropic.
|
27 |
| -JParam = Param(J, 4) |
| 28 | +JParam = Param(J, 4) |
28 | 29 | AddIsotropicBonds!(JParam, UC, 1.0, U, "Heisenberg Interaction")
|
29 |
| -Interactions= [JParam] |
| 30 | +Interactions = [JParam] |
30 | 31 | ##### Brillouin zone
|
31 |
| -const n = 10 |
32 |
| -const kSize = 6 * n + 3 |
33 |
| -bz = BZ([kSize, kSize]) |
| 32 | +const n = 10 |
| 33 | +const kSize = 6 * n + 3 |
| 34 | +bz = BZ([kSize, kSize]) |
34 | 35 | FillBZ!(bz, UC)
|
35 | 36 | ##### Thermodynamic parameters
|
36 |
| -const T = 0.001 |
37 |
| -const filling = 0.5 |
38 |
| -const stat = -1 |
| 37 | +const T = 0.001 |
| 38 | +const filling = 0.5 |
| 39 | +const stat = -1 |
39 | 40 | ##### Mixing alpha used for self-consistency solver
|
40 |
| -const mixingAlpha = 0.5 |
| 41 | +const mixingAlpha = 0.5 |
41 | 42 | ##### Different flux configuration ansatzes to be tested.
|
42 |
| -phis = collect(LinRange(-pi/2, pi/2, 21)) |
| 43 | +phis = collect(LinRange(-pi / 2, pi / 2, 21)) |
43 | 44 |
|
44 | 45 | for ϕ in phis
|
45 | 46 |
|
46 | 47 | ##### Hopping expectation params
|
47 | 48 | ##### Hopping with flux ϕ through the triangles, and π-2ϕ through the hexagons
|
48 |
| - t_flux = Param(1.0, 2) |
49 |
| - AddAnisotropicBond!(t_flux, UC, 1, 2, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
50 |
| - AddAnisotropicBond!(t_flux, UC, 2, 3, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
51 |
| - AddAnisotropicBond!(t_flux, UC, 3, 1, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
52 |
| - AddAnisotropicBond!(t_flux, UC, 4, 5, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
53 |
| - AddAnisotropicBond!(t_flux, UC, 5, 6, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
54 |
| - AddAnisotropicBond!(t_flux, UC, 6, 4, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
55 |
| - AddAnisotropicBond!(t_flux, UC, 4, 2, [ 0, 0], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 49 | + t_flux = Param(1.0, 2) |
| 50 | + AddAnisotropicBond!(t_flux, UC, 1, 2, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 51 | + AddAnisotropicBond!(t_flux, UC, 2, 3, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 52 | + AddAnisotropicBond!(t_flux, UC, 3, 1, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 53 | + AddAnisotropicBond!(t_flux, UC, 4, 5, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 54 | + AddAnisotropicBond!(t_flux, UC, 5, 6, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 55 | + AddAnisotropicBond!(t_flux, UC, 6, 4, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 56 | + AddAnisotropicBond!(t_flux, UC, 4, 2, [0, 0], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
56 | 57 |
|
57 |
| - AddAnisotropicBond!(t_flux, UC, 2, 6, [ 0, -1], exp( im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
58 |
| - AddAnisotropicBond!(t_flux, UC, 4, 6, [ 0, -1], exp(-im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
59 |
| - AddAnisotropicBond!(t_flux, UC, 5, 3, [ 1, -1], exp( im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
60 |
| - AddAnisotropicBond!(t_flux, UC, 5, 1, [ 1, 0], exp(-im * (pi + ϕ/3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
61 |
| - AddAnisotropicBond!(t_flux, UC, 3, 1, [ 0, 1], exp( im * ϕ/3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 58 | + AddAnisotropicBond!(t_flux, UC, 2, 6, [0, -1], exp(im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 59 | + AddAnisotropicBond!(t_flux, UC, 4, 6, [0, -1], exp(-im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 60 | + AddAnisotropicBond!(t_flux, UC, 5, 3, [1, -1], exp(im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 61 | + AddAnisotropicBond!(t_flux, UC, 5, 1, [1, 0], exp(-im * (pi + ϕ / 3)) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
| 62 | + AddAnisotropicBond!(t_flux, UC, 3, 1, [0, 1], exp(im * ϕ / 3) * OnSite[4], 1.0, "[ϕ, ϕ, π - 2*ϕ] Hopping with ϕ/π = $(ϕ/pi)") |
62 | 63 |
|
63 | 64 | ######### Weiss fields for ordering
|
64 |
| - Ordering = Param(1.0, 2) |
65 |
| - AddAnisotropicBond!(Ordering, UC, 1, 1, [ 0, 0], sum([ -sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
66 |
| - AddAnisotropicBond!(Ordering, UC, 2, 2, [ 0, 0], sum([ sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
67 |
| - AddAnisotropicBond!(Ordering, UC, 3, 3, [ 0, 0], sum([ 0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
68 |
| - AddAnisotropicBond!(Ordering, UC, 4, 4, [ 0, 0], sum([ -sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
69 |
| - AddAnisotropicBond!(Ordering, UC, 5, 5, [ 0, 0], sum([ sqrt(3)/2, -1/2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
70 |
| - AddAnisotropicBond!(Ordering, UC, 6, 6, [ 0, 0], sum([ 0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
| 65 | + Ordering = Param(1.0, 2) |
| 66 | + AddAnisotropicBond!(Ordering, UC, 1, 1, [0, 0], sum([-sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
| 67 | + AddAnisotropicBond!(Ordering, UC, 2, 2, [0, 0], sum([sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
| 68 | + AddAnisotropicBond!(Ordering, UC, 3, 3, [0, 0], sum([0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
| 69 | + AddAnisotropicBond!(Ordering, UC, 4, 4, [0, 0], sum([-sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
| 70 | + AddAnisotropicBond!(Ordering, UC, 5, 5, [0, 0], sum([sqrt(3) / 2, -1 / 2, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
| 71 | + AddAnisotropicBond!(Ordering, UC, 6, 6, [0, 0], sum([0.0, 1.0, 0.0] .* OnSite[1:3]), 0.0, "120 degree ordering") |
71 | 72 |
|
72 |
| - HoppingOrders = [t_flux, Ordering] |
| 73 | + HoppingOrders = [t_flux, Ordering] |
73 | 74 | ##### Build the model (trivial since we are looking at the pure spin Heisenberg model)
|
74 |
| - H = Hamiltonian(UC, bz) |
| 75 | + H = Hamiltonian(UC, bz) |
75 | 76 | DiagonalizeHamiltonian!(H)
|
76 |
| - M = Model(UC, bz, H ; T=T, filling=filling, stat=stat) |
| 77 | + M = Model(UC, bz, H; T=T, filling=filling, stat=stat) |
77 | 78 | ##### Build the mean-field theory object, with chosen scaling such that on-site ordering is suppressed.
|
78 |
| - mft = TightBindingMFT(M, HoppingOrders, Interactions, Function[InterQuarticToHopping], Dict{String, Float64}("ij" => 1.0, "ii" => 0.0, "jj" => 0.0) ) |
| 79 | + mft = TightBindingMFT(M, HoppingOrders, Interactions, Function[InterQuarticToHopping], Dict{String,Float64}("ij" => 1.0, "ii" => 0.0, "jj" => 0.0)) |
79 | 80 | ##### File to save data to
|
80 |
| - fileName = "./KagomeHeis_Data/J=$(round(J, digits=3))_phi=$(round(ϕ/pi, digits=3))Pi_New.jld2" |
| 81 | + fileName = "./KagomeHeis_Data/J=$(round(J, digits=3))_phi=$(round(ϕ/pi, digits=3))Pi_New.jld2" |
81 | 82 | ##### Solve the mean-field theory and save the results in fileName
|
82 | 83 | SolveMFT!(mft, fileName)
|
83 | 84 |
|
|
0 commit comments