EFIELD_PEAD E-field response for excited local state

[oscarlb]

Dear Support,

I am trying to calculate the self-consistent response to an electric field of a wide-bandgap semiconductor containing a defect (silicon vacancy in SiC). I am using the EFIELD_PEAD flag for this purpose. However, the calculation I am trying to do skips the PEAD evaluation because the system is classified as a metal for an excited state.

As a defect system, I am interested in the behavior of the defect states localized deep in the bandgap, and want to calculate the E-field response even for an excited defect state. I therefore set the FERWE/FERDO flags to change the occupation within the defect states (no valence/conduction band states are involved). But due to this, there is a gap in occupation of the system and the isMetal function returns TRUE due to this, although there is no change in conduction band states. The system is still sufficiently insulating for the applied field I believe, it is just a different localized wavefunction being occupied.

I realize this is a special case for designed routine, but do you think there is a way to perform the calculation I want, despite this issue? Or do you think the routine is unstable for such a change in the system? I am up for making some tweaks myself to circumvent the issue, but if there is some deeper obstacle I would greatly appreciate your guidance. Including failing run which should run on ~64 cores within a few minutes as an example.

[michael_wolloch]

Dear Oscar,

I made some small test calculations for your system, and I am quite certain that your issues have at least partially to do with insufficient k-point sampling.
When I increased your KPOINTS mesh to 2x2x2 (still Gamma-centered for PEAD), electronic convergence is much improved. I was able to switch to ALGO = Normal and converge in 12 steps. There is also a small but finite gap of about 86 meV, which seems to be enough for PEAD calculations.

Note that so far I have not set FERWE/FERDO, and I also tweaked other settings a bit. E.g., I found no magnetism for the ground state, so I removed ISPIN = 2. Generally, I would recommend keeping your INCAR file as minimal as possible.

Here are my preliminary KPOINTS and INCAR files:

Automatic mesh
0
Gamma
2 2 2

1. 0. 0.

SYSTEM = PBE of T_center in SiC
LASPH=.TRUE.
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1E-06
NSW = 0
LREAL = Auto
ALGO = Normal
PREC = Accurate
KPAR = 4
NCORE=1
EFIELD_PEAD = 0 0 1E-4
SKIP_EDOTP = .TRUE

Please check if your results are converged with respect to the k-point sampling. Just because it works with a 2x2x2 mesh, does not mean that the results are correct.
Please also test if the excited state calculations works as expected and let me know if you need anything else. I am happy to do some additional testing with FERWE/FERDO if needed.

Cheers, Michael

[oscarlb]

Dear Michael,
I appreciate your effort in testing my calculation, but the FERWE/FERDO tags are absolutely vital to the issue. Sorry if I was unclear.

I have no issues running the ground state, for now. I run very small systems and in gamma-point in order to try out the functionality, and convergence tests are planned if I can get the PEAD algorithm to work. My issue is very precisely that the PEAD routine refuses to run when setting FERWE/FERDO to an excited state of the defect. Error/warning:

The calculation of changes in the macroscopic polarization due to the application of a finite electric field, by means of the PEAD method (LCALCEPS=.TRUE. or EFIELD_PEAD/=0), requires your system to be insulating. This does not seem to be the case. VASP will skip this part of your job, sorry...

I tried running with higher K-points (see added input/output), but made no difference. If I look at the source, it looks to me like the criteria for the warning (and skipping PEAD) is having any gaps in occupation. This is an issue for me, who wants to try running precisely when I have a hole and an excited electron inside the bandgap of the semiconductor.

I am looking for insight into how to achieve the PEAD calculation even for such an occupation, or if it is unphysical/unstable to start with due to the internals of the PEAD routine. Would you be able to provide a workaround or guidance to how I could tweak the routines? Or do you know beforehand that the PEAD technique simply cannot handle excited states?

Best regards,
Oscar

PS:
Many system settings I had in the file are vital for my (point-defect) system, which understandably can seem like over-specifying the calculation. But any issues fundamentally arising from these tags, would also defeat the purpose of my calculation, so I would include the ones below in any test.

