◆

プロジェクト報告「機能元素のナノ材料科学」◆

ナノ機能元素のマルチスケールシミュレーション
鶴 田 健 二小 山 敏 幸小林亮1)

尾 形 修 司2 兵 頭 志 明)
岡山大学准教授大学院自然科学研究科(〒7008530 岡山市北区津島中 311)
物質・材料機構計算科学センター主幹研究員
名古屋工業大学大学院工学研究科 1)特任研究員(現・助教) 2)教授
豊田中央研究所計算物理研究室室長
Multiscale Simulations for Atomicscale Modification of Nanomaterials; Kenji Tsuruta, Toshiyuki Koyama, Ryo KobayashiandShuji Ogata, Shiaki Hyodo(Graduate School of Natural Science and Technology, Okayama University, Okayama, Computational Materials Science Center, National Institute for Materials Science, Tsukuba, Graduate School of Engineering, Nagoya Instituteof Technology, Nagoya, Toyota Central R&D Labs, Aichi)
Keywords: multiscale simulation, molecular dynamics, coarsegrained particle, phasefield method, hybrid densityfunctional/classical model,
dislocation, grain boundary, segregation, nanoindentation
2009年 1 月17日受理

.

は じめに

ナノ材料・ナノデバイスの機能発現の解明において，原子レベルの微小領域で起きる現象の詳細と系全体の巨視的な振る舞いとを切り離した理論・シミュレーションでは，その発現メカニズムを十分に捉えられない場合がある．例えば，ナノ多結晶体の力学物性や電気的特性の主要な部分は結晶粒界の方位関係だけでなく，そこに偏在する不純物の元素種，位置，周辺原子との結合様式などの詳細に強く影響を受ける(1)．これらの微視的な様相が材料の特性に如何に影響するかを定量的に解明し，その知見から逆に所望の材料物性値を得るための添加元素(“機能元素”)を，試行錯誤を経ずして選び出すことが本研究課題の目標である．この目標に対し，電子・原子レベルからのボトムアップ的理論・シミュレーショ
図1ンを考える場合，偏析元素周辺の領域にのみ電子状態の変化

マルチスケール計算体系の概念図．

をあらわに取り入れた高精度計算手法を適用し，その他の領域は経験的相互作用模型に基づく古典的分子動力学法を適用

質を粗視化粒子法により記述する，あるいは歪分布などからすれば，サブナノ(Å)～ナノスケールの物理現象の主要要素

直接に系の熱力学的時間発展を PhaseField 法でシミュレーを押さえられると期待できる．一方，材料の疲労や組織変態トする．本稿では，各スケールで試みられている欠陥・転などは，原子レベルからみればずっとマクロかつ超長時間現位・粒界の微視的構造・偏析過程解明への適用と，スケール象であるので，系を熱力学・統計力学的に粗視化した記述が統合へむけた試みを紹介する．

必要である．

.

我々の研究チームは，科研費・特定領域研究「機能元素のナノ材料科学」のメンバーとして，複数のスケールを記述する計算手法をそれぞれ高度化，協調・統合し，ナノ領域に局

方法論

ハイブリッド分子動力学法

在する機能元素がマクロ物性に及ぼす影響を定量的に予測すハイブリッド分子動力学(MD)法は精度の異なる複数の手る新しい計算材料科学手法の枠組みを構築しようとしてい法を階層的に適用し，異なるスケール間でそれぞれに必要なる．図に，本研究で構築を目指すマルチスケール計算体系精度のレベルを段階的に結びつける方法の一つである．例えの模式図を示す．粒界などのナノ領域に偏析する不純物元素

ば，図で模式的に示されるように，材料中の欠陥や不純物がその周囲でどのような結合状態にあるかを大規模な電子状などの局在領域の原子間力・安定構造は電子状態計算を含む態計算法で評価し，それによる周囲の原子変位・歪分布など高精度計算で，その他の領域の原子配置・運動を経験的な古の変化を古典分子動力学法で求め，その原子集団の統計的性典的原子間相互作用で記述する．このハイブリッド系のハミま て り あ 第48巻 第 6 号(2009)Materia Japan図2

図3

ハイブリッド QM/MD 法の模式図．

