Hubbard 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 introduce the python interface of block2
. This interface allows the user to design custom workflow for running DMRG. It is also more convenient and efficient for model Hamiltonians like Hubbard and Heisenberg models, and non-standard quantum chemistry models with non-Hermitian or complex number integrals.
In block2
we support a few different symmetry modes. For fermionic models, if the Hamiltonian conserves the total spin \(S\), the projected spin \(S_z\), and the number of electrons \(N\), we can label the states using quantum numbers \((N, S_z)\) or \((N, S)\) or simply \(N\).
The symmetry group for \((N, S_z)\) is \(\mathrm{U(1)\times U(1)}\), which is Abelian. This symmetry mode is activated using
symm_type=SymmetryTypes.SZ
. We will first introduce the DMRG calculation using this symmetry since it is the most convenient choice.If we label states using \((N, S)\), the symmetry group is \(\mathrm{U(1)\times SU(2)}\). SU(2) is a non-Abelian symmetry, which is slightly difficult to work with, but it can be two to four times faster than the
SZ
mode. This symmetry mode is activated usingsymm_type=SymmetryTypes.SU2
.If we label states using only \(N\), we can remove the
alpha
andbeta
subscripts in the Hamiltonian, and the symmetry group is only \(\mathrm{U(1)}\). This symmetry mode is activated usingsymm_type=SymmetryTypes.SGF
(general spin fermionic). Note that the fermion statistics is considered in this mode.If we label states using only \(N\), we can also transform the Fermionic Hamiltonian to qubit Hamiltonian using JW transformation and then do DMRG on qubits without considering the fermion statistics. This symmetry mode is activated using
symm_type=SymmetryTypes.SGB
(general spin bosonic). This mode has the worst efficiency (most of the time the efficiency ofSGF
andSGB
should be the same), but it may show a clearer connection to other JW-based methods.For 2D Hubbard model with only nearest-neighbor interactions, the particle number is also a pseudo-spin (particle-hole) \(\mathrm{SU(2)}\) symmetry. If we label states using \((N, S_z)\), the symmetry group is \(\mathrm{SU(2)\times U(1)}\). This symmetry mode is activated using
symm_type=SymmetryTypes.SAnyPHSU2
.For 2D Hubbard model with only nearest-neighbor interactions, we can use the spin and pseudo-psin symmetries simultaneously. If we label states using \((N, S)\), the symmetry group is \(\mathrm{SU(2)\times SU(2)~/~Z_2} = \mathrm{SO(4)}\). This symmetry mode is activated using
symm_type=SymmetryTypes.SAnySO4
.
Next, we will explain the settings of the Hubbard Hamiltonian in each of the modes (1) (2) (3) (4) (5) and (6).
[2]:
import numpy as np
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
The SZ
Mode
Initialization
The SZ mode should be the most convenient for the Hubbard model. We first set a driver using this mode. scratch
sets a folder for writing temparary scratch files. symm_type
sets the symmetry mode (and whether the complex number should be used. Here we allow only real numbers). n_threads
sets the number of threads for shared-memory parallelization.
Next, we use the initialize_system
method to set the property of the target wavefunction (namely, which symmetry sector is targeted).
n_sites
is the number of sites \(L\) (or length of the lattice). Here each site has four local states: \(|0\rangle, |\alpha\rangle, |\beta\rangle\), and \(|\alpha\beta\rangle\).n_elec
is the number of electrons \(N\) in the target wavefunction. For half-filling Hubbard model, we should haven_elec == n_sites
.spin
is two times the total spin \(2S\) (inSU2
mode) or two times the projected spin \(2S_z\) (inSZ
mode) quantum number of the target wavefunction. Since we want the final wavefunction to have a equal number of alpha and beta electrons (namely, 4 alpha and 4 beta electrons for \(N=8\)), this number is set to zero here. Note that \(2S_z = N_{\alpha} - N_{\beta}\).
[3]:
L = 8
N = 8
TWOSZ = 0
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=N, spin=TWOSZ)
Build Hamiltonian
Now we set up the Hamiltonian. We first write the nearest-neighbor Hubbard model as the following:
Where we have considered the open boundary condition, and \(\sigma \in \{ \alpha, \beta\}\). In the code we use characters c
, d
, C
, D
to represent \(a^\dagger_\alpha, a_\alpha, a^\dagger_\beta, a_\beta\), respectively.
The method expr_builder
will return a ExprBuilder
object.
The method add_term
of the ExprBuilder
object has three parameters.
The first parameter is a string of characters in
cdCD
.The second parameter is a list of integers indicating the site indices associated with each charactor (operator) in the first parameter
The third parameter is a floating point number (or a list of floating point numbers), indicating the coefficient of the added term.
For example, b.add_term("dCc", [1, 3, 0], 0.7)
will add the term: \(0.7 \ a_{1\alpha} a^\dagger_{3\beta} a^\dagger_{0\alpha}\).
Most of the time, there can be multiple terms differ only in the operator indices and coefficients. For this case, one can invoke add_term
only once. So b.add_term("CD", [1, 3, 3, 1, 2, 2, 2, 4], [0.7, 0.6, 0.5, 0.4])
is equivalent to
b.add_term("CD", [1, 3], 0.7)
b.add_term("CD", [3, 1], 0.6)
b.add_term("CD", [2, 2], 0.5)
b.add_term("CD", [2, 4], 0.4)
b.finalize
will do some necessary operator reordering and simplification and driver.get_mpo
will create an optimal MPO for the given Hamiltonian. So we can set the MPO for the above Hubbard Hamiltonian using the following:
[4]:
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("cd", np.array([[[i, i + 1], [i + 1, i]] for i in range(L - 1)]).flatten(), -t)
b.add_term("CD", np.array([[[i, i + 1], [i + 1, i]] for i in range(L - 1)]).flatten(), -t)
# onsite term
b.add_term("cdCD", np.array([[i, ] * 4 for i in range(L)]).flatten(), U)
mpo = driver.get_mpo(b.finalize(), iprint=2)
Build MPO | Nsites = 8 | Nterms = 36 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 6 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Site = 1 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 2 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 3 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 4 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 5 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 6 / 8 .. Mmpo = 6 DW = 0.00e+00 NNZ = 11 SPT = 0.6944 Tmvc = 0.000 T = 0.002
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 6 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.019 Tmvc-total = 0.000 MPO bond dimension = 6 MaxDW = 0.00e+00
NNZ = 78 SIZE = 228 SPT = 0.6579
Rank = 0 Ttotal = 0.047 MPO method = FastBipartite bond dimension = 6 NNZ = 78 SIZE = 228 SPT = 0.6579
Note that the above method should allow the user to define arbitrary second-quantized fermionic (particle-number conserving) Hamiltonians (with possibly long-range and high-order terms).
Run DMRG
Finally, we can run DMRG to find the ground state energy.
We first use driver.get_random_mps
to create an initial guess for the wavefunction (MPS).
tag
is part of the scratch file name for this MPS. If there are multiple MPS objects, their tags must be all different.bond_dim
is the bond dimension of the initial guess MPS. The final optimized MPS may have a larger bond dimension.nroots
is the number of roots. Ifnroots > 1
, it will compute the ground state andnroots - 1
low-energy excited states.
The DMRG algorithm consists of many “sweeps”. The bond dimension of the MPS can increase after several sweeps. After a few tens of the sweeps, the MPS gradually converges to the ground state if nroots == 1
. To prevent from getting stuck in local minima, we may add noise in most sweeps, but the last sweep will not have any noise (namely, noises[-1] == 0
) so that the accuracy of the final printed energy is not affected. In each sweep, the DMRG algorithm will optimize \(L\) tensors in
the MPS, one by one. For optimizing each tensor, an effective Hamiltonian will be created and the Davidson algorithm is used to solve the eigenvalue problem of this effective Hamiltonian. The eigenvalue (energy) of this effective Hamiltonian is the same as the eigenvalue of the full many-body Hamiltonian (if the bond dimension of the MPS is big enough).
Therefore, to run a DMRG calculation, we need to set up a “sweep schedule”, namely, we need to set for each sweep: (a) the MPS bond dimension used in this sweep; (b) the noise used in this sweep; and (c) the Davidson algorithm convergence threshold for this sweep. Typically, these parameters should be changed every 4 to 5 sweeps (for Hubbard model with several hundred sites they may be changed every 30 to 50 sweeps), and the MPS bond dimension increases, while the noise and Davidson algorithm convergence decrease. The Davidson algorithm convergence cannot be set to zero.
bond_dims
, noises
, and thrds
are three lists. The first number in each list is the parameter used in the first sweep, the second number is for the second sweep, etc.
We call driver.dmrg
to execute the DMRG algorithm.
mpo
is the Hamiltonian represented as an MPO.ket
is the initial guess of the MPS. After DMRG finishes, the content of the MPS will be changed and it will be the optimized state.n_sweeps
is the maximal number of sweeps.bond_dims
is a list of integers indicating the MPS bond dimension for each sweep.noises
is a list of real numbers indicating the noise for each sweep.thrds
is a list of real numbers indicating the Davidson algorithm convergence threshold for each sweep.cutoff
is a real number. Singular values below this number will be discarded during the optimization. Setting this to zero will not discard any small singular value in the MPS unless the MPS bond dimension becomes higher than the set value.iprint
can be 0 (quiet), 1 (print energy after each sweep), or 2 (print energy after each iteration in each sweep).
[5]:
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-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
return driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, cutoff=0, iprint=1)
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.288 | E = -6.2256341447 | DW = 5.89e-16
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.403 | E = -6.2256341447 | DE = 2.22e-14 | DW = 7.51e-16
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.519 | E = -6.2256341447 | DE = -7.99e-15 | DW = 1.19e-16
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.636 | E = -6.2256341447 | DE = 8.88e-16 | DW = 1.18e-16
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.868 | E = -6.2256341447 | DE = 5.33e-15 | DW = 9.04e-29
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.130 | E = -6.2256341447 | DE = -1.78e-15 | DW = 1.10e-28
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.382 | E = -6.2256341447 | DE = -8.88e-16 | DW = 1.27e-28
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.619 | E = -6.2256341447 | DE = 8.88e-16 | DW = 1.32e-28
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 1.844 | E = -6.2256341447 | DE = 2.66e-15 | DW = 2.76e-51
DMRG energy = -6.225634144662398
Note that in the outputs, DE
represents the energy difference between two sweeps. DW
is the maximal discarded weight (sum of discarded eigenvalues of the density matrix) during the current sweep.
The remaining part of this docuementation introduces some other alternative symmetry to solve the same problem.
The SU2
Mode
Initialization
In this section, we try to solve the same problem using the SU(2) symmetry, which has a lower computational cost. Now spin
represents the total spin (instead of the projected spin).
[6]:
TWOS = 0
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=N, spin=TWOS)
Build Hamiltonian
We need to first do some rewriting of the Hamiltonian to use it with the SU(2) symmetry. This is mainly because operators like \(a_\alpha^\dagger\) do not conserve the total spin quantum number. We have to write the Hamiltonian in terms of only the total spin \(S\) preserving elementary operators.
First, note that
We can rewrite the Hubbard Hamiltonian as
Then we introduce the spin tensor operators which are spin-adaptation of elementary operators for alpha and beta spins. The superscript or subscript \([S]\) indicates the total spin quantum number of each spin tensor operator. Define
Because the total spin is a non-Abelian symmetry, we have the following notes:
The tensor product (coupling) of two spins may generate more than one possible resulting spins. For example, the tensor product between \(S=1/2\) (doublet) and \(S=1/2\) (doublet) can generate \(S=0\) (singlet) and \(S=1\) (triplet). Therefore, when we combine two spin operators \(\big(a_p^\dagger\big)^{[1/2]}\) and \(\big(a_q\big)^{[1/2]}\), we have to indicate which resulting total spin is desired. So we need the following notation to diffentiate the two cases:
There is no associative law. Namely, in general, \((S_1\otimes S_2) \otimes S_3 \neq S_1\otimes (S_2 \otimes S_3)\). The two products are connected by the spin recoupling, and the coefficients between the two products can be computed using Clebsch-Gordan factors. So when writing a string of spin tensor operators, we have to use parenthesis to indicate the coupling order of the operators.
Using spin tensor operators, we can rewrite the Hubbard Hamiltonian as (some derivation details can be found in https://block2.readthedocs.io/en/latest/theory/su2.html)
Next, we have to introduce a string notation to represent the above formula. We can use C
and D
for \(a^\dagger\) and \(a\), repsectively. Then we use +
to connect operators (namely, +
represents the tensor product). We use ()
to indicate the coupling (associative) order. We add a number (2 times the total spin) after each )
to indicate the resulting spin number of the tensor product. So, for example, (C+D)0
represents the singlet product between
\(a^\dagger\) and \(a\), and (C+D)2
represents the triplet product.
As a result, we can set the MPO for the above SU(2) Hubbard Hamiltonian using the following:
[7]:
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("(C+D)0", np.array([[[i, i + 1], [i + 1, i]] for i in range(L - 1)]).flatten(),
np.sqrt(2) * -t)
# onsite term
b.add_term("((C+(C+D)0)1+D)0", np.array([[i, ] * 4 for i in range(L)]).flatten(), U)
mpo = driver.get_mpo(b.finalize(), iprint=2)
Build MPO | Nsites = 8 | Nterms = 30 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Site = 1 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 2 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 3 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.001 T = 0.004
Site = 4 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 5 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 6 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.002
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.020 Tmvc-total = 0.001 MPO bond dimension = 4 MaxDW = 0.00e+00
NNZ = 50 SIZE = 104 SPT = 0.5192
Rank = 0 Ttotal = 0.045 MPO method = FastBipartite bond dimension = 4 NNZ = 50 SIZE = 104 SPT = 0.5192
Run DMRG
Finally, we can run DMRG as normal:
[8]:
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.122 | E = -6.2256341447 | DW = 1.62e-31
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.193 | E = -6.2256341447 | DE = -1.78e-15 | DW = 1.34e-31
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.259 | E = -6.2256341447 | DE = 8.88e-16 | DW = 2.12e-32
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.325 | E = -6.2256341447 | DE = 1.78e-15 | DW = 2.46e-31
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.396 | E = -6.2256341447 | DE = -7.99e-15 | DW = 0.00e+00
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.466 | E = -6.2256341447 | DE = 1.78e-15 | DW = 0.00e+00
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.587 | E = -6.2256341447 | DE = 4.44e-15 | DW = 0.00e+00
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.703 | E = -6.2256341447 | DE = -5.33e-15 | DW = 0.00e+00
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.796 | E = -6.2256341447 | DE = 1.78e-15 | DW = 0.00e+00
DMRG energy = -6.225634144676156
The SGF
Mode
In this section, we try to solve the same Hubbard problem using only particle number symmetry. The number of sites will now be the number of spin orbitals (which is equal to the number of qubits used in the next section), which is two times the number of the original (spatial) sites.
[9]:
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.initialize_system(n_sites=L * 2, n_elec=N, heis_twos=1)
Build Hamiltonian
Now we first split every Hubbard site into the alpha (odd) and beta (even) sites (next to each other). We have
We use characters C
and D
to represent \(a^\dagger\) and \(a\), respectively.
So we can set the MPO for the above Hubbard Hamiltonian using the following:
[10]:
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("CD", np.array([[[i, i + 2], [i + 2, i]] for i in range(2 * (L - 1))]).flatten(), -t)
# onsite term
b.add_term("CDCD", np.array([[2 * i, 2 * i, 2 * i + 1, 2 * i + 1] for i in range(L)]).flatten(), U)
mpo = driver.get_mpo(b.finalize(), iprint=2)
Build MPO | Nsites = 16 | Nterms = 36 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Site = 1 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.004
Site = 2 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 3 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.004
Site = 4 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 5 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 6 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 7 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.006
Site = 8 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 9 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 10 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 11 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 12 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 13 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 14 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.002
Site = 15 / 16 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.044 Tmvc-total = 0.000 MPO bond dimension = 7 MaxDW = 0.00e+00
NNZ = 128 SIZE = 560 SPT = 0.7714
Rank = 0 Ttotal = 0.095 MPO method = FastBipartite bond dimension = 7 NNZ = 128 SIZE = 560 SPT = 0.7714
Run DMRG
Finally, we can run DMRG as normal:
[11]:
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.921 | E = -6.2256341447 | DW = 2.51e-15
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.099 | E = -6.2256341447 | DE = -8.88e-15 | DW = 2.67e-15
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.272 | E = -6.2256341447 | DE = 1.78e-15 | DW = 2.33e-15
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.444 | E = -6.2256341447 | DE = 1.07e-14 | DW = 2.38e-15
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.853 | E = -6.2256341447 | DE = -8.88e-16 | DW = 9.93e-33
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.250 | E = -6.2256341447 | DE = -1.78e-15 | DW = 4.03e-32
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.697 | E = -6.2256341447 | DE = -8.88e-16 | DW = 1.85e-32
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.094 | E = -6.2256341447 | DE = -1.78e-15 | DW = 5.42e-32
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.330 | E = -6.2256341447 | DE = 3.55e-15 | DW = 0.00e+00
DMRG energy = -6.225634144669716
The SGB
Mode
Initialization
In this section, we try to solve the same Hubbard problem by first transfroming the model into a qubit (spin) model. The number of sites will now be the number of qubits (spins), which is two times the number of the original sites.
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.
[12]:
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGB, n_threads=4)
driver.initialize_system(n_sites=L * 2, n_elec=N, heis_twos=1)
Build Hamiltonian
Now we first need to do the Jordan-wigner transform. We split every Hubbard site into the alpha (odd) and beta (even) sites (next to each other). We have
We use charactors P
(plus), M
(minus), and Z
to represent operators \(S^+\), \(S^-\) and \(\frac{1}{2}S^z\), respectively. Note that we assume \(S^z|\alpha\rangle = 1\cdot |\alpha\rangle\) and \(S^z|\beta\rangle = -1\cdot |\beta\rangle\) here.
So we can set the MPO for the above Hubbard Hamiltonian using the following:
[13]:
t = 1
U = 2
b = driver.expr_builder()
# hopping term
b.add_term("PZZM",
np.array([[i, i, i + 1, i + 2] for i in range(2 * (L - 1))]).flatten(),
[-t * 4] * 2 * (L - 1))
b.add_term("MZZP",
np.array([[i, i, i + 1, i + 2] for i in range(2 * (L - 1))]).flatten(),
[+t * 4] * 2 * (L - 1))
# onsite term
b.add_term("PMPM",
np.array([[i, i] for i in range(2 * L)]).flatten(),
[U] * 2 * L)
mpo = driver.get_mpo(b.finalize(), iprint=2)
Build MPO | Nsites = 16 | Nterms = 36 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Site = 1 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.003
Site = 2 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 3 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 4 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 5 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 6 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 7 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 8 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.004
Site = 9 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 10 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.003
Site = 11 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 12 / 16 .. Mmpo = 7 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 13 / 16 .. Mmpo = 6 DW = 0.00e+00 NNZ = 9 SPT = 0.7857 Tmvc = 0.000 T = 0.002
Site = 14 / 16 .. Mmpo = 4 DW = 0.00e+00 NNZ = 6 SPT = 0.7500 Tmvc = 0.000 T = 0.002
Site = 15 / 16 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.040 Tmvc-total = 0.000 MPO bond dimension = 7 MaxDW = 0.00e+00
NNZ = 128 SIZE = 560 SPT = 0.7714
Rank = 0 Ttotal = 0.087 MPO method = FastBipartite bond dimension = 7 NNZ = 128 SIZE = 560 SPT = 0.7714
Run DMRG
Finally, we can run DMRG as normal:
[14]:
energies = run_dmrg(driver, mpo)
print('DMRG energy = %20.15f' % energies)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.925 | E = -6.2256341447 | DW = 4.38e-15
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.078 | E = -6.2256341447 | DE = 2.13e-14 | DW = 4.26e-15
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.236 | E = -6.2256341447 | DE = 0.00e+00 | DW = 4.25e-15
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.432 | E = -6.2256341447 | DE = 5.33e-15 | DW = 3.95e-15
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.785 | E = -6.2256341447 | DE = 0.00e+00 | DW = 2.54e-32
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.172 | E = -6.2256341447 | DE = 4.44e-15 | DW = 7.02e-32
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.570 | E = -6.2256341447 | DE = 0.00e+00 | DW = 1.86e-32
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.929 | E = -6.2256341447 | DE = 1.78e-15 | DW = 5.74e-32
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.157 | E = -6.2256341447 | DE = -3.55e-15 | DW = 0.00e+00
DMRG energy = -6.225634144672112
The PHSU2
Mode
Let
We have (assuming \(p\) and \(q\) differ by 1)
So
Now consider the onsite term
So the product of above terms
[15]:
t = 1
U = 2
n_elec = L
driver = DMRGDriver(scratch='./tmp', symm_type=SymmetryTypes.SAnyPHSU2, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=n_elec, spin=0)
b = driver.expr_builder()
b.add_term("(E+F)0", [g for i in range(L - 1) for g in [i, i + 1]], 2 ** 0.5)
b.add_term("(F+E)0", [g for i in range(L - 1) for g in [i, i + 1]], 2 ** 0.5)
b.add_term("((E+E)0+(F+F)0)0", np.array([i for i in range(L) for _ in range(4)]), 0.5 * U)
b.add_const(n_elec * U / 2)
mpo = driver.get_mpo(b.finalize(adjust_order=False), iprint=2)
ket = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
energies = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises, thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energies)
Build MPO | Nsites = 8 | Nterms = 22 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Site = 1 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 2 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 3 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.005
Site = 4 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.005
Site = 5 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 6 / 8 .. Mmpo = 4 DW = 0.00e+00 NNZ = 7 SPT = 0.5625 Tmvc = 0.000 T = 0.003
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.028 Tmvc-total = 0.000 MPO bond dimension = 4 MaxDW = 0.00e+00
NNZ = 50 SIZE = 104 SPT = 0.5192
Rank = 0 Ttotal = 0.069 MPO method = FastBipartite bond dimension = 4 NNZ = 50 SIZE = 104 SPT = 0.5192
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.144 | E = -6.2256341447 | DW = 2.13e-20
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.187 | E = -6.2256341447 | DE = 4.26e-14 | DW = 2.06e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.232 | E = -6.2256341447 | DE = -8.88e-15 | DW = 1.49e-20
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.274 | E = -6.2256341447 | DE = 1.78e-15 | DW = 2.86e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.317 | E = -6.2256341447 | DE = -3.55e-15 | DW = 1.51e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.359 | E = -6.2256341447 | DE = 0.00e+00 | DW = 1.15e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.402 | E = -6.2256341447 | DE = 1.78e-15 | DW = 1.96e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.444 | E = -6.2256341447 | DE = 7.11e-15 | DW = 1.37e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.478 | E = -6.2256341447 | DE = -2.84e-14 | DW = 1.94e-20
DMRG energy = -6.225634144668852
The SO4
Mode
Similarly, we can consider the \(\mathrm{SU(2) \times SU(2)}\) operator \(G\), where the row represents the pseudo spin symmetry and the column represents the spin symmetry.
So the components can be represented as
[16]:
t = 1
U = 2
n_elec = L
driver = DMRGDriver(scratch='./tmp', symm_type=SymmetryTypes.SAnySO4, n_threads=4)
driver.initialize_system(n_sites=L, n_elec=n_elec, spin=0)
b = driver.expr_builder()
b.add_term("(G[1,1]+G[1,1])[0,0]",
[g for i in range(L - 1) for g in [i, i + 1]], [2.0 if i % 2 == 0 else -2.0 for i in range(L - 1)])
b.add_term("((G[1,1]+G[1,1])[0,2]+(G[1,1]+G[1,1])[0,2])[0,0]",
np.array([i for i in range(L) for _ in range(4)]), U / (2 * 3 ** 0.5))
b.add_const(n_elec * U / 2)
mpo = driver.get_mpo(b.finalize(adjust_order=False), iprint=2)
ket = driver.get_random_mps(tag="KET", bond_dim=250, nroots=1)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
energies = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises, thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energies)
Build MPO | Nsites = 8 | Nterms = 15 | Algorithm = FastBIP | Cutoff = 1.00e-14
Site = 0 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 3 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Site = 1 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 2 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 3 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 4 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 5 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 6 / 8 .. Mmpo = 3 DW = 0.00e+00 NNZ = 5 SPT = 0.4444 Tmvc = 0.000 T = 0.003
Site = 7 / 8 .. Mmpo = 1 DW = 0.00e+00 NNZ = 3 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.022 Tmvc-total = 0.000 MPO bond dimension = 3 MaxDW = 0.00e+00
NNZ = 36 SIZE = 60 SPT = 0.4000
Rank = 0 Ttotal = 0.047 MPO method = FastBipartite bond dimension = 3 NNZ = 36 SIZE = 60 SPT = 0.4000
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.066 | E = -6.2256341447 | DW = 1.22e-20
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.113 | E = -6.2256341447 | DE = 3.02e-14 | DW = 2.80e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.159 | E = -6.2256341447 | DE = -3.55e-15 | DW = 7.31e-21
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.192 | E = -6.2256341447 | DE = 8.88e-15 | DW = 9.99e-21
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.229 | E = -6.2256341447 | DE = -1.95e-14 | DW = 1.22e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.261 | E = -6.2256341447 | DE = -1.78e-15 | DW = 5.33e-21
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.295 | E = -6.2256341447 | DE = 3.55e-15 | DW = 8.67e-21
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.333 | E = -6.2256341447 | DE = -1.78e-15 | DW = 7.17e-21
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.367 | E = -6.2256341447 | DE = -1.78e-15 | DW = 1.25e-20
DMRG energy = -6.225634144673016