鉄の相変態および界面諸物性に対する分子論的アプローチ
澁田靖

時間スケールが大きくなる．本稿では，長時間，大規模な計

.

は

じ

め

に

算が可能な古典分子動力学法を用いた計算例を紹介し，その長所と短所について議論したい．

材料プロセスにおける組織形成過程の理解のために，原

古典分子動力学法では，得られる計算精度は原子間相互作

子・分子レベルでの粒界や界面諸現象に関する知識の蓄積が用を表すポテンシャル関数に依存する．一般的にポテンシャ

不可欠である．長年にわたる実験的アプローチに加え，近ル関数は，二体(ペア)ポテンシャルと多体ポテンシャルに分

年，分子レベルでの数値解析手法が開発されつつある．ここ類される．二体ポテンシャル Uij を用いると，以下のようにでは，原子・分子レベルの数値計算手法の一つである古典分

全エネルギー U を二原子間の距離 rij のみで記述でき，取り子動力学法の概要を説明し，組織形成過程に関連した固液相

扱いが容易である．

変化や相変態過程での界面・粒界近傍原子の挙動を取り扱っ

U＝た数値計算例を紹介してみたい．具体的には，固液界面での

1∑∑Uij(rij)
2 i j

(1)

融解・凝固過程，ナノ微粒子での融点のサイズ依存性，対称代表例として， Lennard Jones ポテンシャル(2) や Morse ポ

傾角粒界の粒界エネルギーとその温度依存性，鉄 fccbcc 変

テンシャル(3)が挙げられる．鉄やその他の遷移金属に関して態における界面速度の結晶方位依存性等の諸現象に対する最Morse ポテンシャルのパラメータが提案されているが(4) ，近の研究結果を紹介し，分子動力学法計算が持つ長所・短所

二体ポテンシャルでは，弾性定数に関して Cauchy の関係について考えてみたい．

(C12＝C44)を満たしてしまうことや(遷移金属の実験値は C12≠C44)，表面や欠陥での振舞いを正しく表現出来ないなどの.

分子シミュレーション問題が知られている(5)．そこで，多体効果を記述できるポテ
ンシャルの提案がなされ，遷移金属に関しては，現在では原

原子・分子レベルの数値計算といっても，その計算手法は多岐に渡る(1)．化学反応をより正確の再現するために，電子子挿入法( Embedded Atom Method, EAM )(6) が主流となっている．

状態を厳密に計算するほど，取り扱える原子数および時間ス

EAM では，全エネルギー U は二体項(クーロン相互作用)

ケールが限られるというジレンマがあり，その目的によってVij と電子密度 ri の汎関数からなる埋め込み関数 F の項で記

適切な手法を選ばなくてはいけない．動力学に着目した主な述される．

手法としては，電子の波動関数を時間依存の問題として解く
時間依存密度汎関数法，定常の密度汎関数法に原子の動力学を含めたカー・パリネロ法，第一原理計算の行列要素の一部分を経験的パラメータで与えて簡略化したタイトバインディ

U＝1
F ( r i)∑∑Vij(rij)＋∑
2 i jir＝∑fij(rij)

(2)
(3)

j≠i

ング分子動力学法，原子間相互作用をポテンシャル関数とし電子密度 ri は，原子 i を除く近傍の原子が原子 i の位置に作て与える古典分子動力学法などがある．順に原子間相互作用

る電子密度であり多体的な効果を記述できる．具体的な関数の精度は劣るが，計算負荷は軽くなり，取り扱える原子数と形については，Daw と Baskes によるオリジナルの EAM(6)東京大学助教大学院工学系研究科マテリアル工学専攻(〒1138656 東京都文京区本郷 731)
Numerical Approach to the Phase Transformation of Iron and Related Properties at the Interface by Molecular Simulation; Yasushi
Shibuta(Department of Materials Engineering, The University of Tokyo, Tokyo)
Keywords: molecular dynamics simulation, iron, phase transformation, interface, kinetics
2008年 9 月 2 日受理

ま て り あ 第48巻 第 2 号(2009)Materia Japan

