Many term Hamiltonian to MPO

I am wondering if there are any tricks for speeding up (and possibly reducing the memory usage of) constructing a quantum-number conserving MPO from a Hamiltonian with millions of terms.

I’m trying to run construct the trancorrelated momentum space Fermi-Hubbard Hamiltonian as a MPO as done in Transcorrelated density matrix renormalization group, with the Hamiltonian shown in equation 8. My end goal is to run DMRG with the MPO, and so I when I construct the MPO I set conserve_qns=true so that the number of electrons and the total spin of my input DMRG state is preserved.

The form of the Hamiltonian isn’t particularly important, but for a 4x4 grid, there are almost 2 million terms in the Hamiltonian and the vast majority are three electron terms (so three creation and three annihilation operators). I have been trying in vain to construct the MPO using OpSum. The call to MPO(os, sites; maxdim=100) takes longer than I have been able to run so far. I’ve tried running on my laptop for two days to no avail. I’ve tried running on a cluster with 128 GB of RAM only to run out of memory after eight hours. Note that I’ve requested a maximum bond dimension of 100 so the resulting MPO should be quite small.

I have tried out the DMRG interface that takes a list of MPOs, but on a 3x3 grid this seems to be much slower. Both slower to construct the many MPOs than a single MPO with OpSum and slower to run DMRG. From some simple extrapolation I estimate that it would take half a day running on the cluster just to produce the list of MPOs for the 4x4 grid.

I am using MKL on the cluster and have tried playing around with ITensors.enable_threaded_blocksparse. Not that I can tell if either of them is an improvement because there’s no way for me to know how much progress was made before the program exits.

I’m willing to accept that this calculation is just too big (I’d actually like to repeat it many times, so a day of run time on a cluster is pushing it), but thought I’d ask around here before throwing in the towel. Much appreciated!

julia> versioninfo()
Julia Version 1.9.0
Commit 8e630552924 (2023-05-07 11:25 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin22.4.0)
  CPU: 16 × Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
  Threads: 1 on 16 virtual cores

julia> using Pkg; Pkg.status("ITensors")
Status `~/.julia/environments/v1.9/Project.toml`
  [9136182c] ITensors v0.3.34

It’s a known issue that the OpSum to MPO construction doesn’t scale well for Hamiltonians with O(N^4) terms or more. Miles could comment more on the potential for designing automated methods for improving the scaling of the MPO construction algorithm, but it is a common challenge for quantum chemistry and related Hamiltonians that constructing the MPO of a Hamiltonian with so many terms can become the bottleneck.

We are investigating using methods for splitting the OpSum into multiple groups of terms which each get compiled to MPOs based on papers like [2103.09976] Low communication high performance ab initio density matrix renormalization group algorithms in GitHub - ITensor/ITensorParallel.jl: Parallel tools for ITensors.jl., which can make the MPO construction faster and make the individual MPO bond dimensions smaller than the original MPO bond dimension would have been if you split the Hamiltonian into enough parts in the correct way (which isn’t always possible for a generic Hamiltonian), and then also allows you to parallelize DMRG over the sums of MPOs. Your mileage may vary depending on the particular Hamiltonian, you may have to think hard about how many terms you need to split your OpSum into and how to split your OpSum in order to see a benefit and determine if it is practically feasible to compile your Hamiltonian into MPOs.

In the paper you cite, they mention that you could choose to use a hybrid real/momentum space representation of the Hamiltonian, that could potentially decrease the number of terms needed or make it easier to split into multiple MPOs.

Much appreciated, I’ll take a look at those papers! The real space Hamiltonian (transcorrelated or not) only has a number of terms that scales as O(N) as opposed to O(N^5) for the transcorrelated momentum space hamiltonian, and I can run that no problem. The project I’m working on though is specifically interested in the momentum space representation.

Do you need it to be in momentum space in both directions or is one direction sufficient?

Additionally, it may be possible to analyze spectral properties of the state if you have the real space representation, see for example Phys. Rev. B 104, 115142 (2021) - Efficient matrix product state methods for extracting spectral information on rings and cylinders or Phys. Rev. Lett. 113, 010401 (2014) - Fourier Transform for Fermionic Systems and the Spectral Tensor Network.