Multiorbital electron mixed Hamiltonian

DMRG and Numerical Methods
1

Hi, thanks for the development of this library. I’ve been trying to create a set up similar to the one described in this paper (Phys. Rev. Lett. 126, 197202 (2021) - Quantum Spin Torque Driven Transmutation of an Antiferromagnetic Mott Insulator) were, in a 1D chain, all sites are of the ‘Electron’ type, except for a central region where on top of the electron sites, the authors add another orbital for electrons to interact.
I’ve been trying to use a mixed site setup as described in other examples. Thus, I’ve created a custom sites for the 2 electron orbital central region. I also want to conserve the total number of fermions. When I try to run DMRG on the constructed Hamiltonian I get the following error: Block Block(1, 1, 4) has flux QN((“Nf”,1,-1),(“Nt”,-1,-1)) that is inconsistent with the desired flux QN(“Nf”,0,-1).
The error appears when I consider the hopping elements between the central region and the two other single electron regions.
Here is a minimal code:

First I define the new site type for the central region including the symmetry of total number of fermions.

using ITensors

ITensors.enable_debug_checks()

function ITensors.space(::SiteType"2E";
                        conserve_qns=false, conserve_sz=false)
    if conserve_qns
        return [QN("Nt", 0, -1) => 1, QN("Nt", 1, -1) => 4,
                QN("Nt", 2, -1) => 6, QN("Nt", 3, -1) => 4,
                QN("Nt", 4, -1) => 1]
    end
    return 16
end

Then I label the 16 states

ITensors.val(::ValName"AEmpBEmp", ::SiteType"2E") = 1
ITensors.val(::ValName"AupBEmp", ::SiteType"2E") = 2
ITensors.val(::ValName"AdnBEmp", ::SiteType"2E") = 3
ITensors.val(::ValName"AEmpBup", ::SiteType"2E") = 4
ITensors.val(::ValName"AEmpBdn", ::SiteType"2E") = 5
ITensors.val(::ValName"AupdnBEmp", ::SiteType"2E") = 6
ITensors.val(::ValName"AEmpBupdn", ::SiteType"2E") = 7
ITensors.val(::ValName"AupBup", ::SiteType"2E") = 8
ITensors.val(::ValName"AdnBdn", ::SiteType"2E") = 9
ITensors.val(::ValName"AupBdn", ::SiteType"2E") = 10
ITensors.val(::ValName"AdnBup", ::SiteType"2E") = 11
ITensors.val(::ValName"AupdnBup", ::SiteType"2E") = 12
ITensors.val(::ValName"AupdnBdn", ::SiteType"2E") = 13
ITensors.val(::ValName"AupBupdn", ::SiteType"2E") = 14
ITensors.val(::ValName"AdnBupdn", ::SiteType"2E") = 15
ITensors.val(::ValName"AupdnBupdn", ::SiteType"2E") = 16

ITensors.state(::StateName"AEmpBEmp", ::SiteType"2E") =   [1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AupBEmp", ::SiteType"2E") =    [0, 1.0, 0, 0, 0, 0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AdnBEmp", ::SiteType"2E") =    [0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AEmpBup", ::SiteType"2E") =    [0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AEmpBdn", ::SiteType"2E") =    [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AupdnBEmp", ::SiteType"2E") =  [0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AEmpBupdn", ::SiteType"2E") =  [0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0 ,0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AupBup", ::SiteType"2E") =     [0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AupBdn", ::SiteType"2E") =     [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AdnBup", ::SiteType"2E") =     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0]
ITensors.state(::StateName"AupdnBup", ::SiteType"2E") =   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0, 0]
ITensors.state(::StateName"AupdnBdn", ::SiteType"2E") =   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0, 0]
ITensors.state(::StateName"AupBupdn", ::SiteType"2E") =   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0]
ITensors.state(::StateName"AdnBupdn", ::SiteType"2E") =   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0, 0]
ITensors.state(::StateName"AupdnBupdn", ::SiteType"2E") = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]

Then I create some basic operators. For the example I’m including only the hopping operators of the first orbital, which would be the orbital later connected to the single electron regions. I’m also adding the F fermonic operator

