4E11

Fock 近似における原子構造計算のための高精度数値計算法
石川 英明
Highly accurate, numerical methods for atomic structure calculationsin the Fock approximation
Hideaki Ishikawa【序】
原子構造計算では、原子の電子状態（固有値と固有関数）を計算する。原子は一般に多電子系である。平均場近似では以下の定式化を行う。(1) 系の全波動関数を 1 電子波動関数の積和或いは積和の線形結合で近似する。(2) 全波動関数を用いた全エネルギーの期待値を１電子波動関数で表す際、電子間の Coulomb 反発エネルギーを１次摂動で取り込む。(3) 規格直交性の束縛条件
下で全エネルギーの期待値を 1 電子波動関数に関して変分することにより、固有値、固有関数とポテンシャルの連立非線形微分方程式を導く。これらの方程式はセルフ・コンシステントに解く問題に帰着される [1-3]。
我々はこれまで最も簡単な Hartree 近似で高精度の数値計算法を提示してきた[4]。現実の系を扱うには、近似の精度を更に高める事が必要である。その第１歩として、Pauli の排他原理を満たすよう、全電子波動関数を１電子波動関数から成る１個の行列式（Slater 行列式）で近似する（単一配置）。この全電子波動関数を用いて全エネルギーの表式を求め（Slater による[1,2]）、更にそれを１電子波動関数に関して変分することにより得られる連立非線形微分方程式（Fock 方程式）を扱う。本報告では単一配置の Fock 方程式を解くために必要な高精度数値計算法を述べる。
【Fock 方程式の特徴】
問題の設定：
典型例として、Fock 近似で最も簡単なケースである閉殻の場合の Fock 方程式を以下に示す（文献[1]の p. 54 と p. 50）

 d2 2l (l + 1) 
 2 + Y (nl ; r ) − ε nl ,nl − P(nl ; r ) = X (nl ; r ) + ∑ ε nl ,n 'l P(n ' l ; r ),r
r2 
n '≠ n dr
Y (nl , r ) = N − ∑ q(nl )Y0 (nl , nl; r ) + ∑ α lk Yk (nl , nl ; r ),nlk

α l 0 = 1, α lk = 2 Alk / q(nl ),X (nl ; r ) = −(2 / r )∑ n 'l ' k β ll ' k Yk (nl , n ' l '; r ) P(n ' l ', r ),'β ll ' k = Bll ' k / q (nl ),kk +1∞ r s
Yk (nl , n ' l '; r ) = ∫   P (nl ; s ) P (n ' l ' : s)ds + ∫   P(nl ; s) P(n ' l '; s )ds.
s =0 r
s=r s
 
 
r

ここで、P(nl;r)は 1 電子の動径波動関数、n と l はそれぞれ全動径量子数と方位量子数、Y(nl;r)は全ポテンシャルの係数、εnl,nl はエネルギー・パラメータの対角成分、N は原子番号、q(nl)は状態 nl の占有数、Alk と Bll’k は定数で、閉殻の場合の値は表に示されている。Fock 近似では、交換積分の存在により、Y(nl;r)の中には関数 Yk(nl,n’l’;r)の k = 0 の成分のみならず k > 0の成分も出現する。更に他の状態 P(n’l’;r)と関係した（非対角の）項として、右辺第 1 項のX(nl;r)と第 2 項が現れる。Hartree の方程式は、Y(nl;r)の中で Yk(nl,n’l’;r)の k = 0 のみを残し、かつ右辺の非対角のポテンシャルを含む項を落とし、非対角のパラメータεnl,nl’をゼロに置いたものである。数学的な観点からは、Fock 方程式は非同次項を含む一般化された固有値
方程式（或いは単に境界値問題の微分方程式とも呼ばれる）であり、Hartree 方程式は非同次項がない固有値方程式である。
境界条件は P(nl;r=0) = 0 及び r → ∞で P(nl;r)→0 である。
【数値解法】
Fock 方程式を解くためには、全ポテンシャルと非対角ポテンシャルを計算すること、一般化された固有値方程式を解くこと、及びこれらをセルフ・コンシステントになるまで計算を繰り返すことが必要である。ポテンシャル計算の要点は、構成要素である Yk(nl,n’l’;r)を定義に従って波動関数 P(nl;r)から計算することである。これについては、k > 0 の場合も含めて、これまでに高精度の数値計算法を確立している。また、Yk(nl,n’l’;r)を合成して全ポテン
シャル等を計算するために必要な電子間クーロン反発エネルギーの 2 電子積分でクーロン積
分及び交換積分に現れる Slater の a 係数と b 係数、及びそれらを計算するための c 係数[2]を数式処理で求め、文献の表にある数値を確認した。（これらは係数 Alk, Bll’k に反映される[1]。）
Fock 方程式が非同次項を含むことから、同次方程式である Hartree 方程式で開発した動径固有値問題の数値計算法をそのまま使うことはできない。このため、数値計算法を拡張する必要がある。これまで、種々の数値計算法が提出されてきたが[1,3]、数学的にスッキリと表現された方法は、広義の（或いは一般化された）Green 関数を用いる方法である[5]。この仕事は 1930 年代（手回し計算機の時代）に行われた。現代では、計算機と数値計算法はこの時代のものに比べて格段に進歩している。この一般的な方法に従って数値解を得るため、数値計算法の開発を進めている：原点近傍での線形独立な二つの冪級数解、微分方程式の初期値
問題の数値解、無限遠点近傍での線形独立な二つの漸近展開の解、及びそれらの接続方法、
種々のパラメータの更新法、等。方法と結果は当日発表する。
参考文献[1] D. R. Hartree, The Calculation of Atomic Structures, Wiley, 1957.
[2] J. C. Slater, Quantum Theory of Atomic Structure, 2 vols., McGraw-Hill, 1960.
[3] C. Froese Fischer, The Hartree-Fock Method for Atoms, Wiley, 1977.
[4] 石川英明、第 6 回分子科学討論会 2012 講演要旨 1E05.
[5] V. Fock, M. J. Petrashen, Z. Sowjetunion 6 (1934) 368-414, reissued in Selected Works V. A.
Fock, Quantum Mechanics and Quantum Field Theory, edited by L. D. Faddeev, L. A. Khalfin, and
I. V. Komarov, Chapman & Hall/CRC, 2004.

