# Quantum Chemistry Hamiltonians

```
[1]:
```

```
!pip install block2==0.5.2rc13 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
!pip install pyscf==2.3.0 -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`

).For atom and diatomic molecules, we can perform spin-adapted/non-spin-adapted/spin-orbital DMRG (

`SAnySU2LZ/SAnySZLZ/SAnySGFLZ`

modes in`block2`

) with the \(L_z\) symmetry.

Next, we will explain how to set up the integrals and perform DMRG in each of the modes (1) (2) (3) (4) and (5). 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))
impo = driver.get_identity_mpo()
expt = driver.expectation(ket, mpo, ket) / driver.expectation(ket, impo, ket)
print('Energy from expectation = %20.15f' % expt)
```

```
integral symmetrize error = 5.977257773040786e-14
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | 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.010
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.000 T = 0.009
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.009
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.000 T = 0.010
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.007
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.009
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.004
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.007
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.006
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.010
Ttotal = 0.082 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.197 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.239 | 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.842 | E = -107.6541224475 | DE = -5.40e-12 | DW = 5.66e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.663 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.203 | E = -107.6541224475 | DE = -7.96e-13 | DW = 5.72e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.918 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.92e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.457 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.42e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.188 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.01e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.683 | E = -107.6541224475 | DE = -2.84e-14 | DW = 5.57e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.061 | E = -107.6541224475 | DE = 2.84e-14 | DW = 2.64e-20
DMRG energy = -107.654122447524443
Energy from pdms = -107.654122447524387
Energy from expectation = -107.654122447524344
```

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 = 7.512924107430347e-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.005
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 143 SPT = 0.9167 Tmvc = 0.001 T = 0.006
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 283 SPT = 0.9610 Tmvc = 0.001 T = 0.008
Site = 3 / 10 .. Mmpo = 138 DW = 0.00e+00 NNZ = 1023 SPT = 0.9326 Tmvc = 0.004 T = 0.013
Site = 4 / 10 .. Mmpo = 158 DW = 0.00e+00 NNZ = 535 SPT = 0.9755 Tmvc = 0.001 T = 0.010
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 463 SPT = 0.9842 Tmvc = 0.001 T = 0.010
Site = 6 / 10 .. Mmpo = 106 DW = 0.00e+00 NNZ = 415 SPT = 0.9790 Tmvc = 0.000 T = 0.010
Site = 7 / 10 .. Mmpo = 58 DW = 0.00e+00 NNZ = 163 SPT = 0.9735 Tmvc = 0.000 T = 0.025
Site = 8 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 87 SPT = 0.9423 Tmvc = 0.000 T = 0.010
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Ttotal = 0.103 Tmvc-total = 0.009 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.178 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 = 1.281 | 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 = 2.089 | E = -107.6541224475 | DE = -1.31e-11 | DW = 5.10e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.909 | E = -107.6541224475 | DE = -1.42e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.136 | E = -107.6541224475 | DE = 1.42e-12 | DW = 5.19e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.572 | E = -107.6541224475 | DE = -1.17e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.850 | E = -107.6541224475 | DE = 9.95e-13 | DW = 1.60e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.709 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.581 | E = -107.6541224475 | DE = -1.25e-12 | DW = 1.60e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 9.169 | E = -107.6541224475 | DE = -2.84e-14 | DW = 5.07e-20
DMRG energy = -107.654122447524500
```

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 = 7.512924107430346e-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.005
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.006
Site = 2 / 20 .. Mmpo = 45 DW = 0.00e+00 NNZ = 45 SPT = 0.9500 Tmvc = 0.001 T = 0.005
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 131 SPT = 0.9530 Tmvc = 0.001 T = 0.005
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 159 SPT = 0.9683 Tmvc = 0.001 T = 0.005
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 203 SPT = 0.9759 Tmvc = 0.001 T = 0.006
Site = 6 / 20 .. Mmpo = 125 DW = 0.00e+00 NNZ = 265 SPT = 0.9796 Tmvc = 0.001 T = 0.006
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 974 SPT = 0.9382 Tmvc = 0.001 T = 0.008
Site = 8 / 20 .. Mmpo = 151 DW = 0.00e+00 NNZ = 188 SPT = 0.9901 Tmvc = 0.001 T = 0.010
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 344 SPT = 0.9846 Tmvc = 0.001 T = 0.037
Site = 10 / 20 .. Mmpo = 177 DW = 0.00e+00 NNZ = 246 SPT = 0.9906 Tmvc = 0.001 T = 0.007
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.001 T = 0.005
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 242 SPT = 0.9835 Tmvc = 0.000 T = 0.004
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 150 SPT = 0.9805 Tmvc = 0.000 T = 0.004
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 94 SPT = 0.9774 Tmvc = 0.000 T = 0.003
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 82 SPT = 0.9611 Tmvc = 0.000 T = 0.003
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.002
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.134 Tmvc-total = 0.012 MPO bond dimension = 178 MaxDW = 0.00e+00
NNZ = 3934 SIZE = 200866 SPT = 0.9804
Rank = 0 Ttotal = 0.186 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.861 | E = -107.6541216205 | DW = 7.62e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.471 | E = -107.6541223420 | DE = -7.21e-07 | DW = 6.45e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.030 | E = -107.6541224347 | DE = -9.27e-08 | DW = 7.64e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.597 | E = -107.6541224347 | DE = 8.81e-13 | DW = 6.15e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.707 | E = -107.6541224379 | DE = -3.17e-09 | DW = 6.46e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.003 | E = -107.6541224379 | DE = -3.37e-11 | DW = 7.34e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.306 | 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.088 | E = -107.6541224379 | DE = 7.11e-13 | DW = 7.10e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 7.514 | E = -107.6541224379 | DE = 2.84e-14 | DW = 8.86e-20
DMRG energy = -107.654122437940984
```

## Read and Write FCIDUMP Files

Instead of generating integrals (`h1e`

and `g2e`

) using `pyscf`

, we can also read these integrals from a FCIDUMP file (which can be generated using any of many other quantum chemistry packages) then perform DMRG. Additionally, we also provide methods to write the FCIDUMP file using the data in the `h1e`

and `g2e`

arrays.

After invoking `driver.read_fcidump`

, the integrals and target state infomation can be obtained from `driver.h1e`

, `driver.g2e`

, `driver.n_sites`

, etc.

```
[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.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)
# write integrals to file
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
driver.write_fcidump(h1e, g2e, ecore, filename='N2.STO3G.FCIDUMP', h1e_symm=True, pg='d2h')
# read integrals from file
driver.read_fcidump(filename='N2.STO3G.FCIDUMP', pg='d2h')
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec,
spin=driver.spin, orb_sym=driver.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=driver.h1e, g2e=driver.g2e, ecore=driver.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)
```

```
symmetrize error = 2.3300000000000018e-14
integral symmetrize error = 0.0
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | 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.011
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.000 T = 0.018
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.032
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.001 T = 0.007
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.004
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.012
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.004
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.007
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.004
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Ttotal = 0.103 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.171 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.556 | 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 = 0.903 | E = -107.6541224475 | DE = -3.52e-12 | DW = 4.66e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.295 | 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 = 1.686 | E = -107.6541224475 | DE = -5.97e-13 | DW = 5.90e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.136 | E = -107.6541224475 | DE = -5.68e-14 | DW = 2.02e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.523 | E = -107.6541224475 | DE = 0.00e+00 | DW = 4.86e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.946 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.08e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.342 | E = -107.6541224475 | DE = -2.84e-14 | DW = 5.58e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.607 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.33e-20
DMRG energy = -107.654122447524614
```

## The `SZ`

Mode

Here we use `get_uhf_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.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.5366482610151817e-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.001 T = 0.018
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 143 SPT = 0.9167 Tmvc = 0.001 T = 0.014
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 283 SPT = 0.9610 Tmvc = 0.001 T = 0.012
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.012
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 463 SPT = 0.9842 Tmvc = 0.001 T = 0.007
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.003
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.093 Tmvc-total = 0.007 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.141 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.783 | 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.233 | E = -107.6541224475 | DE = -1.78e-11 | DW = 5.23e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.662 | E = -107.6541224475 | DE = -1.19e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.104 | E = -107.6541224475 | DE = 1.28e-12 | DW = 5.02e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.629 | E = -107.6541224475 | DE = -6.82e-13 | DW = 3.61e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.139 | E = -107.6541224475 | DE = 8.53e-13 | DW = 1.92e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.672 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.173 | E = -107.6541224475 | DE = -1.68e-12 | DW = 2.25e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.529 | E = -107.6541224475 | DE = 2.84e-14 | DW = 7.36e-20
DMRG energy = -107.654122447524529
```

## The `SGF`

Mode

Here we use `get_ghf_integrals`

function to get the integrals.

```
[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.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.1082787907764326e-13
integral cutoff error = 0.0
mpo terms = 5914
Build MPO | Nsites = 20 | Nterms = 5914 | 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.007
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.009
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 49 SPT = 0.9479 Tmvc = 0.002 T = 0.010
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 251 SPT = 0.9139 Tmvc = 0.002 T = 0.013
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 273 SPT = 0.9456 Tmvc = 0.004 T = 0.019
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 357 SPT = 0.9576 Tmvc = 0.003 T = 0.013
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 563 SPT = 0.9580 Tmvc = 0.001 T = 0.012
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 2318 SPT = 0.8574 Tmvc = 0.001 T = 0.014
Site = 8 / 20 .. Mmpo = 155 DW = 0.00e+00 NNZ = 208 SPT = 0.9893 Tmvc = 0.002 T = 0.012
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 686 SPT = 0.9701 Tmvc = 0.001 T = 0.014
Site = 10 / 20 .. Mmpo = 181 DW = 0.00e+00 NNZ = 388 SPT = 0.9855 Tmvc = 0.001 T = 0.010
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 720 SPT = 0.9777 Tmvc = 0.001 T = 0.011
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 496 SPT = 0.9810 Tmvc = 0.001 T = 0.006
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 412 SPT = 0.9720 Tmvc = 0.000 T = 0.004
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 234 SPT = 0.9696 Tmvc = 0.000 T = 0.005
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 118 SPT = 0.9716 Tmvc = 0.000 T = 0.003
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 124 SPT = 0.9411 Tmvc = 0.000 T = 0.003
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 76 SPT = 0.9026 Tmvc = 0.000 T = 0.004
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 22 SPT = 0.8429 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.176 Tmvc-total = 0.023 MPO bond dimension = 181 MaxDW = 0.00e+00
NNZ = 7328 SIZE = 204350 SPT = 0.9641
Rank = 0 Ttotal = 0.269 MPO method = FastBipartite bond dimension = 181 NNZ = 7328 SIZE = 204350 SPT = 0.9641
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.012 | E = -107.6541220013 | DW = 8.06e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.754 | E = -107.6541223288 | DE = -3.27e-07 | DW = 7.52e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.474 | E = -107.6541224307 | DE = -1.02e-07 | DW = 8.06e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.210 | E = -107.6541224307 | DE = 1.14e-12 | DW = 6.96e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.149 | E = -107.6541224313 | DE = -5.22e-10 | DW = 8.65e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.121 | E = -107.6541224313 | DE = 1.65e-12 | DW = 7.88e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.140 | E = -107.6541224313 | DE = -5.68e-14 | DW = 8.65e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.097 | E = -107.6541224313 | DE = -2.47e-12 | DW = 7.25e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 7.676 | E = -107.6541224313 | DE = 0.00e+00 | DW = 8.25e-20
DMRG energy = -107.654122431266543
```

## The `SGB`

Mode

In this section, we try to solve the problem by first transfroming the model into a qubit (spin) model. The code will automatically use Jordan-Wigner transform to change the fermionic operators in the Hamiltonian into spin operators, before constructing the MPO.

To use the `SGB`

mode for ab initio systems, remember to add the `heis_twos=1`

parameter (indicating the 1/2 spin at each site) in `driver.initialize_system`

.

```
[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.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.SGB, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym, heis_twos=1)
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.140864794887064e-13
integral cutoff error = 0.0
mpo terms = 5984
Build MPO | Nsites = 20 | Nterms = 5984 | 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.006
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.009
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 49 SPT = 0.9479 Tmvc = 0.002 T = 0.011
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 251 SPT = 0.9139 Tmvc = 0.002 T = 0.014
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 273 SPT = 0.9456 Tmvc = 0.002 T = 0.014
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 357 SPT = 0.9576 Tmvc = 0.002 T = 0.014
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 563 SPT = 0.9580 Tmvc = 0.002 T = 0.016
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 2316 SPT = 0.8575 Tmvc = 0.002 T = 0.015
Site = 8 / 20 .. Mmpo = 155 DW = 0.00e+00 NNZ = 214 SPT = 0.9890 Tmvc = 0.001 T = 0.008
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 710 SPT = 0.9690 Tmvc = 0.001 T = 0.012
Site = 10 / 20 .. Mmpo = 181 DW = 0.00e+00 NNZ = 398 SPT = 0.9851 Tmvc = 0.001 T = 0.013
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 720 SPT = 0.9777 Tmvc = 0.001 T = 0.013
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 496 SPT = 0.9810 Tmvc = 0.001 T = 0.007
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 412 SPT = 0.9720 Tmvc = 0.000 T = 0.005
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 254 SPT = 0.9670 Tmvc = 0.000 T = 0.006
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 120 SPT = 0.9711 Tmvc = 0.000 T = 0.004
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 124 SPT = 0.9411 Tmvc = 0.000 T = 0.003
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 76 SPT = 0.9026 Tmvc = 0.000 T = 0.003
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 22 SPT = 0.8429 Tmvc = 0.000 T = 0.002
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.177 Tmvc-total = 0.018 MPO bond dimension = 181 MaxDW = 0.00e+00
NNZ = 7388 SIZE = 204350 SPT = 0.9638
Rank = 0 Ttotal = 0.249 MPO method = FastBipartite bond dimension = 181 NNZ = 7388 SIZE = 204350 SPT = 0.9638
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.079 | E = -107.6541211530 | DW = 6.36e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.838 | E = -107.6541223677 | DE = -1.21e-06 | DW = 5.89e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.571 | E = -107.6541224370 | DE = -6.93e-08 | DW = 6.35e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.284 | E = -107.6541224370 | DE = 4.26e-13 | DW = 5.71e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.581 | E = -107.6541224379 | DE = -9.81e-10 | DW = 9.19e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.311 | E = -107.6541224379 | DE = 1.96e-12 | DW = 6.93e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.757 | E = -107.6541224379 | DE = -5.68e-14 | DW = 9.19e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.724 | E = -107.6541224379 | DE = -1.88e-12 | DW = 7.39e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 9.323 | E = -107.6541224379 | DE = -2.84e-14 | DW = 8.26e-20
DMRG energy = -107.654122437940899
```

## Relativistic DMRG

For relativistic DMRG, we use `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.

The `symm_type`

parameter `SymmetryTypes.SGFCPX`

can also be written as `SymmetryTypes.SGF | SymmetryTypes.CPX`

.

```
[10]:
```

```
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 = 1.5501285354427146e-20
mpo terms = 44346
Build MPO | Nsites = 20 | Nterms = 44346 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.004 T = 0.017
Site = 1 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 23 SPT = 0.9087 Tmvc = 0.005 T = 0.025
Site = 2 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 404 SPT = 0.7710 Tmvc = 0.006 T = 0.023
Site = 3 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 592 SPT = 0.8795 Tmvc = 0.006 T = 0.028
Site = 4 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 932 SPT = 0.8768 Tmvc = 0.010 T = 0.035
Site = 5 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 1315 SPT = 0.8870 Tmvc = 0.008 T = 0.030
Site = 6 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 1722 SPT = 0.9024 Tmvc = 0.008 T = 0.033
Site = 7 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 2142 SPT = 0.9181 Tmvc = 0.012 T = 0.038
Site = 8 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 2597 SPT = 0.9315 Tmvc = 0.009 T = 0.038
Site = 9 / 20 .. Mmpo = 252 DW = 0.00e+00 NNZ = 2985 SPT = 0.9444 Tmvc = 0.011 T = 0.040
Site = 10 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 17958 SPT = 0.6654 Tmvc = 0.014 T = 0.068
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 2590 SPT = 0.9317 Tmvc = 0.003 T = 0.020
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 2184 SPT = 0.9165 Tmvc = 0.003 T = 0.016
Site = 13 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 1783 SPT = 0.8989 Tmvc = 0.002 T = 0.019
Site = 14 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 1397 SPT = 0.8800 Tmvc = 0.002 T = 0.013
Site = 15 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 1008 SPT = 0.8668 Tmvc = 0.001 T = 0.010
Site = 16 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 653 SPT = 0.8671 Tmvc = 0.001 T = 0.010
Site = 17 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 494 SPT = 0.7200 Tmvc = 0.001 T = 0.007
Site = 18 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 26 SPT = 0.8968 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.475 Tmvc-total = 0.105 MPO bond dimension = 252 MaxDW = 0.00e+00
NNZ = 40823 SIZE = 323082 SPT = 0.8736
Rank = 0 Ttotal = 0.567 MPO method = FastBipartite bond dimension = 252 NNZ = 40823 SIZE = 323082 SPT = 0.8736
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 23.177 | E = -107.6929206929 | DW = 4.40e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 35.078 | E = -107.6929209450 | DE = -2.52e-07 | DW = 8.35e-10
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 45.366 | E = -107.6929209492 | DE = -4.14e-09 | DW = 2.09e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 53.524 | E = -107.6929209492 | DE = -3.41e-13 | DW = 8.36e-10
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 71.358 | E = -107.6929209492 | DE = -4.26e-13 | DW = 6.57e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 86.911 | E = -107.6929209492 | DE = 2.84e-14 | DW = 1.18e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 103.723 | E = -107.6929209492 | DE = 2.84e-14 | DW = 7.77e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 119.169 | E = -107.6929209492 | DE = -2.84e-14 | DW = 1.19e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 128.099 | E = -107.6929209492 | DE = 2.84e-14 | DW = 8.33e-20
DMRG energy = -107.692920949170755
```

## The `LZ`

Mode

For diatomic molecules, we can set symmetry `dooh`

in `pyscf`

, and then use `itg.lz_symm_adaptation`

to adapt the atomic orbitals for the `LZ`

symmetry before running Hartree-Fock. Then we can use the `LZ`

modes in `block2`

to perform DMRG.

The `LZ`

mode can be combined with `SU2`

, `SZ`

or `SGF`

spin symmetries, and the `SAny`

prefix in `SymmetryTypes`

. To activate the `SAny`

prefix, the `block2`

code needs to be compiled with the `-DUSE_SANY=ON`

option (this option is ON by default in the `pip`

precompiled binaries). Optionally, when `-DUSE_SANY=ON`

, one can also set `-DUSE_SG=OFF -DUSE_SU2SZ=OFF`

to disable the normal `SU2/SZ/SGF`

modes. One can use
`SymmetryTypes.SAnySU2/SymmetryTypes.SAnySZ/SymmetryTypes.SAnySGF`

instead for normal symmetries (with some limitations) when `-DUSE_SG=OFF -DUSE_SU2SZ=OFF`

is used.

With `SU2`

(spin-adapted DMRG):

```
[11]:
```

```
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="dooh", verbose=0)
mol.symm_orb, z_irrep, g_irrep = itg.lz_symm_adaptation(mol)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym_z = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=1, irrep_id=z_irrep)
print(orb_sym_z)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySU2LZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym_z, pg_irrep=0)
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)
```

```
[ 0 0 0 0 -1 1 0 1 -1 0]
integral symmetrize error = 2.0562981879479644e-15
integral cutoff error = 0.0
mpo terms = 1881
Build MPO | Nsites = 10 | Nterms = 1881 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.001 T = 0.031
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 107 SPT = 0.7752 Tmvc = 0.001 T = 0.054
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 189 SPT = 0.9007 Tmvc = 0.001 T = 0.063
Site = 3 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 822 SPT = 0.7776 Tmvc = 0.001 T = 0.063
Site = 4 / 10 .. Mmpo = 88 DW = 0.00e+00 NNZ = 154 SPT = 0.9735 Tmvc = 0.000 T = 0.033
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 318 SPT = 0.9616 Tmvc = 0.000 T = 0.031
Site = 6 / 10 .. Mmpo = 64 DW = 0.00e+00 NNZ = 236 SPT = 0.9608 Tmvc = 0.000 T = 0.023
Site = 7 / 10 .. Mmpo = 40 DW = 0.00e+00 NNZ = 145 SPT = 0.9434 Tmvc = 0.000 T = 0.016
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 49 SPT = 0.9125 Tmvc = 0.000 T = 0.009
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.007
Ttotal = 0.332 Tmvc-total = 0.005 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 2048 SIZE = 29320 SPT = 0.9302
Rank = 0 Ttotal = 0.402 MPO method = FastBipartite bond dimension = 94 NNZ = 2048 SIZE = 29320 SPT = 0.9302
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.679 | E = -107.6541224475 | DW = 2.06e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.030 | E = -107.6541224475 | DE = 1.19e-12 | DW = 2.77e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.383 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.06e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.721 | E = -107.6541224475 | DE = -8.67e-12 | DW = 3.29e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.225 | E = -107.6541224475 | DE = -2.84e-14 | DW = 4.20e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.894 | E = -107.6541224475 | DE = 5.68e-14 | DW = 2.00e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.682 | E = -107.6541224475 | DE = -2.84e-14 | DW = 4.90e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.379 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.69e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.809 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.21e-20
DMRG energy = -107.654122447523761
```

With `SZ`

(non-spin-adapted DMRG):

```
[12]:
```

```
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="dooh", verbose=0)
mol.symm_orb, z_irrep, g_irrep = itg.lz_symm_adaptation(mol)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym_z = itg.get_uhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=1, irrep_id=z_irrep)
print(orb_sym_z)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySZLZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym_z, pg_irrep=0)
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)
```

```
[ 0 0 0 0 -1 1 0 1 -1 0]
integral symmetrize error = 8.030340066584701e-15
integral cutoff error = 0.0
mpo terms = 5231
Build MPO | Nsites = 10 | Nterms = 5231 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 30 SPT = 0.0000 Tmvc = 0.002 T = 0.080
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 271 SPT = 0.8631 Tmvc = 0.001 T = 0.135
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 480 SPT = 0.9339 Tmvc = 0.001 T = 0.136
Site = 3 / 10 .. Mmpo = 124 DW = 0.00e+00 NNZ = 2273 SPT = 0.8334 Tmvc = 0.001 T = 0.135
Site = 4 / 10 .. Mmpo = 171 DW = 0.00e+00 NNZ = 363 SPT = 0.9829 Tmvc = 0.001 T = 0.066
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 811 SPT = 0.9745 Tmvc = 0.001 T = 0.072
Site = 6 / 10 .. Mmpo = 126 DW = 0.00e+00 NNZ = 583 SPT = 0.9751 Tmvc = 0.001 T = 0.045
Site = 7 / 10 .. Mmpo = 82 DW = 0.00e+00 NNZ = 263 SPT = 0.9745 Tmvc = 0.000 T = 0.024
Site = 8 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 182 SPT = 0.9260 Tmvc = 0.000 T = 0.020
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 30 SPT = 0.0000 Tmvc = 0.000 T = 0.009
Ttotal = 0.722 Tmvc-total = 0.009 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 5286 SIZE = 112178 SPT = 0.9529
Rank = 0 Ttotal = 0.773 MPO method = FastBipartite bond dimension = 186 NNZ = 5286 SIZE = 112178 SPT = 0.9529
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.157 | E = -107.6541224475 | DW = 3.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.836 | E = -107.6541224475 | DE = -1.13e-11 | DW = 3.32e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.477 | E = -107.6541224455 | DE = 2.05e-09 | DW = 3.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.168 | E = -107.6541224455 | DE = -2.27e-13 | DW = 3.27e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.974 | E = -107.6541224455 | DE = -1.42e-13 | DW = 3.08e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.811 | E = -107.6541224455 | DE = 2.84e-14 | DW = 1.35e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.270 | E = -107.6541224455 | DE = -2.84e-14 | DW = 3.08e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.758 | E = -107.6541224455 | DE = 2.56e-13 | DW = 1.55e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 8.655 | E = -107.6541224455 | DE = -2.84e-14 | DW = 9.59e-20
DMRG energy = -107.654122445468332
```

With `SGF`

(spin-orbital DMRG):

```
[13]:
```

```
from pyscf.scf.ghf_symm import GHF
from pyscf import symm, lib
from pyscf.scf import hf_symm
import scipy.linalg
import numpy as np
# fix pyscf 2.3.0 bug in ghf_symm for complex orbtials
def ghf_eig(self, h, s, symm_orb=None, irrep_id=None):
if symm_orb is None or irrep_id is None:
mol = self.mol
symm_orb = mol.symm_orb
irrep_id = mol.irrep_id
nirrep = len(symm_orb)
symm_orb = [scipy.linalg.block_diag(c, c) for c in symm_orb]
h = symm.symmetrize_matrix(h, symm_orb)
s = symm.symmetrize_matrix(s, symm_orb)
cs = []
es = []
orbsym = []
for ir in range(nirrep):
e, c = self._eigh(h[ir], s[ir])
cs.append(c)
es.append(e)
orbsym.append([irrep_id[ir]] * e.size)
e = np.hstack(es)
c = hf_symm.so2ao_mo_coeff(symm_orb, cs)
c = lib.tag_array(c, orbsym=np.hstack(orbsym))
return e, c
GHF.eig = ghf_eig
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="dooh", verbose=0)
mol.symm_orb, z_irrep, g_irrep = itg.lz_symm_adaptation(mol)
mf = scf.GHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym_z = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=1, irrep_id=z_irrep)
print(orb_sym_z)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySGFLZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym_z, pg_irrep=0)
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)
```

```
[ 0 0 0 0 0 0 0 0 -1 -1 1 1 0 0 1 1 -1 -1 0 0]
integral symmetrize error = 9.305335353847714e-15
integral cutoff error = 5.457899320407648e-20
mpo terms = 12773
Build MPO | Nsites = 20 | Nterms = 12773 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.002 T = 0.098
Site = 1 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 25 SPT = 0.9008 Tmvc = 0.002 T = 0.170
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 309 SPT = 0.7652 Tmvc = 0.004 T = 0.245
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 357 SPT = 0.8775 Tmvc = 0.003 T = 0.264
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 514 SPT = 0.8977 Tmvc = 0.003 T = 0.280
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 667 SPT = 0.9208 Tmvc = 0.003 T = 0.308
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 2106 SPT = 0.8430 Tmvc = 0.003 T = 0.299
Site = 7 / 20 .. Mmpo = 118 DW = 0.00e+00 NNZ = 3960 SPT = 0.7399 Tmvc = 0.003 T = 0.249
Site = 8 / 20 .. Mmpo = 153 DW = 0.00e+00 NNZ = 250 SPT = 0.9862 Tmvc = 0.001 T = 0.127
Site = 9 / 20 .. Mmpo = 164 DW = 0.00e+00 NNZ = 690 SPT = 0.9725 Tmvc = 0.001 T = 0.134
Site = 10 / 20 .. Mmpo = 185 DW = 0.00e+00 NNZ = 1223 SPT = 0.9597 Tmvc = 0.002 T = 0.131
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 910 SPT = 0.9724 Tmvc = 0.001 T = 0.101
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 814 SPT = 0.9689 Tmvc = 0.001 T = 0.077
Site = 13 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 632 SPT = 0.9642 Tmvc = 0.001 T = 0.052
Site = 14 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 383 SPT = 0.9671 Tmvc = 0.001 T = 0.040
Site = 15 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 249 SPT = 0.9671 Tmvc = 0.000 T = 0.038
Site = 16 / 20 .. Mmpo = 57 DW = 0.00e+00 NNZ = 382 SPT = 0.9141 Tmvc = 0.000 T = 0.022
Site = 17 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 91 SPT = 0.9430 Tmvc = 0.000 T = 0.009
Site = 18 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 28 SPT = 0.8889 Tmvc = 0.000 T = 0.011
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Ttotal = 2.660 Tmvc-total = 0.031 MPO bond dimension = 185 MaxDW = 0.00e+00
NNZ = 13608 SIZE = 222306 SPT = 0.9388
Rank = 0 Ttotal = 2.764 MPO method = FastBipartite bond dimension = 185 NNZ = 13608 SIZE = 222306 SPT = 0.9388
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.655 | E = -107.6541219864 | DW = 1.34e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.673 | E = -107.6541222928 | DE = -3.06e-07 | DW = 4.92e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 5.224 | E = -107.6541224418 | DE = -1.49e-07 | DW = 1.33e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 6.259 | E = -107.6541224418 | DE = -1.02e-12 | DW = 4.94e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.849 | E = -107.6541224449 | DE = -3.13e-09 | DW = 6.76e-12
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 9.282 | E = -107.6541224455 | DE = -5.68e-10 | DW = 5.49e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 10.873 | E = -107.6541224455 | DE = 5.68e-14 | DW = 6.76e-12
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 12.297 | E = -107.6541224455 | DE = 2.84e-14 | DW = 6.04e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 13.165 | E = -107.6541224455 | DE = 2.84e-14 | DW = 7.22e-20
DMRG energy = -107.654122445469270
```

## Expectation and N-Particle Density Matrices

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.

```
[14]:
```

```
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.557 | 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.833 | E = -106.9391328597 | DE = -3.24e-12 | DW = 9.75e-19
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.130 | E = -106.9391328597 | DE = 2.84e-14 | DW = 3.38e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.421 | E = -106.9391328597 | DE = 1.71e-13 | DW = 1.52e-18
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.755 | 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.063 | E = -106.9391328597 | DE = 5.68e-14 | DW = 1.40e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.402 | E = -106.9391328597 | DE = 0.00e+00 | DW = 1.12e-16
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.709 | E = -106.9391328597 | DE = 0.00e+00 | DW = 1.38e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.933 | E = -106.9391328597 | DE = 0.00e+00 | DW = 3.39e-20
DMRG energy = -106.939132859666600
Norm = 1.000000000000000
Energy expectation = -106.939132859666614
<S^2> expectation = 2.000000000000000
```

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

```
[15]:
```

```
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.999995824361892
```

We can then verify this number using 1PDM:

```
[16]:
```

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

```
N0 expectation from 1pdm = 1.999995824361892
```

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`

.

```
[17]:
```

```
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 = 4.622993128790701e-06
```

## Extract CSF and Determinant Coefficients

We can extract CSF (or determinant) coefficients from the spin-adapted MPS (or the non-spin-adapted MPS). The algorithm can compute all CSF or determinant with the absolute value of the coefficient above a threshold (called `cutoff`

). The square of coefficient is the probability (weight) of the CSF or determinant.

Extracting CSF coefficients from spin-adapted MPS in the `SU2`

mode:

```
[18]:
```

```
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)
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)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.639 | 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.016 | E = -107.6541224475 | DE = -3.38e-12 | DW = 3.81e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.407 | E = -107.6541224475 | DE = 5.68e-14 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.777 | E = -107.6541224475 | DE = -5.68e-13 | DW = 4.69e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.209 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.73e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.604 | E = -107.6541224475 | DE = 2.84e-14 | DW = 5.35e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.011 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.66e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.376 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.98e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.647 | E = -107.6541224475 | DE = 2.84e-14 | DW = 5.56e-20
DMRG energy = -107.654122447524614
dtrie finished 0.01697694599999977
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368811359554
CSF 0 2222222000 = 0.957506528669401
CSF 1 2222202200 = -0.131287867807394
CSF 2 2222022020 = -0.131287794595557
CSF 3 2222+-2+-0 = 0.118594301088608
CSF 4 2222++2--0 = -0.084968224982022
CSF 5 2220222020 = -0.054710635444627
CSF 6 2220222200 = -0.054710465272207
CSF 7 2222+2+0-- = 0.053881241407894
CSF 8 22222++-0- = 0.053881215166769
```

Extracting determinant coefficients from spin-adapted MPS in the `SZ`

mode:

```
[19]:
```

```
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=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)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.373 | 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 = 2.089 | E = -107.6541224475 | DE = -2.07e-12 | DW = 5.05e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.609 | E = -107.6541224475 | DE = -1.34e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.051 | E = -107.6541224475 | DE = 1.28e-12 | DW = 5.21e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.572 | 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 = 4.085 | E = -107.6541224475 | DE = 8.81e-13 | DW = 1.37e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.625 | 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 = 5.142 | E = -107.6541224475 | DE = -1.28e-12 | DW = 1.52e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 5.505 | E = -107.6541224475 | DE = 0.00e+00 | DW = 8.34e-20
DMRG energy = -107.654122447524443
dtrie finished 0.02674672699998837
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331638354918
DET 0 2222222000 = 0.957506568620143
DET 1 2222022020 = -0.131287835040785
DET 2 2222202200 = -0.131287814797523
DET 3 2222ab2ab0 = 0.083825279692261
DET 4 2222ba2ba0 = 0.083825279691480
DET 5 2220222200 = -0.054710476471258
DET 6 2220222020 = -0.054710439530678
```

## Construct MPS from CSFs or Determinants

If we know important CSFs or determinants in the state and their coefficients, we can also use this information to construct MPS, and this can be used as an initial guess and further optimized. Note that this initial guess can generate very good initial energies, but if the given CSFs have many doubly occupied and empty orbitals, the MPS will very likely optimize to a local minima.

In the `SU2`

mode:

```
[20]:
```

```
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)
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)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
mps = driver.get_mps_from_csf_coefficients(csfs, coeffs, tag="CMPS", dot=2)
impo = driver.get_identity_mpo()
print(driver.expectation(mps, impo, mps))
print(driver.expectation(mps, mpo, mps) / driver.expectation(mps, impo, mps))
energy = driver.dmrg(mpo, mps, n_sweeps=5, bond_dims=[500] * 5, noises=[1E-5],
thrds=[1E-10] * 5, iprint=1)
print('Ground state energy = %20.15f' % energy)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.496 | 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 = 0.733 | E = -107.6541224475 | DE = -1.24e-11 | DW = 4.89e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.990 | 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 = 1.208 | E = -107.6541224475 | DE = -1.36e-12 | DW = 3.97e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.488 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.48e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.713 | E = -107.6541224475 | DE = -2.84e-14 | DW = 6.11e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.974 | E = -107.6541224475 | DE = 5.68e-14 | DW = 2.41e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.209 | E = -107.6541224475 | DE = -2.84e-14 | DW = 6.14e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.482 | E = -107.6541224475 | DE = 5.68e-14 | DW = 2.91e-20
DMRG energy = -107.654122447524358
dtrie finished 0.016507582999963688
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368809657036
CSF 0 2222222000 = 0.957506521875127
CSF 1 2222202200 = -0.131287813613834
CSF 2 2222022020 = -0.131287676854091
CSF 3 2222+-2+-0 = 0.118594392287595
CSF 4 2222++2--0 = -0.084968279487231
CSF 5 2220222020 = -0.054710673444202
CSF 6 2220222200 = -0.054710596031681
CSF 7 2222+2+0-- = 0.053881289057939
CSF 8 22222++-0- = 0.053881233355398
0.9843688096570361
-107.6155366353171
Sweep = 0 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.054 | E = -107.6285084535 | DW = 3.52e-21
Sweep = 1 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.109 | E = -107.6288109091 | DE = -3.02e-04 | DW = 9.04e-21
Sweep = 2 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.164 | E = -107.6289019710 | DE = -9.11e-05 | DW = 4.57e-21
Sweep = 3 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.217 | E = -107.6289019713 | DE = -3.35e-10 | DW = 2.07e-20
Ground state energy = -107.628901971339317
```

In the `SZ`

mode:

```
[21]:
```

```
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=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)
dets, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
mps = driver.get_mps_from_csf_coefficients(dets, coeffs, tag="CMPS", dot=2)
impo = driver.get_identity_mpo()
print(driver.expectation(mps, impo, mps))
print(driver.expectation(mps, mpo, mps) / driver.expectation(mps, impo, mps))
energy = driver.dmrg(mpo, mps, n_sweeps=5, bond_dims=[500] * 5, noises=[1E-5],
thrds=[1E-10] * 5, iprint=1)
print('Ground state energy = %20.15f' % energy)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.809 | 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.294 | E = -107.6541224475 | DE = -5.80e-12 | DW = 5.17e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.713 | E = -107.6541224475 | DE = -9.38e-13 | DW = 4.15e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.170 | E = -107.6541224475 | DE = 1.36e-12 | DW = 5.37e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.687 | E = -107.6541224475 | DE = -1.36e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.212 | E = -107.6541224475 | DE = 1.05e-12 | DW = 1.85e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.737 | 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.252 | E = -107.6541224475 | DE = -1.14e-12 | DW = 1.77e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.598 | E = -107.6541224475 | DE = 0.00e+00 | DW = 8.68e-20
DMRG energy = -107.654122447524529
dtrie finished 0.03569589099998893
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331637808218
DET 0 2222222000 = -0.957506566009713
DET 1 2222022020 = 0.131287836810967
DET 2 2222202200 = 0.131287818787713
DET 3 2222ab2ab0 = -0.083825284176507
DET 4 2222ba2ba0 = -0.083825284073794
DET 5 2220222200 = 0.054710477807739
DET 6 2220222020 = 0.054710451475846
0.971331637808218
-107.59447325662364
Sweep = 0 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.047 | E = -107.6101211360 | DW = 6.22e-21
Sweep = 1 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.104 | E = -107.6102787809 | DE = -1.58e-04 | DW = 8.84e-21
Sweep = 2 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.163 | E = -107.6104759466 | DE = -1.97e-04 | DW = 1.22e-20
Sweep = 3 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.218 | E = -107.6104759470 | DE = -4.45e-10 | DW = 1.18e-20
Ground state energy = -107.610475947024796
```

## MPS Initial Guess from Occupancies

We can also construct the MPS using an estimate of the occupancy information at each site from a cheaper method (such as CCSD). We can use this information to construct the initial quantum number distribution in the MPS.

Note that this initial guess can generate very good energies in the first few sweeps, but if the given occupancies have many doubly occupied and empty orbitals, the MPS will very likely optimize to a local minima. One can shift the occupancies into the equal probability occupancies (for example, setting `bias=0.4`

in the example below) to randomize the initial guess.

In the `SU2`

mode:

```
[22]:
```

```
from pyscf import gto, scf, cc
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)
bias = 0.1 # make it more random
occs = np.diag(cc.CCSD(mf).run().make_rdm1())
occs = occs + bias * (occs < 1) - bias * (occs > 1)
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, occs=occs, 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)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.044 | E = -107.5033999317 | DW = 2.59e-21
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.135 | E = -107.5830591524 | DE = -7.97e-02 | DW = 3.51e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.217 | E = -107.5836402922 | DE = -5.81e-04 | DW = 9.17e-21
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.294 | E = -107.5867932578 | DE = -3.15e-03 | DW = 1.12e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.371 | E = -107.5868061724 | DE = -1.29e-05 | DW = 1.17e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.443 | E = -107.5868152464 | DE = -9.07e-06 | DW = 1.23e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.517 | E = -107.5868152464 | DE = -5.68e-14 | DW = 1.45e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.589 | E = -107.5868152464 | DE = 2.84e-14 | DW = 8.04e-21
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.655 | E = -107.5868152464 | DE = 2.84e-14 | DW = 1.69e-20
DMRG energy = -107.586815246376730
```

In the `SZ`

mode:

```
[23]:
```

```
from pyscf import gto, scf, cc
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)
bias = 0.1 # make it more random
occs = np.diag(np.sum(cc.UCCSD(mf).run().make_rdm1(), axis=0))
occs = occs + bias * (occs < 1) - bias * (occs > 1)
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=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, occs=occs, 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)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.097 | E = -107.5033999316 | DW = 3.37e-21
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.267 | E = -107.5830606948 | DE = -7.97e-02 | DW = 7.61e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.376 | E = -107.5836402307 | DE = -5.80e-04 | DW = 1.27e-20
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.505 | E = -107.5867932582 | DE = -3.15e-03 | DW = 1.21e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.624 | E = -107.5868061673 | DE = -1.29e-05 | DW = 1.84e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.847 | E = -107.5868152464 | DE = -9.08e-06 | DW = 1.03e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.008 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.88e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.132 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.76e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 1.218 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.33e-20
DMRG energy = -107.586815246376730
```

In the `SGF`

mode:

```
[24]:
```

```
from pyscf import gto, scf, cc
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)
bias = 0.25 # make it more random
occs = np.sum(cc.GCCSD(mf).run().make_rdm1(), axis=0)
occs = occs + bias * (occs < 0.5) - bias * (occs > 0.5)
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=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, occs=occs, 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)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.853 | E = -107.6538893531 | DW = 2.24e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.548 | E = -107.6539226742 | DE = -3.33e-05 | DW = 8.93e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.168 | E = -107.6539462697 | DE = -2.36e-05 | DW = 5.04e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.808 | E = -107.6539506829 | DE = -4.41e-06 | DW = 1.95e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.596 | E = -107.6539506829 | DE = 3.24e-12 | DW = 1.84e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.386 | E = -107.6539506829 | DE = -4.12e-12 | DW = 5.67e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.236 | E = -107.6539506829 | DE = -2.84e-14 | DW = 1.84e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.014 | E = -107.6539506829 | DE = 2.05e-12 | DW = 4.05e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.497 | E = -107.6539506829 | DE = 0.00e+00 | DW = 7.66e-20
DMRG energy = -107.653950682884570
```

## Change from SU2 MPS to SZ MPS

We can also transform the spin-adapted MPS generated in the `SU2`

mode to the non-spin-adapted MPS, which can be used in the `SZ`

mode. After obtaining the non-spin-adapted MPS, one need to redo `driver.initialize_system`

, `driver.get_qc_mpo`

, etc. to make sure every object is now represented in the `SZ`

mode, then you can operate on the `SZ`

non-spin-adapted MPS.

In the following example, we first compute the spin-adapted MPS, then translate it into the non-spin-adapted MPS to extract the determinant coefficients.

```
[25]:
```

```
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)
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)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
zket = driver.mps_change_to_sz(ket, "ZKET")
driver.symm_type = SymmetryTypes.SZ
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)
impo = driver.get_identity_mpo()
print(driver.expectation(zket, mpo, zket) / driver.expectation(zket, impo, zket))
csfs, vals = driver.get_csf_coefficients(zket, cutoff=0.05, iprint=1)
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.902 | 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.281 | E = -107.6541224475 | DE = -2.87e-12 | DW = 4.49e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.682 | 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.068 | E = -107.6541224475 | DE = -5.12e-13 | DW = 5.07e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.610 | E = -107.6541224475 | DE = -5.68e-14 | DW = 2.62e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.019 | E = -107.6541224475 | DE = 2.84e-14 | DW = 4.59e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.435 | E = -107.6541224475 | DE = 2.84e-14 | DW = 2.01e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.808 | E = -107.6541224475 | DE = -5.68e-14 | DW = 3.69e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.095 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.68e-20
DMRG energy = -107.654122447524671
dtrie finished 0.026360811000017748
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368806039259
CSF 0 2222222000 = -0.957506495595307
CSF 1 2222022020 = 0.131287987564377
CSF 2 2222202200 = 0.131287965813458
CSF 3 2222+-2+-0 = -0.118594461577159
CSF 4 2222++2--0 = 0.084967942537817
CSF 5 2220222020 = 0.054710523921704
CSF 6 2220222200 = 0.054710517406346
CSF 7 2222+2+0-- = -0.053881222618422
CSF 8 22222++-0- = -0.053881215803966
-107.65412244752456
dtrie finished 0.0680388100000755
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331619872548
DET 0 2222222000 = -0.957506495595307
DET 1 2222022020 = 0.131287987564377
DET 2 2222202200 = 0.131287965813458
DET 3 2222ab2ab0 = -0.083825363036928
DET 4 2222ba2ba0 = -0.083825363036928
DET 5 2220222020 = 0.054710523921704
DET 6 2220222200 = 0.054710517406346
```

## Change between Real and Complex MPS

We can also change between the MPS with complex numbers and the MPS with real numbers. For complex MPS to real MPS, the imaginary part will be discarded (and the norm of the transformed MPS may decrease). This may be useful when you do a ground state calculation in the real domain and then do the real time evolution in the complex domain.

From real to complex:

```
[26]:
```

```
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)
impo = driver.get_identity_mpo()
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)
print(driver.expectation(ket, impo, ket))
print(driver.expectation(ket, mpo, ket) / driver.expectation(ket, impo, ket))
zket = driver.mps_change_complex(ket, "ZKET")
driver.symm_type = driver.symm_type ^ SymmetryTypes.CPX
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
impo = driver.get_identity_mpo()
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, integral_cutoff=1E-8, iprint=1)
print(driver.expectation(zket, impo, zket))
print(driver.expectation(zket, mpo, zket) / driver.expectation(zket, impo, zket))
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.804 | 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.144 | E = -107.6541224475 | DE = -3.95e-12 | DW = 4.49e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.379 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.601 | E = -107.6541224475 | DE = -7.67e-13 | DW = 5.89e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.868 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.31e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.110 | E = -107.6541224475 | DE = -2.84e-14 | DW = 6.59e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.366 | E = -107.6541224475 | DE = 0.00e+00 | DW = 8.80e-21
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.606 | E = -107.6541224475 | DE = -5.68e-14 | DW = 6.06e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.780 | E = -107.6541224475 | DE = 2.84e-14 | DW = 4.42e-20
DMRG energy = -107.654122447524671
1.0
-107.65412244752456
integral symmetrize error = 1.76819336753967e-14
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 13 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.001 T = 0.007
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.001 T = 0.005
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.006
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.000 T = 0.006
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.007
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.005
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.008
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.003
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.004
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.054 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.108 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
(0.9999999999999999+0j)
(-107.65412244752457+0j)
```

From complex to real:

```
[27]:
```

```
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 | SymmetryTypes.CPX, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
impo = driver.get_identity_mpo()
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)
print(driver.expectation(ket, impo, ket))
print(driver.expectation(ket, mpo, ket) / driver.expectation(ket, impo, ket))
rket = driver.mps_change_complex(ket, "rket")
driver.symm_type = driver.symm_type ^ SymmetryTypes.CPX
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
impo = driver.get_identity_mpo()
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, integral_cutoff=1E-8, iprint=1)
print(driver.expectation(rket, impo, rket))
print(driver.expectation(rket, mpo, rket) / driver.expectation(rket, impo, rket))
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.189 | 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 = 2.240 | E = -107.6541224475 | DE = -1.41e-11 | DW = 2.95e-19
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.105 | 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 = 3.655 | E = -107.6541224475 | DE = -2.05e-12 | DW = 1.13e-19
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.299 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.31e-19
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.784 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.06e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.305 | E = -107.6541224475 | DE = 5.68e-14 | DW = 1.24e-19
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.770 | E = -107.6541224475 | DE = -5.68e-14 | DW = 1.04e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.086 | E = -107.6541224475 | DE = 0.00e+00 | DW = 4.25e-20
DMRG energy = -107.654122447523960
(1.0000000000000002-7.905175135209388e-18j)
(-107.6541224475238-8.610748230937272e-16j)
integral symmetrize error = 1.6922854083945986e-14
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | 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.004
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.001 T = 0.004
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.004
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.000 T = 0.005
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.004
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.004
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.003
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.003
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.002
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.036 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.064 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
0.9626534962131837
-107.65412244752437
```

## MPS Bipartite Entanglement

We can get the the bipartite entanglement \(S_k=-\sum_i \Lambda_k^2 \log \Lambda_k^2\) at each virtual bond (at site \(k\)) in MPS in the `SZ`

mode, where \(\Lambda_k\) are singular values in the bond at site \(k\).

```
[28]:
```

```
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.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=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)
bip_ent = driver.get_bipartite_entanglement()
import matplotlib.pyplot as plt
plt.plot(np.arange(len(bip_ent)), bip_ent, linestyle='-', marker='o',
mfc='white', mec="#7FB685", color="#7FB685")
plt.xlabel("site index $k$")
plt.ylabel("bipartite entanglement $S_k$")
plt.show()
```

```
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.769 | 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.254 | E = -107.6541224475 | DE = -1.43e-11 | DW = 5.08e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.670 | E = -107.6541224475 | DE = -5.12e-13 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.130 | E = -107.6541224475 | DE = 1.31e-12 | DW = 5.21e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.655 | E = -107.6541224475 | DE = -1.14e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.189 | E = -107.6541224475 | DE = 9.38e-13 | DW = 1.75e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.722 | 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.245 | E = -107.6541224475 | DE = -1.22e-12 | DW = 1.54e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.595 | E = -107.6541224475 | DE = 2.84e-14 | DW = 7.71e-20
DMRG energy = -107.654122447524443
```

## Orbital Entropy and Mutual Information

For the optimized MPS in the `SZ`

mode, we can compute the 1- and 2- orbital density matrices and mutual information for pairs of orbitals.

```
[29]:
```

```
odm1 = driver.get_orbital_entropies(ket, orb_type=1)
odm2 = driver.get_orbital_entropies(ket, orb_type=2)
minfo = 0.5 * (odm1[:, None] + odm1[None, :] - odm2) * (1 - np.identity(len(odm1)))
import matplotlib.pyplot as plt
plt.matshow(minfo)
```

```
[29]:
```

```
<matplotlib.image.AxesImage at 0x7d0b6521bd30>
```

## Excited States

To obtain the excited states and their energies, we can perform DMRG for a state-averged MPS, optionally followed by a state-specific refinement.

```
[30]:
```

```
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.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=0)
ket = driver.get_random_mps(tag="KET", bond_dim=100, nroots=3)
energies = driver.dmrg(mpo, ket, n_sweeps=10, bond_dims=[100], noises=[1e-5] * 4 + [0],
thrds=[1e-10] * 8, iprint=1)
print('State-averaged MPS energies = [%s]' % " ".join("%20.15f" % x for x in energies))
kets = [driver.split_mps(ket, ir, tag="KET-%d" % ir) for ir in range(ket.nroots)]
for ir in range(ket.nroots):
energy = driver.dmrg(mpo, kets[ir], n_sweeps=10, bond_dims=[200], noises=[1e-5] * 4 + [0],
thrds=[1e-10] * 8, iprint=0, proj_weights=[5.0] * ir, proj_mpss=kets[:ir])
print('State-specific MPS E[%d] = %20.15f' % (ir, energy))
```

```
Sweep = 0 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.037 | E[ 3] = -107.6541193407 -107.0314407783 -106.9595996333 | DW = 8.51e-05
Sweep = 1 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.220 | E[ 3] = -107.6541193407 -107.0314407783 -106.9595996332 | DE = 7.11e-12 | DW = 1.62e-04
Sweep = 2 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.616 | E[ 3] = -107.6541138874 -107.0314486179 -106.9594564026 | DE = 1.43e-04 | DW = 7.20e-05
Sweep = 3 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.853 | E[ 3] = -107.6541138874 -107.0314486178 -106.9594564026 | DE = -3.97e-11 | DW = 1.47e-04
Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10
Time elapsed = 6.621 | E[ 3] = -107.6540972369 -107.0314487271 -106.9594391421 | DE = 1.73e-05 | DW = 6.81e-05
Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10
Time elapsed = 7.410 | E[ 3] = -107.6540972369 -107.0314487271 -106.9594391422 | DE = -1.72e-11 | DW = 1.43e-04
State-averaged MPS energies = [-107.654097236905088 -107.031448727142006 -106.959439142163390]
State-specific MPS E[0] = -107.654122447496988
State-specific MPS E[1] = -107.031449471612873
State-specific MPS E[2] = -106.959626154678844
```