@yuichirominato 2019.02.19更新 280views

【量子化学】とりあえず水素分子の基底エネルギー求める


はじめに

量子コンピュータで量子化学ということで手順がわからないということで、初心者の僕が解説してみます。基本ツールを使います。。。

古典量子化学のツール

PySCFやPsi4などのオープンソースツールを使って、hogehogeとデータを作るようです。今回はさらにOpenfermionというツールの中に最初から入ってるデータを使ってみたいと思います。

jupyterもしくはgoogle colaboratoryでやってみます。

インストール

今回はBlueqat,OpenFermionBlueqat,OpenFermionをインストールしました。

!pip install blueqat openfermionblueqat openfermion

!マークがついてるのはgoogle colabでやってるからです。インストールできましたか?早速読み込んでみます。

ちなみにBlueqatはこちら、

https://github.com/mdrft/Blueqat

OpenFermion-Blueqatはこちらにあります。

https://github.com/mdrft/OpenFermion-Blueqat

from blueqat import *
from openfermion import *
from openfermionblueqat import*

データの読み込み

早速よみこんでみます。MoluculearDataというのでOpenFermionに最初から入っているデータを使えるようです。H2でSto-3g系を読み込んでみます。

def get_molecule(bond_len):
  geometry = [('H',(0.,0.,0.)),('H',(0.,0.,bond_len))]
  
  description = format(bond_len)
  molecule = MolecularData(geometry, "sto-3g",1,description=description)
  
  molecule.load()
  return molecule

チェックしてみると、原子間距離が0.4のものを読み込んで、

m = get_molecule(0.4)

そのあと、

m.get_molecular_hamiltonian()

これによってハミルトニアンと呼ばれる基本の式が取得できます。

() 1.322943021475
((0, 1), (0, 0)) -1.4820918858979102
((1, 1), (1, 0)) -1.4820918858979102
((2, 1), (2, 0)) -0.1187350527865787
((3, 1), (3, 0)) -0.1187350527865787
((0, 1), (0, 1), (0, 0), (0, 0)) 0.36843967630348756
((0, 1), (0, 1), (2, 0), (2, 0)) 0.08225771204699692
((0, 1), (1, 1), (1, 0), (0, 0)) 0.36843967630348756
((0, 1), (1, 1), (3, 0), (2, 0)) 0.08225771204699692
((0, 1), (2, 1), (0, 0), (2, 0)) 0.082257712046997
((0, 1), (2, 1), (2, 0), (0, 0)) 0.3626667179796745
((0, 1), (3, 1), (1, 0), (2, 0)) 0.082257712046997
((0, 1), (3, 1), (3, 0), (0, 0)) 0.3626667179796745
((1, 1), (0, 1), (0, 0), (1, 0)) 0.36843967630348756
((1, 1), (0, 1), (2, 0), (3, 0)) 0.08225771204699692
((1, 1), (1, 1), (1, 0), (1, 0)) 0.36843967630348756
((1, 1), (1, 1), (3, 0), (3, 0)) 0.08225771204699692
((1, 1), (2, 1), (0, 0), (3, 0)) 0.082257712046997
((1, 1), (2, 1), (2, 0), (1, 0)) 0.3626667179796745
((1, 1), (3, 1), (1, 0), (3, 0)) 0.082257712046997
((1, 1), (3, 1), (3, 0), (1, 0)) 0.3626667179796745
((2, 1), (0, 1), (0, 0), (2, 0)) 0.36266671797967454
((2, 1), (0, 1), (2, 0), (0, 0)) 0.08225771204699726
((2, 1), (1, 1), (1, 0), (2, 0)) 0.36266671797967454
((2, 1), (1, 1), (3, 0), (0, 0)) 0.08225771204699726
((2, 1), (2, 1), (0, 0), (0, 0)) 0.08225771204699728
((2, 1), (2, 1), (2, 0), (2, 0)) 0.38272169831413727
((2, 1), (3, 1), (1, 0), (0, 0)) 0.08225771204699728
((2, 1), (3, 1), (3, 0), (2, 0)) 0.38272169831413727
((3, 1), (0, 1), (0, 0), (3, 0)) 0.36266671797967454
((3, 1), (0, 1), (2, 0), (1, 0)) 0.08225771204699726
((3, 1), (1, 1), (1, 0), (3, 0)) 0.36266671797967454
((3, 1), (1, 1), (3, 0), (1, 0)) 0.08225771204699726
((3, 1), (2, 1), (0, 0), (1, 0)) 0.08225771204699728
((3, 1), (2, 1), (2, 0), (3, 0)) 0.38272169831413727
((3, 1), (3, 1), (1, 0), (1, 0)) 0.08225771204699728
((3, 1), (3, 1), (3, 0), (3, 0)) 0.38272169831413727