4E12
NTChem: 超並列環境に資する分子科学計算ソフトウェアの開発(理研･計算科学研究機構) 中嶋隆人
NTChem: a high-performance software package for molecular electronic structure calculation(RIKEN AICS) Takahito Nakajima
量子化学に基づいた分子科学計算ソフトウェアは物質科学･生命科学･モノづくりなどの多くの分野の共通基盤である．Gaussian (米国) や GAMESS (米国) などの量子化学計算ソフト
ウェアは計算科学者や実験科学者を問わず世界中で幅広く利用されている．コンピュータの高度化･高性能化に伴い，大規模な分子系を高精度に計算する要望は急速に増しつつある．この要望に応えるため，超並列計算が可能で一般ユーザーが利用できる国産分子科学計算ソフ
トウェアの開発が国際的にも急務である．広く用いられている既存ソフトウェアの多くは単一プロセッサの時代に設計･開発されたものであり，その拡張としての単純な並列化は可能であるが，京コンピュータのような超並列スーパーコンピュータにおいて並列化効率が問題となる．全系丸ごと計算に関して言えば，千原子分子程度の第一原理電子状態計算が上限で，新規の機能発現などが期待できるナノスケールサイズの数万原子系に対する計算は不可能である．また，数百原子分子系の第一原理化学反応計算も多数の電子状態計算が必要であり，ほぼ困難である．FMO 法のような領域分割法を使っても並列化効率の点では数万原子系が限度であろう．これらの問題について，分子化学分野をはじめとする多くの領域において京コ
ンピュータを利用した解決の強い要望が出ている．このため，理化学研究所･計算科学研究機
構･量子系分子科学研究チームでは以下に示す研究を行い，幅広い分野の多くのユーザーの利用に資する汎用分子科学計算ソフトウェア「NTChem」の開発を行ってきた．
(1) 他の量子化学計算プログラムでは扱うことができない大規模分子向け分子計算法，高速計算法，高精度分子計算法の理論とそのアルゴリズムを開発する．
(2) それらの理論･アルゴリズムに基づき，我が国独自の分子科学計算ソフトウェアを新たに開発し公開することで，タンパク質の丸ごと計算のような数万原子分子の電子状態計算や，数百~千原子分子系の化学反応の詳細な追跡計算を高速かつ高精度に実現する．
「NTChem」は一から設計をした新しい国産分子科学計算ソフトウェアである．既存ソフト
ウェアの持つ多くの機能をカバーしつつ，われわれが新たに開発してきた理論手法の集大成でもあって，他のプログラムでは利用することのできない多くの量子化学計算法を含んでいる．
「NTChem」の第一版には数千原子分子系に対する第一原理電子状態計算や数百原子分子系の化学反応過程追跡計算を実現するための分子科学理論が実装されている．さらに，京コンピュータなどのマルチコア超並列クラスタ計算システムの性能を引き出すことが可能な並列アルゴリズムが実装されている．
「NTChem」の第一版の主な機能は以下の通りである．
(1) Hartree–Fock 法および密度汎関数 (DFT) 法に基づく基底状態の電子状態計算(2) 時間依存密度汎関数法と遷移ポテンシャル法に基づく励起状態の電子状態計算(3) GFC 法，RIDFT 法，Dual-level DFT 法などによる線形あるいは低スケーリング DFT 計算(4) PDM や擬対角化法などの対角化フリー法による低スケーリング SCF 計算(5) Coupled-cluster 法と量子モンテカルロ法に基づく基底状態および励起状態の高精度電子相
関計算(6) 高並列化 RI-MP2 法による高精度大規模分子電子相関計算(7) 高次 Douglas–Kroll 法や RESC 法などに基づいたスピン–軌道相互作用を含む 2 成分相対論
的電子状態計算(8) QM/MM 法や ONIOM 法による巨大分子のモデル計算(9) 効率良い化学反応経路探索計算(10) ab initio 分子動力学法計算(11) NMR 化学シフトや電子スピン共鳴などの磁気的性質の計算(12) Maximally interaction orbital (MIO) 法や相互作用軌道 (PIO) 法による分子間相互作用解析
本年秋には京コンピュータ上において「NTChem2013」を一般公開する予定である．本講演
では「NTChem」の特徴，機能，性能に関し実例を踏まえながら紹介する．

4E13
漸近展開法による二電子積分の高速計算（イリノイ大学アーバナ・シャンペーン校 1、JST-CREST 2）○山田朋範 1,2、Ryan Brewster 1、平田聡 1,2
Asymptotic expansion for fast computation of two-electron integral(University of Illinois at Urbana-Champaign1, JST-CREST1,2)○Tomonori Yamada 1,2, Ryan Brewster 1, So Hirata1,2【序】二電子積分の計算はハートリー・フォック（ＨＦ）法において最もコストのかかる箇所である。特に大規模分子や周期系（固体）の計算において二電子積分はその計算時間、およびデータを保存するディスク容量の両面において大きなボトルネックとなる。３０年以上前からそのコストを少なくするための多くの研究がなされてきた[1,2]。我々は本研究で二電子積分の新しい近似法、漸近展開法を提案する。この近似法は二つの電子が互いに離れているほど精度が良く、大規模系の計算コストを大幅に削減する。本発表ではこれを用いて１次元、および２次元周期系のエネルギーをＨＦ法により計算した例を紹介する。
【方法】以下では一次元周期系の二電子積分に漸近展開法を適用した場合の説明をする。
一次元周期系の原子基底二電子積分は
0) )|  )  )) =    ∗  )   ) ∗   )   )


（１）

と書ける。ここで  −  )は番目の単位格子に属する番目の基底関数であり、は格子ベクトルを表す。式（１）に漸近展開法を適用すると
0) )|  )  )) ≈
)*)

%)
$ )


%) '$ ) $ ) '$ )
 

 &
&



#$ ) #$ )（２）となる。ここで( ) は重なり積分、  ) は  −  )の中心座標である。式（２）の右辺は、式（１）の右辺を多極子展開し、第一項目のみを考慮した式である。またこれは、仮
,*)

. )

想的な二つの電荷、(+ ) と(-) 、の間のクーロン相互作用とも解釈できる。本手法は、既に開発されている多極子展開に基づいた手法[1,2]と似ているが以下の点で異なる。これらの手法は近似を適用するために空間をボックスで区切り、このボックスごとに近似を適用する。
これに対し、本手法は原子基底ごとに近似を行う。これによりクーロン積分や核－電子間引
力積分だけでなく、交換積分にも同じ精度の近似を用いることができ、また同じアルゴリズ
図１に式（１）と（２）で計算された一次元水素鎖の 1s

Two-Electron integral/ a.u.

ムで容易にプログラムを書くことができる。

図1

2.0
(1)
(2)

1.5

軌道の二電子積分の比較を示す。横軸は二電子間の距離の指標となる m2 を示している。漸近展開法は二電子間の距離が離れているときに良い近似となるので、m2 が大きい場

1.0
0.5
0.0
0

5

10
m2

15

20

合にのみ適用する。これを Coulomb 積分
,*)

- )

/+) = 0.,- 0 0 1.) 0) )|  )  ))

（３）

- )
を計算する場合で説明する。ここで1.) は、電子密度の行列要素である。m2 についての和をとるとき、パラメータ L と X を導入し、| | ≤ 4では（１）を用いて 2 電子積分を計算し、4 < | | ≤ 6 では式（２）を用いる。交換積分や核-電子間引力積分の計算も同様に行う。

