Fermionic OpSum to TTN constructor

src/models.jl

@@ -1,6 +1,90 @@
 _maybe_fill(x, n) = x
 _maybe_fill(x::Number, n) = fill(x, n)
 
+#= Make the free fermion Hamiltonian for the up spins
+os_up = OpSum()
+for n in 1:(N - 1)
+  os_up .+= -t, "Cdagup", n, "Cup", n + 1
+  os_up .+= -t, "Cdagup", n + 1, "Cup", n
+end
+
+# Make the free fermion Hamiltonian for the down spins
+os_dn = OpSum()
+for n in 1:(N - 1)
+  os_dn .+= -t, "Cdagdn", n, "Cdn", n + 1
+  os_dn .+= -t, "Cdagdn", n + 1, "Cdn", n
+end
+=#
+function tight_binding(g::AbstractGraph;t=1.0, tp=0.0, h::Union{<:Real,Vector{<:Real}}=0)
+  h = _maybe_fill(h, nv(g))
+  ℋ = OpSum()
+  if !iszero(t)
+    for e in edges(g)
+      ℋ += -t, "Cdag", maybe_only(src(e)), "C", maybe_only(dst(e))
+      ℋ += -t, "Cdag", maybe_only(dst(e)), "C", maybe_only(src(e))
+    end
+  end
+  if !iszero(t')
+    # TODO, more clever way of looping over next to nearest neighbors?
+    for (i, v) in enumerate(vertices(g))
+      nnn = [neighbors(g, n) for n in neighbors(g, v)]
+      nnn = setdiff(Base.Iterators.flatten(nnn), neighbors(g, v))
+      nnn = setdiff(nnn, vertices(g)[1:i])
+      for nn in nnn
+        ℋ += -tp, "Cdag", maybe_only(v), "C", maybe_only(nn)
+        ℋ += -tp, "Cdag", maybe_only(nn), "C", maybe_only(v)
+      end
+    end
+  end
+  for (i, v) in enumerate(vertices(g))
+    if !iszero(h[i])
+      ℋ -= h[i], "N", maybe_only(v)
+    end
+  end
+  return ℋ
+end
+
+
+
+"""
+t-t' Hubbard Model 
+"""
+function hubbard(g::AbstractGraph; U=0.0, t=1.0, tp=0.0, h::Union{<:Real,Vector{<:Real}}=0)
+  h = _maybe_fill(h, nv(g))
+  ℋ = OpSum()
+  if !iszero(t)
+    for e in edges(g)
+      ℋ += -t, "Cdagup", maybe_only(src(e)), "Cup", maybe_only(dst(e))
+      ℋ += -t, "Cdagup", maybe_only(dst(e)), "Cup", maybe_only(src(e))
+      ℋ += -t, "Cdagdn", maybe_only(src(e)), "Cdn", maybe_only(dst(e))
+      ℋ += -t, "Cdagdn", maybe_only(dst(e)), "Cdn", maybe_only(src(e))
+    end
+  end
+  if !iszero(tp)
+    # TODO, more clever way of looping over next to nearest neighbors?
+    for (i, v) in enumerate(vertices(g))
+      nnn = [neighbors(g, n) for n in neighbors(g, v)]
+      nnn = setdiff(Base.Iterators.flatten(nnn), neighbors(g, v))
+      nnn = setdiff(nnn, vertices(g)[1:i])
+      for nn in nnn
+        ℋ += -tp, "Cdagup", maybe_only(v), "Cup", maybe_only(nn)
+        ℋ += -tp, "Cdagup", maybe_only(nn), "Cup", maybe_only(v)
+        ℋ += -tp, "Cdagdn", maybe_only(v), "Cdn", maybe_only(nn)
+        ℋ += -tp, "Cdagdn", maybe_only(nn), "Cdn", maybe_only(v)
+      end
+    end
+  end
+  for (i, v) in enumerate(vertices(g))
+    if !iszero(h[i])
+      ℋ -= h[i], "Sz", maybe_only(v)
+    end
+    if !iszero(U)
+      ℋ = U, "Nupdn", maybe_only(v)
+    end
+  end
+  return ℋ
+end
+
 """
 Random field J1-J2 Heisenberg model on a general graph
 """
@@ -43,7 +127,7 @@
   ℋ = OpSum()
   if !iszero(J1)
     for e in edges(g)
-      ℋ += J1, "Z", maybe_only(src(e)), "Z", maybe_only(dst(e))
+      ℋ += J1, "Sz", maybe_only(src(e)), "Sz", maybe_only(dst(e))
     end
   end
   if !iszero(J2)
@@ -53,13 +137,13 @@
       nnn = setdiff(Base.Iterators.flatten(nnn), neighbors(g, v))
       nnn = setdiff(nnn, vertices(g)[1:i])
       for nn in nnn
-        ℋ += J2, "Z", maybe_only(v), "Z", maybe_only(nn)
+        ℋ += J2, "Sz", maybe_only(v), "Sz", maybe_only(nn)
       end
     end
   end
   for (i, v) in enumerate(vertices(g))
     if !iszero(h[i])
-      ℋ += h[i], "X", maybe_only(v)
+      ℋ += h[i], "Sx", maybe_only(v)
     end
   end
   return ℋ

src/treetensornetworks/opsum_to_ttn.jl

@@ -4,6 +4,7 @@
 # Utility methods
 # 
 
+
 # linear ordering of vertices in tree graph relative to chosen root, chosen outward from root
 function find_index_in_tree(site, g::AbstractGraph, root_vertex)
   ordering = reverse(post_order_dfs_vertices(g, root_vertex))
@@ -59,8 +60,10 @@
   vertextype_sites = vertextype(sites)
   thishasqns = any(v -> hasqns(sites[v]), vertices(sites))
 
+  linkdir_ref = ITensors.using_auto_fermion() ? ITensors.In : ITensors.Out
   # traverse tree outwards from root vertex
   vs = reverse(post_order_dfs_vertices(sites, root_vertex))                                 # store vertices in fixed ordering relative to root
+  # ToDo: Add check in ttn_svd that the ordering matches that of find_index_in_tree, which is used in sorteachterm #fermion-sign!
   es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex)))                          # store edges in fixed ordering relative to root
   # some things to keep track of
   degrees = Dict(v => degree(sites, v) for v in vs)                                           # rank of every TTN tensor in network
