# Quantum Chemistry Hamiltonians

```
[1]:
```

```
!pip install block2==0.5.2rc5 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
!pip install pyscf==2.1.1 -qq --progress-bar off
```

## Introduction

In this tutorial we explain how to perform quantum chemistry DMRG using the python interface of `block2`

.

The quantum chemistry Hamiltonian in its second quantized form has to be defined in a set of orbitals, such as the Hartree-Fock (or Density Functional Theory) orbitals. The symmetry that can be used in the DMRG calculation is thus has a dependence on the symmetry of the Hartree-Fock orbitals.

For spin-restricted Hartree-Fock (RHF) orbitals, we can perform spin-adapted DMRG (

`SU2`

mode in`block2`

) or non-spin-adapted DMRG with any lower symmetries (`SZ`

or`SGF`

).For spin-unrestricted Hartree-Fock (UHF) orbitals, we can perform non-spin-adapted DMRG (

`SZ`

mode in`block2`

) or DMRG with lower symmetries (such as`SGF`

).For general Hartree-Fock (GHF) orbitals, we can perform DMRG in spin-orbitals (

`SGF`

mode in`block2`

) or first translate the Hamiltonian into the qubit Hamiltonian then do DMRG (`SGB`

mode in`block2`

).For relativistic Dirac Hartree-Fock (DHF) orbitals, we can perform DMRG in complex spin-orbitals (

`SGFCPX`

mode in`block2`

).

Next, we will explain how to set up the integrals and perform DMRG in each of the modes (1) (2) (3) and (4). The quantum chemistry integrals will be generated using `pyscf`

and transformed using funtions defined in `pyblock2._pyscf.ao2mo`

.

```
[2]:
```

```
import numpy as np
from pyblock2._pyscf.ao2mo import integrals as itg
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
```

## Spin-Restricted Integrals

Here we use `get_rhf_integrals`

function to get the integrals. Note that in order to do DMRG in a CASCI space, one can set the `ncore`

(number of core orbitals) and `ncas`

(number of active orbitals) parameters in `get_*_integrals`

. `ncas=None`

will include all orbitals in DMRG.

For medium to large scale DMRG calculations, it is highly recommended to use a scratch space with high IO speed rather than the `./tmp`

used in the following example. One also needs to set a suitable `stack_mem`

in the `DMRGDriver`

constructor to set the memory used for storing renormalized operators (in bytes). The default is `stack_mem=int(1024**3)`

(1 GB). For medium scale calculations 10 to 30 GB might be required.

For the meaning of DMRG parameters, please have a look at the Hubbard - Run DMRG page.

```
[3]:
```

```
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
pdm1 = driver.get_1pdm(ket)
pdm2 = driver.get_2pdm(ket).transpose(0, 3, 1, 2)
print('Energy from pdms = %20.15f' % (np.einsum('ij,ij->', pdm1, h1e)
+ 0.5 * np.einsum('ijkl,ijkl->', pdm2, driver.unpack_g2e(g2e)) + ecore))
```

