Hi,
I’m wondering if it’s possible to model an additional system coupled to the bulk Hamiltonian in ITensor. More concretely, let’s say the Hilbert space of our system takes the form
|\alpha \rangle \otimes |s_1s_2\cdots s_N\rangle
In the simplest case, s_i are fermion sites with QN conservation, and |\alpha\rangle is a two-level system coupled by some interaction hamiltonian
\langle \alpha s_1\cdots s_N| H_{int} | \alpha' s'_1\cdots s'_N \rangle \sim ( 1 - \delta_{\alpha, \alpha'}) f( s_1, s_1'\cdots )
And there is some diagonal bulk Hamiltonian only acting on the s_i sites. In ED this is easily written as
\begin{pmatrix}
H_{bulk} & H_{int} \\
H_{int} & H_{bulk}
\end{pmatrix}
My first implementation is to simply define the sites as
sites = Vector{Index}(undef, L + 1)
#two level system
sites[1] = siteind("Fermion"; addtags="n=1", conserve_qns =false)
for mid = 2: L + 1
sites[mid] = siteind("Fermion"; addtags="n=$mid", conserve_qns =true)
end
It seems even with the simplest bulk/interaction Hamiltonian, it could not be properly transformed into a MPO by calling H = MPO(ampo, sites), as I get the following errors
ERROR: MethodError: no method matching contraction_output(::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, ::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, ::Tuple{Index{Int64}, Index{Int64}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}})
Closest candidates are:
contraction_output(::NDTensors.Tensor, ::Any, ::NDTensors.Tensor, ::Any, ::Any) at ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:602
contraction_output(::NDTensors.Tensor, ::NDTensors.DiagBlockSparseTensor, ::Any) at ~/.julia/packages/NDTensors/c2BpJ/src/blocksparse/diagblocksparse.jl:290
contraction_output(::NDTensors.DiagBlockSparseTensor, ::NDTensors.Tensor, ::Any) at ~/.julia/packages/NDTensors/c2BpJ/src/blocksparse/diagblocksparse.jl:287
...
Stacktrace:
[1] contraction_output(T1::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT1::Tuple{Int64, Int64}, T2::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, labelsT2::Tuple{Int64, Int64}, labelsR::NTuple{4, Int64})
@ NDTensors ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:604
[2] contract(T1::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT1::Tuple{Int64, Int64}, T2::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, labelsT2::Tuple{Int64, Int64}, labelsR::NTuple{4, Int64})
@ NDTensors ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:620
[3] contract(T1::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT1::Tuple{Int64, Int64}, T2::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, labelsT2::Tuple{Int64, Int64})
@ NDTensors ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:611
[4] _contract(A::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, B::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}})
@ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1864
[5] _contract(A::ITensor, B::ITensor)
@ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1870
[6] contract(A::ITensor, B::ITensor)
@ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1974
[7] *
@ ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1961 [inlined]
[8] svdMPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ITensors ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo.jl:119
[9] svdMPO
@ ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo.jl:2 [inlined]
[10] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ITensors ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo_generic.jl:226
[11] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index})
@ ITensors ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo_generic.jl:217
As far as I know, QN conservation is a U(1) symmetry that enforces certain ‘flux’ structure on state tensors, so I’m wondering if this has to do with the specifics of the Hamiltonian or it’s just some general incompatibility.