std::bad_alloc when performing dmrg

I am considering an open quantum system which is described by a non-Hermitian Hamiltonian. The Hamiltonian consists of two bosonic chains coupled to two custom sites (for which I implemented my own class with proper op method). In order to find a steady state (zero value eigenstate) of the Hamiltonian I construct an operator H^\dag H and look for its ground state with a standard dmrg algorithm. For some reason it tries to allocate a lot of memory and fails due to that reason. I’ve even tried to run it on a server with 128Gb available memory but still it fails. I feel stuck because I don’t even know how to debug such an issue, can you give me some advice or recommendation?

Sorry to hear it’s using so much memory - a good starting point to understand what’s going on would be to know what bond dimension the MPO has that you are inputting to DMRG then what bond dimension the MPS is reaching as the code runs. What values of those do you observe?

The bond dimension of the MPO is 25. The bond dimension of the MPS I use for the DMRG sweeps is really small, like 5. It does not proceed beyond the first link though.

Apparently, I have some progress with my problem. This topic helped a bit. Using

auto HdagH = nmultMPO(H, prime(Hdag));

instead of

auto HdagH = nmultMPO(prime(Hdag), H);

(I was absolutely confident that these commands are equivalent) and

HdagH.mapPrime(2, 1);

(also a very counter-intuitive thing to do)
changed an error to “size of initial vector should match linear matrix size”. However, if I consider a standard bosonic chain, then it works. But it does not work for my custom chain.

Another question to make sure I understand: is the MPO you are inputting into DMRG Hermitian or not Hermitian? I think you’re saying it’s supposed to be Hermitian (= H^\dagger H for some H) but I wanted to make sure I understood this part.

The reason I ask is that our C++ DMRG code can only handle Hermitian Hamiltonians and if you put in a non-Hermitian one it can lead to issues like what you are seeing. Our Julia DMRG code can directly work with non-Hermitian MPOs.

H is non-Hermitian, but H^\dag H which I pass to the DMRG algorithm is.

I might need to try out the case you are working on to see what may be causing the memory error, but one guess is that computing H^\dagger H numerically could be subtly causing it to be non-Hermitian even though theoretically it should be. This is because the MPO truncations involved are not necessarily guaranteed to preserve the Hermitian property.

Just mentioning that in case it helps.