CO分子

以下の例では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) 計算

出力ファイルを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
               @@                          |___________________|
                ***** ...                [_______________________]

が表示される。 収束が難しい場合は一旦休憩し、ミキシングパラメーター等を調整し、計算を再開する。

構造最適化(マニュアル構造最適化)

ここでは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を得る。
    etot_CO_fit.png

構造最適化

以下では入力ファイルを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であり、フィッティングによって得られたものと良い一致を示す。

構造最適化の各ステップでの原子位置が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で計算を続ける方がかえって効率が良い場合もある。

添付ファイル: fileetot_CO_fit.pdf 276件 [詳細] fileetot_CO.pdf 295件 [詳細]
トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-01-23 (月) 11:56:58 (461d)