引っ張り変形を行う
はじめに
さて、前回の講座では材料のせん断変形をモデリングするためのMDについて学んだ。材料科学においてせん断変形と同等に重要な変形が引っ張り変形である。これは、材料を上下方向に伸ばすように変形することを意味する。
実験的な研究において、引っ張り変形は材料の延性（伸び易さ）や破壊強度を解析するために頻繁に使われる方法である。破壊は実生活においてあまり望ましい現象ではないが、落とした皿が割れることなど身近な現象でもある。より深刻なのは、地震、台風、津波などの影響で橋や道路など社会インフラにおいて起きる場合である。
このように、破壊は巨視的なすなわち肉眼で見える現象として観察されるが、その物理的本質は原子間結合の解離であり、その現象の本質を原子レベルで実験的に直接観察することは今の科学ではできない。そのため、いままで壊れなかったものがどうしてその時になって初めて壊れたのか、またどうして他の場所でなくその場所で破壊が起きたのかなど、破壊について理解できない単純な疑問はいくつもある。よって破壊のメカニズム研究は、原子レベルの挙動をモデリングするMDの応用が強く望まれる分野の一つである。
すなわち、ここで学ぶMDは実学的に重要なテクニックの基礎となっている。
引っ張り変形MDの実行
引っ張り変形のMDを行う場合、破壊の起点となる切り欠きなどをあらかじめ入れてMDを行うことが多いが、ここでは前回と同じようにこの変形をMDで行う方法のみに注目したいので、鉄の完全結晶を変形する例を示す。
このMDのスクリプトはscriptディレクトリの中のbccFe_tensile.lcmである。実行は以下のコマンドで行う。
$work> lmp_serial -in script/bccFe_tensile.lcm
これによりbccFe_tensile_*.outという形式のファイルが多数出力される。その一つ、例えば最初のファイルbccFe_tensile_0.outをOVITOに読み込ませると自動的に同じ形式のその他のファイルも読み込まれる。OVITO画面下のPlayボタンでアニメーションを実行して、下の動画のように上下に伸びている様子が確認できれば成功である。
Video 1: 引っ張り変形MDのアニメーション
スクリプトの説明
基本結晶を作る
以下で、このMDで使われたスクリプトをステップバイステップで見ていくことにする。
units metal
boundary p p p
atom_style atomic
variable lc equal 2.83
lattice bcc ${lc} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
variable x_box equal 30*${lc}
variable y_box equal 30*${lc}
variable z_box equal 50*${lc}
region box block 0 ${x_box} 0 ${y_box} 0 ${z_box} units box
create_box 1 box
create_atoms 1 box
unitsコマンド、boundaryコマンド、atom_styleコマンドの説明はもう繰り返さないので、わからない受講者は構造緩和を行うの講座へ戻って確認してほしい。
次のvariableコマンドではlcという変数を定義している。後からわかるがこれはここで使われる鉄の原子間ポテンシャルの格子定数である。
次のlatticeコマンドは過去の講座で出てきた。サンプルの初期値を作るときに使う格子系を定義するためのコマンドである。
次の３つのvariableコマンドはx_box、y_box、z_boxの3つの変数をそれぞれ格子定数の30倍として定義する。
次のregionコマンドも既出であるが、ここではboxという名前の直方体領域を定義する。ここで、最後の引数がunits boxになっていることに注目して欲しい。これまではunits boxの代わりにunits latticeが用いられてきており、その場合は長さの単位が格子定数がなっていた。units boxが指定されると、このコマンドで使用された長さの単位がAngstromになる。要するに上の3つのvariableコマンドはunits boxを使うための準備だったわけである。
ボックスを定義するときに長さの単位を前回までの講座のようにlattice（格子定数）としても今回のようにbox（すなわちAngstrom）しても本質的に違いはない。自分の経験ではboxにする方がプログラミングの自由度が増すので、スクリプトが複雑になるにつれboxの方が楽になる傾向がある。
その後の2つのコマンドは既出であるが、create_boxは上で作られたboxと呼ばれる領域から計算に使用するボックスを作っており、最後のcreate_atomsはそのボックスの中に、上のlatticeコマンドで指定した結晶系を作るようにタイプ1の原子を入れるためのコマンドである。
原子間相互作用の定義
次のブロックを示す。
pair_style eam/fs
pair_coeff * * potentials/Fe_mm.eam.fs Fe
neigh_modify every 1 delay 0 check yes
ここでは原子間ポテンシャルの指定および、それを使って原子間力を計算が行われる原子ペアのリスト（Neighborリスト）の更新に関するオプションを指定している。Neighborリストについては前回の講義で詳しく説明しているの。このブロックの詳しい説明は不要と思う。
次に行こう。
構造緩和を行う
variable dump_step equal 200
thermo ${dump_step}
thermo_style custom step pe lx ly lz press
min_style cg
min_modify dmax 0.05
minimize 1.0e-25. 1.0e-12 10000 50000
このブロックでは構造緩和を行っている。まずvariableコマンドでdump_stepという変数を定義した。thermoコマンドではシミュレーションの途中に熱力学変数を計算し出力する頻度を指定する。その次のthermo_styleコマンドではどういう変数を標準出力に出力させるかを指定する。ここでは時間ステップ、全ポテンシャルエネルギー、ブロックのx、y、z方向の長さ、圧力である。
次のmin_styleコマンドは構造緩和のアルゴリズムの指定、ここでは共役勾配法が使われている。
min_modifyコマンドでは構造緩和の速さの上限を指定する。ここでは1イタレーションでそれぞれの原子位置が0.05 Angstrom以上は変化しないとする。構造緩和がうまくいかない場合はより小さい値にし、収束が遅すぎる場合にはより大きい値にすべきである。ディフォルト値が0.1なのでそれ以外の値を指定したい時だけこのコマンドを使う。普段は何もしなくて良いが、変える場合はこのコマンドを使う。
最後のminimizeコマンドで構造緩和が実行される。4つの引数については構造緩和を行うの講座で説明したのでここでは繰り返さない。
これで引っ張り変形をするための初期値ができた。以下では変形を実行していく。
引っ張り変形に必要な変数を定義する
variable lx_current equal lx
variable ly_current equal ly
variable lz_current equal lz
variable lx0 equal ${lx_current}
variable ly0 equal ${ly_current}
variable lz0 equal ${lz_current}
variable dstrain equal 0.001
variable poisson_ratio equal 0.27
variable lz_delta equal ${dstrain}*${lz0}/2
variable lx_delta equal ${lz_delta}*${poisson_ratio}
variable ly_delta equal ${lz_delta}*${poisson_ratio}
lx_current、 ly_current、およびlz_currentはボックスすなわち材料の大きさである。また、lx0、 ly0、およびlz0も同じものとして定義されている。後でわかるが、lx_currentらは引っ張り変形を行っている時の現在値で、lx0らはその初期値、すなわち構造緩和直後の大きさとして使用する。
さて後で詳細を説明するが、この講座のMDでは前講座のせん断変形を行うのように系をアニーリングした後にある部分に強制的な速度を与えて系を変形するような方法を使わない。ここでは絶対零度のまま材料を小さな刻みで無限の速さで変形し変形後に構造緩和をさせる過程を繰り返すような方法を採用する。
次の行で定義されたdstrainは1刻みあたり変形の割合であり、ここではz方向の初期長さlz0の0.1%ずつ引っ張り変形していくことになる。
次のpoisson_ratioは材料のポアソン比である。鉄結晶がz方向に伸ばされれば横方向に縮むことになるがその割合を示すのがポアソン比であり文献値を代入しておく。
後で出てくるが1刻みの微小な変形はボックスの3方向の境界を微小に動かすことによって実現する。lz_delta、lx_delta、ly_deltaは1刻みあたりに動く幅であり、z方向の定義で2で割られているのは上下の境界を両方とも動かすためで、これにより1刻みあたり上下の長さはlz_delta*2だけ伸びることになる。またlx_delta、ly_deltaの定義にはポアソン比が使われる。
引っ張り変形を実行する
それではスクリプトの最後の部分を示そう。
reset_timestep 0
dump 1 all custom ${dump_step} bccFe_tensile_*.out id type xs ys zs
variable final_strain equal 0.04
variable nloop equal ${final_strain}/${dstrain}
variable loop loop ${nloop}
label loop_start
print " "
print " "
print "%%%%% LOOP INCREMENT ${loop} %%%%"
change_box all x delta ${lx_delta} -${lx_delta} y delta ${ly_delta} -${ly_delta} &
z delta -${lz_delta} ${lz_delta} remap units box
minimize 1.0e-25 1.0e-25 1000 10000
run 0
next loop
jump SELF loop_start
最初のreset_timestepはこれまで何度か出てきたが、時刻の初期化である。
dump ... customコマンドもここまで何度か出てきたので詳細は説明しないが、各原子の途中経過の状態をファイルに出力するためのコマンドである。
次の変数final_strainはこの引っ張りシミュレーションによる最終的な伸びを定義する。ここでは、最終的に4%まで伸ばすことにする。
次の変数nloopは最終的な伸びを上で出てきた微小な引っ張りによる伸びで割っているので微小な引っ張りを行う回数になる。
次のvariable ... loopコマンドは初出である。variable ??? loop $$$という構文で繰り返し計算でインクリメントされる変数を定義する。ここで???は変数名でその初期値は1にセットされる。$$$はその繰り返し数または繰り返し過程におけるloop変数の最後の値である。これで以下の繰り返し計算をする準備が整った。
次のlabelコマンドは一番下のjumpコマンドの行き先を指定する。この後の説明でこのlabelコマンドが繰り返し計算においてどのように使われるかが明らかになるので、このまま進もう。
この後インデントされている部分が繰り返される部分である。まず最初の3つのprintコマンドは標準出力へこの繰り返し計算の途中経過を表示するためのものである。
次の2行は&で繋がっていて1つの文である。このchange_boxコマンドはボックスの境界を強制的に移動させるコマンドであり、x delta、y delta、z deltaのキーワードが記述された場合、それぞれボックスの下境界と上境界の移動差分が指定される。ここではxとy方向では下境界が増加し上境界が減少しているので圧縮され、z方向では逆に伸ばされることがわかる。remapはこのコマンドのキーワードで変形の度合いに比例して原子位置を変えることを要求し、その後のunits boxはこのコマンドに使われた長さの単位がAngstromであることを示している。これが上で何度か述べた微小な引っ張り変形の実体である。
minimizeは上の微小な強制的変形後の構造緩和させるコマンドである。
次のrun 0は少し使い方がトリッキーだ。これまでのrunコマンドは与えられたステップ数だけ動的な計算を行うために使われてきたが、ここでは引数が"0"となっていて奇妙だ。これは単に系の熱力学変数を計算させるためだけにこのコマンドが使われている。
next loopではloopという変数を次の値にする。ここではこれまでの値に1が足されることになる。
最後のjumpコマンドはスクリプトのある位置にジャンプさせるためのコマンドであり、最初の引数にSELFが記述された場合にはこのスクリプトの第2引数で指定されたlabelにジャンプする。最初の引数はスクリプトのファイル名でもよく、他のスクリプトが記述された場合には、そのスクリプトに移動する。ここではこの繰り返し計算ループの最初に戻ることになる。ただし、前の>nextコマンドでloopが繰り返しの範囲（ここではnloop）を超えていたらこの文は無視されて、繰り返し計算が終了することになる。このケースではスクリプト全体の最後なのでLAMMPSが終了することになる。
ディスカッション
少し進んだ内容で恐縮だが、この講座を終える前に受講者に知っておいてもらいたいことがある。１つ前の講座ではせん断変形を行うMDを紹介した。そこでは材料のある部分に強制的に速度を与えて変形する方法を用いた。その方法をこの講座で紹介した引っ張り変形に適用することも可能だ。
ただし、速度を与える方法だと、どうしても実験で行う変形よりも変形速度が非現実的に速くなるという問題がある。実験における変形時間は通常秒単位で電子顕微鏡などで観察できる。しかし、MDで追跡できる時間は長くてもnano秒単位である。一つの解決法は今回のように0 Kで微小に変形して緩和させる手続きを繰り返すことである。この場合は微小な変形の後、毎回構造緩和させているので、変形速度が速くて現実的でないという批判からは逃れられる。
ただし、ここで紹介した方法は有限温度では使えないテクニックなので、格子振動が変形に影響を及ぼす場合などには応用できない。実際、材料は温度が高くなるにつれ塑性変形が起き易くなり、破壊が抑制される。よって有限温度で非常にゆっくり変形するMDが理想的だが、最新のスーパーコンピュータでさえそれには限界がある。すなわち、どちらの方法を選んだとしても一長一短があるということだ。
勘違いしないで欲しいが、有限温度の変形MDが全く現実を反映していないということではない。「MD結果は常に正しい」と盲目的に信じるのは危険であり、丁寧な考察が求められるということを忘れてはならない。
目次へ 前はせん断変形を行う 次は照射損傷をモデル化する