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

Pyrite (FeS2) †

In this section, how to perform the electronic structure analysis is described using pyrite (FeS2) as an example.

SCF calculation †

First of all, let us perform an SCF calculation using pw.x to get wave functions and charge density. If necessary, let us perform the geometry optimization. Below is the input file used for the SCF calculation

scf.in

&CONTROL
  calculation   = 'scf'
  etot_conv_thr = 1.2000000000d-04
  forc_conv_thr = 1.0000000000d-04
  outdir        = './out/'
  prefix        = 'fes2'
  pseudo_dir    = '../pseudo/'
  tprnfor       = .true.
  tstress       = .true.
  verbosity     = 'high'
/
&SYSTEM
  degauss     = 2.0000000000d-02
  ecutrho     = 1.0800000000d+03
  ecutwfc     = 9.0000000000d+01
  ibrav       = 0
  nat         = 12
  nosym       = .false.
 !nspin       = 2
  nbnd        = 80
  ntyp        = 2
  occupations = 'smearing'
  smearing    = 'cold'
  !starting_magnetization(1) =   3.1250000000d-01
  !starting_magnetization(2) =   1.0000000000d-01
/
&ELECTRONS
  conv_thr         = 2.4000000000d-09
  electron_maxstep = 80
  mixing_beta      = 4.0000000000d-01
/
ATOMIC_SPECIES
Fe     55.845 Fe.pbe-spn-kjpaw_psl.0.2.1.UPF
S      32.065 s_pbe_v1.4.uspp.F.UPF
ATOMIC_POSITIONS crystal
Fe           0.0000000000       0.0000000000       0.0000000000
Fe           0.5000000000       0.0000000000       0.5000000000
Fe           0.0000000000       0.5000000000       0.5000000000
Fe           0.5000000000       0.5000000000       0.0000000000
S            0.3850400000       0.3850400000       0.3850400000
S            0.6149600000       0.6149600000       0.6149600000
S            0.1149600000       0.6149600000       0.8850400000
S            0.8850400000       0.3850400000       0.1149600000
S            0.6149600000       0.8850400000       0.1149600000
S            0.3850400000       0.1149600000       0.8850400000
S            0.8850400000       0.1149600000       0.6149600000
S            0.1149600000       0.8850400000       0.3850400000
K_POINTS automatic
4 4 4 0 0 0
CELL_PARAMETERS angstrom
      5.4281000000       0.0000000000       0.0000000000
      0.0000000000       5.4281000000       0.0000000000
      0.0000000000       0.0000000000       5.4281000000

↑

Density of states calculation †

After confirming the convergence of SCF, let us perform density of states (DOS) calculation using dos.x. Below is the input file for the DOS calculation.

dos.in

&DOS
  outdir ='./out/'
  prefix ='fes2'
  fildos ='fes2.dos'
  bz_sum = 'tetrahedra'
  Emin   = -85.0
  Emax   = 25.0
  DeltaE = 0.01
/

In this exercise, the SCF calculation was performed with the smearing method, while DOS is calculated using the tetrahedron method.

↑

Refined Density of states calculation †

To get more precise, let us perform a non-SCF (NSCF) calculation. This is not always necessary, but if necessary, perform a NSCF with a finer k-point mesh. Remember this NSCF calculation can take longer than SCF calculation, depending on the number of k-point.
For this purpose set calculation nscf in the &CONTROL namelist as

   calculation   = 'nscf'

and finer k-point grid in the K_POINTS card, for example:

K_POINTS automatic
9 9 9 0 0 0

then the nscf calculation is done, perform the DOS calculation using dos.x.

↑

Projected DOS calculation †

To get DOS projected onto the atomic orbitals (PDOS), let us use projwfc.x. It is important to remember that unlike the DOS calculation, to get "accurate" PDOS, the preceding SCF or NSCF calculation should be done using the tetrahedron method, otherwise PDOS are calculated using the smearing method with the Gaussian function to approximate the delta function. Below is an input file for the PDOS calculation

projwfc.in
```
&PROJWFC
 outdir  = './out/'
 prefix  = 'fes2'
 Emin    = -25.00
 Emax    =  25.00
 DeltaE  =   0.01
/
```
For a better characterization of PDOS, it is useful to use rotated atomic orbitals in such a way that the occupation matrix is diagonalized. This can be done by setting diag_basis .true. as
```
 diag_basis = .true.
```

