Merge branch 'main' of github.com:Anjishnubose/TightBindingToolkit.jl

Project.toml

@@ -1,9 +1,10 @@
 name = "TightBindingToolkit"
 uuid = "2233325a-6eb3-486f-aff0-670e0939fa7e"
 authors = ["Anjishnu Bose", "Sreekar Voleti", "Samuel Vadnais"]
-version = "2.3.0"
+version = "2.3.2"
 
 [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]

Sample/Vortex.jl

@@ -19,13 +19,15 @@ thirdNNdistance     = 2/sqrt(3)
 
 UC      =   UnitCell([a1, a2], 2, 2)
 UC.BC   =   zeros(ComplexF64, 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  =     -10.0
+const tInter  =     -1.5
 
 tIntraParam     =   Param(tIntra, 2)
 AddAnisotropicBond!(tIntraParam, UC, 1, 2, [ 0, 0], SpinVec[4], firstNNdistance, "Intra hopping")
@@ -48,38 +50,77 @@ CreateUnitCell!(UC, [tIntraParam, tInterParam])
 
 UCPlot = 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")
 Plot_Lattice!(lattice ; bond_cmp=:viridis, site_size=10.0)
 
 # M   =   LatticeModel(lattice, H)
 
-######################## Symmetry ##########################
+# ######################## Symmetry ##########################
 
+# 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))
 # mGround = mod.([sum(sum.(m1)), sum(sum.(m2))], Ref(1.0))

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

src/Flux.jl

@@ -1,9 +1,10 @@
 module Flux
 
-    export GetBondPhase, CheckGaugeValidity, ModifyGauge!
+    export GetBondPhase, CheckGaugeValidity, ModifyGauge!, GetStringPhase, InsertMonopolePair!
 
     using LinearAlgebra
 
+    using ..TightBindingToolkit.Useful: SegmentIntersection
     using ..TightBindingToolkit.LatticeStruct: Lattice, FillLattice!
 
 
@@ -41,6 +42,37 @@ module Flux
     end
 
 
+    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 && r>=0.0 && r<1.0
+            return exp(im * flux)
+        else
+            return 1.0
+        end
+    end
+
+
+    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
+
+
 
 
 

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

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

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