Rewrite models based on Square struct, add LineMPO struct

src/ITensorNetworks/ITensorNetworks.jl

@@ -5,9 +5,9 @@ using ITensors
 using ITensors: data
 
 include("subnetwork.jl")
-include("models/models.jl")
 include("ITensors.jl")
 include("lattices.jl")
+include("models/models.jl")
 include("inds_network.jl")
 include("itensor_network.jl")
 include("boundary_mps.jl")

src/ITensorNetworks/chain_rules.jl

@@ -80,15 +80,15 @@ end
 
 # gradient of this function returns nothing.
 @non_differentiable generate_inner_network(
-  peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, Hlocal::Array
+  peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, Hs::Array
 )
 
 @non_differentiable generate_inner_network(
   peps::PEPS,
   peps_prime::PEPS,
   peps_prime_ham::PEPS,
   projectors::Array{<:ITensor,1},
-  Hlocal::Array,
+  Hs::Array,
 )
 
 @non_differentiable insert_projectors(peps::PEPS, center, cutoff, maxdim)

src/ITensorNetworks/lattices.jl

@@ -51,7 +51,6 @@ function neighbor(
   )
 end
 
-#function isboundary(site::dir::Int,
 # Check if the edge connecting to the specified neighbor of the site
 # crosses the boundary of the lattice
 function isboundary(
@@ -61,16 +60,33 @@ function isboundary(
 end
 
 # The neighboring sites of the specified site
-function neighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}) where {N}
+function filterneighbors(
+  f, lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false
+) where {N}
   lattice_size = size(lattice)
   site_neighbors = Vector{NTuple{N,Int}}()
   for dim in 1:N, dir in (-1, 1)
     site_neighbor = neighbor(site, dim, dir; lattice_size=lattice_size)
-    push!(site_neighbors, neighbor)
+    bc_condition = periodic || !(isboundary(site, dim, dir; lattice_size=lattice_size))
+    if f(site, site_neighbor) && bc_condition
+      push!(site_neighbors, site_neighbor)
+    end
   end
   return site_neighbors
 end
 
+function neighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false) where {N}
+  return filterneighbors(≠, lattice, site; periodic=periodic)
+end
+
+function inneighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false) where {N}
+  return filterneighbors(>, lattice, site; periodic=periodic)
+end
+
+function outneighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}; periodic=false) where {N}
+  return filterneighbors(<, lattice, site; periodic=periodic)
+end
+
 # All of the edges connected to the vertex `site`
 function incident_edges(lattice::HyperCubic{N}, site::NTuple{N,Int}) where {N}
   lattice_size = size(lattice)
@@ -90,6 +106,20 @@ function incident_edges(lattice::HyperCubic{N}, site::NTuple{N,Int}) where {N}
   return site_edges
 end
 
-function bonds(lattice::HyperCubic)
-  return ((s, n) for s in sites(lattice) for n in outneighbors(lattice, s))
+function bonds(lattice::HyperCubic; periodic=false)
+  return [
+    (s, n) for s in sites(lattice) for n in outneighbors(lattice, s; periodic=periodic)
+  ]
+end
+
+function bonds(lattice::Square, coord::Tuple{Colon,<:Integer})
+  rowsize = lattice.dims[1]
+  colsites = [(i, coord[2]) for i in 1:(rowsize - 1)]
+  return [(s, (s[1] + 1, s[2])) for s in colsites]
+end
+
+function bonds(lattice::Square, coord::Tuple{<:Integer,Colon})
+  colsize = lattice.dims[2]
+  rowsites = [(coord[1], i) for i in 1:(colsize - 1)]
+  return [(s, (s[1], s[2] + 1)) for s in rowsites]
 end

src/ITensorNetworks/models/hamiltonians.jl

@@ -1,9 +1,15 @@
 using ITensors
+include("../lattices.jl")
 
 struct LocalMPO
   mpo::MPO
-  coord1::Pair{<:Integer,<:Integer}
-  coord2::Pair{<:Integer,<:Integer}
+  coord1::Tuple{<:Integer,<:Integer}
+  coord2::Tuple{<:Integer,<:Integer}
+end
+
+struct LineMPO
+  mpo::MPO
+  coord::Union{Tuple{Colon,<:Integer},Tuple{<:Integer,Colon}}
 end
 
 # Transverse field
