Quantum Chemistry Hamiltonians

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

Introduction

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

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

  1. 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).

  2. 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).

  3. 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).

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

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

[2]:
import numpy as np
from pyblock2._pyscf.ao2mo import integrals as itg
from pyblock2.driver.core import DMRGDriver, SymmetryTypes

Spin-Restricted Integrals

Here we use get_rhf_integrals function to get the integrals. Note that in order to do DMRG in a CASCI space, one can set the ncore (number of core orbitals) and ncas (number of active orbitals) parameters in get_*_integrals. ncas=None will include all orbitals in DMRG.

For medium to large scale DMRG calculations, it is highly recommended to use a scratch space with high IO speed rather than the ./tmp used in the following example. One also needs to set a suitable stack_mem in the DMRGDriver constructor to set the memory used for storing renormalized operators (in bytes). The default is stack_mem=int(1024**3) (1 GB). For medium scale calculations 10 to 30 GB might be required.

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

[3]:
from pyscf import gto, scf

mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
    ncore=0, ncas=None, g2e_symm=8)

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)

bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
    thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)

pdm1 = driver.get_1pdm(ket)
pdm2 = driver.get_2pdm(ket).transpose(0, 3, 1, 2)
print('Energy from pdms = %20.15f' % (np.einsum('ij,ij->', pdm1, h1e)
    + 0.5 * np.einsum('ijkl,ijkl->', pdm2, driver.unpack_g2e(g2e)) + ecore))
integral symmetrize error =  6.220297364062989e-14
integral cutoff error =  0.0
mpo terms =        972

Build MPO | Nsites =    10 | Nterms =        972 | Algorithm = FastBIP | Cutoff = 1.00e-20
 Site =     0 /    10 .. Mmpo =    13 DW = 0.00e+00 NNZ =       13 SPT = 0.0000 Tmvc = 0.000 T = 0.009
 Site =     1 /    10 .. Mmpo =    34 DW = 0.00e+00 NNZ =       59 SPT = 0.8665 Tmvc = 0.001 T = 0.009
 Site =     2 /    10 .. Mmpo =    56 DW = 0.00e+00 NNZ =      117 SPT = 0.9386 Tmvc = 0.000 T = 0.009
 Site =     3 /    10 .. Mmpo =    70 DW = 0.00e+00 NNZ =      363 SPT = 0.9074 Tmvc = 0.001 T = 0.009
 Site =     4 /    10 .. Mmpo =    80 DW = 0.00e+00 NNZ =      217 SPT = 0.9613 Tmvc = 0.000 T = 0.007
 Site =     5 /    10 .. Mmpo =    94 DW = 0.00e+00 NNZ =      201 SPT = 0.9733 Tmvc = 0.000 T = 0.008
 Site =     6 /    10 .. Mmpo =    54 DW = 0.00e+00 NNZ =      169 SPT = 0.9667 Tmvc = 0.000 T = 0.007
 Site =     7 /    10 .. Mmpo =    30 DW = 0.00e+00 NNZ =       73 SPT = 0.9549 Tmvc = 0.000 T = 0.006
 Site =     8 /    10 .. Mmpo =    13 DW = 0.00e+00 NNZ =       41 SPT = 0.8949 Tmvc = 0.000 T = 0.008
 Site =     9 /    10 .. Mmpo =     1 DW = 0.00e+00 NNZ =       13 SPT = 0.0000 Tmvc = 0.000 T = 0.006
Ttotal =      0.078 Tmvc-total = 0.004 MPO bond dimension =    94 MaxDW = 0.00e+00
NNZ =         1266 SIZE =        26498 SPT = 0.9522

