I have used pyblock3 to construct MPS (Matrix Product State) and MPO (Matrix Product Operator) and then found the ground state energy with DMRG (like in the code). I tried to save all the MPS information. Now I want to map this MPS I have into a quantum circuit and measure the density matrix from this quantum MPS to verify if it gives the correct ground state energy.
I'm confused about how to map my MPS into a quantum circuit. Would somebody please explain this to me?
Thank you.
import numpy as np
from scipy.linalg import sqrtm, polar
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
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector, Operator
from qiskit_nature.second_q.operators import FermionicOp
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit.primitives import Estimator
import matplotlib.pyplot as plt
# Part 1: Classical MPS calculation with DMRG
# ------------------------------------------
fd = 'H2O.STO3G.FCIDUMP'
hamil = Hamiltonian(FCIDUMP(pg='d2h').read(fd), flat=True)
mpo = hamil.build_qc_mpo()
mpo, _ = mpopress(cutoff=1E-9, norm_cutoff=1E-9)
print('MPO (compressed) = ', mpo.show_bond_dims())
# Construct MPS
bond_dim = 200
mps = hamil.build_mps(bond_dim)
# Canonicalize MPS
mps = mps.canonicalize(center=0)
mps /= mps.norm()
# DMRG optimization
dmrg = MPE(mps, mpo, mps).dmrg(bdims=[bond_dim], noises=[1E-6, 0],
dav_thrds=[1E-3], iprint=2, n_sweeps=10)
ener = dmrg.energies[-1]
print("Energy(Ground State) = %20.12f" % ener)
print('MPS energy = ', np.dot(mps, mpo @ mps))
print('MPS = ', mps.show_bond_dims())
print('MPS norm = ', mps.norm())
print('DMRG: ', dmrg)
np.save("h2o_energy.npy", ener)
print("MPS after(bond dim): ", mps.show_bond_dims())
print(mps[0])
# Calculate one-particle density matrix (1PDM)
pdm1 = np.zeros((hamil.n_sites, hamil.n_sites))
for i in range(hamil.n_sites):
diop = OpElement(OpNames.D, (i, 0), q_label=SZ(-1, -1, hamil.orb_sym[i]))
di = hamil.build_site_mpo(diop)
for j in range(hamil.n_sites):
djop = OpElement(OpNames.D, (j, 0), q_label=SZ(-1, -1, hamil.orb_sym[j]))
dj = hamil.build_site_mpo(djop)
# factor 2 due to alpha + beta spins
pdm1[i, j] = 2 * np.dot((di @ mps).conj(), dj @ mps)
print("1PDM calculated from classical MPS:")
print(pdm1)
print("MPS after(bond dim): ", mps.show_bond_dims())
np.save("h2o_pdm1.npy", pdm1)
# Save the complete MPS information
mps_data = {
'n_sites': hamil.n_sites,
'bond_dims': [int(dim) for dim in mps.show_bond_dims().split('|')],
'tensors': [t.data.copy() if hasattr(t, 'data') else t.copy() for t in mps.tensors],
'q_labels': [t.q_labels if hasattr(t, 'q_labels') else None for t in mps.tensors],
'energy': ener,
'pdm1': pdm1
}
np.save("h2o_mps_complete.npy", mps_data, allow_pickle=True)
mps_data = np.load("h2o_mps_complete.npy", allow_pickle=True).item()
n_sites = mps_data['n_sites']
tensors = mps_data['tensors']
bond_dims = mps_data['bond_dims']
q_labels = mps_data['q_labels']
pdm1 = mps_data['pdm1']
energy_classical = mps_data['energy']
and this is the result of code:
QC MPO site 0 / 7
QC MPO site 1 / 7
QC MPO site 2 / 7
QC MPO site 3 / 7
QC MPO site 4 / 7
QC MPO site 5 / 7
QC MPO site 6 / 7
MPO (compressed) = 1|16|46|92|60|54|16|1 ( +9.81948)
Sweep = 0 | Direction = forward | BondDim = 200 | Noise = 1.00E-06 | DavThrd = 1.00E-03
--> Site = 0- 1 .. Mmps = 4 Ndav = 12 E = -53.939250566276 DW = 0.00E+00 FLOPS = 4.21E+06 Tdav = 0.011 T = 0.038 MEM = 197 MB
--> Site = 1- 2 .. Mmps = 16 Ndav = 18 E = -54.198125160401 DW = 0.00E+00 FLOPS = 6.47E+07 Tdav = 0.029 T = 0.041 MEM = 197 MB
--> Site = 2- 3 .. Mmps = 64 Ndav = 16 E = -54.845592489463 DW = 0.00E+00 FLOPS = 1.40E+08 Tdav = 0.051 T = 0.076 MEM = 197 MB
--> Site = 3- 4 .. Mmps = 144 Ndav = 73 E = -74.513037376288 DW = 8.31E-14 FLOPS = 1.12E+08 Tdav = 0.252 T = 0.296 MEM = 198 MB
--> Site = 4- 5 .. Mmps = 200 Ndav = 1 E = -74.513037376289 DW = 2.20E-13 FLOPS = 1.17E+08 Tdav = 0.005 T = 0.169 MEM = 212 MB
--> Site = 5- 6 .. Mmps = 57 Ndav = 1 E = -74.513037376289 DW = 4.96E-15 FLOPS = 5.19E+07 Tdav = 0.002 T = 0.295 MEM = 218 MB
Time elapsed = 0.915 | E = -74.513037376288793 | DE = 0.00E+00 | MDW = 2.20E-13 | MEM = 218 MB
Time sweep = 0.915 | Time davidson = 0.350 | Time decomp = 0.176
Sweep = 1 | Direction = backward | BondDim = 200 | Noise = 0.00E+00 | DavThrd = 1.00E-07
<-- Site = 5- 6 .. Mmps = 4 Ndav = 33 E = -74.925851220409 DW = 0.00E+00 FLOPS = 4.78E+07 Tdav = 0.065 T = 0.135 MEM = 208 MB
<-- Site = 4- 5 .. Mmps = 16 Ndav = 28 E = -74.931916324842 DW = 0.00E+00 FLOPS = 2.01E+08 Tdav = 0.081 T = 0.142 MEM = 205 MB
<-- Site = 3- 4 .. Mmps = 19 Ndav = 14 E = -74.931918959699 DW = 5.80E-18 FLOPS = 1.46E+08 Tdav = 0.037 T = 0.057 MEM = 196 MB
<-- Site = 2- 3 .. Mmps = 29 Ndav = 1 E = -74.931918959699 DW = 1.42E-17 FLOPS = 3.34E+07 Tdav = 0.006 T = 0.018 MEM = 196 MB
<-- Site = 1- 2 .. Mmps = 16 Ndav = 1 E = -74.931918959699 DW = 2.42E-17 FLOPS = 2.09E+07 Tdav = 0.004 T = 0.017 MEM = 196 MB
<-- Site = 0- 1 .. Mmps = 4 Ndav = 1 E = -74.931918959699 DW = 5.81E-18 FLOPS = 2.33E+06 Tdav = 0.002 T = 0.012 MEM = 196 MB
Time elapsed = 1.295 | E = -74.931918959698706 | DE = 4.19E-01 | MDW = 2.42E-17 | MEM = 208 MB
Time sweep = 0.380 | Time davidson = 0.196 | Time decomp = 0.003
Sweep = 2 | Direction = forward | BondDim = 200 | Noise = 0.00E+00 | DavThrd = 1.00E-07
--> Site = 0- 1 .. Mmps = 4 Ndav = 1 E = -74.931918959699 DW = 0.00E+00 FLOPS = 2.50E+06 Tdav = 0.002 T = 0.005 MEM = 196 MB
--> Site = 1- 2 .. Mmps = 16 Ndav = 1 E = -74.931918959699 DW = 0.00E+00 FLOPS = 2.28E+07 Tdav = 0.004 T = 0.011 MEM = 196 MB
--> Site = 2- 3 .. Mmps = 29 Ndav = 1 E = -74.931918959699 DW = 6.27E-23 FLOPS = 3.21E+07 Tdav = 0.006 T = 0.017 MEM = 196 MB
--> Site = 3- 4 .. Mmps = 19 Ndav = 1 E = -74.931918959699 DW = 9.49E-19 FLOPS = 4.83E+07 Tdav = 0.006 T = 0.019 MEM = 196 MB
--> Site = 4- 5 .. Mmps = 16 Ndav = 1 E = -74.931918959699 DW = 7.09E-20 FLOPS = 1.48E+07 Tdav = 0.003 T = 0.016 MEM = 196 MB
--> Site = 5- 6 .. Mmps = 4 Ndav = 1 E = -74.931918959699 DW = 8.89E-18 FLOPS = 2.84E+06 Tdav = 0.001 T = 0.009 MEM = 196 MB
Time elapsed = 1.372 | E = -74.931918959698692 | DE = 1.42E-14 | MDW = 8.89E-18 | MEM = 196 MB
Time sweep = 0.076 | Time davidson = 0.022 | Time decomp = 0.003
Energy(Ground State) = -74.931918959699
MPS energy = -74.93191895969869
MPS = 1|4|16|29|19|16|4|1
MPS norm = 1.0
DMRG: DMRG Energy = -74.931918959698692
MPS after(bond dim): 1|4|16|29|19|16|4|1
0 (Q=) (< N=0 SZ=0 PG=0 >, < N=2 SZ=0 PG=0 >, < N=2 SZ=0 PG=0 >) (R=) array([[[1.]]])
1 (Q=) (< N=0 SZ=0 PG=0 >, < N=1 SZ=-1/2 PG=0 >, < N=1 SZ=-1/2 PG=0 >) (R=) array([[[1.]]])
2 (Q=) (< N=0 SZ=0 PG=0 >, < N=0 SZ=0 PG=0 >, < N=0 SZ=0 PG=0 >) (R=) array([[[1.]]])
3 (Q=) (< N=0 SZ=0 PG=0 >, < N=1 SZ=1/2 PG=0 >, < N=1 SZ=1/2 PG=0 >) (R=) array([[[1.]]])
1PDM calculated from classical MPS:
[[ 1.99999624e+00 -2.15948426e-05 -9.24074207e-12 -7.73644238e-05
-7.63509115e-18 8.84393710e-05 -2.06222956e-11]
[-2.15948426e-05 1.99338711e+00 2.59887358e-10 8.79200962e-03
8.14000783e-16 -7.01051772e-03 2.22685962e-09]
[-9.24074207e-12 2.59887358e-10 1.97391091e+00 8.28097926e-11
1.24864262e-15 -6.57973151e-11 -1.57130600e-02]
[-7.73644238e-05 8.79200962e-03 8.28097926e-11 1.98251380e+00
-1.26136347e-15 -1.04953257e-02 2.46189704e-09]
[-7.63509115e-18 8.14000783e-16 1.24864262e-15 -1.26136347e-15
1.99850304e+00 -2.49761104e-14 1.79029966e-14]
[ 8.84393710e-05 -7.01051772e-03 -6.57973151e-11 -1.04953257e-02
-2.49761104e-14 2.32196065e-02 2.46069604e-11]
[-2.06222956e-11 2.22685962e-09 -1.57130600e-02 2.46189704e-09
1.79029966e-14 2.46069604e-11 2.84692985e-02]]
MPS after(bond dim): 1|4|16|29|19|16|4|1