Get the groundstate for a Hamiltonian (not VQE!) using Yao.jl

j1j2_groundstate.jl

using KrylovKit, Yao

"""
    J1J2{D}
    J1J2(size::Int...; J2::Real, periodic::Bool) -> J1J2

D ∈ {1, 2} is the dimension target system.
"""
struct J1J2{D}
    size::NTuple{D, Int}
    periodic::Bool
    J2::Float64
    J1J2(size::Int...; J2::Real, periodic::Bool) = new{length(size)}(size, periodic, Float64(J2))
end

nspin(model::J1J2) = prod(model.size)

"""
    hamiltonian(model::J1J2) -> MatrixBlock

Get the Hamiltonian operator for J1-J2 model.
"""
function hamiltonian(model::J1J2)
    nbit = nspin(model)
    sum(x->x[3]*heisenberg_ij(nbit, x[1], x[2]), get_bonds(model))*0.25
end
heisenberg_ij(nbit::Int, i::Int, j::Int=i+1) = put(nbit, i=>X)*put(nbit, j=>X) + put(nbit, i=>Y)*put(nbit, j=>Y) + put(nbit, i=>Z)*put(nbit, j=>Z)

"""get weighted bonds for J1J2 model"""
function get_bonds(model::J1J2{2})
    m, n = model.size
    cis = LinearIndices(model.size)
    bonds = Tuple{Int, Int, Float64}[]
    for i=1:m, j=1:n
        for (_i, _j) in [(i+1, j), (i, j+1)]
            sites = get_site((_i, _j), (m, n), Val(model.periodic))
            if all(sites .> 0)
                push!(bonds, (cis[i,j], cis[sites...], 1.0))
            end
        end
        for (_i, _j) in [(i-1, j-1), (i-1, j+1)]
            sites = get_site((_i, _j), (m, n), Val(model.periodic))
            if all(sites .> 0)
                push!(bonds, (cis[i,j], cis[sites...], model.J2))
            end
        end
    end
    bonds
end
@inline function get_site(ij, mn, pbc::Val{true})
    Tuple(mod(i-1,m)+1 for (i,m) in zip(ij, mn))
end

@inline function get_site(ij, mn, pbc::Val{false})
    Tuple(i<=m ? i : 0 for (i,m) in zip(ij, mn))
end

function get_bonds(model::J1J2{1})
    nbit, = model.size
    vcat([(i, i%nbit+1, 1.0) for i in 1:(model.periodic ? nbit : nbit-1)], [(i, (i+1)%nbit+1, model.J2) for i in 1:(model.periodic ? nbit : nbit-2)])
end

# construct a J1J2 model, with J2=0.5 and open boundary.
model = J1J2(4, 4; periodic=false, J2=0.5)
# construct the hamiltonian, try to print it!
hami = hamiltonian(model)
# Krylov eigensolver in KrylovKit
E, v = eigsolve(mat(hami), 1, :SR)

using Test
@test E[1] ≈ -7.505556950081073