Dissipative TEBD

README.md

@@ -385,11 +385,13 @@ If you use PastaQ.jl in your work, for now please consider citing the Github pag
 If you used PastaQ.jl and your paper does not appear in this list, please let us know at [[email protected]](mailto:[email protected]).
 
 **2022**   
+- [2204.10061](https://arxiv.org/abs/2204.10061) *Scalable measure of magic for quantum computers*, T Haug and MS Kim.
+- [2204.03454](https://arxiv.org/abs/2204.03454) *Variational dynamics as a ground-state problem on a quantum computer*,S Barison, F Vicentini, I Cirac and G Carleo. 
 - [2203.04948](https://arxiv.org/abs/2203.04948) *Fragile boundaries of tailored surface codes*, O Higgott, TC Bohdanowicz, A Kubica, ST Flammia, ET Campbell.     
 
 **2021**   
 - [2106.12627](https://arxiv.org/abs/2106.12627) *Provably efficient machine learning for quantum many-body problems*, H-Y Huang, R Kueng, G Torlai, VV Albert, J Preskill.     
-- [2106.03769](https://arxiv.org/abs/2106.03769) *Measurement-induced phase transition in trapped-ion circuits*, S Czischek, G Torlai, S Ray, R Islam, RG Melko, *Phys. Rev. A 104, 062405*.       
+- [2106.03769](https://arxiv.org/abs/2106.03769) *Simulating a measurement-induced phase transition in trapped-ion circuits*, S Czischek, G Torlai, S Ray, R Islam, RG Melko, *Phys. Rev. A 104, 062405*.       
 
 **2020**    
 - [2009.01760](https://arxiv.org/abs/2009.01760) *Classical variational simulation of the Quantum Approximation Optimization Algorithm*, M Medvidovic and G Carleo, *Nature Communication, 7, 101*.    

src/circuits/gates.jl

@@ -38,6 +38,24 @@ op(::OpName"a†b", st::SiteType"Qubit") = _op(OpName("a†b"), SiteType("Qudit"
 op(::OpName"ab†", st::SiteType"Qubit") = _op(OpName("ab†"), SiteType("Qudit"); dim=(2, 2))
 op(::OpName"a†b†", st::SiteType"Qubit") = _op(OpName("a†b†"), SiteType("Qudit"); dim=(2, 2))
 
+op(sites::Vector{<:Index}, O⃗::Tuple{AbstractString,Vararg}...) = 
+  reduce(*, [op(sites, O) for O in O⃗])
+
+function op(
+  sites::Vector{<:Index},
+  f::Function,
+  O⃗::Tuple{AbstractString,Vararg}...;
+  kwargs...
+)
+  operator = op(sites, O⃗...; kwargs...)
+  return f(operator)
+end
+
+op(sites::Vector{<:Index}, 
+   fO⃗::Tuple{Function,<:Tuple}) = 
+  op(sites, first(fO⃗), last(fO⃗)...)
+
+
 function phase(v::AbstractVector{ElT}) where {ElT<:Number}
   for x in v
     absx = abs(x)

src/circuits/trottersuzuki.jl

@@ -19,14 +19,82 @@ sort_gates(gates) = sort(gates; by=sort_gates_by, lt=sort_gates_lt)
 WORKING WITH TUPLES  (TEMPORARY)
 """
 
-# simplified version
-function trotter1(H::Vector{<:Tuple}, δτ::Number)
+function trotter1(δτ::Number, H::Vector{<:Tuple}; kwargs...)
   layer = Tuple[]
   for k in 1:length(H)
     length(H[k]) > 3 &&
       error("Only the format (coupling, opname, support) currently allowed")
     coupling, Hdata... = H[k]
-    layer = vcat(layer, [(x -> exp(-δτ * coupling * x), Hdata...)])
+    opname = first(Hdata)
+    layer = vcat(layer, [(x -> exp(-δτ * coupling * x), Hdata)]) 
+  end
+  return layer
+end
+
+function _lindblad_terms(δτ, hilbert, lindblad)
+  rate, opname, site = lindblad
+  !(site isa Int) && error("Only single-body lindblad operators allowed")
+  s = hilbert[site]
+
+  L = array(gate(opname, s))
+  G = -im * δτ * rate * kron(conj(L), L)
+
+  expG = reshape(exp(G),(size(L)..., size(L)...))
+  expG = reshape(permutedims(expG, (1,3,2,4)), size(G))
+  @assert ishermitian(expG)
+
+  λ, U = eigen(expG)
+  λsqrt = diagm(sqrt.(λ .+ atol))
+  K = U * λsqrt 
+  K = reshape(K, (size(L)..., size(K)[2]))
+  krausind = Index(size(K)[3]; tags="kraus")
+
+  T = ITensors.itensor(K, prime(s), ITensors.dag(s), krausind)
+
+  L2 = transpose(L) * conj(L)
+  R = exp(0.5 * im * rate * δτ * op(L2, s))
+  return [T, R]
+end
+
+function trotter1(δτ::Number, hilbert::Vector{<:Index}, H::Vector{<:Tuple}; lindbladians = [], atol = 1e-15, kwargs...)
+  layer = buildcircuit(hilbert, trotter1(δτ, H))
+  if !isempty(lindbladians)
+    for lindblad in lindbladians
+     # layer = vcat(layer, _lindblad_terms(δτ, hilbert, lindblad; atol = atol))
+      rate, opname, site = lindblad
+      !(site isa Int) && error("Only single-body lindblad operators allowed")
+      s = hilbert[site]
+
+      L = array(gate(opname, s))
+      G = -im * δτ * rate * kron(conj(L), L)
+
+      expG = reshape(exp(G),(size(L)..., size(L)...))
+      expG = reshape(permutedims(expG, (1,3,2,4)), size(G))
+      @assert expG ≈ adjoint(expG) atol = 1e-10
+
+
+      #@assert ishermitian(expG)
+      λ, U = eigen(expG)
+      λsqrt = diagm(sqrt.(λ .+ atol))
+      K = U * λsqrt 
+
+      ∑ = zeros(size(K,1), size(K,2))
+      for α in 1:size(K, 3)
+        ∑ += adjoint(K[:,:,α]) * K[:,:,α]
+      end
+      display(∑) 
+
+
+
+      K = reshape(K, (size(L)..., size(K)[2]))
+      krausind = Index(size(K)[3]; tags="kraus")
+      T = ITensors.itensor(K, prime(s), ITensors.dag(s), krausind)
+      layer = vcat(layer, [T])
+
+      R = transpose(L) * conj(L)
+      T = exp(0.5 * im * rate * δτ * op(R, s))
+      layer = vcat(layer, [T])
+    end
   end
   return layer
 end
@@ -35,18 +103,18 @@ end
     trotter2(H::OpSum; δt::Float64=0.1, δτ=im*δt)
 Generate a single layer of gates for one step of 2nd order TEBD.
 """
-function trotter2(H::Vector{<:Tuple}, δτ::Number)
-  tebd1 = trotter1(H, δτ / 2)
+function trotter2(δτ::Number, args...; kwargs...)
+  tebd1 = trotter1(δτ / 2, args...; kwargs...)
   tebd2 = vcat(tebd1, reverse(tebd1))
   return tebd2
 end
 
-function trotter4(H::Vector{<:Tuple}, δτ::Number)
+function trotter4(δτ::Number, args...; kwargs...)
   δτ1 = δτ / (4 - 4^(1 / 3))
   δτ2 = δτ - 4 * δτ1
 
-  tebd2_δ1 = trotter2(H, δτ1)
-  tebd2_δ2 = trotter2(H, δτ2)
+  tebd2_δ1 = trotter2(δτ1, args...; kwargs...)
+  tebd2_δ2 = trotter2(δτ2, args...; kwargs...)
 
   tebd4 = vcat(tebd2_δ1, tebd2_δ1)
   tebd4 = vcat(tebd4, tebd2_δ2)
@@ -58,41 +126,59 @@ end
     trotterlayer(H::OpSum; order::Int = 2, kwargs...)
 Generate a single layer of gates for one step of TEBD.
 """
-function trotterlayer(H::Vector{<:Tuple}, δτ::Number; order::Int=2)
-  order == 1 && return trotter1(H, δτ)
-  order == 2 && return trotter2(H, δτ)
+function trotterlayer(args...; order::Int=2, kwargs...)
+  order == 1 && return trotter1(args...; kwargs...)
+  order == 2 && return trotter2(args...; kwargs...)
   return error("Automated Trotter circuits with order > 2 not yet implemented")
   # TODO: understand weird behaviour of trotter4
   #order == 4 && return trotter4(H, δτ)
   #error("Automated Trotter circuits with order > 2 not yet implemented")
 end
 
-function _trottercircuit(
-  H::Vector{<:Vector{Tuple}}, τs::Vector; order::Int=2, layered::Bool=false, kwargs...
-)
-  @assert length(H) == (length(τs) - 1) || length(H) == length(τs)
+function _trottercircuit(H::Vector{<:Vector{Tuple}}, τs::Vector; layered::Bool = true, lindbladians = [], kwargs...)
+  !isempty(lindbladians) && error("Trotter simulation with Lindblad operators requires a set of indices")
+  @assert length(H) == (length(τs) -1) || length(H) == length(τs)
   δτs = diff(τs)
-  circuit = [trotterlayer(H[t], δτs[t]; order=order) for t in 1:length(δτs)]
+  circuit = [trotterlayer(δτs[t], H[t]; kwargs...) for t in 1:length(δτs)]
   layered && return circuit
   return reduce(vcat, circuit)
 end
 
-#XXX simplified version for Zygote
-function _trottercircuit(
-  H::Vector{<:Tuple}, τs::Vector; order::Int=2, layered::Bool=false, kwargs...
-)
+function _trottercircuit(hilbert::Vector{<:Index}, H::Vector{<:Vector{Tuple}}, τs::Vector; layered::Bool = true, kwargs...)
+  @assert length(H) == (length(τs) -1) || length(H) == length(τs)
+  δτs = diff(τs)
+  circuit = [trotterlayer(δτs[t], hilbert, H[t]; kwargs...) for t in 1:length(δτs)]
+  layered && return circuit
+  return reduce(vcat, circuit)
+end
+
+function _trottercircuit(H::Vector{<:Tuple}, τs::Vector; layered::Bool = true, lindbladians = [], kwargs...)
+  !isempty(lindbladians) && error("Trotter simulation with Lindblad operators requires a set of indices")
+
+  nlayers = length(τs) - 1
+  Δ = τs[2] - τs[1]
+  layer = trotterlayer(Δ, H; kwargs...)
+  layered && return [layer for _ in 1:nlayers]
+  return reduce(vcat, [layer for _ in 1:nlayers])
+end
+
+function _trottercircuit(hilbert::Vector{<:Index}, H::Vector{<:Tuple}, τs::Vector; layered::Bool = true, kwargs...)
   nlayers = length(τs) - 1
-  # XXX: Zygote: this breaks (?)
-  #circuit = [trotterlayer(H, τ; order = order) for τ in τs]
-  #!layered && return reduce(vcat, circuit)
-  #return circuit
   Δ = τs[2] - τs[1]
-  layer = trotterlayer(H, Δ; order=order)
+  layer = trotterlayer(Δ, hilbert, H; kwargs...)
   layered && return [layer for _ in 1:nlayers]
   return reduce(vcat, [layer for _ in 1:nlayers])
 end
 
-trottercircuit(H; kwargs...) = _trottercircuit(H, get_times(; kwargs...); kwargs...)
+function trottercircuit(args...; kwargs...)
+  return _trottercircuit(args..., get_times(; kwargs...); kwargs...)
+end
+
+function get_times(;
+  δt=nothing, δτ=nothing, t=nothing, τ=nothing, ts=nothing, τs=nothing, kwargs...
+)
+  return get_times(δt, δτ, t, τ, ts, τs)
+end
 
 function get_times(;
   δt=nothing, δτ=nothing, t=nothing, τ=nothing, ts=nothing, τs=nothing, kwargs...

src/itensor.jl

@@ -69,38 +69,55 @@ end
 isapprox(x::AbstractMPS, y::ITensor; kwargs...) = isapprox(prod(x), y; kwargs...)
 isapprox(x::ITensor, y::AbstractMPS; kwargs...) = isapprox(y, x; kwargs...)
 
-function expect(T₀::ITensor, ops::AbstractString...; kwargs...)
+function expect(T₀::ITensor, ops; kwargs...)
   T = copy(T₀)
-  N = nsites(T)
-  ElT = real(ITensors.promote_itensor_eltype([T]))
-  Nops = length(ops)
+  s = inds(T; plev=0)
+  N = length(s)
 
-  site_range::UnitRange{Int} = get(kwargs, :site_range, 1:N)
+  ElT = ITensors.promote_itensor_eltype([T])
+
+  is_operator = !isempty(inds(T; plev=1))
+
+  if haskey(kwargs, :site_range)
+    @warn "The `site_range` keyword arg. to `expect` is deprecated: use the keyword `sites` instead"
+    sites = kwargs[:site_range]
+  else
+    sites = get(kwargs, :sites, 1:N)
+  end
+
+  site_range = (sites isa AbstractRange) ? sites : collect(sites)
   Ns = length(site_range)
   start_site = first(site_range)
-  offset = start_site - 1
 
-  normalization = is_operator(T) ? tr(T) : norm(T)^2
+  el_types = map(o -> ishermitian(op(o, s[start_site])) ? real(ElT) : ElT, ops)
+
+  normalization = is_operator ? tr(T) : norm(T)^2
 
-  ex = ntuple(n -> zeros(ElT, Ns), Nops)
-  for j in site_range
-    for n in 1:Nops
-      s = firstind(T; tags="Site, n=$j", plev=0)
-      if is_operator(T)
-        Top = replaceprime(T * op(ops[n], s'), 2 => 1; tags="Site, n=$j")
-        ex[n][j - offset] = real(tr(Top) / normalization)
+  ex = map((o, el_t) -> zeros(el_t, Ns), ops, el_types)
+  for (entry, j) in enumerate(site_range)
+    for (n, opname) in enumerate(ops)
+      if is_operator
+        val = replaceprime(op(opname, s[j])' * T, 2 => 1; inds=s[j]'')
+        val = tr(val) / normalization
       else
-        ex[n][j - offset] =
-          real(scalar(dag(T) * noprime(op(ops[n], s) * T))) / normalization
+        val = scalar(dag(T) * noprime(op(opname, s[j]) * T)) / normalization
       end
+      ex[n][entry] = (el_types[n] <: Real) ? real(val) : val
     end
   end
 
-  if Nops == 1
-    return Ns == 1 ? ex[1][1] : ex[1]
-  else
-    return Ns == 1 ? [x[1] for x in ex] : ex
+  if sites isa Number
+    return map(arr -> arr[1], ex)
   end
+  return ex
+end
+
+function expect(T::ITensor, op::AbstractString; kwargs...)
+  return first(expect(T, (op,); kwargs...))
+end
+
+function expect(T::ITensor, op1::AbstractString, ops::AbstractString...; kwargs...)
+  return expect(T, (op1, ops...); kwargs...)
 end
 
 @non_differentiable ITensors.name(::Any)

test/gates.jl

@@ -673,21 +673,22 @@ end
   @test PastaQ.array(gtest) ≈ G
 end
 
-@testset "function applied to a gate" begin
-  s = siteinds("Qubit", 2)
-
-  θ = 0.1
-  rx = array(gate("Rx", s[1]; θ=0.1))
-  exp_rx = exp(rx)
-  gtest = gate(x -> exp(x), "Rx", s[1]; θ=0.1)
-  @test exp_rx ≈ array(gate(x -> exp(x), "Rx", s[1]; θ=0.1))
-  @test exp_rx ≈ array(gate(x -> exp(x), ("Rx", 1, (θ=0.1,)), s))
-
-  cx = 0.1 * reshape(array(gate("CX", s[1], s[2])), (4, 4))
-  exp_cx = reshape(exp(cx), (2, 2, 2, 2))
-  @test exp_cx ≈ array(gate(x -> exp(0.1 * x), "CX", s[1], s[2]))
-  @test exp_cx ≈ array(gate(x -> exp(0.1 * x), ("CX", (1, 2)), s))
-end
+#@testset "function applied to a gate" begin
+#  s = siteinds("Qubit", 2)
+#  
+#  θ = 0.1
+#  rx = array(gate("Rx", s[1]; θ = 0.1))
+#  exp_rx = exp(rx)
+#  gtest = gate(x -> exp(x), "Rx",s[1]; θ = 0.1)
+#  @test exp_rx ≈ array(gate(x -> exp(x), "Rx",s[1]; θ = 0.1))
+#  @test exp_rx ≈ array(gate(x -> exp(x), ("Rx", 1, (θ = 0.1,)), s))
+#  
+#  cx = 0.1*reshape(array(gate("CX", s[1], s[2])),(4,4))
+#  exp_cx = reshape(exp(cx),(2,2,2,2))
+#  @test exp_cx ≈ array(gate(x -> exp(0.1*x), "CX", s[1], s[2]))
+#  @test exp_cx ≈ array(gate(x -> exp(0.1*x), ("CX", (1,2)), s))
+#end
+#
 
 ##@testset "qudit gates" begin
 ##  dim = 3