X はより大きい方が計算そのものが信頼でき、また L はどの位置から式（２）による近似を適用するかの指標になっていることに注意されたい。
【結果】漸近展開法を用いて計算した結果を以下に示す。表１と表２に 1 次元と 2 次元周期系のエネルギー（E）と計算コストをそれぞれ示す。計算コストには、二電子積分計算とフォック行列の対角化を実行するのに要する計算時間と、二電子積分の保存に必要なディスク容量を示す。前者は SCF のサイクル毎に二電子積分を直接計算する場合を想定した計算時間であり、後者は二電子積分をハードディスクに保存してＳＣＦの各サイクルで読み込む方法を想定した場合に必要なディスク容量である。計算した 1 次元系は水素鎖、およびフッ化水素
鎖で、それぞれ単位格子に H4 と HF を含んでおり、1d-(H4)x、1d-(HF)x と表記する。2 次元系は水素原子シート（2d-(H4)x）とグラフェン（2d-(C2)x）である。
漸近展開法を用いなかった場合（L=X）と用い用いた
用いた場合（L<X）を比較すると、

表1

L

X

きくコストが削減されていることがわ

100

100

かる。また同じ X の値でそれらのエネ

40

どの系においても漸近展開法により大ルギーを比較すると、漸近展開法を用いない場合と近い値を与えることがわかる。
1d-(H4)x の系では電子密度の距離に対する減衰が遅く、L=X=10 の計算では長距離の電子密度を切り捨てているE (a.u.)CPU

Disk (GB)1d-(H4)x [HF/6-31G(d,p)]-1.82822

1475

31

40

-1.82800

598

13

156

3

183

3

10

10

Not-convergent

10

100

-1.82823

100

100-99.99350

1190

9

40

40

-99.99350

482

4

10

10

-99.99339

132

1

10

100

-99.993501401

1d-(HF)x [HF/6-31G(d,p)]2d-(H4)x [HF/STO-3G]

ことが原因でＳＣＦが収束しない。こ40

40-1.31893

26096

103

れは言い換えれば長距離の交換積分を

30

30-1.81890

14788

60

20

20

-1.31889

13081

28

10

40

-1.318944356

8

無視していることが原因である。
1d-(HF)x のようにイオン性が強い系で

2d-(C4)x [HF/3-21G]

は、長距離の静電相互作用を多極子展

1515-75.33388

7473

129

開法[2]で計算するのが非常に有効であ

88-75.33387

4850

87

る。しかし長距離の交換積分が重要に

4

4-75.33372

673

14

8

15

-75.33388

3206

43

なる 1d-(H4)x などの場合にはこれは適用できず、本手法が有効になる。
（参考文献）[1] C.A. White, B.G. Johnson, P.M.W. Gill, M. Head-Gordon, Chem. Phys. Lett. 230 (1994) 8.
[2] J. Delhalle, L. Piela, J.-L. Bredas, J.-M. Andre, Phys. Rev. B 22 (1980) 6254.

DMRG 線形応答理論と post-DMRG 法：

4E14

Tamm-Dancoff 近似による強相関系励起状態計算への応用(プリンストン大学 1，ゲント大学 2) ○中谷直輝 1，Sebastian Wouters2，Garnet K.-L. Chan1

DMRG linear response theory and post-DMRG theory:
Tamm-Dancoff approximation for the excited state problem in strongly correlated systems(Princeton University1, Ghent University2) Naoki Nakatani1, Sebastian Wouters2, and Garnet K.-L. Chan1【諸言】密度行列繰込み群(DMRG)[1]は、多体波動関H[i ] A ni  ESAni数を効率よく記述できる新しい方法論として注目され、

ここで H[i]はサイト i 以外からの寄与が繰込まれた有強相関系のスピン状態や電子状態を計算する上で非常

効 Hamiltonian である。DMRG 法はこの式に用いて行に有効なツールとして用いられている。

列 An を反復的に最適化するための手法である。

DMRG 法はさまざまな強相関系の基底状態を研究するために用いられている一方で、励起状態への応用

同様に、時間依存の変分原理を適用することで、MPS波動関数に対する運動方程式が得られる[4]。

例は少ない。これは励起状態を求める際、状態平均密iS [t i ]度から繰込み行列を計算するため(SA-DMRG)、変分条
件が決まらないことや、各状態に有効な繰込み基底が減少することによって必要な繰込み次元が大きくなる

次に、繰込み行列 An の時間変化を 0 次の行列要素からの調和振動で展開することを考え、Atni  Xni eit  Yni *eitなどの問題があるためである。
本研究では、時間依存の DMRG 方程式から線形応答理論を用いて励起状態を効率的に計算するための手法を開発し、強相関系の励起状態計算への応用を行った。
また、Hartree-Fock 法との相似性から DMRG 波動関数

A tni H [ti ] A tnit

これを運動方程式に代入して 1 次までの項で近似すると RPA 方程式が得られる[3, 5]。

 H
 *
W

W  X 
0  X S   ω
 
* 
* 
H  Y  0  S  Y を参照波動関数とする post-DMRG 法への展開につい

ここで、行列要素 H、W は RPA における有効ハミルて概説する。

トニアン、S は重なり行列で、それぞれ M2k×M2k の【DMRG 線形応答理論】DMRG 波動関数はサイトと呼行列である。また、状態ベクトル X、Y は An の 1 次変

ばれる局所的な Fock 空間の基底に独立な行列の積に化の集合であり、M2k の次元を持つ。本研究では、簡よって与えられ、行列積状態(MPS)と呼ばれる[2]。

単のため RPA 方程式の実部 H、 X のみを考える0   A n1  A ni  A nk n1  nk{n}(Tamm-Dancoff 近似)。
この方程式の解は M2k×M2k の行列を対角化するこ

ここで、k はサイトの数、ni はサイト i における局所とで得られ、その固有値は励起エネルギーを与える。

Fock 空間の基底、An は n に独立な M×M の行列であ

このような巨大次元の行列の対角化には Davidson のる。MPS 波動関数はサイトに独立な平均場近似を与え

部分対角化アルゴリズムが有効である。本研究では、ており、変分原理を適用することで導かれる方程式は、部分空間の有効ハミルトニアンを DMRG の sweep アル

粒子 の運動に 対する平 均場近似 の結果与 えられる

ゴリズムにより計算することで、基底状態と同じく

Hartree-Fock(HF)方程式と相似的な特徴を持つことがO(M3k3+M2k4)のコストで計算するアルゴリズムの開発

知られている[3]。

を行った。

【post-DMRG 法】MPS 波動関数がサイトに関する平均場近似であることは述べたが、ここでは、サイト間の相関効果を後から取り込むことを考える。これにより
MPS の 1 次元格子の entanglement を効率良く改善する
新たな理論の構築ができると考えられ、より小さな繰込み次元で高精度な波動関数を得ることができるよう

