ಢಀϤѐࡌϤ‫״‬ట೎࠙մÕࢗѐߵะ͆‫״‬ԚÖǖȖǷȈ»ǩǆǒǜƛƸƻಔ‫އ‬ǪǙǇȖ¦಍‫ړ‬আҮƘడ಩‫ځ‬ǣǇǭȂǒǜ
಩࣬22୿૭࣬Ϧೝ‫ࠕڜ‬

಩࣬ 23 ୿ 4 ‫מ‬

ะ͆য়౉

λ

‫ܛ‬

ࠃ

É଀Բৡѐৡѐٕ͝ѐֿ‫״‬ԚϤಔกٕѐप٢¥Շ޾Êはしがき文部科学省科学研究費補助金「新学術領域研究」として、「コンピューティクスによる物質デザイン：複合相関と非平衡ダイナミクス」が平成 22 年 7 月より 5 ヶ年の計画
でスタートした（http://computics-material.jp/）。コンピューティクス（computics）という聞き慣れない言葉を訝しむ方々もいらっしゃるかもしれないが、これはわれわれの造語である。理論物理と実験物理は物理学の発展を支える両輪であった。1980 年代より第三のアプローチとして計算物理（computational physics）が出現した。従来からの理論的あるいは実験的アプローチではとらえきれない、自然現象への有力なアプローチとなり、いまや物理学あるいは化学の範疇にとどまらず、自然科学全般さらには工学等の様々な分野において、その有用性が認識されている。しかし、これまでの計算科学（computational science）における研究スタンスは、自然をつかさどる基本方程式をコン
ピュータの助けを借りて解くというものであり、コンピュータそのものはブラックボッ
クスであり、情報科学、計算機科学（computer science）との相互作用は稀であった。
コンピュータ・アーキテクチャは近年劇的に変化している。半導体デバイスの微細化限界、スケーリング則の破綻に伴い、科学の最前線を切り開くスーパー・コンピュータ
は、従来のベクトル・アーキテクチャから超並列マルチコア・アーキテクチャへと変化している。さらに次々世代スーパー・コンピュータでは、基本的計算部分をそれ専用のハード加速器で実行することが視野に入ってきている。
こうした状況で、計算（computation）という第三のアプローチを活用して、先端的科学研究を行うには、コンピュータのアーキテクチャを意識し、そのアーキテクチャに最適の理論手法、数学的アルゴリズムを編み出し、物質科学の最前線にチャレンジすることが不可欠である。そのためには、計算機科学と計算科学という、異なる二つの分野の研究者による、共通の場での共同研究がもっとも有効であろう。われわれはこれをコン
ピューティクスという言葉で表現している。ニュートンのプリンキピア以来、数学（Mathematics）は自然科学を推進する強固な枠組みであったが、今世紀はそれに加えてコンピューティクスという推進エンジンが手に入りつつあると我々は考えている。
そのコンピューティクスによるアプローチで、われわれが狙っている物質科学のター
ゲットは、「複合相関と非平衡ダイナミクス」である。ナノ世界の到来に伴い前世紀に培われた常識は破綻しつつある。そこでは、元素の特質に加え、ナノスケールの形状が電子状態に大きな影響を与え、バルク物質では封印されていた新機能が出現している。
共有結合性、イオン性、電子相関等の競合する因子が、ナノ形状と絡み合い、新たな複合相関現象を生み出している。これら複合因子の相関を、量子論の第一原理に基づき非経験的に解明し、新たな物質設計の強固な方法論的基盤を形成することが第一のターゲットである。また、アト秒スケールの電子励起がピコ秒スケールの原子移動を引き起こし、新たな物質相が形成される現象も多々見られている。そうしたマルチ時間スケールでの複合相関の非平衡ダイナミクスを明らかにし、物質生成の機構を解明・制御することは、時間軸を制御した新たな物質機能の発現には不可欠の要素である。これが第二のターゲットである。
本冊子は、今年度にスタートした１１の計画研究の成果報告書である。研究項目 A01「計算機アーキテクチャと高速計算アルゴリズム」には、稲葉、高橋、張の３つの研究班、A02「密度汎関数法の新展開」には、押山、常行、渡邊、中西、倭の５つの研究班、
A03「密度汎関数法を越えて」には、今田、高田、佐藤の３つの研究班が構成されている。本新学術領域が開始されてまだ間もないが、すでにいくつかの成果が得られている。
ご高覧、ご批判賜れば幸いである。
平成２３年４月
領域代表
押山

淳

目

次

はしがき

研究項目 A01「計算機アーキテクチャと高速計算アルゴリズム」
超高速・超低消費電力物質科学シミュレーション方式の研究開発 ........................

1

大規模並列環境における数値計算アルゴリズム ........................................................

9

計算物質科学の基盤となる超大規模系のための高速解法 ........................................ 15

研究項目 A02「密度汎関数法の新展開」
ナノ構造形成．新機能発現における電子論ダイナミクス ........................................ 21
第一原理分子動力学法による構造サンプリングと非平衡ダイナミクス ................ 28
密度汎関数法理論に基づく非平衡ナノスケール電気伝導ダイナミクス ................ 36
プロトン・ミューオンで探る新物性と量子ダイナミクス ........................................ 42
多自由度・大規模系における反応と構造空間探索 .................................................... 58

研究項目 A02「密度汎関数法の新展開」
第一原理有効模型と相関科学のフロンティア .............................................................. 62
第一原理からの多体理論 ................................................................................................ 67
スピンエレクトロニクス材料の探索 ............................................................................ 73

超高速・超低消費電力物質科学シミュレーション方式の研究開発
High-speed ultra low-power simulation methodology for material science
稲葉真理 1、今井浩 1、須田礼二 1
M. Inaba, H. Imai, R. Suda
東京大学 1The University of Tokyo

平成 22 年度は、研究の第一年度として研究チームを形成するとともに、直接的にハードウェアに問題を写像することにより、超低消費電力かつ超高速に物質科学における第一原理に基づいたシミュレーションを実施するための基盤となる研究を実施した。具体的には、超低消費電力シミュレーションを実現するためのベースとなるアーキテクチャ
技術と消費電力の関係性の研究、FPGA ベース演算アクセラレータの研究、極超高速ネットワークを有効に活用するためのデータ転送方式の研究開発、将来の物質科学シミュ
レーションを容易かつ高速に記述するための高生産性言語（HPC Ruby）の研究、GPGPU
実行環境の研究開発を実施した。
これらの研究開発において、本研究の A02 班、A03 班との有機的連携を実現するた
め、1 月に 1 回以上共同セミナーを実施し、計算システム、数値アルゴリズムと物質科
学シミュレーション間の密な連携を実現した。
１． はじめに「超高速・超低消費電力物質科学シミュレーション方式の研究開発」では、計算流体
力学、有限要素法など、すでに確立した分野と比較して著しく困難である、物質科学における第一原理に基づいたシミュレーションの高速化、大規模化を実現することを大局的研究目的とする。これらの困難の最大原因は、物質科学シミュレーションが FFT、密行列計算、疎行列計算、多体相互作用計算など多くの計算要素を複雑に組み合わせて実現していることにある。この性質から、従来用いられてきたベクトル演算器や SIMD ア
クセラレータでは計算全体の加速が困難であった。
本研究開発では、将来の Exa Flops, Zetta Flops スケールの物質科学シミュレーションを実現するための基礎技術として、汎用性を保ったまま Intel 等の汎用プロセッサ・ク
ラスタを用いたシミュレーションより２桁以上の演算速度当たりの消費電力、設置面積とコストの削減すること、高生産性言語を数値シミュレーションで実用に耐えるレベルまで高速化することにより、物質シミュレーションのプログラム記述を高生産性言語で実現することを具体的目標とする。
本研究項目では、上記目標を、物理レベルから直接的にハードウェアに写像し、オーバーヘッドを極限まで低下させた計算機構を求め、それを実現するためのソフトウェア
層を構築すること、また物質シミュレーションを容易かつ高性能にプログラムするための高生産性プログラミング言語処理系を確立することを目的とする。
平成 22 年度は、研究の第一年度として研究チームを形成するとともに、直接的にハードウェアに問題を写像することにより、超低消費電力かつ超高速に物質科学における第一原理に基づいたシミュレーションを実施するための基盤となる研究を実施した。

－1―

２．超低消費電力アーキテクチャ技術の研究
超低消費電力シミュレーションを実現するためのベースとなるアーキテクチャ技術と消費電力の関係性の研究を実施した。
将来の超高性能スーパーコンピュータのアーキテクチャの候補として、（１）高性能
汎用プロセッサ（例えば Intel 社製 x86 プロセッサ）のマルチコア化、（２）低消費電力
汎用プロセッサ（例えば ARM, Intel ATOM）のマルチコア化、および（３）数値計算用
アクセラレータ（GRAPE-DR や GPGPU のようなもの）の併用による応用指向プロセッ
サが想定される。本研究細目では、過去から現在までの多数のプロセッサを同一の標準的ベンチマークプログラムで計測し、プログラム実行中の電力を計測した。図 1 に整数
ベンチマーク（Dhrystone）,図 2 に整数実行速度に対する消費電力を示す。

図１．Dhrystone MIPS の遷移

スペースの関係で詳細データを省略するが、測定の結果、2018 年に Exa Flops を実現するためには 10 倍以上の演算性能加速が必要であり、応用指向プロセッサ技術が必須であることが示された。また、現在最も標準的である SPEC CPU2006 ベンチマークプ
ログラムにおいて、SPEC CPU2006INT（整数ベンチマーク）と Dhrystone ベンチマーク
プログラム、SPEC CPU2006FLT（数値計算ベンチマーク）と NAS Parallel ベンチマーク
プログラムの性能は、ベクトル計算機を含め、長期間、広い範囲のアーキテクチャに対して高い相関があることを示した。SPEC ベンチマーク実行には、非常に完成し、品質が高いコンパイラが必要であり、x86, Power, Alpha, MIPS など広く使われたプロセッサ－2―

以外では計測が困難であることを考慮すると、プロセッサ性能の測定における
Dhrystone, NAS Parallel ベンチマークプログラムの有効性が示された。

図 2．Dhrystone MIPS あたり消費電力の遷移

３．FPGA ベース演算アクセラレータの研究
FPGA ベース演算アクセラレータの研究では、先行研究である汎用 FPGA ベース演算
アクセラレータ、Master-2 を活用して、様々な研究開発を実施した。平成 22 年度は、
超高速ネットワーク上を通信する多用なトラフィックの測定および成分分析の研究を実施し、平成 23 年度に試作を予定している FPGA ベース演算アクセラレータ設計の基盤を構築するとともに、超高速ネットワークを有効に活用するために必須である、バースト性を持つネットワークトラフィックに対するバンド幅予測技術を確立した。
４．極超高速ネットワークを有効に活用するためのデータ転送方式の研究開発
超高速ネットワーク、特に多くのユーザが共用する Shared Network を高効率でデータ転送するためのトラフィック予測方式およびデータ転送プロトコルの研究開発を実施した。また、今後の極高速ネットワークを用いるデータ転送方式の研究のための基盤として、40Gbps ネットワーク実験環境を構築した。
５．高生産性言語（HPC Ruby）の研究
これまで、プログラムを書きやすいが、実行速度が非常に遅いために数値シミュレーシ
ョンでは用いられなかった高生産言語、Ruby の高速化を目的として静的データフロー
解析に基づく Ruby の最適化の研究開発を実施した。

－3―

開発中の Ruby 最適化コンパイラでは、・クラス再定義・メソッド再定義等に対する投機的な検査命令の挿入
・データフロー解析で変数・メソッド定義使用グラフ(VDUG,MDUG)構築
・副作用の伝播を考慮した抽象解釈アルゴリズム
・型解析によるメソッドの静的バインディング・インライン化
・整数範囲解析による多倍長⇔固定長変換コスト除去
・生存解析によるメモリアロケーション除去
・反復的な部分実行アルゴリズム
・クラス定義・メソッド定義のアップデート ⇔ 解析・部分実行・クラス・メソッドの再定義に対応し、Rails 等を静的に解析可能化
・部分実行不可能命令により発生する副作用を解析・分類
・不要な投機的検査の除去を組み込み、著しい高速化を実現した。
６．GPGPU 実行環境の研究開発
GPGPU実行環境の研究開発では、自動並列化トランスレータによるCからCUDAへのプロ
グラム変換の研究開発およびGPGPUのための最適スケジューリングアルゴリズムに取り組んだ。
CからCUDAへの自動並列化トランスレータの研究開発では、Cで記述されたプログラム
を、指示行なしでCUDAのホストコードとデバイスコードに分割し、変換する。変換の処理は以下の手順で実現される：
１．どこをGPUで処理するか決める
２．配列の添字を解析し、アクセスされる範囲を決定する
３．GPU側のメモリ確保、データ転送、メモリ解放コードを生成
４．計算部分をCUDAのカーネル関数（GPUで実行されるコード）に変換
５．ホストコードを、GPUコードを呼び出すように変換

トランスレータで変換したコードをGPGPU搭載したクラスタで性能評価した。被評価
コードは、N-body, CG法、PDE, K-means, 画像フィルタである。評価の結果、N-bodyでは、問題サイズが100を超えると自動変換したGPGPU使用コードの方が高速となることが示された。今後は、変換後コードの更なる高速化を追求する予定である。
GPGPUにおける最適スケジューリングアルゴリズムの研究では、データ転送とカーネル実行のスケジューリングをヒューリスティックにスケジュールする。具体的には、
１．カーネルの実行性能を測定
２．CPU⇒GPU転送とGPU計算を並列化
３．適切な時間だけ、GPU計算のみ実行
４．GPU計算とGPU⇒CPU転送を並列化という手順でGPUスケジューリングを実現した。本方式を8192×256サイズの行列積に適用して評価を行った。図３は評価結果である。

－4―

図３．GPGPUスケジューリングの評価（行列積）

７．連携研究者・研究協力者
連携研究者：
平木敬（東京大学大学院情報理工学系研究 教授）

８．本研究課題における平成２２年度の発表論文等

発表論文(国際会議)：
1) Tomohiro Sonobe, Mary Inaba, “SAT CNF Encoding with Multi-modeling”, The 3rd Annual
Meeting of the Asian Association for Algorithms and Computation (AAAC2010), pp.45,
Pohang, Korea, Apr.17-19, 2010
2) Yoshiki Iguchi, Naoki TANIDA, Kenichi Koizumi, Mary Inaba, Kei Hiraki, “USADAFOX:
Ultra-High-Speed file-Acquisition-system over Distance with Apache and fireFOX”,
TERENA Networking Conference 2010 (TNC2010), Web-page, Vilnius, Lithuania, May.31Jun.3, 2010
3) Yasuo Ishii, Mary Inaba, Kei Hiraki, “Cache Replacement Policy Using Map-based Adaptive
Insertion”, JWAC-1: 1st JILP Workshop on Computer Architecture Competitions: cachereplacement championship, pp.41-44, Saint-Malo, France, Jun.20, 2010
4) Ryo Wakuta and Tomohiro Sonobe, “SPGF - Search Places by Geographical Features allaround the world”, the 10th Annual International Symposium on Applications and the Internet(SAINT2010), pp.133-136, Seoul, Korea, Jun.19-23, 2010
5) Naoki Tanida, Mary Inaba, Kei Hiraki, “Adaptive auto-tuning of TCP pacing”, The Fifth
International Workshop on Automatic Performance Tuning (iWAPT2010), California, USA,
Jun.22, 2010
6) Hisanobu Tomari, Mary Inaba, Kei Hiraki, “Compressing Floating-Point Number Stream for
Numerical Applications”, The First International Conference on Networking and Computing,
Hiroshima, Japan, Nov.17-19, 20107) Seongsoo Moon, Mary Inaba, “Modeling and simulation for EMS”, 3rd International
Conference on Wireless Information Networks & Business Information System (WINBIS'11),
pp.94-100, Kathmandu, Nepal, Feb.27-Mar.1, 20118) N. Fu, H. Imai and S. Moriyama; Voronoi Diagrams on Periodic Graphs. 2010 International
Symposium on Voronoi Diagrams in Science and Engineering (ISVD '), pp.189-198, Quebec,
2010

－5―

9) T. Tanuma, H. Imai and S. Moriyama: Revisiting Hyperbolic Voronoi Diagrams from
Theoretical, Applied and Generalized Viewpoints. 2010 International Symposium on Voronoi
Diagrams in Science and Engineering (ISVD '10), pp.23-32, Quebec, 201010) Kamil Rocki, Reiji Suda, "MPI-GPU Monte Carlo Tree Search", IEEE 2011 International
Conference on Information and Computer Applications(ICICA), Dubai, UAE, Mar.18-20,
201111) Kamil Rocki, Reiji Suda, “Improving the GPU parallel Monte Carlo Tree Search performanceby the standard deviation based error estimation”, 2011 3rd International Conference on
Machine Learning and Computing (ICMLC 2011) Singapore, February 26-28, 201112) Kamil Rocki and Reiji Suda, "Massively Parallel Monte Carlo Tree Search", Proceedings of
VECPAR 2010, 9th International Meeting High Performance Computing for Computational
Science, June 23-25 (25), 201013) Reiji Suda, "Methods of Parallel Experimental Design of Online Automatic Tuning and their
Application to Parallel Sparse Matrix Data Structure", Proc. iWAPT 2010, June 22, 2010

書籍：
1) Naoki Tanida, Kei Hiraki, Mary Inaba, “Highly efficient data transmission facility throughlong distance high-speed networks with background traffic”, Fusion Engineering and Design(FED), vol.85, pp.553-556, Apr. 2010
2) Y. Nagayama, M. Emoto, Y. Kozaki, H. Nakanishi, S. Sudo, T. Yamamoto, K. Hiraki, S.
Urushidani, “A Proposal for the ITER remote Participation system in Japan”, Fusion
Engineering and Design (FED), vol.85, pp.535-539, Jun. 2010
3) Yasuo Ishii, Mary Inaba, Kei Hiraki, “Access Map Pattern Matching for High Performance
Data Cache Prefetch”, The Journal of Instruction Level Parallelism, Vol.13, (ISSN
1942-9525), pp.1-24, Jan. 2011
4) Editors: Ken Naono, Keita Teranishi, John Cavazos and Reiji Suda, "Software Automatic
Tuning: From Concepts to the State-of-the-Art Results", Springer, Sep. 2010.
5) R. Suda, K. Naono, K. Teranishi and J. Cavazos, "Software Automatic Tuning: Concepts andState-of-the-Art Results", Chapter 1 of "Software Automatic Tuning: From Concepts to theState-of-the-Art Results", Springer, Sep. 2010., Editors: Ken Naono, Keita Teranishi, John
Cavazos and Reiji Suda, pp.3-18.
6) R. Suda, "A Bayesian Method of Online Automatic Tuning", Chapter 16 of "Software
Automatic Tuning: From Concepts to the State-of-the-Art Results", Springer, Sep. 2010,
Editors: Ken Naono, Keita Teranishi, John Cavazos and Reiji Suda, pp.275-294.
7) Y. Sawa and R. Suda, "Autotuning Method for Deciding Block Size Parameters in
Dynamically Load-balanced BLAS", Chapter 3 of "Software Automatic Tuning: From
Concepts to the State-of-the-Art Results", Springer, Sep. 2010., Editors: Ken Naono, Keita
Teranishi, John Cavazos and Reiji Suda, pp.33-48.
8) T. Kajiyama, A. Nukada, R. Suda, H. Hasegawa, and A. Nishida, "Toward Automatic
Performance Tuning for Numerical Simulations in the SILC Matrix Computation
Framework", Chapter 11 of "Software Automatic Tuning: From Concepts to theState-of-the-Art Results", Springer, Sep. 2010., Editors: Ken Naono, Keita Teranishi, John
Cavazos and Reiji Suda, pp.175-192

－6―

その他：
1) 入江 英嗣, 本城 剛毅, 平木 敬, 「動的推定によるプリフェッチ量最適化」, 情報処
理学会論文誌コンピューティングシステム, vol.3, no.3, pp.56-66, Sep. 20102) 小川周吾，入江英嗣，平木敬, 「アクセス履歴の不要なマルチコアCPU向け共有キャ
ッシュ配分方式」, 先進的計算基盤システムシンポジウムSACSIS2010, pp.267-276, 奈
良県新公会堂, May.27-28, 20103) 入江英嗣, 本城剛毅, 平木敬, 「動的推定によるプリフェッチ量最適化」, 先進的計算
基盤システムシンポジウムSACSIS2010, pp.277-285, 奈良県新公会堂, May.27-28, 20104) 山本孝志，中西秀哉，長山好夫, 平木敬，稲葉真理，谷田直輝, 「広帯域ネットワー
クを用いた日仏間高速データ転送実験」, 第8回核融合エネルギー連合講演会、(ポスター), 高山市民文化会館 , Jun.10-11, 20105) 谷田直輝, 稲葉真理, 平木敬, 「パケット喪失履歴に基づいたTCP輻輳制御方式」, 2010
年並列／分散／協調処理に関する『金沢』サマー・ワークショップ(SWoPP金沢2010),
電子情報通信学会技術研究報告, コンピュータシステム(CPSY), vol.110, no.167,
pp.25-30, 金沢市文化ホール, Aug.3-5, 20106) 柴田有, 平木敬, 「オフライン環境における多様性の高い実行時自己改変ソフトウェ
ア」, 2010年並列／分散／協調処理に関する『金沢』サマー・ワークショップ(SWoPP
金沢2010), 電子情報通信学会技術研究報告, コンピュータシステム(CPSY), vol.110,
no.167, pp.85-90, 金沢市文化ホール, Aug.3-5, 20107) 石井 康雄・稲葉 真理・平木 敬, 「マップ型履歴を用いたプリフェッチ方式とキャッシュ置換方式の協調動作」, 2010年並列／分散／協調処理に関する『金沢』サマー・
ワークショップ(SWoPP金沢2010), 情報処理学会研究報告, 計算機アーキテクチャ(ARC), vol.190, no.13, PDF, 金沢市文化ホール, Aug.3-5, 20108) 小川周吾，平木敬, 「メニーコアプロセッサ向き共有キャッシュ配分方式」, 2010年並
列／分散／協調処理に関する『金沢』サマー・ワークショップ(SWoPP金沢2010), 情
報処理学会研究報告, 計算機アーキテクチャ(ARC), vol.190, no.13, PDF, 金沢市文化ホール, Aug.3-5, 20109) 小泉 賢一，稲葉 真理、平木 敬, 「FPGAを用いた広帯域高遅延ネットワーク向けの利用可能帯域推定」, 第9回情報科学技術フォーラム(FIT2010), RL-001, 第四分冊pp.1-5,
九州大学, Sep.7-9, 201010) 山田大地，中野雄介，長谷部雅彦，和久田龍，薗部知大，手塚宏史，稲葉真理, 「衛星データ探索フレームワークのためのData-Process Memoization cache の提案」, 第3回
楽天研究開発シンポジウム，（ポスター）, 楽天本社, Dec.18, 201011) 長谷部雅彦, 山田大地, 中野雄介, 和久田龍, 薗部知大, 手塚宏史, 稲葉真理, 「衛星データ探索フレームワーク GGSⅡ の提案と実装」, 第52回プログラミング・シンポジ
ウム, ニューウェルシティ湯河原, pp.75-82, , Jan.7-9, 201112) 泊久信, 平木敬, “Retrospective Study of Performance and Power Consumption of Computer
Systems”, 先進的計算基盤システムシンポジウムSACSIS2011, to appear, 秋葉原コンベ
ンションホール, May.25-27, 201113) 今井浩, 夫紀恵, 森山園子: 量子グラフ状態による情報処理. 「ナノ量子情報エレクト
ロニクスの進展」シンポジウム, 東京大学, Dec.22, 201014) 今井浩: Feynmannの2つの提唱からグラフと計算量への展開. 電子情報通信学会 第23回量子情報技術研究会(QIT 23), 東京大学, Nov.15-16, 2010－7―15) 今井浩: 幾何から量子そしてその先. 「計算機科学の発展と未来」研究集会, 京都大学,
Mar.11, 201116) 本間咲来，須田礼仁, 「GPGPUにおけるデータ転送とカーネル実行のヒューリスティ
ックスケジューリング」、情報処理学会研究報告 HPC-129-22, 201117) 須田礼仁, 「自動チューニング数理基盤ライブラリ ATMathCoreLib」, 情報処理学会
研究報告 HPC-129-14, 201118) 須田礼仁, 「自動チューニングの数理モデルと最適化」, 第22回 RAMP シンポジウム
論文集, 名古屋大学シンポジオンホール、Oct. 28-29 (28), 2010, pp.63-77. 社団法人日本
オペレーションズ・リサーチ学会常設研究部会数理計画
19) 須田礼仁,「自動チューニングのための数理基盤技術」, 応用数理, Vol.20, No.3, pp.5-14,
Sep. 2010.
20) 須田礼仁,「自動チューニングのための実験計画におけるランダマイズの効用」, 日本
応用数理学会 2010年度年会 講演予稿集 pp.295-296. 明治大学駿河台キャンパス，Sep.
6-9 (9th), 2010.

－8―大規模並列環境における数値計算アルゴリズム
Numerical Computation Algorithm on Large-Scale Parallel Environment
高橋大介 1，今村俊幸 2，多田野寛人 1
D. Takahashi, T. Imamura, H. Tadano
筑波大学 1，電気通信大学 2
University of Tsukuba 1, The University of Electro-Communications 2

ペタフロップスを超える性能を持つ大規模並列環境における数値計算アルゴリズ
ムとして、高速フーリエ変換（FFT）、GPU による 4 倍・8 倍精度 BLAS、GPU における疎行列ベクトル積、CPU+GPU 環境上で動作する実対称固有値ソルバ、そして Block Krylov アルゴリズムによる連立一次方程式の求解高速化について研究を行うと共に、性能チューニング手法について検討した。

１． はじめに2011 年 3 月現在、
1PFlops を超える性能を持つスーパーコンピュータが 7 システム
（そのうち、GPU を搭載したものは 3 システム）出現している。次世代、次々世代のスー
パーコンピュータとしては、マルチコア CPU に加えて GPU などのアクセラレータを搭載した計算ノードを数千～数万台以上接続したものが主流になることが予想される。
このようなマルチコア CPU（+マルチ GPU）から構成されるスーパーコンピュータにおいては、プロセッサコア数の増加や演算性能あたりのメモリバンド幅の不足などにより、高い実行効率を得ることが困難になりつつある。したがって、今後計算科学においてグランドチャレンジを行うためには、これまでに提案されてきた並列数値計算アルゴ
リズムや性能チューニング手法を用いるだけでは不十分である。
そこで、平成 22 年度の本研究課題においては、大規模並列環境における数値計算ア
ルゴリズムとして、・ 次世代スーパーコンピュータに向けた高速フーリエ変換（FFT）のアルゴリズム
・ GPU による 4 倍・8 倍精度 BLAS
・ GPU における疎行列ベクトル積
・ CPU+GPU 環境上で動作する実対称固有値ソルバ
・ Block Krylov アルゴリズムによる連立一次方程式の求解について研究を行うと共に、性能チューニング手法について検討した。
これらの並列数値計算アルゴリズムや性能チューニング手法は、計算物質科学の実
アプリケーションプログラムに反映させ、これまでに不可能とされてきた規模の計算を実現することができると期待される。

２． 次世代スーパーコンピュータに向けた高速フーリエ変換（FFT）のアルゴリズムに関する研究（高橋大介）
これまでに提案されてきた並列三次元 FFT における典型的な配列の分散方法として
は、三次元（x、y、z 方向）のうちの一次元のみ（例えば z 方向）を分割して配列を格納することが多かった。この場合、z 方向のデータ数はプロセッサ数以上となる必要が
ある。

－9－

次世代スーパーコンピュータにおいては、性能を向上させるためにコア数やプロセッ
サ数が増える傾向にあるため、MPI プロセス数は 1 万個以上になることもあると考えられる。したがって z 方向に一次元分割した場合には z 方向のデータ数が 1 万点以上でなければならないことになり、三次元 FFT の問題サイズに制約を受けることになる。
この問題点を改善する方法として、x、y、z 方向に三次元分割する方法[1、2]が提案されている。三次元分割では、x、y、z 方向に FFT を行う際に、その都度全対全通信でデータを交換する必要がある。それに対して、二次元分割においては、x、y、z 方向の
うち、分割されていない方向が一つあるため、全対全通信の回数を減らすことができるという利点がある[3]。本研究課題においては、二次元分割により通信時間を削減すると共に、演算と通信をオーバーラップさせることで、従来の実装に比べてさらに性能を改善した。
図 1 に T2K-Tsukuba における並列三次元 FFT の性能を示す。演算と通信をオーバーラップさせることにより、性能が改善されていることが分かる。
500
No
overlap
Overlap

GFlops

400
300
200
100

2048^31024^3

512^3

256^3

128^3

64^3

0

Length of transform N1xN2xN3

図 1．T2K-Tsukuba（64 ノード, 256MPI プロセス×4 スレッド）における並列三次元 FFTの性能

