アニーリングを行う
はじめに
ここでは鉄結晶をある有限温度に保ちながら時間発展させていく計算法を学ぶ。受講者はすでに構造緩和を行う方法を学んだはずだ。構造緩和は絶対0度で系のエネルギーを最小にしていくプロセスである。アニーリングも有限温度で構造が緩和される場合もあるが、ここでは便宜上絶対0度の場合を構造緩和、有限温度の場合はアニーリングと呼ぶことにする。
アニーリングMDの実行
さっそくアニーリングのMDを実行させてみよう。スクリプト・ファイルbccFe_anneal.lcmがscriptディレクトリーに存在することを確認した上で、以下のようにLAMMPSを実行しよう。
$work> lmp_serial -in script/bccFe_anneal.lcm
プログラムが無事終了したら、bccFe_anneal.outというファイルが作業ディレクトリにできるはずだ。OVITOを起動し、このファイルをOVITOの画面にドラッグすると、そのファイルがOVITOに読み込まれる。その後、OVITOの画面下部の再生ボタンを押せば下のような動画が見られる。300 Kで室温に近い温度なので、よく見ると粒子の位置が揺らいでいるのが確認できるだろう。
Video 1: アニーリングMDから得られた動画
このMDは絶対0度で緩和した鉄の結晶に温度入れて時間発展させることによって実現することができる。動画でわかるようにアニーリング自体はさほど興味深い現象ではない。しかし、これは照射や変形などの興味深い現象を扱う場合の初期値となるのでこのシミュレーション法を確実に理解しておく必要がある。
スクリプトの説明
構造緩和を行う
さて、前回と同様にスクリプトを最初から見ていこう。スクリプトの前半を一気に下に示す。
units metal
boundary p p p
atom_style atomic
lattice bcc 2.83 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region sim_box block 0 25 0 25 0 25 units lattice
create_box 1 sim_box
create_atoms 1 box
pair_style eam/fs
pair_coeff * * ./potentials/Fe_mm.eam.fs Fe
neigh_modify every 1 delay 0 check yes
fix 1 all box/relax iso 0.0 vmax 0.001
min_style cg
minimize 0.0 1.0e-12 1000 10000
これは前回の学んだ鉄結晶の構造緩和のスクリプトである。このスクリプト部の中身はすでに説明したので、次に行く。ここがよくわからない受講者は一つ前の講座に戻って復習していただきたい。
アニーリングを行う準備
unfix 1
variable set_temp equal 300.
velocity all create ${set_temp} 1582775
fix 1 all nve
fix 2 all temp/rescale 10 ${set_temp} ${set_temp} 0.1 1.0
unfixコマンドはすでに実行されたfixコマンドを無効にするためのコマンドである。この例では、上の構造緩和の過程で実行されたfix ... box/relaxが有効であり、この後のアニールの過程では不要になるためそれをここで無効化している。unfixの引数の"1"は無効化したいfixコマンドの識別名である。
variableコマンドは、前の講座ですでに出てきたが、変数を定義するコマンドである。ここでは後でアニール時の温度として使われる値をset_tempという変数で宣言し300に設定している。
velocity ... createコマンドは原子にランダムな速度を与えるためのコマンドである。固体の温度は各原子がどの程度激しく振動しているかで決まる。前の講座で述べたがminimizeコマンドで構造緩和された系は絶対零度になっている。すなわち、原子はまったく振動しておらず、その場に静止している。構造緩和された系を使ってアニールのMDを行うためには温度を引き上げなくてはならない。このコマンドではallで指定されたグループに属する原子（すなわちすべての原子）の速度を温度がset_tempになるように設定する。最後の7桁の整数は乱数のシードである。これを変えれば同じ温度で別のケースを計算できるので結果の再現性をチェックすることができる。
fix ... nveコマンドは、MDの時間積分に関する指定である。ここのnveは原子数、体積、エネルギー一定で時間積分すること（ミクロカノニカルアンサンブル）を意味している。他に、圧力一定で体積を可変とするnptなどを使うこともある。
次のfix ... temp/rescaleコマンドによって、MDが進行している間、ある一定の時間間隔で系の温度を指定された設定温度に強制的に合わせることができる。これを行なうことにより系の温度がMDの時間発展中に設定温度からずれることを防ぐことができる。fixの識別名が"2"になっているが、これはすでに"1"が使われているためである。文法上"1"から順番に使う必要はないが、そうすることによってスクリプトが読みやすくなるだろう。
次の"10"は分子動力学計算における時間ステップを10回繰り返したら強制的な温度合わせを行うことを示している。次の2つの引数はシミュレーション開始時と終了時での温度も設定である。ここでは一定の温度を保つので両方とも${set_temp}となる。系が徐々に温められたり、冷やされたりする場合にはこれらの2つの引数には異なった値を代入することになる。次の引数はウィンドウと呼ばれる値であり、設定温度からどこまで乖離したら強制的な温度調整が行われるかを指定する。この例では0.1°離れたら温度調整されることになる。最後の引数は、温度調整される原子の割合で、ここでは1.0なので調整される場合は全ての原子が調整されることになる。よりコマンドを詳しく知りたい方は、マニュアルを参照してほしい。
アニーリングの実行
さて、最後に残りのコマンドを書き出す。
reset_timestep 0
dump 1 all custom 50 bccFe_anneal.out id type xs ys zs
thermo 200
run 1000
reset_timestepコマンドはMD実行前に時刻を設定するコマンドで、ここでは時刻を0に設定している。このコマンドがなくとも初めは0に設定されるので厳密には不要なのだが、いい機会なのでここでこのコマンド紹介しておく。
すでに構造緩和のコースでdumpが出てきたが、重要なので再度説明する。dumpコマンドは各原子の状態を一定の間隔で出力させるためのコマンドであり、dumpとcustomが組み合わさって1つのコマンドになっていると考えた方がいい。第1引数はdumpコマンドの識別名でここでは"1"を選択している。次のallは系全体を示すグループ名で、ここでは系に存在する全ての原子について出力することを要求している。次の"50"は出力の頻度でありMDの時間ステップが50ステップ進むごとに出力される。指定されない限り、1ステップは1fsになっているので50fsごとに原子の位置などに関する情報がファイルに出力されることになる。
次は出力されるファイル名であり、ここではbccFe_anneal.outである。残りは各原子に対して出力される変数が記されている。説明すると、idは原子1個1個に与えられた通し番号、typeは原子種でありこの例ではFe原子のみがシミュレーションの対象であるので全ての原子で"1"となる。xs、ys、zsは系空間で規格化された原子のx、y、z相対座標値になる。
よって出力されるファイルの形式は
1 1 0.0 0.0 0.0
2 1 0.0 0.0 0.1
...
...
座標位置は全て0.0と1.0の間になる。実際のファイルではボックスの大きさに関する情報が先頭につく。シミュレーションが終了した後にbccFe_anneal.outをOVITOに読み込ませればアニメーションを見れることはすでに述べた。単に、MD結果をOVITOで解析するだけなら、このファイルの中身を詳しく知る必要はない。dump ... customコマンドの詳細はマニュアルを参照してほしい。
MD実行中、系の熱力学的変数は1シミュレーションステップ毎に計算されるわけではない。それら熱力学的変数の計算頻度を指定するのがthermoコマンドである。ここでは時間ステップ200回に１度、熱力学的変数が計算されることになる。また、同じ頻度でMDの途中経過が標準出力に出力される。
さて、ここまでがMDを実行するためのお膳立てであり、やっとrunコマンドで初めてMDが実行されることになる。第１引数は時間発展の総ステップ数である。上でunits metalを宣言されているが、この場合、特別な指定がない限り計算ステップ幅は1fsである。よって今回は1psまで実行されることになる。なお計算ステップ幅は別のコマンドで変えられ、幅を広くすると計算効率が上がるが計算精度は下がる。照射環境下などで出現する高エネルギー状態の原子がなければ通常1fsから変える必要はない。
これでアニーリングのスクリプトの説明は終わりだ。このコースでの時間発展がある分子動力学、すなわちrunコマンドによるMD実行の最初の例である。前回、構造緩和シミュレーションを学んだのはこれを行うためにお膳立てだったわけでる。
目次へ 前は構造緩和を行う 次はせん断変形を行う