JuliaでOpenFermionを呼び出してVQE計算してみた

"Pythonのように書けてCのように動く"というのが売り文句に惹かれて、Juliaでいくつか簡単な量子計算を実装してみました。順次公開予定。まずはVQEです。
VQEの計算部分は湊さんの実装[1]と大きく変わりませんが、今回はハミルトニアンを天下り的に用意するのではなく、電子積分の計算、第二量子化、Jordan-Wigner変換をしてqubitで扱える形に変換するところまでをOpenFermion[2]を用いて行いました。イケてるからJuliaを使っているというのに、結局Pythonのライブラリを読み込んで使っていることになりますが…。使えるものは使っていきましょう。
対象は水素分子です。Quantum Native Dojo[3]を参考にさせていただきました。

PyCallでOpenFermionをインポート

JuliaでPythonのライブラリを使うために、PyCallを使います。@pyimportでopenfermionを読み込みます。ついでにクロネッカー積の演算子も定義します。

using PyCall, LinearAlgebra, Random, Optim, PyPlot
@pyimport openfermion
@pyimport openfermionpyscf
⊗(A,B) = kron(A,B)

ハミルトニアンの計算&変換

Quantum Native Dojo[3]の計算を使わせていただきました。

basis = "sto-3g"
multiplicity = 1
charge = 0
distance  = 0.9
geometry = [["H", [0,0,0]],["H", [0,0,distance]]]
description  = "H2"
molecule = openfermion.hamiltonians.MolecularData(geometry, basis, multiplicity, charge, description)
molecule = openfermionpyscf.run_pyscf(molecule,run_scf=1,run_fci=1)
N = molecule.n_qubits
n_electron = molecule.n_electrons
fermionic_hamiltonian = openfermion.transforms.get_fermion_operator(molecule.get_molecular_hamiltonian())
jw_hamiltonian_py = openfermion.transforms.jordan_wigner(fermionic_hamiltonian)

output:

PyObject (-0.25905412221337204+0j) [] +
(-0.04764292343981623+0j) [X0 X1 Y2 Y3] +
(0.04764292343981623+0j) [X0 Y1 Y2 X3] +
(0.04764292343981623+0j) [Y0 X1 X2 Y3] +
(-0.04764292343981623+0j) [Y0 Y1 X2 X3] +
(0.14907478844731503+0j) [Z0] +
(0.1611381637816487+0j) [Z0 Z1] +
(0.11162723403394147+0j) [Z0 Z2] +
(0.15927015747375767+0j) [Z0 Z3] +
(0.1490747884473151+0j) [Z1] +
(0.15927015747375767+0j) [Z1 Z2] +
(0.11162723403394147+0j) [Z1 Z3] +
(-0.16071249108067315+0j) [Z2] +
(0.167371259483041+0j) [Z2 Z3] +
(-0.16071249108067315+0j) [Z3]

Jordan-Wigner変換後のハミルトニアンが得られました。これで各qubitの|1>と|0>が各スピン軌道の占有・非占有に対応することになります。

VQEで使う関数の定義

(1ビット)回転行列やCNOTのような量子ゲートと、ansatz生成、期待値計算、コスト関数を定義します。ansatz回路は以下のものを使います(参考[4])。

#1量子ビット回転ゲート
function R_x(theta::Float64)
    return [cos(theta/2) -sin(theta/2)im; -sin(theta/2)im cos(theta/2)]
end
function R_z(theta::Float64)
    return [exp(-theta/2*im) 0; 0 exp(theta/2*im)]
end

#CNOT_i,j
function CNOT(control::Int64, target::Int64)
    CNOT_ide = 1
    CNOT_X = 1
    for i=1:N
        if i==control
            CNOT_ide = CNOT_ide ⊗ [1.0 0.0; 0.0 0.0]
            CNOT_X = CNOT_X ⊗ [0.0 0.0; 0.0 1.0]
        elseif i==target
            CNOT_ide = CNOT_ide ⊗ [1.0 0.0; 0.0 1.0]
            CNOT_X = CNOT_X ⊗ [0.0 1.0; 1.0 0.0]
        else
            CNOT_ide = CNOT_ide ⊗ [1.0 0.0; 0.0 1.0]
            CNOT_X = CNOT_X ⊗ [1.0 0.0; 0.0 1.0]
        end
    end
    return CNOT_ide + CNOT_X
end

#初期状態 |000…>を生成
function get_initial_state(nqubits::Int64)
    state = 1
    for i=1:nqubits
        state = state ⊗ [1.0, 0.0]
    end
    return state
