*CO分子 [#u76be93a]
以下の例ではrectangular boxに配置したCO分子の計算例を示す。
ソースおよび擬ポテンシャルは以下を使用すると仮定する。
- ソースディレクトリ: ${HOME}/STATE/src/state-5.6.3/src
- 擬ポテンシャルディレクトリ: ${HOME}/STATE/gncpp
- 擬ポテンシャル: pot_C_pbe1, pot_O_pbe1
-入力ファイル (nfinp_1)
0 0 0 0 0 0 : dummy line (6 integers)
5.50 20.00 2 2 2 : GMAX, GMAXP, NTYP, NATM, NATM2
1 0 : space group number, bravis lattice type
6.00 4.00 4.00 90.00 90.00 90.00 : a, b, c, alpha, beta, gamma
1 1 1 1 1 1 : knx, kny, knz, k-point shift
1 0 : NCORD, NINV
0.0000 0.0000 0.0000 1 1 1 : cps, iwei, imdtyp, ityp
2.2000 0.0000 0.0000 1 1 2 : cps, iwei, imdtyp, ityp
6 0.1500 51577.50 3 1 0.d0 : IATOMN, ALFA, AMION, ILOC, IVAN, ZETA1
8 0.1500 51577.50 3 1 0.d0 : IATOMN, ALFA, AMION, ILOC, IVAN, ZETA1
0 0 0 0 0 : ICOND, INIPOS, INIVEL, ININOSE, INIACC
0 1 : IPRE, IPRI
200 200 0 57200.00 0 : NMD1, NMD2, iter_last, CPUMAX, ifstop
3 1 : way_mix, mix_what
0 8 0.8 : starting mixing, kbxmix,alpha
0.60 0.50 0.60 0.70 1.00 : DTIM1, DTIM2, DTIM3, DTIM4, dtim_last
30.00 2 1 0.10D-08 1.d-06 : DTIO, IMDALG, IEXPL, EDELTA
0.0010 0.10D+02 0 : WIDTH, FORCCR, ISTRESS
ggapbe 1 : XCTYPE, nspin
1.00 3 : destm, n_stm
102 : NBZTYP
0 0 0 : NKX, NKY, NKZ (dummy)
0 0 0 : NKX2, NKY2, NKZ2 (dummy)
8 : NEG (# of bands)
1 : NEXTST (1: G-space, 0: R-space)
0 : 0; random numbers, 1; matrix diagon
2 : imsd (2: Davidson, 1: RMM)
0 : eval. eko diff.: .0 = no ,1 = yes
0 : npdosao
0 0.0 : SM_dopping
以下を実行しSTATEおよび擬ポテンシャルへのシンボリックリンクを作成
ln -s ${HOME}/STATE/src/state-5.6.3/src/STATE STATE
ln -s ${HOME}/STATE/gncpp/pot_C_pbe1 fort.37
ln -s ${HOME}/STATE/gncpp/pot_O_pbe1 fort.38
あるいは
ln -s ${HOME}/STATE/gncpp/C_pbe1/#vnew.data fort.37
ln -s ${HOME}/STATE/gncpp/O_pbe1/#vnew.data fort.38
** Single point (SCF) 計算 [#p8468442]
出力ファイルをnfout_1とし以下を実行する (コマンドやジョブスクリプトは環境に応じて書き換えること)
mpirun -np 2 ./STATE < nfinp_1 > nfout_1
ジョブが開始すると以下のロゴがnfout_1に表示される
***********************************************************************
* *
* *
* *
* ****** ******** ** ******** ******** *
* ******** ******** **** ******** ******** *
* ** ** ** ** ** ** *
* *** ** ******** ** ****** *
* *** ** ********** ** ****** *
* ** ** ** ** ** ** *
* ******** ** ** ** ** ******** *
* ****** ** VERSION 5.6.5 ** ******** *
* RICS-AIST *
* OSAKA UNIVERSITY *
* *
***********************************************************************
SCF計算が始まるところでは以下が出力される。
***********************************************************************
* *
* START SCF *
* *
***********************************************************************
全エネルギーの収束をモニターする場合は
grep ETOT\: nfout_1
を実行する。その結果は以下のようになる
ETOT: 1 -16.71058056 0.1671E+02 0.8965E-01
ETOT: 2 -20.04069483 0.3330E+01 0.6387E-01
ETOT: 3 -21.96017776 0.1919E+01 0.4847E-01
ETOT: 4 -22.11633389 0.1562E+00 0.3198E-01
ETOT: 5 -22.20286500 0.8653E-01 0.1510E-01
ETOT: 6 -22.21912414 0.1626E-01 0.3085E-02
ETOT: 7 -22.21938566 0.2615E-03 0.7750E-03
ETOT: 8 -22.21941988 0.3422E-04 0.2094E-03
ETOT: 9 -22.21942413 0.4249E-05 0.4735E-04
ETOT: 10 -22.21942395 0.1857E-06 0.4811E-04
ETOT: 11 -22.21942422 0.2798E-06 0.1838E-04
ETOT: 12 -22.21942425 0.2761E-07 0.6088E-05
ETOT: 13 -22.21942426 0.3338E-08 0.3279E-06
ETOT: 14 -22.21942426 0.8036E-11 0.8071E-07
ETOT: 15 -22.21942426 0.1084E-11 0.1565E-07
ETOT: 16 -22.21942426 0.3197E-13 0.7047E-08
SCF計算が収束した場合、全エネルギーとその成分が以下のように表示される。
TOTAL ENERGY AND ITS COMPONENTS
TOTAL ENERGY = -22.21942426 A.U.
FREE ENERGY = -22.21942426 A.U.
KINETIC ENERGY = 9.92111448 A.U.
HARTREE ENERGY = 5.12121891 A.U.
XC ENERGY = -5.89585656 A.U.
LOCAL ENERGY = -20.23161767 A.U.
NONLOCAL ENERGY = 6.73686187 A.U.
EWALD ENERGY = -17.87114528 A.U.
PC ENERGY = 0.00000000 A.U.
ENTROPIC ENERGY = 0.00000000 A.U.
原子にかかる力は以下のように表示される。
ATOM COORDINATES FORCES
MD: 1
MD: 1 C 0.000000 0.000000 0.000000 0.01852 0.00000 -0.00000
MD: 2 O 2.200000 0.000000 0.000000 -0.01858 -0.00000 0.00000
また一般にSCFあるいは構造最適化が収束した場合には以下の"victory cat"が表示される。
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
_______________________
__________ _______/______v______v______v___]
D | | |
D A A | | Congratulations! | C( > < )D
-- =(^.^)= | | The calculation has converged. | = o =
| @@@@@ | | | ( )~
/--=O=-+-=O=---+--=O=--+--==O==--+--==O==--+--=O=-+--=O=---=O=-/
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
一方、収束しなかった場合には
Sorry! < < <
The calculation has not converged. < < <
< < <
___________________
@ @ | |
* *** | |XXX
* * * * | Have a break! | X
*** *** ... | | X
| | X
| |XXX
@@ |___________________|
***** ... [_______________________]
が表示される。
収束が難しい場合は一旦休憩し、ミキシングパラメーター等を調整し、計算を再開する。
** 構造最適化（マニュアル構造最適化） [#pdaa4b56]
ここではCO分子の構造をマニュアルで最適化する。以下の入力ファイルを用意する。
- nfinp_d2.12 nfinp_d2.14 nfinp_d2.16 nfinp_d2.18 nfinp_d2.20 nfinp_d2.22
これらの入力ファイルは上記のnfinp_1と基本的に同じであるが、原子位置（CO間距離）のみ異なる。
- nfinp_d2.12
0.0000 0.0000 0.0000 1 1 1 : cps, iwei, imdtyp, ityp
2.1200 0.0000 0.0000 1 1 2 : cps, iwei, imdtyp, ityp
- nfinp_d2.14
0.0000 0.0000 0.0000 1 1 1 : cps, iwei, imdtyp, ityp
2.1400 0.0000 0.0000 1 1 2 : cps, iwei, imdtyp, ityp
...~
- nfinp_d2.22
0.0000 0.0000 0.0000 1 1 1 : cps, iwei, imdtyp, ityp
2.2200 0.0000 0.0000 1 1 2 : cps, iwei, imdtyp, ityp
これらの入力ファイルを用いて順次計算を行う。
mpirun -np 2 ./STATE < nfinp_d2.12 > nfout_d2.12
mpirun -np 2 ./STATE < nfinp_d2.14 > nfout_d2.14
...
mpirun -np 2 ./STATE < nfinp_d2.22 > nfout_d2.22
すべての計算が終了したならば
grep 'TOTAL ENERGY' nfout_d2.* | grep A.U.| awk '{print $5}'
を実行すると
-22.21807672
-22.21900649
-22.21951476
-22.21964087
-22.21942426
-22.21890297
が得られる。これを編集して
#d (Bohr) Etot (Hartree)
2.12 -22.21807672
2.14 -22.21900649
2.16 -22.21951476
2.18 -22.21964087
2.20 -22.21942426
2.22 -22.21890297
を得る。
これを3次関数でフィットすると平衡距離として2.177 Bohrを得る。
#ref(http://www-cp.prec.eng.osaka-u.ac.jp/puki_state/graph/etot_CO_fit.png,center,nowrap,nolink)
** 構造最適化 [#wd5642ab]
以下では入力ファイルをnfinp_2としforccr=1.0D-03 (1.0^-3 Hartree/BOhr)とすることで構造最適化を実行する。ここで構造最適化のアルゴリズムとしてquenched MD (IMDALG=2)、タイムステップ (DTIO)を50 a.u.とした。
- 入力ファイル (nfinp_2)
0 0 0 0 0 0 : I_CTRL(1:6)
5.50 20.00 2 2 2 : GMAX GMAXP NTYP NATM NATM2
1 0 : number of space group, type of bravis lattice
6.00 4.00 4.00 90.00 90.00 90.00 : A B C ALPHA BETA GAMMA
1 1 1 1 1 1 : N1 N2 N3 M1 M2 M3
1 0 : NCORD, NINV
0.00 0.00 0.00 1 1 1 : CPS(1,1:3) IWEI IMDTYP ITYP
2.20 0.00 0.00 1 1 2 : CPS(1,1:3) IWEI IMDTYP ITYP
6 0.15 12.01 3 1 0.d0 : IATOMN ALFA AMION ILOC IVAN
8 0.15 16.00 3 1 0.d0 : IATOMN ALFA AMION ILOC IVAN
0 0 0 0 0 : ICOND INIPOS INIVEL ININOS INIACC
0 1 : IPRE IPRI
200 200 0 57200.00 0 : NMD1 NMD2 ITER_LAST CPUMAX IFSTOP
3 1 : WAY_MIX MIX_WHAT
0 8 0.8 : STARTING_MIXING KBXMIX ALPHA
0.60 0.50 0.60 0.70 1.00 : DTIM1 DTIM2 DTIM3 DTIM4 DTIM_LAST
50.00 2 1 0.10D-08 : DTIO IMDALG IEXPL EDELTA
0.0010 1.0D-03 0 : WIDTH FORCCR ISTRESS
ggapbe 1 : XCTYPE NSPIN
1.00 : DESTM
102 : NBZTYP
0 0 0 : NKX, NKY, NKZ
0 0 0 : NKX2,NKY2,NKZ2
8 : NEG
1 : NEXTST
0 : (DUMMY)
2 : IMSD
0 : EVALUATE_EKO_DIFF
0 : NPDOSAO
0 0.0 : SM_N DOPPING
以下を実行。
mpirun -np ./STATE < nfinp_2 > nfout_2
以下を実行すると構造最適化における各ステップでの全エネルギーが得られる
grep 'TOTAL ENERGY' nfout_opt | grep A.U.
出力結果は以下のようになる。
TOTAL ENERGY = -22.21942426 A.U.
TOTAL ENERGY = -22.21948810 A.U.
TOTAL ENERGY = -22.21957767 A.U.
TOTAL ENERGY = -22.21963740 A.U.
TOTAL ENERGY = -22.21962705 A.U.
TOTAL ENERGY = -22.21962705 A.U.
TOTAL ENERGY = -22.21963215 A.U.
TOTAL ENERGY = -22.21963885 A.U.
TOTAL ENERGY = -22.21964257 A.U.
また以下を実行すると各ステップでの原子にかかる力を取り出すことができる。
grep MD\: nfout_opt
出力結果は以下の通り。
MD: 1
MD: 1 C 0.000000 0.000000 0.000000 0.01852 0.00000 -0.00000
MD: 2 O 2.200000 0.000000 0.000000 -0.01858 -0.00000 0.00000
MD: 2
MD: 1 C 0.002115 0.000000 -0.000000 0.01567 -0.00000 0.00000
MD: 2 O 2.198408 -0.000000 0.000000 -0.01572 0.00000 -0.00000
MD: 3
MD: 1 C 0.006018 -0.000000 -0.000000 0.01025 -0.00000 -0.00000
MD: 2 O 2.195467 -0.000000 0.000000 -0.01030 0.00000 -0.00000
MD: 4
MD: 1 C 0.011092 -0.000000 -0.000000 0.00289 -0.00000 0.00000
MD: 2 O 2.191645 -0.000000 0.000000 -0.00293 -0.00000 -0.00000
MD: 5
MD: 1 C 0.016496 -0.000000 -0.000000 -0.00535 -0.00000 0.00000
MD: 2 O 2.187570 -0.000000 0.000000 0.00531 -0.00000 -0.00000
MD: 6
MD: 1 C 0.016496 -0.000000 -0.000000 -0.00535 0.00000 0.00000
MD: 2 O 2.187570 -0.000000 0.000000 0.00532 -0.00000 0.00000
MD: 7
MD: 1 C 0.015885 -0.000000 -0.000000 -0.00440 0.00000 0.00000
MD: 2 O 2.188026 -0.000000 0.000000 0.00437 -0.00000 -0.00000
MD: 8
MD: 1 C 0.014771 -0.000000 -0.000000 -0.00269 -0.00000 0.00000
MD: 2 O 2.188856 -0.000000 0.000000 0.00265 0.00000 0.00000
MD: 9
MD: 1 C 0.013351 -0.000000 0.000000 -0.00053 -0.00000 0.00000
MD: 2 O 2.189913 -0.000000 0.000000 0.00049 0.00000 0.00000
Hellmann-Feynman力を用いて最適化したCO距離は2.177 Bohrであり、フィッティングによって得られたものと良い一致を示す。
- 構造最適化の際のヒント~
構造最適化を初めて行う際は原子の質量を炭素あるいは珪素のものにし、DTIOを100-200程度、IMDALG=2とする。その後原子にかかる力の最大値（f_max）が10^-2程度になるとDTIO=400程度、IMDALG=4として計算を再開すると効率良く構造最適化ができることが多い。ただし、ポテンシャル面が浅く（相互作用が弱い分子ー表面系など）収束が遅い場合はIMDALG=2で計算を続ける方がかえって効率が良い場合もある。
構造最適化の各ステップでの原子位置がrestart.dataとGEOMETRYファイルに出力される。
restart.dataはバイナリー、GEOMETRYはasciiファイルであり、構造解析にはGEOMETRYファイルを利用するのが便利である。
- GEOMETRYファイル
&SYSTEM
1
2 2 2
6 8
1 1
6.00000 4.00000 4.00000 90.0000 90.0000 90.0000
&END
&COORDINATES
0.013350814786 -0.000000027561 0.000000006935
2.189913192072 -0.000000052830 0.000000019228
&END
&VELOCITIES
-0.000028411950 0.000000000055 0.000000000229
0.000021143930 -0.000000000198 0.000000000143
&END
&FORCES
-0.002688036207 -0.000000009984 0.000000054207
0.002651144409 0.000000007139 0.000000079038
&END
&SYSTEM ... &END
- 1行目: num_space_group (空間群番号)
- 2行目: ktyp, katm, katm2 (原子種の数、反転対称性について非等価な原子数、全原子数)
- 3行目: iatomn(1:ktyp) (原子番号)
- 4行目: iatom(1:ktyp) (各原子種毎の原子数)
- 5行目: a, b, c, alpha, beta, gamma
&COORDINATES ... &END~
カーテシアン座標での原子位置~
&VELOCITIES ... &END~
カーテシアン座標での各原子の速度~
&FORCES ... &END~
カーテシアン座標での原子にかかる力~
- 構造最適化の際のヒント~
構造最適化を初めて行う際は原子の質量を炭素あるいは珪素のものにし、DTIOを100-200程度、IMDALG=2とする。その後原子にかかる力の最大値（f_max）が10^-2程度になるとDTIO=400程度、IMDALG=4として計算を再開すると効率良く構造最適化ができることが多い。ただし、ポテンシャル面が浅く（相互作用が弱い分子ー表面系など）収束が遅い場合はIMDALG=2で計算を続ける方がかえって効率が良い場合もある。
![[PukiWiki]](image/logo_state-senri3_80x80.png "[PukiWiki]"){kind=logo}