３．GPU による 4 倍・8 倍精度 BLAS の実装と評価（高橋大介）
4 倍・8 倍精度演算に対応した BLAS（Basic Linear Algebra Subprograms）関数を GPU（Graphics Processing Unit）向けに実装し評価を行った。4 倍・8 倍精度演算には double
型倍精度数を 2 つ連結して 4 倍精度数を表す double-double（DD）型 4 倍精度演算、および double 型倍精度数を 4 つ連結して 8 倍精度数を表現する quad-double（QD）型 8 倍精度演算を用いた。
NVIDIA Tesla C2050 による性能評価では、Intel Core i7 920 での同一処理と比べ、4 倍
精度 AXPY が約 9.5 倍、8 倍精度 AXPY が約 19 倍高速化された。また 4 倍精度 GEMM
は CPU に比べて約 29 倍、8 倍精度 GEMM は約 24 倍の高速化を達成した。さらに Tesla
C2050 では 4 倍精度 AXPY が倍精度演算の高々2.1 倍の演算時間となり、
GEMV、GEMM
でも倍精度演算に対する計算時間の増大が CPU の場合と比べて大幅に削減された。
一方で PCI-Express（PCIe）によるデータ転送時間を考慮した場合、倍精度 GEMM は
PCIe データ転送性能に律速される傾向が見られたが、4 倍・8 倍精度 GEMM ではこれがほぼ解消されることが分かった。

４．GPU における格納形式自動選択による疎行列ベクトル積の高速化（高橋大介）－10－

科学技術計算に用いられる連立一次方程式では、疎行列を対象とするものが多い。この場合、反復法により解を求めることになり、疎行列ベクトル積が何度も呼び出されることになるため、この疎行列ベクトル積の高速化が求められる。疎行列の格納形式には様々なものが提案されているが、対象となる疎行列によって各格納形式での疎行列ベク
トル積の性能に差が生じる。疎行列ベクトル積は演算量よりもメモリアクセス量が多い演算であるため、GPU の性能を発揮させるためには最適な格納形式を選択することが重要である。
そこで本研究では、疎行列の格納形式を自動選択することで疎行列ベクトル積の高速化を行った。まず、予備評価として 4 種類の格納形式での疎行列ベクトル積を実装し、性能を評価した。その結果、非零要素のばらつきと非零要素率によって各格納形式での疎行列ベクトル積の性能に違いが生じることが分かった。この予備評価の結果を用いて、最適な格納形式を自動選択するためのアルゴリズムを実装し、反復法を用いた連立一次方程式の求解を用いて評価を行った。
NVIDIA Tesla C2050 による性能評価の結果、今回用いた疎行列のうち多くの疎行列において最適な格納形式を選択し、CRS 形式のままで CUSPARSE [4]を用いて実行する場合と同等、またはそれよりも高速であることが分かった。また、特に反復回数が少ない
時でも本研究で提案した自動選択が有効であることが分かった。

５．Multi-{core, socket, GPU}向け数値計算ソフトウェアの最適実装・チューニング技術の研究（今村俊幸）
５．１． はじめに
2011 年 3 月現在、世界最高性能のスパコンは中国の Tianhe-1A（天河 1A）で、 GPU
を 7168 基搭載する。東工大の TSUBAME2.0 も GPU を搭載するスパコンである。現在のスパコンのトレンドはマルチ GPU+マルチコア CPU である。マルチコア CPU を使いこなすことも非常に複雑であるにもかかわらず、利用者は更に GPU の効率利用を求められる。GPU のプログラミングは NVIDIA 社の CUDA [5]に代表されるが、これまでのマルチスレッドプログラミング以上に多数のスレッドを生成するプログラミングが要請される。
本年度の研究主眼は、数値計算ソフト、特に本研究領域で非常に利用頻度の高い「固有値ソルバ」の GPGPU 化である。GPU を利用する数値計算ライブラリとして CULA [6]や MAGMA [7]が先行している。既存の LAPACK を GPGPU 化したものであるが、固有値ソルバを完全サポートしているわけではない。また、逐次計算機向けライブラリの
LAPACK をベースとしていることから、分散メモリ・クラスタへの対応には至っていない。本研究計画ではマルチ GPU クラスタ上での数値計算ライブラリの実装方法の技術確立を目標としており、既存プロジェクトに先駆けての研究開発の第一歩として、単体 CPU+GPU での固有値ソルバを実装し、プログラミング上の問題点・数値計算ソフト
ウェア構成方法の検討、性能評価などを実施する。
５．２． CPU と GPU を併用した固有値ソルバの性能測定結果
固有値計算の過程で最も負荷の高いハウスホルダー三重対角化ならびに逆変換中に現れる行列ベクトル積・行列更新操作を CUDA-BLAS で置き換え、最も簡単な GPGPU 化を実施した。CUDA 上で動作する BLAS は
NVIDIA 社が CUDA 環境に同封する CUBLAS [8]テネシー大学の MAGMA プロジェクト内の MAGMABLAS [7]が知られている。高度に最適化された関数もあるが、固有値ソルバに必要な幾つかの関数はあまり最適化されていない。したがって、最適化された BLAS 関数を作成し、これ－11－

までに別プロジェクトにて開発してきた固有値ソルバ(eigen_s)に組み込んだ。表 1 に
CUDA-BLAS を組み込んだ固有値ソルバ（eigen_sg）の性能測定結果を示す。表の上段に計算時間【秒】、下段に性能【FLOPS】を記録している。注目すべきは、市場で安価に手に入る CPU+GPU の組み合わせ（SandyBridge+GTX580）で、2000 次元の全固有値・全固有ベクトルの計算が 1 秒以下で実現されたことである。CPU+GPU の理論性能は約200GFLOPS であるので、1/6~1/3 程度の効率で計算が達成されたことになる。固有値計
算はメモリバンド幅で性能が頭打ちになる部分があることを考慮すると、この効率は極めて高いものである。
表 1．Eigen_sg の性能測定結果
行列の次元

1000

2000

3000

4000

5000

6000

7000

8000

GTX285

0.379

1.596

4.080

8.267

14.60

23.02

34.37

48.35

8.79G

16.7G

22.0G

25.8G

28.5G

31.2G

33.1G

35.2G

0.349

1.492

3.971

7.989

14.19

22.85

31.85

NA

9.52G

17.8G

22.6G

26.7G

29.3G

31.4G

35.8G

---

0.294

1.160

2.865

5.559

9.871

15.19

22.64

30.79

11.3G

22.9G

31.4G

38.3G

42.2G

47.3G

50.4G

55.4G

0.226

0.864

2.077

4.126

7.333

11.45

17.18

23.75

14.7G

30.8G

43.3G

51.6G

56.8G

62.8G

66.5G

71.8G

/Phenom II X4
GTX460/Corei7
TeslaC2050/Coire7
GTX580/SandyBridge

５．３． まとめ次世代スパコンの原型とも言われている「GPU ク ラスタ」の最小単位である
CPU+GPU 環境上で動作する実対称固有値ソルバ開発を実施した。計測に使用した計算機環境は、多くの場合、市場で安価で購入できるもので 15~40 万円程度である。開発途中ではあるが、予備実験において SandyBridge+GTX580 の組み合わせで、2000 次元の完全対角化（数学的な意味で全固有値固有ベクトルの計算）が 1 秒以下で実現された。今後は、GPU クラスタを指向した GPU 間通信部分のインフラを整備し GPU クラスタ向
け固有値ソルバの開発を進めていきたい。

６．Block Krylov アルゴリズムによる連立一次方程式の求解高速化（多田野寛人）
大規模固有値問題の一部の固有対を求める解法である Sakurai-Sugiura 法（SS 法）では、複数本の右辺ベクトルをもつ連立一次方程式の求解を必要とする。同方程式の求解時間は SS 法の演算の大部分を占めるため、求解時間の高速化が望まれている。
複数本の右辺ベクトルをもつ連立一次方程式の反復法として、Block Krylov 部分空間
反復法がある。同法で複数本のベクトルをまとめて解くことにより、1 本ずつ解いた場
合よりも少ない反復回数で解が得られることがある。しかしながら、右辺ベクトル数 L
が大きくなると、数値的不安定性により残差が発散する場合もある。Block Krylov 部分
空間反復法の 1 つである Block BiCGSTAB 法の相対残差履歴を、図 2 に示す。図 2 の赤線が示すように途中まで相対残差が減少しているものの、数値的不安定性の影響で発散に転じている。同法で現れる L 本のベクトル列に対して直交化を施すことで数値的不安定性を抑制でき、相対残差の収束性を改善できた。

－12－

101

1400.010-21200.0Time [sec]Relative residual norm

図 3 に右辺ベクトル数を変化させたときの Block BiCGSTAB 法の計算時間変化を示す。
図中の J は、前処理として用いた Jacobi 法の反復回数を表す。なお、J=0 は前処理を行っていないことを示す。前処理を行わない場合は、右辺ベクトル数 L の増加に伴い計算時間も大きく増加する。しかしながら、Jacobi 法で前処理を行うことにより、計算時間を短縮することができた。また、J=8 の場合に注目すると、L=1 と L=12 における計算時間がほぼ同じであることが分かる。これは、Block Krylov 部分空間反復法による反復回数の減少と、係数行列と複数本のベクトルの積の効率化の相乗効果によるものである。

10-5
10-8
10-11
10-14

J=0
J=2
J=4
J=6
J=8

1000.0
800.0
600.0

0

1000
2000
Iteration number

400.0

3000

図 2. L=8 の Block BiCGSTAB 法の相対残差

0

2
4
6
8
10
12
Number of right-hand sides, L

図 3. 右辺ベクトル数 L の変化による直交化

履歴。■：直交化なし、■：直交化付き。

付きBlock BiCGSTAB法の計算時間変化。

７．まとめ
平成 22 年度の本研究課題においては、大規模並列環境における数値計算アルゴリズ
ムとして、・ 次世代スーパーコンピュータに向けた高速フーリエ変換（FFT）のアルゴリズム
・ GPU による 4 倍・8 倍精度 BLAS
・ GPU における疎行列ベクトル積
・ CPU+GPU 環境上で動作する実対称固有値ソルバ
・ Block Krylov アルゴリズムによる連立一次方程式の求解について研究を行うと共に、性能チューニング手法について検討した。
来年度は、今年度の研究成果をさらに発展させると共に、計算物質科学の研究者と連携して、実アプリケーションプログラムの高速化を進めていきたい。

８． 参考文献[1] M. Eleftheriou, B. G. Fitch, A. Rayshubskiy, T. J. C. Ward, R. S. Germain: Scalable
framework for 3D FFTs on the Blue Gene/L supercomputer: Implementation and early
performance measurements, IBM Journal of Research and Development, Vol. 49, pp.
457–464, 2005.
[2] B. Fang, Y. Deng, and G. Martyna. Performance of the 3D FFT on the 6D network
torus QCDOC parallel supercomputer, Computer Physics Communications, Vol. 176, pp.
531–538, 2007.
[3] D. Takahashi: An Implementation of Parallel 3-D FFT with 2-D Decomposition on a
Massively Parallel Cluster of Multi-core Processors, Proc. 8th International Conference on
Parallel Processing and Applied Mathematics (PPAM 2009), Part I, Workshop on Memory
Issues on Multi- and Manycore Platforms, Lecture Notes in Computer Science, No. 6067, pp.

－13－

606–614, Springer-Verlag, 2010.
[4] NVIDIA: CUSPARSE User Guide,
http://developer.download.nvidia.com/compute/cuda/3_2/toolkit/docs/CUSPARSE_Library.pdf[5] NVIDIA: CUDA C Programming Guide 3.1.1,
http://developer.download.nvidia.com/compute/cuda/3_1/toolkit/docs/NVIDIA_CUDA_C_Pr
ogrammingGuide_3.1.pdf[6] J. R. Humphrey, D. K. Price, K. E. Spagnoli, A. L. Paolini, E. J. Kelmelis, "CULA:
Hybrid GPU Accelerated Linear Algebra Routines," SPIE Defense and Security Symposium(DSS), April, 2010.
[7] MAGMA プロジェクト, http://icl.cs.utk.edu/magma/[8] NVIDIA : CUDA CUBLAS Library,
http://developer.download.nvidia.com/compute/cuda/2 3/toolkit/docs/CUBLAS Library
2.3.pdf

９． 連携研究者・研究協力者
連携研究者：
佐藤三久（筑波大学大学院システム情報工学研究科教授）
朴泰祐（筑波大学大学院システム情報工学研究科教授）

櫻井鉄也（筑波大学大学院システム情報工学研究科教授）

１０．

本研究課題における平成２２年度の発表論文と招待講演

発表論文：
1)

久保田悠司，高橋大介：GPUにおける格納形式自動選択による疎行列ベクトル積の高速化，情報処理学会研究報告，2010-ARC-192，2010-HPC-128，No. 19，北海道大学，2010年12月17日．

2)椋木大地，高橋大介：GPUによる4倍・8倍精度BLASの実装と評価，2011年ハイパフォーマンスコンピューティングと計算科学シンポジウムHPCS2011論文集，pp.
148–156，産業技術総合研究所，2011年1月19日．

3)

Yuji Kubota and Daisuke Takahashi：Autotuning of Sparse Matrix-Vector Multiplication by
Selecting Storage Schemes on GPU，2011 SIAM Conference on Computational Science andEngineering，Grand Sierra Resort and Casino，Reno，Nevada，USA，2011年3月1日．
今村俊幸，山田進，町田昌彦 : eigen_sg: マルチコア+GPGPU 環境における固有値ソルバ開発, 2011 年ハイパフォーマンスコンピューティングと計算科学シンポジウム,
HPCS2011 論文集（ポスター発表）, 産業技術総合研究所, 2011 年 1 月 18 日.

4)5)

今村俊幸 : GPGPU を用いた固有値計算の可能性, 日本応用数理学会2010年度年会,
明治大学, 2010年9月9日.

招待講演：
なし

－14－

計算物質科学の基盤となる超大規模系のための高速解法
Fast solution for ultra-large scale systems as a basis of computational materials science張紹良 1、阿部邦美 2、曽我部知広 3、今堀慎治 1、宮田考史 1、山本有作 4
S.-L. Zhang, K. Abe, T. Sogabe, I. Imahori, T. Miyata, Y. Yamamoto
名古屋大学 1、岐阜聖徳学園大学 2、愛知県立大学 3、神戸大学 4
Nagoya University, Gifu Shotoku University, Aichi Prefectural University,
Nagoya University, Nagoya University, Kobe University
本領域は様々の異なる計算手法により複合相関と非平衡ダイナミクスに焦点を当てた物質デザインに挑戦するものであり、その研究過程において数多くの超大規模線形方程式と固有値問題を直面しなければならない。この領域に現れる線形方程式と固有値問題を総合的に研究し、より高効率、よりロバストな解法の系を開発・整備するため、本研究は電子構造計算のために蓄積してきたノウハウが領域内での各研究グループの有機的な結合に活用することにより、統一的な観点か
らそれらの数理的特徴を捉えたうえ、数理手法重視の高速解法の誕生を目指す。

１． はじめに既知の基礎法則や支配原理から出発して、極めて大きな自由度をもつ自然系や人工物系などの振る舞いを理解・予測するためには、計算機による超大規模科学技術計算が不可欠であり、近年の計算科学分野のめざましい進歩は高速計算機の性能によるものだけ
でなく、数値解法の驚異的な進歩にもよるものである。物質デザインコンピューティク
スでは、数百万規模以上の系が頻繁に現れて、計算時間の大半がそれを解くことに費やされると予想される。この部分の計算効率の向上は物質デザインコンピューティクスの分野においてきわめて重要である。
本研究では、物質デザインコンピューティクスに現れる様々な超大規模系の数理的諸特徴を研究すると同時に、最新の計算機を高度に駆使するための高速解法に対して総合的開発を行うことを目的とする。本研究の目標としては、解きにくい問題を簡単に、計算時間のかかる問題を高速に、計算精度の不十分な問題を高精度に解けるようにすることである。
２． 新しい高速解法の研究・開発
以下では、線形方程式と固有値問題を分けて、それぞれにおいて初年度の研究成果を紹介する。

2．1 線形方程式に対する部分空間法の開発・新型前処理の研究
応用分野の広い範囲に現れる線形方程式に対してクリロフ部分空間法の開発は８０
年代までに大きな進展があって、Bi-CG法とGMRES法がその代表として欧米の学者によって提案されて、多くの応用分野に利用されてきた。
本研究では近年この二つの解法に対して理論に基づいてその改良を目指している。
最近はGMRES法の実用版であるGMRES(m)法に対して、その収束所要時間を大幅に短縮できる改良版Look-Back GMRES(m) 法の開発に成功した。
一方、応用分野への適用の重要性を鑑み、本研究では計算物質科学のみならず、流体
力学、電気回路設計、線形計画問題、音場計算、プラズマシミュレーション、フォトニ－15－クス等の分野の専門家とも交流を深め、様々なスタイルで共同研究を進めている。とりわけ、東京大学の藤原グループが開発をしている超大規模電子構造計算プログラムパッ
ケージに対して、その問題に特化した高速解法Shifted-COCG法を提案し、数万以上の原子の系の大規模計算に対して、高速解法の数理的な突破口を開いて、さらに計算効率の向上に関する知見を深めた。
また、各反復毎に異なる前処理を適用できる概念が提案され、GMRES 法にその概念を実装した FGMRES 法や GMRESR 法が開発されている。一方、SOR 法を用いた可変的前処理付き GCR 法が提案されている。前処理を実行する際に Krylov 空間法を用いるより SOR 法の方が効果的である場合が報告されている。しかし、対角項に零を持つ行列や長方行列を係数に持つような線形方程式に対しては、通常の SOR 法を使うことはできない。そこで、特異行列や対角項に零を持つような行列に対しても SOR 法を適用できるように、反復行列を生成する際の行列分離方法を改良し、ある種の GSOR 法を提案した。数値実験では、提案した GSOR 法を可変的前処理付き GCR 法のアルゴリ
ズムに適用した際の有効性を示した。さらに、前処理のための反復が計算時間の大半を占める。そこで、残差多項式係数の計算で現れる内積演算の精度に着目し、精度が低下した場合に内部反復を打ち切る方法（エコロジカルに内部反復を実行する方法）を提案
し、計算時間を削減した。
２．２

固有値問題に対する直接法と部分空間法の開発

その１：フォトニクス結晶の固有値問題に対する高速解法
工学における科学技術計算では、大規模行列の固有値問題が様々な応用分野に頻繁に現れ、計算時間の大半が固有値問題を解くことに費やされるという事実がある。大規模固有値問題を高速かつ高精度に解くことは、計算科学の分野において極めて重要である。
固有値問題は、複素平面上に分布した膨大な数の点（固有値）の中から、工学上重要な固有値だけを選別して求める問題である。必要な固有値は工学問題によって異なるため、様々な固有値のニーズに対応した数値解法が求められている。
近年、フォトニクス結晶の応用が盛んに行われており、電子の代わりに光を利用した集積回路が注目を集めている。このような次世代集積回路の設計においては、従来と異なる固有値のニーズが生じている[1]。そのため、既存解法の効率的な適用は困難であり、新たな高速解法の開発が望まれている。本研究班は、新しい固有値のニーズに対して、数理的な側面から高速解法の研究開発を行っている。
従来、応用分野に応じて様々な固有値のニーズがある。例えば、電磁流体力学や電子構造計算の応用分野では、複素平面上において指定された点近傍の少数の固有値が必要とされており、このようなニーズに対して様々な反復法が提案されている[2]。一方、フォトニクス結晶の応用分野から生じた固有値問題では、単位円周上に存在するすべての固有値が必要とされている。そのため、必要な固有値は一点の近傍に存在する保証がなく、既存の反復法はその効率的な適用が困難である。
本研究では、必要なすべての固有値を確実かつ高速に計算することを目指し、数値解法の開発を行っている。具体的には、単位円の内部と外部にそれぞれ円を設置し、二つの円で囲まれた多重連結領域を設定する（図 1 参照、二つの円の半径の差を 0.01 とした）。この領域は単位円を含んでいるため、必要な固有値もすべて領域内部に含まれている。そのため、固有値のニーズは「指定された多重連結領域内部のすべての固有値」として捉えることが可能となる。次に、多重連結領域内部のすべての固有値を求めるた
め、既存の解法（Sakurai-Sugiura 法[3]）を拡張した。
その２：独立成分分析のための LUJ2D アルゴリズムのキャッシュ向け最適化－16－近年注目されている多変量解析の新手法として、独立成分分析がある[4]。独立成分
分析は、n 個の独立な信号源からの信号 s(t)（s(t)は n 次元ベクトル）の線形結合として観測信号 x(t)=As(t)（A は n×n 行列）が与えられたときに、x(t)のみから s(t)を復元する手法である。元々は信号処理の手法として提案されたが、最近ではタンパク質ダイナミ
クスの解析[5]などへも応用が始まっている。
独立成分分析のアルゴリズムとしては、FastICA、JADE をはじめとして様々なアルゴリ
ズムが提案されているが、大規模計算で必須となる実装の高速化については、まだあまり研究が進んでいない。そこで、安定性・汎用性の面で優れたアルゴリズムである LUJ2D[6]を取り上げ、キャッシュを意識した実装の高速化を行った。
LUJ2D は、異なる時間ラグを持つ複数の相関行列を同時対角化することで独立成分を復元する手法であり、計算パターンは固有値計算のためのヤコビ法に類似している。
このアルゴリズムに対して、演算の依存関係の解析を行い、演算順序を変更することで、収束性を損なわずにアルゴリズムのブロック化が行えることを示した 1)。ブロック化により、計算の主要部分を行列乗算の形に書き直すことが可能となる。行列乗算はキャッシュを有効に利用できるため、これにより、メモリアクセスのオーバーヘッドを削減し、高速化が可能となる。
独立成分分析のための LUJ2D アルゴリズムに対してキャッシュ向けの最適化を行い、従来法の約 1.8 倍の高速化を達成した。今後は、他の独立成分分析アルゴリズムに対しても最適化を行うとともに、タンパク質ダイナミクスの解析など、物性分野への応用も行っていく予定である。
３． 結果
香港科技大学の Prof. Lu から提供されたフォトニクス結晶の固有値問題に対して、本解法の有効性を検証した。数値実験の結果、本解法を用いて、単位円周上のすべての固有値が得られた。

図 1.

複素平面上の固有値分布（●必要な固有値、○不要な固有値）。

－17－

LUJ2D アルゴリズムのキャッシュ向け最適化の計算結果は図２に結果を示す。信号源の数が n=240 本で同時対角化を行う相関行列が 240 個の大規模問題の場合、従来法では Opteron 1.8GHz プロセッサで 2,000 秒程度の時間がかかっていた。提案手法ではそれが約 1.8 倍高速化されている。

図 2. 提案手法による高速化の効果（信号源 240 本、Opteron プロセッサの場合）
４． まとめ本研究で開発を行っている高速解法は、実問題に対する数値実験の結果、その有効性が確認された。
今後は、東京大学渡邉グループと鳥取大学星研との非線形シフト方程式・非線形固有値問題に関する共同研究を強化させながら、より大規模な実問題の解析に役立てるため、数理的な側面から解法の高速化に取り組む。
総合的開発の成果としては、我々は、専門家と利用者のインタラクションを計るため、様々なサービスを提供する計算サーバを構築し始めている。我が国に今現在そのようなシステムが存在しない現状に対して、我々はウェブサーバの構築を通して、専門家と利用者の双方にとって使いやすいインターフェスを有するものを目指す。
５． 参考文献[1] Y. Huang, Y. Y. Lu, and S. Li, “Analyzing photonic crystal waveguides by Dirichletto-Neumann maps”, J. Opt. Soc. Am. B, 24 (2007), 2860-2867.
[2] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe and H. van der Vorst, eds., “Templates for thesolution of Algebraic Eigenvalue Problems: A Practical Guide”, SIAM, Philadelphia, 2000.
[3] T. Sakurai and H. Sugiura, “A projection method for generalized eigenvalue problems using
numerical integration”, J. Comput. Appl. Math., 159 (2003), 119-128.
[4] A. Hyvärinen, J. Karhunen and E. Oja, “Independent Component Analysis”,
Wiley-Interscience, 2001.
[5] 成富 佑輔 他、 “独立成分分析によるタンパク質ダイナミクスの解析: 長時間
スケールの揺らぎ”、 第 48 回日本生物物理学会年会講演、 2010。

－18－

[6] B. Afsari, “Simple LU and QR based non-orthogonal matrix joint diagonalization”, in:Proc. of the 6th Int. Conf. on Independent Component Analysis and Blind Source
Separation (J. Rosca et al. ed.), Lecture Notes in Computer Science, Vol. 3889, pp. 1-7,
Springer-Verlag, 2006.
[7] S. Imahori, Y. Chien, Y. Tanaka, M. Yagiura, "Enumerating bottom-left stable
positions for rectangles with overlap" , 第 9 回情報科学技術フォーラム（FIT2010）、
講演論文集(第 1 分冊)、 pp. 25-30.
[8] 川島大貴、田中勇真、今堀慎治、柳浦睦憲、 "3次元箱詰め問題に対する構築型
解法の効率的実現法"、 第9回情報科学技術フォーラム（FIT2010）、講演論文集(第1分冊)、 pp. 31-38.
６． 連携研究者・研究協力者
連携研究者：杉原正顕教授（東京大学大学院情報理工系学研究科）

７． 本研究課題における平成２２年度の発表論文と招待講演
発表論文：
1) Y. Mizuno, K. Ohi, T. Sogabe, Y. Yamamoto and Y. Kaneda, “Four-point correlation
function of a passive scalar field in rapidly fluctuating turbulence: Numerical analysisof an exact closure equation”, Phys. Rev. E 82, 036316 (2010), pp. 1-9.
2) S. Imahori, Y. Chien, Y. Tanaka and M. Yagiura, "Enumerating bottom-left stable positions
for rectangles with overlap", 第9回情報科学技術フォーラム（FIT2010）、審査あり、講演論文集(第1分冊)、 pp. 25-30.
3)

川島大貴、田中勇真、今堀慎治、柳浦睦憲、 "3次元箱詰め問題に対する構築型解法の効率的実現法、" 第9回情報科学技術フォーラム（FIT2010）、審査あり、 講演
論文集(第1分冊)、 pp, 31-38.

4)

M. E. A. El-Mikkawy and T. Sogabe, “A new family of k-Fibonacci numbers”, Appl.
Math. Comput., 215 (2010), pp. 4456-4461.
M. E. A. El-Mikkawy and T. Sogabe, “Notes on particular symmetric polynomials
with applications”, Appl. Math. Comput., 215 (2010), pp. 3311-3317.
D. Aoto, E. Ishiwata and K. Abe, “A variable preconditioned GCR method using the
GSOR method for singular linear systems”, J. of Comput. and Appl. Math., Vol. 234,
No.2(2010), pp. 703-712.

5)
6)

7)

8)

9)

T. Miyata, Y. Yamamoto and S.-L. Zhang, “Performance Modeling of Multishift QR
Algorithms for the Parallel Solution of Symmetric Tridiagonal Eigenvalue Problems”,(Proceeding) Lecture Notes in Computer Science, 6082 (2010), pp. 401-412.
K. Aihara, E. Ishiwata and K. Abe, “A strategy for reducing the inner iteration counts
for the variable preconditioning GCR(m) method”, JSIAM Letters, Vol. 2(2010), pp.
77-80.
Y. Hirota, Y. Yamamoto and S.-L. Zhang, “Cache optimization of a non-orthogonaljoint diagonalization method”, JSIAM Letters, Vol. 3(2011), pp. 9-12.

10) A. Fukuda, Y. Yamamoto, M. Iwasaki, E. Ishiwata and Y. Nakamura, “A Bäcklundtransformation between two integrable discrete hungry systems”, Physics Letters A,
Vol. 375, No. 3 (2011), pp. 303-308.

－19－

招待講演：
1)

T. Miyata, T. Sogabe and S.-L. Zhang, “A Krylov Subspace Method for the Correction
Equation in the Jacobi-Davidson Method”, Numerical Linear Algebra Seminar, Gifu,
Japan (May 7, 2010)

2)

S.-L. Zhang, A. Imakura and T. Sogabe, “GMRES(m) method with Look-Back-type restartfor solving nonsymmetric linear systems”, 7th International Conference on Scientific
Computing and Applications, Dalian, China (Jun. 13-15, 2010)

3)

S.-L. Zhang, T. Miyata, L. Du, T. Sogabe and Y. Yamamoto, “A Projection Approach
Based on the Residue Theorem for Eigenvalues within a Multiply Connected Domain”,
The 2010 NIMS Conference & The Third China-Japan-Korea Joint Conference on
Numerical Mathematics, Gangneung, Korea(Aug. 19-23, 2010)4)

宮田考史、曽我部知広、張紹良、”固有値問題に対するJacobi-Davidson法について”、
線形計算研究会（東京、2010年9月21日）

5)

張紹良、宮田考史、”固有値計算分野においての実用的なニーズに対する数値解法”、第25回信号処理シンポジウム、奈良(2010年11月25-26日)

6)

S.-L. Zhang and T. Miyata, “A Numerical Method for Solving the Eigenvalues in a
Multiple Connected Domain”, 2011 LHD Workshop, Beijing, China(Jan. 15, 2011)－20－ナノ構造形成．新機能発現における電子論ダイナミクス
Electronic-Structure-Theory-Based Dynamics in Formation and Properties ofNano-Structures
押山淳 1、宮崎剛 2、尾崎泰助 3、岩田潤一 4、土田英二 5
A. Oshiyama, T. Miyazaki, T. Ozaki, J.-I. Iwata, E. Tsuchida
東京大学 1、物質・材料研究機構 2、北陸先端科学技術大学院大学 3、筑波大学 4、産業技術総合研究所 5
The University of Tokyo, NIMS, JAIST, University of Tsukuba, AIST

