%run init.ipynb
from qiskit import *

def Hxxx():
    H = tp(tp(pauli(1),pauli(1)),id(2)) + tp(id(2),tp(pauli(1),pauli(1)))\
    + tp(tp(pauli(2),pauli(2)),id(2)) + tp(id(2),tp(pauli(2),pauli(2)))\
    + tp(tp(pauli(3),pauli(3)),id(2)) + tp(id(2),tp(pauli(3),pauli(3)))
    return H
H = Hxxx(); H

#eig

#eig = H.eigenvects(); eig[0], eig[0][0], eig[0][2][0]

def mat_exp(A):
    d = A.shape[0]; eig = A.eigenvects(); Aexp = zeros(d,d); ne = 0; j = 0; lk = 0
    while ne < d:
        mult = eig[j][1]; ne += mult
        for k in range(0,mult):
            Proj = proj(eig[j][2][k]); Aexp += exp(eig[j][0])*(Proj/trace(Proj))
        j += 1
    return Aexp
simplify(mat_exp(pauli(1)))

def UHxxx(t):
    H = Hxxx()
    return mat_exp(-1j*t*H)
U = UHxxx(0); U

psi0 = tp(tp(cb(2,1),cb(2,1)),cb(2,0)); #psi0
dt = 0.05; t = np.arange(0, 2*pi+dt, dt); d = len(t); F = np.zeros(d);  psit = psi0
for j in range(0, d):
    U = UHxxx(t[j]); psit = U*psi0
    ip = inner_product(psi0, psit); F[j] = ip*conjugate(ip)
    #print(t[j], F[j])    

plt.figure(figsize = (6,4), dpi = 100)
plt.plot(t, F); plt.xlabel(r'$t$'); plt.ylabel(r'$F$'); plt.grid(); plt.show()

t = symbols('t'); U = UHxxx(t); psit = U*psi0; simplify(psit)

tp(tp(cb(2,0),cb(2,0)), cb(2,0)).T, tp(tp(cb(2,0),cb(2,0)), cb(2,1)).T, tp(tp(cb(2,0),cb(2,1)), cb(2,0)).T,\
tp(tp(cb(2,1),cb(2,0)), cb(2,0)).T, tp(tp(cb(2,0),cb(2,1)), cb(2,1)).T, tp(tp(cb(2,1),cb(2,0)), cb(2,1)).T,\
tp(tp(cb(2,1),cb(2,1)), cb(2,0)).T, tp(tp(cb(2,1),cb(2,1)), cb(2,1)).T

1/3, 1/6

H12 = tp(tp(pauli(1),pauli(1)),id(2)) + tp(tp(pauli(2),pauli(2)),id(2)) + tp(tp(pauli(3),pauli(3)),id(2))
H23 = tp(id(2),tp(pauli(1),pauli(1))) + tp(id(2),tp(pauli(2),pauli(2))) + tp(id(2),tp(pauli(3),pauli(3)))
dH = comm(H12,H23); dH, dH.norm()

H1 = tp(tp(pauli(1),pauli(1)),id(2)) + tp(id(2),tp(pauli(1),pauli(1))) + tp(id(2),tp(pauli(3),pauli(3)))
H2 = tp(tp(pauli(2),pauli(2)),id(2)) + tp(id(2),tp(pauli(2),pauli(2))) + tp(tp(pauli(3),pauli(3)),id(2))
dH = comm(H1,H2); dH, dH.norm()

comm(tp(pauli(1),pauli(1)),tp(pauli(2),pauli(2)))

XX = tp(pauli(1), pauli(1)); YY = tp(pauli(2), pauli(2)); ZZ = tp(pauli(3), pauli(3)); XX, YY, ZZ

XX.eigenvects(), YY.eigenvects()#, ZZ.eigenvects()

def Hxxx12():
    H = tp(pauli(1),pauli(1)) + tp(pauli(2),pauli(2)) + tp(pauli(3),pauli(3))
    return tp(H,id(2))
