The lower bond-dim. of AKLT model

ITensor Julia Questions
1

I’m doing DMRG simulation for AKLT model, and Heisenberg model. Here is my Hamiltonian for AKLT model, for example:

ITensors.space(::SiteType"S=1") = 3

ITensors.op(::OpName"Sx", ::SiteType"S=1") = (1/sqrt(2)) * [0 1 0; 1 0 1; 0 1 0]
ITensors.op(::OpName"Sy", ::SiteType"S=1") = (1/sqrt(2)) * [0 -1im 0; 1im 0 -1im; 0 1im 0]
ITensors.op(::OpName"Sz", ::SiteType"S=1") = [1 0 0; 0 0 0; 0 0 -1]

os = OpSum()
# OBC
for j = 1:N-1
    os += 1.0, "Sx", j, "Sx", j + 1
    os += 1.0, "Sy", j, "Sy", j + 1
    os += 1.0, "Sz", j, "Sz", j + 1
    os += 1 / 3, "Sx", j, "Sx", j + 1, "Sx", j, "Sx", j + 1
    os += 1 / 3, "Sy", j, "Sy", j + 1, "Sy", j, "Sy", j + 1
    os += 1 / 3, "Sz", j, "Sz", j + 1, "Sz", j, "Sz", j + 1
end

# PBC
os += 1.0, "Sx", 1, "Sx", N
os += 1.0, "Sy", 1, "Sy", N
os += 1.0, "Sz", 1, "Sz", N
os += 1 / 3, "Sx", 1, "Sx", N, "Sx", 1, "Sx", N
os += 1 / 3, "Sy", 1, "Sy", N, "Sy", 1, "Sy", N
os += 1 / 3, "Sz", 1, "Sz", N, "Sz", 1, "Sz", N

H = MPO(os, sites)

Now, I’m setting a dims-array, like dims=[200, 300, 400, 600], and then the simulation code is

for k in 1:length(dims)
    dim = dims[k]
    println("dim = ", dim)

    nsweeps = 10 # number of sweeps
    maxdim = [dim] # gradually increase states kept
    cutoff = [1E-15] # desired truncation error
    # if k == 1
    #     nsweeps = 20 # number of sweeps
    # end

    energy = 0.00
    if k == 1
        psi0 = random_mps(sites; linkdims=dim)
        energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)
    else
        energy, psi = dmrg(H, psi; nsweeps, maxdim, cutoff)
    end
end

I use dims=[300, 400, 600, 800] for L=128 to 192, then, I get the output