１． はじめに
本計画研究では、量子論の第一原理に立脚した、実空間アプローチとオーダーＮ法を計算手法の軸にすえ、計算機科学分野および実験科学分野の研究者との連携により、ナノメートル・スケールの複合構造体のナノ形状と電子機能の複合相関を明らかにし、また物質創成の場における非平衡ダイナミクスの解明により、新機能を有するナノ構造体の提唱を行うことを目的としている。
２． 概要
実空間アプローチは、マルチコア超並列アーキテクチャに適した計算手法であり、密度汎関数理論に基づく全エネルギー・電子構造計算手法として、すでに我々は、RSDFT（Real Space Density Functional Theory）コードの開発とその高度化を、計算機科学分野との共同で進めている[Iwata et al., J. Comp. Phys. 229, 2339 (2010) ]。今年度は、岩田を中心として、RSDFT コードの多機能化、Car-Parrinello Molecular Dynamics（RSDFT）手法との結合、さらには固有値問題の新たな解法である櫻井－杉浦法のコード化を行った。
またこの RSDFT の応用として、シリコン・ナノワイヤーの電子物性解明、シリコン表面での炭素ナノワイヤーの制御可能性探索を行った。
密度行列の最適化に基づく、オーダーN 手法は超大規模系に対する密度汎関数理論計算を可能にする有望手法である。そのひとつである CONQUEST コードは、宮崎を中心に開発されてきたが、今年度は、それをマルチコア超並列アーキテクチャにおいて高度化するとともに、バイオ系での有用性を調べる計算を実行した。
擬ポテンシャル・局在基底法は計算精度と計算コストのバランスのとれた手法であり、また様々な物理量を計算する上で拡張性が高いという利点を有している。そのひとつである OpenMX コードは、尾崎を中心に開発されてきたが、本年度はその手法の拡張を行った。ひとつの例は、新しい固有値解法の提唱であり、従来手法と比べて、サイズ N
に対するより低次のスケーリングを示すことを実証した。
第一原理分子動力学計算における、位相空間の探索を効率化するために、土田は原子の質量を非対角行列の形に一般化する手法の検討を行っている。本年度はその端緒的成果が得られた。
密度汎関数理論の適用物質範囲を広げる試みは重要である。それには、より量子効果を取り入れた交換相関エネルギー汎関数を構築することがひとつのアプローチである。
本年度、我々はハイブリッド汎関数を導入し、半導体、絶縁体のバンドギャップが、従来近似に比して、より実験値に近い値を生み出すことを見出した。

－21－

３． RSDFT コードの高度化と応用計算
RSDFT の機能拡張については、本年度、(1) LDA、GGA の近似の範囲でのスピン密度汎関数法の実装（押山研究室古家真之介研究員担当）、(2) Car-Parinello 分子動力学法の実装（重田、Boero）、(3) ウルトラソフト擬ポテンシャル法の実装（内田）、(4) 対称性の利用、の４点が行われた。いずれも一通りの実装を完了し、更なるチューニング
および実際的な応用計算を今後進めていく予定である。
RSDFT の中で重要な部分を占める演算として、部分空間の対角化がある。これは大規模計算において計算量のネックとなる要素の一つである。その固有値問題数値解法の新しいアルゴリズムの探索を本年度開始し、A01 高橋班との共同により、「櫻井‐杉浦法」と呼ばれる新固有値解法アルゴリズムの RSDFT への適用可能性について研究を行った。櫻井‐杉浦法は、固有値を Green 関数の極として、固有ベクトルをその留数として取り出す方法であり、複素平面上でのエネルギーに関する周回積分を数値的に実行することでこれを実現する。現在、RSDFT 上での実装を完了し、1000 原子程度の系のKohn-Sham 方程式を自己無憧着に解くテスト計算に成功している。これらのテストを通じ、櫻井‐杉浦法を用いれば、注目するエネルギー領域(例えばフェルミレベル近傍)だ
けを抜き出してバンド計算、状態密度計算が可能であること、複素周回積分の積分路を
少数の固有値のみを囲むように取ることで、各積分路について完全に並列化が可能であることなどを見出し、今後の超高並列計算機に適したアルゴリズムとして非常に有望であることがわかった。
シリコン・ナノワイヤー型電界効果トランジスター（SiNW-FET）は、スケーリング
則が行き詰った半導体テクノロジーにおいて、次世代デバイスへの突破口と目されているナノ構造である。それはワイヤー状のゲート構造によるオフ時のリーク電流の軽減と、ナノスケールのゲート長形成によるバリスティック伝導・トランジスター作用、の２点に起因している。そうした SiNW-FET のデバイス設計においては、ワイヤー軸方向の選択、ワイヤー断面形状の最適化が重要であり、予測能力の高い計算が望まれている。
本年度我々は、RSDFT コードにより、断面スケール 3nm から 6nm 程度の様々な SiNW
について系統的な計算を実行し、(100)軸 SiNW では有効質量が増加するのに対し、(110)軸 SiNW では逆に減少することを見出した。またバリスティック・レジームで重要になる伝導チャネル数を系統的に計算した。

図１：様々な断面形状の[110]-SiNW。濃灰色、薄灰色のボールが Si 原子、水素原子を表している。 (a) 直径 4.3 nm の円形断面ワイヤー、(b) [001]-楕円ワイヤー、(c) [001]ダンベル・ワイヤー、(d) [1-10]-楕円ワイヤー、(e) [1-10]-ダンベル・ワイヤー。

図１は、一例として、4 – 6 nm スケールの様々な断面形状の、(110)方向を軸としてもつ
SiNW が示してある。ワイヤーの側面は、ある種の実験に対応して水素で終端された場合を考えている。これはワイヤー原子構造として理想的な状況に対応し、量子サイズ
効果に的を絞った計算とみなせる。図１(b)と(d)は楕円断面であるが、その長軸方向は－22－非等価な結晶軸を向いているので、異なるワイヤーである。また図１(c)と(e)も同様に、
互いに非等価なダンベル断面形状をもつワイヤーである。

図 2：図１の(a)-(e)で示した[110]-SiNW のエネルギー帯。上パネルが伝導帯、下パネル
が価電子帯。伝導帯の最低点と第二最低点とのエネルギー差を図中に示す。

図２は、図１で示した各 SiNW のエネルギー帯構造である。伝導帯の最低点はいずれのワイヤーでも、ブリルアン域の中心（Γ点）に位置している。第二最低点は、域境界との途中に位置し、最低点よりのエネルギー差∆Eshift は量子閉じ込め効果によって決まり、ワイヤー形状、サイズによって異なる。Γ点での最低点は、もともとの Si 結晶の伝導帯バレーの内、[001]、[00-1]方向に位置する２つのバレーに起因し、一方第二最低点
は、[010]および[100]方向に位置する２つのバレーに起因している。従っていずれの最低点も二重に縮退している。∆Eshift は図１(b)の[001]-楕円ワイヤーが最も小さいことがわかった。これは、量子閉じ込め効果による、ワイヤー垂直面内の定在波の波数が断面形状によって異なることに起因している。これより、[001]-楕円断面形状の[110]-SiNW
が、その低バイアス下でのチャネル数の多さにより、FET に最適であることが予想される。また今回の計算により、[110]-SiNW では有効質量がバルクの値より軽くなることもわかった。
今後さらに、実デバイスサイズ（～10 nm）で、かつ実プロセス等によるラフな側面をもつナノワイヤーの電子状態を明らかにし、得られた結果とコンパクトモデルを利用して、シリコンナノワイヤのバリスティックなトランジスタ特性(電流‐電圧特性)の評価を行う。
４． 密度行列最適化法によるオーダーN 法（CONQUEST）の開発とバイオ系への応用
生体系に対するこれまでの理論研究は、古典力場を用いた研究がほとんどであるが、計算量が少ないという利点の代償として、結果が力場のパラメータに強く依存し、また化学結合の切断や再結合を表すことは不可能に近い。このため、化学結合の状態が変化するような重要な領域に対してのみ量子論（QM）にもとづいた手法を用い、それを囲む領域は古典力場で扱う（MM）という QM/MM 法なども適用されているが、これらの手法には領域の設定（QM 領域のサイズ、境界の扱い、境界の更新等）に問題がある場合がある。オーダーN 法第一原理計算手法は巨大な領域に対して全原子第一原理計算を可能とするので、これらの問題を解決することが期待され、生体系はそのオーダーN 法
第一原理計算の適用が望まれる重要な分野の一つである。
本年度は、生体系に対するオーダーN 法第一原理計算として、ジヒドロ葉酸還元酵素（DHFR）とイオンチャネル gramicidin A（GA）の系を扱った。この二つの系に対して、オーダーN 法第一原理計算による本格的な理論研究を行うための準備計算と、生体系へ－23－

図 3：水中の DHFR の系（16, 042 原子）に対してオーダーN 法第一原
理計算によって得られた電荷密度分布図

図 4：オーダーN 法第一原理計算でセルフコンシステン
ト計算を行った、イオンチャ
ネル gramicidin A の系（約15,500 原子）の応用計算に必要な手法の整備、開発を行った。具体的には、１）実験結果と古典分子動力学により、計算する系を構築、２）生体系のような複雑系での我々のプログラムの並列化効率を調べ、実際に大規模計算を実現すること、３）セルフコンシステント計算
を安定、効率的に行う手法の確立、４）計算された電荷密度の解析、可視化方法の開発、を行った。DHFR に関しては、実験結果から得られる構造の周辺に水分子を加えた後に
古典分子動力学を行うことによって、いくつかの構造を作成した。この際に、水分子の配置を均等にすると同時に DHFR の構造が実験構造と大きく異ならないようにするためには、DHFR（と dihydrofolate や NADPH 分子）を構成する原子の動きに適切な束縛条件を考慮しなければならないことが分かった。得られた構造に対してオーダーN 法第一原理計算を適用した結果、セルフコンシステント計算が安定して行える事も分かった（図３）。さらに、密度行列最適化のステップとセルフコンシステントな電荷密度を得るためのステップという二つのステップを同時に行うことにより、より高速かつ安定した計算が可能であることが分かった。また、DHFR と Folate の合体系と孤立系との差電荷密度を求めることによって、相互作用によって変化を受けやすい部位の解析が行える
事を確認した。DHFR による酵素反応の様々な過程に対する計算は次年度以降に行う事を予定している。
また、gramicidin A の系に関しても、イオンチャネル GA を囲む脂質二重膜、さらに、この膜の上下の水の領域を含んだ複雑系（図 4、約１万６千原子系）に対して、セルフ
コンシステント第一原理計算を実現することに成功した。この際に、密度行列を表すための局在軌道（サポート関数）は部分系に対する第一原理計算を用いて最適化すること－24－によって得ている。この手法の精度を調べ、全体系の計算において、各原子に働く力を計算した。本格的な計算や詳細な解析は今後行う。
５． 低次の計算オーダーを有する新しい固有値解法の開発
我々が開発してきた OpenMX コードは、擬ポテンシャル・局在基底法であり、計算
精度と計算コストのバランスがとれた手法である。平成 22 年度は、OpenMX のさらなる高精度大規模計算への拡張に向けて、近似を導入することなしに低次の計算オーダーを持つ新しい対角化手法の開発に成功した。
2

4

10

1D
2D
3D
1D (fitted)
2D (fitted)
3D (fitted)

0

10

O(N

)

Conventional (-4130.938589871644)
New method (-4130.938589871899)

0

10
1.90

O(N

−2

10

2.35

0.90

O(N

Norm of residual

Elapsed Time (sec.)

10

)

2

(log2N) )

-4

10

-8

10

-12

10

−4

10

10

1

2

10

10

3

4

10

Number of Atoms

10

20
SCF step3040

図６：DNA 分子の SCF 計算

図５: 逆行列の計算時間

Car-Parrinello 法の発展を契機に、第一原理密度汎関数計算の大規模系への適用限界を拡張するために非常に多くの対角化計算手法がこれまでに開発されている。これらの対角化手法は(i) 共役勾配法や Davidson 法などの繰り返し計算に基づくものと、(ii) 近似を導入して計算量を O(N)に低減したオーダーN 法に大別される。前者においては与えられた Kohn-Sham ハミルトニアンに対して近似を導入することなしに、固有状態の計算を行うが、占有軌道間の直交条件を保持する計算に対して O(N3)の計算量が必要となる。その中間に位置する O(N2～)法の開発はこれまであまり行われていないが、近似を導入することなしに計算オーダーが低減できるならば、O(N3)法に代わる手法として興味の持たれる所である。今回、擬原子局在基底や有限要素法などの局在基底関数の場合に対して、近似を導入することなしに低い計算オーダーを持った対角化手法の開発に取り組んだ。そのアイデアは留数積分によって直接、密度行列を計算するというものであり、収束に必要な極の数が系のサイズに依存しないことを示すことができる。主要な計算は疎行列の逆行列計算(グリーン関数の計算)となり、ここで計算オーダーが支配される。我々は逆行列要素を計算するために、行列の Nested Dissection(ND)法と Shur 補元を用いて漸化式を導出し、これを用いて必要な逆行列要素を求める新しい計算手法を開発した。この漸化式に基づく逆行列計算法の計算オーダーは 3 次元系で O(N7/3)、2 次元系で O(N2)、1 次元系で O(N(log2N)2)であり、近似を導入することなしに計算オーダーを低減することが可能となった。
理論的解析から得られた計算オーダーが実際の計算でも成り立っていることを確認するため、隣接強結合モデルで記述された一次元系(直鎖分子)、二次元系(正方格子)、
三次元系(立方格子)に対して逆行列の計算時間の実測を行った(図 5)。実線は実測値へのフィッティングから得られたものであり、実測の計算オーダーは 3 次元系で O(N2.35)、
2 次元系で O(N1.90)、1 次元系で O(N0.90(log2N)2)と見積もることが出来る。この結果は理論値と良く一致しており、実測値からも低次スケーリングが確認された。次に本手法が数値厳密であることを実証するために、有限数の極(80 極)で DNA の SCF 計算を行ったものが図６である。従来法と比較し、収束特性は同一であり、また得られた全エネル－25－

ギーは 13 桁一致している。本手法で用いられている留数積分は 100 極程度で倍精度限界まで収束する手法であり、例え、有限数の極を用いたとしても本手法は数値厳密な手法であると言える。各極に対する計算は独立に実行できるため、本手法は並列計算に適した手法である。我々は本手法を OpenMP と MPI を用いてハイブリッド並列化した。
速度向上率を図７に示す。160 プロセス、４スレッドを用いた速度向上率はおよそ 350倍である。通常の対角化法と比較し、大幅な並列効率の向上が得られているが、今後、さらに並列化法を工夫することで、並列化効率の向上が必要である。
1 thread2 threads
4 threads
Conventional (1 thread)

Speed−up Ratio

400

300

3.76 sec.

200

100
7.09 sec.
0
0

20

40

60

80

100 120 140 160 180

Number of Processes

図 7：OpenMP/MPI ハイブリッド並列による速度向上比

図 8：本手法（MTMD）と従来の方法（REF）
による、酸素間の動径分布関数の比較

６． 一般化された質量を用いた第一原理ＭＤの高速化
第一原理分子動力学計算を行う場合には、常に計算の負荷が非常に高いという難点があり、高速化するような手法が強く望まれている。本年度は、原子の質量を一般化し、有限温度における位相空間の探索を効率化するようなアルゴリズムについて研究を進めた。通常、原子の質量は原子の座標に依存しない対角型の行列で与えられるが、この方法では原子座標に依存するような非対角行列に拡張する。即ち、ハミルトニアンは以下の式で表せる：H (q, p) =1
pi M ij−1 (q) p j + U (q)∑2 i, j

この行列Ｍの成分を古典的な力場で近似することで、位相空間が等方的になるような変換を行うことができる。この際、力場のパラメータについては、第一原理計算の結果から Force Matching 法を用いてフィッティングを行った。また、運動方程式の数値積分についても効率の良い方法を開発した。液体の水（６４分子系）について検証を行い、アンサンブル平均を不変に保ちつつ３倍程度の高速化を実現できることを示した。この
他にも、座標に依存しない範囲で最適な原子質量を見積もる方法について提案した。また、元のハミルトニアンに対応するダイナミクスを復元する可能性についても幾つか議論した。
本手法は特に平衡状態における平均構造や自由エネルギー計算において有用であると期待できる。また、オーダーＮ法のような電子状態計算を高速化するような手法とも両立することが可能である。
７． 連携研究者・研究協力者
連携研究者：
重田育照（大阪大学基礎工学研究科准教授）、内田和之（東京大学工学系研究科助教）、－26－

平山博之（東京工業大学総合理工学研究科教授）、研究協力者：
D. R. Bowler （University College London, Reader）、M. Boero （University of Strasbourg,
Professor）
８． 本研究課題における平成２２年度の発表論文と招待講演
発表論文
1)

J.-I. Iwata, A. Oshiyama adn K. Shiraishi, ``Large-Scale First-Principles Electronic Structure
Calculations for Silicon Nanostructures" Proc. 30th Int. Conf. Physics on Semiconductros(Seoul, July 26 - 30, 2010).

2)

Y. Gohda and A. Oshiyama, ``First-principles calculations on spin polarization of vacancies
in nitride semiconductors" Proc. 30th Int. Conf. Physics on Semiconductros (Seoul, July 26 30, 2010).

3)

T. Kurita, K. Uchida and A. Oshiyama, ``Atomic and Electronic Structures of α-Al2O3
Surfaces", Phys. Rev. B 82, 155319 (2010).

4)

K. Uchida and A. Oshiyama, ``Structure, spin, and charge of decavacancy V10 in crystalline
silicon" Proc. 30th Int. Conf. Physics on Semiconductros (Seoul, July 26 - 30, 2010).

5)

D. Sung, J.-I. Iwata, ``Structural stability and energy bands of Si nanowires along [110]direction" Proc. 30th Int. Conf. Physics on Semiconductros (Seoul, July 26 - 30, 2010).

6)

M. P. Sena, T. Miyazaki and D. R. Bowler, “Linear Scaling Constrained Density Functional
Theory in CONQUEST”, Journal of Chemical Theory and Computation, 印刷中。(Mar 7,
2011, on web)7)

M. Toyoda and T. Ozaki, "Fast spherical Bessel transform via fast Fourier transform andrecurrence formula", Comp. Phys. Comm. 181, 277 (6 pages) (2010).

8)

T. Ozaki, K. Nishio, H. Weng, and H. Kino, "Dual spin filter effect in a zigzag grapheme
nanoribbon", Phys. Rev. B 81, 075422 (5 pages) (2010).

9)

M.J. Han, T. Ozaki, and J. Yu, "Effect of on-site Coulomb interactions on the electronic
structure and magnetic property of Gd2 cluster", Chem. Phys. Lett. 492, 89-92 (2010).

10) M. Toyoda and T. Ozaki, "LIBERI: Library for numerical evaluation of electron-repulsion
integrals", Comp. Phys. Commun. 181, 1455-1463 (2010).
11) K. Nishio, T. Ozaki, T. Morishita, and M. Mikami, "Tunable electronic transport properties ofsilicon-fullerene-linked nanowires: Semiconductor, conducting wire, and tunnel diode", Phys.
Rev. B 81, 115444 (11 pages) (2010).
12) T. Ozaki, "Efficient low-order scaling method for large-scale electronic structure calculations
with localized basis functions", Phys. Rev. B 82, 075131 (17 pages) (2010).
13) G. Otero, G. Biddau, T. Ozaki, B. Gomez-Lor, J. Mendez, R. Perez, and J.A. Martin-Gago,"Spontaneous Discrimination of Polycyclic Aromatic Hydrocarbon (PAH) Enantiomers on a
Metal Surface", Chem. -A. Euro. J. 16, 13920-13924 (2010).
14) E. Tsuchida, “Ab initio mass tensor molecular dynamics”, J.Chem.Phys. 134, 044112 (2011).

招待講演
1) A. Oshiyama, “Current Status of Density-Functional-Theory-Based Calculations for Nanoand Bio-Materials”, Int. Sympo. “Nanoscience and Quantum Physics 2011” (Tokyo,
January 26-28, 2011).
2) T. Ozaki, "Numerically exact low-order scaling method for large-scale electronic structure
calculations", Psi-k 2010 Conference, Henry Ford Building, Berlin, Sept. 12 - 16, 2010.
3) 尾崎泰助, "Localized basis method in OpenMX", 第１７回コンピュテーショナル・マテ
リアルズ・デザイン （CMD®）ワークショップ, 大阪大学サイバーメディアセンター,
2010年9月6日-10日.

－27－

第一原理分子動力学法による構造サンプリングと非平衡ダイナミクス
First-Principles Calculation of Structure Sampling and Non-Equilibrium Dynamics
常行真司 1，吉本芳英 2，山内 淳 3，大谷 実 4
S. Tsuneyuki, Y. Yoshimoto, J. Yamauchi, M. Otani
東京大学 1，鳥取大学 2，慶應大学 3，産業技術総合研究所 4
The University of Tokyo1, Tottori University2, Keio University3, AIST4

原子間相互作用の非調和性が本質的に重要となる大きな原子変位を伴う非平衡物
理現象の予測とダイナミクスの解明を目指し、平面波基底第一原理分子動力学法
プログラムの高速化を目指したアクセラレータ（Grape-DR）の利用、第一原理非調和格子モデルの開発と熱伝導計算への応用、有効遮蔽媒質法（ESM 法）の拡張を行った。またシリコン結晶中のホウ素欠陥の XPS スペクトル計算、窒化物／硼
化物界面の構造と電子物性計算を行った。

１． はじめに本研究グループでは、次世代半導体デバイスや熱電素子、電池等エネルギー変換素子への応用を念頭に、第一原理分子動力学法を用いてナノ構造体や新材料の熱科学の解明を目的とする。具体的には、材料およびナノ構造体の熱伝導度、熱膨張率、熱破壊の前駆現象、固液相変化とナノスケールでの相関や揺らぎ、分子固体中や分子/電極界面での電子移動による再配置エネルギーと電子移動度など、原子間相互作用の非調和性が本質的に重要となる大きな原子変位を伴う非平衡物理現象の予測と、ダイナミクスの解明を目指す。
上記のような物理量を意味のある統計量として計算し、物理現象を正しく理解・予測するためには、これまでにない大規模かつ長時間のシミュレーションと統計的なサンプ
リングが必要であり、これを第一原理分子動力学法だけで達成することは、ペタフロッ
プス級の次世代スパコンをもってしても不可能である。そこで本研究では、平面波基底関数を用いた第一原理計算コードを高速化して長時間シミュレーションを実現し、各種
構造計算、動力学計算に適用する。また比較的短時間の第一原理分子動力学法シミュ
レーションを用いて原子間相互作用の有効モデルを導出し、それを高速な古典分子動力学法に適用することによって、必要とされる長時間シミュレーションを達成する手法を開発する。
２． 平面波基底第一原理計算プログラム xTAPP の開発（吉本）
物性物理学の父母は量子力学と統計力学である。そのため第一原理電子状態計算による構造サンプリングが自由に行えるようにすることは大変重要である。しかしながら第一原理電子状態計算のコストは大変大きいため、構造の多様性が重要になってくる数
百原子以上の系を経済的な問題なく自由に扱うことは困難である。計算機の経済性を向上させる要素としてアクセラレータなどが開発されており、したがってこれらの技術への対応を進める必要がある。そこで吉本は A02 の稲葉班の協力を得て、物性における第一原理電子状態計算の標準的枠組み（平面波基底＋密度汎関数法＋擬ポテンシャル）を実行するプログラムパッケージ xTAPP を稲葉班で開発されたアクセラレータ－28－

Grape-DR を使えるように改変した。
この開発の狙いは、その能力がよくわかっている平面波基底第一原理電子状態計算プ
ログラムをたたき台として、計算科学と計算機科学の間の対話を実際に行い、課題を認識することで、A03 今田班が開発している他の第一原理的手法などさらに計算機資源を必要とするその他の手法を経済的に実行できるようにする土台作りをすることにある。
吉本は、4 chip モデルの Grape-DR を稲葉班から借受け、まずは Grape-DR システムを構築することから始めた。（図 1）その結果判明したことは、計算機科学の研究の成果である Grape-DR を計算科学側でシステムが正常に動作するようにするだけでも、想定以上の人的コストが必要になることである。
具体的に発生したトラブルの中で大きな物に、Grape-DR のボードの firm ware が稲葉班の認識とは異なるものであること判明するのに時間がかかったことがある。結局計算科学側で状況証拠の積み上げを行った上でボードを鳥取から東京に持ち込み検証することで問題の存在を証明して解決することができた。
問題解決に必要であったことをまとめると、（１）正しい情報を入手すること（２）
計算科学側でのシステム動作の解析（３）現物の確認、である。Grape-DR のようにそれ自体が研究の対象物であるものを用いて研究をする場合、（A）連絡を密にすること、（B）
実物を共有し動くところまで指導すること、の二点が重要であるとまとめ直すこともできる。このことは計算機科学の成果物を計算科学が利用する場合だけではなく、その逆に計算科学の課題を計算機科学に持ち込む場合にも同様に成り立っているであろう。
正常に動作するようになった Grape-DR システムへの xTAPP の移植においては、そのプログラムの構成において計算負荷が大きい部分を数理的な構造（行列積など）に着目して取り替え可能なモジュールとして抽出する作業が役立った。
ここでいうところの行列積構造を取り出すときに具体的に必要となるのは、角運動量など短い添字をその親である原子の番号と合併して一つの長い添字を構成することである。このようにコードの構造をまとめ直すことで、アクセラレータなどを活用する上で有利な数理構造である行列積構造を取り出すことができる。
Grape-DR には付属の DGEMM ライブラリが用意されているので、これを用いてこの行列積の数理構造を担当するモジュールを作成し、xTAPP の Grape-DR 対応化を行った。
火成岩の一種 Pyrope の 640 原子系を用いたベンチマークの結果を表１に示す。実行
環境は Core i7 960 である。Core i7 とあるのは行列積相当の部分の実行に CPU を用いたもの、そのうち MKL とあるのは数値演算ライブラリ MKL の DGEMM を用いてこれを実行させたもの、hand opt とあるものは行列積相当の部分を吉本が手で書いたものである。なお MKL や Grape-DR を使うためにはデータ形状の変更のためメモリが余分に必要になるが、直接書く場合には必要ない。
表

行列積[s]比FFT [s]
固有系 [s]その他 [s]
合計 [s]比1

図1

Grape-DR Core i7 MKL Core i7
DGEMM DGEMM
hand opt
3346
9989
33872
1
2. 99
10. 12
1335
1333
1339
338
335
337
1770
1774
1764
6788
13431
37313
1
1. 98
5. 50

hand opt と MKL DGEMM の結果の比較から、行列積構造を取り出した上でそれを－29－

DGEMM で処理することで相当に高速化することがわかる。そして Grape-DR を用いるとそこからさらに２倍高速化する。
このように、コードの構造を計算機科学的な意味での DGEMM にしておくことは、
Grape-DR が使えるようになるのみならず、ホストの CPU そのものをうまく扱うことにつながっている。また GPGPU においても高性能な DGEMM ライブラリは入手可能であるから GPGPU を活用できるようにもなっているはずである。
一方で Grape-DR と Core i7 の最大演算能力の差は 10 倍以上あるはずであるが実際の行列積部分の差は３倍にとどまっており、Grape-DR の実効性能に課題があることがわかる。また Grape-DR のデータにおいて行列積以外の部分に必要時間な時間は行列積にかかる時間に比べて小さくはない。このことはプログラム全体のさらなる高速化には
その他の部分の高速化を併せて行う必要があることを示している。
３． 第一原理非調和格子モデルの導出と熱伝導計算（常行）
格子熱伝導率の第一原理計算は、新規材料設計を行う上でその重要性を増している。
例えば、高効率の熱電変換材料を探索する際にはその性能指数 ZT が指針となるが、この値を理論的に評価するためには格子熱伝導率の高精度計算が必要になる。格子欠陥
や表面・界面および複雑なナノ構造を持つ物質の格子熱伝導率を計算するには、広く用いられてきた Boltzmann 方程式よりも非平衡分子動力学(NEMD) 法が適している。この手法では、系の両端に温度差を設け、その際に系に生じる温度勾配と熱流から Fourierの法則に従って熱伝導率を評価する。ただし、NEMD と第一原理分子動力学との組み合わせは計算コスト的に難しく、特定の原子種に対して良い近似となる古典的なモデル
ポテンシャルを用いた計算に限られていた。
この原子種に関する制限をなくし、任意のナノ構造を持つ任意の物質の熱伝導度を決定するため、我々は格子熱伝導率を第一原理的に評価する新たな計算手法の開発を行った。この手法では、NEMD を行う際に用いる原子間ポテンシャルとして、第一原理に基づく非調和格子モデル[1]を用いる。これは、全系のポテンシャルエネルギーV を各原子の平衡位置からの変位 {ui } に関して Taylor 展開を行い、熱伝導率に本質的な 4 次の非調和項まで含めたものである。

V ≅ V0 +1
11
Φ ij ui u j + ∑ Φ ijk ui u j uk + ∑ Φ ijkl ui u j uk ul∑
2 i, j
3! i , j ,k4! i , j , k ,l

