Create a new version of paste:

! this program written in fortran 90 analyzes the DOSCAR output file of VASP


! depending on user answers, it may compute the total density of states, or projected density of states of any atom labeled as in the POSCAR file.


! the program askes the user what it needs during execution. no input file is needed. for simplicity.


! the energies may be translated so that Fermi energy is the reference (abscissa is E-Efermi) if user wants.


! it has been written by Maximilien Levesque, while in postdoc @ Ecole Normale Superieure, Paris


! in the group of theoretical physical-chemistry of Daniel Borgis


! it is free of any copyright. just send me an email max.....en.lev....e at gmail.com for thanks :) or bug reports.


! it's based on VASP 4.6 documentation written found at http://cms.mpi.univie.ac.at/vasp/vasp/DOSCAR_file.html


! I recommand to use the GNU fortran compiler (gfortran)


! for compilation, just type in almost any terminal :    gfortran doscar_analysis_spd.f90 -o doscar_analysis_spd


! and then execute it :    ./doscar_analysis_spd


! the POSCAR file should be in the directory from which you call doscar_analysis_spd


! for instance, if you have compilde doscar_analysis in  /home/max/bin , and your DOSCAR is in /home/max/vaspoutput/DOSCAR


! type :            cd /home/max/vaspoutput/


! and then :        ../bin/doscar_analysis_spd


! Maximilien Levesque 201108040017

! Maximilien Levesque 201109011504 debug of number of steps in case of local density of states


program doscar_analysis


implicit none


integer :: number_of_atoms ! total number of atoms in POSCAR


integer :: i ! dummy


double precision :: t ! dummy


double precision :: Emax , Emin ! max and min value of the energy (absciss limits)


integer :: number_of_steps ! discretization of the abscissa


double precision :: Efermi ! Fermi energy of the system


integer :: choice ! choice made by user


integer :: choice_translate ! choice made by user : translate or not to fermi energy


integer :: choice_polarization ! choice made by user : polarization enabled or not


integer :: choice_atom ! choice made by user : label of the atom the user wants the PDOS of


integer :: choice_spdf ! choice made by user : whether to specify orbital contributions separately


double precision :: dos_up , dos_down ! local value of dos spin up and spin down


double precision :: dos_up_s, dos_down_s ! local values of dos spin up and spin down for s-shell


double precision :: dos_up_p, dos_down_p ! local values of dos spin up and spin down for p-shell


double precision :: dos_up_d, dos_down_d ! local values of dos spin up and spin down for d-shell


double precision :: n1=0, n2=0, n3=0, n4=0, n5=0, n6=0, n7=0, n8=0, n9=0


double precision :: n10=0, n11=0, n12=0, n13=0, n14=0, n15=0, n16=0, n17=0, n18=0


! open DOSCAR


open ( 10 , file = 'DOSCAR' , form = 'formatted' )


! first line is the total number of atoms and 3 numbers i don't know what they mean


read ( 10 , * ) number_of_atoms , i , i , i


! line 2 to 5 are useless


do i = 2 , 5 ; read ( 10 , * ) ; end do


! line 6 contains Emax , Emin , number_of_steps , Efermi and ?


read ( 10 , * ) Emax , Emin , number_of_steps , Efermi , t


! ask user what he wants of me


77 write ( * , * ) "Do you want to compute the total density of states (press 1) or a projected one (press 2) ?"


read ( * , * ) choice


! test if answer is correct


if ( choice /= 1 .and. choice /= 2 ) then


  write ( * , * ) "Only 1 or 2 is possible."


  goto 77 


end if


! ask user if he wants to translate to Fermi level or not


88 write ( * , * ) "Do you want to translate all the energies to Fermi energy (press 1 for yes , press 2 for no) ?"


read ( * , * ) choice_translate


! test if answer is correct


if ( choice_translate /= 1 .and. choice_translate /= 2 ) then


  write ( * , * ) "Only 1 or 2 is possible."


  goto 88 


end if


! ask user if the calculation is polarized or not


99 write ( * , * ) "Is system polarized (press 1 for yes , press 2 for no) ?"


read ( * , * ) choice_polarization


! test if answer is correct


if ( choice_polarization /= 1 .and. choice_polarization /= 2 ) then


  write ( * , * ) "Only 1 or 2 is possible."


  goto 99 


end if


! open file for writting output


open ( unit = 11 , file = 'doscar_analysis_spd.dat' , form = 'formatted' )


! write it Fermi energy


write ( 11 , * ) "# E fermi = " , Efermi


! if total DOS is wanted


if ( choice == 1 ) then


  ! read the number_of_steps next lines.


  if ( choice_polarization == 1 ) then


    do i = 1 , number_of_steps


      read ( 10 , * ) t , dos_up , dos_down


      write ( 11 , * ) "Energy - DOS_UP - DOS_DOWN - DOS_TOT"


      if ( choice_translate == 1 ) write ( 11 , * ) t - Efermi , dos_up , dos_down , dos_up + dos_down


      if ( choice_translate == 2 ) write ( 11 , * ) t , dos_up , dos_down , dos_up + dos_down


    end do 


  else if ( choice_polarization == 2 ) then


    do i = 1 , number_of_steps


      read ( 10 , * ) t , dos_up

    
      write ( 11 , * ) "Energy - DOS_TOT"


      if ( choice_translate == 1 ) write ( 11 , * ) t - Efermi , dos_up


      if ( choice_translate == 2 ) write ( 11 , * ) t , dos_up


    end do 


  end if


