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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 第一性原理 » Vasp&MedeA » 【求助】怎么根据PROCAR作DOS和PDOS图
221/3返回列表123下一页
查看: 2822  |  回复: 21
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。

huangyc

木虫 (正式写手)

[交流] 【求助】怎么根据PROCAR作DOS和PDOS图

大家好，计算DOS时得到PROCAR，准备学习画DOS和PDOS，共126个K点，每个K点100 bands，用lev00处理时提示缺少PROCAR或PROCAR错误，用p4vasp也打不开，请帮我看看，谢谢
PROCAR new format
# of k-points:  126         # of bands: 100         # of ions:  24

k-point    1 :    0.04166667 0.04166667 0.04166667     weight = 0.00462963

band   1 # energy  -29.57997136 # occ.  2.00000000

ion      s      p      d      f    tot
  1  0.233  0.000  0.000  0.000  0.233
  2  0.233  0.000  0.000  0.000  0.233
  3  0.233  0.000  0.000  0.000  0.233
  4  0.233  0.000  0.000  0.000  0.233
  5  0.000  0.000  0.000  0.000  0.000
  6  0.000  0.000  0.000  0.000  0.000
  7  0.000  0.000  0.000  0.000  0.000
  8  0.000  0.000  0.000  0.000  0.000
  9  0.000  0.000  0.000  0.000  0.000
10  0.000  0.000  0.000  0.000  0.000
11  0.000  0.000  0.000  0.000  0.000
12  0.000  0.000  0.000  0.000  0.000
13  0.000  0.000  0.000  0.000  0.000
14  0.000  0.000  0.000  0.000  0.000
15  0.000  0.000  0.000  0.000  0.000
16  0.000  0.000  0.000  0.000  0.000
17  0.000  0.000  0.000  0.000  0.000
18  0.000  0.000  0.000  0.000  0.000
19  0.000  0.000  0.000  0.000  0.000
20  0.000  0.000  0.000  0.000  0.000
21  0.000  0.000  0.000  0.000  0.000
22  0.000  0.000  0.000  0.000  0.000
23  0.000  0.000  0.000  0.000  0.000
24  0.000  0.000  0.000  0.000  0.000
tot  0.937  0.001  0.000  0.000  0.938
这是其中一部分，下面还有band2-100，接着就是k-points2-126，很长的文件
不知道大家都怎么处理DOS？学习VASP时间不长，还希望虫友们指点

[ Last edited by SHY31 on 2008-12-4 at 20:41 ]
回复此楼

» 收录本帖的淘帖专辑推荐

VASP

» 猜你喜欢
1楼 2008-11-27 22:17:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

intro987

铁杆木虫 (小有名气)
★ ★
jghe(金币+2,VIP+0):感谢你参与讨论！
我使用DOSCAR处理的，应该和PROCAR是一起出来的结果，里面很多列数据，第一列是能量，后面对应s, p_y, p_z, p_x, d_xy, d_yz, d_z2, d_xz, d_x2-y2轨道（VASP出来的结果应该顺序都是一样的）。数据上下顺序对应POSCAR里面的原子顺序
lev00是处理电荷密度的，p4vasp我只用来打开过POSCAR，没试过PROCAR，可能不支持这种类型吧~
赞一下(10人) 回复此楼
2楼2008-12-04 18:43:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hym5701

木虫 (小有名气)
★ ★
jghe(金币+2,VIP+0):谢谢你提供的资源！
引用回帖:
Originally posted by huangyc at 2008-11-27 22:17:
大家好，计算DOS时得到PROCAR，准备学习画DOS和PDOS，共126个K点，每个K点100 bands，用lev00处理时提示缺少PROCAR或PROCAR错误，用p4vasp也打不开，请帮我看看，谢谢
PROCAR new format
# of k-points:  126   ...