ルトニアンを次のように 3 項の和として定義する(2)(3)
H ＝Hsystem
MD

粗視化することにより得られ，これを通して MD 領域からの波動を CGP 領域に伝え，CGP 領域から入射する波動には({r… i, r…_ i; i∈total })＋E cluster… i; i∈QM＋Handshake})
QM ({r－Ecluster
CL({r… i; i∈QM＋Handshake}).

原子と粗視化粒子との接続の概念図(各シンボ
ル，記号の意味は本文参照のこと)．

無反射条件を実現する．一方，仮想原子の変位は，ニュート(1)ン方程式に従う時間ステップ毎に，右辺第 1 項は古典相互作用で計算される全系のハミルトニ

NCG

ãu(i e)＝ai u(i e)＋∑NIiアン，第 2 項は高精度に評価したい“量子(QM)クラスター領域”の電子状態計算によるポテンシャルエネルギー，第 3I(U －a ∑ f u )NatomIi

Ij j(2)jの変更を受ける．ここで Natom, NCG はそれぞれ原子数およ

項は第 1 項と同じ古典モデルで計算する QM 領域のポテン

び CG 粒子数， fIi は粗視化の度合いを決定する重み関数，シャルエネルギーである．次章で例示するアルミナ結晶中のNIi は CG 粒子変位から原子変位を決定する内挿関数であ転位芯構造の解析では，QM 計算には実空間表式の密度汎関る．式( 2 )において，仮想原子の変位の長波長成分は CGP数法を適用し，その他の領域は Vashishta らの経験的モデ系から内挿のものに置き換えられるが，短波長成分は減衰因

ル(4)を用いる．この多階層化手法はクラスター計算に妥当性子 ai を掛けて徐々に減衰する．これにより，CGP 領域からがある対象においてユニバーサルな方法論であり，次節で概の波動は MD 領域へ伝えられ， MD 領域からの波動に対し説する粗視化粒子法と MD 法の接合にも適用される．計算ては仮想原子の領域で長波長成分のみ透過，短波長成分は減の詳細は文献(2)(3)を参照されたい．

衰させる無反射条件を実現する．



粗視化粒子法PhaseField 法

粗視化粒子( CGP )法(5) の基本的考え方は，原子集団のあPhase Field ( PF )法(9)(10) は，現象論的な連続体モデルにる平均変位を粗視化粒子の変位で表し，フォノン近似のもと基づき全自由エネルギーを基礎に複数の発展方程式を数値解での熱統計平均を通じて系のエネルギー期待値を求め，析することによって，組織形成過程を計算する手法である．

CGP 間の実効バネ係数を得ることにある(つまり，バネ定数計算対象の組織形態を表現する秩序変数を定義し，この秩序は方向および粒子ペア毎に異なる)．CGP 間ポテンシャルを変数にて，組織の有する全自由エネルギーを書き下し，その用いることで，少ない計算コストで精度良く材料の弾性挙動全自由エネルギーが最も効率的に減少するように秩序変数のを再現できる．フォノン近似の範囲ではあるが，低温極限の時間空間変化を記述する発展方程式が定義される

固体に関して， CGP 系のエネルギーが原子系の基底エネル&ci(r, t)dGsys& s j (r , t )dGsys＝;･ Mc;＝－Ls,.
&tdc i&tds j

ギーと正確に一致する(6)．Kobayashi らはこれまで，再帰的{}(3)操作による粗視化計算の高効率化(6)，さらに局所的な回転もci(r, t)と sj(r, t)はそれぞれ位置 r および時間 t における濃度考慮して複雑構造体の大規模変形にも適用可能にし(6)，流体場および phase field を表す秩序変数，i と j はそれぞれの秩中の微細材料の振動問題に応用した(7)．

序変数の番号(成分および相の種類に対応)である．Mc と LsMD

は，各々の秩序変数の時間変化に対する易動度である．Gsys領域と CGP 領域を適切に接続する方法の開発が必須であは相変態組織全体の全自由エネルギーで，化学的自由エネルる．特に，波動伝播などを伴う過程の場合，二つの領域で表ギー，勾配エネルギー，弾性歪エネルギー，および電磁気エ現できる最小の波長が異なるので，MD 領域で発生した短波ネルギーの総和にて与えられる．上式を数値計算することに長の波は，そのままでは CG との境界で非物理的な反射をおより，組織形成過程をシミュレートする．この PF 法には，こす．これを抑制する方法がこれまでにいくつか提案されて材料における各種の格子不整合領域(転位，粒界，相境界ないる(8)が，それぞれに長短所がある．我々はそれらの欠点をど)の形態形成，および相変態・析出における元素の拡散挙改善し，MD 領域形状の選択自由度が高く，計算コストが低動などを，総合的かつ定量的にモデル化できる利点がある．

本研究の目標であるマルチスケール化においては，い接続方法を提案している．

.

図に提案接続法における MD 領域と CGP 領域の境界設定の模式図を示す． MD 系， CGP 系共に，境界から奥にも仮想の原子，粒子が存在している．実在の原子および粒子は，それぞれ白丸および白四角，変位には上付きの(r)を付

適用例

ハイブリッド QM / MD 法による a Al2O3 転位芯構造と電子状態

け，同様に仮想の原子および粒子は黒丸および黒四角，変位アルミナ(Al2O3 )は構造材料としても電子材料としても用には上付きの(e )で表す．仮想 CG 粒子の変位は， MD 系をいられる代表的な酸化物セラミックスであるが，その構造の

プロジェクト報告複雑さゆえに，金属・半導体に比べて転位構造・挙動に不明な ら び に MD と CGP の ハ イ ブ リ ッ ド 計 算 を 行 っ た 場 合な点が多い．図のような basal 転位芯(バーガーズベクト(Hybrid MD/CGP)についてそれぞれプロットしたものであ〈 112̃0 〉)において，余剰原子半面が Al 原子列で終端
ル＝ 1/ 3る． Full CGP ではインデンター周辺の大きな変形による非するか，O列かによって転位線の有効電荷が変わり，その安調和効果やボンドの切断を表現できないのに対し，インデン定性と運動のメカニズムが大きく異なってくる(11) ．高分解

ター周辺に MD 領域を採用した Hybrid MD / CGP の結果は

能電子顕微鏡の進歩により，最近になってそれぞれに荷電すFull MD の結果とほぼ一致する．

る部分転位がペアとなって安定に存在することが初めて観察

実際の実験ではグラフェン・ナノドラムの直径 1 mm 程度され(11)(12)，現在，その安定性を裏付ける原子レベルの解析で，単原子膜であっても原子数は 3000 万程度にもなり膨大が必要とされている．我々は，ハイブリッド QM / MD 法にであるが，例えば100原子を 1 粒子に粗視化するハイブリッより，aAl2O3 結晶中転位芯の安定構造ならびに，転位芯付

ド CGP/MD 計算を行えば，現実的な計算量にまで軽減できる．

近の局所電子状態を求めた(13)．図 4(a)は本手法で求めた転

また，CG 粒子の運動を特徴付ける時間スケールが原子スケ位芯ペアの安定構造図と，ハイブリッド計算を行った QMールよりも長い事から，系のダイナミクスを比較的長い時間領域の拡大図である．各転位芯は正負両方の部分転位の近接ステップで行える．上記のインデンテーション破壊の場合，ペアを想定し， MD 緩和において転位芯表面付近の O 原子

Full MD に対し約30倍早く計算結果を得ることが出来た．

が大きく変位する表面再構成が起きることがわかった．図4 ( b )は転位芯中心を原点とした各距離にある原子の有効電荷分布のバルク値からの変化を表す．原子位置の変化は O



転位・粒界析出の PF シミュレーション

図(a)は PF 法を用いた Fe10 atCu の 823 K 等温時効原子が大きかったのに対し，電荷分布の変動はむしろ Al 原における Cu の転位上析出過程のシミュレーションである．

子が大きい．これは，Al 原子の O への結合ボンドの方位依

黒さの度合いが局所的な Cu 濃度を表し，中央の白色の横線存性が強く，表面再構成によるボンド歪みの局所電子状態への部分は，原子面を 1 枚余分に挟み込んだ位置を示していの影響が大きいことを反映していると考えられる．現在，こる．つまり刃状転位が白線の両端に存在し，転位線は紙面にの転位芯への不純物添加効果を解析中である．

垂直で，バーガースベクトルは図の上下方向を向いている．



ハイブリッド MD/CGP 法によるグラフェン・ナノドラムのインデンテーション破壊

初期状態は bcc 相の過飽和固溶体で，相分解の初期に刃状転位の外側に Cu は濃縮し始め，逆に刃状転位側の Cu が枯渇する．その後，時効の進行に伴い Cu 相が成長していく．

図( a )は，前章で概説したハイブリッド MD / CGP 法にCu の原子サイズが母相の Fe の原子サイズよりも大きいのよるグラフェン・ナノドラムのインデンテーション破壊におで，以上の挙動は刃状転位の応力場を緩和するように Cu がける原子と CG 粒子の初期配置である．インデンター先端を転位の内側から外側方向に向けて拡散したためと理解でき球で模擬し，押込まれる部分の近傍でのみ炭素原子を明示的る．したがって，この析出過程は転位および析出相の弾性応に扱い，その周囲は CG 粒子で表現している．図 5 ( b )は，力場に支配されていることがわかる．析出の駆動力の観点か

直径 50 nm 程度のナノドラムのシミュレーションを行ったらは，転位の応力場が Cu 析出の駆動力をより大きくしてい際の押込み荷重と押込み深さの関係を，全系を原子のみ(Full MD)，全系を CG 粒子のみで扱った場合(Full CGP)，

図5
図4

アルミナ転位芯のハイブリッド QM / MD 計算
( a ) 安定構造(緑球 Al 原子，赤球 O 原子)，(b) 有効電荷変位分布．

ま て り あ 第48巻 第 6 号(2009)Materia Japan

ハイブリッド MD/CGP 法によるグラフェン・ナノドラムのインデンテーションシミュレーション(a) 初期配置図(b) 押込み荷重と押込み深さの関係．



び不純物偏析過程の解析，ならびに異なるスケール間の接続・マルチスケール化の試みについて報告した．現状では，各スケールの計算手法が得意とする物質に対する有用性の検証の段階であるが，今後，本特定研究領域内で実験的に創成・観察される共通試料分析への全手法の協調適用を試みる．また，ハイブリッド化に関しては主に QM→MD→CGPまでの接続法開発に重点をおいてきたが，今後は PF 法へのスケールアップ，あるいは PF 法からミクロスケールへのトッ プ ダ ウ ン 手 法 の 開 発 が 必 要 と な る ． こ れ ま で に ， Ginzburg Landau 方程式や Cahn Hilliard 方程式を微視的な運動方程式を基に再構成する試み(15) や，第一原理的計算によるクラスター変分法との PF 法の連成(16) ，あるいは第一原理計算から微視的な界面エネルギー汎関数を近似計算する手法(17) などが報告されているが，精密かつ広範に蓄積された実験データに基づく現 PF 法の高い汎用性を凌駕していな
い．本稿で概説した CGP 法の統計力学的粗視化過程を介した PF 法への接続など，ミクロ・メゾスケールからの熱力学的記述への新たな接続法の開発も，本研究における今後の重要な課題の一つである．
本研究は，文部科学省科学研究費補助金特定領域研究「機図6PF シミュレーションの計算結果( a ) 刃状転位上における Cu ( bcc )析出( b ) Fe 0.4 mass  C 合金の 1023 K 等温時効における g 相の粒界析出過程．

能元素のナノ材料科学(領域番号 474)」(マルチスケール手法によるナノ機能元素材料解析(科学研究費番 19053007 ))により行われた．ここに謝意を表する．

ることになるので， Cu の母相における均一析出に比較して，転位上への不均一析出は相変態のより初期から始まることもわかる．この計算は，転位周辺における溶質原子の拡散挙動の点から見れば，ちょうどコットレル雰囲気の形成過程の計算にも対応する．
次に，FeC 系における粒界析出の解析例について説明する(14) ．計算対象は a 相(多結晶)からの g 相の析出である．
この多結晶 a 相の初期組織形態には，通常の結晶粒成長のPF 計算から導き，初期組織における炭素の濃度場は均一と設 定 し た ． 図 6 ( b ) は 1023 K 等 温 時 効 に お け る 組 織 変 化
で，図中の数値は無次元化された時効時間である．上段が多結晶 a 相および g 相の phase field qi であり，灰色の部分がa 相，白い部分が g 相である(黒は結晶粒界)．下段が炭素濃度場xC である．明暗，(白純 Fe)および(黒Fe1 massC )として，この間を炭素濃度に合わせて連続的にグレイスケールで結んでいる．まず相分解初期(図 6(b)の a))において，結晶粒界に炭素が濃化するとともに，濃化の顕著な領域において g 相が形成される．時効の進行に伴い g 相の体積分率が増加する(同図の b )～ c))．後期では g 相のオストワルド成長段階へと移行していく(同図 c )～ d ))． g 相の全体的な配置を眺めると明確な方向性はなく，a 相内で空間的にほぼ均一に g 相は形成・成長していることがわかる．

