STATE_tutorial/H2_jp - 第一原理分子動力学プログラム STATE Senri Wiki

水素分子 †

このチュートリアルでは分子系について平面波カットオフとユニットセルの大きさについての収束性を調べる手順を説明します。周期境界条件の下で孤立分子のシミュレーションを行う場合、分子は周期イメージとの相互作用が無視できる程度に大きな箱（ユニットセル）を用います。ここでは水素分子を取り上げ、局所密度近似を用いて計算を行います。使用するポテンシャルはpot.H_lda1です。

↑

SCF 計算 †

入力ファイル (nfinp_1)

      0     0     0     0     0     0       : H2 molecule in a box 
    5.0  10.0     1     2     2             : GMAX GMAXP NTYP NATM NATM2
      1     0                               : NUM_SPACE_GROUP TYPE_OF_BRAVIS_LATTICE 
   10.0  10.0  10.0  90.0  90.0  90.0       : A B C ALPHA BETA GAMMA
      1     1     1     1     1     1       : N1 N2 N3 M1 M2 M3 (k-point mesh and shift)        
      1     0                               : NCORD NINV
 -0.710 0.000 0.000     1     1     1       : CPS(1:3) IWEI IMDTYP ITYP
  0.710 0.000 0.000     1     1     1       : CPS(1:3) IWEI IMDTYP ITYP
      1  0.15     1.00794     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                               : MIX_WAY MIX_WHAT
      0     8   0.7                         : STARTING_MIXING KBXMIX ALPHA
   0.60  0.50  0.60  0.70  1.00             : DTIM1 DTIM2 DTIM3 DTIM4 DTIM_LAST
  100.0     4     1    0.10D-08             : DTIO IMDALG IEXPL EDELTA
 0.0010    0.50D+03     0                   : WIDTH FORCCR ISTRESS
ldapw91           1                         : XCTYPE NSPIN
   1.00                                     : DESTM       
    101                                     : NBZTYP
      0     0     0                         : DUMMY         
      0     0     0                         : DUMMY         
      2                                     : NEG
      1                                     : NEXTST
      0                                     : DUMMY
      2                                     : IMSD
      0                                     : EVALUATE_EKO_DIFF
      0                                     : NPDOSAO
      0    0.0                              : SM_N dopping

実行

$ mpirun -np 4 ./STATE < nfinp_1 > nfout_1

結果
STATEを実行すると以下のメッセージが表示されます。

 ***********************************************************************
 *                                                                     *
 *                                                                     *
 *                                                                     *
 *              ******  ********    **    ******** ********            *
 *             ******** ********   ****   ******** ********            *
 *             **          **     **  **     **    **                  *
 *              ***        **    ********    **    ******              *
 *                ***      **   **********   **    ******              *
 *                  **     **  **        **  **    **                  *
 *             ********    ** **          ** **    ********            *
 *              ******     ** VERSION 5.6.5  **    ********            *
 *                               RICS-AIST                             *
 *                           OSAKA UNIVERSITY                          *
 *                                                                     *
 ***********************************************************************

SCF計算が収束すると以下の全エネルギーおよび運動エネルギーなどの各寄与が表示されます。

                     TOTAL ENERGY AND ITS COMPONENTS 
                  TOTAL ENERGY     =          -1.13758782 A.U.
                   FREE ENERGY     =          -1.13758782 A.U.
                KINETIC ENERGY     =           0.96824571 A.U.
                HARTREE ENERGY     =           0.73570075 A.U.
                     XC ENERGY     =          -0.64555732 A.U.
                  LOCAL ENERGY     =          -2.31372346 A.U.
               NONLOCAL ENERGY     =          -0.02336940 A.U.
                  EWALD ENERGY     =           0.14111590 A.U.
                     PC ENERGY     =           0.00000000 A.U.
               ENTROPIC ENERGY     =           0.00000000 A.U.

