Fermi energy too high?

[Phil G.]

Dear people and experts of CP2K,

after the geometry optimization of the lithium niobate (LiNbO3) unit cell I would like to obtain pdos in order to determine the band gap and Fermi energy of the bulk system.

After the calculation with ENERGY_FORCE I got pdos files of the three atoms (indexing depends on the z-position of the atoms) and I'm wondering about the value of Fermi energy: E_F = 0.300174 a.u. which is 8.168 eV. Is that not too high? And which energy has the value 0 and what is the reference? What is the Fermi energy defined in the language of CP2K?

The energy band gap (HOMO-LUMO gap) of 3.62 eV agrees well with experimental values of 3.7 to 3.9 eV. But I cannot imagine that Fermi level has too high energy values.

Has anyone an idea what is the reason for such high Fermi energy values?

Here the input and output files are attached here.

Kind regards,

Phil

[Matt W]

In a periodic system the zero of the one electron levels is arbitrary. If you need a reference you need to run a slab calculation with vacuum or try to align semi-core states to something.

Matt

[Phil G.]

Dear Matt,

thank you for your reply and good suggestions. Now I have let different LiNbO3 slab systems to be calculated:

a) 14 trilayer system as from Sanna et al., Appl. Surf. Sci. 301 (2014), 70-78 with Nb-O3-Li2 surface termination on the one side of the slab and Li-O surface termination on the other side. Vacuum space of at least 40 Angstroms was included. The bulk region was already geometry-optimized and bulk atoms were fixed in the inner 6 trilayers. Geometry optimization on the whole slab system was performed and then the pdos of the system was calculated and plotted for every atom layers.

Result: E_F = 0.1552 eV  (fermi energy is overall constant, in every atom layers)

b) the same as a), but the bulk region was not already geometry-optimized before. Geometry optimization was performed and calculation of pdos.

Result: E_F = - 0.8516 eV

c) the same as b), but 26 trilayers instead of 14 trilayers. Geometry optimization and calculation of pdos were performed.

Result: E_F = 2.3372 eV

So, I am wondering why these values differ so much. Should I need band structure calculation of the bulk LiNbO3 in order to find the global valence band edge maximum (with KPOINT calculation)?

Kind regards,

Phil

[Matt W]

Hello again,

did you find the potential in the vacuum and align to that? You need to set a reference to get absolute values.

You could also try using the wavelet solver 

&POISSON

   PSOLVER WAVELET

   PERIODIC XZ
&END

and PERIODIC XZ  in the &CELL section. The Y direction must be the non-periodic one. That gives an absolute reference (if there is no dipole in the cell otherwise you need the  dipole correction switched on).

Matt

[Phil G.]

Dear Matt,

how can find the potential in the vacuum (which type of potential? potential energy or electric/electrostatic potential?) ?

For the case of electric/electrostatic potential, there is a flat curve with a step near the vacuum center as a consequence of dipole correction in Z direction, while in the bulk slab there is a periodic curve.

I will try to use the wavelet solver with PERIODIC XY.

Phil

[Phil G.]

Hello again,

have tried some attempts to start calculation with WAVELET poisson solver, but all attempts have failed due to following error messages:

1)  the FFT in the x direction is not allowed
     n01 dimension         154

     (pw/ps_wavelet_util.F:358)

     ===== Routine Calling Stack =====

           13 S_FFT_dimensions
           12 RS_z_slice_distribution
           11 ps_wavelet_create
     the FFT in the x direction is not allowed
     n01 dimension         154
           10 pw_poisson_rebuild
            9 pw_poisson_solve
            8 qs_ks_build_kohn_sham_matrix
            7 rebuild_ks_matrix
     the FFT in the x direction is not allowed
     n01 dimension         154
     the FFT in the x direction is not allowed
     n01 dimension         154
     the FFT in the x direction is not allowed
     n01 dimension         154
     the FFT in the x direction is not allowed
     n01 dimension         154
            6 qs_ks_update_qs_env
            5 scf_env_do_scf_inner_loop
            4 scf_env_do_scf
            3 qs_energies
            2 qs_forces
            1 CP2K

2)  after that I turn off the command EXTENDED_FFT_LENGTHS, then:

    Index to radix array not found.

    (pw/fft_tools.F:293)

     ===== Routine Calling Stack =====

            6 pw_grid_setup
            5 pw_env_rebuild
            4 qs_env_rebuild_pw_env
            3 qs_env_setup
            2 qs_init_subsys
            1 CP2K

