水素の先端計算による水素機能の高精度解析
常行

真司

することは，原子数が増えるにつれて可能性の数が爆発的に.

は じめに増えるため，非常に困難になる．そこで可能性を絞るために実験データを利用しようというのが，我々の提案するデータ

ユビキタスな元素である水素は，様々な物質中に不純物と同化手法である．

して存在し，半導体の不動態化(注キャリアドープした半導体の電気的活性が失われること)を引き起こす例のように，原子周りの環境に応じて異なる荷電状態をとることで，物性を大きく左右することが知られている．また近年，様々.


水素を含む系の第一原理シミュレーション
電子状態計算な酸化物中の酸素をヒドリドで大量に置換してキャリアドー物質・材料の結晶構造，電子状態，物性のシミュレーショピングできることや，新しいヒドリド伝導体，高温超伝導をン手法のうち，周期系，液体，表面・界面といった凝縮系一示す超高圧下の硫化水素などが見つかり，固体中の水素が改般で経験パラメータを使わないシミュレーションを可能にすめて脚光を浴びている．一方，電極反応や各種触媒反応，合るのが，密度汎関数理論(Density Functional Theory，以下金の水素吸蔵・放出など，水素社会実現のために必要な化学DFT)に基づく第一原理電子状態計算と，そこから計算され反応も，近年さらに重要性が増している．こういった水素をる原子間力を使って原子を動かす第一原理分子動力学法

含む物質・材料の研究では，水素量や水素原子位置，ダイナ(FirstPrinciples Molecular Dynamics，以下 FPMD)である．

ミクスの計測の難しさが，常に大きな問題となっている．

DFT は1960年台に W.

Kohn らが作り上げた理論で，物新学術領域研究「ハイドロジェノミクス高次水素機能に質の基底状態エネルギーが電子密度分布だけで決まり，そのよる革新的材料・デバイス・反応プロセスの創成」の計画研計算に必要なエネルギー汎関数は物質によらないユニバーサ

究 A05 2 「水素の先端計算による水素機能の高精度予測」

ル な も の で あ る と い う こ と 示 し た Hohenberg Kohn の 定

では，この水素の問題を計算科学の立場で解決するため，第理(1)と，電子密度の計算を相互作用しない仮想的な電子系の一原理シミュレーション手法の開発を行っている．中でも，問題にマップすることで，電子密度とエネルギーの精密な計計測とシミュレーションを直接つなぐデータ同化手法の開発

算 手 法 を 与 え た Kohn Sham 理 論(2) か ら な る ． DFT に よは，本領域に特徴的な試みである．

り，波動関数理論ではできなかったような凝縮系の電子状態

本稿では，まず第 2 章で物質・材料の研究を支える第一原理シミュレーション手法の現状を紹介する．とくに水素で計算が可能となり，開発者の Kohn は 1998 年にノーベル化学賞を受賞した(3)．

顕著な原子核の量子効果，有機物のシミュレーションで重要

DFT の実際の計算では，電子密度からエネルギーを計算なファンデルワールス相互作用，水素化合物の超伝導転移温する際に，何らかの近似に基づく交換相関エネルギー汎関数度の定量計算は，ハイドロジェノミクス研究を進める上で重を利用する．有名なのは一様電子ガスの結果を利用する局所要な手法である．第 3 章では，実験的観測が困難な水素を密度近似( LDA )や密度勾配の補正を加えた一般化密度勾配含む物質・材料のデータ同化構造探索手法について述べる．

近似(GGA)であるが(文献(3)およびその引用文献を参照)，計算機シミュレーションだけを使って物質の結晶構造を予測

最近では，波動関数理論であるハートリー・フォック法の交東京大学大学院理学系研究科教授(〒3020124 東京都文京区本郷 731)
Hydrogen Function Analysis by Advanced Computer Simulations; Shinji Tsuneyuki(School of Science, The University of Toyo, Tokyo)
Keywords: hydrogenomics, data assimilation, firstprinciples simulation, density functional theory (DFT), molecular dynamics (MD)2021年 1 月25日受理[doi:10.2320/materia.60.176]

