Rename tdvp to time_evolve. Add tests.

Author: emstoudenmire

src/solvers/applyexp.jl

@@ -7,16 +7,16 @@ import ConstructionBase: setproperties
   current_time::Number = 0.0
 end
 
-ITensorNetworks.state(tdvp::ApplyExpProblem) = tdvp.state
-operator(tdvp::ApplyExpProblem) = tdvp.operator
-current_time(tdvp::ApplyExpProblem) = tdvp.current_time
+ITensorNetworks.state(A::ApplyExpProblem) = A.state
+operator(A::ApplyExpProblem) = A.operator
+current_time(A::ApplyExpProblem) = A.current_time
 
 function region_plan(tdvp::ApplyExpProblem; nsites, time_step, sweep_kwargs...)
   return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...)
 end
 
 function updater(
-  T::ApplyExpProblem,
+  A::ApplyExpProblem,
   local_state,
   region_iterator;
   nsites,
@@ -25,26 +25,26 @@ function updater(
   outputlevel,
   kws...,
 )
-  local_state, info = solver(x->optimal_map(operator(T), x), time_step, local_state; kws...)
+  local_state, info = solver(x->optimal_map(operator(A), x), time_step, local_state; kws...)
 
   if nsites==1
     curr_reg = current_region(region_iterator)
     next_reg = next_region(region_iterator)
     if !isnothing(next_reg) && next_reg != curr_reg
-      next_edge = first(edge_sequence_between_regions(state(T), curr_reg, next_reg))
+      next_edge = first(edge_sequence_between_regions(state(A), curr_reg, next_reg))
       v1, v2 = src(next_edge), dst(next_edge)
-      psi = copy(state(T))
+      psi = copy(state(A))
       psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2]))
-      shifted_operator = position(operator(T), psi, NamedEdge(v1=>v2))
+      shifted_operator = position(operator(A), psi, NamedEdge(v1=>v2))
       R_t, _ = solver(x->optimal_map(shifted_operator, x), -time_step, R; kws...)
       local_state = psi[v1]*R_t
     end
   end
 
-  curr_time = current_time(T) + time_step
-  T = setproperties(T; current_time=curr_time)
+  curr_time = current_time(A) + time_step
+  A = setproperties(A; current_time=curr_time)
 
-  return T, local_state
+  return A, local_state
 end
 
 function applyexp_sweep_printer(
@@ -90,7 +90,7 @@ end
 
 process_real_times(z) = round(-imag(z); digits=10)
 
-function tdvp(
+function time_evolve(
   H,
   init_state,
   time_points;

test/solvers/test_tree_tdvp.jl

@@ -0,0 +1,76 @@
+using Test: @test, @testset
+using ITensors
+using TensorOperations # Needed to use contraction order finding
+import ITensorNetworks: dmrg, maxlinkdim, siteinds, time_evolve, ttn
+import Graphs: add_vertex!, add_edge!, vertices
+import NamedGraphs: NamedGraph
+import ITensorMPS: OpSum
+
+function chain_plus_ancilla(; nchain)
+  g = NamedGraph()
+  for j in 1:nchain
+    add_vertex!(g, j)
+  end
+  for j in 1:(nchain - 1)
+    add_edge!(g, j=>j+1)
+  end
+  # Add ancilla vertex near middle of chain
+  add_vertex!(g, 0)
+  add_edge!(g, 0=>nchain÷2)
+  return g
+end
+
+@testset "Tree TDVP on chain plus ancilla" begin
+  outputlevel = 1
+
+  N = 10
+  g = chain_plus_ancilla(; nchain=N)
+
+  sites = siteinds("S=1/2", g)
+
+  # Make Heisenberg model Hamiltonian
+  h = OpSum()
+  for j in 1:(N - 1)
+    h += "Sz", j, "Sz", j+1
+    h += 1/2, "S+", j, "S-", j+1
+    h += 1/2, "S-", j, "S+", j+1
+  end
+  H = ttn(h, sites)
+
+  # Make initial product state
+  state = Dict{Int,String}()
+  for (j, v) in enumerate(vertices(sites))
+    state[v] = iseven(j) ? "Up" : "Dn"
+  end
+  psi0 = ttn(state, sites)
+
+  cutoff = 1E-10
+  maxdim = 100
+  nsweeps = 5
+
+  nsites = 2
+  trunc = (; cutoff, maxdim)
+  E, gs_psi = dmrg(H, psi0; inserter_kwargs=(; trunc), nsites, nsweeps, outputlevel)
+  (outputlevel >= 1) && println("2-site DMRG energy = ", E)
+
+  inserter_kwargs=(; trunc)
+  nsites = 1
+  tmax = 0.10
+  time_range = 0.0:0.02:tmax
+  psi1_t = time_evolve(H, gs_psi, time_range; inserter_kwargs, nsites, outputlevel)
+  (outputlevel >= 1) && println("Done with $nsites-site TDVP")
+
+  @test norm(psi1_t) > 0.999
+
+  nsites = 2
+  psi2_t = time_evolve(H, gs_psi, time_range; inserter_kwargs, nsites, outputlevel)
+  (outputlevel >= 1) && println("Done with $nsites-site TDVP")
+  @test norm(psi2_t) > 0.999
+
+  @test abs(inner(psi1_t, gs_psi)) > 0.99
+  @test abs(inner(psi1_t, psi2_t)) > 0.99
+
+  # Test that accumulated phase angle is E*tmax
+  z = inner(psi1_t, gs_psi)
+  @test abs(atan(imag(z)/real(z)) - E*tmax) < 1E-4
+end