Wigner-Seitz radii and CHGCAR

[alexabr]

Hello community,

I am calculating partial and total charges using certain set of Wigner-Seitz radii with VASP. At the same time I perform similar calculations with the same set of radii using electronic density and positions of ions provided in CHGCAR, independently of VASP. I am obtaining discrepancy in the total charge around 2 electrons.
Does VASP use the same electronic density in these calculations as it is printed in CHGCAR? If so how it is possible to subdivide density into s, p, and d on the basis of information provided in CHGCAR? Where is the place in the source code at which it would be possible to print the density which is used for charges calculations in VASP, if the density is different from CHGCAR density?

Alex

[alexabr]

Some new data about the issue. My test system is cubic boron-nitride. There are two types of atoms B and N. I varied radii for each type in the range 0.1-1.9, step was 0.1. At each step I enforced Wigner-Seitz condition, i.e. 4/3*pi*SUMi(Ri^3)=Vcell. I used these Wigner-Seitz enforced radii to calculate total charge with VASP. My findings showed that the maximum charge is recovered when boron radius is bigger than nitrogen radius. This is in contrast with physical intuition as the system is ionic and electronegativity of nitrogen is larger. This qualitative conclusion is supported by visualisation of the electronic density in CHGCAR with VMD, indeed nitrogen has much more density around.

This might possibly be numerical artefact due to double counting of electronic density in VASP calculation. I suspect that VASP double counts density in overlapping regions of Wigner-Seitz spheres as some recovered charges are bigger than the number of valence electrons in the system.

Any way it would be better to get some clarifications from developers as well as answers to my previous questions.

Alex

[admin]

please note that local charges in all PW and mixed basis set codes are obtained by integrating the charge density over the volume given in spheres (r=RWIGS) around each ion. Hence all local charges are somewhat arbitrary (some space may be accounted for twice if the spheres overlap, some 'interstital' may be left void), depending on the choice of RWIGS. This behavior is inherent to all these methods, not just vasp.

[alexabr]

Thank you for reply.
There should definitely be the way to recover maximum amount of electronic density by appropriate selection of the radii. I am doing the integration you have mentioned using information in CHGCAR and I am obtaining discrepancy with VASP results in both case when I double count charge and when I do not. Moreover the VASP results are in contrast with qualitative understanding and with visual representation of CHGCAR provided by VMD (I have pointed that out in previous post using BN system). That is the reason I have doubts is VASP doing the integration over the density provided in CHGCAR or over some other density. If latter is the case then where is the place in the source code at which it would be possible to print this density out.

Alex

[admin]

I am sorry, I fear I don't understand correctly what you mean. You cannot subdivide any polyhedron into spheres such that the volumes of the spheres sum up to the same volume as the cell , without leaving space void and counting space (at least) double. Therefore any splitting into local charges is somewhat arbitrary, even if you weight the radii by the ratio of the covalent or ionic shpere radii of the atoms, or the ratio you can estimate from CHGCAR plots...
VASP does the integration correctly, but it does not integrate over local wavefunctions, but over l,m projected charge densities, obtained from the PW to obtain the local charges. Please mind the difference!

[alexabr]

The space is double counted but the charge in the overlapping region can be divided among different ions sharing this region (in several ways actually). Even in this case the local charges are arbitrary, I agree.


---
VASP does the integration correctly, but it does not integrate over local wavefunctions, but over l,m projected charge densities, obtained from the PW to obtain the local charges. Please mind the difference!
---

Thank you for this clarification.

[alexabr]

By the way how the density elements (parallelepipeds) of CHGCAR are centred in the unit cell? The very first one has coordinates 0, 0, 0 or
LatticeVectorX / NGXF / 2; LatticeVectorY / NGYF / 2; LatticeVectorZ / NGZF / 2

Where the dimensions of each density parallelepiped are
LatticeVectorX / NGXF; LatticeVectorY / NGYF; LatticeVectorZ / NGZF

Thanks in advance.

[alexabr]

Perhaps I need to clarify what I mean. Let us take a cubic unit cell AxAxA for simplicity. Then, taking into account the number of points of FFT grid on which the charge is given, the dimensions of a small cube around each grid point are ax=A/NGXF; ay=A/NGYF; az=A/NGZF. For the cube we also should have NGXF=NGYF=NGZF and ax=ay=az=a. It is also reasonable to assume that the charge density within the small cube is constant and equal to the charge density given on the corresponding grid point (this is important for the charge integration only and is not relevant for the question).

Let us now calculate coordinates of the center of the each small cube along X axis only. The coordinate system is centered at the vertex of the cubic unit cell and X, Y, Z are directed along the lattice vectors. The very first element along X 0; the second one 0+a; the third one 0+2a; ... the very last one 0+( NGXF-1)a.

Alternatively we may think that the coordinates along X are: the very first element a/2; the second one a/2+a; the third one a/2+2a; ... the very last one a/2+( NGXF-1)a.

In the last case the charge density is somehow centered within the cubic unit cell. In the former case it is shifted in negative direction by a/2.

Hopefully it is clear. I just want to be sure that I understand CHGCAR format correctly. Could you please correct me if I am wrong or tell me which way of thinking is proper?