Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_tetraef
8 : ! -----------------------------------------------------------------------
9 : ! This subroutine evaluates the density of states with the tetrahedron method
10 : ! and sets the weight factors needed for the charge density for bulk systems.
11 : ! Adapted for the FLEUR code GB 2000
12 : ! -----------------------------------------------------------------------
13 : USE m_types
14 : USE m_constants
15 : USE m_juDFT
16 :
17 : IMPLICIT NONE
18 :
19 : CONTAINS
20 :
21 2 : SUBROUTINE tetra_ef(kpts,jspins,lb,ub,eig,zc,xfac,efermi,w,l_output)
22 :
23 : TYPE(t_kpts), INTENT(IN) :: kpts
24 : INTEGER, INTENT(IN) :: jspins
25 : REAL, INTENT(IN) :: lb,ub,zc,xfac
26 : REAL, INTENT(INOUT) :: eig(:,:,:) !(neigd,nkptd,jspd)
27 : REAL, INTENT(OUT) :: w(:,:,:) !(neigd,nkptd,jspd)
28 : REAL, INTENT(OUT) :: efermi
29 : LOGICAL, INTENT(IN) :: l_output
30 :
31 : INTEGER :: i,j,jspin,iBand,ikpt,nelec,ncr,itet,it,icorn,jcorn
32 : REAL :: elow,dlow,eup,dup,ttt,dfermi,wgs
33 2 : REAL :: weight(4),ecmax(2,size(w,1))
34 2 : REAL :: wght(2,kpts%nkpt,size(w,1)),eval(4), tetra_weight(kpts%nkpt)
35 :
36 764 : w=0.0
37 2198 : wght=0.0
38 38 : DO iBand = 1,size(w,1)
39 74 : DO jspin = 1,jspins
40 36 : ecmax(jspin,iBand) = -1.0e25
41 792 : DO ikpt = 1,kpts%nkpt
42 720 : w(iBand,ikpt,jspin) = 0.0
43 756 : IF(eig(iBand,ikpt,jspin).GT.ecmax(jspin,iBand)) ecmax(jspin,iBand) = eig(iBand,ikpt,jspin)
44 : ENDDO
45 : ENDDO
46 : ENDDO
47 : !
48 : ! check for energy degeneracies in tetrahedrons
49 : !
50 4 : DO jspin = 1,jspins
51 92 : DO itet = 1,kpts%ntet
52 1674 : DO iBand = 1,size(w,1)
53 6424 : DO i = 1,3
54 4752 : icorn = kpts%ntetra(i,itet)
55 15840 : DO j = i+1,4
56 9504 : jcorn = kpts%ntetra(j,itet)
57 14256 : IF(abs(eig(iBand,icorn,jspin)-eig(iBand,jcorn,jspin)).LT.1.0e-7) THEN
58 16 : eig(iBand,icorn,jspin) = eig(iBand,icorn,jspin) + 1.0e-7
59 16 : eig(iBand,jcorn,jspin) = eig(iBand,jcorn,jspin) - 1.0e-7
60 : ENDIF
61 : ENDDO
62 : ENDDO
63 : ENDDO
64 : ENDDO
65 : ENDDO
66 : !
67 : ! calculate weight factors
68 : !
69 42 : tetra_weight = 0.0
70 90 : DO itet=1,kpts%ntet
71 440 : DO i=1, 4
72 440 : tetra_weight(kpts%ntetra(i,itet)) = tetra_weight(kpts%ntetra(i,itet)) + kpts%voltet(itet)/(4.0*kpts%ntet)
73 : ENDDO
74 1674 : DO iBand=1,size(w,1)
75 3256 : DO jspin=1,jspins
76 :
77 15840 : eval = eig(iBand,kpts%ntetra(:,itet),jspin)
78 :
79 7920 : IF(ANY(eval.GE.9999.9)) CYCLE
80 :
81 7920 : DO i=1,4
82 6336 : weight(i) = 1.0
83 33264 : DO j=1,4
84 31680 : IF(i.NE.j) weight(i) = weight(i)*(eval(j)-eval(i))
85 : ENDDO
86 : ENDDO
87 9504 : DO i=1,4
88 6336 : icorn = kpts%ntetra(i,itet)
89 6336 : weight(i) = 6.0*kpts%voltet(itet)/(weight(i)*kpts%ntet)
90 7920 : wght(jspin,icorn,iBand) = wght(jspin,icorn,iBand) + weight(i)
91 : ENDDO
92 :
93 : ENDDO
94 : ENDDO
95 : ENDDO
96 : !
97 : !xfac = 2.0/jspins
98 42 : tetra_weight = xfac*tetra_weight
99 2198 : wght = xfac*wght
100 : !
101 : !---------------------------------------------------
102 : ! determine fermi energy
103 : !---------------------------------------------------
104 : !
105 2 : nelec = zc ! determine # of electrons
106 2 : ncr = 0 ! (modified gb2000)
107 :
108 2 : elow = lb ! determine lower bound
109 2 : dlow = ncr
110 42 : DO ikpt = 1,kpts%nkpt
111 762 : DO iBand = 1,size(w,1)
112 1480 : DO jspin = 1,jspins
113 720 : ttt = elow - eig(iBand,ikpt,jspin)
114 720 : IF ( elow.GT.ecmax(jspin,iBand) ) THEN
115 0 : dlow = dlow + tetra_weight(ikpt)
116 0 : cycle
117 : endif
118 720 : IF (ttt.LT.0.0e0) cycle
119 1440 : dlow = dlow + wght(jspin,ikpt,iBand)*ttt*ttt*ttt/6
120 : ENDDO
121 : ENDDO
122 : ENDDO
123 2 : IF (dlow.GT.nelec) THEN
124 0 : WRITE (oUnit,180) elow,dlow,nelec
125 : 180 FORMAT (' valence band too high ',/,' elow ',f10.5,' dlow ',f10.5,' nelec ',i5)
126 0 : CALL juDFT_error("dos: valence band too high ",calledby="tetra_ef")
127 : ENDIF
128 :
129 :
130 2 : it = 0
131 2 : eup = ub
132 2 : dup = 0.0 ! determine upper bound
133 4 : DO WHILE ((dup-nelec).LT.0.00001)
134 2 : dup = ncr
135 42 : DO ikpt = 1,kpts%nkpt
136 762 : DO iBand = 1,size(w,1)
137 1480 : DO jspin = 1,jspins
138 720 : ttt = eup - eig(iBand,ikpt,jspin)
139 720 : IF ( eup.GT.ecmax(jspin,iBand) ) THEN
140 400 : dup = dup + tetra_weight(ikpt)
141 400 : cycle
142 : endif
143 320 : IF (ttt.LT.0.0e0) cycle
144 1040 : dup = dup + wght(jspin,ikpt,iBand)*ttt*ttt*ttt/6
145 : ENDDO
146 : ENDDO
147 : ENDDO
148 :
149 4 : IF((dup-nelec).LT.0.00001) THEN
150 0 : eup = eup + 0.2
151 0 : it = it + 1
152 0 : IF(it .gt. 10) THEN
153 0 : WRITE (oUnit,200) eup,dup,nelec
154 : 200 FORMAT (' valence band too low ',/,' eup ',f10.5,' dup ',f10.5,' nelec ',i5)
155 0 : CALL juDFT_error("dos: valence band too low ",calledby ="tetra_ef")
156 : END IF
157 : ENDIF
158 : ENDDO
159 :
160 :
161 :
162 68 : DO WHILE ((eup-elow).GT.1.0e-10) ! iterate for fermi-energy
163 66 : efermi = 0.5*(elow+eup)
164 66 : dfermi = real(ncr)
165 1386 : DO ikpt = 1,kpts%nkpt
166 25146 : DO iBand = 1,size(w,1)
167 48840 : DO jspin = 1,jspins
168 23760 : ttt =efermi-eig(iBand,ikpt,jspin)
169 23760 : IF ( efermi.GT.ecmax(jspin,iBand) ) THEN
170 12640 : dfermi = dfermi + tetra_weight(ikpt)
171 12640 : cycle
172 : endif
173 11120 : IF (ttt.LT.0.0e0) cycle
174 34880 : dfermi = dfermi + wght(jspin,ikpt,iBand)*ttt*ttt*ttt/6
175 : ENDDO
176 : ENDDO
177 : ENDDO
178 68 : IF (dfermi.GT.nelec) THEN
179 36 : eup = efermi
180 : ELSE
181 30 : elow = efermi
182 : ENDIF
183 : ENDDO
184 2 : if (l_output) then
185 2 : WRITE (oUnit,220) efermi,dfermi,nelec
186 : 220 FORMAT (//,'>>> D O S <<<',//,' fermi energy =',f10.5,' dtot =',f10.5,' nelec =',i5)
187 : endif
188 : !
189 : !---------------------------------------------------
190 : ! calculate weight factors for charge density integration
191 : !---------------------------------------------------
192 : !
193 90 : DO itet = 1,kpts%ntet
194 1674 : DO iBand = 1,size(w,1)
195 3256 : DO jspin = 1,jspins
196 15840 : eval = eig(iBand,kpts%ntetra(:,itet),jspin)
197 :
198 7920 : IF(ANY(eval.GE.9999.9)) CYCLE
199 :
200 7920 : DO i = 1,4
201 6336 : weight(i) = 1.0
202 31680 : DO j = 1,4
203 31680 : IF(i.NE.j) THEN
204 19008 : weight(i) = weight(i) * (eval(j) - eval(i))
205 : ENDIF
206 : ENDDO
207 7920 : weight(i) = 6.0 * kpts%voltet(itet)/(weight(i)*kpts%ntet)
208 : ENDDO
209 :
210 1584 : IF ( efermi.GT.ecmax(jspin,iBand) ) THEN
211 880 : wgs = kpts%voltet(itet)/(4.0* kpts%ntet)
212 : else
213 : wgs = 0.0e0
214 3520 : DO i = 1,4
215 2816 : ttt = efermi - eval(i)
216 2816 : IF( ttt.LT.0.0e0 ) cycle
217 3520 : wgs = wgs + ttt**3*weight(i)
218 : ENDDO
219 704 : wgs = wgs / 24.0
220 : endif
221 :
222 30096 : w(iBand,kpts%ntetra(:,itet),jspin) = w(iBand,kpts%ntetra(:,itet),jspin) + wgs
223 :
224 : ENDDO
225 : ENDDO
226 : ENDDO
227 :
228 : ! DO jspin = 1,jspins
229 : ! DO ikpt = 1,kpts%nkpt
230 : ! DO iBand = 1,size(w,1)
231 : ! w(iBand,ikpt,jspin) = xfac * w(iBand,ikpt,jspin)
232 : ! ENDDO
233 : ! ENDDO
234 : ! ENDDO
235 :
236 2 : END SUBROUTINE tetra_ef
237 : END MODULE m_tetraef
|