@@ -81,7 +84,7 @@
         op_cache[ITensors.which_op(st) => ITensors.site(st)] = op_tensor
       end
       if !isnothing(flux(op_tensor))
-        q -= flux(op_tensor)
+        q += flux(op_tensor)
       end
     end
     return q
@@ -93,13 +96,13 @@
   for v in vs
     is_internal[v] = isempty(sites[v])
     if isempty(sites[v])
-      push!(sites[v], Index(Hflux => 1)) ###FIXME: using Hflux may be a problem with AutoFermion, the dummy index has to be bosonic (but the QN of H should be too?) ...
+      push!(sites[v], Index(Hflux => 1))
     end
   end
   inbond_coefs = Dict(
     e => Dict{QN,Vector{ITensors.MatElem{coefficient_type}}}() for e in es
   )                         # bond coefficients for incoming edge channels
-  site_coef_done = Prod{Op}[]                                                               # list of terms for which the coefficient has been added to a site factor
+  site_coef_done = Prod{Op}[] # list of terms for which the coefficient has been added to a site factor
   # temporary symbolic representation of TTN Hamiltonian
   tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs)
 
@@ -139,10 +142,10 @@
 
       # compute QNs
       incoming_qn = calc_qn(incoming)
-      not_incoming_qn = calc_qn(not_incoming)  ###does this work?
+      not_incoming_qn = calc_qn(not_incoming)
       outgoing_qns = Dict(e => calc_qn(outgoing[e]) for e in edges_out)
       site_qn = calc_qn(onsite)
-
+      
       # initialize QNArrayElement indices and quantum numbers 
       T_inds = MVector{degrees[v]}(fill(-1, degrees[v]))
       T_qns = MVector{degrees[v]}(fill(QN(), degrees[v]))
@@ -151,9 +154,10 @@
       bond_col = -1
       if !isempty(incoming)
         # get the correct map from edge=>QN to term and channel
