化学気相反応における触媒に関する大規模シミュレーション研究
プロジェクト責任者
藤澤著

義和

株式会社本田技術研究所 基礎技術研究センター者

徳根敏生＊ 1、手島正吾＊ 2、古田照実＊ 1、市川政夫＊ 1、中村

＊1

株式会社本田技術研究所 基礎技術研究センター＊2

財団法人 高度情報科学技術研究機構＊3

独立行政法人海洋研究開発機構 計算システム計画・運用部

壽＊ 2、数納

広哉＊ 3

利用施設 ： 独立行政法人海洋研究開発機構 地球シミュレータ
利用期間 ： 平成 19 年 10 月 31 日〜平成 20 年 3 月 31 日

アブストラクト
我が国を支える先端技術・製品には希少金属が多用されている。これら資源の供給不安に備え、我が国はその使用量削減を目指した代替材料開発が求められている。本研究は、触媒材料に着目し、その特性解析と新たな代替触媒材料を開発することを目指す。本年度は、単層カーボンナノチューブの化学気相反応合成における Fe 触媒の挙動を明らかにすることに焦点をあてた。大規模古典分子動力学・第一原理計算を組み合わせたマルチスケールモデリングを用い、複雑な触媒反応過程の追跡とその反応過程における詳細なメカニズム解明、その主要な制御因子の検討を行った。
Fe 触媒上の炭化水素分解解析には第一原理分子動力学コード PWscf（Ref. 3-1）を、カーボンナノチューブの生成・成長過程解析には古典分子動力学コード DL̲POLY（Ref. 3-2）を用いた。
DL̲POLY 対しては、粒子配置に依存した粒子間ポテンシャル関数の切り替え機能を追加し、さらに地球シミュレータ向きにベクトル最適化を行い大規模計算に対応させた。この結果、我々が提唱する
VLS 反応過程モデル（Ref. 2-1,2-2）に示される鉄触媒表面での CH4 分解過程と鉄触媒内部への炭素拡散を確認することができた。

キーワード： 触媒反応、古典分子動力学、第一原理分子動力学、VLS 反応過程モデル

1.

背景及び目的
我が国は先端科学技術を駆使した自動車や電気・電子製品等のハイテク製品を生産し、世界をリードしている。一方、これらの先端技術を応用した製品には多種類の希少金属（レアメタル）が不可欠であるが、これらは主要産出国である中国、モンゴル等への依存度が高く、必ずしも必要量が安定に確保される状態ではない。希少金属の供給が不足すれば、先端技術製品の生産量は不安定化し、さらには経済にも大きな影響を与えることになる。この様な背景から、希少金属を利用している材料に対し、その使用量削減を目的とした代替材料開発が求められている。本研究では代替触媒材料開発に着目し、巨大な系・長時間計算に適応可能な地球シミュレータを用いた大規模計算による既存触媒の特― 77 ―

性解析と新たな代替触媒の開発を目指す。
今年度は構造が単純で工業的にも期待が高い単層カーボンナノチューブの化学気相反応（以下、CVD）合成における Fe 触媒の挙動を明らかにすることを重点課題とし、大規模古典分子動力学・第一原理計算を組み合わせたマルチスケールモデリングを用いた、複雑な触媒反応過程の追跡とその反応過程における詳細なメカニズム解明、その主要な制御因子の検討を行った。

2.

鉄ナノ触媒反応モデル
カーボンナノチューブの CVD 合成は操作が比較的簡単であるが、その制御因子は、導入ガス種、濃度、流速、温度等と多く複雑である。
（一般的に、炭化水素原料ガスの濃度、圧力、基盤の種類・構造、触媒の温度等が、チューブの品質を決定するパラメータと言われている。）我々は、Fe 触媒の融点解
析、メタン同位体を用いた成長の時差解析などの実験データから、カーボンナノチューブの場合にもVLS 成長モデルが適用されることを提唱している（Ref. 2-1, 2-2）。VLS 成長モデルとは、（1）気相の原料が液相の触媒に溶け込み、（2）原料が過飽和状態となり、（3）固相が析出し、（4）ナノワイヤやチューブが成長するというモデルである（Fig.2-1）。このモデルの詳細な振る舞いの解明は、実験に多くの示唆を与え、新規の機能物質の創生の大きなヒントとなりうると考えられる。我々の反応過程モデルを以下に示す。
（1）炭化水素原料ガスは、触媒表面へ吸着し、分解（2）炭素のみが表面に残り、水素はガスとなり放出（3）反応の進行に伴い、炭素は鉄触媒内部に溶解（4）吸着部の非一様性により、密度勾配、温度勾配が生じ、炭素はさらに鉄触媒内部に拡散（5）鉄触媒の融点は低下し、拡散粒子は飽和状態に達し、炭素は表面に析出
（6）触媒の表面との複雑な相互作用の影響により、析出した炭素は結晶化、又は独自の構造体を形成