@@ -24,47 +30,80 @@ function mpo(::Model"tfim", sites::Matrix{<:Index}; h::Float64)
   return MPO(opsum, sites_vec)
 end
 
-# Get the local hamiltonian term of a 2D grid
-function localham_term(::Model"tfim", sites::Matrix{<:Index}, bond; h::Float64)
+function localham_term(
+  ::Model"tfim",
+  sites::Matrix{<:Index},
+  bond::Tuple{Tuple{<:Integer,<:Integer},Tuple{<:Integer,<:Integer}};
+  h::Float64,
+)
   Ny, Nx = size(sites)
-  sites_vec = vec(sites)
-  n1, n2 = bond.s1, bond.s2
+  coord1, coord2 = bond
   opsum = OpSum()
   opsum += -1, "X", 1, "X", 2
-  if n2 == n1 + 1
+  if coord2[1] == coord1[1] + 1
     opsum += h, "Z", 1
   end
-  if n2 == n1 + 1 && n2 % Ny == 0
+  if coord2[1] == coord1[1] + 1 && coord2[1] == Ny
     opsum += h, "Z", 2
   end
-  mpo = MPO(opsum, sites_vec[[n1, n2]])
-  coord1 = ((n1 - 1) % Ny + 1) => trunc(Int, (n1 - 1) / Ny) + 1
-  coord2 = ((n2 - 1) % Ny + 1) => trunc(Int, (n2 - 1) / Ny) + 1
+  mpo = MPO(opsum, [sites[coord1...], sites[coord2...]])
   return LocalMPO(mpo, coord1, coord2)
 end
 
-# Return a list of LocalMPO
 function localham(m::Model, sites; kwargs...)
   Ny, Nx = size(sites)
-  lattice = square_lattice(Nx, Ny; yperiodic=false)
-  return [localham_term(m, sites, bond; kwargs...) for bond in lattice]
+  lattice = Square((Ny, Nx))
+  bds = bonds(lattice; periodic=false)
+  return [localham_term(m, sites, bond; kwargs...) for bond in bds]
+end
+
+function localham_term(
+  ::Model"tfim", sites::Matrix{<:Index}, bond::Tuple{Colon,<:Integer}; h::Float64
+)
+  Ny, Nx = size(sites)
+  opsum = OpSum()
+  for i in 1:(Ny - 1)
+    opsum += -1, "X", i, "X", i + 1
+    opsum += h, "Z", i
+  end
+  opsum += h, "Z", Ny
+  return LineMPO(MPO(opsum, sites[bond...]), bond)
+end
+
+function localham_term(
+  ::Model"tfim", sites::Matrix{<:Index}, bond::Tuple{<:Integer,Colon}; h::Float64
+)
+  Ny, Nx = size(sites)
+  opsum = OpSum()
+  for i in 1:(Nx - 1)
+    opsum += -1, "X", i, "X", i + 1
+  end
+  return LineMPO(MPO(opsum, sites[bond...]), bond)
+end
+
+function lineham(m::Model, sites; kwargs...)
+  Ny, Nx = size(sites)
+  lattice = Square((Ny, Nx))
+  bonds_row = [(i, :) for i in 1:Ny]
+  bonds_column = [(:, i) for i in 1:Nx]
+  bds = vcat(bonds_row, bonds_column)
+  return [localham_term(m, sites, bond; kwargs...) for bond in bds]
 end
 
 # Check that the local Hamiltonian is the same as the MPO
-function checklocalham(Hlocal, H, sites)
+function checkham(Hlocal::Array{LocalMPO}, H, sites)
   @disable_warn_order begin
     Ny, Nx = size(sites)
-    sites_vec = reshape(sites, Nx * Ny)
-    lattice = square_lattice(Nx, Ny; yperiodic=false)
-
-    # This scales exponentially
+    lattice = Square((Ny, Nx))
+    bds = bonds(lattice; periodic=false)
     Hlocal_full = ITensor()
-    for (i, bond) in enumerate(lattice)
+    for (i, bond) in enumerate(bds)
       Hlocalterm_full = prod(Hlocal[i].mpo)