function ITensors.op!(Op::ITensor, ::OpName"Nt", ::SiteType"2E", s::Index)
    Op[s' => 2, s => 2] = 1
    Op[s' => 3, s => 3] = 1
    Op[s' => 4, s => 4] = 1
    Op[s' => 5, s => 5] = 1
    Op[s' => 6, s => 6] = 2
    Op[s' => 7, s => 7] = 2
    Op[s' => 8, s => 8] = 2
    Op[s' => 9, s => 9] = 2
    Op[s' => 10, s => 10] = 2
    Op[s' => 11, s => 11] = 2
    Op[s' => 12, s => 12] = 3
    Op[s' => 13, s => 13] = 3
    Op[s' => 14, s => 14] = 3
    Op[s' => 15, s => 15] = 3
    return Op[s' => 16, s => 16] = 4
end

function ITensors.op!(Op::ITensor, ::OpName"Caup", ::SiteType"2E", s::Index)
    Op[s' => 1, s => 2] = 1
    Op[s' => 3, s => 6] = -1
    Op[s' => 4, s => 8] = -1
    Op[s' => 5, s => 10] = -1
    Op[s' => 7, s => 14] = 1
    Op[s' => 9, s => 13] = 1
    Op[s' => 11, s => 12] = 1
    return Op[s' => 15, s => 16] = -1
end

function ITensors.op!(Op::ITensor, ::OpName"Cadagup", ::SiteType"2E", s::Index)
    Op[s' => 2, s => 1] = 1
    Op[s' => 6, s => 3] = -1
    Op[s' => 8, s => 4] = -1
    Op[s' => 10, s => 5] = -1
    Op[s' => 14, s => 7] = 1
    Op[s' => 13, s => 9] = 1
    Op[s' => 12, s => 11] = 1
    return Op[s' => 16, s => 15] = -1
end

function ITensors.op!(Op::ITensor, ::OpName"F", ::SiteType"2E", s::Index)
    Op[s' => 1, s => 1] = 1
    Op[s' => 2, s => 2] = -1
    Op[s' => 3, s => 3] = -1
    Op[s' => 4, s => 4] = -1
    Op[s' => 5, s => 5] = -1
    Op[s' => 6, s => 6] = 1
    Op[s' => 7, s => 7] = 1
    Op[s' => 8, s => 8] = 1
    Op[s' => 9, s => 9] = 1
    Op[s' => 10, s => 10] = 1
    Op[s' => 11, s => 11] = 1
    Op[s' => 12, s => 12] = -1
    Op[s' => 13, s => 13] = -1
    Op[s' => 14, s => 14] = -1
    Op[s' => 15, s => 15] = -1
    return Op[s' => 16, s => 16] = 1
end

I make sure to include the has_fermion_string for the necessary operators

ITensors.has_fermion_string(::OpName"Cadagup", ::SiteType"2E") = true
ITensors.has_fermion_string(::OpName"Caup", ::SiteType"2E") = true

Then I construct the Hamiltonian (for the example I only include one ‘Electron’ site to each side of the central ‘2E’ site)

N = 3
sites = siteinds(n -> n!=2 ? "Electron" : "2E", N; conserve_qns=true, conserve_sz=false)
println(sites)
os = OpSum()
os += -1, "Cdagup", 1, "Caup",  2
os += -1, "Cadagup", 2,  "Cup", 1
os += -1, "Cadagup", 2, "Cup",  1
os += -1, "Cdagup", 1,  "Caup", 2

H = MPO(os, sites)
state = ["Up", "AupBup", "Emp"]
ψ = randomMPS(sites, state, 2)
sweeps = Sweeps(30)
maxdim!(sweeps, 20, 20, 30, 50, 100, 100, 100, 200, 200, 200, 200, 200, 400)
noise!(sweeps, 1E-6, 1E-7, 1E-9)
 cutoff!(sweeps, 1E-9)
E, ϕ = dmrg(H, ψ, sweeps)

Thank you.

2

Hi, thanks for your patience with our slow reply. It looks like you are using the library correctly and way to go about delving into some of the more advanced features. My best guess about the error you are getting is that your Hamiltonian may not actually conserve the set of quantum numbers you are wishing to either due to an error in typing it in, or a slight conceptual error somewhere.

One thing you can do is call the functions checkflux and flux on your initial wavefunction to check that it has a well-defined flux and it is the one you want.

Another thing to definitely do is to put

ITensor.enable_debug_checks()

at the top of your code which will turn on additional, possibly more informative error messages.

If none of these reveal the problem, please send me an email with a short, fully working example code that reproduces the issue and I can probably find it for you.

3

Hello, Miles. Thank you for the answer. Yes, the problem was spotted when I used the flux option. The error was very simple, namely I had named the QN for the new site different than the one for the electron site and thus two QN were defined and the Hamiltonian did not conserve them.

4

Glad to hear you found the issue!