<jupyter_start><jupyter_text>Caution !!! I heve not checked bugs yet.  Two-leg ladder spin-$1/2$ Heisenberg modelThis program uses the single-site unit cell VUMPS algorithm [1].This program corresponds to Fig. 5 in Ref. [3].  Reference[1] Phys. Rev. B 97, 045145 (2018)https://journals.aps.org/prb/abstract/10.1103/PhysRevB.97.045145https://arxiv.org/abs/1701.07035[2] Phys. Rev. B 86, 125441 (2012)https://journals.aps.org/prb/abstract/10.1103/PhysRevB.86.125441https://arxiv.org/abs/1204.0704[3] Phys. Rev. B 85, 075125 (2012)https://journals.aps.org/prb/abstract/10.1103/PhysRevB.85.075125https://arxiv.org/abs/0909.4059<jupyter_code>import numpy as np
import scipy as sp
from scipy.sparse.linalg import LinearOperator

import sys
sys.path.append('Library')
import MathFunctions as MF
import MPSOperators as MO
import SingleVUMPS as SV
dtype = np.dtype("float") 
#dtype = np.dtype("complex")

D = 2 # physical bond D = 2 corresponds to spin-1/2
M = 30 # virtual bond  
max_iter = 500 # the number of the VUMPS steps
D2 = D * D 

Sx,Sy,Sz,Su,Sd = MF.Spin(D)
Sz2 = np.einsum("ab,bc -> ac",Sz,Sz)

Jleg = 1.0; Jrung = 1.0;
h_xxx = ( np.kron(Sx,Sx) + np.kron(Sy,Sy) + np.kron(Sz,Sz) ).real
I = np.eye(D*D)
h = 0.5 * Jrung * np.kron(h_xxx,I).reshape(D2,D2,D2,D2) \
    + 0.5 * Jrung * np.kron(I,h_xxx).reshape(D2,D2,D2,D2) \
    + Jleg * np.kron(h_xxx,I).reshape(D,D,D,D,D,D,D,D).transpose(0,2,1,3,4,6,5,7).reshape(D**4,D**4).reshape(D2,D2,D2,D2) \
    + Jleg * np.kron(I,h_xxx).reshape(D,D,D,D,D,D,D,D).transpose(0,2,1,3,4,6,5,7).reshape(D**4,D**4).reshape(D2,D2,D2,D2)
# initial state
np.random.seed(3162277)
A = np.random.rand(M,D2,M) 
HR = np.random.rand(M,M); HL = np.random.rand(M,M)
tol = 1e-10

AC,C,AR,AL = MO.MixedCanonicalForm(A,dtype)

for i in range (max_iter):
    HR,er = SV.Calc_HR(AR,HR,h,dtype,tol=tol)
    HL,el = SV.Calc_HL(AL,HL,h,dtype,tol=tol)
    AC = SV.Next_AC(AC,AR,AL,HR,HL,h,dtype,tol=tol)
    C = SV.Next_C(C,AR,AL,HR,HL,h,dtype,tol=tol)
    #AR = SV.Next_AR_SVD(AC,C)
    AR = SV.Next_AR_PolarDecomposition(AC,C)
    #AL = SV.Next_AL_SVD(AC,C)
    AL = SV.Next_AL_PolarDecomposition(AC,C)
    B = SV.Calc_B(AC,C,AR,AL)
    tol = B / 100
    if B < 1e-12: 
        print ("Converged!")
        break
    E = 0.5 * ( er + el ).real
    print ("step {}, E {}, |B| {}".format(i,E/2,B))

#Tensors = np.empty(4, dtype=object)
#Tensors[0] = AC; Tensors[1] = C; Tensors[2] = AR; Tensors[3] = AL
#np.save("Data/Tensors_Heisenberg_{}_{}_{}_{}_{}.npy".format(D,M,JJ,BB,DD),Tensors)<jupyter_output>step 0, E 0.3747957383451404, |B| 0.6901806392842745
step 1, E -0.4612725447683895, |B| 0.06960191599807383
step 2, E -0.577781576589631, |B| 0.0009100835238643438
step 3, E -0.5780400087111095, |B| 0.000246792654496597
step 4, E -0.5780420727775445, |B| 0.0001562412785329729
step 5, E -0.5780421957134088, |B| 0.00016645682931763095
step 6, E -0.5780421952743988, |B| 0.00018238592750027382
step 7, E -0.5780421870690309, |B| 0.00018652472074347367
step 8, E -0.5780421977346266, |B| 0.00018773207064901848
step 9, E -0.5780422158046312, |B| 0.00017708586219277044
step 10, E -0.5780422368773837, |B| 0.00016986518790324093
step 11, E -0.5780422548967421, |B| 0.00016569172879399196
step 12, E -0.5780422695940132, |B| 0.0001646029161318995
step 13, E -0.5780422811436825, |B| 0.0001668291579029103
step 14, E -0.578042289845023, |B| 0.00017163751353287062
step 15, E -0.5780422962500299, |B| 0.00016765130351065126
step 16, E -0.5780423329285013, |B| 0.0001534460008433905
step 17, E -0.5780423571[...]