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)]})
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)
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()
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
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()
[1]https://doi.org/10.48550/arXiv.1804.09796