Problem with dipole moment correction for a slab

[ymantz]

Dear Support,

Recently, I computed the energy of an unrelaxed 12-layer slab or surface model for the LaMnO3 (001) surface as a function of vacuum region size, without and with a dipole moment correction. I believe that the results with the dipole moment correction are not correct, and I would like your help to obtain the correct results.

The surface model can be represented as (LaO-MnO_2)_6. Thus, it has a net dipole moment along the surface normal. The energy of the surface model should converge very quickly as a function of vacuum region size with the dipole moment correction (https://www.vasp.at/wiki/index.php/Mono ... orrections). However, it actually converges more slowly with the correction at vacuum regions below 42 Å and at a similar rate at vacuum regions above 42 Å (see the attached .jpg file).

According to the VASP manual, for a surface model oriented perpendicularly to the third lattice vector, the dipole moment correction is made by setting IDIPOL = 3
(https://www.vasp.at/wiki/index.php/Mono ... orrections). I believe that my computational setup is correct, and the fact that the energy of the 12-layer surface model is slow to converge with the dipole moment correction is puzzling.

The VASP input and output files, the .jpg file of results referred to above, and a file READ_ME providing a description of all the files are included in the attached .zip file.
The slowly converging energy of the surface model with the dipole moment correction is a technical issue that is not covered in the VASP manual and, I believe,
is not addressed in previous VASP forum posts. Thus, I am asking the VASP team to help me to resolve this issue. Thanks for considering this issue, and I’ll wait to hear from you.

Sincerely,

Yves

[henrique_miranda]

Dear Yves,

The system you are studying is rather complicated.
There are many things that can go wrong not necessarily related to the implementation of the dipole corrections in VASP.
Here are a couple of suggestions that might help you to trace down the problem:
1. Check if the dipole corrections work as you expect on a simpler system ex: ZnO
2. Try to compute a thinner layer of your system
3. Try to visualize the potential LOCPOT file (using LVTOT=.TRUE. in the input file)

Here a couple of references for further reading with a description of the problems involved in these calculations:
https://doi.org/10.1088/0953-8984/12/31/201
https://doi.org/10.1103/PhysRevB.68.245409
https://doi.org/10.1088/0022-3719/12/22/036

Do not hesitate in sharing your findings here in case you find something definitely wrong in the implementation of the dipole corrections as done in VASP.

[ymantz]

Dear Support,

Thanks for your detailed reply. I would like to ask a few follow-up questions:

First, I want to make sure that we agree that the problem I’m describing is a problem. In my previous post, I was asking for help in determining why VASP is giving the “wrong” energies of the 12-layer surface model for different vacuum region sizes when the dipole moment correction option is used (IDIPOL = 3). I claim that the energies are wrong, because the energy of the surface model is slow to converge as a function of vacuum region size with the dipole moment correction. In the VASP manual, it’s suggested, but not stated explicitly, that the energy of a surface model with a net dipole moment along the surface normal should converge very quickly as a function of vacuum region size with the dipole moment correction
(https://www.vasp.at/wiki/index.php/Mono ... orrections). Do you agree that the energy of a surface model with a net dipole moment along the surface normal should converge very quickly as a function of vacuum region size with the dipole moment correction? Also, referring to the plot that I attached to my original post and the energies of the surface model without and with the dipole moment correction as a function of vacuum region size (red and blue curves, respectively), do you agree that the energy of the surface model is converging too slowly as a function of vacuum region size with the dipole moment correction?

Also, I would like to ask about your approach or suggestions #1-3 to address the problem of the slowly converging energy. Because the energy of the 12-layer surface model is slow to converge as a function of vacuum region size with the dipole moment correction, I believe that the dipole corrections to the energies of the surface model are wrong. If they are wrong for the 12-layer surface model, then they are likely to be wrong for a thinner surface model or a polar surface model of a simpler material such as ZnO. Don’t you agree? Also, if the energies of a thinner surface model or a polar surface model of a simpler material were found to converge slowly as a function of vacuum region size with the dipole moment correction, what would be the next step to address the problem? Lastly, what is the reason to look at the potentials of the 12-layer surface model for different vacuum region sizes using LVTOT = .TRUE. when the energies without dipole corrections are fine (red curve in plot)?

On a different note, what do you think of the idea of someone modifying the VASP code to print the intermediate steps of the procedure to obtain the dipole correction to the energy to determine what part or parts of the procedure are slow to converge as a function of vacuum region size? Could the intermediate steps of the procedure be printed or examined in a practical way, and do you think that printing or examining them would reveal why the energy is slow to converge and how to address the problem?

Thanks,

Yves

[henrique_miranda]

Yes, it is a problem of course.
However, at this point, we only know that for this material the energies converge slowly with vacuum size, we don't know why.
It might be an issue in how the dipole corrections are implemented in VASP but it might also not be.
We need to find out.

I don't know your level of experience with using the dipole corrections in VASP that is why I suggested that you perform some calculations in ZnO. Notice that one of the references deals with ZnO so you can check if in that case, you obtain what you would expect. But you can ignore this advice in case you are already experienced with using the dipole corrections in VASP.

Because the energy of the 12-layer surface model is slow to converge as a function of vacuum region size with the dipole moment correction, I believe that the dipole corrections to the energies of the surface model are wrong. If they are wrong for the 12-layer surface model, then they are likely to be wrong for a thinner surface model or a polar surface model of a simpler material such as ZnO. Don’t you agree?

Yes, but what if the energies converge quickly with vacuum size for a ZnO or a thinner layer of your system? Then you might come to the conclusion that the corrections are well implemented but have some limitations of applicability.
In one of the references, I sent you there a discussion related to the dependence on the thickness of the slab and charge transfer.

The dipole corrections change the potential so visualizing it might help you in understanding what the problem is.
Of course, it might also just be completely inconclusive.
You can already have access to more information related to the dipole corrections, just look for the DIPCOR tag in your OUTCAR file.

[ymantz]

Dear Support,

Thanks for the additional feedback. Based on one of your comments, I decided to look more closely at the information associated with the DIPCOR tags in the OUTCAR files for the 12-layer surface model of the LaMnO3 (001) surface. Surprisingly, I found that the dipole moment of the surface model falls off as a function of vacuum region size, from 0.93 e x A at a vacuum region size of 13.9 A to 0.36 e x A at a vacuum region size of 65.9 A (see the attached plot).

Thus, the question that I think needs to be addressed now is, why is the dipole moment of the system changing as a function of vacuum region size? The dipole moment is defined with respect to Rcenter in the VASP manual and could depend on the choice of Rcenter if the charge density does not fall to zero at some distance from Rcenter
(https://www.vasp.at/wiki/index.php/DIPOL). In my computations, the DIPOL tag is not set, and Rcenter is determined by VASP as described in the manual, that is, where the charge density averaged over one plane drops to a minimum plus half a lattice vector perpendicular to the plane where the charge density has a minimum. Thus, Rcenter will change as a function of vacuum region size. Do you have any idea why the dipole moment of the 12-layer surface model is changing as a function of vacuum region size? If not, how do you suggest proceeding to resolve this issue?

Thanks,

Yves

[henrique_miranda]

Did you try setting the DIPOL tag to a fixed position in the center of your slab that does not change with the vacuum?
Also, you can check what are the coordinates for DIPOL that VASP determines in the xml file.
Look for "electronic dipolcorrection" and "DIPOL".
Hope this helps.

[ymantz]

Dear Support,

Thanks for your most recent post and suggestions. I have very good news. I can now explain the slowly converging energy and dipole moment of the 12-layer surface model.

After receiving your post, I decided to look at a simpler system, a two-dimensional layer of H2O molecules. This system can be modeled as a single H2O molecule in a box with the dipole moment of the molecule oriented along the z-axis. It was examined previously (Phys. Rev. B 59, 12301). First, I computed the energy of the molecule as a function of vacuum region size and found a slow convergence from below, as expected. Next, I performed computations with IDIPOL = 3. I found that the energy converged faster, but still too slowly, from above. Also, the dipole moment of the system converged slowly. After obtaining these results, I came to the conclusion that the results were wrong because the electron density of the system was not accurate. Because the potential of the system is forced to be periodic, both sides of the layer or model are effectively acted on by an electric field, and the electrons are perturbed by this field. Stated differently, because of the dipole-dipole coupling between periodic images, the electron density of the system is polarized. Thus, the absolute dipole moment is overestimated as well as the dipole correction to the energy, resulting in a convergence of the energy from above. A dipole correction to the potential is needed to produce a step in the potential and remove the fictitious electric field. In fact, using IDIPOL = 3, LDIPOL = .True., and setting DIPOL to be equal to the center of mass of the H2O molecule, I found that the energy and dipole moment of the system converged very fast. Also, by calculating the x,y-averaged potential from the LOCPOT file obtained by setting LVHAR = .True., I saw that a step in the potential had been introduced. Please see the attached file with plots for the H2O layer.

Next, I returned to the LaMnO3 (001) 12-layer surface model. Upon applying a dipole moment correction with IDIPOL = 3, LDIPOL = .True., and setting DIPOL to be equal to the center of mass of the surface model, I found that the energy and dipole moment of the system converged essentially right away starting at a vacuum region of 13.9 Å. Also, by setting LVHAR = .True., I saw that a step in the potential had been introduced. Please see the attached file with plots for the 12-layer surface model.

Thus, the problem was that the flag LDIPOL = .True. needed to be set. To conclude, I would like to summarize my understanding of when to use IDIPOL versus IDIPOL, LDIPOL, and DIPOL. As stated on the VASP wiki page Monopole, Dipole, and Quadrupole corrections
(https://www.vasp.at/wiki/index.php/Mono ... orrections), the flag IDIPOL should be used for slab calculations. Preferably, it should be used with LDIPOL and DIPOL, even for slabs whose atoms are fixed. Consider a slab with fixed atoms and a large vacuum region. If LDIPOL = .False., then the electron density that is perturbed by the interaction with the periodic images is used to calculate the dipole moment of the slab and correction to the energy. However, if LDIPOL = .True., then a correction to the potential is added, and the electron density is solved for self-consistently to obtain a new density that is not perturbed and is used to obtain a more accurate dipole moment and correction to the energy. Also, I know that the paper by Neugebauer and Scheffler is cited on the page Monopole, Dipole, and Quadrupole corrections (Phys. Rev. B 46, 16067). However, an additional paper that is helpful because of the comparison that is made between dipole moment corrections without and with the correction to the potential is Phys. Rev. B 59, 12301 [see Eqs. (9) and (13) and discussion at the end of Section II.B]. Do you agree with the above summary of when to use IDIPOL versus IDIPOL, LDIPOL, and DIPOL?

Also, the VASP wiki page Monopole, Dipole, and Quadrupole corrections could be changed to make it more clear when to use IDIPOL versus IDIPOL, LDIPOL, and DIPOL. What do you think?

Lastly, I would like to point out that a typo needs to be corrected on the VASP wiki page LOCPOT (https://www.vasp.at/wiki/index.php/LOCPOT):
Depending on LVHAR the total local Kohn-Sham potential (LVHAR=.FALSE.) or the sum of the Hartree and ionic potential (LVHAR=.FALSE. ---> .TRUE.) is written to the LOCPOT file. Would you please make this change?

Thanks,

Yves

[henrique_miranda]

Dear Yves,

What an adventure
I am happy that you figured out what the problem was finally.
And thank you for taking the time to explain your conclusions.
I believe this might be useful for other users.
I agree with your summary of the use of these variables.

We will go through the documentation in
wiki/index.php/Monopole_Dipole_and_Quad ... orrections
once again with your personal experience in mind and check if we can improve it.

Thanks for pointing out the mistake in our documentation.
I did some overhauling in the pages referring to LOCPOT:
wiki/index.php/LOCPOT
wiki/index.php/LVHAR
wiki/index.php/LVTOT

[ymantz]

Dear Support,

Thanks for your reply and letting me know that you’ll review the page Monopole, Dipole, and Quadrupole corrections
(wiki/index.php/Monopole_Dipole_and_Quad ... orrections) and decide whether to change it. My suggestion would be not to edit the page. Instead, I would add a note at the end of the page, such as the following:

Note for slab calculations:
As stated above, the flag IDIPOL should be used for slab calculations. Preferably, it should be used in conjunction with LDIPOL and DIPOL, even for a slab with fixed atoms. Consider a slab with fixed atoms, oriented perpendicularly to the z-axis, and with a dipole moment along the z-axis. Because of the dipole-dipole coupling between periodic images, the energy of the system is lowered and will converge from below as the size of the vacuum region is increased. Thus, a dipole correction to the energy is needed. If this correction is made with IDIPOL = 3 only, then the electron density that is perturbed or polarized by the interaction with the periodic images is used to calculate the dipole moment of the slab and dipole correction to the energy. However, if the dipole correction is made with IDIPOL = 3 and LDIPOL = .True., then a correction to the potential is added to cancel the dipole-dipole coupling, and the electron density is solved for self-consistently to obtain a new density that is not perturbed. This density is then used to obtain a more accurate dipole moment and correction to the energy, and the energy of the slab will converge rapidly with increasing vacuum region size. In VASP, the dipole corrections to the energy without and with the correction to the potential are similar / identical to those in Phys. Rev. B 59, 12301, Eqs. (9) and (13), respectively.

If you keep the last sentence of the above note, you need to indicate whether the corrections to the energy are similar or identical to those in the cited paper.

Would you please let me know what you decide to do about the wiki page by replying to this post?

Thanks,

Yves

[ymantz]

Dear Support,

I’m writing to ask whether you decided to modify the page Monopole, Dipole, and Quadrupole corrections (https://www.vasp.at/wiki/index.php/Mono ... orrections) and, more specifically, add a note to the page such as the one that I suggested in my last post. Because the page is unchanged since my last post, I think that you decided not to modify it. Am I correct? If I am, could you explain why you decided not to modify it? Is it because the forum is the preferred way to resolve issues such as how to perform the dipole correction for a slab calculation?

Thanks,

Yves

[henrique_miranda]

Dear Yves,

We still did not find an occasion to discuss internally in detail how and if to modify this wiki page.
I assure you that this issue is not forgotten and we will try to make this wiki page clearer.
We really appreciate your suggestion of modification and will consider using it

The forum is the preferred way to answer user questions.
But we try to improve the documentation whenever possible.

Thanks,
Henrique Miranda

[ymantz]

Dear Henrique,

Thanks for letting me know that the issue of whether and how to modify the wiki page still needs to be addressed. I’m glad to hear that the issue’s considered to be an important one to address and that my suggestion to improve the page is helpful.

I’ll wait for the issue to be addressed and to hear from you or someone whether Support decided to modify the page and, if they decided not to modify it, why they decided not to modify it.

Thanks,

Yves