Dynamical correlation function for thermal states

ITensor Julia Questions
1

Hi, I am computing dynamical correlation function using TEBD of the type \langle S^z(t)_i S^z(0)_j \rangle where the average \langle \rangle is with respect to thermal state (computed using ancilla based expansion algorithm) at inverse temperature \beta i.e. we have \rho(\beta). The notations (i,j) above are spin sites and t denotes at a particular time where S^z_i is evolved in the Heisenberg picture. Since I need to track the said correlation at long times, I am using fourth-order Trotter expansion, so that the error doesnt accummulate a lot. I find that the code is exceptionally slow for spin chains (S=1/2) even for N >= 10 sites.

Is there anything inbuilt in ITensor that I am not using properly, or even missing altogether that can fasten the code? Any help, will be greatly appreciated. I am attaching my code block below and the resources I used

using ITensors

function time_evolution_rho_based_driver(time, no_spins, sites, coeff, del_t, rho0, probe_index, order, cut_off=1E-8)
  
  function trotter_gates_construct(lambda=-1.0, order="second-order", bond_type="continuous")   
      
       start_index = 1
       steplength = 1
       for site_index in start_index:steplength:(no_spins-1)  
            push!(index_array, (site_index, site_index+1))
       end 
    
       gates = ITensor[]  
       for item in index_array
           s1 = sites[item[1]]
           s2 = sites[item[2]]
    
           #println(s1, s2)
       if order == "first-order"
             Gj = op("expτH", s1, s2, tau=-1.0*del_t * im, coeff=coeff) 
       elseif order == "second-order"   
             Gj = op("expτH", s1, s2, tau=-1.0*del_t * im / 2, coeff=coeff)
       else order == "fourth-order" 
             Gj = op("expτH", s1, s2, tau=lambda*del_t * im / 2, coeff=coeff)   
       end        
       push!(gates, Gj)
       end  
      # Include gates in reverse order too
      # (N,N-1),(N-1,N-2),...
      if order in ["second-order", "fourth-order"]        
         append!(gates, reverse(gates));  # This is for second order Trotter
      end
      return gates 
  end  
  
  # Compute the number of steps to do    
  Nsteps = Int(time / del_t) 
    
  # Compute and print initial <Sz> value
  t = 0.0

  # Do the time evolution by applying the gates
  # for Nsteps steps
  rho = rho0
  cut_off=1E-4
  rho = apply(one_body_op(sites, 1, 1,  "Sz"), rho; cut_off)  
  
  if order in ["second-order", "first-order"]  
      gate_array = trotter_gates_construct(-1.0, order, "continuous")   # this is second order , for first and second order lambda=-1.0
      for step in 1:Nsteps
         #println(step)
          rho = apply(gate_array, rho; cut_off, apply_dag=true)
          #println(psi)
          #t += del_t
      end
      #@show inner(psi', one_body_op(sites, 1, probe_index,  "Sz"), psi)     
      #psi = apply(op("Sz", sites, probe_index), psi; cut_off)
      @show abs(inner(one_body_op(sites, 1, probe_index,  "Sz"), rho))
        
    
  elseif order in ["fourth-order"]
      
      factor_1 = 1.0/(4-4^(1/3))
      factor_2 = 1-4*factor_1  
      gate_array_fact_1 = trotter_gates_construct(-factor_1, order, "continuous") #lambda=-1.0*factor
      gate_array_fact_2 = trotter_gates_construct(-factor_2, order, "continuous")  
      
      for step in 1:2*Nsteps
        rho = apply(gate_array_fact_1, rho; cut_off, apply_dag=true)
        #println(psi)
        #t += del_t
      end
        
      #t=0  
      for step in 1:Nsteps
        rho = apply(gate_array_fact_2, rho; cut_off, apply_dag=true)
        #println(psi)
        #t += del_t
      end
      
      for step in 1:2*Nsteps
        rho = apply(gate_array_fact_1, rho; cut_off, apply_dag=true)
        #println(psi)
        #t += del_t
      end 
      @show abs(inner(one_body_op(sites, 1, probe_index,  "Sz"), rho))  
        
  end  
end

The resources which I am using are the following

  1. For computing thermal states I am using- https://github.com/ITensor/ITensors.jl/blob/main/examples/finite_temperature/purification.jl
    Thermal states using this is computed correctly. I have verified that against exact diagonalization for small sizes

  2. For implementing TEBD
    https://github.com/ITensor/ITensors.jl/blob/main/examples/gate_evolution/mpo_gate_evolution.jl

I cannot use wave-function based MPS formulation (the reason being I am not really at low T limit, where the thermal state is the ground state, neither am I at the infinite T limit where the state is maximally mixed, I am somewhere in between and the probability distribution for thermal state has support over many eigenstates of H , so computing dynamical correlation using MPS would mean accessing many eigenvectors of H.

I feel I am probably missing something so any help to fasten the code, will be helpful, right now it takes hours even for N=10 spins.

2

Hi Manas,
I think the problem is that you are passing a keyword named cut_off to the apply function, while the keyword that function accepts is spelled cutoff without any underscore. If I’m correct about that issue, a good way to catch an issue like that is to print out or monitor the bond dimension of the tensor network state as you are evolving it. In this case I suspect that it is growing exponentially quickly with the number of time steps because no truncation is occurring.

3

@miles this is a good motivating use case for why we should fix up our keyword argument system so it catches misspelled keyword arguments. In this case, it may be sufficient to list all of the keyword arguments explicitly in a low level function like svd.

1 Like
4

Good point, Matt!

5

Thank you so much Miles and Matthew. That was indeed the issue. I kept overlooking those lines while debugging all the time thinking them to be de-facto correct and suspecting the issue to be elsewhere. While much faster than before, it is slightly slower than I would like it to be for N>=20 but that is probably the fourth-order Trotterization that I am using coupled with thermal state MPO.