```
integral symmetrize error = 6.220297364062989e-14
integral cutoff error = 0.0
mpo terms = 972
Build MPO | Nsites = 10 | Nterms = 972 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 13 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.000 T = 0.009
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 59 SPT = 0.8665 Tmvc = 0.001 T = 0.009
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 117 SPT = 0.9386 Tmvc = 0.000 T = 0.009
Site = 3 / 10 .. Mmpo = 70 DW = 0.00e+00 NNZ = 363 SPT = 0.9074 Tmvc = 0.001 T = 0.009
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 217 SPT = 0.9613 Tmvc = 0.000 T = 0.007
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 201 SPT = 0.9733 Tmvc = 0.000 T = 0.008
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 169 SPT = 0.9667 Tmvc = 0.000 T = 0.007
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.006
Site = 8 / 10 .. Mmpo = 13 DW = 0.00e+00 NNZ = 41 SPT = 0.8949 Tmvc = 0.000 T = 0.008
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.000 T = 0.006
Ttotal = 0.078 Tmvc-total = 0.004 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1266 SIZE = 26498 SPT = 0.9522
Rank = 0 Ttotal = 0.154 MPO method = FastBipartite bond dimension = 94 NNZ = 1266 SIZE = 26498 SPT = 0.9522
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.107 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.640 | E = -107.6541224475 | DE = -1.07e-11 | DW = 4.78e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.154 | E = -107.6541224475 | DE = 2.84e-14 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.552 | E = -107.6541224475 | DE = -1.17e-12 | DW = 5.48e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.113 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.86e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.668 | E = -107.6541224475 | DE = 0.00e+00 | DW = 4.95e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.254 | E = -107.6541224475 | DE = -2.84e-14 | DW = 3.20e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.842 | E = -107.6541224475 | DE = 8.53e-14 | DW = 6.54e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 5.168 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.97e-20
DMRG energy = -107.654122447523932
Energy from pdms = -107.654122447523847
```

We can also run non-spin-adapted DMRG (`SZ`

mode) using the restricted integrals.

```
[4]:
```

```
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
```

```
integral symmetrize error = 9.016818189452533e-14
integral cutoff error = 0.0
mpo terms = 2778
Build MPO | Nsites = 10 | Nterms = 2778 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.001 T = 0.009
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 143 SPT = 0.9167 Tmvc = 0.001 T = 0.007
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 283 SPT = 0.9610 Tmvc = 0.001 T = 0.015
Site = 3 / 10 .. Mmpo = 138 DW = 0.00e+00 NNZ = 1023 SPT = 0.9326 Tmvc = 0.002 T = 0.015
Site = 4 / 10 .. Mmpo = 158 DW = 0.00e+00 NNZ = 535 SPT = 0.9755 Tmvc = 0.001 T = 0.011
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 463 SPT = 0.9842 Tmvc = 0.001 T = 0.011
Site = 6 / 10 .. Mmpo = 106 DW = 0.00e+00 NNZ = 415 SPT = 0.9790 Tmvc = 0.000 T = 0.006
Site = 7 / 10 .. Mmpo = 58 DW = 0.00e+00 NNZ = 163 SPT = 0.9735 Tmvc = 0.000 T = 0.004
Site = 8 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 87 SPT = 0.9423 Tmvc = 0.000 T = 0.004
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Ttotal = 0.087 Tmvc-total = 0.008 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.130 MPO method = FastBipartite bond dimension = 186 NNZ = 3164 SIZE = 102772 SPT = 0.9692
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.831 | E = -107.6541224475 | DW = 4.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.323 | E = -107.6541224475 | DE = -4.46e-12 | DW = 5.27e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.792 | E = -107.6541224475 | DE = -1.02e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.389 | E = -107.6541224475 | DE = 1.42e-12 | DW = 5.36e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.364 | E = -107.6541224475 | DE = -1.11e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.355 | E = -107.6541224475 | DE = 9.09e-13 | DW = 1.51e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.296 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.108 | E = -107.6541224475 | DE = -1.14e-12 | DW = 1.68e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.506 | E = -107.6541224475 | DE = -2.84e-14 | DW = 8.66e-20
DMRG energy = -107.654122447524529
```

We can also run DMRG in spin orbitals (`SGF`

mode) using the restricted integrals, which will be much slower (for more realistic systems).

```
[5]:
```

```
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.n_sites = ncas
g2e = driver.unpack_g2e(g2e)
orb_sym = [orb_sym[i // 2] for i in range(len(orb_sym) * 2)]
n_sites = ncas * 2
driver.initialize_system(n_sites=n_sites, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
```