Fig.2-1

3.

カーボンナノチューブの VLS 成長モデル

計算手法と地球シミュレータへの最適化
前述した反応過程モデルは触媒が炭化水素を分解する反応過程［（1）-（2）］と炭素からチューブを生成する反応過程［（3）-（6）］の二つから成り立っている。前半を第一原理分子動力学コード
PWscf（Ref. 3-1）にて、後半を古典分子動力学コード DL̲POLY（Ref. 3-2）にてシミュレーションを実施した。水素が介在する分解反応は古典分子動力学で正しく表現することが困難であるため、第一原理手法を採用した。PWscf は既に地球シミュレータに最適化済みである（Ref. 3-3）
。以下、DL̲POLY の機能拡張と地球シミュレータへのベクトル最適化について述べる。

― 78 ―

3-1.

目的とする触媒シミュレーションに対応するための DL̲POLY の機能拡張触媒反応では結合の型が大きく変化する化学反応が生じる。DL̲POLY コードはこの様な化学反応を想定しておらず、原子間ポテンシャルが固定されている。このため近接粒子種と距離に応じ、ポテンシャル関数を変更できるようにコードを改良した。尚、ポテンシャル関数の変更時、エネルギーのジャンプが生じるが（Ref. 3-4）、このエネルギー差は化学反応熱（生成熱）として考え、非常に大きなエネルギージャンプに制限を加えた以外は、現時点ではそのままとした。ポテンシャルとその適用条件を Table.3-1 に示した。
Table.3-1.

DL̲POLY に適用したポテンシャルとその適用条件

結合種

ポテンシャル

適用条件

Fe-Fe

Gupta (Ref.3-6)

距離 2.553 Å以内（カットオフ値）

Fe-C（粒子内部）

Johnson (Ref.3-7)（結合力 0.35eV）

炭素周囲（2.3 Å内）の鉄 5 個以上

Fe-C（表面近接）Johnson（結合力 1.50eV）

炭素周囲（2.3 Å内）の鉄 4 個以下、炭素 2 個以下

Fe-C（表面遠方）

Johnson（結合力 0.14eV）

炭素周囲（2.3 Å内）の鉄 4 個以下、炭素 3 個C-C（粒子内部）

LJ (Ref.3-8)炭素周囲（3 Å内）の鉄 5 個以上

C-C（表面）

Tersoﬀ (Ref.3-9)炭素周囲（3 Å内）の鉄 4 個以下

3-2.

ベクトル最適化

オリジナルコードでは、炭素の Tersoﬀ ポテンシャルの 3 体問題計算がボトルネックであることが判明した。
「近接粒子の抽出（2 重）と相互作用計算（3 重）のための 5 重の do ループ構造」とその中の「炭素粒子抽出のための 7 つの IF 文制御構造」という複雑な流れである為に、ベクトル最適化が阻害されていた。炭素約 5 万粒子を用いた場合、この部分の計算性能は 130Mﬂops, ベクトル化率92％、ベクトル長 142 であった。そこで、（1）IF 文による制御部分を演算部分と分離し、（2）5 重の do ループをベクトル化可能な参照関係の無い 1 重の do ループ構造に書き換えてベクトル長を伸ばす最適化を行った。その結果、1200Mﬂops、ベクトル化率 99.52%、ベクトル長 216 が得られた。ベクトル長は並列数に依存しない状態にすることができた。1000 ステップの MD 計算を実施した場合、40 プロセッサーで 417 秒、80 プロセッサーで 242 秒となり、高速性と共に高並列性も確認された。
オリジナルコードの約 10 倍の高速化を実現し、大規模計算に向けた最適化を完了した。

4.

計算結果

4-1.

第一原理分子動力学による鉄触媒の炭化水素分解シミュレーション

カーボンナノチューブ生成初期段階である、Fe 触媒上での炭化水素の分解・拡散過程を詳細に追跡するため、第一原理分子動力学計算を行った。Fe 触媒粒子はγ相（fcc 構造）を想定した、Fe 55原子のクラスターとした。炭化水素は CH4 を用いた。CH4 の分解、拡散過程を追跡するため、Fe クラスター表面の（001）、（111）面上に、CH4 を 4 個配置した（Fig.4-1（a））。第一原理計算コード
は PWscf を用い、計算条件は、Γ点、セルサイズ 12.5 Å× 12.5 Å× 12.5 Å、エネルギーカットオ
フ 470eV、収束条件 1.6 × 10 − 7eV 以下、T ＝ 1200K、時間刻み 0.13fsec、ポテンシャルは Perdew-― 79 ―