dim = 300
After sweep 1 energy=-175.6700890962815 maxlinkdim=300 maxerr=3.82E-03 time=537.875
After sweep 2 energy=-175.7735744967423 maxlinkdim=300 maxerr=1.44E-10 time=522.379
After sweep 3 energy=-175.78601502265704 maxlinkdim=300 maxerr=3.35E-12 time=521.770
After sweep 4 energy=-175.78880911048893 maxlinkdim=300 maxerr=3.68E-13 time=426.888
After sweep 5 energy=-175.78980642882556 maxlinkdim=287 maxerr=2.42E-15 time=233.428
After sweep 6 energy=-175.79032944824485 maxlinkdim=255 maxerr=1.00E-15 time=142.528
After sweep 7 energy=-175.79070715694272 maxlinkdim=261 maxerr=1.00E-15 time=110.557
After sweep 8 energy=-175.79097864921798 maxlinkdim=231 maxerr=1.00E-15 time=93.435
After sweep 9 energy=-175.79120153292308 maxlinkdim=215 maxerr=1.00E-15 time=89.942
After sweep 10 energy=-175.7913800158873 maxlinkdim=200 maxerr=1.00E-15 time=96.384
M = 1.8748952508287908
dim = 400
After sweep 1 energy=-175.79152353841388 maxlinkdim=257 maxerr=9.99E-16 time=102.785
After sweep 2 energy=-175.7916364828066 maxlinkdim=282 maxerr=1.00E-15 time=116.808
After sweep 3 energy=-175.7917243558656 maxlinkdim=241 maxerr=1.00E-15 time=121.927
After sweep 4 energy=-175.7917937179097 maxlinkdim=258 maxerr=1.00E-15 time=124.382
After sweep 5 energy=-175.79184741320583 maxlinkdim=222 maxerr=9.99E-16 time=119.985
After sweep 6 energy=-175.79188901694073 maxlinkdim=201 maxerr=1.00E-15 time=117.388
After sweep 7 energy=-175.79192096191525 maxlinkdim=199 maxerr=1.00E-15 time=116.917
After sweep 8 energy=-175.79194529746766 maxlinkdim=225 maxerr=1.00E-15 time=117.195
After sweep 9 energy=-175.79196390679022 maxlinkdim=261 maxerr=1.00E-15 time=120.797
After sweep 10 energy=-175.79197853005198 maxlinkdim=295 maxerr=9.99E-16 time=120.388
M = 1.6199579500487467
dim = 600
After sweep 1 energy=-175.7919895910903 maxlinkdim=200 maxerr=1.00E-15 time=119.369
After sweep 2 energy=-175.79199808350663 maxlinkdim=204 maxerr=9.99E-16 time=116.197
After sweep 3 energy=-175.79200452952952 maxlinkdim=200 maxerr=1.00E-15 time=114.137
After sweep 4 energy=-175.7920094574526 maxlinkdim=196 maxerr=1.00E-15 time=114.450
After sweep 5 energy=-175.79201323282263 maxlinkdim=200 maxerr=1.00E-15 time=115.410
After sweep 6 energy=-175.79201611908252 maxlinkdim=199 maxerr=1.00E-15 time=115.764
After sweep 7 energy=-175.79201830913698 maxlinkdim=202 maxerr=1.00E-15 time=114.662
After sweep 8 energy=-175.79202000072257 maxlinkdim=199 maxerr=1.00E-15 time=114.225
After sweep 9 energy=-175.79202129106199 maxlinkdim=198 maxerr=1.00E-15 time=114.275
After sweep 10 energy=-175.79202227249195 maxlinkdim=198 maxerr=9.99E-16 time=114.096
M = 1.5496199741177734
dim = 800
After sweep 1 energy=-175.79202301984716 maxlinkdim=199 maxerr=1.00E-15 time=116.524
After sweep 2 energy=-175.79202359423184 maxlinkdim=197 maxerr=1.00E-15 time=113.040
After sweep 3 energy=-175.79202403002986 maxlinkdim=197 maxerr=1.00E-15 time=113.935
After sweep 4 energy=-175.79202436428295 maxlinkdim=197 maxerr=1.00E-15 time=113.379
After sweep 5 energy=-175.79202461923285 maxlinkdim=198 maxerr=1.00E-15 time=113.582
After sweep 6 energy=-175.79202481307132 maxlinkdim=197 maxerr=1.00E-15 time=113.391
After sweep 7 energy=-175.79202496037252 maxlinkdim=198 maxerr=9.99E-16 time=113.416
After sweep 8 energy=-175.79202507216422 maxlinkdim=196 maxerr=1.00E-15 time=112.726
After sweep 9 energy=-175.79202515921943 maxlinkdim=198 maxerr=9.99E-16 time=112.785
After sweep 10 energy=-175.79202522452573 maxlinkdim=195 maxerr=1.00E-15 time=113.054
M = 1.5312736167728396

It seems that it gives no convergence, and the bond-dim is reducing step by step.

Why it appears? and how to solve it?

2

Hi @souhang,
I think you are starting with a too big bond dimension, at the beginning I found it better to start on a lower one (you can see from the examples in documentation to start with a bond dimension O(10)).

To me it looks like it just needs more iterations to converge, but it is going toward a convergence, you can see that from the fact that the energy is still decreasing (you should worry when it assumes values in a loop without decreasing, which means it is stuck in a local minimum).

It is normal for the bond dimension to decrease, remember that nearest neighbor 1D hamiltonians, when far from criticality, have ground states that are low entangled. Since you are starting very high in bond dimension it then needs to lower.

To check the convergence you can also check the relative variation of the energy:

\varepsilon_n = \left | \frac{E_n-E_{n-1}}{E_{n-1}} \right |
2 Likes
3

Thanks for your nice suggestion @VinceNeede

Here I have a naive question about that, the ground state of AKLT model is known VBS state, which the bond dim. should be only 2.

In my output file, the bond dim. is smaller than the setting bond dim, which seems converges the real bond dim of the VBS states? And the energy seems converges… Are there any explanation about it?

But I decide to receive your suggestion, to set smaller bond dim firstly, and I would improve the sweeping step to 15 or 20. Thanks.

4

