Line data Source code
1 : MODULE m_tetrahedronInit
2 :
3 : !------------------------------------------------------------------------------------
4 : !This module provides the weights for the linear tetrahedron/triangular method
5 : !for Brillouin-zone integration with the step function as the weight
6 : !
7 : !This module should be used by calling the interface tetrahedronInit
8 : !
9 : !Structure:
10 : ! tetrahedronInit:
11 : ! getWeightKpoints: gives weight at one (fermi) energy for all kpoints
12 : ! getWeightEnergyMesh: gives weight on an energy grid for one kpoint.
13 : ! By providing the dos switch as .TRUE. the weights are
14 : ! numerically differentiated for use for DOS calculations
15 : !
16 : !These subroutines all call getWeightSingleBand which handles the actual calculation
17 : !of weights. These differ in the order of the loops and the arguments provided.
18 : !------------------------------------------------------------------------------------
19 :
20 : USE m_types_kpts
21 : USE m_types_input
22 : USE m_juDFT
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 : PUBLIC :: tetrahedronInit
28 :
29 : INTERFACE tetrahedronInit
30 : PROCEDURE getWeightKpoints, getWeightEnergyMesh
31 : END INTERFACE
32 :
33 : REAL, PARAMETER :: weightCutoff=1e-14
34 : real, parameter :: tol_degeneracy=1e-9
35 :
36 : CONTAINS
37 :
38 1794 : SUBROUTINE getWeightKpoints(kpts,input,eig,neig,efermi,weightSum,weights)
39 :
40 : TYPE(t_kpts), INTENT(IN) :: kpts
41 : TYPE(t_input), INTENT(IN) :: input
42 : REAL, INTENT(IN) :: eig(:,:)
43 : INTEGER, INTENT(IN) :: neig
44 : REAL, INTENT(IN) :: efermi
45 : REAL, OPTIONAL,INTENT(INOUT) :: weightSum
46 : REAL, OPTIONAL,INTENT(INOUT) :: weights(:,:)
47 :
48 : INTEGER :: ikpt,ncorn,itet,icorn,iband,num_degenerate_states,last_band
49 1794 : REAL :: w(1),etetra(SIZE(kpts%ntetra,1)),vol(kpts%ntet)
50 : logical :: l_weights_pres, l_weightsum_pres
51 :
52 1794 : IF(.NOT.PRESENT(weightSum).AND..NOT.PRESENT(weights)) THEN
53 : CALL juDFT_error("No output variable provided (either weightSum or weights)",&
54 0 : calledby="getWeightKpoints")
55 : ENDIF
56 :
57 1794 : l_weights_pres = PRESENT(weights)
58 1794 : l_weightsum_pres = PRESENT(weightSum)
59 : !Tetrahedra or Triangles?
60 1794 : ncorn = MERGE(3,4,input%film)
61 :
62 37604 : IF(PRESENT(weights)) weights = 0.0
63 1794 : IF(PRESENT(weightSum)) weightSum = 0.0
64 :
65 503334 : vol = kpts%voltet(:)/kpts%ntet
66 :
67 : !More efficient to just loop through all tetrahedra
68 : !$OMP parallel do default(none) &
69 : !$OMP shared(neig,ncorn,vol,l_weights_pres,l_weightsum_pres) &
70 : !$OMP shared(kpts,input,eig,weights,efermi) &
71 : !$OMP private(itet,ikpt,icorn,iband,etetra,w) &
72 1794 : !$OMP reduction(+:weightSum) collapse(3)
73 : DO itet = 1, kpts%ntet
74 : DO icorn = 1, ncorn
75 : DO iband = 1, neig
76 :
77 : ikpt = kpts%ntetra(icorn,itet)
78 : etetra = eig(iband,kpts%ntetra(:,itet))
79 :
80 : IF( ALL(etetra>efermi + 1E-8) .AND. .NOT.input%l_bloechl ) CYCLE
81 :
82 : w = getWeightSingleBand([efermi],etetra,icorn,vol(itet),input%film,input%l_bloechl)
83 :
84 : IF(l_weightsum_pres) weightSum = weightSum + w(1)
85 :
86 : IF(l_weights_pres) THEN
87 : !$OMP critical
88 : weights(iband,ikpt) = weights(iband,ikpt) + w(1)
89 : !$OMP end critical
90 : ENDIF
91 : ENDDO
92 : ENDDO
93 : ENDDO
94 : !$OMP end parallel do
95 :
96 1794 : if (l_weights_pres) then
97 : !find degenerate states
98 2110 : do ikpt=1,kpts%nkpt
99 : iband = 1
100 31258 : do while(iband<neig)
101 : num_degenerate_states=1
102 30116 : do while (abs(eig(iband,ikpt)-eig(iband+num_degenerate_states,ikpt))<tol_degeneracy)
103 986 : num_degenerate_states=num_degenerate_states+1
104 30116 : if (iband+num_degenerate_states>neig) exit
105 : enddo
106 :
107 29148 : if (num_degenerate_states>1) THEN
108 850 : last_band=iband+num_degenerate_states-1
109 : !Make sure all weights are equal
110 4522 : weights(iband:last_band,ikpt)=sum(weights(iband:last_band,ikpt))/num_degenerate_states
111 : iband=last_band
112 : endif
113 29148 : iband=iband+1
114 : enddo
115 : enddo
116 : endif
117 :
118 1794 : END SUBROUTINE getWeightKpoints
119 :
120 144 : SUBROUTINE getWeightEnergyMesh(kpts,input,ikpt,eig,neig,eMesh,weights,resWeights,bounds,dos)
121 :
122 : USE m_differentiate
123 :
124 : TYPE(t_kpts), INTENT(IN) :: kpts
125 : TYPE(t_input), INTENT(IN) :: input
126 : REAL, INTENT(IN) :: eig(:,:)
127 : REAL, INTENT(INOUT) :: weights(:,:)
128 : INTEGER,OPTIONAL, INTENT(INOUT) :: bounds(:,:)
129 : REAL, OPTIONAL, INTENT(INOUT) :: resWeights(:,:)
130 :
131 : INTEGER, INTENT(IN) :: neig
132 : INTEGER, INTENT(IN) :: ikpt
133 : REAL, INTENT(IN) :: eMesh(:)
134 : LOGICAL,OPTIONAL, INTENT(IN) :: dos
135 :
136 : INTEGER :: itet,iband,ncorn,ie,icorn,ntet,num_degenerate_states, last_band,i
137 : LOGICAL :: l_dos, l_bounds_pres, l_resWeights_pres
138 144 : REAL :: etetra(SIZE(kpts%ntetra,1),neig,SIZE(kpts%tetraList,1)),del,vol
139 144 : REAL, ALLOCATABLE :: dos_weights(:)
140 : !Temporary Arrays to include end points
141 : !to avoid numerical trouble with differentiation
142 144 : REAL, ALLOCATABLE :: calc_weights(:,:)
143 144 : REAL, ALLOCATABLE :: calc_weights_thread(:,:),resWeights_thread(:,:)
144 144 : REAL, ALLOCATABLE :: calc_eMesh(:)
145 144 : real, allocatable :: mean_weight(:)
146 :
147 : !Extract the right eigenvalues
148 144 : IF(kpts%tetraList(1,ikpt)==0) CALL juDFT_error("No tetrahedrons in tetraList", calledby="tetrahedronInit")
149 : ntet = 1
150 : do
151 3936 : itet = kpts%tetraList(ntet,ikpt)
152 :
153 178486 : do iband = 1, neig
154 1723036 : etetra(:,iband,ntet) = eig(iband,kpts%ntetra(:,itet))
155 : enddo
156 4080 : if(ntet < SIZE(kpts%tetraList,1)) THEN
157 3936 : if(kpts%tetraList(ntet+1,ikpt)>0) THEN
158 : ntet = ntet + 1
159 : else
160 : exit
161 : endif
162 : else
163 : exit
164 : endif
165 : enddo
166 :
167 144 : l_dos = PRESENT(dos)
168 144 : IF(PRESENT(dos))THEN
169 144 : l_dos = dos.AND.SIZE(eMesh)>1
170 : ENDIF
171 :
172 144 : IF(PRESENT(resWeights)) resWeights = 0.0
173 :
174 144 : IF(l_dos) THEN
175 34540152 : ALLOCATE(calc_weights(SIZE(eMesh)+2,neig),source=0.0)
176 713520 : ALLOCATE(calc_eMesh(SIZE(eMesh)+2),source=0.0)
177 713520 : ALLOCATE(mean_weight(SIZE(eMesh)+2),source=0.0)
178 144 : del = eMesh(2)-eMesh(1)
179 1426320 : calc_eMesh = [eMesh(1)-del,eMesh(:),eMesh(SIZE(eMesh))+del]
180 : ELSE
181 0 : ALLOCATE(calc_weights(SIZE(eMesh),neig),source=0.0)
182 0 : ALLOCATE(calc_eMesh(SIZE(eMesh)),source=0.0)
183 0 : ALLOCATE(mean_weight(SIZE(eMesh)),source=0.0)
184 0 : calc_eMesh = eMesh
185 : ENDIF
186 :
187 144 : l_resWeights_pres = PRESENT(resWeights)
188 : !$OMP parallel default(none)&
189 : !$OMP shared(neig,ntet,l_resWeights_pres,kpts,input,ikpt,resWeights,calc_weights) &
190 : !$OMP shared(etetra,calc_eMesh,eMesh) &
191 144 : !$OMP private(itet,iband,vol,icorn,resWeights_thread,calc_weights_thread)
192 : IF(l_resWeights_pres) ALLOCATE(resWeights_thread(SIZE(resWeights,1),SIZE(resWeights,2)),source = 0.0)
193 : ALLOCATE(calc_weights_thread(SIZE(calc_weights,1),SIZE(calc_weights,2)),source = 0.0)
194 : !$OMP do collapse(3)
195 : DO itet = 1, ntet
196 : DO iband = 1, neig
197 : DO icorn = 1, SIZE(kpts%ntetra,1)
198 :
199 : vol = kpts%voltet(kpts%tetraList(itet,ikpt))/kpts%ntet
200 : IF(kpts%ntetra(icorn,kpts%tetraList(itet,ikpt)).NE.ikpt) CYCLE
201 :
202 : IF(l_resWeights_pres) THEN
203 : resWeights_thread(:,iband) = resWeights_thread(:,iband) &
204 : + getWeightSingleBand(eMesh,etetra(:,iband,itet),icorn,vol,&
205 : input%film,input%l_bloechl,l_res=.TRUE.)
206 : ENDIF
207 :
208 : IF( ALL(etetra(:,iband,itet)>MAXVAL(calc_eMesh) + 1E-8) .AND. .NOT.input%l_bloechl ) CYCLE
209 :
210 : calc_weights_thread(:,iband) = calc_weights_thread(:,iband) &
211 : + getWeightSingleBand(calc_eMesh,etetra(:,iband,itet),icorn,vol,&
212 : input%film,input%l_bloechl)
213 : ENDDO
214 : ENDDO
215 : ENDDO
216 : !$OMP end do
217 : !$OMP critical
218 : IF(l_resWeights_pres) resWeights = resWeights + resWeights_thread
219 : calc_weights = calc_weights + calc_weights_thread
220 : !$OMP end critical
221 : DEALLOCATE(calc_weights_thread)
222 : IF(l_resWeights_pres) DEALLOCATE(resWeights_thread)
223 : !$OMP end parallel
224 :
225 : !find degenerate states
226 144 : iband = 1
227 6604 : do while(iband<neig)
228 : num_degenerate_states=1
229 7448 : do while (abs(eig(iband,ikpt)-eig(iband+num_degenerate_states,ikpt))<tol_degeneracy)
230 992 : num_degenerate_states=num_degenerate_states+1
231 7448 : if (iband+num_degenerate_states>neig) exit
232 : enddo
233 :
234 6460 : if (num_degenerate_states>1) THEN
235 854 : last_band=iband+num_degenerate_states-1
236 : !Make sure all weights are equal
237 13841908 : mean_weight = sum(calc_weights(:,iband:last_band), dim=2)/num_degenerate_states
238 2700 : do i = iband, last_band
239 9477992 : calc_weights(:,i) = mean_weight
240 : enddo
241 : iband=last_band
242 : endif
243 6604 : iband=iband+1
244 : enddo
245 :
246 34524536 : weights = 0.0
247 144 : l_bounds_pres = PRESENT(bounds)
248 : !-------------------------------------
249 : ! PostProcess weights
250 : !-------------------------------------
251 : !$OMP parallel default(none) &
252 : !$OMP shared(neig,l_dos,del,eMesh) &
253 : !$OMP shared(calc_weights,weights,bounds,l_bounds_pres,l_resWeights_pres) &
254 144 : !$OMP private(iband,dos_weights,ie)
255 : IF(l_dos) ALLOCATE(dos_weights(SIZE(eMesh)+2),source=0.0)
256 : !$OMP do schedule(dynamic,1)
257 : DO iband = 1, neig
258 : !---------------------------------------------------
259 : ! Weights for DOS -> differentiate with respect to E
260 : !---------------------------------------------------
261 : IF(l_dos) THEN
262 : CALL diff3(calc_weights(:,iband),del,dos_weights)
263 : weights(1:SIZE(eMesh),iband) = dos_weights(2:SIZE(eMesh)+1)
264 : ELSE
265 : weights(:,:neig) = calc_weights
266 : ENDIF
267 : !--------------------------------------------------------------
268 : ! Find the range where the weights are bigger than weightCutoff
269 : !--------------------------------------------------------------
270 : IF(l_bounds_pres.AND..NOT.l_resWeights_pres) THEN
271 : !--------------------
272 : ! Lower bound
273 : !--------------------
274 : ie = 1
275 : DO
276 : IF(ABS(weights(ie,iband)).GT.weightCutoff) THEN
277 : bounds(iband,1) = ie
278 : EXIT
279 : ELSE
280 : weights(ie,iband) = 0.0
281 : ie = ie + 1
282 : IF(ie.EQ.SIZE(eMesh)) THEN
283 : bounds(iband,1) = SIZE(eMesh)
284 : EXIT
285 : ENDIF
286 : ENDIF
287 : ENDDO
288 : !--------------------
289 : ! Upper bound
290 : !--------------------
291 : ie = SIZE(eMesh)
292 : DO
293 : IF(ABS(weights(ie,iband)).GT.weightCutoff) THEN
294 : bounds(iband,2) = ie
295 : EXIT
296 : ELSE
297 : weights(ie,iband) = 0.0
298 : ie = ie - 1
299 : IF(ie.EQ.0) THEN
300 : bounds(iband,2) = 1
301 : EXIT
302 : ENDIF
303 : ENDIF
304 : ENDDO
305 : !For safety
306 : IF(bounds(iband,1).GT.bounds(iband,2)) THEN
307 : bounds(iband,1) = 1
308 : bounds(iband,2) = 1
309 : ENDIF
310 : ELSE IF(l_bounds_pres) THEN
311 : bounds(iband,1) = 1
312 : bounds(iband,2) = SIZE(eMesh)
313 : ENDIF
314 : ENDDO
315 : !$OMP end do
316 : IF(l_dos) DEALLOCATE(dos_weights)
317 : !$OMP end parallel
318 :
319 34524536 : IF(ANY(weights(:,:neig)<0.0).AND..NOT.input%l_bloechl) THEN
320 0 : CALL juDFT_error("TetraWeight error: Unexpected negative weight", calledby="getWeightEnergyMesh")
321 : ENDIF
322 :
323 144 : END SUBROUTINE getWeightEnergyMesh
324 :
325 :
326 25068210 : PURE FUNCTION getWeightSingleBand(eMesh,etetra,icorn,vol,film,l_bloechl,l_res) Result(weight)
327 :
328 : !--------------------------------------------------------------
329 : ! This is the core routine calculating the integration
330 : ! weight for one kpoint in one tetrahedron for one band index
331 : ! for the provided energy points
332 : ! All other routines wrap this procedure so that it is most
333 : ! efficient in all cases
334 : !--------------------------------------------------------------
335 :
336 : USE m_tetsrt
337 : USE m_tetraWeight
338 : USE m_resWeight
339 : USE m_bloechl
340 :
341 : REAL, INTENT(IN) :: eMesh(:) !Energy points, where the weights are calculated
342 : REAL, INTENT(IN) :: etetra(:) !Eigenvalues at the corners of the tetrahedron
343 : INTEGER, INTENT(IN) :: icorn !Current k-point index (We calculate the weights for this corner)
344 : REAL, INTENT(IN) :: vol !Volume of the tetrahedron
345 : LOGICAL, INTENT(IN) :: film !Switch controls wether tetrahedron/triangular method is used
346 : LOGICAL, INTENT(IN) :: l_bloechl !Controls bloechl corrections
347 : LOGICAL, OPTIONAL,INTENT(IN) :: l_res
348 :
349 : REAL :: weight(SIZE(eMesh))
350 25068210 : INTEGER :: i,ie,nstart,nend,ind(SIZE(etetra)),icornSorted
351 : REAL :: eb,et,del
352 : LOGICAL :: l_calcres
353 :
354 :
355 : !Sort the eigenvalues at the corners (ascending order in ind)
356 25068210 : ind = tetsrt(etetra)
357 :
358 : !Find out where the corner went
359 125281488 : DO i = 1, SIZE(ind)
360 125281488 : IF(ind(i)==icorn) icornSorted = i
361 : ENDDO
362 :
363 25068210 : l_calcres = .FALSE.
364 25068210 : IF(PRESENT(l_res)) l_calcres = l_res
365 :
366 0 : IF(l_calcres) THEN
367 0 : weight = resWeight(eMesh,etetra(ind),icornSorted,vol,film)
368 : ELSE
369 : !Find the range in the (equidistant) energy mesh where the weights are changing
370 25068210 : IF( SIZE(eMesh)>1 .AND. .NOT.l_bloechl) THEN !With bloechl corrections we cannot cut as clean
371 : !Extract basic parameters of the equidistant eMesh
372 400933320 : eb = MINVAL(eMesh)
373 : et = MAXVAL(eMesh)
374 77240 : del = eMesh(2)-eMesh(1)
375 :
376 : !Get last point to the left of the lowest eigenvalue in the tetrahedron
377 77240 : nstart = INT((etetra(ind(1))-eb)/del)+1
378 77240 : nstart = MAX(1,nstart)
379 :
380 : !Get first point to the right of the highest eigenvalue in the tetrahedron
381 77240 : nend = INT((etetra(ind(SIZE(etetra)))-eb)/del)+2
382 77240 : nend = MIN(SIZE(eMesh),nend)
383 77240 : nend = MAX(1,nend)
384 : ELSE
385 25068210 : nstart = 1
386 25068210 : nend = SIZE(eMesh)
387 : ENDIF
388 :
389 194421342 : IF(nstart /= 1) weight(:nstart-1) = 0.0
390 : !Calculate the weights
391 62097536 : DO ie = nstart, nend
392 221880378 : weight(ie) = tetraWeight(eMesh(ie),etetra(ind),icornSorted,vol,film)
393 155547776 : IF(l_bloechl) weight(ie) = weight(ie) + bloechl(eMesh(ie),etetra(ind),icornSorted,vol,film)
394 : ENDDO
395 :
396 : !The loop terminates if the energy is larger than
397 : !all eigenvalues at the tetrahedron/triangle corners (nend)
398 : !For all consecutive values the weight is constant
399 219541832 : IF(nend /= SIZE(eMesh)) weight(nend+1:) = weight(nend)
400 : ENDIF
401 :
402 25068210 : END FUNCTION getWeightSingleBand
403 :
404 : END MODULE m_tetrahedronInit
|
