24小时热门版块排行榜

返回列表

【奖励】本帖被评价64次，作者jpchou增加金币 50.8 个

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

jpchou

专家顾问 (著名写手)

[资源] dos-procar.f90

???f???~С?????????]?????w????
??????~????S????Щ???~???п???????~???????????~?x?x

???????^ PROCAR ???? (http://muchong.com/bbs/viewthread.php?tid=6657343)
??? PROCAR ??????????? PROCAR ?D?????????D??n??~????ò??????N?????y

????????X??~?????T????W??~?]?д????M??????
???????????ò??????N?????
?????????????x????~??F?????????W??
????????????W?I??~?????κγ??????????δ??????????
fortran, C, perl, python, shell script, ..... ???????????κ?????????
??????????????.... ???... ????????????????^?

???????????f???
???????~ dos-procar.f ?@???????^?????????????W????????
??????????@????~???Р??h???? ^_^"

????
???????????? dos-procar.f90 ????a?N???? (????? dos-procar.f)
?M?????f???e?????
(ps. ???е??_?^?????????~?????????o?P)
????????:
  1) ?x? PROCAR (LORBIT=10) and OUTCAR
  2) ?a?? total density of state (TDOS) and partial density of state (PDOS) of s-, p-, d-orbital
  3) ??????????? DOS (LDOS)

  1 !===========================================================================
  2 !  Program dos-procar.f90 : > ?_?^????е???????????? "!"?????????????"?f??"
  3 !    INPUT FILES : OUTCAR, PROCAR
  4 !
  5 ! Modified by JPCHOU, Nov. 2013
  6 !===========================================================================
  7    implicit real*8(a-h,o-z)  > FORTRAN?????????r???????棬???????f?????_?^?? a, b, c, ...h, o, p, q, ....z ?????????~???? real*8
  8    character*64 tmp
  9    dimension ddosP(4),ddosT(4)  > ddosP and ddosT ??e?? "???????DOS" ?c "??????y?? DOS "
                                                   > ???????~??e???? s-, p-, d-orbital, and total