def Hxxx23():
    H = tp(pauli(1),pauli(1)) + tp(pauli(2),pauli(2)) + tp(pauli(3),pauli(3))
    return tp(id(2),H)
def UHxxx12(t):
    H = Hxxx12()
    return mat_exp(-1j*t*H)
def UHxxx23(t):
    H = Hxxx23()
    return mat_exp(-1j*t*H)

dt = 0.05; t = np.arange(0, 2*pi+dt, dt); d = len(t)
F1 = np.zeros(d); F2 = np.zeros(d); F3 = np.zeros(d); F4 = np.zeros(d); F5 = np.zeros(d)
psit1 = psi0; psit2 = psi0; psit3 = psi0; psit4 = psi0; psit5 = psi0
for j in range(0, d):
    U12 = UHxxx12(t[j]); U23 = UHxxx23(t[j]); U1 = U12*U23; psit1 = U1*psi0
    ip = inner_product(psi0,psit1); F1[j] = ip*conjugate(ip)
    #U12 = UHxxx12(t[j]/2); U23 = UHxxx23(t[j]/2); U2 = U12*U23*U12*U23; psit2 = U2*psi0
    #ip = inner_product(psi0,psit2); F2[j] = ip*conjugate(ip)
    #U12 = UHxxx12(t[j]/3); U23 = UHxxx23(t[j]/3); U3 = U12*U23*U12*U23*U12*U23; psit3 = U3*psi0
    #ip = inner_product(psi0,psit3); F3[j] = ip*conjugate(ip)
    #U12 = UHxxx12(t[j]/4); U23 = UHxxx23(t[j]/4); U4 = U12*U23*U12*U23*U12*U23*U12*U23; psit4 = U4*psi0
    #ip = inner_product(psi0,psit4); F4[j] = ip*conjugate(ip)
    #U12 = UHxxx12(t[j]/5); U23 = UHxxx23(t[j]/5); U5 = U12*U23*U12*U23*U12*U23*U12*U23*U12*U23; psit5 = U5*psi0
    #ip = inner_product(psi0,psit5); F5[j] = ip*conjugate(ip)

plt.figure(figsize = (6,3), dpi = 100)
plt.plot(t, F, '^', label = r'$F$'); plt.plot(t, F1, '--', label = r'$F_{1}$')
#plt.plot(t, F2, '-.', label = r'$F_{2}$'); plt.plot(t, F3, ':', label = r'$F_{3}$')
#plt.plot(t, F4, '.', label = r'$F_{4}$'); plt.plot(t, F5, '*', label = r'$F_{5}$')
plt.xlabel(r'$t$'); plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.grid(); plt.show()
# esse cálculo "simbólico" demora muito. Melhor usar o numpy!

# matrizes de Pauli
I = np.array([[1,0],[0,1]]); #print(I)
X = np.array([[0,1],[1,0]]); #print(X)
Y = np.array([[0,-1j],[1j,0]]);# print(Y)
Z = np.array([[1,0],[0,-1]]); #print(Z)

def Hxxx_num():
    H = np.kron(np.kron(X,X),I) + np.kron(I,np.kron(X,X))\
    + np.kron(np.kron(Y,Y),I) + np.kron(I,np.kron(Y,Y))\
    + np.kron(np.kron(Z,Z),I) + np.kron(I,np.kron(Z,Z))
    return H
H = Hxxx_num(); print(H.real)

# para calcular a função de matriz (teve problemas de pequenas oscilações)
def proj_num(psi): 
    #return np.outer(psi, np.ndarray.conjugate(psi))
    d = len(psi); proj = np.zeros((d, d), dtype = complex)
    for j in range(0, d):
        for k in range(0, d):
            proj[j,k] = psi[j]*np.conj(psi[k])
    return proj
proj_num(np.array([[1],[1j]]))

import cmath
#from math import e
# para calcular a função de matriz (teve problemas de pequenas oscilações)
def mat_exp_num(A): 
    d = A.shape[0]; Aexp = np.zeros((d,d), dtype = complex); 
    eva, eve = np.linalg.eig(A);
    # eva, eve = scipy.linalg.eig(A)
    for j in range(0, d):
        Aexp += cmath.exp(eva[j])*proj_num(eve[:,j])
        #Aexp += e**(eva[j])*proj_num(eve[:,j])
    return Aexp