ここで i,j,k,l は原子の番号と x,y,z 座標を表すインデックスである。展開係数 Φ ij , Φ ijk , Φ ijklには系の対称性によってさまざまな制限が課せられるが、欠陥構造などにも使える汎用性を重視して、最小限の（自明な）制限だけにとどめる。和は原子間距離があるカット
オフ以下のものについてとる。
原子が結晶格子の平衡点近傍で動く限りにおいて、上記のモデルポテンシャルは極めて汎用性の高いものである。モデルパラメータである展開係数 Φ ij , Φ ijk , Φ ijkl を決定するに当たっては、まず有限温度の第一原理分子動力学計算を短時間（数 ps 程度）実行し、原子変位の時間発展 {ui (t )} と原子に働く力 {Fi (t )} の時系列を得る. 次に、第一原理計算で得られた力が、上記のモデルポテンシャル V の座標微分を使って最も良く再現される
ように、最小自乗法で展開係数を決定する。その際、系の対称操作を用いて変位と力のデータを複製して利用しすることで、系の対称性を反映した展開係数を得ることができる。モデルポテンシャルは展開係数について線形であるため、これは線形最小自乗法の問題となり、安定して解くことができる。

－30－

シリコンとダイヤモンドに対して本手法を適用し、格子熱伝導率の評価を行った。その結果、システムサイズの増加に伴って熱伝導率が非線形に増加する、準バリスティッ
ク領域に特徴的な振る舞いが見られた（図 2）。またシリコンのほうが非調和効果が大きく、熱伝導度が桁違いに低いことが確認された。ただし現状の非調和格子モデルは高温シミュレーションで不安定になるため、パラメータ決定法やモデルの高精度化などを検討中である。

図 2. 非平衡分子動力学法で計算されたシリコンとダイヤモンドの熱伝導度のシステム
サイズ Lz に対する依存性。
４. シリコン結晶中のホウ素欠陥の XPS スペクトル計算（山内）シリコン結晶中の不純物元素の原子形態に関する知見は、基礎科学上の興味だけにとどまらず、半導体デバイス微細化におけるスケーリング則を維持するという応用面でも重要である。一方で、表面系における走査型プローブ顕微鏡のような直接的観測手段がないために、半導体中の不純物欠陥構造の同定は困難であり、種々の間接的な観測手段を総合して原子構造を決定している。
このような間接的な手段の一つとして内殻電子による X 線光電子スペクトル（XPS）
測定があげられるが、格子欠陥構造に関しては、これまで信頼性のある結果はほとんど得られていなかった。それには幾つかの要因が絡んでおり、その中でも重要なものは計算上の境界条件の問題である。XPS 観測値は電子の励起エネルギーであるが、そのエネ
ルギーの基準値を互いに比較する個々の欠陥系で統一する必要がある。しかしながら、
欠陥系は周りの結晶に対して歪などで弾性的な影響を及ぼすに留まらず、欠陥種により余分な電子、ホール等を放出して静電ポテンシャル的な意味での境界条件に大きな影響を与える。境界条件を揃え、理想的にはバルクのシリコン結晶と同等の条件を実現する
ためには、基本的には計算するモデルを大きくすることが必要である。大きな系を取り扱うためには平面波基底の擬ポテンシャル法が計算精度並びに効率から定評があるが、XPS の中心的役割を果たす内殻電子を naive な擬ポテンシャルでは取り扱えず工夫が必要となる。
本報告では、シリコン中の B を含む欠陥系について、これらの計算手法の有効性を確認し、XPS 計算法として知られる frozen orbital 近似、Slater の遷移状態（STS）法、
ΔSCF 法について比較した結果について述べる。
XPS 計算用に内殻原子の効果を擬ポテンシャルで取り扱う方法としては、Schefflerのグループにより提案された screened core-hole 付きの擬ポテンシャルが最もよく知ら－31－

れている[2]。この擬ポテンシャルは XPS に関与する内殻電子を hole にして、valenceの最外殻軌道に１電子付け加えた配置で作成するもので、ΔSCF 計算を行うことで相対
XPS シフト量を計算する。本報告では、これとは異なり、Vanderbilt により提案された
超ソフト擬ポテンシャル（USP）[3]の多重参照機能をコア電子にも適用することによって計算を行う。この方法は、前者に比べて hole となる状態を凍結していないので精度の上がること、spin を露わに取り入れることができるので、強磁性状態計算などにも適用できると期待される。一方で、作成方法からわかるように原子番号の大きな元素に適用することには困難が予想される。
USP の多重参照機能を用いた方法の有効性について、テスト系として C と B を含む分子に関して frozen orbital 近似、STS 法、ΔSCF 法について調べてみた。結果は、frozen orbital
近似は誤差が非常に大きかったのに比べて、STS 法、ΔSCF 法は共に 10%の精度で実験を再現した。これは XPS 計算における電子緩和の重要性を反映している。境界条件に関しては、立方体のスーパーセルモデルを用いて、セル境界付近に位置する Si の 2s 準位を probe
としてプロットした結果を図 3 に示す。青で描かれているのが純粋な Si 結晶のもので、その他の記号は B を含む欠陥を意味する。これから 512Si 原子セルを用いると 0. 1eV 程度の精度で XPS スペクトルが議論できるということがわかる。この計算では XPS の絶対値は計算できず、相対値のみ得られるので、置換配置 B（図中 Bs）からの相対値を使って、イオン注入によって作られた実験結果[4]と比較してみると、B12 クラスター（図中 ICO，S4）、<001>B-Si クラスター等の存在を強く示唆する結果が得られた。

Si2s BINDING ENERGY
(eV)

ICO
S4
T4
Si only
splitB2
O6
BCH
Bs
B+VH3
Bs-Bs
B+V
64
216
SUPERCEL L SIZE

512

図 3 スーパーセル中の Si 2s 軌道エネルギーのセルサイズ依存性. ICO, S4, T4, splitB2,
O6, BCH, Bs は、スーパーセル中の欠陥構造モデルの違いを表す。
５. 有効遮蔽媒質法の拡張（大谷）有効遮蔽媒質法（ESM 法）は電池などの電極表面で起こる電気化学反応のダイナミクスを追うために開発された手法である[5]. ESM 法では、系に 3 次元の周期境界条件を課して行う計算手法において、表面垂直方向への周期境界条件を外し、代わりに任意の境界条件を課す事を可能とする。これにより、従来は難しかった電極表面に電圧を印加した状況下でのシミュレーションが可能になった。ESM 法を用いて実際に水の電気分解反応の一部の

－32－

(a)

(b)

図 4．(a)従来の ESM 法での境界条件。(b)拡張された ESM 法での境界条件。 は ESMの誘電率を表し、 は灰色の ESM が置かれた位置を示す。
シミュレーションが実際に行われ、電気化学反応の様子や電圧印加による界面構造の変化などが第一原理計算で初めて明らかにされた[6]。
このように ESM 法では従来は不可能だった計算を可能にしたが、まだ改善の余地がある。
その中でも最も重要な問題は、“計算対象の表面/界面と ESM の間には真空領域が必要”ということである。
従来の ESM 法では図 4(a)のような誘電率を用いており、z=z1 で誘電率は不連続に変化する。これにより、電子が z>z1 の領域に進入すると静電ポテンシャルの 1 次微分が不連続なる。このような不連続性は計算の安定性を損なう恐れがあるために、従来は電子が z>z1 の領域に進入しないように十分な真空領域を z<z1 に導入する必要があった。この真空領域に印加されている電圧の定義が難しいために、電極表面に印加している正味の電圧を定義できなかった. また、真空領域に原子が進入しないようにバリアポテンシャルを z<z1 の領域に導入する必要があり、モデリングの自由度を損なっていた。
今回は図 4(b)のような滑らかな誘電率の変化を表現できる関数を導入することによって、電子が z>z1 の領域に進入しても静電ポテンシャルの 1 次微分の不連続性を除去することができた。これにより、より自然な電気化学系のモデリングが可能になった. 新たな関数に関しても解析的にグリーン関数を計算することが可能であり、従来のプログラムへ加える修正は僅かである。
z=z1 におけるグリーン関数の連続性より、静電ポテンシャルの値と 1 次微分の連続性は保証されている。新たな境界条件の有効性を確かめるために、水分子を用いて ESM との間に働く相互作用を計算した。図 5 に示すように、水分子を向きは固定したままで ESM に近づ
けていき全エネルギーの変化を計算した。ESM の遠方から徐々に水分子を近づけていくと、

図 5．ESM と水分子の相互作用エネルギー。インセットは ESM と水分子の相対的位置関係のイメージ。

－33－

引力的相互作用が働き、1. 5Å 付近に最安定となる領域がある。ESM は古典的な金属のように振る舞うために、水分子に鏡像力が働き、引力的な相互作用となる。
さらに、水分子を近づけると逆に斥力的相互作用が働き ESM 領域には水分子は進入しな
いことがわかる。金属的な ESM が z<z1 の領域に存在するイオンの静電力を完全に遮蔽するために、z>z1 の領域に進入した電子にはイオンからの束縛力が働かず、エネルギーが上昇し斥力が働いていると考えられる。
これらの ESM が電子及びイオンに及ぼす影響は、あくまでも第一原理計算で扱う領域に対する有効的な効果であり、最安定位置やエネルギーの深さに物理的な意味は無いことに注意する必要がある。
６. 窒化物／硼化物界面の強磁性（常行）
平面波基底第一原理計算を用いて、仮想的な AlN/MgB2(0001)界面の安定構造と電子状態を調べ、最安定状態で強磁性が発現することを見出した[7]。スピン分極は界面に局在しているという意味で 2 次元的であり、強磁性状態の起源は窒素原子の 2p 状態のスピン分極によるものである。界面窒素 2p 状態の局所状態密度を計算し鉄の状態密度と比較した結果、Fermi 準位近傍においてスピン分極によりピーク位置が Fermi 準位からずれているという意
味において類似している事が分かった。また窒素 2p 状態の局所状態密度は Bloch 状態となった後においても完全に 2 重縮退していた。このため、この界面における界面強磁性は Hund
結合と Fermi 準位近傍の高い局所状態密度によるバンド強磁性であると考えられる。これに対して先行研究のある GaN/ZrB2(0001)界面では、金属誘起ギャップ状態(MIGS)の効果によりスピン分極が得られない事も分かった。さらに界面構造のコンダクタンスを計算した結果、スピン分極した電子輸送特性が得られ、スピントロニクスへの応用の可能性が示された。
７. まとめ
本年度の成果は以下のとおりである。
(1) 動力学の第一原理計算に用いる平面波基底第一原理計算プログラム xTAPP を、稲葉班で開発されたアクセラレータ Grape-DR を使えるように改変し、問題点の洗い出しを行った。
(2) 物質の熱伝導度の第一原理計算を実現するため、短時間の第一原理分子動力学シミュ
レーションから非調和格子モデルを導出する手法を開発した。また得られたモデルを使った非平衡分子動力学法シミュレーションを行い、シリコンとダイヤモンドの熱伝導特性を調べた。
(3) 超ソフト擬ポテンシャルの多重参照機能を用いて、結晶中の欠陥による内殻 XPS を
高精度で計算する手法を提案し、Si 結晶中の B 欠陥について、0. 1eV 程度の精度で実験を再現できることを示した。
(4) 新たな境界条件を導入することにより ESM 法を拡張した。これにより、従来の ESM 法の問題であった“真空領域の導入”と“原子に対するバリアポテンシャルの導入”の二つの問題が同時に解決された。また、真空領域の除去により電極表面への印加電圧が明瞭に定義することが可能となった。
(5) 平面波基底第一原理計算を用いて、仮想的な AlN/MgB2(0001)界面の安定構造と電子状態を調べ、最安定状態で強磁性が発現することを見出した。

８. 参考文献[1] K. Esfarjani and H. T. Stokes, Phys. Rev. B. 77, 144112 (2008).
[2] E. Pehlke and M. Scheffler, Phys. Rev. Lett. 71 2338 (1993).
[3] D. Vanderbilt, Phys. Rev. B 41 7892 (1990).
[4] I. Mizushima et al. , Appl. Phys. Lett. 63, 373 (1993)－34－

[5] M. Otani and O. Sugino, Phys. Rev. B 73, 115407 (2006).
[6] M. Otani, I. Hamada, O. Sugino, Y. Morikawa, Y. Okamoto and T. Ikeshoji, Phys. Chem.
Chem. Phys. 10, 3609 (2008) and references therein.
[7] Y. Gohda and S. Tsuneyuki, Phys. Rev. Lett. 106, 047201 (2011).
９. 連携研究者・研究協力者
連携研究者：
中山隆史（千葉大学）、杉野 修（東京大学）、森川良忠（大阪大学）、赤木和人（東北大学）、館山佳尚（物質・材料研究機構）、諏訪雄二（日立基礎研究所）、合田義弘（東京大学）

１０. 本研究課題における平成２２年度の発表論文と招待講演
発表論文：
1)

J. Yamauchi and Y. Yoshimoto, ``X-ray Photoelectron Spectroscopy for the Boron Impuritiesin Silicon: a First-Principles Study" Proc. 30th Int. Conf. Physics on Semiconductors(Seoul, July 26 - 30, 2010).

2)

Y. Gohda and S. Tsuneyuki, ``Two-Dimensional Intrinsic Ferromagnetism at Nitride-Boride
Interfaces", Phys. Rev. Lett. 106, 047201 (2011).

招待講演：
1)

常行真司「次世代スパコン戦略機関と第一原理電子状態計算の挑戦的課題」（2011. 9.
10 CMDワークショップ，大阪大学）.

－35－

密度汎関数法理論に基づく非平衡ナノスケール
電気伝導ダイナミクス
Nanoscale Non-Equilibrium Electric Transport Dynamics Based on
Density Functional Theory
渡邉聡 1、渡辺一之 2、相馬聡文 3、小野倫也 4
S. Watanabe1, K. Watanabe2, S. Souma3, T. Ono4
東京大学 1、東京理科大学 2、神戸大学 3、大阪大学 4
1

The University of Tokyo, 2Tokyo University of Science, 3Kobe University,
4
Osaka University

ナノスケール電気伝導現象について、様々な面から計算科学的に研究した結果を報告する。交流応答特性、過渡応答特性、光刺激電子励起といった電子ダイナミ
クスが関係した現象、単分子架橋や原子スイッチ等の原子ダイナミクスに係わる現象、および最近注目を集めているナノリボンの電気特性等に関する解析で興味深い知見を得た他、方法論・プログラム開発にも進展があった。

１． はじめに
ナノスケール電気伝導は、1990 年代から活発に研究されているトピックである。計算科学の面からも、非平衡開放系の電気伝導特性を計算する手法の開発が進み、様々なナノ構造の電気伝導特性が解析されてきた。しかし、界面ラフネスをはじめとした原子配列の詳細、局所高電界場の印加による非平衡過程の出現、高クロック動作に伴う過渡現象の重要性の増大等が絡み合ったダイナミックな過程がナノデバイスやその実験的研究の場において発現していることに鑑みると、我々のナノスケール電気伝導に対する理解はまだ不十分である。そこで本計画班では、様々な意味でのダイナミクスに特に重点を置き、電気伝導とフォノン・熱、イオン伝導、スピン、電子励起等との絡み合いやより実際に近いモデル化にも注意を払いながら、ナノスケール電気伝導のより深い理解を目指して計算科学研究を進めている。以下に、その平成 22 年度における成果を述べる。
２．交流応答特性および過渡応答特性に関する研究
ナノスケール電気伝導に関する従来の計算科学研究は大部分が定常状態に関する物であったが、デバイスの動作においてはスイッチング時の過渡応答や交流電圧印加時の特性が重要である。この点に鑑み、交流応答特性、過渡応答特性の検討を進めた。
交流応答特性については、サブ THz～THz 領域での高速動作を狙う次世代ナノ電子デバイスの材料としてナノカーボン材料が有力であることから、金属カーボンナノチューブ（CNT）のサブ THz 交流電気伝導を解析している。非平衡グリーン関数（NEGF）法と強結合法と組み合わせ、電極の取り扱いにおいてはワイドバンド極限近似を、交流輸送係数であるアドミッタンスの評価には線形応答近似を用いた計算により、交流応答特性の電極‐CNT 接触強度依存性を解析した結果、接触強度が強くなると共に容量性サセプタンスから誘導性サセプタンスへの転移が起きること、そしてこの転移がちょうど量子化コンダクタンス g = 2e2/h で起こっていることがわかった（図１参照）。さらに、
この結論は様々なカイラリティーの金属 CNT に対して成立する普遍的なものである－36－

ことがわかった 1)。また、NEGF 法と密度汎関数法（DFT）計算とを組み合わせたプログ
ラムも開発し、単原子空孔欠陥や Stone-Wales
欠陥を含む CNT の交流応答のシミュレーショ
ンを進めた。結果の詳細は解析中である。
過渡応答特性については、NEGF 法、強結合法、ワイドバンド極限近似および線形応答
近似に基づく計算で量子ドット‐電極結合系にステップ関数的に電圧を印加した際の応答特性について以前に解析した結果[1]を踏まえ、
水素分子またはオクタテトラエン分子が金属電極と結合した系に対して解析した。これら 図 1. ジグザグ型金属ナノチューブ架橋のの系における過渡応答特性は量子ドット系の エミッタンス E と直流コンダクタンス g の場合より複雑であったが、電極と分子との結 関係。E0 は接触強度パラメータγ=1 の場合合が弱い場合には、量子ドット系の場合と定 のエミッタンス。挿入図は E とγの関係。
性的に同じ振舞いを示す独立な成分に分解できることがわかった[2]。
３．光誘起電気伝導に関する研究
ナノスケール物質の光応答と電子励起ダイナミクスを多体電子論と DFT によって明らかにすることを目的として、平成 22 年度は（1）CNT スピン三重項励起子発光メカ
ニズムの解明と（2）炭素系ナノリボンからのレーザー刺激電界電子放射のシミュレー
ションを行った。以下に概要を述べる。
CNT の光物性の中で、励起子からの発光は最も重要なテーマの一つである。最近の実験で、原子状水素が吸着した CNT のフォトルミネッセンスの発光スペクトルに、そのエネルギー位置からスピン三重項励起子からの発光とみられるサブピークが観測された[3]。これは、スピン三重項励起子は光学禁制で発光しないという今までの常識を覆すものである。この観測結果についてベーテ・サルピータ方程式を解いて解析することにより、原子状水素吸着によってCNT 局所電子状態が SP2 から SP3 結合に変わることでスピン軌道相互作用が大きくなり、スピン一重項・三重項励起子間混成によって光学禁制であった暗いスピ
ン三重項励起子が光学許容になることを明らかにした 2)。吸収スペクトルの計算結果（図２）からわかるように、低エネル
ギー位置のスピン三重項励起子ピークが
図２．水素吸着(20,0)CNT の光吸収スペクトル
吸着水素量に比例して成長している 2)。これは実験結果をよく説明している。
レーザー刺激電界電子放射に関しては、最近、フェムト秒パルスレーザーを電界電子源に照射することによって電界放射特性を制御する実験が注目を集めている[4]。そこで、フェムト秒パルスレーザー刺激を受けた炭素系ナノリボンから電界放射する電子のダ
イナミクスを追跡するために、時間依存密度汎関数を使ったシミュレーションに着手した。得られた結果の一つは、グラフェンリボンに 1V/Åの静電界と 0.1V/Åの強度で 4eVの振動数のレーザーを照射すると、レーザーを照射しない場合に比較して 2 倍の電流が放射されたことである。これは、Optical Field Emission と定義される特性と理解できる。

－37－

４．原子ダイナミクスと電気伝導との相関に関する研究
熱振動や電圧・電場による原子移動は電気伝導特性にも大きな影響を及ぼす。電子‐
フォノン散乱のナノスケール電気伝導への影響については既に理論研究が進んでいるが、本計画班ではこれに限らず原子ダイナミクスとナノスケール電気伝導との相関について広く研究を進めている。ここでは、単分子架橋、原子スイッチ、4 端子抵抗測定の3 つのトピックについて成果を報告する。
金属電極間の単分子架橋の電気伝導特性は、ナノスケール電気伝導の中でも理論・実験両面から特に盛んに研究されているトピックの一つである。しかし、室温・溶液中の分子架橋については、実験が盛んであるにもかかわらず理論解析は少なかった。そこで、本班メンバーは以前に第一原
理分子動力学法による原子配置計算と NEGF 法による電気伝導計算とを組み合わせた解析を行った[5]が、その方法論は整合性が不十分であり、またバイアス電圧印加下での分子動力学（MD）計算もできなかった。そこで、NEGF 法の枠組みの中で MD 計算も行って、金電極間ベンゼンジチオール架橋の電気伝導度と構造の時間変化を計算した[6]。伝導度ヒストグラムの計算結果を図３に示す。バイアス電圧ゼロの場合については、伝導度揺らぎ幅の溶媒水分子の有無による変化や電気伝導度と C‐S 結合長の負の相関といった、先行研究[5]に現れた重要な特徴をおおむね定性的に再現する結果が得られた。バイアス電圧を 1V 印 図３．( a ) 水中および( b ) 真空中の単加した場合には、これらの特徴は変わらないもの 分子架橋の伝導度ヒストグラム。低伝の、伝導度の平均値が大きくなった。さらに、伝 導度ピークはバイアス電圧 0 V 、高伝
導度および C‐S 結合長のフーリエスペクトルの 導度ピークは 1 V の場合。
解析を進めている[6]。伝導度のフーリエスペクト
ルの形状に水中の場合と真空中の場合とで明らかに異なること、スペクトルのピークが、
C‐S 結合に関する対称振動モードと対応すること等を見出した。
原子スイッチ（固体電解質を酸化可能な金属電極と不活性な金属電極の間に挟んだ構造を持ち、電圧の印加により低抵抗状態と高抵抗状態との間をスイッチングする素子）に関しては、Cu‐Ta2O5‐Pt 接合系関する DFT 計算により Ta2O5 層中に添加 Cu 原子が 1列に並んで電極間を架橋した構造が高い伝導性を示すことを既に明らかにしている[7]。
これを踏まえ、低抵抗状態の形成過程を解明するための予備検討として、Cu‐Ta2O5 接合系に対して有効遮蔽媒体法[8]を用いた電場印加計算を行った[9]。その結果、電場印加
時にも Ta2O5‐Cu 界面近傍の実効的な電場は非常に小さいことを見出した。これは、低抵抗状態の形成過程がバイアス電圧印加でまず Cu 電極中の界面付近の Cu が動き出すという単純なものではないことを示唆している。
4 端子抵抗測定のシミュレーションに関しては、開発を進めてきた多端子電気伝導シ
ミュレータ[10]に、有限バイアス電圧印加の効果および電子‐フォノン散乱の効果を陽に取り込むことを試みた。有限バイアス電圧に関しては、電圧プローブに流れる電流がゼロとなるような電圧プローブ電圧を透過関数スペクトルがバイアス電圧により変化しないとの仮定の下で求める関係式を導出することができた他、有限バイアス電圧による
静電ポテンシャル変化と透過関数スペクトル変化を陽に考慮した計算にも炭素原子鎖系で成功し、上述の関係式の妥当性を確認することができた[11]。電子‐フォノン散乱については、試料、プローブとも水素原子鎖から成る系について中央の 1 原子の振動のみ－38－

考慮した計算に成功し、振動エネルギー（42meV）に対応した電流プローブ間バイアス
電圧において 4 端子抵抗値が増大することの他、電流プローブ間電圧が 75meV のとこ
ろでも 4 端子抵抗値が増大することを見出した[12]。さらに、75meV での抵抗値増大は
電流プローブ‐電圧プローブ間の電圧が振動エネルギーと等しくなるためであることを明らかにした。以上に加え、多端子電気伝導シミュレータの高速化のための計算アルゴリズムについて、A01-3 張紹良班と議論した。現在は克服すべき課題が明らかになった段階であるが、今後連携を進めることにより、高速化が進む可能性がある。
５．ナノリボン構造の電気伝導特性
グラフェンは、シリコンに代わる半導体材料として近年注目を集めている材料であり、
特に，それをリボン状に加工したグラフェンナノリボンは、結晶に対するリボン軸の方向やリボンの幅によって異なる電気的、光学的性質を持つ事から応用上の興味が高まっている。これに関連して、BN ないし BNC 複合ナノリボン構造も興味を持たれている。本計画
班では、神戸大と阪大の 2 つのサブグループが、異なる面からナノリボンの電気特性に関する研究を進めた。
神戸大グループは、ジグザグ端グラフェンナノリボンによる n-i-n 構造（n:電子ドープソー
ス／ドレイン電極，i：ノンドープチャネル）の特性を解析した。グラフェンナノリボンの
うちアームチェア端のリボンはその幅によって金属的な場合と半導体的な場合に分かれる
が、ジグザグ端を持ったリボンは幅によらず端状態と呼ばれる端に局在した状態がフェルミ面近傍に存在し、そのため金属的な性質を持つと言われている。しかし実はジグザグ端
グラフェンナノリボンによる n-i-n 構造（n:電子ドープソース／ドレイン電極，i：ノンドー
プチャネル）は、その幅が偶数個のジグザグ状炭素原子列から構成される場合と奇数個のジグザグ状炭素原子列から構成される場合とで異なった電流電圧特性を持つことがわかっている。特に幅が偶数個のジグザグ状炭素原子列から構成される場合、価電子帯頂上のバンドと伝導帯底のバンドがそれぞれ幅方向に反対称、対称な波動関数を持つことから、それらの間での電子の透過がブロックされ、通常の電界効果型トランジスタに類似した電流の飽和特性が現れる[13]。
この現象は対称性に起因した現象であるため、格子欠陥が存在する場合にはその影響を強く受ける事が予想される。ところが、リボン端に原子欠陥がある場合のシミュレー
ションを行ったところ、欠陥がドレイン電極近くにある場合にはたしかにその影響が強く現れたのに対し、ソース電極近く 図４．偶数幅ジグザグ端グラフェンナノリボンの電流電圧特性ににある場合にはほとんど おける格子欠陥の影響。欠陥位置が異なる３つの場合の比較。
影響が現れない事がわかった 3)。
阪大グループは、まず方法論の改良に取り組んだ。グリーン関数を用いて輸送特性を計算する従来の Overbridging boundary-matching 法は、グリーン関数を計算するためにほとんどの計算時間を費やしていた。そこで、グリーン関数を計算する代わりに、入射波に対する連立方程式を解くことにより散乱の波動関数を計算する方法を開発した 4)。
この方法は、Jellium 電極を用いることを強いられるが、従来法に比べ 5 倍以上の計算
速度で輸送特性を計算することができる。この方法を用いて、BNC ナノリボン構造の－39－

スピン輸送特性を調べた。
2 次元 BNC
構造において、BN(Boron-Nitride)領域と C(Carbon)領域を分けるストライ
プ構造を導入すると、C 領域がスピン分極することが報告されている。
しかしながら、2 次元 BNC 構造の両端にグラフェンを接合した際に、ス
ピン分極を示すのか、あるいは、電流を流した際にスピン偏流が得られるのかについては、わかっていない。
電極で挟んだ状況でもスピン分極を維持し、より多くのスピン偏流を得ることは、スピン伝導性を用いた応 図５．グラフェンに接合された BNC 構造(挿入用的な側面からも重要であると考え 図)と入射電子のエネルギーに対する輸送特性。
られる。そこで、両端をグラフェン 挿入図では、白丸が C、灰丸が B、黒丸が N 原に接合されたストライプ構造をもつ 子を示す。
BNC 構造について、電子構造と輸送特性の計算を行い、ストライプ構造が BNC 構造のスピン偏流に与える影響について調
べた。ストライプ構造の変化や C 領域間距離の変化がスピン分極に与える影響を調べたところ、より小さな C 領域が BN 領域に挟まれて孤立した構造を持つ程、C 領域でのスピン分極は大きくなることがわかった。また輸送特性は、図５に示すようにハーフメ
タリック金属のような完全スピン偏極輸送特性を示すことを発見した。
６．参考文献[1] K. Sasaoka, T. Yamamoto, and S. Watanabe, Appl. Phys. Lett. 96, 102105 (2010).
[2] W. Liu, K. Sasaoka, T. Yamamoto, and S. Watanabe, submitted to J. Appl. Phys.
[3] K. Nagatsu, S. Chiashi, S. Konabe, and Y. Homma, Phys. Rev. Lett. 105, 157403 (2010).
[4] H. Yanagisawa, et al., Phys. Rev. Lett. 103, 257603 (2009).
[5] A. Tawara, T.Tada, and S. Watanabe, Phys. Rev. B 80, 073409 (2009).
[6] 稲塚大氣、東京大学大学院工学系研究科マテリアル工学専攻修士論文(2011).
[7] T. K. Gu, T. Tada, and S. Watanabe, ACS Nano 4, 6477 (2010).
[8] M. Otani and O. Sugino, Phys. Rev. B 7 3, 115407 (2006).
[9] 松本拓也、東京大学工学部マテリアル工学科卒業論文(2011).
[10] A. Terasawa, T. Tada, and S. Watanabe, Phys. Rev. B 79, 195436 (2009).
[11] 寺澤麻子、東京大学大学院工学系研究科マテリアル工学専攻博士論文(2011).
[12] 飛松啓司、東京大学大学院工学系研究科マテリアル工学専攻修士論文(2011).
[13] S. Souma, M. Ogawa, T. Yamamoto, and K. Watanabe, J. Comp. Electronics 7, 390(2008).
７．連携研究者・研究協力者
連携研究者：
多田朋史（東京大学大学院工学系研究科特任講師）、山本貴博（東京大学工学系研究科助教（平成23年4月より東京理科大学工学部講師））、小鍋哲（筑波大学数理科学研究科研究員）、酒井明（京都大学工学研究科教授）－40－