NELECT: sets the charge state of the system to -e. The added charge will be localized to the defect, and add a monopole to the system.
ISPIN: The defect actually has an effective spin of 3/2 in the intended charge state, and the ISPIN=2 is vital to capture this, and it is indeed found in my example (maybe 1/2 if done in the neutral state). The excitation I need is only in the spin-down channel.
LPEAD_RELAX: I am actually trying to relax the structure accounting for the electric field. But I am stuck on the step before that. It also limits the ALGO (https://www.vasp.at/forum/viewtopic.php?t=19899).
ISYM: The structural relaxation of the defect can reduce the symmetry of the system, so as standard is not determined/constrained in the calculation (ISYM=0).

There is also a small but finite gap of about 86 meV, which seems to be enough for PEAD calculations

Just to avoid confusion: The gap of 86 meV is not the bandgap of the system. The true bandgap of the system is 2.2 eV (with PBE) because the first conduction band state is the first non-localized unoccupied states and there may be several unoccupied defect states. This is what I think is creating incompatibility with the metallic checks currently implemented. What you specify is simply the difference between localized, dispersionless, defect levels, as could be seen in the bandplot if I had it available. It would be fine if the occupation between these states were to change due to an applied electric field. It would simply reflect charge transport within the defect and it would be interesting, but most importantly it would not imply any metallic behavior. But I will keep this in mind if the PEAD algorithm will behave unphysically if I breach that defect state gap.

[michael_wolloch]

Dear Oscar,

I am sorry that I misunderstood your initial post. I am usually trying to start with the simplest possible version of a problem, and even though I did not have much insight last evening, I wanted to get back to you asap. Thank you very much for the clarification.
I now understand what you are trying to do, but unfortunately, this is not possible with the PEAD routines in VASP.

While the system you are trying to simulate is indeed not metallic, the function we use for checking that (isMetal in wave.F) assumes that the occupied bands are grouped together, as you noticed yourself. It does NOT deal correctly with the mid-gap state you want to introduce.

Fixing this function would not be a major issue, but the problem is that the PEAD routines themselves rely on the order of the occupied bands. They compute the derivative of the Bloch orbitals with respect to the Bloch vectors by finite differences by summing over occupied bands (see equations in the LPEAD documentation on the wiki). It is assumed that the occupied bands are bunched together. So, even if one would change the isMetal function to correctly take occupied gap states into account, the PEAD routines would not work with your system.

Reworking the PEAD routines so that they could handle non-continuously occupied states might be possible in theory, but would require significant effort from a developer familiar with those routines (which unfortunately does not include me). This effort is very unlikely to be undertaken soon.

If you would be willing and able to do this development work yourself, I can ask a colleague if they would be willing to jump on a Zoom call with you and discuss some general ideas.

Sorry that there is no better news,
Michael

[oscarlb]

Dear Michael,

Thank you for telling me. I was just looking for clarity, and I understood it is a special case I was asking for.

However, I am still somewhat hopeful it could be done in my special case, and I would like to try. But I should familiarize myself more with the methods applied in the PEAD algorithm and I will update this topic if I think I grasp it to the level a discussion with your developers would be productive.

Naively speaking, I think it should work as long as the energy ordering of the state being occupied is well defined across the Brillouin zone (mostly okay for defect states as they are dispersionless up to the supercell convergence, maybe degenerate states are more of an issue though). Although there issue with nearly degenerate valence bands switching should be addressed similarly already. The procedure is supposed to use the occupied states weighted with the occupation numbers according to Resta, but I understand the current implementation assumes these are only the valence band states. Finite-difference integration over k-space or whatever trick is used should be fine if we just perform the calculation over occupied states, defined as such, perhaps accomplished by a remapping of the wavefunction indices somewhere or (conceptually speaking) gathering the wavefunctions needed separately. However, I should make sure no assumption in the tricks employed is broken because of this (the way of transforming the 1D k-space integration into a logarithm of products of determinants (eq 15 of King-Smith and Vanderbilt) is especially unclear to me).

Thank you for your help. I hope we'll keep each other posted.

Best regards,
Oscar