# VaspUtilities

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3     // Use modern global access method and strict wave access.
function vaspbands(path1,fname1,path2,fname2,fermilevel)
string path1,fname1,path2,fname2
variable fermilevel  //fermilevel
variable a,i,t,j,k
string str,fword,nkpts,nbnds
Open/R a path1 + ":" + fname1
for(i=0;i<5;i++)
endfor
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==2)
nkpts=fword
endif
if(t==3)
nbnds=fword
endif
endif
i=i+1
while(t<3)
make/d/n=(str2num(nkpts),str2num(nbnds)) w
for(i=0;i<str2num(nkpts);i++)
for(j=0;j<str2num(nbnds);j++)
k=0
t=0
do
fword=stringfromlist(k,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==2)
w[i][j]=str2num(fword)
endif
endif
k=k+1
while(t<2)
endfor
endfor
for(i=0;i<str2num(nbnds);i++)
make/d/n=(str2num(nkpts)) tmp
tmp=w[p][i]-fermilevel
duplicate tmp \$"wave"+num2str(i)
killwaves tmp
endfor
killwaves w
close a
Open/R a path2 + ":" + fname2
make/d/n=3 kold,knew
make/d/n=(str2num(nkpts)) kpt
do
if(strlen(str)==0)
break
endif
if(stringmatch(str,"*k-points in units of 2pi/SCALE*")==1)
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
kold[0]=str2num(fword)
elseif(t==2)
kold[1]=str2num(fword)
elseif(t==3)
kold[2]=str2num(fword)
endif
endif
i=i+1
while(t<3)
kpt[0]=0
for(i=1;i<str2num(nkpts);i++)
k=0
t=0
do
fword=stringfromlist(k,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
knew[0]=str2num(fword)
elseif(t==2)
knew[1]=str2num(fword)
elseif(t==3)
knew[2]=str2num(fword)
endif
endif
k=k+1
while(t<3)
if(mod(i,100)==0)
kpt[i]=kpt[i-1]
else
kpt[i]=kpt[i-1]+sqrt((knew[0]-kold[0])^2+(knew[1]-kold[1])^2+(knew[2]-kold[2])^2)
endif
kold=knew
endfor
endif
while(1)
killwaves kold,knew
close a
for(i=0;i<str2num(nbnds);i++)
if(i==0)
display \$"wave"+num2str(i) vs kpt
else
appendtograph \$"wave"+num2str(i) vs kpt
endif
endfor
end

function vaspbands_mbj(path1,fname1,path2,fname2,fermilevel)
string path1,fname1,path2,fname2
variable fermilevel  //fermilevel
variable a,i,t,j,k
variable cc=120
string str,fword,nkpts,nbnds
Open/R a path1 + ":" + fname1
for(i=0;i<5;i++)
endfor
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==2)
nkpts=fword
endif
if(t==3)
nbnds=fword
endif
endif
i=i+1
while(t<3)
make/d/n=(str2num(nkpts),str2num(nbnds)) w
for(i=0;i<str2num(nkpts);i++)
for(j=0;j<str2num(nbnds);j++)
k=0
t=0
do
fword=stringfromlist(k,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==2)
w[i][j]=str2num(fword)
endif
endif
k=k+1
while(t<2)
endfor
endfor
for(i=0;i<str2num(nbnds);i++)
make/d/n=(str2num(nkpts)) tmp
tmp=w[p][i]-fermilevel
duplicate tmp \$"wave"+num2str(i)
killwaves tmp
endfor
killwaves w
close a
Open/R a path2 + ":" + fname2
make/d/n=3 kold,knew
make/d/n=(str2num(nkpts)) kpt
do
if(strlen(str)==0)
break
endif
if(stringmatch(str,"*k-points in units of 2pi/SCALE*")==1)
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
kold[0]=str2num(fword)
elseif(t==2)
kold[1]=str2num(fword)
elseif(t==3)
kold[2]=str2num(fword)
endif
endif
i=i+1
while(t<3)
kpt[0]=0
for(i=1;i<str2num(nkpts);i++)
k=0
t=0
do
fword=stringfromlist(k,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
knew[0]=str2num(fword)
elseif(t==2)
knew[1]=str2num(fword)
elseif(t==3)
knew[2]=str2num(fword)
endif
endif
k=k+1
while(t<3)
if(mod(i,100)==0)
kpt[i]=kpt[i-1]
else
kpt[i]=kpt[i-1]+sqrt((knew[0]-kold[0])^2+(knew[1]-kold[1])^2+(knew[2]-kold[2])^2)
endif
kold=knew
endfor
endif
while(1)
killwaves kold,knew
close a
for(i=0;i<str2num(nbnds);i++)
DeletePoints 0,cc, \$"wave"+num2str(i)
endfor
DeletePoints 0,cc,kpt
for(i=0;i<str2num(nbnds);i++)
if(i==0)
display \$"wave"+num2str(i) vs kpt
else
appendtograph \$"wave"+num2str(i) vs kpt
endif
endfor
end

function splitdos(path1,fname1,path2,fname2) //fname1=DOSCAR,fname2=POSCAR
string path1,fname1,path2,fname2
variable a,b,i,t,j,k
string str,fword,natms,nspin
Open/R a path1 + ":" + fname1
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
natms=fword
endif
if(t==3)
nspin=fword
endif
endif
i=i+1
while(t<3)
for(i=0;i<4;i++)
endfor
printf "Number of atoms: %d\r  Number of spin: %d\r" str2num(natms),str2num(nspin)
Open/R b path2 + ":" + fname2
for(i=0;i<5;i++)
endfor
i=0
variable/g nspec=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
nspec=nspec+1
endif
i=i+1
while(i<strlen(str))
make/t/n=(nspec)/o atmnam
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==nspec)
variable memory=strlen(fword)
string lastnm
sscanf fword[0,memory-1],"%s",lastnm
atmnam[t-1]=lastnm
else
atmnam[t-1]=cleanupname(fword,0)
endif
endif
i=i+1
while(t<nspec)
make/d/n=(nspec)/o atmnum
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
atmnum[t-1]=str2num(fword)
endif
i=i+1
while(t<nspec)
close b
printf "Number of species: %d\r" nspec
variable temp=0
for(i=0;i<nspec;i++)
temp+=atmnum[i]
endfor
if(str2num(natms)!=temp)
printf "DOSCAR info & POSCAR does not match!\r"
close a
return -1
endif
i=0
t=0
do
fword=stringfromlist(i,str," ",32)
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
variable/g nedos=str2num(fword)
endif
if(t==2)
variable/g efermi=str2num(fword)
endif
endif
i=i+1
while(t<2)
for(i=0;i<nedos;i++)
if(i==0)
j=0
variable tl=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
tl=tl+1
endif
j=j+1
while(j<strlen(str))
make/o/n=(nedos,tl)/d pdos_tot
endif
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
pdos_tot[i][t-1]=str2num(fword)
endif
j=j+1
while(t<tl)
endfor
for(i=0;i<nspec;i++)
for(j=0;j<atmnum[i];j++)
for(k=0;k<nedos;k++)
if(k==0)
variable in=0
tl=0
do
fword=stringfromlist(in,str," ")
if (stringmatch(fword,"")!=1)
tl=tl+1
endif
in=in+1
while(in<strlen(str))
make/o/n=(nedos,tl)/d tempwave
endif
in=0
t=0
do
fword=stringfromlist(in,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
tempwave[k][t-1]=str2num(fword)
endif
in=in+1
while(t<tl)
endfor
duplicate/o tempwave, \$"pdos_"+atmnam[i]+"_"+num2str(j+1)
endfor
endfor
killwaves tempwave
end

function splitdosif(path1,fname1,path2,fname2)
string path1,fname1,path2,fname2
variable a,b,i,t,j,k
string str,fword,natms,nspin
Open/R a path1 + ":" + fname1
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
natms=fword
endif
if(t==3)
nspin=fword
endif
endif
i=i+1
while(t<3)
for(i=0;i<4;i++)
endfor
printf "Number of atoms: %d\r  Number of spin: %d\r" str2num(natms),str2num(nspin)
Open/R b path2 + ":" + fname2
for(i=0;i<5;i++)
endfor
i=0
variable/g nspec=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
nspec=nspec+1
endif
i=i+1
while(i<strlen(str))
make/t/n=(nspec)/o atmnam
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==nspec)
variable memory=strlen(fword)
string lastnm
sscanf fword[0,memory-1],"%s",lastnm
atmnam[t-1]=lastnm
else
atmnam[t-1]=cleanupname(fword,0)
endif
endif
i=i+1
while(t<nspec)
make/d/n=(nspec)/o atmnum
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
atmnum[t-1]=str2num(fword)
endif
i=i+1
while(t<nspec)
close b
printf "Number of species: %d\r" nspec
variable temp=0
for(i=0;i<nspec;i++)
temp+=atmnum[i]
endfor
if(str2num(natms)!=temp)
printf "DOSCAR info & POSCAR does not match!\r"
close a
return -1
endif
i=0
t=0
do
fword=stringfromlist(i,str," ",32)
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
variable/g nedos=str2num(fword)
endif
if(t==2)
variable/g efermi=str2num(fword)
endif
endif
i=i+1
while(t<2)
for(i=0;i<nedos;i++)
if(i==0)
j=0
variable tl=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
tl=tl+1
endif
j=j+1
while(j<strlen(str))
make/o/n=(nedos,tl)/d pdos_tot
endif
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
pdos_tot[i][t-1]=str2num(fword)
endif
j=j+1
while(t<tl)
endfor
for(i=0;i<nspec;i++)
for(j=0;j<atmnum[i];j++)
for(k=0;k<(2*nedos);k++)
if(k==0)
variable in=0
tl=0
do
fword=stringfromlist(in,str," ")
if (stringmatch(fword,"")!=1)
tl=tl+1
endif
in=in+1
while(in<strlen(str))
make/o/n=(nedos,tl)/d tempwave
endif
if(k==1)
in=0
variable tll=0
do
fword=stringfromlist(in,str," ")
if (stringmatch(fword,"")!=1)
tll=tll+1
endif
in=in+1
while(in<strlen(str))
redimension/n=(nedos,(tl+tll))/d tempwave
endif
in=0
t=0
if(mod(k,2)==0)
do
fword=stringfromlist(in,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
tempwave[(k/2)][t-1]=str2num(fword)
endif
in=in+1
while(t<tl)
else
do
fword=stringfromlist(in,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
tempwave[floor(k/2)][tl+t-1]=str2num(fword)
endif
in=in+1
while(t<tll)
endif
endfor
duplicate/o tempwave, \$"pdos_"+atmnam[i]+"_"+num2str(j+1)
endfor
endfor
killwaves tempwave
end

function sumdos()
variable i,j
variable/g nspec,nedos
wave atmnum
wave/t atmnam
for(i=0;i<nspec;i++)
for(j=0;j<atmnum[i];j++)
if(j==0)
duplicate/o \$"pdos_"+atmnam[i]+"_"+num2str(j+1),tempwave
duplicate/o/r=[][0] \$"pdos_"+atmnam[i]+"_"+num2str(j+1),energywave
else
duplicate/o \$"pdos_"+atmnam[i]+"_"+num2str(j+1),tempwave2
tempwave+=tempwave2
killwaves tempwave2
endif
endfor
variable k
for(k=0;k<nedos;k++)
tempwave[k][0]=energywave[k]
endfor
duplicate/o tempwave \$"pdos_"+atmnam[i]+"_tot"
endfor
killwaves tempwave,energywave
end

function plotdosif()
variable i
variable/g nspec
wave/t atmnam
for(i=0;i<nspec;i++)
variable j
for(j=0;j<65;j++)
if(j==1)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" \$"pdos_"+atmnam[i]+"_tot_s"
endif
if(j==5)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_p
endif
if(j==9)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_p_t
temp_p+=temp_p_t
endif
if(j==13)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_p_t
temp_p+=temp_p_t
duplicate/o temp_p,\$"pdos_"+atmnam[i]+"_tot_p"
killwaves temp_p_t,temp_p
endif
if(j==17)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d
endif
if(j==21)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
endif
if(j==25)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
endif
if(j==29)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
endif
if(j==33)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
duplicate/o temp_d,\$"pdos_"+atmnam[i]+"_tot_d"
killwaves temp_d_t,temp_d
endif
if(j==37)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_f
endif
if(j==41)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_f_t
temp_f+=temp_f_t
endif
if(j==45)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_f_t
temp_f+=temp_f_t
endif
if(j==49)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_f_t
temp_f+=temp_f_t
endif
if(j==53)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_f_t
temp_f+=temp_f_t
endif
if(j==57)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_f_t
temp_f+=temp_f_t
endif
if(j==61)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_f_t
temp_f+=temp_f_t
duplicate/o temp_f,\$"pdos_"+atmnam[i]+"_tot_f"
killwaves temp_f_t,temp_f
endif
endfor
endfor
for(i=0;i<nspec;i++)
if(i==0)
display \$"pdos_"+atmnam[i]+"_tot_s" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
else
appendtograph \$"pdos_"+atmnam[i]+"_tot_s" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
endif
appendtograph \$"pdos_"+atmnam[i]+"_tot_p" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
appendtograph \$"pdos_"+atmnam[i]+"_tot_d" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
appendtograph \$"pdos_"+atmnam[i]+"_tot_f" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
endfor
end

function plotdoslf()
variable i
variable/g nspec
wave/t atmnam
for(i=0;i<nspec;i++)
variable j
for(j=0;j<65;j++)
if(j==1)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" \$"pdos_"+atmnam[i]+"_tot_s"
endif
if(j==5)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_p
endif
if(j==9)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_p_t
temp_p+=temp_p_t
endif
if(j==13)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_p_t
temp_p+=temp_p_t
duplicate/o temp_p,\$"pdos_"+atmnam[i]+"_tot_p"
killwaves temp_p_t,temp_p
endif
if(j==17)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d
endif
if(j==21)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
endif
if(j==25)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
endif
if(j==29)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
endif
if(j==33)
duplicate/o/r=[][j] \$"pdos_"+atmnam[i]+"_tot" temp_d_t
temp_d+=temp_d_t
duplicate/o temp_d,\$"pdos_"+atmnam[i]+"_tot_d"
killwaves temp_d_t,temp_d
endif
endfor
endfor
for(i=0;i<nspec;i++)
if(i==0)
display \$"pdos_"+atmnam[i]+"_tot_s" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
else
appendtograph \$"pdos_"+atmnam[i]+"_tot_s" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
endif
appendtograph \$"pdos_"+atmnam[i]+"_tot_p" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
appendtograph \$"pdos_"+atmnam[i]+"_tot_d" vs \$"pdos_"+atmnam[i]+"_tot"[][0]
endfor
end

function Findbands(numbnds,efermi)
variable numbnds,efermi
variable i
make/o/n=(numbnds) mini,maxm
for(i=0;i<numbnds;i++)
endfor
for(i=0;i<numbnds;i++)
wavestats \$"wave"+num2str(i)
mini[i]=V_min
maxm[i]=V_max
killwaves \$"wave"+num2str(i)
endfor
for(i=0;i<numbnds;i++)
if((efermi>mini[i])&&(efermi<maxm[i]))
print i+1
endif
endfor
end

string path1,fname1,path2,fname2
variable nstart,nend
variable a,b,i,t,j,k
string str,temp,fword
variable nx,ny,nz
Open/R a path1 + ":" + fname1
Open b path2 + ":" + fname2
for(i=0;i<14;i++)
temp=str[0,strlen(str)-2]
fprintf b,temp
fprintf b,"\n"
endfor

i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
temp=fword[0,strlen(fword)-2]
endif
endif
i=i+1
while(t<1)

if((nend<nstart)||(str2num(temp)<nend))
fprintf b,"Wrong indices,exit...\n"
close a
close b
return -1
endif

fprintf b,"%15d\n",nend-nstart+1

temp=str[0,strlen(str)-2]
fprintf b,temp
fprintf b,"\n"
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==1)
temp=fword[0,strlen(fword)-1]
nx=str2num(temp)
endif
if(t==2)
temp=fword[0,strlen(fword)-1]
ny=str2num(temp)
endif
if(t==3)
temp=fword[0,strlen(fword)-2]
nz=str2num(temp)
endif
endif
i=i+1
while(t<3)

for(i=0;i<4;i++)
temp=str[0,strlen(str)-2]
fprintf b,temp
fprintf b,"\n"
endfor

for(i=1;i<nstart;i++)
for(j=0;j<(nx*ny*nz+1);j++)
endfor
endfor

for(i=0;i<(nend-nstart+1);i++)
fprintf b,"BAND:%12d\n",i+1
for(j=0;j<(nx*ny*nz);j++)
temp=str[0,strlen(str)-2]
fprintf b,temp
fprintf b,"\n"
endfor
endfor

fprintf b," END_BANDGRID_3D\n"
fprintf b,"  END_BLOCK_BANDGRID_3D\n"

close a
close b
end

wave atoms
wave orbitals
string path1,fname1
variable a,i,t,j,k,u
variable s
variable inner,np,ccc
string str,fword,nkpts,nbnds,nions
Open/R a path1 + ":" + fname1
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==4)
nkpts=fword
endif
if(t==8)
nbnds=fword
endif
if(t==12)
nions=fword
endif
endif
i=i+1
while(t<12)
make/o/d/n=(str2num(nkpts),str2num(nbnds)) tem
for(u=0;u<str2num(nkpts);u++)
for(i=0;i<3;i++)
endfor
for(k=0;k<str2num(nbnds);k++)
for(i=0;i<3;i++)
endfor
s=0