网上不是用vasp的一些tools吗，里面就有一个分析PROCAR的文件，文件名是dos-procar.f，你现将这个文件编译一下，然后就运行这个文件，他会提示你要干什么，运行完后，就可以得到totdos。dat 和ldos。dat，接着用origin等画图就可以了。
赞一下(10人) 回复此楼
3楼2008-12-04 20:22:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

huangyc

木虫 (正式写手)
楼上的这个东东编译不起来啊，有语法错误
哪位再指点，谢谢
赞一下 回复此楼
4楼2008-12-05 11:15:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hym5701

木虫 (小有名气)
★ ★
wuli8(金币+2,VIP+0):感谢您的交流！(*^__^*) ……
引用回帖:
Originally posted by huangyc at 2008-12-5 11:15:
楼上的这个东东编译不起来啊，有语法错误
哪位再指点，谢谢

不可能吧我都编译过啊，你是这么编译的，命令贴出来看看
赞一下(10人) 回复此楼
5楼2008-12-05 14:15:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cuihang

铁杆木虫 (正式写手)

VASPer
★
wuli8(金币+1,VIP+0):Good idea!
不如楼上把网址贴出来
赞一下(10人) 回复此楼
聪明人知道自己是傻子，傻子以为自己是聪明人
6楼2008-12-05 18:21:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wyl04

金虫 (小有名气)
★ ★ ★
wuchenwf(金币+3,VIP+0):xiexie
implicit real*8(a-h,o-z)
      parameter (nkd =900)
      parameter (nbd = 50 )
      parameter (natmd = 20)
      parameter (ned = 1001)
      dimension dump(20),oc(nkd,nbd,natmd,4),eig_ev(nkd,nbd),wt(nkd)
      dimension ee(ned),gpdos(ned),gpdost(ned)
      dimension gpdos_s(ned),gpdos_p(ned),gpdos_d(ned)
      open(7,file='PROCAR')
      pi = 3.141592654

      write(*,*) 'Spin polarized calculation? (no=1,yes=2):'
      read (*,*) ispin
      if ((ispin.ne.1).and.(ispin.ne.2)) then
      write(*,*) ' INPUT ERROR, ispin must equal to 1 or 2 '
      stop
      endif
      if (ispin.eq.1) then
      open(31,file='ldos.dat')
      open(32,file='dos-tot.dat')
      elseif (ispin.eq.2) then
      open(31,file='ldos-up.dat')
      open(32,file='dos-tot-up.dat')
      open(41,file='ldos-dn.dat')
      open(42,file='dos-tot-dn.dat')
      endif

      write(6,*) ' Enter fermi energy: '
      read(5,*) fermi
c     write(6,102) fermi
      read(7,103) dump                    
c     write(6,103) dump                    

      write(6,*) 'which atom (na) you want to plot LDOS:'
      read(5,*) na
      write(6,*) 'enter the gaussian smearing factor:'
      read(5,*) gaussian

      do 9000 isp=1,ispin
      read(7,104) nk,nband,nion
c     write(6,104) nk,nband,nion
      if (nk .gt. nkd) stop ' nk too large '
      emin = 1000.0
      emax = -1000.0
      do 1000 k = 1,nk
      read(7,103) dump
      read(7,105) kp,pt1,pt2,pt3,wt(k)
c     write(6,105) kp,pt1,pt2,pt3,wt(k)
      read(7,103) dump
      do  nb = 1,nband
      read(7,106) nb1,eig_ev(k,nb),occ
c     write(6,106) nb1,eig_ev(k,nb),occ
      eig_ev(k,nb) = eig_ev(k,nb)-fermi
      if (eig_ev(k,nb) .gt. emax) emax = eig_ev(k,nb)
      if (eig_ev(k,nb) .lt. emin) emin = eig_ev(k,nb)
      read(7,103) dump
      read(7,103) dump
c     write(6,*) 'nion=',nion
      niont = nion +1
      if (nion .eq. 1) niont = 1
      do  ion = 1,niont
      read(7,107) (oc(k,nb,ion,j),j=1,4)
