Fix error in determining exponents. Change many variable names to more generic and correct exponent terminology.

Author: emstoudenmire

src/solvers/applyexp.jl

@@ -3,39 +3,46 @@ using Printf: @printf
 @kwdef mutable struct ApplyExpProblem{State}
   operator
   state::State
-  current_time::Number = 0.0
+  current_exponent::Number = 0.0
 end
 
 ITensorNetworks.state(A::ApplyExpProblem) = A.state
 operator(A::ApplyExpProblem) = A.operator
-current_time(A::ApplyExpProblem) = A.current_time
+current_exponent(A::ApplyExpProblem) = A.current_exponent
+function current_time(A::ApplyExpProblem)
+  t = im*A.current_exponent
+  return iszero(imag(t)) ? real(t) : t
+end
 
 function set_operator(A::ApplyExpProblem, operator)
-  ApplyExpProblem(operator, A.state, A.current_time)
+  ApplyExpProblem(operator, A.state, A.current_exponent)
+end
+function set_state(A::ApplyExpProblem, state)
+  ApplyExpProblem(A.operator, state, A.current_exponent)
 end
-set_state(A::ApplyExpProblem, state) = ApplyExpProblem(A.operator, state, A.current_time)
-function set_current_time(A::ApplyExpProblem, current_time)
-  ApplyExpProblem(A.operator, A.state, current_time)
+function set_current_exponent(A::ApplyExpProblem, current_exponent)
+  ApplyExpProblem(A.operator, A.state, current_exponent)
 end
 
 function region_plan(tdvp::ApplyExpProblem; nsites, time_step, sweep_kwargs...)
-  return tdvp_regions(state(tdvp), time_step; nsites, sweep_kwargs...)
+  return applyexp_regions(state(tdvp), time_step; nsites, sweep_kwargs...)
 end
 
 function update(
   prob::ApplyExpProblem,
   local_state,
   region_iterator;
   nsites,
-  time_step,
+  exponent_step,
   solver=runge_kutta_solver,
   outputlevel,
   kws...,
 )
+  iszero(abs(exponent_step)) && return prob, local_state
+
   local_state, info = solver(
-    x->optimal_map(operator(prob), x), time_step, local_state; kws...
+    x->optimal_map(operator(prob), x), exponent_step, local_state; kws...
   )
-
   if nsites==1
     curr_reg = current_region(region_iterator)
     next_reg = next_region(region_iterator)