でました。これは第二量子化されたものですので、このまま量子コンピュータの演算子に変換できます。

量子コンピュータ向け変換

ブラビキタエフ変換やジョルダンウィグナー変換がありますが、今回は前者のBK変換をしてみると、

h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))

確認してみると、

print(h)

(0.7407724940116754+0j)*I + (0.23528824284103544+0j)*Z[0] + (0.23528824284103542+0j)*Z[0]*Z[1] + (-0.45353118471995524+0j)*Z[2] + (-0.45353118471995524+0j)*Z[1]*Z[2]*Z[3] + (0.18421983815174378+0j)*Z[1] + (0.041128856023498556+0j)*Y[0]*Z[1]*Y[2]*Z[3] + (0.041128856023498556+0j)*X[0]*Z[1]*X[2] + (0.041128856023498556+0j)*X[0]*Z[1]*X[2]*Z[3] + (0.041128856023498556+0j)*Y[0]*Z[1]*Y[2] + (0.14020450296633868+0j)*Z[0]*Z[2] + (0.18133335898983727+0j)*Z[0]*Z[1]*Z[2]*Z[3] + (0.18133335898983727+0j)*Z[0]*Z[1]*Z[2] + (0.14020450296633868+0j)*Z[0]*Z[2]*Z[3] + (0.19136084915706864+0j)*Z[1]*Z[3]

このようにきちんと変換されました。そして、このhをVQEにかけます。

runner = vqe.Vqe(UCCAnsatz(h,2,Circuit().x[0]))
result = runner.run(verbose = True)

そうすると、

params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [1.00803084 0.88310007] val: 0.10544430487391798
params: [-1.61000316  0.88310007] val: 0.48593858308326665
params: [0.00803084 0.88310007] val: -0.3685674394653079
params: [-0.61000314  0.88310007] val: -0.24619628358384726
params: [0.38999684 0.88310007] val: -0.27563514593504546
params: [-0.07664849  0.88310007] val: -0.371744879738929
params: [-0.07664849  0.88310007] val: -0.371744879738929
params: [-0.07664849  1.88310007] val: 1.0541144286439048
params: [-0.07664849 -0.73493393] val: -0.44860542849648677
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.16132782 -0.82515946] val: -0.32853423118712466
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [0.92335151 0.0289703 ] val: -0.3234253290431143
params: [-1.69468249  0.0289703 ] val: 0.7219628224155241
params: [-0.07664849  0.0289703 ] val: -0.8993554487133357
params: [-0.69468246  0.0289703 ] val: -0.5631470828007092
params: [0.30531751 0.0289703 ] val: -0.8360012977529356
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  1.0289703] val: -0.19293944513988606
params: [-0.0024937 -1.5890637] val: 0.5730186565494813
params: [-0.0024937  0.0289703] val: -0.903736928997239
params: [-0.0024937  -0.58906367] val: -0.6556036819150325
params: [-0.0024937  0.4109363] val: -0.7807595396151468
params: [-0.0024937  -0.00256402] val: -0.9043519838381495
params: [-0.0024937  -0.16023561] val: -0.885812805585636

出てきました。パラメータが最適化されながら最小基底状態が探索されています。中で実際にどのような回路が格納されているかはIBMのqasmに変換できます。

result.circuit.to_qasm()

