Heisenberg Model
[1]:
!pip install block2==0.5.2rc13 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
Introduction
In this tutorial we explain how to solve the Heisenberg model using the python interface of block2
.
First, we have to define the “site-spin” of the model. The parameter heis_twos
represents two times the spin in each site. Namely, for \(S = 1/2\) Heisenberg model, heis_twos = 1
, for \(S = 1\) Heisenberg model, heis_twos = 2
, etc. Note that arbitrary non-negative half-integer and integer \(S\) can be supported in block2
.
Second, we can solve the model using the SU2 symmetry (SU2
mode) or U1 symmetry (SGB
mode). The SU2 symmetry can be more efficient (for large S
) and can generate states with well-defined total spin symmetry, but requires some additional rearrangement of the Hamiltonian terms.
The SGB
Mode
The Hamiltonian is
We can solve this model using the following code:
[2]:
import numpy as np
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
L = 100
heis_twos = 1
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGB, n_threads=4)
driver.initialize_system(n_sites=L, heis_twos=heis_twos, heis_twosz=0)
b = driver.expr_builder()
for i in range(L - 1):
b.add_term("PM", [i, i + 1], 0.5)
b.add_term("MP", [i, i + 1], 0.5)
b.add_term("ZZ", [i, i + 1], 1.0)
heis_mpo = driver.get_mpo(b.finalize(), iprint=0)
def run_dmrg(driver, mpo):
ket = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-5] * 4 + [1e-6] * 2 + [0]
thrds = [1e-6] * 4 + [1e-8] * 4
return driver.dmrg(
mpo,
ket,
n_sweeps=8,
bond_dims=bond_dims,
noises=noises,
thrds=thrds,
cutoff=1E-24,
iprint=1,
)
energies = run_dmrg(driver, heis_mpo)
print('DMRG energy = %20.15f' % energies)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 32.844 | E = -44.1201487994 | DW = 1.12e-09
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 38.582 | E = -44.1277383773 | DE = -7.59e-03 | DW = 1.21e-13
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 41.882 | E = -44.1277383773 | DE = -4.01e-12 | DW = 7.28e-14
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 44.918 | E = -44.1277383773 | DE = -1.34e-12 | DW = 2.79e-14
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-06 | Dav threshold = 1.00e-08
Time elapsed = 60.141 | E = -44.1277398826 | DE = -1.51e-06 | DW = 4.18e-17
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-06 | Dav threshold = 1.00e-08
Time elapsed = 79.821 | E = -44.1277398826 | DE = -1.92e-13 | DW = 3.84e-17
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-08
Time elapsed = 94.075 | E = -44.1277398826 | DE = -1.35e-13 | DW = 1.61e-20
DMRG energy = -44.127739882610513
The SU2
Mode
To solve the model in the SU2
mode, we define the following (triplet) spin tensor operator:
Then we have
Note that in the above calculation, we have used the following CG factors:
[3]:
from block2 import SU2CG
print('<Jp=1,Jzp=+1; Jq=1,Jzq=-1|J=0,Jz=0> = ', SU2CG().cg(2, 2, 0, 2, -2, 0))
print('<Jp=1,Jzp=-1; Jq=1,Jzq=+1|J=0,Jz=0> = ', SU2CG().cg(2, 2, 0, -2, 2, 0))
print('<Jp=1,Jzp= 0; Jq=1,Jzq= 0|J=0,Jz=0> = ', SU2CG().cg(2, 2, 0, 0, 0, 0))
print(1 / 3 ** 0.5)
<Jp=1,Jzp=+1; Jq=1,Jzq=-1|J=0,Jz=0> = 0.5773502691896257
<Jp=1,Jzp=-1; Jq=1,Jzq=+1|J=0,Jz=0> = 0.5773502691896257
<Jp=1,Jzp= 0; Jq=1,Jzq= 0|J=0,Jz=0> = -0.5773502691896257
0.5773502691896258
So the Hamiltonian in SU2
notation is
We can solve this model using the following code:
[5]:
import numpy as np
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
L = 100
heis_twos = 1
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=L, heis_twos=heis_twos, spin=0)
b = driver.expr_builder()
for i in range(L - 1):
b.add_term("(T+T)0", [i, i + 1], -np.sqrt(3) / 2)
heis_mpo = driver.get_mpo(b.finalize(adjust_order=False), iprint=0)
energies = run_dmrg(driver, heis_mpo)
print('DMRG energy = %20.15f' % energies)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 48.587 | E = -44.1275743858 | DW = 7.23e-12
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 53.973 | E = -44.1277393752 | DE = -1.65e-04 | DW = 6.90e-14
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 59.107 | E = -44.1277393752 | DE = -3.84e-13 | DW = 3.59e-16
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-05 | Dav threshold = 1.00e-06
Time elapsed = 63.144 | E = -44.1277393752 | DE = -1.28e-13 | DW = 8.81e-17
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-06 | Dav threshold = 1.00e-08
Time elapsed = 75.269 | E = -44.1277398850 | DE = -5.10e-07 | DW = 3.06e-22
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-06 | Dav threshold = 1.00e-08
Time elapsed = 94.287 | E = -44.1277398850 | DE = -1.99e-13 | DW = 6.38e-23
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-08
Time elapsed = 103.865 | E = -44.1277398850 | DE = 0.00e+00 | DW = 1.42e-23
DMRG energy = -44.127739885005937