広島工業大学紀要研究編
第 57 巻（2023）23-25

報

告Tersoff 型原子間ポテンシャルを用いた分子動力学法による引張変形下における単結晶シリコンの不安定挙動
金舛育実＊・大村

訓史＊（令和 4 年10月28日受付）

Molecular dynamics study of instability behavior of single-crystal silicon
under uniaxial tension using Tersoff potential
Ikumi KANEMASU and Satoshi OHMURA（Received Oct. 28, 2022）

Abstract
Instability behavior of single-crystal silicon under tensile deformation has been investigated bymolecular-dynamics simulations based on Tersoff-type interatomic potential for the two proposed
potential parameters. From the stress-strain curves obtained from the simulations, it is confirmed thatthe maximum tensile stress decreases with increasing temperature for both parameters. However,there are significant differences in the behavior after the stress release depending on the parameter
type, indicating that defects such as dislocations or voids are formed differently.
Key Words: molecular dynamics, uniaxial tension, stress-strain curve, silicon, Tersoff-type potential

陥に適しているという特徴がある．結晶シリコンの強度特1 ．はじめに性を原子レベルで理解する上では，このような特徴に加シリコンは安定した構造を持つことや加工が容易なこと

え，引張変形下のような非平衡状態においてポテンシャルから半導体の基板として広く使用されており，その強度特パラメーターの違いが計算結果にどのような影響を与える性に関する研究がさかんに行われている．さらに近年の超かを知る必要がある．そこで，今回 Si（B），Si（C）の 2 種

微細加工技術の進展によって，強度特性を原子レベルで理類のポテンシャルパラメーターを用いた分子動力学シミュ解することが必要とされており，分子動力学法がその強力レーションを行い，引張変形下における単結晶シリコンの1, 2な手法として使われている ．分子動力学法とは，物質を不安定挙動の違いを調べた．

構成する各原子の運動方程式を数値的に解くことによりそ2 ．計算の詳細の物質の性質を調べる手法である．分子動力学法から得られる結果は原子に作用する力を求めるための原子間ポテン

結晶 Si の単位格子を 5 × 5 × 5 並べた総粒子数1000個シャルに依存することが分かっている．結晶シリコンの分の系を用い，x，y，z のいずれの方向にも周期境界条件を子動力学シミュレーションにおいては，Tersoff 型の原子適用することで単結晶のバルクとした．分子動力学計算は

34間ポテンシャルが用いられることが多く，Si（B），Si（C）

大規模計算用の LAMMPS5 を使用した．原子に作用するの 2 種類のパラメーターセットが用意されている．Si（B）

力は Tersoff 型原子間ポテンシャルを用い， 2 種類のポテは表面の構造を精度よく記述でき，Si（C）はバルクや点欠

ンシャルパラメーター Si（B），Si（C）を使用した．具体的＊

広島工業大学工学部環境土木工学科― 23 ―

金舛育実・大村訓史

なパラメーターの値は表 1 に示している．温度は 1 K，いる．ともにグラフの概形は一致しており，応力が解放さ300 K，600 K の 3 通 り と し， そ れ ぞ れ の 温 度 に お い てれた後はひずみを加えても応力の増加はほとんど見られな

Maxwell 分布に従うように各原子の初速度の与え方を変い．

え， 8 セット用意した．
まず，それぞれの温度に対して20 ps の温度・圧力一定の分子動力学シミュレーションを行い，その後，等温・等圧アンサンブルを用いて x，y 方向の応力が 0 になるように緩和させながら，ひずみ速度1.0×109 /s で z 方向［001］に引張荷重を加え続けた．
表1

本研究で使用した Tersoff 型ポテンシャルの2, 3, 4
パラメーター

Si (B)Si (C)

A [eV]

3

3.2647×10

1.8308×103

B [eV]

9.5373×101

4.7118×102

λ 1 [Å-1]

3.23942.4799-1

λ 2 [Å ]1.32581.7322-1

λ 3 [Å ]1.32581.7322β3.3675×10

n-11.1000×10-6

1

7.8734×10c2.2956×10
4.8381d2.04171.6217×101

h0.0000-5.9825×10-1

R [Å]

3.0

2.85

D [Å]0.20.15-1

図 1 （上部）パラメーター Si（B）を用いた計算から得られた応力-ひずみ曲線．赤，青，緑の線はそれぞれ，
1 K，300 K，600 K の 結 果 を 示 す．（ 下 部 ）600 Kの引張過程（a）ε＝0.1，（b）ε＝0.27，（c）ε＝0.4における原子配置51.0039×10

3. 2

Si（B），Si（C）のどちらのパラメーターを使用しても，3 ．結果と考察
3. 1

Si（B），Si（C）の応力-ひずみ曲線の違い

温度が高くなるにつれ，応力開放が起こるまでに与えられ応力-ひずみ関係の温度依存性るひずみは小さくなり，最大応力も低下する．しかし，Siパラメーター Si（B），Si（C）を用いた分子動力学シミュ（C）の方が300 K と600 K の差が大きい．また，600 K ではレーションによって得られた応力-ひずみ曲線の温度依存どちらもε＝0.27付近で応力開放が起こるが，その後の挙性をそれぞれ図 1 ，図 2 に示す．ここで横軸は工学ひずみ動が Si（C）に比べ Si（B）の方が初速度の与え方による影響

