My Community

I have a question/sanity check about treatment of electrostatics in VASP. I think my confusion may just stem from not thinking properly about the physics, but nevertheless I want to ask about this.

I am performing calculations of a heterostructure consisting of an antiferromagnetic insulator (Cr2O3) and a heavy metal (Pt), terminated with vacuum on both sides. Due to the difference in electrostatic potential between the heavy metal and the insulator , I would expect an electric dipole to be formed at the heavy metal-insulating interface, and also separately, an electric dipole formed across the entire slab due to the difference in the delta E between vacuum and the insulating interface, and the delta E between vacuum and the heavy metal interface. I attach the INCAR and structure file here.

If we look at the LOCPOT plotted along the c-axis of the heterostructure (generated with LVHAR=.TRUE.), we do see a finite slope in the vacuum region, consistent with the second dipole source I described. Also, there is a clear "crossover region" between the insulator and heavy metal at their interface, before the potential on either side look like that of the isolated consitutents.

What I'm confused by is why the average electrostatic potential across the insulator (between about .5 and .7 of the fractional c-coordinate) appears completely flat, and does not have a slope. I guess maybe I can understand for the heavy-metal region (0.7-1) as the electric dipole might be screened here, but for the insulating compound I don't see how the dipole would be screened. Clearly, the slab has a finite electric dipole through the structure as manifested in the vacuum area. But I don't understand why this doesn't also seem to show up in the electrostatic potential profile within the slab region. As an aside, this calculation was done without dipole corrections but the same thing happens if I include them (i.e. slope in vacuum, no slope across structure ).

Does this simply have to do with the details of electrostatic boundary conditions in VASP? Or is this profile consistent with what you would expect based on physics (in which case if you could clarify that would be great.)

Thanks in advance!

Thank you for your detailed question.

I am not sure that I understood it completely.
In particular, I don't understand what is the red/blue lines in the plot.
Could you explain how you obtained them?

Hello,

Sorry for not explaining better. the red lines are the "average" electrostatic potential (where the electrostatic potential is plotted in dark blue and is simply the output of the locpot for LVHAR=TRUE) where the average is separately taken in the different regions divided by the vertical blue lines. basically I take the average in the region of the insulator in the heterostructure, and in the matallic region of the heterostructure (and I separate out the interface region where the electrostatic potential clearly transitions).

The point I am struggling with is why there is not a linear slope superimposed on the periodic electrostatic potential variation due to atoms in the slab, throughout the who slab, rather than just in the vacuum region. Clearly their is an electric dipole generated at the interface so why doesn't this manifest throughout the entire slab?

Best,
Sophie

Ok, thank you for the clarification!
I think I understand your question now: because the Cr O compound is an insulator, you would expect to see a linear slope superimposed in the periodic potential, indicating a finite electric field through the slab similarly to what one sees in the vaccum region.
Unfortunately, I don't know how to explain it either.

The first explanation that comes to my mind would be that the layer structure of the Cr O compound is not an insulator in this heterostructure, but I don't know.
I guess you could try plotting the bandstructure projected onto the atomic sites of Cr O and see how they are located w.r.t the fermi level.
One way to do it would be using py4vasp reading the vaspout.h5 file:

import py4vasp
c = py4vasp.Calculation.from_path('.')
graph = c.band.to_graph('Cr O')
graph.show()

Hi,

Thanks a lot for your response. I haven't calculated the band structure, but the same information should show up in the layer-projected density of states, which I have. It is indeed the case that the layer of the insulator directory interfaced with the metal has a very small finite dos at the fermi level (dashed black line, first slide of pdf). But the layer dos becomes insulating one layer below and is insulating for all further layers until vacuum. Do you really think that the single metallic layer at the insulator-metal interface would be enough to completely screen the electric dipole? Also, if this was the case, why would there be any finite slope in the vacuum?

As a sanity check I looked at a calculation where I used an applied electric field (EFIELD) on a completely insulating structure, and as we both expected, there is a linear slope superimposed on the bare crystal electrostatic potential (see thin black lines for guidance on slide 2), which has about the same slope as the vacuum. So indeed the "normal" behavior for an insulator is that a slope should persist throughout the entire structure. But in slide 3, it's quite clear that the average potential in the insulating region is really flat.

I have an additional confusion about a seeming descrepancy between the LOCPOT and the layer-projected density of states on slide 1. For the layer dos, if you look at "layer 5" (second layer of the insulator away from the insulator-metal interface), the highest occupied states are below the fermi level of the metal. And clearly you get an upward band bending as you go from the insulator to the metal. I would interpret this to mean that your effective electric field at the interface (with the normal convention where electric fields point in the direction of motion of POSITIVE charges), which points from the insulator (lower electron energy/higher positive charge energy) to the metal (higher e energy /lower positive charge energy).

However if you instead look at the LOCPOT, the electrostatic potential of the metal, for electrons, is LOWER than for the insulator. This would imply an electric field with the same convention would point in the opposite direction, from the metal to the insulator. I would think that you ought to be able to determine the sign of the electric dipole/effective electric field at the interface independently with either of these two methods. But they seem to give opposite results. Can you tell me what I am missing?

Thanks very much!

It is hard for me to answer to this question, "Do you really think that the single metallic layer at the insulator-metal interface would be enough to completely screen the electric dipole?" without performing additional calculations.
I think you have more experience than me in this type of calculations, and you are very careful with your analysis.

One thing that you could maybe try is to increase the number of layers of the insulator and check if with sufficient layers you get a slope in the potential as you expect.

The projected density of states as you plotted in the slides is perfect to check whether you still have an insulator!
From the looks of it, I would say no. The fact that you have some density of states at the fermi level means that your insulator is not so much an insulator anymore, I think.
Hence, my suggestion to try checking this dependence with more layers of insulator or maybe a different type of insulator?
Perhaps with increasing number of layers, the insulator is able to screen the top from the bottom layer and retain its insulating character.
But I am not sure about this, sorry.

In the upcoming release of VASP you will be able to use Coulomb truncation which might be helpful for these calculations, but it is really hard for me to say if it will help to answer your questions without performing tests and looking at the data.