DDMRG++ for Green’s Function
References:
Ronca, E., Li, Z., Jimenez-Hoyos, C. A., Chan, G. K. L. Time-step targeting time-dependent and dynamical density matrix renormalization group algorithms with ab initio Hamiltonians. Journal of Chemical Theory and Computation 2017, 13, 5560-5571.
The following example calculate the Green’s function for H10 (STO6G):
\[G_{ij}(\omega) = \langle \Psi_0 | a_j^\dagger \frac{1}{\omega + \hat{H}_0 - E_0 + i \eta} a_i |\Psi_0\rangle\]
where \(|\Psi_0\rangle\) is the ground-state, \(i = j = 4\) (isite
), \(\omega = -0.17, \eta = 0.05\).
import time
import numpy as np
from pyblock3.algebra.mpe import MPE
from pyblock3.hamiltonian import Hamiltonian
from pyblock3.fcidump import FCIDUMP
from pyblock3.symbolic.expr import OpElement, OpNames
from pyblock3.algebra.symmetry import SZ
np.random.seed(1000)
fd = 'data/H10.STO6G.R1.8.FCIDUMP'
ket_bond_dim = 500
bra_bond_dim = 750
hamil = Hamiltonian(FCIDUMP(pg='d2h').read(fd), flat=True)
tx = time.perf_counter()
mpo = hamil.build_qc_mpo()
print('MPO (NC) = ', mpo.show_bond_dims())
print('build mpo time = ', time.perf_counter() - tx)
tx = time.perf_counter()
mpo, _ = mpo.compress(left=True, cutoff=1E-9, norm_cutoff=1E-9)
print('MPO (compressed) = ', mpo.show_bond_dims())
print('compress mpo time = ', time.perf_counter() - tx)
mps = hamil.build_mps(ket_bond_dim, occ=occ)
print('MPS = ', mps.show_bond_dims())
bdims = [500]
noises = [1E-4, 1E-5, 1E-6, 0]
davthrds = None
dmrg = MPE(mps, mpo, mps).dmrg(bdims=bdims, noises=noises,
dav_thrds=davthrds, iprint=2, n_sweeps=20, tol=1E-12)
ener = dmrg.energies[-1]
print("Energy = %20.12f" % ener)
isite = 4
mpo.const -= ener
omega, eta = -0.17, 0.05
dop = OpElement(OpNames.D, (isite, 0), q_label=SZ(-1, -1, hamil.orb_sym[isite]))
bra = hamil.build_mps(bra_bond_dim, target=SZ.to_flat(
dop.q_label + SZ.from_flat(hamil.target)))
dmpo = hamil.build_site_mpo(dop)
print('DMPO = ', dmpo.show_bond_dims())
MPE(bra, dmpo, mps).linear(bdims=[bra_bond_dim], noises=noises,
cg_thrds=davthrds, iprint=2, n_sweeps=20, tol=1E-12)
np.random.seed(0)
gbra = hamil.build_mps(bra_bond_dim, target=SZ.to_flat(
dop.q_label + SZ.from_flat(hamil.target)))
print('GFMPO = ', mpo.show_bond_dims())
print(MPE(bra, dmpo, mps).greens_function(mpo, omega, eta, bdims=[bra_bond_dim], noises=noises,
cg_thrds=[1E-4] * 10, iprint=2, n_sweeps=10, tol=1E-4))