特集換エネルギーの表式を取り入れたハイブリッド汎関数もよく(MD)を使って温度効果も含めてシミュレーションすること用いられている(4)．また局在性の強い d 軌道や f 軌道に，電ができる．この手法を経路積分分子動力学法(PathIntegral

子相関効果を格子モデル的に取り入れる DFT＋ U 法も使わ

Molecular Dynamics，以下 PIMD)と言い，とくに第一原理れている(5)．に基づくものを FPPIMD と言う．FPPIMD は D.

Marx

これらの汎関数では，極性のない分子間の相互作用としてらによって初めて水素を含む分子に適用され(10) ，その後，重要なファンデルワールス( vdW )相互作用の記述が不十分

著者のグループが初めて凝縮系(結晶シリコン中の不純物水で，分子性の固体，液体，表面分子吸着系などの計算結果の素(11) ，超高圧下の固体水素(12) )に適用して，原子核の量子精度が悪い．そこで近年，簡便で正確に vdW 相互作用を計効果を論じた．

算するための手法が研究されている．濱田により提案されたPIMD では，量子系を古典系にマップするにあたって，vdW 密度汎関数は，その中でも成功を収め，広く使われ始計算量が仮想粒子の数の分だけ増える．量子効果が大きい場めた手法である(6)．

合ほど，近似の精度を担保するのに必要な仮想粒子の数を増

動力学計算やす必要があり，水素を含む系の典型的な計算量は通常のMD の100倍のオーダーとなる．そのため，もともと計算量1985 年になると，原子の動きとそれに伴う電子状態変化の多い DFT を使った FPPIMD の応用はなかなか広がらなを同時に解くことで計算量や計算に必要な計算機のメモリーかった．幸い FP PIMD は並列計算しやすいため，近年の量を大幅に減らすことのできる CarParrinello 法が登場し，並列計算機の発展に伴って利用しやすくなり，日本では志賀

液体や表面の動力学シミュレーションも DFT で行えるようらが非常に汎用性の高いプログラムを開発・公開したことになった(7)．この手法を電子状態の逐次解法へと発展させ，で，さまざまな系に応用が広がっている(13)．

より精密な電子状態計算と分子動力学計算を実現したのが，現在の FPMD である(8)．FPMD は動力学計算だけでなく，ポテンシャルエネルギー面をたどる構造最適化にも利用され

超伝導第一原理に基づく超伝導転移温度の計算には様々な報告例があるが，多くの場合， DFT による第一原理電子状態計算る．
近年，電極での水素発生のように，電子のリザーバーであと電子格子相互作用を表す経験パラメータを組み合わせて用る固体表面と外部のイオンの間で電子移動が起きる化学反応いており，予言力には限界があった．近年，ドイツのを， FPMD で調べたいというニーズが高まってきた．その

E.K.U. Gross らが，DFT を拡張して超伝導転移を記述するような電子移動を伴うシミュレーションを有限サイズの系で理論(超伝導密度汎関数理論，以下，SCDFT )を提案し，少行 う と ，電 子 1 個 の 移 動で 表 面 が無 視 でき な い ほど 帯 電なくともフォノンを媒介とする超伝導については，真に非経し，電位が変化するという問題が生じる．大谷と杉野による験的な転移温度予測の道が開かれた(14) ．その後フォノン以有効遮蔽媒質法は，この問題を解決して電子数一定ではなく外，たとえばプラズモン(15)やスピン揺らぎ(16)の効果を取り電位一定条件という開放系のシミュレーョンを実現する手法入れる試みもなされている．

であり， FPMD と組み合わせて広く利用されるようになっている(9)．


原子核の量子効果

固体水素は以前から圧力により金属化すると予想され，金属化すれば高温超伝導体になるであろうと考えられてきた．
残念ながら固体水素の金属化が起きる圧力はあまりにも高く，水素単体での超伝導は実現されていないが，水素を多量FPMD では，与えられた原子位置を使って DFT に基づに含む系で高温超伝導の発現がいくつか実証されている．そく電子状態計算を行い，そこから原子に働く力を計算し，ニのきっかけとなったのは，ドイツのマックスプランク研究所

ュートンの運動方程式を使って原子を動かす．すなわち原子のグループが発見した，硫化水素高圧相での200 K を超えるは古典的な粒子として扱っているので，しばしば用いられる温度での超伝導転移である(17)．このような水素系の超伝導「量子 MD」という言い方は誤解を生みやすく推奨できない．

