UT theoretical chemistry code forum

Hi,

I used Quantum Espresso to calculate the charge density distribution for graphene (hexagonal). For the calculations four valence electrons per carbon atom have been considered, the remaining two electrons have been treated within the framework of pseudopotentials. Although I expect a symmetric charge distribution among both atoms, the Bader analysis gives a significant asymmetric result. I also found a related topic in the forum (viewtopic.php?f=1&t=562&p=2482&hilit=asymmetric#p2482), but it is rather old and gives no answer to the problem.
Quantum Espresso uses the following basis vectors for hexagonal Bravais lattices (a...lattice constant):
a1=a(1,0,0)
a2=a(-1/2,1/2*sqrt(3),0)
a3=a(0,0,c/a)
with lattice parameters
a=4.63 bohr and
c/a=6.05 (corresponds to a vacuum layer of about 28 bohr).
For graphene the basis vectors for the carbon atoms are
s1=1/3*a1+2/3*a2
s2=2/3*a1+1/3*a2.

The content of ACF.dat tells me, that there are +0.11 electrons more on one carbon atom, while -0.11 electrons are missing on the other atom. I also tries a higher FFT grid (5x more points in x and y direction than usual), but the asymmetry survives.

You can download the .cube file of the charge density from the following link, if anybody is interested in the result.
https://webmail.uni-jena.de/upload-data ... hg_pp.cube


Furthermore there seems to be a bug. When one calculates the Bader volumes around each atom using
./bader -p sel_atom 1 2 chg_pp.cube
the header of the written Bader charge densities contains wrong lattice vectors. For camparison:

------original header in chg_pp.cube:
Cubfile created from PWScf calculation
Total SCF Density
2 0.000000 0.000000 0.000000
24 0.193000 0.000000 0.000000 <---please compare these lines
24 -0.096500 0.167143 0.000000 <--- with the corr. ones given below
144 0.000000 0.000000 0.194448
6 6.000000 0.000000 2.674286 14.000220
6 6.000000 2.316000 1.337143 14.000220
....
....

------header in e.g. BvAt0001.cube
Gaussian cube file
Bader charge
2 0.000000 0.000000 0.000000
24 0.193000 -0.096500 0.000000 <---
24 0.000000 0.167143 0.000000 <---
144 0.000000 0.000000 0.194448
6 6.000000 0.000000 2.674286 14.000220
6 6.000000 2.316000 1.337143 14.000220
....
....

Thank you in advance
Lars

Thank you for the bug post. The cell coordinates for the written cube files will be fixed asap.

The asymmetry in the charge density could be more tricky. It is most likely in the charge density file -- if there is any indication that there is a problem with the analysis, we'll try to track it down. Check, for example, if the asymmetry is present for the Voronoi calculation. If so, it is a property of the charge density and not the Bader calculation.

Thank you for your reply.

The asymmetry is preserved even for Voronoi decomposition of the density.
I'm still worried about the order of magnitude of the deviations between expected charge (4 electrons per atom) and obtained charge with Bader/Voronoi analysis (3.89 and 4.11 electrons) because it is more than just numerical noise.

Concerning the Voronoi decomposition I think in my case (highly symmetric unit cell) the Voronoi cells should be of equal size, because they don't depend on the charge density, right? This is the output obtained with ./bader -c voronoi chg_pp.cube :

VORONOI ANALYSIS RESULT
# X Y Z CHARGE ATOMIC VOL
----------------------------------------------------------------------
1 0.0000 2.6741 13.9999 4.1072 269.1173
2 2.3158 1.3370 13.9999 3.8928 251.0557
-----------------------------------------------------------------------
NUMBER OF ELECTRONS: 7.99998

But why are the atomic volumes different? As I explained before, in my opinion they should be equal. By the way, the positions of both atoms are correct.
In total the cell has a volume of about 520.031 bohr^3 and the sum of both Voronoi volumes is 520.173 bohr^3. Maybe some problems with charge density at the Voronoi boundaries or at the boundaries of the periodic unit cell where charge density is not well seperated?


With best regards
Lars

It does seem fishy and worth looking into.

Perhaps look to see if the volume difference changes with grid spacing. I see that you have a course grid of 24x24xZZ in the graphene plane. Depending upon how the grid points fall between the atoms, you could get errors on the order of ~1/24 in the estimated voronoi volumes, which is not far from what you see. If you use a 48x48xZZ grid and do not see the error drop by a factor of 2, then I would agree that something is not right.

As you suggested I used a 48x48x(144) FFT grid and now I obtained the following expected (correct) results with Voronoi decomposition:
VORONOI ANALYSIS RESULT
# X Y Z CHARGE ATOMIC VOL
----------------------------------------------------------------------
1 0.0000 2.6741 13.9999 4.0003 259.4105
2 2.3158 1.3370 13.9999 3.9997 260.7652
-----------------------------------------------------------------------
NUMBER OF ELECTRONS: 8.00002

But the problem with asymmetric charges in the Bader analysis remains. First I calculated the charges for the same grid as above (48x48 points) and also for 96x96 points with the results:

BADER:
FFT-grid: 48 x 48 x 144
# X Y Z CHARGE MIN DIST ATOMIC VOL
--------------------------------------------------------------------------------
1 0.0000 2.6741 13.9999 4.0555 1.2544 260.2431
2 2.3158 1.3370 13.9999 3.9445 1.1699 259.9326
--------------------------------------------------------------------------------
VACUUM CHARGE: 0.0000
VACUUM VOLUME: 0.0000
NUMBER OF ELECTRONS: 8.0000


BADER:
FFT-grid: 96 x 96 x 144
# X Y Z CHARGE MIN DIST ATOMIC VOL
--------------------------------------------------------------------------------
1 0.0000 2.6741 13.9999 4.0279 1.2955 258.3895
2 2.3158 1.3370 13.9999 3.9721 1.2535 261.7862
--------------------------------------------------------------------------------
VACUUM CHARGE: 0.0000
VACUUM VOLUME: 0.0000
NUMBER OF ELECTRONS: 8.0000


Is the convergence with respect to the grid points in general really as slow as in this example, because for larger systems I cannot easily increase the number of points in each (2D) direction by a factor of 4 or 5 in order to get the correct charges and, in particular, to get reliable values for charge transfer complex structures.

These highly symmetric systems are probably the worst for convergence because there is no cancelation of errors along surfaces dividing the Bader volumes. The error should go like 1/N where N is the number of points along an axis.

For higher accuracy, we are putting in a weight method that was developed by Dallas Trinkle. This scheme partitions boundary points fractionally to the Bader volumes. I'll post when this is working. Here is a link to the method, including nice plots of convergence:

https://dtrinkle.matse.illinois.edu/_me ... cp2011.pdf

Thanks a lot for this information, I didn't know that in case of highly symmetric system the method scales worse then for 'noisy' systems. In the end I want to compare charge transfer of molecules adsorbed on graphene. Therefore, I need supercells of graphene containing more then 150 carbon atoms. While calculating the charge transfer I realized (by means of Bader analysis) a periodic accumulation and depletion of electrons between neighboring carbon atoms. Unfortunately I thought it was a feature of the system and not a problem of Bader analysis. Anyway, for the total charge transfer to graphene it gives reasonable values compared to a method based on projections of quantum mechanical solid state wavefunctions onto single atomic orbitals, but of course for this method no sum rule is fulfilled.

If possible, can you please inform me as soon as the new version of 'bader' with the weight method is released? I could also provide you with data sets for testing purposes if needed.

With best regards
Lars

dear admin,

any updates here? is the latest update (Version 0.27d Released) supposed to fix this problem or was this something different?

i'm asking because i always got the same results (different charges on the two Carbon atoms) and just tested the latest version and still got different charge contributions on the two Carbon atoms. (but only on a coarse grid yet).

The new version fixed the printing of the lattice in cube files.

The weight-based interpolation algorithm to help with accuracy for course grids will take a while to implement. For now, you must increase the grid density until the charges are converged.