に 始 ま り ， 多 元 系 へ の 拡 張 が 容 易 な generalized GEAM
( GEAM )(7) や ， 角 度 依 存 を 考 慮 し た Modified EAM(MEAM )(8)などが開発されている．より精度が高いポテン
シャル関数になるにつれ，関数系が複雑になり，計算機負荷も大きくなるというジレンマがここでも存在する．またFinnis Sinclair(FS )ポテンシャル(9)は，EAM と同様な全エ
ネルギーの記述を取り，簡便な関数形と，短距離相互作用のみで bcc 金属の物性を良好に再現できることから，頻繁に用いられてきた．本稿で紹介する計算も，簡便ながら bcc
金属の物性を良好に再現できる FS ポテンシャルを採用している．その他，ポテンシャル関数についての詳細は文献(1 )を参考にしていただきたい．
以下の章において，FS ポテンシャルを用いた古典分子動力学法による固液界面での融解・凝固過程，ナノ微粒子での融点のサイズ依存性，対称傾角粒界の粒界エネルギーとその温度依存性，鉄 fccbcc 変態における界面速度の結晶方位依存性を解析した例を紹介する．本稿で紹介する分子動力学法の具体的な計算手法としては，運動方程式の数値積分には蛙
飛び法を用い，時間刻みは 2.0 fs とした．系の温度制御にはBerendsen 法(10) を用い，圧力制御には Andersen 法(11) を用いた．その他の詳細は，以下の章で適宜説明する．
図1

.

FS ポテンシャルの融点の見積りと固液界面の挙動(12)固液界面を有するバルク系における融解・凝固過程．(a )初期条件．(b )融点より低い場合．(c)融点より高い場合．

分子動力学計算で用いられるポテンシャルパラメータは，第一原理計算(もしくは実験値)で求めた格子定数，体積弾性る．本稿では粒界近傍原子の挙動の温度依存性に関しては，率，エネルギー曲線といった力学的特性にフィッティングす

上記で求めた融点で規格化した温度，TN＝T/2400で表現しる形で構築されるのが一般的である．そのため融点等の熱力た．また， FS ポテンシャルの A3 および A4 変態点は明確

学的特性までカバーしていないことが多く，温度の定量的議に定義できないため fccbcc 変態過程についての温度依存性論には注意を要する．そこで始めに，FS ポテンシャルのバなどの議論はしない．これらの問題点は，個々の原子間の相

ルクの融点を以下の方法で算出した(12) ．鉄原子 20500 個を互作用のみをパラメータとする分子論的手法に特有の問題で体心立方格子( bcc )構造に配置した周期境界条件セルを，あり，多数の原子が集まって初めて定義できる物性値(融点1000 K お よ び 4000 K で ， NVT 一 定 の 条 件 で 100 ps 熱 緩や変態点)を陽的なパラメータとする現象論的数値解法と大和し，固相および液相の構造を得た．これを図(a )に示す

きく異なる点である．

ように(100)bcc 面を介して結合させたものを初期構造とし，様々な制 御温度で NPT 一定の条 件で 50 ps の 間熱緩和し.

鉄ナノ微粒子の融解・凝固過程(13)(14)た．なお，周期境界条件の為，固液界面はセルの中心と端の二箇所に存在している．さらに温度制御を止めて( NPH 一

次に，半径 1.5 5.0 nm の鉄ナノ微粒子について，その融定の条件)，950 ps の間緩和計算を行った．初期制御温度が解・凝固過程の振舞いや融点のサイズ依存性について検討し融点より低い場合，凝固が起こり固液界面は液相側に移動すた．具体的な計算方法(14)として，まず原子数 2n3 個(n＝10

る(図 1 ( b ))．その際，潜熱の放出により系の温度が上が25 )の鉄 bcc 結晶を真空下で， 1 ns の間 3000 K で加熱し，る．一方，初期制御温度が融点より高い場合，融解が起こり様々な大きさの液滴を準備する．これを 3000 K から 0 K ま

固液界面は固相側に進行する(図 1 ( c ))．その際，潜熱の吸で 2 K/10 ps の割合で冷却し，液滴からの凝固過程を観察し収により系の温度が下がる．最終的にある一定の温度で収束た ． さ ら に 得 ら れ た 鉄 ナ ノ 微 粒 子 を 2 K / 10 ps の 割 合 でし，これを融点とみなすことができる．FS ポテンシャルを

3000 K まで再加熱し，鉄ナノ微粒子の融解過程を観察し用いた bcc 鉄の場合，いずれの初期温度からでも 2400 ± 10た．図に凝固・融解過程の 31250 原子(n ＝ 25 )からなる鉄K に収束した(12)．よって，FS ポテンシャルの bcc 鉄のバル

