LAMMPSで構造最適化
原子間ポテンシャルに従ってMD計算をするだけでなく構造最適化もできるらしい
原子位置はMD計算中に緩和されるとしても、格子定数の最適化は別で行うしかなさそう？
キーとなるコマンド
minimize: 構造最適化を実行fix box/relax: 外圧を加え、計算モデルの大きさと形も最適化するように変更
Wurtzite-GaNを構造最適化(失敗)
から取得できるw-GaN構造を初期構造として使用
LAMMPSのpotentialsディレクトリにいくつかプリセットがあるのでその中でGaN.swを使用して最適化してみます
## 基本設定 units metal dimension 3 boundary p p p atom_style atomic newton on ## 計算モデル作成 variable a equal 3.216 variable c equal 5.240 variable b_x equal $a*(-0.5) variable b_y equal $a*sqrt(3.0)/2.0 variable 1_3 equal 1.0/3.0 variable 2_3 equal 2.0/3.0 lattice custom 1.0 & a1 $a 0.0 0.0 & a2 ${b_x} ${b_y} 0.0 & a3 0.0 0.0 $c & basis ${1_3} ${2_3} 0.0 & basis ${2_3} ${1_3} 0.5 & basis ${1_3} ${2_3} 0.375 & basis ${2_3} ${1_3} 0.875 region myreg block 0 5 0 5 0 5 create_box 2 myreg create_atoms 2 box & basis 1 1 & basis 2 1 & basis 3 2 & basis 4 2 ## 原子をGa,Nに設定 mass 1 69.72 mass 2 14.01 ## 原子間ポテンシャルをStillinger-Weberに設定 pair_style sw pair_coeff * * GaN.sw Ga N ## 構造最適化実行 timestep 0.0001 # 0.1fs thermo_style custom step etotal temp press lx ly lz vol thermo 1000 # aniso: xyzそれぞれに等しく外圧を加える&長さ比は保存しない # couple: 指定した軸の長さ比を保存 # 外圧0, 1回のiterationで変化する体積は最大1% fix 1 all box/relax aniso 0.0 couple xy vmax 0.01 # エネルギーによる収束条件なし, 原子間力他の収束条件を使用 # 最大で10,000step, 1,000,000eval minimize 0.0 1.0e-20 10000 1000000
実行ログより一部抜粋
LAMMPS (16 Mar 2018)
Lattice spacing in x,y,z = 4.824 2.78514 5.24
Created orthogonal box = (0 0 0) to (24.12 13.9257 26.2)
...
Step TotEng Temp Press Lx Ly Lz Volume
0 -2877.0956 0 73125.066 24.12 13.925688 26.2 8800.2553
1000 -3100.5642 0 212.02536 24.548929 14.173331 26.081804 9074.9054
2000 -3100.7013 0 0.40261942 24.562407 14.181112 26.076651 9083.0776
3000 -3100.7028 0 -0.00062158384 24.562064 14.180914 26.077135 9082.9924
4000 -3100.7028 0 -3.3787783e-06 24.562064 14.180914 26.077135 9082.9924
5000 -3100.7028 0 5.4301294e-10 24.562064 14.180914 26.077135 9082.9924
6000 -3100.7028 0 3.1547809e-09 24.562064 14.180914 26.077135 9082.9924
7000 -3100.7028 0 8.1352718e-10 24.562064 14.180914 26.077135 9082.9924
8000 -3100.7028 0 1.5203109e-09 24.562064 14.180914 26.077135 9082.9924
9000 -3100.7028 0 3.0510664e-09 24.562064 14.180914 26.077135 9082.9924
10000 -3100.7028 0 3.1167627e-09 24.562064 14.180914 26.077135 9082.9924
Loop time of 2.77113 on 144 procs for 10000 steps with 750 atoms
99.5% CPU use with 144 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = max iterations
Energy initial, next-to-last, final =
-2877.0955666 -3100.70275883 -3100.70275883
Force two-norm initial, final = 727.066 4.09501e-11
Force max component initial, final = 450.072 3.9289e-11
Final line search alpha, max atom move = 1 3.9289e-11
Iterations, force evaluations = 10000 19987
なお、LAMMPSでは実行時に計算モデルが直方体となるように自動で変換されます
今回は初期構造として底面の辺=a, 高さ=cの菱形柱の単位格子を設定したが、
x=3a/2, y=√3a/2, z=cの直方体に変換されてます
最適化前後の構造を可視化したところ、以下のようになりました(OVITOで可視化)
不自然な曲がり方をしている。。
上の曲がった構造を周期的境界条件(PBC)に従って拡張すると原子同士が重なったりぽっかり穴が開いたり……もしかしてPBCが動作してない？
うまくいかない原因は…？
サンプルコードがどこかにないかよく探しました
なんだかんだでこれが一番大事で近道だったりします！！
今回はLAMMPSのexamplesディレクトリの中から以下の条件に当てはまるサンプルを探しました。
minimizeコマンドを使用pair_style swコマンドを使用
最終的な目的は違いますが、neb/in.neb.sivacが当てはまったのでこれを参考にしました
minimizeと合わせるとPBCが機能しない？
=>サンプルコードを動かしたところ、物質はGaNではなくSiだが構造が崩れることなくちゃんと動作しているので違う初期状態で速度を持ってしまっている？
=>velocityコマンドで明示的に速度0にしたが変わらずfix box/relaxコマンドで与えた外圧が悪さしている？
=>minimizeコマンドだけでも構造が崩れるので違うStillinger-WeberポテンシャルではWurtziteを再現できない？
=>むしろWurtziteのような4配位構造を再現するためのポテンシャルtimestepやmaxstep、
minimizeコマンドの収束条件などが不適？ =>いくつか試したが変わらず初期構造がオカシイ
=>自分で直方体の単位格子を取り直し、初期構造を設定したところちゃんと動作した！！
Wurtzite-GaNを構造最適化(今度こそ成功)
LAMMPSスクリプト
# 基本設定 units metal atom_style atomic boundary p p p # 単位格子 variable a equal 3.216 variable b equal $a*sqrt(3.0) variable c equal 5.240 variable 1_2 equal 1.0/2.0 variable 2_3 equal 2.0/3.0 variable 1_6 equal 1.0/6.0 variable 3_8 equal 3.0/8.0 variable 7_8 equal 7.0/8.0 lattice custom 1.0 & a1 $a 0.0 0.0 & a2 0.0 $b 0.0 & a3 0.0 0.0 $c & basis 0.0 0.0 0.0 & basis ${1_2} ${1_2} 0.0 & basis 0.0 ${2_3} ${1_2} & basis ${1_2} ${1_6} ${1_2} & basis 0.0 0.0 ${3_8} & basis ${1_2} ${1_2} ${3_8} & basis 0.0 ${2_3} ${7_8} & basis ${1_2} ${1_6} ${7_8} region myreg block 0 5 0 5 0 5 create_box 2 myreg create_atoms 2 box & basis 1 1 & basis 2 1 & basis 3 1 & basis 4 1 & basis 5 2 & basis 6 2 & basis 7 2 & basis 8 2 mass 1 69.72 mass 2 14.01 # ペアポテンシャルを設定 pair_style sw pair_coeff * * GaN.sw Ga N # 構造最適化 thermo_style custom step etotal temp press lx ly lz vol thermo 10 fix 1 all box/relax aniso 0.0 vmax 0.01 minimize 0.0 1.0e-20 1000 10000
最適化前後の構造は以下のようになりました(格子定数が変わっただけなのでぱっと見は変化なし
SWポテンシャルの下で構造最適化されたw-GaNの格子定数はそれぞれ以下のようになりました
a: 3.216 -> 3.190597
c: 5.240 -> 5.2102232
以上です！！