T able 1. Energy errors for lowest 3 excited states of C16 H18from converged calc. DMRG M = 1000 (in mH)

DMRG-TDA
state
M = 25

50

100

150

200

1

X Ag

0.7

0.1

0.0

0.0

0.0

21 A g

18.4

6.4

1.4

0.9

0.7

1

13.4

4.3

0.8

0.5

0.3

1 Bu

になると期待される。
DMRG の繰込み行列 An は、左側の繰込み基底{l}とサイトに独立な Fock 空間の基底{n}の積空間から右側の繰込み基底{r}への射影を与える。ここで、{r}の補空間への射影を Bn として、サイトに関する励起演算子
B+を、0 次の繰込み行列 An を Bn へと変化させるように定義する。

8SA-DMRGstate
M = 25

50

100

150

200

X1 A g

17.4

2.2

0.6

0.2

0.1

1

2 Ag

73.0

10.8

1.6

0.5

0.2

11 Bu

64.2

9.4

1.4

0.5

0.2

T able 2. Energy errors of low-lying states for Hubbard model(L = 8, U/t = 10) from FCI (in mH)

Bˆ i 0   A n1 B ni  A nk n1  nk{n}

このような励起演算子を用いることで、MPS 波動関数を参照する多体相関理論を容易に導出することができる。例えば、CISD や CCSD 波動関数は次のような形で与えられる[5]。

1CISD    Ci Bˆ i   Dij Bˆ i Bˆ j  02CCSD  e1  Ci Bˆ i   Dij Bˆ i Bˆ j20

【応用】C16H18 の励起状態について SA-DMRG とDMRG-TDA による π-valence CI (16e, 32o)/cc-pVDZ の計算結果を Table 1 に示す。

S

N

TDA M = 3 TDA M = 9 CISD M = 3

0

4

104.2

0.5

0.3

1/2

5

84.9

0.8

2.4

1

4

117.8

1.5

3.4

1/2

5

120.7

9.3

7.2

以上の結果から、DMRG 波動関数を参照波動関数とする post-DMRG 法が多体相関理論の新たな展開として非常に期待できると考えられる。
今後は、DMRG-CISD や DMRG-CCSD を行うための実用的なアルゴリズムの開発と分子系への応用を目指す予定である。

繰込み次元 M が小さい範囲では DMRG-TDA の方が良い結果を与えている。しかし、M が大きい範囲では
SA-DMRG の方が良い結果を与えており、DMRG-TDA【参考文献】[1] S. R. White, Phys. Rev. Lett. 69, 2863 (1992)[2] U. Schollwöck, Ann. Phys. 326, 96 (2011)の M に対するエネルギーの収束が遅いことが分かる。
これは、DMRG-TDA では 2 サイト以上の相関を取り[3] J. M. Kinder, C. C. Ralph, G. K.-L. Chan, Adv. Chem.
Phys. in press (2013)

込むことが出来ないためであると考えられる。一方で[4] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pizorn, H.

2 サイト以上の相関を取り込むことで M を小さく抑え

Verschelde, F. Verstraete, Phys. Rev. Lett. 107, 070601

つつ、エネルギーを高精度に計算できる可能性も示唆(2011)される。
Table 2 には Hubbard モデルを用いた DMRG-TDA とDMRG-CISD のテスト計算の結果を示す。この結果から 2 サイトの相関を取り込むことによって TDA と比較して結果が劇的に改善していることが分かる。

[5] S. Wouters, N. Nakatani, G. K.-L. Chan, Phys. Rev. B.
in press (2013)4E15
モデル空 間 量 子 モンテカルロ法 I 基本的なアイデアと予備的計算（神 戸 大 ）天 能 精 一 郎
Model space quantum Monte Carlo method I: Basic ideas and pilot applications(Kobe univ.) Seiichiro Ten-no	
 最近、スレーター行列式を用いた配位空間での量子モンテカルロ法が注目を集めている[1,2]。これらの手法は、通常の完全 CI 計算ではメモリーの制限で取り扱えない大次元の CI 係数を離散的なウォーカーの統計分布で表現し、(1.1)の虚時間発展に従って全ての励起状態の寄与を指数関数的に減衰させ、厳密な基底状態の波動関数とエネルギーを得るものである。	
 	 (1.1)	 Ψ(τ ) = e−τ ( Ĥ − E ) Ψ(0) 	
 完全基底関数極限での完全 CI 解を用いるための F12 法[3,4]や、周期系への拡張による固体のスピンギャップの厳密解の計算[5]が行われており、従来手法では定性的な記述すら困難な電子状態の計算が可能になってきている。一方、これらの量子
モンテカルロ法では、完全 CI に近いエネルギーが現実の計算時間内で得られる反面、擬縮重を伴う相関の強い電子系では符号問題が顕著になる事や、励起状態の計算に適用出来ないという問題がある。本研究では、CI 係数の代わりに有効ハミルト
ニアンを確率的に求めるモデル空間量子モンテカルロ法(MSQMC 法)[6]を提案し、擬縮重を含む基底状態と励起状態を取り扱い可能にする。	
 	
  擬縮重電子状態を取り扱う有効ハミルトニアンの構築には、Löwdin 分割法と呼ばれるエネルギー依存分割法(EDP)と、行列型 Bloch 方程式を用いたエネルギー非依存分割法(EIP)とがある。EIP は侵入状態の問題(Intruder	
 s tate	
 p roblem)を伴う事が多いため、ここでは EDP を用いた量子モンテカルロ法の定式化を行う。EDP では、与えられたエネルギー近傍の電子状態について、	
 
	
 

⎛ H
PP
⎜
⎜⎝ HQP⎛ C ⎞H PQ ⎞ ⎛ C P ⎞P⎟⎜⎟ = E⎜⎟ 	 
HQQ ⎟⎠ ⎜⎝ CQ ⎟⎠⎜⎝ CQ ⎟⎠

(1.2)	
 モデル空間(P-空間)の CI 係数とその外側(Q-空間)の CI 係数がトランスファー行列によって CQ = TQPC P のように関係付けられ、モデル空間での有効ハミルトニアンの固有値問題に帰着する。	
 	 (1.3)	 Heff
C = EC P 	
 
PP P
eff
H PP = H PP + H PQ TQP 	 	 (1.4)	
 トランスファー行列は線形方程式により求める事が出来る。	
 −(HQQ − I QQ E)TQP − HQP = 0 	 	 (1.5)	 一方、(1.1)から得られる CQ の虚時間発展	 	 dCQ (τ )dτ= −(HQQ − I QQ E)CQ − HQPC P 	 (1.6)	
 に CQ = TQPC P を代入し、 C P の和を分解すれば、十分条件としてトランスファー行列

の虚時間発展が得られる。	
 	
 dTQP (τ )dτ= −(HQQ − I QQ E)TQP − HQP 	 (1.7)	
 この定常状態 (dTQP (τ ) / dτ = 0) は(1.5)に他ならず、EDP の解である事が分かる。MSQMC