Furthermore, Lowdin population analysis is performed during the PDOS calculation. See the output file and search the word Lowdin Charges, which looks like:

Lowdin Charges:

     Atom #   1: total charge =  16.8478, s =  2.4992,
     Atom #   1: total charge =  16.8478, p =  7.3552, p1=  2.4500, p2=  2.4720, p3=  2.4332,
     Atom #   1: total charge =  16.8478, d =  6.9934, d1=  0.9573, d2=  0.8688, d3=  1.7250, d4=  1.6888, d5=  1.7535,
     Atom #   2: total charge =  16.8460, s =  2.4917,
     Atom #   2: total charge =  16.8460, p =  7.3512, p1=  2.4353, p2=  2.4631, p3=  2.4528,
     Atom #   2: total charge =  16.8460, d =  7.0031, d1=  0.9289, d2=  0.8879, d3=  1.7192, d4=  1.6805, d5=  1.7866,
     Atom #   3: total charge =  16.8204, s =  2.4916,
     Atom #   3: total charge =  16.8204, p =  7.3495, p1=  2.4327, p2=  2.4857, p3=  2.4311,
     Atom #   3: total charge =  16.8204, d =  6.9793, d1=  0.9570, d2=  0.8475, d3=  1.7116, d4=  1.6795, d5=  1.7836,
     Atom #   4: total charge =  16.8502, s =  2.4957,
     Atom #   4: total charge =  16.8502, p =  7.3523, p1=  2.4398, p2=  2.4296, p3=  2.4828,
     Atom #   4: total charge =  16.8502, d =  7.0022, d1=  0.9437, d2=  0.8749, d3=  1.7397, d4=  1.6725, d5=  1.7714,
     Atom #   5: total charge =   5.4643, s =  1.4831,
     Atom #   5: total charge =   5.4643, p =  3.9813, p1=  1.1366, p2=  1.4327, p3=  1.4119,
     Atom #   5: total charge =   5.4643, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Atom #   6: total charge =   5.4643, s =  1.4831,
     Atom #   6: total charge =   5.4643, p =  3.9813, p1=  1.1366, p2=  1.4320, p3=  1.4126,
     Atom #   6: total charge =   5.4643, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Atom #   7: total charge =   5.4600, s =  1.4921,
     Atom #   7: total charge =   5.4600, p =  3.9679, p1=  1.1348, p2=  1.4123, p3=  1.4208,
     Atom #   7: total charge =   5.4600, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Atom #   8: total charge =   5.4600, s =  1.4921,
     Atom #   8: total charge =   5.4600, p =  3.9679, p1=  1.1348, p2=  1.4122, p3=  1.4210,
     Atom #   8: total charge =   5.4600, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Atom #   9: total charge =   5.4591, s =  1.4921,
     Atom #   9: total charge =   5.4591, p =  3.9670, p1=  1.1314, p2=  1.4170, p3=  1.4186,
     Atom #   9: total charge =   5.4591, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Atom #  10: total charge =   5.4591, s =  1.4921,
     Atom #  10: total charge =   5.4591, p =  3.9670, p1=  1.1314, p2=  1.4188, p3=  1.4168,
     Atom #  10: total charge =   5.4591, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Atom #  11: total charge =   5.4606, s =  1.4873,
     Atom #  11: total charge =   5.4606, p =  3.9732, p1=  1.1304, p2=  1.3943, p3=  1.4486,
     Atom #  11: total charge =   5.4606, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Atom #  12: total charge =   5.4606, s =  1.4873,
     Atom #  12: total charge =   5.4606, p =  3.9732, p1=  1.1304, p2=  1.3957, p3=  1.4471,
     Atom #  12: total charge =   5.4606, d =  0.0000, d1=  0.0000, d2=  0.0000, d3=  0.0000, d4=  0.0000, d5=  0.0000,
     Spilling Parameter:   0.0085

This may be helpful to get an insight into the charge (state). However, the extreme care is required to make a conclusion on the atomic charge, as one can may see from the above result.

↑

Bader charge analysis †

To further gain an insight into the charge state/atomic charge, let us perform the Bader charge analysis.

