## Why ExpMPOEvolution updates my MPDO incorrectly.

How do I use this algorithm? What does that parameter do?
Jirawat S
Posts: 8
Joined: 01 Mar 2022, 14:12

### Why ExpMPOEvolution updates my MPDO incorrectly.

EDIT: Problem solved, see in comment.

Hi,

I am trying to simulate time evolution of matrix product density operator (MPDO) using ExpMPOEvolution. My system is given by MPDOSite which is modified from SpinSite to has local dimension d^2. Here I define charge to be trivial for simplicity. Operators in MPO are given by kronecker product of local operators.

Code: Select all

class MPDOSite(Site):
def __init__(self, S=0.5, conserve='None'):
self.S = S = float(S)
d = int(2*S + 1)
#Rho doubles Hibert space.
leg = npc.LegCharge.from_trivial(d**2)
self.conserve = conserve
#Save index of eigenstates in vectorized rho. Will need for creating projectors.
self.DiagIdx = {str(S):idx*(d+1) for idx,S in enumerate(np.arange(-S,S+1,1))}

Iden_d = np.eye(d) #Local identity operator, not the same as Id from Site.
Sz_diag = -S + np.arange(d)
Sz = np.diag(Sz_diag)
Sp = np.zeros([d, d])
for n in np.arange(d - 1):
# Sp |m> =sqrt( S(S+1)-m(m+1)) |m+1>
m = n - S
Sp[n + 1, n] = np.sqrt(S * (S + 1) - m * (m + 1))
Sm = np.transpose(Sp)
# Sp = Sx + i Sy, Sm = Sx - i Sy
Sx = (Sp + Sm) * 0.5
Sy = (Sm - Sp) * 0.5j

Sx_Sx =  np.kron(Sx,Sx)
Sy_Sy =  np.kron(Sy,Sy)
Sz_Sz =  np.kron(Sz,Sz)
Sz_I = np.kron(Sz,Iden_d)
I_Sz = np.kron(Iden_d,Sz.T)
Sx_I = np.kron(Sx,Iden_d)
I_Sx = np.kron(Iden_d,Sx.T)
Sy_I = np.kron(Sy,Iden_d)
I_Sy = np.kron(Iden_d,Sy.T)
ops = dict(Sz_I=Sz_I, Sx_I=Sx_I, Sy_I=Sy_I,
I_Sz=I_Sz, I_Sx=I_Sx, I_Sy=I_Sy, Sx_Sx=Sx_Sx, Sy_Sy=Sy_Sy, Sz_Sz=Sz_Sz)
Site.__init__(self, leg, **ops)
self.state_labels = self.DiagIdx
self.state_labels.update({'down':self.DiagIdx[str(-S)]})
self.state_labels.update({'up':self.DiagIdx[str(S)]})

Calling init_sites to create MPDO

Code: Select all


def init_sites(S, MPDO=False):
"""
Define projectors for calculating expectation values.
"""
# Physical dimensions.
d = 2*S + 1
if MPDO == True:
site = MPDOSite(S=S, conserve='None')
d = d**2
else:
site = SpinSite(S=S, conserve='None')

if S == 1.0:
proj_op = zip(['-1.0', '0.0', '1.0'], ['Pm', 'P0', 'Pp'])
elif S == 0.5:
proj_op = zip(['-0.5', '0.5'], ['P0', 'Pm'])
else:
raise ValueError("Wrong S!!!")

for label, opname in proj_op:
diag = np.zeros(int(d)) #  diagonal
diag[site.state_index(label)] = 1.
op = np.diag(diag)  # 3x3 numpy array
return site

#Define MPDO
L = 2
S = 0.5
site = init_sites(S=S, MPDO=True)
p_leg = site.leg
chinfo = p_leg.chinfo
lat = Chain(L, site, bc_MPS='finite')
state = ["down", "up"]*(L//2) + ["down"]*(L%2)
psi_MPDO = MPS.from_product_state(lat.mps_sites(), state, lat.bc_MPS)


Now I implemented time evolution by define MPO from unitary terms in Hamiltonian $$\mathcal{L} = -\frac{i}{\hbar} H\otimes I + \frac{i}{\hbar} I\otimes H^{T}$$, where $$H = S_{x}$$.
So basically what I did is

Code: Select all

  #Define extended MPO.
J = 1e-4
Model = TenPy_MPO(lat)
Model.add_coupling(J, 0, 'Sz_I', 0, 'Sz_I', 1, plus_hc=True)
Model.MPOModel_init()

*coupling term is there to avoid dimensional issues.
And Called ExpMPOEvolution to update the state.

Code: Select all

  #Measure observables.
times = []
m_x, m_y, m_z = [], [], [] #expentation values
Pg, Pm, Pp = [], [], [] #Population in state |0>np.array([]), np.array([]), np.array([])
Op_proj = ['P0', 'Pm', 'Pp'] if S==1.0 else ['P0', 'Pm'] if S==0.5 else 0
legend =  ['P(|0>)', 'P(|-1>)', 'P(|+1>)'] if S == 1.0 else ['P(|-1/2>)', 'P(|+1/2>)'] if S == 0.5 else 0
assert legend != 0, "Wrong S!"

for i in tqdm(range(steps)):
#Getting probability in |0> state of spin 1.
Pg.append(psi_MPDO.expectation_value(Op_proj))  # two values, one for each site!
Pm.append(psi_MPDO.expectation_value(Op_proj))
if S == 1.0:
Pp.append(psi_MPDO.expectation_value(Op_proj))
engine = mpo_evolution.ExpMPOEvolution(psi_MPDO, Model, options=evo_params)
engine.run()
times.append(dt*i)
psi_MPDO = engine.psi

However, the results I got is not consistent to one from MPS (See in the plot). And I notice that the dynamics will be correct if I define MPO as $$\mathcal{L} = -\frac{i}{\hbar} H\otimes H^{T}$$

Code: Select all

  Model = TenPy_MPO(lat)
Model.add_coupling(J, 0, 'Sz_I', 0, 'Sz_I', 1, plus_hc=True)
Model.MPOModel_init() MPDO_fig.png (70.19 KiB) Viewed 640 times
Did I do something wrong here?

https://doi.org/10.48550/arXiv.1804.09796
Last edited by Jirawat S on 12 Sep 2022, 12:54, edited 1 time in total.
Jirawat S
Posts: 8
Joined: 01 Mar 2022, 14:12

### Re: Why ExpMPOEvolution updates my MPDO incorrectly.

Okay, I've found the bug. Since the function expectation_value(Op_proj) is defined to get a probability $$p = <O_{proj}>$$ from MPS, when $$O_{proj}$$ is a projection operator, by computing $$\frac{<\psi|O_{proj}|\psi>}{<\psi|\psi>}$$.
When I use it with MPDO, what the function really compute is $$p^2$$ (because it is $$\frac{<\rho|O_{proj}|\rho>}{<\rho|\rho>}$$).

After adding a square root to arrays of probability, I got the expected dynamics.
Hope this is useful to anyone who get the same problem. MPDO_fig_sqrt.png (80.37 KiB) Viewed 584 times