mat_exp_num(Y)

def UHxxx_num(t):
    H = Hxxx_num()
    #return mat_exp_num(-1j*t*H)
    return scipy.linalg.expm(-1j*t*H) # resolveu o problema das pequenas flutuações
U = UHxxx_num(0*math.pi); U.real

def ip_num(a, b):
    d = len(a); ip = 0
    for j in range(0, d):
        ip += np.conj(a[j])*b[j]
    return ip
a = np.array([[1],[1j]]); b = np.array([[1],[1j]]); ip_num(a, b)

ket0 = np.array([[1],[0]]); ket1 = np.array([[0],[1]]); #ket0, ket1
dt = 0.05; t = np.arange(0, math.pi+dt, dt); d = len(t); F = np.zeros(d);  psit = np.zeros(d, dtype = complex)
psi0_ = np.zeros(d, dtype = complex); psi0_ = np.kron(ket1, np.kron(ket1, ket0)); #psi0
for j in range(0, d):
    U = UHxxx_num(t[j]); psit = np.dot(U, psi0_)
    #ip = ip_num(psi0_, psit); F[j] = np.abs(ip)**2
    F[j] = qiskit.quantum_info.state_fidelity(psi0_, psit)
plt.figure(figsize = (6,4), dpi = 100)
plt.plot(t, F); plt.xlabel(r'$t$'); plt.ylabel(r'$F$'); plt.grid(); plt.show()

H2 = np.kron(X, X) + np.kron(Y, Y) + np.kron(Z, Z)
def UHxxx12_num(t):
    H = np.kron(H2, I)
    return scipy.linalg.expm(-1j*t*H)
def UHxxx23_num(t):
    H = np.kron(I, H2)
    return scipy.linalg.expm(-1j*t*H)