10    allocatable wt(

,oc(:,:,:,

,eigen(:,

,ee(

,na(

> wt: weight, oc: occupations, eigen: eigenvalue, ee: energy, na: number of atom
11 allocatable gdosP(:,

,gdosT(:,

12    pi=4.0d0*atan(1.0d0)  ;  netot = 1001 > ?A?????????? 1000 ?c
13    allocate(ee(netot),gdosP(4,netot),gdosT(4,netot))
14    open(1,file='OUTCAR')
15 !=================================================================
16    nline = 0 ;  open(99,status='scratch')  > ?@???????x? OUTCAR???K??? E-fermi ?@???
17  111  read(1,'(a64)',end=222) tmp          > ????????? E-fermi ????????? energy ?????? E-fermi
18    if(tmp(2:8) .eq. 'E-fermi') nline=nline+1  > open 99 ?@???n???????n~??????????~???????????w(memory, ????)??
19    if(tmp(2:8) .eq. 'E-fermi') write(99,'(a64)') tmp
20    go to 111
21  222  rewind(99)
22    do i=1,nline ; read(99,*) tmp,tmp,fermi ; end do
23    close(99)
24 !=================================================================
25    write(6,'("fermi energy = ",f10.5)') fermi    > 25~30 ?????c????g????
26    write(6,'("enter the gaussian "

') ; read(5,*) gaussian > ?o??? Gaussian function ?? width
27 write(6,'("ISPIN(1/2) = ? "

') ; read(5,*) ispin > ???V???????y???????? spin, ????@?Y????????? OUTCAR ??x???
28 write(6,'("How many atoms are integrated? "

') ;  read(5,*) Nna >?????? "??????" ????????w
29    allocate(na(Nna))
30    write(6,*) '  atom number :'                ;  read(5,*) na  >?????? "??????" ???ǎ??w (?o???~?????POSCAR??????)
31
32    open( 7,file='PROCAR')
33    if(ispin .eq. 1) then
34       open(11,file='gdosP-s.dat') ; open(51,file='gdosT-s.dat')  > ?????y?]?????SPIN
35       open(12,file='gdosP-p.dat') ; open(52,file='gdosT-p.dat')  > ?t????n?? gdosP-s, p, d, A.dat  (?@????????? DOS)
36       open(13,file='gdosP-d.dat') ; open(53,file='gdosT-d.dat')  >                gdosT-s, p, d, A.dat (?@????????y?? DOS)
37       open(14,file='gdosP-A.dat') ; open(54,file='gdosT-A.dat') >    A: all
38    else if(ispin .eq. 2) then
39       open(11,file='gdosP-up-s.dat') ; open(51,file='gdosT-up-s.dat')  > ??????y?п??] spin-polarization???t????n??? up and down
40       open(12,file='gdosP-up-p.dat') ; open(52,file='gdosT-up-p.dat')
41       open(13,file='gdosP-up-d.dat') ; open(53,file='gdosT-up-d.dat')
42       open(14,file='gdosP-up-A.dat') ; open(54,file='gdosT-up-A.dat')
43
44       open(21,file='gdosP-dn-s.dat') ; open(61,file='gdosT-dn-s.dat')
45       open(22,file='gdosP-dn-p.dat') ; open(62,file='gdosT-dn-p.dat')
46       open(23,file='gdosP-dn-d.dat') ; open(63,file='gdosT-dn-d.dat')
47       open(24,file='gdosP-dn-A.dat') ; open(64,file='gdosT-dn-A.dat')
48    else
49       write(6,*) 'Input ERROR, you must input 1 or 2 ' ; stop
50    end if
51 !===========================================================================
52    do is=1,ispin                                                 >  ?_??x? PROCAR?? ?????x?? spin-up and spin-down
53       if(is .eq. 1) read(7,'(a64)') tmp                   > ??]????????????? tmp ??x???
54       read(7,'(16x,i3,20x,i4,19x,i4)') nk,nband,nion  > ???K?c????(nk)~BAND????(nband)~?c??????(nion)
55       if(is .eq. 1) then
56          write(6,'("number of kpoints: ",i4)') nk
57          write(6,'("number of bands: ",i4)') nband
58          write(6,'("number of ions: ",i4)') nion
59       end if
60       emin = 9999.0  ;  emax= -9999.0
61       allocate(wt(nk),eigen(nk,nband),oc(nk,nband,nion+1,4))
62       do k=1,nk
63          read(7,'(a64)') tmp
64          read(7,'(10x,i3,5x,3f11.8,13x,f11.8)') kp,pt1,pt2,pt3,wt(k)  > ?x?????K?c??????g?????c????
65          read(7,'(a64)') tmp
66          do nb=1,nband
67          read(7,'(8x,9x,f14.8,7x,f12.8)') eigen(k,nb),occ  > ?x???? eigenvalue ?c occupation
68          eigen(k,nb) = eigen(k,nb) - fermi                   > ??????? E-fermi??? fermi level ???? 0 ?c
69          if(eigen(k,nb) .gt. emax) emax=eigen(k,nb)    > ????????????
70          if(eigen(k,nb) .lt. emin) emin=eigen(k,nb)
71          read(7,'(a64)') tmp ; read(7,'(a64)') tmp
72          niont = nion +1
73          if(nion .eq. 1) niont = 1
74          do ion = 1,niont
75             read(7,'(3x,4f7.3)') (oc(k,nb,ion,j),j=1,4)
76          end do
77          read(7,'(a64)') tmp
78          end do
79       end do
80 ! ##############################################################
81       weight = 0.0 ; gdosP=0.0d0 ; gdosT=0.0d0
82       do k=1,nk ; weight = weight + wt(k) ; end do
83       do k=1,nk ; wt(k) = wt(k) / weight  ; end do
84       estart_ev = int(emin -5.0)
85       eend_ev = int(emax +5.0)
86       de_ev    = (eend_ev - estart_ev)/(netot-1)
87       do ne = 1,netot ; ee(ne) = estart_ev + (ne-1) * de_ev ; end do
88       ascal = 1.0/(gaussian*sqrt(pi))
89 ! ###############################################################
90       do k=1,nk ; do nb=1,nband
91          do i=1,4 ; ddosT(i) =    oc(k,nb,niont,i)*wt(k)  ; end do
92          do i=1,4 ; ddosP(i) = sum(oc(k,nb,na(

,i)*wt(k)) ; end do
93          do ne=1,netot
94          dij = ( eigen(k,nb)-ee(ne) )**2 / (gaussian**2)
95          do i=1,4
96             gdosT(i,ne) = gdosT(i,ne) + ascal*ddosT(i)*exp(-dij)
97             gdosP(i,ne) = gdosP(i,ne) + ascal*ddosP(i)*exp(-dij)
98          end do
99          end do
100       end do ; end do
101 ! ###############################################################
102       if(is .eq. 1) ns=+1 ; if(is .eq. 2) ns=-1
103       do i=1,4 ; write(is*10+i ,1) (ee(j),ns*gdosP(i,j),j=1,netot) ; end do
104       do i=1,4 ; write(is*10+i+40,1) (ee(j),ns*gdosT(i,j),j=1,netot) ; end do
105       deallocate(wt,eigen,oc)
106    end do
107 !===========================================================================
108 1 format(f10.5,f12.4)
109    end

??????
??????..............................................................
??l?X??????U?????~??????ж?????fЩ??N
???@?????]Ч???????
??????????y
???????X??.......

????????????d?????????????
???.... ???????????
???ж??????
???????????}?????????????f????
?????????
?@?????????^??Ч???????

?????????? dos-procar.f ???????
??????????~???????

回复此楼