That's strange and I don't know what to do.

In my input file there are some info about commands:

[...]

    SURFACE_DIPOLE_CORRECTION .TRUE.
    SURF_DIP_DIR Y

[...]

    &MGRID
      CUTOFF 600
      NGRIDS 5
      REL_CUTOFF 50
    &END MGRID

[...]

    &POISSON
      POISSON_SOLVER WAVELET
      PERIODIC XZ
    &END POISSON

[...]

[...]

    &CELL
      A    5.148      0.0         0.0
      B    0.000      100.0      0.0
      C    0.0          0.0         8.9166        
      PERIODIC XZ
    &END CELL

[...]

Phil

[Matt W]

Ah, OK. The extended FFT lengths only works with FFTW not with the wavelet FFT.

You do not need such a large vacuum with the wavelet solver as it is genuinely non-periodic. Place the slab in the center and 5A of vacuum either side should be sufficient - allow 10 A either side to get a clear decay to vacuum level(s). You will get two vacuum levels if you have  a dipole. I don't know if wavelet will work with the dipole correction, I'd turn it off to start with.

Matt

[Phil G.]

Dear Matt,

thank you for the suggestions and after centering the slab in y-direction and turning off the surface dipole correction, the program finally runs and I get the result, but there are some error messages such as:

 *** WARNING in pw/ps_wavelet_methods.F:236 :: Density non-zero on the ***
 *** edges of the unit cell: wrong results in WAVELET solver           ***

    89 Broy./Diag. 0.10E+00    2.3     0.00000880     -3020.8753805067  4.85E-05

  *** SCF run converged in    89 steps ***

I have chosen the cell length of 40 angstroms in Y direction (slab length in y-direction is about 27.6 angstroms).

For the LiNbO3 slab consisting of 14 trilayers as stated in the message of 9th January, I obtain the result of the Fermi energy:

  E_F = 11.174 eV    (in comparison to the -0.8516 eV in case b) )

This is unrealistic...should I have to enlarge the cell length in y-direction or should I turn on the dipole correction?

Phil

[Matt W]

Is your slab in the centre of the cell (Y direction)? -  the cell runs from 0 to L, so the slab must be centred at L/2.

Matt

[Phil G.]

yes, the slab was centered in y-direction (xyz file attached here), but in x- and z- direction from 0 to the corresponding lengths of the slab.

Maybe I should extend the size of vacuum space or turn on the dipole correction?

For searching for reason, the input file is also attached here.

Phil

[Matt W]

Hi,

the xyz file you sent doesn't have the y coords centred in the cell. The cell goes from 0 to L. You centre as if it goes from -L/2 to L/2.

If you print V_HARTREE_CUBE you can work out the potential profile across the slab. 

BTW your slab is too small in X and Z for reliable results without k-point sampling...

Matt

[Matt W]

I added

&CENTER_COORDS

&END

in the topology section. I get the 'fermi level' reported at -4.61 eV. And the vacuum levels at +- 0.75 eV

So the work functions of the two faces are ~5.4 and 3.8 eV respectively. I got the profile of the potential by adding &V_HARTREE_CUBE and using the cubecruncher utility (bundled in the tools directory) to extract the averaged y profile to find the vacuum levels.

Matt

[Phil G.]

Dear Matt,

thank you very much for the suggestions and after some tests it also worked ! :-)

To check the effect of dipole correction on the calculation, I started different slab systems and used cell length in y-direction at 50 Angstrom (much more than this value will lead to error (index to radix array not found)):

The slab systems are defined as in the first message above (9th January)

 

 a)  without dipole correction:  E_F = -3.31475 eV

      with dipole correction:  E_F = -3.2983 eV

 b)  without dipole correction:  E_F = -4.642 eV

      with dipole correction:  E_F = -4.72 eV

It seems that the dipole correction leads to a marginal reduction of the Fermi energy.

Now I want to check, if the vacuum level has another value that the Fermi energy as you showed (+/- 0.75 eV) by plotting pdos as a function of atom coordinates, i.e. in slab y-direction.

Thanks again for your helpful suggestions! :-)

Kind regards,

Phil

[Phil G.]

Dear Matt (or someone in the CP2K group),

