diff --git a/Sample/KagomeHeis.jl b/Sample/KagomeHeis.jl
index b53e0e3..8bbd876 100644
--- a/Sample/KagomeHeis.jl
+++ b/Sample/KagomeHeis.jl
@@ -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)
 
diff --git a/Sample/RenormalizedSquaretJ.jl b/Sample/RenormalizedSquaretJ.jl
index 272ff62..b5187b1 100644
--- a/Sample/RenormalizedSquaretJ.jl
+++ b/Sample/RenormalizedSquaretJ.jl
@@ -1,18 +1,16 @@
-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)
@@ -20,31 +18,31 @@ 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
 
@@ -52,39 +50,39 @@ 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
\ No newline at end of file
diff --git a/Sample/SquareAltermagnet.jl b/Sample/SquareAltermagnet.jl
index 4e5bc83..4e58c03 100644
--- a/Sample/SquareAltermagnet.jl
+++ b/Sample/SquareAltermagnet.jl
@@ -1,98 +1,94 @@
-using MKL
+using MeanFieldToolkit, TightBindingToolkit, FixedPointToolkit
+using LinearAlgebra, Distributions, Distributions, LaTeXStrings, Base.Threads
+##### YOU NEED TO CREATE THE FOLDER Sample/Altermagnetism BEFORE RUNNING THIS SCRIPT
+########## Defining square lattice with two orbitals
+##### Primitives
+const a1 = [1.0, 1.0]
+const a2 = [1.0, -1.0]
 
-include("../src/MeanFieldToolkit.jl")
-using .MeanFieldToolkit
+const b1 = [0.0, 0.0]
+const b2 = [1.0, 0.0]
 
