====== 概要 ======
結晶構造の最適化を行います。
* 内部座標の最適化 ... 格子定数は変化させず、単位胞の原子位置のみ動かす
* 格子定数の最適化 ... 結晶座標での原子位置は変化させず、格子定数のみ動かす
の2通りがあります。
どちらもエネルギーが最小になるような原子位置あるいは格子定数を探します。
スピン軌道相互作用を入れた場合は構造最適化することができないので注意してください。
====== 内部座標の最適化 ======
内部座標の最適化には、MSR1aという方法を使うものと、min_lapw という名前のプログラムを使うものの2通りがあります。
min_lapw は過去に使用されていたもので、現在は使用を推奨されていませんので、MSR1aを使った方法を紹介します。
- 内部座標最適化用のオプション -min をつけてSCF計算を実行します。原子に加わる力を収束させるためにオプション -fc をつけて実行してください。
$ init_lapw -b -numk 1000 -rkmax 7.0
$ run_lapw -ec 0.00001 -cc 0.0001 -fc 0.1 -p -min
* 構造緩和の途中でマフィンチン半径が隣の原子と重なるために、エラーが出て止まることがあります。
Mixer - Error. no feasible step small enough, check RMT and model
この場合は
$ reduce_rmt_lapw -r 3
$ run_lapw -ec 0.00001 -cc 0.0001 -fc 0.1 -p -min
としてください（reduce_rmt_lapw は、オプション -r X でマフィンチン半径をX%減らします。適当に増やしていって、エラーが出なくなったらこの設定を外して再計算してください）
* 内部座標最適化の計算が終われば case.struct は新しい構造のものに自動的に置き換わります
* もとのstructファイルは init_lapw 実行時に生成された case.struct_orig で見ることができます。
* 構造最適化を中断したい場合は .minstop という空のファイルをおいてください。原子位置の変更が行われず、通常のSCF計算のみが続行されます
$ touch .minstop
- setrmt_lapw case でマフィンチン半径を再計算します。違う値になるようであれば、新しい半径で再びその構造で最初から（init_lapwから）最適化を行ってください
- k点数を変えて最適化を行い、収束を確認してください。
各原子に加わる力のベクトルは、以下のようにして確認できます。FGL001は1番目の原子に加わる力（単位はmRy）です
$ grep ":FGL001:" case.scf
====== 格子定数の最適化 ======
必要に応じて、格子定数の最適化もできます。
格子定数を少しずつ変化させたstructファイルを複数用意して、内部座標の最適化を行いつつエネルギーを計算し、最後にエネルギーを2次関数でフィッティングすることで、もっとも安定となる構造を推定します。
* 計算には非常に時間がかかる割にそんなに精度が出ないので、単純な格子構造を除いてあまり使用は推奨されません。
* 事前に内部座標の最適化を行ってください。
* ブラベー格子が三斜晶系 (Triclinic) となるものには対応していません。
* WIEN2kがデフォルトで使うPBE汎関数の特徴として、体積を固定せずに計算する場合は実験値よりも格子定数が大きな値になってしまう傾向があります。
* init_lapw の実行時にオプション -vxc PBESOL をつけることで、格子定数の値を改善してくれるPBEsol汎関数を使うことができます。$ init_lapw -b -numk 1000 -rkmax 7.0 -vxc PBESOLただしPBEsolは凝集エネルギー、解離エネルギーの値があまり正しく出せないので、格子定数を決める以外での使用は推奨されていません。
- 構造を少しずつ変化させたときのstructファイルを作成します。
$ x_lapw optimize
********************************************
GENERATES STRUCT-FILES AND optimize.job
PLEASE CHOOSE ONE OF THE FOLLOWING FEATURES:
[1] VARY VOLUME with CONSTANT RATIO A:B:C
[2] VARY C/A RATIO with CONSTANT VOLUME (tetr and hex lattices)
[3] VARY C/A RATIO with CONSTANT VOLUME and B/A (orthorh lattice)
[4] VARY B/A RATIO with CONSTANT VOLUME and C/A (orthorh lattice)
[5] VARY A and C (2D-case) (tetragonal or hexagonal lattice)
[6] VARY A, B and C (3D-case) (orthorhombic lattice)
[7] VARY A, B, C and Gamma (4D-case) (monoclinic lattice)
[8] VARY C/A RATIO and VOLUME (2D-case) (tetr and hex lattices)
********************************************
* 扱っている物質のブラベー格子や実験結果に合わせて、動かしたいパラメタを含んでいるものを選んでください
* [1] A:B:Cを保ちながら、体積を変化させます。
* [2] 体積を保ちながら、C/Aを変化させます（正方晶系と六方晶系）
* [3] 体積とB/Aを保ちながら、C/Aを変化させます（斜方晶系）
* [4] 体積とC/Aを保ちながら、B/Aを変化させます（斜方晶系）
* [5] AとCを変化させます（正方晶系と六方晶系、フィッテイング変数:2）
* [6] A, B, Cを変化させます（斜方晶系、フィッテイング変数:3）
* [7] A, B, C, γを変化させます（単斜晶系、フィッテイング変数:4）
* [8] C/Aと体積を変化させます（正方晶系と六方晶系、フィッテイング変数:2）
* A, B, C, γはブラベー格子の格子定数です
* 各項目ごとに対話形式で入力していきます
* 格子定数を動かすパーセンテージは最大でも±10%以内に収まるようにします。
* エネルギー的に安定な解が、この格子定数を動かした範囲内に入っている必要があります。
* 元のstructファイルは case_initial.struct という名前で複製されます。
- 設定が終わると、格子定数を変更させた複数のstructファイルと、最適化計算用のスクリプト optimize.job ができます。optimize.job を開いて run_lapw の行を自分の収束条件に合わせて変更します。
- 前回の計算履歴を削除し、optimize.job を実行します
$ rm *.broyd*
$ ./optimize.job
* 構造を変化させたときにマフィンチン半径が隣の原子とがぶってしまい、エラーが出ることがあります。その場合は
$ rm -f case_*.struct
$ reduce_rmt_lapw -r 3
でマフィンチン半径の減らして再度 x_lapw optimize を実行します。あるいは x_lapw optimize を再度実行して格子定数の変化のパーセンテージを小さくしてください。
- （[1]-[4]を選んだ場合）eplot_lapw を使うことで、もっとも安定になる格子定数を推測します
$ eplot_lapw -a " " -t vol/coa/boa
* オプション-tで変化させたパラメタを指定します。
* -t vol ... 体積を変化させた場合。[1]
* -t coa ... C/Aを変化させた場合。[2], [3]
* -t boa ... B/Aを変化させた場合。[4]
- （[5]-[8]を選んだ場合）parabolfit_lapw を使うことで、もっとも安定になる格子定数を推測します
$ parabolfit_lapw -t 2/3/4
* オプション-tでフィッテイングする変数の数を指定します。
* -t 2 ... [5], [8]
* -t 3 ... [6]
* -t 4 ... [7]
- 得られた格子定数を使って、再度内部座標の最適化を行います
$ vim case_initial.struct
（structファイルを編集して、新しい格子定数に置き換えます ）
$ setrmt_lapw case_initial
$ cp -f case_initial.struct_setrmt case.struct
$ init_lapw -b -numk 1000 -rkmax 7.0
$ rm *.broyd*
$ run_lapw -ec 0.0001 -cc 0.001 -fc 0.5 -p -min
===== 例：TiO2（ルチル構造） =====
例としてルチル（正方晶系）でやってみました。
CIFファイルは https://staff.aist.go.jp/nomura-k/japanese/itscgallary.htm から入手できますので、ここからstructファイルを作成してください。
格子定数は A = 8.679512 bohr = 4.5932 Ang、C = 5.591700 bohr = 2.9592 Ang です。
プロジェクト名は Rutile としました。
- 内部構造の最適化
$ init_lapw -b -numk 1000 -rkmax 7.5
$ run_lapw -ec 0.0001 -cc 0.001 -fc 0.5 -p -NI -min
- マフィンチン半径を小さくして再設定
$ setrmt -r 4 Rutile
$ cp -f Rutile.struct_setrmt Rutile.struct
==== [2] VARY C/A RATIO with CONSTANT VOLUME (tetr and hex lattices) を選ぶ場合 ====
- 体積を保ちつつC/Aを変化させます。計算は-2%から+2%までの5つで行うことにしました。
$ x_lapw optimize
********************************************
GENERATES STRUCT-FILES AND optimize.job
PLEASE CHOOSE ONE OF THE FOLLOWING FEATURES:
[1] VARY VOLUME with CONSTANT RATIO A:B:C
[2] VARY C/A RATIO with CONSTANT VOLUME (tetr and hex lattices)
[3] VARY C/A RATIO with CONSTANT VOLUME and B/A (orthorh lattice)
[4] VARY B/A RATIO with CONSTANT VOLUME and C/A (orthorh lattice)
[5] VARY A and C (2D-case) (tetragonal or hexagonal lattice)
[6] VARY A, B and C (3D-case) (orthorhombic lattice)
[7] VARY A, B, C and Gamma (4D-case) (monoclinic lattice)
[8] VARY C/A RATIO and VOLUME (2D-case) (tetr and hex lattices)
********************************************
2
***************************************************
Generating Rutile_initial.struct
next time this file will be used as template unless you remove it explicitly.
***************************************************
NUMBER OF STRUCTURE CHANGES ?
5
PLEASE ENTER VALUE 1 (IN %)
-2
PLEASE ENTER VALUE 2 (IN %)
-1
PLEASE ENTER VALUE 3 (IN %)
0
PLEASE ENTER VALUE 4 (IN %)
1
PLEASE ENTER VALUE 5 (IN %)
2
Rutile_coa___-2.00.struct
8.738160 8.738160 5.516893 90.000000
Rutile_coa___-1.00.struct
8.708639 8.708639 5.554360 90.000000
Rutile_coa____0.00.struct
8.679513 8.679513 5.591700 90.000000
Rutile_coa____1.00.struct
8.650772 8.650772 5.628916 90.000000
Rutile_coa____2.00.struct
8.622409 8.622409 5.666010 90.000000
Now run optimize.job
0.007u 0.000s 0:11.49 0.0% 0+0k 0+88io 0pf+0w
$ vim optimize.job
（run_lapwで始まる行のコメントを外して次のようにします ... run_lapw -ec 0.0001 -p -it -cc 0.01 -fc 1 -min）
$ rm *.broyd*
$ ./optimize.job
* 計算にはめっちゃ時間かかります
- 結果を整理します
$ eplot_lapw -a "" -t coa
####################################
# #
# E - PLOT #
# Extension by Morteza Jamal #
# (2014) #
####################################
iter chisq delta/lim lambda a1 a2 a3 a4 a5
0 8.1119622037e+07 0.00e+00 5.26e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
1 3.4613698959e+07 -1.34e+05 5.26e-01 -3.624361e+02 7.155063e-01 -2.197764e+02 2.153182e-02 -1.766469e+02
2 4.4453663198e+05 -7.69e+06 5.26e-02 -3.366452e+03 1.402766e-01 -6.812561e+02 -3.653918e-02 1.290336e+02
3 2.1224532734e+01 -2.09e+09 5.26e-03 -4.013859e+03 4.997830e-04 -7.021329e+00 -1.067682e-04 1.497482e+00
4 1.7489236145e-07 -1.21e+13 5.26e-04 -4.018078e+03 1.261558e-04 -3.565157e-04 -8.273419e-06 1.303460e-04
5 1.7153761350e-19 -1.02e+17 5.26e-05 -4.018078e+03 1.261458e-04 2.897735e-04 -8.270833e-06 -8.678723e-06
6 2.0679515314e-25 -8.30e+10 5.26e-06 -4.018078e+03 1.261458e-04 2.897738e-04 -8.270833e-06 -8.678750e-06
7 2.0679515314e-25 0.00e+00 5.26e-07 -4.018078e+03 1.261458e-04 2.897738e-04 -8.270833e-06 -8.678750e-06
iter chisq delta/lim lambda a1 a2 a3 a4 a5
After 7 iterations the fit converged.
final sum of squares of residuals : 2.06795e-25
rel. change during last iteration : 0
Exactly as many data points as there are parameters.
In this degenerate case, all errors are zero by definition.
Final set of parameters
=======================
a1 = -4018.08
a2 = 0.000126146
a3 = 0.000289774
a4 = -8.27083e-06
a5 = -8.67875e-06
Do you want a hardcopy? (Y/n)y
Specify a filename (default is Rutile.eplot_coa_.ps)
Printing hardcopy
Rutile.eplot_coa_.ps generated
******************************
==========================================================
Lowest point on X-axis is = -0.216300
Minimum value of C/A is = 0.642848
Minimum value of ENE is = -4018.078221 Ry
==========================================================
Value of A is = 8.68578 bohr ; 4.59632 Ang
Value of B is = 8.68578 bohr ; 4.59632 Ang
Value of C is = 5.58363 bohr ; 2.95473 Ang
==========================================================
Initial value of C/A is = 0.644241
Conv. Unit Cell Vol = 421.2447 bohr^3 ; 62.4220 Ang^3
==========================================================
******************************
- フィッテングの結果、A = 8.68578 bohr = 4.59632 Ang、C = 5.58363 bohr = 2.95473 Ang となることがわかりました。
==== [5] VARY A and C (2D-case) (tetragonal or hexagonal lattice) を選ぶ場合 ====
- 3 × 3 = 9 のメッシュで格子定数AとCを変化させます。
$ x_lapw optimize
********************************************
GENERATES STRUCT-FILES AND optimize.job
PLEASE CHOOSE ONE OF THE FOLLOWING FEATURES:
[1] VARY VOLUME with CONSTANT RATIO A:B:C
[2] VARY C/A RATIO with CONSTANT VOLUME (tetr and hex lattices)
[3] VARY C/A RATIO with CONSTANT VOLUME and B/A (orthorh lattice)
[4] VARY B/A RATIO with CONSTANT VOLUME and C/A (orthorh lattice)
[5] VARY A and C (2D-case) (tetragonal or hexagonal lattice)
[6] VARY A, B and C (3D-case) (orthorhombic lattice)
[7] VARY A, B, C and Gamma (4D-case) (monoclinic lattice)
[8] VARY C/A RATIO and VOLUME (2D-case) (tetr and hex lattices)
********************************************
5
***************************************************
Generating Rutile_initial.struct
next time this file will be used as template unless you remove it explicitly.
***************************************************
number of structures: 6, 9 (3x3), 16 (4x4), 25 (5x5), 36
9
PLEASE enter a percentage change of a
2
Rutile_a+c____1.00.struct
8.505922 8.505922 5.479866 90.000000
Rutile_a+c____2.00.struct
8.679512 8.679512 5.479866 90.000000
Rutile_a+c____3.00.struct
8.853102 8.853102 5.479866 90.000000
Rutile_a+c____4.00.struct
8.505922 8.505922 5.591700 90.000000
Rutile_a+c____5.00.struct
8.679512 8.679512 5.591700 90.000000
Rutile_a+c____6.00.struct
8.853102 8.853102 5.591700 90.000000
Rutile_a+c____7.00.struct
8.505922 8.505922 5.703534 90.000000
Rutile_a+c____8.00.struct
8.679512 8.679512 5.703534 90.000000
Rutile_a+c____9.00.struct
8.853102 8.853102 5.703534 90.000000
Now run optimize.job
0.000u 0.007s 0:10.77 0.0% 0+0k 0+168io 0pf+0w
$ vim optimize.job
（run_lapwで始まる行のコメントを外して次のようにします ... run_lapw -ec 0.0001 -p -it -cc 0.01 -fc 1 -min）
$ rm *.broyd*
$ ./optimize.job
- フィッテイングします
$ parabolfit_lapw -t 2
The following scf files were used for analysis:
Rutile_a+c____1.00_default.scf
Rutile_a+c____2.00_default.scf
Rutile_a+c____3.00_default.scf
Rutile_a+c____4.00_default.scf
Rutile_a+c____5.00_default.scf
Rutile_a+c____6.00_default.scf
Rutile_a+c____7.00_default.scf
Rutile_a+c____8.00_default.scf
Rutile_a+c____9.00_default.scf
Rutile.ene and Rutile.latparam generated
Enter dimension of fit (number of variable lattice parameters, 1-6):
2 fitcase 6 parameter
lowest data point: -4018.08034234000 8.85310000000000
5.59170000000000
9 -4018.05375612000 -4018.07178379000
-4018.07753850000 -4018.06454837000 -4018.07809085000
-4018.08034234000 -4018.06932397000 -4018.07887390000
-4018.07702213000
I INITIAL X(I) D(I)
1 -0.401808D+04 0.100D+01
2 0.100000D+00 0.600D+00
3 0.885310D+01 0.600D+00
4 0.100000D+00 0.600D+00
5 0.559170D+01 0.600D+00
6 0.100000D+00 0.600D+00
IT NF F RELDF PRELDF RELDX MODEL STPPAR D*STEP NPRELDF
0 1 0.592D-04
1 2 0.370D-04 0.38D+00 0.10D+01 0.1D-04 G 0.0D+00 0.1D+00 0.10D+01
2 3 0.192D-05 0.95D+00 0.10D+01 0.3D-05 G 0.0D+00 0.3D-01 0.10D+01
3 4 0.213D-06 0.89D+00 0.11D+01 0.3D-06 S 0.0D+00 0.3D-02 0.11D+01
4 5 0.487D-07 0.77D+00 0.77D+00 0.2D-06 G 0.0D+00 0.2D-02 0.77D+00
5 6 0.487D-07 0.27D-03 0.27D-03 0.8D-09 G 0.0D+00 0.9D-05 0.27D-03
***** X-CONVERGENCE *****
FUNCTION 0.486922D-07 RELDX 0.817D-09
FUNC. EVALS 6 GRAD. EVALS 36
PRELDF 0.267D-03 NPRELDF 0.267D-03
I FINAL X(I) D(I) G(I)
1 -0.401808D+04 0.100D+01 0.534D-08
2 0.193393D+00 0.121D+00 0.172D-09
3 0.878393D+01 0.138D+00 0.246D-09
4 0.235417D+00 0.467D-01 0.476D-10
5 0.561585D+01 0.122D+00 0.178D-09
6 0.207137D+00 0.467D-01 0.137D-10
Parabolic equation of state: info 3
E = x1 + x2(a-x3)^2
+ x4(b-x5)^2 + x6(a-x3)(b-x5)
Fitparameter are
-4018.080980 0.193393 8.783926 0.235417
5.615850 0.207137
lattic parameters energy de(EOS)
8.505920 5.479870 -4018.053756 -0.000094
8.679510 5.479870 -4018.071784 0.000206
8.853100 5.479870 -4018.077539 -0.000112
8.505920 5.591700 -4018.064548 0.000043
8.679510 5.591700 -4018.078091 -0.000121
8.853100 5.591700 -4018.080342 0.000079
8.505920 5.703530 -4018.069324 0.000051
8.679510 5.703530 -4018.078874 -0.000085
8.853100 5.703530 -4018.077022 0.000033
Sigma: 0.000110
Optionally create data points from fit function
Enter number of datapoints for your 2 dimensional Energy surface
NI=0 terminates; NI=1 will use 1 specific value in I-th component and allows to
generate 2D-cuts
0.016u 0.020s 0:00.03 100.0% 0+0k 0+8io 0pf+0w
* フィッティングの結果、A = 8.783926 bohr, C = 5.615850 bohr となりました。
==== で、どっちが正しいの？ ====
この質問はあまり意味がないかもしれません。
なぜなら自然科学というのは「実験で得られたデータが正しい」ものであり、我々は単に密度汎関数理論の枠組みでPBE汎関数を使ったときのエネルギー的に安定な格子構造を求めたにすぎないからです。
実験で確実にわかっている量をもとにて、不確実な部分だけ第一原理で決める、というのが基本的な運用方法です。