８．本研究課題における平成２２年度の発表論文と招待講演
発表論文：
1)

T. Yamamoto, K. Sasaoka, and S. Watanabe, “Universal transition between inductive andcapacitive admittance of metallic single-walled carbon nanotubes”, Phys. Rev. B 82, 205404(2010).

2)

S. Konabe and K. Watanabe, “Mechanism for optical activation of dark spin-triplet excitonsin hydrogenated single-walled carbon nanotubes”, Phys. Rev. B 83, 045407 (2011).

3)

S. Souma, M. Ogawa, T. Yamamoto, and K. Watanabe, “Parity induced edge-current
saturation and current distribution in zigzag-edged graphene nano-ribbon devices”, acceptedfor publication in J. Comp. Electronics (2011).

4)

Y. Egami, K. Hirose, and T. Ono, “Time-saving first-principles calculation method forelectron transport between jellium electrodes”, Phys. Rev. E 82, 056706 (2010).

5)

T. Ono, M. Heide, N. Atodiresei, P. Baumeister, S. Tsukamoto, and S. Blügel, “Real-spaceelectronic-structure calculations with full-potential all-electron precision for transitionmetals”, Phys. Rev. B 82, 205115 (2010).

6)

T. Ono and S. Saito, “First-Principles Study on Electronic Structure of Dangling Bond at
Ge/GeO2 Interfaces”, Appl. Phys. Express 4, 021303 (2011).

7)

A. Terasawa, K. Tobimatsu, T. Tada, T. Yamamoto, and S. Watanabe, “Effects of Resonant
Scattering by Probe Contacts on Nanoscale Four-Probe Resistance Measurements”, New. J.
Phys. 12, 083107 (2010).

8)

寺澤麻子，飛松啓司，多田朋史，山本貴博, 渡邉聡, “ナノスケールの4端子抵抗測定
における量子干渉効果の理論研究”, 表面科学 31, 374 (2010).

9)

S. Kasamatsu, T. Tada, and S. Watanabe, “Theoretical analysis of space charge layer
formation at metal/ionic conductor interfaces”, Solid State Ionics 183, 20 (2011).

10) 渡邉聡, 多田朋史, 山本貴博, “ナノスケール電気伝導ダイナミクスのシミュレーショ
ン”, ナノ学会会報 9 (2011), 印刷中.
招待講演：（下線を付したのは登壇者）
1)

S. Watanabe, T. Tada, S. Kasamatsu and T.K. Gu, “Nanoscale simulations on electronic andatomic conduction in solid electrolyte/metal junction systems”, IUMRS 11th International
Conference in Asia (Qingdao, China, September 25-28, 2010).

2)

S. Watanabe, W. Liu, D. Hirai, K. Sasaoka, and T. Yamamoto, “Simulations on time-varying
nanoscale electronic transport”, Asian Consortium for Computational Materials Science thethird Working Group Meeting (ACCMS-WGM3) on Advances in Nano-device Simulation(Jeju Island, Korea, March 31 - April 2, 2011).

3)

T. Yamamoto, K. Sasaoka, and S. Watanabe, “Coherent AC Transport in Metallic
Carbon-Nanotube Bridges”, Workshop on Simulation and Modeling of Emerging Electronics(SMEE 2010) (Hong Kong, China, December 6-10, 2010).

4)

T. Ono, “First-Principles Study on Transport Property of Nanostructures”, International
Conference on Core Research and Engineering Science of Advanced Materials & Third
International Conference on Nanospintronics Design and Realization (3rd-ICNDR) (Osaka,
May 30- June 4, 2010).

5)

小野倫也, “第一原理計算によるSiおよびGe MOS界面原子構造とリーク電流特性の予測”, ゲートスタック研究会-材料・プロセス・評価の物理, (Tokyo, January 21-23, 2011).

－41－

プロトン・ミューオンで探る新物性と量子ダイナミクス
New properties of materials probed by proton and muon, and their quantum dynamics
中西寛 1、後藤英和 1、下司雅章 1、Markus Wilde2、Wilson Dino1、福谷克之 2*、笠井秀明 1*
H. Nakanishi 1, H. Goto 1, M. Geshi 1, M. Wilde 2, W. Dino 1,
K. Fukutani 2*, H. Kasai 1*
大阪大学 1、東京大学 2Osaka University 1, The University of Tokyo 2
本研究班は、物質環境下におけるプロトン・ミューオン等の粒子の振る舞いに関する理論的取り扱い方法、およびその第一原理計算コードを開発し、その量子ダイナミクスを探る。また、それら粒子の関わる新規物性を探査する。第 1 章では、物質環境下におけるプロトン・ミューオンの第一原理計算手法について、第 2 章では、プロトン・ミュ
ーオン計算手法に供する為、より簡便で高精度な電子状態計算手法の開発状況について報告する。 新規物性の探査として、第 3 章で高圧下の金属水素化物への取り組み、第4 章で金属表面での水素吸収ダイナミクスと触媒反応性に関する実験グループの取り組
みについて報告する。なお、実験グループで見いだされた結果は、理論グループの開発するシミュレーションコードの評価にも用いる予定である。

第１章「物質環境下におけるプロトン・ミューオンの第一原理計算手法の開発」
物質環境下における水素等の質量の小さな原子の核の運動に対して量子力学を適応する第一原理量子ダイナミクスコード(Naniwa-series)を開発している[1-6]。今年度は、
Naniwa-Static[3-6]を用いて Pd(111)表面における水素の同位体（μ＋、H、D）効果を調査した。その結果、μ＋は H、D に比べ全く異なる量子状態を Pd (111) 表面上で示すことがわかった。
1.1 はじめに密度汎関数理論に基づく電子状態計算法は、広義の局所密度近似を用いることにより現在の計算機の能力で演算量・データ量を妥当な範囲に収めることができ、様々な物質
系に適応された、最も成功した固体物性論の計算手法の一つである。さらに原子と原子の間の相互作用にこの計算手法を援用した第一原理分子動力学法は、様々な物質の動的過程（化学反応を含む）に適用され成功を収めつつある。しかしながら小さな質量の水素原子の振る舞いはトンネル効果、干渉効果、束縛状態のエネルギー離散化効果（零点運動を含む）、非局在効果等の量子力学的効果が顕著になり、明らかに分子動力学法の適応範囲外にある。
また、近年ミューオン（正確には反ミューオン：μ＋）が物性のプローブ粒子として活用されている。μ＋は、陽子と同じ電荷とスピンをもつ安定な素粒子である。質量は、陽子の１／９で、物質中では水素の同位体としてふるまう。μ＋のスピンは、パイ中間
子の自然崩壊で生成される際、運動量とは反対方向に偏極しており、物質中に放たれたμ＋のスピンは物質内の局所的な磁場を受けて歳差運動を始める。μ＋の平均寿命は 2.2
μ秒と十分長いため、物質中の安定な位置に落ち着いた後に崩壊する。その際スピン－42－方向に選択的に放出される陽電子を観測することにより、μ＋の滞留した位置での局所的な磁場を評価することができる（μSR）。
我々はこれまで電子と同じく水素原子核にも量子力学を適用するための計算方法を模索してきた。その中で、様々な元素からなる様々な構造の固体表面に対して一様の近似精度で評価することが可能な第一原理（電子状態）計算手法の現在の利点を生かしつ
つ、その固体表面における水素の量子力学的振る舞いを記述する方法として電子系－水素原子核－環境格子系の二段階の断熱近似を用いる手法を実践してきた。 この方法で、金属表面上の吸着水素原子の基底状態の運動量分布や、振動励起エネルギーが、定量的に実験結果（前者は共鳴核反応法、後者は電子線エネルギー損失分光法）と一致することも示してきた。
今年度は、表面での水素の量子状態を計算する
Naniwa-Static を用いて、Pd (111) 表面におけるμ＋の吸着量子状態を調査し、H および D の吸着量子状態[6]と比較検討を行った。
1.2 第一原理量子状態計算コード：Naniwa-Static
原子・分子の吸着量子状態を調べるため、まず固体表面の上に、吸着分子・原子を置き、原子・分子の位置座標の関数として系のトータルエネル
ギー（原子・分子の運動エネルギーを除く）を、密度汎関数法をもとにした第一原理計算によりもとめ、固体表面で原子・分子が受ける断熱ポテン
シャルエネルギー超曲面を導出する。得られた断熱ポテンシャル超曲面で原子・分子のシュレーデ
ィンガー方程式を解くことにより吸着量子状態の固有エネルギーと波動関数を得ることができる[3-6]。
1.3 計算結果
図 1.1(a)は、Naniwa で計算した Pd(111)表面上でのμ ＋の基底状態の波動関数である。図1(b)に示す軽水素 H(の原子核 proton)の波動関数が fcc hollow サイトに局在しているのに対し、μ＋は fcc hollow サイトでの頻度が高い
ことは変わりないが top サイトを除く表面全
体に広がり、局在性が非常に弱いことがわかる。図 1.2 にμ＋、H、D(の原子核 deuteron)の各エネルギー準位を示す。H と D の零点振動の基底状態および低エネルギーの励起状態
では、波動関数の対称性は変化せず、エネル
ギー準位のシフトと局在性の若干の強弱が生じていたが、μ＋では、基底状態から著しく局在性が失われる為、低エネルギー準位の対称性も異なる。H、D では、基底状態と第一励
起状態は僅かにエネルギー分裂をおこし、それぞれ fcc hollow サイトと hcp hollow サイトに局在したものであった。そこから数百 meVのギャップを挟んだ第二励起状態以上で振動－43－

図 1.1. Pd(111)表面上での反ミューオン(a)と軽水素(b)の基底状態の波動関数。

図 1.2. Pd(111)表面上での反ミューオン、軽水素、重水素のエネルギー準位(a)。

図(a)内シアン色

枠内の拡大図(b)。各準位右上の数字は、縮重度若
しくは僅かに分裂した近接準位数。

励起モードが現れていた。それに対しμ＋では、基底状態と第一励起状態の間で既に数十 meV のエネルギー差がある。また、第一励起状態から振動励起モードを表すノード
が波動関数に生じ始める。このように Pd(111)表面でのμ＋と H、D では、定性的に異なる様相を示すことがわかった。
1.4 まとめ
Pd(111)表面上での反ミューオンの吸着量子状態を調査し、水素原子(H、D)と比較を行った。以前に調査した水素(H、D)の基底状態および低エネルギー励起状態への同位体効果とは、質的に異なる変化が見いだされた。現在のところ、ミューオンビームの強度
は表面系の探査のためには弱く、これを直接実験で検証することは現時点では難しい。
そこで、実験がなされているバルク中での振る舞いも取り扱えるように計算コードを拡張し、調査を行うことを計画している。
1.5

参考文献

[1] H. Kasai, A.Okiji, Progress in Surface Science, 44 (1993) 101.
[2] W.A. Diño, H. Kasai, A. Okiji, Progress in Surface Science 63 (2000) 63.
[3] K. Nobuhara, H. Kasai, H. Nakanishi, A. Okiji, Surface Science, 507 (2002) 82.
[4] N. Ozawa, N. B. Arboleda Jr., H. Nakanishi, N. Shimoji, H. Kasai, Surface and Interface Analysis, 40(2008) 1108.
[5] N. Ozawa, M. Sakaue, H. Kasai, Journal of Vacuum Society of Japan, 53 (2010) 592.
[6] N. Ozawa, T. Roman, N. B. Arboleda Jr., W. A. Diño, H. Nakanishi, H. Kasai, Journal of Physics:
Condensed Matter 19 (2007) 365214.

第２章「多体系量子状態計算手法の開発」
多体効果を高精度に取り入れることが可能な計算手法によるプロトン・ミューオン系の量子状態シミュレーションを実現することを目的として、多電子状態の高精度・高効率計算手法の開発を行った。多電子波動関数としてスレーター行列式の線形結合を用い、変分原理に基づいて 1 電子波動関数の更新を繰り返すことで非直交基底関数系を生成し、できるだけ少ない数のスレーター行列式で基底状態へと収束させる方法を提案した。
コード開発を行い、ガウス関数基底セット(6-31G*)を用いて炭素原子の基底エネルギー計算に適用したところ、厳密解を与える full CI (Configuration Interaction) 法の 1/2000 分の数の Slater 行列式で、誤差 0.003%の基底エネルギーが得られた。また、系の電子数の増加に対し、必要な Slater 行列式の数は full CI 法に比較して穏やかに増加する傾向
であることがわかった。
2.1 はじめに
当該計画研究では、プロトンやミューオンの量子トンネル効果や干渉効果などの量子力学的ダイナミクスを理論的・実験的に解明し、新規物性を探索することを目的としている。電子と原子のみから成る系の量子状態計算においては、密度汎関数理論による電子状態シミュレーションが現在最も用いられているツールである。この方法は数多くの成功を収めている信頼性かつ実用性の高い手法ではあるが、交換・相関エネルギー項には近似を用いており、計算対象によっては定量的・定性的精度に限界があることが多くの実例により指摘されていることもまた事実である。加えて本プロジェクトが対象とするプロトンやミューオンを含む系の量子ダイナミクスシミュレーションにそのまま用いることはできない。以上のような背景から、本分担研究では密度汎関数理論に依らずに電子相関エネルギーを高精度かつ高効率に計算することが可能な手法の開発を目的－44－

とした研究を行う。さらには、開発した手法とこれまでに研究・開発を行ってきたイン
パルス・レスポンス法 [1-3] と呼ばれる電子輸送特性シミュレーション手法を基に、電子相関を正確に取り入れた電荷・スピン輸送特性シミュレーションを行うことも計画している。
そこで、本年度の研究では電子相関エネルギーの高精度・高効率計算手法の提案と開発を行った。分担者はこれまでの研究で、基底関数を用いずに空間に一定の間隔で設定
したグリッド点上の物理量のみを扱う実空間差分法 [4, 5] に基づき、多電子状態計算手法である Direct Energy Minimization (DEM) 法 [6-12] の提案と開発を行ってきた。この手法では、多体波動関数を Slater 行列式の線形結合で表わし、変分原理に基づき各 1
電子波動関数のグリッド点上の値を基底状態に向かって繰り返し更新してゆく。重要なポイントは、波動関数には直交・規格化条件を課さずに更新を行うことで、非直交な 1
電子波動関数系を生成する点である。非直交基底系を用いることで、1 個の Slater 行列式で複数の配置関数を効果的に取り込むことが可能であり、直交基底を用いて全ての配置関数を用いる配置間相互作用 (CI: Configuration Interaction) 法よりも少ない数のSlater 行列式で基底状態を表現できることが指摘されている [13-16]。本年度の研究では、この手法における計算時間の短縮を目的としてガウス関数基底セット [17] を導入した。ガウス関数基底セットを用いれば、グリッド点を利用した DEM 法 (Grid DEM) で多大な時間を要した数値積分にガウス関数の積分公式を使用することが可能となり、大幅な計算時間の短縮が可能となる。さらに、1 電子波動関数の更新方法についても効率の良い方法を提案した。シミュレーションコード作成を行うとともに簡単な原子・分子に対して適用し、主に基底状態への収束性と必要な Slater 行列式の数に着目し本手法の性能評価を行った。
次節以下では、今回提案した計算手法の概略を説明した後に計算結果例を紹介し、最後にまとめを述べる。
2.2 計算手法の概要
N 電子系の波動関数 Ψ (τ 1 ,τ 2 , Λ , τ N ) を次式のように Slater 行列式の線形結合で表わす。

LΨ (τ 1 ,τ 2 ,Λ ,τ N ) = ∑ C A
A=1

ψ 1A (τ 1 ) ψ 2A (τ 1 ) Κψ 1A (τ 2 ) ψ 2A (τ 2 ) Κψ NA (τ 1 )ψ NA (τ 2 )ΜΜΟΜψ (τ N ) ψ 2A (τ N ) Κ ψ NA (τ N )

(2.1)

A
1

L
L
ρ ρ ρ
ρ
≡ ∑ C A ψ 1A ψ 2A ψ 3A Λ ψ NA ≡ ∑ C A Φ A
A=1

A=1

ここに ψ (τ j ) ≡ φ i (r j )γ i (σ j ) は 1 電子波動関数であり、 r j , σ j はそれぞれ j 番目の電子の位置座標とスピン座標、 φiA (r ), γ i (σ ) は、それぞれ i 番目の 1 電子波動関数の空間部分とスピン部分である。
1 電子波動関数の空間部分 φiA (r ) は、次式のようにガウス基底関数 χ s (r ) の線形結合
Aで表わす。ここに Di , s は展開係数、M は基底関数の数である。
A
i

A

M

φiA ( ri ) = ∑ DiA, s χ s ( ri )

(2.2)

s =1

ここで、1 電子波動関数を基底状態へ向かって更新するための修正関数
N baseξ µ (r ) = ∑ Gµ , s χ s (r )

(2.3)

s =1

を N C 個導入し、次式のように A 番目の Slater 行列式の i 番目の 1 電子軌道関数を修正

－45－

する。

φi

A( new)( r ) = CiφiA( old )NC( r ) + ∑ CL + µ ξ µ ( r )

(2.4)

µ =1

この時、N 電子波動関数は次式のように L + N C 個のスレーター行列式の線形結合で表わされる。

Ψ (r1 , r2 ,Λ , rN )( new )

L
ρ
ρ
ρ
ρ
ρ
ρ
ρ ρ ρ
ρ
= C L +1 φ1A Λ ξ1A Λ φNA + Λ + C L + N C φ1A Λ ξ NAC Λ φNA + ∑ Ck φ1k φ2k φ3k Λ φNk

(2.5)

k =1

NC +L

≡ ∑ Ck Φ k
k =1

係数 C i は、全エネルギーE に変分原理を適用することで得られる一般化固
有値方程式により求めることができ、A 番目の Slater 行列式の i 番目の 1 電子
波動関数は式(2.4)により修正される。
この修正を全ての Slater 行列式の全ての 1 電子波動関数について行うことで1 回の更新作業が終了する。そして、この更新を繰り返すことで基底状態を表わす互いに非直交な 1 電子波動関数が形成される。実際の計算では、最初に Hartree-Fock (HF) 方程式を解き、得られた HF 解を初期波動関数として採用している。
2.3 計算結果例
炭素原子の計算例を図 2.1 に示す。
厳密解を与える full CI (Configuration
Interaction) 法の 1/2000 分の数(100 個以下)のスレーター行列式で、誤差0.003%の基底エネルギーが得られた。
今回提案した方法で作成された非直交基底が、効率よく電子励起配置を取り込んでいるものと思われる。
以上のような計算を種々の原子に対して行い、基底エネルギーを得るた
めに必要であった Slater 行列式の数を
図 2.2 に示した。Full CI 法では原子番号の増加とともに Slater 行列式の数が爆発的に増大するが、本手法では激しい増加は見られないことがわかる。図2.3 はいくつかの簡単な分子についての結果を示しており、図 2.2 と同様穏やかな増加傾向を示している。
次に、1 電子波動関数の更新プロセ
スにおけるパラメーターである修正

図 2.1
関係

C 原子における Slater 行列式の数とエネルギーとのガウス基底セット：6-31G*

修正関数の数 Nc=8

図 2.2 各種原子における基底状態に必要な Slater 行列式の数

図 2.3 各種分子における基底状態に必要な Slater 行列式の数－46－

関数の数 Nc (式(2.4)参照) と全エネ
ルギーの収束性との関係について検討した。ガウス基底セットを
6-31G*として Be 原子の基底状態を求めた結果を図 2.4 に示す。修正関数の数 Nc が多いほど探索空間の次元が大きくなるため、基底状態への収束性が向上することがわかる。
2.4 まとめ
非直交基底を用いて多電子系の基底状態を高精度かつ高効率に求図 2.4 Be 原子における更新回数とエネルギーとの関係
める手法を提案した。また、コード作成を行うとともに簡単な元素や分子に適用し性能を評価した。一例として 6-31G*ガウス関数基底セットを用いて炭素原子の基底エネルギー計算に適用した場合、厳密解を与える full CI (Configuration
Interaction) 法の 1/2000 分の数のスレーター行列式で、誤差 0.003%の基底エネルギーが得られた。また、基底状態に必要な Slater 行列式の数は、原子番号の増加に対して full CI法のように激しく増加しないという優れた特徴が見出された。以上のように、非直交基底による配置関数の効果的な取り込みが実証できたが、本手法は繰り返し計算が必要であり、計算時間については実用的なレベルにまで至っていない。今後は非直交基底の作成方法や 1 電子波動関数の更新方法も含めて、より短時間で基底状態に到達するための工夫が必要である。さらには、周期境界条件やプロトン・ミューオンを含む系への適用試験を行っていきたい。
2.5

参考文献

[1] H. Goto, T. Ono and K. Hirose, J. Phys.: Condens. Matter, 19, 365205 (2007).
[2] T. Suzuki, H. Goto and K. Hirose, Surf. Interface Anal. 40, 1113 (2008).
[3] H. Goto, S. Aiba, T. Suzuki and K. Hirose, J. Comput. Theor. Nanosci. 6, 2656 (2009)[4] J. R. Chelikowsky, N. Troullier, K. Wu and Y. Saad, Phys. Rev. B 50, 11355 (1994)[5] K. Hirose, T. Ono, Y. Fujimoto and S. Tsukamoto, First-Principles Calculations in Real-SpaceFormalism, Imperial College Press, London (2005)[6] H. Goto and K. Hirose, J. Phys.: Condens. Matter 21, 064231 (2009)[7] F. Mauri, G. Galli and R. Car, Phys. Rev. B 47, 9973 (1993)[8] F. Mauri and G. Galli, Phys. Rev. B 50, 4316 (1994)[9] K. Hirose and T. Ono, Phys. Rev. B 64, 085105 (2001)[10] T. Sasaki, T. Ono and K. Hirose, Phys. Rev. E 74, 056704 (2006)[11] H. Goto, T. Yamashiki, S. Saito and K. Hirose, J. Comput. Theor. Nanosci. 6, 2576 (2009)[12] H. Goto and K. Hirose, J. Nanosci. Nanotechnol. 11, 2997 (2011)[13] H. Fukutome, Prog. Theor. Phys. 80, 417 (1988)[14] Y. Noda and M. Imada, Phys. Rev. Lett., 89, 176803 (2002)[15] T. Kashima and M. Imada, J. Phys. Soc. Jpn. 70, 2287 (2001)[16] 渡辺真仁、水崎高浩、今田正俊、固体物理、39, 1 (2004)[17] 例えば、A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced
Electronic Structure Theory, Macmillan, London (1982)－47－

第３章「高圧下における金属水素化物にみる水素効果」
3.1 はじめに長年の電子ガス及び金属水素の理論的研究に関連して、近年水素を組成として多く含む金属化合物の物性に興味が持たれている。その理由の一つは、圧力をかけた時金属元素があることで水素が測定値よりも高圧状態になっていると考えられるためである。水素は 400GPa くらいの圧力でようやく金属化すると理論的には予測されているが、それに相当する状態がずっと低圧で実現していると考えられ、金属化や超伝導状態の出現、しかもその転移温度が高いことが期待されている。Ashcroft によるこのような提案[1]がなされてすぐに理論計算が行われ、例えば Yao らは SiH4 では 125GPa での超伝導転移
温度は 50K 程度という高いものになると予測された[2]。そのほかにも CH4、GeH4、SnH4
等について理論的予測が報告されている。
（これらを MH4(M は金属元素)と表す。
）SiH4については実験が行われたが、最高 17K であって、理論予測に比べてかなり低い値である。他にも MH3 という組成の金属化合物についても理論予測が行われている。これらは MH4 に比べて超伝導転移温度は低いと予想されているが、MH4 よりも低圧で出現すると予想されている。
例えば Goncharenko らは、理論予測では 110GPa で 24K とした。
彼らで実験も行ったが、160GPa を超えても超伝導状態は確認されなかった[3]。理論予
測は電子格子相互作用を起源とする BCS 理論の範囲で、Allen-Dynes の公式で見積もっているものばかりである。この方法では考慮されない効果が起因しているのかもしれな
いが、理論予測と実験との不一致の原因はまだ分かっていない。
3.2 今年度の準備と今後の展望
今年度は、このような高圧下での金属水化物の研究の現状を調査し、これらのいくつかの計算を行うことにより、高圧下での新しい現象の可能性の探索を行った。また、いくつかの新規物質の力学的構造安定性を弾性定数の計算やフォノン計算により行った。
水素については比較的低圧で水素自身はより高圧での状態と同等の状態にあることは確かで、実験可能な圧力範囲で興味深い物性が現れる可能性は十分にある。ただ、先行研究では金属化や超伝導に焦点が当てられているが、それらは高圧下での興味深い現象の一例である。Ashcroft の提案に啓発された研究は、金属水素など水素に興味が偏った面があるので、金属元素と水素との相互作用での興味深い物性についても注目していく。
3.3

参考文献

[1] N. W. Ashcroft, Phys. Rev. Lett. 92, 187002 (2004).
[2] Y. Yao, J. S. Tse, Y. Ma, and K. Tanaka, Eur. Phys. Lett. 78, 37003(2007).
[3] I. Goncharenko, M. I. Elemets, M. Hanfland, J. S. Tse, M. Amboage, Y. Yao, and I. A. Trojan, Phys.
Rev. Lett. 100, 045504 (2008).

第４章「Pd(110)単結晶表面における水素吸収ダイナミクスと触媒反応性」
The hydrogen absorption dynamics at a Pd(110) single crystal surface are investigated by
combined LEED, NRA, and isotopic labeled TDS measurements. We consistently identify allfeatures in the complex H2 TDS spectrum of Pd(110) as surface-adsorbed and subsurfaceabsorbed H-species, respectively. A distinction of two forms of absorbed hydrogen in Pd(110) (aconcentrated near-surface hydride phase (α1) and bulk-dissolved H atoms (α3)) is achieved forthe first time. The kinetics of the α1, α3 state populations in the 85-160 K range reveal unusuallysmall activation energies, which are at odds with the conventional concept that hydrogen－48－

absorption proceeds as a surface-to-subsurface transition of chemisorbed H atoms. Isotopic
exchange between absorbed and surface-adsorbed hydrogen occurs only for the bulk-absorbed
α3 species, but not for the α1 hydride, suggesting that more than one microscopic pathway exists
for H atoms between gas phase H2 and the Pd metal interior. The reactivity of the Pd-catalyzedhydrogenation of butene (C4H8) to butane (C4H10) is found to proceed only with α3-hydrogen,but not with α1, which furthers our suspect that the recombination/resurfacing of these twoabsorbed H species may occur on different pathways.
4.1

はじめに
Hydrogen absorption by palladium is of paramount importance for industrial catalysis ofolefin hydrogenation, in fuel cells, and hydrogen storage applications. The exact pathway for thereversible transition of H atoms between molecular H2 in the gas phase and atomic H in themetal interior however, remains elusive, despite extensive research on the hydrogen interaction
with Pd surfaces [1]. The Pd(110) surface is a prototype open surface known to be particularly
prone to hydrogen absorption, for which its reactivity is enhanced by a H-induced pairing-row
reconstruction upon contact with hydrogen [2]. A variety of ordered H adsorption phases has
been identified [3] but the correlation to the rather complex TDS spectrum of Pd(110), whichfeatures multiple H2 desorption peaks of debated origin [4], is not yet clear. We therefore aim
here to achieve a clear assignment of the numerous TDS signals to the surface adsorption phasesof hydrogen, and especially to discriminate these features from those originating from hydrogen
absorbed below the surface. To this end we complement TDS measurements by LEED and NRA
data acquisitions, to obtain information on the surface structures and possible subsurface
absorption. The temperature range of formation, thermal stability and depth extension of the‘subsurface’ hydrogen state are of further interest in this study.
Having identified the TDS trace of ‘subsurface’ hydrogen, our second step is to clarify
details of the microscopic hydrogen absorption mechanism. The conventional concept ofhydrogen absorption regards this process as a surface-to-subsurface transition of chemisorbed H
atoms. Since the surface adsorption energy is profoundly more exothermic (-0.5 eV/H) than theheat of H solution in bulk Pd (-0.1 eV/H) and subsurface sites are only slightly stabilized (-0.2eV), this model leads to the expectation that the activation energy for the subsurface/bulk site
population should amount to at least 0.3-0.4 eV. Experimental values for H-uptake into
subsurface sites, e.g. for Pd(100), are much smaller (0.05 eV), and special sites such as defects
are suggested to be the responsible H-uptake channels [5]. In order to test this hypothesis wechose to measure the H-absorption kinetics at Pd(110), which should structurally resemble, e.g.,
step edge defects on Pd(100).
We further perform isotope-labeled (H2, D2) exposure sequences in order to investigate theisotopic exchange between surface-adsorbed and subsurface-absorbed hydrogen species, and thepossible involvement of chemisorbed H atoms in the H-absorption process in the presence ofgas phase H2. Such isotopic scrambling may be interpreted as an indication for an absorption
reaction coordinate between the gas phase H2 and H in bulk sites, on which interaction with thesurface-adsorbed H species may take place.
This question is intimately related to a long-standing problem in hydrogenation catalysis.
Subsurface-absorbed hydrogen has recently been demonstrated to be a necessary requirement
for the catalytic activity of Pd in olefin hydrogenation reactions, i.e., the H-addition to the C=C
bond of organic molecules [6]. It is unclear so far, however, how the subsurface hydrogenachieves this reactivity. Two scenarios are conceivable: a) absorbed H directly attacks theadsorbed olefin molecule, or b) the absorbed H modifies the surface electronic structure, thereby
rendering co-adsorbed surface H atoms more reactive to attach to the olefin. In order to help－49－

