feram コードによる強誘電体薄膜キャパシタのヒステリシス・ループの超高速分子動力学シミュレーション
岩本昌也, 西松毅, 川添良幸
東北大学

金属材料研究所

計算材料学研究部門

iwamoto@imr.edu
概要: スーパーコンピュータの活用により, 強誘電体キャパシタの分子動力学シミュレーションが実際的なシステムサイズ（〜 100 nm）と時間（〜 1 ns）で可能になってきた. feram コード(http://loto.sourceforge.net/feram/) は我々が独自に開発を進めている強誘電体のための分子動力学計算プログラムである. 第一原理計算から得られた有効ハミルトニアンに基づいて, ペロブスカイト型強誘電体のバルクやキャパシタの高速なシミュレーションが可能である. 今回, feram コー
ドで用いられている物理と並列化等の高速化手法とを紹介するとともに, 計算結果の例として薄膜キャパシタの相転移温度 TC と coercive ﬁeld Ec の膜厚依存性（ナノサイズ効果）やエピタキシャル成長に伴う圧縮応力依存性について報告する. さらに, 次世代スーパーコンンピュータを使うことで計算可能になる強誘電体の物理現象について述べる.

1

はじめに強誘電体薄膜は積層セラミックコンデンサや不揮発性ランダムアクセスメモリ (FeRAM) [1], ナノアクチュエー
タなどに幅広く応用され, 現在, それらの微細化（ダウンサイジング）が強く求められているところである. 分
極特性やヒステリシスの膜厚依存性（ナノサイズ効果）が強誘電体の基礎物性として実験と理論の両面から盛んに研究されているところである. 近年, 分子線エピタキ
シ (MBE) やパルスレーザーデポジション (PLD) などの低温非平衡膜成長技術により品質の高い酸化物強誘電体薄膜の作製が可能になりつつあり, より正確で詳細な物性が明かになりつつある [2].
一方, 強誘電体キャパシタの特性は強誘電体と電極との間の界面の状態に大きく影響を受けるということがよく知られている. 例えば, 強誘電体キャパシタの疲労は,電極近くの不活性層が原因となっている [3, 4]. また, 電極による遮蔽が不完全な場合には, 強誘電体薄膜中に反分極場が部分的に残ってしまい, 膜厚が薄くなるにした
がって coercive ﬁeld Ec が小くなるであろうということが予想されている [5].
1 しかしながら, 強誘電体キャパシタのヒステリシス・ループ, 分極反転, ドメイン壁の動力学的研究は, それらの温度依存性やナノサイズ効果の理解もふくめて困難であった. その理由は, 実験的には「その場観察」が難しく, 理論的には長距離の双極子−双極子相互作用が大規模・長時間の分子動力学シミュレーションを不可能にしていたためである. 今までは, また, 強誘電体−電極間界面構造によって引き起こされる表面効果や反分極場をどのようにシミュレーションに取り込むかが明確ではなかった.

2

計算手法と高速化手法

上述の問題を解決するために, 我々は分子動力学シミュ
レーションコード feram の開発を進めている. feram
コードは, 周期的境界条件 (Fig. 1) の下, ペロプスカイ
ト型強誘電体 ABO3 の第一原理計算で求められた有効ハミルトニアン Heﬀ [6, 7] に基づいた分子動力学計算を行う. 長距離の双極子−双極子相互作用と非一様歪みは,その k-局所性 (k-locality) を用いて逆空間で計算され,その他の r-局所的な項は実空間で計算される. 実空間と逆空間との間の変換は高速フーリエ変換 (FFT) を使うため, 今までは N 2 に比例（O(N 2 ), ここで N はスーパーセル中のユニットセルの数）していた双極子−双極子相互作用の計算時間が大幅に短縮され O(N log N )となり, 強誘電体薄膜キャパシタの実際的な規模（システムサイズ, 〜 100 nm）と時間（〜 1 ns）でのシミュレーションが可能となった. その他の計算手法の詳細は(a)ind = PzEd(d)EEPEE(b)free =(c)P

indlEdlEdEaEFigure 1: 厚さ l ユニットセル（ここでは l = 2）の強誘電体薄膜の模式図. (a) 真空中に置かれた強誘電体薄膜.
(b) 短絡した「遮蔽の完全な」電極間に挟まれた強誘電体薄膜. 強誘電体薄膜のシミュレーションのための「２倍」周期境界条件の (c) 短絡した遮蔽の完全な電極間に挟まれている場合と (d) 短絡した遮蔽の不完全な電極間に挟まれている場合. “E” とマークされた太い水平な線
は, 電極をシミュレートするための静電的な「鏡」である. 強誘電体の表面と電極（静電的な「鏡」）との距離は d/2 ((c) d = 0, (d) d = 1) である. 細い矢印は, ペロブスカイト型結晶のユニットセル毎の局所双極子を表す. 太い破線はシミュレーション用のスーパーセルの周期を表している.
文献 [8] に譲る. 現在まで feram コードは強誘電体薄膜キャパシタの分極特性の温度特性の研究 [9] や疲労によりダメージを受けたキャパシタのヒステリシス・ループの再現 [8] などに成功している.
プログラムの高速化手法としては FFT に FFTW ラ
イブラリ (http://www.FFTW.org/) やスーパーコンピ
ュータ用に最適化されたライブラリを用いている. OepnMP (http://OpenMP.org/) によって並列化されている.

Table 1: 誘電率の周波数依存性を計算する場合に必要とされるであろう計算時間. 32×32×32 ユニットセルのスーパーセルを用い、∆t=2fs の場合.
周波数f1THz
1GHz
1MHz

3

周期
T
1ps
1ns1µsMD ステップ
500
500,000
500,000,000

計算結果の例

Fig. 2 に温度 100 〜 475 K での BaTiO3 強誘電体薄膜キャパシタ（膜厚 l = 4 ユニットセル）のヒステリシス・ループの計算結果を示す. 同様の計算を膜厚と圧縮歪みとを変えて行うことで, 強誘電体−電極間に不活性層がなくても, 膜厚が薄くなるほど相転移温度 TC が低くなることが判った. また, 圧縮歪みの存在は coercive ﬁeld
Ec を増大させることがわかった.

60Pz [µC / cm2]

4020

SR11000
1node(16 cores)
8.4 sec
140 min97 days

動力学シミュレーションコード feram の紹介とそれによる BaTiO3 強誘電体薄膜キャパシタのヒステリシス・ループのシミュレーション計算を示した. feram により分極特性やヒステリシス・ループなどの温度, 膜厚, 圧縮応力, 強誘電体−電極間界面構造依存性がシミュレートできるようになった.
これを基に今後, (1)PbTiO3 など BaTiO3 以外の様々なペロプスカイト型強誘電体の有効ハミルトニアンを用意しそれらの分子動力学シミュレーションを可能にする. (2) 次世代スーパーコンンピュータを用いることでさらに大規模, 長時間の計算を可能にして, たとえばリラクサ（ある特殊な強誘電体）の誘電率の周波数依存性の原理の理解を進める. (Table1 に予想される計算時間を示す. さらに大規模, 長時間の計算を行うにはより一層の高速化・並列化を必要とされるだろう.) (3) 非鉛(Pb-free) の圧電材料の探索など, 新しい強誘電体の物質設計への応用. などに繋げてゆきたいと考えている.

References

0[1] J. F. Scott: Ferroelectric Memories (SpringerVerlag, Berlin Heidelberg, 2000).

-20[2] K. J. Choi, M. Biegalski, Y. L. Li, A. Sharan, J.
Schubert, R. Uecker, P. Reiche, Y. B. Chen, X. Q.
Pan, V. Gopalan, L. Q. Chen, D. G. Schlom andC. B. Eom: Science 306 (2004) 1005.

-40-60-5000-2500

T = 100 K
T = 175 K
T = 250 K
T = 325 K02500
E [kV/cm]

5000

T = 400 K
T = 475 K
T = 550 K

Figure 2: 膜厚 l = 4 ユニットセルの BaTiO3 強誘電体
薄膜キャパシタのヒステリシス・ループの温度依存性のシミュレーション. 面内の 1% の圧縮歪みを導入することでエピタキシャル成長に伴う圧縮応力下を仮定している. TC ≈ 450 K で常誘電相（ヒステリシスなし）から強誘電相（ヒステリシスあり）に相転移している.

4AMD64
1.8GHz dual core
54 sec
900 min
620 days

まとめと今後の課題

第一原理計算から求められた有効ハミルトニアンに基づいたペロプスカイト型強誘電体 ABO3 のための分子[3] V. V. Lemanov and V. K. Yarmarkin: Phys. Solid
State 38 (1996) 1363 [Fiz. Tverd. Tela 38, 2482
(1996)].
[4] H. Z. Jin and J. Zhu: J. Appl. Phys. 92 (2002)
4594.
[5] M. Dawber, P. Chandra, P. B. Littlewood and J. F.
Scott: J. Phys.-Condes. Matter 15 (2003) L393.
[6] R. D. King-Smith and D. Vanderbilt: Phys. Rev.
B 49 (1994) 5828.
[7] W. Zhong, D. Vanderbilt and K. M. Rabe: Phys.
Rev. Lett. 73 (1994) 1861.
[8] T. Nishimatsu, U. V. Waghmare, Y. Kawazoe andD. Vanderbilt, arXiv:0804.1853, to be published inPRB.
[9] J. Paul, T. Nishimatsu, Y. Kawazoe and U. V.
Waghmare: Phys. Rev. Lett. 99 (2007) 077601.