Rank =     0 Ttotal =      0.154 MPO method = FastBipartite bond dimension =      94 NNZ =         1266 SIZE =        26498 SPT = 0.9522

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.107 | E =    -107.6541224475 | DW = 1.87e-10

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.640 | E =    -107.6541224475 | DE = -1.07e-11 | DW = 4.78e-20

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      2.154 | E =    -107.6541224475 | DE = 2.84e-14 | DW = 1.87e-10

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      2.552 | E =    -107.6541224475 | DE = -1.17e-12 | DW = 5.48e-20

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      3.113 | E =    -107.6541224475 | DE = -2.84e-14 | DW = 2.86e-20

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      3.668 | E =    -107.6541224475 | DE = 0.00e+00 | DW = 4.95e-20

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.254 | E =    -107.6541224475 | DE = -2.84e-14 | DW = 3.20e-20

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.842 | E =    -107.6541224475 | DE = 8.53e-14 | DW = 6.54e-20

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      5.168 | E =    -107.6541224475 | DE = 0.00e+00 | DW = 2.97e-20

DMRG energy = -107.654122447523932
Energy from pdms = -107.654122447523847

We can also run non-spin-adapted DMRG (SZ mode) using the restricted integrals.

[4]:
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
    thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error =  9.016818189452533e-14
integral cutoff error =  0.0
mpo terms =       2778

Build MPO | Nsites =    10 | Nterms =       2778 | Algorithm = FastBIP | Cutoff = 1.00e-20
 Site =     0 /    10 .. Mmpo =    26 DW = 0.00e+00 NNZ =       26 SPT = 0.0000 Tmvc = 0.001 T = 0.009
 Site =     1 /    10 .. Mmpo =    66 DW = 0.00e+00 NNZ =      143 SPT = 0.9167 Tmvc = 0.001 T = 0.007
 Site =     2 /    10 .. Mmpo =   110 DW = 0.00e+00 NNZ =      283 SPT = 0.9610 Tmvc = 0.001 T = 0.015
 Site =     3 /    10 .. Mmpo =   138 DW = 0.00e+00 NNZ =     1023 SPT = 0.9326 Tmvc = 0.002 T = 0.015
 Site =     4 /    10 .. Mmpo =   158 DW = 0.00e+00 NNZ =      535 SPT = 0.9755 Tmvc = 0.001 T = 0.011
 Site =     5 /    10 .. Mmpo =   186 DW = 0.00e+00 NNZ =      463 SPT = 0.9842 Tmvc = 0.001 T = 0.011
 Site =     6 /    10 .. Mmpo =   106 DW = 0.00e+00 NNZ =      415 SPT = 0.9790 Tmvc = 0.000 T = 0.006
 Site =     7 /    10 .. Mmpo =    58 DW = 0.00e+00 NNZ =      163 SPT = 0.9735 Tmvc = 0.000 T = 0.004
 Site =     8 /    10 .. Mmpo =    26 DW = 0.00e+00 NNZ =       87 SPT = 0.9423 Tmvc = 0.000 T = 0.004
 Site =     9 /    10 .. Mmpo =     1 DW = 0.00e+00 NNZ =       26 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Ttotal =      0.087 Tmvc-total = 0.008 MPO bond dimension =   186 MaxDW = 0.00e+00
NNZ =         3164 SIZE =       102772 SPT = 0.9692

Rank =     0 Ttotal =      0.130 MPO method = FastBipartite bond dimension =     186 NNZ =         3164 SIZE =       102772 SPT = 0.9692

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.831 | E =    -107.6541224475 | DW = 4.14e-08

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.323 | E =    -107.6541224475 | DE = -4.46e-12 | DW = 5.27e-09

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.792 | E =    -107.6541224475 | DE = -1.02e-12 | DW = 4.14e-08

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      2.389 | E =    -107.6541224475 | DE = 1.42e-12 | DW = 5.36e-09

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      3.364 | E =    -107.6541224475 | DE = -1.11e-12 | DW = 3.60e-11

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.355 | E =    -107.6541224475 | DE = 9.09e-13 | DW = 1.51e-19

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      5.296 | E =    -107.6541224475 | DE = 2.84e-14 | DW = 3.60e-11

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      6.108 | E =    -107.6541224475 | DE = -1.14e-12 | DW = 1.68e-19

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      6.506 | E =    -107.6541224475 | DE = -2.84e-14 | DW = 8.66e-20