c     write(6,107) (oc(k,nb,ion,j),j=1,4)
      enddo
      read(7,103) dump
c     write(6,103) dump
      enddo
1000 continue
      weight = 0.0
      do k = 1, nk
      weight = weight + wt(k)
      enddo
      do k = 1,nk
      wt(k) =  wt(k) / weight
      enddo
c     estart_ev = - 5.0
c     eend_ev = 2.0
      estart_ev = int(emin -1.0)
      eend_ev = int(emax +1.0)
      netot =  701
      de_ev = (eend_ev - estart_ev)/(netot-1)
      do 2000 ne = 1, netot
      ee(ne) = estart_ev + (ne-1) * de_ev
      gpdos(ne) = 0.0d0
      gpdos_s(ne) = 0.0d0
      gpdos_p(ne) = 0.0d0
      gpdos_d(ne) = 0.0d0
      gpdost(ne) = 0.0d0
2000 continue
      ascal = de_ev   / ( gaussian * sqrt(pi) )
      ascal = 1.0     / ( gaussian * sqrt(pi) )
      scal_spin = 2.0
      do 5000 k=1,nk
      do 4000 nb=1,nband
      ddos = scal_spin * oc(k,nb,na,4)*wt(k)
      ddos_s = scal_spin * oc(k,nb,na,1)*wt(k)
      ddos_p = scal_spin * oc(k,nb,na,2)*wt(k)
      ddos_d = scal_spin * oc(k,nb,na,3)*wt(k)
      ddost = scal_spin *wt(k)
      do ne = 1,netot
      dij = ( eig_ev(k,nb) - ee(ne) )**2 / (gaussian**2)
      gpdos(ne) = gpdos(ne) + ascal * ddos * exp(-dij)
      gpdos_s(ne) = gpdos_s(ne) + ascal * ddos_s * exp(-dij)
      gpdos_p(ne) = gpdos_p(ne) + ascal * ddos_p * exp(-dij)
      gpdos_d(ne) = gpdos_d(ne) + ascal * ddos_d * exp(-dij)
      gpdost(ne) = gpdost(ne) + ascal * ddost * exp(-dij)
      enddo
4000 continue
5000 continue
      do 6000 ne = 1,netot
      if (gpdos(ne) .lt. 1.0E-29) gpdos(ne)=0.0
      if (gpdos_s(ne) .lt. 1.0E-29) gpdos_s(ne)=0.0
      if (gpdos_p(ne) .lt. 1.0E-29) gpdos_p(ne)=0.0
      if (gpdos_d(ne) .lt. 1.0E-29) gpdos_d(ne)=0.0
      if (gpdost(ne) .lt. 1.0E-29) gpdost(ne)=0.0
      write(21+10*isp,108) ee(ne),gpdos(ne),gpdos_s(ne),gpdos_p(ne)
     &       ,gpdos_d(ne)  
      write(22+10*isp,108) ee(ne),gpdost(ne)
6000 continue
9000 continue
101  format(10x,f9.5)
102  format(f10.5)
103  format(20a4)
104  format(16x,i3,20x,i5,19x,i4)
105  format(10x,i3,5x,3f11.8,13x,f11.8)
106  format(4x,i4,9x,f14.8,7x,f12.8)
107  format(3x,4f7.3)
108  format(f10.5,4e12.4)
      stop
      end
赞一下(10人) 回复此楼
7楼2008-12-05 18:57:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wyl04

金虫 (小有名气)
★ ★
fegg7502(金币+2,VIP+0):thanks
procar->bands

      implicit real*8(a-h,o-z)
      parameter (nbd = 200)
      parameter (nkd = 500)
      parameter (nxd = 300)
      parameter (natmd = 20)
      dimension a(3,3),b(3,3),c(3),e(nkd,nbd),sk(nkd,3)
      dimension xx(nxd) ,wei(nkd)
      dimension dump(20),oc(nkd,nbd,natmd,4)
      open(7,file='PROCAR',form='FORMATTED',status='OLD')
      pi = 3.141592654
      read(7,103) dump                    

      write(*,*) 'Spin polarized calculation? (no=1,yes=2):'
      read (*,*) ispin
      if ((ispin.ne.1).and.(ispin.ne.2)) then
      write(*,*) ' INPUT ERROR, ispin must equal to 1 or 2 '
      stop
      endif

