Dear All,
I have encountered an error whilst using the new expand function included in the ITensorTDVP package.
I am using a Hamiltonian which includes a boson-fermion coupling term:
H =-t_h \sum_{i=1}^{N} (e^{{i\frac{g}{\sqrt{L}}}(a + a^\dag)}c^\dag_jc_{j+1} + h.c.)
My system is modelled with one boson site and a chain of length N of spinless fermion sites. I am evolving the initial state in imaginary time to obtain a thermal state. I am using N physical sites and N ancillary sites to achieve this.
When I use this Hamiltonian with the coupling term (the exponential term) included I run into an error. If I alter the Hamiltonian and remove the coupling, the expand() function works fine for me.
Here is the code that I used to generate the Hamiltonian:
function H_kinetic(L, sites, N_max, t_h, g, mu)
op = OpSum()
for j in 3:(2*L) # fermion sites are from 3:2L+2
if j % 2 == 0 # if it's an even site ignore it (ancilla site)
else
op += -t_h, "Cdag", j, "C", j+2 # if it's a physical site add the term
end
end
op += -t_h, "Cdag", ((2*L)+1), "C", 3 # periodic boundary conditions
Kinetic = MPO(op, sites)
op_conj = OpSum()
for j in 3:(2*L)
if j % 2 == 0
else
op_conj += -t_h, "Cdag", j+2, "C", j
end
end
op_conj += -t_h, "Cdag", 3, "C", ((2*L)+1) # Periodic boundary conditions
Kinetic_Conjugate = MPO(op_conj, sites)
# now need to add boson-fermion coupling
a = boson_creation_op(N_max) # Create the boson creation operator
adag = a' # Hermitian Conjugate
exponent = 1im * g/sqrt(L) * (a + adag)
exponent_conj = -1im * g/sqrt(L) * (a + adag)
exp_matrix = exp(exponent)
exp_matrix_conj = exp(exponent_conj)
s = sites[1]
spp = prime(prime(s))
Boson_ITensor = ITensor(exp_matrix, s, spp)
Boson_ITensor_Conj = ITensor(exp_matrix_conj, s, spp)
Bosonic_Coupling = Kinetic[1] * Boson_ITensor
Kinetic[1] = Bosonic_Coupling
Bosonic_Coupling_Conj = Kinetic_Conjugate[1] * Boson_ITensor_Conj
Kinetic_Conjugate[1] = Bosonic_Coupling_Conj
Kinetic = mapprime(Kinetic, 2 => 0)
Kinetic_Conjugate = mapprime(Kinetic_Conjugate, 2 => 0)
return Kinetic, Kinetic_Conjugate
end
I then give this Hamiltonian to expand:
H1, H2 = H_kinetic(L, sites, N_max, t_h, g, mu)
ψ_0 = inital_thermal_state(sites, L, N_max, cutoff)
expanded_ψ = expand(ψ_0, H; alg="global_krylov", krylovdim=2)
Which then results in this error:
InexactError: Float64(2.000000000000001 - 4.3774855293491566e-24im)
Stacktrace:
[1] Float64(z::ComplexF64)
@ Base ./complex.jl:44
[2] expand(::NDTensors.BackendSelection.Algorithm{:orthogonalize, NamedTuple{(), Tuple{}}}, state::MPS, references::Vector{MPS}; cutoff::Float64)
@ ITensorTDVP ~/.julia/packages/ITensorTDVP/YqF4o/src/expand.jl:94
[3] expand
@ ~/.julia/packages/ITensorTDVP/YqF4o/src/expand.jl:73 [inlined]
[4] expand(state::MPS, reference::Vector{MPS}; alg::String, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:cutoff,), Tuple{Float64}}})
@ ITensorTDVP ~/.julia/packages/ITensorTDVP/YqF4o/src/expand.jl:29
[5] expand
@ ~/.julia/packages/ITensorTDVP/YqF4o/src/expand.jl:28 [inlined]
[6] expand(::NDTensors.BackendSelection.Algorithm{:global_krylov, NamedTuple{(), Tuple{}}}, state::MPS, operator::MPO; krylovdim::Int64, cutoff::Float64, apply_kwargs::NamedTuple{(:maxdim,), Tuple{Int64}})
@ ITensorTDVP ~/.julia/packages/ITensorTDVP/YqF4o/src/expand.jl:158
[7] expand(state::MPS, reference::MPO; alg::String, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:krylovdim,), Tuple{Int64}}})
@ ITensorTDVP ~/.julia/packages/ITensorTDVP/YqF4o/src/expand.jl:29
[8] top-level scope
@ ~/Desktop/Summer_Placement/Julia_DMRG/Quantum Variance/test.ipynb:171
I am hoping that there is an easy fix. If not I figured it was still worth reporting to you anyway.
I will put a full example code block below which reproduces this error:
using ITensors, ITensorMPS
# create identity matrix for given dimension (dimxdim)
function identity(dim)
Id = zeros(dim, dim)
for i in 1:dim
Id[i, i] = 1
end
return Id
end
# create matrix rep of boson creation operator
function boson_creation_op(N_max)
dim = N_max + 1
a = zeros(dim, dim)
for i in 1:dim
for j in 1:dim
if i == j-1
a[i, j] = sqrt(i)
else
a[i, j] = 0
end
end
end
return a
end
function H_kinetic(L, sites, N_max, t_h, g, mu)
op = OpSum()
for j in 3:(2*L) # fermion sites are from 3:2L+2
if j % 2 == 0 # if it's an even site ignore it (ancilla site)
else
op += -t_h, "Cdag", j, "C", j+2 # if it's a physical site add the term
end
end
op += -t_h, "Cdag", ((2*L)+1), "C", 3 # periodic boundary conditions
Kinetic = MPO(op, sites)
op_conj = OpSum()
for j in 3:(2*L)
if j % 2 == 0
else
op_conj += -t_h, "Cdag", j+2, "C", j
end
end
op_conj += -t_h, "Cdag", 3, "C", ((2*L)+1) # Periodic boundary conditions
Kinetic_Conjugate = MPO(op_conj, sites)
# now need to add boson-fermion coupling
a = boson_creation_op(N_max) # Create the boson creation operator
adag = a' # Hermitian Conjugate
exponent = 1im * g/sqrt(L) * (a + adag)
exponent_conj = -1im * g/sqrt(L) * (a + adag)
exp_matrix = exp(exponent)
exp_matrix_conj = exp(exponent_conj)
s = sites[1]
spp = prime(prime(s))
Boson_ITensor = ITensor(exp_matrix, s, spp)
Boson_ITensor_Conj = ITensor(exp_matrix_conj, s, spp)
Bosonic_Coupling = Kinetic[1] * Boson_ITensor
Kinetic[1] = Bosonic_Coupling
Bosonic_Coupling_Conj = Kinetic_Conjugate[1] * Boson_ITensor_Conj
Kinetic_Conjugate[1] = Bosonic_Coupling_Conj
Kinetic = mapprime(Kinetic, 2 => 0)
Kinetic_Conjugate = mapprime(Kinetic_Conjugate, 2 => 0)
return Kinetic, Kinetic_Conjugate
end
# create the sites for the chain + ancillas
function Ancillary_Chain(L, N_max)
sites = Vector{Index}(undef, ((2*L)+2)) # add on boson and ancillary boson sites
sites[1] = siteind("Boson", dim=(N_max+1)) # give the first site a boson index
sites[2] = siteind("Boson", dim=(N_max+1)) # give the ancillary site a boson index
for j in 3:((2*L)+2)
sites[j] = siteind("Fermion") # and give the rest of them fermion indices
end
return sites
end
# generate the initial, maximally entangled state for use in TDVP
function inital_thermal_state(sites, L, N_max, cutoff)
list = 1:(2L+2) # list of site numbers
cdag = ITensors.ops(sites, [(n == 1 || n == 2 ? ("Id", n) : ("Cdag", n)) for n in list]) # identity on boson sites, cdag on fermion sites
adag = ITensors.ops(sites, [(n == 1 || n == 2 ? ("Adag", n) : ("Id", n)) for n in list]) # adag on boson sites, identity on fermion sites
Bos_Id = identity(N_max+1) # boson identity matrix - identity function just constructs this
Ferm_Id = identity(2) # fermion identity matrix
Id = Vector{ITensor}(undef, (2L+2)) # +2 is for the bosonic site and ancilla. identity MPO
for i = 1:(2L+2)
if i == 1 || i == 2 # if it's a boson site
Id[i] = ITensor(Bos_Id, sites[i]', sites[i]) # use the boson identity
else # else if it's a fermion site
Id[i] = ITensor(Ferm_Id, sites[i]', sites[i]) # use fermion identity
end
end
vacstate = productMPS(sites, "0")
fermgates = [(cdag[n]*cdag[n+1] + Id[n]*Id[n+1])/sqrt(2) for n in 3:2:(2L+2)] # gates for fermion chain
ferm_maxentangled = apply(fermgates, vacstate; cutoff=1E-30) # fermions now in correct state
ops_array = [Id[1]*Id[2], (adag[1]*adag[2])] # Initialize ops_array with identity and |11>
loop1_array = [adag[1]]
loop2_array = [adag[2]]
for i in 2:N_max
adag1_copy1 = deepcopy(adag[1])
previous_op1 = deepcopy(loop1_array[i-1])
prime!(adag1_copy1, 2, plev=0)
prime!(previous_op1, 1, plev=1)
newop1 = adag1_copy1 * previous_op1
adag2_copy1 = deepcopy(adag[2])
previous_op2 = deepcopy(loop2_array[i-1])
prime!(adag2_copy1, 2, plev=0)
prime!(previous_op2, 1, plev=1)
newop2 = adag2_copy1 * previous_op2
finalop = (newop1 * newop2)/factorial(i)
push!(loop1_array, newop1)
push!(loop2_array, newop2)
push!(ops_array, finalop)
end
bosgates = [sum(ops_array)]
thermalstate = apply(bosgates, ferm_maxentangled; cutoff=1E-30)
return thermalstate
end
L = 6
N_max = 16
t_h = 1
g = 1
mu = 0
sites = Ancillary_Chain(L, N_max)
H1, H2 = H_kinetic(L, sites, N_max, t_h, g, mu)
H = +(H1, H2)
ψ_0 = inital_thermal_state(sites, L, N_max, 1E-20)
expanded_ψ = expand(ψ_0, H; alg="global_krylov", krylovdim=2)
I appreciate any and all advice you could offer.
Many thanks!