Charge in the dipole correction layer

[yujia_teng]

LDIPOL = .TRUE. puts a virtual dipole layer in the vacuum to correct the potential. Is it possible to obtain the charge of this virtual dipole from VASP directly? I'm trying some other ways to obtain dipole moment, such as from the potential profile (from LVHAR = .True.), which should match direct charge integration result as implemented in dipole correction.

I tried it myself but the result doesn't seem reliable. So in the figure below, I can find potential difference ΔV across the dipole layer and distance d of the potential drop region. Then I view the dipole as a parallel plate capacitor and use ΔV = σd/ε₀ to obtain surface charge density of the dipole. The charge is then obtained by multiplying the in-plane area.

But both ΔV and d change a small amount of value based on different choices, and approximate dipole layer as parallel plate capacitor seems very rough.

[andreas.singraber]

Hello!

The magnitude of the dipole moment which is introduced via LDIPOL=.TRUE. is listed in the OUTCAR file, look for the word dipolmoment (see also here):

Code: Select all

grep dipolmoment OUTCAR | tail -1

On a side note: Unfortunately this is probably one of the occasions where early VASP versions used an incorrect spelling of the word (actually it is the German spelling of "dipole moment") which later could not be changed any more because user scripts, third-party software, etc. rely on it. Sorry for that...

Regarding the plot you attached: I discussed this also with my colleagues and we think that the position of the potential step may be too close to the slab which could also cause the unreliable results you describe. Did you also set DIPOL explicitly in your INCAR file or did you let VASP determine this automatically?

It would be very helpful if you could upload the relevant VASP input and output files for your system! Thank you!

All the best,
Andreas Singraber

[yujia_teng]

Dear Andreas,
I attach the two kinds of input files here. The output files and the potential figure are also attached. One is for just ferroelectric material, and the other is heterostructure where FE material is combined with metal. I set DIPOL = 0.5 0.5 0.5 in heterostructure case while comment it out in freestanding layer case. The vacuum region should be 15-20 Å, so the dipole layer should be in the middle of vacuum region.

What I mean not reliable is my result extracted from VASP output data, not VASP result itself. Because I choose the region based on potential gradient. If the gradient is within a threshold, then it's considered as the dipole layer region. These values in my post processing script are are user-defined. So the result could change a little bit based on different parameters. The output is the 'find-kink.png' file.

Hope these make the confusing problem more clear.

[andreas.singraber]

Hello!

Ok, I think I understand that you want to compute the charge of a model capacitor which you imagine sitting at the position of the potential kink. Furthermore you assume that the thickness of the capacitor corresponds to the kink region and you determine this region by setting some thresholds on the potential (or potential change) and search from both sides. I agree that you can then compute the charge on this capacitor with the formula \(Q = \frac{\Delta V \varepsilon_0 A}{d}\) but of course it is dependent on your definition of the capacitor thickness \(d\).

I am not an expert in this kind of calculations but at this point I cannot follow your intentions any more: isn't the thickness of the kink region somewhat arbitrary because it depends on the implementation of the dipole correction scheme? As far as I understand the LDIPOL tag is used to remove the potential difference induced by the slab which is a problem in the presence of periodic boundary conditions. So I imagine that the "correction potential drop" \(\Delta V\) could have different extent (a.k.a. capacitor thickness) in z-direction, as long as it compensates exactly the \(\Delta V\) from the slab. However, I would assume there are some restrictions: it should not be too narrow because I guess then we will get convergence issues. On the other hand, it should not be too wide in order to stay isolated in the vacuum and prevent any interaction with the slab. Anyway, within some boundaries the model capacitor thickness would be arbitrary and therefore also the computed charge. Does it then still make sense to compute this charge? Wouldn't it make more sense to model the slab itself as the capacitor because its thickness is known (to some degree) and we can then compute the model capacitor charge?

Sorry if I am missing something here because of lack of experience... maybe you can elaborate some more why you need to compute this hypothetical capacitor charge and how the inaccuracies come into your calculation? Thank you!

All the best,
Andreas Singraber

[yujia_teng]

Dear Andreas,
Yes your 1st paragraph is what I'm trying to do.

I also agree that the kink region could vary, because LDIPOL just compensates ΔV. And as you said, the thickness d could vary, then the computed charge Q from this also varies. So I feel the result is unreliable. For this, I also tried SIESTA https://docs.siesta-project.org/project ... ric-fields, which is another first-principle package but based on atomic orbital instead of plane wave, and its dipole correction result (from tag Slab.DipoleCorrection) also gives the electric field from dipole correction:

siesta: Electric dipole (a.u.) = -0.000000 -0.000000 -0.182917
siesta: Electric dipole (Debye) = -0.000000 -0.000000 -0.464929
Dipole moment in unit cell = -0.0000 -0.0000 -0.4649 D
Electric field for dipole correction = 0.000000 0.000000 0.041922 eV/Ang/e

Since it gives electric field, we can use formula \(E=\dfrac{\sigma}{\varepsilon_0}=\dfrac{Q}{A\varepsilon_0}\), so \(Q=\varepsilon_0 AE\). All quantities are known, and there's no worry about the problem of varying of \(d\) here.

Looks like VASP doesn't output electric field from dipole correction, is it possible to add this feature in the future? For now, I think we can't get it from post-processing due to the varying \(d\) you mentioned.

Anyway, I did all of this is because the following: the screening of graphene is not perfect, so in the heterostructure I attached, there's still ΔV in vacuum arise from the dipole moment. In the perfect screening case, ΔV should be 0 in vacuum. So I think the imperfect screening is because the induced charge on top & bottom graphene is not enough. If I add the charge from the virtual dipole layer to graphene, the screening should be perfect. In other words, the dipole moment computed from graphene + virtual dipole layer charge (model the slab itself as the capacitor here) should equal to the result from integration \(p=\int \rho(r) (r-r_{center})d^3r\), where integration region is only around FE slab.

The root cause that drives me to do this is because it looks like VASP dipole moment output in OUTCAR is for the whole system only, and there's no way to know the dipole moment of a FE slab in heterostructure case directly from VASP output.

In post processing, I found the total charge around integration around FE slab is nonzero due to charge transfer between graphene, so the dipole moment obtained from this integration depends on origin. Because of these, I'm trying to find another way to obtain dipole moment, which is described above.

Hope it make things more clear.