doscar_analysis_spd

! 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