new hamiltonian solvers

Author: GiggleLiu

src/QuAlgorithmZoo.jl

@@ -27,6 +27,7 @@ include("RotBasis.jl")
 include("Grover.jl")
 include("PhaseEstimation.jl")
 include("HHL.jl")
+include("hamiltonian_solvers.jl")
 
 
 end # module

src/hamiltonian_solvers.jl

@@ -0,0 +1,46 @@
+export heisenberg, iter_groundstate!, itime_groundstate!
+
+"""
+    heisenberg(nbit::Int; periodic::Bool=true)
+
+heisenberg hamiltonian, for its ground state, refer `PRB 48, 6141`.
+"""
+function heisenberg(nbit::Int; periodic::Bool=true)
+    sx = i->put(nbit, i=>X)
+    sy = i->put(nbit, i=>Y)
+    sz = i->put(nbit, i=>Z)
+    mapreduce(i->(j=i%nbit+1; sx(i)*sx(j)+sy(i)*sy(j)+sz(i)*sz(j)), +, 1:(periodic ? nbit : nbit-1))
+end
+
+"""
+    iter_groundstate!({reg::AbstractRegister}, h::MatrixBlock; niter::Int=100) -> AbstractRegister
+
+project wave function to ground state by iteratively apply -h.
+"""
+iter_groundstate!(h::MatrixBlock; niter::Int=100) = reg -> iter_groundstate!(reg, h, niter=niter)
+function iter_groundstate!(reg::AbstractRegister, h::MatrixBlock; niter::Int=100)
+    for i = 1:niter
+        reg |> h
+        i%5 == 0 && reg |> normalize!
+    end
+    reg |> normalize!
+end
+
+"""
+    itime_groundstate!({reg::AbstractRegister}, h::MatrixBlock; τ::Int=20, tol=1e-4) -> AbstractRegister
+
+project wave function to ground state by exp(-hτ). `tol` is for `expmv`.
+"""
+itime_groundstate!(h::MatrixBlock; τ::Real=20, tol=1e-4) = reg -> itime_groundstate!(reg, h; τ=τ, tol=tol)
+function itime_groundstate!(reg::AbstractRegister, h::MatrixBlock; τ::Int=20, tol=1e-4)
+    span = 1
+    te = timeevolve(h, -im*span)
+    for i = 1:τ÷span
+        reg |> te |> normalize!
+    end
+    if τ%span != 0
+        reg |> timeevolve(h, τ%span) |> normalize!
+    end
+    reg
+end
+

test/hamiltonian_solvers.jl

@@ -0,0 +1,22 @@
+using Yao, Yao.Blocks
+using LinearAlgebra
+using Test
+using QuAlgorithmZoo
+
+@testset "solving hamiltonian" begin
+    nbit = 8
+    h = heisenberg(nbit) |> cache
+    @test ishermitian(h)
+    reg = rand_state(nqubits(h))
+    E0 = expect(h, reg)/nbit
+    reg |> iter_groundstate!(h, niter=1000)
+    EG = expect(h, reg)/nbit/4
+    @test isapprox(EG, -0.4564, atol=1e-4)
+
+    # using Time Evolution
+    reg = rand_state(nqubits(h))
+    reg |> itime_groundstate!(h, τ=20)
+    EG = expect(h, reg)/nbit/4
+    @test isapprox(EG, -0.4564, atol=1e-4)
+end
+

test/runtests.jl

@@ -32,3 +32,7 @@ end
 @testset "QCOptProblem" begin
     include("QCOptProblem.jl")
 end
+
+@testset "hamiltonian solvers" begin
+    include("hamiltonian_solvers.jl")
+end