ナノ微粒子の断面図を示す．図中，赤色の原子は bcc 配位

ク融点はおよそ 2400 K と見積られ，実験値 1808 K と 30 を持つ原子，白色の原子はそれ以外を示す．凝固過程では，程度のずれがあることが分かった．温度に関して直接的に実液滴の表面近傍で初期核が発生した後，凝固が液滴全体に伝験と比較することは FS ポテンシャルを用いる限り困難であ播していく様子が観察された．一方，融解過程では，鉄ナノ




最近の研究

図2

図3

31250原子からなる鉄ナノ微粒子の凝固および融解過程．赤色の原子は bcc 配位をとる原子，白色の原子はそれ以外を表す．

鉄ナノ微粒子の冷却および加熱過程における結合揺らぎの二乗平均．

図4

鉄ナノ微粒子の融点および核生成温度のサイズ依存性．

微粒子は表面から均一に融解していく様子が分かる．また融解の直前では，固液共存のいわゆる表面融解(surface melt-

加熱過程では障壁がないため過加熱が起こらず平衡融点で融ing)の状態が確認された．

解するという実験的経験と一致している．なお，結合揺らぎ

ナノ微粒子の凝固・融解温度は，それぞれ冷却・加熱過程の二乗平均の他に，比熱からも融点を見積ることができにおける結合揺らぎの二乗平均( Lindemann Index )(15) d をる(14) ．凝固・融解は潜熱を伴う一次の相転移であり，相転調べることで確認できる．結合揺らぎの二乗平均は，移温度において比熱はピークを示す．比熱がピークを示す温

d＝2∑n(n＝1) i＜j2〈r 2ij〉－〈rij〉

〈rij〉

(4)

で表される．〈 〉は時間平均を表す．図に冷却・加熱過程

度は結合揺らぎの二乗平均が不連続になる温度と対応しており，これら二つの指標より，鉄ナノ微粒子の融点を定義できる．

における原子数 2n3 個( n ＝ 10 25 )の鉄ナノ微粒子の結合揺

図に，鉄ナノ微粒子の融点および核生成温度のサイズ依

らぎの二乗平均を示す．結合揺らぎの二乗平均は温度が上昇

存性を示す．融点は微粒子のサイズが小さくなるにつれ単調するにつれて大きくなるが，固相から液相に相転移する際に減少するが，核生成温度については，基本的にはサイズの減

不連続的に変化する．よって図 3 において d が不連続に変少とともに低下するが，過冷却液滴中での初期核発生確率と化する温度から，各サイズの鉄ナノ微粒子の凝固および融解の関連から，温度依存性に多少偏差が生じているのが特徴で温度(融点)が分かる．この結果より，鉄ナノ微粒子の凝固温ある．なお，融点の低下は，半径の逆数に比例しており，

度は融点より 5001000 K も低いことが分かった．これは，Pawlow モデル(16)と合致している．

冷却過程では核生成障壁を考慮するため過冷却が起こるが，

ま て り あ 第48巻 第 2 号(2009)Materia Japanた．また(221)S9, (310)S5, (210)S5 面などに小さなエネル

.

鉄対称傾角粒界の粒界エネルギーと予融解(12)(17)

ギーカスプが確認された．これらの傾向は，0 K 下での静的なエネルギー計算に導出された粒界エネルギー(18)(19)と合致

鉄対称傾角粒界に関して，その粒界エネルギーの傾角依存している．

性についての計算例は数多く存在する(18)(20) が，それらの有限温度下での分子動力学法を用いる利点は，高温下でのほとんどは 0 K で静的に評価したものである．一方，実験界面近傍の原子の挙動を直接観察できることである．図においては，高温下で界面の移動が促進された状態で，複数

は，(111)S3 および(112)S3 面の粒界を含む bcc 結晶を，の界面の熱力学的な釣り合い状態から粒界自由エネルギーを

規格化温度0.416および0.916で 1000 ps 緩和した際の粒界近算出している．静的なエネルギー計算による粒界エネルギー

傍原子の軌跡を示したものである．図中の軌跡は， 500 psにはエントロピー項の効果が含まれておらず，厳密には実験から 1000 ps 間の 10 ps 毎の原子位置をたどったものであ

値と直接対応していない．しかし，数値計算によりエントロ