-      n1, n2 = bond.s1, bond.s2
-      for m in 1:(Nx * Ny)
-        if !(m in (n1, n2))
-          Hlocalterm_full *= op("Id", sites_vec, m)
+      for y in 1:Ny
+        for x in 1:Nx
+          if !((y, x) in bond)
+            Hlocalterm_full *= op("Id", vec(sites), (x - 1) * Ny + y)
+          end
         end
       end
       Hlocal_full += Hlocalterm_full
@@ -73,3 +112,33 @@ function checklocalham(Hlocal, H, sites)
   end
   return isapprox(norm(Hlocal_full), norm(prod(H)))
 end
+
+function checkham(Hline::Array{LineMPO}, H, sites)
+  @disable_warn_order begin
+    Ny, Nx = size(sites)
+    Hlocal_full = ITensor()
+    for h in Hline
+      h_full = prod(h.mpo)
+      if h.coord[1] isa Colon
+        for y in 1:Ny
+          for x in 1:Nx
+            if x != h.coord[2]
+              h_full *= op("Id", vec(sites), (x - 1) * Ny + y)
+            end
+          end
+        end
+      else
+        for y in 1:Ny
+          for x in 1:Nx
+            if y != h.coord[1]
+              h_full *= op("Id", vec(sites), (x - 1) * Ny + y)
+            end
+          end
+        end
+      end
+      Hlocal_full += h_full
+    end
+    @show norm(Hlocal_full - prod(H))
+  end
+  return isapprox(norm(Hlocal_full), norm(prod(H)))
+end

src/ITensorNetworks/peps.jl

@@ -1,4 +1,5 @@
 using Random
+using .Models
 
 """
 A finite size PEPS type.
@@ -126,8 +127,8 @@ function inner_network(
   dimy, dimx = size(peps.data)
   for ii in 1:dimx
     for jj in 1:dimy
-      if (jj => ii) in coordinates
-        index = findall(x -> x == (jj => ii), coordinates)
+      if (jj, ii) in coordinates
+        index = findall(x -> x == (jj, ii), coordinates)
         @assert(length(index) == 1)
         network = vcat(network, [mpo.data[index[1]]])
         network = vcat(network, [peps_prime_ham.data[jj, ii]])
@@ -160,19 +161,26 @@ Parameters
 peps: a peps network with datatype PEPS
 peps_prime: prime of peps used for inner products
 peps_prime_ham: prime of peps used for calculating expectation values
-Hlocal: An array of MPO operators with datatype LocalMPO
+Hs: An array of MPO operators with datatype LocalMPO
 Returns
 -------
 An array of networks.
 """
 function generate_inner_network(
-  peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, Hlocal::Array
+  peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, Hs::Array
 )
   network_list = Vector{Vector{ITensor}}()
-  for H_term in Hlocal
-    inner = inner_network(
-      peps, peps_prime, peps_prime_ham, H_term.mpo, [H_term.coord1, H_term.coord2]
-    )
+  for H_term in Hs
+    if H_term isa Models.LocalMPO
+      coords = [H_term.coord1, H_term.coord2]
+    elseif H_term isa Models.LineMPO
+      if H_term.coord[1] isa Colon
+        coords = [(i, H_term.coord[2]) for i in 1:length(H_term.mpo)]
+      else
+        coords = [(H_term.coord[1], i) for i in 1:length(H_term.mpo)]
+      end
+    end
+    inner = inner_network(peps, peps_prime, peps_prime_ham, H_term.mpo, coords)
     network_list = vcat(network_list, [inner])
   end
   inner = inner_network(peps, peps_prime)
@@ -185,9 +193,9 @@ function generate_inner_network(
   peps_prime::PEPS,
   peps_prime_ham::PEPS,
   projectors::Array{<:ITensor,1},
-  Hlocal::Array,
+  Hs::Array,
 )
-  network_list = generate_inner_network(peps, peps_prime, peps_prime_ham, Hlocal)
+  network_list = generate_inner_network(peps, peps_prime, peps_prime_ham, Hs)
   return map(network -> vcat(network, projectors), network_list)
 end