QUANTUM ESPRESSO bands

QUANTUM ESPRESSO：bands

Posted by Xiehua on December 15, 2021

QUANTUM ESPRESSO：bands.x

采用scf计算得到电荷密度和波函数进行bands计算，非自洽（bands）计算能带

设置nscf-bands.in文件,需要修改的参数包括：calculation,nbnd,K_POINTS

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

&CONTROL
    calculation   = "bands"  
    restart_mode  = "from_scratch"
    prefix        = "carbyne"
    outdir        = "./outdir/"
    pseudo_dir    = "./pseudo/"
    verbosity     = "high"
    tprnfor       = .true.  
    tstress       = .true.
    etot_conv_thr =  1.0d-6
/

&SYSTEM
    ibrav       = 6
    celldm(1)   = 18.89726125
    celldm(3)   = 0.256551169
    nat         = 2
    ntyp        = 1
    nbnd        = 22
    occupations = 'fixed'
    ecutwfc     =  50
    ecutrho     =  400
/

&ELECTRONS
    conv_thr         =  1.000e-9
    electron_maxstep =  200
    mixing_beta      =  0.7
    startingpot      = "atomic"
    startingwfc      = "atomic+random"
/

K_POINTS crystal_b
2
0.0 0.0 -0.5  400
0.0 0.0  0.5  1   

ATOMIC_SPECIES
C      12.01070  C.pbe-n-kjpaw_psl.1.0.0.UPF

ATOMIC_POSITIONS (angstrom)
C             5.0000000000       5.0000000000       0.0000000000    0   0   0
C             5.0000000000       5.0000000000       1.2640744811    0   0   1

NSCF-bands计算完成后，使用bands.x程序来处理得到态密度，输入文件bands.in如下

1
2
3
4
5
6

&bands
prefix = 'carbyne'
outdir = './outdir'
filband = 'carbyne.bands.dat'
lsym = .true.
/

对文件carbyne.bands.dat使用Originlab作图，并将能量减去VBM值，使得VBM=0 ，能带结构如下：
bands

计算完能带后使用projwfc.x计算投影能带(Pbands)进行fatband处理，用来得到Wannier拟合的初始函数指定

依次运行pw.x、projwfc.x，在能带计算之后，用projwfc.x生成的fatband.projwfs_up文件，可以和bands.x生成的文件carbyne.bands.dat结合画出各个能带的原子轨道投影，画图脚本如下：

输入文件projwfc.in如下：

1
2
3
4
5
6

&projwfc
prefix = 'carbyne'
outdir = './outdir'
lsym = .false.
filproj = 'fatband'
/