ε，縦軸は応力σ（GPa）である．また，600 K における引が大きく，ばらつきが出ている．

張過程の原子配置を図の下部に示す．
図 1 に示すように，パラメーター Si（B）を用いた場合，3. 3

Si（B），Si（C）の引張変形下の原子配置の違い

1 K ではε＝0.4付近で一度応力が減少したものの応力開600 K における原子配置（a）ε＝0.1，（b）ε＝0.27，（c）

放は起こらず再び応力は増加している．300 K においてε＝0.4について考察する．図 1 ，パラメーター Si（B）を使は，ε＝0.3付近で応力開放が起こる．また，600 K におい用した場合の（a）ε＝0.1では，構造はほとんど変化していては，ε＝0.25付近で応力開放が起こっている． 1 K，ない．しかし（b）ε＝0.27においては原子構造が乱れ，数300 K では初速度の与え方による影響は小さく，ほぼグラ原子サイズのボイドが形成された（図 1 －（b）の左上部）．

フが一致している．一方，600 K では応力開放が起こるま（c）ε＝0.4では，ε＝0.27で形成された微小なボイドが成でに与えるひずみにばらつきがあり，解放後も初速度の与長し数十原子サイズのボイドとなっており，上部に集中しえ方による影響が他の温度に比べて大きい．て発生している．一方，図 2 のパラメーター Si（C）を使用

パラメーター Si（C）を用いた場合（図 2 ）， 1 K ではεした場合において，（a）ε＝0.1では，Si（B）の場合と同様，＝0.42付近で応力開放が起こり，応力が12 GPa ほど低下構造はほとんど変化しておらず図 1 －（a）とあまり違いは見し，その後応力は増加している．300 K と600 K においてられない．（b）ε＝0.27においては，原子構造の乱れが生はε＝0.33，ε＝0.26付近でそれぞれ応力開放が起こってじている．（c）ε＝0.4では系全体の原子配置が崩れている― 24 ―

Tersoff 型原子間ポテンシャルを用いた分子動力学法による引張変形下における単結晶シリコンの不安定挙動が，Si（B）のようなボイド（図 1 －（c））は形成されていな

1 K における応力-ひずみ曲線を見ると，図 1 に示すよい．この違いは，応力開放後の挙動にも影響を与えていうに，パラメーター Si（B）を用いた場合，ε＝0.4付近で応る．図 1 の応力-ひずみ曲線を見ると，600 K の応力開放力が 1 GPa ほど減少しているが，再び応力は増加してい後の挙動が初速度の与え方によって大きくばらついているる．これは図 3 －（a）に示すように構造全体で転位が発生しことがわかる．

た後すぐに構造が回復したためだと考えられる．一方，パラメーター Si（C）を用いた場合，このような構造の回復は見られない．しかしながら，応力開放後に応力が再び増加している．これはε＝0.41で構造が変化したあとも，ある程度引張強度を持っているためだと考える．

4 ．まとめ今回，Si（B），Si（C）の 2 種類のポテンシャルパラメー
ターを用いた分子動力学シミュレーションによって引張変形下における単結晶シリコンの不安定挙動について検討し，以下の結果が得られた．
どちらのパラメーターにおいても温度が高くなるにつれ応力開放が起こるまでに与えられるひずみは小さくなり，最大応力も低下する．しかしながらパラメーターの種類によって応力開放後の挙動に違いがあり転位やボイドなどの欠陥の形成状況が異なる． 1 K での引張変形の場合，Si（B）を用いると応力開放は起こらないが Si（C）を用いたシ
図 2 （上部）パラメーター Si（C）を用いた計算から得られた応力-ひずみ曲線．赤，青，緑の線はそれぞれ，
1 K，300 K，600 K の 結 果 を 示 す．（ 下 部 ）600 Kの引張過程（a）ε＝0.1，（b）ε＝0.27，（c）ε＝0.4における原子配置

ミュレーションではε＝0.4付近で構造が乱れ応力開放が起こる．また600 K では Si（B）を用いると初速度の与え方によって応力開放が起こるまでに与えられるひずみと応力開放後の挙動にばらつきが出ている．一方 Si（C）を用いた
場合，初速度の与え方によるばらつきは小さい．

これは図 1 （－ c）のようなボイドが形成された場合は破断に伴って応力が解放されやすいが，初速度の与え方によっ謝てはボイドが形成される場所や大きさが異なるなどして，ε＝0.5までに破断に至らないことがあるためである．一方，パラメーター Si（C）を用いた場合はボイドが形成されないため，応力開放後の挙動は初速度の与え方に依存せ辞本研究は公益財団法人サタケ技術振興団の助成を受けたものです。

文献

ず，Si（B）の場合（図 1 ）に比べグラフの概形にばらつき1がない．

K. Shimamura et al J. Phys.: Conf. Ser. 40２ 012044（2012）

2大石直樹，神戸大学大学院博士課程前期課程機械工学専攻

修士論文（2010）

3J. Tersoff Phys. Rev. B 37 6991 （1988）

4J. Tersoff Phys. Rev. B 38 9902 （1988）

5A. P. Thompson et al Comp Phys Comm, ２7１ 10817（2022）

図3

1 K における引張変形下での原子配置．（a）パラメー
ター Si（B）を用いた場合 （b）パラメーター Si（C）を用いた場合．

― 25 ―