c     write(*,*) 'Enter # of interval (npoints) and division (ndiv):'
c     read (*,*) npoints,ndiv
      open(9,file='KPOINTS',form='FORMATTED',status='OLD')
      read(9,100) temp
      read(9,*) ndiv

      write(*,*) 'Enter the range of energy to plot:'
      read (*,*) er1,er2
      emin=min(er1,er2)
      emax=max(er1,er2)
      write(*,*) 'Enter the value of fermi energy:'
      read(*,*)  ef

      if (ispin.eq.1) then
      open(11,file='band.dat')
      elseif (ispin.eq.2) then
      open(11,file='band-up.dat')
      open(12,file='band-dn.dat')
      endif

      open(8,file='POSCAR',form='FORMATTED',status='OLD')
      read(8,100) temp
c      write(6,100) temp
100   format(20a4)
      read (8,*) aa
c      WRITE(6,*) aa
c   
c      *** read lattice constant from POSCAR**
c            
      do i=1,3
         read (8,*) (a(i,j),j=1,3)
c       WRITE(6,500) (a(i,j),j=1,3)
500   format (3f12.8)
      enddo
      do i=1,3
         do j=1,3
         a(i,j)=aa*a(i,j)
         enddo
c       WRITE(6,500) (a(i,j),j=1,3)
      enddo
c
c     *** read lattice vector from POSCAR***
c
      volume=a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)
     &+a(1,3)*a(2,1)*a(3,2)-a(1,1)*a(2,3)*a(3,2)
     &-a(1,2)*a(2,1)*a(3,3)-a(1,3)*a(2,2)*a(3,1)
      do i=1,3
          if (i .eq. 1) then
            j=2
            k=3
          else if (i .eq. 2) then
            j=3
            k=1
          else
            j=1
            k=2
          endif
        c(1)=a(j,2)*a(k,3)-a(j,3)*a(k,2)
        c(2)=a(j,3)*a(k,1)-a(j,1)*a(k,3)
        c(3)=a(j,1)*a(k,2)-a(j,2)*a(k,1)
        do j=1,3
           b(i,j)=2*pi*c(j)/volume
c         WRITE (6,*) b(i,j)
        enddo
       enddo


      do 9000 isp=1,ispin

      read(7,104) nk,nband,nion
      do 1000 k = 1,nk
      read(7,103) dump
      read(7,105) kp,(sk(k,j),j=1,3),wei(k)
c     write(6,105) kp,(sk(k,j),j=1,3),wei(k)
      read(7,103) dump
      do  nb = 1,nband
      read(7,106) nb1,e(k,nb),occ
c     write(6,106) nb1,e(k,nb),occ
      read(7,103) dump
      read(7,103) dump
c     write(6,*) 'nion=',nion
      niont = nion +1
      if (nion .eq. 1) niont = 1
      do  ion = 1,niont
      read(7,107) (oc(k,nb,ion,j),j=1,4)
c     write(6,107) (oc(k,nb,ion,j),j=1,4)
      enddo
      read(7,103) dump
c     write(6,103) dump
      enddo
1000 continue

      weight = 0.0
      do k = 1, nk
      weight = weight + wei(k)
      enddo

      do k = 1,nk
      wei(k) =  wei(k) / weight
      enddo

101  format(10x,f9.5)
102  format(f10.5)
103  format(20a4)
104  format(16x,i3,20x,i5,19x,i4)
105  format(10x,i3,5x,3f11.8,13x,f11.8)
106  format(4x,i4,9x,f14.8,7x,f12.8)
107  format(3x,4f7.3)