after I get the 'reliable' Fermi level and the 'vacuum levels' which are the same as Matt got them (-4.61 eV and +- 0.75 eV, see message on this topic on 16th January). But I am bit confusing, why the vacuum levels are of the same value (but with opposite signs) on both surfaces of the slab. Should they do depend on the electronic properties of the surface? What is the vacuum level exactly mean? (I mean that their values are not zero and what correlates to the zero level?) I want to plot the PDOS against the certain energy level scale as defined in the electrochemistry, e.g. NHE (normal hydrgen electrode), SHE (standard hydrogen electrode). Is the zero level in the potential plot defined as the vacuum level at infinity distance and the both vacuum levels (+- 0.75 eV) the local vacuum levels?

Important question: How was the Fermi level or the Fermi energy of a slab (with two surfaces) in the CP2K calculated? Is it the same as the energy where the valence band maximum is (globally?) or is it defined as the value where the probability for an electron to occupy between conductance and valence band edges is exactly 0.5 (according to the Fermi-Dirac distribution function)?

If I get the pdos file of the slab, I get only a value of the Fermi energy. Is that the average value for the Fermi energies throughout the whole slab or is that rather the bulk property? (Because I would like to plot PDOS as a function of atom layer coordinates (I have inset integer indices to the atom types referring to the atom layers along the axis perpendicular to their surfaces), and you can see that the energy positions of the band edges change with the atom layer coordinate. It seems to be like a graphical plot of band bending on both surface sites, but the energy position change is overall linear throughout the slab (see attached image, the Fermi level was shifted to zero in the energy scale). In my opinion, this is not the realistic/reproducible band bending diagram as we know it from the semiconductor heterostructures (junction between semiconductor and a medium (solid, electrolyte or something)). Is it rather the result of constant potential difference between the slab surfaces (at position z = 0 and z = - 54 Angstrom), and it seems to be the same irrespective of the slab lengths (26 Angstroms). The potential difference is conceived to be determined by the polarization of the slab material. How can I change the constant potential model to the constant (surface) charge model and how to implement it into the CP2K code?

Kind regards,

Phil

[Matt W]

Hi Phil, 

a lot of questions / good points here, mainly physics rather than CP2K, I guess. Quick reply, I'll try and do a bit better later (if I remember).

The fermi level depends on the algorithm. If there is no smearing, the reported fermi level is the HOMO of the system. If smearing, it is really the fermi level at the given electronic temperature.

Maybe try looking at papers by Jacek Goniakowski (metal insulator systems) or  Emilio Artacho (2d metals) that discuss some issues with semi-infinite slab models / dipoles / electostatics etc.

If you can clearly define the boundary conditions that you actually want then we can probably suggest how to get them in CP2K.

Matt

[Phil G.]

Dear Matt,

The Fermi smearing was turned on during the calculations, so probably the Fermi level will be the electronic level of 50% occupation probability of electrons (lies in the band gap between the band edges).

The slab used in the calculation has to mimic the ferroelectric single crystal with the two polar surfaces exposed to vacuum. The both polar surfaces have the opposite surface charges due to the spontaneous polarization of the bulk crystal (in the middle region of the slab), so for example the top surface (in the calculation using wavelet poisson solver it is the positively y-directed surface) is positively polarized and the bottom surface (negatively y-directed surface in this calculation) is negatively polarized. The polarization is aligned along the polar axis of the crystal. The value of the surface charge depends on the magnitude of the spontaneous polarization in the bulk region (so: sigma = P_s * normalvector ). The slab system has to fulfill the 2D periodic boundary (surface problem, in this calculation x-z periodicity).

I would like to calculate the electronic (structure) properties of the single crystal surface exposed to vacuum (air) and later to water in order to get the Fermi level, work function, electron affinity, ionization energy of both the polar surfaces.

Kind regards,

Phil

[Harsha]

Hi Matt or Phill,

I have a question based on your conversion above. How do I extract vacuum level?

I printed &V_HARTREE_CUBE file used cubecruncher to extract the averaged Y profile.

But how do I know the vacuum level? Is it printed somewhere like the Fermi level which is printed in the output file or is that something that we have to get from the plot? 

Best regards

Harsha

[Matt W]

You need to look at the values on the cell edges of the non-periodic direction.  If you get an planar avereraged potential from cube cruncher you should see it plateus either side of system (depending on charge / dipole etc in your cell).