-        coutmap = get!(outmaps, edge_in => -not_incoming_qn, Dict{Vector{Op},Int}())
-        cinmap = get!(inmaps, edge_in => incoming_qn, Dict{Vector{Op},Int}())
         # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index
+        coutmap = get!(outmaps, edge_in => not_incoming_qn , Dict{Vector{Op},Int}())
+        cinmap = get!(inmaps, edge_in => -incoming_qn, Dict{Vector{Op},Int}())
+
         bond_row = ITensors.posInLink!(cinmap, incoming)
         bond_col = ITensors.posInLink!(coutmap, not_incoming) # get incoming channel
         bond_coef = convert(coefficient_type, ITensors.coefficient(term))
@@ -162,14 +166,14 @@
         )
         push!(q_inbond_coefs, ITensors.MatElem(bond_row, bond_col, bond_coef))
         T_inds[dim_in] = bond_col
-        T_qns[dim_in] = incoming_qn
+        T_qns[dim_in] = -incoming_qn 
       end
       for dout in dims_out
         coutmap = get!(
-          outmaps, edges[dout] => -outgoing_qns[edges[dout]], Dict{Vector{Op},Int}()
+          outmaps, edges[dout] => outgoing_qns[edges[dout]], Dict{Vector{Op},Int}()
         )
         T_inds[dout] = ITensors.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel
-        T_qns[dout] = -outgoing_qns[edges[dout]]  # minus sign to account for outgoing arrow direction
+        T_qns[dout] = outgoing_qns[edges[dout]]
       end
       # if term starts at this site, add its coefficient as a site factor
       site_coef = one(coefficient_type)
@@ -192,6 +196,7 @@
 
       push!(tempTTN[v], el)
     end
+
     ITensors.remove_dups!(tempTTN[v])
     # manual truncation: isometry on incoming edge
     if !isnothing(dim_in) && !isempty(inbond_coefs[edges[dim_in]])
@@ -210,20 +215,15 @@
   # compress this tempTTN representation into dense form
 
   link_space = Dict{edgetype_sites,Index}()
-  linkdir = ITensors.using_auto_fermion() ? ITensors.In : ITensors.Out ###FIXME: ??
   for e in es
     qi = Vector{Pair{QN,Int}}()
     push!(qi, QN() => 1)
     for (q, Vq) in Vs[e]
       cols = size(Vq, 2)
-      if ITensors.using_auto_fermion() # <fermions>
-        push!(qi, (-q) => cols)
-      else
-        push!(qi, q => cols)
-      end
+      push!(qi, q => cols)
     end
     push!(qi, Hflux => 1)
-    link_space[e] = Index(qi...; tags=edge_tag(e), dir=linkdir)
+    link_space[e] = Index(qi...; tags=edge_tag(e), dir=linkdir_ref)
   end
 
   H = TTN(sites0)   # initialize TTN without the dummy indices added
@@ -240,12 +240,9 @@
     edges = filter(e -> dst(e) == v || src(e) == v, es)
     dim_in = findfirst(e -> dst(e) == v, edges)
     dims_out = findall(e -> src(e) == v, edges)
-
     # slice isometries at this vertex
     Vv = [Vs[e] for e in edges]
-
     linkinds = [link_space[e] for e in edges]
-    linkdims = dim.(linkinds)
 
     # construct blocks
     blocks = Dict{Tuple{Block{degrees[v]},Vector{Op}},Array{coefficient_type,degrees[v]}}()
@@ -272,32 +269,24 @@
 
       # set non-trivial helper inds
       for d in normal_dims
-        block_helper_inds[d] = qnblock(linkinds[d], T_qns[d])
+        block_helper_inds[d] = qnblock(linkinds[d],  T_qns[d])
       end
-
       @assert all(≠(-1), block_helper_inds)# check that all block indices are set
-      # make Block 
+
+      # make and fill Block 
       theblock = Block(Tuple(block_helper_inds))
-
-      # T_qns takes into account incoming/outgoing arrow direction, while Vs are stored for incoming arrow
-      # flip sign of QN for outgoing arrows, to match incoming arrow definition of Vs
-      iT_qns = deepcopy(T_qns)
-      for d in dims_out
-        iT_qns[d] = -iT_qns[d]
-      end
-
       if isempty(normal_dims)
         M = get!(blocks, (theblock, terms(t)), zero_arr())
         @assert isone(length(M))
         M[] += ct
       else
         M = get!(blocks, (theblock, terms(t)), zero_arr())