DMRG energy = -107.654122447524529

We can also run DMRG in spin orbitals (SGF mode) using the restricted integrals, which will be much slower (for more realistic systems).

[5]:
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)

driver.n_sites = ncas
g2e = driver.unpack_g2e(g2e)
orb_sym = [orb_sym[i // 2] for i in range(len(orb_sym) * 2)]
n_sites = ncas * 2

driver.initialize_system(n_sites=n_sites, n_elec=n_elec, spin=spin, orb_sym=orb_sym)

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
    thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error =  9.016818189452533e-14
integral cutoff error =  0.0
mpo terms =       2438

Build MPO | Nsites =    20 | Nterms =       2438 | Algorithm = FastBIP | Cutoff = 1.00e-20
 Site =     0 /    20 .. Mmpo =     7 DW = 0.00e+00 NNZ =        7 SPT = 0.0000 Tmvc = 0.000 T = 0.010
 Site =     1 /    20 .. Mmpo =    20 DW = 0.00e+00 NNZ =       19 SPT = 0.8643 Tmvc = 0.001 T = 0.009
 Site =     2 /    20 .. Mmpo =    45 DW = 0.00e+00 NNZ =       45 SPT = 0.9500 Tmvc = 0.002 T = 0.008
 Site =     3 /    20 .. Mmpo =    62 DW = 0.00e+00 NNZ =      131 SPT = 0.9530 Tmvc = 0.001 T = 0.010
 Site =     4 /    20 .. Mmpo =    81 DW = 0.00e+00 NNZ =      159 SPT = 0.9683 Tmvc = 0.001 T = 0.008
 Site =     5 /    20 .. Mmpo =   104 DW = 0.00e+00 NNZ =      203 SPT = 0.9759 Tmvc = 0.001 T = 0.007
 Site =     6 /    20 .. Mmpo =   125 DW = 0.00e+00 NNZ =      265 SPT = 0.9796 Tmvc = 0.002 T = 0.008
 Site =     7 /    20 .. Mmpo =   126 DW = 0.00e+00 NNZ =      974 SPT = 0.9382 Tmvc = 0.001 T = 0.009
 Site =     8 /    20 .. Mmpo =   151 DW = 0.00e+00 NNZ =      188 SPT = 0.9901 Tmvc = 0.001 T = 0.006
 Site =     9 /    20 .. Mmpo =   148 DW = 0.00e+00 NNZ =      344 SPT = 0.9846 Tmvc = 0.001 T = 0.010
 Site =    10 /    20 .. Mmpo =   177 DW = 0.00e+00 NNZ =      246 SPT = 0.9906 Tmvc = 0.001 T = 0.006
 Site =    11 /    20 .. Mmpo =   178 DW = 0.00e+00 NNZ =      402 SPT = 0.9872 Tmvc = 0.001 T = 0.006
 Site =    12 /    20 .. Mmpo =   147 DW = 0.00e+00 NNZ =      306 SPT = 0.9883 Tmvc = 0.000 T = 0.006
 Site =    13 /    20 .. Mmpo =   100 DW = 0.00e+00 NNZ =      242 SPT = 0.9835 Tmvc = 0.000 T = 0.009
 Site =    14 /    20 .. Mmpo =    77 DW = 0.00e+00 NNZ =      150 SPT = 0.9805 Tmvc = 0.000 T = 0.006
 Site =    15 /    20 .. Mmpo =    54 DW = 0.00e+00 NNZ =       94 SPT = 0.9774 Tmvc = 0.000 T = 0.004
 Site =    16 /    20 .. Mmpo =    39 DW = 0.00e+00 NNZ =       82 SPT = 0.9611 Tmvc = 0.000 T = 0.008
 Site =    17 /    20 .. Mmpo =    20 DW = 0.00e+00 NNZ =       50 SPT = 0.9359 Tmvc = 0.000 T = 0.004
 Site =    18 /    20 .. Mmpo =     7 DW = 0.00e+00 NNZ =       20 SPT = 0.8571 Tmvc = 0.000 T = 0.003
 Site =    19 /    20 .. Mmpo =     1 DW = 0.00e+00 NNZ =        7 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal =      0.140 Tmvc-total = 0.014 MPO bond dimension =   178 MaxDW = 0.00e+00
NNZ =         3934 SIZE =       200866 SPT = 0.9804

Rank =     0 Ttotal =      0.204 MPO method = FastBipartite bond dimension =     178 NNZ =         3934 SIZE =       200866 SPT = 0.9804

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.935 | E =    -107.6541214061 | DW = 7.63e-08

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.649 | E =    -107.6541223348 | DE = -9.29e-07 | DW = 6.45e-08

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      2.289 | E =    -107.6541224347 | DE = -9.99e-08 | DW = 7.64e-08

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      2.976 | E =    -107.6541224347 | DE = 7.39e-13 | DW = 6.15e-08

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      3.938 | E =    -107.6541224379 | DE = -3.18e-09 | DW = 6.46e-11

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.837 | E =    -107.6541224379 | DE = -3.39e-11 | DW = 6.30e-20

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      5.788 | E =    -107.6541224379 | DE = 0.00e+00 | DW = 6.46e-11

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      7.237 | E =    -107.6541224379 | DE = 6.82e-13 | DW = 7.32e-20

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      8.093 | E =    -107.6541224379 | DE = 2.84e-14 | DW = 6.78e-20

DMRG energy = -107.654122437940956

The SZ Mode

Here we use the get_uhf_integrals function to get the integrals.

[6]:
from pyscf import gto, scf

mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_uhf_integrals(mf,
    ncore=0, ncas=None, g2e_symm=8)

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
    thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error =  1.6638583155958613e-13
integral cutoff error =  0.0
mpo terms =       2778

Build MPO | Nsites =    10 | Nterms =       2778 | Algorithm = FastBIP | Cutoff = 1.00e-20
 Site =     0 /    10 .. Mmpo =    26 DW = 0.00e+00 NNZ =       26 SPT = 0.0000 Tmvc = 0.002 T = 0.009
 Site =     1 /    10 .. Mmpo =    66 DW = 0.00e+00 NNZ =      143 SPT = 0.9167 Tmvc = 0.001 T = 0.015
 Site =     2 /    10 .. Mmpo =   110 DW = 0.00e+00 NNZ =      283 SPT = 0.9610 Tmvc = 0.002 T = 0.014
 Site =     3 /    10 .. Mmpo =   138 DW = 0.00e+00 NNZ =     1023 SPT = 0.9326 Tmvc = 0.001 T = 0.012
 Site =     4 /    10 .. Mmpo =   158 DW = 0.00e+00 NNZ =      535 SPT = 0.9755 Tmvc = 0.001 T = 0.011
 Site =     5 /    10 .. Mmpo =   186 DW = 0.00e+00 NNZ =      463 SPT = 0.9842 Tmvc = 0.001 T = 0.008
 Site =     6 /    10 .. Mmpo =   106 DW = 0.00e+00 NNZ =      415 SPT = 0.9790 Tmvc = 0.001 T = 0.007
 Site =     7 /    10 .. Mmpo =    58 DW = 0.00e+00 NNZ =      163 SPT = 0.9735 Tmvc = 0.000 T = 0.004
 Site =     8 /    10 .. Mmpo =    26 DW = 0.00e+00 NNZ =       87 SPT = 0.9423 Tmvc = 0.000 T = 0.003
 Site =     9 /    10 .. Mmpo =     1 DW = 0.00e+00 NNZ =       26 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal =      0.086 Tmvc-total = 0.009 MPO bond dimension =   186 MaxDW = 0.00e+00
NNZ =         3164 SIZE =       102772 SPT = 0.9692

Rank =     0 Ttotal =      0.126 MPO method = FastBipartite bond dimension =     186 NNZ =         3164 SIZE =       102772 SPT = 0.9692

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.823 | E =    -107.6541224475 | DW = 4.14e-08

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.336 | E =    -107.6541224475 | DE = -1.26e-11 | DW = 5.18e-09

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.794 | E =    -107.6541224475 | DE = -9.66e-13 | DW = 4.14e-08

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      2.282 | E =    -107.6541224475 | DE = 1.22e-12 | DW = 5.37e-09

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      2.835 | E =    -107.6541224475 | DE = -9.95e-13 | DW = 3.60e-11

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      3.405 | E =    -107.6541224475 | DE = 9.09e-13 | DW = 1.69e-19

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      3.992 | E =    -107.6541224475 | DE = 2.84e-14 | DW = 3.60e-11

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      4.568 | E =    -107.6541224475 | DE = -1.22e-12 | DW = 1.63e-19

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      4.956 | E =    -107.6541224475 | DE = -2.84e-14 | DW = 9.42e-20

DMRG energy = -107.654122447524443

The SGF Mode

Here we use the get_ghf_integrals function to get the integrals.

[7]:
from pyscf import gto, scf

mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.GHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_ghf_integrals(mf,
    ncore=0, ncas=None, g2e_symm=8)

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
    thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error =  2.195446592452528e-13
integral cutoff error =  0.0
mpo terms =       5948

Build MPO | Nsites =    20 | Nterms =       5948 | Algorithm = FastBIP | Cutoff = 1.00e-20
 Site =     0 /    20 .. Mmpo =     7 DW = 0.00e+00 NNZ =        7 SPT = 0.0000 Tmvc = 0.001 T = 0.024
 Site =     1 /    20 .. Mmpo =    20 DW = 0.00e+00 NNZ =       19 SPT = 0.8643 Tmvc = 0.003 T = 0.016
 Site =     2 /    20 .. Mmpo =    47 DW = 0.00e+00 NNZ =       49 SPT = 0.9479 Tmvc = 0.003 T = 0.012
 Site =     3 /    20 .. Mmpo =    62 DW = 0.00e+00 NNZ =      251 SPT = 0.9139 Tmvc = 0.003 T = 0.015
 Site =     4 /    20 .. Mmpo =    81 DW = 0.00e+00 NNZ =      269 SPT = 0.9464 Tmvc = 0.003 T = 0.015
 Site =     5 /    20 .. Mmpo =   104 DW = 0.00e+00 NNZ =      355 SPT = 0.9579 Tmvc = 0.003 T = 0.013
 Site =     6 /    20 .. Mmpo =   129 DW = 0.00e+00 NNZ =      563 SPT = 0.9580 Tmvc = 0.004 T = 0.021
 Site =     7 /    20 .. Mmpo =   126 DW = 0.00e+00 NNZ =     2316 SPT = 0.8575 Tmvc = 0.002 T = 0.021
 Site =     8 /    20 .. Mmpo =   155 DW = 0.00e+00 NNZ =      214 SPT = 0.9890 Tmvc = 0.004 T = 0.027
 Site =     9 /    20 .. Mmpo =   148 DW = 0.00e+00 NNZ =      708 SPT = 0.9691 Tmvc = 0.003 T = 0.013
 Site =    10 /    20 .. Mmpo =   181 DW = 0.00e+00 NNZ =      398 SPT = 0.9851 Tmvc = 0.001 T = 0.021
 Site =    11 /    20 .. Mmpo =   178 DW = 0.00e+00 NNZ =      718 SPT = 0.9777 Tmvc = 0.004 T = 0.023
 Site =    12 /    20 .. Mmpo =   147 DW = 0.00e+00 NNZ =      496 SPT = 0.9810 Tmvc = 0.001 T = 0.031
 Site =    13 /    20 .. Mmpo =   100 DW = 0.00e+00 NNZ =      410 SPT = 0.9721 Tmvc = 0.001 T = 0.006
 Site =    14 /    20 .. Mmpo =    77 DW = 0.00e+00 NNZ =      254 SPT = 0.9670 Tmvc = 0.000 T = 0.005
 Site =    15 /    20 .. Mmpo =    54 DW = 0.00e+00 NNZ =      118 SPT = 0.9716 Tmvc = 0.001 T = 0.004
 Site =    16 /    20 .. Mmpo =    39 DW = 0.00e+00 NNZ =      124 SPT = 0.9411 Tmvc = 0.000 T = 0.004
 Site =    17 /    20 .. Mmpo =    20 DW = 0.00e+00 NNZ =       74 SPT = 0.9051 Tmvc = 0.000 T = 0.003
 Site =    18 /    20 .. Mmpo =     7 DW = 0.00e+00 NNZ =       20 SPT = 0.8571 Tmvc = 0.000 T = 0.003
 Site =    19 /    20 .. Mmpo =     1 DW = 0.00e+00 NNZ =        7 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal =      0.282 Tmvc-total = 0.037 MPO bond dimension =   181 MaxDW = 0.00e+00
NNZ =         7370 SIZE =       204350 SPT = 0.9639

Rank =     0 Ttotal =      0.397 MPO method = FastBipartite bond dimension =     181 NNZ =         7370 SIZE =       204350 SPT = 0.9639

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      2.170 | E =    -107.6541220534 | DW = 7.49e-08

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      4.014 | E =    -107.6541223446 | DE = -2.91e-07 | DW = 6.92e-08

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      5.890 | E =    -107.6541224348 | DE = -9.02e-08 | DW = 7.50e-08

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      7.570 | E =    -107.6541224348 | DE = 7.96e-13 | DW = 6.60e-08

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      9.612 | E =    -107.6541224379 | DE = -3.16e-09 | DW = 6.35e-11

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     11.483 | E =    -107.6541224379 | DE = 4.26e-13 | DW = 7.23e-20

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     13.573 | E =    -107.6541224379 | DE = 0.00e+00 | DW = 6.35e-11

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     15.436 | E =    -107.6541224379 | DE = -2.47e-12 | DW = 7.77e-20

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =     16.106 | E =    -107.6541224379 | DE = -2.84e-14 | DW = 6.27e-20

DMRG energy = -107.654122437941297

Relativistic DMRG

For relativistic DMRG, we use the get_dhf_integrals function to get the integrals. We use the SGFCPX Mode in block2 to execute DMRG. Note that the integrals, MPO, and MPS will all contain complex numbers in this mode.

[8]:
from pyscf import gto, scf

mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.DHF(mol).set(with_gaunt=True, with_breit=True).run(conv_tol=1E-12)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_dhf_integrals(mf,
    ncore=0, ncas=None, pg_symm=False)

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGFCPX, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
    thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error =  0.0
integral cutoff error =  2.9170571933900254e-20
mpo terms =      26566

Build MPO | Nsites =    20 | Nterms =      26566 | Algorithm = FastBIP | Cutoff = 1.00e-20
 Site =     0 /    20 .. Mmpo =     9 DW = 0.00e+00 NNZ =        9 SPT = 0.0000 Tmvc = 0.003 T = 0.013
 Site =     1 /    20 .. Mmpo =    28 DW = 0.00e+00 NNZ =       23 SPT = 0.9087 Tmvc = 0.003 T = 0.018
 Site =     2 /    20 .. Mmpo =    63 DW = 0.00e+00 NNZ =      320 SPT = 0.8186 Tmvc = 0.006 T = 0.026
 Site =     3 /    20 .. Mmpo =    78 DW = 0.00e+00 NNZ =      440 SPT = 0.9105 Tmvc = 0.005 T = 0.024
 Site =     4 /    20 .. Mmpo =    97 DW = 0.00e+00 NNZ =      605 SPT = 0.9200 Tmvc = 0.005 T = 0.025
 Site =     5 /    20 .. Mmpo =   120 DW = 0.00e+00 NNZ =      728 SPT = 0.9375 Tmvc = 0.007 T = 0.029
 Site =     6 /    20 .. Mmpo =   147 DW = 0.00e+00 NNZ =      992 SPT = 0.9438 Tmvc = 0.012 T = 0.041
 Site =     7 /    20 .. Mmpo =   178 DW = 0.00e+00 NNZ =     1163 SPT = 0.9556 Tmvc = 0.008 T = 0.033
 Site =     8 /    20 .. Mmpo =   213 DW = 0.00e+00 NNZ =     1971 SPT = 0.9480 Tmvc = 0.007 T = 0.031
 Site =     9 /    20 .. Mmpo =   252 DW = 0.00e+00 NNZ =     2107 SPT = 0.9607 Tmvc = 0.007 T = 0.025
 Site =    10 /    20 .. Mmpo =   213 DW = 0.00e+00 NNZ =    11022 SPT = 0.7947 Tmvc = 0.012 T = 0.055
 Site =    11 /    20 .. Mmpo =   178 DW = 0.00e+00 NNZ =     1717 SPT = 0.9547 Tmvc = 0.003 T = 0.018
 Site =    12 /    20 .. Mmpo =   147 DW = 0.00e+00 NNZ =     1641 SPT = 0.9373 Tmvc = 0.002 T = 0.015
 Site =    13 /    20 .. Mmpo =   120 DW = 0.00e+00 NNZ =     1197 SPT = 0.9321 Tmvc = 0.002 T = 0.014
 Site =    14 /    20 .. Mmpo =    97 DW = 0.00e+00 NNZ =      859 SPT = 0.9262 Tmvc = 0.002 T = 0.011
 Site =    15 /    20 .. Mmpo =    78 DW = 0.00e+00 NNZ =      682 SPT = 0.9099 Tmvc = 0.001 T = 0.008
 Site =    16 /    20 .. Mmpo =    63 DW = 0.00e+00 NNZ =      328 SPT = 0.9333 Tmvc = 0.001 T = 0.006
 Site =    17 /    20 .. Mmpo =    28 DW = 0.00e+00 NNZ =      238 SPT = 0.8651 Tmvc = 0.000 T = 0.004
 Site =    18 /    20 .. Mmpo =     9 DW = 0.00e+00 NNZ =       34 SPT = 0.8651 Tmvc = 0.000 T = 0.003
 Site =    19 /    20 .. Mmpo =     1 DW = 0.00e+00 NNZ =        9 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal =      0.403 Tmvc-total = 0.086 MPO bond dimension =   252 MaxDW = 0.00e+00
NNZ =        26085 SIZE =       323082 SPT = 0.9193

Rank =     0 Ttotal =      0.500 MPO method = FastBipartite bond dimension =     252 NNZ =        26085 SIZE =       323082 SPT = 0.9193

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     27.037 | E =    -107.6929206223 | DW = 6.66e-10

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     39.215 | E =    -107.6929209453 | DE = -3.23e-07 | DW = 7.27e-10

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     49.219 | E =    -107.6929209492 | DE = -3.91e-09 | DW = 2.32e-10

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =     57.774 | E =    -107.6929209492 | DE = -3.98e-13 | DW = 7.32e-10

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     75.576 | E =    -107.6929209492 | DE = -5.12e-13 | DW = 6.34e-20

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =     91.974 | E =    -107.6929209492 | DE = 2.84e-14 | DW = 1.05e-19

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =    112.150 | E =    -107.6929209492 | DE = -2.84e-14 | DW = 7.11e-20

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =    130.356 | E =    -107.6929209492 | DE = -2.84e-14 | DW = 1.12e-19

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =    138.078 | E =    -107.6929209492 | DE = 2.84e-14 | DW = 7.25e-20

DMRG energy = -107.692920949169959

Expectation and N-Particle Density Matrix

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

In this example, we compute the triplet state.

[9]:
from pyscf import gto, scf

mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
    ncore=0, ncas=None, g2e_symm=8)

spin = 2

driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)

mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)

ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
    thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)

impo = driver.get_identity_mpo()

norm = driver.expectation(ket, impo, ket)
ener = driver.expectation(ket, mpo, ket)

print('Norm = %20.15f' % norm)
print('Energy expectation = %20.15f' % (ener / norm))

# <S^2> [ in spin-adapted mode this is always S(S+1) ]
ssq_mpo = driver.get_spin_square_mpo(iprint=0)
ssq = driver.expectation(ket, ssq_mpo, ket)
print('<S^2> expectation = %20.15f' % (ssq / norm))

Sweep =    0 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.632 | E =    -106.9391328597 | DW = 3.38e-10

Sweep =    1 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      0.963 | E =    -106.9391328597 | DE = -4.92e-12 | DW = 8.92e-19

Sweep =    2 | Direction =  forward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.267 | E =    -106.9391328597 | DE = 0.00e+00 | DW = 3.38e-10

Sweep =    3 | Direction = backward | Bond dimension =  250 | Noise =  1.00e-04 | Dav threshold =  1.00e-10
Time elapsed =      1.599 | E =    -106.9391328597 | DE = 1.71e-13 | DW = 1.20e-18

Sweep =    4 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      1.964 | E =    -106.9391328597 | DE = 0.00e+00 | DW = 1.12e-16

