BHZ Model DMRG using ITensors.jl

Simulate the interacting 2D BHZ model using the Density Matrix Renormalization Group (DMRG) algorithm with the ITensors.jl package. Hamiltonian $H= H_0 + H_{\rm int}$:

$$H_0 = M \sum_{\mathbf{r},\sigma=\uparrow,\downarrow} (n_{\mathbf{r},s,\sigma} - n_{\mathbf{r},p,\sigma}) + B \sum_{\mathbf{r},\sigma=\uparrow,\downarrow}\sum_{\alpha=s,p} \left[ c^{\dagger}_{\mathbf{r},\alpha,\sigma}(\sigma^z)_{\alpha\alpha}c_{\mathbf{r}+\hat{x},\alpha,\sigma} + c^{\dagger}_{\mathbf{r},\alpha,\sigma}(\sigma^z)_{\alpha\alpha}c_{\mathbf{r} + \hat{y},\alpha,\sigma} + H.c. \right] + \frac{A}{2}\sum_{\mathbf{r},\sigma=\uparrow,\downarrow}\sum_{\alpha,\beta=s,p\\ \alpha\neq\beta}\left[(-1)^{\sigma}c^{\dagger}_{\mathbf{r},\alpha,\sigma}(\iota\sigma^{x})_{\alpha\beta}c_{\mathbf{r}+\hat{x},\beta,\sigma} + c^{\dagger}_{\mathbf{r},\alpha,\sigma}(-\iota\sigma^{y})_{\alpha\beta}c_{\mathbf{r}+\hat{y},\beta,\sigma} + h.c. \right]\ ,$$ $$H_{\rm int} = U \sum_{\mathbf{r},\alpha} n_{\mathbf{r},\alpha,\uparrow} n_{\mathbf{r},\alpha,\downarrow} + \left(U'-\frac{J_H}{2}\right) \sum_{\mathbf{r}, \alpha<\alpha'} n_{\mathbf{r},\alpha} n_{\mathbf{r},\alpha'} - 2J_H \sum_{\mathbf{r},\alpha <\alpha'} \mathbf{S}_{\mathbf{r},\alpha} \cdot \mathbf{S}_{\mathbf{r},\alpha'} + J_H\sum_{\mathbf{r},\alpha<\alpha'}\left(P^\dagger_{\mathbf{r},\alpha}P_{\mathbf{r},\alpha'}+ H.c.\right) \ .$$

DMRG can run for the canonical setup with fixed particle number $N_\uparrow$ and $N_\downarrow$ or for the grand canonical ensemble with additional chemical potential term $$\sum_{i,\sigma} \mu_{i,\sigma} \hat{n}_{i,\sigma} \ .$$ In the grand canonical case, where ITensors uses dense tensors, the code may be accelerated by using a GPU with CUDA.jl (--gpu flag).

Install dependencies

This package requires Julia 1.11.1 or later. To install the required packages, run the following command in the Julia REPL:

using Pkg
Pkg.add(["NPZ","JLD2","MKL","CUDA","ITensors","ArgParse","Printf","Dates","Random","HDF5","Logging","TimerOutputs","Adapt","ITensorMPS","KrylovKit","ProgressBars","Observers"])

Add chemical potential profile

To add a chemical potential profile function, edit the file ./BHZitensorsDMRG/chemical_potentials/mu_defs.jl by adding a function with signature (x,y, σ, α, L, W, JH, U, Up) -> Float64. Then add the function to the μ_funcs dictionary in the same file and give it a unique key. Here an additional paramater can be used and handed via the command line argument --mu-param in the run_bhz.jl script.

Example usage

To run a 1D ($W=1$) canonical simulation with $L=44$ sites, 42 spin-up and 42 spin-down electrons in the initial state, $U=1.0$, $m=-1.0$, $A=0.3$, $B=0.5$, $J_H=0.25U$, $U'=U-2J_H$ (measured in units of m=-M), and chemical potential according to the stepHFw6 function, use the following command:

julia -t 12 ./BHZitensorsDMRG/run_bhz.jl 44 42 42 --canonical --W 1 --U 1.0 --M -1.0 --A 0.3 --B 0.5 --out ./out --states-out-folder "./out/states" --save-states --mu-func-name stepHFw6 --cutoff 1e-7

This saves results to the folder ./out. In addition, the state is saved after each sweep to ./out/states.

To run the same simulation with a grand canonical ensemble on gpu, use:

julia -t 12 ./BHZitensorsDMRG/run_bhz.jl 44 42 42  --W 1 --U 1.0 --M -1.0 --A 0.3 --B 0.5 --out ./out --states-out-folder "./out/states" --save-states --mu-func-name stepHFw6 --cutoff 1e-7 --gpu

Usage

usage: run_bhz.jl [--W W] [--out FOLDER] [--save-states]
                  [--auto-restart] [--first-sweep FIRST_SWEEP]
                  [--write-when-maxdim-exceeds WRITE_WHEN_MAXDIM_EXCEEDS]
                  [--states-out-folder FOLDER] [--no-flush]
                  [--threaded-blocksparse] [--symmetrize-initial]
                  [--symmetrize-after-each-sweep] [--skip-dmrg]
                  [--skip-S2] [--cutoff CUTOFF] [--noise NOISE]
                  [--nmax NMAX] [--random-initial] [--U U] [--M M]
                  [--A A] [--B B] [--H H] [--Hp Hp] [--Utilde Utilde]
                  [--mu-func-name MU_FUNC_NAME] [--canonical]
                  [--y-periodic] [--BH-only] [--JH JH] [--Up Up]
                  [--SIA SIA] [--edge-potential EDGE_POTENTIAL]
                  [--half-edge] [--mu-param MU_PARAM] [--seed seed]
                  [--gpu] [--unified-memory] [-h] L Nup Ndn

positional arguments:
  L                     length of the lattice (should be ≧ W) (type:
                        Int64)
  Nup                   number of spin-up electrons (type: Int64)
  Ndn                   number of spin-down electrons (type: Int64)

optional arguments:
  --W W                 width of the lattice (type: Int64, default: 1)
  --out FOLDER          path to output folder
  --save-states         save state after each iteration (only keeps
                        last state). Also uses an observer with early
                        convergence checks.
  --auto-restart        tries to load a previously saved state as
                        initial state to dmrg. This can be used to
                        restart a non-finished dmrg computation that
                        used --save-states.
  --first-sweep FIRST_SWEEP
                        Set when restarting dmrg to start using dmrg
                        parameters from a later sweep. Set to the last
                        completet sweep for which the state is loaded
                        to restart with the paramters from the
                        previous run. (type: Int64, default: 1)
  --write-when-maxdim-exceeds WRITE_WHEN_MAXDIM_EXCEEDS
                        ITensorMPS.dmrg keyword argument
                        write_when_maxdim_exceeds:
                        `write_when_maxdim_exceeds::Int` - when the
                        allowed maxdim exceeds this value, begin
                        saving tensors to disk to free RAM memory in
                        large calculations. Default is -1, meaning
                        never. If set to valid integer, tensors will
                        be saved to --states-out-folder which defaults
                        to --out if not set. (type: Int64, default:
                        -1)
  --states-out-folder FOLDER
                        folder to save states to, defaults to --out 
  --threaded-blocksparse
                        use threaded blocksparse for ITensors 
  --skip-dmrg           load state but skip dmrg computation
  --skip-S2             skips measuring S^2
  --cutoff CUTOFF       cutoff for dmrg (type: Float64, default:
                        1.0e-7)
  --noise NOISE         noise for dmrg (type: Float64, default: 0.0)
  --nmax NMAX           maximum number of states to keep (type: Int64,
                        default: 10000)
  --random-initial      not use idmrg but a random intial state
  -h, --help            show this help message and exit

BHZ parameters:
  --U U                 Interaction strength U in units of t. (type:
                        Float64, default: 2.0)
  --M M                 Onsite energy M. (type: Float64, default: 1.0)
  --A A                 Interorbital coupling. (type: Float64,
                        default: 1.0)
  --B B                 Intraorbital coupling. (type: Float64,
                        default: 1.0)
  --H H                 Additional field at rx=1. (type: Float64,
                        default: 0.0)
  --Hp Hp               Add a term: -Hp S^z_{r_h,p}; Hp = 2JH *
                        <S^z_{r_h,s}>^{N0+2} (type: Float64, default:
                        0.0)
  --Utilde Utilde       Add a term: Utilde n_{r_h,p}; Utilde = (U' -
                        JH/2) * <n_{r_h,s}>^{N0+2} (type: Float64,
                        default: 0.0)
  --mu-func-name MU_FUNC_NAME
                        Dictionary key for chemical potential function
                        in ./chemical_potentials/mu_defs.jl:μ_funcs OR
                        a float used as chemical potential value.
                        (default: "zero")
  --canonical           Use fixed particle number Nup and Ndown.
  --y-periodic          If W>2, use periodic boundary conditions in
                        y-direction.
  --BH-only             Only use Bose-Hubbard part of the interaction
                        term.
  --JH JH               Interaction parameter JH. Default: 0.25 U
                        (type: Float64)
  --Up Up               Interaction parameter U'. Default: U - 2JH
                        (type: Float64)  
  --half-edge           Add sia and edge potential only to right side.
  --mu-param MU_PARAM   Potential parameter (type: Float64)

General parameters:
  --seed seed           seed for random number generator (type: Int64,
                        default: 12345)
  --gpu                 use gpu for dmrg, (cannot be used with
                        --canonical)
  --unified-memory      use cuda unified memory