statevector__simulator
Qiskitでトロッター分解を用いた量子シミュレーションの記事を以前に掲載しておりますが、
qasm_simulatorを用いたものでした。
https://qiita.com/notori48/items/dbac7bb0226bcfbfbf40
qasm_simulatorは実機と同様に状態作成後にサンプリングを行います。
対して、statevector__simulatorは状態作成をして、その状態を直接配列として取り出します。
（もちろん、実機ではこのようなことは通常出来ません。状態を推定する別のアルゴリズム、例えば量子状態トモグラフィなどが必要となります。この場合も推定誤差や雑音の影響が問題になります。）
今回はstatevector__simulatorを使います。
スクリプト
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
from qiskit.quantum_info.operators import Operator, Pauli
# Observable 作成
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
print(obs)
# ダイナミクスをシミュレーションする時間
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
from qiskit.quantum_info import Statevector
backend = Aer.get_backend('statevector_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()
result_statevec = Statevector(results.get_statevector())
expectation = result_statevec.expectation_value(obs_op)
print(expectation)
y.append(expectation)
imaginary part が出てきていますが、シミュレーターの計算誤差に起因するものです。
np.real() で排除できます。
ちゃんと出来ています。
https://dojo.qulacs.org/ja/latest/notebooks/4.2_trotter_decomposition.html
の結果と一致してますね
解説
statevector_simulatorを使うには、
backend = Aer.get_backend('statevector_simulator')
results = execute(circuit_trotter_Ising, backend, shots=1024).result()
でよいです。shotsは特に意味がない（1でもいい）かと思います。
なお、回路の最後にはmeasureを入れないことをおすすめします。
measureがあると、回路実行後に状態が確定（射影）してしまって、それが読み出されてしまいます。
シミュレーション後は、
results.get_statevector()
で状態ベクトルを読み出します。ただ、今回はここからオブザーバブルの期待値を計算したいです。
残念ながら直接はできないので、qiskit.quantum_infoにあるクラス Statevector を作ります。
そして、Statevector のメソッド（機能）としてある .expectation_value を使います。
https://qiskit.org/documentation/locale/ja_JP/stubs/qiskit.quantum_info.Statevector.html
result_statevec = Statevector(results.get_statevector())
expectation = result_statevec.expectation_value(obs_op)
引数はオブザーバブルと、オブザーバブルが作用する量子ビットのようです。
今回、オブザーバブルは、例えば $ZIIIII$ のように6量子ビットフルの表示をしてますので、量子ビットの指定は必要ないです。もし$Z$と書いた場合は、.expectation_value(obs_Z, 0) のように書くのではないでしょうか。
また、オブザーバブルとしては Operator のクラスをそのまま渡すか、Operator.data で配列に持ってきてから渡すかです。
このあたり、qiskitの関数によってまちまちです。
Operator(U) とすると、行列UをOperatorに変換できます。
結論
statevector_simulator でも計算できた。
確かにDojoと一致することがわかった。
なお、他のバックエンドとしてunitary_simulatorというのものあるようです。
これは回路全体を記述するユニタリ行列を取るsimulatorのようです。