-using LinearAlgebra, TightBindingToolkit, Distributions, LaTeXStrings, Base.Threads
+SpinVec = SpinMats(1 // 2)
 
-##### Square lattice
-const a1  =   [1.0, 1.0]
-const a2  =   [1.0, -1.0]
-
-const b1  =   [0.0, 0.0]
-const b2  =   [1.0, 0.0] 
-
-SpinVec     =   SpinMats(1//2)
-
-const t     =   1.0
-const t1    =   0.4 * t
-const t2    =   0.2 * t
+const t = 1.0
+const t1 = 0.4 * t
+const t2 = 0.2 * t
 
 ########## Brillouin Zone
-const n       =   20
-const kSize   =   6 * n + 3  
+const n = 20
+const kSize = 6 * n + 3
 
 ##### Thermodynamic parameters
-const T         =   0.001
-const filling   =   0.5
-const stat      =   -1
+const T = 0.001
+const filling = 0.5
+const stat = -1
 
-const mixingAlpha    =    0.5   
+const mixingAlpha = 0.5
 
-Us              =   collect(LinRange(1.0 * t, 6.0 * t, 41))
+Us = collect(LinRange(1.0 * t, 6.0 * t, 41))
 
 # Threads.@threads for U in Us
-U               =   4.0
+U = 4.0
 
-HoppingUC       =   UnitCell([a1, a2], 2, 2)
-ChiUC           =   UnitCell([a1, a2], 2, 2)
-InteractionUC   =   UnitCell([a1, a2], 2, 4)
+HoppingUC = UnitCell([a1, a2], 2, 2)
+ChiUC = UnitCell([a1, a2], 2, 2)
+InteractionUC = UnitCell([a1, a2], 2, 4)
 
-InitialField    =   zeros(Float64, HoppingUC.localDim^2)
+InitialField = zeros(Float64, HoppingUC.localDim^2)
 
-AddBasisSite!(HoppingUC, b1, InitialField)
-AddBasisSite!(HoppingUC, b2, InitialField)
+AddBasisSite!(HoppingUC, b1)
+AddBasisSite!(HoppingUC, b2)
 
-AddBasisSite!(ChiUC, b1, InitialField)
-AddBasisSite!(ChiUC, b2, InitialField)
+AddBasisSite!(ChiUC, b1)
+AddBasisSite!(ChiUC, b2)
 
-AddBasisSite!(InteractionUC, b1, InitialField)
-AddBasisSite!(InteractionUC, b2, InitialField)
+AddBasisSite!(InteractionUC, b1)
+AddBasisSite!(InteractionUC, b2)
 
 ######################## HoppingParams
-tNNParam    =   Param(-t, 2)
+tNNParam = Param(-t, 2)
 AddIsotropicBonds!(tNNParam, HoppingUC, 1.0, SpinVec[4], "NN Hopping")
 
-t1Param     =   Param(-t1, 2)
-AddAnisotropicBond!(t1Param, HoppingUC, 1, 1, [ 1,  1], SpinVec[4], 2.0, "3NN Hopping t1")
-AddAnisotropicBond!(t1Param, HoppingUC, 2, 2, [ 1, -1], SpinVec[4], 2.0, "3NN Hopping t1")
+t1Param = Param(-t1, 2)
+AddAnisotropicBond!(t1Param, HoppingUC, 1, 1, [1, 1], SpinVec[4], 2.0, "3NN Hopping t1")
+AddAnisotropicBond!(t1Param, HoppingUC, 2, 2, [1, -1], SpinVec[4], 2.0, "3NN Hopping t1")
 
-t2Param     =   Param(-t2, 2)
-AddAnisotropicBond!(t2Param, HoppingUC, 1, 1, [ 1, -1], SpinVec[4], 2.0, "3NN Hopping t2")
-AddAnisotropicBond!(t2Param, HoppingUC, 2, 2, [ 1,  1], SpinVec[4], 2.0, "3NN Hopping t2")
+t2Param = Param(-t2, 2)
+AddAnisotropicBond!(t2Param, HoppingUC, 1, 1, [1, -1], SpinVec[4], 2.0, "3NN Hopping t2")
+AddAnisotropicBond!(t2Param, HoppingUC, 2, 2, [1, 1], SpinVec[4], 2.0, "3NN Hopping t2")
 
-HoppingParams   =   [tNNParam, t1Param, t2Param]
+HoppingParams = [tNNParam, t1Param, t2Param]
 CreateUnitCell!(HoppingUC, HoppingParams)
 
 ######################## Hopping expectation params
-Neel    =   Param(1.0, 2)
-AddAnisotropicBond!(Neel, ChiUC, 1, 1, [ 0, 0],  SpinVec[3], 0.0, "Neel order")
-AddAnisotropicBond!(Neel, ChiUC, 2, 2, [ 0, 0], -SpinVec[3], 0.0, "Neel order")
+Neel = Param(1.0, 2)
+AddAnisotropicBond!(Neel, ChiUC, 1, 1, [0, 0], SpinVec[3], 0.0, "Neel order")
+AddAnisotropicBond!(Neel, ChiUC, 2, 2, [0, 0], -SpinVec[3], 0.0, "Neel order")
 
-ChiParams    =  [Neel]
-HoppingBlock =  ParamBlock(ChiParams, ChiUC)
+ChiParams = [Neel]
 
 ####################### Interaction Block                                  
-n_up        =   [1.0 0.0; 0.0 0.0]
-n_down      =   [0.0 0.0; 0.0 1.0]
-Hubbard     =   DensityToPartonCoupling(n_up, n_down)
+n_up = [1.0 0.0; 0.0 0.0]
+n_down = [0.0 0.0; 0.0 1.0]
+Hubbard = DensityToPartonCoupling(n_up, n_down)
 
-UParam     =   Param(U, 4)
+UParam = Param(U, 4)
 AddIsotropicBonds!(UParam, InteractionUC, 0.0, Hubbard, "Hubbard Interaction")
-InteractionBlock   =  ParamBlock([UParam], InteractionUC)
 
-bz            =   BZ([kSize, kSize])
+
+bz = BZ([kSize, kSize])
 FillBZ!(bz, HoppingUC)
 
-H               =   Hamiltonian(HoppingUC, bz)
+H = Hamiltonian(HoppingUC, bz)
 DiagonalizeHamiltonian!(H)
 
-M    =   Model(HoppingUC, bz, H ; T=T, filling=filling, stat=stat)
+M = Model(HoppingUC, bz, H; T=T, filling=filling, stat=stat)
 SolveModel!(M)
 
-mft         =   TightBindingMFT(M, HoppingBlock, InteractionBlock, IntraQuarticToHopping)
+mft = TightBindingMFT(M, ChiParams, [UParam], IntraQuarticToHopping)
 
-fileName    =   "./Sample/Altermagnetism/U=$(round(U, digits=2))_t1=$(round(t1, digits=2))_t2=$(round(t2, digits=2))_new.jld2"
+fileName = "./Altermagnetism/U=$(round(U, digits=2))_t1=$(round(t1, digits=2))_t2=$(round(t2, digits=2))_new.jld2"
 # SolveMFT!(mft, fileName)
 SolveMFT!(mft)
 
diff --git a/Sample/SquareHubbard.jl b/Sample/SquareHubbard.jl
new file mode 100644
index 0000000..2551d1e
--- /dev/null
+++ b/Sample/SquareHubbard.jl
@@ -0,0 +1,68 @@
+using MeanFieldToolkit, TightBindingToolkit, FixedPointToolkit
+using LinearAlgebra, Distributions, Base.Threads
+##### YOU NEED TO CREATE THE FOLDER Sample/SquareHubbard_Data BEFORE RUNNING THIS SCRIPT
+########## Defining square lattice 
+##### Primitives
+const a1 = [0.0, 1.0]
+const a2 = [1.0, 0.0]
+
+const b1 = [0.0, 0.0]
+
+const InitialField = zeros(Float64, 4)
+
+const n = 10
+const kSize = 6 * n + 3
+
+SpinVec = SpinMats(1 // 2)
+const t = 1.0
+const U = 4.0
+
+##### Thermodynamic parameters
+const T = 0.001
+const stat = -1
+const mixingAlpha = 0.5
+
+
+fillings = collect(LinRange(0.1, 0.5, 20))
+UC = UnitCell([a1, a2], 2, 2)
+
+AddBasisSite!(UC, b1)
+##### HoppingParams
+t1 = -t
+t1Param = Param(t1, 2)
+AddIsotropicBonds!(t1Param, UC, 1.0, SpinVec[4], "t1")
+
+HoppingParams = [t1Param]
+n_up = [1.0 0.0; 0.0 0.0]
+n_down = [0.0 0.0; 0.0 1.0]
+Hubbard = DensityToPartonCoupling(n_up, n_down)
+
+UParam = Param(U, 4)
+AddIsotropicBonds!(UParam, UC, 0.0, Hubbard, "Hubbard Interaction")
+CreateUnitCell!(UC, HoppingParams)
+
+bz = BZ([kSize, kSize])
+FillBZ!(bz, UC)
+for filling in fillings
+    ##### Hopping expectation params
+    t_s = Param(1.0, 2)
+    AddIsotropicBonds!(t_s, UC, 1.0, SpinVec[4], "s Hopping")
+    Neel = Param(1.0, 2)
+    AddAnisotropicBond!(Neel, UC, 1, 1, [0, 0], SpinVec[3], 0.0, "Neel order")
+    AddAnisotropicBond!(Neel, UC, 2, 2, [0, 0], -SpinVec[3], 0.0, "Neel order")
+
+    ChiParams = [t_s, Neel]
+
+    ####################### Interaction Block                                  
+
+
+    H = Hamiltonian(UC, bz)
+    DiagonalizeHamiltonian!(H)
+
+    M = Model(UC, bz, H; T=T, filling=filling, stat=stat)
+    SolveModel!(M)
+    mft = TightBindingMFT(M, ChiParams, [UParam], IntraQuarticToHopping)
+    fileName = "./SquareHubbard/filling=$(round(filling, digits=3))_U=$(round(U, digits=2))_t1=$(round(t1, digits=2)).jld2"
+    SolveMFT!(mft, fileName)
+
+end
diff --git a/Sample/VortexSpins.jl b/Sample/VortexSpins.jl
index 219d7af..b701b2d 100644
--- a/Sample/VortexSpins.jl
+++ b/Sample/VortexSpins.jl
@@ -1,60 +1,62 @@
-using TightBindingToolkit, MeanFieldToolkit, LinearAlgebra
-
-############# Dual honeycomb lattice
-a1   =  [ 3/2, sqrt(3)/2]
-a2   =  [-3/2, sqrt(3)/2]
+using MeanFieldToolkit, TightBindingToolkit, FixedPointToolkit
+using LinearAlgebra, Distributions, Distributions, LaTeXStrings, Base.Threads
+##### YOU NEED TO CREATE THE FOLDER Sample/DualHoneycomb_Data BEFORE RUNNING THIS SCRIPT
+########## Defining dual honeycomb lattic
+##### Primitives
+a1 = [3 / 2, sqrt(3) / 2]
+a2 = [-3 / 2, sqrt(3) / 2]
 ##### 6 sublattices
-b1   =  [ 0.0 , 0.0 ]
-b2   =  [ 1/2 , -1/(2 * sqrt(3)) ]
-b3   =  [ 1/2 , -sqrt(3)/(2) ]
-b4   =  [ 0.0 , -2/sqrt(3)   ]
-b5   =  [-1/2 , -sqrt(3)/(2) ]
-b6   =  [-1/2 , -1/(2 * sqrt(3))]
+b1 = [0.0, 0.0]
+b2 = [1 / 2, -1 / (2 * sqrt(3))]
+b3 = [1 / 2, -sqrt(3) / (2)]
+b4 = [0.0, -2 / sqrt(3)]
+b5 = [-1 / 2, -sqrt(3) / (2)]
+b6 = [-1 / 2, -1 / (2 * sqrt(3))]
 
-firstNNdistance     = 1.0/sqrt(3)
-secondNNdistance    = 1.0
-thirdNNdistance     = 2/sqrt(3)
+firstNNdistance = 1.0 / sqrt(3)
+secondNNdistance = 1.0
+thirdNNdistance = 2 / sqrt(3)
 
-UC      =   UnitCell([a1, a2], 2, 2)
-UC.BC   =   [0.0, 0.0]
+UC = UnitCell([a1, a2], 2, 2)
+UC.BC = [0.0, 0.0]
 
 AddBasisSite!.(Ref(UC), [b1, b2, b3, b4, b5, b6])
 
-SpinVec     =   SpinMats(1//2)
+SpinVec = SpinMats(1 // 2)
 ##### XY Spin Model exchange matrix
-Jmatrix     =   [1.0 0.0 0.0 ; 0.0 1.0 0.0 ; 0.0 0.0 0.0]
+Jmatrix = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0]
 ##### XY Spin interaction written in terms of 4-parton interactions --> rank-4 array.
-U           =   SpinToPartonCoupling(Jmatrix, 1//2)
+U = SpinToPartonCoupling(Jmatrix, 1 // 2)
 
-JIntraParam     =   Param(1.0,  4)
-AddAnisotropicBond!(JIntraParam, UC, 1, 2, [ 0, 0], U, firstNNdistance, "Intra Spin Interaction")
-AddAnisotropicBond!(JIntraParam, UC, 2, 3, [ 0, 0], U, firstNNdistance, "Intra Spin Interaction")
-AddAnisotropicBond!(JIntraParam, UC, 3, 4, [ 0, 0], U, firstNNdistance, "Intra Spin Interaction")
-AddAnisotropicBond!(JIntraParam, UC, 4, 5, [ 0, 0], U, firstNNdistance, "Intra Spin Interaction")
-AddAnisotropicBond!(JIntraParam, UC, 5, 6, [ 0, 0], U, firstNNdistance, "Intra Spin Interaction")
-AddAnisotropicBond!(JIntraParam, UC, 6, 1, [ 0, 0], U, firstNNdistance, "Intra Spin Interaction")
+JIntraParam = Param(1.0, 4)
+AddAnisotropicBond!(JIntraParam, UC, 1, 2, [0, 0], U, firstNNdistance, "Intra Spin Interaction")
+AddAnisotropicBond!(JIntraParam, UC, 2, 3, [0, 0], U, firstNNdistance, "Intra Spin Interaction")
+AddAnisotropicBond!(JIntraParam, UC, 3, 4, [0, 0], U, firstNNdistance, "Intra Spin Interaction")
+AddAnisotropicBond!(JIntraParam, UC, 4, 5, [0, 0], U, firstNNdistance, "Intra Spin Interaction")
+AddAnisotropicBond!(JIntraParam, UC, 5, 6, [0, 0], U, firstNNdistance, "Intra Spin Interaction")
+AddAnisotropicBond!(JIntraParam, UC, 6, 1, [0, 0], U, firstNNdistance, "Intra Spin Interaction")
 
-JInterParam     =   Param(1.0, 4)
-AddAnisotropicBond!(JInterParam, UC, 1, 4, [ 1,  1], U, firstNNdistance, "Inter Spin Interaction")
-AddAnisotropicBond!(JInterParam, UC, 3, 6, [ 0, -1], U, firstNNdistance, "Inter Spin Interaction")
-AddAnisotropicBond!(JInterParam, UC, 5, 2, [-1,  0], U, firstNNdistance, "Inter Spin Interaction")
+JInterParam = Param(1.0, 4)
+AddAnisotropicBond!(JInterParam, UC, 1, 4, [1, 1], U, firstNNdistance, "Inter Spin Interaction")
+AddAnisotropicBond!(JInterParam, UC, 3, 6, [0, -1], U, firstNNdistance, "Inter Spin Interaction")
+AddAnisotropicBond!(JInterParam, UC, 5, 2, [-1, 0], U, firstNNdistance, "Inter Spin Interaction")
 
-Interactions   =   [JIntraParam, JInterParam]
+Interactions = [JIntraParam, JInterParam]
 
-JInters        =   collect(range(-3.0, 0.0, 16))
+JInters = collect(range(-3.0, 0.0, 16))
 ##### 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
 
 for JInter in JInters
 
@@ -62,40 +64,40 @@ for JInter in JInters
 
 
     ##### Order parameters
-    tIntraParam     =   Param(1.0, 2)
-    AddAnisotropicBond!(tIntraParam, UC, 1, 2, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
-    AddAnisotropicBond!(tIntraParam, UC, 2, 3, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
-    AddAnisotropicBond!(tIntraParam, UC, 3, 4, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
-    AddAnisotropicBond!(tIntraParam, UC, 4, 5, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
-    AddAnisotropicBond!(tIntraParam, UC, 5, 6, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
-    AddAnisotropicBond!(tIntraParam, UC, 6, 1, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
-
-    tInterParam     =   Param(1.0, 2)
-    AddAnisotropicBond!(tInterParam, UC, 1, 4, [ 1,  1], SpinVec[4], firstNNdistance, "Inter hopping")
-    AddAnisotropicBond!(tInterParam, UC, 3, 6, [ 0, -1], SpinVec[4], firstNNdistance, "Inter hopping")
-    AddAnisotropicBond!(tInterParam, UC, 5, 2, [-1,  0], SpinVec[4], firstNNdistance, "Inter hopping")
-
-    tIntraSzParam     =   Param(1.0, 2)
-    AddAnisotropicBond!(tIntraSzParam, UC, 1, 2, [ 0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
-    AddAnisotropicBond!(tIntraSzParam, UC, 2, 3, [ 0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
-    AddAnisotropicBond!(tIntraSzParam, UC, 3, 4, [ 0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
-    AddAnisotropicBond!(tIntraSzParam, UC, 4, 5, [ 0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
-    AddAnisotropicBond!(tIntraSzParam, UC, 5, 6, [ 0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
-    AddAnisotropicBond!(tIntraSzParam, UC, 6, 1, [ 0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
-
-    tInterSzParam     =   Param(1.0, 2)
-    AddAnisotropicBond!(tInterSzParam, UC, 1, 4, [ 1,  1], SpinVec[3], firstNNdistance, "Inter hopping Sz")
-    AddAnisotropicBond!(tInterSzParam, UC, 3, 6, [ 0, -1], SpinVec[3], firstNNdistance, "Inter hopping Sz")
-    AddAnisotropicBond!(tInterSzParam, UC, 5, 2, [-1,  0], SpinVec[3], firstNNdistance, "Inter hopping Sz")
-
-    HoppingOrders   =   [tIntraParam, tInterParam, tIntraSzParam, tInterSzParam]
-
-    H             =   Hamiltonian(UC, bz)
+    tIntraParam = Param(1.0, 2)
+    AddAnisotropicBond!(tIntraParam, UC, 1, 2, [0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
+    AddAnisotropicBond!(tIntraParam, UC, 2, 3, [0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
+    AddAnisotropicBond!(tIntraParam, UC, 3, 4, [0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
+    AddAnisotropicBond!(tIntraParam, UC, 4, 5, [0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
+    AddAnisotropicBond!(tIntraParam, UC, 5, 6, [0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
+    AddAnisotropicBond!(tIntraParam, UC, 6, 1, [0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
+
+    tInterParam = Param(1.0, 2)
+    AddAnisotropicBond!(tInterParam, UC, 1, 4, [1, 1], SpinVec[4], firstNNdistance, "Inter hopping")
+    AddAnisotropicBond!(tInterParam, UC, 3, 6, [0, -1], SpinVec[4], firstNNdistance, "Inter hopping")
+    AddAnisotropicBond!(tInterParam, UC, 5, 2, [-1, 0], SpinVec[4], firstNNdistance, "Inter hopping")
+
+    tIntraSzParam = Param(1.0, 2)
+    AddAnisotropicBond!(tIntraSzParam, UC, 1, 2, [0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
+    AddAnisotropicBond!(tIntraSzParam, UC, 2, 3, [0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
+    AddAnisotropicBond!(tIntraSzParam, UC, 3, 4, [0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
+    AddAnisotropicBond!(tIntraSzParam, UC, 4, 5, [0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
+    AddAnisotropicBond!(tIntraSzParam, UC, 5, 6, [0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
+    AddAnisotropicBond!(tIntraSzParam, UC, 6, 1, [0, 0], SpinVec[3], firstNNdistance, "Intra hopping Sz")
+
+    tInterSzParam = Param(1.0, 2)
+    AddAnisotropicBond!(tInterSzParam, UC, 1, 4, [1, 1], SpinVec[3], firstNNdistance, "Inter hopping Sz")
+    AddAnisotropicBond!(tInterSzParam, UC, 3, 6, [0, -1], SpinVec[3], firstNNdistance, "Inter hopping Sz")
+    AddAnisotropicBond!(tInterSzParam, UC, 5, 2, [-1, 0], SpinVec[3], firstNNdistance, "Inter hopping Sz")
+
+    HoppingOrders = [tIntraParam, tInterParam, tIntraSzParam, tInterSzParam]
+
+    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, InterQuarticToHopping, Dict{String, Float64}("ij" => 1.0, "ii" => 0.0, "jj" => 0.0) )
-    fileName      =   "./Sample/DualHoneycomb_Data/JInter=$(round(JInter, digits=2))_JIntra=$(round(JIntraParam.value[end], digits=2)).jld2"
+    mft = TightBindingMFT(M, HoppingOrders, Interactions, InterQuarticToHopping, Dict{String,Float64}("ij" => 1.0, "ii" => 0.0, "jj" => 0.0))
+    fileName = "./DualHoneycomb_Data/JInter=$(round(JInter, digits=2))_JIntra=$(round(JIntraParam.value[end], digits=2)).jld2"
     SolveMFT!(mft, fileName)
 
 end