Qiskitでトロッター分解を用いた量子シミュレーション

Dojoの以下記事をqulacsからqiskitに焼き直してみます。
https://dojo.qulacs.org/ja/latest/notebooks/4.2_trotter_decomposition.html

スクリプト

pip install qiskit pylatexenc
#必要なモジュールのインポート
from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute, Aer
from qiskit.qasm import pi
from qiskit.tools.visualization import plot_histogram, circuit_drawer
import numpy as np
# ダイナミクスをシミュレーションする時間
t = 3.0
## 横磁場の強さ
h = 3.
# トロッター分解の分割数
M = 100
# 時間の刻み幅
delta = t/M

## 時間と磁化を記録するリスト
x = [i*delta for i in range(M+1)]
y = []

import matplotlib.pyplot as plt

nqubits = 6

from qiskit.aqua.utils import tensorproduct
from qiskit.quantum_info.analysis import average_data
backend = Aer.get_backend('qasm_simulator')

for p in range(M+1):
  circuit_trotter_Ising = QuantumCircuit(nqubits) 
  for k in range(p):
    for i in range(nqubits):
      circuit_trotter_Ising.cx(i,(i+1)%(nqubits))
      circuit_trotter_Ising.rz(2*delta,(i+1)%nqubits) ## RZ(a)=exp(-i*a/2*Z)に注意. rz(rotation_angle, qbit_label)
      circuit_trotter_Ising.cx(i,(i+1)%(nqubits))
      circuit_trotter_Ising.rx(2*delta*h,i) ## RX(a)=exp(-i*a/2*X)に注意
  circuit_trotter_Ising.measure_all()           
  results = execute(circuit_trotter_Ising, backend, shots=1024).result()
  counts = results.get_counts()
  I_mat = np.array([[1, 0], [0, 1]])
  Z_mat = np.array([[1, 0], [0, -1]])
  obs = tensorproduct(Z_mat,I_mat,I_mat,I_mat,I_mat,I_mat) + tensorproduct(I_mat,Z_mat,I_mat,I_mat,I_mat,I_mat) + tensorproduct(I_mat,I_mat,Z_mat,I_mat,I_mat,I_mat) + + tensorproduct(I_mat,I_mat,I_mat,Z_mat,I_mat,I_mat) + tensorproduct(I_mat,I_mat,I_mat,I_mat,Z_mat,I_mat) + tensorproduct(I_mat,I_mat,Z_mat,I_mat,I_mat,Z_mat)
  obs = obs / nqubits
  expectation = average_data(counts, obs)
  print(expectation)
  y.append(expectation )

print(y)
print(x)
plt.plot(x,y)

Dojoの結果（ガタツキがないのは、状態ベクトルから解析的に物理量を計算しているため。）

まぁ、だいたいあっていそうです。
なぜか t= 0.5 のあたりがずれているけど・・・？

20201227追記
typoが見つかりました。

  obs = tensorproduct(Z_mat,I_mat,I_mat,I_mat,I_mat,I_mat) + tensorproduct(I_mat,Z_mat,I_mat,I_mat,I_mat,I_mat) + tensorproduct(I_mat,I_mat,Z_mat,I_mat,I_mat,I_mat) + + tensorproduct(I_mat,I_mat,I_mat,Z_mat,I_mat,I_mat) + tensorproduct(I_mat,I_mat,I_mat,I_mat,Z_mat,I_mat) + tensorproduct(I_mat,I_mat,Z_mat,I_mat,I_mat,Z_mat)

++もおかしいですが、最後の項に$Z$が二個あります・・・

これを正すと、

一致しました。
ちなみに、オブザーバブルを作るにはqiskit.quantum_info.operatorsを使うことも出来ます。
https://qiskit.org/documentation/locale/ja_JP/tutorials/circuits_advanced/02_operators_overview.html
若干、可読性が上がるかもしれません。

from qiskit.quantum_info.operators import Operator, Pauli
pauli0 = Pauli(label='ZIIIII')
pauli1 = Pauli(label='IZIIII')
pauli2 = Pauli(label='IIZIII')
pauli3 = Pauli(label='IIIZII')
pauli4 = Pauli(label='IIIIZI')
pauli5 = Pauli(label='IIIIIZ')

obs_op = (1/6)*(Operator(pauli0)+Operator(pauli1)+Operator(pauli2)+Operator(pauli3)+Operator(pauli4)+Operator(pauli5))
obs = obs_op.data

ちなみに回路の深さは、

という絶望的な値です。

解説

基本的にはDojo (Qulacs)のコードを移植すればよいのですが、
- Rzゲート、Rxゲートの位相の定義がプラスマイナス逆
- Qiskitでは任意のオブザーバブルの期待値を計算する機能がない

のが注意点です。後者はqiskit.aquaのライブラリを使うか、あるいはこの記事のように
qiskit.quantum_info.analysis.average_data
を使うかになります。
average_data は、測定結果の頻度と行列形式のオブザーバブルを引数に取ります。

Z基底で測定しているのでZの線形和の期待値しか計算できないように思うのですが、
なぜかオブザーバブルにXが含まれていても値が返ってきます。謎だ。。

2020/12/26追記
この部分は別記事で解説します。結構重要なところです。

結論

qiskitでシュレディンガー方程式が解けた。