QN Conservation by Site Property vs Added Operator

DMRG and Numerical Methods
1

I’m still pretty inexperienced with DMRG, so wanted to share my current experience with this issue and see what more experienced people think.

I’m investigating a system where charge mod 4 is conserved, and have been trying to figure out how to find the ground state/energy in each symmetry sector.

The first way I attempted was by defining an operator to have large energy penalties for the sectors I wanted to omit, which I then add to the regular Hamiltonian. The runtime per sweep is quite fast, and the bond dimension never gets beyond a few hundred even though I set the max dimension to 10K. The energy also seems to converge fairly well to a value, typically decreasing monotonically, and I consistently verify the resulting state is in the sector I want. I’m still working on ensuring it’s an energy eigenstate by looking at the variance. As for whether it’s the lowest energy or not is another concern which I’ll worry about later.

My second way which I’ve now been trying as I read more into the documentation is defining custom sites on top of the built-in “Electron” site class, with a charge mod 4 quantum number. I then define a random MPS from a state with well defined charge mod 4. I’ve been fiddling around with noise a bit alternating between low (10E-10) and zero noise by the end, and have rather low cutoff (10E-10). Despite that, for the same systems as I tested my first method on, the bond dimension seems to get into the 1000s at times, and the energies often fluctuate a lot, going up and down instead of converging nicely, and so on.

Is there any ideas what might be going on with my second approach of built-in conservation? And is there a general consensus for which of these two approaches is generally better?

That wasn’t the most coherent question, so I’d be happy to elaborate in words or provide code snippets of mine if that’d be helpful. Thank you!

2

It sounds like there’s a chance your Hamiltonian (as coded into ITensor I mean) isn’t Hermitian in the second approach, because an accidentally non-Hermitian H can lead to the effects you describe. Could you post the piece of your code where you define H?

3

Based on your other question about QN’s, another possibility is that you are defining your Hamiltonian correctly but not registering certain operators as fermionic for your custom site type. Please let us know a bit more about how you are making your site indices (e.g. what tags are you giving them?) and we can see if your setup is likely to give correct results.

4

Thanks so much for the quick replies on this and my other post. Here’s my Hamiltonian first:

H = OpSum()

for i=1:N-1
H += -1, “Cdagup”, i, “Cup”, i+1
H += -1, “Cdagdn”, i, “Cdn”, i+1
H += -1, “Cdagup”, i+1, “Cup”, i
H += -1, “Cdagdn”, i+1, “Cdn”, i
end

for i=1:N
H += -mu, “Ntot”, i
H += -U/2, “Ntot”, i
end

for i=1:N
H += U, “Nup”, i, “Ndn”, i (Planning to shorten this to “Nupdn”, but shouldn’t change things)
end

for i=1:N-1
H += delta, “Cdagup”, i, “Cdagdn”, i, “Cdagup”, i+1, “Cdagdn”, i+1
H += delta, “Cdn”, i+1, “Cup”, i+1, “Cdn”, i, “Cup”, i
end

H = MPO(H, sites)

5

Here’s some more code snippets regarding my custom sites to try to enforce charge mod 4 conservation:

function ITensors.space(::SiteType"E4"; conserve_q4=false)
if conserve_q4
return [
QN(“Q4”, 0, -4) => 1,
QN(“Q4”, 1, -4) => 2,
QN(“Q4”, 2, -4) => 1
]
end
return 4
end

And in my code, I’m setting conserve_q4 = true. Then for the operators, I pretty just much just paste everything directly from the electron site class like so:

function ITensors.op(::OpName"Nupdn", ::SiteType"E4")
return [
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0
]
end

function ITensors.op(::OpName"Ntot", ::SiteType"E4")
return [
0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0
0.0 0.0 1.0 0.0
0.0 0.0 0.0 2.0
]
end

function ITensors.op(::OpName"Cup", ::SiteType"E4")
return [
0.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0
]
end

function ITensors.op(::OpName"Cdagdn", ::SiteType"E4")
return [
0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0
0.0 -1.0 0.0 0.0
]
end