dt = 0.05; t = np.arange(0, math.pi+dt, dt); d = len(t)
#F = np.zeros(d); 
F1 = np.zeros(d); F2 = np.zeros(d); F3 = np.zeros(d); F4 = np.zeros(d); F5 = np.zeros(d)
F6 = np.zeros(d); F7 = np.zeros(d); F8 = np.zeros(d); F9 = np.zeros(d); F10 = np.zeros(d); F11 = np.zeros(d)
#psi0_ = np.zeros(d, dtype = complex); 
psi0_ = np.kron(ket1, np.kron(ket1, ket0))
for j in range(0, d):
    U12 = UHxxx12_num(t[j]); U23 = UHxxx23_num(t[j])
    B = np.dot(U12, U23); U1 = np.linalg.matrix_power(B, 1); psit1 = np.dot(U1, psi0_)
    #ip = ip_num(psi0_, psit1); F1[j] = np.abs(ip)**2
    F1[j] = qiskit.quantum_info.state_fidelity(psi0_, psit1)
    U12 = UHxxx12_num(t[j]/2); U23 = UHxxx23_num(t[j]/2)
    B = np.dot(U12, U23); U2 = np.linalg.matrix_power(B, 2); psit2 = np.dot(U2, psi0_)
    #ip = ip_num(psi0_, psit2); F2[j] = np.abs(ip)**2
    F2[j] = qiskit.quantum_info.state_fidelity(psi0_, psit2)
    U12 = UHxxx12_num(t[j]/3); U23 = UHxxx23_num(t[j]/3)
    B = np.dot(U12, U23); U3 = np.linalg.matrix_power(B, 3); psit3 = np.dot(U3, psi0_)
    #ip = ip_num(psi0_, psit3); F3[j] = np.abs(ip)**2
    F3[j] = qiskit.quantum_info.state_fidelity(psi0_, psit3)
    U12 = UHxxx12_num(t[j]/4); U23 = UHxxx23_num(t[j]/4)
    B = np.dot(U12, U23); U4 = np.linalg.matrix_power(B, 4); psit4 = np.dot(U4, psi0_)
    #ip = ip_num(psi0_, psit4); F4[j] = np.abs(ip)**2
    F4[j] = qiskit.quantum_info.state_fidelity(psi0_, psit4)
    U12 = UHxxx12_num(t[j]/5); U23 = UHxxx23_num(t[j]/5)
    B = np.dot(U12, U23); U5 = np.linalg.matrix_power(B, 5); psit5 = np.dot(U5, psi0_)
    #ip = ip_num(psi0_, psit5); F5[j] = np.abs(ip)**2
    F5[j] = qiskit.quantum_info.state_fidelity(psi0_, psit5)
    U12 = UHxxx12_num(t[j]/6); U23 = UHxxx23_num(t[j]/6)
    B = np.dot(U12, U23); U6 = np.linalg.matrix_power(B, 6); psit6 = np.dot(U6, psi0_)
    #ip = ip_num(psi0_, psit6); F6[j] = np.abs(ip)**2
    F6[j] = qiskit.quantum_info.state_fidelity(psi0_, psit6)
    U12 = UHxxx12_num(t[j]/7); U23 = UHxxx23_num(t[j]/7)
    B = np.dot(U12, U23); U7 = np.linalg.matrix_power(B, 7); psit7 = np.dot(U7, psi0_)
    #ip = ip_num(psi0_, psit7); F7[j] = np.abs(ip)**2
    F7[j] = qiskit.quantum_info.state_fidelity(psi0_, psit7)
    U12 = UHxxx12_num(t[j]/8); U23 = UHxxx23_num(t[j]/8)
    B = np.dot(U12, U23); U8 = np.linalg.matrix_power(B, 8); psit8 = np.dot(U8, psi0_)
    #ip = ip_num(psi0_, psit8); F8[j] = np.abs(ip)**2
    F8[j] = qiskit.quantum_info.state_fidelity(psi0_, psit8)
    U12 = UHxxx12_num(t[j]/9); U23 = UHxxx23_num(t[j]/9)
    B = np.dot(U12, U23); U9 = np.linalg.matrix_power(B, 9); psit9 = np.dot(U9, psi0_)
    #ip = ip_num(psi0_, psit9); F9[j] = np.abs(ip)**2
    F9[j] = qiskit.quantum_info.state_fidelity(psi0_, psit9)
    U12 = UHxxx12_num(t[j]/10); U23 = UHxxx23_num(t[j]/10)
    B = np.dot(U12, U23); U10 = np.linalg.matrix_power(B, 10); psit10 = np.dot(U10, psi0_)
    #ip = ip_num(psi0_, psit10); F10[j] = np.abs(ip)**2
    F10[j] = qiskit.quantum_info.state_fidelity(psi0_, psit10)
    U12 = UHxxx12_num(t[j]/11); U23 = UHxxx23_num(t[j]/11)
    B = np.dot(U12, U23); U11 = np.linalg.matrix_power(B, 11); psit11 = np.dot(U11, psi0_)
    #ip = ip_num(psi0_, psit11); F11[j] = np.abs(ip)**2
    F11[j] = qiskit.quantum_info.state_fidelity(psi0_, psit11)
    if j == d-1:
        print('F5 = ', F5[j], ', F6 = ', F6[j])

plt.figure(figsize = (8,5), dpi = 100)
plt.plot(t, F, '^', label = r'$F$'); 
plt.plot(t, F1, '-.', label = r'$F_{1}$')
plt.plot(t, F2, '', label = r'$F_{2}$'); plt.plot(t, F3, '-.', label = r'$F_{3}$')
plt.plot(t, F4, '', label = r'$F_{4}$'); plt.plot(t, F5, '-.', label = r'$F_{5}$')
plt.plot(t, F6, '', label = r'$F_{6}$'); plt.plot(t, F7, '-.', label = r'$F_{7}$')
plt.plot(t, F8, '', label = r'$F_{8}$'); plt.plot(t, F9, '-.', label = r'$F_{9}$')
plt.plot(t, F10, '*', label = r'$F_{10}$'); plt.plot(t, F11, 'o', label = r'$F_{11}$')
plt.xlabel(r'$t$'); plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.grid(); plt.show()