else if ( choice == 2 ) then


  ! if projected DOS is wanted


  ! first ask user which atom he is interested in


  1010 write ( * , * ) "Which atom (same label as in POSCAR and first label is 1 (not 0)) ?"


  read ( * , * ) choice_atom


  ! test if answer is correct


  if ( choice_atom < 1 .or. choice_atom > number_of_atoms ) then


    write ( * , * ) "Wrong answer"


    goto 1010 


  end if


  ! Next ask if the specific orbital DOS are to be printed as well


  1111 write ( * , * ) "Also specify s-p-d contributions separately (1 for yes or 2 for no) ?"


  read ( * , * ) choice_spdf


  ! test if the answer is a valid option


  if ( choice_spdf /= 1 .and. choice_spdf /= 2) then

   
     write ( * , * ) "Wrong answer"

   
     goto 1111


  end if 


  ! go to line corresponding to good label


  ! ignore first total dos lines


  do i = 1 , number_of_steps ; read ( 10 , * ) ; end do


  ! ignore non wanted atoms


  if ( choice_atom > 1 ) then


    do i = 1 , ( choice_atom - 1) * ( number_of_steps + 1 ) ; read ( 10 , * ) ; end do


  end if


  ! read useless line


  read ( 10 , * )


  ! for spin polarised s+ s- px+ px- py+ py- pz+ pz- dxy+ dxy- dyz+ dyz- dxz+ dxz- dx2-y2+ dx2-y2- dz2+ dz2-

  ! for non-spin polarized s px py pz dxy dyz dxz dx2-y2 dz2


    if ( choice_spdf == 1 ) then

	
      if ( choice_polarization == 1 ) then

	
	write ( 11 , * ) "Energy - DOS_UP - DOS_DOWN - DOS_TOT - DOS_UP_s - DOS_DOWN_s - DOS_TOT_s - &

                          &DOS_UP_p - DOS_DOWN_p - DOS_TOT_p - DOS_UP_d - DOS_DOWN_d - DOS_TOT_d"


	do i = 1 , number_of_steps      	

	
	  read(10,*) t, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16, n17, n18

	  
	  dos_up = n1+n3+n5+n7+n9+n11+n13+n15+n17


	  dos_down = n2+n4+n6+n8+n10+n12+n14+n16+n18


	  dos_up_s = n1


	  dos_down_s = n2


	  dos_up_p = n3+n5+n7


	  dos_down_p = n4+n6+n8


	  dos_up_d = n9+n11+n13+n15+n17


	  dos_down_d = n10+n12+n14+n16+n18


	  if ( choice_translate == 1 ) write ( 11 , * ) t - Efermi, dos_up, dos_down, dos_up+dos_down,&

	  & dos_up_s, dos_down_s, dos_up_s+dos_down_s, dos_up_p, dos_down_p, dos_up_p+dos_down_p,&

	  & dos_up_d, dos_down_d, dos_up_d+dos_down_d


	  if ( choice_translate == 2 ) write ( 11 , * ) t, dos_up, dos_down, dos_up+dos_down,&

	  & dos_up_s, dos_down_s, dos_up_s+dos_down_s, dos_up_p, dos_down_p, dos_up_p+dos_down_p,&

	  & dos_up_d, dos_down_d, dos_up_d+dos_down_d


	end do

      
      else if ( choice_polarization == 2 ) then


	write ( 11 , * ) "Energy - DOS_TOT - DOS_s - DOS_p - DOS_d"

	
	do i = 1 , number_of_steps      	

	
	  read(10,*) t, n1, n2, n3, n4, n5, n6, n7, n8, n9


	  dos_up = n1+n2+n3+n4+n5+n6+n7+n8+n9


	  dos_up_s = n1


	  dos_up_p = n2+n3+n4


	  dos_up_d = n5+n6+n7+n8+n9


	  if ( choice_translate == 1 ) write ( 11 , * ) t - Efermi, dos_up, dos_up_s, dos_up_p, dos_up_d


	  if ( choice_translate == 2 ) write ( 11 , * ) t, dos_up, dos_up_s, dos_up_p, dos_up_d


	end do


      end if  

    
    else if ( choice_spdf == 2 ) then

    
      if ( choice_polarization == 1 ) then


	write ( 11 , * ) "Energy - DOS_UP - DOS_DOWN - DOS_TOT"


	do i = 1 , number_of_steps

	
	  read(10,*) t, n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12, n13, n14, n15, n16, n17, n18


	  dos_up = n1+n3+n5+n7+n9+n11+n13+n15+n17


	  dos_down = n2+n4+n6+n8+n10+n12+n14+n16+n18


	  if ( choice_translate == 1 ) write ( 11 , * ) t - Efermi, dos_up, dos_down, dos_up+dos_down


	  if ( choice_translate == 2 ) write ( 11 , * ) t, dos_up, dos_down, dos_up+dos_down


	end do


      else if ( choice_polarization == 2 ) then


	write ( 11 , * ) "Energy - DOS_TOT"


	do i = 1 , number_of_steps


      	  read(10,*) t, n1, n2, n3, n4, n5, n6, n7, n8, n9


	  dos_up = n1+n2+n3+n4+n5+n6+n7+n8+n9


	  if ( choice_translate == 1 ) write ( 11 , * ) t - Efermi , dos_up


	  if ( choice_translate == 2 ) write ( 11 , * ) t , dos_up


	end do

      
      end if


    end if


end if


! close DOSCAR


close ( 10 )


! close output file doscar_analysis_spd.dat


close ( 11 )


write ( * , * ) 'done. look at => doscar_analysis_spd.dat'


write ( * , * ) 'To read it, you may use => xmgrace -nxy doscar_analysis_spd.dat'


end program doscar_analysis