end

#ansatz生成
function ansatz(theta_list::Array{Float64,3})
    state = get_initial_state(N)
    for i=1:size(theta_list,3), j=1:size(theta_list,2)
        gate = 1
        for k=1:size(theta_list,1)
            if j == 2
                gate = gate ⊗ R_x(theta_list[k,j,i])
            else
                gate = gate ⊗ R_z(theta_list[k,j,i])
            end
        end
        state = gate * state
        for j=1:N-1
            state = CNOT(j,j+1) * state
        end
    end
    return state
end

#期待値の計算
function get_expectation_value(state::Array{Complex{Float64},1}, observable::Array{Float64,2})
    return real(dot(state, observable*state))
end

#最適化計算用のコスト関数
function cost(theta_list::Array{Float64,3})
    return get_expectation_value(ansatz(theta_list), Hamiltonian)
end

計算実行

ハミルトニアンに使う$X_i, Y_i, Z_i;(\mathrm{dim}(X_i=2^N$))ゲートの定義して、OpenFermionで得たハミルトニアンを取り出して行列形式への変換を行います。attributesについて詳しくはOpenFemionのドキュメント[5]を確認してください(SymbolicOperatorクラス)。
最適化計算はOptim.jl[6]を使います。　

layer = 3
itr = 3

#Pauliゲートを定義(X[i]=X_i)
X = [Matrix(1.0I,1,1) for i in 1:N]
Y = [Matrix{Complex}(1.0I,1,1) for i in 1:N]
Z = [Matrix(1.0I,1,1) for i in 1:N]
for i=1:N, j=1:N
    if i==j
        X[i] = X[i] ⊗ [0.0 1.0; 1.0 0.0]
        Y[i] = Y[i] ⊗ [0.0 -1.0im; 1.0im 0.0]
        Z[i] = Z[i] ⊗ [1.0 0.0; 0.0 -1.0]
    else
        X[i] = X[i] ⊗ Matrix(1.0I,2,2)
        Y[i] = Y[i] ⊗ Matrix(1.0I,2,2)
        Z[i] = Z[i] ⊗ Matrix(1.0I,2,2)
    end
end

Hamiltonian = zeros(2^N,2^N)

#OpenFermion形式のハミルトニアンを行列に変換
for (pauli,coefficient) in (jw_hamiltonian_py[:terms])
    if length(pauli)==0
        Hamiltonian += real(coefficient) * Matrix(1.0I,2^N,2^N)
        continue
    end
    pauli_gates = 1
    for i=1:length(pauli)
        if pauli[i][2] == "X"
            pauli_gates = pauli_gates * X[(pauli[i][1]+1)]
        elseif pauli[i][2] == "Y"
            pauli_gates = pauli_gates * Y[(pauli[i][1]+1)]
        elseif pauli[i][2] == "Z"
            pauli_gates = pauli_gates * Z[(pauli[i][1]+1)]
        else
            println("something went wrong…")
        end
    end
    pauli_gates = real(coefficient) * pauli_gates
    Hamiltonian += pauli_gates
end
Hamiltonian = real(Hamiltonian)

#θの初期値をランダムに決定し、最適化計算を実行
rng = MersenneTwister(4838)
theta_list = rand(rng, N,layer,itr)*2pi
res = optimize(cost, theta_list, iterations=10000, store_trace=true)
plot(Optim.f_trace(res))
plt.xlabel("Iteration")
plt.ylabel("Energy (Hartree)")

output:

println(minimum(res))

output:

-1.1205602587271528

エネルギーが求められました。ハミルトニアンを直接対角化して、エネルギーが一致するか見てみましょう。juliaでは簡単にできますね。

eigvals(Hamiltonian)

output:

16-element Array{Float64,1}:
 -1.120560281299988
 -0.6828493923576944
 -0.6828493923576942
 -0.6828493923576939
 -0.5742460086733262
 -0.5742460086733259
 -0.5634367948093946
 -0.5634367948093945
 -0.30170600483916443
  0.032862358979865766
  0.03286235897986578
  0.0686039556493665
  0.0686039556493665
  0.1758813173718267
  0.5879746787999998
  0.634525489333432

VQEの結果とよく一致しています。

参考文献

[1]Juliaでお手製VQEを実行
[2]OpenFermion: The Electronic Structure Package for Quantum Computers, arXiv:1710.07629 [quant-ph]
[3]6-1. OpenFermionの使い方 (Quantum Native Dojo)
[4]量子コンピュータによる量子化学計算入門
[5] Optim.jl