法 で は 、 右 辺 の 寄 与 を ３ つ に 分 割 し 、 1)diagonal	  death/cloning,	
  2)Q-space	 spawning,	
 3 )P-space	 s pawning の各ステップによってトランスファー行列の振幅に対するウォーカーの分布密度ダイナミ
クスを行い、マクロステップ毎の有効
ハミルトニアンの対角化により状態エ
ネルギーの更新を行う。任意の励起状
態が状態エネルギーを変える事により計算可能である。	
 	
  図１に DZP 基底を用いた H4 モデルに対するエネルギーの虚時間発展を示す。
線 形 の 構 造 で は | (1a1 )2 (1b2 )2 | の 単 一 配 置が主配置であるため、全ての方法でエネルギーの分散が小さい。一方、正方構造では | (1a1 )2 (1b2 )2 | と | (1a1 )2 (2a2 )2 | の二つの配置が競合する擬縮重状態となる
図 1.	
 各 種 QMC 法 で の エ ネ ル ギ ー 虚 時 間 発 展 	
 ため、従来の FCIQMC 法やモデル空間の次元を１に制限した MSQMC/1 では、符号問題が顕著になり、統計誤差が大きくなる事が分かる。MSQMC/SPD では、44次元のモデル空間が動的に発生され、非常に速やかに完全 CI エネルギーに収束している。同様に N 2 分子の核間距離を変えたエネルギー誤差を図２に示す。
MSQMC/SPD 法 で は 全 て の 核 間 距 離 で 基
底 状 態 と 5 Σ +g へ の 励 起 エ ネ ル ギ ー を 定
図 2.	
 N 2 分 子 の 誤 差 と 原 子 間 距 離 の 関 係 	
 量的に計算可能である。	
 	
  MSQMC 法は、擬縮重を含むあらゆる励起状態についてほぼ厳密解を得る事が可能な手法であり、これまでの分子軌道計算で取り扱いが困難であった多くの現象に切り込めるポテンシャルを持っている。又、伝播はモデル空間の行列式に対して独立である事を利用して、大塚らによる超並列実装も行われている。	
 【 参 考 文 献 】 	
 1. Y. Ohtsuka and S. Nagase, Chem. Phys. Lett., 463, 431 (2008).
2. G. H. Booth, A. J. W. Thom and A. Alavi, J. Chem. Phys., 131, 054106 (2009).
3. Y. Ohtsuka and S. Ten-no, AIP Conf. Proc., 1456, 97 (2012).
4. G. H. Booth, D. Cleland, A. Alavi, and D. P. Tew, J. Chem. Phys., 137, 164112 (2012).
5. G. H. Booth, A. Grueneis, G. Kresse, and A. Alavi, Nature, 493, 365 (2013).
6. S. Ten-no, J. Chem. Phys., 138, 164126 (2013).

4E16
モデル空間量子モンテカルロ法 II hybrid 並列実装といくつかの応用例（神戸大システム情報） ○大塚 勇起, 天能 精一郎

Model space quantum Monte Carlo method II Hybrid parallel
implementation and some applications（Kobe Univ.） Yuhki Ohtsuka and Seiichiro Ten-no【序】 モデル空間モンテカルロ法（MSQMC 法）[1]は、前講演(4E15)で紹介された通り、スレーター行列式を使用した量子モンテカルロ法(PMC-SD 法[2]や FCIQMC 法[3])の結合
解離領域における統計誤差の増加や、励起状態計算の困難さを解決するために提案された。テスト計算によって、結合解離や励起状態を精度良く安定に計算できることが確かめられているが、計算労力はモデル空間（P-space）に含まれる電子配置の数に比例して増加するため、並列実装が不可欠である。今回、MPI/OpenMP hybrid 並列実装を行い、MSQMC法の高速な計算を可能にした。
【アルゴリズム】 今回行った MSQMC 法の並列実装の概要を図 1 に示す。P-space に含まれる電子配置は、MPI によって各ノードに分配される。モンテカルロ法による P-space からの新しい電子配置の生成（P-space spawning）は、OpenMP を使用して並列化を行った。このス
テップでは、それぞれのウォーカー（サンプル）において独立に生成する確率を計算できるので、並列化は容易である。P-space 以外（Q-space）の電子配置からの spawning（Q-space
spawning）に関しても、同様に OpenMP を使用して並列化した。
Do

I = 1,Do

Macro step

IP =
Do

1,

J=

NP
1,

(P-space determinants) (Distribute using MPI)

Micro stepP-space spawning
Q-space spawning
Death/Cloning step (Q-space only)

OpenMP

Annihilation step
Sigma step for IP(Construct effective Hamiltonian)

Enddo
Enddo
Reduce effective Hamiltonian(MPI)

Diagonalization of effective Hamiltonian
Broadcast Energy

(MPI)

Enddo

図 1.

MSQMC 法の MPI/OpenMP hybrid 並列実装の概要

Death/Cloning step では、Q-space にあるウォーカーを、電子配置のエネルギーに従って、増加や消滅させるが、Q-space spawning と同時に行われるため、自動的に並列化されている。Annihilation step は、Transfer matrix を計算するために、それぞれの電子配置に存在す
るウォーカー数を数える過程であるが、電子配置のラベルをソートすることによって効率的に行うことができる。今回は、OpenMP を使用して並列ソートを実装することによって高速化した。
effective Hamiltonian の行列要素は、各ノードで並列計算され(Sigma step)、Macro step 毎に MPI によってマスターノードに集められる。対角化によって計算された目的の状態のエネ
ルギーは、全てのノードに Broadcast され、次のステップで使用される。MPI によってノード間で通信されるのは、effective Hamiltonian の行列要素とエネルギーのみである。
【計算結果と考察】 開発した並列プログラムを使用して、C2 分子の基底状態と励起状態ポテンシャルカーブを計算し、Full-CI 法による結果[5]との比較を行った。P-space に含まれる
電子配置は、CASCI 計算を行い、目的の状態の波動関数の中で、係数が閾値（0.1）以上のものとした。必要な配置が P-space に選択されない場合も、モンテカルロシミュレーション中にその配置のウォーカー数が増加し、重要だとわかるので容易に P-space を改善することができる。また、サンプリングの効率を向上させるために、initiator 法[4]も利用した。図 2 に示すよ
うに、1～2.8Å にわたって X1Σg+ , B1∆g , B’1Σg+の 3 つの状態の計算を行ったが、誤差は最大でも 4mHartree 以下であった。誤差の主な原因は、initiator 法によるウォーカーのカットオフ
から来るものである。当日は、アルゴリズムの詳しい説明と他の系への応用例を紹介する予
定である。
X1Σg+

全エネルギー (a.u.)

B1∆g
B’1Σg+