画图脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: yyyu200@163.com
"""

import numpy as np

# plot range for y-axis
ymin=-23
ymax=10
dline=30 # vertical line intervals
lw=0.5 # line width
nband=26 # highest valence band

feig=open('bd.dat')
l=feig.readline()
nbnd=int(l.split(',')[0].split('=')[1])
nks=int(l.split(',')[1].split('=')[1].split('/')[0])

eig=np.zeros((nks,nbnd),dtype=float)
for i in range(nks):
    l=feig.readline()
    count=0
    if nbnd%10==0:
        n=nbnd//10
    else:
        n=nbnd//10+1
    for j in range(n):
        l=feig.readline()
        for k in range(len(l.split())):
            eig[i][count]=l.split()[k]
            count=count+1

feig.close()

import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt

F=plt.gcf()
F.set_size_inches([4,5])
p1=plt.subplot(1, 1, 1)

eig_vbm=max(eig[:, nband-1])
#eig_vbm= 0.00 # fermi energy level in scf output for metals
for i in range(nbnd):
    line1=plt.plot(np.arange(0,nks), eig[:,i]-eig_vbm,color='grey',linewidth=lw )

vline=dline
while vline<nks-1:
    plt.axvline(x=vline, ymin=ymin, ymax=ymax,linewidth=lw,color='black')
    vline=vline+dline

elem=['Sn','O']
ielem=np.array([2,4],dtype=np.int32) # number of atoms for each element
orb=[['s','p','d'],['s','p']]  # projectors for each element

N=len(elem)
iorb=np.zeros([N,],dtype=np.int32)  # number of projectors for each element
for i in range(N):
    iorb[i]=len(orb[i])

lorb=np.zeros([N,],dtype=np.int32) # number of local orbital for each element
for i in range(N):
    for j in orb[i]:
        if j is 's':
            lorb[i]+=1
        elif j is 'p':
            lorb[i]+=3
        elif j is 'd':
            lorb[i]+=5
        elif j is 'f':
            lorb[i]+=7
        else:
            print("unexpect: ",j)
            assert False

nlorb=np.dot(ielem,lorb)

pjsum=np.zeros([nlorb, nks, nbnd], dtype=np.float32)

filproj=open('sno.projwfc_up')
nline_io_header=14 # line number at '    F    F'

for i in range(nline_io_header):
    filproj.readline()

for i in range(nlorb):
    filproj.readline()
    for j in range(nks):
        for k in range(nbnd):
            pjsum[i,j,k]=float(filproj.readline().split()[2])

nplotline=np.sum(iorb)
# oo, orbital index for each kind of color, oo can be generated by the following commands
#grep '[a-zA-Z]' sno.projwfc_up |grep 2P|awk '{printf( $1-1",")}'
oo=[[0,9],
[1,2,3,10,11,12],
[4,5,6,7,8,13,14,15,16,17],
[18,22,26,30],
[19,20,21,23,24,25,27,28,29,31,32,33]]

s_of_o=np.zeros([nks,],dtype=np.float32)
color=['r','g','orange','cyan','blue']
label=['Sn s','Sn p','Sn d','O s','O p']

scale=90.0
st=[]
for i in range(len(oo)):
    for k in range(nbnd):
        s_of_o=np.zeros([nks,],dtype=np.float32)
        for j in oo[i]:
            s_of_o[:]+=pjsum[j,:,k]
        if k == 0:
            st.append(plt.scatter(-1, ymin-1, 20, c=color[i], alpha=0.5, label=label[i],marker='.',edgecolor='none'))
        st.append(plt.scatter(np.arange(0,nks), eig[:,k]-eig_vbm, s=scale*s_of_o, c=color[i], alpha=0.5, marker='.',edgecolor='none'))

plt.xlim([0,nks-1]) # 201 points
plt.ylim([ymin,ymax])
plt.ylabel(r'E (eV)',fontsize=16)
plt.xticks((0,30,60), ('M', r'${\Gamma}$', 'Z') )

plt.subplots_adjust(left=0.20, right=0.75, top=0.95, bottom=0.1)
p1.legend(scatterpoints =1, numpoints=1,markerscale=2.0, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)

plt.savefig('pjband.png',dpi=400)

得到投影能带（Pjbands）如下：

pdos

脚本使用说明，见Link

(1) 首先在终端运行：grep ‘[a-zA-Z]’ sno.projwfc_up

输出了文件中所有包含字母的行，可以看出，Sn有5S，5P，4D轨道，O有2S，2P轨道，轨道的种类5类（当然还可以细分为P的3个小类Pz,Px,Py，和D轨道的5类，顺序见这里[2]），这里选择5类颜色作投影能带。

(2) 再运行grep '[a-zA-Z]' sno.projwfc_up |grep 'Sn 5S'|awk '{printf( $1-1",")}END{print("");}'

得到Sn 5S是哪些序号的，依次找到Sn的5S，5P，4D轨道，O的2S，2P轨道共5种，所以oo设置成：

1
2
3
4
5

oo=[[0,9],
[1,2,3,10,11,12],
[4,5,6,7,8,13,14,15,16,17],
[18,22,26,30],
[19,20,21,23,24,25,27,28,29,31,32,33]]

(3) 导体根据scf输出费米能级

eig_vbm = 0.0

绝缘体设置nband为最高价带，以下两行保留一行注释另一行

1
2

eig_vbm=max(eig[:, nband-1])
#eig_vbm= 0.00 # fermi energy level in scf output for metals

ymin,ymax，画图中y坐标范围

(4) 设置元素名，顺序不限

1

elem=['Sn','O']

每种元素的原子个数，顺序与下面orb保持一致

ielem=np.array([32,1],dtype=np.int32) # number of atoms for each element

每种原子的轨道类型，限s、p、d、f，每个元素一个list，要确保赝势含有对应轨道。

orb=[['s','p'],['s','p','s','d']] # projectors for each element

(5) 指定文件名，projwfc.x中由filproj设置的文件

filproj=open('sno.projwfc_up')

(6) 指定数据开始位置

运行grep -n ' F F' sno.projwfc_up

得到行号是41。设置如下：

1

nline_io_header=41

(7) 投影能带的颜色

1
2

color=['r','g','orange','cyan','blue']
label=['Sn s','Sn p','Sn d','O s','O p']

，更多颜色名称见(https://matplotlib.org/3.1.0/gallery/color/named_colors.html)。相应轨道名称，在画图图例中显示。

(8) 投影能带点的大小

1

scale=30.0

(9) k点路径名称

1

plt.xticks(vline, ('M'，r'${\Gamma}$', 'Z') )

k点横坐标

1

vline=(0,30,60)

(10) 输出图片文件名，分辨率

1

plt.savefig('pjband.png',dpi=400)

使用自写fortran脚本进行fatband画图

脚本下载1(需要能上GitHub)：https://github.com/xh125/QE-changecode/tree/main/bandfat

脚本下载2（Onedrive共享）：bandfat.f90 bandfat.in

脚本使用说明： 依次运行pw.x、projwfc.x，在能带计算之后，用projwfc.x生成的fatband.projwfs_up文件，可以和bands.x生成的文件carbyne.bands.dat和carbyne.bands.dat.gnu结合画出各个能带的原子轨道投影.将fortran代码编译生成可执行文件bandfat.x.输入文件使用bandfat.in,其输入文件如下：

1
2
3
4
5

&fatinput
    filband  = "Si.bands.dat"
    filproj  = "fatband"
    bandfatfile  = "fatband.dat.gnu"
/

其中输入参数需要满足 filband = Si.bands.dat filband need to set as the input of bands.x

filproj = fatband filproj need to set as the input of projwfc.x

bandfatfile = fatband.dat.gnu the fatband are write to file: fatband.dat.gnu

能带投影文件名为：

1
2
3
4
5
6
7
8
9
10
11
12
13

  DO is=1,nspin
     IF (nspin==2) THEN
        IF (is==1) filename=trim(filproj)//'.projwfc_up'
        IF (is==2) filename=trim(filproj)//'.projwfc_down'
        nksinit=(nkstot/2)*(is-1)+1
        nkslast=(nkstot/2)*is
        nk_ = nkstot/2
     ELSE
        filename=trim(filproj)//'.projwfc_up'
        nksinit=1
        nkslast=nkstot
        nk_ = nkstot
     ENDIF

投影能带中的内容：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

   WRITE (iunplot, '(a)') title
   WRITE (iunplot, '(8i8)') nr1x, nr2x, nr3x, nr1, nr2, nr3, nat, ntyp
   WRITE (iunplot, '(i6,6f12.8)') ibrav, celldm
   IF  (ibrav == 0) THEN
       WRITE ( iunplot, * ) at(:,1)
       WRITE ( iunplot, * ) at(:,2)
       WRITE ( iunplot, * ) at(:,3)
   ENDIF
   WRITE (iunplot, '(3f20.10,i6)') gcutm, dual, ecutwfc, 9
   WRITE (iunplot, '(i4,3x,a2,3x,f5.2)') &
         (nt, atm (nt), zv (nt), nt=1, ntyp)
   WRITE (iunplot, '(i4,3x,3f15.9,3x,i2)') (na, &
         (tau (i, na), i = 1, 3), ityp (na), na = 1, nat)
   WRITE (iunplot, '(3i8)') natomwfc, nkstot, nbnd
   WRITE (iunplot, '(2l5)') noncolin, lspinorb

     DO nwfc = 1, natomwfc
        IF (lspinorb) THEN
           WRITE(iunproj,'(2i5,1x,a4,1x,a2,1x,2i5,f5.1,1x,f5.1)') &
                nwfc, nlmchi(nwfc)%na, atm(ityp(nlmchi(nwfc)%na)), &
                nlmchi(nwfc)%els, nlmchi(nwfc)%n, nlmchi(nwfc)%l, &
                nlmchi(nwfc)%jj, &
                compute_mj(nlmchi(nwfc)%jj,nlmchi(nwfc)%l, nlmchi(nwfc)%m)
        ELSE IF (noncolin) THEN
           WRITE(iunproj,'(2i5,1x,a4,1x,a2,1x,3i5,1x,f4.1)') &
                nwfc, nlmchi(nwfc)%na, atm(ityp(nlmchi(nwfc)%na)), &
                nlmchi(nwfc)%els, nlmchi(nwfc)%n, nlmchi(nwfc)%l, &
                nlmchi(nwfc)%m, &
                0.5d0-int(nlmchi(nwfc)%ind/(2*nlmchi(nwfc)%l+2))
        ELSE
           WRITE(iunproj,'(2i5,1x,a4,1x,a2,1x,3i5)') &
                nwfc, nlmchi(nwfc)%na, atm(ityp(nlmchi(nwfc)%na)), &
                nlmchi(nwfc)%els, nlmchi(nwfc)%n, nlmchi(nwfc)%l,  &
                nlmchi(nwfc)%m
        END IF
        DO ik=nksinit,nkslast
           DO ibnd=1,nbnd
              WRITE(iunproj,'(2i8,f20.10)') ik,ibnd, abs(proj(nwfc,ibnd,ik))
           ENDDO
        ENDDO
     ENDDO

首先运行：grep "[a-zA-Z]" fatband.projwfc_up，得到的每一个投影轨道的信息：

fatband-pwf

具体投影轨道的波函数可以参考Wannier90 User Guide中53页的projector的相关波函数。分别对应于轨道n,l,m好量子数。投影轨道，与晶格的取向有关。 angular function radialfunc

输出了文件中所有包含字母的行，可以看出，C有2S，2P 轨道，轨道的种类4类，这里选择4类颜色作投影能带。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228

module kinds
  implicit none
  integer,parameter :: dp= kind(1.0d0)
end module 

module constants
  implicit none
  integer,parameter :: maxlen = 86
end module 

module io
  use kinds
  use constants
  implicit none
  integer,public,save :: stdout,stdin
  character(len=maxlen) :: stdout_name="bandfat.out"
  character(len=maxlen) :: stdin_name = "bandfat.in"
  character(len=maxlen) :: msg
  integer :: ierr
  contains
  
  function io_file_unit() !得到一个当前未使用的unit，用于打开文件
  !===========================================                                     
  !! Returns an unused unit number
  !! so we can later open a file on that unit.                                       
  !===========================================
  implicit none

    integer :: io_file_unit,unit_index
    logical :: file_open

    unit_index = 9
    file_open  = .true.
    do while ( file_open )
      unit_index = unit_index + 1
      inquire( unit=unit_index, OPENED = file_open ) !用于检查文件状态，参考P.536
    end do
    
    io_file_unit = unit_index

    return
  
  end function io_file_unit
  
  subroutine open_file(file_name,file_unit)
    implicit none
    
    character(len=*),intent(in) :: file_name
    integer,intent(in)          :: file_unit
    open(unit=file_unit, file=file_name,iostat=ierr,iomsg=msg)
    if(ierr /= 0 ) then
      call io_error('Error: Problem opening "'//trim(adjustl(file_name))//' " file')
      call io_error(msg)
    endif
  end subroutine open_file
  
  subroutine close_file(file_name,file_unit)
    implicit none
    
    integer,intent(in)          :: file_unit
    character(len=*),intent(in) :: file_name
    close(file_unit,iostat=ierr,iomsg=msg)
    if(ierr /= 0 ) then
      call io_error('Error: Problem close "'//trim(adjustl(file_name))//' " file')
      call io_error(msg)
    endif   
  end subroutine close_file  
  
  !========================================
  subroutine io_error ( error_msg )
  !========================================
  !! Abort the code giving an error message 
  !========================================

    implicit none
    character(len=*), intent(in) :: error_msg

    write(stdout,*)  'Exiting.......' 
    write(stdout, '(1x,a)') trim(error_msg)    
    close(stdout)    
    write(*, '(1x,a)') trim(error_msg)
    write(*,'(A)') "Error: examine the output/error file for details" 
    STOP
         
  end subroutine io_error
  
  subroutine findkline(funit,kline,indexi,indexe)
    implicit none
      integer,intent(in):: funit
      character(len=*),intent(in) :: kline
      integer,intent(in)::indexi,indexe
      logical :: lfindkline
      character(len=maxlen) :: ctmp
      
      lfindkline = .FAlSE.
      do while(.NOT. lfindkline)
        read(funit,"(A)") ctmp
        if (ctmp(indexi:indexe)==kline) then
          lfindkline = .TRUE.
          backspace(unit=funit)
        endif
      enddo  
  end subroutine findkline  
  
end module 

module parameters
  use kinds
  use constants
  use io
  implicit none
  character(len=maxlen) :: filband,filband0
  character(len=maxlen) :: filproj
  character(len=maxlen) :: bandfatfile
  integer :: bandsunit,projsunit,bandfatunit,bands0unit
  
  namelist /fatinput/ filband,filproj,bandfatfile
  
  contains
  
  subroutine readinput()
    implicit none
    
    stdin = io_file_unit()
    call open_file(stdin_name,stdin)
    
    filband="bands.dat"
    filproj="fatband"
    bandfatfile="fatband.dat.gnu"

    read(UNIT=stdin,nml=fatinput,iostat=ierr,iomsg=msg)
    if(ierr /= 0) then
      call io_error('Error: Problem reading namelist file fatinput')
      call io_error(msg)
    endif  
    close(stdin)    
    
    write(stdout,"(1X,A13,1X,A)") "filband     =", trim(filband)
    write(stdout,*) "filband need to set as the input of bands.x  "
    write(stdout,"(/,1X,A13,1X,A)") "filproj     =", trim(filproj)
    write(stdout,*) "filproj need to set as the input of projwfc.x "
    write(stdout,"(/,1X,A13,1X,A)") "bandfatfile =", trim(bandfatfile)
    write(stdout,*) "the fatband are write to file: ",trim(bandfatfile)
    
  end subroutine readinput
  
end module parameters  

program main
  use kinds
  use constants
  use parameters
  use io
  implicit none
  
  integer :: nbnd,nks,natomwfc
  integer :: ik,iband,ipol,iwfc
  character(len=maxlen) :: ctmp,cformat
  real(kind=dp),allocatable :: kpoints(:,:),Enk(:,:),proj(:,:,:)
  character(len=13),allocatable :: nameatomwfc(:)
  real(kind=dp),allocatable   :: lkline(:)
  
  stdout = io_file_unit()
  call open_file(stdout_name,stdout)
  call readinput()
  bandfatunit = io_file_unit()
  call open_file(bandfatfile,bandfatunit)
  
  
  bandsunit = io_file_unit()
  call open_file(filband,bandsunit)
  read(bandsunit,"(12X,I4,6X,I6)") nbnd,nks
  allocate(kpoints(3,nks),Enk(nbnd,nks))
  write(ctmp,*) nbnd
  do ik=1,nks
    read(bandsunit,'(10x,3f10.6)') (kpoints(ipol,ik),ipol=1,3)
    read(bandsunit,"(10f9.3)") (Enk(iband,ik),iband=1,nbnd)
  enddo
  
  allocate(lkline(nks))
  lkline = 0.0
  filband0 = trim(filband)//".gnu"
  bands0unit = io_file_unit()
  call open_file(filband0,bands0unit)
  do ik =1,nks
    read(bands0unit,"(f10.4)") lkline(ik)
  enddo
  call close_file(filband0,bands0unit)
  
  
  
  
  !WRITE (iunplot, '(3i8)') natomwfc, nkstot, nbnd
  !WRITE (iunplot, '(2l5)') noncolin, lspinorb
  projsunit = io_file_unit()
  filproj = trim(filproj)//".projwfc_up"
  call open_file(filproj,projsunit)
  call findkline(projsunit,"    F    F",1,10)
  backspace(projsunit)
  read(projsunit,"(3i8)") natomwfc,nks,nbnd
  allocate(proj(natomwfc,nbnd,nks))
  allocate(nameatomwfc(natomwfc))
  read(projsunit,*)

  
  do iwfc=1,natomwfc
    read(projsunit,"(5X,A13)") nameatomwfc(iwfc)
    do ik=1,nks
      do iband=1,nbnd
        read(projsunit,"(16X,F20.10)") proj(iwfc,iband,ik)
      enddo
    enddo
  enddo
  
  
  
  write(ctmp,*) natomwfc
  cformat = "(A10,1X,A10,"//trim(adjustl(ctmp))//"(1X,A13))"
  write(bandfatunit,cformat) "longkpoint","Energy(eV)",(nameatomwfc(iwfc),iwfc=1,natomwfc)
  cformat = "(F10.5,1X,F10.5,"//trim(adjustl(ctmp))//"(1X,F13.5))"
  do iband=1,nbnd
    do ik=1,nks
      write(bandfatunit,cformat) lkline(ik),Enk(iband,ik),(proj(iwfc,iband,ik),iwfc=1,natomwfc) 
    enddo
    write(bandfatunit,*)
  enddo

end program

采用Originlab作图得到的fatband的图如下所示：其中红色为pz轨道，蓝色为sp2轨道。使用命令：grep high bands.out得到高对称点在线性模式下的坐标：

1
2
3
4

     high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   0.0000
     high-symmetry point:  0.5000 0.2887 0.0000   x coordinate   0.5774
     high-symmetry point:  0.6667 0.0000 0.0000   x coordinate   0.9107
     high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   1.5774

fatband

CATALOG
    FEATURED TAGS
    LVCSH QUANTUM ESPRESSO Wannier Fortran QE 激活秘钥 DFPT 安装 EPW Endnote X20 Visual Studio WSL phonon
    FRIENDS