24小时热门版块排行榜

返回列表

dx0620

木虫 (正式写手)

应助: 14 (小学生)
金币: 4053.6
散金: 5
红花: 6
帖子: 393
在线: 910小时
虫号: 1303311
注册: 2011-05-22
性别: GG
专业: 凝聚态物性 II ：电子结构

[求助] 用dos_procar处理dos数据不准已有1人参与

我的vasp版本是5.3，用的dos_procar处理的PROCAR，想分别得到s，p，d的态密度，但是最后得到的态密度曲线不对，总态密度跟s，p，d的对不上，如下图

我的dos_procar.f90的代码如下

!===========================================================================
!    INPUT FILES : OUTCAR , PROCAR
!===========================================================================
   implicit real*8(a-h,o-z)
   character*64 tmp
   dimension ddosP(4),ddosT(4)
   allocatable wt(

,oc(:,:,:,

,eigen(:,

,ee(

,na(

allocatable gdosP(:,

,gdosT(:,

   pi=4.0d0*atan(1.0d0)  ;  netot = 1001
   allocate(ee(netot),gdosP(4,netot),gdosT(4,netot))
   open(1,file='OUTCAR')
!=================================================================
   nline = 0 ;  open(99,status='scratch')
111  read(1,'(a64)',end=222) tmp
   if(tmp(2:8) .eq. 'E-fermi') nline=nline+1
   if(tmp(2:8) .eq. 'E-fermi') write(99,'(a64)') tmp
   go to 111
222  rewind(99)
   do i=1,nline ; read(99,*) tmp,tmp,fermi ; end do
   close(99)
!=================================================================
   write(6,'("fermi energy = ",f10.5)') fermi
   write(6,'("enter the gaussian "

') ; read(5,*) gaussian
write(6,'("ISPIN(1/2) = ? "

') ; read(5,*) ispin
write(6,'("How many atoms are integrated? "

') ;  read(5,*) Nna
   allocate(na(Nna))
   write(6,*) '  atom number :'                ;  read(5,*) na

   open( 7,file='PROCAR')
   if(ispin .eq. 1) then
      open(11,file='gdosP-s.dat') ; open(51,file='gdosT-s.dat')
      open(12,file='gdosP-p.dat') ; open(52,file='gdosT-p.dat')
      open(13,file='gdosP-d.dat') ; open(53,file='gdosT-d.dat')
      open(14,file='gdosP-A.dat') ; open(54,file='gdosT-A.dat')
   else if(ispin .eq. 2) then
      open(11,file='gdosP-up-s.dat') ; open(51,file='gdosT-up-s.dat')
      open(12,file='gdosP-up-p.dat') ; open(52,file='gdosT-up-p.dat')
      open(13,file='gdosP-up-d.dat') ; open(53,file='gdosT-up-d.dat')
      open(14,file='gdosP-up-A.dat') ; open(54,file='gdosT-up-A.dat')

      open(21,file='gdosP-dn-s.dat') ; open(61,file='gdosT-dn-s.dat')
      open(22,file='gdosP-dn-p.dat') ; open(62,file='gdosT-dn-p.dat')
      open(23,file='gdosP-dn-d.dat') ; open(63,file='gdosT-dn-d.dat')
      open(24,file='gdosP-dn-A.dat') ; open(64,file='gdosT-dn-A.dat')
   else
      write(6,*) 'Input ERROR, you must input 1 or 2 ' ; stop
   end if
!===========================================================================
   do is=1,ispin
      if(is .eq. 1) read(7,'(a64)') tmp
      read(7,'(16x,i3,20x,i4,19x,i4)') nk,nband,nion
      if(is .eq. 1) then
      write(6,'("number of kpoints: ",i4)') nk
      write(6,'("number of bands: ",i4)') nband
      write(6,'("number of ions: ",i4)') nion
      end if
      emin = 9999.0  ;  emax= -9999.0
      allocate(wt(nk),eigen(nk,nband),oc(nk,nband,nion+1,4))
      do k=1,nk
      read(7,'(a64)') tmp
      read(7,'(10x,i3,5x,3f11.8,13x,f11.8)') kp,pt1,pt2,pt3,wt(k)
      read(7,'(a64)') tmp
      do nb=1,nband
         read(7,'(8x,9x,f14.8,7x,f12.8)') eigen(k,nb),occ
         eigen(k,nb) = eigen(k,nb) - fermi
         if(eigen(k,nb) .gt. emax) emax=eigen(k,nb)
         if(eigen(k,nb) .lt. emin) emin=eigen(k,nb)
         read(7,'(a64)') tmp ; read(7,'(a64)') tmp
         niont = nion +1
         if(nion .eq. 1) niont = 1
         do ion = 1,niont
            read(7,'(3x,4f7.3)') (oc(k,nb,ion,j),j=1,4)
         end do
         read(7,'(a64)') tmp
      end do
      end do
! ##############################################################
      weight = 0.0 ; gdosP=0.0d0 ; gdosT=0.0d0
      do k=1,nk ; weight = weight + wt(k) ; end do
      do k=1,nk ; wt(k) = wt(k) / weight  ; end do
      estart_ev = int(emin -5.0)
      eend_ev = int(emax +5.0)
      de_ev    = (eend_ev - estart_ev)/(netot-1)
      do ne = 1,netot ; ee(ne) = estart_ev + (ne-1) * de_ev ; end do
      ascal = 1.0/(gaussian*sqrt(pi))
! ###############################################################
      do k=1,nk ; do nb=1,nband
      do i=1,4 ; ddosT(i) =    oc(k,nb,niont,i)*wt(k)  ; end do
      do i=1,4 ; ddosP(i) = sum(oc(k,nb,na(

,i)*wt(k)) ; end do
      do ne=1,netot
         dij = ( eigen(k,nb)-ee(ne) )**2 / (gaussian**2)
         do i=1,4
            gdosT(i,ne) = gdosT(i,ne) + ascal*ddosT(i)*exp(-dij)
            gdosP(i,ne) = gdosP(i,ne) + ascal*ddosP(i)*exp(-dij)
         end do
      end do
      end do ; end do
! ###############################################################
      if(is .eq. 1) ns=+1 ; if(is .eq. 2) ns=-1
      do i=1,4 ; write(is*10+i ,1) (ee(j),ns*gdosP(i,j),j=1,netot) ; end do
      do i=1,4 ; write(is*10+i+40,1) (ee(j),ns*gdosT(i,j),j=1,netot) ; end do
      deallocate(wt,eigen,oc)
   end do
!===========================================================================
  1 format(f10.5,f12.4)
   end
############################################################################

INCAR 如下：

ENCUT = 400
ISTART = 1
ICHARG = 11
ISMEAR = 1
LORBIT = 11
#EMIN = -25
#EMAX = 22.5
NBANDS = 64
#SIGMA = 0.1
#EDIFF = 1E-6
#EDIFFG = -1E-3
#NSW = 200
#IBRION = 2
#POTIM = 0.2
#ISIF = 3
#ALGO = Normal
#LREAL = Auto
PREC = Accurate
PSTRESS = 500.00

回复此楼