転移は，基本的にはフォノン媒介であると考えられ，実際これに対して，本稿の主題である水素原子は質量が軽いた

SCDFT を 使う と転 移温 度が よく 再現 され る(18)(19) ．そ のめ，固体中にある時の典型的な零点振動エネルギーが0.1 eV時，水素が軽く量子効果による分布の広がりが大きいことを程度となり，これは温度に換算すると 1000 K 以上の大きな反映し，格子振動の非調和効果が大きいことも指摘されてい値である．つまり室温程度の温度では，分布関数の広がりのる(20)．

主な要因は熱揺らぎではなく量子効果ということになる．このような，エネルギーや分布に現れる量子効果を，我々は通.

データ同化構造探索

常無視することが多いが，本当に無視して良いかどうかは物質や現象ごとによく検討しなければならない．

これまで紹介した第一原理シミュレーション手法では，い量子効果が重要と思われる場合には，ファインマンの経路ずれの場合も初期原子座標を入力として与え，そこからエネ積分法を使うことで，量子力学的な原子をバネでつながれたルギーが下がるように原子を動かしたり，温度に見合った初複数の仮想古典粒子として近似的に記述し，分子動力学法速度を与えて動力学計算をしたりする．では適切な初期原子

ま て り あ 第60巻 第 3 号(2021)
Materia Japan座標がわからない場合，どうすれば良いだろうか．というベイズの定理が成り立つ．有限温度 T でのカノニカ