Basis set: 6-31G* (1s core frozen)
Full-CI dimension: 5.24 ×107

C-C 距離 (Å)
図 2. MSQMC 法と Full-CI 法による C2 分子のポテンシャルカーブの比較[1]

S. Ten-no, J. Chem. Phys., 138, 164126, (2013).

[2]

Y. Ohtsuka and S. Nagase, Chem. Phys. Lett., 463, 431, (2008).

[3]

G.H. Booth, A.J.W. Thom and A. Alavi, J. Chem. Phys., 131 , 054106, (2009).

[4]

D. Cleland, G.H. Booth, and A. Alavi, J. Chem. Phys., 132 , 041103, (2010).

[5]

M. L. Abrams and D. Sherrill, J. Chem. Phys., 121, 9211, (2004).

4E17

量子モンテカルロ法による多原子分子の振動励起状態への陽電子吸着機構の理論的解析（横市大院・生命ナノ）北 幸海，立川 仁典

Theoretical investigation of the binding of a positron to vibrational excitedstates of polyatomic molecules with quantum Monte Carlo method(Yokohama City Univeristy) Yukiumi KITA, Masanori TACHIKAWA【 緒 言 】 陽電子 (e+) は電子と同質量・同スピンおよび電荷+1 を持つ電子の反粒子であり、電子との衝突により、２〜３個の光子を放出しながら対消滅を起こす。物質中に入射された陽電子は対消滅を起こす前に、原子・分子のイオン化や励起、電子と陽電子から成る水素様原子であるポジトロニウム形成、そして原子・分子との一時的な束縛状態である陽電子化合物形成など、様々な反応を起こす事が示唆されている[1]。
近年、カルフォルニア大学サンディエゴ校の Surko らは、入射エネルギーを制御した低速
陽電子を用いた対消滅率測定実験を行い、アルカン、アルデヒド、ケトン、ニトリルなど様々な分子の陽電子束縛エネルギー（陽電子親和力, PA）を報告している[2,3]。振動 Feshbach 共鳴を利用した彼らの実験では、分子の振動励起状態に吸着した陽電子の束縛エネルギーを測定していると考えられているが、対応する理論計算はほとんど報告されておらず、その吸着機構等の詳細は十分明らかにはなっていない。
そこで本研究では、多原子分子の振動励起状態への陽電子吸着機構の詳細を明らかにすることを目的に、幾つかの小分子に対して振動励起状態における陽電子親和力の理論的解析を行った。具体的には、量子モンテカルロ (QMC) 法に基づいた分子の非調和振動計算、および第一原理多成分分子軌道 (MC_MO) 法に基づいた陽電子親和力計算を組み合わせること
で、様々な振動状態における振動平均陽電子親和力の解析を行った。対象分子は最も単純なニトリルであるシアン化水素 (HCN) である。
【 方 法 】 本研究で解析した振動平均陽電子親和力は以下で定義される：

PAνA ≡ ∫ PA A (Q) × ψ νA (Q) dQ,
2+

PA A (Q) = E A (Q) − E [A;e ] (Q).

ここで Q は基準振動座標、 ψ νA は分子 A の規格化された振動の波動関数（振動量子数 v）、+

PA A (Q) は座標 Q における分子 A の陽電子親和力、 E A と E [A;e ] は分子 A とその陽電子化合A[A;e + ]物[A;e+]の全エネルギーである。本研究では、 E を Hartree-Fock (HF) 法により、 Eを1
電子、1 陽電子、および 1 電子-1 陽電子励起配置のみを考慮した配置間相互作用法 (CISD) [4]により計算した。
振動の波動関数 ψ νA の解析には、変分モンテカルロ (VMC) 法と Reptation Monte Carlo(RMC) 法[5]という２つの QMC 法を用いた。VMC 法は、変分パラメータを含む試行波動関数を用いて物理量の期待値をモンテカルロ積分によって算出する方法である。本研究では平均場近似に対応する試行波動関数を仮定し、Umrigar らによって提案された励起状態にも適用可能なエネルギー最小化法[6]を用いて、振動基底および励起状態における試行波動関数の最適化を行った。RMC 法は、時間依存 Schrödinger 方程式の虚時間発展を利用した QMC 法の１つであり、任意の初期試行波動関数に虚時間プロパゲータを作用させる事で、系の正確な固有状態へと射影する方法である。一般的に用いられる (importance-sampled) 拡散モンテカル
2ロ(DMC) 法とは異なり、RMC 法では系の正確な確率密度分布 Ψ 0 を直接得ることができる。
【 計 算 の 詳 細 】 HCN 分 子 の 陽 電 子 化 合 物 [HCN;e+] の CISD 計 算 で は 、 電 子 基 底 に
6-31++G(2df,2pd)を、陽電子基底に 15s15p3d2f1g の Gauss 型関数を用いた。この計算レベルでの HCN の平衡構造における PA は 39.6meV となり、DMC 法による高精度計算の値 38(5) meV[7]を定量的に良く再現している。非調和振動解析における HCN 分子の多次元ポテンシャル
エネルギー曲面 (PES) は、Born-Oppenheimer 近似の下、CCSD(T)/ aug-cc-pVTZ レベルで計算したグリッド点をスプライン補完することで生成した（全グリッド数は 14692 点）。
【 結 果 と 考 察 】 Figure 1 に HCN 分子の振動基底状態、および基音・２倍音・結合音
準位における振動平均 PA を示す。各振動
状態は平均場近似の試行波動関数を用いた
VMC 法によって得られたものである。振動
基底状態における PA は、平衡構造における値(39.6meV)とほぼ同一の値となった。一方、基音準位では、変角モード(v1,2)において PA は減少し、対称伸縮モード(v3)と逆対
称伸縮モード(v4)において PA が増大することがわかった。これら PA の増大／減少は、Fig. 1 Vibrational averaged positron affinities (PA) for

各基音準位における双極子モーメントの期 vibrational ground and excited states of HCN molecule obtainedwith VMC calculations. The v1,2 is vibrational quantum number

待値と相関しており、 v1,2, v3, v4 の基音準位 for two equivalent bending modes, v3 the symmetric stretchingmode, and v4 the anti-symmetric stretching mode.
における振動平均双極子モーメントはそれ
ぞれ 3.25D, 3.29D, 3.32D である（振動基底状態では 3.28D）。また、２倍音・結合音準位における PA の変化に関しても、各振動準位における振動平均双極子モーメントと強い相関があることがわかった。全ての振動順位における PA と双極子モーメントに対する線形回帰分析
から、その決定係数は R2=0.977 と大きな値となり、HCN 分子の振動励起状態への陽電子吸着
は、各振動準位における双極子モーメントの寄与が支配的であることがわかった。RMC 法による結果、および他分子の振動平均 PA の解析結果については、当日報告を行う。
【参考文献】[1] 陽電子計測の科学 (日本アイソトープ協会, 1993). [2] J. R. Danielson, J. J. Gosselin, and C. M. Surko,
Phys. Rev. Lett. 104, 233201 (2010). [3] J. R. Danielson, A.C.L. Jones, M.R. Natisin, and C. M. Surko, Phys.
Rev. Lett. 109, 113201 (2012). [4] M. Tachikawa et al., Phys. Chem. Chem. Phys., 13, 2701 (2011). [5] S. Baroni
and S. Moroni, Phys. Rev. Lett., 82, 4745 (1999). [6] J. Toulouse and C. J. Umrigar, J. Chem. Phys. 126, 084102(2007).	