def qc_ezz(t):
    qc = QuantumCircuit(2, name = 'e^(-itZZ)')
    qc.cx(0, 1); qc.rz(2*t, 1); qc.cx(0, 1) 
    return qc
qc_ezz_ = qc_ezz(math.pi); qc_ezz_.draw(output='mpl')

def qc_exx(t):
    qc = QuantumCircuit(2, name = 'e^(-itXX)')
    qc.h([0,1]); qc.cx(0, 1); qc.rz(2*t, 1); qc.cx(0, 1); qc.h([0,1])
    #qc.h([0,1]); qc_ezz_ = qc_ezz(t); qc.append(qc_ezz_, [0, 1]); qc.h([0,1])
    return qc
qc_exx_ = qc_exx(math.pi); qc_exx_.draw(output='mpl')

#qc_exx_.decompose().draw(output='mpl')

def qc_eyy(t):
    qc = QuantumCircuit(2, name = 'e^(-itYY)')
    qc.sdg([0,1]); qc.h([0,1]); qc.cx(0, 1); qc.rz(2*t, 1); qc.cx(0, 1); qc.h([0,1]); qc.s([0,1])
    #qc.sdg([0,1]); qc_exx_ = qc_exx(t); qc.append(qc_exx_, [0, 1]); qc.s([0,1])
    return qc
qc_eyy_ = qc_eyy(math.pi); qc_eyy_.draw(output='mpl')

#qc_eyy_.decompose().draw(output='mpl')

def qc_Bj(t):
    qc = QuantumCircuit(3, name = 'B_j')
    qc_ezz_ = qc_ezz(t); qc_eyy_ = qc_eyy(t); qc_exx_ = qc_exx(t)
    qc.append(qc_ezz_, [1, 2]); qc.append(qc_eyy_, [1, 2]); qc.append(qc_exx_, [1, 2])
    qc.append(qc_ezz_, [0, 1]); qc.append(qc_eyy_, [0, 1]); qc.append(qc_exx_, [0, 1])
    return qc
qc_Bj_ = qc_Bj(math.pi/2); qc_Bj_.draw(output='mpl')

qc_Bj_.decompose().draw(output='mpl')

def F_psi_rho_num(psi, rho):
    d = rho.shape[0]; F = 0
    for j in range(0, d):
        for k in range(0, d):
            F += np.conj(psi[j])*rho[j,k]*psi[k]
    return F.real
psi = [1/math.sqrt(2),1j/math.sqrt(2)]; rho = 0.5*np.array([[1,-1j],[1j,1]])
F = F_psi_rho_num(psi, rho); print(F) # ok!

nshots = 8192
qiskit.IBMQ.load_account()
#provider = IBMQ.get_provider(hub='ibm-q-community', group='ibmquantumawards', project='open-science-22')
#device = provider.get_backend('ibmq_jakarta')
provider = qiskit.IBMQ.get_provider(hub = 'ibm-q-research-2', group = 'federal-uni-sant-1', project = 'main')
device = provider.get_backend('ibmq_bogota')
simulator = Aer.get_backend('qasm_simulator')
from qiskit.tools.monitor import job_monitor
from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
from qiskit.tools.monitor import backend_overview, backend_monitor

backend_overview()

def qc_psi0():
    qc = QuantumCircuit(3, name = 'psi0')
    qc.x([0,1])
    return qc
qc_psi0_ = qc_psi0(); qc_psi0_.draw(output = 'mpl')

psi0_ = np.kron(ket1, np.kron(ket1, ket0)); 
psi0__ = np.kron(ket0, np.kron(ket1, ket1)) # para calcular a fidelidade com o estado da tomografia
psi0_.T, psi0__.T

