Out of Memory error when applying swap operator

ITensor Julia Questions
1

Hi,

The following code involves a swap operation between particles of two regions R1, R3. I am trying to implement Eq. 5 (https://journals.aps.org/prb/abstract/10.1103/PhysRevB.103.075102) in which V, W are defined in Eq. 10. By setting the angle to 0 in the code, one effectively just performs the swap operation S.
I used the ground state of the hofstadter model(Eq. 79) for a 5x6 system. However, the code gave an out of memory error even with a random mps. The linkdims of random_mps() is changed to match values for the hofstadter model(It was noted that having a linkdims around 250 leads to the error). Reducing the linkdims leads to poor results for the eigenspectrum when implementing the hofstadter model.

struct Region_parameters
    Nx::Int 
    Ny::Int
    R1_x_from::Int
    R1_x_to::Int
    R2_x_from::Int
    R2_x_to::Int
    R3_x_from::Int
    R3_x_to::Int
    y_from::Int
    y_to::Int
end

function xy_to_label(x, y, Nx, Ny)  
    return Ny*(x-1) + (y) 
end

#Define the number operator for the S=1/2 site. 
ITensors.op(::OpName"N",::SiteType"S=1/2")=
    [0.75 0 
     0    0.25]

function swap_one_copy(Rp::Region_parameters, sites::Any, psi::MPS, theta::Float64, prefactor::Any)
    
    N = length(psi) 
    L = Rp.R3_x_from - Rp.R1_x_from

    #REGION 1
    x = Rp.R1_x_from
    y = 1                                
    p = xy_to_label(x, y, Rp.Nx, Rp.Ny)  
    pL = xy_to_label(x + L, y, Rp.Nx, Rp.Ny)
    angle=0
    
    T1=psi[p]  
    T1*=exp(1im*angle*op(sites, "N", p)) 
    T1*=delta(prime(dag(sites[p])), sites[pL])  
    T1*=dag(prime(psi[pL], "Link"))  
    
    for y in 2:(Rp.Ny) 
        p=xy_to_label(x, y, Rp.Nx, Rp.Ny)
        pL=xy_to_label(x+L, y, Rp.Nx, Rp.Ny)
        angle=0
        
        T1*=psi[p]          
        T1*=exp(1im*angle*op(sites, "N", p))         
        T1*=delta(prime(dag(sites[p])), sites[pL])        
        T1*=dag(prime(psi[pL], "Link"))
            
    end

    
    for x in (Rp.R1_x_from+1):Rp.R1_x_to 
        for y in 1:(Rp.Ny) 
            p=xy_to_label(x, y, Rp.Nx, Rp.Ny)
            pL=xy_to_label(x+L, y, Rp.Nx, Rp.Ny)
            angle=0

            T1*=psi[p]  
            T1*=exp(1im*angle*op(sites, "N", p))
            T1*=delta(prime(dag(sites[p])), sites[pL])
            T1*=dag(prime(psi[pL], "Link"))
        end   
    end
    
    #REGION 2
    x = Rp.R2_x_from
    y = 1 
    p = xy_to_label(x, y, Rp.Nx, Rp.Ny)
    angle=0

    T2=psi[p]
    T2*=exp(1im*angle*op(sites, "N", p))
    T2*=dag(prime(psi[p]))  
    
    for y in 2:(Rp.Ny) 
        p = xy_to_label(x, y, Rp.Nx, Rp.Ny)
        angle=0
        
        T2*=psi[p]
        T2*=exp(1im*angle*op(sites, "N", p))
        T2*=dag(prime(psi[p])) 
        
    end

    for x in (Rp.R2_x_from+1):Rp.R2_x_to  
        
        for y=1: (Rp.Ny) 
            p=xy_to_label(x, y, Rp.Nx, Rp.Ny)
            angle=0

            T2*=psi[p]
            T2*=exp(1im*angle*op(sites, "N", p))
            T2*=dag(prime(psi[p])) 
        end
    end
    
    #REGION 3
    x = Rp.R3_x_from
    y = 1  
    p = xy_to_label(x, y, Rp.Nx, Rp.Ny)
    pL = xy_to_label(x-L, y, Rp.Nx, Rp.Ny)
    angle=0
    
    T3=psi[p]
    T3*=exp(1im*angle*op(sites, "N", p))
    T3*=delta(prime(dag(sites[p])), sites[pL])
    T3*=dag(prime(psi[pL], "Link"))
    
    for y in 2:(Rp.Ny)  
        p=xy_to_label(x, y, Rp.Nx, Rp.Ny)
        pL=xy_to_label(x-L, y, Rp.Nx, Rp.Ny)
        angle=0
        
        T3*=psi[p]
        T3*=exp(1im*angle*op(sites, "N", p))
        T3*=delta(prime(dag(sites[p])), sites[pL])
        T3*=dag(prime(psi[pL], "Link"))
    end
    
    for x in (Rp.R3_x_from+1):Rp.R3_x_to
            
        for y in 1:(Rp.Ny) 
            p=xy_to_label(x, y, Rp.Nx, Rp.Ny)
            pL=xy_to_label(x-L, y, Rp.Nx, Rp.Ny)
            angle=0
            
            T3*=psi[p]
            T3*=exp(1im*angle*op(sites, "N", p))
            T3*=delta(prime(dag(sites[p])), sites[pL])
            T3*=dag(prime(psi[pL], "Link"))
        end
    end 

    T=T1*T2*T3
    Ledge = xy_to_label(Rp.R1_x_from, 1, Rp.Nx, Rp.Ny)-1
    Redge = xy_to_label(Rp.R3_x_to, Rp.Ny, Rp.Nx, Rp.Ny)+1  

    for i in Ledge:-1:1  
        T*=psi[i]  
        T*=dag(prime(psi[i], "Link"))
    end

    for i in Redge:N  
        T*=psi[i]  
        T*=dag(prime(psi[i], "Link")) 
    end
           
    return complex(T)          
         
end

Nx=5
Ny=6
N=Nx*Ny 

Rp = Region_parameters(Nx,Ny,1,2,3,3,4,5,1,6) 

sites_spin=siteinds("S=1/2", N; conserve_qns=true) 
init_state=[isodd(n) ? "Up" : "Dn" for n in 1:N]
psi_GS=random_mps(sites_spin, init_state, linkdims=250) 


tensor=swap_one_copy(Rp, sites_spin, psi_GS, 0.1, 1.0)

Here are my doubts-

  1. Is there anything important to keep in mind when implementing swap operations?
  2. I wonder if my implementation is responsible for leading to the out of memory error. Or is it the value of linkdims?
    Any inputs about improving the code will be helpful. Thanks!
2

Sorry to hear this. I would suggest printing out information about your state or your tensors at each step (such as printing inds(T) for each ITensor), to find out what sizes they are just before you get the out of memory error. If you see bond dimensions reaching into the thousands, say, then it will let you back-track to what operations or parameters are leading to that happening.

3

I haven’t looked at your code in detail, but note that in general swap gates generate entanglement and therefore increase the bond dimension of your state, which could lead to you running out of memory if you don’t truncate enough (but truncating too much will eventually lead to inaccurate answers).

4

Thanks for your reply, I’ll try this out and give an update.

5

Makes sense, thanks for your reply.

6

I checked how different steps change the tensors. The out of memory error is caused at the line in Region 2: T2*=dag(prime(psi[p])) which creates an itensor with 4 link indices, each of which have a dimension of 250.
I’ll be trying to change the order of contraction to resolve this issue. The other option I think would be to run the code on a system with more memory.