Sweep =    5 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      2.288 | E =    -106.9391328597 | DE = 8.53e-14 | DW = 1.23e-19

Sweep =    6 | Direction =  forward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      2.664 | E =    -106.9391328597 | DE = 2.84e-14 | DW = 1.12e-16

Sweep =    7 | Direction = backward | Bond dimension =  500 | Noise =  1.00e-05 | Dav threshold =  1.00e-10
Time elapsed =      2.984 | E =    -106.9391328597 | DE = -2.84e-14 | DW = 1.23e-19

Sweep =    8 | Direction =  forward | Bond dimension =  500 | Noise =  0.00e+00 | Dav threshold =  1.00e-09
Time elapsed =      3.236 | E =    -106.9391328597 | DE = 0.00e+00 | DW = 1.72e-20

DMRG energy = -106.939132859666685
Norm =    0.999999999999999
Energy expectation = -106.939132859666685
<S^2> expectation =    2.000000000000000

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

\[\hat{N}_0 = a^\dagger_{0\alpha} a_{0\alpha} + a^\dagger_{0\beta} a_{0\beta} = \sqrt{2} \big(a_0^\dagger\big)^{[1/2]} \otimes_{[0]} \big(a_0\big)^{[1/2]}\]

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

[10]:
b = driver.expr_builder()
b.add_term("(C+D)0", [0, 0], np.sqrt(2))
n_mpo = driver.get_mpo(b.finalize(), iprint=0)

n_0 = driver.expectation(ket, n_mpo, ket)
print('N0 expectation = %20.15f' % (n_0 / norm))
N0 expectation =    1.999995824360303

We can then verify this number using 1PDM:

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

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

\[\mathrm{DM}_{ijklmn} := \langle E_{ij} E_{kl} E_{mn} \rangle\]

So we have to use the same convention in block2 by setting the npdm_expr parameter in block2 to ((C+D)0+((C+D)0+(C+D)0)0)0.

[12]:
pdm3_b2 = driver.get_3pdm(ket, iprint=0, npdm_expr="((C+D)0+((C+D)0+(C+D)0)0)0")

from pyscf import fci

mx = fci.addons.fix_spin_(fci.FCI(mf), ss=2)
mx.kernel(h1e, g2e, ncas, nelec=n_elec, nroots=3, tol=1E-12)
print(mx.e_tot)
pdm3_fci = fci.rdm.make_dm123('FCI3pdm_kern_sf', mx.ci[0], mx.ci[0], ncas, n_elec)[2]

print('diff = ', np.linalg.norm(pdm3_fci - pdm3_b2))
[-106.93913286 -106.85412245 -106.70055113]
diff =  8.328760079421614e-06