```
integral symmetrize error = 9.016818189452533e-14
integral cutoff error = 0.0
mpo terms = 2438
Build MPO | Nsites = 20 | Nterms = 2438 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.010
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.009
Site = 2 / 20 .. Mmpo = 45 DW = 0.00e+00 NNZ = 45 SPT = 0.9500 Tmvc = 0.002 T = 0.008
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 131 SPT = 0.9530 Tmvc = 0.001 T = 0.010
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 159 SPT = 0.9683 Tmvc = 0.001 T = 0.008
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 203 SPT = 0.9759 Tmvc = 0.001 T = 0.007
Site = 6 / 20 .. Mmpo = 125 DW = 0.00e+00 NNZ = 265 SPT = 0.9796 Tmvc = 0.002 T = 0.008
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 974 SPT = 0.9382 Tmvc = 0.001 T = 0.009
Site = 8 / 20 .. Mmpo = 151 DW = 0.00e+00 NNZ = 188 SPT = 0.9901 Tmvc = 0.001 T = 0.006
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 344 SPT = 0.9846 Tmvc = 0.001 T = 0.010
Site = 10 / 20 .. Mmpo = 177 DW = 0.00e+00 NNZ = 246 SPT = 0.9906 Tmvc = 0.001 T = 0.006
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 402 SPT = 0.9872 Tmvc = 0.001 T = 0.006
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 306 SPT = 0.9883 Tmvc = 0.000 T = 0.006
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 242 SPT = 0.9835 Tmvc = 0.000 T = 0.009
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 150 SPT = 0.9805 Tmvc = 0.000 T = 0.006
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 94 SPT = 0.9774 Tmvc = 0.000 T = 0.004
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 82 SPT = 0.9611 Tmvc = 0.000 T = 0.008
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 50 SPT = 0.9359 Tmvc = 0.000 T = 0.004
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 20 SPT = 0.8571 Tmvc = 0.000 T = 0.003
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.140 Tmvc-total = 0.014 MPO bond dimension = 178 MaxDW = 0.00e+00
NNZ = 3934 SIZE = 200866 SPT = 0.9804
Rank = 0 Ttotal = 0.204 MPO method = FastBipartite bond dimension = 178 NNZ = 3934 SIZE = 200866 SPT = 0.9804
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.935 | E = -107.6541214061 | DW = 7.63e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.649 | E = -107.6541223348 | DE = -9.29e-07 | DW = 6.45e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.289 | E = -107.6541224347 | DE = -9.99e-08 | DW = 7.64e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.976 | E = -107.6541224347 | DE = 7.39e-13 | DW = 6.15e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.938 | E = -107.6541224379 | DE = -3.18e-09 | DW = 6.46e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.837 | E = -107.6541224379 | DE = -3.39e-11 | DW = 6.30e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.788 | E = -107.6541224379 | DE = 0.00e+00 | DW = 6.46e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.237 | E = -107.6541224379 | DE = 6.82e-13 | DW = 7.32e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 8.093 | E = -107.6541224379 | DE = 2.84e-14 | DW = 6.78e-20
DMRG energy = -107.654122437940956
```

## The `SZ`

Mode

Here we use the `get_uhf_integrals`

function to get the integrals.

```
[6]:
```

```
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_uhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
```

