UT theoretical chemistry code forum

Dear all,

I am trying to compute the Bader charges of cubic LaMnO3 from VASP files. However, I am having a problem with the convergence of the charges versus numbers of grid points NG(X,Y,Z)F. When I compute the charges, I see a slight drift in the charges versus NG(X,Y,Z)F for large NG(X,Y,Z)F. Thus, I am hoping that you can suggest how to eliminate the drift.

Representing cubic LaMnO3 using the model that is given at the end of this post, I use the procedure that is given at the HenkelmanGroup website to compute the charges of the La, Mn, and O atoms (https://theory.cm.utexas.edu/henkelman/code/bader/). Specifically, I generate AECCAR0 and AECCAR2 files, add them using the chgsum.pl script, and then run bader CHGCAR -ref CHGCAR_sum. To generate the VASP files, I use a stringent wave function convergence criterion, EDIFF = 1e-6.

The charges of the La, Mn, and O atoms are given below as average charges qLa, qMn, and qO +/- sigmas for different grids. The different grids are specified by values of NG(X,Y,Z)F that are 1x, 2x, 3x, 4x, and 5x greater than the default values given by VASP and are labeled “original”, “2x greater”, “3x greater”, “4x greater”, and “5x greater”, respectively:

grid………..….qLa…..………….....….…qMn…….…....……………qO…….…....……………
original…….…2.0653 +/- 0.0003..….1.6211 +/- 0.0005…...-1.2288 +/- 0.0034
2x greater……2.0780 +/- 0.0002...…1.6453 +/- 0.0007.….-1.2411 +/- 0.0020
3x greater……2.0816 +/- 0……….....1.6542 +/- 0……..…....-1.2453 +/- 0.0030
4x greater……2.0720 +/- 0.0001...…1.6591 +/- 0.0001..….-1.2437 +/- 0.0017
5x greater……2.0756 +/- 0……...…..1.6623 +/- 0.0001…...-1.2460 +/- 0.0002

Examining the charges given above, qLa appears to level off versus grid. However, qMn changes from 1.6211 (original) to 1.6453 (2x greater) and then drifts from 1.6453 (2x greater) to 1.6623 (5x greater), while qO also seems to drift.

In addition to the above charges that are obtained using the default method (-b neargrid), the charges that are obtained using the -b weight flag are given below:

grid………..….qLa…..………….....….…qMn…….…....……………qO…….…....……………
original…….…2.0727 +/- 0……….....1.6394 +/- 0.0002…...-1.2374 +/- 0.0001
2x greater……2.0754 +/- 0……….....1.6524 +/- 0.0001…...-1.2426 +/- 0.0001
3x greater……2.0759 +/- 0……….....1.6581 +/- 0……..…....-1.2447 +/- 0.0001
4x greater……2.0762 +/- 0……….....1.6614 +/- 0.0001…...-1.2459 +/- 0.0001
5x greater……2.0764 +/- 0……….....1.6636 +/- 0……..…....-1.2467 +/- 0.0001

Examining these charges, the charges obtained using “original” and “2x greater” grids are better converged than the corresponding values in the first table. However, qLa, qMn, and qO are all drifting for grids of “2x greater” or more. For example, qMn changes from 1.6524 (2x greater) to 1.6636 (5x greater).

I would appreciate any suggestions or help to eliminate the drift in the charges versus grid for grids of “2x greater” or more and resolve the problem of converging the charges versus NG(X,Y,Z)F.

Thanks,

Yves

VASP POSCAR file of LaMnO3 model
---------------------------------------------------

LaMnO3 model - atoms are given in the order La, Mn, and O, and atoms of a certain type are symmetrically equivalent
1.00
6.75153404790348 0.00000000000000 0.00000000000000
0.00000000000000 5.51260446613033 0.00000000000000
0.00000000000000 0.00000000000000 9.54811101736883
6 6 18
Cartesian
0.00000000000000 2.75630223306517 0.00000000000000
4.50102269860232 0.00000000000000 1.59135183622814
2.25051134930116 2.75630223306517 3.18270367245627
0.00000000000000 0.00000000000000 4.77405550868442
4.50102269860232 2.75630223306517 6.36540734491256
2.25051134930116 0.00000000000000 7.95675918114069
3.37576702395174 2.75630223306517 0.00000000000000
1.12525567465058 0.00000000000000 1.59135183622814
5.62627837325290 2.75630223306517 3.18270367245627
3.37576702395174 0.00000000000000 4.77405550868442
1.12525567465058 2.75630223306517 6.36540734491256
5.62627837325290 0.00000000000000 7.95675918114069
0.00000000000000 0.00000000000000 0.00000000000000
2.25051134930116 1.37815111653258 0.79567591811407
2.25051134930116 4.13445334959775 0.79567591811407
4.50102269860232 2.75630223306517 1.59135183622814
0.00000000000000 1.37815111653258 2.38702775434221
0.00000000000000 4.13445334959775 2.38702775434221
2.25051134930116 0.00000000000000 3.18270367245627
4.50102269860232 1.37815111653258 3.97837959057035
4.50102269860232 4.13445334959775 3.97837959057035
0.00000000000000 2.75630223306517 4.77405550868442
2.25051134930116 1.37815111653258 5.56973142679848
2.25051134930116 4.13445334959775 5.56973142679848
4.50102269860232 0.00000000000000 6.36540734491256
0.00000000000000 1.37815111653258 7.16108326302662
0.00000000000000 4.13445334959775 7.16108326302662
2.25051134930116 2.75630223306517 7.95675918114069
4.50102269860232 1.37815111653258 8.75243509925476
4.50102269860232 4.13445334959775 8.75243509925476

Looking at your data, I think that the convergence rates are what you should expect from the methods. I focused on the -b weight method, and in the attached figure (below) I've plotted Log error vs Log grid density (points) and a line with a slope of -2. The agreement indicates that the charges are converging quadratically with the grid density, as expected for the -b weight method.

I do not understand why the charges systematically over- or under-estimate for La/Mn and O respectively but I don't think it is fair to describe it as drift because of the quadratic convergence. It's also not clear to me that this is a property of the Bader algorithm as opposed to being a property of the charge density data from vasp. I think that one would need to take the fine grid and average points to make a coarser grid to distinguish between these cases. Either way, the convergence rate is expected.

Dear Professor Henkelman,

Thanks very much for your reply. I’d like to ask a few follow-up questions:

First, I’d like to understand better what you’re plotting. You say that you’re plotting log error versus log grid density, but I’m not sure how you’re determining log error and log grid density. I know that you’re not looking at the qMn obtained using the -b flag and plotting log [qMn(“5x greater” grid) - qMn] versus log (relative number of grid points), for then I would expect to obtain data points with y-coordinates of log (1.6636 - 1.6394), log (1.6636 - 1.6524), log (1.6636 - 1.6581), and log (1.6636 - 1.6614) and x-coordinates of log (1), log (2^3), log (3^3), and log (4^4), respectively. However, the data points in the plot do not have these x- and y-coordinates. Would you please explain what you’re plotting in more detail?

Also, you say, “I do not understand why the charges systematically over- or under-estimate for La/Mn and O respectively but I don't think it is fair to describe it as drift because of the quadratic convergence.” Indeed, I thought that I might not have been obtaining accurate charges of La, Mn, and O, because the charges of La and Mn were slowly increasing or becoming more positive versus grid while the charge of O was slowly decreasing or becoming more negative versus grid. However, you’re saying that the charges are not incorrect due to the slow changes of the charges versus grid, because they are converging quadratically with the grid density, as expected for the -b weight method. Am I correct?

The other issue that I believe you’re raising is to explain why the charges are changing in the directions that they are changing, that is, the charges of La and Mn are systematically becoming more positive while the charge of O is systematically becoming more negative versus grid as opposed to the charges of La and Mn systematically becoming more negative and the charge of O systematically becoming more positive versus grid. You say, “It’s… not clear to me that this is a property of the Bader algorithm as opposed to being a property of the charge density data from vasp.” To explain why the charges are changing in the directions that they are changing, are you suggesting to take a CHGCAR file with a fine grid, average points to make a CHGCAR file with the same number of points and located in similar positions as in another CHGCAR file, and then compare the charges obtained using the two CHGCAR files? Would you please explain how it would help?

Thanks,

Yves

For the plot, I was looking at the -b weight La data, so the y-axis values are ln(-(2.0653-2.0756)), ln(-(2.0780-2.0756)), .. and the x-axis values are ln(1), ln(2), ... The -b weight algorithm should converge quadratically with the spacing between points, not the total number of points.

I am not sure why the charges are changing monotonically with grid density. I would have expected the error to be more statistical in nature. However, while I don't understand that, I see quadratic convergence and that is all we can expect from the algorithm.

My final small point was just a technical one, and that is that these DFT calculations of different CHGCAR files with different grids do not nessecarily represent the exact same underlying function of charge density. We are basically assuming that the only difference between the charge densities is the grid size, but the DFT calculation and the underlying FFTs are also being done on those different grids. So I was saying that (only in principle) one could take the finest grid as 'exact' and then generate courser grids based upon averaging the finer grid. That would be a fair test for the convergence of the weight method with grid size using the same underlying data. For that, I would be very surprised to systematic changes in the charges with grid size; rather, I would expect random statistical errors that drop quadratically with the grid spacing.

Another way of saying this is that if you consider some other DFT parameter, such as ENCUT, you might well find monotonic convergence of Bader charges with a fixed grid. In the same way, increasing the grid density is also increasing the precision of the calculation which could systematically change the Bader charges.

Dear Professor Henkelman,

Thanks for your additional comments. Now I understand the points that you were making in your previous post.

Following the procedure that you used to examine the -b weight La data, I plotted the -b weight La, Mn, and O data and fitted linear functions to the data points. The functions had slopes of -2.03 +/- 0.17, -1.65 +/- 0.28 and -1.69 +/- 0.27, respectively. Because the slopes are all about equal to -2 within error, I see why you’re saying that the charges are quadratically convergent versus grid spacing.

A few more questions:

To determine the rate of convergence for the charges, I believe that the values of NG(X,Y,Z)F need to be increased by the same constant factor when recomputing the charges using new values of NG(X,Y,Z)F. For example, I increased NG(X,Y,Z)F by a factor of two to determine the charges labeled “2x greater”. I didn’t determine a new set of charges by, for example, increasing NGXF and NGYF by one factor and NGZF by a different factor. The NG(X,Y,Z)F are simply the numbers of grid points along the different lattice vectors. Thus, I assume that the rate of convergence for the charges in a non-orthorhombic cell could be determined by plotting ln (y) versus ln (x) for each of the charges as I did above, as long as the values of NG(X,Y,Z)F are changed by the same constant factor when recomputing the charges using new values of NG(X,Y,Z)F?

Also, as you said, the charges are changing monotonically with grid spacing and not exhibiting a random error that drops quadratically with grid spacing. You suggested that the monotonic change of a charge could be due to the fact that the charge densities obtained using CHGCAR files with different grid densities will be slightly different. This explanation is not specific to the system being examined here, cubic LaMnO3. Thus, it would seem that the charges of other systems will change monotonically versus grid density when they are computed versus grid density using VASP. In fact, in the paper describing the -b weight method [M. Yu and D. R. Trinkle, Accurate and efficient algorithm for Bader charge integration, J. Chem. Phys. 134, 064111 (2011)], the charge of Na is shown to decrease monotonically versus grid density (Fig. 7). Do you believe that the monotonic change of a charge versus grid density could be the norm rather than the exception when the charge of a system is computed versus grid density using VASP? Alternatively, do you not expect to see a monotonic change of a charge versus grid density for all systems because the changes in charge density as the grid is changed are not necessarily going to result in the monotonic change of a charge?

Lastly, because the charges are changing monotonically versus grid density, I’d like to identify the best way to determine the final values of the charges. One way would be to plot the value of a charge versus relative number of grid points along a lattice vector and fit the function f(x) = a - b/x, where a and b are adjustable constants, to the data points to determine the charge’s asymptotic limit, equal to a. Although the value of a will depend on the number of fitted data points, I think that the error will be relatively small. For example, if I plot qMn of 1.6394, 1.6524, …, 1.6636 at values of 1, 2, …, 5 and fit f(x) to the data points #1-5, #2-5, #3-5, and #4-5, I obtain a = 1.6685, 1.6707, 1.6717, and 1.6717, respectively. Thus, the value of qMn could be taken as 1.67e. Do you agree that this procedure is a reasonable way to determine the final values of the charges? Do you have a better way to determine the final values of the charges?

Thanks,

Yves

Yves,
I think I understand and agree with everything that you say. Unfortunately, I have very little to add to answer your questions. I am surprised by the monotonic change in charges with grid density and I don't know if that would be typical of other systems. Following that, I can also not suggest if an extrapolation to an infinite grid density makes sense, or how you would do that. Sorry that I'm not more help.
Graeme

Dear Professor Henkelman,

Thanks for your reply. The monotonic change of the charges versus grid density is odd. However, even if it’s attributed to a systematic error in the determination of the charges such as the error that you suggested, that is, slightly different charge densities in CHGCAR files with different grids, at least the charges can be accurately extrapolated to an infinite grid density to obtain final values of the charges. Notably, the charges obtained with reasonable grid densities of “2x greater” or “3x greater” are close to the extrapolated values. Similarly, the qMn and qO obtained using the -b neargrid method can be extrapolated to an infinite grid density.

Thanks again for your comment about the charges obtained using the -b weight method following the expected convergence behavior and a possible explanation for the monotonic change of the charges versus grid density.

Yves