for(np=0;np<numpnts(atoms);np++)
if(np==0)
if(atoms[0]==1)
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
if(atoms[0]>1)
for(inner=0;inner<(atoms[0]-1);inner++)
endfor
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
else
if((atoms[np]-atoms[np-1])==1)
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
if((atoms[np]-atoms[np-1])>1)
for(inner=0;inner<(atoms[np]-atoms[np-1]-1);inner++)
endfor
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
endif
endfor

if(str2num(nions)==atoms[numpnts(atoms)-1])
else
for(i=0;i<(str2num(nions)-atoms[numpnts(atoms)-1]);i++)
endfor
endif

tem[u][k]=s

endfor
endfor
for(i=0;i<str2num(nbnds);i++)
duplicate/o/r=[0,str2num(nkpts)-1][i,i] tem \$"weight"+num2str(i)
endfor
close a
end

wave atoms
wave orbitals
string path1,fname1
variable a,i,t,j,k,u
variable s
variable inner,np,ccc
string str,fword,nkpts,nbnds,nions
Open/R a path1 + ":" + fname1
i=0
t=0
do
fword=stringfromlist(i,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
if(t==4)
nkpts=fword
endif
if(t==8)
nbnds=fword
endif
if(t==12)
nions=fword
endif
endif
i=i+1
while(t<12)
make/o/d/n=(str2num(nkpts),str2num(nbnds)) tem
for(u=0;u<str2num(nkpts);u++)
for(i=0;i<3;i++)
endfor
for(k=0;k<str2num(nbnds);k++)
for(i=0;i<3;i++)
endfor
s=0

for(np=0;np<numpnts(atoms);np++)
if(np==0)
if(atoms[0]==1)
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
if(atoms[0]>1)
for(inner=0;inner<(atoms[0]-1);inner++)
endfor
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
else
if((atoms[np]-atoms[np-1])==1)
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
if((atoms[np]-atoms[np-1])>1)
for(inner=0;inner<(atoms[np]-atoms[np-1]-1);inner++)
endfor
j=0
t=0
do
fword=stringfromlist(j,str," ")
if (stringmatch(fword,"")!=1)
t=t+1
for(ccc=0;ccc<numpnts(orbitals);ccc++)
if(t==orbitals[ccc])
s+=str2num(fword)
endif
endfor
endif
j=j+1
while(t<orbitals[numpnts(orbitals)-1])
endif
endif
endfor

variable jjj
for(jjj=0;jjj<(3*(str2num(nions)+1+1));jjj++)
endfor

if(str2num(nions)==atoms[numpnts(atoms)-1])
else
for(i=0;i<(str2num(nions)-atoms[numpnts(atoms)-1]);i++)
endfor
endif

tem[u][k]=s

endfor
endfor
for(i=0;i<str2num(nbnds);i++)
duplicate/o/r=[0,str2num(nkpts)-1][i,i] tem \$"weight"+num2str(i)
endfor
close a
end

function projbulkbnds(m,n,p) //projected bulk bands plot.m:the number of points between two kpoints;
//n:number of different kzs;p:numberof bands.
variable m,n,p
variable i,j,k
wave kpt
make/n=(2*m) refer
for(i=0;i<(2*m);i++)
refer[i]=kpt[i]
endfor
make/n=(2*m) tem1
for(i=0;i<p;i++)
duplicate \$"wave"+num2str(i),tem2
for(j=i;j<i+p*(n-1)+1;j=j+p)
for(k=0;k<(2*m);k++)
tem1[k]=tem2[k]
endfor
duplicate tem1,\$"w"+num2str(j)
if(j==0)
display \$"w"+num2str(j) vs refer
else
appendtograph \$"w"+num2str(j) vs refer
endif
DeletePoints 0,3*m, tem2
endfor
killwaves tem2
endfor
killwaves tem1
end

function difkzplot(nbnds,nkz,ymin,ymax,j)
variable nbnds
variable nkz
variable ymin,ymax //display range
variable j  //j can take value 0-20,corresponding to kz=j/(nkz-1)*(gamma-z)
variable kz=j/(nkz-1)
wave xwave
variable i
for(i=0;i<nbnds;i++)
if(i==0)
display \$"w"+num2str(j*nbnds) vs xwave
TextBox/C/N=text0/F=0/A=MC num2str(kz)+"*G-Z"
SetAxis left ymin,ymax
else
appendtograph \$"w"+num2str(j*nbnds+i) vs xwave
endif
endfor
end

function fcc111()  //for writing fcc 111 KPOINTS
variable a,b,c
a=    5.429
make/o/d/n=3 a1={0,a/2,a/2}
make/o/d/n=3 a2={a/2,0,a/2}
make/o/d/n=3 a3={a/2,a/2,0}
cross a2,a3
variable vprim
wave W_Cross
vprim=matrixDot(a1,W_Cross)
make/o/d/n=3 b1,b2,b3
b1=2*pi*W_Cross/vprim
cross a3,a1
b2=2*pi*W_Cross/vprim
cross a1,a2
b3=2*pi*W_Cross/vprim
Concatenate/O {b1,b2,b3},primbasis

make/o/d/n=3 aa1={-a/2,a/2,0}
make/o/d/n=3 aa2={0,-a/2,a/2}
make/o/d/n=3 aa3={a,a,a}
cross aa2,aa3
variable vconv
vconv=matrixDot(aa1,W_Cross)
make/o/d/n=3 bb1,bb2,bb3
bb1=2*pi*W_Cross/vconv
cross aa3,aa1
bb2=2*pi*W_Cross/vconv

bb3=1/2*(b1+b2+b3)
Concatenate/O {bb1,bb2,bb3},convbasis

matrixinverse primbasis
wave M_inverse
matrixop/o matC=M_inverse x convbasis

print vprim
print vprim/vconv
cross b2,b3
variable kvprim
kvprim=matrixDot(b1,W_Cross)
cross bb2,bb3
variable kvconv
kvconv=matrixDot(bb1,W_Cross)
print kvprim/kvconv
end

function simplecubic111()  //for writing simple cubic 111 KPOINTS
variable a,b,c
a=    2.3537117117709037*2
make/o/d/n=3 a1={a,0,0}
make/o/d/n=3 a2={0,a,0}
make/o/d/n=3 a3={0,0,a}
cross a2,a3
variable vprim
wave W_Cross
vprim=matrixDot(a1,W_Cross)
make/o/d/n=3 b1,b2,b3
b1=2*pi*W_Cross/vprim
cross a3,a1
b2=2*pi*W_Cross/vprim
cross a1,a2
b3=2*pi*W_Cross/vprim
Concatenate/O {b1,b2,b3},primbasis

make/o/d/n=3 aa1={a,-a,0}
make/o/d/n=3 aa2={0,a,-a}
make/o/d/n=3 aa3={a,a,a}
cross aa2,aa3
variable vconv
vconv=matrixDot(aa1,W_Cross)
make/o/d/n=3 bb1,bb2,bb3
bb1=2*pi*W_Cross/vconv
cross aa3,aa1
bb2=2*pi*W_Cross/vconv

bb3=1/2*(b1+b2+b3)
Concatenate/O {bb1,bb2,bb3},convbasis

matrixinverse primbasis
wave M_inverse
matrixop/o matC=M_inverse x convbasis

print vprim/vconv
cross b2,b3
variable kvprim
kvprim=matrixDot(b1,W_Cross)
cross bb2,bb3
variable kvconv
kvconv=matrixDot(bb1,W_Cross)
print kvprim/kvconv
end

function hr111()  //for writing hr 111 KPOINTS
variable a,c
a=    3.6281/sqrt(3)
c=26.233/3
make/o/d/n=3 a1={-a/2,sqrt(3)/2*a,c}
make/o/d/n=3 a2={-a/2,-sqrt(3)/2*a,c}
make/o/d/n=3 a3={a,0,c}
cross a2,a3
variable vprim
wave W_Cross
vprim=matrixDot(a1,W_Cross)
make/o/d/n=3 b1,b2,b3
b1=2*pi*W_Cross/vprim
cross a3,a1
b2=2*pi*W_Cross/vprim
cross a1,a2
b3=2*pi*W_Cross/vprim
Concatenate/O {b1,b2,b3},primbasis

make/o/d/n=3 aa1=a2-a3
make/o/d/n=3 aa2=a3-a1
make/o/d/n=3 aa3=a1+a2+a3
cross aa2,aa3
variable vconv
wave W_cross
vconv=matrixDot(aa1,W_Cross)
make/o/d/n=3 bb1,bb2,bb3
bb1=2*pi*W_Cross/vconv
cross aa3,aa1
bb2=2*pi*W_Cross/vconv

bb3=1/2*(b1+b2+b3)
Concatenate/O {bb1,bb2,bb3},convbasis

matrixinverse primbasis
wave M_inverse
matrixop/o matC=M_inverse x convbasis

print vprim
print vprim/vconv
cross b2,b3
variable kvprim
kvprim=matrixDot(b1,W_Cross)
cross bb2,bb3
variable kvconv
kvconv=matrixDot(bb1,W_Cross)
print kvprim/kvconv
end

function changecoord() //for KPOINTS coordinate transformation
variable x,y,z
variable m,n,l
wave matC
variable i
variable k=0
for(i=0;i<21;i++)
x=1/2;y=0;z=i/20 // M bar point
m=matC[0][0]*x+matC[0][1]*y+matC[0][2]*z
n=matC[1][0]*x+matC[1][1]*y+matC[1][2]*z
l=matC[2][0]*x+matC[2][1]*y+matC[2][2]*z
if(i==0)
printf "%.12f %.12f %.12f\r" m, n, l
else
printf "%.12f %.12f %.12f\r" m, n, l
printf "%.12f %.12f %.12f\r" m, n, l
endif
x=0;y=0;z=i/20 //gamma bar point
m=matC[0][0]*x+matC[0][1]*y+matC[0][2]*z
n=matC[1][0]*x+matC[1][1]*y+matC[1][2]*z
l=matC[2][0]*x+matC[2][1]*y+matC[2][2]*z
printf "%.12f %.12f %.12f\r" m, n, l
printf "%.12f %.12f %.12f\r" m, n, l
x=1/3;y=1/3;z=i/20 //K bar point
m=matC[0][0]*x+matC[0][1]*y+matC[0][2]*z
n=matC[1][0]*x+matC[1][1]*y+matC[1][2]*z
l=matC[2][0]*x+matC[2][1]*y+matC[2][2]*z
if(i==20)
printf "%.12f %.12f %.12f\r" m, n, l
else
printf "%.12f %.12f %.12f\r" m, n, l
printf "%.12f %.12f %.12f\r" m, n, l
endif
endfor
end

function fcc111slab(nlayer,a) //for writing fcc 111 slab POSCAR
variable nlayer,a
make/o/d matA={{-1/2,1/2,0},{0,-1/2,1/2},{1,1,1}}
matrixinverse matA
wave M_inverse
make/o/d matB={{0,1/2,1/2},{1/2,0,1/2},{1/2,1/2,0}}
matrixop/o matC=M_inverse x matB
make/o/d/n=(3,3) Cefilm
variable x,y,z,m,n,l,i
for(i=0;i<3;i++)
x=0;y=0;z=i
m=matC[0][0]*x+matC[0][1]*y+matC[0][2]*z
n=matC[1][0]*x+matC[1][1]*y+matC[1][2]*z
l=matC[2][0]*x+matC[2][1]*y+matC[2][2]*z
do
if(m<0)
m=m+1
elseif((m>1)||(m==1))
m=m-1
endif
while((m<0)||(m>1)||(m==1))
do
if(n<0)
n=n+1
elseif((n>1)||(n==1))
n=n-1
endif
while((n<0)||(n>1)||(n==1))
do
if(l<0)
l=l+1
elseif((l>1)||(l==1))
l=l-1
endif
while((l<0)||(l>1)||(l==1))
Cefilm[i][0]=m ;Cefilm[i][1]=n ;Cefilm[i][2]=l
endfor
bubble(Cefilm,3)
redimension/n=(nlayer,3) Cefilm
for(i=3;i<nlayer;i++)
Cefilm[i][0]=Cefilm[i-3][0]
Cefilm[i][1]=Cefilm[i-3][1]
Cefilm[i][2]=Cefilm[i-3][2]+1
endfor
for(i=0;i<nlayer;i++)
printf "%.12f %.12f %.12f\r",Cefilm[i][0],Cefilm[i][1], (Cefilm[i][2]*a*sqrt(3)/((Cefilm[nlayer-1][2]-Cefilm[0][2])*a*sqrt(3)+18))
endfor
print "\r"
printf "%.12f\r" (Cefilm[nlayer-1][2]-Cefilm[0][2])*a*sqrt(3)+18
end

function bubble(w,n)
wave w
variable n
variable i,j,k,t1,t2,t3
for(i=1;i<n;i++)
k=0
for(j=0;j<(n-i);j++)
if(w[j][2]>w[j+1][2])
t1=w[j][0];w[j][0]=w[j+1][0];w[j+1][0]=t1
t2=w[j][1];w[j][1]=w[j+1][1];w[j+1][1]=t2
t3=w[j][2];w[j][2]=w[j+1][2];w[j+1][2]=t3
else
k++
endif
endfor
if(k==(n-i))
break
endif
endfor
end

function xcoordi()
wave b1,b2
make/o/d/n=200 xwave
variable i
for(i=0;i<100;i++)
xwave[i]=1/2*sqrt(b1[0]^2+b1[1]^2+b1[2]^2)/99*(i-99)
endfor
for(i=100;i<200;i++)
xwave[i]=1/3*sqrt((b1[0]+b2[0])^2+(b1[1]+b2[1])^2+(b1[2]+b2[2])^2)/99*(i-100)
endfor
end

Vasp utilities written in Igor. It can be used to plot bands, projected dos, weighted bands and find bands cross the fermi level. It can also be used to write POSCAR and KPOINTS for several crystal structures.

Would you be so kind to do two things: Package this as an Igor Pro pxp file and provide an example of the plots that it can produce.

I don't quite understand what is meant by "Package this as an Igor Pro pxp file". Could you explain it to me? If you want to do some plot, you can just run one of the functions, e.g. vaspbands("D","EIGENVAL","D","OUTCAR",1.2). As for the plots it can give, I can show some graphs I have done later.

In reply to by Yuan Fang

Yuan Fang wrote:

If you want to do some plot, you can just run one of the functions, e.g. vaspbands("D","EIGENVAL","D","OUTCAR",1.2).

Well, it seems that some files are required which many of the forum users will not have readily available.

> I don't quite understand what is meant by "Package this as an Igor Pro pxp file".

* Open Igor Pro.
* Copy the code into the procedure window.
* Save the Igor Pro experiment (e.g. as VASP Utitlities Code.pxp).
* Post the pxp file.

The better approach for long code that does something specific such as your example is to create an Igor Pro project rather than post the raw code as forum discussion. See the brief information and longer list here. https://www.wavemetrics.com/projects

Finally, it would be helpful when the pxp project file contained sample raw data so that the functions could be run to create example plots directly. Then, it could be called VASP Utilities Demo.pxp by example.

Here is a synopsis of the steps to make what you have posted into a project on Igor Exchange.

Create an Igor Pro experiment (pxp file) with the code. Include example data in the experiment. Make plots using the example data. Save that file as VASP Demo. Compress that file as a ZIP archive. Create a NEW PROJECT on this site. Upload the ZIP archive.