After all the operator definitions, I also ensure that my creation/annihilation operators have fermionic statistics with the following lines:

ITensors.has_fermion_string(::OpName"Cup", ::SiteType"E4") = true
ITensors.has_fermion_string(::OpName"Cdagup", ::SiteType"E4") = true
ITensors.has_fermion_string(::OpName"Cdn", ::SiteType"E4") = true
ITensors.has_fermion_string(::OpName"Cdagdn", ::SiteType"E4") = true

The site states are just from the electron site type too:

ITensors.val(::ValName"Emp", ::SiteType"E4") = 1
ITensors.val(::ValName"Up", ::SiteType"E4") = 2
ITensors.val(::ValName"Dn", ::SiteType"E4") = 3
ITensors.val(::ValName"UpDn", ::SiteType"E4") = 4

ITensors.state(::StateName"Emp", ::SiteType"E4") = [1.0, 0, 0, 0]
ITensors.state(::StateName"Up", ::SiteType"E4") = [0.0, 1, 0, 0]
ITensors.state(::StateName"Dn", ::SiteType"E4") = [0.0, 0, 1, 0]
ITensors.state(::StateName"UpDn", ::SiteType"E4") = [0.0, 0, 0, 1]

6

I feel like I made my Hamiltonian to be Hermitian, and with fermionic operators, but maybe I’m wrong about one of those claims.

I still have my second approach if somehow this doesn’t pan out, but even those states often have order 1 variance which isn’t the best. However, I’d imagine this approach to be the most robust if I can get it to work.

One of the few times when this approach did work, it hit what seemed to be the ground state and stayed there for the remaining sweeps, ending with variance on the order of 10E-7 which was awesome. That was for charge mod 4 of 0, but then for the same Hamiltonian parameters, other charge sectors fluctuate like so (maybe this will give some more insight):

After sweep 47 energy=-3.2530958041802336 maxlinkdim=724 maxerr=9.99E-11 time=38.223
After sweep 48 energy=-3.337560639050551 maxlinkdim=736 maxerr=1.00E-10 time=41.093
After sweep 49 energy=-4.129265930762855 maxlinkdim=737 maxerr=1.00E-10 time=37.960
After sweep 50 energy=-3.3465975257097753 maxlinkdim=740 maxerr=1.00E-10 time=36.979
After sweep 51 energy=-3.3648561468674174 maxlinkdim=741 maxerr=1.00E-10 time=38.588

Or like the following:

After sweep 194 energy=-2.4222737145329774 maxlinkdim=496 maxerr=1.00E-10 time=14.694
After sweep 195 energy=-2.521500237317497 maxlinkdim=499 maxerr=9.99E-11 time=14.477
After sweep 196 energy=-2.4635004631312967 maxlinkdim=473 maxerr=9.97E-11 time=14.468
After sweep 197 energy=-2.5302800029978982 maxlinkdim=504 maxerr=9.99E-11 time=13.877
After sweep 198 energy=-2.4391579528858642 maxlinkdim=447 maxerr=9.98E-11 time=13.886
After sweep 199 energy=-2.502675046675442 maxlinkdim=493 maxerr=9.99E-11 time=14.450

7

This all looks reasonable. Thanks for sharing these helpful details (and also the right details).

One question I have is, did you explicitly define the “F” operator? I think you may need to do that for a new fermionic site type like yours. I think the definition you’ll need would be:

function ITensors.op(::OpName"F", ::SiteType"E4")
  return [
    +1.0 0.0 0.0 0.0
    0.0 -1.0 0.0 0.0
    0.0 0.0  +1.0 0.0
    0.0 0.0 0.0  -1.0
  ]
end
8

I should add the energy oscillating up and down like that in DMRG definitely “feels” like a kind of bug in the Hamiltonian. But from the code you’ve showed it may be a pretty subtle one, quite possibly related to how fermions get handled in our OpSum system. But I think you are quite close to having it all working.

9

Thanks so much! That seems to have fixed my issue. Initially I didn’t include those fermion operators because I didn’t think they were needed, but now I see why they are.

10

Glad it’s working now!