分子動力学法による Fe とCu の融解の指標に関する検討 - 山口大学

23
分子動力学法による:FeとCuの融解の指標に関する検討
藤田和孝＊
う
服部勝己＋
渡邊匡史＊＊
弘中良† 大木順司††
Mahadevan Khanthast Va. sek Vitekst
Study on Melting Index Using Molecular Dynamics Simulation
of Cu and Fe Single Component Systems
by
Kazutaka FuJITA＊， Masash三WATANABE＊＊， Ryou HIRoNAKA†， Jyunji OHGI††，
Katsumi HATToRI‡， Mahadevan KHANTHA『， Vasek VITEK¶
＊ 1)epartment qノ『ルIecゐanical Engineering， σわεNational College qプTecんnology，
Tokiwadai， Ube 755-8555， Japan
＊＊Advαneεd Course studenちUbe Nationat Oo〃ege of Teeんnology，
コP()kiωadai， Ube 755-8555，」αpαη
†A・αη・h・gurum・k・U・aku・y・鋸，
Higasiikej'iri， Osaka＿Sayama 589-0004， Japan
††Fa・吻・f・E・ginee伽9，ひ痂・・吻・ノY・m・9U・ゐi，
コ監)kiωa(lai2557;σわε 755-8611，」αpαη
＋Department・f Gεneral Education， Ube Nαti・nαt College of Technology
¶0・P・伽翻砿M罐短・1・S・i・nce. ・and・E・9慨而・9，σ珈・rs吻・ゾP・nn・yt・ani・，
Phぬdψゐiα， Pennsylvania 19104，乙1:SA
Abstract: Molecular dynamic simulations were carried out to study the melting indexes using
single component atomic systems (Cu and Fe) under the constant pressure condition.
As all the
indexes(the three dimensional view of atomic structures， the values of the enthalpy and the volume，
the slope of the mean square displacement， the distribution of the radial distribution function， and
the rate of Voronoi polyhedrons)， which were studied on meltings， showed clear changes just before
and after meltings， these indexes could be used for the decision of meltings or not.
Key words: Moteeular Dynamies， Simulation， Amorphous， Melting， Glass， Voronoi Polyhedron，
Enthalpy， Mean Square Displacement， Radial Distre'bution Function.
1緒
言
近年の計算機性能の飛躍的な向上に伴い，コンピュー
タシミュレーションを用いた現象の解明はある程度の
成果を収めつつある. 特に，原子・分子レベルのミク
ロな世界においては，実験および解析困難な対象が数
多く存在するため，計算機の中で現実の世界をある程
ついて，平均二乗変位，3次元観察，エンタルピ，原
子体積，動画分布関数，ボロノイ多面体率を対象に検
討した.
2. 分子動力学法の概要
分子動力学法とは
度モデル化した計算モデルを駆使し，現象を解釈説明
する事で理解を深める方法が有効と言える.
コンピュータシミュレーションの一つの手法として
分子動力学法では，系を構成する全ての原子や分子
についてニュートンの運動方程式をたて，これを数値
分子動力学(Molecular Dynamics;MD)法1)が知られて
積分し，時々刻々の原子や分子の位置や速度を求め，
いる. ここではparrinello-Rahmanのアルゴリズム2＞
これを解析して所要の量を計算して，系のマクロな性
質，系の起こす現象の特徴を調べることができるシミュ
を用いた圧カー定のMD法を使い，432個のFe原子
と500個のCu原子の各単元素結晶を融解させる過程
の原子レベルシミュレーションを行い，融解の指標に
レーション方法である.
2. 1 原子の遅動の求め方
(1998年9月24日受理)
＊宇部工業高等専門学校機械工学科
F部工業高等専門学校専攻科学生
t株式会社浅野歯車工作所
††山口大学工学部機械工学科
‡宇部工業高等専門学校一般科
¶ペンシルバニア大学材料科学工学科
＊＊
原子Zに作用する力:Fiに基づき，ニュートンの運動
の第2法則により加速度aiを求め，これを解いて，原
子の運動を求める.
他の原子から力を受けながら運動するN個の原子か
宇部工業高等専門学校研究報告第45号平成11年3月
24
らなる質点系，このような体系が本研究で取り扱う対
象である. 原子は同種でなくてもよいが，式を複雑に
湘
Z一一一1，2r'. N (1)
創
で表される. Fiは原子iのポテンシャルをriで微分する
ことにより計算される. ここでポテンシャルはFinnis-
湘
Sincl・i・(F-S)ポテンシャノレ3)4)を用いた・F-Sポテン
シャルで用いられているパラメータは，実験で得られ
ている弾性率，擬集エネルギーなどと適合するように
､
も研めト
dt 2 一 m'
ノ。＼○研
モン)P
動方程式は
d2ri ］堅
も研
するだけなので，質量mの同種原子とする. 原子Zの
受ける力をFiとし，原子zの位置をriとすると，運
ンρ ＼⊃ ノ)'b
､､
湘卵め幽
m
Cト
L
決定されている.
数値積分法は，時間刻みを十分小さくとると，極め
て高い精度がえられるとされるギアの予測子一修正子
Fig. 1 Schematic illustration of atom movement
in unit cells under periodic boundary condition.
法(4次精度)5)を用いた.
2. 3 Parinello-Ra㎞anのアルゴリズム
2. 2 周期境界条件
熱力学的平衡状態にある均質な物質のシュミレー
ションを実施し，統計熱力学的関係からその性質を求
めるには，原子数100万といえども立体にすれば1辺
100原子にしかすぎず，これをシュミレーションして
も，表面の存在効果が大きすぎるので，バルクの性質
を計算するには十分ではない. そこでこの問題に対し
て，周期境界条件が適用される.
N個の原子が一辺しの立方体の容器内で運動して
いると考える. 空間に固定した座標系に関する各原子
は弾性的にはね返されるので，原子の運動範囲は
O 一く Xi，Yi，Z'i 一くし
ParinelloとRahmanは，基本セルを表わす3方向
のベクトルを変化させて異方性のある物質に対しても
圧カー定のシミュレーションが行える方法を開発して
いる2). すなわち，Fig2に示した基本セルを表わすマ
トリックスll＝・(a，b，c)と，原子の実座標位置rをマ
トリックスhで規格化した内部相対座標Si(ri ＝ hsi)
に対する運動に関して，下式(2)のうグランジェの方
程式をたて，それから，Siとhの運動方程式(3)，(4)
を導いた.
L一廻Σ野嘱Gξ・Σ野眺Φ(T'iゴ)
十四丁ア［h‘h］一∫)V
(2)
である. この場合，原子は他の原子以外に，この壁か
らも力を受ける. このように，弾性壁によって原子の
9・＝一Σゴ≠、mi(錦):il(・r・ゴ)一G-1GSi
(3)
運動を容器内に制限する条件を固体弾性壁の境界条件
と呼ぶ. 周期境界条件では，原子は壁から力を受ける
Wh ＝ (T-P) a
(4)
ことなく，そのまま容器の外に出てしまい，対面の壁
の相対する位置から同じ速度で同種の原子が入ってく
るとする. この状況を2次元平面内の運動について示
ここで，mは原子質量， Gはhth，Φは原子間ポテン
すとFig. 1のようになる. 中央の立方体を，注目する
シャルである. 式(2)の右辺二幅暢Σ野mβ1磁
体系とし，この立方体を基本セルと呼ぶ. また他の前
は原子の運動エネルギーs第二項のΣ野ΣみΦ(吻
後，左右，上下に積み重ねた基本セル以外のセルをイ
は原子のポテンシャルエネルギーを表している. ま
た，第三項の讐丁仙‘坤ま基本セルの運動エネル
メージセルとよぶことにする. Fig. 1に示したように，
基本セル内の原子の運動がすべてのイメージセルで繰
り返されているので，基本セルからとび出した原子は
ギー，第四項の」)Vは基本セルのポテンシャルエネル
相対する壁から入ってくることになる.
圧力，σはV(ht)一1， Wは質量:の次元をもったパラ
ギーを表している. 式(4)のπは内部応力，Pは外部
メータ(基本セルの質量で，内部応力πと設定圧力値
Pがつり合った緩は物理的意味はない)を表している.
分子動力学法によるFeとCuの融解の指標に関する検討(藤田和孝、渡邊匡史、他5名)
この基本セルの運動方程式は，式(3)のG-iG6‘に影
25
2. 6 温度の制御
響している. そして，このπがPより小さくなれば右
辺は負になり，iiも負となって基本セルの体積が小さ
くなって，πが大きくなる. またπがPより大きけれ
ば右辺は正になり，hも正となって基本セルの体積が
大きくなり，πが小さくなる. このようにπとPがつ
り合うように基本セルの大きさを変化させて圧力が一
温度の制御は速度スケーリング法により行った. こ
れは，原子速度Viをiv/17iii7T viに変更し，現在の温度
Tから設定したい温度TRに調節するものである.
統計熱力学から現在の温度Tは
9NkBT＝Σ野卑1(5)
定に制御される. ここでは静水圧下で行った.
また，今回は液体状態において妥当と考えられるど
の方向にも同じ割合で膨張，収縮する一様膨張6)を仮
で求められる. ここでkBはボルツマン定数， mは原
子質量である. ViにVπをかけてスケーリングし
た、v/iZii7ii;viを式(5)のViに代入すると
定した.
号争Σ∬v3＝㍗×iNkBT ・ gNkB TR(6)
となり，スケーリング後の温度はTRとなる.
2. 7 昇温過程
速度スケーリング法による温度制御は，200ステッ
プに一回行った.
FeはOKから700Kずつ上昇させ，4200Kまで
行った. 各昇温は2万ステップ温度制御を行って，そ
の後2万ステップ緩和を行った. そして原子体積が急
a
激に変化していた2800Kとその前の2100Kの間では，
融点を探すため，温度制御の幅を細かくした. 融点に
Fig. 2 lllustration of unit cell and vector
近ずくにつれて急激な変化がより遅く起こったため長
時間(50万ステップ程度)温度制御させた. そして，
representing unit cell.
その後に2万ステップ緩和させた.
2. 4 初期条件
Cuは， OI〈から2000Kまでは10万ステップ，それ
原子の無次元化初速度島は一様乱数によって与えた
以降は2万ステップ温度制御し，その後に2万ステッ
プ緩和を行った. Cuの場合も融点直前では急激な変
(一〇. 5∼十〇. 5). また初期位置はFeでは体心立方格子
点上に，Cuでは面心立方格子点上に置いた.
化が起こっていた2000Kとその前の1500Kの間で，
2. 5 基本セルにおける条件
温度制御の幅を狭くし，融点を探した.
FeもCuも基本セルは， table Iに示す直方体であ
り・各辺はx，y， z座標と一致させ，それぞれ［100］，
［010］，［001］方向を有する.
Table I Model of analysis
Iron
Number of atoms
Size of unit ce11
Copper
N＝432
O
P7. 199、4×
o
o
17. 199、4×17. 1994
N＝・500
o
o
o
P8。075/4× 18. 075/1× 18。075ノ室
Lattice constant
6
Q. 8665/1
R. 0615A
Mass of atom
9・272×10-26kg
1. 055×10｝25kg
Time step
2. 0×10 15sec
Extemal pressure
0. 1MPa
Mass of box
9. 272×10-28
2. 0×10-15sec
0. 1MPa
1. 055×10-27
宇部工業高等専門学校研究報告第45号平成11年3月
26
3.
FeとCuにおける融解の指標の検討
3. 2 3次元観察
今回の研究では平均二乗変位，3次元観察，エンタ
ルピ，原子体積，動径分布関数，ボロノイ多面体解析
物質には，気相，液相，固相の3つの相がある. 固
相の結晶では原子は規則的に並んで運動するが，気相
によりFeとCuにおける融解の指標を検討した.
と液相では原子は不規則な配置をとって運動する. こ
の相の変化を3次元的に原子位置を表示できるソフト
ウェアを用いて観察した.
3. 1 平均二乗変位
Feについて融解直前の原子配置をx方向からみた
t＝0から時刻t間における各原子の変位の二乗和を
図をFig. 4に，融解直後をx方向からみた図をFig・5
に示す. 融解直前では結晶格子点近傍に配置してい
求め，これを全原子数Nで平均した
て，ここでも固体(結晶)であることが確認できる.
ノ＞
融解直後では，原子はランダムな構造を示していて液
2(ri(t)一ri(o))2/N (7)
体状態であることが確認できる.
i＝o
また，Cuでも同じ結果が得られた.
を平均二乗変位と呼ぶ.
平均二乗変位は，原子がどれだけ移動したかを示す
oψ︾書診O
6噸●9②幽ヤ
06轡9電6
り‘ρも壱魯
6畠㌔壁蓼9。ρ'㍉
2686Kを融解直後とした.
やひ・び轡亀ら
Fig. 3に1414K，2650K，2686K，4317Kの平均二乗
変位の値を示す. Feでは2650Kまでは，平均二乗変
位の値がほぼ一定しており2650Kまでは固体だと分
かる. そして2686K以上では，平均二乗変位の値が
増加しているので液体だと分かる. 融解は2650Kと
2686Kの間で生じており，2650Kを融解の直前，
.
固体状態では増加しない.
''。
轟098亀昏
系が液体や気体状態にあれば時間と共に増加するが，
ﾗ論轡簡＆
轟魅も噌ξ6
◎◇義ず。禽豊
。◎蒐8讐轡δ
目安となる. この量は，原子の拡散と関係しており，
Fig. 4 Atomic arrangement of Fe just before
5
melting(at 2650K， atom diameter＝O. 7 A).
1414K
一一一一一一一一・
/
@26501〈
/
//
2686K
@43171＜:
/
/
/
/
/
肪
く(引掛v
一一一一一・
ノ/
/
/
/
/
ダ
/
/
/
1''
ノ
/
，
1！
1！
ノ
ノ1
！ノノ
ノ
ノ
/！〆:'噂一一一一一一曹一一一一一一噂一一一一一一一一一一一一一一一一
。
〃ノ
100
Time steps
200
Fig. 3 Mean square displacement vesus timesteps.
Cuは，18141〈までは平均二乗変位が一定していて
固体だと分かった. そして1841K以上は平均二乗変
位の値が増加しているので液体だと分かる. そこで
18141〈を融解直前，1841Kを融解直後とした.
Fig. 5 Atomic arrangement of Fe just after
む diameter＝0. 7 A).
melting(at 2686K， atom
3. 3 原子体積
ここでいう原子体積とは，全体の体積を原子数で
割ったもので密度の逆数である.
分子動力学法によるFeとCuの融解の指標に関する検討(藤田和孝、渡邊匡史、他5名)
MD計算を温度を変えて行い，基本セルの体積を
基本セル内の原子数で除して平均原子体積(以下原子
体積とする)を求めた. Fig。6に原子体積と温度の関
係を示す. 実線は温度上昇過程を示していて，その他
の線は冷却過程を示している. Fig. 6(a)では， A点と
B点の間で不連続が起こっており，このA点は融解直
前の点と一致しており，B点も融解直後の点と一致し
ていた. よって原子体積でも固体と液体の区別がはっ
きり分かることが示された. Cu(Fig. 6(b))も同じ結果
27
4. 4 エンタルピ
エンタルピは∬＝U十PVで求められる. ここで，
U(内部エネルギー)は基本セル内の全原子の運動エネ
ルギーの和とポテンシャルエネルギーの和を加えたも
のであり，Pは外部静水圧， Vは基本セルの体積であ
る. 評価はHを全原子数で除して原子1個当りの平
均エンタルピを求め行った.
平均エンタルピと温度の関係をFig. 7に示す. ここ
でも原子体積と同様にA点とB点の間で不連続が生
が得られた.
じており，A点は融解直前と一致し， B点も融解直後
と一致していた. これよりエンタルピでも固体と液体
の区別がよく分かることが示された.
Fe(432At，oms at， 0. IMPa)
＋ Uniforrn dilat，ion
e General dilation /'
一一一一一一一
14
@一loi3 K/sec . /！
@一loi4 1〈7sec if
一一一一一
@一loi'' K/sec ，V'
一2. 5
Fe(432Atoms at O. IMPa)
，∠グ
一一一一
一一一一一
﹁
ジ〆
r一一一一一一一一loi3 K/sec . 一
，
?10i5 K/se. c. . ・一！''''
/
5
1000 2000 3000 4000
3
12
黷撃盾堰eK/sec . ‘」
P
A
o General dilation
-/1
P
〃夕'
一一一！t. r' :
一
・//IB
，！汐ク
，ラノ
A
一一. 一一 Uniform dilation
3
8ε︿＼﹀①眺旧邸ヨq国
一一一一
13
§のくこ・。q＜舞三￡
15
A
ノ〆
. 一/一一2」.
竃
Temperature K
，4タ/
一4.
Relation between temperature and
Fig. 6 (a)
volume/atom in Fe.
1000 2000 3000 4000
Temperature K
Relation between temperature and
Fig. 7 (a)
Enthalpy/atom in Fe.
Cu(500 A乱oms at oユMPa)，/'
薯
18
＋Ulliform dilation
/
ノ
/
o
Cll(500Atoms at O. 1MPa) /'
/
P0i''' K/sec
ノ
∠
撃盾奄S 1〈/sec
/
Ce
C
e
＜＜
I
11
/S
/S
暮
善
一一一一
一'一一
4rり
0
10
1
ー∴蹴三下㎞/ノ
fi 6
＋Uniform dirar，ion
鋒ε︿＼﹀・
'
/
'
/
！/
蹄
、〃/
4
8
ノ
！！
碁14
A
o looo 200armt 60
Temperature K
/'
一一一
/2
/ '
P013K/sec ，・'
一『雪一
P014K/sec z'！:
一一一一
P015K/sec，'
！〃''
3
眺3邸調帽国
α斥∫二/
''''i
12
Fig. 6 (b)
。 General dilation
〃γ A
ノノ
. /！
ノク'
読/
Relation between temperature and
一3. ，r)
volume/atom in Cu.
竃
O 1000 2000 3000
Temperature 1〈
ここで不連続が生じる理由は物質が固体から液体に
変わると密度は小さくなると考えられ，その逆数であ
る原子体積は大きくなり，その変化が不連続になって
現れている.
Fig. 7 (b)
Relation between temperature and
Enthalpy/atom in Cu.
なお，このHの不連続の量は，潜熱に相当する7).
宇部工業高等専門学校研究報告第45号平成11年3月
28
3. 5 動径分布関数
n5， n6は，それぞれ三，四，五，六角形の数である.
ボロノイ多面体解析は，局所的な原子構造を与える.
動径分布関数とは，着目する原子から，ある距離r離
れた位置に原子が存在する確率のことである. その求
Fig. 9にFeの融解が起こるまでに多かったボロノ
イ多面体(0608)，(0446)(多い方から2つ)と融解
め方は，MD計算であらかじめ求められた原子の座標
した後に多かったボロノイ多面体(0284)，(034
データからnk(r， t)(時刻tに原子kを中心とした半径
r一△r/2 r十△r/2の2球面ではさまれた球殻中
の原子数)を計算し，全原子数Nで割って平均したも
のを
6)(多い方から2つ)を示した・(0608)と(0446)
は，OKから融解する前までは多く存在しているが，
融解直後には，激減していることが分かる((0608)
は0％・(0446)はほぼ2％)・また，(0284)と(03
64)は融解するまでははぼ2％であるのに融解直後か
〈n(・)〉＝三幅蘭，t)(8)
ら急には5％，6％と多くなっているのが分かる.
k＝1
によって求め，動径分布関数g(r)は
から得られる. ここでρは数密度(≡N/V)であり，
△rはO. 05a(a:格子定数)とした. 動径分布関数は基
本セル中の平均的な原子構造を示す.
Fig. 8は， FeにおけるOK，融解直前および融解直
後の原子位置の動径分布関数である. 完全結晶である
OKと，融解直前を比較すると，ピークを示す位置は
ほとんど変わらない. 融解直後では，第2ピーク以降
のも曾ゐ。食。q。﹄。﹀も￡娼属
g(r) ＝〈 n(r) 〉 /4Tr2pAr (9)
100
x
Fe
N
N
@
@
N
XNyX
@
xR
@
T0
@［9698) :r
一'一
''［9446) il
'一
Z''
一. 一
@［9？8‘) i
@(0364) 1
〆届
器
2000 ' '4000
のピークは平滑になっており，融解直前と融解直後の
区別がはっきり表われている. よって動径分布関数で
も融解の指標になることが示された.
N
x
Temperature K
Fig. 9 Rate of voronoi polyhedron under melting
from crystalline Fe.
OI〈500step平圭勾
2s s
Cuのボロノイ多面体率をFig. 10に示す. 融解する
前に多く存在していた(0364)と(0446)と(0447)と
(0365)を示した. この4つの多面体は融解する前は多
5
﹂Q属
2660K」「OOstep平均
く存在しているが融解してからは急激に減っている.
)4)
))
ゐρ0 7︻U
T0
噂魅嚇
黄，囚
(
﹁＼
ら(n3，n4，n5，n6)多面体と表示する. ここで， n3，n4，
ロ
鷲
の
ロ
^︻
定義される. 個々の多面体を構成する面の形状と数か
鴫ー
@
o
ボロノイ多面体は，ある特定の原子と他の任意の原
子とを結ぶ線分の二等分面が作る最も内側の多面体で
コ
凶
3. 6 ボ日ノイ多面体解析
一
@
≦
Radial distribution function.
@
Cu
4ρO
Fig. 8
@
nδ4轟﹂4り0
。
1一 2 一'3
Dist，ance(a)
三
@
ρ04
餌
、∼ノ！'''''一一'一一. 一＿一一一＿
㌫d七①毫で臼89￡申。＄邸国
5
﹂O.
26801〈500g. teps平均
0(0
(
((00
100
1000 2000 V
Temperature K
R00b
Fig. 10 Rate of voronoi polyhedron under melting
from crystalline Cu.
分子動力学法によるFeとCuの融解の指標に関する検討(藤田和孝、渡邊匡史、他5名)
このようにボロノイ多面体解析は融解の指標に使
えることが分かる. 特にFeでは体心立方構造を示す
参考文献
(0608)多面体が，Cuでは面心立方構造を示す(0446)
1)上田顕著，''コンピュータシミュレーション''，
多面体が最も有効である.
(1990)，朝倉書店.
2)M. Parinenello and A. Rahman， J. Appl. Phys. ， 52，
7182(1981).
6. 結
言
3)M. W. Finnis and J. E. Sinclair， Philosophical
Magazine A. ， 50，45(1983).
4)G. J. Acland， G. Tichy， V. Vitek， and M. W. Finnis，
Parrinello 一RahmanのアルゴリズムとFinnis 一
Philosophical Magazin A. ， 56，735(1987).
Sinclairポテンシャルを使った分子動力学法を用いて，
5)M. P. Allren and D. J. Tildesley， ''Computer
CuとFeの各単結晶における融解の指標の検討を
Simulation of Liquids''， CLARENDON PRESS，
行った.
その結果，3次元観察，エンタルピ，原子体積，動
OXFORD， 341(1987).
6) S. NoSe and M. L. Klein， Molecular Phys. ， 50， 1055
径分布関数，体心立方と面心立方の結晶構造の特徴を
(1983).
示すボロノイ多面体の割合は，平均二乗変位で確認し
7)D. A. Porter and K. E. Easterling， ''Phase
た融解の直前と直後の温度において急激な変化を生じ
Transformations in Metals and Al｝oys'' ， CHPMAN
ており，いずれも融解の指標とすることができた.
＆ HALL， 8(1980).
29

Download Report