Bader電荷解析

Henkelman Groupが配布しているプログラムを利用してBader電荷の計算を行う。

グラフェン 1x1

gra1x1.png

以下ではノルム保存擬ポテンシャルC_pbe6TMを用いる。

SCF計算

入力ファイル

0 0 0 0 0 0
8.0000  20.0000   1   2   2   : GMAX, GMAXP, NTYP, NATM, NATM2
191  0  : number of space group, type of bravis lattice
4.655538364 4.655538364 20.000  90.0   90.0  120.0  : a,b,c,alpha,beta,gamma
24    24       1       1      1     1     : knx,kny,knz, k-point shift
0   0   : NCORD, NINV,   : IWEI, IMDTYP, ITYP
0.6666666666      0.3333333333      0.0000000000  1   1  1
0.3333333333      0.6666666666      0.0000000000  1   1  1
6  0.1500  12.0107 1 1 0.d0  : TYPE 1IATOMN,ALFA,AMION,ILOC,IVAN
0 0 0 0 0 : ICOND 0-MD, 1-CONT.MD, 2-WAVE FN,, 3-WAVE FN CONT.,  iconstpw
    0   1   : IPRE, IPRI
    200   200     0    3600.00  0 : NMD1, NMD2, iter_last, CPUMAX,ifstop
    3   1  : Simple=1,Broyd2=3,Blugel=6,    1:charge, 2:potential mix.
    0   20  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
 101  : NBZTYP 0-SF, 1-BK, 2-SC, 3-BCC, 4-FCC, 5-DIA, 6-HEX
    0   0   0   : NKX, NKY, NKZ
    0   0   0   : NKX2,NKY2,NKZ2
    8           : NEG
        1       : NEXTST(MB)
        0       : 0; random numbers, 1; matrix diagon
        2      0      0      0(MB) : imsd, i_2lm, i_sd2another, wksz for phase
        0       : evaluation of eko difference.0 = no ,1 = yes
        0                       : npdosao
        0    0.0     : SM_N, DOPPING

計算結果

$ grep ETOT your_output_file_name|grep "[0-9]"
ETOT:   1     15.49339609  0.1549E+02  0.4996E-01
ETOT:   2      0.34620940  0.1515E+02  0.1209E-01
ETOT:   3     -5.83837125  0.6185E+01  0.1717E-01
ETOT:   4     -8.60080202  0.2762E+01  0.1845E-01
ETOT:   5    -10.11782100  0.1517E+01  0.2536E-01
ETOT:   6    -11.03537609  0.9176E+00  0.2397E-01
ETOT:   7    -11.40509046  0.3697E+00  0.2251E-01
ETOT:   8    -11.45549051  0.5040E-01  0.2050E-01
ETOT:   9    -11.48179862  0.2631E-01  0.1827E-01
ETOT:  10    -11.47035224  0.1145E-01  0.1958E-01
ETOT:  11    -11.53852896  0.6818E-01  0.1092E-01
ETOT:  12    -11.57400237  0.3547E-01  0.4760E-02
ETOT:  13    -11.58514681  0.1114E-01  0.1420E-02
ETOT:  14    -11.58579130  0.6445E-03  0.5748E-03
ETOT:  15    -11.58587175  0.8045E-04  0.1339E-03
ETOT:  16    -11.58587362  0.1875E-05  0.4189E-04
ETOT:  17    -11.58587413  0.5126E-06  0.1486E-04
ETOT:  18    -11.58587420  0.6892E-07  0.3361E-05
ETOT:  19    -11.58587421  0.4461E-08  0.1102E-05
ETOT:  20    -11.58587421  0.1802E-08  0.3413E-06
ETOT:  21    -11.58587421  0.2176E-09  0.1103E-06
ETOT:  22    -11.58587421  0.3718E-10  0.4886E-07
   ETOT(Q) + SM_energy =  -11.5858742100388        8.35680549919333

実空間の電荷密度

icond=9として再計算して実空間の電荷密度nfchgt_r.dataを作成する。

Gaussian cube fileの作成

/home/hamada/STATE/tools/ChargeUtil/chg2cubeを用いて電荷密度をGaussian cube形式に変換する。

$ chg2cube
Enter the name of the STATE input file> your_input_file_name
Enter the name of the charge density file> nfchgt_r.data
Enter the prefix for output file(s)> charge
reading nfchgt_r.data...done.
# of data in nfchgt_r.data : 1
data #1  : ispin : 1  : ilevel: 1
Creating charge.cube...done.
Program successfully ended

chg2cubeの実行後、対話的に

  • STATEの入力ファイルの名前
  • 電荷密度ファイルの名前(nfchgt_r.data)
  • cubeファイルの名前(上ではcharge)

を入力するとcharge.cubeが生成される。

Gaussian cube形式に関しては例えばここを参照。

Bader電荷の計算

上で得られたGaussian cube形式の電荷密度ファイルcharge.cubeを用いてBader電荷を計算する。

