Add Model*s submodules

src/ITensorNetworks.jl

@@ -14,7 +14,6 @@ include("sitetype.jl")
 include("abstractitensornetwork.jl")
 include("contraction_sequences.jl")
 include("expect.jl")
-include("models.jl")
 include("tebd.jl")
 include("itensornetwork.jl")
 include("mincut.jl")
@@ -67,6 +66,8 @@ include("solvers/sweep_plans/sweep_plans.jl")
 include("apply.jl")
 include("environment.jl")
 include("exports.jl")
+include("ModelHamiltonians/ModelHamiltonians.jl")
+include("ModelNetworks/ModelNetworks.jl")
 
 using PackageExtensionCompat: @require_extensions
 using Requires: @require

src/ModelHamiltonians/ModelHamiltonians.jl

@@ -0,0 +1,127 @@
+module ModelHamiltonians
+using Dictionaries: Dictionary
+using Graphs: AbstractGraph, dst, edges, edgetype, neighborhood, path_graph, src, vertices
+using ITensors.Ops: OpSum
+
+to_callable(value::Type) = value
+to_callable(value::Function) = value
+to_callable(value::Dict) = Base.Fix1(getindex, value)
+to_callable(value::Dictionary) = Base.Fix1(getindex, value)
+to_callable(value) = Returns(value)
+
+# TODO: Move to `NamedGraphs.jl` or `GraphsExtensions.jl`.
+# TODO: Add a tet for this.
+function nth_nearest_neighbors(g, v, n::Int)
+  isone(n) && return neighborhood(g, v, 1)
+  return setdiff(neighborhood(g, v, n), neighborhood(g, v, n - 1))
+end
+
+# TODO: Move to `NamedGraphs.jl` or `GraphsExtensions.jl`.
+# TODO: Add a tet for this.
+next_nearest_neighbors(g, v) = nth_nearest_neighbors(g, v, 2)
+
+function tight_binding(g::AbstractGraph; t=1, tp=0, h=0)
+  (; t, tp, h) = map(to_callable, (; t, tp, h))
+  h = to_callable(h)
+  ℋ = OpSum()
+  for e in edges(g)
+    ℋ -= t(e), "Cdag", src(e), "C", dst(e)
+    ℋ -= t(e), "Cdag", dst(e), "C", src(e)
+  end
+  for v in vertices(g)
+    for nn in next_nearest_neighbors(g, v)
+      e = edgetype(g)(v, nn)
+      ℋ -= tp(e), "Cdag", v, "C", nn
+      ℋ -= tp(e), "Cdag", nn, "C", v
+    end
+  end
+  for v in vertices(g)
+    ℋ -= h(v), "N", v
+  end
+  return ℋ
+end
+
+"""
+t-t' Hubbard Model g,i,v
+"""
+function hubbard(g::AbstractGraph; U=0, t=1, tp=0, h=0)
+  (; U, t, tp, h) = map(to_callable, (; U, t, tp, h))
+  ℋ = OpSum()
+  for e in edges(g)
+    ℋ -= t(e), "Cdagup", src(e), "Cup", dst(e)
+    ℋ -= t(e), "Cdagup", dst(e), "Cup", src(e)
+    ℋ -= t(e), "Cdagdn", src(e), "Cdn", dst(e)
+    ℋ -= t(e), "Cdagdn", dst(e), "Cdn", src(e)
+  end
+  for v in vertices(g)
+    for nn in next_nearest_neighbors(g, v)
+      e = edgetype(g)(v, nn)
+      ℋ -= tp(e), "Cdagup", v, "Cup", nn
+      ℋ -= tp(e), "Cdagup", nn, "Cup", v
+      ℋ -= tp(e), "Cdagdn", v, "Cdn", nn
+      ℋ -= tp(e), "Cdagdn", nn, "Cdn", v
+    end
+  end
+  ℋ -= h(v), "Sz", v
+  ℋ += U(v), "Nupdn", v
+  return ℋ
+end
+
+"""
+Random field J1-J2 Heisenberg model on a general graph
+"""
+function heisenberg(g::AbstractGraph; J1=1, J2=0, h=0)
+  (; J1, J2, h) = map(to_callable, (; J1, J2, h))
+  ℋ = OpSum()
+  for e in edges(g)
+    ℋ += J1(e) / 2, "S+", src(e), "S-", dst(e)
+    ℋ += J1(e) / 2, "S-", src(e), "S+", dst(e)
+    ℋ += J1(e), "Sz", src(e), "Sz", dst(e)
+  end
+  for v in vertices(g)
+    for nn in next_nearest_neighbors(g, v)
+      e = edgetype(g)(v, nn)
+      ℋ += J2(e) / 2, "S+", v, "S-", nn
+      ℋ += J2(e) / 2, "S-", v, "S+", nn
+      ℋ += J2(e), "Sz", v, "Sz", nn
+    end
+  end
+  for (i, v) in enumerate(vertices(g))
+    if !iszero(h[i])
+      ℋ += h[i], "Sz", v
+    end
+  end
+  return ℋ
+end
+
+"""
+Next-to-nearest-neighbor Ising model (ZZX) on a general graph
+"""
+function ising(g::AbstractGraph; J1=-1, J2=0, h=0)
+  (; J1, J2, h) = map(to_callable, (; J1, J2, h))
+  ℋ = OpSum()
+  for e in edges(g)
+    ℋ += J1(e), "Sz", src(e), "Sz", dst(e)
+  end
+  for v in vertices(g)
+    for nn in next_nearest_neighbors(g, v)
+      e = edgetype(g)(v, nn)
+      ℋ += J2(e), "Sz", v, "Sz", nn
+    end
+  end
+  for v in vertices(g)
+    ℋ += h(v), "Sx", v
+  end
+  return ℋ
+end
+
+"""
+Random field J1-J2 Heisenberg model on a chain of length N
+"""
+heisenberg(N::Integer; kwargs...) = heisenberg(path_graph(N); kwargs...)
+
+"""
+Next-to-nearest-neighbor Ising model (ZZX) on a chain of length N
+"""
+ising(N::Integer; kwargs...) = ising(path_graph(N); kwargs...)
+end

