Quantum Chemistry Hamiltonians
[1]:
!pip install block2==0.5.2rc13 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
!pip install pyscf==2.3.0 -qq --progress-bar off
Introduction
In this tutorial we explain how to perform quantum chemistry DMRG using the python interface of block2.
The quantum chemistry Hamiltonian in its second quantized form has to be defined in a set of orbitals, such as the Hartree-Fock (or Density Functional Theory) orbitals. The symmetry that can be used in the DMRG calculation is thus has a dependence on the symmetry of the Hartree-Fock orbitals.
For spin-restricted Hartree-Fock (RHF) orbitals, we can perform spin-adapted DMRG (
SU2mode inblock2) or non-spin-adapted DMRG with any lower symmetries (SZorSGF).For spin-unrestricted Hartree-Fock (UHF) orbitals, we can perform non-spin-adapted DMRG (
SZmode inblock2) or DMRG with lower symmetries (such asSGF).For general Hartree-Fock (GHF) orbitals, we can perform DMRG in spin-orbitals (
SGFmode inblock2) or first translate the Hamiltonian into the qubit Hamiltonian then do DMRG (SGBmode inblock2).For relativistic Dirac Hartree-Fock (DHF) orbitals, we can perform DMRG in complex spin-orbitals (
SGFCPXmode inblock2).For atom and diatomic molecules, we can perform spin-adapted/non-spin-adapted/spin-orbital DMRG (
SAnySU2LZ/SAnySZLZ/SAnySGFLZmodes 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
Spin-Restricted Integrals
Here we use get_rhf_integrals function to get the integrals. Note that in order to do DMRG in a CASCI space, one can set the ncore (number of core orbitals) and ncas (number of active orbitals) parameters in get_*_integrals. ncas=None will include all orbitals in DMRG.
For medium to large scale DMRG calculations, it is highly recommended to use a scratch space with high IO speed rather than the ./tmp used in the following example. One also needs to set a suitable stack_mem in the DMRGDriver constructor to set the memory used for storing renormalized operators (in bytes). The default is stack_mem=int(1024**3) (1 GB). For medium scale calculations 10 to 30 GB might be required.
For the meaning of DMRG parameters, please have a look at the Hubbard - Run DMRG page.
[3]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
pdm1 = driver.get_1pdm(ket)
pdm2 = driver.get_2pdm(ket).transpose(0, 3, 1, 2)
print('Energy from pdms = %20.15f' % (np.einsum('ij,ij->', pdm1, h1e)
+ 0.5 * np.einsum('ijkl,ijkl->', pdm2, driver.unpack_g2e(g2e)) + ecore))
impo = driver.get_identity_mpo()
expt = driver.expectation(ket, mpo, ket) / driver.expectation(ket, impo, ket)
print('Energy from expectation = %20.15f' % expt)
integral symmetrize error = 5.977257773040786e-14
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 13 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.000 T = 0.010
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.000 T = 0.009
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.009
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.000 T = 0.010
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.007
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.009
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.004
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.007
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.006
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.010
Ttotal = 0.082 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.197 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.239 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.842 | E = -107.6541224475 | DE = -5.40e-12 | DW = 5.66e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.663 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.203 | E = -107.6541224475 | DE = -7.96e-13 | DW = 5.72e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.918 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.92e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.457 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.42e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.188 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.01e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.683 | E = -107.6541224475 | DE = -2.84e-14 | DW = 5.57e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.061 | E = -107.6541224475 | DE = 2.84e-14 | DW = 2.64e-20
DMRG energy = -107.654122447524443
Energy from pdms = -107.654122447524387
Energy from expectation = -107.654122447524344
We can also run non-spin-adapted DMRG (SZ mode) using the restricted integrals.
[4]:
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error = 7.512924107430347e-14
integral cutoff error = 0.0
mpo terms = 2778
Build MPO | Nsites = 10 | Nterms = 2778 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.001 T = 0.005
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 143 SPT = 0.9167 Tmvc = 0.001 T = 0.006
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 283 SPT = 0.9610 Tmvc = 0.001 T = 0.008
Site = 3 / 10 .. Mmpo = 138 DW = 0.00e+00 NNZ = 1023 SPT = 0.9326 Tmvc = 0.004 T = 0.013
Site = 4 / 10 .. Mmpo = 158 DW = 0.00e+00 NNZ = 535 SPT = 0.9755 Tmvc = 0.001 T = 0.010
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 463 SPT = 0.9842 Tmvc = 0.001 T = 0.010
Site = 6 / 10 .. Mmpo = 106 DW = 0.00e+00 NNZ = 415 SPT = 0.9790 Tmvc = 0.000 T = 0.010
Site = 7 / 10 .. Mmpo = 58 DW = 0.00e+00 NNZ = 163 SPT = 0.9735 Tmvc = 0.000 T = 0.025
Site = 8 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 87 SPT = 0.9423 Tmvc = 0.000 T = 0.010
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Ttotal = 0.103 Tmvc-total = 0.009 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.178 MPO method = FastBipartite bond dimension = 186 NNZ = 3164 SIZE = 102772 SPT = 0.9692
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.281 | E = -107.6541224475 | DW = 4.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.089 | E = -107.6541224475 | DE = -1.31e-11 | DW = 5.10e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.909 | E = -107.6541224475 | DE = -1.42e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 4.136 | E = -107.6541224475 | DE = 1.42e-12 | DW = 5.19e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.572 | E = -107.6541224475 | DE = -1.17e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.850 | E = -107.6541224475 | DE = 9.95e-13 | DW = 1.60e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.709 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.581 | E = -107.6541224475 | DE = -1.25e-12 | DW = 1.60e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 9.169 | E = -107.6541224475 | DE = -2.84e-14 | DW = 5.07e-20
DMRG energy = -107.654122447524500
We can also run DMRG in spin orbitals (SGF mode) using the restricted integrals, which will be much slower (for more realistic systems).
[5]:
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.n_sites = ncas
g2e = driver.unpack_g2e(g2e)
orb_sym = [orb_sym[i // 2] for i in range(len(orb_sym) * 2)]
n_sites = ncas * 2
driver.initialize_system(n_sites=n_sites, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error = 7.512924107430346e-14
integral cutoff error = 0.0
mpo terms = 2438
Build MPO | Nsites = 20 | Nterms = 2438 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.005
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.006
Site = 2 / 20 .. Mmpo = 45 DW = 0.00e+00 NNZ = 45 SPT = 0.9500 Tmvc = 0.001 T = 0.005
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 131 SPT = 0.9530 Tmvc = 0.001 T = 0.005
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 159 SPT = 0.9683 Tmvc = 0.001 T = 0.005
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 203 SPT = 0.9759 Tmvc = 0.001 T = 0.006
Site = 6 / 20 .. Mmpo = 125 DW = 0.00e+00 NNZ = 265 SPT = 0.9796 Tmvc = 0.001 T = 0.006
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 974 SPT = 0.9382 Tmvc = 0.001 T = 0.008
Site = 8 / 20 .. Mmpo = 151 DW = 0.00e+00 NNZ = 188 SPT = 0.9901 Tmvc = 0.001 T = 0.010
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 344 SPT = 0.9846 Tmvc = 0.001 T = 0.037
Site = 10 / 20 .. Mmpo = 177 DW = 0.00e+00 NNZ = 246 SPT = 0.9906 Tmvc = 0.001 T = 0.007
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 402 SPT = 0.9872 Tmvc = 0.001 T = 0.006
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 306 SPT = 0.9883 Tmvc = 0.001 T = 0.005
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 242 SPT = 0.9835 Tmvc = 0.000 T = 0.004
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 150 SPT = 0.9805 Tmvc = 0.000 T = 0.004
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 94 SPT = 0.9774 Tmvc = 0.000 T = 0.003
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 82 SPT = 0.9611 Tmvc = 0.000 T = 0.003
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 50 SPT = 0.9359 Tmvc = 0.000 T = 0.004
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 20 SPT = 0.8571 Tmvc = 0.000 T = 0.002
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.134 Tmvc-total = 0.012 MPO bond dimension = 178 MaxDW = 0.00e+00
NNZ = 3934 SIZE = 200866 SPT = 0.9804
Rank = 0 Ttotal = 0.186 MPO method = FastBipartite bond dimension = 178 NNZ = 3934 SIZE = 200866 SPT = 0.9804
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.861 | E = -107.6541216205 | DW = 7.62e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.471 | E = -107.6541223420 | DE = -7.21e-07 | DW = 6.45e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.030 | E = -107.6541224347 | DE = -9.27e-08 | DW = 7.64e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.597 | E = -107.6541224347 | DE = 8.81e-13 | DW = 6.15e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.707 | E = -107.6541224379 | DE = -3.17e-09 | DW = 6.46e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.003 | E = -107.6541224379 | DE = -3.37e-11 | DW = 7.34e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.306 | E = -107.6541224379 | DE = 0.00e+00 | DW = 6.46e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.088 | E = -107.6541224379 | DE = 7.11e-13 | DW = 7.10e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 7.514 | E = -107.6541224379 | DE = 2.84e-14 | DW = 8.86e-20
DMRG energy = -107.654122437940984
Read and Write FCIDUMP Files
Instead of generating integrals (h1e and g2e) using pyscf, we can also read these integrals from a FCIDUMP file (which can be generated using any of many other quantum chemistry packages) then perform DMRG. Additionally, we also provide methods to write the FCIDUMP file using the data in the h1e and g2e arrays.
After invoking driver.read_fcidump, the integrals and target state infomation can be obtained from driver.h1e, driver.g2e, driver.n_sites, etc.
[6]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
# write integrals to file
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
driver.write_fcidump(h1e, g2e, ecore, filename='N2.STO3G.FCIDUMP', h1e_symm=True, pg='d2h')
# read integrals from file
driver.read_fcidump(filename='N2.STO3G.FCIDUMP', pg='d2h')
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec,
spin=driver.spin, orb_sym=driver.orb_sym)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
mpo = driver.get_qc_mpo(h1e=driver.h1e, g2e=driver.g2e, ecore=driver.ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
symmetrize error = 2.3300000000000018e-14
integral symmetrize error = 0.0
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 13 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.000 T = 0.011
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.000 T = 0.018
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.032
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.001 T = 0.007
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.004
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.012
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.004
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.007
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.004
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Ttotal = 0.103 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.171 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.556 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.903 | E = -107.6541224475 | DE = -3.52e-12 | DW = 4.66e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.295 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.686 | E = -107.6541224475 | DE = -5.97e-13 | DW = 5.90e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.136 | E = -107.6541224475 | DE = -5.68e-14 | DW = 2.02e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.523 | E = -107.6541224475 | DE = 0.00e+00 | DW = 4.86e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.946 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.08e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.342 | E = -107.6541224475 | DE = -2.84e-14 | DW = 5.58e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.607 | E = -107.6541224475 | DE = 0.00e+00 | DW = 2.33e-20
DMRG energy = -107.654122447524614
The SZ Mode
Here we use get_uhf_integrals function to get the integrals.
[7]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_uhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error = 1.5366482610151817e-13
integral cutoff error = 0.0
mpo terms = 2778
Build MPO | Nsites = 10 | Nterms = 2778 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.001 T = 0.018
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 143 SPT = 0.9167 Tmvc = 0.001 T = 0.014
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 283 SPT = 0.9610 Tmvc = 0.001 T = 0.012
Site = 3 / 10 .. Mmpo = 138 DW = 0.00e+00 NNZ = 1023 SPT = 0.9326 Tmvc = 0.002 T = 0.015
Site = 4 / 10 .. Mmpo = 158 DW = 0.00e+00 NNZ = 535 SPT = 0.9755 Tmvc = 0.001 T = 0.012
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 463 SPT = 0.9842 Tmvc = 0.001 T = 0.007
Site = 6 / 10 .. Mmpo = 106 DW = 0.00e+00 NNZ = 415 SPT = 0.9790 Tmvc = 0.000 T = 0.006
Site = 7 / 10 .. Mmpo = 58 DW = 0.00e+00 NNZ = 163 SPT = 0.9735 Tmvc = 0.000 T = 0.004
Site = 8 / 10 .. Mmpo = 26 DW = 0.00e+00 NNZ = 87 SPT = 0.9423 Tmvc = 0.000 T = 0.003
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.093 Tmvc-total = 0.007 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 3164 SIZE = 102772 SPT = 0.9692
Rank = 0 Ttotal = 0.141 MPO method = FastBipartite bond dimension = 186 NNZ = 3164 SIZE = 102772 SPT = 0.9692
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.783 | E = -107.6541224475 | DW = 4.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.233 | E = -107.6541224475 | DE = -1.78e-11 | DW = 5.23e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.662 | E = -107.6541224475 | DE = -1.19e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.104 | E = -107.6541224475 | DE = 1.28e-12 | DW = 5.02e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.629 | E = -107.6541224475 | DE = -6.82e-13 | DW = 3.61e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.139 | E = -107.6541224475 | DE = 8.53e-13 | DW = 1.92e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.672 | E = -107.6541224475 | DE = 0.00e+00 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.173 | E = -107.6541224475 | DE = -1.68e-12 | DW = 2.25e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.529 | E = -107.6541224475 | DE = 2.84e-14 | DW = 7.36e-20
DMRG energy = -107.654122447524529
The SGF Mode
Here we use get_ghf_integrals function to get the integrals.
[8]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.GHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error = 2.1082787907764326e-13
integral cutoff error = 0.0
mpo terms = 5914
Build MPO | Nsites = 20 | Nterms = 5914 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.001 T = 0.007
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.009
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 49 SPT = 0.9479 Tmvc = 0.002 T = 0.010
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 251 SPT = 0.9139 Tmvc = 0.002 T = 0.013
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 273 SPT = 0.9456 Tmvc = 0.004 T = 0.019
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 357 SPT = 0.9576 Tmvc = 0.003 T = 0.013
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 563 SPT = 0.9580 Tmvc = 0.001 T = 0.012
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 2318 SPT = 0.8574 Tmvc = 0.001 T = 0.014
Site = 8 / 20 .. Mmpo = 155 DW = 0.00e+00 NNZ = 208 SPT = 0.9893 Tmvc = 0.002 T = 0.012
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 686 SPT = 0.9701 Tmvc = 0.001 T = 0.014
Site = 10 / 20 .. Mmpo = 181 DW = 0.00e+00 NNZ = 388 SPT = 0.9855 Tmvc = 0.001 T = 0.010
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 720 SPT = 0.9777 Tmvc = 0.001 T = 0.011
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 496 SPT = 0.9810 Tmvc = 0.001 T = 0.006
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 412 SPT = 0.9720 Tmvc = 0.000 T = 0.004
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 234 SPT = 0.9696 Tmvc = 0.000 T = 0.005
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 118 SPT = 0.9716 Tmvc = 0.000 T = 0.003
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 124 SPT = 0.9411 Tmvc = 0.000 T = 0.003
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 76 SPT = 0.9026 Tmvc = 0.000 T = 0.004
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 22 SPT = 0.8429 Tmvc = 0.000 T = 0.003
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Ttotal = 0.176 Tmvc-total = 0.023 MPO bond dimension = 181 MaxDW = 0.00e+00
NNZ = 7328 SIZE = 204350 SPT = 0.9641
Rank = 0 Ttotal = 0.269 MPO method = FastBipartite bond dimension = 181 NNZ = 7328 SIZE = 204350 SPT = 0.9641
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.012 | E = -107.6541220013 | DW = 8.06e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.754 | E = -107.6541223288 | DE = -3.27e-07 | DW = 7.52e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.474 | E = -107.6541224307 | DE = -1.02e-07 | DW = 8.06e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.210 | E = -107.6541224307 | DE = 1.14e-12 | DW = 6.96e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.149 | E = -107.6541224313 | DE = -5.22e-10 | DW = 8.65e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.121 | E = -107.6541224313 | DE = 1.65e-12 | DW = 7.88e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.140 | E = -107.6541224313 | DE = -5.68e-14 | DW = 8.65e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.097 | E = -107.6541224313 | DE = -2.47e-12 | DW = 7.25e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 7.676 | E = -107.6541224313 | DE = 0.00e+00 | DW = 8.25e-20
DMRG energy = -107.654122431266543
The SGB Mode
In this section, we try to solve the problem by first transfroming the model into a qubit (spin) model. The code will automatically use Jordan-Wigner transform to change the fermionic operators in the Hamiltonian into spin operators, before constructing the MPO.
To use the SGB mode for ab initio systems, remember to add the heis_twos=1 parameter (indicating the 1/2 spin at each site) in driver.initialize_system.
[9]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.GHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGB, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym, heis_twos=1)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error = 2.140864794887064e-13
integral cutoff error = 0.0
mpo terms = 5984
Build MPO | Nsites = 20 | Nterms = 5984 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.001 T = 0.006
Site = 1 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 19 SPT = 0.8643 Tmvc = 0.001 T = 0.009
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 49 SPT = 0.9479 Tmvc = 0.002 T = 0.011
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 251 SPT = 0.9139 Tmvc = 0.002 T = 0.014
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 273 SPT = 0.9456 Tmvc = 0.002 T = 0.014
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 357 SPT = 0.9576 Tmvc = 0.002 T = 0.014
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 563 SPT = 0.9580 Tmvc = 0.002 T = 0.016
Site = 7 / 20 .. Mmpo = 126 DW = 0.00e+00 NNZ = 2316 SPT = 0.8575 Tmvc = 0.002 T = 0.015
Site = 8 / 20 .. Mmpo = 155 DW = 0.00e+00 NNZ = 214 SPT = 0.9890 Tmvc = 0.001 T = 0.008
Site = 9 / 20 .. Mmpo = 148 DW = 0.00e+00 NNZ = 710 SPT = 0.9690 Tmvc = 0.001 T = 0.012
Site = 10 / 20 .. Mmpo = 181 DW = 0.00e+00 NNZ = 398 SPT = 0.9851 Tmvc = 0.001 T = 0.013
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 720 SPT = 0.9777 Tmvc = 0.001 T = 0.013
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 496 SPT = 0.9810 Tmvc = 0.001 T = 0.007
Site = 13 / 20 .. Mmpo = 100 DW = 0.00e+00 NNZ = 412 SPT = 0.9720 Tmvc = 0.000 T = 0.005
Site = 14 / 20 .. Mmpo = 77 DW = 0.00e+00 NNZ = 254 SPT = 0.9670 Tmvc = 0.000 T = 0.006
Site = 15 / 20 .. Mmpo = 54 DW = 0.00e+00 NNZ = 120 SPT = 0.9711 Tmvc = 0.000 T = 0.004
Site = 16 / 20 .. Mmpo = 39 DW = 0.00e+00 NNZ = 124 SPT = 0.9411 Tmvc = 0.000 T = 0.003
Site = 17 / 20 .. Mmpo = 20 DW = 0.00e+00 NNZ = 76 SPT = 0.9026 Tmvc = 0.000 T = 0.003
Site = 18 / 20 .. Mmpo = 7 DW = 0.00e+00 NNZ = 22 SPT = 0.8429 Tmvc = 0.000 T = 0.002
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.177 Tmvc-total = 0.018 MPO bond dimension = 181 MaxDW = 0.00e+00
NNZ = 7388 SIZE = 204350 SPT = 0.9638
Rank = 0 Ttotal = 0.249 MPO method = FastBipartite bond dimension = 181 NNZ = 7388 SIZE = 204350 SPT = 0.9638
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.079 | E = -107.6541211530 | DW = 6.36e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.838 | E = -107.6541223677 | DE = -1.21e-06 | DW = 5.89e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.571 | E = -107.6541224370 | DE = -6.93e-08 | DW = 6.35e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.284 | E = -107.6541224370 | DE = 4.26e-13 | DW = 5.71e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.581 | E = -107.6541224379 | DE = -9.81e-10 | DW = 9.19e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.311 | E = -107.6541224379 | DE = 1.96e-12 | DW = 6.93e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.757 | E = -107.6541224379 | DE = -5.68e-14 | DW = 9.19e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 8.724 | E = -107.6541224379 | DE = -1.88e-12 | DW = 7.39e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 9.323 | E = -107.6541224379 | DE = -2.84e-14 | DW = 8.26e-20
DMRG energy = -107.654122437940899
Relativistic DMRG
For relativistic DMRG, we use get_dhf_integrals function to get the integrals. We use the SGFCPX Mode in block2 to execute DMRG. Note that the integrals, MPO, and MPS will all contain complex numbers in this mode.
The symm_type parameter SymmetryTypes.SGFCPX can also be written as SymmetryTypes.SGF | SymmetryTypes.CPX.
[10]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.DHF(mol).set(with_gaunt=True, with_breit=True).run(conv_tol=1E-12)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_dhf_integrals(mf,
ncore=0, ncas=None, pg_symm=False)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGFCPX, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
integral symmetrize error = 0.0
integral cutoff error = 1.5501285354427146e-20
mpo terms = 44346
Build MPO | Nsites = 20 | Nterms = 44346 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.004 T = 0.017
Site = 1 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 23 SPT = 0.9087 Tmvc = 0.005 T = 0.025
Site = 2 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 404 SPT = 0.7710 Tmvc = 0.006 T = 0.023
Site = 3 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 592 SPT = 0.8795 Tmvc = 0.006 T = 0.028
Site = 4 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 932 SPT = 0.8768 Tmvc = 0.010 T = 0.035
Site = 5 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 1315 SPT = 0.8870 Tmvc = 0.008 T = 0.030
Site = 6 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 1722 SPT = 0.9024 Tmvc = 0.008 T = 0.033
Site = 7 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 2142 SPT = 0.9181 Tmvc = 0.012 T = 0.038
Site = 8 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 2597 SPT = 0.9315 Tmvc = 0.009 T = 0.038
Site = 9 / 20 .. Mmpo = 252 DW = 0.00e+00 NNZ = 2985 SPT = 0.9444 Tmvc = 0.011 T = 0.040
Site = 10 / 20 .. Mmpo = 213 DW = 0.00e+00 NNZ = 17958 SPT = 0.6654 Tmvc = 0.014 T = 0.068
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 2590 SPT = 0.9317 Tmvc = 0.003 T = 0.020
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 2184 SPT = 0.9165 Tmvc = 0.003 T = 0.016
Site = 13 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 1783 SPT = 0.8989 Tmvc = 0.002 T = 0.019
Site = 14 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 1397 SPT = 0.8800 Tmvc = 0.002 T = 0.013
Site = 15 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 1008 SPT = 0.8668 Tmvc = 0.001 T = 0.010
Site = 16 / 20 .. Mmpo = 63 DW = 0.00e+00 NNZ = 653 SPT = 0.8671 Tmvc = 0.001 T = 0.010
Site = 17 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 494 SPT = 0.7200 Tmvc = 0.001 T = 0.007
Site = 18 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 26 SPT = 0.8968 Tmvc = 0.000 T = 0.003
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.475 Tmvc-total = 0.105 MPO bond dimension = 252 MaxDW = 0.00e+00
NNZ = 40823 SIZE = 323082 SPT = 0.8736
Rank = 0 Ttotal = 0.567 MPO method = FastBipartite bond dimension = 252 NNZ = 40823 SIZE = 323082 SPT = 0.8736
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 23.177 | E = -107.6929206929 | DW = 4.40e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 35.078 | E = -107.6929209450 | DE = -2.52e-07 | DW = 8.35e-10
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 45.366 | E = -107.6929209492 | DE = -4.14e-09 | DW = 2.09e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 53.524 | E = -107.6929209492 | DE = -3.41e-13 | DW = 8.36e-10
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 71.358 | E = -107.6929209492 | DE = -4.26e-13 | DW = 6.57e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 86.911 | E = -107.6929209492 | DE = 2.84e-14 | DW = 1.18e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 103.723 | E = -107.6929209492 | DE = 2.84e-14 | DW = 7.77e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 119.169 | E = -107.6929209492 | DE = -2.84e-14 | DW = 1.19e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 128.099 | E = -107.6929209492 | DE = 2.84e-14 | DW = 8.33e-20
DMRG energy = -107.692920949170755
The LZ Mode
For diatomic molecules, we can set symmetry dooh in pyscf, and then use itg.lz_symm_adaptation to adapt the atomic orbitals for the LZ symmetry before running Hartree-Fock. Then we can use the LZ modes in block2 to perform DMRG.
The LZ mode can be combined with SU2, SZ or SGF spin symmetries, and the SAny prefix in SymmetryTypes. To activate the SAny prefix, the block2 code needs to be compiled with the -DUSE_SANY=ON option (this option is ON by default in the pip precompiled binaries). Optionally, when -DUSE_SANY=ON, one can also set -DUSE_SG=OFF -DUSE_SU2SZ=OFF to disable the normal SU2/SZ/SGF modes. One can use
SymmetryTypes.SAnySU2/SymmetryTypes.SAnySZ/SymmetryTypes.SAnySGF instead for normal symmetries (with some limitations) when -DUSE_SG=OFF -DUSE_SU2SZ=OFF is used.
With SU2 (spin-adapted DMRG):
[11]:
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="dooh", verbose=0)
mol.symm_orb, z_irrep, g_irrep = itg.lz_symm_adaptation(mol)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym_z = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=1, irrep_id=z_irrep)
print(orb_sym_z)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySU2LZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym_z, pg_irrep=0)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
[ 0 0 0 0 -1 1 0 1 -1 0]
integral symmetrize error = 2.0562981879479644e-15
integral cutoff error = 0.0
mpo terms = 1881
Build MPO | Nsites = 10 | Nterms = 1881 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.001 T = 0.031
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 107 SPT = 0.7752 Tmvc = 0.001 T = 0.054
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 189 SPT = 0.9007 Tmvc = 0.001 T = 0.063
Site = 3 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 822 SPT = 0.7776 Tmvc = 0.001 T = 0.063
Site = 4 / 10 .. Mmpo = 88 DW = 0.00e+00 NNZ = 154 SPT = 0.9735 Tmvc = 0.000 T = 0.033
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 318 SPT = 0.9616 Tmvc = 0.000 T = 0.031
Site = 6 / 10 .. Mmpo = 64 DW = 0.00e+00 NNZ = 236 SPT = 0.9608 Tmvc = 0.000 T = 0.023
Site = 7 / 10 .. Mmpo = 40 DW = 0.00e+00 NNZ = 145 SPT = 0.9434 Tmvc = 0.000 T = 0.016
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 49 SPT = 0.9125 Tmvc = 0.000 T = 0.009
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.007
Ttotal = 0.332 Tmvc-total = 0.005 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 2048 SIZE = 29320 SPT = 0.9302
Rank = 0 Ttotal = 0.402 MPO method = FastBipartite bond dimension = 94 NNZ = 2048 SIZE = 29320 SPT = 0.9302
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.679 | E = -107.6541224475 | DW = 2.06e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.030 | E = -107.6541224475 | DE = 1.19e-12 | DW = 2.77e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.383 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.06e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.721 | E = -107.6541224475 | DE = -8.67e-12 | DW = 3.29e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.225 | E = -107.6541224475 | DE = -2.84e-14 | DW = 4.20e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.894 | E = -107.6541224475 | DE = 5.68e-14 | DW = 2.00e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.682 | E = -107.6541224475 | DE = -2.84e-14 | DW = 4.90e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.379 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.69e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.809 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.21e-20
DMRG energy = -107.654122447523761
With SZ (non-spin-adapted DMRG):
[12]:
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="dooh", verbose=0)
mol.symm_orb, z_irrep, g_irrep = itg.lz_symm_adaptation(mol)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym_z = itg.get_uhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=1, irrep_id=z_irrep)
print(orb_sym_z)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySZLZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym_z, pg_irrep=0)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
[ 0 0 0 0 -1 1 0 1 -1 0]
integral symmetrize error = 8.030340066584701e-15
integral cutoff error = 0.0
mpo terms = 5231
Build MPO | Nsites = 10 | Nterms = 5231 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 30 SPT = 0.0000 Tmvc = 0.002 T = 0.080
Site = 1 / 10 .. Mmpo = 66 DW = 0.00e+00 NNZ = 271 SPT = 0.8631 Tmvc = 0.001 T = 0.135
Site = 2 / 10 .. Mmpo = 110 DW = 0.00e+00 NNZ = 480 SPT = 0.9339 Tmvc = 0.001 T = 0.136
Site = 3 / 10 .. Mmpo = 124 DW = 0.00e+00 NNZ = 2273 SPT = 0.8334 Tmvc = 0.001 T = 0.135
Site = 4 / 10 .. Mmpo = 171 DW = 0.00e+00 NNZ = 363 SPT = 0.9829 Tmvc = 0.001 T = 0.066
Site = 5 / 10 .. Mmpo = 186 DW = 0.00e+00 NNZ = 811 SPT = 0.9745 Tmvc = 0.001 T = 0.072
Site = 6 / 10 .. Mmpo = 126 DW = 0.00e+00 NNZ = 583 SPT = 0.9751 Tmvc = 0.001 T = 0.045
Site = 7 / 10 .. Mmpo = 82 DW = 0.00e+00 NNZ = 263 SPT = 0.9745 Tmvc = 0.000 T = 0.024
Site = 8 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 182 SPT = 0.9260 Tmvc = 0.000 T = 0.020
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 30 SPT = 0.0000 Tmvc = 0.000 T = 0.009
Ttotal = 0.722 Tmvc-total = 0.009 MPO bond dimension = 186 MaxDW = 0.00e+00
NNZ = 5286 SIZE = 112178 SPT = 0.9529
Rank = 0 Ttotal = 0.773 MPO method = FastBipartite bond dimension = 186 NNZ = 5286 SIZE = 112178 SPT = 0.9529
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.157 | E = -107.6541224475 | DW = 3.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.836 | E = -107.6541224475 | DE = -1.13e-11 | DW = 3.32e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.477 | E = -107.6541224455 | DE = 2.05e-09 | DW = 3.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.168 | E = -107.6541224455 | DE = -2.27e-13 | DW = 3.27e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.974 | E = -107.6541224455 | DE = -1.42e-13 | DW = 3.08e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.811 | E = -107.6541224455 | DE = 2.84e-14 | DW = 1.35e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.270 | E = -107.6541224455 | DE = -2.84e-14 | DW = 3.08e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.758 | E = -107.6541224455 | DE = 2.56e-13 | DW = 1.55e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 8.655 | E = -107.6541224455 | DE = -2.84e-14 | DW = 9.59e-20
DMRG energy = -107.654122445468332
With SGF (spin-orbital DMRG):
[13]:
from pyscf.scf.ghf_symm import GHF
from pyscf import symm, lib
from pyscf.scf import hf_symm
import scipy.linalg
import numpy as np
# fix pyscf 2.3.0 bug in ghf_symm for complex orbtials
def ghf_eig(self, h, s, symm_orb=None, irrep_id=None):
if symm_orb is None or irrep_id is None:
mol = self.mol
symm_orb = mol.symm_orb
irrep_id = mol.irrep_id
nirrep = len(symm_orb)
symm_orb = [scipy.linalg.block_diag(c, c) for c in symm_orb]
h = symm.symmetrize_matrix(h, symm_orb)
s = symm.symmetrize_matrix(s, symm_orb)
cs = []
es = []
orbsym = []
for ir in range(nirrep):
e, c = self._eigh(h[ir], s[ir])
cs.append(c)
es.append(e)
orbsym.append([irrep_id[ir]] * e.size)
e = np.hstack(es)
c = hf_symm.so2ao_mo_coeff(symm_orb, cs)
c = lib.tag_array(c, orbsym=np.hstack(orbsym))
return e, c
GHF.eig = ghf_eig
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="dooh", verbose=0)
mol.symm_orb, z_irrep, g_irrep = itg.lz_symm_adaptation(mol)
mf = scf.GHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym_z = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=1, irrep_id=z_irrep)
print(orb_sym_z)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SAnySGFLZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym_z, pg_irrep=0)
bond_dims = [250] * 4 + [500] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
[ 0 0 0 0 0 0 0 0 -1 -1 1 1 0 0 1 1 -1 -1 0 0]
integral symmetrize error = 9.305335353847714e-15
integral cutoff error = 5.457899320407648e-20
mpo terms = 12773
Build MPO | Nsites = 20 | Nterms = 12773 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.002 T = 0.098
Site = 1 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 25 SPT = 0.9008 Tmvc = 0.002 T = 0.170
Site = 2 / 20 .. Mmpo = 47 DW = 0.00e+00 NNZ = 309 SPT = 0.7652 Tmvc = 0.004 T = 0.245
Site = 3 / 20 .. Mmpo = 62 DW = 0.00e+00 NNZ = 357 SPT = 0.8775 Tmvc = 0.003 T = 0.264
Site = 4 / 20 .. Mmpo = 81 DW = 0.00e+00 NNZ = 514 SPT = 0.8977 Tmvc = 0.003 T = 0.280
Site = 5 / 20 .. Mmpo = 104 DW = 0.00e+00 NNZ = 667 SPT = 0.9208 Tmvc = 0.003 T = 0.308
Site = 6 / 20 .. Mmpo = 129 DW = 0.00e+00 NNZ = 2106 SPT = 0.8430 Tmvc = 0.003 T = 0.299
Site = 7 / 20 .. Mmpo = 118 DW = 0.00e+00 NNZ = 3960 SPT = 0.7399 Tmvc = 0.003 T = 0.249
Site = 8 / 20 .. Mmpo = 153 DW = 0.00e+00 NNZ = 250 SPT = 0.9862 Tmvc = 0.001 T = 0.127
Site = 9 / 20 .. Mmpo = 164 DW = 0.00e+00 NNZ = 690 SPT = 0.9725 Tmvc = 0.001 T = 0.134
Site = 10 / 20 .. Mmpo = 185 DW = 0.00e+00 NNZ = 1223 SPT = 0.9597 Tmvc = 0.002 T = 0.131
Site = 11 / 20 .. Mmpo = 178 DW = 0.00e+00 NNZ = 910 SPT = 0.9724 Tmvc = 0.001 T = 0.101
Site = 12 / 20 .. Mmpo = 147 DW = 0.00e+00 NNZ = 814 SPT = 0.9689 Tmvc = 0.001 T = 0.077
Site = 13 / 20 .. Mmpo = 120 DW = 0.00e+00 NNZ = 632 SPT = 0.9642 Tmvc = 0.001 T = 0.052
Site = 14 / 20 .. Mmpo = 97 DW = 0.00e+00 NNZ = 383 SPT = 0.9671 Tmvc = 0.001 T = 0.040
Site = 15 / 20 .. Mmpo = 78 DW = 0.00e+00 NNZ = 249 SPT = 0.9671 Tmvc = 0.000 T = 0.038
Site = 16 / 20 .. Mmpo = 57 DW = 0.00e+00 NNZ = 382 SPT = 0.9141 Tmvc = 0.000 T = 0.022
Site = 17 / 20 .. Mmpo = 28 DW = 0.00e+00 NNZ = 91 SPT = 0.9430 Tmvc = 0.000 T = 0.009
Site = 18 / 20 .. Mmpo = 9 DW = 0.00e+00 NNZ = 28 SPT = 0.8889 Tmvc = 0.000 T = 0.011
Site = 19 / 20 .. Mmpo = 1 DW = 0.00e+00 NNZ = 9 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Ttotal = 2.660 Tmvc-total = 0.031 MPO bond dimension = 185 MaxDW = 0.00e+00
NNZ = 13608 SIZE = 222306 SPT = 0.9388
Rank = 0 Ttotal = 2.764 MPO method = FastBipartite bond dimension = 185 NNZ = 13608 SIZE = 222306 SPT = 0.9388
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.655 | E = -107.6541219864 | DW = 1.34e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.673 | E = -107.6541222928 | DE = -3.06e-07 | DW = 4.92e-08
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 5.224 | E = -107.6541224418 | DE = -1.49e-07 | DW = 1.33e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 6.259 | E = -107.6541224418 | DE = -1.02e-12 | DW = 4.94e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 7.849 | E = -107.6541224449 | DE = -3.13e-09 | DW = 6.76e-12
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 9.282 | E = -107.6541224455 | DE = -5.68e-10 | DW = 5.49e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 10.873 | E = -107.6541224455 | DE = 5.68e-14 | DW = 6.76e-12
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 12.297 | E = -107.6541224455 | DE = 2.84e-14 | DW = 6.04e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 13.165 | E = -107.6541224455 | DE = 2.84e-14 | DW = 7.22e-20
DMRG energy = -107.654122445469270
Expectation and N-Particle Density Matrices
Once the optimized MPS is obtained, we can compute the expectation value on it, including its norm, the energy expectation, \(\langle S^2 \rangle\), N-particle density matrix, or any operator that can be constructed as an MPO.
In this example, we compute the triplet state.
[14]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
spin = 2
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
impo = driver.get_identity_mpo()
norm = driver.expectation(ket, impo, ket)
ener = driver.expectation(ket, mpo, ket)
print('Norm = %20.15f' % norm)
print('Energy expectation = %20.15f' % (ener / norm))
# <S^2> [ in spin-adapted mode this is always S(S+1) ]
ssq_mpo = driver.get_spin_square_mpo(iprint=0)
ssq = driver.expectation(ket, ssq_mpo, ket)
print('<S^2> expectation = %20.15f' % (ssq / norm))
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.557 | E = -106.9391328597 | DW = 3.38e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.833 | E = -106.9391328597 | DE = -3.24e-12 | DW = 9.75e-19
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.130 | E = -106.9391328597 | DE = 2.84e-14 | DW = 3.38e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.421 | E = -106.9391328597 | DE = 1.71e-13 | DW = 1.52e-18
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.755 | E = -106.9391328597 | DE = 0.00e+00 | DW = 1.12e-16
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.063 | E = -106.9391328597 | DE = 5.68e-14 | DW = 1.40e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.402 | E = -106.9391328597 | DE = 0.00e+00 | DW = 1.12e-16
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.709 | E = -106.9391328597 | DE = 0.00e+00 | DW = 1.38e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.933 | E = -106.9391328597 | DE = 0.00e+00 | DW = 3.39e-20
DMRG energy = -106.939132859666600
Norm = 1.000000000000000
Energy expectation = -106.939132859666614
<S^2> expectation = 2.000000000000000
We can also evaluate expectation of arbitray operator such as the occupancy in the first orbital
\[\hat{N}_0 = a^\dagger_{0\alpha} a_{0\alpha} + a^\dagger_{0\beta} a_{0\beta}
= \sqrt{2} \big(a_0^\dagger\big)^{[1/2]} \otimes_{[0]} \big(a_0\big)^{[1/2]}\]
[15]:
b = driver.expr_builder()
b.add_term("(C+D)0", [0, 0], np.sqrt(2))
n_mpo = driver.get_mpo(b.finalize(), iprint=0)
n_0 = driver.expectation(ket, n_mpo, ket)
print('N0 expectation = %20.15f' % (n_0 / norm))
N0 expectation = 1.999995824361892
We can then verify this number using 1PDM:
[16]:
pdm1 = driver.get_1pdm(ket)
print('N0 expectation from 1pdm = %20.15f' % pdm1[0, 0])
N0 expectation from 1pdm = 1.999995824361892
We can compute the 3PDM and compare the result with the FCI 3PDM. Note that in pyscf the 3PDM is defined as
\[\mathrm{DM}_{ijklmn} := \langle E_{ij} E_{kl} E_{mn} \rangle\]
So we have to use the same convention in block2 by setting the npdm_expr parameter in block2 to ((C+D)0+((C+D)0+(C+D)0)0)0.
[17]:
pdm3_b2 = driver.get_3pdm(ket, iprint=0, npdm_expr="((C+D)0+((C+D)0+(C+D)0)0)0")
from pyscf import fci
mx = fci.addons.fix_spin_(fci.FCI(mf), ss=2)
mx.kernel(h1e, g2e, ncas, nelec=n_elec, nroots=3, tol=1E-12)
print(mx.e_tot)
pdm3_fci = fci.rdm.make_dm123('FCI3pdm_kern_sf', mx.ci[0], mx.ci[0], ncas, n_elec)[2]
print('diff = ', np.linalg.norm(pdm3_fci - pdm3_b2))
[-106.93913286 -106.85412245 -106.70055113]
diff = 4.622993128790701e-06
Extract CSF and Determinant Coefficients
We can extract CSF (or determinant) coefficients from the spin-adapted MPS (or the non-spin-adapted MPS). The algorithm can compute all CSF or determinant with the absolute value of the coefficient above a threshold (called cutoff). The square of coefficient is the probability (weight) of the CSF or determinant.
Extracting CSF coefficients from spin-adapted MPS in the SU2 mode:
[18]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.639 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.016 | E = -107.6541224475 | DE = -3.38e-12 | DW = 3.81e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.407 | E = -107.6541224475 | DE = 5.68e-14 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.777 | E = -107.6541224475 | DE = -5.68e-13 | DW = 4.69e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.209 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.73e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.604 | E = -107.6541224475 | DE = 2.84e-14 | DW = 5.35e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.011 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.66e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.376 | E = -107.6541224475 | DE = 0.00e+00 | DW = 5.98e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 3.647 | E = -107.6541224475 | DE = 2.84e-14 | DW = 5.56e-20
DMRG energy = -107.654122447524614
dtrie finished 0.01697694599999977
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368811359554
CSF 0 2222222000 = 0.957506528669401
CSF 1 2222202200 = -0.131287867807394
CSF 2 2222022020 = -0.131287794595557
CSF 3 2222+-2+-0 = 0.118594301088608
CSF 4 2222++2--0 = -0.084968224982022
CSF 5 2220222020 = -0.054710635444627
CSF 6 2220222200 = -0.054710465272207
CSF 7 2222+2+0-- = 0.053881241407894
CSF 8 22222++-0- = 0.053881215166769
Extracting determinant coefficients from spin-adapted MPS in the SZ mode:
[19]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_uhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.373 | E = -107.6541224475 | DW = 4.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.089 | E = -107.6541224475 | DE = -2.07e-12 | DW = 5.05e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.609 | E = -107.6541224475 | DE = -1.34e-12 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.051 | E = -107.6541224475 | DE = 1.28e-12 | DW = 5.21e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.572 | E = -107.6541224475 | DE = -9.95e-13 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.085 | E = -107.6541224475 | DE = 8.81e-13 | DW = 1.37e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.625 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.142 | E = -107.6541224475 | DE = -1.28e-12 | DW = 1.52e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 5.505 | E = -107.6541224475 | DE = 0.00e+00 | DW = 8.34e-20
DMRG energy = -107.654122447524443
dtrie finished 0.02674672699998837
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331638354918
DET 0 2222222000 = 0.957506568620143
DET 1 2222022020 = -0.131287835040785
DET 2 2222202200 = -0.131287814797523
DET 3 2222ab2ab0 = 0.083825279692261
DET 4 2222ba2ba0 = 0.083825279691480
DET 5 2220222200 = -0.054710476471258
DET 6 2220222020 = -0.054710439530678
Construct MPS from CSFs or Determinants
If we know important CSFs or determinants in the state and their coefficients, we can also use this information to construct MPS, and this can be used as an initial guess and further optimized. Note that this initial guess can generate very good initial energies, but if the given CSFs have many doubly occupied and empty orbitals, the MPS will very likely optimize to a local minima.
In the SU2 mode:
[20]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
mps = driver.get_mps_from_csf_coefficients(csfs, coeffs, tag="CMPS", dot=2)
impo = driver.get_identity_mpo()
print(driver.expectation(mps, impo, mps))
print(driver.expectation(mps, mpo, mps) / driver.expectation(mps, impo, mps))
energy = driver.dmrg(mpo, mps, n_sweeps=5, bond_dims=[500] * 5, noises=[1E-5],
thrds=[1E-10] * 5, iprint=1)
print('Ground state energy = %20.15f' % energy)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.496 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.733 | E = -107.6541224475 | DE = -1.24e-11 | DW = 4.89e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.990 | E = -107.6541224475 | DE = 2.84e-14 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.208 | E = -107.6541224475 | DE = -1.36e-12 | DW = 3.97e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.488 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.48e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.713 | E = -107.6541224475 | DE = -2.84e-14 | DW = 6.11e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.974 | E = -107.6541224475 | DE = 5.68e-14 | DW = 2.41e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.209 | E = -107.6541224475 | DE = -2.84e-14 | DW = 6.14e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.482 | E = -107.6541224475 | DE = 5.68e-14 | DW = 2.91e-20
DMRG energy = -107.654122447524358
dtrie finished 0.016507582999963688
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368809657036
CSF 0 2222222000 = 0.957506521875127
CSF 1 2222202200 = -0.131287813613834
CSF 2 2222022020 = -0.131287676854091
CSF 3 2222+-2+-0 = 0.118594392287595
CSF 4 2222++2--0 = -0.084968279487231
CSF 5 2220222020 = -0.054710673444202
CSF 6 2220222200 = -0.054710596031681
CSF 7 2222+2+0-- = 0.053881289057939
CSF 8 22222++-0- = 0.053881233355398
0.9843688096570361
-107.6155366353171
Sweep = 0 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.054 | E = -107.6285084535 | DW = 3.52e-21
Sweep = 1 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.109 | E = -107.6288109091 | DE = -3.02e-04 | DW = 9.04e-21
Sweep = 2 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.164 | E = -107.6289019710 | DE = -9.11e-05 | DW = 4.57e-21
Sweep = 3 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.217 | E = -107.6289019713 | DE = -3.35e-10 | DW = 2.07e-20
Ground state energy = -107.628901971339317
In the SZ mode:
[21]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_uhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
dets, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
mps = driver.get_mps_from_csf_coefficients(dets, coeffs, tag="CMPS", dot=2)
impo = driver.get_identity_mpo()
print(driver.expectation(mps, impo, mps))
print(driver.expectation(mps, mpo, mps) / driver.expectation(mps, impo, mps))
energy = driver.dmrg(mpo, mps, n_sweeps=5, bond_dims=[500] * 5, noises=[1E-5],
thrds=[1E-10] * 5, iprint=1)
print('Ground state energy = %20.15f' % energy)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.809 | E = -107.6541224475 | DW = 4.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.294 | E = -107.6541224475 | DE = -5.80e-12 | DW = 5.17e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.713 | E = -107.6541224475 | DE = -9.38e-13 | DW = 4.15e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.170 | E = -107.6541224475 | DE = 1.36e-12 | DW = 5.37e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.687 | E = -107.6541224475 | DE = -1.36e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.212 | E = -107.6541224475 | DE = 1.05e-12 | DW = 1.85e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.737 | E = -107.6541224475 | DE = -2.84e-14 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.252 | E = -107.6541224475 | DE = -1.14e-12 | DW = 1.77e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.598 | E = -107.6541224475 | DE = 0.00e+00 | DW = 8.68e-20
DMRG energy = -107.654122447524529
dtrie finished 0.03569589099998893
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331637808218
DET 0 2222222000 = -0.957506566009713
DET 1 2222022020 = 0.131287836810967
DET 2 2222202200 = 0.131287818787713
DET 3 2222ab2ab0 = -0.083825284176507
DET 4 2222ba2ba0 = -0.083825284073794
DET 5 2220222200 = 0.054710477807739
DET 6 2220222020 = 0.054710451475846
0.971331637808218
-107.59447325662364
Sweep = 0 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.047 | E = -107.6101211360 | DW = 6.22e-21
Sweep = 1 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.104 | E = -107.6102787809 | DE = -1.58e-04 | DW = 8.84e-21
Sweep = 2 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.163 | E = -107.6104759466 | DE = -1.97e-04 | DW = 1.22e-20
Sweep = 3 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.218 | E = -107.6104759470 | DE = -4.45e-10 | DW = 1.18e-20
Ground state energy = -107.610475947024796
MPS Initial Guess from Occupancies
We can also construct the MPS using an estimate of the occupancy information at each site from a cheaper method (such as CCSD). We can use this information to construct the initial quantum number distribution in the MPS.
Note that this initial guess can generate very good energies in the first few sweeps, but if the given occupancies have many doubly occupied and empty orbitals, the MPS will very likely optimize to a local minima. One can shift the occupancies into the equal probability occupancies (for example, setting bias=0.4 in the example below) to randomize the initial guess.
In the SU2 mode:
[22]:
from pyscf import gto, scf, cc
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
bias = 0.1 # make it more random
occs = np.diag(cc.CCSD(mf).run().make_rdm1())
occs = occs + bias * (occs < 1) - bias * (occs > 1)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, occs=occs, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.044 | E = -107.5033999317 | DW = 2.59e-21
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.135 | E = -107.5830591524 | DE = -7.97e-02 | DW = 3.51e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.217 | E = -107.5836402922 | DE = -5.81e-04 | DW = 9.17e-21
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.294 | E = -107.5867932578 | DE = -3.15e-03 | DW = 1.12e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.371 | E = -107.5868061724 | DE = -1.29e-05 | DW = 1.17e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.443 | E = -107.5868152464 | DE = -9.07e-06 | DW = 1.23e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.517 | E = -107.5868152464 | DE = -5.68e-14 | DW = 1.45e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.589 | E = -107.5868152464 | DE = 2.84e-14 | DW = 8.04e-21
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 0.655 | E = -107.5868152464 | DE = 2.84e-14 | DW = 1.69e-20
DMRG energy = -107.586815246376730
In the SZ mode:
[23]:
from pyscf import gto, scf, cc
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.UHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_uhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
bias = 0.1 # make it more random
occs = np.diag(np.sum(cc.UCCSD(mf).run().make_rdm1(), axis=0))
occs = occs + bias * (occs < 1) - bias * (occs > 1)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, occs=occs, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.097 | E = -107.5033999316 | DW = 3.37e-21
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.267 | E = -107.5830606948 | DE = -7.97e-02 | DW = 7.61e-21
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.376 | E = -107.5836402307 | DE = -5.80e-04 | DW = 1.27e-20
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.505 | E = -107.5867932582 | DE = -3.15e-03 | DW = 1.21e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.624 | E = -107.5868061673 | DE = -1.29e-05 | DW = 1.84e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.847 | E = -107.5868152464 | DE = -9.08e-06 | DW = 1.03e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.008 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.88e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.132 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.76e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 1.218 | E = -107.5868152464 | DE = 0.00e+00 | DW = 1.33e-20
DMRG energy = -107.586815246376730
In the SGF mode:
[24]:
from pyscf import gto, scf, cc
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.GHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_ghf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
bias = 0.25 # make it more random
occs = np.sum(cc.GCCSD(mf).run().make_rdm1(), axis=0)
occs = occs + bias * (occs < 0.5) - bias * (occs > 0.5)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SGF, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, occs=occs, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.853 | E = -107.6538893531 | DW = 2.24e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.548 | E = -107.6539226742 | DE = -3.33e-05 | DW = 8.93e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.168 | E = -107.6539462697 | DE = -2.36e-05 | DW = 5.04e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.808 | E = -107.6539506829 | DE = -4.41e-06 | DW = 1.95e-08
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.596 | E = -107.6539506829 | DE = 3.24e-12 | DW = 1.84e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.386 | E = -107.6539506829 | DE = -4.12e-12 | DW = 5.67e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.236 | E = -107.6539506829 | DE = -2.84e-14 | DW = 1.84e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 6.014 | E = -107.6539506829 | DE = 2.05e-12 | DW = 4.05e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.497 | E = -107.6539506829 | DE = 0.00e+00 | DW = 7.66e-20
DMRG energy = -107.653950682884570
Change from SU2 MPS to SZ MPS
We can also transform the spin-adapted MPS generated in the SU2 mode to the non-spin-adapted MPS, which can be used in the SZ mode. After obtaining the non-spin-adapted MPS, one need to redo driver.initialize_system, driver.get_qc_mpo, etc. to make sure every object is now represented in the SZ mode, then you can operate on the SZ non-spin-adapted MPS.
In the following example, we first compute the spin-adapted MPS, then translate it into the non-spin-adapted MPS to extract the determinant coefficients.
[25]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
csfs, coeffs = driver.get_csf_coefficients(ket, cutoff=0.05, iprint=1)
zket = driver.mps_change_to_sz(ket, "ZKET")
driver.symm_type = SymmetryTypes.SZ
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
impo = driver.get_identity_mpo()
print(driver.expectation(zket, mpo, zket) / driver.expectation(zket, impo, zket))
csfs, vals = driver.get_csf_coefficients(zket, cutoff=0.05, iprint=1)
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.902 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.281 | E = -107.6541224475 | DE = -2.87e-12 | DW = 4.49e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.682 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.068 | E = -107.6541224475 | DE = -5.12e-13 | DW = 5.07e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.610 | E = -107.6541224475 | DE = -5.68e-14 | DW = 2.62e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.019 | E = -107.6541224475 | DE = 2.84e-14 | DW = 4.59e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.435 | E = -107.6541224475 | DE = 2.84e-14 | DW = 2.01e-20
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.808 | E = -107.6541224475 | DE = -5.68e-14 | DW = 3.69e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.095 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.68e-20
DMRG energy = -107.654122447524671
dtrie finished 0.026360811000017748
Number of CSF = 9 (cutoff = 0.05)
Sum of weights of included CSF = 0.984368806039259
CSF 0 2222222000 = -0.957506495595307
CSF 1 2222022020 = 0.131287987564377
CSF 2 2222202200 = 0.131287965813458
CSF 3 2222+-2+-0 = -0.118594461577159
CSF 4 2222++2--0 = 0.084967942537817
CSF 5 2220222020 = 0.054710523921704
CSF 6 2220222200 = 0.054710517406346
CSF 7 2222+2+0-- = -0.053881222618422
CSF 8 22222++-0- = -0.053881215803966
-107.65412244752456
dtrie finished 0.0680388100000755
Number of DET = 7 (cutoff = 0.05)
Sum of weights of included DET = 0.971331619872548
DET 0 2222222000 = -0.957506495595307
DET 1 2222022020 = 0.131287987564377
DET 2 2222202200 = 0.131287965813458
DET 3 2222ab2ab0 = -0.083825363036928
DET 4 2222ba2ba0 = -0.083825363036928
DET 5 2220222020 = 0.054710523921704
DET 6 2220222200 = 0.054710517406346
Change between Real and Complex MPS
We can also change between the MPS with complex numbers and the MPS with real numbers. For complex MPS to real MPS, the imaginary part will be discarded (and the norm of the transformed MPS may decrease). This may be useful when you do a ground state calculation in the real domain and then do the real time evolution in the complex domain.
From real to complex:
[26]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
impo = driver.get_identity_mpo()
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
print(driver.expectation(ket, impo, ket))
print(driver.expectation(ket, mpo, ket) / driver.expectation(ket, impo, ket))
zket = driver.mps_change_complex(ket, "ZKET")
driver.symm_type = driver.symm_type ^ SymmetryTypes.CPX
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
impo = driver.get_identity_mpo()
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, integral_cutoff=1E-8, iprint=1)
print(driver.expectation(zket, impo, zket))
print(driver.expectation(zket, mpo, zket) / driver.expectation(zket, impo, zket))
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.804 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.144 | E = -107.6541224475 | DE = -3.95e-12 | DW = 4.49e-20
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.379 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.601 | E = -107.6541224475 | DE = -7.67e-13 | DW = 5.89e-20
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.868 | E = -107.6541224475 | DE = -2.84e-14 | DW = 2.31e-20
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.110 | E = -107.6541224475 | DE = -2.84e-14 | DW = 6.59e-20
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.366 | E = -107.6541224475 | DE = 0.00e+00 | DW = 8.80e-21
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.606 | E = -107.6541224475 | DE = -5.68e-14 | DW = 6.06e-20
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 2.780 | E = -107.6541224475 | DE = 2.84e-14 | DW = 4.42e-20
DMRG energy = -107.654122447524671
1.0
-107.65412244752456
integral symmetrize error = 1.76819336753967e-14
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 13 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.001 T = 0.007
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.001 T = 0.005
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.006
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.000 T = 0.006
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.007
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.005
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.008
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.003
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.004
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.054 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.108 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
(0.9999999999999999+0j)
(-107.65412244752457+0j)
From complex to real:
[27]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2 | SymmetryTypes.CPX, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
impo = driver.get_identity_mpo()
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
print(driver.expectation(ket, impo, ket))
print(driver.expectation(ket, mpo, ket) / driver.expectation(ket, impo, ket))
rket = driver.mps_change_complex(ket, "rket")
driver.symm_type = driver.symm_type ^ SymmetryTypes.CPX
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
impo = driver.get_identity_mpo()
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, integral_cutoff=1E-8, iprint=1)
print(driver.expectation(rket, impo, rket))
print(driver.expectation(rket, mpo, rket) / driver.expectation(rket, impo, rket))
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.189 | E = -107.6541224475 | DW = 1.87e-10
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.240 | E = -107.6541224475 | DE = -1.41e-11 | DW = 2.95e-19
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.105 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.87e-10
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 3.655 | E = -107.6541224475 | DE = -2.05e-12 | DW = 1.13e-19
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.299 | E = -107.6541224475 | DE = 0.00e+00 | DW = 1.31e-19
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.784 | E = -107.6541224475 | DE = -2.84e-14 | DW = 1.06e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.305 | E = -107.6541224475 | DE = 5.68e-14 | DW = 1.24e-19
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.770 | E = -107.6541224475 | DE = -5.68e-14 | DW = 1.04e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 6.086 | E = -107.6541224475 | DE = 0.00e+00 | DW = 4.25e-20
DMRG energy = -107.654122447523960
(1.0000000000000002-7.905175135209388e-18j)
(-107.6541224475238-8.610748230937272e-16j)
integral symmetrize error = 1.6922854083945986e-14
integral cutoff error = 0.0
mpo terms = 1030
Build MPO | Nsites = 10 | Nterms = 1030 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 10 .. Mmpo = 13 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.000 T = 0.004
Site = 1 / 10 .. Mmpo = 34 DW = 0.00e+00 NNZ = 63 SPT = 0.8575 Tmvc = 0.001 T = 0.004
Site = 2 / 10 .. Mmpo = 56 DW = 0.00e+00 NNZ = 121 SPT = 0.9364 Tmvc = 0.000 T = 0.004
Site = 3 / 10 .. Mmpo = 74 DW = 0.00e+00 NNZ = 373 SPT = 0.9100 Tmvc = 0.000 T = 0.005
Site = 4 / 10 .. Mmpo = 80 DW = 0.00e+00 NNZ = 269 SPT = 0.9546 Tmvc = 0.000 T = 0.004
Site = 5 / 10 .. Mmpo = 94 DW = 0.00e+00 NNZ = 169 SPT = 0.9775 Tmvc = 0.000 T = 0.004
Site = 6 / 10 .. Mmpo = 54 DW = 0.00e+00 NNZ = 181 SPT = 0.9643 Tmvc = 0.000 T = 0.003
Site = 7 / 10 .. Mmpo = 30 DW = 0.00e+00 NNZ = 73 SPT = 0.9549 Tmvc = 0.000 T = 0.003
Site = 8 / 10 .. Mmpo = 14 DW = 0.00e+00 NNZ = 41 SPT = 0.9024 Tmvc = 0.000 T = 0.002
Site = 9 / 10 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.036 Tmvc-total = 0.003 MPO bond dimension = 94 MaxDW = 0.00e+00
NNZ = 1317 SIZE = 27073 SPT = 0.9514
Rank = 0 Ttotal = 0.064 MPO method = FastBipartite bond dimension = 94 NNZ = 1317 SIZE = 27073 SPT = 0.9514
0.9626534962131837
-107.65412244752437
MPS Bipartite Entanglement
We can get the the bipartite entanglement \(S_k=-\sum_i \Lambda_k^2 \log \Lambda_k^2\) at each virtual bond (at site \(k\)) in MPS in the SZ mode, where \(\Lambda_k\) are singular values in the bond at site \(k\).
[28]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="GS", bond_dim=250, nroots=1)
energy = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises,
thrds=thrds, iprint=1)
print('DMRG energy = %20.15f' % energy)
bip_ent = driver.get_bipartite_entanglement()
import matplotlib.pyplot as plt
plt.plot(np.arange(len(bip_ent)), bip_ent, linestyle='-', marker='o',
mfc='white', mec="#7FB685", color="#7FB685")
plt.xlabel("site index $k$")
plt.ylabel("bipartite entanglement $S_k$")
plt.show()
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.769 | E = -107.6541224475 | DW = 4.14e-08
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.254 | E = -107.6541224475 | DE = -1.43e-11 | DW = 5.08e-09
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 1.670 | E = -107.6541224475 | DE = -5.12e-13 | DW = 4.14e-08
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 2.130 | E = -107.6541224475 | DE = 1.31e-12 | DW = 5.21e-09
Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.655 | E = -107.6541224475 | DE = -1.14e-12 | DW = 3.60e-11
Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.189 | E = -107.6541224475 | DE = 9.38e-13 | DW = 1.75e-19
Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.722 | E = -107.6541224475 | DE = 2.84e-14 | DW = 3.60e-11
Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.245 | E = -107.6541224475 | DE = -1.22e-12 | DW = 1.54e-19
Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 4.595 | E = -107.6541224475 | DE = 2.84e-14 | DW = 7.71e-20
DMRG energy = -107.654122447524443
Orbital Entropy and Mutual Information
For the optimized MPS in the SZ mode, we can compute the 1- and 2- orbital density matrices and mutual information for pairs of orbitals.
[29]:
odm1 = driver.get_orbital_entropies(ket, orb_type=1)
odm2 = driver.get_orbital_entropies(ket, orb_type=2)
minfo = 0.5 * (odm1[:, None] + odm1[None, :] - odm2) * (1 - np.identity(len(odm1)))
import matplotlib.pyplot as plt
plt.matshow(minfo)
[29]:
<matplotlib.image.AxesImage at 0x7d0b6521bd30>
Excited States
To obtain the excited states and their energies, we can perform DMRG for a state-averged MPS, optionally followed by a state-specific refinement.
[30]:
from pyscf import gto, scf
mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="sto3g", symmetry="d2h", verbose=0)
mf = scf.RHF(mol).run(conv_tol=1E-14)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf,
ncore=0, ncas=None, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
mpo = driver.get_qc_mpo(h1e=h1e, g2e=g2e, ecore=ecore, iprint=0)
ket = driver.get_random_mps(tag="KET", bond_dim=100, nroots=3)
energies = driver.dmrg(mpo, ket, n_sweeps=10, bond_dims=[100], noises=[1e-5] * 4 + [0],
thrds=[1e-10] * 8, iprint=1)
print('State-averaged MPS energies = [%s]' % " ".join("%20.15f" % x for x in energies))
kets = [driver.split_mps(ket, ir, tag="KET-%d" % ir) for ir in range(ket.nroots)]
for ir in range(ket.nroots):
energy = driver.dmrg(mpo, kets[ir], n_sweeps=10, bond_dims=[200], noises=[1e-5] * 4 + [0],
thrds=[1e-10] * 8, iprint=0, proj_weights=[5.0] * ir, proj_mpss=kets[:ir])
print('State-specific MPS E[%d] = %20.15f' % (ir, energy))
Sweep = 0 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 2.037 | E[ 3] = -107.6541193407 -107.0314407783 -106.9595996333 | DW = 8.51e-05
Sweep = 1 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 3.220 | E[ 3] = -107.6541193407 -107.0314407783 -106.9595996332 | DE = 7.11e-12 | DW = 1.62e-04
Sweep = 2 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 4.616 | E[ 3] = -107.6541138874 -107.0314486179 -106.9594564026 | DE = 1.43e-04 | DW = 7.20e-05
Sweep = 3 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 5.853 | E[ 3] = -107.6541138874 -107.0314486178 -106.9594564026 | DE = -3.97e-11 | DW = 1.47e-04
Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10
Time elapsed = 6.621 | E[ 3] = -107.6540972369 -107.0314487271 -106.9594391421 | DE = 1.73e-05 | DW = 6.81e-05
Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10
Time elapsed = 7.410 | E[ 3] = -107.6540972369 -107.0314487271 -106.9594391422 | DE = -1.72e-11 | DW = 1.43e-04
State-averaged MPS energies = [-107.654097236905088 -107.031448727142006 -106.959439142163390]
State-specific MPS E[0] = -107.654122447496988
State-specific MPS E[1] = -107.031449471612873
State-specific MPS E[2] = -106.959626154678844