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++)
         FReadLine a, str
     endfor
     FReadLine a, str
     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++)
         FReadLine a, str
         FReadLine a, str
         for(j=0;j<str2num(nbnds);j++)
             FReadLine a, str
              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
         FReadLine a, str
         if(strlen(str)==0)
             break
         endif
         if(stringmatch(str,"*k-points in units of 2pi/SCALE*")==1)
             FReadLine a, str
             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++)
                 FReadLine a, str
                 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++)
         FReadLine a, str
     endfor
     FReadLine a, str
     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++)
         FReadLine a, str
         FReadLine a, str
         for(j=0;j<str2num(nbnds);j++)
             FReadLine a, str
              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
         FReadLine a, str
         if(strlen(str)==0)
             break
         endif
         if(stringmatch(str,"*k-points in units of 2pi/SCALE*")==1)
             FReadLine a, str
             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++)
                 FReadLine a, str
                 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
     FReadLine a, str
     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++)
         FReadLine a, str
     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++)
         FReadLine b, str
     endfor
     FReadLine b, str
     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
     FReadLine b, str
     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
    FReadLine a, str
    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++)
         FReadLine a, str
         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++)
             FReadLine a, str
             for(k=0;k<nedos;k++)
                 FReadLine a, str
                 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
     FReadLine a, str
     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++)
         FReadLine a, str
     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++)
         FReadLine b, str
     endfor
     FReadLine b, str
     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
     FReadLine b, str
     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
    FReadLine a, str
    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++)
         FReadLine a, str
         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++)
             FReadLine a, str
             for(k=0;k<(2*nedos);k++)
                 FReadLine a, str
                 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++)
loadwave/a/o/d/g/L={0,21+1030302*i,1030301,0,0} "C:\Users\Lenovo\Desktop\wannier90.bxsf"
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


function extract_fs(path1,fname1,path2,fname2,nstart,nend) //fname1 read;fname2 write
    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++)
         FReadLine a, str
         temp=str[0,strlen(str)-2]
         fprintf b,temp
         fprintf b,"\n"
     endfor
     
     FReadLine a, str
     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
     
     FReadLine a, str
     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++)
         FReadLine a, str
         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++)
             FReadLine a, str
         endfor
     endfor
     
     for(i=0;i<(nend-nstart+1);i++)
         FReadLine a, str
         fprintf b,"BAND:%12d\n",i+1
         for(j=0;j<(nx*ny*nz);j++)
             FReadLine a, str
             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


function vaspbands_proj(path1,fname1,atoms,orbitals) //read procar non-soc
    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
     FReadLine a, str
     FReadLine a, str
     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++)
             FReadLine a, str
         endfor
         for(k=0;k<str2num(nbnds);k++)
             for(i=0;i<3;i++)
                 FReadLine a, str
             endfor
             s=0
             
             for(np=0;np<numpnts(atoms);np++)
               if(np==0)
                 if(atoms[0]==1)
                     FReadLine a, str
                     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++)
                         FReadLine a, str
                     endfor  
                     FReadLine a, str
                     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)
                     FReadLine a, str
                     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++)
                         FReadLine a, str
                     endfor  
                     FReadLine a, str
                     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])
                 FReadLine a, str
                 FReadLine a, str
             else
                 for(i=0;i<(str2num(nions)-atoms[numpnts(atoms)-1]);i++)
                     FReadLine a, str
                 endfor
                 FReadLine a, str
                 FReadLine a, str
             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 vaspbands_projsoc(path1,fname1,atoms,orbitals) //read procar soc
    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
     FReadLine a, str
     FReadLine a, str
     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++)
             FReadLine a, str
         endfor
         for(k=0;k<str2num(nbnds);k++)
             for(i=0;i<3;i++)
                 FReadLine a, str
             endfor
             s=0
             
             for(np=0;np<numpnts(atoms);np++)
               if(np==0)
                 if(atoms[0]==1)
                     FReadLine a, str
                     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++)
                         FReadLine a, str
                     endfor  
                     FReadLine a, str
                     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)
                     FReadLine a, str
                     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++)
                         FReadLine a, str
                     endfor  
                     FReadLine a, str
                     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++)
               FReadLine a, str
            endfor
           
            if(str2num(nions)==atoms[numpnts(atoms)-1])
                 FReadLine a, str
                 FReadLine a, str
                 FReadLine a, str
             else
                 for(i=0;i<(str2num(nions)-atoms[numpnts(atoms)-1]);i++)
                     FReadLine a, str
                 endfor
                 FReadLine a, str
                 FReadLine a, str
                 FReadLine a, str
             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.

 

 

In reply to by jjweimer

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.