MPS Import (Renormalizer)
[ ]:
!pip install block2==0.5.4rc6 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
!pip install pyscf==2.8.0 -qq --progress-bar off
!pip install renormalizer==0.0.11 -qq --progress-bar off
Preparing metadata (setup.py) ... done
Building wheel for renormalizer (setup.py) ... done
Introduction
In this tutorial we explain how to import an MPS from Renormalizer to block2. Renormalizer (https://github.com/shuaigroup/Renormalizer) is a Python tensor network package focused on electron-phonon quantum dynamics. This can be used for computing the overlap between a block2 MPS (in the SZ or SU2 mode) and an external MPS, which may be needed in quantum information or in Quantum Monte Carlo calculations (https://arxiv.org/abs/2405.05440).
The external MPS
We use the hydrogen chain as a test system. The required integrals for quantum chemistry Hamiltonian can be generated using pyscf.
[2]:
from pyscf import gto, scf, lo, tools
import numpy as np
from pyblock2._pyscf.ao2mo import integrals as itg
N, R = 6, 3.0
mf = scf.RHF(gto.M(atom=[['H', (0, 0, i * R)] for i in range(N)], basis='sto3g', verbose=3, symmetry='c1', unit='bohr'))
mf.mo_coeff = lo.orth.lowdin(mf.mol.intor('cint1e_ovlp_sph'))
ncas, n_elec, spin, ecore, h1e, h2e, orb_sym = itg.get_rhf_integrals(mf, g2e_symm=8)
We can find the ground state using the Renormalizer package with SZ symmetry, which is used as an external MPS.
[4]:
import tempfile
from renormalizer import Model, Mps, Mpo, optimize_mps
from renormalizer.model import h_qc
from renormalizer.utils import log
from logging import INFO
np.product = np.prod
logger = log.logging.getLogger("renormalizer")
logger.setLevel(INFO)
ftmp = tempfile.NamedTemporaryFile()
tools.fcidump.from_mo(mf.mol, ftmp.name, mf.mo_coeff)
h1e_spin, h2e_spin, nuc = h_qc.read_fcidump(ftmp.name, ncas)
basis, ham_terms = h_qc.qc_model(h1e_spin, h2e_spin)
model = Model(basis, ham_terms)
mpo = Mpo(model)
print(f"mpo_bond_dims:{mpo.bond_dims}")
M = 30
procedure = [[M, 0.4], [M, 0.2], [M, 0.1], [M, 0], [M, 0], [M, 0], [M, 0]]
mps = Mps.random(model, mf.mol.nelec, M, percent=1.0)
mps.optimize_config.procedure = procedure
mps.optimize_config.method = "2site"
energies, mps = optimize_mps(mps.copy(), mpo)
gs_e = min(energies)+nuc
print(f'lowest energy from Renormalizer MPS: {gs_e}')
2025-06-14 20:16:30,306[INFO] nuclear repulsion: 2.899999999999999
INFO:renormalizer.model.h_qc:nuclear repulsion: 2.899999999999999
2025-06-14 20:16:30,309[INFO] spin norbs: 12
INFO:renormalizer.model.h_qc:spin norbs: 12
2025-06-14 20:16:32,492[INFO] optimization method: 2site
INFO:renormalizer.mps.gs:optimization method: 2site
2025-06-14 20:16:32,495[INFO] e_rtol: 1e-06
INFO:renormalizer.mps.gs:e_rtol: 1e-06
2025-06-14 20:16:32,497[INFO] e_atol: 1e-08
INFO:renormalizer.mps.gs:e_atol: 1e-08
2025-06-14 20:16:32,499[INFO] procedure: [[30, 0.4], [30, 0.2], [30, 0.1], [30, 0], [30, 0], [30, 0], [30, 0]]
INFO:renormalizer.mps.gs:procedure: [[30, 0.4], [30, 0.2], [30, 0.1], [30, 0], [30, 0], [30, 0], [30, 0]]
mpo_bond_dims:[1, 4, 16, 37, 54, 68, 81, 68, 54, 37, 16, 4, 1]
2025-06-14 20:16:35,500[INFO] DMRG has converged!
INFO:renormalizer.mps.gs:DMRG has converged!
2025-06-14 20:16:35,575[INFO] mps current size: 48.4KiB, Matrix product bond dim:[1, 2, 4, 8, 16, 30, 30, 30, 16, 8, 4, 2, 1]
INFO:renormalizer.mps.gs:mps current size: 48.4KiB, Matrix product bond dim:[1, 2, 4, 8, 16, 30, 30, 30, 16, 8, 4, 2, 1]
lowest energy from Renormalizer MPS: -2.957646039005119
The following code is used to convert the external MPS from a spin-orbital basis to fermionic spatial orbtial basis with the left-canonicalized format.
[5]:
mps.ensure_left_canonical()
def spin_to_fermion_mps(mps_spin, qn):
nsites_spin = len(mps_spin)
assert nsites_spin % 2 == 0
nsites_fermion = nsites_spin // 2
qnl_fermion = []
mps_fermion = []
for i in range(nsites_fermion):
ml = mps_spin[2*i].shape[0]
mr = mps_spin[2*i+1].shape[-1]
merged_array = np.einsum("ipj, jql->ipql", mps_spin[2*i], mps_spin[2*i+1])
merged_array = np.transpose(merged_array, (0, 2, 1, 3))
merged_array = merged_array.reshape(ml, 4, mr)
mps_fermion.append(merged_array)
qnl_fermion.append(qn[2*i])
qnl_fermion.append(qn[-1])
return mps_fermion, qnl_fermion
qn = []
for i_mps in range(len(mps)):
iqn = [x.tolist() for x in mps.qn[i_mps]]
qn.append(iqn)
qn[0] = [[0, 0]]
qn.append([[mf.mol.nelec[0], mf.mol.nelec[1]]])
mps_fermion, qn_fermion = spin_to_fermion_mps(mps, qn)
mps_tensors = [mps_fermion[i_mps] for i_mps in range(len(mps_fermion))]
print("bond dimensions:", [ts.shape[0] for ts in mps_tensors])
print("quantum numbers:", qn_fermion)
bond dimensions: [1, 4, 16, 30, 16, 4]
quantum numbers: [[[0, 0]], [[1, 0], [0, 1], [1, 1], [0, 0]], [[0, 1], [0, 1], [1, 2], [1, 2], [2, 1], [2, 1], [0, 0], [1, 1], [1, 1], [1, 1], [1, 1], [2, 0], [0, 2], [2, 2], [1, 0], [1, 0]], [[0, 1], [1, 2], [1, 2], [1, 2], [1, 2], [2, 1], [2, 1], [2, 1], [2, 1], [3, 1], [3, 1], [1, 1], [1, 1], [1, 1], [1, 1], [0, 3], [2, 0], [2, 0], [3, 0], [2, 3], [0, 2], [0, 2], [2, 2], [2, 2], [2, 2], [2, 2], [1, 0], [3, 2], [1, 3], [1, 3]], [[1, 2], [1, 2], [2, 1], [2, 1], [3, 1], [1, 1], [2, 3], [2, 3], [3, 3], [2, 2], [2, 2], [2, 2], [2, 2], [3, 2], [3, 2], [1, 3]], [[2, 3], [3, 3], [2, 2], [3, 2]], [[3, 3]]]
MPS overlap
Here we show how to transform this external MPS to a block2 MPS and compute the overlap between this MPS and the ground state MPS optimized in block2 for the same system.
[6]:
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
from pyblock2.algebra.io import MPSTools
bond_dims = [20] * 4 + [30] * 4
noises = [1e-4] * 4 + [1e-5] * 4 + [0]
thrds = [1e-10] * 8
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
pybra_sz = MPSTools.from_renormalizer_dense_sz(mps_tensors, qn_fermion)
pybra_su2 = MPSTools.trans_sz_to_su2(pybra_sz, driver.basis, driver.target, target_twos=0)
The overlap computation can be done in either the SZ or SU2 mode in block2.
In the SZ mode:
[7]:
driver.symm_type = SymmetryTypes.SZ
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
zbra = MPSTools.to_block2(pybra_sz, driver.basis, tag='ZBRA')
mpo = driver.get_qc_mpo(h1e=h1e, g2e=h2e, ecore=ecore, iprint=1)
zket = driver.get_random_mps(tag="ZKET", bond_dim=20, nroots=1)
energies = driver.dmrg(mpo, zket, n_sweeps=20, bond_dims=bond_dims, noises=noises, thrds=thrds, iprint=1)
zket = driver.adjust_mps(zket, dot=1)[0]
impo = driver.get_identity_mpo()
print('Overlap (SZ):', driver.expectation(zbra, impo, zket))
integral symmetrize error = 0.0
integral cutoff error = 0.0
mpo terms = 2286
Build MPO | Nsites = 6 | Nterms = 2286 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 6 .. Mmpo = 26 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.001 T = 0.009
Site = 1 / 6 .. Mmpo = 66 DW = 0.00e+00 NNZ = 243 SPT = 0.8584 Tmvc = 0.001 T = 0.010
Site = 2 / 6 .. Mmpo = 110 DW = 0.00e+00 NNZ = 459 SPT = 0.9368 Tmvc = 0.001 T = 0.011
Site = 3 / 6 .. Mmpo = 66 DW = 0.00e+00 NNZ = 1147 SPT = 0.8420 Tmvc = 0.001 T = 0.012
Site = 4 / 6 .. Mmpo = 26 DW = 0.00e+00 NNZ = 243 SPT = 0.8584 Tmvc = 0.000 T = 0.008
Site = 5 / 6 .. Mmpo = 1 DW = 0.00e+00 NNZ = 26 SPT = 0.0000 Tmvc = 0.000 T = 0.006
Ttotal = 0.057 Tmvc-total = 0.004 MPO bond dimension = 110 MaxDW = 0.00e+00
NNZ = 2144 SIZE = 18004 SPT = 0.8809
Rank = 0 Ttotal = 0.120 MPO method = FastBipartite bond dimension = 110 NNZ = 2144 SIZE = 18004 SPT = 0.8809
Sweep = 0 | Direction = forward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.663 | E = -2.9576460853 | DW = 6.66567e-05
Sweep = 1 | Direction = backward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.724 | E = -2.9576460853 | DE = 1.04e-11 | DW = 5.92236e-05
Sweep = 2 | Direction = forward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.779 | E = -2.9576460853 | DE = -3.64e-11 | DW = 6.66526e-05
Sweep = 3 | Direction = backward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.837 | E = -2.9576460853 | DE = -2.38e-12 | DW = 5.92239e-05
Sweep = 4 | Direction = forward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.896 | E = -2.9576460853 | DE = 2.36e-12 | DW = 5.02379e-06
Sweep = 5 | Direction = backward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.959 | E = -2.9576460853 | DE = -3.69e-11 | DW = 2.73774e-06
Sweep = 6 | Direction = forward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.024 | E = -2.9576460853 | DE = -1.95e-11 | DW = 5.02112e-06
Sweep = 7 | Direction = backward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.084 | E = -2.9576460853 | DE = 1.93e-11 | DW = 2.73775e-06
Sweep = 8 | Direction = forward | Bond dimension = 30 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 1.144 | E = -2.9576460852 | DE = 1.66e-10 | DW = 1.91390e-06
Overlap (SZ): 0.9999999659573686
In the SU2 mode:
[8]:
driver.symm_type = SymmetryTypes.SU2
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
bra = MPSTools.to_block2(pybra_su2, driver.basis, tag='BRA')
mpo = driver.get_qc_mpo(h1e=h1e, g2e=h2e, ecore=ecore, iprint=1)
ket = driver.get_random_mps(tag="KET", bond_dim=20, nroots=1)
energies = driver.dmrg(mpo, ket, n_sweeps=20, bond_dims=bond_dims, noises=noises, thrds=thrds, iprint=1)
ket = driver.adjust_mps(ket, dot=1)[0]
impo = driver.get_identity_mpo()
print('Overlap (SU2):', driver.expectation(bra, impo, ket))
integral symmetrize error = 0.0
integral cutoff error = 0.0
mpo terms = 863
Build MPO | Nsites = 6 | Nterms = 863 | Algorithm = FastBIP | Cutoff = 1.00e-20
Site = 0 / 6 .. Mmpo = 13 DW = 0.00e+00 NNZ = 13 SPT = 0.0000 Tmvc = 0.000 T = 0.003
Site = 1 / 6 .. Mmpo = 34 DW = 0.00e+00 NNZ = 97 SPT = 0.7805 Tmvc = 0.000 T = 0.003
Site = 2 / 6 .. Mmpo = 56 DW = 0.00e+00 NNZ = 180 SPT = 0.9055 Tmvc = 0.001 T = 0.005
Site = 3 / 6 .. Mmpo = 34 DW = 0.00e+00 NNZ = 405 SPT = 0.7873 Tmvc = 0.000 T = 0.008
Site = 4 / 6 .. Mmpo = 14 DW = 0.00e+00 NNZ = 100 SPT = 0.7899 Tmvc = 0.000 T = 0.003
Site = 5 / 6 .. Mmpo = 1 DW = 0.00e+00 NNZ = 14 SPT = 0.0000 Tmvc = 0.000 T = 0.002
Ttotal = 0.024 Tmvc-total = 0.002 MPO bond dimension = 56 MaxDW = 0.00e+00
NNZ = 809 SIZE = 4753 SPT = 0.8298
Rank = 0 Ttotal = 0.038 MPO method = FastBipartite bond dimension = 56 NNZ = 809 SIZE = 4753 SPT = 0.8298
Sweep = 0 | Direction = forward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.046 | E = -2.9576460853 | DW = 1.81931e-06
Sweep = 1 | Direction = backward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.081 | E = -2.9576460854 | DE = -5.92e-11 | DW = 9.59991e-07
Sweep = 2 | Direction = forward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.116 | E = -2.9576460854 | DE = -3.03e-12 | DW = 1.81905e-06
Sweep = 3 | Direction = backward | Bond dimension = 20 | Noise = 1.00e-04 | Dav threshold = 1.00e-10
Time elapsed = 0.150 | E = -2.9576460854 | DE = 4.31e-12 | DW = 9.60002e-07
Sweep = 4 | Direction = forward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.597 | E = -2.9576460854 | DE = -4.28e-12 | DW = 4.15116e-09
Sweep = 5 | Direction = backward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.718 | E = -2.9576460854 | DE = -3.61e-12 | DW = 7.70782e-10
Sweep = 6 | Direction = forward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 0.865 | E = -2.9576460854 | DE = -4.10e-12 | DW = 4.13558e-09
Sweep = 7 | Direction = backward | Bond dimension = 30 | Noise = 1.00e-05 | Dav threshold = 1.00e-10
Time elapsed = 1.129 | E = -2.9576460854 | DE = 1.87e-12 | DW = 7.70814e-10
Sweep = 8 | Direction = forward | Bond dimension = 30 | Noise = 0.00e+00 | Dav threshold = 1.00e-09
Time elapsed = 1.398 | E = -2.9576460851 | DE = 2.61e-10 | DW = 6.52891e-18
Overlap (SU2): 0.9999990300129581