24小时热门版块排行榜    

查看: 1603  |  回复: 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的回帖

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的回帖
查看全部 4 个回答

mywai520

铁杆木虫 (著名写手)


【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +1
dx0620: 金币+5, 有帮助 2014-05-04 09:29:41
换一个处理脚本吧,最好用p4vasp,直观图形化。若是不想也可以用vaspkit,所有的处理打包,是个不错的选择。
2楼2014-05-02 12:24:18
已阅   关注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的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 材料与化工304求B区调剂 +6 邱gl 2026-03-11 7/350 2026-03-17 01:27 by 学员FnSWZj
[考研] [复试调剂]西南科技大学国防/材料导师推荐 +3 尖角小荷 2026-03-16 6/300 2026-03-16 23:21 by 尖角小荷
[考研] 302求调剂 +3 小贾同学123 2026-03-15 5/250 2026-03-16 20:39 by zhq0425
[考研] 328求调剂,英语六级551,有科研经历 +3 生物工程调剂 2026-03-16 4/200 2026-03-16 20:13 by Wangjingyue
[考研] 333求调剂 +3 文思客 2026-03-16 7/350 2026-03-16 18:21 by 文思客
[考研] 0854控制工程 359求调剂 可跨专业 +3 626776879 2026-03-14 9/450 2026-03-16 17:42 by 626776879
[考研] 304求调剂 +5 素年祭语 2026-03-15 5/250 2026-03-16 17:00 by 我的船我的海
[考研] 321求调剂 +5 大米饭! 2026-03-15 5/250 2026-03-16 16:33 by houyaoxu
[考研] 0703一志愿211 285分求调剂 +5 ly3471z 2026-03-13 5/250 2026-03-16 16:16 by 哦哦123
[考研] 277材料科学与工程080500求调剂 +3 自由煎饼果子 2026-03-16 3/150 2026-03-16 14:10 by 运气yunqi
[考研] 求老师收留调剂 +4 jiang姜66 2026-03-14 5/250 2026-03-15 20:11 by Winj1e
[考研] 0856专硕279求调剂 +5 加油加油!? 2026-03-15 5/250 2026-03-15 11:58 by 2020015
[基金申请] 现在如何回避去年的某一个专家,不知道名字 +3 zk200107 2026-03-12 6/300 2026-03-14 17:13 by zk200107
[考研] 331求调剂(0703有机化学 +5 ZY-05 2026-03-13 6/300 2026-03-14 10:51 by Jy?
[考研] 求材料调剂 +5 隔壁陈先生 2026-03-12 5/250 2026-03-13 22:03 by 星空星月
[考研] 314求调剂 +7 无懈可击的巨人 2026-03-12 7/350 2026-03-13 15:40 by JourneyLucky
[考研] 289求调剂 +3 李政莹 2026-03-12 3/150 2026-03-13 11:02 by 求调剂zz
[考研] 274求调剂0856材料化工 +12 z2839474511 2026-03-11 13/650 2026-03-13 10:39 by peike
[考研] 工科0856专硕化学工程269能调剂吗 +10 我想读研11 2026-03-10 10/500 2026-03-13 10:14 by Yuyi.
[考研] 420求调剂 +4 莫向外求11 2026-03-10 6/300 2026-03-12 14:41 by ruiyingmiao
信息提示
请填处理意见