る．規格化温度0.416 では，どちらの粒界面においても，原ピー項の効果を検討することは，手法，計算機能力の両面に子が平衡位置近傍で振動し，対称傾角粒界の形を保っている

おいて困難であるのが現状である．本稿では，有限温度下が，規格化温度0.916 では，(111)S3 面近傍で構造の乱れがで，対称傾角粒界を有する系を緩和し，粒界を有しない同原観察された．これは一種の予融解(premelting)状態と考えら子数の系との内部エネルギー差から，粒界エネルギーを算出

れ，高温下において，粒界エネルギーが上昇する因子である

した例(12) を紹介する．この方法を用いると，粒界エネルギと考えられる(12) ．予融解は図 5 で計算されたほとんどの粒ーの温度依存性，および高温下における界面近傍原子の挙動等が可能である．しかし，エントロピー項の効果を検討する
ためのツールとしての必要十分性は満たしておらず，さらなる物理的考察が必要な事を明記しておく．
図に，鉄 bcc 構造において，〈 100 〉，〈 110 〉，〈 111 〉，〈112〉方位を回転軸とした対称傾角粒界の粒界エネルギーの傾角依存性を示す．図中，粒界エネルギーはステレオ極点図に投影された面指数上に三次元的に表現している．粒界エネ
ルギーは有限温度下で導出したが，FS ポテンシャルの融点が実験値とずれているため，ここでは 3 章で求めたバルクの融点で規格化した温度，TN＝T/2400で表現している．その他，詳細な計算条件は文献( 12 )をご参考いただきたい．
粒界エネルギーは傾角に大きく依存し，特に( 112 ) S3 面で大きなエネルギーカスプが確認された．一方，同じ S 値を持つ( 111 ) S3 面ではエネルギーカスプは確認されなかっ

図6




図5

a ( bcc )鉄対称傾角粒界の粒界エネルギーの傾角
依存性．

(111) S3 および(112 )S3 面の粒界を含む bcc 結晶を，規格化温度0.416および0.916 で緩和した際の粒界近傍原子の軌跡．

最近の研究

界面で観察されたが，大きなエネルギーカスプを示した( 112 ) S3 面でのみ観察されなかった．また面心立方格子(fcc)構造についても同様に，〈100〉, 〈110〉, 〈111〉, 〈112〉方位を回転軸とした対称傾角粒界の粒界エネルギーを求めたが，大きなエネルギーカスプを示した( 111 ) S3 面でのみ，
高温下で予融解が確認されなかった(12)．このような，融点近傍での粒界近傍原子の振る舞いの傾角依存性は，実験で優先的に観察される方位や，複数の界面の釣り合い状態での方位関係などと密接に関連しており，興味深いトピックの一つ
である．

.

鉄 fccbcc 変態における界面のキネティクス(21)鉄 fccbcc 変態過程における界面のキネティクスやその結晶方位依存性については，長年にわたる実験的研究およびそれを裏付ける理論により，多くの知識が蓄積されてきた(22)(23)．そのような背景のもと，分子スケールの数値解析

図7(a )西山Wasserman 方位関係，(b )Kurdjumov 
Sachs 方位関係にある fccbcc 異相界面を有する系の熱緩和過程．

図8

西 山 Wasserman ( N W ) お よ び Kurdjumov 
Sachs (K S )方位関係における整合領域の分布の違い．

が新たに貢献できる方法を考えてみたい．これまでに何度か述べてきたように分子動力学法の利点は，界面や粒界近傍の原子の挙動を直接観察できる点であり， in situ で見ることが困難な実験に対し有効な方法である．固相変態過程に関し，実際に界面のキネティクスを分子スケールで議論することは実験的には未だ困難である．そこで，ここでは様々な方
位の鉄 fcc bcc 界面を有する系を熱緩和し， fcc bcc 異相界面近傍の原子の挙動を直接観察することで，界面のキネティ
クスを検討してみた．

bcc と fcc の異相界面としては，(110)bcc//(111)fcc の関係を満たす界面のうち，[001]bcc//[11̃0]fcc を満たす西山Was-