典型的な例は，化学組成から未知の結晶構造をシミュレーションで探索したいという場合である．そのような場合にル分布を考えると，右辺の r(R)は，ポテンシャルエネルギー E(R)と温度の逆数 b(＝1/kBT)を用いて，r(R) ∝ exp [－bE(R)］は，たとえばランダムな初期原子配置をたくさん用意してそ(2)れぞれ構造最適化を行い，最も安定なものを選ぶ( Randomと書ける．R が決まれば回折ピークの角度や大きさはすべてSampling)，液体のように原子が動き回れる高い温度の MD計算できるので，その計算値を Icalc(u; R)とする．r(I|R)をから始めて徐々に温度を下げて安定構造を見つける( Simu-

考える際には， u の関数である Icalc ( u; R )が I ( u )からずれる

Annealing ，以下 SA )，準安定構造を一つ見つけたにつれてその確からしさが減ると考え，ずれが大きいほど大ら，その近傍でポテンシャルエネルギー障壁を越えるようにきなペナルティーを与えるようなエネルギーの次元をもった原子を動かして次の安定点を見つける(Basin Hopping, Mini-汎関数 D[I, Icalc]を導入して，latedr(I|R) ∝ exp (－bD[I, Icalc])ma Hopping)，といった様々な手法が提案されている．また最近成功を収めている手法には，ランダムな初期配置の構造と書く．これらをベイズの定理の式に代入すると，最適化に，遺伝子のような交叉と突然変異を組み合わせる遺伝アルゴリズム(Genetic(3)r(R|u) ∝ exp {－b(E(R)＋D[I, Icalc])}(4)

Algorithm)や，鳥が群れになってが得られる．構造シミュレーションの各ステップでの R が餌場を探すかのように多次元空間を探索する粒子群最適化法

わかればその構造に対する Icalc が計算され，そこからペナ(Particle Swarm Optimization)がある［レビュー論文としルティー D も決まるので， D [ I, Icalc ]を明示的に D ( R )と書て，文献(21)をあげておく］．これらの結晶構造探索は，シくことにしよう．結果的に我々が欲しい r(R|u)は，仮想的ミュレーションセルに含まれる原子の座標からなる多次元空なポテンシャルエネルギー関数間での最適化問題であり，正しい解を得ることは原子数(次F(R)≡E(R)＋D(R)(5)

元)が増えるにつれて急激に困難になる．現状では 100 原子のもとでの原子分布関数であることがわかる．構造探索のもくらいになると，もはや絶望的であろう．

っとも簡単な手法の一つである SA では， E ( R )で相互作用

そこで我々が注目したのは，データ同化の考え方である．

しながら運動する原子のシミュレーションを融点より高い温気象データは離散的な点でしか得られないため，気象モデル度で開始し，絶対零度まで徐々に温度を下げていくことで最(運動方程式)を解くための初期値を内挿で決めるしかない．

安定構造を探す．回折ピークに関して何らかの情報がある場

また気象モデルも完全なものではない．そのため，時間が経合には，E(R)に替えて F(R)を使った SA を行えば，その情てばシミュレーション結果と現実のずれが大きくなるので，報を利用できるわけである．

時間と共に次々と入ってくる観測データを使って，シミュレこのことをもう少し直感的に示したのが，図である．第ーションを改善しようと考えた．これが気象予報の分野で長一原理シミュレーションで得られる E(R)は，一般に非常に年にわたり様々な手法が開発されてきた，データ同化の基本

複雑な多谷構造をもち，SA などの手法で最安定点を見つけ的なアイデアである．

ようとしても準安定点にトラップされて構造探索に失敗す

結晶構造予測の場合には，時系列の実験データではないる．ペナルティー D(R)はこれとは異なる多谷構造を持つとが，手に入れやすい X 線や中性子線の粉末回折実験データ考えられるが，最安定構造(図 1 の星印)は E (R )と D (R)にを使う．もちろん良質のサンプルがあってきれいな回折実験データが手に入るならば，第一原理に基づく計算機シミュレーションは必要ない．データ同化が役に立つのは，試料の質が悪かったり量が少なかったり，あるいはダイヤモンドアンビルセルを用いた超高圧実験のように実験装置の制限で回折角が限られるなどの理由から，回折実験データだけでは構造が決まらない場合である．水素のように X 線ではほとんど見えない元素の位置を決めたい場合にも有用である．以下では回折実験データを用いたデータ同化構造探索手法について，少し詳しく説明しよう．
全原子の位置をまとめて R ，粉末回折実験データを I ( u )とシンボリックに書くことにする．I は回折角 u の関数である．有限温度での原子の分布関数を r ( R )， I の分布関数をr(I)，実験から I が与えられたときの原子の分布関数を r(R|I)，逆に R がわかったとき I が観測される確率分布を r(I|R)とすると，これらの間にはr(R|I)＝

r(I|R)r(R)r( I )(1)

図1

データ同化構造探索手法の模式図．(オンラインカラー)

特集共通の安定点なので，両者をたしあわせた F(R)ではその安ッケージ LAMMPS(24)を使って初期温度5000 K から 0 K ま定性が強調され，探索が容易になる(22)．

で 5 ps かけて SA を行った．左下図の横軸は結晶化度 l でペナルティー関数は人為的なものなので，様々な選び方がある．SA の時間発展は概ね左から右に進み，この例では最ある．素直に考えれば，観測された回折データ I ( u )と計算終的に l＝1 となって完全結晶にたどりついている．データされた回折データ Icalc ( u; R )の相関係数(を 1 から減じたも同化に使ったのは，回折角 60 度以下に見られる 3 本のピーの)を利用すると良いが，たとえば回折ピークの現れる回折ク位置のみである．初期構造はランダムなので，回折ピーク角だけで定義される結晶化度 l を使うこともできる(22)．

は周期的境界条件が許すすべての角度に現れる．時間が経つ

D(R)＝D[I(u), Icalc(u; R)]＝aN(1－l)

fl＝∑
uobs

uobs＋duobs－df

Icalc(u; R)du/

umax
umin

Icalc(u; R)du

(6)
(7)

式( 6 )の a はエネルギーの次元を持つ定数で，データ同化につれて，データ同化に使った 3 箇所の正しい回折ピーク(×印のついた線)と同じ角度のピークが成長し，正しい結晶構造が得られる．この過程で， E(R)の変化を見ると活性障壁を超えているように見えるが，計算に使った F(R)の方はの重みを決めるコントロールパラメータ， N はシミュレー単調に減少しており，これにより E(R)を使う SA よりも短ションに用いる原子数である．式( 7 )で定義される l の分時間で結晶構造に到達したことがわかる．

母は，実験で計測された角度範囲 umin  u  umax で Icalc ( u 

シリコンは結晶構造が簡単だが，それでも通常の SA でこR)を角度積分した規格化因子，分子は同じく Icalc(u; R)の積れほど急冷してしまうと，結晶構造を得るのは容易ではな分であるが，積分範囲は実験で回折ピークが観測される角度い．対してデータ同化で正しい結晶構造が得られる確率(成uobs の近傍(幅 2d)にのみ限定する．こうすれば常に 0l1功率)は，ほぼ 100 であった．もちろん，このようなデーであり， umin  u  umax に観測されたすべての回折ピークのタ同化がうまくいくかどうかは，コントロールパラメータ a

回折角が(強度を別にして)計算値と一致すれば， l＝ 1(D ＝の大きさ，SA の初期温度やどのくらいゆっくりと温度を下0 )となる．ポテンシャルエネルギーに対するペナルティー

げるかによって変わる．化学組成や凝集機構，構造の複雑の重みを決めるパラメータ a は，何らかの原理により決めさ，実験データの質と量にも依存してしる．したがって，未ることができないが，経験的には融点の熱エネルギーの 10

知物質の構造探索を行う場合，条件設定のための試行錯誤は倍程度でうまくいくようである．

避けられない．しかしながら我々がこれまで様々なタイプの実際にこのデータ同化手法を用いてシリコンの構造探索を結晶(酸化物，水酸化物，水素化物，合金，金属間化合物，行った時の構造およびそこから計算される X 線回折データイオン伝導体など)でテストした限りでは，いずれの場合もの 変 化 の 様 子 と ， そ の 過 程 で 計 算 さ れ た E ( R ) (破 線 ) とデータ同化によって構造探索の成功率が大きく向上すること

F(R)(実線)を，図に示す．この計算は原子数64個の周期が確認された．たとえば高温超伝導を示す硫化水素の高圧相セルと原子間力モデル(23)を用い，古典 MD のプログラムパの場合， X 線回折実験で明瞭に観測されたピークは 4 本し

図2

シリコンのデータ同化構造探索シミュレーションで得られた原子配列の時間変化と，そこから計算される角度範囲 0＜2u＜60°の粉末 X 線回折データ(赤)．データ同化に用いた既知の回折データを，×印つきの線(黒)で表示した．結晶化度 l は時間と共に増加し，完全結晶の値である 1 に近づく．左下図では，その過程で得られた E(R)(破線)と F(R)(実線)を，結晶化度 l に対してプロットした．(オンラインカラー)

ま て り あ 第60巻 第 3 号(2021)
Materia Japanかないが(25)，この情報があるだけで実験を非常に良く説明し，かつ超伝導転移温度の計算値もよく一致する結晶構造を，容易に見つけることができた．多少なりとも回折実験データがあるならば，データ同化に利用してみる価値は大いにあると言って良いだろう．
実は本手法と関連する手法は，過去にもいくつか提案されている．粉末回折実験データから結晶構造を決めるという意
味で最も近いのは， Putz らによる Combined method である(26)．彼らは理論的に結晶構造を探査するというより，粉末回折実験データから逆問題を解いて結晶構造を決める際に，粗い原子間相互作用モデルを使って探索範囲を絞り，構造探索を加速することを目的としていた．我々の手法では，第一原理計算等で得られる信頼性の高いポテンシャルエネルギーを利用し，そこに粉末回折実験データを同化することで，情報不足のために構造決定に役立たないとして捨てられていた実験データでさえ有効利用し，シミュレーションによる構造探索を加速する．つまりデータ同化により，理論と実験のコラボレーションが可能になるのである．

.

お わりに

本稿では，まず第 2 章で，第一原理シミュレーション手法の現状を概観した．その中で述べた vdW 相互作用の精密計算，電極反応のような開放系を取り扱う手法，原子核の量子効果を取り扱う手法は，ハイドロジェノミクス研究にとって重要であり，それぞれの開発者が科研費新学術領域研究「ハイドロジェノミクス」に参加して，領域内の実験グループと様々な共同研究を行っている．
また第 3 章では，粉末回折実験データを利用して第一原理シミュレーションによる結晶構造探索を加速する，データ同化手法を紹介した．この手法は SA 以外の手法と組み合わせたり，E(R)と D(R)を足し合わせずに利用したりと(27)，さまざまなバリエーションが考えられ，いまだ発展途上であるが，すでに実験で見えない水素を計算機で「見る」ための手法として，上記新学術領域研究で未知構造の探索に使われ始めている．現在は粉末回折実験のデータを利用している( 2 ) W. Kohn and L. J. Sham: Phys. Rev., 140(1965), A1133
A1138.
( 3 ) W. Kohn: Rev. Mod. Phys., 71(1999), 12531266.
( 4 ) M. Marsman, J. Paier, A. Stroppa and G. Kresse: J. Phys.
Cond. Mat., 20(2008), 064201.
( 5 ) B. Himmetoglu, A. Floris, S. de Gironcoli and M. Cococcioni:Int. J. Quantum Chem., 114(2014), 1449.
( 6 ) I. Hamada: Phys. Rev. B, 89(2014), 121103(R), Phys. Rev. B,91(2015), 119902 (Erratum).
( 7 ) R. Car and M. Parrinello: Phys. Rev. Lett., 55(1985), 2471
2474.
( 8 ) M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D.
Joannopoulos: Rev. Mod. Phys., 64(1992), 10451096.
( 9 ) M. Otani and O. Sugino: Phys. Rev. B, 73(2006), 115407.
(10) D. Marx and M. Parrinello: Z. Phys. B, 95(1994), 143144.
(11) T. Miyake, T. Ogitsu and S. Tsuneyuki: Phys. Rev. Lett., 81(1998), 18731876.
(12) H. Kitamura, S. Tsuneyuki, T. Ogitsu and T. Miyake: Nature,404(2000), 259262.
(13) M. Shiga: PIMD version 2.4.0(2020), M. Shiga, M. Tachikawa
and S. Miura: J. Chem. Phys., 115(2001), 91499159, M.
Shiga, M. Tachikawa and S. Miura: Chem. Phys. Lett., 332(2000), 396
402.
(14) L. N. Oliveira, E. K. U. Gross and W. Kohn: Phys. Rev. Lett.,
60(1988), 24302433.
(15) R. Akashi and R. Arita: Phys. Rev. Lett., 111(2013), 057006.
(16) K. Tsutsumi, Y. Hizume, M. Kawamura, R. Akashi and S.
Tsuneyuki: Phys. Rev. B, 102(2020), 214515.
(17) A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov
and S. I. Shylin: Nature, 525(2015), 7376.
(18) R. Akashi, M. Kawamura, S. Tsuneyuki, Y. Nomura and R.
Arita: Phys. Rev. B, 91(2015), 224513.
(19) R. Akashi, W. Sano, R. Arita and S. Tsuneyuki: Phys. Rev.
Lett., 117(2016), 075503.
(20) I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y.
Li, H. Liu, Y. Zhang, Y. Ma and F. Mauri: Phys. Rev. Lett.,
114(2015), 157004.
(21) Y. Wang and Y. Ma: J. Chem. Phys., 140(2014), 040901.
(22) N. Tsujimoto, D. Adachi, R. Akashi, S. Todo and S.
Tsuneyuki: Phys. Rev. Materials, 2(2018), 053801.
(23) T. Kumagai, S. Izumi, S. Hara and S. Sakai: Comput. Mater.
Sci., 39(2007), 457464.
(24) S. Plimpton: J. Comp. Phys., 117(1995), 1
19.
(25) M. Einaga, M. Sakata, T. Ishikawa, K. Shimizu, M. I.
Eremets, A. P. Drozdov, I. A. Troyan, N. Hirao and Y. Ohishi:
Nature Physics, 12(2016), 835838.
(26) H. Putz, J. C. Sch äon and M. Jansen: J. Appl. Crystallogr., 32(1999), 864
870.
(27) D. Adachi, N. Tsujimoto, R. Akashi, S. Todo and S.
Tsuneyuki: Comp. Phys. Comm., 241(2019), 9297.

が，透過型電子顕微鏡や走査プローブ顕微鏡の画像など，原★★★★★★★★★★★★★★★★★★★★★

子位置から容易にシミュレーションできる実験データであれ

1986年 東京大学大学院理学研究科修士課程修了
1991年 東京大学物性研教授
2002年 東京大学大学院理学系研究科助教授
2007年 8 月現職
専門分野物性理論，計算物質科学◎第一原理計算を用いた物性研究やその手法開発に従事．

ば，同じようにデータ同化できると考えられる．

文献( 1 ) P. Hohenberg and W. Kohn: Phys. Rev., 136(1964), B864
B871.



★★★★★★★★★★★★★★★★★★★★★
常行真司

特集