コンテンツにスキップ

分子動力学シミュレーションを用いた水素添加時における鉄の劣化解析#

分子動力学シミュレーションでは，与えられた原子構造からある関数を使ってポテンシャルエネルギーを計算する．この関数は「力場」と呼ばれ，通常，経験的に決められた関数形を使用する．一方，ニューラルネットワーク力場では経験的な知識を使わずに，学習によってゼロから力場を構成する． 本事例では，水素添加時に鉄が劣化する現象の分子動力学シミュレーションの解析事例について紹介する．金属に適する力場としてはEAM（Embedded Atom Method，原子挿入法）力場が知られているが，今回は既存の力場を使わずに，Advance/NeuralMDを使ってニューラルネットワーク力場を作り，水素添加量が異なる鉄にひずみを加えたときの応力の応答から，剛性の指標として用いられるヤング率を調査した．

教師データの作成とニューラルネットワークの学習#

ニューラルネットワーク力場では，系の中の1つひとつの原子について，「近傍の原子構造」から「その原子のエネルギー，およびその原子に働く力」を求める．そのためニューラルネットワークの学習に使う教師データとして，本番の計算で現れ得る「近傍の原子構造」を万遍なく用意する必要がある．

本解析事例では，まずMaterials Projectより取得した体心立方構造のユニットセル(mp-13)1を3×3×3で拡張したスーパーセルモデルを作成した．鉄原子数は54個であり，その内部に3個の水素原子を添加することで水素添加された鉄をモデル化した．Advance/NanoLaboでは，ユニットセルを決定した後に各軸方向の繰り返し回数および添加原子種および数を指定するだけでこのようなモデルを容易に作成することができる．力場作成に使用したモデルを以下の図に示した． training data

本モデルを使用して，自己学習ハイブリッドモンテカルロ（SLHMC）法によってニューラルネットワーク力場を作成した．SLHMC計算の詳細を以下に示す．

  • QE計算はNVTアンサンブルで行った．温度は実際にMD計算にてヤング率の解析を行う300 Kに設定した.
  • SCF計算における交換相関汎関数はPseudopotentialを採用し，鉄にはfe_pbe_v1.5.uspp.F，水素にはh_pbe_v1.4.uspp.Fを使用した．エネルギーのカットオフは40.0 Ry，電荷のカットオフは225.0 Ryとした．
  • NNPにおける対称関数にはChebyshev関数を用いた．カットオフ関数にはtanhを用いており，半径関数の数は50，角度関数の数は30に設定した．またカットオフ半径はいずれも6.0 Åとした．Neural Networkは40ノード×2層とし，活性関数はtwisted tanhを用いた．
  • 初期力場は200ステップのFPMDを行って作成し，NNP-MDはSCF計算毎に2.5 fs～ 500.0 fsで行った．教師データが100個増える度に教師データおよびNNPの更新を行った．教師データには採択されたものと棄却されたものの両方を含み，総数は2600個になった．

MD計算モデルとスキームの設定#

MD計算モデルは体心立方構造のユニットセル(mp-13)1を10×10×10で拡張したスーパーセルモデルに水素を添加することで作成した．鉄原子数は2000個であり，その内部に20個の水素を添加した水素濃度1.0%, 40個添加した水素濃度2.0%の場合と未添加時のヤング率を比較することで水素脆化の解析を行った．計算に使用したモデルを以下に示した（左から，水素濃度0.0%, 1.0%, 2.0%）．

calculation model 0 calculation model 1 calculation model 2

計算スキームは3つのステップから構成した．

  1. まず系のエネルギーを緩和するためにエネルギー安定化処理（Optimization）を行った．
  2. その後，300K，1barの条件下で10万ステップ（50 ps）のシミュレーションを行い，系を平衡化した．タイムステップは0.5 fsとし，セルには等方的な変形のみを許した．
  3. 次に，各モデルについてX方向に0.001Å/psのひずみを与えながら，ビリアル応力の計算を50万ステップ（500 ps）行った．タイムステップは1.0 fsとした．本シミュレーションでは一軸応力の測定を目的としているため，ひずみを与えている方向と垂直な方向の圧力は0barで一定となるように設定を行った．

Advance/NanoLaboではOption画面から，任意の方向へのセルの変形を容易に設定することが可能である．また，User's画面からは任意の物理量の定義，計算及びグラフによる可視化の設定を行うことができる．本シミュレーションではUser's画面からビリアル応力の計算設定を行った．以下に計算スキームの設定画面およびOption画面を示す．

scheme option

Tips 計算を実行する直前の画面でinput fileを直接編集することができる．スキームの設定部分にfix ensemble all npt temp 300.0 300.0 0.1 y 0 0 1.0 z 0 0 1.0と入力することで，y,z方向の圧力を0barで一定に保つように設定することができる．

計算結果#

Advance/NanoLaboでは，計算結果のGUI上での可視化や，CSVファイルの出力を容易に行うことができる．ヤング率はひずみと応力の比として定義されるため，ひずみに対する応力の応答を一次関数でフィッティングすることでヤング率を算出した．応力―ひずみ曲線および，その結果得られた各水素濃度におけるヤング率を表に示した． 我々の解析結果は水素濃度0.0%の純鉄のヤング率は218.98 Gpaであり，実験結果201-216 Gpa2の範囲におおよそ一致していることから，本解析は妥当であると考えられる．水素濃度1.0%のヤング率は208.51 Gpa，2.0%においては195.42 Gpaとなっており，水素濃度の増加によってヤング率が低下していることから，水素脆化現象を適切に再現することが出来た． 本手法を応用することで，鉄に限らず様々な物質の水素脆化現象について網羅的に解析が可能であることが示された． result

水素濃度 (%) ヤング率 (GPa) (計算値) ヤング率 (GPa) (実験値)
0.0 218.98 201-216
1.0 208.51 -
2.0 195.42 -

関連ページ#

  1. mp-13: Fe (Cubic, Im-3m, 229). Materials Project. https://next-gen.materialsproject.org/materials/mp-13 (accessed 2024-03-06). ↩

  2. Takaki, S.; Masumura, T.; Tsuchiyama, T. Young’s Modulus of Single Crystalline Iron and Elastic Stiffness. Tetsu To Hagane 2020, 106 (9), 679–682. ↩