[CP2K-user] Fermi energy too high?
Matt W
mattwa... at gmail.com
Fri Apr 5 15:39:23 UTC 2019
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
On Friday, April 5, 2019 at 1:28:56 PM UTC+1, Phil G. wrote:
>
> 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
>
>
> Am Montag, 21. Januar 2019 12:25:54 UTC+1 schrieb 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
>>
>>
>> Am Mittwoch, 16. Januar 2019 14:20:18 UTC+1 schrieb 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
>>>
>>>
>>>
>>> On Wednesday, January 16, 2019 at 9:26:32 AM UTC, Matt W wrote:
>>>>
>>>> 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
>>>>
>>>> On Wednesday, January 16, 2019 at 8:54:56 AM UTC, Phil G. wrote:
>>>>>
>>>>> 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
>>>>>
>>>>> Am Dienstag, 15. Januar 2019 16:07:55 UTC+1 schrieb 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
>>>>>>
>>>>>> On Tuesday, January 15, 2019 at 2:21:26 PM UTC, Phil G. wrote:
>>>>>>>
>>>>>>> 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
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> Am Dienstag, 15. Januar 2019 12:31:22 UTC+1 schrieb 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
>>>>>>>>
>>>>>>>> On Tuesday, January 15, 2019 at 8:52:35 AM UTC, Phil G. wrote:
>>>>>>>>>
>>>>>>>>> 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
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Am Donnerstag, 10. Januar 2019 10:56:24 UTC+1 schrieb 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
>>>>>>>>>>
>>>>>>>>>> Am Mittwoch, 9. Januar 2019 14:05:04 UTC+1 schrieb 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
>>>>>>>>>>>
>>>>>>>>>>> On Wednesday, January 9, 2019 at 8:18:58 AM UTC, Phil G. wrote:
>>>>>>>>>>>>
>>>>>>>>>>>> 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
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> Am Freitag, 14. Dezember 2018 17:41:03 UTC+1 schrieb 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
>>>>>>>>>>>>>
>>>>>>>>>>>>> On Friday, December 14, 2018 at 4:33:13 PM UTC, Phil G. wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> 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
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20190405/37718e52/attachment.htm>
More information about the CP2K-user
mailing list