|
| 1 | +include("../src/MeanFieldToolkit.jl") |
| 2 | +using .MeanFieldToolkit |
| 3 | +using TightBindingToolkit, FixedPointToolkit, JLD2, Plots, LaTeXStrings |
| 4 | + |
| 5 | +##### primitive vectors |
| 6 | +const a1 = [1/2, sqrt(3)/2] |
| 7 | +const a2 = [-1/2, sqrt(3)/2] |
| 8 | + |
| 9 | +const SpinVec = SpinMats(1//2) |
| 10 | +##### sublattices |
| 11 | +const b1 = [0.0, 0.0] |
| 12 | +const b2 = [0.0, 1/sqrt(3)] |
| 13 | + |
| 14 | +function honeycomb(t1::Float64, t3::Float64, inPlaneField::Float64, outPlaneField::Float64) |
| 15 | + HoppingUC = UnitCell([a1, a2], 2, 2) |
| 16 | + AddBasisSite!.(Ref(HoppingUC), [b1, b2]) |
| 17 | + |
| 18 | + tNN = Param(t1, 2) |
| 19 | + t3NN = Param(t3, 2) |
| 20 | + InPlaneField = Param(inPlaneField, 2) |
| 21 | + OutPlaneField = Param(outPlaneField, 2) |
| 22 | + |
| 23 | + AddIsotropicBonds!(tNN, HoppingUC, 1/sqrt(3), 2*SpinVec[3], "NN Hopping") |
| 24 | + AddIsotropicBonds!(t3NN, HoppingUC, 2/sqrt(3), 2*SpinVec[3], "3NN Hopping") |
| 25 | + AddIsotropicBonds!(InPlaneField, HoppingUC, 0.0, SpinVec[1], "InPlane Field") |
| 26 | + AddIsotropicBonds!(OutPlaneField, HoppingUC, 0.0, SpinVec[3], "OutPlane Field") |
| 27 | + |
| 28 | + CreateUnitCell!(HoppingUC,[tNN,t3NN,InPlaneField,OutPlaneField]) |
| 29 | + |
| 30 | + return HoppingUC |
| 31 | +end |
| 32 | + |
| 33 | +##### Thermodynamic parameters |
| 34 | +const T = 0.05 |
| 35 | +const stat = -1 |
| 36 | +const filling = 0.5 |
| 37 | + |
| 38 | +function honeycombMFT(t1::Float64, t3::Float64, inPlaneField::Float64, outPlaneField::Float64, |
| 39 | + J1::Float64, J3::Float64, |
| 40 | + fileName::String, |
| 41 | + scalings::Dict{String, Any} = Dict{String, Any}("ij" => 1.0, "ii" => 1.0, "jj" => 1.0)) |
| 42 | + HoppingUC = honeycomb(t1, t3, inPlaneField, outPlaneField) |
| 43 | + bz = BZ([33, 33]) |
| 44 | + FillBZ!(bz, HoppingUC) |
| 45 | + |
| 46 | + Jmatrix = [[1.0 0.0 0.0];[0.0 1.0 0.0];[0.0 0.0 0.0]] |
| 47 | + U = SpinToPartonCoupling(Jmatrix, 1//2) |
| 48 | + JParam3 = Param(J3, 4) |
| 49 | + JParam1 = Param(J1, 4) |
| 50 | + AddIsotropicBonds!(JParam3,HoppingUC, 2/sqrt(3), U, "J3 Interaction") |
| 51 | + AddIsotropicBonds!(JParam1,HoppingUC, 1/sqrt(3), U, "J1 Interaction") |
| 52 | + InteractionParams = [JParam3,JParam1] |
| 53 | + |
| 54 | + ##### Order parameters |
| 55 | + ferro = Param(1.0, 2) |
| 56 | + direction = 1 |
| 57 | + ##### Ferromagnetic order |
| 58 | + ferroSigns = [1, 1] |
| 59 | + for (b, basis) in enumerate(HoppingUC.basis) |
| 60 | + AddAnisotropicBond!(ferro, HoppingUC, b, b, [0, 0], ferroSigns[b] * SpinVec[direction], 0.0, "ferro order along $(direction)") |
| 61 | + end |
| 62 | + |
| 63 | + t1Chi = Param(1.0, 2) |
| 64 | + AddIsotropicBonds!(t1Chi, HoppingUC, 1/sqrt(3), 2*SpinVec[3], "t1 expectation") |
| 65 | + t3Chi = Param(1.0, 2) |
| 66 | + AddIsotropicBonds!(t3Chi, HoppingUC, 2/sqrt(3), 2*SpinVec[3], "t3 expectation") |
| 67 | + |
| 68 | + expectations = [ferro, t1Chi, t3Chi] |
| 69 | + |
| 70 | + H = Hamiltonian(HoppingUC, bz) |
| 71 | + DiagonalizeHamiltonian!(H) |
| 72 | + |
| 73 | + model = Model(HoppingUC, bz, H ; T=T, filling=filling, stat=stat) |
| 74 | + mft = TightBindingMFT(model, expectations, InteractionParams, InterQuarticToHopping, scalings) |
| 75 | + |
| 76 | + if fileName != "" |
| 77 | + sc = SolveMFT!(mft, fileName; max_iter=200, tol=1e-4); |
| 78 | + else |
| 79 | + sc = SolveMFT!(mft; max_iter=200, tol=1e-4) |
| 80 | + end |
| 81 | + # save(fileName, Dict("ferro"=> ferro.value[end], "t1" => t1Chi.value[end], "t3" => t3Chi.value[end], |
| 82 | + # "energy" => mft.MFTEnergy[end])) |
| 83 | + |
| 84 | + return mft |
| 85 | + |
| 86 | +end |
| 87 | + |
| 88 | +input = load("./J1=-1.0_J3=0.3_T=0.05_wBx_Scaling=0.0.jld2") |
| 89 | + |
| 90 | + |
| 91 | +##################### parameters |
| 92 | +const J = 1.0 |
| 93 | +const g = 1.0 |
| 94 | +const alpha1 = 0.75 |
| 95 | +const alpha3 = 0.7*alpha1 |
| 96 | + |
| 97 | +# const t1as = -0.5 * J * alpha1 * abs.(input["t1s"]) * g |
| 98 | +# const t3as = 0.2 * t1as |
| 99 | +const t1 = 0.0 |
| 100 | +const t3 = 0.0 |
| 101 | +const outPlaneField = 0.0 |
| 102 | + |
| 103 | +# J = 3.0 |
| 104 | +# theta = 0.8*pi |
| 105 | +# const J1 = J*cos(theta) |
| 106 | +# const J3 = J*sin(theta) |
| 107 | + |
| 108 | +# const J1 = -(1-alpha1)*J * g |
| 109 | +# const J3 = (1-alpha3)*J*0.25 * g |
| 110 | +const J1 = -1.0*J*g |
| 111 | +const J3 = 0.3*J*g |
| 112 | + |
| 113 | +const inPlanes = collect(range(0.0, 0.5, length=51))*J |
| 114 | +# const Bx = -1.0 |
| 115 | +# const theta = 1.0*pi |
| 116 | +# const Js = collect(LinRange(0.0, 2.0, 41)) |
| 117 | + |
| 118 | +ferros = Float64[] |
| 119 | +t1s = Float64[] |
| 120 | +t3s = Float64[] |
| 121 | +energies = Float64[] |
| 122 | + |
| 123 | +# alpha1 = 1.0 |
| 124 | +# OnSiteScaling = 1-alpha1 |
| 125 | +# scalings = Dict{String, Float64}("ij" => alpha1, "ii" => OnSiteScaling, "jj" => OnSiteScaling) |
| 126 | +scalings1 = Dict{String, Float64}("ij" => alpha1, "ii" => 1.0 - alpha1, "jj" => 1.0 - alpha1) |
| 127 | +scalings3 = Dict{String, Float64}("ij" => alpha3, "ii" => 1.0 - alpha3, "jj" => 1.0 - alpha3) |
| 128 | +scalings = Dict{String, Any}("J1 Interaction" => scalings1, "J3 Interaction" => scalings3) |
| 129 | + |
| 130 | +# for J in Js |
| 131 | +# J1 = J*cos(theta) |
| 132 | +# J3 = J*sin(theta) |
| 133 | +# fileName = "./MeanFieldToolkit.jl/Sample/HoneycombXXZ_Data/J1=$(round(J1, digits=2))_J3=$(round(J3, digits=2))_Bx=$(Bx)_T=$(T)_wBx_Scaling=$(OnSiteScaling).jld2" |
| 134 | +# sc = honeycombMFT(t1, t3, Bx, outPlaneField, J1, J3, fileName, scalings) |
| 135 | +# push!(ferros, sc.HoppingOrders[1].value[end]) |
| 136 | +# push!(t1s, sc.HoppingOrders[2].value[end]) |
| 137 | +# push!(t3s, sc.HoppingOrders[3].value[end]) |
| 138 | +# push!(energies, sc.MFTEnergy[end]) |
| 139 | +# println("J = $(J) done") |
| 140 | +# end |
| 141 | + |
| 142 | +for (b, Bx) in enumerate(inPlanes) |
| 143 | + # fileName = "./Sample/HoneycombXXZ/J1=$(round(J1, digits=2))_J3=$(round(J3, digits=2))_Bx=$(Bx)_T=$(T)_kSpace.jld2" |
| 144 | + fileName = "" |
| 145 | + sc = honeycombMFT(0.0, 0.0, Bx, outPlaneField, J1, J3, fileName, scalings) |
| 146 | + push!(ferros, sc.HoppingOrders[1].value[end]) |
| 147 | + push!(t1s, sc.HoppingOrders[2].value[end]) |
| 148 | + push!(t3s, sc.HoppingOrders[3].value[end]) |
| 149 | + push!(energies, sc.MFTEnergy[end]) |
| 150 | + println("Bx = $(Bx) done") |
| 151 | +end |
| 152 | +# save("./J1=$(round(J1, digits=2))_J3=$(round(J3, digits=2))_T=$(T)_wBx_Scaling=$(round(alpha1, digits=2))_ratio=0.5.jld2", |
| 153 | +# Dict("ferros" => ferros, "energies" => energies, "inPlanes" => inPlanes, |
| 154 | +# "t1s" => t1s, "t3s" => t3s)) |
| 155 | + |
| 156 | +# save("./J1=$(round(J1, digits=2))_J3=$(round(J3, digits=2))_T=$(T)_wBx_Scaling=$(OnSiteScaling).jld2", |
| 157 | +# Dict("ferros" => ferros, "energies" => energies, "inPlanes" => inPlanes, |
| 158 | +# "t1s" => t1s, "t3s" => t3s)) |
| 159 | + |
| 160 | +# save("./Sample/HoneycombXXZ/Dirac_J1=$(round(J1, digits=2))_J3=$(round(J3, digits=2))_T=$(T)_wBx_Scaling=$(OnSiteScaling).jld2", |
| 161 | +# Dict("ferros" => ferros, "energies" => energies, "inPlanes" => inPlanes, |
| 162 | +# "t1s" => t1s, "t3s" => t3s)) |
| 163 | + |
| 164 | +# save("./Sample/HoneycombXXZ/Dirac_J=$(J)_theta=$(round(theta/pi, digits=2))pi_T=$(T)_wBx_Scaling=$(OnSiteScaling).jld2", |
| 165 | +# Dict("ferros" => ferros, "energies" => energies, "inPlanes" => inPlanes, |
| 166 | +# "t1s" => t1s, "t3s" => t3s)) |
| 167 | + |
| 168 | +# Js = [0.0, 1.0, 2.0, 3.0] |
| 169 | + |
| 170 | +# datas = [] |
| 171 | +# for J in Js |
| 172 | +# push!(datas, load("./Sample/HoneycombXXZ/Dirac_J=$(J)_theta=$(round(theta/pi, digits=2))pi_T=$(T)_wBx_Scaling=$(OnSiteScaling).jld2")) |
| 173 | +# end |
| 174 | + |
| 175 | +# p = plot(framestyle=:box, grid=false, |
| 176 | +# guidefont = "Computer Modern", legendfont = "Computer Modern", tickfont = "Computer Modern", |
| 177 | +# guidefontsize = 14, legendfontsize = 12, tickfontsize = 12, |
| 178 | +# xlabel = L"B/t_1", |
| 179 | +# ylabel = L"M_x", |
| 180 | +# title = L"\theta=0.8\pi") |
| 181 | + |
| 182 | +# for (j, J) in enumerate(Js) |
| 183 | +# plot!(p, inPlanes, abs.(datas[j]["ferros"]), label = L"J=%$(J)", lw = 2) |
| 184 | +# end |
| 185 | + |
| 186 | + |
| 187 | + |
| 188 | +# scalings = [0.0, 0.1, 0.2, 0.3] |
| 189 | +# datas = [] |
| 190 | + |
| 191 | +# for scaling in scalings |
| 192 | +# push!(datas, load("./Sample/HoneycombXXZ/J1=$(round(J1, digits=2))_J3=$(round(J3, digits=2))_T=$(T)_wBx_Scaling=$(scaling).jld2")) |
| 193 | +# end |
| 194 | + |
| 195 | +# pMx = plot(framestyle=:box, grid=false, |
| 196 | +# guidefont = "Computer Modern", legendfont = "Computer Modern", tickfont = "Computer Modern", |
| 197 | +# guidefontsize = 14, legendfontsize = 12, tickfontsize = 12, |
| 198 | +# xlabel = L"B/J_1", |
| 199 | +# ylabel = L"M_x", |
| 200 | +# xlims = (-0.025, 1.025), |
| 201 | +# title = L"J_1=-1.0\,,J_3=0.32") |
| 202 | + |
| 203 | +# for (s, scaling) in enumerate(scalings) |
| 204 | +# plot!(pMx, inPlanes, abs.(datas[s]["ferros"]), label = L"\alpha=%$(1-scaling)", lw = 2, marker=:o) |
| 205 | +# end |
| 206 | + |
| 207 | + |
| 208 | +# pt1 = plot(framestyle=:box, grid=false, |
| 209 | +# guidefont = "Computer Modern", legendfont = "Computer Modern", tickfont = "Computer Modern", |
| 210 | +# guidefontsize = 14, legendfontsize = 12, tickfontsize = 12, |
| 211 | +# xlabel = L"B/J_1", |
| 212 | +# ylabel = L"t_1", |
| 213 | +# xlims = (-0.025, 1.025), |
| 214 | +# title = L"J_1=-1.0\,,J_3=0.32") |
| 215 | + |
| 216 | +# for (s, scaling) in enumerate(scalings) |
| 217 | +# plot!(pt1, inPlanes, abs.(datas[s]["t1s"])/2, label = L"\alpha=%$(1-scaling)", lw = 2, marker=:o) |
| 218 | +# end |
0 commit comments