24小时热门版块排行榜

返回列表

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

jpchou

专家顾问 (著名写手)

[资源] dos-procar.f90

先說一下~小弟台灣人，沒法打簡體中文
用字遣詞上也許會有些不同~若有看不懂之處~請不用客氣地提出~謝謝

之前有討論過 PROCAR 內容 (http://muchong.com/bbs/viewthread.php?tid=6657343)
了解 PROCAR 之後要寫程式將 PROCAR 轉化為可以畫圖的檔案~就變得不是那麼的困難

我個人是覺得~如果是碩士班學生~沒有打算進修博士班
如何寫程式就變得不是那麼的重要
可是如果打算要讀博班~或現在就是博班的學生
在計算物理化學領域~不會任何程式會讓自己的未?碜兊暮苄量�
fortran, C, perl, python, shell script, ..... 能夠熟悉其中任何一個都是好
但一個都不會的話.... 嘿嘿... 以後的日子會很不好過喔

我還是想要再說一下
就我所知~ dos-procar.f 這程式最源頭是台灣的中正大學梁�?全老師
但我不能肯定這件事~若有爭議請提出 ^_^"

以下
我會將我自己的 dos-procar.f90 程式碼貼出?� (修改自 dos-procar.f)
盡可能說明裡面內容
(ps. 每一行的開頭數字是行號~跟程式本身無關)
程式用途是:
  1) 讀取 PROCAR (LORBIT=10) and OUTCAR
  2) 產生 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 : > 開頭前六行的第一個字元都是 "!"，代表前六行是"說明"
  3 !    INPUT FILES : OUTCAR, PROCAR
  4 !
  5 ! Modified by JPCHOU, Nov. 2013
  6 !===========================================================================
  7    implicit real*8(a-h,o-z)  > FORTRAN在使用變數時必須要宣告，此行是說所有開頭是 a, b, c, ...h, o, p, q, ....z 的單字之變數~都是 real*8
  8    character*64 tmp
  9    dimension ddosP(4),ddosT(4)  > ddosP and ddosT 分別是 "特定原子的DOS" 與 "整個系統的 DOS "
                                                   > 各有四個~分別代表 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 > 預期將能量分成 1000 點
13    allocate(ee(netot),gdosP(4,netot),gdosT(4,netot))
14    open(1,file='OUTCAR')
15 !=================================================================
16    nline = 0 ;  open(99,status='scratch')  > 這一段是去讀取 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 這個檔是個暫存檔~不會真的存在~只暫存在記憶體(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 是程式與人之間的對話
26    write(6,'("enter the gaussian "

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

') ; read(5,*) ispin > 告訴程式此系統是否有計算 spin, 其�?這資�?也可以?� OUTCAR 去讀就好
28 write(6,'("How many atoms are integrated? "

') ;  read(5,*) Nna >告知程式 "特定原子" 總共有幾顆
29    allocate(na(Nna))
30    write(6,*) '  atom number :'                ;  read(5,*) na  >告知程式 "特定原子" 是那幾顆 (給編號~對應於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')  > 如果系統沒有計算SPIN
35       open(12,file='gdosP-p.dat') ; open(52,file='gdosT-p.dat')  > 則輸出檔為 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 (這是整個系統的 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')  > 若是系統有考慮 spin-polarization，則輸出檔就分 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                                                 >  開始讀取 PROCAR，分兩次讀， spin-up and spin-down
53       if(is .eq. 1) read(7,'(a64)') tmp                   > 把沒有用到的內容都用 tmp 去讀取掉
54       read(7,'(16x,i3,20x,i4,19x,i4)') nk,nband,nion  > 得知K點總數(nk)~BAND總數(nband)~與原子總數(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)  > 讀取每一個K點在倒空間座標與權重
65          read(7,'(a64)') tmp
66          do nb=1,nband
67          read(7,'(8x,9x,f14.8,7x,f12.8)') eigen(k,nb),occ  > 讀取每個 eigenvalue 與 occupation
68          eigen(k,nb) = eigen(k,nb) - fermi                   > 能量平移 E-fermi，使 fermi level 落在 0 點
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

對不起
我累了..............................................................
我發覺我個人?U言太多~程式每一行都想要說些什麼
但這是個很沒效率的做法
且版面會變很?y
所以我覺得.......

對程式內容有興趣的人請自行提出
不然.... 就直接用就好了
每一行都有行號
對任一行有問題的人就把指定行號說出?�
我再一一回應
這樣似乎是比較有效率的做法

或是對於原程式 dos-procar.f 有疑問者
也可以提出?韣大家討論討論

24小时热门版块排行榜

[资源] dos-procar.f90

» 本帖附件资源列表

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

» 本帖已获得的红花（最新10朵）

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