c
c     *** find reciprocal lattice vector ***
      xx(1) = 0.0
      nn = 1
      do k = 1,nk-1
      dkx=(sk(k+1,1)-sk(k,1))*b(1,1) + (sk(k+1,2)-sk(k,2))*b(2,1)
     &   + (sk(k+1,3)- sk(k,3))*b(3,1)
      dky=(sk(k+1,1)-sk(k,1))*b(1,2) + (sk(k+1,2)-sk(k,2))*b(2,2)
     &   + (sk(k+1,3)- sk(k,3))*b(3,2)
      dkz=(sk(k+1,1)-sk(k,1))*b(1,3) + (sk(k+1,2)-sk(k,2))*b(2,3)
     &   + (sk(k+1,3)- sk(k,3))*b(3,3)
      del =  sqrt ( dkx**2 + dky**2 + dkz**2 )
      nn = nn +1
      xx(nn) = xx(nn-1) + del
      enddo

      do n=1,nband
        if (mod(n,2).ne.0) then
         do k=1,nk
          ee = e(k,n) - ef
          if ( ee .gt. emax ) ee = emax
          if ( ee .lt. emin ) ee = emin
         write (10+isp,300) xx(k),ee
         enddo
        elseif (mod(n,2).eq.0) then
         do i=nk,1,-1
          ee = e(i,n) - ef
          if ( ee .gt. emax ) ee = emax
          if ( ee .lt. emin ) ee = emin
          write (10+isp,300) xx(i),ee
         enddo
        endif
      enddo
300   format (f12.8,2x,f12.8)

      if (mod(nband,2) .ne. 0) then
          write (10+isp,300) xx(nk),emin
          write (10+isp,300) xx(1),emin
      else
          write (10+isp,300) xx(1),emin
      endif
c
c     *** write xx-ee ***
c
        npoints=nk/ndiv
       do n=2,npoints
         kk=(n-1)*ndiv
        write (10+isp,300) xx(kk),emin
        write (10+isp,300) xx(kk),emax
        write (10+isp,300) xx(kk),emin
      enddo
        write (10+isp,300) xx(nk),emin
        write (10+isp,300) xx(nk),emax
        write (10+isp,300) xx(1),emax
        write (10+isp,300) xx(1),emin
        zero=0.0
        write (10+isp,300) xx(1),zero
        write (10+isp,300) xx(nk),zero

9000  continue

      stop
      end
赞一下(10人) 回复此楼
8楼2008-12-05 19:03:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

huangyc

木虫 (正式写手)
我用
>ifort -o dos-procar dos-procar.f90编译的。楼上的两个版本我都编译不通啊
赞一下 回复此楼
9楼2008-12-07 19:38:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hym5701

木虫 (小有名气)
★ ★
zdhlover(金币+2,VIP+0):非常感谢，欢迎常来
引用回帖:
Originally posted by huangyc at 2008-12-7 19:38:
我用
>ifort -o dos-procar dos-procar.f90编译的。楼上的两个版本我都编译不通啊

奇怪了，你为什么要将f改成f90呢，改了当然就错了下面是我编译的信息，没有任何错误信息你看看吧
[root@localhost ~]# ls
anaconda-ks.cfg  Desktop  dos-procar.f  install.log  install.log.syslog  xcrys_tmp
[root@localhost ~]# ifort -o dos-procar dos-procar.f
[root@localhost ~]# ls
anaconda-ks.cfg  Desktop  dos-procar  dos-procar.f  install.log  install.log.syslog  xcrys_tmp
[root@localhost ~]# ./dos-procar
Spin polarized calculation? (no=1,yes=2):
1
  Enter fermi energy:
.......................................
.......

[ Last edited by hym5701 on 2008-12-7 at 21:35 ]
赞一下(10人) 回复此楼
10楼2008-12-07 21:34:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 huangyc 的主题更新
221/3返回列表123下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