应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 第一性原理 » Vasp&MedeA » 用dos_procar处理dos数据不准
北京石油化工学院2026年研究生招生接收调剂公告
41/1返回列表
查看: 1620  |  回复: 3
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2014-05-02 10:38:18
已阅   关注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 的主题更新
41/1返回列表
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 287分求调剂 有专利国奖一志愿哈工大085406 +3 白易辰 2026-04-06 4/200 2026-04-06 15:24 by going home
[考研] 302分求调剂 一志愿安徽大学085601 +9 zyx上岸！ 2026-04-04 9/450 2026-04-06 13:05 by 猪会飞
[考研] 319分085702安全工程求调剂 +4 rious 2026-04-05 4/200 2026-04-06 11:08 by 静静静静静静静�
[考研] 一志愿北京交通大学材料工程总分358求调剂 +4 cs0106 2026-04-04 4/200 2026-04-05 18:46 by imissbao
[考研] 一志愿同济大学323分（080500）求调剂 +8 yikeniu 2026-04-01 8/400 2026-04-05 18:15 by cql1109
[考研] 一志愿上海海洋大学083200食品学硕，求调剂，接受其他专业083200 +4 what张 2026-04-04 5/250 2026-04-05 14:07 by chw1980_0
[考研] 26考研调剂0710 0860 +9 补补不补 2026-04-03 14/700 2026-04-04 23:32 by 果冻大王
[考研] 085400电子信息319求调剂（接受跨专业调剂） +5 星星不眨眼喽 2026-04-03 6/300 2026-04-04 21:50 by hemengdong
[考研] 349求调剂 +11 zwjjjjjj 2026-03-31 11/550 2026-04-04 19:52 by 蓝云思雨
[考研] 调剂 +4 是可乐不是可乐 2026-04-04 4/200 2026-04-04 19:41 by 唐沐儿
[考研] 309分085801求调剂 +11 MY_angel 2026-03-31 11/550 2026-04-04 19:11 by 蓝云思雨
[考研] 322求调剂 +4 FZAC123 2026-04-03 4/200 2026-04-03 20:55 by zhq0425
[考研] 071000生物学调剂 +8 知昭蔓 2026-04-02 8/400 2026-04-03 10:36 by macy2011
[考研] 309求调剂 +14 呆菇不是戴夫 2026-04-02 14/700 2026-04-03 09:42 by 蓝云思雨
[考研] 学硕化学工程与技术，一志愿中国海洋大学320+求调剂 +8 披星河 2026-04-02 8/400 2026-04-02 14:12 by oooqiao
[考研] 材料专硕322分 +11 哈哈哈吼吼吼哈 2026-04-01 11/550 2026-04-02 10:52 by lnilvy
[考研] 省双一流重点一本大学招收调剂 +4 wwwwffffff 2026-03-31 7/350 2026-04-01 15:23 by wwwwffffff
[考研] 070300化学专业279调剂 +10 哈哈哈^_^ 2026-03-31 10/500 2026-03-31 23:13 by liu823948201
[考研] 本2一志愿C9-333分，材料科学与工程，求调剂 +9 升升不降 2026-03-31 9/450 2026-03-31 18:01 by 无际的草原
[考研] 物理学调剂 +4 小羊36 2026-03-30 4/200 2026-03-31 16:16 by lishahe
信息提示
关闭
请填处理意见
关闭