nshots = 500
dt = 0.2; t = np.arange(0, math.pi+dt, dt); d = len(t)
F_teo = np.zeros(d); F_sim = np.zeros(d); F_exp = np.zeros(d)
for j in range(0, d):fidelity
    # teórico
    U12 = UHxxx12_num(t[j]/2); U23 = UHxxx23_num(t[j]/2)
    B = np.dot(U12, U23); U = np.linalg.matrix_power(B, 2); psit = np.dot(U, psi0_)
    F_teo[j] = qiskit.quantum_info.state_fidelity(psi0_, psit)
    # circuito quântico
    qc = QuantumCircuit(3)
    qc_psi0_ = qc_psi0(); qc.append(qc_psi0_, [0, 1, 2])
    qc_Bj_ = qc_Bj(t[j]/2); qc.append(qc_Bj_, [0, 1, 2]); qc.append(qc_Bj_, [0, 1, 2]) # 2 passos de Trotter
    qstc = state_tomography_circuits(qc, [0, 1, 2])
    # simulação
    job = qiskit.execute(qstc, backend = simulator, shots = nshots)
    qstf = StateTomographyFitter(job.result(), qstc)
    rho = qstf.fit(method = 'lstsq')
    F_sim[j] = qiskit.quantum_info.state_fidelity(psi0__, rho)
    # experimento
    jobE = qiskit.execute(qstc, backend = device, shots = nshots)
    print(jobE.job_id()); job_monitor(jobE)
    qstfE = StateTomographyFitter(jobE.result(), qstc)
    rhoE = qstfE.fit(method = 'lstsq')
    F_exp[j] = qiskit.quantum_info.state_fidelity(psi0__, rhoE)
qc.draw(output = 'mpl')

plt.figure(figsize = (8,5), dpi = 100)
plt.plot(t, F_teo, '-.', label = r'$F_{teo}$'); plt.plot(t, F_sim, '-.', label = r'$F_{sim}$')
plt.plot(t, F_exp, '-.', label = r'$F_{exp}$')
plt.xlabel(r'$t$'); plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')
plt.grid(); plt.show()

device = provider.get_backend('ibmq_jakarta')
jakarta_info = backend_monitor(device); print(jakarta_info)

nshots = 500
t = math.pi
for j in range(0, 7): # muda o No. de passos de Trotter
    # teórico
    U12 = UHxxx12_num(t/(j+1)); U23 = UHxxx23_num(t/(j+1))
    B = np.dot(U12, U23); U = np.linalg.matrix_power(B, j+1); psit = np.dot(U, psi0_)
    F_teo = qiskit.quantum_info.state_fidelity(psi0_, psit)
    # circuit q
    qc = QuantumCircuit(3)
    qc_psi0_ = qc_psi0(); qc.append(qc_psi0_, [0, 1, 2])
    qc_Bj_ = qc_Bj(t/(j+1))
    for k in range(0, j+1):
        qc.append(qc_Bj_, [0, 1, 2])
    qstc = state_tomography_circuits(qc, [0, 1, 2])
    # simulação
    job = qiskit.execute(qstc, backend = simulator, shots = nshots)
    qstf = StateTomographyFitter(job.result(), qstc)
    rho = qstf.fit(method = 'lstsq')
    F_sim = qiskit.quantum_info.state_fidelity(psi0__, rho)
    # experimento
    jobE = qiskit.execute(qstc, backend = device, shots = nshots)
    print(jobE.job_id()); job_monitor(jobE)
    qstfE = StateTomographyFitter(jobE.result(), qstc)
    rhoE = qstfE.fit(method = 'lstsq')
    F_exp = qiskit.quantum_info.state_fidelity(psi0__, rhoE)
    print('No. passos=', j+1, ',F_teo=', F_teo, ',F_sim=', F_sim, ',F_exp=', F_exp)

qc.draw(output = 'mpl')

