Quantum Chemistry Hamiltonians
[1]:
!pip install block2==0.5.3rc5 -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 symmetries that can be used in the DMRG calculation thus have 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 inblock2
) or non-spin-adapted DMRG with any lower symmetries (SZ
orSGF
).For spin-unrestricted Hartree-Fock (UHF) orbitals, we can perform non-spin-adapted DMRG (
SZ
mode inblock2
) or DMRG with lower symmetries (such asSGF
).For general Hartree-Fock (GHF) orbitals, we can perform DMRG in spin-orbitals (
SGF
mode inblock2
) or first translate the Hamiltonian into the qubit Hamiltonian then do DMRG (SGB
mode inblock2
).For relativistic Dirac Hartree-Fock (DHF) orbitals, we can perform DMRG in complex spin-orbitals (
SGFCPX
mode inblock2
).For atom and diatomic molecules, we can perform spin-adapted/non-spin-adapted/spin-orbital DMRG (
SAnySU2LZ/SAnySZLZ/SAnySGFLZ
modes inblock2
) 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
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
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)
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.5986859646912075e-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.008
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.000 T = 0.004
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.001 T = 0.004
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.000 T = 0.004
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.014
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.005
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Ttotal = 0.055 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.125 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.856 | 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.252 | E = -107.6541224475 | DE = -5.37e-12 | DW = 4.75e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.706 | 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.132 | E = -107.6541224475 | DE = -9.09e-13 | DW = 5.74e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.609 | E = -107.6541224475 | DE = 2.84e-14 | DW = 1.65e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.111 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.79e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.623 | E = -107.6541224475 | DE = -5.68e-14 | DW = 2.13e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.076 | E = -107.6541224475 | DE = 2.84e-14 | DW = 6.59e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.452 | E = -107.6541224475 | DE = 2.84e-14 | DW = 2.76e-20
DMRG energy = -107.654122447524415
Energy from pdms = -107.654122447524443
Energy from expectation = -107.654122447524287
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.225136772727681e-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.034
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 143 SPT = 0.9167 Tmvc = 0.001 T = 0.034
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 283 SPT = 0.9610 Tmvc = 0.011 T = 0.042
Site = 3 / 10 .. Mmpo = 138 DW = 0.00e+00 NNZ = 1023 SPT = 0.9326 Tmvc = 0.004 T = 0.018
Site = 4 / 10 .. Mmpo = 158 DW = 0.00e+00 NNZ = 535 SPT = 0.9755 Tmvc = 0.001 T = 0.046
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 463 SPT = 0.9842 Tmvc = 0.001 T = 0.015
Site = 6 / 10 .. Mmpo = 106 DW = 0.00e+00 NNZ = 415 SPT = 0.9790 Tmvc = 0.001 T = 0.021
Site = 7 / 10 .. Mmpo = 58 DW = 0.00e+00 NNZ = 163 SPT = 0.9735 Tmvc = 0.000 T = 0.015
Site = 8 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 87 SPT = 0.9423 Tmvc = 0.000 T = 0.022
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.017
Ttotal = 0.263 Tmvc-total = 0.021 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.504 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.658 | 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.407 | E = -107.6541224475 | DE = -1.32e-11 | DW = 5.09e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.150 | E = -107.6541224475 | DE = -1.31e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.998 | E = -107.6541224475 | DE = 1.42e-12 | DW = 5.21e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.947 | E = -107.6541224475 | DE = -1.19e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.944 | E = -107.6541224475 | DE = 9.09e-13 | DW = 1.83e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.888 | 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 = 7.889 | E = -107.6541224475 | DE = -1.25e-12 | DW = 1.43e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 8.448 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.39e-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 = 7.225136772727683e-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.010
Site = 2 / 20 .. Mmpo = 45 DW = 0.00e+00 NNZ = 45 SPT = 0.9500 Tmvc = 0.001 T = 0.012
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 131 SPT = 0.9530 Tmvc = 0.001 T = 0.020
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 159 SPT = 0.9683 Tmvc = 0.002 T = 0.008
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 203 SPT = 0.9759 Tmvc = 0.001 T = 0.009
Site = 6 / 20 .. Mmpo = 125 DW = 0.00e+00 NNZ = 265 SPT = 0.9796 Tmvc = 0.002 T = 0.019
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 974 SPT = 0.9382 Tmvc = 0.001 T = 0.029
Site = 8 / 20 .. Mmpo = 151 DW = 0.00e+00 NNZ = 188 SPT = 0.9901 Tmvc = 0.001 T = 0.011
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 344 SPT = 0.9846 Tmvc = 0.001 T = 0.015
Site = 10 / 20 .. Mmpo = 177 DW = 0.00e+00 NNZ = 246 SPT = 0.9906 Tmvc = 0.001 T = 0.010
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 402 SPT = 0.9872 Tmvc = 0.001 T = 0.009
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 306 SPT = 0.9883 Tmvc = 0.000 T = 0.012
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 242 SPT = 0.9835 Tmvc = 0.000 T = 0.006
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.010
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.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.004
Ttotal = 0.204 Tmvc-total = 0.014 MPO bond dimension = 178 MaxDW = 0.00e+00
NNZ = 3934 SIZE = 200866 SPT = 0.9804
Rank = 0 Ttotal = 0.341 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 = 1.578 | 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 = 2.716 | 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 = 3.877 | 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 = 4.516 | E = -107.6541224347 | DE = 9.38e-13 | DW = 6.15e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.397 | 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 = 6.224 | E = -107.6541224379 | DE = -3.38e-11 | DW = 7.42e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.123 | 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.919 | E = -107.6541224379 | DE = 6.82e-13 | DW = 7.25e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 8.532 | E = -107.6541224379 | DE = 0.00e+00 | DW = 6.05e-20
DMRG energy = -107.654122437941069
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)
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.4200000000000016e-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.001 T = 0.015
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.002 T = 0.008
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.001 T = 0.008
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.001 T = 0.008
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.004 T = 0.019
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.004
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.007
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.079 Tmvc-total = 0.009 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.140 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.513 | 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.765 | E = -107.6541224475 | DE = -3.13e-12 | DW = 4.68e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.010 | 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.246 | E = -107.6541224475 | DE = -5.97e-13 | DW = 6.02e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.513 | E = -107.6541224475 | DE = -5.68e-14 | DW = 1.68e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.770 | E = -107.6541224475 | DE = 8.53e-14 | DW = 4.69e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.045 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.14e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.294 | E = -107.6541224475 | DE = 0.00e+00 | DW = 4.87e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.472 | E = -107.6541224475 | DE = 2.84e-14 | DW = 4.03e-20
DMRG energy = -107.654122447524557
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.5510746027850855e-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.014
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.013
Site = 3 / 10 .. Mmpo = 138 DW = 0.00e+00 NNZ = 1023 SPT = 0.9326 Tmvc = 0.001 T = 0.011
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.008
Site = 6 / 10 .. Mmpo = 106 DW = 0.00e+00 NNZ = 415 SPT = 0.9790 Tmvc = 0.000 T = 0.008
Site = 7 / 10 .. Mmpo = 58 DW = 0.00e+00 NNZ = 163 SPT = 0.9735 Tmvc = 0.000 T = 0.006
Site = 8 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 87 SPT = 0.9423 Tmvc = 0.000 T = 0.005
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Ttotal = 0.095 Tmvc-total = 0.007 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.170 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.401 | 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.190 | E = -107.6541224475 | DE = -1.80e-11 | DW = 5.07e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.846 | E = -107.6541224475 | DE = -1.11e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.333 | E = -107.6541224475 | DE = 1.11e-12 | DW = 5.10e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.881 | E = -107.6541224475 | DE = -9.38e-13 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.435 | E = -107.6541224475 | DE = 8.81e-13 | DW = 1.28e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.997 | 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 = 5.553 | E = -107.6541224475 | DE = -1.34e-12 | DW = 1.64e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 5.919 | E = -107.6541224475 | DE = -5.68e-14 | DW = 9.83e-20
DMRG energy = -107.654122447524585
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.2056297680002223e-13
integral cutoff error = 0.0
mpo terms = 5992
Build MPO | Nsites = 20 | Nterms = 5992 | 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.017
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.014
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 49 SPT = 0.9479 Tmvc = 0.002 T = 0.015
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 251 SPT = 0.9139 Tmvc = 0.002 T = 0.011
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 273 SPT = 0.9456 Tmvc = 0.003 T = 0.020
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 357 SPT = 0.9576 Tmvc = 0.001 T = 0.013
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 563 SPT = 0.9580 Tmvc = 0.002 T = 0.011
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 2318 SPT = 0.8574 Tmvc = 0.003 T = 0.021
Site = 8 / 20 .. Mmpo = 155 DW = 0.00e+00 NNZ = 212 SPT = 0.9891 Tmvc = 0.001 T = 0.008
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 702 SPT = 0.9694 Tmvc = 0.001 T = 0.010
Site = 10 / 20 .. Mmpo = 181 DW = 0.00e+00 NNZ = 396 SPT = 0.9852 Tmvc = 0.001 T = 0.009
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 728 SPT = 0.9774 Tmvc = 0.001 T = 0.010
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 496 SPT = 0.9810 Tmvc = 0.001 T = 0.008
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 412 SPT = 0.9720 Tmvc = 0.000 T = 0.011
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 = 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.005
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.199 Tmvc-total = 0.021 MPO bond dimension = 181 MaxDW = 0.00e+00
NNZ = 7386 SIZE = 204350 SPT = 0.9639
Rank = 0 Ttotal = 0.296 MPO method = FastBipartite bond dimension = 181 NNZ = 7386 SIZE = 204350 SPT = 0.9639
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.091 | E = -107.6541220194 | DW = 8.01e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.230 | E = -107.6541223360 | DE = -3.17e-07 | DW = 7.44e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.557 | E = -107.6541224309 | DE = -9.49e-08 | DW = 8.01e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.886 | E = -107.6541224309 | DE = 1.79e-12 | DW = 6.88e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.128 | E = -107.6541224313 | DE = -3.64e-10 | DW = 8.01e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.178 | E = -107.6541224313 | DE = 1.31e-12 | DW = 5.93e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.296 | E = -107.6541224313 | DE = -5.68e-14 | DW = 8.01e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 9.355 | E = -107.6541224313 | DE = -2.30e-12 | DW = 8.16e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 10.016 | E = -107.6541224313 | DE = 2.84e-14 | DW = 8.22e-20
DMRG energy = -107.654122431266515
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.03655637126935e-13
integral cutoff error = 0.0
mpo terms = 5906
Build MPO | Nsites = 20 | Nterms = 5906 | 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.009
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.011
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 49 SPT = 0.9479 Tmvc = 0.002 T = 0.012
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 251 SPT = 0.9139 Tmvc = 0.002 T = 0.017
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 273 SPT = 0.9456 Tmvc = 0.003 T = 0.018
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.003 T = 0.017
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 2318 SPT = 0.8574 Tmvc = 0.004 T = 0.019
Site = 8 / 20 .. Mmpo = 155 DW = 0.00e+00 NNZ = 208 SPT = 0.9893 Tmvc = 0.001 T = 0.011
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 684 SPT = 0.9702 Tmvc = 0.001 T = 0.012
Site = 10 / 20 .. Mmpo = 181 DW = 0.00e+00 NNZ = 384 SPT = 0.9857 Tmvc = 0.001 T = 0.011
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 718 SPT = 0.9777 Tmvc = 0.001 T = 0.009
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 496 SPT = 0.9810 Tmvc = 0.001 T = 0.014
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 412 SPT = 0.9720 Tmvc = 0.001 T = 0.010
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.006
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 124 SPT = 0.9411 Tmvc = 0.000 T = 0.005
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 76 SPT = 0.9026 Tmvc = 0.000 T = 0.006
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 22 SPT = 0.8429 Tmvc = 0.000 T = 0.005
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Ttotal = 0.214 Tmvc-total = 0.024 MPO bond dimension = 181 MaxDW = 0.00e+00
NNZ = 7320 SIZE = 204350 SPT = 0.9642
Rank = 0 Ttotal = 0.317 MPO method = FastBipartite bond dimension = 181 NNZ = 7320 SIZE = 204350 SPT = 0.9642
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.172 | E = -107.6541214496 | DW = 5.29e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.994 | E = -107.6541223419 | DE = -8.92e-07 | DW = 6.70e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.754 | E = -107.6541224354 | DE = -9.35e-08 | DW = 5.28e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.536 | E = -107.6541224354 | DE = -1.14e-13 | DW = 6.50e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.606 | E = -107.6541224379 | DE = -2.55e-09 | DW = 6.45e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.653 | E = -107.6541224379 | DE = 2.27e-13 | DW = 5.96e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.786 | E = -107.6541224379 | DE = 5.68e-14 | DW = 6.45e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.843 | E = -107.6541224379 | DE = 1.39e-12 | DW = 6.44e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 8.501 | E = -107.6541224379 | DE = -2.84e-14 | DW = 6.36e-20
DMRG energy = -107.654122437938938
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 = 0.0
mpo terms = 44348
Build MPO | Nsites = 20 | Nterms = 44348 | 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.020
Site = 2 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 404 SPT = 0.7710 Tmvc = 0.006 T = 0.029
Site = 3 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 589 SPT = 0.8801 Tmvc = 0.007 T = 0.029
Site = 4 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 933 SPT = 0.8767 Tmvc = 0.007 T = 0.030
Site = 5 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 1315 SPT = 0.8870 Tmvc = 0.008 T = 0.032
Site = 6 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 1722 SPT = 0.9024 Tmvc = 0.010 T = 0.038
Site = 7 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 2141 SPT = 0.9182 Tmvc = 0.009 T = 0.034
Site = 8 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 2577 SPT = 0.9320 Tmvc = 0.009 T = 0.037
Site = 9 / 20 .. Mmpo = 252 DW = 0.00e+00 NNZ = 2967 SPT = 0.9447 Tmvc = 0.010 T = 0.040
Site = 10 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 17974 SPT = 0.6651 Tmvc = 0.010 T = 0.056
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 2599 SPT = 0.9315 Tmvc = 0.004 T = 0.023
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 2195 SPT = 0.9161 Tmvc = 0.003 T = 0.021
Site = 13 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 1787 SPT = 0.8987 Tmvc = 0.002 T = 0.017
Site = 14 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 1403 SPT = 0.8795 Tmvc = 0.005 T = 0.021
Site = 15 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 1013 SPT = 0.8661 Tmvc = 0.002 T = 0.014
Site = 16 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 652 SPT = 0.8673 Tmvc = 0.001 T = 0.011
Site = 17 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 496 SPT = 0.7188 Tmvc = 0.001 T = 0.012
Site = 18 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 26 SPT = 0.8968 Tmvc = 0.000 T = 0.005
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.488 Tmvc-total = 0.103 MPO bond dimension = 252 MaxDW = 0.00e+00
NNZ = 40834 SIZE = 323082 SPT = 0.8736
Rank = 0 Ttotal = 0.589 MPO method = FastBipartite bond dimension = 252 NNZ = 40834 SIZE = 323082 SPT = 0.8736
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 27.182 | E = -107.6929206708 | DW = 6.28e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 40.137 | E = -107.6929209452 | DE = -2.74e-07 | DW = 7.49e-10
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 51.567 | E = -107.6929209492 | DE = -3.94e-09 | DW = 2.15e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 60.941 | E = -107.6929209492 | DE = -5.40e-13 | DW = 7.52e-10
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 80.487 | E = -107.6929209492 | DE = -4.55e-13 | DW = 7.91e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 98.052 | E = -107.6929209492 | DE = -5.68e-14 | DW = 1.15e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 121.383 | E = -107.6929209492 | DE = 0.00e+00 | DW = 6.27e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 138.658 | E = -107.6929209492 | DE = -2.84e-14 | DW = 1.15e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 147.981 | E = -107.6929209492 | DE = 2.84e-14 | DW = 7.89e-20
DMRG energy = -107.692920949170187
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)
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 = 6.776263578034403e-21
mpo terms = 1882
Build MPO | Nsites = 10 | Nterms = 1882 | 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.036
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 106 SPT = 0.7773 Tmvc = 0.001 T = 0.049
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 189 SPT = 0.9007 Tmvc = 0.001 T = 0.054
Site = 3 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 824 SPT = 0.7771 Tmvc = 0.001 T = 0.058
Site = 4 / 10 .. Mmpo = 88 DW = 0.00e+00 NNZ = 154 SPT = 0.9735 Tmvc = 0.001 T = 0.035
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 318 SPT = 0.9616 Tmvc = 0.000 T = 0.029
Site = 6 / 10 .. Mmpo = 64 DW = 0.00e+00 NNZ = 237 SPT = 0.9606 Tmvc = 0.000 T = 0.024
Site = 7 / 10 .. Mmpo = 40 DW = 0.00e+00 NNZ = 145 SPT = 0.9434 Tmvc = 0.000 T = 0.020
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 49 SPT = 0.9125 Tmvc = 0.000 T = 0.010
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.007
Ttotal = 0.321 Tmvc-total = 0.005 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 2050 SIZE = 29320 SPT = 0.9301
Rank = 0 Ttotal = 0.379 MPO method = FastBipartite bond dimension = 94 NNZ = 2050 SIZE = 29320 SPT = 0.9301
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.704 | 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.071 | E = -107.6541224475 | DE = 1.19e-12 | DW = 3.19e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.438 | 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.808 | E = -107.6541224475 | DE = -8.64e-12 | DW = 3.50e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.291 | E = -107.6541224475 | DE = -2.84e-14 | DW = 5.14e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.713 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.17e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.196 | E = -107.6541224475 | DE = 2.84e-14 | DW = 5.84e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.622 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.86e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.889 | E = -107.6541224475 | DE = 2.84e-14 | DW = 1.94e-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)
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 = 1.249373597200093e-20
mpo terms = 5237
Build MPO | Nsites = 10 | Nterms = 5237 | 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.077
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 278 SPT = 0.8596 Tmvc = 0.001 T = 0.117
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 478 SPT = 0.9342 Tmvc = 0.002 T = 0.126
Site = 3 / 10 .. Mmpo = 123 DW = 0.00e+00 NNZ = 2285 SPT = 0.8311 Tmvc = 0.001 T = 0.131
Site = 4 / 10 .. Mmpo = 170 DW = 0.00e+00 NNZ = 360 SPT = 0.9828 Tmvc = 0.001 T = 0.062
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 814 SPT = 0.9743 Tmvc = 0.001 T = 0.065
Site = 6 / 10 .. Mmpo = 126 DW = 0.00e+00 NNZ = 579 SPT = 0.9753 Tmvc = 0.001 T = 0.046
Site = 7 / 10 .. Mmpo = 82 DW = 0.00e+00 NNZ = 264 SPT = 0.9744 Tmvc = 0.000 T = 0.025
Site = 8 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 187 SPT = 0.9240 Tmvc = 0.000 T = 0.017
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 30 SPT = 0.0000 Tmvc = 0.000 T = 0.011
Ttotal = 0.678 Tmvc-total = 0.009 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 5305 SIZE = 111588 SPT = 0.9525
Rank = 0 Ttotal = 0.740 MPO method = FastBipartite bond dimension = 186 NNZ = 5305 SIZE = 111588 SPT = 0.9525
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.233 | E = -107.6541224475 | DW = 3.06e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.979 | E = -107.6541224475 | DE = -1.18e-11 | DW = 3.32e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.687 | E = -107.6541224455 | DE = 2.05e-09 | DW = 3.06e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.429 | E = -107.6541224455 | DE = -3.41e-13 | DW = 3.27e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.297 | E = -107.6541224455 | DE = -1.14e-13 | DW = 3.09e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.229 | E = -107.6541224455 | DE = 1.14e-13 | DW = 1.40e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.354 | E = -107.6541224455 | DE = -2.84e-14 | DW = 3.09e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.963 | E = -107.6541224455 | DE = 1.71e-13 | DW = 1.58e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 9.102 | E = -107.6541224455 | DE = -2.84e-14 | DW = 1.23e-19
DMRG energy = -107.654122445468445
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)
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 = 1.457638075028426e-14
integral cutoff error = 6.527447649778452e-20
mpo terms = 13645
Build MPO | Nsites = 20 | Nterms = 13645 | 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.100
Site = 1 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 23 SPT = 0.9087 Tmvc = 0.002 T = 0.171
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 307 SPT = 0.7667 Tmvc = 0.003 T = 0.226
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 358 SPT = 0.8771 Tmvc = 0.003 T = 0.276
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 515 SPT = 0.8975 Tmvc = 0.003 T = 0.292
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 667 SPT = 0.9208 Tmvc = 0.003 T = 0.305
Site = 6 / 20 .. Mmpo = 131 DW = 0.00e+00 NNZ = 796 SPT = 0.9416 Tmvc = 0.003 T = 0.309
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 5583 SPT = 0.6618 Tmvc = 0.003 T = 0.305
Site = 8 / 20 .. Mmpo = 161 DW = 0.00e+00 NNZ = 266 SPT = 0.9869 Tmvc = 0.002 T = 0.142
Site = 9 / 20 .. Mmpo = 168 DW = 0.00e+00 NNZ = 732 SPT = 0.9729 Tmvc = 0.002 T = 0.138
Site = 10 / 20 .. Mmpo = 187 DW = 0.00e+00 NNZ = 1400 SPT = 0.9554 Tmvc = 0.002 T = 0.130
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 964 SPT = 0.9710 Tmvc = 0.001 T = 0.110
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 850 SPT = 0.9675 Tmvc = 0.001 T = 0.080
Site = 13 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 667 SPT = 0.9622 Tmvc = 0.001 T = 0.066
Site = 14 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 441 SPT = 0.9621 Tmvc = 0.001 T = 0.041
Site = 15 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 300 SPT = 0.9603 Tmvc = 0.000 T = 0.031
Site = 16 / 20 .. Mmpo = 57 DW = 0.00e+00 NNZ = 417 SPT = 0.9062 Tmvc = 0.000 T = 0.025
Site = 17 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 84 SPT = 0.9474 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.007
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Ttotal = 2.767 Tmvc-total = 0.032 MPO bond dimension = 187 MaxDW = 0.00e+00
NNZ = 14416 SIZE = 229418 SPT = 0.9372
Rank = 0 Ttotal = 2.864 MPO method = FastBipartite bond dimension = 187 NNZ = 14416 SIZE = 229418 SPT = 0.9372
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.314 | E = -107.6541218464 | DW = 1.76e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.490 | E = -107.6541222920 | DE = -4.46e-07 | DW = 5.23e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 5.869 | E = -107.6541224411 | DE = -1.49e-07 | DW = 1.75e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 7.083 | E = -107.6541224411 | DE = -5.34e-12 | DW = 5.24e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.880 | E = -107.6541224448 | DE = -3.71e-09 | DW = 9.90e-12
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 10.506 | E = -107.6541224455 | DE = -6.31e-10 | DW = 5.25e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 12.313 | E = -107.6541224455 | DE = -2.84e-14 | DW = 9.88e-12
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 13.931 | E = -107.6541224455 | DE = 8.53e-14 | DW = 5.12e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 14.954 | E = -107.6541224455 | DE = 0.00e+00 | DW = 8.00e-20
DMRG energy = -107.654122445469611
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.580 | 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.880 | E = -106.9391328597 | DE = -3.04e-12 | DW = 9.74e-19
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.193 | E = -106.9391328597 | DE = -5.68e-14 | DW = 3.38e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.481 | E = -106.9391328597 | DE = 2.56e-13 | DW = 1.53e-18
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.833 | E = -106.9391328597 | DE = 5.68e-14 | DW = 1.12e-16
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.168 | E = -106.9391328597 | DE = -1.42e-13 | DW = 1.35e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.507 | 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.825 | E = -106.9391328597 | DE = 1.14e-13 | DW = 1.38e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.045 | E = -106.9391328597 | DE = 0.00e+00 | DW = 2.47e-20
DMRG energy = -106.939132859666600
Norm = 1.000000000000000
Energy expectation = -106.939132859666586
<S^2> expectation = 2.000000000000001
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.999995824361891
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 = 5.3418848100122e-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.497 | 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.747 | E = -107.6541224475 | DE = -3.33e-12 | DW = 4.05e-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.224 | E = -107.6541224475 | DE = -4.83e-13 | DW = 5.95e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.498 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.70e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.755 | E = -107.6541224475 | DE = 0.00e+00 | DW = 6.84e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.121 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.50e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.617 | E = -107.6541224475 | DE = 2.84e-14 | DW = 6.27e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.911 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.40e-20
DMRG energy = -107.654122447524642
mps center changed (temporarily)
dtrie finished 0.004684378000092693
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368811359562
CSF 0 2222222000 = 0.957506528669410
CSF 1 2222202200 = -0.131287867807398
CSF 2 2222022020 = -0.131287794595548
CSF 3 2222+-2+-0 = 0.118594301088579
CSF 4 2222++2--0 = -0.084968224982026
CSF 5 2220222020 = -0.054710635444609
CSF 6 2220222200 = -0.054710465272209
CSF 7 2222+2+0-- = 0.053881241407867
CSF 8 22222++-0- = 0.053881215166800
Extracting determinant coefficients from non-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 = 0.826 | 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.306 | E = -107.6541224475 | DE = -1.90e-12 | DW = 5.67e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.783 | 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 = 2.259 | E = -107.6541224475 | DE = 1.42e-12 | DW = 5.13e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.822 | E = -107.6541224475 | DE = 2.07e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.380 | E = -107.6541224475 | DE = -3.24e-12 | DW = 1.71e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.960 | 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.519 | E = -107.6541224475 | DE = 8.81e-13 | DW = 1.76e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.910 | E = -107.6541224475 | DE = -2.84e-14 | DW = 7.05e-20
DMRG energy = -107.654122447523420
mps center changed (temporarily)
dtrie finished 0.0024573000000600587
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331592192659
DET 0 2222222000 = -0.957506477722952
DET 1 2222202200 = 0.131287949363695
DET 2 2222022020 = 0.131287934434239
DET 3 2222ab2ab0 = -0.083825435733035
DET 4 2222ba2ba0 = -0.083825435117335
DET 5 2220222020 = 0.054710529942414
DET 6 2220222200 = 0.054710516357351
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 common 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.517 | 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.760 | E = -107.6541224475 | DE = -1.25e-11 | DW = 5.82e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.007 | 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.256 | E = -107.6541224475 | DE = -1.39e-12 | DW = 4.80e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.533 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.06e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.779 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.63e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.052 | E = -107.6541224475 | DE = 5.68e-14 | DW = 4.42e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.328 | E = -107.6541224475 | DE = -5.68e-14 | DW = 5.11e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.508 | E = -107.6541224475 | DE = 2.84e-14 | DW = 5.29e-20
DMRG energy = -107.654122447524387
mps center changed (temporarily)
dtrie finished 0.0019679109999515276
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368809657108
CSF 0 2222222000 = 0.957506521875156
CSF 1 2222202200 = -0.131287813613838
CSF 2 2222022020 = -0.131287676854108
CSF 3 2222+-2+-0 = 0.118594392287611
CSF 4 2222++2--0 = -0.084968279487238
CSF 5 2220222020 = -0.054710673444325
CSF 6 2220222200 = -0.054710596031564
CSF 7 2222+2+0-- = 0.053881289057961
CSF 8 22222++-0- = 0.053881233355426
0.9843688096571082
-107.61553663531718
Sweep = 0 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.052 | E = -107.6285064669 | DW = 6.01e-21
Sweep = 1 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.102 | E = -107.6288048128 | DE = -2.98e-04 | DW = 1.43e-21
Sweep = 2 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.155 | E = -107.6289019710 | DE = -9.72e-05 | DW = 1.63e-21
Sweep = 3 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.206 | E = -107.6289019713 | DE = -3.35e-10 | DW = 1.14e-20
Ground state energy = -107.628901971339630
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.819 | 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.305 | E = -107.6541224475 | DE = -5.66e-12 | DW = 5.12e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.770 | E = -107.6541224475 | DE = -6.25e-13 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.241 | E = -107.6541224475 | DE = 1.14e-12 | DW = 4.98e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.812 | 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.601 | E = -107.6541224475 | DE = 8.53e-13 | DW = 1.77e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.558 | 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 = 5.488 | E = -107.6541224475 | DE = -1.14e-12 | DW = 1.67e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.143 | E = -107.6541224475 | DE = -2.84e-14 | DW = 4.96e-20
DMRG energy = -107.654122447524500
mps center changed (temporarily)
dtrie finished 0.004142805999890697
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331639061928
DET 0 2222222000 = -0.957506569121758
DET 1 2222022020 = 0.131287834673405
DET 2 2222202200 = 0.131287815095743
DET 3 2222ba2ba0 = -0.083825277849152
DET 4 2222ab2ab0 = -0.083825277218030
DET 5 2220222200 = 0.054710476468495
DET 6 2220222020 = 0.054710443995518
0.9713316390619284
-107.5944732555996
Sweep = 0 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.090 | E = -107.6101011814 | DW = 8.24e-21
Sweep = 1 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.185 | E = -107.6102811404 | DE = -1.80e-04 | DW = 1.03e-20
Sweep = 2 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.282 | E = -107.6104759466 | DE = -1.95e-04 | DW = 1.28e-20
Sweep = 3 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.375 | E = -107.6104759470 | DE = -4.56e-10 | DW = 1.14e-20
Ground state energy = -107.610475947023915
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.053 | E = -107.5033999317 | DW = 2.74e-21
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.107 | E = -107.5830578672 | DE = -7.97e-02 | DW = 5.53e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.175 | E = -107.5836402840 | DE = -5.82e-04 | DW = 1.44e-20
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.238 | E = -107.5867932582 | DE = -3.15e-03 | DW = 3.21e-21
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.292 | 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.355 | E = -107.5868152464 | DE = -9.07e-06 | DW = 9.11e-21
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.412 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.51e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.466 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.40e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.521 | E = -107.5868152464 | DE = 0.00e+00 | DW = 2.11e-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.060 | E = -107.5033999316 | DW = 1.96e-21
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.131 | E = -107.5830508082 | DE = -7.97e-02 | DW = 5.41e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.204 | E = -107.5836402367 | DE = -5.89e-04 | DW = 3.59e-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.80e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.382 | E = -107.5868061676 | DE = -1.29e-05 | DW = 1.53e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.469 | E = -107.5868152464 | DE = -9.08e-06 | DW = 1.02e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.557 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.19e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.641 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.39e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.717 | E = -107.5868152464 | DE = 2.84e-14 | DW = 1.67e-20
DMRG energy = -107.586815246376702
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 = 1.097 | E = -107.6538892686 | DW = 2.59e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.444 | E = -107.6539226748 | DE = -3.34e-05 | DW = 9.28e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.623 | E = -107.6539461908 | DE = -2.35e-05 | DW = 5.62e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.728 | E = -107.6539506829 | DE = -4.49e-06 | DW = 1.96e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.627 | E = -107.6539506829 | DE = 3.07e-12 | DW = 2.73e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.502 | E = -107.6539506829 | DE = 7.39e-13 | DW = 4.70e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.419 | E = -107.6539506829 | DE = 0.00e+00 | DW = 2.73e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.281 | E = -107.6539506829 | DE = -3.84e-12 | DW = 4.67e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 8.825 | E = -107.6539506829 | DE = -2.84e-14 | DW = 5.57e-20
DMRG energy = -107.653950682886418
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.499 | 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.747 | E = -107.6541224475 | DE = -4.09e-12 | DW = 5.36e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.058 | 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.439 | E = -107.6541224475 | DE = -5.97e-13 | DW = 6.52e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.891 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.75e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.294 | E = -107.6541224475 | DE = -2.84e-14 | DW = 4.50e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.744 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.05e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.162 | E = -107.6541224475 | DE = 2.84e-14 | DW = 4.93e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.455 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.44e-20
DMRG energy = -107.654122447524585
mps center changed (temporarily)
dtrie finished 0.002923151000004509
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368801370672
CSF 0 2222222000 = -0.957506527302991
CSF 1 2222202200 = 0.131287890948123
CSF 2 2222022020 = 0.131287890439278
CSF 3 2222+-2+-0 = -0.118594401232681
CSF 4 2222++2--0 = 0.084968039984560
CSF 5 2220222020 = 0.054710469983632
CSF 6 2220222200 = 0.054710326229694
CSF 7 2222+2+0-- = -0.053881243065345
CSF 8 22222++-0- = -0.053881235680181
-107.65412244752449
mps center changed (temporarily)
dtrie finished 0.004040642000063599
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331607927357
DET 0 2222222000 = -0.957506527302991
DET 1 2222202200 = 0.131287890948123
DET 2 2222022020 = 0.131287890439278
DET 3 2222ba2ba0 = -0.083825360995141
DET 4 2222ab2ab0 = -0.083825360995141
DET 5 2220222020 = 0.054710469983632
DET 6 2220222200 = 0.054710326229694
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.521 | 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.758 | E = -107.6541224475 | DE = -1.76e-12 | DW = 6.23e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.002 | 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.260 | E = -107.6541224475 | DE = -2.84e-13 | DW = 5.93e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.552 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.65e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.792 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.19e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.053 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.83e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.301 | E = -107.6541224475 | DE = 5.68e-14 | DW = 6.71e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.491 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.80e-20
DMRG energy = -107.654122447524671
1.0
-107.65412244752463
integral symmetrize error = 1.758383454541612e-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.006
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.001 T = 0.006
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.001 T = 0.005
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.001 T = 0.006
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.005
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.004
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.003
Ttotal = 0.044 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.081 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
(0.9999999999999999+0j)
(-107.65412244752464+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 = 0.593 | 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.095 | E = -107.6541224475 | DE = -7.25e-12 | DW = 1.53e-19
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.625 | 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.922 | E = -107.6541224475 | DE = -1.22e-12 | DW = 1.47e-19
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.467 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.49e-19
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.856 | E = -107.6541224475 | DE = 5.68e-14 | DW = 5.30e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.193 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.30e-19
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.511 | E = -107.6541224475 | DE = 0.00e+00 | DW = 6.00e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.733 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.78e-20
DMRG energy = -107.654122447524557
(1+3.326128828347501e-17j)
(-107.65412244752451+5.487212746008987e-15j)
integral symmetrize error = 1.9716021141012156e-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.006
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.000 T = 0.011
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.001 T = 0.005
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.005
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.004
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.003
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.047 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.079 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
0.8950611475659227
-107.6541224475247
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
/SGF
mode, where \(\Lambda_k\) are singular values in the bond at site \(k\).
In the SZ
mode:
[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_sz = driver.get_bipartite_entanglement()
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.789 | E = -107.6541224475 | DW = 4.15e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.279 | E = -107.6541224475 | DE = -4.43e-12 | DW = 5.20e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.734 | E = -107.6541224475 | DE = 1.51e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.209 | E = -107.6541224475 | DE = -1.62e-12 | DW = 5.14e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.793 | E = -107.6541224475 | DE = -1.25e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.342 | E = -107.6541224475 | DE = 9.95e-13 | DW = 2.04e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.991 | 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.931 | E = -107.6541224475 | DE = -1.14e-12 | DW = 1.69e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 5.595 | E = -107.6541224475 | DE = 5.68e-14 | DW = 7.96e-20
DMRG energy = -107.654122447524585
In the SGF
mode:
[29]:
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=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_sgf = driver.get_bipartite_entanglement()
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.182 | E = -107.6541211614 | DW = 8.43e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.011 | E = -107.6541223352 | DE = -1.17e-06 | DW = 7.39e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.806 | E = -107.6541224309 | DE = -9.56e-08 | DW = 8.43e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.602 | E = -107.6541224309 | DE = 3.58e-12 | DW = 6.84e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.655 | E = -107.6541224313 | DE = -4.10e-10 | DW = 8.63e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.697 | E = -107.6541224313 | DE = 1.90e-12 | DW = 8.46e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.823 | E = -107.6541224313 | DE = 8.53e-14 | DW = 8.62e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.152 | E = -107.6541224313 | DE = -2.73e-12 | DW = 7.37e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 9.272 | E = -107.6541224313 | DE = 0.00e+00 | DW = 7.39e-20
DMRG energy = -107.654122431266543
We can then plot the bipartite entanglement in MPS with spatial orbital (SZ
) and spin orbital (SGF
) sites:
[36]:
import matplotlib.pyplot as plt
plt.plot(np.arange(len(bip_ent_sz)) * 2, bip_ent_sz, linestyle='-', marker='o',
mfc='white', mec="#7FB685", color="#7FB685", label='SZ')
plt.plot(np.arange(len(bip_ent_sgf)), bip_ent_sgf, linestyle='-', marker='s',
mfc='white', mec="#3C5B39", color="#3C5B39", label='SGF')
plt.xlim((-1, ncas - 1))
plt.legend(loc='upper left')
plt.xlabel("site index $k$")
plt.ylabel("bipartite entanglement $S_k$")
plt.gcf().set_dpi(150)
plt.show()
Orbital Entropy and Mutual Information
For the optimized MPS in the SZ
/SGF
mode, we can compute the 1- and 2- orbital density matrices and mutual information for pairs of orbitals.
In the SZ
mode:
[37]:
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)
ordm1 = driver.get_orbital_entropies(ket, orb_type=1)
ordm2 = driver.get_orbital_entropies(ket, orb_type=2)
minfo = 0.5 * (ordm1[:, None] + ordm1[None, :] - ordm2) * (1 - np.identity(len(ordm1)))
import matplotlib.pyplot as plt
plt.matshow(minfo, cmap='ocean_r')
plt.gcf().set_dpi(150)
plt.show()
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.404 | 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.188 | E = -107.6541224475 | DE = -1.58e-11 | DW = 5.10e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.906 | 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 = 3.666 | E = -107.6541224475 | DE = 1.48e-12 | DW = 5.15e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.563 | E = -107.6541224475 | DE = -7.11e-13 | DW = 3.61e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.463 | E = -107.6541224475 | DE = 9.09e-13 | DW = 1.59e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.332 | 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 = 7.195 | E = -107.6541224475 | DE = -1.68e-12 | DW = 1.90e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 7.761 | E = -107.6541224475 | DE = 0.00e+00 | DW = 6.47e-20
DMRG energy = -107.654122447524557
In the SGF
mode:
[38]:
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=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)
ordm1 = driver.get_orbital_entropies(ket, orb_type=1)
ordm2 = driver.get_orbital_entropies(ket, orb_type=2)
minfo = 0.5 * (ordm1[:, None] + ordm1[None, :] - ordm2) * (1 - np.identity(len(ordm1)))
import matplotlib.pyplot as plt
plt.matshow(minfo, cmap='ocean_r')
plt.gcf().set_dpi(150)
plt.show()
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.852 | E = -107.6541217988 | DW = 7.63e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.490 | E = -107.6541223313 | DE = -5.33e-07 | DW = 7.48e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 5.184 | E = -107.6541224348 | DE = -1.04e-07 | DW = 7.65e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 6.857 | E = -107.6541224348 | DE = 6.82e-12 | DW = 7.21e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.588 | E = -107.6541224379 | DE = -3.11e-09 | DW = 8.57e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 10.540 | E = -107.6541224379 | DE = -9.01e-12 | DW = 6.84e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 12.280 | E = -107.6541224379 | DE = -2.84e-14 | DW = 8.57e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 13.677 | E = -107.6541224379 | DE = 2.84e-12 | DW = 9.28e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 14.361 | E = -107.6541224379 | DE = -2.84e-14 | DW = 7.96e-20
DMRG energy = -107.654122437939137
Orbital Reordering
We support several algorithms (such as “fiedler” and “gaopt”) to reorder the orbitals used in DMRG, so that the entanglement in the MPS can be decreased. With an optimal orbital reordering, one can expect a better energy expectation with a low bond dimension MPS.
The orbital reordering can be applied directly on the integrals h1e
, g2e
(unpacked), and orb_sym
(if there are point group symmetries), before they are used by pyblock2
. This is the recommended way, and the user is responsible for tracking the orbital ordering.
In the SZ
mode (using the exchange integral as the cost function):
[39]:
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_orig, g2e_orig, orb_sym_orig = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
g2e_orig = driver.unpack_g2e(g2e_orig, n_sites=ncas)
minfos = {}
for method in ["fiedler", "gaopt"]:
idx = driver.orbital_reordering(h1e_orig, g2e_orig, method=method)
h1e = h1e_orig[idx][:, idx]
g2e = g2e_orig[idx][:, idx][:, :, idx][:, :, :, idx]
orb_sym = np.array(orb_sym_orig)[idx]
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('Method = %s DMRG energy = %20.15f' % (method, energy))
minfos[method] = driver.get_orbital_interaction_matrix(ket)
import matplotlib.pyplot as plt
f, axs = plt.subplots(1, 2, sharey=True)
f.set_size_inches(w=14 / 2, h=7 / 2)
for ax, (method, minfo) in zip(axs, minfos.items()):
ax.matshow(minfo, cmap='ocean_r')
ax.set_title(method)
plt.gcf().set_dpi(150)
plt.show()
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.154 | E = -107.6541224475 | DW = 3.94e-14
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.779 | E = -107.6541224475 | DE = -1.71e-13 | DW = 5.09e-14
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.519 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.53e-14
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.267 | E = -107.6541224475 | DE = -1.71e-13 | DW = 5.99e-14
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.325 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.60e-19
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.298 | E = -107.6541224475 | DE = 0.00e+00 | DW = 9.10e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.249 | E = -107.6541224475 | DE = 1.42e-13 | DW = 1.68e-19
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.029 | E = -107.6541224475 | DE = 5.68e-14 | DW = 8.36e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 7.552 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.65e-20
Method = fiedler DMRG energy = -107.654122447502473
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.165 | E = -107.6541224475 | DW = 4.14e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.551 | E = -107.6541224475 | DE = -1.21e-11 | DW = 4.20e-11
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.929 | E = -107.6541224475 | DE = -2.30e-12 | DW = 4.13e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.314 | E = -107.6541224475 | DE = -2.27e-13 | DW = 4.41e-11
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.783 | E = -107.6541224475 | DE = 5.68e-14 | DW = 8.96e-21
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.298 | E = -107.6541224475 | DE = -1.14e-13 | DW = 1.20e-14
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.812 | E = -107.6541224475 | DE = 5.68e-14 | DW = 1.14e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.337 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.21e-14
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.638 | E = -107.6541224475 | DE = -1.42e-13 | DW = 4.79e-20
Method = gaopt DMRG energy = -107.654122447524131
In the SGF
mode (using the exchange integral as the cost function):
[40]:
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_orig, g2e_orig, orb_sym_orig = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
g2e_orig = driver.unpack_g2e(g2e_orig, n_sites=ncas)
minfos = {}
for method in ["fiedler", "gaopt"]:
idx = driver.orbital_reordering(h1e_orig, g2e_orig, method=method)
h1e = h1e_orig[idx][:, idx]
g2e = g2e_orig[idx][:, idx][:, :, idx][:, :, :, idx]
orb_sym = np.array(orb_sym_orig)[idx]
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('Method = %s DMRG energy = %20.15f' % (method, energy))
minfos[method] = driver.get_orbital_interaction_matrix(ket)
import matplotlib.pyplot as plt
f, axs = plt.subplots(1, 2, sharey=True)
f.set_size_inches(w=14 / 2, h=7 / 2)
for ax, (method, minfo) in zip(axs, minfos.items()):
ax.matshow(minfo, cmap='ocean_r')
ax.set_title(method)
plt.gcf().set_dpi(150)
plt.show()
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.089 | E = -107.6421167128 | DW = 1.94e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 6.269 | E = -107.6541224404 | DE = -1.20e-02 | DW = 2.02e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 7.843 | E = -107.6541224449 | DE = -4.57e-09 | DW = 2.63e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 9.347 | E = -107.6541224419 | DE = 3.04e-09 | DW = 2.48e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 11.393 | E = -107.6541224475 | DE = -5.62e-09 | DW = 2.71e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 13.951 | E = -107.6541224475 | DE = -9.66e-13 | DW = 2.95e-11
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 16.841 | E = -107.6541224475 | DE = -3.13e-13 | DW = 2.74e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 18.893 | E = -107.6541224475 | DE = 1.99e-13 | DW = 3.67e-11
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 19.757 | E = -107.6541224475 | DE = 7.11e-12 | DW = 3.64e-20
Method = fiedler DMRG energy = -107.654122447517537
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.646 | E = -107.6540666856 | DW = 5.18e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.397 | E = -107.6541222892 | DE = -5.56e-05 | DW = 4.78e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 5.321 | E = -107.6541224442 | DE = -1.55e-07 | DW = 5.35e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 6.216 | E = -107.6541224442 | DE = -8.87e-12 | DW = 4.55e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.531 | E = -107.6541224475 | DE = -3.35e-09 | DW = 1.02e-10
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.955 | E = -107.6541224475 | DE = -6.82e-13 | DW = 2.97e-11
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 10.368 | E = -107.6541224475 | DE = 4.55e-13 | DW = 1.01e-10
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 11.805 | E = -107.6541224475 | DE = -7.67e-13 | DW = 2.97e-11
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 12.680 | E = -107.6541224475 | DE = 3.69e-13 | DW = 2.59e-20
Method = gaopt DMRG energy = -107.654122447524415
The orbital reordering can also be done implicitly using the reorder
argument in DMRGDriver.get_qc_mpo
. pyblock2
will automatically recover the original ordering for observables, whenever possible, such as \(N\)-PDMs. But ordering for quantities like bipartite/orbital entropy are not recovered (kept as the computational ordering since they are indicators of the orbital topology). In order to avoid any confusion, this implicit reordering approach is not recommended unless only the
ground state energy and \(N\)-PDMs are needed.
Implicit orbital reordering using gaopt
in the SZ
mode:
[41]:
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, reorder="gaopt", 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)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.869 | E = -107.6541224475 | DW = 4.11e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.480 | E = -107.6541224475 | DE = -1.18e-11 | DW = 4.28e-11
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.116 | E = -107.6541224475 | DE = -2.10e-12 | DW = 4.13e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.707 | E = -107.6541224475 | DE = -1.42e-13 | DW = 4.76e-11
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.456 | E = -107.6541224475 | DE = -8.53e-14 | DW = 1.46e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.258 | E = -107.6541224475 | DE = 5.68e-14 | DW = 1.22e-14
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.046 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.52e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.088 | E = -107.6541224475 | DE = -5.68e-14 | DW = 1.21e-14
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.794 | E = -107.6541224475 | DE = 0.00e+00 | DW = 4.16e-20
DMRG energy = -107.654122447523960
One can also run a small bond dimension DMRG with the default orbital ordering to get the mutual information for orbital correlation, then use this matrix as the cost function to get the orbital ordering for doing more efficient DMRG (this is supported in the SZ
and SGF
modes). See the following example.
In the SGF
mode (using the mutual information as the cost function):
[43]:
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_orig, g2e_orig, orb_sym_orig = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
g2e_orig = driver.unpack_g2e(g2e_orig, n_sites=ncas)
# approx DMRG to get orbital_interaction_matrix
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym_orig)
mpo = driver.get_qc_mpo(h1e=h1e_orig, g2e=g2e_orig, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=50, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=10, bond_dims=[50] * 8, noises=noises,
thrds=thrds, iprint=1)
print('Approx DMRG energy = %20.15f' % energy)
minfo_orig = driver.get_orbital_interaction_matrix(ket)
minfos = {}
for method in ["fiedler", "gaopt"]:
idx = driver.orbital_reordering_interaction_matrix(minfo_orig, method=method)
h1e = h1e_orig[idx][:, idx]
g2e = g2e_orig[idx][:, idx][:, :, idx][:, :, :, idx]
orb_sym = np.array(orb_sym_orig)[idx]
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('Method = %s DMRG energy = %20.15f' % (method, energy))
minfos[method] = driver.get_orbital_interaction_matrix(ket)
import matplotlib.pyplot as plt
f, axs = plt.subplots(1, 2, sharey=True)
f.set_size_inches(w=14 / 2, h=7 / 2)
for ax, (method, minfo) in zip(axs, minfos.items()):
ax.matshow(minfo, cmap='ocean_r')
ax.set_title(method)
plt.gcf().set_dpi(150)
plt.show()
Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.790 | E = -107.6461019596 | DW = 1.12e-04
Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.199 | E = -107.6518646319 | DE = -5.76e-03 | DW = 8.45e-05
Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.658 | E = -107.6524131412 | DE = -5.49e-04 | DW = 4.76e-05
Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.123 | E = -107.6524673390 | DE = -5.42e-05 | DW = 4.58e-05
Sweep = 4 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.589 | E = -107.6525301747 | DE = -6.28e-05 | DW = 3.18e-05
Sweep = 5 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.076 | E = -107.6525427750 | DE = -1.26e-05 | DW = 3.17e-05
Sweep = 6 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.514 | E = -107.6525457470 | DE = -2.97e-06 | DW = 3.25e-05
Sweep = 7 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.945 | E = -107.6525492373 | DE = -3.49e-06 | DW = 3.25e-05
Sweep = 8 | Direction = forward | Bond dimension = 50 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.315 | E = -107.6525424769 | DE = 6.76e-06 | DW = 2.77e-05
Sweep = 9 | Direction = backward | Bond dimension = 50 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.642 | E = -107.6525448992 | DE = -2.42e-06 | DW = 2.79e-05
ATTENTION: DMRG is not converged to desired tolerance of 1.00e-08
Approx DMRG energy = -107.652544899236361
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.878 | E = -107.6541143612 | DW = 4.90e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.658 | E = -107.6541224474 | DE = -8.09e-06 | DW = 2.78e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.499 | E = -107.6541224475 | DE = -1.10e-10 | DW = 4.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.253 | E = -107.6541224475 | DE = 1.76e-12 | DW = 2.78e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.260 | E = -107.6541224475 | DE = -3.67e-12 | DW = 7.63e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.300 | E = -107.6541224475 | DE = -5.68e-14 | DW = 9.52e-12
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.397 | E = -107.6541224475 | DE = 1.42e-12 | DW = 8.05e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.455 | E = -107.6541224475 | DE = 1.14e-13 | DW = 9.44e-12
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 9.077 | E = -107.6541224475 | DE = 3.27e-12 | DW = 2.36e-20
Method = fiedler DMRG energy = -107.654122447518702
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.205 | E = -107.6541223034 | DW = 1.54e-09
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.032 | E = -107.6541224434 | DE = -1.40e-07 | DW = 2.01e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.921 | E = -107.6541224475 | DE = -4.10e-09 | DW = 1.56e-09
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.728 | E = -107.6541224475 | DE = 2.64e-12 | DW = 2.01e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.899 | E = -107.6541224475 | DE = -5.54e-12 | DW = 1.05e-13
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.131 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.63e-14
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.302 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.05e-13
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 11.248 | E = -107.6541224475 | DE = -5.68e-14 | DW = 3.58e-14
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 12.502 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.66e-20
Method = gaopt DMRG energy = -107.654122447523818
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.
[44]:
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 = 1.873 | E[ 3] = -107.6541184424 -107.0314355209 -106.9595839138 | DW = 8.49e-05
Sweep = 1 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.807 | E[ 3] = -107.6541184424 -107.0314355209 -106.9595839138 | DE = -1.51e-11 | DW = 1.61e-04
Sweep = 2 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.640 | E[ 3] = -107.6541138875 -107.0314486179 -106.9594564029 | DE = 1.28e-04 | DW = 7.20e-05
Sweep = 3 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.520 | E[ 3] = -107.6541138875 -107.0314486178 -106.9594564029 | DE = -9.18e-12 | DW = 1.47e-04
Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10
Time elapsed = 5.353 | E[ 3] = -107.6540972369 -107.0314487272 -106.9594391430 | 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 = 6.315 | E[ 3] = -107.6540972369 -107.0314487272 -106.9594391430 | DE = -1.78e-11 | DW = 1.43e-04
State-averaged MPS energies = [-107.654097236858902 -107.031448727151187 -106.959439143016382]
State-specific MPS E[0] = -107.654122447521317
State-specific MPS E[1] = -107.031449471624697
State-specific MPS E[2] = -106.959626154671795