How does one parse the chg file to plot the planar average charge density for a surface slab using VASP?? I am trying to calculate the planar averaged charge density difference.
Hi Eric,
You need to subtract the two CHGCAR files to get the charge density difference file and then do a planar averaged plot of that new CHGCAR file.
You can either write your own code for performing these tasks or you can use available tools such as the charge density scripts in the VTST Tools .
I hope that helps.
Hi Eric,
The plane-averaged charge density difference can be written as ∆σ=σ_A+B −σ_A −σ_B . Here, σ_A+B is the plane-averaged density of the combined A+B system, and σ_A and σ_B are the plane-averaged densities of the A and B, respectively, which will calculated by freezing the atomic positions of the respective components in the combined system.
If you are familiar with FORTRAN, the following code calculates the plane-averaged charge density difference (and bunch of other things) for a Cu+P system. Make sure to reformat/rename the CHGCARs accordingly.
***********************************************************************************
PROGRAM readsum
IMPLICIT NONE
INTEGER :: n,ii,nr,nrx,nry,nrz,maxx,maxy,maxz,maxatom
REAL :: totC, x,y,z, lx,ly,lz,vol,scale,k1x,k1y,k1z,&
& k2x,k2y,k2z,k3x,k3y,k3z,ax,ay,az,bx,by,bz,cx,cy,cz,charge,&
&totCW,totCP,totCCu, bcharge,ncharge,Sxy,bchargesum,nchargesum,&
&tchargesum ,d,lbcharge,d1,d2,lchargesum, xycharge,xychargesum,nrz1,nrz2
REAL,allocatable ::CHgridW(:,:,:),CHgridP(:,:,:),CHgridCu(:,:,:) ! dynamic 3D matrix
open(unit=20,file='CHGCARCu',action='read')
open(unit=40,file='CHGCARP',action='read')
open(unit=60,file='CHGCARW',action='read')
open(100,file='calc',FORM='formatted',STATUS='unknown')
open(200,file='denz',FORM='formatted',STATUS='unknown')
! lattice vectors
k1x=
k1y=
k1z=
k2x=
k2y=
k2z=
k3x=
k3y=
k3z=
! number of grid-points
maxx=
maxy=
maxz=
d1= ! n layer below the Cu surface
d2= ! above the top of P
Sxy=(k1x**2+k1y**2)**0.5*(k2y*k2y+k2z*k2z)**0.5
allocate(CHGridW(1:maxx,1:maxy,1:maxz))
read(60,*) (((CHGridW(nrx,nry,nrz),nrx=1,maxx),nry=1,maxy),nrz=1,maxz) !this reads the entire dataset of charges.
allocate(CHGridP(1:maxx,1:maxy,1:maxz))
read(40,*) (((CHGridP(nrx,nry,nrz),nrx=1,maxx),nry=1,maxy),nrz=1,maxz) !this reads the entire dataset of charges.
allocate(CHGridCu(1:maxx,1:maxy,1:maxz))
read(20,*) (((CHGridCu(nrx,nry,nrz),nrx=1,maxx),nry=1,maxy),nrz=1,maxz) !this reads the entire dataset of charges.
bcharge=0.0
ncharge=0.0
totC=0.0
totCW=0.0
totCCu=0.0
totCP=0.0
xycharge=0.0
lbcharge=0.0
charge=0.0
Vol=Sxy*(k3z**2+k3y**2)**0.5
do nrz=1, maxz
do nry=1,maxy
do nrx=1,maxx
ax=nrx*k1x/maxx
ay=nrx*k1y/maxy
az=nrx*k1z/maxz
bx=nry*k2x/maxy
by=nry*k2y/maxy
bz=nry*k2z/maxy
cx=nrz*k1x/maxz
cy=nrz*k3y/maxz
cz=nrz*k3z/maxz
x= ax+bx+cx
y= ay+by+cy
z= az+bz+cz
totCW=totCW+CHgridW(nrx,nry,nrz)
totCCu=totCCu+CHgridCu(nrx,nry,nrz)
totCP=totCP+CHgridP(nrx,nry,nrz)
charge=CHGridW(nrx,nry,nrz)-CHGridP(nrx,nry,nrz)-CHGridCu(nrx,nry,nrz)
xycharge=charge+xycharge
d=(((k3y**2+k3z**2)**0.5)*nrz)/maxz
xycharge=xycharge/(maxx*maxy*maxz)
charge=charge/(maxx*maxy*maxz)
IF (charge > 0 ) THEN
bcharge= charge+bcharge
Else
ncharge= charge+ncharge
endif
enddo
enddo
write(200,*)d,bcharge/Sxy,ncharge/Sxy
bchargesum=bcharge+bchargesum
nchargesum=ncharge+nchargesum
xychargesum=xycharge+xychargesum
lchargesum=lbcharge+lchargesum
bcharge=0.0
ncharge=0.0
lbcharge=0.0
xycharge=0.0
charge=0.0
d=(((k3y**2+k3z**2)**0.5)*nrz)/maxz
nrz1=(d1*maxz)/((k3y**2+k3z**2)**0.5)
nrz2=(d2*maxz)/((k3y**2+k3z**2)**0.5)
enddo
write(100,*) 'The volume of the cell is: ', Vol
write(100,*) 'The charge of P+Cu system =', totCW/(maxx*maxy*maxz)
write(100,*) 'The charge of Cu system=', totCCu/(maxx*maxy*maxz)
write(100,*) 'The charge of P system=', totCP/(maxx*maxy*maxz)
write(100,*) 'The surface charge density gain=', bchargesum/Sxy
write(100,*) 'The charge gain=', bchargesum
write(100,*) 'The charge loss=', nchargesum
write(100,*) 'The total charge=', xychargesum
END PROGRAM readsum
ChGCAR is typical kind of a matrix, i think matlab is most suitable to do this
Hi
You can use vaspkit to generate output file and use vesta to open.