MPO Construction
References:
Knowles, P. J., Handy, N. C. A determinant based full configuration interaction program. Computer Physics Communications 1989, 54, 75-83.
Chan, G. K.-L.; Keselman, A.; Nakatani, N.; Li, Z.; White, S. R. Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms. The Journal of Chemical Physics 2016, 145, 014102.
Ren, J., Li, W., Jiang, T., & Shuai, Z. A general automatic method for optimal construction of matrix product operators using bipartite graph theory. The Journal of Chemical Physics 2020, 153, 084118.
Hubig, C., McCulloch, I. P., & Schollwöck, U. Generic construction of efficient matrix product operators. Physical Review B 2017, 95, 035129.
From FCIDUMP
MPO can be constructed from a FCIDUMP file:
from pyblock3.hamiltonian import Hamiltonian
from pyblock3.fcidump import FCIDUMP
fd = 'data/H8.STO6G.R1.8.FCIDUMP'
hamil = Hamiltonian(FCIDUMP(pg='d2h').read(fd), flat=False)
mpo = hamil.build_qc_mpo()
This will build an MPS object (representing MPO) using a list of FermionTensor.
If flat parameter is set to True in Hamiltonian, the code will use more efficient C++ code for building
MPO, and the resulting MPO is an MPS object with a list of FlatFermionTensor included.
from pyblock3.hamiltonian import Hamiltonian
from pyblock3.fcidump import FCIDUMP
fd = 'data/H8.STO6G.R1.8.FCIDUMP'
hamil = Hamiltonian(FCIDUMP(pg='d2h').read(fd), flat=True)
mpo = hamil.build_qc_mpo()
One can also use mpo.to_flat() to transform a FermionTensor MPO to a FlatFermionTensor MPO.
From Hamiltonian expression (pure python)
A slower but more general way of build MPO is from Hamiltonian expression.
One can set the explicit expression for Hamiltonian for Hubbard model and quantum chemistry:
def build_hubbard(u=2, t=1, n=8, cutoff=1E-9):
fcidump = FCIDUMP(pg='c1', n_sites=n, n_elec=n, twos=0, ipg=0, orb_sym=[0] * n)
hamil = Hamiltonian(fcidump, flat=False)
def generate_terms(n_sites, c, d):
for i in range(0, n_sites):
for s in [0, 1]:
if i - 1 >= 0:
yield t * c[i, s] * d[i - 1, s]
if i + 1 < n_sites:
yield t * c[i, s] * d[i + 1, s]
yield u * (c[i, 0] * c[i, 1] * d[i, 1] * d[i, 0])
return hamil, hamil.build_mpo(generate_terms, cutoff=cutoff).to_sparse()
def build_qc(filename, pg='d2h', cutoff=1E-9):
fcidump = FCIDUMP(pg=pg).read(fd)
hamil = Hamiltonian(fcidump, flat=False)
def generate_terms(n_sites, c, d):
for i in range(0, n_sites):
for j in range(0, n_sites):
for s in [0, 1]:
t = fcidump.t(s, i, j)
if abs(t) > cutoff:
yield t * c[i, s] * d[j, s]
for i in range(0, n_sites):
for j in range(0, n_sites):
for k in range(0, n_sites):
for l in range(0, n_sites):
for sij in [0, 1]:
for skl in [0, 1]:
v = fcidump.v(sij, skl, i, j, k, l)
if abs(v) > cutoff:
yield (0.5 * v) * (c[i, sij] * c[k, skl] * d[l, skl] * d[j, sij])
return hamil, hamil.build_mpo(generate_terms, cutoff=cutoff, const=hamil.fcidump.const_e).to_sparse()
Then the MPO can be built by:
hamil, mpo = build_hubbard(n=4)
or
fd = 'data/H8.STO6G.R1.8.FCIDUMP'
hamil, mpo = build_qc(fd, cutoff=1E-12)
From Hamiltonian expression (fast)
If C++ optimized code and numba are available, when there are very large number of terms in Hamiltonian,
the MPO building process can be accelerated:
First, we can use numba optimized functions to set the Hamiltonian terms:
import numpy as np
import numba as nb
flat = True
SPIN, SITE, OP = 1, 2, 16384
@nb.njit(nb.types.Tuple((nb.float64[:], nb.int32[:, :]))(nb.int32, nb.float64, nb.float64))
def generate_hubbard_terms(n_sites, u, t):
OP_C, OP_D = 0 * OP, 1 * OP
h_values = []
h_terms = []
for i in range(0, n_sites):
for s in [0, 1]:
if i - 1 >= 0:
h_values.append(t)
h_terms.append([OP_C + i * SITE + s * SPIN, OP_D + (i - 1) * SITE + s * SPIN, -1, -1])
if i + 1 < n_sites:
h_values.append(t)
h_terms.append([OP_C + i * SITE + s * SPIN, OP_D + (i + 1) * SITE + s * SPIN, -1, -1])
h_values.append(0.5 * u)
h_terms.append([OP_C + i * SITE + s * SPIN, OP_C + i * SITE + (1 - s) * SPIN,
OP_D + i * SITE + (1 - s) * SPIN, OP_D + i * SITE + s * SPIN])
return np.array(h_values, dtype=np.float64), np.array(h_terms, dtype=np.int32)
@nb.njit(nb.types.Tuple((nb.float64[:], nb.int32[:, :]))
(nb.int32, nb.float64[:, :], nb.float64[:, :, :, :], nb.float64))
def generate_qc_terms(n_sites, h1e, g2e, cutoff=1E-9):
OP_C, OP_D = 0 * OP, 1 * OP
h_values = []
h_terms = []
for i in range(0, n_sites):
for j in range(0, n_sites):
t = h1e[i, j]
if abs(t) > cutoff:
for s in [0, 1]:
h_values.append(t)
h_terms.append([OP_C + i * SITE + s * SPIN,
OP_D + j * SITE + s * SPIN, -1, -1])
for i in range(0, n_sites):
for j in range(0, n_sites):
for k in range(0, n_sites):
for l in range(0, n_sites):
v = g2e[i, j, k, l]
if abs(v) > cutoff:
for sij in [0, 1]:
for skl in [0, 1]:
h_values.append(0.5 * v)
h_terms.append([OP_C + i * SITE + sij * SPIN,
OP_C + k * SITE + skl * SPIN,
OP_D + l * SITE + skl * SPIN,
OP_D + j * SITE + sij * SPIN])
return np.array(h_values, dtype=np.float64), np.array(h_terms, dtype=np.int32)
def build_hubbard(u=2, t=1, n=8, cutoff=1E-9):
fcidump = FCIDUMP(pg='c1', n_sites=n, n_elec=n,
twos=0, ipg=0, orb_sym=[0] * n)
hamil = Hamiltonian(fcidump, flat=flat)
terms = generate_hubbard_terms(n, u, t)
return hamil, hamil.build_mpo(terms, cutoff=cutoff).to_sparse()
def build_qc(filename, pg='d2h', cutoff=1E-9, max_bond_dim=-1):
fcidump = FCIDUMP(pg=pg).read(filename)
hamil = Hamiltonian(fcidump, flat=flat)
terms = generate_qc_terms(
fcidump.n_sites, fcidump.h1e, fcidump.g2e, cutoff)
return hamil, hamil.build_mpo(terms, cutoff=cutoff, max_bond_dim=max_bond_dim, const=hamil.fcidump.const_e).to_sparse()
Then the MPO can be built by:
hamil, mpo = build_hubbard(n=16, cutoff=cutoff)
or
fd = 'data/H8.STO6G.R1.8.FCIDUMP'
hamil, mpo = build_qc(fd, cutoff=cutoff, max_bond_dim=-1)
From pyscf
The FCIDUMP can also be initialized using integral arrays, such as those obtained from pyscf.
Here is an example for H10 (STO6G).
Note that running pyblock3 and pyscf in the same python script with openMP activated may cause some conflicts in parallel MKL library, in some cases. One need to check number of threads used by pyblock3 during DMRG, to make sure that number of openMP threads is correct.
from pyscf import gto, scf, lo, symm, ao2mo
from pyblock3.fcidump import PointGroup
# H chain
N = 10
BOHR = 0.52917721092 # Angstroms
mpg = 'c1' # point group: d2h or c1
R = 1.8 * BOHR
mol = gto.M(atom=[['H', (0, 0, i * R)] for i in range(N)],
basis='sto6g', verbose=0, symmetry=mpg)
pg = mol.symmetry.lower()
mf = scf.RHF(mol)
ener = mf.kernel()
print("SCF Energy = %20.15f" % ener)
if pg == 'd2h':
fcidump_sym = ["Ag", "B3u", "B2u", "B1g", "B1u", "B2g", "B3g", "Au"]
elif pg == 'c1':
fcidump_sym = ["A"]
mo_coeff = mf.mo_coeff
n_mo = mo_coeff.shape[1]
orb_sym_str = symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mo_coeff)
xorb_sym = np.array([fcidump_sym.index(i) + 1 for i in orb_sym_str])
h1e = mo_coeff.T @ mf.get_hcore() @ mo_coeff
g2e = ao2mo.restore(1, ao2mo.kernel(mol, mo_coeff), n_mo)
ecore = mol.energy_nuc()
na = nb = mol.nelectron // 2
orb_sym = [PointGroup[mpg][i] for i in xorb_sym]
fd = FCIDUMP(pg='c1', n_sites=n_mo, n_elec=na + nb, twos=na - nb, ipg=0, uhf=False,
h1e=h1e, g2e=g2e, orb_sym=orb_sym, const_e=ecore, mu=0)
hamil = Hamiltonian(fd, flat=True)
mpo = hamil.build_qc_mpo()