文( 1 ) 例えば，J. P. Buban, K. Matsunaga, J. Chen, N. Shibata, W. Y.
215.
Ching, T. Yamamoto, Y. Ikuhara: Science, 311(2006), 212
( 2 ) K. Tsuruta, C. Totsuji, H. Totsuji, S. Ogata: Trans. MRSJ,29(2004), 3669
3672.
( 3 ) S. Ogata: Phys. Rev. B, 72(2005), 045348.
( 4 ) P. Vashishta, R. K. Kalia, A. Nakano, J. P. Rino: J. Appl. Phys.,103(2008), 083504.
( 5 ) R. E. Rudd and J. Q. Broughton: Phys. Rev. B, 58(1998), R5893
R5896.
( 6 ) R. Kobayashi, T. Nakamura and S. Ogata: Mater. Trans.,
49(2008), 2541
2549.
( 7 ) Y. Inoue, J. Tanaka, R. Kobayashi, S. Ogata and T. Gotoh: Mater.
Trans., 49(2008), 2550
2558.
( 8 ) 例えば，S. Tang, T. Y. Hou, W. K. Liu: Int. J. Numer. Meth. Eng.,
65(2008) 16881713.
( 9 ) 小山敏幸まてりあ， 42 ( 2003 ) , 397, 470; ふぇらむ， 9 ( 2004 ) ,240, 301, 376, 497, 905; ibid 11(2006), 647; まぐね，3(2008),
564.
(10) T. Koyama: Chapter 21 in Springer Handbook of Materials Measurement Methods, H. Czichos, T. Saito and L. Smith (Eds), Springer
Verlag, (2007), 1031.
(11) A. Nakamura, T. Yamamoto, Y. Ikuhara: Acta Mater., 50(2002),101
108.
(12) N. Shibata, M. F. Chisholm, A. Nakamura, S. J. Pennycook, T.
85.
Yamamoto and Y. Ikuhara: Science, 316(2007), 82(13) K. Tsuruta, T. Koyama and S. Ogata: Mater. Trans., 50(2009),1018.
1015(14) T. Koyama and H. Onodera: ISIJ International, 46(2006), 1277.
(15) Q. Bronchart, Y. Le Bouar and A. Finel: Phys. Rev. Lett.,
100(2008), 015702.
(16) T. Mohri, M. Ohno and Y. Chen: J. Phase Equilibria and Diffusion,27(2006), 47
53.
(17) S. Yamakawa, K. OkazakiMaeda, M. Kohyama, S. Hyodo: J.
Phys: Conf., 100(2008), 072042.

★★★★★★★★★★★★★★★★★★★★★

このように転位・粒界などの格子不整合領域と相変態を組

1994年 東京大学大学院理学系研究科博士課程修了
1994
1998年 ルイジアナ州立大学博士研究員
1998年 講師，助教授を経て2007年 4 月より現職
専門分野計算材料科学
ナノ構造半導体・セラミックスの大規模シミュレーシ
ョンに従事．

み合わせる考え方は，材料設計に大きな自由度およびプロセスの簡素化を与える新しい方法論として重要である．

.

献課題と今後の計画★★★★★★★★★★★★★★★★★★★★★

多スケールにわたる計算科学手法を用いた欠陥・界面およ


鶴田健二

プロジェクト報告