$ bader charge.cube

  GRID BASED BADER ANALYSIS  (Version 0.95a 02/26/16)

  OPEN ... charge.cube
  GAUSSIAN-STYLE INPUT FILE
  FFT-grid:    30 x  30 x 128
  CLOSE ... charge.cube
  RUN TIME:    0.03 SECONDS

  CALCULATING BADER CHARGE DISTRIBUTION
                 0  10  25  50  75  100
  PERCENT DONE:  **********************

  REFINING AUTOMATICALLY
  ITERATION: 1
  EDGE POINTS:         23177
  REASSIGNED POINTS:     969
  ITERATION: 2
  CHECKED POINTS:       9850
  REASSIGNED POINTS:       0

  RUN TIME:       0.15 SECONDS

  CALCULATING MINIMUM DISTANCES TO ATOMS
                 0  10  25  50  75  100
  PERCENT DONE:  **********************
  RUN TIME:    0.01 SECONDS

  WRITING BADER ATOMIC CHARGES TO ACF.dat
  WRITING BADER VOLUME CHARGES TO BCF.dat

  NUMBER OF BADER MAXIMA FOUND:              6
      SIGNIFICANT MAXIMA FOUND:              6
                 VACUUM CHARGE:         0.0205
           NUMBER OF ELECTRONS:        8.00003

オプション

プログラムbaderのオプションはbader -hで確認できる:

$ bader -h

  GRID BASED BADER ANALYSIS  (Version 0.95a 02/26/16)

 Description of flags

    -c | -n  < bader | voronoi >
         Turn on [-c] or off [-n] the following calculations
            bader: Bader atoms in molecules (default)
            voronoi: population analysis based on distance

    -b < neargrid | ongrid | weight >
         Use the default near-grid bader partitioning, the
         original on-grid based algorithm, or the weight method
         of Yu and Trinkle

    -r < refine_edge_method >
         By default (-r -1) , only the points around reassigned
         points are checked during refinements. The old method,
         which checks every edge point during each refinement, can
         be enabled using the -r -2 switch:
            bader -r -2 CHGCAR

    -ref < reference_charge >
         Use a reference charge file to do the Bader partitioning.
         This is the recommended way to analyze vasp output files:
            bader CHGCAR -ref CHGCAR_total
         where CHGCAR_total is the sum of AECCAR0 and AECCAR2.

    -vac < off | auto | vacuum_density >
         Assign low density points to vacuum.
           auto: vacuum density cutoff is 1E-3 e/Ang^3 by default
           off: do not assign low density points to a vacuum volume
           vacuum_density: maximum density assigned to a vacuum volume

    -m < known | max >
         Determines how trajectories terminate
            known: stop when a point is surrounded by known points
            max: stop only when a charge density maximum is reached

    -p < all_atom | all_bader >
         Print calculated Bader volumes
            all_atom: all atomic volumes
            all_bader: all Bader volumes

    -p < sel_atom | sel_bader > <volume list or range>
         Print calculated Bader volumes
            sel_atom: atomic volume(s) around the selected atom(s)
            sel_bader: selected Bader volumes

    -p < sum_atom | sum_bader > <volume list or range>
         Print calculated Bader volumes
            sum_atom: sum of atomic volume(s) around the selected atom(s)
            sum_bader: sum of selected Bader volumes

    -p < atom_index | bader_index >
         Print index of atomic or Bader volumes
            atom_index: print atomic volume indicies
            bader_index: print Bader volume indicies

    -i < cube | chgcar >
         Input charge density file type.  If not specified, the
         program will try to determine the charge density file type
         automatically.

    -h
         Help.

    -v
         Verbose output.

    -cp
         Find critical points of the charge density
         Calculate eigenvalues and eigenvectors at those points
         Store results in file named CPF.dat

結果

Bader電荷の計算結果はACF.datに出力される:

$ cat ACF.dat
   #         X           Y           Z        CHARGE     MIN DIST    ATOMIC VOL
--------------------------------------------------------------------------------
   1   2.3277690   1.3439380   0.0000000   3.9020033   1.0751426  65.2790699
   2   0.0000000   2.6878760   0.0000000   4.0774900   1.2120244  69.3785643
--------------------------------------------------------------------------------
   VACUUM CHARGE:               0.0205
   VACUUM VOLUME:             240.7492
   NUMBER OF ELECTRONS:         8.0000

CHARGEに各原子の電荷が表示される。単位胞内の2つの炭素原子は等価なのに電荷に偏りが生じているのは、電荷密度のグリッドが粗いためと考えられる。

カットオフ依存性

エネルギーカットオフを増加させて(電荷密度のグリッドを細かくして)電荷分布の収束性を調べる。

gmaxp2030405060708090100
grid-X3048608090108120144150
grid-Y3048608090108120144150
grid-Y128192256320384450512576640
炭素原子1の電荷3.9023.9343.9453.9563.9603.9653.9673.9713.972
炭素原子2の電荷4.0774.0444.0344.0224.0194.0144.0114.0084.007
電荷の偏り0.1750.1100.0880.0660.0590.0490.0440.0370.035
電荷の偏り(%)4.3752.7502.2001.6501.4751.2251.1000.9250.875