tackling this delicate question we investigate also the isotopic exchange between the surfaceadsorbed and subsurface-absorbed H species and their interactions with co-adsorbed olefin
molecules.
4.2

実験手法の概要

A Pd(110) specimen was precision-cut from a Laue-oriented single crystal rod, mechanicallypolished to mirror finish, and mounted onto a UHV sample manipulator that allowed for liquid
nitrogen cooling to 80 K and electron bombardment heating up to 1300 K. The experiments
were performed in an ultra-high vacuum apparatus, which has a base pressure < 1×10-8 Pa and isequipped for surface structure characterization by low energy electron diffraction (LEED) andfor thermal desorption spectroscopy (TDS). This system is connected to the 5 MV
van-de-Graaff tandem accelerator (MALT, University of Tokyo), which delivers the ~6.4 MeV
15 2+
N ion beam for surface-near hydrogen depth profiling and absolute surface-H coveragemeasurements by 1H(15N,αγ)12C nuclear reaction analysis (NRA) [7,8]. After careful in-situ
surface cleaning by Ar+ sputtering/annealing (1000 K) cycles followed by subsequent O2/H2
exposures that eliminate residual carbon impurities, the hydrogen ad-/absorption states andabsorption kinetics of the Pd(110) single crystal were characterized systematically by controlled
H2 (D2) exposures at 85-160 K and combined TDS, LEED, and NRA measurements.
4.3 実験結果a) TDS peak assignments
Fig. 4.1. shows the adsorption phases of the H2 exposed Pd(110)
surface at 100 K [9]. The (1×1) LEED pattern of the clean surface converts into a faint (2×1)
structure at minute H2 exposure (0.3 L), after which the (β1, β2) TDS peaks at 310 and 340 K are

0L

0.3 L

(1×1)

_

[110]

θ=0 ML

0.5 L

θ=1.5 ML

θ=? ML

α2 表面

表面
[001]

β1

streaky
(1×2)

(1×2)

(2×1)

θ=1.0 ML

50 L

α3
α1

β2

Fig. 4.1 Evolution of the LEED pattern and H2 TDS traces (1 K/s) of Pd(110) for increasing H2
exposures at 100 K. Small H2 dosages (< 1 L) suffice to saturate the surface-adsorbed H species.

－50－

observed. Further dosage to a total of 0.5 L turns the surface into a (1×2) pairing-row
reconstruction, and an additional peak (2) at 230 K appears in the TDS spectrum. The small H2
exposures and changes in the LEED pattern unambiguously identify the (1, 2, 2) TDS
features as surface-adsorbed H species. H coverages follow from the indicated structure models.
Subsequent changes take place on a slower time scale. Larger H2 exposures induce increasing
streaks in the (1×2) LEED pattern that can be rationalized by a loss of periodicity in the [001]surface direction through evolution of increasing numbers of missing rows in the reconstruction.
Concomitantly, two low temperature features (3, 1) slowly grow in the TDS spectrum. To
identify their origin, TDS spectra and NRA H depth profiles were recorded for larger H2dosages (1000-2000 L, 1 L = 1.33×10-4 Pa·s), after which these peaks become more prominent(Fig. 4.2).

a) NRA: H Depth Distributionb) TDS (1000 L H2)
α 2

α 1

PdHx

X=0.14 (水素化物)
α 3

1+3
3

X=0.007
(固溶体)

β 1

β 2

Fig. 4.2. a) NRA depth profiles of H in Pd(110)exposed to 2000 L H2 (in 100 s) at 130 K (blue) and→ Complete Peak Assignment:
145 K (red). The profile of the H-saturated surface(1
is shown for comparison (black).
L2,H2 1exposure), 2: Surface-adsorbed
Hb) TDS spectrum of 1000 L H2 dosed on Pd(110) at the1 : temperatures,
Near-surface
hydridesame
in comparison to the H2-saturated
3 : TDS trace
Bulk-dissolvedH
surface’s(dotted line).

Compared to the only surface-H saturated Pd(110) surface (black symbols), the NRA H
profiles after the large H2 exposure of 2000 L are asymmetrically broadened on the high energyside, which indicates H uptake into the subsurface region of the crystal. The extension of andthe H-concentration in this ‘subsurface region’ differs strongly depending on the H2 exposure
temperature (Te). At 130 K, a region of about 20 nm is affected. The topmost 3 nm include ~14
at.%, which largely exceeds the solubility of H in the low-concentration solid solution Pd
-phase at this temperature [10]. This condition has thus to be ascribed to a near-surface hydride
phase (-PdHx). It corresponds to the 1-feature in the TDS spectrum. At 145 K, the depth
profile indicates a much lower H-concentration below the surface, which, however, extends
rather far into the bulk of the crystal (data above 20 nm not shown). The corresponding TDS
trace indicates that this solid solution H desorbs from the Pd bulk in the 3 peak near 200 K. It
shows a tailing at high temperatures as indicative of H diffusion before desorption. The strong
Te dependence of the subsurface H distribution may be attributed to the fact that the H-diffusion
length in Pd increases rapidly with Te (1960 nm in 100 s at 130 K vs. 3570 nm at 145 K).
b) H-absorption kinetics
We determined the population rates of the (1, 3) absorbedhydrogen states in the H2 (and D2) exposure temperature (Te) range of 90-160 K from extensive－51－

TDS measurements. An exemplary uptake series for Te=90 K is shown in Fig. 4.3 a). Peakfitting analysis of the spectra is performed to determine the respective intensities of the α1 andα3 components in each spectrum. By graphing these populations against the H2 (D2) exposure,the absorption rate is obtained for a given Te. The Arrhenius plot of the (α1, α3) population
rates is finally shown in Fig. 4.3 b).

α2

a)
α1

α3

b)
β1

β2

Fig. 4.3. a) Uptake series of H2/Pd(110) TDS spectra for H2 exposures (t=100 s, H2 pressure varied) at
Te=90 K. b) Arrhenius plot of (α1, α3) population rates for the Te range of 90-160 K.

The analysis yields activation energies (Ea) of the α1 and α3 state populations of 0.13 eV and
0.06 eV for H2 and 0.17 eV for D2 (α3). These hydrogen absorption activation energies are
considerably smaller than expected by the most widely considered kinetic model, which
imagines that hydrogen absorption proceeds as a surface-to-subsurface transition ofchemisorbed H atoms and hence falsely predicts Ea≥0.3 eV. We present further evidence insection c) that this intuitive yet simple concept may need a revision.
The α1 feature was found to be much smaller in the D2 TDS spectra and hence not analyzed
further. We interpret this as a consequence of the enhanced D diffusion speed in the Pd bulk,
which is famously showing a rare ‘inverse’ isotope effect (D diffuses faster than H). The
absorbed D may thus more easily diffuse from the subsurface into the deeper crystal with lesschance to accumulate in the near-surface region to form a hydride phase. It is further worth
noting in this context that the activation energy for H2 absorption is smaller than that for D2, i.e.,the surface penetration step proceeds faster for H than for D. This is a ‘normal’ isotope effect,
remarkably opposite to the diffusion behavior of the hydrogen isotopes in the Pd bulk. A similar
trend has been reported for the absorption of H and D by Pd(100) [6].
c) Absorption mechanism
In order to investigate the possible involvement of
surface-adsorbed H atoms in the subsurface state population we performed sequential isotope
exposure experiments, where the surface was first saturated with dissociated D atoms by a small
D2 dosage (1.25 L) followed by a large (1000 L) H2 exposure to populate the subsurface states.
The result of this experiment (at 115 K) is shown in Fig. 4.4 a). All surface H species (β1, β2, α2)and the bulk-absorbed H (α3) show complete isotopic scrambling of H and D. Only the α1feature emerges purely as the post-dosed isotope (here H). This behavior is consistent against
reversal of the H2 and D2 dosage order (α1 selectively contains D when D2 is post-dosed).
This observation indicates that isotopic exchange with surface-hydrogen apparently takes
place for the bulk-absorbed H species (α3), but not for the near-surface hydride phase (α1). The
reason for this striking behavior has yet to be clarified, but it hints a the possibility that there
exists a separate reaction coordinate, on which H from the near-surface hydride can desorb
without interaction with the layer of surface-adsorbed hydrogen atoms. A similar bypassing of
surface hydrogen in the desorption from the subsurface region is known to occur also on
Pd(100) [5]. Microscopic reversibility lets us predict that the same reaction coordinate should

－52－

α1

b) 140x10

α2
α3

β1

QMS ion current [A]

a)

β2

-12

α2

120
100

0.8 L H2 at 130 K
quenched to 85 K
kept 80 sec at 130 K

β1 β

2

80
60
40
20
0
150

200

250

300

350

400

Temperature [K]Fig. 4.4. a) TDS spectrum of H2 and D2 after exposing Pd(110) at 115 to 1.25 L D2 followed by 1000 L
H2. b) Two TDS spectra of surface-H saturated Pd(110) (0.8 L H2 at 130 K): One tempered at the
exposure temperature for 80 s, the other recorded immediately after the H2 exposure.

also enable gaseous H2 to enter into the subsurface without the need to accommodate a single H
atom in the surface chemisorption wells (which are pre-occupied by the saturated surface
adsorption layer in our experiments). This raises interesting questions concerning the actual
point on the reaction coordinate where the H-H bond cleaves in the absorption step on the
H-saturated surface, which should in turn correspond to the location where recombination
occurs during H2 desorption from the near-surface hydride phase (α1).
We strongly suggest that this novel absorption reaction pathway involving molecular H2 anda H-saturated Pd surface be given theoretical attention, because in practice it will be the only
relevant absorption mechanism since at moderate temperatures the surface is always H-saturated
in the presence of any sizable H2 gas phase. To underline this point further, we present the resultof a control experiment in Fig. 4.4 b). Here, only surface-H is prepared on the Pd(110) surface at
130 K by a small H2 exposure (0.8 L). In UHV, the sample is then tempered at the same Te for aperiod of time (80 s) during which subsequent dosage of H2 or D2 would lead to substantial
population of the absorbed H states (similar to Fig. 4.4 a). Although all surface-H species (β1, β2,α2) are present at their maximum coverage in this experiment, not the faintest H uptake into thesubsurface states is discernible (every trace of α1 or α3 is missing in Fig. 4.4 b). This result
emphasizes that gaseous molecular H2 (or D2) is required in addition to surface-H in order toinduce hydrogen absorption. Surface-adsorbed H atoms obviously do not penetrate into thesubsurface by themselves; assistance (a portion of the H2 dissociation energy?) by gaseous H2 isnecessary to enable this process. Again, any surface-to-subsurface transition of H atoms doesapparently not start at the bottom of the chemisorption well. An alternative mechanism (a novel
reaction coordinate) is waiting to be discovered.
d) Catalytic reactivity
To address the question of the catalytic reactivity of subsurfaceabsorbed hydrogen in olefin hydrogenation, we prepared a Pd(110) surface with populated
subsurface states (α1 + α3) at 120 K and subsequently co-adsorbed cis-2-butene (C4H8). To test
this system for the hydrogenation reaction to butane (C4H10), the sample was warmed in astandard TDS run to monitor the reaction products by the mass spectrometer. The initial surface
condition before butene co-adsorption and the resulting reaction product are shown in Fig. 4.5.
Fig. 4.5 b) clearly indicates that a hydrogenation reaction of the butene to butane has
occurred during the TDS run. We confirmed that no such reaction occurs without subsurfaceabsorbed hydrogen, neither with surface-hydrogen only [11]. Apparently, desorption of the
butane reaction product peaks at 190 K and shows a nearly perfect overlap with the shape of the
α3-H2 desorption signal. On the contrary, almost no butane formation is seen in the temperature
range of the α1 desorption peak. This result indicates that although both α1 and α3 states－53－

originate below the Pd(110) surface, only one has the potential to hydrogenate thesurface-adsorbed olefin molecule.

a)b)Fig. 4.5. a) H2 TDS of 1000 L H2 at 115 K. b) TDS run of the identically H2-precharged sample after
co-adsorption of butene (C4H8), monitoring the hydrogenation reaction product butane (C4H10)
simultaneously with the H2 desorption. Only the H2 α3 component is shown for clarity.

One primitive explanation might be that the α1 near-surface hydride phase alreadydecomposes and desorbs at 170 K, a temperature possibly below the activation threshold of thehydrogenation reaction. A more interesting interpretation involves the possibly different
resurfacing/recombination reaction coordinates of the two absorbed H species that we discussed
in section c). We saw that α1-H desorbs independently and bypasses surface-adsorbed H atoms,
whereas α3-H emerging from the bulk does interact with surface species leading to isotopic
exchange. The evolution of the two H species from the subsurface on different spatial/energetic
reaction coordinates might also explain their different reactivity towards surface-adsorbed olefin
molecules. While our data are not conclusive on this interesting point at the present stage, theyat least do allow excluding one of the hypothetical origins of the catalytic hydrogenation
reactivity of ‘subsurface’ hydrogen. The maximum rate of butane formation is at 190 K, wherethe α1 near-surface hydride phase is almost completely decomposed. This makes the possibility
extremely unlikely that butene would react with surface hydrogen that has received its reactivity
through a modification of the surface electronic structure due to the presence of a concentrated
near-surface hydride layer.
4.4

まとめ
We achieved the first complete assignment of all H states in the complex H2 TDS spectrumof Pd(110), by identifying the (β1, β2, α2) features as surface-adsorbed H-species, (α1) as a
subsurface hydride phase (which forms only below ~140 K), and (α3) as bulk-dissolved H
atoms. The latter two species were distinguished for the first time. We determined the activationenergies (Ea) of the α1 and α3 state populations (0.13 eV, 0.06 eV) for H2 and 0.17 eV (α3) for
D2, which demonstrate with unprecedented clarity the inconsistency of conventional kinetic
models that envisage the hydrogen absorption process as a surface-to-subsurface transition ofchemisorbed H atoms (predicting Ea≥0.3 eV). Isotopic scrambling with surface-adsorbed H only
occurs during bulk-absorbed H (α3) desorption, but not during decomposition of thenear-surface hydride phase (α1), suggesting that the respective reaction coordinates to gas phase
H2 may differ for these two ‘subsurface’ species. Their reactivity towards catalytic butene(C4H8) hydrogenation shows a similar disparity: only the bulk-absorbed H (α3) affects H
addition to the olefinic C=C double bond, but not the H from the near-surface hydride phase(α1). The latter (α1) decomposes at 180 K before the butane (C4H10) reaction product is formed
at 190 K, which makes electronic structure changes induced by near-surface hydride an unlikely
scenario for the catalytic hydrogenation reactivity of Pd.

－54－

4.5

参考文献

[1] K. Christmann, Surf. Sci. Rep. 9 (1988) 1.
[2] K. H. Rieder, M. Baumberger, W. Stocker, Phys. Rev. Lett. 51 (1983) 1799.
[3] W. Dong, V. Ledentu, Ph. Sautet, A. Eichler, J. Hafner, Surf. Sci. 411 (1998) 123.
[4] M. G. Cattania, V. Penka, R. J. Behm, K. Christmann, G. Ertl, Surf. Sci. 126 (1983) 382.
[5] H. Okuyama, W. Siga, N. Takagi, M. Nishijima, T. Aruga, Surf. Sci. 401 (1998) 344.
[6] M. Wilde, K. Fukutani, W. Ludwig, B. Brandt, J.-H. Fischer, S. Schauermann, and H.-J.
Freund, Angew. Chem. Int. Ed. 47 (2008), 9289-9293.
[7] M. Wilde, M. Matsumoto, K. Fukutani, T. Aruga, Surf. Sci. 346 (2001) 482.
[8] M. Wilde, K. Fukutani, Phys. Rev. B 78 (2008) 115411.
[9] S. Ohno, M. Wilde, K. Fukutani, in preparation.
[10] Y. Fukai, The Metal Hydrogen System, 2nd Ed., Springer, Berlin, 2005.
[11] S. Katano et al., J. Phys. Chem. C 112 (2008) 17219.

5. 連携研究者・研究協力者
連携研究者*：
福谷克之(東京大学生産技術研究所教授)
笠井秀明（大阪大学大学院工学研究科教授）
研究協力者：
広瀬喜久治(大阪大学大学院工学研究科特任教授)
6. 本研究課題における平成２２年度の発表論文と講演
発表論文：
1)

2)
3)

4)

5)

H. K. Dipojono, A. A. B. Padama, N. Ozawa, H. Nakanishi, H. Kasai, "A First Principles
Study on Dissociation and Adsorption Processes of H2 on Pd3Ag(111) Surface", Japanese
Journal of Applied Physics , Vol. 49, pp. 115702 (2010).
N. Ozawa, M. Sakaue, H. Kasai, "Interactions of Hydrogen with Solid Surfaces", Journalof Vacuum Society of Japan, Vol. 53, No. 10, pp. 592-601 (2010).
A. A. B. Padama, M. C. S. Escaño, M. Y. David, H. Kasai, H. Kawai, "First principles
investigation for H2 dissociative adsorption on Ni and Cr-decorated Ni surfaces - Anapplication to alkaline polymer electrolyte fuel cell", e-Journal of Surface Science andNanotechnology, Vol. 8, pp. 325 (2010).
S. Aspera, M. David, H. Kasai, "First principles study of the Adsorption of Water
Molecule on Tri-s-triazine based Graphitic Carbon Nitride (g-C3N4)", Japanese Journalof Applied Physics, Vol. 49, pp. 115703 (2010).
D. N. Son, B. T. Cong, H. Kasai, "Hydronium Adsorption on OOH Precovered Pt(111)
Surface: Effects of Electrode Potentiale", Journal of Nanoscience and Nanotechnology,
Vol. 11, pp. 2983 (2011).

－55－

6)

7)

8)

9)

H. S. Kuncoro, M. Sakaue, H. Nakanishi, H. K. Dipojono, H. Kasai, "First-principles
Investigation on Ionization Strength, Volume Expansion, and Water Rotational Rigidity
of Small Water Cluster Systems Formed around Sodium(I), Calcium(II), and Iron(II)Ions", Journal of the Physical Society of Japan, Vol. 80, pp.024601 (2011).
Hidekazu Goto, Kikuji Hirose, “Electron-Electron Correlations in Square-Well
Quantum Dots: Direct Energy Minimization Approach”, J. Nanosci. Nanotechnol., 11(2011) pp.2997-3004.
H. Nagara, K. Mukose, T. Ishikawa, M. Geshi, and N. Suzuki, “Origin of the simple
modulated structures and the pressure induced superconductivity”, Journal of Physics:
Conference Series, 215, 012107(2010).
M. Geshi, “Structural properties of ferromagnetic CaN in a CsCl-type and rock salt
structure: A first-principles study”, Physica B, 405, 517-523(2010) .

講演：
1. 中西寛,“磁性制御による表面反応デザインと量子ダイナミクス”, 山田科学振興財
団 新規プログラム（若手研究者の分野交錯研究推進プログラム）研究会「ナノスケールの原子構造および電子状態とその変化の理論的設計」, 2010 年 8 月 31 日-9
月 1 日, 国際高等研究所，日本【招待】2. 中西寛，笠井秀明，“プロトン,ミューオンの量子ダイナミクス”，中間子科学研究会「超低速ミューオン：その無限の可能性を探る」，2010 年 9 月 25 日，山梨大学，
日本
3. Allan Abraham Padama, Melanie Yadao David, 笠井秀明, Kawai Hiroyuki, "Density
functional theory investigation on H2 dissociative adsorption on Ni(111) and Cr-decorated
Ni(111) surfaces", 第 51 回真空に関する連合講演会, 2010 年 11 月 4-6 日, 大阪大学,
日本
4. 中西寛,“プロトン・ミューオンのための次世代第一原理計算手法”, 東北大通研共同
プロジェクト研究「次世代第一原理計算手法の開発•応用」, 2010 年 11 月 10-11 日,
東北大学電気通信研究所，日本 【招待】5. 笠井秀明，Wilson Agerico Diño, "燃料電池の課題と展望", 2011 年春季 第 58 回応用物理学関係連合講演会, シンポジウム「エネルギー・環境研究会企画「新エネルギ
ーとそのネットワーク化 スマートグリッドにおける応用物理の果たすべき役割と新しい可能性 」, 神奈川工科大学, 2011 年 3 月 24-27 日【招待】6. 中西寛、笠井秀明，“水素・ミューオンの固体表面近傍での量子ダイナミクスと電子状態”, 日本物理学会第 66 回年次大会, 2011 年 3 月 25-28 日, 新潟大学，日本
7. 笠井秀明，“表面・界面・固体中のプロトン伝導とその量子ダイナミクス計算”、
日本物理学会第 66 回年次大会 領域４シンポジウム「ナノスケール量子輸送の計
算科学的研究の現状・展望と次世代スパコンへの期待」,2011 年 3 月 25-28 日, 新潟大学，日本【招待】8. 國貞雄治，中西寛，Wilson Agerico Diño，笠井秀明, "銀表面上のステップ構造近傍での水素分子の吸着状態",日本物理学会第 66 回年次大会, 2011 年 3 月 25-28 日, 新潟大学，日本9. A. Sasaki, K. Hirose, H. Goto, “Electron-Electron Interactions in Few-Electron Quantum
Dots: Direct Energy Minimization Approach”, Third International Symposium on
Atomically Controlled Fabrication Technology (Suita, November 24-26, 2010)
10. 佐々木晃、広瀬喜久治、後藤英和、”非直交基底による電子相関エネルギーの高精度計算手法”,日本物理学会第 66 回年次大会, 2011 年 3 月 25-28 日, 新潟大学，日本－56－

11. 下司雅章，長柄一誠，向瀬紀一郎，“ハロゲン incommensurate 相のフォノン計算”
日本物理学会第 66 回年次大会，2011 年 3 月 25 日-28 日，新潟大学，日本
12. Markus Wilde, “Surface/subsurface transitions of hydrogen on metallic single crystals andnanoparticles”, International Conference on Diffusion in Solids and Liquids (DSL-2010)(Paris, France, July 5-7, 2010) 【招待】13. 大野哲, Markus Wilde, 福谷克之, “核反応法による Pd(110)サブサーフェス水素深さ分布の曝露条件依存性に関する研究”, 日本物理学会 2010 年秋季大会, 2010 年 9
月 23-26 日, 大阪府立大学, 日本(24pWS-9)
14. Markus Wilde, “Elucidating Hydrogen Absorption at Metal Surfaces and Nanocrystals by
1
H(15N,)12C Nuclear Reaction Analysis”, 第１１回「イオンビームによる表面・界面解析」特別研究会, 2010 年 12 月 3-4 日, 名城大学, 日本【招待】15. 大野哲、Markus Wilde, 福谷克之,“共鳴核反応法を用いた Pd(110)への水素吸蔵機構の研究”, 日本物理学会第 66 回年次大会, 2011 年 3 月 25-28 日, 新潟大学、日本（25aTG-1）－57－

多自由度・大規模系における反応と構造空間探索
Reaction Analysis and Conformational Sampling forMultidimensional Large Scale Systems
倭

剛久 1、光武亜代理 2
T. Yamato, A. Mitsutake

名古屋大学 1、慶応義塾大学 2
Nagoya University, Keio University
視物質ロドプシンの吸収波長制御メカニズム：
ロドプシンは超高感度の光センサータンパク質である。その光吸収極大波長は棲息域の光学的環境に適応している。海中の魚類は、深度に応じてロドプシンのアミノ酸配列を変化させ、光吸収極大波長を調節している。我々はアナゴのロドプシンとその祖先型ロドプシンの立体構造をホモロジーモデリング法で予測し、長時間分子シミュレーションで統計アンサンブル
を生成、さらに高精度の分子軌道法で吸収スペクトルを解析した。そして、吸収波長制御の分子機構をあきらかにした。

１． はじめに
地球上の生物にとって太陽の光は生命活動に欠くことのできない存在である。光を検知したり、光エネルギーを化学エネルギーに変換する役割は光受容蛋白質や光合成蛋白質などが担っている。これらの蛋白質の機能メカニズムを解き明かすには、分子内部でおこっている物理的なプロセスを解析し、蛋白質を構成する個々のアミノ酸残基の役割を調べることが肝要である。そのために、分子シミュレーションと電子状態計算を活用したアプローチが有効となる。
光情報伝達に関して、生物は進化の過程で視物質ロドプシンという優れたナノマシンを獲得した。例えばほ乳類のロドプシンは光子を１つ吸収すると、約７０％の確率で光反応を起こすことが知られている。また、様々な生物が持つ種々のロドプシン分子は360nm から 570nm にわたる幅広い範囲の吸収極大波長(λmax)を持つ。これはロドプシン中に含まれるレチナールを取り囲むアミノ酸残基の顔ぶれが異なることで、レチナール
近傍の構造が変化したり、レチナールと周囲との相互作用が変化したりすることに起因する。
ロドプシンの電子状態計算はコンピューティクスにとってチャレンジングな問題を提供する。熱揺らぎをする蛋白質環境の中でレチナールの構造や電子状態を考慮し、さらにその励起状態を精密に評価することで、ロドプシンの吸収波長制御の問題にアプローチできる。
2008 年、Yokoyama らはアナゴの祖先型ロドプシン(p501、図１)に２つのアミノ酸（N195, A292）を同時に置換することでλmax が 13nm ブルーシフトすることを発見した[1-2]。深度の深い海では赤い光は吸収されて見えないため、棲息する深度に応じてロド
プシンのλmax が異なっていると考えられる。興味深いことにこれら２つのアミノ酸の１つはレチナール近傍に、もう１つはレチナールから遠く離れた場所に位置しており、いずれの一方をアミノ酸置換してもλmax は変化しない。すなわち、２つのアミノ酸が共同的に働いてλmax を制御していると考えられる。

－58－

２． 手法 [3]アナゴロドプシンにおける吸収波長制御機構をしらべるため、まず、立体構造未知のアナゴ祖先型ロドプシン(p501)、アミノ酸残基の単置換体(SM: A292S)、および２重置換体(DM: N195A/A292S)の構造をホモロジーモデル法を用いて作成した [4]。
ロドプシンの構造空間上に存在する多くのエネルギー極小点の分布により、レチナー
ルの可視吸収スペクトルにはブロードな不均一広がりが生じていると考えられる。本研
究では、この不均一広がりを考慮するため、上記３種類のロドプシンの構造サンプリン
グを各々, 実効的に 48 ns の分子動力学シミュレーションを用いて実行した。各々100個の瞬間構造を抽出し、古典力場関数を用いたエネルギー最適化と QM/MM 法によるエネルギー最適化をつづけて施した。つぎに、これらのエネルギー極小構造のすべてについて CASPT2 法を用いてエネルギー計算を実行し、励起エネルギーの分布を解析した。

３． 結果
表１にλmax の計算値を示す。励起エネルギーの絶対値は過大評価になっているが、２重置換体のブルーシフトを正しく再現することができた。CASPT2 法レベルの計算をこれほど大規模に行った計算は著者の知りうる限り初めての例である。
次に、励起エネルギーの分布を制御するロドプシンの構造因子の抽出である（図２）。
そのため、われわれは、主成分分析法を改変して用いた。100 個の構造の各々の励起エネルギーの値を拡張次元として考慮し、「立体構造の自由度＋１」次元の拡張空間で主成分解析を実行した。レチナールの共役鎖の一重結合／二重結合の結合長の交代の度合
いを bond length alternation と呼び, シッフ塩基レチナール（正電荷）と近傍の脱プロト
ン化アミノ酸（負電荷）の間の距離変化を counterion displacement と呼ぶ。図２の第一
主成分は励起エネルギー変化と bond length alternation, conterion displacement の成分を多く含んでおり、これら２つの構造要因がλmax の制御に密接に関係していると考えられる。また、第二主成分はレチナールの C6-C7 ボンドの捻れに対応していた。

5

0

-5

-10

-15

図１．ロドプシンの構造

-10

-5

0

5

10

15

図２. 主成分分析
横軸：第一主成分、縦軸：第二主成分

－59－

p501
SM
DM

Calculation
2.91 (426) / 0.12
2.90 (427) / 0.15
2.96 (418) / 0.11

Experiment2.47 (501)/ (0.53)*
2.47 (502) / —2.54 (488) / —表１
励起エネルギーのピーク値 [eV] ( λmax ([nm]) / 半値幅 [eV]*) p501 の半値幅の実験値はウシロドプシンの実測値による参考値 [5].

４． まとめホモロジーモデリング、分子動力学シミュレーションによる構造サンプリング、および蛋白質の電子状態計算により、魚類ロドプシンの進化過程でおこる吸収波長変化の機構を調べた。
アナゴの祖先型ロドプシン、およびそのアミノ酸一重置換体、二重置換体について励起エネルギーの 100 点計算を各々実行した結果、二重置換体におけるブルーシフトを再現することができた。
拡張した主成分分析をもちいて、蛋白質構造変化と吸収極大波長の相関を解析したと
ころ、bond length alternation と counterion displacement との強い相関をみた。