-        dim_ranges = Tuple(size(Vv[d][iT_qns[d]], 2) for d in normal_dims)
-        for c in CartesianIndices(dim_ranges)
+        dim_ranges = Tuple(size(Vv[d][T_qns[d]], 2) for d in normal_dims)        
+        for c in CartesianIndices(dim_ranges) # applies isometries in a element-wise manner
           z = ct
           temp_inds = copy(T_inds)
           for (i, d) in enumerate(normal_dims)
-            V_factor = Vv[d][iT_qns[d]][T_inds[d], c[i]]
+            V_factor = Vv[d][T_qns[d]][T_inds[d], c[i]]
             z *= (d == dim_in ? conj(V_factor) : V_factor) # conjugate incoming isometry factor
             temp_inds[d] = c[i]
           end
@@ -308,32 +297,28 @@
 
     H[v] = ITensor()
 
+    # make the final arrow directions
     _linkinds = copy(linkinds)
     if !isnothing(dim_in)
       _linkinds[dim_in] = dag(_linkinds[dim_in])
     end
 
     for ((b, q_op), m) in blocks
       Op = computeSiteProd(sites, Prod(q_op))
-      if hasqns(Op) ###FIXME: this may not be safe, we may want to check for the equivalent (zero tensor?) case in the dense case as well
+      if hasqns(Op) # FIXME: this may not be safe, we may want to check for the equivalent (zero tensor?) case in the dense case as well
         iszero(nnzblocks(Op)) && continue
       end
       sq = flux(Op)
       if !isnothing(sq)
+        rq =(b[1]==1 ? Hflux : first(space(linkinds[1])[b[1]])) # get row (dim_in) QN
+        cq = rq - sq # get column (out_dims) QN
         if ITensors.using_auto_fermion()
-          @assert false
-          ###assuming we want index ordering dim_in, phys..., dims_out...
-          ###but are assembling the tensor by tensor product (dim_in,dims_out...) x (phys...)
-          ###we need to undo the autofermion logic.
-          ###this requires the perm and the quantum numbers
-          ###perm is 1,3,2 (switching phys and dims_out dimensions)
-          ###quantum numbers need to be retrieved from block and linkinds
-          ###basically first(space(linkinds[i])[j]) for (i,j) in block
-          ###this logic requires dim_in to be 1 always, i don't think this is guaranteed
-          ###however, maybe it is required when using autofermion?
-          ###so we would have to figure out a perm to bring it into that format
-          ##compute perm factor to bring into this format, combine with perm factor that results from MPS-like logic
-          #perm=(1,3,2)
+          # we need to account for the direct product below ordering the physical indices as the last indices
+          # although they are in between incoming and outgoing indices in the canonical site-ordering
+          perm=(1,3,2)
+          if ITensors.compute_permfactor(perm,rq,sq,cq) == -1
+            Op .*= -1
+          end
         end
       end
       T = ITensors.BlockSparseTensor(coefficient_type, [b], _linkinds)
@@ -342,9 +327,13 @@
       if !thishasqns
         iT = removeqns(iT)
       end
+
       if is_internal[v]
         H[v] += iT
       else
+        if hasqns(iT)
+          @assert flux(iT * Op) == Hflux
+        end
         H[v] += (iT * Op)
       end
     end
@@ -411,7 +400,8 @@
   return T
 end
 
-###CHECK/ToDo: Understand the fermion logic on trees...
+# This is almost an exact copy of ITensors/src/opsum_to_mpo_generic:sorteachterm except for the site ordering being
+# given via find_index_in_tree
 # changed `isless_site` to use tree vertex ordering relative to root
 function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {V}
   os = copy(os)
@@ -458,7 +448,7 @@
       prevsite = currsite
 
       if fermionic
-        error("No verified fermion support for automatic TTN constructor!") # no verified support, just throw error
+        #error("No verified fermion support for automatic TTN constructor!") # no verified support, just throw error
         parity = -parity
       else
         # Ignore bosonic operators in perm