'OPENQASM 2.0;\ninclude "qelib1.inc";\nqreg q[4];\ncreg c[4];\nx q[0];\nrz(0.0036865832224720592) q[0];\ncx q[0],q[1];\nrz(0.003686583222472059) q[1];\ncx q[0],q[1];\nrz(-0.007106094364375361) q[2];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(-0.007106094364375361) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\nrz(0.0028864245674849803) q[1];\nrx(-1.5707963267948966) q[0];\nrx(-1.5707963267948966) q[2];\ncx q[0],q[1];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(0.0006444221298305115) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\ncx q[0],q[1];\nrx(1.5707963267948966) q[0];\nrx(1.5707963267948966) q[2];\nh q[0];\nh q[2];\ncx q[0],q[1];\ncx q[1],q[2];\nrz(0.0006444221298305115) q[2];\ncx q[1],q[2];\ncx q[0],q[1];\nh q[0];\nh q[2];\nh q[0];\nh q[2];\ncx q[0],q[1];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(0.0006444221298305115) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\ncx q[0],q[1];\nh q[0];\nh q[2];\nrx(-1.5707963267948966) q[0];\nrx(-1.5707963267948966) q[2];\ncx q[0],q[1];\ncx q[1],q[2];\nrz(0.0006444221298305115) q[2];\ncx q[1],q[2];\ncx q[0],q[1];\nrx(1.5707963267948966) q[0];\nrx(1.5707963267948966) q[2];\ncx q[0],q[2];\nrz(0.002196776014430724) q[2];\ncx q[0],q[2];\ncx q[0],q[1];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(0.0028411981442612363) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\ncx q[0],q[1];\ncx q[0],q[1];\ncx q[1],q[2];\nrz(0.0028411981442612363) q[2];\ncx q[1],q[2];\ncx q[0],q[1];\ncx q[0],q[2];\ncx q[2],q[3];\nrz(0.002196776014430724) q[3];\ncx q[2],q[3];\ncx q[0],q[2];\ncx q[1],q[3];\nrz(0.002998312569392093) q[3];\ncx q[1],q[3];\nrz(0.0037905403357534834) q[0];\ncx q[0],q[1];\nrz(0.003790540335753483) q[1];\ncx q[0],q[1];\nrz(-0.007306477486699262) q[2];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(-0.007306477486699262) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\nrz(0.0029678181906944716) q[1];\nrx(-1.5707963267948966) q[0];\nrx(-1.5707963267948966) q[2];\ncx q[0],q[1];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(0.0006625940414107214) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\ncx q[0],q[1];\nrx(1.5707963267948966) q[0];\nrx(1.5707963267948966) q[2];\nh q[0];\nh q[2];\ncx q[0],q[1];\ncx q[1],q[2];\nrz(0.0006625940414107214) q[2];\ncx q[1],q[2];\ncx q[0],q[1];\nh q[0];\nh q[2];\nh q[0];\nh q[2];\ncx q[0],q[1];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(0.0006625940414107214) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\ncx q[0],q[1];\nh q[0];\nh q[2];\nrx(-1.5707963267948966) q[0];\nrx(-1.5707963267948966) q[2];\ncx q[0],q[1];\ncx q[1],q[2];\nrz(0.0006625940414107214) q[2];\ncx q[1],q[2];\ncx q[0],q[1];\nrx(1.5707963267948966) q[0];\nrx(1.5707963267948966) q[2];\ncx q[0],q[2];\nrz(0.002258722396542493) q[2];\ncx q[0],q[2];\ncx q[0],q[1];\ncx q[1],q[2];\ncx q[2],q[3];\nrz(0.002921316437953215) q[3];\ncx q[2],q[3];\ncx q[1],q[2];\ncx q[0],q[1];\ncx q[0],q[1];\ncx q[1],q[2];\nrz(0.002921316437953215) q[2];\ncx q[1],q[2];\ncx q[0],q[1];\ncx q[0],q[2];\ncx q[2],q[3];\nrz(0.002258722396542493) q[3];\ncx q[2],q[3];\ncx q[0],q[2];\ncx q[1],q[3];\nrz(0.0030828612966605916) q[3];\ncx q[1],q[3];'

最小基底エネルギーを確認してみると、

runner.ansatz.get_energy(result.circuit,runner.sampler)

-0.9043519838381495

求まりました。このように0.4でやったものを全体でやってみます。

グラフを作る

最初からの手順をループにするだけです。下記の手順で計算ができます。

from blueqat import *
from openfermion import *
from openfermionblueqat import*
import numpy as np

x = [];e=[];fullci=[]
for bond_len in np.arange(0.2,2.5,0.1):
  m = get_molecule("{:.2}".format(bond_len))
  h = bravyi_kitaev(get_fermion_operator(m.get_molecular_hamiltonian()))
  runner = vqe.Vqe(UCCAnsatz(h,6,Circuit().x[0]))
  result = runner.run()
  x.append(bond_len)
  e.append(runner.ansatz.get_energy(result.circuit,runner.sampler))
  fullci.append(m.fci_energy)

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(x,fullci)
plt.plot(x,e,"o")

このようにFullCIと比較した結果が簡単に取得できました。計算時間はかかりますが、このようにして少しずつ大きな分子で計算をすることができます。水素はすぐ終わります。

ファイルは下記にあげましたのでjupyterなどで実行できます。詳しい話は僕はあまりよくわからないので、開発者のgyu-donさんにきいてみてください。

https://github.com/minatoyuichiro/openfermionblueqat_h2/blob/master/h2.ipynb

以上です。

Recommended


Wikiへ移動