今回は、例題を通してLAMMPSのスクリプトの骨組みを紹介していきます。各コマンドの簡単な解説も行います。
LAMMPSコマンドについて
LAMMPSのコマンドの説明は公式がとても丁寧です。困ったらまず参照するといいでしょう。この記事に登場する各コマンドについても詳しい使い方はリンク先を参照してみてください。
スクリプト全体の説明として、改行がコマンドの区切りになります。コマンドは上から順に読みながら実行されていきます。(途中に破綻した文字列があってもそこまでは到達します。)
コマンドは単語区切りになっており、必ず1語目にコマンドの種類、2語目以降にそのコマンドに対応したオプションを書いていきます。
スクリプト内の#記号以降はコメント扱いです。
LAMMPSには様々な種類のコマンドがありますが、全体の流れは
- 初期化、シミュレーション全体に関わる設定
- シミュレーション格子の作成と原子の配置
- MD計算実行中に行う各種制御の追加
- MD計算を実行
という順序になっています。
まずは書いてみる
MD計算実行はrunコマンドですが、runコマンドだけを流すと、シミュレーションする対象も作っていないのに実行できないよ、とエラーが出てきます。
エラーなくrunコマンドを流すための最小限のスクリプトは以下のようになります。まだ原子を配置していないのでほんとうに最小限です。
region region1 block 0 1 0 1 0 1 units box
create_box 1 region1
pair_style lj/cut 1.0
pair_coeff * * 1 1
mass 1 1.0
run 0
ざっくりと説明をしておくと、regionとcreate_boxでシミュレーション格子の作成、pair_*とmassでポテンシャルと質量の定義をしています。
上記のスクリプトを実行すると様々な情報が出力されます。ほとんどはrun 0実行開始時と終了時に出力されるものです。
とはいえ、今回は原子を1つも配置していないため、おもしろい結果はないと思います。
構造緩和計算: Fe結晶の安定構造
スクリプト
次はもう少し内容のある例をやってみましょう。今回の例題はFe結晶の安定構造の探索と、物性値の解析です。
MDを行う上で(局所)安定構造を考える場面は数多く出てきます。結晶の多形構造のエネルギーの差を見積もったり、NEB計算を行う際の出発点と終点に設定したりなどです。
ここではよくある問題設定として、適当に作成した初期構造からスタートして緩和計算を行い、安定構造を探ってみましょう。
atom_style atomic
units metal
boundary p p p
lattice bcc 2.9
region region1 block 0 3 0 3 0 3 units lattice
create_box 1 region1
create_atoms 1 box
pair_style eam/fs
pair_coeff * * Fe_mm.eam.fs Fe
mass 1 55.845
thermo_style custom step etotal temp lx vol press pxx pyy pxy
thermo 10
fix f1 all box/relax iso 0.0
minimize 1.0e-10 0.0 1000 100000
MD計算を実行するrunがないじゃないかと思われるかもしれませんが、ここではminimizeというコマンドがrunの亜種に相当すると理解してください。
はじめのブロックは初期設定に関わるものです。原子の取扱方法、単位系、境界条件を設定しています。
エネルギーや長さの次元がポテンシャルによって決まっているため、単位系は使うポテンシャルの指示に合わせましょう。
境界条件はMDでは基本的に周期境界条件を適用します。原子スケールの現象では表面を出すことそのものが物性に関わってきてしまうため、周期境界条件によってこれをうまく避けている形になります。一方でこのままでは本来は外部から影響を受けるべき何か(温度や圧力など)が入れられなくなってしまうので、MDではそれらはアンサンブルという形で組み込んでいます。アンサンブルについては次回解説します。次のブロックの
create_boxでシミュレーション格子の作成を行っています。create_boxの後の1という数字は系に存在する原子の種類の数を指示しています。LAMMPS内部の原子のデータの扱いに関わるため、この段階で決めておきます。create_atomsで原子を追加します。後ろのオプションで1個ずつ置くこともできますが、latticeコマンドで結晶格子を指定しておくとこのようにシミュレーション格子全体を一気に埋めることができます。
create_boxと異なり、このコマンドは原子を追加したければ何度でも呼ぶことができます。
read_*コマンドを使って外部のファイルを読み込むことで複雑な構造を展開することもできます。後の回でまた詳しく書きます。4番目のブロックはポテンシャルと質量の設定です。今回はEAM/FSポテンシャル(Finnis-Sinclairポテンシャル)を使っています。
このように、ポテンシャルの種類を決めるpair_styleとポテンシャルパラメータを決めるpair_coeffはいつもペアで使います。pair_coeffはポテンシャルの種類によって書き方が変わり、パラメータを直書きしたりファイルを指示したりします。今回は後者です。5番目のブロックの
thermo_styleとthermoは画面に出力する情報を決めています。MD実行中、10ステップおきに温度や圧力など様々な統計量を出力するようにしています。6番目のブロックの
fixでminimize計算実行中の格子サイズの変形を指示しています。
正確に表現すると、ここで定義したfixがminimize計算中に呼ばれて格子の変形を行います。最後に
minimizeでエネルギーの緩和計算を行います。minimizeはrunコマンドと似ていますが、動力学計算を行うのではなく、よりエネルギーが低くなるように原子を動かすコマンドです。
なので、minimizeが終わったときはエネルギーが(局所)安定点に達した構造が得られることが期待されます。つまり0ケルビンのときの構造です。デフォルトはCG法が使われています。
fixコマンドについてもう少し詳しく説明しておきます。
fixコマンドは、コマンド中3番目のオプションで指定する何らかの「制御」を系に追加するものです。この「制御」はMD計算実行中にどんな計算をするかを決めるもので、とても重要なコマンドです。
具体的には、特定の原子に外力を与えたい・固定したい・移動を制限したいといった原子の動きに関する制御、NVTアンサンブルやNPTアンサンブルといった系全体の温度や体積を変化させる制御、また一定時間おきに何らかの値の標準出力やファイル出力を行うものなど多様なものが含まれており、fixだけで副コマンド群を形成しています。
今回のfix box/relaxコマンドはminimize計算中に系が安定になるように(→圧力が0になるように)体積を変化させるというものです。
実行結果
なお、今回のスクリプトを実行する際にはポテンシャルファイル"Fe_mm.eam.fs"が同じディレクトリに必要になります。これはLAMMPSのソースファイルのpotentialsディレクトリ以下に置いてありますので、コピーして持ってきましょう。
実行中の結果の出力は以下のようになると思います。
Step TotEng Temp Lx Volume Press Pxx Pyy Pxy
0 -221.85018 0 8.7 658.503 -79218.494 -79218.494 -79218.494 -1.9628292e-11
10 -221.94489 0 8.6913 656.52947 -74541.381 -74541.381 -74541.381 -4.4582695e-11
20 -222.0336 0 8.6826 654.55988 -69760.564 -69760.564 -69760.564 2.2597797e-11
30 -222.1162 0 8.6739 652.59423 -64874.803 -64874.803 -64874.803 -2.2481665e-11
40 -222.19259 0 8.6652 650.63253 -59882.527 -59882.527 -59882.527 2.3569327e-12
50 -222.26266 0 8.6565 648.67476 -54791.317 -54791.317 -54791.317 -3.1896949e-11
60 -222.32635 0 8.6478 646.72092 -49660.404 -49660.404 -49660.404 9.6674003e-12
70 -222.38366 0 8.6391 644.77101 -44502.086 -44502.086 -44502.086 -2.218868e-11
80 -222.43456 0 8.6304 642.82502 -39314.184 -39314.184 -39314.184 -5.633681e-12
90 -222.47906 0 8.6217 640.88296 -34097.681 -34097.681 -34097.681 2.9908225e-11
100 -222.51714 0 8.613 638.9448 -28853.033 -28853.033 -28853.033 -2.1703726e-11
110 -222.54879 0 8.6043 637.01056 -23578.707 -23578.707 -23578.707 1.5543501e-11
120 -222.57401 0 8.5956 635.08023 -18276.072 -18276.072 -18276.072 2.0280005e-11
130 -222.59279 0 8.5869 633.1538 -12945.259 -12945.259 -12945.259 -1.8541946e-11
140 -222.60511 0 8.5782 631.23127 -7584.2829 -7584.2829 -7584.2829 -9.7042885e-12
150 -222.61096 0 8.5695 629.31263 -2192.8485 -2192.8485 -2192.8485 3.1661592e-11
156 -222.6115 0 8.5659746 628.53627 -1.7710042e-05 -1.7710045e-05 -1.7710055e-05 1.0170402e-12
縦の位置がずれていますが、上の"Step TotEng Temp..."が下の数字と対応しています。
この最後の値は別の初期値(別の格子サイズなど)から始めても同じ結果になることが期待されます。
実行すると、エネルギーが低下していくこと、格子が小さくなっていくこと、圧力が0に近づいていくことなどがわかります。
単位はunitsコマンドで指定したものが出てくるので、単位系で迷ったらマニュアルを見ましょう。
圧力は序盤でとてつもなく大きな値(GPaのオーダー)が出ていますが、これは正常な振る舞いです。マクロな世界で鉄を1%縮める操作がどれほど大変かを考えてみればわかります。
これはまだ原子を安定位置に動かしただけのいわゆるMolecular Static計算で、動力学計算を行ってはいないのですが、この状態でもわかる物性値はたくさんあります。
例えばFeの凝集エネルギーや格子定数などが求まっています。体積と原子数がわかっているので密度も計算することができます。(おおよそ7.9g/cm^3と求まり、実験結果と対応します。)
次回は結晶を溶かす例を題材にダイナミックな分子動力学計算を行ってみます。可視化にも挑戦していきます。
← LAMMPS入門の手引き 0:環境の設定
→ LAMMPS入門の手引き 2:動的な計算
応用課題
例題に関連した話題を載せておきます。
ここから原子を1つ除去すると、結晶中に空孔ができ、格子欠陥となる。Feの格子欠陥のエネルギーを求めてみよう。
得られた構造から原子を動かさず格子を1方向にだけ伸ばすと真空領域ができる。周期境界条件なので、これは結晶と真空がウエハース上に重なった構造となる。ここから表面エネルギーを求めてみよう。
系全体をx軸方向に圧縮させるなど微妙に変形させると、ひずみが生じ、圧力が0でなくなる。ここからFeの弾性定数(C11, C12, C44)を求めてみよう。