Why ExpMPOEvolution updates my MPDO incorrectly.

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

Why ExpMPOEvolution updates my MPDO incorrectly.

Post by Jirawat S »

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
        site.add_op(opname, op)
    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[1] \(\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_onsite_with_entry(1e2,0,'Sx_I')
  Model.add_onsite_with_entry(-1e2,0,'I_Sx')
  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[0]))  # two values, one for each site!
      Pm.append(psi_MPDO.expectation_value(Op_proj[1]))
      if S == 1.0:
          Pp.append(psi_MPDO.expectation_value(Op_proj[2]))
      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_onsite_with_entry(2e2,0,'Sx_Sx')
  Model.add_coupling(J, 0, 'Sz_I', 0, 'Sz_I', 1, plus_hc=True)
  Model.MPOModel_init()
MPDO_fig.png
MPDO_fig.png (70.19 KiB) Viewed 1535 times
Did I do something wrong here?

[1]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.

Post by Jirawat S »

Okay, I've found the bug. :D
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
MPDO_fig_sqrt.png (80.37 KiB) Viewed 1479 times
Post Reply