6-1. OpenFermionの使い方¶
この節では量子化学計算用のPythonライブラリである、OpenFermion [1] を用いて、相互作用する電子系のハミルトニアンを、量子コンピュータ上で扱いやすい形に変換する方法を紹介する。OpenFermion には量子化学計算のオープンソースライブラリである Psi4 および PySCF との接続が用意されており、これらのライブラリの詳細な使い方を理解しなくても、分子の構造を入力するだけで、量子化学計算において現れる電子系のハミルトニアンを得られるようになっている。ここでは PySCF を使用する。
[ ]:
## 各種ライブラリがインストールされていない場合は実行してください
## Google Colaboratory上で実行する場合'You must restart the runtime in order to use newly installed versions.'と出ますが無視してください。
## runtimeを再開するとクラッシュします。
!pip install qulacs pyscf openfermion openfermionpyscf
## Google Colaboratory / (Linux or Mac)のjupyter notebook 環境の場合にのみ実行してください。
## Qulacsのエラーが正常に出力されるようになります。
!pip3 install wurlitzer
%load_ext wurlitzer
[1]:
#必要なライブラリのインポート
#エラーが出る場合は openfermion を v.1.0.0 以上にしてみてください
import numpy as np
import matplotlib.pyplot as plt
from openfermion.chem import MolecularData
from openfermion.transforms import get_fermion_operator, jordan_wigner, bravyi_kitaev
from openfermion.linalg import get_sparse_operator
from openfermion.ops import FermionOperator
from openfermionpyscf import run_pyscf
from pyscf import fci
水素分子を計算してみる¶
openfermion では、分子を記述するデータを MolecularData というクラスに入力する。
[2]:
#define constants
basis = "sto-3g" #basis set
multiplicity = 1 #spin multiplicity
charge = 0 #total charge for the molecule
distance = 0.65
geometry = [("H",(0,0,0)),("H", (0,0,distance))] #xyz coordinates for atoms
description = str(distance) #description for the psi4 output file
molecule = MolecularData(geometry, basis, multiplicity, charge, description)
変数の説明¶
以下で上記のコード内で現れている変数の意味を説明する。
basis: 基底関数¶
分子軌道を表現するための基底関数を設定する。sto-3g, 6-31G などいろいろな基底関数系がある。
ここで使った sto-3g (Slater Type Orbital - 3 gaussian) は Slater type orbital を 3つのgaussianで近似した基底関数である。
Slater type orbital とは、水素原子の解に似せた軌道であり、動径方向の関数として
を使用し、角度方向は球面調和関数\(Y_{lm}(\theta,\phi)\)を使用するものである。sto-3g では、この動径方向の波動関数\(R_{nl}(r)\)を、3つのgaussianで近似した関数を用いる。
multiplicity: スピン多重度¶
電子はスピン1/2を持っているので、1つの電子が孤立して存在しているときスピン多重度は2である。しかし水素分子の場合、基底状態では電子はsingletを組み、全体ではスピン0になっていると考えられる。スピン0は1状態のみなので、この場合ではスピン多重度は1とする。
charge: 全電荷¶
全体の電荷を入力する。イオンを考える場合は + になったり − になったりする。
geometry: 原子核配置¶
原子種とその座標を x,y,z で指定する。
description¶
pyscf が計算した出力結果は openfermion のライブラリが保存されているディレクトリ内に保存される。そのファイルの名前を決めるための変数である。
PySCF による計算¶
上記で設定した MolecularData を関数 run_pyscf に投げて PySCFによる量子化学計算を行ってみよう。数秒で終わるはずである。
[3]:
molecule = run_pyscf(molecule,run_scf=1,run_fci=1)
HF & Full-CI energy¶
PySCF の計算によって求まった Hartree-Fock エネルギーと Full-CI エネルギー (=厳密な基底エネルギー) を見てみよう。(1 Hartree = 27.2116 eV)
[4]:
print("HF energy: {} (Hartree)".format(molecule.hf_energy))
print("FCI energy: {} (Hartree)".format(molecule.fci_energy))
HF energy: -1.1129965456691684 (Hartree)
FCI energy: -1.1299047843229135 (Hartree)
1 電子積分 \(h_{ij}\)・2電子積分 \(h_{ijkl}\)¶
1 電子積分や 2 電子積分といった量も MolecularData クラスに保存されている。
[5]:
print(molecule.one_body_integrals)
[[-1.30950987e+00 1.98461056e-17]
[ 8.14840166e-17 -4.10026381e-01]]
[6]:
print(molecule.two_body_integrals)
[[[[ 6.91904405e-01 -4.16333634e-17]
[-2.77555756e-17 1.76318452e-01]]
[[-2.77555756e-17 1.76318452e-01]
[ 6.79683914e-01 0.00000000e+00]]]
[[[-4.16333634e-17 6.79683914e-01]
[ 1.76318452e-01 8.32667268e-17]]
[[ 1.76318452e-01 8.32667268e-17]
[ 0.00000000e+00 7.14671111e-01]]]]
第二量子化形式のハミルトニアン¶
openfermionはこれらの積分値から第二量子化形式のハミルトニアン
を計算してくれる（第二量子化については、例えば参考文献[2]を参照）。 get_molecular_hamiltonian メソッドを呼ぶことでハミルトニアンが計算できる。
表示は (3,1)が \(c_3^\dagger\), (1,0)が \(c_1\) といった具合。
[7]:
print(molecule.get_molecular_hamiltonian())
() 0.8141187860307693
((0, 1), (0, 0)) -1.309509868464871
((1, 1), (1, 0)) -1.309509868464871
((2, 1), (2, 0)) -0.4100263808117837
((3, 1), (3, 0)) -0.4100263808117837
((0, 1), (0, 1), (0, 0), (0, 0)) 0.3459522026149022
((0, 1), (0, 1), (2, 0), (2, 0)) 0.0881592258051036
((0, 1), (1, 1), (1, 0), (0, 0)) 0.3459522026149022
((0, 1), (1, 1), (3, 0), (2, 0)) 0.0881592258051036
((0, 1), (2, 1), (0, 0), (2, 0)) 0.0881592258051036
((0, 1), (2, 1), (2, 0), (0, 0)) 0.3398419569652305
((0, 1), (3, 1), (1, 0), (2, 0)) 0.0881592258051036
((0, 1), (3, 1), (3, 0), (0, 0)) 0.3398419569652305
((1, 1), (0, 1), (0, 0), (1, 0)) 0.3459522026149022
((1, 1), (0, 1), (2, 0), (3, 0)) 0.0881592258051036
((1, 1), (1, 1), (1, 0), (1, 0)) 0.3459522026149022
((1, 1), (1, 1), (3, 0), (3, 0)) 0.0881592258051036
((1, 1), (2, 1), (0, 0), (3, 0)) 0.0881592258051036
((1, 1), (2, 1), (2, 0), (1, 0)) 0.3398419569652305
((1, 1), (3, 1), (1, 0), (3, 0)) 0.0881592258051036
((1, 1), (3, 1), (3, 0), (1, 0)) 0.3398419569652305
((2, 1), (0, 1), (0, 0), (2, 0)) 0.33984195696523034
((2, 1), (0, 1), (2, 0), (0, 0)) 0.0881592258051036
((2, 1), (1, 1), (1, 0), (2, 0)) 0.33984195696523034
((2, 1), (1, 1), (3, 0), (0, 0)) 0.0881592258051036
((2, 1), (2, 1), (0, 0), (0, 0)) 0.0881592258051036
((2, 1), (2, 1), (2, 0), (2, 0)) 0.35733555551906837
((2, 1), (3, 1), (1, 0), (0, 0)) 0.0881592258051036
((2, 1), (3, 1), (3, 0), (2, 0)) 0.35733555551906837
((3, 1), (0, 1), (0, 0), (3, 0)) 0.33984195696523034
((3, 1), (0, 1), (2, 0), (1, 0)) 0.0881592258051036
((3, 1), (1, 1), (1, 0), (3, 0)) 0.33984195696523034
((3, 1), (1, 1), (3, 0), (1, 0)) 0.0881592258051036
((3, 1), (2, 1), (0, 0), (1, 0)) 0.0881592258051036
((3, 1), (2, 1), (2, 0), (3, 0)) 0.35733555551906837
((3, 1), (3, 1), (1, 0), (1, 0)) 0.0881592258051036
((3, 1), (3, 1), (3, 0), (3, 0)) 0.35733555551906837
量子コンピュータの扱いやすい演算子に変換する¶
量子コンピュータ上で一番扱いやすいのは、 Pauli 演算子 \(I, X, Y, Z\) とそのテンソル積である。そこで、普通電子のハミルトニアンを量子コンピュータで扱うには、第二量子化形式のハミルトニアン
を、
の形に変換する。様々な変換方法が提案されているが、ここでは Jordan-Wigner 変換と呼ばれている一番簡単なものを使う。Jordan-Wigner 変換では、分子軌道 \(i\) を \(i\) 番目の qubit に対応させ、その分子軌道を電子が占有しているという状況を \(|1\rangle\), そうでないときには \(|0\rangle\) で表すという約束をする。
このような約束の下で、fermion の生成消滅演算子の反交換関係
を満たすようにパウリ演算子を構成すると、
という対応関係を得る。
Jordan-Wigner 変換以外の変換方式については、[2][3] などを参照されたい。
openfermion では Jordan-Wigner 変換が実装されている。jordan_wigner 関数に FermionOperator を渡すことで、その演算子の Jordan-Wigner 変換に対応する QubitOperator を返してくれる。以下では、上で作り出した水素分子の MolecularData から FermionOperator を作り出し、Jordan-Wigner 変換することで水素分子のハミルトニアンを量子コンピュータの扱いやすい形に変換している。
[8]:
jw_hamiltonian = jordan_wigner(get_fermion_operator(molecule.get_molecular_hamiltonian()))
print(jw_hamiltonian)
(0.0377511039464572+0j) [] +
(-0.0440796129025518+0j) [X0 X1 Y2 Y3] +
(0.0440796129025518+0j) [X0 Y1 Y2 X3] +
(0.0440796129025518+0j) [Y0 X1 X2 Y3] +
(-0.0440796129025518+0j) [Y0 Y1 X2 X3] +
(0.18601648886230573+0j) [Z0] +
(0.1729761013074511+0j) [Z0 Z1] +
(0.12584136558006342+0j) [Z0 Z2] +
(0.16992097848261523+0j) [Z0 Z3] +
(0.18601648886230576+0j) [Z1] +
(0.16992097848261523+0j) [Z1 Z2] +
(0.12584136558006342+0j) [Z1 Z3] +
(-0.26941693141632095+0j) [Z2] +
(0.17866777775953419+0j) [Z2 Z3] +
(-0.26941693141632095+0j) [Z3]
このハミルトニアンから、Hartree-Fock (HF) エネルギーを計算してみよう。Jordan-Wigner 変換では、qubitの\(\left|0\right\rangle, \left|1\right\rangle\)と軌道の占有数が 1対1 対応していることから、HFエネルギーを計算するには、下から電子数分だけを詰めていった \(\left|1100\right\rangle\) に対する期待値をとれば良い。
[9]:
#テンソル積を計算するための関数
def kron_N(*ops):
tmp = ops[0]
for op in ops[1:]:
tmp = np.kron(tmp,op)
return tmp
bra0 = np.array([[1,0]])
bra1 = np.array([[0,1]])
HFbra = kron_N(bra1, bra1, bra0, bra0)
HFket = HFbra.T
print(HFbra)
jw_matrix = get_sparse_operator(jw_hamiltonian)
print(np.real(HFbra.dot(jw_matrix.dot(HFket))), molecule.hf_energy)
[[0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0]]
[[-1.11299655]] -1.1129965456691684
pyscf の計算と殆ど一致していることが確認できる。
次にハミルトニアンを対角化して、その結果がFull-CI (厳密解) エネルギーと一致することを確かめてみよう。
[10]:
eigenenergies, eigenvecs = np.linalg.eigh(jw_matrix.toarray())
print(eigenenergies[0], molecule.fci_energy)
-1.129904784322913 -1.1299047843229135
これも殆ど一致していることが確認できる。基底状態の波動関数 \(\left|\psi_g\right\rangle\) は
[11]:
print(eigenvecs[:,0])
[ 0. +0.j 0. +0.j 0. +0.j 0.09545811+0.j
0. +0.j 0. +0.j 0. +0.j 0. +0.j
0. +0.j 0. +0.j 0. +0.j 0. +0.j
-0.99543345+0.j 0. +0.j 0. +0.j 0. +0.j]
すなわち \(\left|\psi_g \right\rangle \approx - 0.995\left|1100\right\rangle + 0.0955 \left|0011\right\rangle\)。Hartree-Fock解 \(\left|1100\right\rangle\) に対して少し二電子励起状態 \(\left|0011\right\rangle\) が混ざっていることがわかる。