24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 1621  |  回复: 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

铁杆木虫 (著名写手)


【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +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的回帖

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的回帖
相关版块跳转 我要订阅楼主 dx0620 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 材料考研求调剂总分280 +12 mkjlz1 2026-04-06 14/700 2026-04-06 22:48 by mkjlz1
[考研] 材料334求调剂 +19 Eecho# 2026-04-03 19/950 2026-04-06 08:37 by 小小树2024
[考研] 考研调剂 +5 美丽的youth_ 2026-04-04 6/300 2026-04-06 06:57 by houyaoxu
[考研] 0703求调剂383分 +9 W55j 2026-04-03 9/450 2026-04-06 06:50 by houyaoxu
[考研] 调剂 +3 李广火 2026-04-05 3/150 2026-04-05 18:57 by 蓝云思雨
[考研] 312求调剂 +3 Say Never 2026-04-04 3/150 2026-04-05 11:19 by guoweigw
[考研] 材料调剂 +11 一样YWY 2026-04-02 13/650 2026-04-04 23:10 by 无际的草原
[考研] 可跨专业调剂 +3 周的得地 2026-04-04 6/300 2026-04-04 22:21 by barlinike
[考研] 278求调剂 +14 范婷娜 2026-04-04 15/750 2026-04-04 22:15 by lqwchd
[考研] 333求调剂 +12 wfh030413@ 2026-04-03 13/650 2026-04-04 21:02 by jj987
[考研] 297求调剂 +11 ljy20040718! 2026-04-03 13/650 2026-04-04 09:23 by 来看流星雨10
[考研] 400分求调剂 +3 尴尬且挠头 2026-04-04 3/150 2026-04-04 08:41 by jp9609
[考研] 求调剂,一志愿北京中医药大学 +3 小小达不溜 2026-04-02 3/150 2026-04-03 22:55 by 冲矢昴星团
[考研] 301求调剂 +14 A_JiXing 2026-04-01 14/700 2026-04-03 18:31 by ls刘帅
[考研] 303求调剂 +3 一色清羽 2026-04-02 4/200 2026-04-03 10:22 by 蓝云思雨
[考研] 260求调剂 +3 朱芷琳 2026-04-02 3/150 2026-04-03 08:44 by yulian1987
[考研] 348求调剂 +11 zzzzyk123 2026-04-01 11/550 2026-04-02 16:52 by Wang200018
[考研] 调剂 +3 好好读书。 2026-04-01 6/300 2026-04-02 15:49 by liumengping
[考研] 321求调剂 一志愿 浙江工业大学生物医药 +5 嘿嘿HC 2026-04-01 6/300 2026-04-02 15:23 by sophie2180
[考研] 本2一志愿C9-333分,材料科学与工程,求调剂 +9 升升不降 2026-03-31 9/450 2026-03-31 18:01 by 无际的草原
信息提示
请填处理意见