```
integral symmetrize error = 1.6638583155958613e-13
integral cutoff error = 0.0
mpo terms = 2778
Build MPO | Nsites = 10 | Nterms = 2778 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.002 T = 0.009
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 143 SPT = 0.9167 Tmvc = 0.001 T = 0.015
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 283 SPT = 0.9610 Tmvc = 0.002 T = 0.014
Site = 3 / 10 .. Mmpo = 138 DW = 0.00e+00 NNZ = 1023 SPT = 0.9326 Tmvc = 0.001 T = 0.012
Site = 4 / 10 .. Mmpo = 158 DW = 0.00e+00 NNZ = 535 SPT = 0.9755 Tmvc = 0.001 T = 0.011
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 463 SPT = 0.9842 Tmvc = 0.001 T = 0.008
Site = 6 / 10 .. Mmpo = 106 DW = 0.00e+00 NNZ = 415 SPT = 0.9790 Tmvc = 0.001 T = 0.007
Site = 7 / 10 .. Mmpo = 58 DW = 0.00e+00 NNZ = 163 SPT = 0.9735 Tmvc = 0.000 T = 0.004
Site = 8 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 87 SPT = 0.9423 Tmvc = 0.000 T = 0.003
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.086 Tmvc-total = 0.009 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.126 MPO method = FastBipartite bond dimension = 186 NNZ = 3164 SIZE = 102772 SPT = 0.9692
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.823 | E = -107.6541224475 | DW = 4.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.336 | E = -107.6541224475 | DE = -1.26e-11 | DW = 5.18e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.794 | E = -107.6541224475 | DE = -9.66e-13 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.282 | E = -107.6541224475 | DE = 1.22e-12 | DW = 5.37e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.835 | E = -107.6541224475 | DE = -9.95e-13 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.405 | E = -107.6541224475 | DE = 9.09e-13 | DW = 1.69e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.992 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.568 | E = -107.6541224475 | DE = -1.22e-12 | DW = 1.63e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.956 | E = -107.6541224475 | DE = -2.84e-14 | DW = 9.42e-20
DMRG energy = -107.654122447524443
```

## The `SGF`

Mode

Here we use the `get_ghf_integrals`

function to get the integrals.

```
[7]:
```

```
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.GHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
```

```
integral symmetrize error = 2.195446592452528e-13
integral cutoff error = 0.0
mpo terms = 5948
Build MPO | Nsites = 20 | Nterms = 5948 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.001 T = 0.024
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.003 T = 0.016
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 49 SPT = 0.9479 Tmvc = 0.003 T = 0.012
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 251 SPT = 0.9139 Tmvc = 0.003 T = 0.015
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 269 SPT = 0.9464 Tmvc = 0.003 T = 0.015
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 355 SPT = 0.9579 Tmvc = 0.003 T = 0.013
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 563 SPT = 0.9580 Tmvc = 0.004 T = 0.021
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 2316 SPT = 0.8575 Tmvc = 0.002 T = 0.021
Site = 8 / 20 .. Mmpo = 155 DW = 0.00e+00 NNZ = 214 SPT = 0.9890 Tmvc = 0.004 T = 0.027
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 708 SPT = 0.9691 Tmvc = 0.003 T = 0.013
Site = 10 / 20 .. Mmpo = 181 DW = 0.00e+00 NNZ = 398 SPT = 0.9851 Tmvc = 0.001 T = 0.021
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 718 SPT = 0.9777 Tmvc = 0.004 T = 0.023
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 496 SPT = 0.9810 Tmvc = 0.001 T = 0.031
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 410 SPT = 0.9721 Tmvc = 0.001 T = 0.006
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 254 SPT = 0.9670 Tmvc = 0.000 T = 0.005
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 118 SPT = 0.9716 Tmvc = 0.001 T = 0.004
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 124 SPT = 0.9411 Tmvc = 0.000 T = 0.004
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 74 SPT = 0.9051 Tmvc = 0.000 T = 0.003
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 20 SPT = 0.8571 Tmvc = 0.000 T = 0.003
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.282 Tmvc-total = 0.037 MPO bond dimension = 181 MaxDW = 0.00e+00
NNZ = 7370 SIZE = 204350 SPT = 0.9639
Rank = 0 Ttotal = 0.397 MPO method = FastBipartite bond dimension = 181 NNZ = 7370 SIZE = 204350 SPT = 0.9639
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.170 | E = -107.6541220534 | DW = 7.49e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.014 | E = -107.6541223446 | DE = -2.91e-07 | DW = 6.92e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 5.890 | E = -107.6541224348 | DE = -9.02e-08 | DW = 7.50e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 7.570 | E = -107.6541224348 | DE = 7.96e-13 | DW = 6.60e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 9.612 | E = -107.6541224379 | DE = -3.16e-09 | DW = 6.35e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 11.483 | E = -107.6541224379 | DE = 4.26e-13 | DW = 7.23e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 13.573 | E = -107.6541224379 | DE = 0.00e+00 | DW = 6.35e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 15.436 | E = -107.6541224379 | DE = -2.47e-12 | DW = 7.77e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 16.106 | E = -107.6541224379 | DE = -2.84e-14 | DW = 6.27e-20
DMRG energy = -107.654122437941297
```

