added additional helper functions in base as well as a J1J2 example

examples/J1J2Model.jl

@@ -1,16 +1,48 @@
-using fermions, Serialization
+using Fermions, Serialization, Plots
 include("../src/modelHamiltonians.jl")
 
-totalSites = 33
-@assert totalSites % 3 == 0
-partitions = collect(3:3:totalSites)
-J1byJ2Values = 0:0.1:2
-gapValues = zeros(length(J1byJ2Values))
-@time for (i, J1byJ2) in enumerate(J1byJ2Values)
-    hamiltonian = J1J2Model(J1byJ2, totalSites)
+global maxSize = 2000
+
+function GetGap(
+        J1J2Ratio::Float64,
+        totalSites::Int64,
+    )
+    partitions = collect(2:2:totalSites)
+    hamiltonian = J1J2Model(J1J2Ratio, totalSites)
     hamiltonianFlow = MinceHamiltonian(hamiltonian, partitions)
-    savePaths, results = IterDiag(hamiltonianFlow, 2000)
-    eigVals = deserialize(savePaths[end-1])["eigVals"]
-    gapValues[i] = eigVals[2] - eigVals[1]
+    savePaths, results = IterDiag(hamiltonianFlow, maxSize; symmetries=['N'], silent=true)
+    energies = deserialize(savePaths[end-1])["eigVals"]
+    quantumNos = first.(deserialize(savePaths[end-1])["quantumNos"]) .- totalSites/2
+    display(quantumNos[1:5])
+    return energies[findfirst(==(0.), quantumNos)] - energies[findfirst(==(1.), quantumNos)]
+end
+
+totalSitesList = [80] # round.(Int, 10 .^ (1:0.3:2))
+for (i, totalSites) in enumerate(totalSitesList)
+    if totalSites % 2 ≠ 0
+        totalSitesList[i] += 1
+    end
+end
+display(totalSitesList)
+J1J2RatioCrit = zeros(length(totalSitesList))
+for (i, totalSites) in enumerate(totalSitesList)
+    J1J2RatioBounds = [0., 4.]
+    tolerance = 1e-2
+    gapBounds = [GetGap(J1J2RatioBounds[1], totalSites), GetGap(J1J2RatioBounds[2], totalSites)]
+    println(gapBounds)
+    @assert gapBounds[1] * gapBounds[2] < 0
+    while abs(J1J2RatioBounds[2] - J1J2RatioBounds[1]) > tolerance
+        midPoint = 0.5 * sum(J1J2RatioBounds)
+        gap = GetGap(midPoint, totalSites)
+        if gap * gapBounds[1] > 0
+            J1J2RatioBounds[1] = midPoint
+            gapBounds[1] = gap
+        else
+            J1J2RatioBounds[2] = midPoint
+            gapBounds[2] = gap
+        end
+    end
+    J1J2RatioCrit[i] = sum(J1J2RatioBounds)/2
+    println(J1J2RatioCrit[i])
 end
-display(gapValues)
+plot(1 ./ totalSitesList, J1J2RatioCrit; xscale=:log10)

src/base.jl

@@ -562,3 +562,36 @@ function Dagger(
     return operator
 end
 export Dagger
+
+
+"""
+    VacuumState(basisStates)
+
+Returns the vacuum (completely unoccupied) state for the 
+
+# Examples
+```jldoctest
+julia> operator = [("++", [1, 2], 1.), ("+n", [3, 4], 2.)]
+2-element Vector{Tuple{String, Vector{Int64}, Float64}}:
+ ("++", [1, 2], 1.0)
+ ("+n", [3, 4], 2.0)
+
+julia> Dagger(operator)
+2-element Vector{Tuple{String, Vector{Int64}, Float64}}:
+ ("--", [2, 1], 1.0)
+ ("n-", [4, 3], 2.0)
+ ```
+"""
+function VacuumState(
+        basisStates::Vector{Dict{BitVector,Float64}};
+        tolerance::Float64=1e-14,
+    )
+    numSites = basisStates |> first |> keys |> first |> length
+    numberOperators = [("n", [i], 1.) for i in 1:numSites]
+    for state in basisStates
+        if GenCorrelation(state, numberOperators) < tolerance
+            return state
+        end
+    end
+    @assert false "Vacuum state not found!"
+end

src/modelHamiltonians.jl

@@ -382,11 +382,11 @@ export Dispersion
 
 
 function J1J2Model(
-        J1byJ2::Float64,
+        J2J1Ratio::Float64,
         numSites::Int64
     )
-    J2 = 1.
-    J1 = J2 * J1byJ2
+    J1 = 1.
+    J2 = J1 * J2J1Ratio
     hamiltonian = Tuple{String, Vector{Int64}, Float64}[]
     for site in 1:numSites-1
         # J1 terms
@@ -409,3 +409,17 @@ function J1J2Model(
     end
     return hamiltonian
 end
+
+
+function CollatzModel(
+        numSites::Int64,
+    )
+    hamiltonian = Tuple{String, Vector{Int64}, Float64}[]
+    for site in 2:2:numSites
+        push!(hamiltonian, ("+-", [div(site, 2), site], 1.))
+    end
+    for site in 1:2:div(numSites-1, 3)
+        push!(hamiltonian, ("+-", [site * 3 + 1, site], 1.))
+    end
+    return hamiltonian
+end

src/thermalisation.jl

@@ -40,3 +40,22 @@ function OTOC(
     return otoc
 end
 export OTOC
+
+
+function StateEvolution(
+        initState::Vector{ComplexF64}, 
+        hamiltonian::Matrix{Float64},
+        timeSteps::Vector{Float64},
+    )
+    stateEvolution = Vector{Complex}[initState]
+    for step in eachindex(timeSteps)[2:end]
+        deltaTime = timeSteps[step] - timeSteps[step-1]
+        newState = inv(I + 0.5im * deltaTime * hamiltonian) * (I - 0.5im * deltaTime * hamiltonian) * stateEvolution[end]
+        push!(stateEvolution, newState)
+    end
+    for (i, v) in enumerate(stateEvolution)
+        stateEvolution[i] = v/norm(v)
+    end
+    return stateEvolution
+end
+export StateEvolution