serman ( N W ) 方 位 関 係 ， [ 11̃1 ]bcc // [ 11̃0 ]fcc を 満 た す
KurdjumovSachs(KS)方位関係が良く知られている．KS
方位関係は NW 方位関係から[111]fcc 軸に対し5.26°
回転させた関係にあたる．ここでは， N W, K S 方位関係および，これらの中間である K S 方位関係から[ 111 ]fcc 軸に対だけ回 転させ た方位 関係の 4 つの条 件に対し 2.00 °, 4.00 °し，これらの界面を有する系を作成し，系全体を規格化温度0.5 で熱緩和させた場合の界面の移動速度およびその優先方位について比較した．
図(a )に，鉄47640原子で構成される NW 方位関係の界面を有する系の熱緩和過程を示す．図中，界面の伝播過程を見やすくするために bcc 配位の原子のみを表示している．
茶色の原子がバルクの bcc 配位(考慮している原子の隣もbcc 配位を満たす 場合)を，橙色の原子が bcc 以外の原子( fcc 等)と接する bcc 配位の原子である． N W 方位関係の場合，界面が初期配置の方位関係を保ったまま，平面的にfcc 側に進行する様子が確認された．界面ではレッジとよばれる階段状の乱れはあるものの， N W 方位関係を持った界面は基本的に平面状に伝播していくことが分かる．図 7 (b )に 99910 原子で構成される K S 方位関係の界面を有する系の熱緩和過程を示す． K S 方位関係の場合，最初のうちはN W 関係同様，界面は平面的に fcc 側に伝播するが，だん
だん界面の乱れが大きくなり，途中で特定の方位に向かい針

ま て り あ 第48巻 第 2 号(2009)Materia Japan

状に優先的に成長した．このスケールで針状成長の速度を正ここで紹介した結果は，東京大学の鈴木俊夫先生ご指導の確に算出することは難しいが，針状成長の先端位置から見積下，大学院生の鷹本真弥君(現トヨタ自動車株式会社)，立ったおおよその速度は，平面成長の速度に比べ約 10 倍早く山真司君のご尽力により得られたものであり，紙面を借りてなっている．針状成長は，KS 関係方位，および NW から深く御礼申し上げます．また，本研究に関してご指導やご助回転させた界面にのみ観察され，NW および NW か
4.00°

言をくださった，東京大学の長汐晃輔先生，熊谷知久先生，

回転させた界面では確認されなかった(21)．また，各ら2.00°

原祥太郎先生，ケンブリッジ大学ジェームズ・エリオット先方位関係における平面成長中の界面伝播速度を比較すると生に厚く御礼申し上げます．

N W 方位関係から[ 111 ]fcc 軸に対する回転角度が大きくなるほど，界面伝播速度が遅くなることが分かった(21)．
図に N W および K S 方位関係を満たす fcc bcc 界面の原子配列を示す．NW 方位関係では，界面原子の整合性

本研究の遂行にあたり，科学研究費補助金・若手研究(A )( No. 18686017 )，萌芽研究( No. 20656036 )および文部科学
省「ナノテクノロジー・材料を中心とした融合新興分野研究
開発研究」複層鋼板プロジェクトの支援を受けました．

の高い部分(図中青点線で囲んだ部分)が網目状に， K S 方

文

位関係では整合性の高い部分が縦長に現れ，NW と KS 方

献位関係の間で，整合領域の分布に違いがあることが確認できる．上記の熱緩和過程の計算では，整合性の高い部分から優先的に，bcc 配位から fcc 配位へと変化していることが確認されており(21)，K S では縦長な整合性の高い部分が他の部分に比べ優先的に変態が進行するため，途中で針状成長に変化したと考えられる．一方，NW 方位関係では整合性の高い部分が網目状に分布しているため，これがレッジ構造を含む平面成長の原因になっていることが考えられる．
以上，鉄 fccbcc 変態過程のシミュレーションは，使用したポテンシャルの信頼性等の問題から現段階では定量的な議論をすることは難しいが，実験的に直接観察することが難し
い領域に新しいアイデアを提供してくれることは間違いな
い．現在，固相変態過程を正しく再現できるポテンシャル関
数を，ベインの格子変形パスや各相の熱膨張率等に着目しながら開発中であり，近いうちに定量的な議論が可能になると信じている．

.

お

わ

り

に