５． 参考文献
1. S. Yokoyama, T. Tada, H. Zhang and L. Britt, Proc Natl Acad Sci U S A 105 (36), 13480-13485(2008).
2. C. Day, Physics Today 61 (10), 20-22 (2008).
3. H. C. Watanabe, Y. Mori, T. Tada, S. Yokoyama and T. Yamato, Biophysics (Oxf) 6, 67-68(2010).
4. A. Sali and T. L. Blundell, J Mol Biol 234 (3), 779-815 (1993).
5. R. R. Birge, Biochimica Et Biophysica Acta 1016 (3), 293-327 (1990).

６． 連携研究者・研究協力者
連携研究者：
岡本祐幸（名古屋大学理学研究科教授）、足立伸一（高エネルギー加速器研究機構准教授）
研究協力者：
杉田有治（理化学研究所准主任研究員）、渡辺宙志（名古屋大学理学研究科学振特別研究員）

７． 本研究課題における平成２２年度の発表論文と招待講演
1)

発表論文：
H. C. Watanabe, Y. Mori, T. Tada, S. Yokoyama, and T. Yamato
Molecular mechanism of long-range synergetic color tuning between multiple amino acid residues inconger rhodopsin Biophysics 6, 67-78 (2010).

－60－

2)

3)

4)

A. Mitsutake, Yoshiharu Mori, and Yuko Okamoto, Multi-dimensional multicanonicalalgorithm, simulated tempering, replica-exchange method, and all that. Physics Procedia 4,89-105 (2010).
倭 剛久
はたらくタンパク質：エネルギー、電子、プロトン、リガンド移動の計算生体高分子科学 高
分子論文集（特集「計算高分子科学」） 67, 179-186 (2010)
分子シミュレーション研究会会誌 多変数拡張アンサンブル法（特集） 光武亜代理、分子シミュレーション研究会会誌「アンサンブル」12(2),29—35 (2010)

招待講演：
1) T. Yamato
Intramolecular communication chart of proteins: Computational analysis of energy transfer
pathways, Thermal transport at the nanoscale, Telluride Science Research Center Workshop,
Telluride, CO, USA, Jun. 21-25, 2010.
2) T. Yamato
Proteins at work: Computational biophysics of energy, electron, proton transfer and ligand
migration, COE Start-up International Workshop, Chiba University, Chiba, Japan, Nov. 11-12,
2010.
3) A. Mitsutake
The Fourth Shanghai International Conference on Biophysics and Molecular Biology(2010SICBM), Shanghai, Chaina, August 2010－61－

第一原理有効模型と相関科学のフロンティア
Frontiers of Ab Initio Effective Models and Correlation Science
今田正俊 1、三宅隆 2、中村和磨 1、佐久間怜 3
有田亮太郎 1、石橋章司 2、小口多美夫 4、藤森淳 1
M. Imada, T. Miyake, K. Nakamura, R. Arita, S. Ishibashi, T. Oguchi, A. Fujimori
東京大学 1、産業技術総合研究所 2、千葉大学 3、大阪大学 4、
The University of Tokyo, AIST, Chiba University, Osaka University
第一原理的に強相関電子系の電子状態を明らかにする計算手法を開発し、広範に応用していくとともに、超伝導体やスピン軌道相互作用の強い系など特徴ある物性を示す物質群の解析を進めている。

１． はじめに
我々の日常目にする物質の物性の多くが多数の電子の挙動によって決まる。この多電子の中には、原子核付近のコア電子から、フェルミエネルギー付近の電子に至るまで、数十から百 eV におよぶエネルギー差があり、これをまたぐマルチスケールの相互作用
が多電子の階層構造を形成する。この階層構造の果てに 3 桁以上小さな常温以下のエネ
ルギースケールの物性が決まる。階層性を利用して自由度を縮減するダウンフォールディ
ング法により、着目するエネルギースケールのための少数自由度のみの有効模型が非経験的に導出され、有効模型を精緻な低エネルギーソルバーで解くことにより、従来までの密度汎関数法では困難な強相関電子系を高精度、低負荷で取り扱える新手法を確立した。我々はこの密度汎関数法と強相関模型解法を融合した手法のさらなる展開を図り、現実の物質群への広範な実証研究を展開し始めている。
２．手法の開発
２－１．次元縮約法（中村、今田）
従来までの「制限乱雑位相近似 (cRPA) 法」は常に三次元模型のパラメータを導出するのに対し、擬一次元・二次元系の物性解析では、純粋一次元・二次元模型が解析対象となる。導出模型の次元と解析模型の次元は異なっており、次元縮約された有効模型の導出理論およびプログラム開発が必須となる。cRPA 法で導入される「制限」をバンド
自由度だけでなく、空間自由度にも拡張することにより、上記の問題を解決した。擬二次元系鉄系超伝導体 LaFeAsO の純粋二次元模型の導出を行った。FeAs 層間の遮蔽効果を取り込むことにより、(i) 三次元相互作用特有の長距離型クーロン相互作用が短距離型に変わること、(ii) オンサイト相互作用に対してさえ、10-15%の値減少効果があることを見出した。
２－２．応力テンソル計算法（小口）
構造最適化や弾性定数の第一原理計算で重要となる応力テンソルの全電子フルポテ
ンシャル線形化補強平面波（FLAPW）法での定式化を提案した。応力テンソルの計算
は擬ポテンシャル法によるものがほとんどで、全電子 FLAPW 法による報告は１報あるのみでその応用も Al と Si の圧力計算だけである。ここでは、Soler-Williams の LAPW
基底による応力テンソルの式は、従来の LAPW 基底のものよりも間便な表式を与える
ことを示した。種々の遷移金属系や化合物系の応力テンソル計算を実行し、圧力（応力－62－テンソルのトレース）の評価からいずれの系でも 0.1GPa 以内の高い精度を与えることが検証された。また、非対角要素から弾性定数を見積り、エネルギーから求めたものとコンシステントであることが示された。
３．鉄系超伝導体の電子状態
３－１．磁気秩序と磁気モーメント（中村、今田）
第一原理的に導出された鉄系超伝導体の有効模型に対して、多変数変分モンテカルロ
法を適用し、磁性秩序構造、秩序磁気モーメントの両方において定量的に LaFeAsO の実験結果を再現することを示し、またそのほかの化合物の磁気モーメントについてもその多様な変化を説明できることを示した。実験的には小さな磁気秩序モーメントであるのに対して、従来の LSDA 計算で大きなモーメントを予測してしまっていた謎を解き、
LaFeAsO が磁気量子臨界点近傍にあることを第一原理的に立証した。また多くの異なる
鉄系超伝導体が相互作用の強さを 1 つのパラメタとする普遍的な模型で磁性を良く説明できることを示し、電子相関がこの物質群の多様な物性を決める鍵を握っていることを示した。
３－２．Ca4Al2O6Fe2As2 の電子状態 (三宅、石橋)
鉄ヒ素超伝導体に対する精力的な研究が続いている。研究の初期より構造と超伝導の相関が指摘されてきた。様々な鉄ヒ素物質の超伝導転移温度をヒ素-鉄-ヒ素の結合角αに対してプロットするとユニバーサルな線が描け、正四面体に対応する 109 度で極大点をもつことが知られている。昨年、鉄ヒ素超伝導体の中で最小のαをもつ
Ca4Al2O6Fe2As2 が合成された。この物質の電子構造を密度汎関数法で調べたところ、
LaFeAsO よりΓ点の周りでフェルミ面が１枚少ないことが分かった。仮想的にαを変えた一連の計算から、結合角が電子構造に大きく影響することがわかった。Γ点でのエネ
ルギー準位を見ると、αの増加とともにフェルミ面の下-0.5eV 付近のバンドが上昇し、フェルミ面を横切る。やがて約+0.2eV に位置する２重縮退点をすり抜けると共にバンドの組み替えがおこり、バンド構造のトポロジーが変わる。最局在ワニエ関数を用いた解析から、このふるまいは xy 軌道のオンサイトエネルギーと最隣接鉄サイトへのホッ
ピングの変化に起因することがわかった。超伝導転移温度との関連に興味が持たれる。
４．スピン軌道相互作用の強い物質群の電子状態
４－１．Ir 酸化物系の LDA+DMFT 計算（有田、今田）
3d 電子系、4d 電子系、5d 電子系を比べた場合、周期律表で下にいけばいくほど軌道の広がりが大きくなって弱相関系になっていくことが期待される。しかしながら、
Sr2RhO4 と Sr2IrO4 においてはこの単純な予想と逆の振る舞いが見られることが知られている。すなわち、前者が金属であるのに対して、後者はモット絶縁体である。この事実は、後者において、強いスピン軌道相互作用の存在によって、t2g 軌道の組み換えがおこり、J=1/2 のバンドが half-filled になった 2 次元電子系が実現した、という具合に理解できる。
2 次元単一軌道系のモット絶縁体としては、銅酸化物がまず念頭に浮かぶが、Sr2IrO4が銅酸化物同様高い転移温度をもつ超伝導体になり得るかという可能性を考察することは非常に興味深い。両者の結晶構造を比べたときに、Sr2IrO4 に見られる特徴として考えられるのが、O-Ir-O の角度が 180 度から大きくずれていることである。この角度を180 度に近づけるためには、Sr をイオン半径の大きなものにかえることが考えられる。
実際、最近、青山学院大学・秋光グループにおいてこの試みがなされた。そこで、相対論的第一原理計算によって構造最適化を行い、Ba イオンに変えることによってどの程度角度が広がるかについて検討を行った。さらに電子状態計算の結果からワニエ軌道を作って、低エネルギー有効模型を構築し、動的平均場理論を用いて解析した。

－63－

O-Ir-O の角度は Ba の導入により広がる傾向があるが、完全に 180 度になることはな
く、Sr2IrO4 と比べて劇的に銅酸化物の状況に近づくことはできないという結果が得られた。今後外部圧力などの効果を検討し、銅酸化物の 5d analog を作る可能性を検討する予定である。
４－２．QMAS 上でのコードの開発 （石橋、三宅）
トポロジカル絶縁体や電気磁気効果などスピン軌道相互作用が引き起こす物理現象が注目を集めている。これらの効果はエネルギースケールが小さく、表面や界面など複雑な物質系を舞台にすることが多い。私たちは、信頼できる高精度な計算プログラムとして平面波基底の PAW コードである QMAS（http://www.qmas.jp/）上で相対論効果を考慮した第一原理計算プログラムの開発を行っている。初年度はプログラムの整備と簡単な系への応用を行った。fcc 構造の金の電子構造が過去の計算例とよく一致することを確認した。ついで金(111)面の薄膜の計算を行った。膜厚 N が小さいとき、両側の表面状態が混成して結合状態と反結合状態を形成し、Γ点で分裂が見られる。N を増やすとこの分裂は急激に小さくなり、また N>14 ではΓ点からずれたところにバンド分散曲線の極小点の波数 k0 が移動する。さらに N を増やすと k0 は増加するが、そこでの分裂幅は N に対して鈍感であることが分かった。これらの結果に対する自由電子模型とタイ
トバインディング模型による解析も行った。
４－３．スピン軌道相互作用を取り入れた GW 近似（佐久間）近年スピントロニクスの分野でトポロジカル絶縁体やイリジウム酸化物 Sr2IrO4 などスピン軌道相互作用が電子状態の理解に本質的に重要である系が注目されている。これらの系の一電子スペクトルを精密に計算する第一原理電子状態計算手法として、GW 近似をスピンに依存する系に拡張したスピン依存 GW 近似が近年 Aryasetiawan らによって提案されている 。 我 々 は 今 年 度 、 ユ ー リ ッ ヒ ( ド イ ツ ) の 理 論 グ ル ー プ と の 共 同 研 究 に よ り 彼 らの
FPLAPW-GW コードにスピン依存 GW 近似の実装を行った。典型半導体やトポロジカル絶縁体への適用を行い、得られた結果は一般的手法である局所密度近似の結果を改善する傾向があることがわかった。また、スピン軌道相互作用のある系において再局在ワニエ関数と制限 RPA 法によるハバードパラメタ U を評価する手法の実装を行い、第一原理計算と低エネルギー模型を組み合わせた電子状態計算のためのコードを整備した。
５．光電子分光５－１．鉄系超伝導体の光電子分光（藤森）鉄系超伝導体のバンド構造とフェルミ面を角度分解光電子分光法により調べ、第一原理計算と比較した。その結果、次に述べるような系統的な変化を見出した。
（1）PrFeAsO1-y、 BaFe2(As1-xPx)2 の光電子分光測定と第一原理計算から、鉄原子層に対するニクトゲン原子層の高さが低下すると、dz2 バンドが上昇してフェルミ面への寄与が増し、フェルミ面の 3 次元性が増すことが示された。
（2）Ba1-xKxFe2As2、Ba(Fe1-xCox)2As2、Ba(Fe1-xNix)2As2 の光電子分光測定と第一原理計算から、ホールドーピングによりフェルミ面がより 2 次元的に、電子ドープによりフェルミ面の 3 次元性が増すことが示された。
５－２．第一原理計算との比較（小口）角度分解光電子分光で観測されたバンド構造に関して、局所密度近似(LDA)の範囲での第一原理計算に基づき詳細な解析を実行した。YNi2B2C については、ｋ依存超伝導ギャップのバンド毎の異方性を明らかにした。

－64－

６．そのほか （小口）強磁性 Fe に関しては、電子間相互作用及び電子格子相互作用によるカップリングパラ
メータを評価し、有効質量の増大は主として前者に起因することを定量的に示した。
７．連携研究者
辛 埴（東京大学物性研究所）

８．研究協力者
三澤貴宏（東京大学工学系研究科助教）

９．本研究課題における平成２２年度の発表論文と招待講演
発表論文：
1) Takahiro Misawa, Kazuma Nakamura, Masatoshi Imada, "Magnetic Properties of Ab initio
Model for Iron-Based Superconductors LaFeAsO", J. Phys. Soc. Jpn. 80 (2011) 023704(1-4).
2) Masatoshi Imada, Takashi Miyake, "Electronic Structure Calculation by First Principles forStrongly Correlated Electron Systems", J. Phys. Soc. Jpn. 79 (2010) 112001(1-42).
Takashi Miyake, Taichi Kosugi, Shoji Ishibashi and Kiyoyuki Terakura, “Electronic Structureof Novel Superconductor Ca4Al2O6Fe2As2” J. Phys. Soc. Japan 79, 123713 (2010). (4 pages)3) S. Sakai, Y. Motome, M. Imada, "Doped high-Tc cuprate superconductors elucidated in thelight of zeros and poles of electronic Green's function", Phys. Rev. B 82 (2010) 134505(1-16).
4) Kazuma Nakamura, Yoshihide Yoshimoto, Yoshiro Nohara, and Masatoshi Imada, “Ab initio
Low-Dimensional Physics Opened Up by Dimensional Downfolding: Application to
LaFeAsO”, J. Phys. Soc. Jpn. 79, 123708 (2010).
5) Despina Louca, K. Horigane, A. Llobet, R. Arita, S. Ji, N. Katayama, S. Konbu, K. Nakamura,
T.-Y. Koo, P. Tong, and K. Yamada, “Local atomic structure of superconducting FeSe1−xTex”,
Phys. Rev. B 81, 134524 (2010).
6) Kazuma Nakamura, Ryotaro Arita, and Hiroaki Ikeda, “First-principles calculation oftransition-metal impurities in LaFeAsO”, Phys. Rev. B (to be published).
7) Yu-ichiro Matsushita, Kazuma Nakamura, and Atsushi Oshiyama, “Comparative study ofhybrid functionals applied to structural and electronic properties of semiconductors andinsulators”, submitted to Phys. Rev. B.
8) T. Yoshida, I. Nishi, S. Ideta, A. Fujimori, M. Kubota, K. Ono, S. Kasahara, T. Shibauchi, T.
Terashima, Y. Matsuda, H. Ikeda, R. Arita: Two- and three-dimensional Fermi surfaces andtheir nesting properties in superconducting BaFe2(As1-xPx)2, Phys. Rev. Lett., in press－65－9) T. Baba, T. Yokoya, S. Tsuda, T. Watanabe, M. Nohara, H. Takagi, T. Oguchi, and S. Shin,“Angle-resolved photoemission observation of the superconducting-gap minimum and its
relation to the nesting vector in the phonon-mediated superconductor YNi2B2C”, Phys. Rev. B
81, 180509/1-4(R) (2010).
10) X. Y. Cui, K. Shimada, Y. Sakisaka, H. Kato, M. Hoesch, T. Oguchi, Y. Aiura, H. Namatame,and M. Taniguchi, “Evaluation of the coupling parameters of many-body interactions in Fe(110)”, Phys. Rev. B 82, 195132/1-9 (2010).
11) N. Nagasako and T. Oguchi, “Stress formulation in the all-electron full-potential linearlized
augmented plane wave method”, J. Phys. Soc. Jpn. 80, 024701/1-13 (2011).

招待講演：
1) M. Imada: “Renewed Understanding on Doped Mott Insulators”. American Physical Society
March Meeting (Dallas, USA, March 21-25, 2011)

依頼講演:
1) 石橋章司、“QMAS紹介(概要と使い方)”および“QMASのより進んだ使い方”、HPC
産業応用スクール：ナノテクコース、東京、2011年3月9-10日－66－

第一原理からの多体理論
Many-body theory from first principles
髙田康民 1、白井光雲 2、前園涼 3、前橋英明 1、吉澤香奈子 4
Y. Takada, K. Shirai, R. Maezono, H. Maebashi, K. Yoshizawa
東京大学 1、大阪大学 2、北陸先端科学技術大学院大学 3、上智大学 4
University of Tokyo1, Osaka University2, JAIST3, Sophia University4
第一原理系における近藤効果の物理の検証という観点から電子ガス中の１原子系に注目し、拡散量子モンテカルロ法や LDA を用いた詳細な研究に着手した。また、
GWΓ法を高度化して強相関電子系での自己エネルギーも計算可能なように改良を進めている。更に、ホウ素・炭素系における高温超伝導出現の可能性を調べた。

１． はじめにA03 髙田班の研究課題は、「第一原理からの多体問題」という大テーマの中からいく
つかの課題を抽出して、「第一原理系励起状態の多体論と高転移温度超伝導物質デザイン」とした。初年度（平成２２年度）の研究期間は約半年と短いことも考慮して、根本になる大テーマの研究に取り組む際の視点や姿勢を明確にすることから始めた。
一般に、物質科学のためにスーパーコンピュータを駆使して実行される第一原理計算
は第一原理のハミルトニアン HFP を出発点にして既存手法（あるいは、それを漸次改良した手法）を用いたものであるのに対し、我々の視点は(1)計算手法の根本的な開発、
たとえば、密度汎関数理論（DFT）における交換相関エネルギー汎関数 Exc[n(r)]や時間依存密度汎関数理論（TDDFT）における交換相関核 fxc(rt,r’t’)の探索・改良、に直接的に結びつく研究、(2)モデルハミルトニアンを通して確立されてきた物理概念、たとえば、不純物アンダーソン模型における近藤問題の物理やハバード模型におけるモット転
移、を HFP から出発して直接的に検証し、新概念の発見を目指す研究、の２つである。
ところで、「①電子ガスに埋め込まれた１原子系」はこれら両視点のいずれから見ても大変に興味深い研究テーマである。すなわち、(1)の立場からは、一様密度の電子ガス系の研究に立脚して導かれた局所密度近似（LDA）を超えて、原子核のポテンシャルに由来する非一様性を正しく取り込む Exc[n(r)]を研究する際に欠かせないテーマであると同時に、(2)の立場からは、伝導電子間には相互作用のない不純物アンダーソン模型に立脚して導かれている物理を伝導電子間にも長距離クーロン相互作用のある HFP で再検討するものである。このテーマは前園と吉澤が担当している。また、「②GWΓ法：自己無撞着な自己エネルギーの計算手法の確立と遂行」というテーマも、(1)の視点から
は fxc(rt,r’t’)の開発・改良に直接的に関連しており、また、(2)の視点からは朝永・ラッティンジャー流体やモット転移を HFP から出発しての直接的に記述するという目標に関連していて、大変重要なものである。これは前橋と髙田が担当している。
この他、「③高温超伝導体の探索」に関連して超伝導転移温度 Tc の問題がある。今年
度は超伝導密度汎関数理論（SCDFT）における対相互作用を検討し、弱および強結合の両極限におけるその汎関数形のあるべき姿を髙田が議論した。これとは別に、実験との協力関係を全面に押し出した超伝導の研究を白井が行っている。以下、①から③の各テーマについて、それぞれ、もう少し詳しく報告し、得られた成果をとりまとめる。
２．電子ガスに埋め込まれた１原子系
電子ガス中の１原子系は原子核の電荷 Z と電子密度を表す rs の２つのパラメータで記述される。今、Z≪1 では一様密度の電子ガス系から出発した線形応答理論で厳密解が－67－

Number of electrons得られる。そして、遍歴電子による遮蔽効果で原子核の周りに電子は局在しない。反対に、Z≫1 では、電子が原子核に束縛され、局在化してくる。そして、rs も大きくなると、原子が希薄電子ガス中に存在している描像になる。しかし、中間領域の Z や rs では、近藤効果で問題になるスピン偏極の可能性も含めて、その詳細は不明である。
前園はこの系の基底状態の電子密度 n(r)を Z と rs の関数として拡散量子モンテカルロ（DMC）法で正確に得ることに着手した。そして、今年度は少なくとも Z≦2 で rs≦23
では LDA がほぼ正確な n(r)を与えることを確認した。
（なお、rsÆ∞のヘリウム原子では LDA が正しい結果を与えないことは分かっている。）
rs = 15
そこで、Z≦2 に限定するものの、吉澤は線形応答理論の適用
2
範囲を越えた領域での電子状態
を LDA で詳しく研究した。
特に、
1.5
コーン・シャム軌道の「位相のずれδ」から束縛状態の有無を調1べた。その結果、Z=2 の場合、rsが 0.5 を越えるとδが正（引力的）0.5から負（斥力的）に突然に変化
Z=2
することを見いだしたが、これZ=10
は束縛状態が１つできたことを
0
5
10
15
20
25
30
35
40
意味する。この束縛状態ができR [a.u.]
ると、遍歴電子の散乱状態はこ
図１：rs=15 に対する Z=1 と Z=2 の原子核をれに直交するための斥力を感じ中心とした半径 R の球内にある総過剰電子数。
てδが負になるのである。さらに
rs を大きくしていくと、斥力が序々に弱められ連続的に引力に変わり、ついに rs>7 でもう一度引力から斥力に急に変わる。これは２つ目の束縛状態出現を意味する。Z が 1 の場合でも同様で、rs が 1.8 以
上で１つ目の、そして、rs が 7 以上で２つ目の束縛状態ができる。このように、δの解析により、Z が１と２の両方とも系が低密度になると２つの束縛状態ができることが分かった。ただ、いずれの状況でもスピン偏極は見られない。
さらに、束縛状態２つが出現している状況（rs=15）での電子密度を調べ、原子核を中心とした半径 R の球内に存在する総過剰電子数を R の関数としてプロットしたのが図１である。もちろん、この総過剰電荷数は RÆ∞では電荷の完全遮蔽条件から原子核の電荷 Z に収束するが、その極限値には振動しながら近づく。この際、Z=1 の場合の大きな振幅が興味深い。これは共鳴的な束縛状態とはいえ、H-が出来る状況
では、遍歴電子からみると、原子核に近いところでは H+を遮蔽するように動くものが、原子
図 4: シリコン固体の構造相転移(ダイヤモンド構
核から離れると H-のような逆造→βスズ構造)に関する量子モンテカルロ計算。
電荷の遮蔽に動かなければな
曲線の極小点から平衡体積が、 共通接線から転移を生らないからと考えられる。
じる圧力やその際の体積変化などが算定される [8)]。

－68－なお、LDA や GGA による Nazarov らの過去の計算[1]によれば、Z が 5-8 で rs>4 では
スピン偏極した解が得られる。これは近藤シングレットの考え方と矛盾するものである
が、電子ガスには本質的に SDW の不安定性があるとのオーバーハウザー理論とは整合的ともいえる。いずれにしても、この領域の Z について、DMC の結果が待たれる。
最後に、前園は非一様系第一原理電子状態を DMC で直接計算することも行っている。
物質デザインにおける複合相関の様相を記述することは古くから理論への大きな挑戦
課題とされてきた。多体電子論の数値的手法である DMC はこの課題に対して最も客観的情報を与える強力なツールであり、大規模高速計算の進歩により、近年、その重要性を益々増しつつある。電子物性における高次の電子相関効果の解明に資する前園の計算結果の一例は図 2 に示されている。
３．GWΓ法：一次元電子系の一粒子励起スペクトルの解析
多電子系の一粒子励起スペクトルの計算手法開発は固体物理学の重要な課題の一つである。現在用いられている GW 近似は一粒子励起を自己無撞着に決定する近似手法
であるが、自己無撞着性を捨てた G0W0 近似の方が半導体のエネルギーギャップの実験
値など高エネルギー励起を良く再現することが経験的に知られている。また、GW 近似の大きな欠点の一つとして、一次元電子系の低エネルギー励起である朝永ラッティンジ
ャー液体を記述できないという原理的な問題がある。本研究の中心課題である（自己無撞着な）GWΓ法は、朝永ラッティンジャー液体を含む金属の低エネルギー励起を記述できるよう定式化されており[2,3]、超伝導状態への応用も考察されている[10)]。
本年度は、GWΓ法が一次元電子系の低エネルギー励起だけでなく、高エネルギ
ー励起も良く記述することを示すために、1/2 充填ハバード模型の高温における一粒子スペクトル関数の計算を行った。図 3 にその結果を示す。G0W0 近似
では、上部・下部ハバードバンドに対応する高エネルギー励起構造が見られる
が、GW 近似ではこの構造が消失している。しかしながら、GWΓ法では (自己エネルギー補正の高次の項とΓの相殺効果
図 3：1/2 充填一次元ハバード模型のスペクトルのために) 高エネルギー励起の構造が適関数 A(pF,ω)のエネルギーω依存性。 pF はフェルミ運動量である。点線、一点鎖線、実線はそれぞ
切に再現されることが明らかになった。
れ G0W0 近似、GW 近似、GWΓ法に対応する。
なお、連携研究者の大野かおる（横浜国立大）も、まだ論文にはなっていない
ものの、別の立場からの GWΓ法の定式化を試みている。また、TDDFT 電子励起ダイナ
ミクス・シミュレーションにより C60 と Zn フタロシアニンの電荷分離機構を調べた。
この他、詳細な第一原理計算により、金属触媒粒子によるナノチューブの規則成長メカ
ニズムを明らかにした。さらに、GW+T 行列法によりスピン偏極した Al や Na クラス
ターのダブルイオン化エネルギースペクトルの計算を行い、同様の方法により Na クラ
スターの光吸収スペクトルの計算を行った。さらに、Si、C、SiC 結晶の全電子 GW 計算を行い、芯電子 XPS スペクトルを求めた。
４．高温超伝導体の探索：ホウ素、硬い半導体による高 Tc 超伝導物質探索
今世紀に入り、銅酸化物超伝導体以外に高 Tc を予感させる超伝導体の発見が相次いでいる。高圧でのボロン、ボロンドープダイヤモンド、あるいはグラファイトインター
カレーション、ドープされたフラーレン族などである。これらの特徴は、古典的なフォ－69－

ノン媒介でありながら、従来の予測よりも高い Tc が出現していること、および、従来は超伝導体とはみなされていなかった半導体でも超伝導が出現したことである。軽い半
導体は固いものが多く、それはフォノンの立場からは強い電子フォノン相互作用を意味
し、したがってフォノン媒介の超伝導でも 50K を上回る超伝導出現の可能性が出てくる。本研究では、そのような硬い半導体を中心にした新たな超伝導体の物質探索を行う
ことを目的とする。
具体的なターゲットしては多くの結晶構造を持つホウ素・ホウ素系半導体、グラフ
ァイト、C-B 化合物である。特にホウ素に関しては、これまでホウ素の相図をはじめて予測し、その上で a-ホウ素の高圧での安定性を予測し、かつ実験家との共同研究によりその超伝導が実証されてきた。これに関しては、レビュー論文を著し[5),6)]、この分野への重要な貢献をなした。
もう一点はグラファイトである。従来、グラファイトはインターカレーションの形
で超伝導探索が行われていたが、本研究ではそれとは全く異なる視点、硬い物質という観点から超伝導の可能性を探る。ドープせずに高圧をかけることでキャリアの密度を大きくし、超伝導出現を図ろうとするものである。そのための第一歩として、高圧での電子構造を研究した。図 4 に示されているが、p=20-30 GPa くらいになると、元のフェルミ面に加えて、新たなフェルミ面が出現し、キャリア密度が急激に増加することが計算から示された。これであれば，B ドープのダイヤモンドでのキャリア密度程度となるので、超伝導出現の可能性がある。

図 4： 高圧下グラファイトの電子構造[2)]。

この他、連携研究員の是常隆（東京工業大学）は、第一原理からのTc評価として標準的であるマクミラン公式を使った手法の改良を目指し、エリアシュバーグ方程式に基づいたTc計算プログラムを開発し、それをフラーレン超伝導体に適用して電子格子相互作用の効果について定量的な評価をした。特に、状態密度のエネルギー依存性などを正しく取り込むことによりマクミラン公式からの改善も見られた。また、ナノチューブにおけるキャリアドープの可能性、さらには超伝導の可能性についても議論した。