カットオフの増加とともに電荷の偏りが減少し精度が向上していることが確認できる。

グラフェン 1x2/√3

gra1x2_r3.png

入力ファイル

0 0 0 0 0 0
8.0000  20.0000   1   4   4   : GMAX, GMAXP, NTYP, NATM, NATM2
1  0  : number of space group, type of bravis lattice
8.063628983  4.655538364  20.000  90.0   90.0  90.0  : a,b,c,alpha,beta,gamma
12    24       1       1      1     1     : knx,kny,knz, k-point shift
1   0   : NCORD, NINV,   : IWEI, IMDTYP, ITYP
 1.3439381638      0.0000000000      0.0000000000  1   1  1
-1.3439381638      0.0000000000      0.0000000000  1   1  1
 2.6878763276      2.3277691820      0.0000000000  1   1  1
-2.6878763276      2.3277691820      0.0000000000  1   1  1
6  0.1500  12.0107 1 1 0.d0  : TYPE 1IATOMN,ALFA,AMION,ILOC,IVAN
0 0 0 0 0 : ICOND 0-MD, 1-CONT.MD, 2-WAVE FN,, 3-WAVE FN CONT.,  iconstpw
    0   1   : IPRE, IPRI
    200   200     0    3600.00  0 : NMD1, NMD2, iter_last, CPUMAX,ifstop
    3   1  : Simple=1,Broyd2=3,Blugel=6,    1:charge, 2:potential mix.
    0   20  0.8 : starting mixing, kbxmix,alpha
   0.60  0.50  0.60  0.70  1.00    : DTIM1, DTIM2, DTIM3, DTIM4, dtim_last
 100.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
 101  : NBZTYP 0-SF, 1-BK, 2-SC, 3-BCC, 4-FCC, 5-DIA, 6-HEX
    0   0   0   : NKX, NKY, NKZ
    0   0   0   : NKX2,NKY2,NKZ2
   24           : NEG
        1       : NEXTST(MB)
        0       : 0; random numbers, 1; matrix diagon
        2      0      0      0(MB) : imsd, i_2lm, i_sd2another, wksz for phase
        0       : evaluation of eko difference.0 = no ,1 = yes
        0                       : npdosao
        0    0.0     : SM_N, DOPPING

計算結果

$ grep ETOT nfout_1|grep "[0-9]"
ETOT:   1     26.21341351  0.2621E+02  0.3937E-01
ETOT:   2     -5.07267018  0.3129E+02  0.1212E-01
ETOT:   3    -15.33757551  0.1026E+02  0.1318E-01
ETOT:   4    -19.66071003  0.4323E+01  0.1333E-01
ETOT:   5    -22.00674384  0.2346E+01  0.1766E-01
ETOT:   6    -22.88457509  0.8778E+00  0.1587E-01
ETOT:   7    -22.92511351  0.4054E-01  0.1501E-01
ETOT:   8    -22.97540886  0.5030E-01  0.1372E-01
ETOT:   9    -22.95949325  0.1592E-01  0.1391E-01
ETOT:  10    -22.97376104  0.1427E-01  0.1353E-01
ETOT:  11    -23.04325978  0.6950E-01  0.9979E-02
ETOT:  12    -23.15572141  0.1125E+00  0.3768E-02
ETOT:  13    -23.16991809  0.1420E-01  0.1487E-02
ETOT:  14    -23.17160528  0.1687E-02  0.4518E-03
ETOT:  15    -23.17172043  0.1151E-03  0.1618E-03
ETOT:  16    -23.17173717  0.1673E-04  0.7242E-04
ETOT:  17    -23.17174331  0.6144E-05  0.2139E-04
ETOT:  18    -23.17174360  0.2898E-06  0.7831E-05
ETOT:  19    -23.17174357  0.3402E-07  0.3759E-05
ETOT:  20    -23.17174366  0.8909E-07  0.2707E-05
ETOT:  21    -23.17174376  0.1035E-06  0.1720E-05
ETOT:  22    -23.17174377  0.1106E-07  0.1275E-05
ETOT:  23    -23.17174379  0.1527E-07  0.1788E-05
ETOT:  24    -23.17174380  0.1140E-07  0.6630E-06
ETOT:  25    -23.17174380  0.3104E-08  0.3894E-06
ETOT:  26    -23.17174380  0.7199E-09  0.2579E-06
ETOT:  27    -23.17174380  0.6648E-10  0.5746E-07
   ETOT(Q) + SM_energy =  -23.1717438004519        16.7136110440348

カットオフ依存性

gmaxp20406080100
grid-X54108160216270
grid-Y306090120150
grid-Y128256384512640
炭素原子1の電荷3.8903.9673.9563.9753.969
炭素原子2の電荷4.0884.0114.0234.0234.009
炭素原子3の電荷4.0884.0114.0234.0234.009
炭素原子4の電荷3.8903.9673.9563.9753.969
電荷の偏り0.1980.0430.0660.0270.039
電荷の偏り(%)4.9501.0751.6500.6750.975
トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2023-01-23 (月) 11:57:01 (462d)