Electron Density as a Function of Position from CHGCAR

[JonM]

I understand that this may be a super rudimentary question but how can I get a grid of the electron density from a CHGCAR file? Basically, I would like to know (x, y, z, rho) for the system (to within the accuracy of the calculation).

Looking at the CHGCAR, I can clearly see where it gives me the lattice vectors and atomic coordinates. Beyond those, though, I am at a loss as to how to interpret the data. I've looked at the VASP documentation but honestly I'm not exactly clear on what the all of the values mean. Based on what I've read: somewhere in there it should have rho(r)*V-cell, but the position of ?r? or if the V-cell is referring to the total volume of the system or the volume of some tinier, grid-sized, cell remains nebulous. And I assume that I might just not be reading the documentation correctly or overlooking something therein.

If you're curious, I'm looking to get the CHGCAR data into a more manageable format so I can look at various topological properties of electron density and see how the relate to material properties. But first, I'd like to be able to get the electron density data into something I can more easily play around with.

Thanks.

[maartendft]

CHGCAR will give charge density, scaled by the cell volume.

In Matlab, do something like:

[CHG, lat] = importCHGCAR('CHGCAR');
size(CHG); --> a b c

Now you can see how the division of the first cell vector in a steps, the second in b steps and the third in c steps. From this, you can back out the positions.

Alternatively, look at the Electron Localization Function to obtain actual electron "position probabilities" (ELFCAR) in VASP.

<span class='smallblacktext'>[ Edited Fri Sep 07 2012, 06:54PM ]</span>

[maartendft]

Where the importCHGCAR function looks like:

function [CHG,lat] = importCHGCAR(fname)

lat = zeros(3,3);

count = 0;
fid=fopen(fname);
N = 10^16;
while 1
tline = fgetl(fid);
if ~ischar(tline), break, end
count = count + 1;
% Lattice vectors
if count >=3 && count <= 5
lat(count-2,:) = sscanf(tline, '%f %f %f');
end
if count == 7
atomNum = sscanf(tline, '%f %f %f');
N = sum(atomNum);
end
if count == N+10;
CHGdim = sscanf(tline, '%f %f %f');
break
end
end

% Now get all CHG data
data = textscan(fid,'%f');
Nentry = CHGdim(1)*CHGdim(2)*CHGdim(3);
CHG = reshape(data{1}(1:Nentry),CHGdim');

fclose(fid);
end

[JonM]

Thank you very much, this is just what I was looking for.

[askhetan]

dear maartendft,
thank you very much for the script. it works perfectly and i was able to import the data. but i still have some questions to ask. would be very kind if you could please throw some light on these.specifically - what i read the vasp manual is that the CHGCAR is given in the format of
(((C(NX,NY,NZ),NX=1,NGXC),NY=1,NGYZ),NZ=1,NGZC) which makes it 9 columns. but the CHGCAR car file has (after the initial coordinnates and lattice positions) 5 columns and I cannot understand what do these columns stand for. The CHG file has (after the initial coordinnates and lattice positions) 10 columns and I cannot understand that either. When i used your matlab script i got an 84 x 84 x 384 size matrix with several values. i believe each of these values represent the number of electrons in that differential volume dx*dy*dz.... that is basically (a/84)x(b/84)x(c/384) in my case......is that correct ?