全エネルギー以降にフェルミエネルギー、Kohn-Sham固有値と占有数が表示されます。

                                         FERMI ENERGY =      -0.19018506
 NKP=    1  NGP=      2103  K=(  0.00000  0.00000  0.00000)  WKP= 1.0000
                              EIGEN VALUE 
 -0.37074 -0.00963
                               OCCUPATION 
  1.00000  0.00000

さらに原子にかかる力が以下のように表示されます。

    ATOM              COORDINATES                        FORCES
MD:   1
MD:     1  H   -0.710000    0.000000    0.000000  -0.01200 -0.00000  0.00000
MD:     2  H    0.710000    0.000000    0.000000   0.01200 -0.00000  0.00000

↑

収束性の検証 †

平面波基底を用いた計算の場合、波動関数の平面波展開のためのカットオフエネルギーが計算の精度を決定する最も重要なパラメーターの一つとなります。ノルム保存擬ポテンシャルを用い、なおかつ非線形内殻補正を考慮しない場合、電荷密度の平面波カットオフエネルギーは波動関数のそれの4倍になります。しかしながらウルトラソフト擬ポテンシャルの場合は局在した補強電荷を正確に記述するためにしばしば波動関数についての4倍以上のカットオフエネルギーが必要となる場合があります。 STATEでは波動関数を展開する平面波のカットオフ波数ベクトル (GMAX) を設定しています。GMAXの二乗がカットオフエネルギー（Ry）に対応します。電荷密度のカットオフ波数ベクトルはGMAXPで指定され、GMAXPの二乗が電荷密度のカットオフエネルギーになります。

GMAXとGMAXPを変えながら全エネルギー計算を実行することで以下のデータが得られます。 [GMAXP=15 (E_cut^den=225 Ry)]

#GMAX ETOT (Ha)
4.0 -1.12781490
4.5 -1.13418820
5.0 -1.13697688
5.5 -1.13816998
6.0 -1.13847300
6.5 -1.13850127
7.0 -1.13851867
7.5 -1.13857102

[GMAXP=20 (E_cut^den=400 Ry)]

#GMAX ETOT (Ha)
4.0 -1.12780183
4.5 -1.13417130
5.0 -1.13695809
5.5 -1.13815036
6.0 -1.13845313
6.5 -1.13848136
7.0 -1.13849876
7.5 -1.13855109

[GMAXP=25 (E_cut^den=625 Ry)]

#GMAX ETOT (Ha)
4.0 -1.12780130
4.5 -1.13417074
5.0 -1.13695751
5.5 -1.13814979
6.0 -1.13845255
6.5 -1.13848078
7.0 -1.13849818
7.5 -1.13855052

[GMAXP=30 (E_cut^den=900 Ry)]

#GMAX ETOT (Ha)
4.0 -1.12780182
4.5 -1.13417124
5.0 -1.13695801
5.5 -1.13815027
6.0 -1.13845303
6.5 -1.13848127
7.0 -1.13849867
7.5 -1.13855100

GMAX=6とGMAXP=20で全エネルギーは0.1 mRyの精度で収束していることが分かります（エネルギーの変化が0.1 mRy以下）。

孤立分子のシミュレーションを実行する場合、カットオフエネルギーに加えてスーパーセルのサイズについての収束性を検証する必要もあります。 GMAX=6とGMAXP=20を用いて全エネルギーをユニットセル（立方体）のサイズ（辺の長さ）の関数として計算すると以下が得られます。

#A(Bohr) ETOT(Ha)
 5.00 -1.31378028
10.00 -1.13845313
15.00 -1.13662806
20.00 -1.13661187

上の結果から立方体のユニットセルの辺の長さが15 Bohrとなったところでエネルギーが変化が0.1 mRy以下と小さくなり、収束しているといえることが分かります。大きい分子の場合はより大きいユニットセルを用いる必要があります。

↑

構造最適化 (1) †