I’m sorry I didn’t notice before, but are you sure the hamiltonian is correct? The term should be (\vec S_i \vec S_{i+1})^2 = (S^x_iS^x_{i+1}+S^y_iS^y_{i+1}+S^z_iS^z_{i+1})^2 which I don’t think it is equivalent to the one you wrote, the general writing should be this, I’m not sure if the additional terms cancel out.

Also, consider using S^+ and S^- instead of S^x and S^y, in this way the hamiltonian would be real and the computation should be faster.

Furthermore, the S=1 spin with the spin operators is already implemented in ITensors SiteTypes Included with ITensor · ITensorMPS.jl

Edit

I confirm that the error is the implementation of the Hamiltonian, This is the code I have used:

Code Used 
using ITensors, ITensorMPS

os = OpSum()
# OBC
for j = 1:N-1
    global os
    os += 1.0, "Sx", j, "Sx", j + 1
    os -= 1.0, "iSy", j, "iSy", j + 1
    os += 1.0, "Sz", j, "Sz", j + 1
    os += 1 / 3, "Sx", j, "Sx", j + 1, "Sx", j, "Sx", j + 1
    os += 1 / 3, "iSy", j, "iSy", j + 1, "iSy", j, "iSy", j + 1
    os += 1 / 3, "Sz", j, "Sz", j + 1, "Sz", j, "Sz", j + 1
    os -= 1 / 3, "Sx", j, "Sx", j + 1, "iSy", j, "iSy", j + 1
    os += 1 / 3, "Sx", j, "Sx", j + 1, "Sz", j, "Sz", j + 1
    os -= 1 / 3, "iSy", j, "iSy", j + 1, "Sx", j, "Sx", j + 1
    os -= 1 / 3, "iSy", j, "iSy", j + 1, "Sz", j, "Sz", j + 1
    os += 1 / 3, "Sz", j, "Sz", j + 1, "Sx", j, "Sx", j + 1
    os -= 1 / 3, "Sz", j, "Sz", j + 1, "iSy", j, "iSy", j + 1
end

# PBC
# os += 1.0, "Sx", 1, "Sx", N
# os -= 1.0, "iSy", 1, "iSy", N
# os += 1.0, "Sz", 1, "Sz", N
# os += 1 / 3, "Sx", 1, "Sx", N, "Sx", 1, "Sx", N
# os += 1 / 3, "iSy", 1, "iSy", N, "iSy", 1, "iSy", N
# os += 1 / 3, "Sz", 1, "Sz", N, "Sz", 1, "Sz", N
# os -= 1 / 3, "Sx", 1, "Sx", N, "iSy", 1, "iSy", N
# os += 1 / 3, "Sx", 1, "Sx", N, "Sz", 1, "Sz", N
# os -= 1 / 3, "iSy", 1, "iSy", N, "Sx", 1, "Sx", N
# os -= 1 / 3, "iSy", 1, "iSy", N, "Sz", 1, "Sz", N
# os += 1 / 3, "Sz", 1, "Sz", N, "Sx", 1, "Sx", N
# os -= 1 / 3, "Sz", 1, "Sz", N, "iSy", 1, "iSy", N

sites = siteinds("S=1", N)

H = MPO(os, sites)
psi0 = random_mps(sites, linkdims=2)
energy, psi = dmrg(H, psi0; nsweeps=10, maxdim=[10, 10, 20, 40, 80, 100]);

And these are the results obtained with N=128:

julia> include("exp.jl");
After sweep 1 energy=-84.0199504996155  maxlinkdim=10 maxerr=1.01E-04 time=0.383
After sweep 2 energy=-84.65892269585657  maxlinkdim=10 maxerr=5.88E-06 time=0.326
After sweep 3 energy=-84.66601973022019  maxlinkdim=20 maxerr=1.21E-10 time=0.745
After sweep 4 energy=-84.66650798857829  maxlinkdim=40 maxerr=3.71E-13 time=2.431
After sweep 5 energy=-84.66661612521872  maxlinkdim=80 maxerr=1.66E-14 time=6.331
After sweep 6 energy=-84.66665106329236  maxlinkdim=100 maxerr=3.63E-15 time=11.978
After sweep 7 energy=-84.6666626385208  maxlinkdim=100 maxerr=1.23E-15 time=14.806
After sweep 8 energy=-84.66666560968093  maxlinkdim=100 maxerr=1.66E-15 time=13.052
After sweep 9 energy=-84.66666639004092  maxlinkdim=100 maxerr=1.00E-15 time=9.586
After sweep 10 energy=-84.6666666048294  maxlinkdim=69 maxerr=2.22E-16 time=5.896

