| 查看: 2480 | 回复: 21 | ||||
| 当前主题已经存档。 | ||||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||||
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 |
» 猜你喜欢
博士读完未来一定会好吗
已经有27人回复
-
投稿精细化工
已经有4人回复
-
高职单位投计算机相关的北核或SCI四区期刊推荐，求支招！
已经有4人回复
-
导师想让我从独立一作变成了共一第一
已经有9人回复
-
到新单位后，换了新的研究方向，没有团队，持续积累2区以上论文，能申请到面上吗
已经有11人回复
-
读博
已经有4人回复
-
JMPT 期刊投稿流程
已经有4人回复
-
心脉受损
已经有5人回复
-
Springer期刊投稿求助
已经有4人回复
-
小论文投稿
已经有3人回复
★ ★ ★
wuchenwf(金币+3,VIP+0):xiexie
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 |
7楼2008-12-05 18:57:53
★ ★
fegg7502(金币+2,VIP+0):thanks
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 |
8楼2008-12-05 19:03:18
11楼2008-12-08 14:25:29