src/ModelNetworks/ModelNetworks.jl

@@ -0,0 +1,80 @@
+module ModelNetworks
+using Graphs: edges
+using ..ITensorNetworks: IndsNetwork
+using NamedGraphs: NamedGraph
+
+"""
+BUILD Z OF CLASSICAL ISING MODEL ON A GIVEN GRAPH AT INVERSE TEMP BETA
+H = -\\sum_{(v,v') \\in edges}\\sigma^{z}_{v}\\sigma^{z}_{v'}
+OPTIONAL ARGUMENT:
+  h: EXTERNAL MAGNETIC FIELD
+  szverts: A LIST OF VERTICES OVER WHICH TO APPLY A SZ.
+    THE RESULTANT NETWORK CAN THEN BE CONTRACTED AND DIVIDED BY THE ACTUAL PARTITION FUNCTION TO GET THAT OBSERVABLE
+    INDSNETWORK IS ASSUMED TO BE BUILT FROM A GRAPH (NO SITE INDS) AND OF LINK SPACE 2
+"""
+function ising_network(
+  eltype::Type, s::IndsNetwork, beta::Number; h::Number=0.0, szverts=nothing
+)
+  s = insert_missing_internal_inds(s, edges(s); internal_inds_space=2)
+  tn = delta_network(eltype, s)
+  if (szverts != nothing)
+    for v in szverts
+      tn[v] = diagITensor(eltype[1, -1], inds(tn[v]))
+    end
+  end
+  for edge in edges(tn)
+    v1 = src(edge)
+    v2 = dst(edge)
+    i = commoninds(tn[v1], tn[v2])[1]
+    deg_v1 = degree(tn, v1)
+    deg_v2 = degree(tn, v2)
+    f11 = exp(beta * (1 + h / deg_v1 + h / deg_v2))
+    f12 = exp(beta * (-1 + h / deg_v1 - h / deg_v2))
+    f21 = exp(beta * (-1 - h / deg_v1 + h / deg_v2))
+    f22 = exp(beta * (1 - h / deg_v1 - h / deg_v2))
+    q = eltype[f11 f12; f21 f22]
+    w, V = eigen(q)
+    w = map(sqrt, w)
+    sqrt_q = V * ITensors.Diagonal(w) * inv(V)
+    t = itensor(sqrt_q, i, i')
+    tn[v1] = tn[v1] * t
+    tn[v1] = noprime!(tn[v1])
+    t = itensor(sqrt_q, i', i)
+    tn[v2] = tn[v2] * t
+    tn[v2] = noprime!(tn[v2])
+  end
+  return tn
+end
+
+function ising_network(s::IndsNetwork, beta::Number; h::Number=0.0, szverts=nothing)
+  return ising_network(typeof(beta), s, beta; h, szverts)
+end
+
+function ising_network(
+  eltype::Type, g::NamedGraph, beta::Number; h::Number=0.0, szverts=nothing
+)
+  return ising_network(eltype, IndsNetwork(g; link_space=2), beta; h, szverts)
+end
+
+function ising_network(g::NamedGraph, beta::Number; h::Number=0.0, szverts=nothing)
+  return ising_network(eltype(beta), g, beta; h, szverts)
+end
+
+"""Build the wavefunction whose norm is equal to Z of the classical ising model
+s needs to have site indices in this case!"""
+function ising_network_state(eltype::Type, s::IndsNetwork, beta::Number; h::Number=0.0)
+  return ising_network(eltype, s, 0.5 * beta; h)
+end
+
+function ising_network_state(eltype::Type, g::NamedGraph, beta::Number; h::Number=0.0)
+  return ising_network(eltype, IndsNetwork(g, 2, 2), 0.5 * beta; h)
+end
+
+function ising_network_state(s::IndsNetwork, beta::Number; h::Number=0.0)
+  return ising_network_state(typeof(beta), s, beta; h)
+end
+
+function ising_network_state(g::NamedGraph, beta::Number; h::Number=0.0)
+  return ising_network(typeof(beta), IndsNetwork(g, 2, 2), 0.5 * beta; h)
+end
+end