From da403eb2ac45a0cdcabfb498a63f6a8cbc9e6891 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Wed, 27 Sep 2023 18:23:12 -0400
Subject: [PATCH 01/14] added function to check whether given line intersects a
 given line segment

---
 src/Useful.jl | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/src/Useful.jl b/src/Useful.jl
index 3d74ce5..e3fbab3 100644
--- a/src/Useful.jl
+++ b/src/Useful.jl
@@ -1,5 +1,5 @@
 module Useful
-    export GetAllOffsets, VecAngle, Meshgrid, BinarySearch, DistFunction, DeriDistFunction , GetIndexPath, FFTArrayofMatrix, Central_Diff, Arrayfy, DeArrayfy, GetPhase
+    export GetAllOffsets, VecAngle, Meshgrid, BinarySearch, DistFunction, DeriDistFunction , GetIndexPath, FFTArrayofMatrix, Central_Diff, Arrayfy, DeArrayfy, GetPhase, SegmentIntersection
 
     using LinearAlgebra, Statistics, FFTW, TensorCast
 
@@ -255,4 +255,20 @@ PBC   : vector of boolean with length = spatial dimensions = rank of f ---> if t
     end
 
 
+    function SegmentIntersection(SegmentPoints::Vector{Vector{Float64}}, LinePoints::Vector{Vector{Float64}})::Float64
+
+        M1, M2  =   SegmentPoints
+        r1, r2  =   LinePoints
+
+        dx  =   (r2[1] - r1[1])
+        dy  =   (r2[2] - r1[2])
+        dMx =   (M2[1] - M1[1])
+        dMy =   (M2[2] - M1[2])
+
+        t       =   (dx * (r1[2] - M1[2]) - dy * (r1[1] - M1[1])) / (dx * dMy - dy * dMx)
+
+        return t
+    end
+
+
 end
\ No newline at end of file

From 131e5306f2ea3bf4febaa70722cfc7a73a167b59 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Wed, 27 Sep 2023 18:23:34 -0400
Subject: [PATCH 02/14] finding rotation eigenvalues

---
 Sample/Vortex.jl | 73 +++++++++++++++++++++++++++++++++++-------------
 1 file changed, 54 insertions(+), 19 deletions(-)

diff --git a/Sample/Vortex.jl b/Sample/Vortex.jl
index 481dd07..972de33 100644
--- a/Sample/Vortex.jl
+++ b/Sample/Vortex.jl
@@ -18,13 +18,14 @@ secondNNdistance    = 1.0
 thirdNNdistance     = 2/sqrt(3)
 
 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)
 
 const tIntra  =     1.0