二原子分子の場合、原子間距離（ボンド長）を変えながら全エネルギーを計算することで最適な構造（エネルギーが最小となる構造）を求めることができます。ここではGMAX=6、GMAXP=20、辺の長さが10 Bohrの立方体のユニットセルを用い、全エネルギーをボンド長の関数として計算を行います。

# Bond length (Bohr) Total energy (Hartree)
1.28000000 -1.13291015
1.30000000 -1.13434083
1.32000000 -1.13555229
1.34000000 -1.13655937
1.36000000 -1.13737569
1.38000000 -1.13801370
1.40000000 -1.13848494
1.42000000 -1.13880026
1.44000000 -1.13896985
1.46000000 -1.13900337
1.48000000 -1.13890994
1.50000000 -1.13869810
1.52000000 -1.13837576
1.54000000 -1.13795024
1.56000000 -1.13742831
1.58000000 -1.13681622
1.60000000 -1.13611994
1.62000000 -1.13534513
1.64000000 -1.13449712
1.66000000 -1.13358088

得られた全エネルギーを6次の多項式にフィットすることで平衡ボンド長 1.4551 Bohr (0.77 Angstrom)が得られます。

↑

構造最適化 (2) †

より複雑な系では手動で上記のような構造最適化を行うのは大変困難です。そのような場合にはHellmann-Feynman力を求めて構造最適化を行うのが一般的です。 Hellmann-Feynmann力を用いて構造最適化を行う場合、力の閾値（FORCCR）を1.e-3 (1.d-3)から1.e-4 (1.e-4) Hartree/Bohrに設定します（構造最適化を行わない場合FORCCRは1.d+3などの大きい値に設定しています）。

この例では最適化のアルゴリムとしてgeneralized direct inversion of iterative subspace (GDIIS) 法 (IMDALG=4) を採用し、原子位置を更新する時間刻み (DTIO)として100 atomic unitを採用しています。原子の質量 (AMION) は1.00794 a.m.u.としていますが、構造最適化の場合はより重い質量を使う方が計算が安定する場合が多いです。

      0     0     0     0     0     0       : H2 molecule in a box 
    6.0  20.0     1     2     2             : GMAX GMAXP NTYP NATM NATM2
      1     0                               : NUM_SPACE_GROUP TYPE_OF_BRAVIS_LATTICE 
   10.0  10.0  10.0  90.0  90.0  90.0       : A B C ALPHA BETA GAMMA
      1     1     1     1     1     1       : N1 N2 N3 M1 M2 M3 (k-point mesh and shift) 
      1     0                               : NCORD NINV
-0.710  0.000 0.000     1     1     1       : CPS(1:3) IWEI IMDTYP ITYP
 0.710  0.000 0.000     1     1     1       : CPS(1:3) IWEI IMDTYP ITYP
      1  0.15     1.00794     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                               : MIX_WAY MIX_WHAT
      0     8   0.7                         : STARTING_MIXING KBXMIX ALPHA
   0.60  0.50  0.60  0.70  1.00             : DTIM1 DTIM2 DTIM3 DTIM4 DTIM_LAST
  100.0     4     1    0.10D-08             : DTIO IMDALG IEXPL EDELTA
 0.0010    0.50D-03     0                   : WIDTH FORCCR ISTRESS
ldapw91           1                         : XCTYPE NSPIN
   1.00                                     : DESTM 
    101                                     : NBZTYP
      0     0     0                         : DUMMY 
      0     0     0                         : DUMMY 
      2                                     : NEG
      1                                     : NEXTST
      0                                     : DUMMY
      2                                     : IMSD
      0                                     : EVALUATE_EKO_DIFF
      0                                     : NPDOSAO
      0    0.0                              : SM_N dopping

得られた平衡ボンド長は1.455 Bohr (0.77 Angstrom)で手動最適化によって得られた値と良い一致を示していることが分かります。

↑

演習問題 †

ボンド長のユニットセルサイズおよびカットオフ依存生を調べてみましょう。

第一原理分子動力学プログラム STATE Senri Wiki