以上，駆け足ではあるが，古典分子動力学法を用いた鉄の界面諸現象に対する計算例として，固液界面での融解・凝固過程，ナノ微粒子での融点のサイズ依存性，対称傾角粒界の粒界エネルギーとその温度依存性，鉄 fccbcc 変態における界面速度の結晶方位依存性を紹介した．紙面の都合上，計算手法や議論の細部まで述べることができなかったので，詳細は個々の文献を参考にしていただきたい．
一般に，古典分子動力学法は，電子物性論に基づく第一原
理計算に比べ厳密性が劣る反面，いわゆる現象論的数値解法のように，実験と直接比較できるほど時空間スケールは大き( 1 ) S. Yip (Ed.): Handbook of Materials Modeling, Springer,
1068.
Dordrecht, (2005), 449( 2 ) J. E. LennardJones: Proc. Phys. Soc., 43(1931), 461
482.
( 3 ) P. M. Morse: Phys. Rev., 34(1929), 5764.
( 4 ) L. A. Girifalco and V. G. Weizer: Phys. Rev., 114(1959), 687
690.
( 5 ) 川添良幸，三上益弘，大野かおるコンピュータ・シミュレーションによる物質化学，共立出版，(1996), 65
66.
( 6 ) S. Daw and M. I. Baskes: Phys. Rev. B, 29(1984), 64436453.
( 7 ) X. W. Zhou, H. N. G. Wadley, R. A. Johnson, D. J. Larson, N.
Long, G. D. W. Smith, P. H.
Tabat, A. Cerezo, A. K. Petford
Clifton, R. L. Martens and T. F. Kelly: Acta Mater., 49(2001),
4015.
4005
( 8 ) M. I. Baskes: Phys. Rev. B, 46(1992), 2727
2742.
( 9 ) M. W. Finnis and J. E. Sinclair: Philo. Mag. A, 50(1984), 4555.
(10) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A.
DiNola and J. R. Haak: J. Chem. Phys., 81(1984), 36843690.
(11) H. C. Andersen: J. Chem. Phys., 72(1980), 2384
2393.
(12) Y. Shibuta, S. Takamoto and T. Suzuki: ISIJ Int., 48(2008),
1591.
1582
(13) Y. Shibuta and T. Suzuki: Chem. Phys. Lett., 445(2007), 265
270.
(14) Y. Shibuta and T. Suzuki: J. Chem. Phys., 129(2008), 144102.
(15) F. A. Lindemann: Phys. Z., 11(1910), 609.
(16) P. Pawlow: Stoechiom. Verwandtschaftsl., 65(1909), 1.
(17) Y. Shibuta, S. Takamoto and T. Suzuki: Comp. Mater. Sci.,
doi: 10.1016/j.commatsci.2008.07.013.
(18) D. Wolf: Philo. Mag. A, 62(1990), 447.
(19) 中島英治・竹内宗孝鉄と鋼，86(2000), 357362.
(20) A. P. Sutton and R. W. Balluffi: Interfaces in Crystalline
386.
Materials, Oxford University Press, Oxford, (1995), 384(21) S. Tateyama, Y. Shibuta and T. Suzuki: Scripta Mater.,
59(2008), 971974.
( 22 ) 谷野 満，鈴木 茂鉄鋼材料の科学 鉄に凝縮されたテクノロジー，内田老鶴圃，(2001).
( 23 ) 榎本正人金属の相変態 材料組織の科学入門，内田老鶴
圃，(2000).
★★★★★★★★★★★★★★★★★★★★★

くなく，その「帯に短し襷に長し」的な側面を欠点として指2004 年 3 月

摘する声は多い．しかし，私個人的には，これを逆手にとり積極的にその側面を利用することで，各計算手法が学問的に独立に発展を続けることから生じる「ニッチ市場」を新たに開拓できると考えている．ここで紹介した例はそれらの試みの一環であり，固液相変化や固相変態過程における界面や粒界近傍の原子の挙動を数値解として提供することで，キネテ
ィクスに関連した実験や理論の解釈のお手伝いができれば幸
甚である．




澁田

靖

東京大学大学院工学系研究科博士課程修了(機械工学専攻)
2004年 4 月 日本学術振興会特別研究員(PD)2005 年 4 月 東京大学大学院工学系研究科助手(2007
年 4 月より助教)
専門分野分子動力学，計算材料学，相変化を伴う諸現象の数値解析
鉄固相変態からカーボンナノチューブ生成過程まで，
原子分子レベルでの解釈が求められる幅広い現象について，理論( theoretical )でなく数値的( numerical )にアプローチする手法を模索する日々．
URL: http://www.mse.t.u
tokyo.ac.jp/shibuta★★★★★★★★★★★★★★★★★★★★★

最近の研究