## Relativistic DMRG

For relativistic DMRG, we use the `get_dhf_integrals`

function to get the integrals. We use the `SGFCPX`

Mode in `block2`

to execute DMRG. Note that the integrals, MPO, and MPS will all contain complex numbers in this mode.

```
[8]:
```

```
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.DHF(mol).set(with_gaunt=True, with_breit=True).run(conv_tol=1E-12)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_dhf_integrals(mf,
ncore=0, ncas=None, pg_symm=False)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGFCPX, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
```

```
integral symmetrize error = 0.0
integral cutoff error = 2.9170571933900254e-20
mpo terms = 26566
Build MPO | Nsites = 20 | Nterms = 26566 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.003 T = 0.013
Site = 1 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 23 SPT = 0.9087 Tmvc = 0.003 T = 0.018
Site = 2 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 320 SPT = 0.8186 Tmvc = 0.006 T = 0.026
Site = 3 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 440 SPT = 0.9105 Tmvc = 0.005 T = 0.024
Site = 4 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 605 SPT = 0.9200 Tmvc = 0.005 T = 0.025
Site = 5 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 728 SPT = 0.9375 Tmvc = 0.007 T = 0.029
Site = 6 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 992 SPT = 0.9438 Tmvc = 0.012 T = 0.041
Site = 7 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 1163 SPT = 0.9556 Tmvc = 0.008 T = 0.033
Site = 8 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 1971 SPT = 0.9480 Tmvc = 0.007 T = 0.031
Site = 9 / 20 .. Mmpo = 252 DW = 0.00e+00 NNZ = 2107 SPT = 0.9607 Tmvc = 0.007 T = 0.025
Site = 10 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 11022 SPT = 0.7947 Tmvc = 0.012 T = 0.055
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 1717 SPT = 0.9547 Tmvc = 0.003 T = 0.018
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 1641 SPT = 0.9373 Tmvc = 0.002 T = 0.015
Site = 13 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 1197 SPT = 0.9321 Tmvc = 0.002 T = 0.014
Site = 14 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 859 SPT = 0.9262 Tmvc = 0.002 T = 0.011
Site = 15 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 682 SPT = 0.9099 Tmvc = 0.001 T = 0.008
Site = 16 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 328 SPT = 0.9333 Tmvc = 0.001 T = 0.006
Site = 17 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 238 SPT = 0.8651 Tmvc = 0.000 T = 0.004
Site = 18 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 34 SPT = 0.8651 Tmvc = 0.000 T = 0.003
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.403 Tmvc-total = 0.086 MPO bond dimension = 252 MaxDW = 0.00e+00
NNZ = 26085 SIZE = 323082 SPT = 0.9193
Rank = 0 Ttotal = 0.500 MPO method = FastBipartite bond dimension = 252 NNZ = 26085 SIZE = 323082 SPT = 0.9193
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 27.037 | E = -107.6929206223 | DW = 6.66e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 39.215 | E = -107.6929209453 | DE = -3.23e-07 | DW = 7.27e-10
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 49.219 | E = -107.6929209492 | DE = -3.91e-09 | DW = 2.32e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 57.774 | E = -107.6929209492 | DE = -3.98e-13 | DW = 7.32e-10
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 75.576 | E = -107.6929209492 | DE = -5.12e-13 | DW = 6.34e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 91.974 | E = -107.6929209492 | DE = 2.84e-14 | DW = 1.05e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 112.150 | E = -107.6929209492 | DE = -2.84e-14 | DW = 7.11e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 130.356 | E = -107.6929209492 | DE = -2.84e-14 | DW = 1.12e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 138.078 | E = -107.6929209492 | DE = 2.84e-14 | DW = 7.25e-20
DMRG energy = -107.692920949169959
```