nshots = 500; 
Npt = 7 # Número de passos de Trotter 
ph = 2*math.pi # ângulo para a transformação de simetria
t = math.pi
for j in range(0, Npt): # muda o No. de passos de Trotter
    # teórico
    U12 = UHxxx12_num(t/(j+1)); U23 = UHxxx23_num(t/(j+1))
    B = np.dot(U12, U23); U = np.linalg.matrix_power(B, j+1); psit = np.dot(U, psi0_)
    F_teo = qiskit.quantum_info.state_fidelity(psi0_, psit)
    # circuit q
    qc = QuantumCircuit(3)
    qc_psi0_ = qc_psi0(); qc.append(qc_psi0_, [0, 1, 2]) # estado inicial 
    qc_Bj_ = qc_Bj(t/(j+1))
    for k in range(0, j+1):
        qc.rz(ph/(j+1), [0, 1, 2]) # transformação de simetria
        qc.append(qc_Bj_, [0, 1, 2])
    qc.rz(-ph, [0, 1, 2]) # revertendo a transformação de simetria
    if j == 1:
        qc.draw(output = 'mpl')
    qstc = state_tomography_circuits(qc, [0, 1, 2])
    # simulação
    job = qiskit.execute(qstc, backend = simulator, shots = nshots)
    qstf = StateTomographyFitter(job.result(), qstc)
    rho = qstf.fit(method = 'lstsq')
    F_sim = qiskit.quantum_info.state_fidelity(psi0__, rho)
    # experimento
    jobE = qiskit.execute(qstc, backend = device, shots = nshots)
    print(jobE.job_id()); job_monitor(jobE)
    qstfE = StateTomographyFitter(jobE.result(), qstc)
    rhoE = qstfE.fit(method = 'lstsq')
    F_exp = qiskit.quantum_info.state_fidelity(psi0__, rhoE)
    print('No. passos=', j+1, ',F_teo=', F_teo, ',F_sim=', F_sim, ',F_exp=', F_exp)

qc.draw(output = 'mpl')

nshots = 1000; 
Npt = 7 # Número de passos de Trotter 
ph = 2*math.pi # ângulo para a transformação de simetria RZ
th = math.pi # ângulo para a transformação de simetria RX
t = math.pi
for j in range(0, Npt): # muda o No. de passos de Trotter
    # teórico
    U12 = UHxxx12_num(t/(j+1)); U23 = UHxxx23_num(t/(j+1))
    B = np.dot(U12, U23); U = np.linalg.matrix_power(B, j+1); psit = np.dot(U, psi0_)
    F_teo = qiskit.quantum_info.state_fidelity(psi0_, psit)
    # circuit q
    qc = QuantumCircuit(3)
    qc_psi0_ = qc_psi0(); qc.append(qc_psi0_, [0, 1, 2]) # estado inicial 
    qc_Bj_ = qc_Bj(t/(j+1))
    for k in range(0, j+1):
        qc.rx(th/(j+1), [0, 1, 2]) # transformação de simetria RX
        qc.rz(ph/(j+1), [0, 1, 2]) # transformação de simetria RZ
        qc.append(qc_Bj_, [0, 1, 2])
    qc.rz(-ph, [0, 1, 2]) # revertendo a transformação de simetria RZ
    qc.rx(-th, [0, 1, 2]) # revertendo a transformação de simetria RX
    if j == 1:
        qc.draw(output = 'mpl')
    qstc = state_tomography_circuits(qc, [0, 1, 2])
    # simulação
    job = qiskit.execute(qstc, backend = simulator, shots = nshots)
    qstf = StateTomographyFitter(job.result(), qstc)
    rho = qstf.fit(method = 'lstsq')
    F_sim = qiskit.quantum_info.state_fidelity(psi0__, rho)
    # experimento
    jobE = qiskit.execute(qstc, backend = device, shots = nshots)
    print(jobE.job_id()); job_monitor(jobE)
    qstfE = StateTomographyFitter(jobE.result(), qstc)
    rhoE = qstfE.fit(method = 'lstsq')
    F_exp = qiskit.quantum_info.state_fidelity(psi0__, rhoE)
    print('No. passos=', j+1, ',F_teo=', F_teo, ',F_sim=', F_sim, ',F_exp=', F_exp)

qc.draw(output = 'mpl')