# The lower bond-dim. of AKLT model

**URL:** <https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346>\
**Category:** ITensor Julia Questions\
**Created:** [April 16, 2025, 9:50am UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346 "2025-04-16T09:50:52Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![souhang](https://avatars.discourse-cdn.com/v4/letter/s/ba9def/32.png) [@souhang](https://itensor.discourse.group/u/souhang)\
**Post date:** [April 16, 2025, 9:50am UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/1 "2025-04-16T09:50:52Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [April 16, 2025, 1:56pm UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/2 "2025-04-16T13:56:47Z")

</div>

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](https://docs.itensor.org/ITensorMPS/stable/examples/DMRG.html) 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 |

---

<div class="post-metadata">

**Author:** ![souhang](https://avatars.discourse-cdn.com/v4/letter/s/ba9def/32.png) [@souhang](https://itensor.discourse.group/u/souhang)\
**Post date:** [April 16, 2025, 5:18pm UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/3 "2025-04-16T17:18:07Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [April 16, 2025, 7:35pm UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/4 "2025-04-16T19:35:48Z")

</div>

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](https://itensor.discourse.group/t/applying-dmrg-to-aklt-hamiltonian/2075/), 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](https://docs.itensor.org/ITensorMPS/stable/IncludedSiteTypes.html#%22S1%22-States)

# Edit

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

> **Code Used**
>
> ```julia
> 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-repl
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.

---

<div class="post-metadata">

**Author:** ![souhang](https://avatars.discourse-cdn.com/v4/letter/s/ba9def/32.png) [@souhang](https://itensor.discourse.group/u/souhang)\
**Post date:** [April 17, 2025, 4:30am UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/5 "2025-04-17T04:30:30Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![souhang](https://avatars.discourse-cdn.com/v4/letter/s/ba9def/32.png) [@souhang](https://itensor.discourse.group/u/souhang)\
**Post date:** [April 17, 2025, 5:21am UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/6 "2025-04-17T05:21:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [April 17, 2025, 6:15am UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/7 "2025-04-17T06:15:39Z")

</div>

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.

```julia
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

---

<div class="post-metadata">

**Author:** ![souhang](https://avatars.discourse-cdn.com/v4/letter/s/ba9def/32.png) [@souhang](https://itensor.discourse.group/u/souhang)\
**Post date:** [April 17, 2025, 8:04am UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/8 "2025-04-17T08:04:06Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![souhang](https://avatars.discourse-cdn.com/v4/letter/s/ba9def/32.png) [@souhang](https://itensor.discourse.group/u/souhang)\
**Post date:** [April 20, 2025, 10:36am UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/9 "2025-04-20T10:36:18Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![VinceNeede](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/vinceneede/32/731_2.png) [@VinceNeede](https://itensor.discourse.group/u/VinceNeede)\
**Post date:** [April 23, 2025, 5:37pm UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/10 "2025-04-23T17:37:09Z")

</div>

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](https://docs.itensor.org/ITensorMPS/stable/faq/DMRG.html#What-boundary-conditions-should-I-choose:-open,-periodic,-or-infinite?))

---

<div class="post-metadata">

**Author:** ![system](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/d3072b13e047cd06df7f3981547c0917d940dfa4.png) [@system](https://itensor.discourse.group/u/system)\
**Post date:** [May 3, 2025, 5:37pm UTC](https://itensor.discourse.group/t/the-lower-bond-dim-of-aklt-model/2346/11 "2025-05-03T17:37:49Z")

</div>

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