## Expectation and N-Particle Density Matrix

Once the optimized MPS is obtained, we can compute the expectation value on it, including its norm, the energy expectation, \(\langle S^2 \rangle\), N-particle density matrix, or any operator that can be constructed as an MPO.

In this example, we compute the triplet state.

```
[9]:
```

```
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
spin = 2
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
impo = driver.get_identity_mpo()
norm = driver.expectation(ket, impo, ket)
ener = driver.expectation(ket, mpo, ket)
print('Norm = %20.15f' % norm)
print('Energy expectation = %20.15f' % (ener / norm))
# <S^2> [ in spin-adapted mode this is always S(S+1) ]
ssq_mpo = driver.get_spin_square_mpo(iprint=0)
ssq = driver.expectation(ket, ssq_mpo, ket)
print('<S^2> expectation = %20.15f' % (ssq / norm))
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.632 | E = -106.9391328597 | DW = 3.38e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.963 | E = -106.9391328597 | DE = -4.92e-12 | DW = 8.92e-19
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.267 | E = -106.9391328597 | DE = 0.00e+00 | DW = 3.38e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.599 | E = -106.9391328597 | DE = 1.71e-13 | DW = 1.20e-18
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.964 | E = -106.9391328597 | DE = 0.00e+00 | DW = 1.12e-16
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.288 | E = -106.9391328597 | DE = 8.53e-14 | DW = 1.23e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.664 | E = -106.9391328597 | DE = 2.84e-14 | DW = 1.12e-16
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.984 | E = -106.9391328597 | DE = -2.84e-14 | DW = 1.23e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.236 | E = -106.9391328597 | DE = 0.00e+00 | DW = 1.72e-20
DMRG energy = -106.939132859666685
Norm = 0.999999999999999
Energy expectation = -106.939132859666685
<S^2> expectation = 2.000000000000000
```

We can also evaluate expectation of arbitray operator such as the occupancy in the first orbital

For the usage of `add_term`

function, please have a look at the Hubbard - Build Hamiltonian (SZ) and Hubbard - Build Hamiltonian (SU2) page.

```
[10]:
```

```
b = driver.expr_builder()
b.add_term("(C+D)0", [0, 0], np.sqrt(2))
n_mpo = driver.get_mpo(b.finalize(), iprint=0)
n_0 = driver.expectation(ket, n_mpo, ket)
print('N0 expectation = %20.15f' % (n_0 / norm))
```

```
N0 expectation = 1.999995824360303
```

We can then verify this number using 1PDM:

```
[11]:
```

```
pdm1 = driver.get_1pdm(ket)
print('N0 expectation from 1pdm = %20.15f' % pdm1[0, 0])
```

```
N0 expectation from 1pdm = 1.999995824360300
```

We can compute the 3PDM and compare the result with the FCI 3PDM. Note that in `pyscf`

the 3PDM is defined as

So we have to use the same convention in `block2`

by setting the `npdm_expr`

parameter in `block2`

to `((C+D)0+((C+D)0+(C+D)0)0)0`

.

```
[12]:
```

```
pdm3_b2 = driver.get_3pdm(ket, iprint=0, npdm_expr="((C+D)0+((C+D)0+(C+D)0)0)0")
from pyscf import fci
mx = fci.addons.fix_spin_(fci.FCI(mf), ss=2)
mx.kernel(h1e, g2e, ncas, nelec=n_elec, nroots=3, tol=1E-12)
print(mx.e_tot)
pdm3_fci = fci.rdm.make_dm123('FCI3pdm_kern_sf', mx.ci[0], mx.ci[0], ncas, n_elec)[2]
print('diff = ', np.linalg.norm(pdm3_fci - pdm3_b2))
```

```
[-106.93913286 -106.85412245 -106.70055113]
diff = 8.328760079421614e-06
```