LAMMPSでlatticeコマンドのorientオプションを用いて特定の結晶方位を持つセルを作成し、格子緩和を行った際に、せん断応力成分が適切に緩和されないという問題に直面した。本記事では、その現象と解決策について共有する。
結論から言うと、最初に最小単位のセルで格子緩和を完全に行い、その後replicateコマンドで目的のサイズに拡大することで、この問題を回避できた。初期構造の応力が十分に緩和されていることは、その後の計算の信頼性を担保する上で極めて重要である。ちなみに、わたしはこの確認を怠ったために、多くの計算をやり直すこととなった。
問題の現象
BCC鉄の格子緩和を例に説明する。ここでは、一般的な(100), (010), (001)の結晶方位ではなく、X, Y, Z軸がそれぞれ(111), (-110), (-1-12)方向を向くようなシミュレーションセルを準備する。計算すればわかるがこれらは直行している。
以下に、問題が発生した際のLAMMPSの入力スクリプトを示す。
# initializationunits metal # 単位系定義dimension 3 # 次元boundary p p p # p (periodic) or s or f or matom_style atomic # Define what style of atoms to use in a simulation
# create geometrylattice bcc 2.85 orient x 1 1 1 orient y -1 1 0 orient z -1 -1 2region box block 0 20 0 20 0 20create_box 1 boxlabelmap atom 1 Fecreate_atoms Fe box
# potentialspair_style eam/fs # ポテンシャルタイプ、Finnis Sinclairpair_coeff * * ./Fe_Ackland04_ZBL.eam.fs Fe # ポテンシャルファイル読み込み
# calctimestep 0.001 #
dump 1 all custom 1000 relax_*.dump id type x y z xs ys zs vx vy vz fx fy fz # 座標出力
thermo 500 # output thermodynamics every N timestepsthermo_style custom step temp pe pxx pyy pzz pxy pxz pyz # the style and content for printing thermodynamic data to the screen and log file.
# 原子座標緩和neigh_modify delay 0 every 1 check yes # the building and use of pairwise neighbor lists.min_style cg # cg:共役勾配法min_modify dmax 0.001 # The dmax parameter is how far any atom can move in a single line search in any dimension (x, y, or z).minimize 0.0 1.0e-12 10000 10000 # 緩和 (minimize etol ftol maxiter maxeval)
# 原子座標＋体積緩和fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 1.0e-5 nreset 100 # セルの大きさを変える。min_modify line quadratic # the quadratic line search algorithm tries to reduce the forces to zero, while guaranteeing that the energy changes is not positive (uphill).minimize 0.0 1.0e-6 10000 10000 # 緩和 (minimize etol ftol maxiter maxeval)8行目でorientを用いて結晶方位を、9行目でシミュレーションセルのサイズを一度に定義している。また、24行目のthermo_styleで静水圧成分（Pxx, Pyy, Pzz）だけでなく、せん断応力成分（Pxy, Pxz, Pyz）も出力している点がポイントである。
このスクリプトを実行（lmp_serial -in relax.cmd）すると、以下のようなログが出力される。
Step Temp PotEng Pxx Pyy Pzz Pxy Pxz Pyz 1056 0 -256805.04 9934.46 7092.2933 12886.431 4.0878753e-12 -6326.5485 -3.1600352e-12 1392 0 -256819.73 -1.5004335e-07 -6.7607257e-06 6.55549e-06 -5.7846837e-12 -5449.1859 -2.4248676e-11緩和計算は1392ステップで収束しているが、Pxzに-5449.18 barという非常に大きな応力が残存していることがわかる。これでは、適切に緩和された初期構造とは言えない。
原因は完全には特定できていないが、orientで軸を指定した際に、大きなシミュレーション領域を直接作成すると、原子座標の丸め誤差などが累積し、結果としてせん断応力が残留するのではないかと推測している。
解決策
この問題を解決するために、以下の2段階のプロセスを踏むアプローチを試みた。
regionコマンドでまず最小単位のセル（0 1 0 1 0 1）を作成し格子緩和を行う。- 緩和済みの単位セルを
replicateコマンドで目的のサイズ（20 20 20）に複製する。
このアプローチに基づいた最終的な入力スクリプトは以下。
# initializationunits metal # 単位系定義dimension 3 # 次元boundary p p p # p (periodic) or s or f or matom_style atomic # Define what style of atoms to use in a simulation
# create geometrylattice bcc 2.85 orient x 1 1 1 orient y -1 1 0 orient z -1 -1 2region box block 0 1 0 1 0 1create_box 1 boxlabelmap atom 1 Fecreate_atoms Fe box
# potentialspair_style eam/fs # ポテンシャルタイプ、Finnis Sinclairpair_coeff * * ./Fe_Ackland04_ZBL.eam.fs Fe # ポテンシャルファイル読み込み
# calctimestep 0.001 #
dump 1 all custom 1000 relax_*.dump id type x y z xs ys zs vx vy vz fx fy fz # 座標出力
thermo 500 # output thermodynamics every N timestepsthermo_style custom step temp pe pxx pyy pzz pxy pxz pyz # the style and content for printing thermodynamic data to the screen and log file.
# 原子座標緩和neigh_modify delay 0 every 1 check yes # the building and use of pairwise neighbor lists.min_style cg # cg:共役勾配法min_modify dmax 0.001 # The dmax parameter is how far any atom can move in a single line search in any dimension (x, y, or z).minimize 0.0 1.0e-12 10000 10000 # 緩和 (minimize etol ftol maxiter maxeval)
# 原子座標＋体積緩和fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 1.0e-5 nreset 100 # セルの大きさを変える。min_modify line quadratic # the quadratic line search algorithm tries to reduce the forces to zero, while guaranteeing that the energy changes is not positive (uphill).minimize 0.0 1.0e-6 10000 10000 # 緩和 (minimize etol ftol maxiter maxeval)
# cell 拡大replicate 20 20 20
# 原子座標緩和neigh_modify delay 0 every 1 check yes # the building and use of pairwise neighbor lists.min_style cg # cg:共役勾配法min_modify dmax 0.001 # The dmax parameter is how far any atom can move in a single line search in any dimension (x, y, or z).minimize 0.0 1.0e-12 10000 10000 # 緩和 (minimize etol ftol maxiter maxeval)
# 原子座標＋体積緩和fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 1.0e-5 nreset 100 # セルの大きさを変える。min_modify line quadratic # the quadratic line search algorithm tries to reduce the forces to zero, while guaranteeing that the energy changes is not positive (uphill).minimize 0.0 1.0e-6 10000 10000 # 緩和 (minimize etol ftol maxiter maxeval)このスクリプトによる緩和計算の結果は以下の通り。
Step Temp PotEng Pxx Pyy Pzz Pxy Pxz Pyz 2135 0 -252931.74 -0.0006457296 -0.0034065106 0.0067741329 -3.3083e-11 1.3272091e-09 -1.2592623e-09 2136 0 -252931.74 -0.0006457296 -0.0034065106 0.0067741329 -3.3083e-11 1.3272091e-09 -1.2592623e-09静水圧成分、せん断応力成分ともに、ほぼゼロに近い値まで緩和されていることが確認できる。
まとめ
LAMMPSでorientオプションを使って複雑な結晶方位のセルを扱う場合、単純な緩和計算では意図せず応力が残存するケースがあることがわかった。このような場合でも、まずは最小単位で構造緩和を行い、それを元に系を拡大するという手順を踏むことで、解決することができた。
シミュレーションを行う際は、「緩和計算が収束した」という事実だけでなく、応力テンソルなどの物理量が期待される値になっているかを具体的に確認することが、結果の妥当性を保証する上で不可欠であることを再認識することとなった。