Wang 91（Ref. 4-1）とした。
68fsec（540 MD steps）後、一部の CH4 について分解する過程が確認された（Fig.4-1（b））。分解後、C は Fe 表面に残り、H 原子は Fe から離れるものと吸着したものとに分かれた。同時に、Fe クラスターの部分的な構造の崩れが見られた。これらのことから、Fe クラスター構造が崩れても、その表面では CH4 の自発的な分解反応が生じることが予想される。しかしながら、原子状となった C は Feの 3-hole サイトに保持され、内部への拡散は現状、認められていない。今後、十分な長時間計算を行い、この拡散挙動と炭素飽和状態まで追跡する。

Fig. 4-1

鉄触媒表面のメタンの分解（（a）は初期配置、（b）は 68fsec 後。ピンク：鉄、黄色：炭素、水色：水素）4-2.

古典分子動力学による鉄触媒中の炭素拡散シミュレーション

カーボンナノチューブの生成・成長の条件（温度、炭素源供給量等）を明らかにするため、「炭素の鉄触媒内拡散→炭素粒子飽和→炭素粒子の触媒表面露出→チューブ生成」過程について、古典分子動力学よるシミュレーションを試みた。周期境界を持った一辺 50 Å四方のスーパーセル内に、Fe55 原子クラスターの半球面上に C 単独原子 17 個をランダムに配置したものを初期構造モデルとした（Fig.4-2（a））。アンサンブルは NVT、時間刻み 0.5fsec、T=1070K、温度制御アルゴリズムはBerendsen の方法（Ref.4-2）、各原子間のポテンシャルは前述 Table.3-1 のとおりとした。
10000MD steps（5 ピコ秒）までの計算において、Fe クラスターの表面に配置された C 原子は、クラスター内部に拡散していくことが確認された（Fig.4-2（b））。しかしながら、Fe クラスター表面におけるナノチューブの核生成は現状、認められていない。今後、第一原理計算結果を適用したポテンシャル関数型の見直しと、各種環境条件の探索を行い、成長条件を明らかにする。

― 80 ―

Fig.4-2鉄触媒内部への炭素の拡散（（a）は初期配置、（b）は 5 ピコ秒後。灰色：鉄、青：炭素）

5.

まとめ地球シミュレータを用いた大規模計算による既存触媒の特性解析と新たな代替触媒の開発を目指し、その基盤となる古典分子動力学モデルの構築と第一原理計算を組み合わせたマルチスケールモデリングによる材料計算研究を実施した。その結果、単層カーボンナノチューブ合成用 Fe 触媒における、（1）Fe 表面での CH4 分解反応、（2）Fe 内部への炭素拡散挙動を追うことができた。今後はさらに複雑な触媒反応過程を追跡していく。

謝辞
本研究にあたり、多大なご助言およびご指導をいただきました独立行政法人海洋開発機構計算シス

テム計画・運用部の関係各位ならびに信州大学・遠藤守信教授に深謝いたします。

参考文献
Ref.2-1 F. Ding, K. Bolton, and A. Rosen; J. Phys. Chem. B 108（2004）17369-17377
Ref.2-2 A. R. Harutyunyan,E. Mora and T. Tokune, Appl. Phys. Lett., 90（2007）163120
Ref.3-1 URL: http://www.pwscf.org
Ref.3-2 URL: http://www.cse.scitech.ac.uk/ccg/software/DL̲POLY/
Ref.3-3 T.Furuta, Y.Ito, S.Tejima, et al; Annual Report of the Earth Simulator Center （April2005-March 2006）391-395
Ref.3-4 鉄表面で Tersoﬀ ポテンシャルにより安定な結合状態であった炭素粒子対が、熱揺らぎにより鉄内部に入り込むと、その瞬間に LJ ポテンシャルに切り替えが行われる。この際にエネルギーのジャンプが生じる。これにより、粒子の運動エネルギーが急激に変化し系を不安定にさせる。鉄表面から内部に向かい滑らかな関数でポテンシャルを補間するのが望ましいが、その補間関数が構築されていない。ここでは大きすぎるエネルギージャンプで生じる巨大な斥力値に対し、系が破綻しない一定の斥力値で置き換えるという制限を加えている。エネル
ギージャンプにより一時的に系は不安定になるが、次第に緩和され、長時間の計算が可能となっている。
Ref.3-6 Cleri, F. and Rosato, F., Phys. Rev. B. 48（1993）22
Ref.3-7 Johnson, R.A., Phys. Rev. A , 134（1964）1329-1336
Ref.3-8 Lennard-Jones, J. E. Cohesion. Proceedings of the Physical Society , 43（1931）461-482― 81 ―

Ref.3-9 Tersoﬀ, J., Phys. Rev. B , 39（1989）5566
Ref.4-1 J. P. Perdew, et al., Phys. Rev. B , 46（1992）6671
Ref.4-2 Berendsen, H. J. C., et al, J. Chem. Phys., 81（1984）3684-3690― 82 ―