↑

Preparation †

To perform the Bader charge analysis, the charge density of the system in the Gaussian cube formate is required. To do this, one needs to use pp.x. Below is an example for pp.x to generate a charge density file.

pp.in
```
&INPUTPP
  prefix     = 'fes2'
  outdir     = './out/'
  filplot    = 'fes2'
  plot_num   = 0
/
&PLOT
  iflag         = 3
  output_format = 6
  fileout       = 'fes2_val.cube'
/
```
plot_num is used to output the charge density in real space, iflag is used to specify that the dimension of the charge density, and output_format is used to specify the file format.

↑

Execution †

Supposing the path to the program bader is set, one can simply execute the following command

bader fes2_val.cube

Output may look like:

  GRID BASED BADER ANALYSIS  (Version 1.03 11/13/17)

  OPEN ... fes2_val.cube
  GAUSSIAN-STYLE INPUT FILE
  DENSITY-GRID:  108 x 108 x 108
  CLOSE ... fes2_val.cube
  RUN TIME:    0.27 SECONDS

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

  REFINING AUTOMATICALLY
  ITERATION: 1
  EDGE POINTS:        428375
  REASSIGNED POINTS:   28371

  RUN TIME:       2.29 SECONDS

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

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

  NUMBER OF BADER MAXIMA FOUND:             72
      SIGNIFICANT MAXIMA FOUND:             72
                 VACUUM CHARGE:         0.0000
           NUMBER OF ELECTRONS:      112.00003

The calculated Bader charges are written to ACF.dat as

    #         X           Y           Z       CHARGE      MIN DIST   ATOMIC VOL
 --------------------------------------------------------------------------------
    1   10.257622   10.257622   10.257622   15.308219     1.792039    65.277172
    2    5.128811   10.257622    5.128811   15.308219     1.792040    65.277172
    3   10.257622    5.128811    5.128811   15.308219     1.792040    65.277172
    4    5.128811    5.128811   10.257622   15.308219     1.792040    65.277172
    5    3.949595    3.949595    3.949595    6.354199     1.892313   102.381722
    6    6.308027    6.308027    6.308027    6.337589     1.877954   102.164956
    7    1.179216    6.308027    9.078406    6.337589     1.877956   102.164956
    8    9.078406    3.949595    1.179216    6.354199     1.892313   102.381722
    9    6.308027    9.078406    1.179216    6.337589     1.877956   102.164956
   10    3.949595    1.179216    9.078406    6.354199     1.892313   102.381722
   11    9.078406    1.179216    6.308027    6.354199     1.892311   102.381722
   12    1.179216    9.078406    3.949595    6.337589     1.877957   102.164956
 --------------------------------------------------------------------------------
    VACUUM CHARGE:               0.0000
    VACUUM VOLUME:               0.0000
    NUMBER OF ELECTRONS:       112.0000

Note the charge is based on the valence charge defined by the pseudopotentials used. In this example, we use the following pseudopotentials:

Fe.pbe-spn-kjpaw_psl.0.2.1.UPF
s_pbe_v1.4.uspp.F.UPF More specifically, the respective valence configurations can be found in the header of these potential files as

Fe.pbe-spn-kjpaw_psl.0.2.1.UPF

    nl pn  l   occ       Rcut    Rcut US       E pseu
    3S  1  0  2.00      1.100      1.300    -6.910119
    4S  2  0  2.00      0.800      1.400    -0.388933
    3P  2  1  6.00      1.000      1.300    -4.413015
    4P  3  1  0.00      1.000      1.600    -0.097407
    3D  3  2  6.00      1.400      2.000    -0.551558

s_pbe_v1.4.uspp.F.UPF
```
    nl pn  l   occ               Rcut            Rcut US             E pseu
    3S  3  0  2.00      0.00000000000      1.50000000000     -1.26833444300
    3P  3  1  4.00      0.00000000000      1.50000000000     -0.51513600900
```
It can be seen that the numbers of valence electros for Fe and S are 16 and 6, respectively. Based on simple mathematics, it seems Fe is positively charged, whereas S, negatively charged. Compare with the Lowdin charges. For more precise discussion on the actual charge state, one may need further analysis of PDOS (and others).

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

Menu