julia> energy, psi = dmrg(H, psi; nsweeps = 10, maxdim=100);
After sweep 1 energy=-84.66666665655427  maxlinkdim=51 maxerr=2.22E-16 time=2.880
After sweep 2 energy=-84.66666666534802  maxlinkdim=37 maxerr=2.21E-16 time=1.291
After sweep 3 energy=-84.66666666653006  maxlinkdim=25 maxerr=2.22E-16 time=0.592
After sweep 4 energy=-84.66666666665955  maxlinkdim=14 maxerr=2.20E-16 time=0.495
After sweep 5 energy=-84.66666666666681  maxlinkdim=6 maxerr=2.15E-16 time=0.200
After sweep 6 energy=-84.66666666666667  maxlinkdim=2 maxerr=1.50E-17 time=0.203
After sweep 7 energy=-84.66666666666609  maxlinkdim=2 maxerr=3.48E-22 time=0.149
After sweep 8 energy=-84.66666666666649  maxlinkdim=2 maxerr=1.27E-26 time=0.158
After sweep 9 energy=-84.66666666666688  maxlinkdim=2 maxerr=8.51E-29 time=0.154
After sweep 10 energy=-84.6666666666664  maxlinkdim=2 maxerr=9.59E-29 time=0.135

As you can see in less than 20 sweeps it has reached the convergence with a maxlinkdim of 2 (which corresponds to the VBS state)

For the future, I suggest you to use exact diagonalization for small systems and compare the results and the hamiltonian, so that you can spot these mistakes previously.

5

Thanks for your nice solution way @VinceNeede , I think my Hamiltonian is really incorrect.

But in your code, I’m not sure if the ladder operator should be defined as S^+=(1/2)(S^x+iS^y), but it seems that your code does not show the coefficients (1/2).
Maybe my notication is wrong…

6

I guess it may be written as this form:

for j in 1:N-1
# Add S_j ⋅ S_{j+1}
add!(ampo, “Sx”, j, “Sx”, j+1)
add!(ampo, “Sy”, j, “Sy”, j+1)
add!(ampo, “Sz”, j, “Sz”, j+1)
#Add (S_j ⋅ S_{j+1})^2 / 3
add!(ampo, 1/3, “Sx”, j, “Sx”, j+1, “Sx”, j, “Sx”, j+1)
add!(ampo, 1/3, “Sy”, j, “Sy”, j+1, “Sy”, j, “Sy”, j+1)
add!(ampo, 1/3, “Sz”, j, “Sz”, j+1, “Sz”, j, “Sz”, j+1)
add!(ampo, 1/3, “Sx”, j, “Sy”, j+1, “Sy”, j, “Sx”, j+1)
add!(ampo, 1/3, “Sx”, j, “Sz”, j+1, “Sz”, j, “Sx”, j+1)
add!(ampo, 1/3, “Sy”, j, “Sz”, j+1, “Sz”, j, “Sy”, j+1)
end

7

I’m sorry, I should have specified it. I’ve not used S^+ and S^- in my code at the end, but I have used S^x and iS^y, since S^y is purely imaginary, iS^y is real, the price is that if there are 2 S^y terms you have to subtract instead of adding. In fact, you can see the minus in the code.

os -= 1 / 3, "Sx", j, "Sx", j + 1, "iSy", j, "iSy", j + 1

About your last code, I’m not sure to understand how you obtained it

1 Like
8

OK, I think I get your point @VinceNeede , thanks for your nice answer again!

You use a pure matrix iS^y replacing the original Pauli matrix, and the coefficient is changed into -1/3 because each term has “double” S^y, and the minus comes from i*i.

I think I have understood this totally, thanks again!

1 Like
9

Moreover, if I want to simulate other models, like to replacing the factor 1/3 to \beta, where \beta can be equal to 0.70~1.20.

I found that when I set system size from 64 to 96 (PBCs), and bond dimon ranging from 600 to 1200, it consumes much time.

Are there any ways to improve it? Or I just needs more computation resources…

10

PBCs do require much more time, since the scaling goes from m to m^2 (where m is the maximum number of states kept during dmrg calculation). Consider if you actually need PBCs for what you want to study, you may also want to take a look at this (DMRG FAQs · ITensorMPS.jl)

11

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.