For surfaces that have a non-vanishing dipole density quantum-espresso (QE) offers a dipole correction in combination with the sawtooth potential. The relevant variables in the pw.x input are dipfield and tefield, which both have to be set to true. Unfortunately this is not all there is to do – we also have to define where the dipole is located and the location (and size) of the dipole play an important role. Since this confuses me regularly I decided to write it down here in the hope it helps me and others.
Lets assume we have centered the slab in the center of our cell and the cell extends from 0 (the bottom) to 1 (the top). As QE uses periodic boundary conditions everything that exceeds the length of the cell will re-enter the cell at the bottom, i.e. (1+x)=x in relative z-coordinates. As the manual tells us one has to take care of 3 parameters: emaxpos, eopreg and eamp. Since we only want to apply a dipole correction the amplitude eamp=0. The field linearly increases until emaxpos and then decreases to reach 0 at emaxpos+eopreg, as such they represent the positive and negative charges that create the dipole. It is important that the change of slope, that is emaxpos, lies in the vacuum region and so must emaxpos+eopreg (the opposite charge) – otherwise the correction will not work. A very useful setup is to center the slab in the cell (~0.5 in relative z-coordinates) and the dipole around 0 of the slab. One can then slowly increase the length of the dipole and see if the resulting electrostatic potential (plot_num=11 in pp.x) becomes flat in the vacuum region (thanks for Dr. Thomas Brumme for the hint!). As long as the dipole and the charge density of the electrons do not overlap the correction works and the necessary amplitude of the dipole field decreases linearly with the dipole length – as to be expected since p=q*z (see the right side of Figure 1). When the dipole gets too close the correction fails and the potential is no longer flat in the vacuum region. Please note that this does not mean that SCF convergence cannot be reached – in fact, the last case converges faster than all the cases were the results are physically meaningful in this particular case.
Figure 1. Potential as obtained by pp.x (plot_num=11) and averaged in xy-plane (average.x) as a function of the dipole length with the dipole centered at 0 of the cell. On the right side the amplitude of the dipole field is shown as function of the dipole cell.
In recent versions of QE a different correction has been implemented by Sohier et al. as well. This is activated with the input-flag assume_isolated =’2d’ and gives basically the same result. As an advantage the new correction seems to be working better with phonons of 2-dimensional systems but not all codes that build on QE output are able to read it yet (for example the resulting electron-phonon matrix elements in EPW are too large by a factor of at least 100).
Attached shell script will run a calculation for dipoles of varying length – by collecting the results one can reproduce above figure.
Hope this was helpful – if not, let me know in the comments!
Christoph Wolf 
Hi Christoph, thanks for this useful information. But I am still not fully clear. If the bottom of my slab is at the bottom of the cell, and the top of my slab is near 0.3 of the cell z direction, the electrostatic potential without dipol correction clearly increases from 0.35 to 0.95 which actually is the vacuum region. Then for this case, how should I set these two parameters of emaxpos and eopreg?
LikeLike
Dear Lesheng,
Sorry for the late reply! In your case you want to locate the dipole probably at emaxpos =0.6 and eopreg=0.05. The “physical region” is then from 0.6+0.05-1=-0.35 to 0.6 of your slab and contains the entire vacuum. Post-process the potential and sawtooth potential (plot_num=11 and 12, respectively) to see if the dipole is located where you would expect it to be!
LikeLike
.
Dear Christoph,
I am having a difficulty implementing your explanation on dipole correction. I have a slab centered in the supercell. The atoms are distributed between 0.21 and 0.8 fraction of the slab. Please explain to me how to chose emaxpos and eopreg in this case.
Regards.
Timothy Chibueze
University of Nigeria
LikeLike
Hi Timothy,
Your arrangement seems very tight with respect to vacuum; are you sure you put enough vacuum space between the slabs (~15-20 Angstrom is usually a good choice); in any case you could chose emaxpos=0.95 and eopreg=0.05 (then the sawtooth potential is continuous across the region from 0 to 0.95). Don’t forget to plot the potential (num_plot=12 in pp.x) to check that is is situated as you would expect!
HTH!
Chris
LikeLike
first of all, I am very grateful for your interesting information and tools. what about strain calculations. did you have any script already done to automate calculations? what about HSE06 in QE? thank you very much. Best regards
LikeLike
Hi kika,
I would avoid employing hybrid functionals on QE since they do not lead you far within this package (you can’t plot band structures or run post-processing with HSE, for example). For strain calculations you could simply edit Christoph’s script to write the new respective lattice parameters on each scf input.
LikeLike
Hello Christoph,
What would one do if they have vacuum along 2 directions? Is it possible to do dipole corrections along both these directions simultaneously? what should be done in such cases?
LikeLike
Dear Lekshmi,
that does not work
and I am not sure how that would be set up considering periodic boundary conditions. Could you give me an example?
LikeLike
I have a question about metallic systems. If one applies electric field (tefield = true and eamp > 0) to a metallic slab, then the negative and positive charges will accumulate on top and on the bottom of the slab. This will create a Coulomb interaction between neighboring cells. I mean, the negative charges accumulated (lets say) on top of the system within the original supercell will interact with the positive charges accumulated at the bottom of the system within the next image of the supercell (so the charges at each side of the vacuum region). Is there a way to avoid this Coulombic interaction, other than enlarging the vacuum region as much as possible to minimize it?
LikeLike
hi
is that like to monolayer . when applying electric field . for monolayer do I most us dipole correction
LikeLike
Which software did you use to plot such a beautiful graph? Do you have the plotting script available?
LikeLike
I think this was just plotted in gnuplot (which is free and open source); for a publication grade plot I would recommend matplotlib these days, even though I am a big fan of gnuplot 🙂
http://www.gnuplot.info/
LikeLike