@@ -45,21 +52,29 @@ function update(
       psi = copy(state(prob))
       psi[v1], R = qr(local_state, uniqueinds(local_state, psi[v2]))
       shifted_operator = position(operator(prob), psi, NamedEdge(v1=>v2))
-      R_t, _ = solver(x->optimal_map(shifted_operator, x), -time_step, R; kws...)
+      R_t, _ = solver(x->optimal_map(shifted_operator, x), -exponent_step, R; kws...)
       local_state = psi[v1]*R_t
     end
   end
 
-  prob = set_current_time(prob, current_time(prob)+time_step)
+  prob = set_current_exponent(prob, current_exponent(prob)+exponent_step)
 
   return prob, local_state
 end
 
 function sweep_printer(
-  problem::ApplyExpProblem; outputlevel, sweep, nsweeps, process_time=identity, kws...
+  problem::ApplyExpProblem;
+  exponent_description="exponent",
+  outputlevel,
+  sweep,
+  nsweeps,
+  process_time=identity,
+  kws...,
 )
   if outputlevel >= 1
-    @printf("  Current time = %s, ", process_time(current_time(problem)))
+    @printf(
+      "  Current %s = %s, ", exponent_description, process_time(current_exponent(problem))
+    )
     @printf("maxlinkdim=%d", maxlinkdim(state(problem)))
     println()
     flush(stdout)
@@ -75,30 +90,30 @@ function applyexp(
   insert_kwargs=(;),
   outputlevel=0,
   nsites=1,
-  tdvp_order=4,
+  order=4,
   kws...,
 )
   init_prob = ApplyExpProblem(;
     state=align_indices(init_state), operator=ProjTTN(align_indices(operator))
   )
-  time_steps = diff([zero(eltype(exponents)); exponents])[2:end]
-  sweep_kws = (;
-    outputlevel, extract_kwargs, insert_kwargs, nsites, tdvp_order, update_kwargs
-  )
-  kws_array = [(; sweep_kws..., time_step=t) for t in time_steps]
+  exponent_steps = diff([zero(eltype(exponents)); exponents])
+  sweep_kws = (; outputlevel, extract_kwargs, insert_kwargs, nsites, order, update_kwargs)
+  kws_array = [(; sweep_kws..., time_step=t) for t in exponent_steps]
   sweep_iter = sweep_iterator(init_prob, kws_array)
   converged_prob = sweep_solve(sweep_iter; outputlevel, kws...)
   return state(converged_prob)
 end
 
-process_real_times(z) = round(-imag(z); digits=10)
+process_real_times(z) = iszero(abs(z)) ? 0.0 : round(-imag(z); digits=10)
 
 function time_evolve(
   operator,
   time_points,
   init_state;
   process_time=process_real_times,
-  sweep_printer=(a...; k...)->sweep_printer(a...; process_time, k...),
+  sweep_printer=(
+    a...; k...
+  )->sweep_printer(a...; exponent_description="time", process_time, k...),
   kws...,
 )
   exponents = [-im*t for t in time_points]

src/solvers/region_plans/tdvp_region_plans.jl

@@ -1,26 +1,26 @@
-function tdvp_sub_time_steps(tdvp_order)
-  if tdvp_order == 1
+function applyexp_sub_steps(order)
+  if order == 1
     return [1.0]
-  elseif tdvp_order == 2
+  elseif order == 2
     return [1 / 2, 1 / 2]
-  elseif tdvp_order == 4
+  elseif order == 4
     s = (2 - 2^(1 / 3))^(-1)
     return [s/2, s/2, 1/2 - s, 1/2 - s, s/2, s/2]
   else
-    error("TDVP order of $tdvp_order not supported")
+    error("Applyexp order of $order not supported")
   end
 end
 
 function first_order_sweep(
-  graph, time_step, dir=Base.Forward; update_kwargs, nsites, kws...
+  graph, exponent_step, dir=Base.Forward; update_kwargs, nsites, kws...
 )
   basic_fwd_sweep = post_order_dfs_plan(graph; nsites, kws...)
-  update_kwargs = (; nsites, time_step, update_kwargs...)
+  update_kwargs = (; nsites, exponent_step, update_kwargs...)
   sweep = []
   for (j, (region, region_kws)) in enumerate(basic_fwd_sweep)
     push!(sweep, (region, (; nsites, update_kwargs, region_kws...)))
     if length(region) == 2 && j < length(basic_fwd_sweep)
-      rev_kwargs = (; update_kwargs..., time_step=(-update_kwargs.time_step))
+      rev_kwargs = (; update_kwargs..., exponent_step=(-update_kwargs.exponent_step))
       push!(sweep, ([last(region)], (; update_kwargs=rev_kwargs, region_kws...)))
     end
   end
@@ -31,13 +31,13 @@ function first_order_sweep(
   return sweep
 end
 
-function tdvp_regions(graph, time_step; update_kwargs, tdvp_order, nsites, kws...)
+function applyexp_regions(graph, exponent_step; update_kwargs, order, nsites, kws...)
   sweep_plan = []
-  for (step, weight) in enumerate(tdvp_sub_time_steps(tdvp_order))
+  for (step, weight) in enumerate(applyexp_sub_steps(order))
     dir = isodd(step) ? Base.Forward : Base.Reverse
     append!(
       sweep_plan,
-      first_order_sweep(graph, weight*time_step, dir; update_kwargs, nsites, kws...),
+      first_order_sweep(graph, weight*exponent_step, dir; update_kwargs, nsites, kws...),
     )
   end
   return sweep_plan

test/solvers/test_applyexp.jl

@@ -1,8 +1,9 @@
 using Test: @test, @testset
 using ITensors
-using ITensorNetworks: dmrg, maxlinkdim, siteinds, time_evolve, ttn
+using ITensorNetworks: ITensorNetworks, applyexp, dmrg, maxlinkdim, siteinds, time_evolve, ttn
 using Graphs: add_vertex!, add_edge!, vertices
 using NamedGraphs: NamedGraph
+using NamedGraphs.NamedGraphGenerators: named_path_graph
 using ITensorMPS: OpSum
 using TensorOperations: TensorOperations #For contraction order finding
 
@@ -20,7 +21,7 @@ function chain_plus_ancilla(; nchain)
   return g
 end
 
-@testset "Test Applyexp" begin
+@testset "Test Tree Time Evolution" begin
   outputlevel = 0
 
   N = 10
@@ -74,3 +75,48 @@ end
   z = inner(psi1_t, gs_psi)
   @test abs(atan(imag(z)/real(z)) - E*tmax) < 1E-4
 end
+
+@testset "Applyexp Time Point Handling" begin
+  N = 10
+  g = named_path_graph(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)
+
+  # 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)
+
+  nsites = 2
+  trunc = (; cutoff=1E-8, maxdim=100)
+  insert_kwargs=(; trunc)
+
+  # Test that all time points are reached and reported correctly
+  time_points = [0.0,0.1,0.25,0.32,0.4]
+  times = Real[]
+  function collect_times(problem; kws...)
+    push!(times, ITensorNetworks.current_time(problem))
+  end
+  time_evolve(H, time_points, psi0; insert_kwargs, nsites, sweep_callback=collect_times,outputlevel=1)
+  @test norm(times - time_points) < 10*eps(Float64)
+
+  # Test that all exponents are reached and reported correctly
+  exponent_points = [-0.0,-0.1,-0.25,-0.32,-0.4]
+  exponents = Real[]
+  function collect_exponents(problem; kws...)
+    push!(exponents, ITensorNetworks.current_exponent(problem))
+  end
+  applyexp(H, exponent_points, psi0; insert_kwargs, nsites, sweep_callback=collect_exponents,outputlevel=1)
+  @test norm(exponents - exponent_points) < 10*eps(Float64)
+end
+