構造緩和を行う
はじめに
このコースでは実際にLAMMPSを使って簡単な計算を行う。LAMMPSにどんな計算をやらせるかは入力スクリプト（以下簡単にスクリプトと呼ぶ）に記述することになる。スクリプトはある文法に従って作成する必要があり。要するにLAMMPSにやらせる作業は、LAMMPSの実行中に選択するのではなく、実行する前にLAMMPS特有の言語でスクリプトを作っておき、LAMMPS開始時にそのスクリプトを読み込み、それを順番に実行していく、ということだ。
この講座では、その最初の例、鉄の結晶を作り凝集エネルギー(Cohesive energy)を求めるシミュレーションを紹介する。凝集エネルギーとは原子1個を結晶から剥ぎ取って真空に置くのに必要なエネルギーである。鉄は低温から常温では体心立法格子(BCC)構造になる。以下ではLAMMPSに、シミュレーションを行う空間を作って、BCC構造の結晶をそこに入れて、鉄の原子間ポテンシャルで構造緩和する、という一連の作業をさせてみよう。
構造緩和計算の実行
前の講座でLAMMPSインストールの後にscriptディレクトリをダウンロードしたはずだ。そこにはこのコースで使ういくつかのスクリプトが格納されている。その中にbccFe_relax.lcmというファイルがあるが、これが今回使うスクリプトである。さっそく、このファイルを使ってLAMMPSを実行してみよう。LAMMPSの実行は、使っているOSで標準入出力のアプリを実行させ、インストール時に作成したworkディレクトリに入って、標準入力に
$work> lmp_serial -in script/bccFe_relax.lcm
とタイプすればよい。-inはスクリプトを指定する場合に使われるコマンドライン・フラグ（以下簡単にフラグと呼ぶ）である。ファイルの名前の形式は何でもよい。LAMMPS開発者側が提供するスクリプトの例ではin.???というファイル名形式になっているが、これだとWindows環境における拡張子による選別ができない。筆者はLAMMPS commandの略でlcmを使って、???.lcmというファイル名形式を使っている。このコース全体を通じて、このlcm<という拡張子のファイルはのスクリプトを意味することとする。
さて、上の例のように、一般的にLAMMPSは実行形式のファイル（ここではlmp_serial）、スクリプト（ここではbccFe_relax.lcm）、ポテンシャルファイル（ここではFe_mm.eam.fs）の３つのファイルがあればLAMMPSでMDを実行させることができる。
上の例では実行が終了すると最後に以下のように出力されるはずである。
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
-128682.636737 -128826.096944 -128826.096944
Force two-norm initial, final = 32116.9 0.0014228
Force max component initial, final = 32116.9 0.0014228
Final line search alpha, max atom move = 0.00549093 7.8125e-06
Iterations, force evaluations = 11 23
...
...
**** Cohesive energy : 4.12243510220317. *****
Total wall time: 0:00:03
最初にMinimization stats:という部分がある思われるが、その次の文でエネルギー変化が0で最適化計算がこれ以上できなくなってしまったことが記述されている。色々面倒なことが書いてあるが、とりあえず今は構造緩和計算に成功したことがわかれば良い。そのあと、求まった凝集エネルギー(約4.12eV）が出力されている。Total wall timeは計算にかかった時間で、環境によって変化する。
どうだろう、これまでMDは面倒だと思っていた受講者もいるかもしれないが、必ずしもそうでは無いとお分かり頂けただろうか。
エラーが起きて上のような標準出力が得られなかった場合は、スクリプトまたはポテンシャルのファイルが所定の場所にない可能性が高いのでそれを確認しよう。
さて、上の例では絶対0度で鉄結晶がどのような構造に落ち着くかというエネルギー最適化計算を行ったことになる。そのため、この計算には時間の概念もない。このような計算法は厳密には分子動力学(Molecular dynamics)ではなく、分子静力学(Molecular statics)と呼ばれる。この例をあえて紹介したのは、後で説明するように初心者にとって重要なコマンドが凝縮されているとともに、今後の講座でより複雑な計算法を学ぶときに、その初期状態を作るために必要だからである。
スクリプトの説明
環境の設定
プログラムが正常に終了したら、今回使ったスクリプトであるbccFe_relax.lcmの内容を1行ずつ見ていくことにする。まずは最初の4行を以下に表示する。
# Relaxing iron crystal
units metal
boundary p p p
atom_style atomic
#から始まる行はすべてコメントになる。すなわち最初の行はなくてもプログラムは正しく実行される。
unitsコマンドは単位の設定を行う。固体の材料を扱っている場合はmetalと呼ばれる単位のセットを使うことをお勧めする。例えば、エネルギーの単位はeV、距離の単位はAngstrom、時間の単位はピコ秒、圧力の単位はbarという具合である。
boundaryコマンドは境界条件を指定する。パラメータはp p pとなっているが、これは(x,y,z)の3方向に対して周期的境界条件が設定されていることを示す。pのほかにfやsなどの設定があり詳しくはマニュアルに書いてあるが今は気にする必要はない。
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
latticeコマンドは結晶系を定義するのに使われる。ここでは結晶のミラー指数が[100]、[010]、[001]である3方向をそれぞれx、y、z方向の単純な直交座標系に置き、格子定数2.83AのBCC（体心立法）結晶系を定義している。このコマンドでは、この後使用する結晶系を定義しただけで、まだ計算に使用する原子を定義したわけではないことに注意しよう。
regionコマンドは領域を定義する。ここではsim_boxという名前（名前はなんでもいい）の領域を定義する。blockはこのコマンドの詳細を規定するキーワードで、ここでは領域の形が直方体であることを示し、その後に続くx_min、x_max、y_min、y_max、z_min、z_maxの6つの境界によって定義される。unitsもこのコマンドのキーワードで、x_min等に代入された値の単位を規定するための属性でここではlatticeが指定されて格子定数が単位に選択されている。すなわちここで定義された領域は一辺が25a(aは格子定数)の立方体である。
create_boxコマンドでは上で定義したsim_boxという名前のregionを計算空間（ボックス）として定義する。ちょっと2度手間な感じがするが私が知らない理由があるのだろう。このコマンドの第1引数の"1"は原子種の数を表している。ここでは、鉄原子のみを扱うのでこの数値になる。
create_atomsコマンドでは上で定義したボックスに原子を入れる。上のコマンドと同様に、第1引数の"1"は原子種の数を表している。第2引数のboxが選択されると、create_boxですでに定義された計算空間全体に、latticeコマンドで定義された結晶系を作るように原子が配置される。すなわち、このコマンドの前にlatticeコマンドとcreate_boxコマンドが実行されていなくてはならない。
これで一辺が25aの立法体の計算空間の中にBCC構造の原子群が定義されたことになる。上で説明したように、その空間には周期的境界条件が課せられているので擬似的に無限の空間となる。これで計算開始時の初期状態が出来上がったことになる。
原子間ポテンシャル
pair_style eam/fs
pair_coeff * * ./potentials/Fe_mm.eam.fs Fe
pair_styleコマンドは原子間ポテンシャルのタイプを選択する。MDでは、与えられた初期状態が原子間の相互作用によって時間発展していく。この原子間の相互作用を定義するのが原子間ポテンシャルである。この例では引数がeam/fsとなっており、Embedded Atom Method（EAM）ポテンシャルグループの中のFinnis-Sinclare（FS)タイプのポテンシャルが選択されている。原子間ポテンシャルの詳細については後のコースで学ぶので、今は金属においてFSタイプがよく使われることを述べるにとどめる。
pair_coeffコマンドでは原子間ポテンシャルのファイルを指定する。このファイルはpair_styleコマンドで設定されたタイプのものでなければならない。このファイルの中身は数値テーブルである。最初の2つの引数はポテンシャルを適用する原子種の組み合わせを定義している。ここで* *は、すべての組み合わせについて適用することを示している。通常、これらは変える必要はない。3つ目の引数はポテンシャルファイルの位置、4つ目は原子種の定義である。ここでは原子種がFe（すなわち鉄）であるとしている。
その他の準備
dump 1 all custom 100 bccFe_relax.out id type xs ys zs
neigh_modify every 1 check yes
fix 1 all box/relax iso 0.0 vmax 0.001
dumpコマンドは4ワード目（すなわち第3引数）とペアになったコマンドとして認識してほしい。すなわち、ここではdump ... customコマンドと考える。これはMD進行中に各原子の状態をカスタマイズされたフォーマットで出力させるために使われる。
注意して欲しいのは、dump ... customコマンドはあくまでデータの出力の仕方を定義するコマンドで、実際の出力は初期状態から原子運動の計算が開始されてから行われることになる。最初の引数"1"はこのdumpコマンドの識別名である。dumpコマンドが2つ以上ある場合に役に立つが、この例では1つしかないので大きな意味はない。次の引数allはグループの識別名でこの場合はボックス内のすべての原子をひとまとめにしたグループを意味する。一般にグループはスクリプト内のgroupコマンドによって定義されるが、allはすべての原子を意味するグループでありスクリプトで定義しなくても常に存在する。
customは自分で以下に続く複数の引数を用いて出力フォーマットをカスタマイズすることを意味する。ここでは、計算サイクル100回ごとに出力し、その出力先はbccFe_relax.outというファイルで、出力の内容は原子のx、y、z座標で"xs ys zs"となっているのは計算空間のサイズで正規化した数値（すなわち0から1の間の数値）で出力することを求めている。
なお、出力されたファイルは人間が読むというよりも主に計算結果を計算機で可視化する場合に使われる。ここではbccFe_relax.outをOVITO等のプログラムに読み込ませれば可視化できる。可視化については後の講座で少しずつ触れていくことにする。
上で少し触れたが、dump ... customコマンドは2つ以上定義することができて、別のファイルに別の形式や内容の原子データを出力させることもできる。dump ... customコマンドは非常に多機能なコマンドなので一度で全てを理解する必要はない。今はこの例のようにすれば各原子位置のデータがファイルに出力するためのコマンドと理解しておけば良い。
neigh_modifyコマンドは原子間ポテンシャルを適用する原子ペア・リスト（Neighborリスト）の更新の設定である。このコマンドは初心者にはちょっと難しいので、気にしないで次に進んでもらいたい。後の講座で説明する。
fixコマンドは第3引数と組み合わせてMDにおける様々なオプションを指定する。ここではfix ... box/relaxコマンドと理解すれば良く、緩和計算のエネルギー最小化の過程に置いて系の体積緩和も行うことが要求されている。第1引数の"1"はこのfixコマンドの識別名であり、上のdumpの識別名と同様の役割を果たす。次のallは全ての原子をひとまとめしたグループで、dumpコマンドで出てきたallと同様にこのfix ... box/relaxコマンドが適用される範囲を指定していて、この場合は全ての原子にこのコマンドが適用されることになる。
次のisoは等方的な体積緩和をさせることを意味し、それに続く"0.0"は圧力が0.0 barになるように緩和させることを意味する。単位のbarは最初のunitsコマンドで選択されている。その後のvmaxは一度の体積緩和計算ステップでどの程度の体積変化までを許すのかを指定するキーワードで、ここでは一度にできる体積変化は0.001までということになる。経験上、vmaxは0.001にしておけば問題ない。dumpコマンドと同様にfixコマンドは非常に多機能であり、いますぐこのコマンドの全てを理解する必要はなく、いろいろ経験しながら少しずつ慣れていくことにしよう。
構造緩和計算を行う
min_style cg
minimize 0.0 1.0e-12 1000 10000
min_styleコマンドは緩和計算におけるエネルギー最小化のアルゴリズムを指定するコマンドである。ここではcg、すなわち共役勾配法(conjugate gradient法)が選択されている。
ここまで準備が整ったらいよいよ目的の構造緩和計算を行うことができる。それを行うのがminimizeコマンドである。第1引数はエネルギーによる収束条件を示す。エネルギー変化の相対値（全エネルギーに対する割合）がこの値よりも小さいと収束する。ただし、この例では値が"0.0"となっている。第2引数は力による収束条件である。任意の原子のx、y、z方向のすべての力成分が設定値を下回れば収束する。ここで力の単位は上のunitsコマンドで定義されたものである。この例では"1.0e-12" eV/Aが設定値として定義されている。第3引数は最大の繰り返し計算の数である。この値を超えて繰り返し計算を行わない。第4引数は最大の力の評価回数である。通常、力の評価回数は、通常繰り返し計算よりも大きくなる。第3または第4引数によって計算が止まった場合は、系は収束しておらず構造緩和していない状態であることに注意しよう。
結果の出力
variable e_coh equal (-1)*pe/count(all)
print "**** Cohesive enegy : ${e_coh}. *****"
variableコマンドはスクリプト内で使用する変数を定義するコマンドであると考えればよい。ここではe_cohという変数を全エネルギーpeを全原子数で割った値として定義している。peは定義しなくても常に存在する変数で系の全ポテンシャルエネルギーを表す。また、count()は原子数を返す関数であり、"all"はすでに説明したように系全体を表すグループ識別名であり、count(all)は系全体の原子数になる。よってこの文では、e_cohという変数を凝集エネルギーとして定義していることになる。
printコマンドは"..."で囲まれた部分を標準出力に出力するコマンドである。ただし${e_coh}の部分は上で定義した変数e_cohを参照していて、その値が出力される。ここで注意して欲しいことは、LAMMPSのスクリプトでは変数aを参照する場合には、単にaではなく${a}と記述する必要がある、ということである。
これでスクリプトの説明を終える。いかがだっただろうか？各コマンドの詳細な引数は気にする必要はないが、もう一度このコマンドの説明のトップに戻って大まかな流れを復習することをお勧めする。
目次へ 前はインストール 次はアニーリングを行う