Quantum Chemistry Hamiltonians
[1]:
!pip install block2==0.5.2rc5 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
!pip install pyscf==2.1.1 -qq --progress-bar off
Introduction
In this tutorial we explain how to perform quantum chemistry DMRG using the python interface of block2.
The quantum chemistry Hamiltonian in its second quantized form has to be defined in a set of orbitals, such as the Hartree-Fock (or Density Functional Theory) orbitals. The symmetry that can be used in the DMRG calculation is thus has a dependence on the symmetry of the Hartree-Fock orbitals.
For spin-restricted Hartree-Fock (RHF) orbitals, we can perform spin-adapted DMRG (
SU2mode inblock2) or non-spin-adapted DMRG with any lower symmetries (SZorSGF).For spin-unrestricted Hartree-Fock (UHF) orbitals, we can perform non-spin-adapted DMRG (
SZmode inblock2) or DMRG with lower symmetries (such asSGF).For general Hartree-Fock (GHF) orbitals, we can perform DMRG in spin-orbitals (
SGFmode inblock2) or first translate the Hamiltonian into the qubit Hamiltonian then do DMRG (SGBmode inblock2).For relativistic Dirac Hartree-Fock (DHF) orbitals, we can perform DMRG in complex spin-orbitals (
SGFCPXmode inblock2).
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