[7] Y. Kita, et al., J. Chem. Phys. 131, 134310 (2009).

4E18
量子モンテカルロ法による 1,4-ジヨードベンゼン分子結晶多形予測	
 （北陸先端科学技術大学院大学*）	
  ○本郷研太*，前園涼*	
 Quantum Monte Carlo study of predicting 1,4-diiodobenzene polymorphism(JAIST*) ○Kenta Hongo*, Ryo Maezono*【緒言】	
 分子結晶はアスピリンなどの薬理物質や、近年では、クリーンエネルギー実現のための有機半導体
機能性物質として注目されている。その特徴は、系の安定性が分子間力(水素結合・分散力)で決まる
ため、比較的小さいエネルギー(数〜数十 kcal/mol)で容易に構相転移を起こし、その立体配列の違い(結晶多形)から、物理化学的性質が大きく変化する点にある[1]。シミュレーションによる正確な結晶
構予測が可能となれば、ドラッグ・デザインやマテリアル・デザインにおける多大な組	
 合せ探索に劇的な効率化がもたらされる。しかしながら、現行の枠組みにおいては、その構造を決定づけている分散力自体の理論記述が不十分なため、この方向での実用化が阻害されている。分散力は高次の電子相
関効果であり、分子間力の中でも特に、電子間相互作用の精緻な取扱いが必要であり、ドラッグデザ
インや物質デザインにおける挑戦課題となっている。	
 	
 本研究で対象とする 1,4-ジヨードベンゼン(para-diiodobenzene; p-DIB)分子結晶は、高効率太陽電池材料の一つとして注目を集め、その特性理解のために、密度汎関数法による先行研究が行われたが、標準的な密度汎関数(DFT)法による理論予測では結晶多形の相対的安定性について、実験結果を説明できないことが報告されている[2]。実験によれば、p-DIB は 326K で α 相から β 相への構造相転移を示し、従って、絶対零度においては、α 相の方が β 相よりも安定であることが知られている(図 1)。他方、CCSD(T)法などの高精度分子軌道法を直接分子結晶系に適用しようとすると、現状の超並列計算機であっても、計算資源不足から計算が走らないといった問題に直面する。そのため、計算精度と計算コストの両面でバランスのとれた代替法が必要となる。計算精度の観点から、第一原理量子モンテカルロ(QMC)法[3]は、分散力の取扱
い信頼性の高い手法として知られており、最近では、生体分子系において分散力を高精度に再現している[4]。高い並列効率(98%超)故に超並列計算機との親和性が高く、近年、実用性を伸ばしつつある。申請者らは	
 最近、前述の p-DIB 結晶多形の問題に本手法を適用し、その結晶多形の理論予測を行い、p-DIB 図	
 1 :	
 p -DI B 分子結晶の２つの 結晶多形：
α相（上図）およ びβ相（下図）	
 分子結晶多形の相対的安定性を正しく再現した[5,6]。	
 	
 当該先行研究で実施した QMC シミュレーションでは、しかしながら、計算機資源の制限から、
1×1×1 シミュレーション・セルでの QMC 計算を行い、それに対して、Kwee らによって提案されたDFT に基づく手法（KZK 補正スキーム[7]）によって有限サイズ効果を見積もっていた。結晶の周期性を再現するためには、本来、異なるシミュレーション・セルサイズで計算を行い、無限大への外挿が必要となるが、QMC 法は、DFT 法と異なり、非常に大きなセルサイズの計算をすることは実際上不可能であり、信頼性の高い外挿値を得ることは難しく、依然として、有限サイズ効果[8]に起因する誤差の見積もりが必要となる。そこで、本研究では、先行研究よりも大きなシミュレーション・セル
サイズでのQMC 計算を行って、有限サイズ効果を見積もる。	
 	
 【計算方法】	
 	
 本研究では、p-DIB の多形である α と β の 2 つの相の相対的安定性を調べる。これら 2 相に対して、
本研究では、QMC 法の中でも、特に分散力の扱いに高い信頼性を持つ、拡散モンテカルロ(DMC)法を適用し、2 相のエネルギー差を評価する。本研究の DMC 計算は CASINO プログラムを用いる[3]。
DMC 計算で考慮するシミュレーション・セルサイズは、1×1×1 と 1×3×3 の 2 つであり、そこでの有限サイズ効果を調べるとともに、また、1×1×1 セルサイズに対しては、節固定近似の影響を調べるた
めに、3 つの試行節を考慮する。本研究で採用した試行関数（試行節）は、スレーター・ジャストロー
型関数である。スレーター部分については、1×1×1 セルサイズに対して、LDA/PBE/B3LYP 汎関数に

よる DFT 計算を行い、1×3×3 セルサイズに対して、PBE 汎関数のみを考慮した。（これ以降、DMC
計算スペックを表す略記号として「DMC/SCF/1×n×n」を用いる。ここで、SCF = LDA/PBE/B3LYP、
n = 1/3 である。）結晶軌道は平面波で展開し、そのカットオフエネルギー値は 40 Hartree とした。
試行関数生成のための DFT 計算は Qunatum Esrepsso を利用したが、1×1×1 セルの LDA 節生成計算
は先行研究と同様に ABINIT を利用した[4]。DMC 計算では、この結晶軌道をスプライン(blip)基底により再度、基底展開している。ジャストロー部分については、CASINO コードに実装されている多項式表式を採用し、合計で 36 個の変分パラメータを含む。これらの変分パラメータは、分散最小化法
によって、最適化した。p-DIB を構成する炭素・ヨウ素・水素の各元素につき、内殻電子部分をTrail-Needs 型擬ポテンシャルを用いて置き換えた。局所近似に基づく擬ポテンシャル評価で生じるバ
イアスの軽減を目的として、T-move スキームを採用した。固体周期系での電子間相互作用の取扱いとしては、Ewald と MPC (Model Periodic Coulumb)	
 相互作用の 2 つを考慮した。ただし、DMC 虚時
間発展における分岐過程では、Ewald 相互作用のみを考慮して評価した、局所エネルギー値を用いている。Ewald と MPC で得られる結果の相違は、有限サイズ効果の大小を判断する上での良い指標となることが知られている[8]。DMC の計算条件は次の通りである：目標ポピュレーション数は、1280(1×1×1)、20480 (1×3×3)である。平衡化ステップ数は 1500 として、統計蓄積ステップ数は、1×105(1×1×1)、および、7×103 (1×3×3)	
 とした。1×1×1 計算は PC クラスター計算機の 32 コアを利用し、1×3×3
