24小时热门版块排行榜    

查看: 1747  |  回复: 3
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

dx0620

木虫 (正式写手)

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

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

用dos_procar处理dos数据不准

我的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
回复此楼
已阅   关注TA 给TA发消息 送TA红花 TA的回帖

dx0620

木虫 (正式写手)

引用回帖:
2楼: Originally posted by mywai520 at 2014-05-02 12:24:18
换一个处理脚本吧,最好用p4vasp,直观图形化。若是不想也可以用vaspkit,所有的处理打包,是个不错的选择。

您好,用p4vasp的话是怎么处理啊?我用vaspkit处理的话,ipdos是什么啊?为什么感觉tdos跟pdos里面s、p、d加起来差好多,这几个都要除以原子数吗?!谢谢~
3楼2014-05-02 15:03:31
已阅   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 4 个回答

mywai520

铁杆木虫 (著名写手)


【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +1
dx0620: 金币+5, 有帮助 2014-05-04 09:29:41
换一个处理脚本吧,最好用p4vasp,直观图形化。若是不想也可以用vaspkit,所有的处理打包,是个不错的选择。
2楼2014-05-02 12:24:18
已阅   关注TA 给TA发消息 送TA红花 TA的回帖

mywai520

铁杆木虫 (著名写手)


引用回帖:
3楼: Originally posted by dx0620 at 2014-05-02 15:03:31
您好,用p4vasp的话是怎么处理啊?我用vaspkit处理的话,ipdos是什么啊?为什么感觉tdos跟pdos里面s、p、d加起来差好多,这几个都要除以原子数吗?!谢谢~...

你给出的是每一个原子的spd态,要全部原子的加起来才是总的
4楼2014-05-02 16:50:24
已阅   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 2026国自然函评费到账 +3 羊腰板 2026-08-21 3/150 2026-08-21 19:42 by kulium
[基金申请] 人气不行了 +7 fansofjerry 2026-08-21 7/350 2026-08-21 14:13 by Vivilian
[基金申请] 感觉是下周放榜了 +7 angus9576 2026-08-17 12/600 2026-08-21 13:38 by weiyin
[基金申请] 只有每年这种时候来逛逛小木虫 +23 yaoyewhu2008 2026-08-20 25/1250 2026-08-21 12:19 by Poppy1104
[基金申请] 什么时候开奖? +7 CrisMessi 2026-08-18 7/350 2026-08-21 12:15 by weiyin
[基金申请] 今天基金会出结果吗?20260819 +15 kkkl_v 2026-08-19 16/800 2026-08-21 11:22 by 365372687
[基金申请] 应该是下周三26日公布了吧? +4 哈哈蛤? 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[论文投稿] 投稿咨询 +5 wwm09 2026-08-17 7/350 2026-08-21 10:11 by 期刊论文帮手
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
[基金申请] 放榜前的不淡定 40+4 snowwithsea 2026-08-19 13/650 2026-08-21 08:30 by 北京莱茵编辑
[基金申请] 时间戳今天,20号变了 +4 archvillain 2026-08-20 4/200 2026-08-21 08:10 by gatelove
[基金申请] 估计是周四 +3 archvillain 2026-08-18 3/150 2026-08-21 01:48 by jnhyjjm
[基金申请] 今天系统多次维护,明天很可能放榜! +10 zju2000 2026-08-16 11/550 2026-08-20 20:02 by cl479861084
[基金申请] 时间戳变了,能看出什么问题? +18 基诺咪客 2026-08-17 23/1150 2026-08-20 17:19 by Godzela
[教师之家] 为什么余额宝的年化利率越来越低?主要原因有哪些? +5 瞬息宇宙 2026-08-15 5/250 2026-08-19 20:23 by super2002521
[基金申请] 朋友圈看到的 +6 wangzilk 2026-08-18 8/400 2026-08-19 10:55 by Haru815
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
信息提示
请填处理意见