-const tInter  =     -2.5
+const tInter  =     -1.5
 
 tIntraParam     =   Param(tIntra, 2)
 AddAnisotropicBond!(tIntraParam, UC, 1, 2, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
@@ -47,35 +48,69 @@ CreateUnitCell!(UC, [tIntraParam, tInterParam])
 
 # Plot_UnitCell!(UC ; bond_cmp=:viridis, site_size=10.0, plot_conjugate=false, plot_lattice=true)
 
-# kSize   =   6*50 + 3
-# bz      =   BZ([kSize, kSize])
-# FillBZ!(bz, UC)
+kSize   =   6*50 + 3
+bz      =   BZ([kSize, kSize])
+FillBZ!(bz, UC)
+
+H       =   Hamiltonian(UC, bz)
+DiagonalizeHamiltonian!(H)
+
+VortexModel     =   Model(UC, bz, H ; T=0.001, mu=0.0)
+SolveModel!(VortexModel ; get_gap = true)
+
+R3      =   zeros(Float64, 6, 6)
+R3[3, 1]=   1.0
+R3[4, 2]=   1.0
+R3[5, 3]=   1.0
+R3[6, 4]=   1.0
+R3[1, 5]=   1.0
+R3[2, 6]=   1.0
+R3      =   kron(R3, Matrix(I, 2, 2))
+
+R2      =   zeros(Float64, 6, 6)
+R2[4, 1]=   1.0
+R2[5, 2]=   1.0
+R2[6, 3]=   1.0
+R2[1, 4]=   1.0
+R2[2, 5]=   1.0
+R2[3, 6]=   1.0
+R2      =   kron(R2, Matrix(I, 2, 2))
+
+Gammaind=   GetQIndex(bz.HighSymPoints["G"], bz)
+K1ind   =   GetQIndex(bz.HighSymPoints["K1"], bz)
+K2ind   =   GetQIndex(bz.HighSymPoints["K2"], bz)
+
+OrbitalsGamma   =   H.states[Gammaind...]
+OrbitalsK1  =   H.states[K1ind...]
+OrbitalsK2  =   H.states[K2ind...]
+
+REigGamma   =   eigen((OrbitalsGamma' * R3 * OrbitalsGamma)[1:6, 1:6]).values
+REigK1      =   eigen((OrbitalsK1' * R3 * OrbitalsK1)[1:6, 1:6]).values
+REigK2      =   eigen((OrbitalsK2' * R3 * OrbitalsK2)[1:6, 1:6]).values
 
-# H       =   Hamiltonian(UC, bz)
-# DiagonalizeHamiltonian!(H)
 
-# VortexModel     =   Model(UC, bz, H ; T=0.001, mu=0.0)
-# SolveModel!(VortexModel ; get_gap = true)
 
 # Plot_FS!(H, bz, collect(0.05:0.05:0.5), [7]; cbar=true)
 # Plot_Band_Structure!(VortexModel , [bz.HighSymPoints["G"], bz.HighSymPoints["K1"], bz.HighSymPoints["M2"]] ; labels = [L"\Gamma", L"K_1", L"M_2"] )
 
 ###################### Real Space ######################################
 
-lattice     =   Lattice(UC, [6, 6])
-FillLattice!(lattice)
+# lattice     =   Lattice(UC, [6, 6])
+# FillLattice!(lattice)
 
-H = LatticeHamiltonian(lattice)
-DiagonalizeHamiltonian!(H)
+# # Plot_Lattice!(lattice ; bond_cmp=:viridis, site_size=8.0)
+
+# H = LatticeHamiltonian(lattice)
+# DiagonalizeHamiltonian!(H)
 
-plot(H.bands, marker = :circle, label = "Energies")
+# plot(H.bands, marker = :circle, label = "Energies")
 
-M   =   LatticeModel(lattice, H)
+# M   =   LatticeModel(lattice, H)
 
-Ta1 = Translations(lattice ; primitive  = 1, with_local = true)
-Ta2 = Translations(lattice ; primitive  = 2, with_local = true)
+# Ta1 = Translations(lattice ; primitive  = 1, with_local = true)
+# Ta2 = Translations(lattice ; primitive  = 2, with_local = true)
 
-m1 = GetPhase.(FindQuantumNumbers(M, Ta1 ; till_band = length(H.bands)÷2))
-m2 = GetPhase.(FindQuantumNumbers(M, Ta2 ; till_band = length(H.bands)÷2))
+# m1 = GetPhase.(FindQuantumNumbers(M, Ta1 ; till_band = length(H.bands)÷2))
+# m2 = GetPhase.(FindQuantumNumbers(M, Ta2 ; till_band = length(H.bands)÷2))
 
-mGround = mod.([sum(sum.(m1)), sum(sum.(m2))], Ref(1.0))
\ No newline at end of file
+# mGround = mod.([sum(sum.(m1)), sum(sum.(m2))], Ref(1.0))
\ No newline at end of file

From 1a8cc16bc666ce83dfb6f9667fe50e9d99701986 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Wed, 27 Sep 2023 18:23:47 -0400
Subject: [PATCH 03/14] function to add monopole strings

---
 src/Flux.jl | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/Flux.jl b/src/Flux.jl
index df165b6..86d638b 100644
--- a/src/Flux.jl
+++ b/src/Flux.jl
@@ -4,6 +4,7 @@ module Flux
 
     using LinearAlgebra
 
+    using ..TightBindingToolkit.Useful: SegmentIntersection
     using ..TightBindingToolkit.LatticeStruct: Lattice, FillLattice!
 
 
@@ -41,6 +42,30 @@ module Flux
     end
 
 
+    function GetStringPhase(Monopoles::Vector{Vector{Float64}}, r1::Vector{Float64}, r2::Vector{Float64} ;  flux::Float64 = pi) :: ComplexF64
+
+        t = SegmentIntersection(Monopoles, [r1, r2])
+
+        if t>=0.0 && t<1.0
+            return exp(im * flux)
+        else
+            return 1.0
+        end
+    end
+
+
+    function InsertMonopolePair(lat::Lattice{T}, Monopoles::Vector{Vector{Float64}} ; flux::Float64 = pi )
+
+        positions   =   getindex.(getindex.(Ref(lat.positions), collect(1:lat.length)) , Ref(1)) 
+        r1s         =   repeat(positions, 1, size(lat.BondSites, 2))
+        r2s         =   getindex.(getindex.(Ref(lat.positions), lat.BondSites) , Ref(1))
+
+        phases      =   GetStringPhase.(Ref(Monopoles), r1s, r2s ; flux = flux)
+        lat.BondMats=   lat.BondMats .* phases
+
+    end
+
+
 
 
 

From 3ca34dc5f5a2074c53a1ea7c7a8eac8ee90a7297 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 2 Oct 2023 16:39:24 -0400
Subject: [PATCH 04/14] Added an initial guess kwarg to binary search

---
 src/Useful.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/Useful.jl b/src/Useful.jl
index e3fbab3..ed00623 100644
--- a/src/Useful.jl
+++ b/src/Useful.jl
@@ -57,21 +57,19 @@ BinarySearch(target::Float64, xRange::Tuple{Float64, Float64}, f::T, args::Tuple
 General function implementing Binary search on a monotonic function f(x, args...)=target, in the range x ∈ xRange, with tolerance tol. 
 
 """
-    function BinarySearch(target::Float64, xRange::Tuple{Float64, Float64}, f::T, args::Tuple ; x_tol::Float64 = 0.001, target_tol::Float64 = 1e-6) where T<:Function
+    function BinarySearch(target::Float64, xRange::Tuple{Float64, Float64}, f::T, args::Tuple ; initial::Float64 = mean(xRange), x_tol::Float64 = 1e-3, target_tol::Float64 = 1e-6)::Float64 where T<:Function
+        ##### TODO pass an initial guess 
         xExt = collect(xRange)
-        current = nothing
+        @assert xExt[1] <= initial <= xExt[2] "xRange must be an increasing tuple of Float64, and initial guess must be in between."
+        current = initial
         ##### ///TODO : Fix bug when steps = 0
         if xExt[end]!= xExt[begin]
-            steps       =   Int(ceil(log2((xExt[end]-xExt[begin])/(x_tol)))) 
-        end
-
-        if steps == 0
-            steps   =   1
+            steps       =   max(Int(ceil(log2((xExt[end]-xExt[begin])/(x_tol)))), 1)
         end
         
         for i in 1:steps
-            current = mean(xExt)
             check = f(current, args...)
+
             if check - target > target_tol
                 xExt[end] = current
             elseif check - target < -target_tol
@@ -79,6 +77,8 @@ General function implementing Binary search on a monotonic function f(x, args...
             else
                 break
             end
+
+            current     =   mean(xExt)
         end
     
         return current

From 4286a4ddf584591a527de9ba68b5be1d128635ea Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 2 Oct 2023 16:40:11 -0400
Subject: [PATCH 05/14] SolveModel now has an initial guess for mu as a kwarg

---
 src/BdGModel.jl | 10 +++++-----
 src/Model.jl    | 14 +++++++-------
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/BdGModel.jl b/src/BdGModel.jl
index d780a8a..cdb46bc 100644
--- a/src/BdGModel.jl
+++ b/src/BdGModel.jl
@@ -7,7 +7,7 @@ module BdG
     using ..TightBindingToolkit.BZone:BZ, MomentumPhaseFFT
     using ..TightBindingToolkit.Hams:Hamiltonian, ModifyHamiltonianField!
     
-    using LinearAlgebra, Tullio, TensorCast, Logging
+    using LinearAlgebra, Tullio, TensorCast, Logging, Statistics
 
     import ..TightBindingToolkit.TBModel:FindFilling, GetMu!, GetFilling!, GetGk!, GetGr!, SolveModel!, GetGap!, FreeEnergy
 
@@ -108,8 +108,8 @@ GetMu!(M::BdGModel ;  tol::Float64=0.001)
 ```
 Function to get the correct chemical potential given a filling.
 """
-    function GetMu!(M::BdGModel ;  mu_tol::Float64 = 0.001, filling_tol::Float64 = 1e-6)
-        M.mu    =   BinarySearch(M.filling, M.Ham.bandwidth, FindFilling, (M,) ; x_tol = mu_tol, target_tol = filling_tol)
+    function GetMu!(M::BdGModel ;  guess::Float64 = 0.0, mu_tol::Float64 = 0.001, filling_tol::Float64 = 1e-6)
+        M.mu    =   BinarySearch(M.filling, M.Ham.bandwidth, FindFilling, (M,) ; initial = guess, x_tol = mu_tol, target_tol = filling_tol)
         @info "Found chemical potential μ = $(M.mu) for given filling = $(M.filling)."
     
     end
@@ -169,13 +169,13 @@ SolveModel!(M::BdGModel)
 ```
 one-step function to find all the attributes in  BdGModel after it has been initialized.
 """
-    function SolveModel!(M::BdGModel ; get_correlations::Bool = true, get_gap::Bool = false, verbose::Bool = true, mu_tol::Float64 = 1e-3, filling_tol::Float64 = 1e-6)
+    function SolveModel!(M::BdGModel ; mu_guess::Float64 = M.Ham.bandwidth[1] + M.filling * (M.Ham.bandwidth[2] - M.Ham.bandwidth[1]), get_correlations::Bool = true, get_gap::Bool = false, verbose::Bool = true, mu_tol::Float64 = 1e-3, filling_tol::Float64 = 1e-6)
         @assert M.Ham.is_BdG==true "Use other format for pure hopping Hamiltonian"
 
         if M.filling<0    ##### Must imply that filling was not provided by user and hence needs to be calculated from given mu
             GetFilling!(M)
         else
-            GetMu!(M ; mu_tol = mu_tol, filling_tol = filling_tol)
+            GetMu!(M ; guess = mu_guess, mu_tol = mu_tol, filling_tol = filling_tol)
         end
 
         if get_gap
diff --git a/src/Model.jl b/src/Model.jl
index 0236c6c..7e009e6 100644
--- a/src/Model.jl
+++ b/src/Model.jl
@@ -8,7 +8,7 @@ module TBModel
     using ..TightBindingToolkit.Parameters: Param
     using ..TightBindingToolkit.Hams:Hamiltonian
 
-    using LinearAlgebra, Logging
+    using LinearAlgebra, Logging, Statistics
 
 
 @doc """
@@ -77,13 +77,13 @@ GetMu!(M::Model; tol::Float64=0.001)
 Function to get chemical potential for a given `Model`, within a tolerance.
 
 """
-    function GetMu!(M::Model ;  mu_tol::Float64 = 0.001, filling_tol::Float64 = 1e-6)
-    
+    function GetMu!(M::Model ;  guess::Float64 = 0.0, mu_tol::Float64 = 0.001, filling_tol::Float64 = 1e-6)
+        ##### TODO Get an initial guess and pass to BinarySearch
         if M.T≈0.0
             energies  =   sort(reduce(vcat, M.Ham.bands))
             M.mu      =   (energies[floor(Int64, length(energies)*M.filling)] + energies[floor(Int64, length(energies)*M.filling)+1])/2
         else
-            M.mu      =   BinarySearch(M.filling, M.Ham.bandwidth, FindFilling, (M.Ham.bands, M.T, M.stat) ; x_tol=mu_tol, target_tol = filling_tol)
+            M.mu      =   BinarySearch(M.filling, M.Ham.bandwidth, FindFilling, (M.Ham.bands, M.T, M.stat) ; initial = guess, x_tol=mu_tol, target_tol = filling_tol)
             @info "Found chemical potential μ = $(M.mu) for given filling = $(M.filling)."
         end
     end
@@ -153,13 +153,13 @@ SolveModel!(M::Model)
 one-step function to find all the attributes in Model after it has been initialized.
 
 """
-    function SolveModel!(M::Model ; get_correlations::Bool = true, get_gap::Bool = false, verbose::Bool = true, mu_tol::Float64 = 1e-3, filling_tol::Float64 = 1e-6)
+    function SolveModel!(M::Model ; mu_guess::Float64 = M.Ham.bandwidth[1] + M.filling * (M.Ham.bandwidth[2] - M.Ham.bandwidth[1]), get_correlations::Bool = true, get_gap::Bool = false, verbose::Bool = true, mu_tol::Float64 = 1e-3, filling_tol::Float64 = 1e-6)
         @assert M.Ham.is_BdG==false "BdG Hamiltonians should be solved using a BdGModel"
-
+        ##### TODO Pass initial guess of chemical potential to GetMu! if provided by user
         if M.filling<0    ##### Must imply that filling was not provided by user and hence needs to be calculated from given mu
             GetFilling!(M)
         else
-            GetMu!(M ; mu_tol = mu_tol, filling_tol = filling_tol)
+            GetMu!(M ; guess = mu_guess, mu_tol = mu_tol, filling_tol = filling_tol)
         end
 
         if get_gap

From dc81ca9c45e7ac211773cfcdce7c753ed088445e Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 2 Oct 2023 16:40:52 -0400
Subject: [PATCH 06/14] Monopole string now takes into account bonds which go
 across the lattice

---
 src/Flux.jl | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/Flux.jl b/src/Flux.jl
index 86d638b..ea0d642 100644
--- a/src/Flux.jl
+++ b/src/Flux.jl
@@ -1,6 +1,6 @@
 module Flux
 
-    export GetBondPhase, CheckGaugeValidity, ModifyGauge!
+    export GetBondPhase, CheckGaugeValidity, ModifyGauge!, GetStringPhase, InsertMonopolePair!
 
     using LinearAlgebra
 
@@ -42,11 +42,12 @@ module Flux
     end
 
 
-    function GetStringPhase(Monopoles::Vector{Vector{Float64}}, r1::Vector{Float64}, r2::Vector{Float64} ;  flux::Float64 = pi) :: ComplexF64
+    function GetStringPhase(Monopoles::Vector{Vector{Float64}}, r1::Vector{Float64}, r2::Vector{Float64} ;  flux::Float64 = Float64(pi)) :: ComplexF64
 
         t = SegmentIntersection(Monopoles, [r1, r2])
+        r = SegmentIntersection([r1, r2], Monopoles)
 
-        if t>=0.0 && t<1.0
+        if t>=0.0 && t<1.0 && r>=0.0 && r<1.0
             return exp(im * flux)
         else
             return 1.0
@@ -54,15 +55,21 @@ module Flux
     end
 
 
-    function InsertMonopolePair(lat::Lattice{T}, Monopoles::Vector{Vector{Float64}} ; flux::Float64 = pi )
+    function InsertMonopolePair!(lat::Lattice{T}, Monopoles::Vector{Vector{Float64}} ; flux::Float64 = Float64(pi) ) where{T}
 
         positions   =   getindex.(getindex.(Ref(lat.positions), collect(1:lat.length)) , Ref(1)) 
         r1s         =   repeat(positions, 1, size(lat.BondSites, 2))
         r2s         =   getindex.(getindex.(Ref(lat.positions), lat.BondSites) , Ref(1))
 
+        dists       =   norm.(r2s .- r1s)
+        check       =   isapprox.(dists, lat.BondDists, atol = 1e-6, rtol = 1e-6)
+
         phases      =   GetStringPhase.(Ref(Monopoles), r1s, r2s ; flux = flux)
+        phases      =   phases .^ check
         lat.BondMats=   lat.BondMats .* phases
 
+        return findall(==(false), isapprox.(phases, Ref(1.0)))
+
     end
 
 

From ff29b3dce29f136eb9acfcf642ff5b1274549af8 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 2 Oct 2023 16:41:16 -0400
Subject: [PATCH 07/14] added function for gauge fixing wavefunctions : needs
 to be tested

---
 src/LatticeHamiltonian.jl | 34 ++++++++++++++++++++++++++++------
 1 file changed, 28 insertions(+), 6 deletions(-)

diff --git a/src/LatticeHamiltonian.jl b/src/LatticeHamiltonian.jl
index f944f99..39dfba0 100644
--- a/src/LatticeHamiltonian.jl
+++ b/src/LatticeHamiltonian.jl
@@ -1,5 +1,5 @@
 module LatHam
-    export FillHamiltonian, LatticeHamiltonian, DiagonalizeHamiltonian!, Slater, SingleParticleFidelity, SlaterOverlap
+    export FillHamiltonian, LatticeHamiltonian, DiagonalizeHamiltonian!, Slater, SingleParticleFidelity, SlaterOverlap, GaugeFix!
 
     using ..TightBindingToolkit.LatticeStruct: Lattice
     import ..TightBindingToolkit.Hams: FillHamiltonian, DiagonalizeHamiltonian!
@@ -62,11 +62,33 @@ module LatHam
     end
 
 
-    function SingleParticleFidelity(H1::LatticeHamiltonian, H2::LatticeHamiltonian) :: Matrix{ComplexF64}
+    function GaugeFix!(H::LatticeHamiltonian, gauge::Tuple{Int64, Float64} = (1, 0.0))
 
-        @assert size(H1) == size(H2) "The two hamiltonians must be the same size!"
+        site, DesiredPhase = gauge[1], gauge[2]
 
-        fidelity = adjoint(H1.states) * H2.states
+        NormCheck = (abs2.(H.states[site:end, :]) .> 0.0)
+        NonTrivialNormSites = findfirst.(==(1), eachcol(NormCheck)) .+ (site - 1)     
+
+        phases = angle.(getindex.(Ref(H.states), NonTrivialNormSites, 1:length(H.bands)))
+        phaseShift = exp.(im .* (DesiredPhase .- phases))
+
+        H.states = H.states .* phaseShift'
+
+        return (NonTrivialNormSites, phaseShift)
+
+    end
+
+
+    function SingleParticleFidelity(H1::LatticeHamiltonian, H2::LatticeHamiltonian, states::Vector{Int64} = collect(1:length(H1.bands)) ; fixGauge::Bool = true, gauge::Tuple{Int64, Float64} = (1, 0.0)) :: Matrix{ComplexF64}
+
+        @assert size(H1.H) == size(H2.H) "The two hamiltonians must be the same size!"
+
+        if fixGauge
+            GaugeFix!(H1, gauge)
+            GaugeFix!(H2, gauge)
+        end
+
+        fidelity = adjoint(H1.states[:, states]) * H2.states[:, states]
 
         return fidelity
     end
@@ -83,9 +105,9 @@ module LatHam
     end
 
 
-    function SlaterOverlap(H1::LatticeHamiltonian, H2::LatticeHamiltonian) :: ComplexF64
+    function SlaterOverlap(H1::LatticeHamiltonian, H2::LatticeHamiltonian, states::Vector{Int64} = collect(1:Int64(length(H1.bands)//2)); fixGauge::Bool = true, gauge::Tuple{Int64, Float64} = (1, 0.0)) :: ComplexF64
 
-        return det(SingleParticleFidelity(H1, H2))
+        return det((SingleParticleFidelity(H1, H2, states ; fixGauge = fixGauge, gauge = gauge)))
     end
 
 end
\ No newline at end of file

From 0ba62b7b2ad19a33bcfb537f63be640715fb2697 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 2 Oct 2023 16:41:26 -0400
Subject: [PATCH 08/14] updated module code

---
 src/TightBindingToolkit.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/TightBindingToolkit.jl b/src/TightBindingToolkit.jl
index b8bf84b..8e95c6d 100644
--- a/src/TightBindingToolkit.jl
+++ b/src/TightBindingToolkit.jl
@@ -40,7 +40,7 @@ export CreateLattice, ModifyLattice!, ScaleLatticeBonds!, ScaleLattice!, RemoveL
 
 include("Flux.jl")
 using .Flux
-export GetBondPhase, CheckGaugeValidity, ModifyGauge!
+export GetBondPhase, CheckGaugeValidity, ModifyGauge!, GetStringPhase, InsertMonopolePair!
 
 ##### Module to define Brillouin Zone structure and some functions to manipulate it.
 include("BZ.jl")
@@ -54,7 +54,7 @@ export Hamiltonian , FillHoppingHamiltonian, FillPairingHamiltonian, FillHamilto
 
 include("LatticeHamiltonian.jl")
 using .LatHam
-export FillHamiltonian, LatticeHamiltonian, DiagonalizeHamiltonian!, Slater, SingleParticleFidelity, SlaterOverlap
+export FillHamiltonian, LatticeHamiltonian, DiagonalizeHamiltonian!, Slater, SingleParticleFidelity, SlaterOverlap, GaugeFix!
 
 ##### Module to define a Tight-Binding Model structure which takes into account thermodynamical parameters such as temperature and filling etc.
 include("Model.jl")

From ded353520cec19e75cb0ec69ad23b34a57aa2814 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <39477277+Anjishnubose@users.noreply.github.com>
Date: Mon, 2 Oct 2023 18:23:21 -0400
Subject: [PATCH 09/14] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index bee4873..6d95667 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TightBindingToolkit"
 uuid = "2233325a-6eb3-486f-aff0-670e0939fa7e"
 authors = ["Anjishnu Bose", "Sreekar Voleti", "Samuel Vadnais"]
-version = "2.3.0"
+version = "2.3.1"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

From c3723ae58dffec72f75cc5b649904543e878ceba Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Tue, 3 Oct 2023 13:41:45 -0400
Subject: [PATCH 10/14] can use UC of different rank than param when adding
 bonds to param

---
 src/Params.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/Params.jl b/src/Params.jl
index 076a584..ea4b8e9 100644
--- a/src/Params.jl
+++ b/src/Params.jl
@@ -50,7 +50,7 @@ Add a bond with the given attributes to `param`.
 If given `mat` attribute is a number, it is converted into a 1x1 matrix when entered into the bond.
 
 """
-    function AddAnisotropicBond!( param::Param{T, R}, uc::UnitCell{T} , base::Int64 , target::Int64 , offset::Vector{Int64} , mat::Array{<:Number, T} , dist::Float64, label::String ) where {T, R}
+    function AddAnisotropicBond!( param::Param{T, R}, uc::UnitCell{T2} , base::Int64 , target::Int64 , offset::Vector{Int64} , mat::Array{<:Number, T} , dist::Float64, label::String ) where {T, R, T2}
 
         dims 	=	repeat([uc.localDim], T)
         @assert size(mat) == Tuple(dims) "Interaction matrix has the wrong dimension!"
@@ -73,7 +73,7 @@ If given `mat` attribute is a number, it is converted into a 1x1 matrix when ent
         end
     end
     
-    function AddAnisotropicBond!( param::Param{T, R}, uc::UnitCell{T} , base::Int64 , target::Int64 , offset::Vector{Int64} , mat::Number , dist::Float64, label::String ) where {T, R}
+    function AddAnisotropicBond!( param::Param{T, R}, uc::UnitCell{T2} , base::Int64 , target::Int64 , offset::Vector{Int64} , mat::Number , dist::Float64, label::String ) where {T, R, T2}
 	
         @assert uc.localDim == 1
         dims 	=	repeat([uc.localDim], T)
@@ -93,7 +93,7 @@ The input `checkOffsetRange` must be adjusted depending on the input distance.
 The optional input `subs` is meant for isotropic bonds when only a subset of sublattices are involved.
 
 """    
-    function AddIsotropicBonds!( param::Param{T, R}, uc::UnitCell{T} , dist::Float64 , mat::Array{<:Number, T} , label::String; checkOffsetRange::Int64=2 , subs::Vector{Int64}=collect(1:length(uc.basis)) ) where {T, R}
+    function AddIsotropicBonds!( param::Param{T, R}, uc::UnitCell{T2} , dist::Float64 , mat::Array{<:Number, T} , label::String; checkOffsetRange::Int64=2 , subs::Vector{Int64}=collect(1:length(uc.basis)) ) where {T, R, T2}
     
         dims 	=	repeat([uc.localDim], T)
         @assert size(mat) == Tuple(dims) "Interaction matrix has the wrong dimension!"
@@ -123,7 +123,7 @@ The optional input `subs` is meant for isotropic bonds when only a subset of sub
         end
     end
 
-    function AddIsotropicBonds!( param::Param{T, R}, uc::UnitCell{T} , dist::Float64 , mat::Number , label::String; checkOffsetRange::Int64=2 , subs::Vector{Int64}=collect(1:length(uc.basis))) where {T, R}
+    function AddIsotropicBonds!( param::Param{T, R}, uc::UnitCell{T2} , dist::Float64 , mat::Number , label::String; checkOffsetRange::Int64=2 , subs::Vector{Int64}=collect(1:length(uc.basis))) where {T, R, T2}
 
         @assert uc.localDim == 1
         dims 	=	repeat([uc.localDim], T)
@@ -134,13 +134,13 @@ The optional input `subs` is meant for isotropic bonds when only a subset of sub
 
 @doc """
 ```julia
-AddSimilarBond!(param::Param{T, R}, uc::UnitCell{T}, bond::Bond{T} ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R}
-AddSimilarBond!(param::Param{T, R}, uc::UnitCell{T}, base::Int64, target::Int64, offset::Vector{Int64}, mat::Array{ComplexF64, T}, dist::Float64, label::String ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R}
+AddSimilarBond!(param::Param{T, R}, uc::UnitCell{T2}, bond::Bond{T} ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R}
+AddSimilarBond!(param::Param{T, R}, uc::UnitCell{T2}, base::Int64, target::Int64, offset::Vector{Int64}, mat::Array{ComplexF64, T}, dist::Float64, label::String ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R}
 ```
 Function to add bonds which are not completely isotropic, but are still related by translation (not by the unit cell primitives but by the underlying lattice primitives).
 
 """
-    function AddSimilarBonds!(param::Param{T, R}, uc::UnitCell{T}, bond::Bond{T} ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R}
+    function AddSimilarBonds!(param::Param{T, R}, uc::UnitCell{T2}, bond::Bond{T} ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R, T2}
 
         MotherParam     =   Param(1.0, 2)
         AddIsotropicBonds!(MotherParam, uc, bond.dist, bond.mat, bond.label ; checkOffsetRange = checkOffsetRange , subs = subs)
@@ -164,7 +164,7 @@ Function to add bonds which are not completely isotropic, but are still related
 
     end
 
-    function AddSimilarBonds!(param::Param{T, R}, uc::UnitCell{T}, base::Int64, target::Int64, offset::Vector{Int64}, mat::Array{ComplexF64, T}, dist::Float64, label::String ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R}
+    function AddSimilarBonds!(param::Param{T, R}, uc::UnitCell{T2}, base::Int64, target::Int64, offset::Vector{Int64}, mat::Array{ComplexF64, T}, dist::Float64, label::String ;  subs::Vector{Int64}=collect(1:length(uc.basis)), checkOffsetRange::Int64=2) where {T, R, T2}
 
         AddSimilarBonds!(param, uc, Bond(base, target, offset, mat, dist, label) ; subs = subs, checkOffsetRange = checkOffsetRange)
     end

From bd862ab528f1e2bd3ecd65a96aab76305bad2355 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <39477277+Anjishnubose@users.noreply.github.com>
Date: Tue, 3 Oct 2023 13:43:23 -0400
Subject: [PATCH 11/14] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 6d95667..9a61b02 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "TightBindingToolkit"
 uuid = "2233325a-6eb3-486f-aff0-670e0939fa7e"
 authors = ["Anjishnu Bose", "Sreekar Voleti", "Samuel Vadnais"]
-version = "2.3.1"
+version = "2.3.2"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

From 43305ff453bf5397e90a19364e5198011320dcef Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 23 Oct 2023 13:21:57 -0400
Subject: [PATCH 12/14] added bijections

---
 Project.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 6d95667..b8c34c7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Anjishnu Bose", "Sreekar Voleti", "Samuel Vadnais"]
 version = "2.3.1"
 
 [deps]
+Bijections = "e2ed5e7c-b2de-5872-ae92-c73ca462fb04"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -24,9 +25,9 @@ MappedArrays = "^0.4"
 Plots = "^1.38"
 PyPlot = "2.11"
 Statistics = "1.0 - 2.0"
+StatsBase = "^0.34"
 TensorCast = "^0.4"
 Tullio = "^0.3"
-StatsBase = "^0.34"
 julia = "^1.7, ^1.9"
 
 [extras]

From de51b0dd2f63ec4dbeade514e4acfc49ae1b0101 Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 23 Oct 2023 19:25:48 -0400
Subject: [PATCH 13/14] added a secant method for root finding

---
 src/Useful.jl | 77 ++++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 73 insertions(+), 4 deletions(-)

diff --git a/src/Useful.jl b/src/Useful.jl
index ed00623..c34e53d 100644
--- a/src/Useful.jl
+++ b/src/Useful.jl
@@ -1,5 +1,5 @@
 module Useful
-    export GetAllOffsets, VecAngle, Meshgrid, BinarySearch, DistFunction, DeriDistFunction , GetIndexPath, FFTArrayofMatrix, Central_Diff, Arrayfy, DeArrayfy, GetPhase, SegmentIntersection
+    export GetAllOffsets, VecAngle, Meshgrid, BinarySearch, DistFunction, DeriDistFunction , GetIndexPath, FFTArrayofMatrix, Central_Diff, Arrayfy, DeArrayfy, GetPhase, SegmentIntersection, SwitchKroneckerBasis, SecantSearch, ImpIndices
 
     using LinearAlgebra, Statistics, FFTW, TensorCast
 
@@ -58,10 +58,17 @@ General function implementing Binary search on a monotonic function f(x, args...
 
 """
     function BinarySearch(target::Float64, xRange::Tuple{Float64, Float64}, f::T, args::Tuple ; initial::Float64 = mean(xRange), x_tol::Float64 = 1e-3, target_tol::Float64 = 1e-6)::Float64 where T<:Function
-        ##### TODO pass an initial guess 
+        ##### ///TODO pass an initial guess 
         xExt = collect(xRange)
-        @assert xExt[1] <= initial <= xExt[2] "xRange must be an increasing tuple of Float64, and initial guess must be in between."
-        current = initial
+        ##### Choosing initial point to start Binary search with. 
+        if initial < xExt[begin]
+            current = xExt[begin]
+        elseif xExt[end] < initial
+            current = xExt[end]
+        else
+            current = initial
+        end 
+
         ##### ///TODO : Fix bug when steps = 0
         if xExt[end]!= xExt[begin]
             steps       =   max(Int(ceil(log2((xExt[end]-xExt[begin])/(x_tol)))), 1)
@@ -85,6 +92,27 @@ General function implementing Binary search on a monotonic function f(x, args...
     end
 
 
+    function SecantSearch(target::Float64, xInitials::Vector{Float64}, yInitials::Vector{Float64}, f::T, args::Tuple ; x_tol::Float64 = 1e-3, target_tol::Float64 = 1e-4, max_iter::Int64 = 10)::Float64 where T<:Function
+
+        x0, x1  =   xInitials
+        y0, y1  =   yInitials
+        
+        for i in 1:max_iter
+            x  =   x1 - (x1 - x0) * (y1 - target) / (y1 - y0)
+            check   =   f(x, args...)
+
+            if abs(check - target) < target_tol || abs(x - x1) < x_tol
+                break
+            end
+
+            x0, x1  =   x1, x
+            y0, y1  =   y1, check
+        end
+    
+        return x
+    end
+
+
 @doc """
 ```julia
 DistFunction(E ; T::Float64, mu::Float64=0.0, stat::Int64=-1)
@@ -255,6 +283,13 @@ PBC   : vector of boolean with length = spatial dimensions = rank of f ---> if t
     end
 
 
+@doc """
+```julia
+SegmentIntersection(SegmentPoints::Vector{Vector{Float64}}, LinePoints::Vector{Vector{Float64}})-->Float64
+```
+Given two points on a line and two points on a segment, returns the parameter t such that the intersection point is given by (1-t) * SegmentPoints[1] + t * SegmentPoints[2]. 
+If the intersection point is not on the segment, returns a value outside the range [0, 1].
+"""
     function SegmentIntersection(SegmentPoints::Vector{Vector{Float64}}, LinePoints::Vector{Vector{Float64}})::Float64
 
         M1, M2  =   SegmentPoints
@@ -270,5 +305,39 @@ PBC   : vector of boolean with length = spatial dimensions = rank of f ---> if t
         return t
     end
 
+    function Linearize(Ns::Tuple{Int64, Int64}, index::Tuple{Int64})::Int64
+
+        return Ns[2] * (index[1] - 1) + index[2]
+    end
+
+    function SwitchKroneckerBasis(Ns::Tuple{Int64, Int64})::Vector{Int64}
+
+        indices     =   Meshgrid(collect(Ns))
+        indices     =   vcat(indices...)
+        order       =   Linearize.(Ref(Ns), indices)
+
+        return order
+    end
+
+
+    function ImpIndices(A::Matrix{Float64}, sizeB::Tuple{Int64, Int64})::Vector{Tuple{Int64, Int64}}
+        inds = Tuple{Int64, Int64}[]
+
+        for (i, j) in Meshgrid(collect(sizeB))
+
+            B = zeros(Float64, sizeB...)
+            B[i, j] = 1.0
+
+            C = B * A
+
+            if sum(.!(isapprox.(C, zeros(Float64, size(C))))) > 0   ##### if the matrix product is not zero => this element of B can lead to non-zero result
+                push!(inds, (i, j))
+            end
+        end
+
+        return inds
+    end
+
+
 
 end
\ No newline at end of file

From d1cf98be46fbaa6bab4178c241045f8c3289758e Mon Sep 17 00:00:00 2001
From: Anjishnu Bose <anjishnubose98@gmail.com>
Date: Mon, 23 Oct 2023 19:25:59 -0400
Subject: [PATCH 14/14] added RPA calculation

---
 src/Susceptibility.jl | 174 +++++++++++++++++++++++++++++++-----------
 1 file changed, 129 insertions(+), 45 deletions(-)

diff --git a/src/Susceptibility.jl b/src/Susceptibility.jl
index ce73d47..7a98b1d 100644
--- a/src/Susceptibility.jl
+++ b/src/Susceptibility.jl
@@ -2,9 +2,11 @@ module suscep
     export Susceptibility , FindChi , FillChis!
 
     using ..TightBindingToolkit.SpinMatrices:SpinMats
-    using ..TightBindingToolkit.Useful: DistFunction
+    using ..TightBindingToolkit.Useful: DistFunction, SwitchKroneckerBasis, ImpIndices
     using ..TightBindingToolkit.BZone:BZ, GetQIndex
     using ..TightBindingToolkit.TBModel:Model
+    using ..TightBindingToolkit.Parameters:Param
+    using ..TightBindingToolkit.UCell:UnitCell
 
 
     using LinearAlgebra, TensorCast, Tullio, Logging
@@ -12,6 +14,28 @@ module suscep
     directions = ["x" , "y" , "z"]
 
 
+    function ResolvedQInteraction(Q::Vector{Float64}, interaction::Param{T, Float64}, uc::UnitCell{T2})::Matrix{Array{ComplexF64, T}} where {T, T2}
+
+        resolvedInt     =   repeat([zeros(ComplexF64, size(interaction.unitBond[begin].bondMat))], length(uc.basis), length(uc.basis))     ##### Initializing the sublattice resolved interaction matrix
+        flippedIndices 	=	collect(T:-1:1)
+
+
+        for bond in interaction.unitBonds
+
+            contribution    =   bond.mat * exp.(im .* dot.(Ref(Q) , bond.offset .* uc.primitives))     ##### Adding the bond contribution to the resolved interaction matrix
+            contribution    +=  collect(conj.(permutedims(bond.mat, flippedIndices))) * exp.(-im .* dot.(Ref(Q) , bond.offset .* uc.primitives))    ##### Adding the hermitian conjugate bond contribution to the resolved interaction matrix
+
+            resolvedInt[bond.base, bond.target]     +=  contribution / 2
+            resolvedInt[bond.target, bond.base]     +=  contribution / 2
+        end
+
+        return resolvedInt
+    end
+
+
+
+
+
 @doc """
 `Susceptibility` is a data type representing the magnetic response, ``χ^{ab}(Q , Ω)`` for a general tight-binding `Model`.
 
@@ -29,14 +53,42 @@ Susceptibility(M::Model , Qs::Vector{Vector{Float64}}, Omegas::Vector{Float64} ;
 ```
 """
     mutable struct Susceptibility
-        M       ::  Model   
+        
         Qs      ::  Vector{Vector{Float64}}
         Omegas  ::  Vector{Float64}
         Spread  ::  Float64
-        chis    ::  Dict
+        Bare        ::  Matrix{Matrix{ComplexF64}}  ##### The bare susceptibility
+        BareReduced ::  Matrix{Matrix{ComplexF64}}  ##### The reduced bare susceptibility
+        RPA         ::  Matrix{Matrix{ComplexF64}}  ##### The RPA susceptibility
+        RPAReduced  ::  Matrix{Matrix{ComplexF64}}  ##### The reduced RPA susceptibility
+
+        function Susceptibility(Omegas::Vector{Float64} , M::Model;  fill_bz::Bool = true , eta::Float64 = 1e-2)
+
+            if fill_bz
+                Qs      =   reshape(M.bz.ks, length(M.bz.ks))
+            else
+                Qs      =   [zeros(Float64, length(M.uc.primitives))]
+            end
+
+            Bare    =   repeat([zeros(ComplexF64, length(M.uc.basis) * length(M.uc.OnSiteMats), length(M.uc.basis) * length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
+            BareReduced     =   repeat([zeros(ComplexF64, length(M.uc.OnSiteMats), length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
+
+            RPA    =   repeat([zeros(ComplexF64, length(M.uc.basis) * length(M.uc.OnSiteMats), length(M.uc.basis) * length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
+            RPAReduced  =   repeat([zeros(ComplexF64, length(M.uc.OnSiteMats), length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
 
-        Susceptibility(M::Model , Omegas::Vector{Float64} ;  eta::Float64 = 1e-2) = new{}(M, [], Omegas, eta, Dict())
-        Susceptibility(M::Model , Qs::Vector{Vector{Float64}}, Omegas::Vector{Float64} ;  eta::Float64 = 1e-2) = new{}(M, Qs, Omegas, eta, Dict())
+            return new{}(Qs, Omegas, eta, Bare, BareReduced, RPA, RPAReduced)
+        end
+
+        function Susceptibility(Qs::Vector{Vector{Float64}}, Omegas::Vector{Float64},  M::Model;  eta::Float64 = 1e-2) 
+            
+            Bare    =   repeat([zeros(ComplexF64, length(M.uc.basis) * length(M.uc.OnSiteMats), length(M.uc.basis) * length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
+            BareReduced     =   repeat([zeros(ComplexF64, length(M.uc.OnSiteMats), length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
+
+            RPA    =   repeat([zeros(ComplexF64, length(M.uc.basis) * length(M.uc.OnSiteMats), length(M.uc.basis) * length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
+            RPAReduced  =   repeat([zeros(ComplexF64, length(M.uc.OnSiteMats), length(M.uc.OnSiteMats))], length(Omegas), length(Qs))
+
+            return new{}(Qs, Omegas, eta, Bare, BareReduced, RPA, RPAReduced)
+        end
     end
 
 
@@ -47,78 +99,96 @@ nF_factor(Q_index::Vector{Int64}, Omega::Float64 , M::Model; eta::Float64=1e-2)
 function to calculate ``(nF(E(k)) - nF(E(k+Q))) / (Ω + i * η - (E(k+Q) - E(k)))`` for a general multi-band system, at all k-points.
 
 """
-    function nF_factor(Q_index::Vector{Int64}, Omega::Float64 , M::Model; eta::Float64=1e-2) :: Matrix{Matrix{ComplexF64}}
+    function nF_factor(Q_index::Vector{Int64}, Omega::Float64 , M::Model; eta::Float64=1e-2) :: Vector{Matrix{ComplexF64}}
         
         Es_k    =   M.Ham.bands
         Es_kpQ  =   circshift(Es_k, -Q_index)  ##### The bands with momentum indices rolled to k+Q instead of k.
+        dEs     =   reshape(map((x, y) -> x .- y, Es_kpQ, Es_k), length(M.Ham.bands))  ##### The energy differences between the bands at k+Q and k. Flattened momentum index.
 
         ##### The fermi distribution functions for each band at each k-point
         nf_k    =   DistFunction.(Es_k ; T=M.T, mu=M.mu, stat=M.stat) 
         nf_kpQ  =   circshift(nf_k, -Q_index)
+        dnFs    =   reshape(map((x, y) -> x .- y, nf_k, nf_kpQ), length(M.Ham.bands))  ##### The fermi distribution function differences between the bands at k+Q and k. Flatened momentum index.
 
-        @cast nF[k1, k2][i, j] |= (nf_k[k1, k2][i] - nf_kpQ[k1, k2][j]) / (Omega + im * eta - (Es_kpQ[k1, k2][j] - Es_k[k1, k2][i]));
-        ##### Returning a matrix of matrix type
+        @cast nF[k][i, j] |= dnFs[k][i, j] / (Omega + im * eta - dEs[k][j, i]);
         return nF
     end 
 
 
 @doc """
 ```julia
-FindChi(Q::Vector{Float64}, Omega::Float64 , M::Model; a::Int64=3, b::Int64=3, eta::Float64=1e-2) --> ComplexF64
+FindReducedChi(Q::Vector{Float64}, Omega::Float64 , M::Model; a::Int64=3, b::Int64=3, eta::Float64=1e-2) --> ComplexF64
 ```
 function to calculate susceptibility at a fixed Ω=`Omega`, and `Q`, and along a fixed direction given by `a` and `b`.
 
 """
-    function FindChi(Q::Vector{Float64}, Omega::Float64 , M::Model; a::Int64=3, b::Int64=3, eta::Float64=1e-2, SpinMatrix::Vector{Matrix{ComplexF64}} = SpinMats((M.uc.localDim-1)//2)) ::ComplexF64
+    function FindReducedChi(Omega::Float64, Q::Vector{Float64}, M::Model; a::Int64=3, b::Int64=3, eta::Float64=1e-2, Gamma_index::Vector{Int64} = ((M.bz.gridSize .+ 1) .÷ 2)) ::ComplexF64
 
         expQ    =   diagm(exp.(-im .* dot.(Ref(Q) , M.uc.basis)))
-        Qa      =   kron(expQ , SpinMatrix[a])
-        Qb      =   kron(expQ , SpinMatrix[b])
-        @assert size(Qa) == size(Qb) == (length(M.uc.basis) * M.uc.localDim, length(M.uc.basis) * M.uc.localDim)
+        Qa      =   kron(expQ , M.uc.OnSiteMats[a])
+        Qb      =   kron(expQ , M.uc.OnSiteMats[b])
 
-        Q_index =   GetQIndex(Q, M.bz) .- ((M.bz.gridSize .+ 1) .÷ 2)
+        Q_index =   GetQIndex(Q, M.bz) .- Gamma_index
         nF      =   nF_factor(Q_index, Omega , M; eta=eta)
 
         Uk      =   M.Ham.states
         UkQ     =   circshift(M.Ham.states, -Q_index)
 
         if a!=b
-            Ma      =   adjoint.(UkQ) .* Ref(Qa) .* Uk
-            Mb      =   adjoint.(UkQ) .* Ref(Qb) .* Uk
-            @tullio chi     :=   - Ma[k1, k2][i, j] * conj(Mb[k1, k2][i, j]) * nF[k1, k2][j, i]
+            Ma      =   reshape(adjoint.(UkQ) .* Ref(Qa) .* Uk, length(M.Ham.bands))
+            Mb      =   reshape(adjoint.(UkQ) .* Ref(Qb) .* Uk, length(M.Ham.bands))
+            @tullio chi     :=   - Ma[k][i, j] * conj(Mb[k][i, j]) * nF[k][j, i]
         else
-            Ma      =   adjoint.(UkQ) .* Ref(Qa) .* Uk
-            @tullio chi     :=   - Ma[k1, k2][i, j] * conj(Ma[k1, k2][i, j]) * nF[k1, k2][j, i]
+            Ma      =   reshape(adjoint.(UkQ) .* Ref(Qa) .* Uk, length(M.Ham.bands))
+            @tullio chi     :=   - Ma[k][i, j] * conj(Ma[k][i, j]) * nF[k][j, i]
         end
 
         return chi 
     end
 
 
-@doc """
-```julia
-FindChi_FullBZ(Omega::Float64 , M::Model ; a::Int64=3, b::Int64=3, eta::Float64=1e-2) --> Matrix{ComplexF64}
-```
-function to calculate susceptibility at a fixed Ω=`Omega`, but for all `Q` present in the `BZ`, and along a fixed direction given by `a` and `b`.
+    function FindChi(Omega::Float64, Q::Vector{Float64}, M::Model; a::Int64=3, b::Int64=3, eta::Float64=1e-2, Gamma_index::Vector{Int64} = ((M.bz.gridSize .+ 1) .÷ 2)) ::Matrix{ComplexF64}
 
-"""
-    function FindChi_FullBZ(Omega::Float64 , M::Model ; a::Int64=3, b::Int64=3, eta::Float64=1e-2, SpinMatrix::Vector{Matrix{ComplexF64}} = SpinMats((M.uc.localDim-1)//2)) :: Matrix{ComplexF64}
-        return FindChi.(M.bz.ks, Ref(Omega), Ref(M) ; a=a, b=b, eta=eta, SpinMatrix = SpinMatrix)
-    end
+        Q_index =   GetQIndex(Q, M.bz) .- Gamma_index
+        nF      =   nF_factor(Q_index, Omega , M; eta=eta)
+
+        Uk      =   M.Ham.states
+        UkQ     =   circshift(M.Ham.states, -Q_index)
 
+        chiMat     =   zeros(ComplexF64, length(M.uc.basis), length(M.uc.basis))
 
-@doc """
-```julia
-FindChi_Path(Omega::Float64 , M::Model, path::Vector{Vector{Float64}} ; a::Int64=3, b::Int64=3, eta::Float64=1e-2) --> Vector{ComplexF64}
-```
-function to calculate susceptibility at a fixed Ω=`Omega`, but for all `Q` present in the given path, and along a fixed direction given by `a` and `b`.
+        for (i, base) in enumerate(M.uc.basis)
+            for (j, target) in enumerate(M.uc.basis)
 
-"""
-    function FindChi_Path(Omega::Float64 , M::Model, path::Vector{Vector{Float64}} ; a::Int64=3, b::Int64=3, eta::Float64=1e-2, SpinMatrix::Vector{Matrix{ComplexF64}} = SpinMats((M.uc.localDim-1)//2)) :: Vector{ComplexF64}
-        return FindChi.(path, Ref(Omega), Ref(M) ; a=a, b=b, eta=eta, SpinMatrix = SpinMatrix)
+                Uis     =   selectdim.(Uk, 1, (i - 1) * M.uc.localDim + 1 : i * M.uc.localDim)
+                Ujs     =   selectdim.(adjoint.(Uk), 2, (j - 1) * M.uc.localDim + 1 : j * M.uc.localDim)
+
+                UisQ    =   selectdim.(adjoint.(UkQ), 2, (i - 1) * M.uc.localDim + 1 : i * M.uc.localDim)
+                UjsQ    =   selectdim.(UkQ, 1, (j - 1) * M.uc.localDim + 1 : j * M.uc.localDim)
+
+                Ma      =   reshape(UisQ .* Ref(M.uc.OnSiteMats[a]) .* Uis, length(M.Ham.bands))
+                Mb      =   reshape(Ujs .* Ref(M.uc.OnSiteMats[b]) .* UjsQ, length(M.Ham.bands))
+
+                @tullio chi     :=   - Ma[k][i, j] * Mb[k][i, j] * nF[k][j, i]
+                chiMat[i, j]    +=  chi
+            end
+        end
+
+        return chiMat 
     end
 
 
+    function PerformRPA!(chi::Susceptibility, GeneratorInteraction::Param{2, Float64}, uc::UnitCell{T}) where {T}
+        ##### GeneratorInteraction refers to the interaction written in the basis of all the generators of the local Hilbert space (including the identity matrix). Eg. Id + SU(2) spins -> GeneratorInteraction is a 4x4 matrix.
+        resolvedInts    =   ResolvedQInteraction.(chi.Qs, Ref(interaction), Ref(chi.M.uc))
+        resolvedInts    =   map(x -> hvcat(length(uc.basis), permutedims(resolvedInts, [2, 1])...), resolvedInts)  ##### Reshaping the resolved interactions to be a vector of matrices, with the first index being the momentum index, and the second being the sublattice x Generator index.
+        permutation     =   SwitchKroneckerBasis((length(uc.OnSiteMats), length(uc.basis)))
+        map!(x -> x[permutation, permutation], resolvedInts, resolvedInts)  ##### Switching the basis of the resolved interactions to be a vector of matrices into the generator x sublattice index
+        resolvedInts    =   permutedims(reshape(repeat(resolvedInts, length(chi.Omegas)), (length(chi.Qs), length(chi.Omegas))), [2, 1])    ##### Repeating the resolved interactions to be a matrix of matrices, with the outermatrix having the (energy, momentum) index, and the inner being sublattice x Generator index.
+        corrections     =   inv(Ref(I) .- resolvedInts .* chi.Bare)  ##### The RPA corrections to the bare susceptibility coming from resummation
+        chi.RPA         =   chi.Bare .* corrections  ##### The RPA susceptibility
+    end
+
 @doc """
 ```julia
 FillChis!(chi::Susceptibility; fill_BZ::Bool=false, a::Int64=3, b::Int64=3)
@@ -126,21 +196,35 @@ FillChis!(chi::Susceptibility; fill_BZ::Bool=false, a::Int64=3, b::Int64=3)
 function to calculate susceptibility at a all given Ω=`Omegas`, but for all `Q` present in the given path, and along a fixed direction given by `a` and `b`.
 
 """
-    function FillChis!(chi::Susceptibility; fill_BZ::Bool=false, a::Int64=3, b::Int64=3, SpinMatrix::Vector{Matrix{ComplexF64}} = SpinMats((chi.M.uc.localDim-1)//2))
-        @assert prod(chi.M.bz.gridSize .% 2) == 1  "Only works if the Gamma point is included in the grid ---> grid size must be odd!"
+    function FillChis!(chi::Susceptibility, M::Model; a::Int64=3, b::Int64=3, reduce::Bool = true)
+        Gamma_index = GetQIndex(zeros(Float64, length(M.uc.primitives)), M.bz ; nearest = false)    ##### The index of the Gamma point in the BZ grid.
+        @assert !isempty(Gamma_index) "The Gamma point needs to be in the BZ grid for this calculation to work!" 
 
         if fill_BZ
-            chis    =   FindChi_FullBZ.(chi.Omegas, Ref(chi.M) ; a=a, b=b, eta=chi.Spread, SpinMatrix = SpinMatrix)
-            @cast chisNew[i, j, k] |= chis[i][j, k] / length(chi.M.bz.ks);
-        else
-            chis    =   FindChi_Path.(chi.Omegas, Ref(chi.M), Ref(chi.Qs) ; a=a, b=b, eta=chi.Spread, SpinMatrix = SpinMatrix)
-            @cast chisNew[i, j] |= chis[i][j] / length(chi.M.bz.ks);
+            chi.Qs  =   reshape(M.bz.ks, length(M.bz.ks))
         end
 
-        chi.chis[directions[a] * directions[b]]     =   chisNew
-        @info "Chis filled along $(directions[a])-$(directions[b])"
+        if reduce   ##### Calculate only the reduced susceptibility
+            chis    =   FindReducedChi.(chi.Omegas, adjoint(chi.Qs), Ref(chi.M), ; a=a, b=b, eta=chi.Spread)
+            chis    =   chis ./ length(M.bz.ks)
+
+            setindex!.(chi.BareReduced, chis, Ref(CartesianIndex(a, b)))
+        else    ##### Calculate the sublattice resolved susceptibility
+            chis    =   FindChi.(chi.Omegas, adjoint(chi.Qs), Ref(chi.M), ; a=a, b=b, eta=chi.Spread)
+            chis    =   chis ./ length(M.bz.ks)
+            ##### Setting the appropriate indices of the susceptibility matrix which is generator x sublattice.
+            for (i, base) in enumerate(M.uc.basis)
+                for (j, target) in enumerate(M.uc.basis)
+
+                    setindex!.(chi.Bare, Ref(getindex.(chis, i, j)), Ref(CartesianIndex(length(M.uc.basis) * (a - 1) + i, length(M.uc.basis) * (b - 1) + j)))
+                end
+            end
+        end
+        
+        @info "Chis filled along $(a)-$(b)."
     end
 
+
 end