計算はスーパーコンピュータ京の2048 コアを利用して、並列計算を行った。	
 	
 【結果と考察】	
 
	
 
表	
 1:	
 DMC による p-DIB の 2 相 α および β のエネルギー差 ΔE (定義については本文を参照のこと)。
エネルギーの単位は、kcal/mol/cell である。各項目中の記号+/-の右側の数値はエラーバーを表す。	
 

	
 
	
 

LDA/1x1x1

PBE/1x1x1

B3LYP/1x1x1

PBE/1x3x3

Ewald-3.8 +/- 1.6-0.5 +/- 1.5-0.3 +/- 1.4-1.8 +/- 1.3

MPC-3.2 +/- 1.6

2.8 +/- 1.5-6.5 +/- 1.4-1.5 +/- 1.3

表 1 に、DMC 計算によって得られた、2 つの結晶多形のエネルギー差	
 ΔE = E(α)- E(β)	
 の評価値を示す。ここで、E(α/β)は α/β の全エネルギーであり、ΔE < 0 は、α 相の方が β 相よりも安定である
ことを示し、すなわち実験結果との一致を意味する。表１から分かるように、1×1×1 セルサイズの場
合、DMC/LDA、DMC/PBE、DMC/B3LYP の 3 つの結果を比較すると、相互作用にかかわらず、結果が固定節に強く依存することが分かる。特に、DMC/PBE と DMC/B3LYP では、Ewald と MPC の結果に大きな差が存在するが、それに対して、DMC/LDA では、Ewald と MPC の差は小さく、両者
はエラーバーの範囲内で一致しいる先行研究では、DMC/LDA/1×1×1 計算を行って、実験と一致する
結果が得られていたが、本研究の結果から、先行研究は相対的安定性の再現に偶々成功したものであると結論づけられる。次に、DMC/PBE/1×3×3 の結果に移ると、Ewald と MPC の結果が非常に良く一致しており、2 相間のエネルギー差に関する限り、有限サイズ効果の影響は極めて小さいと結論づけられる。また、DMC/PBE/1×1×1 では、ΔE の符号を正しく再現できなかったが、シミュレーション・
セルを 1×3×3 まで大きく取ることで、DMC/PBE/1×3×3 では、ΔE < 0 を正しく再現していることが分かった。当該結果を得るために、いくつかの準備計算を行っており、それらについての詳細は当日に発表する。
	 【参考文献】	
 
	
 
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]

Bernstein, J. Polymorphism in Molecular Crystals; Oxford University Press: New York, 2002.
Brillante, A.; Della Valle, R. G.; Farina, L.; Venuti, E.; Cavazzoni, C.; Emerson, A. P. J.; Syassen, K. J.
Am. Chem. Soc. 2005, 127, 3038–43.
R. J. Needs, M. D. Towler, N. D. Drummond, P. Lopez Rvos, J. Phys.: Condens. Matter 2010, 22,023201.
Hongo, K.; Cuong, N. T.; Maezono, R. J. Chem. Theory Comput. 2013, 9, 1081–1086.
Hongo, K.; Watson, M. A.; Sánchez-Carrera, R. S.; Iitaka, T.; Aspuru-Guzik, A. J. Phys. Chem. Lett.
2010, 1, 1789–1794.
Watson, M. A.; Hongo, K.; Iitaka, T.; Aspuru-Guzik, A. Advances in Quantum Monte Carlo; Chapter 10,
pp 101–117.
Kwee,H.;Zhang,S.;Krakauer,H. Phys. Rev. Lett. 2008, 100, 126404.
Drummond, N. D.; Needs, R. J.; Sorouri, A.; Foulkes, W. M. C. Phys. Rev. B 2008, 78, 125106.

4E19

量子拡散としての電子ダイナミックス(豊橋技術科学大学) 関野秀男・浜田信二

Electron Dynamics as Quantum Diffusion(Toyohashi University of Technology) Hideo Sekino, Shinji Hamada【序論】
分子結合に深く関わる電子ダイナミックスは量子力学の法則に従っているが、その基本方
程式である Schrödinger 方程式の解法にも解くべき電子状態の特性により種々の方法論が考案されている。本稿では確率保存と内部自由度の考慮を行うだけで時間依存 Schrödinger 方
程式(TDSE) の帰結をその近似として得られる量子ウォーク法[1]に注目し、電子ダイナミッ
クス数値解法の新しい方法論としての可能性を論じる。
光速よりずっと低い速度を持つ自由電子の挙動は運動エネルギー項のみをハミルトニアンとする TDSE の解として与えられるが、同様の帰結は内部
自由度を考慮した
粒子の自由拡散の長波長限界として与えられる[2,3]。
図１に実ガウシア
ン分布を初期関数として左右等確率で進む量子ウォークによる時間発展を示す。

図１シンメトリック Schrödinger Quantum Walk (30000 時間ステップ)

隣接格子のみへの移動を許すウォークでは格子点の全体は奇数時間ステップと偶数時間ステップの２グループに分かれる。エネルギーの変化分と運動量の変化分の間の２次分散関係がある Schrödinger Quantum Walk では偶数あるいは奇数時間ステップのみを実際のステップとみなしてもよいが、数値算定に使う格子点は実際の物理空間の２倍となる。現実的な高次元空間への QW の拡張のためには内部自由度に対応するパリティーの物理的考察がより重要となる。一般にｍ次元空間において単純な隣接格子へのウォークを考えれば１時間ステップ
ごとに 2m 格子点への移動を考慮することになるため、意味のある格子点の数にくらべ数値算
定に使う格子点の数は非常に大きくなる。

図２ ３次元空間での量子ウォーク（それぞれ、xy 平面（左上）、xz 平面（右上）、yz 平
面（左）への射影）

我々は量子ウォークを一般的な拡散係数を有する拡散と解釈し、その連続極限でSchrödinger 方程式の解を与えることを示した。更に拡散係数が空間依存をする非均一系でのウォークについて考察し[3]、一般的な場での Schrödinger 方程式の数値解への QW 法の適用有効性を示した。
【参考文献】[1]Quantum Potential Theory: Structure and Applications to Physics at the Alfried Krupp
Kolleg Greifdwald, N.Konno, Quantum Walks (Sangyo Tosho 2008) Feb. 26- March 10,(2007)[2] A Solution of Time Dependent Schrödinger Equation by Quantum Walk, Hideo Sekino,
Masayuki Kawahata and Shinji Hamada, Physics Journal of Physics Conference Series,
352 012013, (2012)[3] Quantum walk with variable parameters, Shinji Hamada, Masayuki Kawahata and Hideo
Sekino, Journal of Quantum Information Science (2013)