－70－

実験研究を進めている 連携研究者の秋光純（青山学院大学）は次のような成果を挙げた：①スピン軌道相互作用が誘起する新しいタイプのモット絶縁体Ba2IrO4の合成に成功した。更にこの系にホールをドープすることにより金属・絶縁体転移を観測した。②
Y5T6Sn18(T=Rh, Ir)の比熱測定の結果からこの系の超伝導ギャップが異方的であることを初めて明らかにした。③新超伝導体W5SiB2（Tc～5.8K)を発見し、その超伝導特性を明らかにした。なお、Ba2IrO4が超伝導にならない理由は理論家にとっても今後の課題である。また、上田寛（東京大学）は新奇な強相関物性を示す物質の開発の一環として、強磁性金属・絶縁体転移を示すK2Cr8O16について、放射光X線回折により低温絶縁体相の構造を明らかにした。その結果、低温絶縁体相では電荷秩序 はなく、代わりにCr-Oの結合交替があり、Crの４量体化が観測された。また、バンド 計算からK2Cr8O16はパ
イエルス不安定性を持ち、25％の過剰電子を４つのCrで共有して絶縁体化するというパイエルス機構が明らかになった。
５．参考文献[1] V. U. Nazarov, C. S. Kim and Y. Takada, Phys. Rev. B 72, 233205 (2005).
[2] Y. Takada, Phys. Rev. Lett. 87, 226402 (2001).
[3] S. Ishii, H. Maebashi, and Y. Takada, arXiv:1003.3342.
６．連携研究者
秋光純（青山学院大学理工学部教授）、春山純志（青山学院大学理工学部准教授）、
上田寛（東京大学物性研究所教授）、大野かおる（横浜国立大学工学研究院教授）
斎藤晋（東京工業大学理工学部教授）、是常隆（東京工業大学理工学部助教）、
７．本研究課題における平成２２年度の発表論文と招待講演
発表論文：
1)

Y. Takada, “Theory of superconductivity in graphite intercalation compounds”, inComprehensive Semiconductor Science and technology, Vol.1, ed. by P. Bhattacharya, R.
Fornari, and H. Kamimura (Elsevier, 2011) pp. 410-426.

2)

N. Nakae, J. Ishisada, K. Shirai, A. Yanase, "Fermi surface and metallic properties of graphite
at high pressures", J. Phys. Chem. Solids 71 (2010) pp. 418–422.

3)

N. Nakae, J. Ishisada, H. Dekura and K. Shirai, "Theoretical investigation on synthesizing
BC5 crystal", J. Phys.: Conf. Ser. 215 (2010) 012116 (7 p).

4)

H. Dekura, K. Shirai and A. Yanase, "Metallicity of boron carbides at high pressure", J.
Phys.: Conf. Ser. 215 (2010) 012117 (6 p).

5)

K. Shirai, "Electronic Structures and Its Mechanical Properties of Boron and Boron-rich
Crystals (Part I)", J. Superhard Materials 32 (2010) pp. 205-225.

6)

K. Shirai, "Electronic Structures and Its Mechanical Properties of Boron and Boron-rich
Crystals (Part II)", J. Superhard Materials 32 (2010) pp. 336-345.

7)

H. Werheit, V. Filipov, K. Shirai, H Dekura, U. Schwarz, and M. Armbrüster, "Raman
scattering and isotopic phonon effects in dodecaborides", J. Phys.: Condens. Matter 23(2011) 065403 (26 pp).

8)

Ryo Maezono, N.D. Drummond, A. Ma, and R.J. Needs, “Diamond to beta-tin phasetransition in Si within quantum Monte Carlo”, Phys. Rev. B. 82, 184108 (2010).

9)

J.R. Trail and Ryo Maezono, “Optimum and efficient sampling for variational quantum
Monte Carlo”, J. Chem. Phys. 133, 174120 (2010).

－71－10) H. Maebashi and Y. Takada, Physica C: Superconductivity 470, S975 (2010).
11) R. Cherian, C. Gerard, P. Mahadevan, N.T. Cuong, and Ryo Maezono, “Size dependence ofthe bulk modulus of semiconductors nanocrystals”, Phys. Rev. B. 82, 235321 (2010).
12) Y. Kodama and K. Ohno, “Charge separation dynamics at molecular heterojunction of C60
and zinc phthalocyanine”, Appl. Phys. Lett. 96, 034101;1-3 (2010).
13) J.-T. Wang, C. Chen, K. Ohno, E. Wang, X.-L. Chen, D.-S. Wang, H. Mizuseki, and Y.
Kawazoe, “Atomistic nucleation and growth mechanism for single-wall carbon nanotubes on
catalytic nanopartile”, Nanotechnology 21, 115602;1-5 (2010).
14) Y. Noguchi, K. Ohno, I. Solovyev, and T. Sasaki, “Cluster size dependence of double
ionization energy spectra of spin-polarized aluminum and sodium clusters: All-electron
spin-polarized GW + T-matrix method”, Phys. Rev. B 81, 165411;1-9 (2010).
15) Y. Noguchi and K. Ohno, “All-electron first-principles GW + Bethe-Salpeter calculation for
optical absorption spectra of sodium clusters”, Phys. Rev. A 81, 045201;1-4 (2010).
16) S. Ishii, S. Iwata, and K. Ohno, “All-Electron GW Calculations of Silicon, Diamond, andSilicon Carbide”, Mater. Trans. 51, 2150-2156 (2010).
招待講演：
1)

Y. Takada, “On the first-principles determination of Tc”, Superstripes 2010 Quantum
Phenomena in Complex Matter (Erice, 19-25 July, 2010).

2)

Y. Takada, “Theory for reliable first-principles prediction of the superconducting Tc”, 13th
Asian Workshop on First-Principles Electronic Structure Calculations (Pohang, 1-3
November, 2010).

3)

K. Shirai and H. Dekura, “Superconductivity research on semiconducting boron”, 14th
International Conference on High Pressure Semiconductor Physics (Jilin, China, August 1-4
2010).

4)

Ryo Maezono, “Binding of chromium dimer studied by QMC”, PACIFICHEM 2010 (Hawaii,
December 15-20, 2010).

5)

Ryo Maezono, “Electronic structure calculation using Quantum Monte Carlo technique”,
JAIST-CNSI Workshop 2011 (Los Angeles, January 13-14, 2011).

6)

J. Akimitsu, “New superconductors recently developed by our group”, The 9th Asia PacificWorkshop on Materials Physics (Hanoi, December 2010).

7)

J. Akimitsu, “MgB2, its application and new superconductors recently developed by ourgroup”, 2011 APCTP Winter Workshop on Frontiers in Electronic Quantum Matter (Pohang,
February 2011).

8)

上田寛、“Hollandites: novel metal-insulator transitions”, 日本物理学会2010年秋季大会シンポジウム招待講演 (大阪府立大学中百舌鳥キャ ンパ, 2010/9/23-26).

9)

大野かおる、“材料シミュレーション”、神戸大学 文科省科学技術振興調整費事業「企業を牽引する計算科学高度技術者の育成」第１回シミュレーションスクール（2010年10月25日、神戸ポートアイランド）.

10) T. Koretsune and S. Saito “Electronic Structures of Boron- or Nitrogen-doped Carbon
Nanotubes”, Int. Conf. on Superconductivity and Magnetism (ICSM) 2010 (Antalya, April
25-30, 2010).

－72－

スピンエレクトロニクス材料の探索
Design and Realization of Spin-electronics Materials
佐藤和則 1、小田竜樹 2、小倉昌子 1、野崎隆行 1
1

K. Sato, 2T. Oda, 1M. Ogura and 1T. Nozaki
大阪大学 1、金沢大学 2

Osaka University, Kanazawa University

半導体スピントロニクス材料として重要な高いキュリー温度を持つ磁性半導体のデザインを行い同時ドーピング法を提案した。また、磁性半導体の精密な物性評
価のために磁気光学効果の第一原理計算法を開発した。電圧による磁気異方性制
御は超低消費電力な磁化書き込み法の基盤技術として、金属磁性薄膜における磁気異方性エネルギーとその電界効果に関する研究を行い Pd/Fe/Pd(001)における実験結果を定性的に説明した。また、実証実験として Ag / 超薄膜 Fe / MgO をベースとした接合構造において電圧印加による超薄膜 Fe の磁気異方性変化を調べ、Au バッファーを用いた従来構造の約 3 倍の電圧効果の増大に成功した。

１． はじめに	  超高速•超高集積•超省エネルギーの次世代エレクトロニクスとしてスピンエレクト
ロニクスが注目されている。本研究ではスピントロニクス実現への大きなブレークスルーを目指して半導体および金属系スピントロニクス材料の第一原理マテリアルデザイ
ンとその実証実験の連携研究を行う。H22 年度は、(1)相分離を抑制し高濃度添加を可能にする同時ドーピング法のデザイン、(2)磁気光学効果の第一原理計算法の開発とそれを用いた磁性半導体の精密な物性評価、(3)Fe 磁性層や MgO/Fe 接合界面における磁気異方性制御および電界効果増大のための指針提案、(4) Ag / 超薄膜 Fe / MgO / Polyimide/ ITO 接合での垂直磁気異方性、およびその電圧印加効果増大の実証実験を実施した。	
 
２． (Ga, Mn)As における同時ドーピング法のデザイン
	  半導体材料をベースとする半導体スピントロニクスでは磁性半導体が有望とされているが室温より高いキュリー温度をもつ磁性半導体の開発とが不可欠である。過去のマテリアルデザインによるとキュリー温度上昇には磁性不純物の高濃度添加が必要であるが相分離のために実現は難しい。例えば(Ga, Mn)As の場合、熱平衡状態では GaAs と
MnAs に分離した状態が安定である。均一に混合した状態と分離した状態の安定性は、混合
エ ネ ル ギ ー ΔE を 計 算 す る こ と で 予 測 で き る 。 例 え ば (Ga, Mn)As に つ い て ΔE は 、
ΔE=TE(Ga1-xMnxAs)-[x TE(MnAs)+(1-x) TE(GaAs)]と定義される。ここで TE は全エネ
ルギーであり計算には赤井の開発した MACHIKANEYAMA2002 パッケージを用いた。
計算の結果(Ga, Mn)As の場合 Mn 濃度の全域にわたって混合エネルギーは正となり均一に固溶した状態は不安定であることがわかった。しかしここに Li を同時に添加すると状況は一変し、Mn の低濃度領域に混合エネルギーが負になる均一添加が可能な領域が現れる。図1 は計算で得た混合エネルギーに配置のエントロピーの効果を加え自由エネルギーの計算から相図を予測した結果である。Li の添加により Mn の固溶限が大きく上昇しているのがわかる。
	  磁性半導体はキャリア誘起型の強磁性を示すので、格子間 Li の同時ドーピングにより正孔が補償されるため、高濃度に Mn を添加できたとしてもこのままでは強磁性転移温度は低－73－いと予想される。しかし格子間 Li の拡散係数が非常に大きいため低温でのアニーリングにより取り除くことができ、その結果高い強磁性転移温度が実現されると期待できる。このことを定量的に示すため、PAW 法により Li の拡散障壁や Li-Mnの結合エネルギーを見積もり、これらのパラメータを使い、動的モンテカルロシミュレーショ
ンから Li の(Ga, Mn)As 中での拡散係数を計算した。Li の場合は、Mn との引力がそれほど大きくなく、現実的な温度とアニーリング時間で薄膜表面まで拡散できることを示した。

図 1（右図）: Li を同時ドープした(Ga, Mn)As の相図。赤線はバイノーダル線、緑線はスピノーダル線を示す。
３． KKR-CPA 法と久保公式による磁性半導体の光吸収 MCD の計算
	  磁気円二色性(MCD)や磁気 Kerr 効果などの磁気光学測定は希薄磁性半導体の電子状態と磁性を調べるのに強力な手法である。しかし、これまでの磁性半導体の磁気光学効果に関する議論
は主に s-d モデルに基づく定性的なもので、必ずしも実験結果を説明できるものではな
かった[1]。本研究では KKR グリーン関数法で求められる 1 粒子グリーン関数から久保公式を用いて光学伝導を第一原理から計算し、磁性半導体について数 eV 領域の光吸収
MCD を調べた。光学伝導の吸収部分は Butler による DC 電気伝導の表式[2]と同様の形で書き表すことができる。
absσ µν(ω ) = lim

η →0

ε
1
∫ F dε [σ˜ µν (ε + iη,ε + ω + iη) + σ˜ µν (ε − iη,ε + ω − iη)4ω ε F −ω
− σ˜ µν (ε + iη,ε + ω − iη) − σ˜ µν (ε − iη,ε + ω + iη)]σ˜ µν (z1,z2 ) = −Tr j µ G(z1 ) jν G(z2 )πV

j µ = −i

e 
∇µ
m

	
  計算の結果得られた MCD スペクトルは実験で得られたスペクトルの特徴をよく再現し、MCD スペクトルに関して磁性半導体は以下の３つに分類することができることがわかった。(1)磁性不純物がホールを作らない、つまりフェルミレベルがバンドギャッ€
プ中にある。(2)片方のスピン状態にだけホールが存在する。(3)両方のスピン状態にホールが存在する。
	  吸収係数は主に光学伝導 σxx(1)で決まるのに対して、MCD は σxy(2)で決まる。だから
MCD の起源を議論するためには σxy(2)のふるまいを調べればよい。磁性半導体で光吸収に寄与するのは主に価電子帯である p 状態から伝導帯である s 状態への励起である。p状態はスピン軌道相互作用によって up スピン状態ではエネルギーの低い方から m = +1,
0, -1 のように分裂する（down スピンでは逆）。m = -1 状態から s 状態への励起は負の、
m = +1 状態からの励起は正の σxy(2)を作るので、up スピン電子の励起は吸収端付近で負の σxy(2)を作り光子エネルギーが上がるにつれて σxy(2)は正に変化する。down スピンでは逆である。このままでは up と down の寄与が打ち消し合って MCD はゼロになってしま
うが、磁性半導体では up スピンの p 状態は磁性不純物の d 状態と強く混成して down
スピンより広いバンド幅を持つ。その結果、MCD スペクトルは正と負が交互に現れるような形状を持つ。以上が(Cd1-xMnx)Te などの(1)の磁性半導体での MCD が生じるメカ－74－

ニズムである。
	  (2)の場合では、(1)の場合に加えて価電子帯からホールへの励起を考えなければならない。up スピン電子の価電子帯からホールへの励起は低光子エネルギー領域に正のMCD を作る。そのために MCD スペクトルは(1)に比べて全体的に底上げされたような形を持つ。このようなスペクトルは(Zn1-xCrx)Te などの強磁性半導体で見られる。(3)の場合は down スピン電子のホールへの励起もあるのでさらに複雑になる。
	  (Ga1-xMnx)As の MCD スペクトルも同様に up スピン状態のホールの変化から説明が可能である。(Ga1-xMnx)As は up スピン状態にだけホールが存在するハーフメタルであるが、通常 GaAs は antisite As などの不純物により n 型になっているため少量の Mn では
ホールは形成されない。そのため、Mn 濃度の低いサンプルでは(1)の MCD スペクトルが得られ、Mn 濃度が増加するにつれて(2)のスペクトルに変化する。理論計算では(1)の場合は強磁性が発現しないが、実験でも(1)の MCD スペクトルを示すサンプルは常磁性である。このことから、磁性半導体の強磁性発現にはホールの存在が関係していることが強く裏付けられる。
４． 金属磁性薄膜における磁気異方性とその電界効果	  金属系スピントロニクスのデザインについては、外部電界下での磁性薄膜の磁気状態を、密度汎関数理論に基づいた相対論的 2 成分擬ポテンシャル第一原理電子状態計算法
を用いて研究した。薄膜に印加する電界毎に、 (001)と(100)の 2 方向の磁化に対する全エネルギー計算の差を求めることで磁気異方性エネルギー（MAE）を見積もった。金属 薄 膜 ( 真 空 /M/Fe/M(001), M=Pt, Pd) 、 接 合 薄 膜 (MgO/Fe/M(001), M=Au, Pt) 、
MgO/Fe/Au(001)について計算を行った。
（A）金属薄膜(真空/M/Fe/M(001), M=Pt, Pd)
FePt や FePd 規則合金薄膜を用いることにより、強い面直な表面磁気異方性が得られる
ことはよく知られているが、電界効果について、これら２つの系は逆の効果を示す[3]。そこで理論計算で Pd 系と Pt 系の電界効果がどの程度再現されるかを確認した。
図 2(a)は計算で用いたスラブモデルを示している。スラブの下３層の Pd はバルクモデルの格子定数に固定し最適化を行った。層間距離は 0.176, 0.179, 0.206 nm(それぞれ Pd(c)–Fe,
Fe–Pd(1), Pd(1)–Pd(2) )となった。これに対して Pt 系では 0.177, 0.184, 0.213 nm となり、よく似た傾向を示している。図 2(b)は、静電ポテンシャルと電界誘起の電子密度およびスピン分
解電子密度を示している。外部電界に対するス
クリーニング効果が良く再現されており、電界誘起されるスピン密度が表面３原子層まで振動しながら減衰している。
図 3 は Pd 系および Pt 系の MAE 電界効果を示している。電場 0 における MAE は Pd 系で‐0.43
mJ/m2、Pt 系においては 7.6mJ/ m2 であった。それ ぞ れ の 電 界 に 対 す る 傾 き は γ Pd(Pt)= ‐21(78)fJ/Vm と得られ、γPd/γPt=‐0.27 と見積もられた。後者の傾きの比については、実験[3]と半定量的に一致した。MAE 密度[4]を計算した結果、やはり電界効果は表面付近で見られた。Pd
図 2：(a)モデル、(b)電界ε=‐9.6(緑),系と Pt 系に対する電界効果γPd(Pt)の比が、実験
0.0(赤), 9.6(青)V/nm の時の静電ポテンシ
結果を再現していることは、MAE 電界効果のデャル(b 上部)と誘起される電子密度(b 中
ザインに明るい見通しを与える。
部)、電場 9.6V/nm の時の誘起される多数
Pd 系と Pt 系で電界効果γPd(Pt)の符号が異なっ
スピン(マゼンダ)と少数スピン(シアていることは、フェルミ準位近傍での電子状態
ン)(b 下部)の詳細に起因するが、決定的な原因を明確にす－75－

ることはできなかった。Pt 系では、非占有のFe d(xz,yz)状態が、内向き電界の増大に伴い成長し、占有状態 Fe d(x2-y2)とのスピン軌道結合の寄与が増大し、その結果 MAE が減少することが明らかとなっている。一方、Pd 系では、
Fe と Pd の軌道混成が弱く、フェルミ準位近傍の電子構造がより複雑であることから、Pt 系で得られた原因以外に、いくつかのものが挙げられている。これらが逆の電界効果を生み出し
図 3：Pd/Fe/Pd(GGA)(赤)、Pt/Fe/Pt(GGA)ている。
(青)および Pt/Fe/Pt(LDA)(緑)における外（B）接合薄膜(MgO/Fe/M(001), M=Au, Pt)
部電界に対する磁気異方性エネルギー変	  強い面直磁気異方性を持つ Pt/Fe/Pt(001)薄化
膜においては、これまでに電界効果についても Fe 層だけの薄膜に比較して数倍程度大きくなる可能性が示唆されている[4]。Pt 層をもつ薄膜での大きな電界効果を明確に示す
ために、MgO/Fe/M(001) (M=Au, Pt)の２種類の薄膜についての比較研究を行った。Pt 基板の系では Au 基板の系に比べ電界効果が増大しており、電界に対する MAE 変化の傾きは
63fJ/Vm と得られた。この値は、Pt/Fe/Pt(001)薄膜(72fJ/Vm)と同程度であり、孤立 Fe 薄膜等に比べ数倍程度大きな値である。Pt 基板の系における電界効果の増大は、電界によりFe3d 軌道と Pt5d 軌道の混成が変化することによると考えられる。	
 ５． Ag/超薄膜 Fe/MgO 接合における電圧誘起磁気異方性変化
	  超薄膜 Fe 層の垂直磁気異方性制御の実証
実 験 と し て 、 Au / 超 薄 膜 Fe / MgO /Polyimide / ITO 接合における電圧誘起磁気
異方性変化の効果増大を目指し、Ag バッフ
ァーを用いた構造における超薄膜 Fe の垂直
磁気異方性とその電圧印加効果を調べた。
	  分子線エピタクシー法を用いて MgO(001)
基板上に Cr / Ag バッファー (50 nm) / 超薄膜 Fe (tFe: 0.5~1 nm) / MgO (10 nm)からなる多層膜構造を成長させた。全て室温成長である
が、Ag バッファー層のみ成長後に①アニール無し、②200℃、③400℃の３種の条件
図 4:	
  Ag バッファーアニール 400oC、tFe=0.9 nm
で処理を行った。試料を大気暴露した後、 試料における垂直磁気異方性への電圧印加効果
スピンコートにより有機絶縁体である
Polyimide 層(1.5 µm)を塗布し、200℃ベー
クで硬化させた。最後に ITO 透明電極(直径 1 mm)をメタルマスク法により作製した。
超薄膜 Fe の磁化曲線は垂直磁場印加下での極カー効果測定により評価した。
Ag バッファー層試料においても、Au バッファー層の場合と同様に界面磁気異方性に起因する垂直磁気異方性が観察された。Au バッファーの場合は 0.5 nm 以下の超薄膜領域においても磁化容易軸は面内であったのに対し、アニール無しの Ag バッファー試料
では 0.7 nm 程度の比較的厚い膜厚領域においても垂直磁化膜が得られ、980 mJ/m2 に及ぶ大きな界面磁気異方性を有することが確認された。同試料において電圧印加効果を調べた。図 4 は今回作成した試料の内、最も大きな異方性変化が観察された実験例である。
Ag バッファー層は 400℃アニール処理しており、Fe 膜厚は 0.9 nm である。 200 V の電圧印加により、明瞭な垂直磁気異方性の変化が観察された。磁化曲線から見積もられる異方性変化量は 31.3 kJ/m3 であり、これは電界印加に対する異方性変化率としては 312－76－

fJ/Vm に相当する。この値は Au / 超薄膜 Fe / MgO 接合系と比較して約 3 倍の増大である。増大の原因には超薄膜 Fe の成長形態が関係していると考えているが、STM 観察などを用いたより詳細な検討が必要である。
６． まとめ	  半導体スピントロニクス材料のデザインとして、磁性半導体へ高濃度に磁性不純物を添加する方法を提案した。格子間 Li を含む(Ga, Mn)As の相図を第一原理から計算し Liの同時ドーピングで Mn の固溶度が上昇することを示した。また結晶成長後に Li を取り除く方法として低温アニーリング法をデザインした。KKR-CPA 法と久保公式を用いて磁性半導体の MCD スペクトルの第一原理計算法を開発し CdMnTe, ZnCrTe や
GaMnAs 等の代表的な磁性半導体の MCD スペクトルの定量的な予測に成功した。
	
  金属系スピントロニクスのデザインについては、外部電界下での磁性薄膜の磁気状態を、密度汎関数理論に基づいた相対論的 2 成分擬ポテンシャル第一原理電子状態計算法
を用いて研究した。金属薄膜の電界効果においては、実験結果を定性的に説明し、接合薄膜を用いて MEA 電界効果を大きくするための指針を検討した。Ag / 超薄膜 Fe / MgO
接合系における電圧印加磁気異方性変化の実証実験を行った。Au バッファー構造と比較して垂直磁気異方性が増大し、電圧効果も約 3 倍にエンハンスされることを見出した。
７． 参考文献[1] K. Ando, H. Saito, K. C. Agarwal, et al., Phys. Rev. Lett. 100, 067204 (2008).
[2] W. H. Butler, Phys. Rev. B 31, 3260 (1985).
[3] M. Weisheit et al., Science, 315, 349 (2007)[4] M. Tsujikawa and T. Oda, Phys. Rev. Lett., 102, 247203 (2009)[5] T. Maruyama et al. Nature Nanotechnol. 4, 158 (2009).
[6] Y. Shiota et al. Appl. Phys. Express 2, 063001 (2009).
８． 連携研究者
福島鉄也（大阪大学基礎工学研究科特任研究員）、黒田眞司（筑波大学物質工学系教授）、吉田博（大阪大学基礎工学研究科教授）、朝日一（大阪大学産業科学研究所教授）、鈴木義茂（大阪大学基礎工学研究科教授）、赤井久純（大阪大学理学研究科教授）、下司雅章（大阪大学ナノサイエンスデザイン教育センター特任講師）

９． 本研究課題における平成２２年度の発表論文と招待講演
発表論文：
1) K. Sato, L. Bergqvist, J. Kudrnovsky, P. H. Dederichs, O. Eriksson, I. Turek, B. Sanyal,
G. Bouzerar, H. Katayama-Yoshida, V. A. Dinh, T. Fukushima, H. Kizaki, R. Zeller,“First-principles theory of dilute magnetic semiconductors”, Rev. Mod. Phys. 82 (2010)
1633-1690.
2) N. D. Vu, K. Sato and H. Katayama-Yoshida, “Giant Peltier effect in self-organizedquasi-one-dimensional nanostructure in CuNi-alloy”, Appl. Phys. Express 4 (2011)
015203 (3pages).
3) L. Bergqvist, K. Sato, H. Katayama-Yoshida and P. H. Dederichs, “Computational
Materials design for high-TC (Ga, Mn)As with Li-codoping”, Phys. Rev. B in printing.
4) H. Fujii, K. Sato, L. Bergqvist, P. H. Dederichs and H. Katayama-Yoshida, “Interstitialdonor codoping method in (Ga, Mn)As to increases solubility limit of Mn and Curie
temperature”, Appl. Phys. Express in pringing.

－77－

5)

M. Tsujikawa, S. Haraguchi, T. Oda et al., “A Comparative ab initio study onelectric-field dependence of magnetic anisotropy in MgO/Fe/Pt and MgO/Fe/Au films”,
J. Appl. Phys. in printing.

6)

S. Haraguchi, M. Tsujikawa, J. Gotou and T. Oda, “Electric-field effects on magnetic
anisotropy in Pd/Fe/Pd(001) surface”, J. Phys. D: Appl. Phys. 44 (2011) 064005(8pages).

7)

M. Ogura and H. Akai, “First-principles calculation of magnetic circular dichroism
spectra of magnetic semiconductors”, Phys. Rev. B 82 (2010) 184426 (9 pages).

招待講演：
1) K. Sato, “First-principles theory of dilute magnetic semiconductors”, International
Conference on Core Research and Engineering Science of Advanced Materials (Osaka,
Japan, 30 May – 4 Jun. 2010).
2) K. Sato et al., “Materials design of dilute magnetic semiconductors by controlling the
spinodal decomposition”, International Conference on Core Research and Engineering
Science of Advanced Materials (Osaka, Japan, 30 May – 4 Jun. 2010).
3) K. Sato et al.,“Computatilnal Materials Design of Spinodal Nanotechnology”,
International Union Materials of Research Societies – International Conference on
Electronic Materials (IUMRS-ICEM) (Seoul, Korea, 22-27 Aug., 2010).
4) K. Sato, “Spinodal decomposition and magnetic properties of dilute magnetic
semiconductors”, European Materials Research Society (EMRS) fall meeting, Warsaw,(Poland, 13-17 Sep., 2010).

5)

M. Ogura, “First-principles calculation of the magneto-optical effects of magnetic
semiconductors”, International Conference on Core Research and Engineering Scienceof Advanced Materials (Osaka, Japan, 30 May – 4 Jun. 2010).

6)

T. Oda, “Electric-Field Effect on the Magnetic Anisotropy of the
Ferromagnetic/Dielectric Films: A First-principles Study”, International Symposium on
Advanced Magnetic Materials and Applications (ISAMMA) (Sendai, Japan, 12-16, Jul.
2010

7)

T. Oda, “Toward a Computer Modeling in Magnetic Anisotropy and its
Electric-Field-Control for Nano-Structures”, European Materials Research Society(EMRS) fall meeting, Warsaw, (Poland, 13-17 Sep., 2010).

8)

T. Oda, “Magnetic anisotropy and its electric field effect in the nano-structures”, The
13th Asian Workshop on First-Principles Electronic structure Calculations (Asian13)(Pohang, Korea, 1-3 Nov. 2010).

9)

小田竜樹，``ナノ構造の磁気異方性と電界効果"，PF研究会「磁性薄膜・多層膜を究める：キャラクタリゼーションから新奇材料の創製へ」（KEK研究本館小林ホール, つ
くば市，2011年3月11-12日)

10) T. Nozaki, “Voltage-control of perpendicular magnetic anisotropy in fully epitaxial magnetic
tunnel junctions”, SPIE 2010 NanoScience+Engineering (San Diego, 1~5 Aug. 2010)11) 野崎隆行	
  “Fe/MgO接合系における電界誘起磁気異方性制御の最近の進展”, 東北大学電気通信研究所	
  共同プロジェクト研究会 (電気磁気および磁気弾性効果の計算
機物質設計とデバイス応用) (仙台, 2010年12月17日).
12) 野崎隆行	
  “電界による磁気異方性制御の試み”磁気記録・情報ストレージ研究会(MR) (大阪, 2011年1月20日)－78－