Line data Source code
1 : MODULE m_resWeight
2 :
3 : IMPLICIT NONE
4 :
5 : PRIVATE
6 : PUBLIC :: resWeight
7 :
8 : CONTAINS
9 :
10 0 : PURE FUNCTION resWeight(eMesh,etetra,ind,vol,film) Result(weights)
11 :
12 : REAL, INTENT(IN) :: eMesh(:)
13 : REAL, INTENT(IN) :: etetra(:)
14 : INTEGER, INTENT(IN) :: ind
15 : REAL, INTENT(IN) :: vol
16 : LOGICAL, INTENT(IN) :: film
17 :
18 : REAL :: weights(SIZE(eMesh))
19 : INTEGER ie
20 :
21 0 : IF(.not.film) then
22 0 : weights = resWeightBulk(eMesh,etetra,ind) * vol
23 : ELSE
24 0 : DO ie = 1, SIZE(eMesh)
25 0 : IF(film) THEN
26 0 : weights(ie) = resWeightFilm(eMesh(ie),etetra,ind) * vol
27 : ENDIF
28 : ENDDO
29 : ENDIF
30 :
31 0 : END FUNCTION resWeight
32 :
33 0 : PURE FUNCTION resWeightBulk(eMesh,e,ind) Result(weights)
34 :
35 : !Calculates resolvent weights for point in tetrahedron (Compare PhysRevB.29.3430)
36 :
37 : REAL, INTENT(IN) :: eMesh(:)
38 : REAL, INTENT(IN) :: e(:)
39 : INTEGER, INTENT(IN) :: ind
40 :
41 : REAL :: weights(SIZE(eMesh))
42 :
43 : REAL, PARAMETER :: tol = 1e-7!tolerance for degeneracy
44 : REAL, PARAMETER :: fac = 5.0
45 : REAL, PARAMETER :: highEnergyCut = 0.1
46 :
47 : REAL a,b,cut,min,z,denom(4)
48 : REAL eShift
49 : REAL eDeg(4)
50 : INTEGER lowerCut,upperCut,ie
51 : INTEGER ndeg,i,j,k,l,m
52 : INTEGER ideg(6,2)
53 :
54 0 : REAL, ALLOCATABLE :: eMesh_shifted(:)
55 : REAL :: eig_shifted(4)
56 :
57 0 : eShift = 0.0
58 0 : if (MAXVAL(eMesh)*MINVAL(eMesh) < 0.0) then
59 : !0 is inside the interval
60 : !The asymptotic relations for energies far from the
61 : !eigenvalues at the corners has a pole for energies at 0
62 : !To prevent these we shift everything (eignevalues and energy grid so that zero)
63 : !is far away. This does not change the results for all other cases
64 : !Since there only energy differences contribute
65 0 : eShift = 2.0*MAX(ABS(MINVAL(eMesh)),ABS(MAXVAL(eMesh)))
66 : endif
67 :
68 0 : weights = 0.0
69 : cut = 0.0
70 : !
71 0 : DO i = 1, 3
72 0 : DO j=i+1,4
73 0 : IF(ABS(e(i)-e(j)).GT.cut) cut = MAX(ABS(e(i)-e(j)),tol)
74 : ENDDO
75 : ENDDO
76 : !Determine the cutoffs
77 0 : upperCut = SIZE(eMesh)
78 0 : DO
79 0 : IF(MINVAL(ABS(eMesh(upperCut)-e)).LT.fac*cut.OR.upperCut.EQ.1) THEN
80 : EXIT
81 0 : ELSE IF(MAXVAL(eMesh(upperCut)-e) < 0.0) THEN
82 : EXIT
83 : ELSE
84 0 : upperCut = upperCut - 1
85 : ENDIF
86 : ENDDO
87 : lowerCut = 1
88 0 : DO
89 0 : IF(MINVAL(ABS(eMesh(lowerCut)-e)).LT.fac*cut.OR.lowerCut.EQ.upperCut) THEN
90 : EXIT
91 0 : ELSE IF(MINVAL(eMesh(lowerCut)-e) > 0.0) THEN
92 : EXIT
93 : ELSE
94 0 : lowerCut = lowerCut + 1
95 : ENDIF
96 : ENDDO
97 0 : eig_shifted = e(:) - eShift
98 :
99 0 : ALLOCATE(eMesh_shifted,mold=eMesh)
100 0 : eMesh_shifted = eMesh(:) - eShift
101 :
102 0 : ndeg = 0
103 0 : ideg = 0
104 : !Search for degenerate eigenvalues
105 0 : eDeg = eig_shifted
106 0 : DO i = 1, 3
107 0 : DO j = i+1,4
108 0 : IF(ABS(eig_shifted(i)-eig_shifted(j)).LT.tol) THEN
109 0 : ndeg = ndeg + 1
110 0 : ideg(ndeg,1) = i
111 0 : ideg(ndeg,2) = j
112 : !Set the two values equal
113 0 : eDeg(i) = (eig_shifted(i)+eig_shifted(j))/2.0
114 0 : eDeg(j) = eDeg(i)
115 0 : DO k = 1, ndeg-1
116 0 : IF(ideg(k,1).EQ.i.OR.ideg(k,1).EQ.j) THEN
117 0 : eDeg(ideg(k,1)) = (eig_shifted(i)+eig_shifted(j)+eig_shifted(ideg(k,1)))/3.0
118 0 : eDeg(j) = eDeg(ideg(k,1))
119 0 : eDeg(i) = eDeg(ideg(k,1))
120 0 : ELSE IF(ideg(k,2).EQ.i.OR.ideg(k,2).EQ.j) THEN
121 0 : eDeg(ideg(k,2)) = (eig_shifted(i)+eig_shifted(j)+eig_shifted(ideg(k,2)))/3.0
122 0 : eDeg(j) = eDeg(ideg(k,2))
123 0 : eDeg(i) = eDeg(ideg(k,2))
124 : ENDIF
125 : ENDDO
126 : ENDIF
127 : ENDDO
128 : ENDDO
129 0 : ndeg = 0
130 0 : ideg = 0
131 0 : DO i = 1, 3
132 0 : DO j = i+1,4
133 0 : IF(ABS(eDeg(i)-eDeg(j)).LT.tol) THEN
134 0 : ndeg = ndeg + 1
135 0 : ideg(ndeg,1) = i
136 0 : ideg(ndeg,2) = j
137 : ENDIF
138 : ENDDO
139 : ENDDO
140 :
141 0 : IF(MAXVAL(MAXVAL(eMesh)-e)<-highEnergyCut) THEN
142 0 : lowerCut=SIZE(eMesh)
143 0 : upperCut=SIZE(eMesh)
144 : ENDIF
145 :
146 0 : IF(MINVAL(MINVAL(eMesh)-e)>highEnergyCut) THEN
147 0 : lowerCut=1
148 0 : upperCut=1
149 : ENDIF
150 :
151 : !asymptotic relation Eqs. 9-11
152 : !the formulas below are unstable for big z arguments
153 0 : a = (eDeg(ind) + SUM(eDeg(:)))/5.0
154 0 : b = 0.0
155 0 : DO j = 1, 4
156 0 : IF(j.EQ.ind) CYCLE
157 0 : b = b + 3*(eDeg(j)-eDeg(ind))**2
158 0 : DO k = 1, 4
159 0 : IF(k.EQ.j) CYCLE
160 0 : b = b + (eDeg(k)-eDeg(j))**2
161 : ENDDO
162 : ENDDO
163 0 : b = b/300.0
164 :
165 0 : weights( :lowerCut) = 0.25/(eMesh_shifted( :lowerCut)-a-b/eMesh_shifted( :lowerCut))
166 0 : weights(upperCut:) = 0.25/(eMesh_shifted(upperCut:)-a-b/eMesh_shifted(upperCut:))
167 :
168 0 : IF(lowerCut==upperCut) RETURN
169 :
170 : IF(ndeg.EQ.0) THEN
171 0 : denom = 1.0
172 0 : DO i = 1, 4
173 0 : DO j = 1, 4
174 0 : IF(i.EQ.j) CYCLE
175 0 : denom(i) = denom(i)*(eDeg(j)-eDeg(i))
176 : ENDDO
177 0 : IF(i.EQ.ind) CYCLE
178 0 : denom(i) = denom(i)*(eDeg(ind)-eDeg(i))
179 : ENDDO
180 0 : DO ie = lowerCut, upperCut
181 0 : z = eMesh_shifted(ie)
182 : !all eigenvalues are non degenerate (EQ.7 from the paper)
183 : !First term
184 0 : weights(ie) = (z-eDeg(ind))**2/denom(ind)
185 0 : DO j = 1, 4
186 0 : IF(j.EQ.ind) CYCLE
187 0 : weights(ie) = weights(ie) + (z-eDeg(j))**3/denom(j)*LOG(ABS((z-eDeg(j)))/ABS((z-eDeg(ind))))
188 : ENDDO
189 : ENDDO
190 :
191 : ELSE IF(ndeg.EQ.1) THEN
192 :
193 : !Two eigenvalues are degenerate (EQs. A1 A2)
194 :
195 0 : l = ideg(1,1)
196 0 : m = ideg(1,2)
197 :
198 0 : IF(ind.EQ.l.OR.ind.EQ.m) THEN
199 : !EQ. A2
200 :
201 : !Get the two indices different from l,m and ind
202 : j = 0
203 0 : DO i = 1, 4
204 0 : IF(i.EQ.l.OR.i.EQ.m) CYCLE
205 0 : j = i
206 : ENDDO
207 : k = 0
208 0 : DO i = 1, 4
209 0 : IF(i.EQ.l.OR.i.EQ.m.OR.i.EQ.j) CYCLE
210 0 : k = i
211 : ENDDO
212 0 : DO ie = lowerCut, upperCut
213 0 : z = eMesh_shifted(ie)
214 : weights(ie) = (z-eDeg(k))**3/((eDeg(k)-eDeg(j))*(eDeg(k)-eDeg(m))**3)*LOG(ABS(z-eDeg(k))) &
215 : +(z-eDeg(j))**3/((eDeg(j)-eDeg(k))*(eDeg(j)-eDeg(m))**3)*LOG(ABS(z-eDeg(j))) &
216 : +(z-eDeg(m))/((eDeg(m)-eDeg(j))*(eDeg(m)-eDeg(k))) * (0.5 + (z-eDeg(j))/(eDeg(m)-eDeg(j))&
217 : +(z-eDeg(k))/(eDeg(m)-eDeg(k)) + ((z-eDeg(j))**2/(eDeg(m)-eDeg(j))**2 &
218 : +(z-eDeg(k))**2/(eDeg(m)-eDeg(k))**2 +(z-eDeg(j))/(eDeg(m)-eDeg(j))*(z-eDeg(k))/(eDeg(m)-eDeg(k))) &
219 0 : *LOG(ABS(z-eDeg(m))))
220 : ENDDO
221 : ELSE
222 : !k is the one site not equal to ind, l or m
223 : k = 0
224 0 : DO i = 1, 4
225 0 : IF(i.EQ.ind.OR.i.EQ.l.OR.i.EQ.m) CYCLE
226 0 : k = i
227 : ENDDO
228 : !EQ. A1
229 0 : DO ie = lowerCut, upperCut
230 0 : z = eMesh_shifted(ie)
231 : weights(ie) = (z-eDeg(ind))**2/((eDeg(ind)-eDeg(m))**2*(eDeg(k)-eDeg(ind))) *&
232 : (1 + (2*(z-eDeg(m))/(eDeg(ind)-eDeg(m))+(z-eDeg(k))/(eDeg(ind)-eDeg(k)))*LOG(ABS(z-eDeg(ind)))) &
233 : +(z-eDeg(m))**2/((eDeg(m)-eDeg(ind))**2*(eDeg(k)-eDeg(m))) *&
234 : (1 + (2*(z-eDeg(ind))/(eDeg(m)-eDeg(ind))+(z-eDeg(k))/(eDeg(m)-eDeg(k)))*LOG(ABS(z-eDeg(m)))) &
235 0 : +(z-eDeg(k))**3/((eDeg(k)-eDeg(ind))**2*(eDeg(k)-eDeg(m))**2)*LOG(ABS(z-eDeg(k)))
236 : ENDDO
237 : ENDIF
238 : ELSE IF(ndeg.EQ.2) THEN
239 : !This is the case E1=E2<E3=E4 => A4
240 0 : IF(ind.LE.2) THEN
241 0 : DO ie = lowerCut, upperCut
242 0 : z = eMesh_shifted(ie)
243 : weights(ie) = 3.0*(z-eDeg(3))**2*(z-eDeg(2))/(eDeg(3)-eDeg(2))**4*LOG(ABS((z-eDeg(2)))/ABS((z-eDeg(3)))) &
244 0 : - 3.0/2.0*(z-eDeg(2))*(2*(z-eDeg(3))-eDeg(3)+eDeg(2))/(eDeg(3)-eDeg(2))**3-1.0/(eDeg(3)-eDeg(2))
245 : ENDDO
246 : ELSE
247 0 : DO ie = lowerCut, upperCut
248 0 : z = eMesh_shifted(ie)
249 : weights(ie) = 3.0*(z-eDeg(2))**2*(z-eDeg(3))/(eDeg(3)-eDeg(2))**4*LOG(ABS((z-eDeg(3)))/ABS((z-eDeg(2)))) &
250 0 : + 3.0/2.0*(z-eDeg(3))*(2*(z-eDeg(2))+eDeg(3)-eDeg(2))/(eDeg(3)-eDeg(2))**3+1.0/(eDeg(3)-eDeg(2))
251 : ENDDO
252 : ENDIF
253 : ELSE IF(ndeg.GE.3 .AND. ndeg.LT.6) THEN
254 : !EQs A3/A5 (here we explicitly write each weight)
255 0 : IF(ALL(ideg(:,:).NE.4)) THEN
256 : !A3
257 0 : IF(ind.NE.4) THEN
258 0 : DO ie = lowerCut, upperCut
259 0 : z = eMesh_shifted(ie)
260 : weights(ie) = (z-eDeg(4))**3/(eDeg(4)-eDeg(3))**4*LOG(ABS((z-eDeg(4))/(z-eDeg(3)))) &
261 0 : + (6*(z-eDeg(4))**2-3*(eDeg(4)-eDeg(3))*(z-eDeg(4))+2*(eDeg(4)-eDeg(3))**2)/(6*(eDeg(4)-eDeg(3))**3)
262 : ENDDO
263 : ELSE
264 0 : DO ie = lowerCut, upperCut
265 0 : z = eMesh_shifted(ie)
266 : weights(ie) = 3.0*(z-eDeg(4))**2*(z-eDeg(3))/(eDeg(4)-eDeg(3))**4*LOG(ABS((z-eDeg(3))/(z-eDeg(4)))) &
267 0 : - 3.0/2.0*(z-eDeg(3))*(2*(z-eDeg(4))-eDeg(4)+eDeg(3))/(eDeg(4)-eDeg(3))**3-1.0/(eDeg(4)-eDeg(3))
268 : ENDDO
269 : ENDIF
270 0 : ELSE IF(ALL(ideg(:,:).NE.1)) THEN
271 : !A5
272 0 : IF(ind.EQ.1) THEN
273 0 : DO ie = lowerCut, upperCut
274 0 : z = eMesh_shifted(ie)
275 : weights(ie) = 3.0*(z-eDeg(1))**2*(z-eDeg(2))/(eDeg(2)-eDeg(1))**4*LOG(ABS((z-eDeg(2))/(z-eDeg(1)))) &
276 0 : + 3.0/2.0*(z-eDeg(2))*(2*(z-eDeg(1))-eDeg(1)+eDeg(2))/(eDeg(2)-eDeg(1))**3+1.0/(eDeg(2)-eDeg(1))
277 : ENDDO
278 : ELSE
279 0 : DO ie = lowerCut, upperCut
280 0 : z = eMesh_shifted(ie)
281 : weights(ie) = (z-eDeg(1))**3/(eDeg(2)-eDeg(1))**4*LOG(ABS((z-eDeg(1))/(z-eDeg(2)))) &
282 0 : - (6*(z-eDeg(1))**2+3*(z-eDeg(1))*(eDeg(2)-eDeg(1))+2*(eDeg(2)-eDeg(1))**2)/(6*(eDeg(2)-eDeg(1))**3)
283 : ENDDO
284 : ENDIF
285 : ENDIF
286 : ELSE IF(ndeg.EQ.6) THEN
287 : !Eq. A6
288 0 : weights(lowerCut:upperCut) = 0.25/(eMesh_shifted(lowerCut:upperCut)-eDeg(1))
289 : ENDIF
290 0 : END FUNCTION resWeightBulk
291 :
292 0 : PURE REAL FUNCTION resWeightFilm(z,e,ind)
293 :
294 : !Calculates resolvent weights for point in triangle (Compare Solid State Phys. 20 (1987) 4823-4831.)
295 :
296 : REAL, INTENT(IN) :: z
297 : REAL, INTENT(IN) :: e(:)
298 : INTEGER, INTENT(IN) :: ind
299 :
300 : REAL, PARAMETER :: tol = 1e-8!tolerance for degeneracy
301 : REAL, PARAMETER :: fac = 4.0
302 :
303 : REAL denom,a,b,cut,min,prod
304 : INTEGER ndeg,i,j,k,l,m
305 : INTEGER ideg(2,2)
306 :
307 0 : min = 9e+20
308 0 : cut = 0.0
309 0 : DO i = 1, 2
310 0 : DO j=i+1,3
311 0 : IF(ABS(e(i)-e(j)).GT.cut) cut = MAX(ABS(e(i)-e(j)),tol)
312 : ENDDO
313 : ENDDO
314 :
315 0 : DO i =1, 3
316 0 : IF(ABS(z-e(i)).LT.min) min = ABS(z-e(i))
317 : ENDDO
318 0 : ndeg = 0
319 0 : ideg = 0
320 0 : IF(min.GT.fac*cut) THEN
321 : !asymptotic relation Eqs. 9-11
322 : !the formulas below are unstable for big z arguments
323 0 : a = (e(ind) + SUM(e(:)))/4.0
324 0 : b = 3*e(ind)**2
325 0 : prod = 1.0
326 0 : DO j = 1,3
327 0 : IF(j.EQ.ind) CYCLE
328 0 : b = b + 2*e(ind)*e(j)+e(j)**2
329 0 : prod = prod * e(j)
330 : ENDDO
331 0 : b = (b+prod)/300.0
332 0 : resWeightFilm = 1.0/3.0/(z-a-b/z)
333 : !The asymptotic relationship has a divergence at z=0
334 : !For now we drop the 1/z term in the denominator here
335 0 : IF(ABS(z).LT.2e-2) resWeightFilm = 1.0/3.0/(z-a)
336 : ELSE
337 : !Search for degenerate eigenvalues
338 0 : DO i = 1, 2
339 0 : DO j = i+1,3
340 0 : IF(ABS(e(i)-e(j)).LT.tol) THEN
341 0 : ndeg = ndeg + 1
342 0 : ideg(ndeg,1) = i
343 0 : ideg(ndeg,2) = j
344 : ENDIF
345 : ENDDO
346 : ENDDO
347 0 : resWeightFilm = 0.0
348 0 : IF(ndeg.EQ.0) THEN
349 :
350 : !all eigenvalues are non degenerate
351 : denom = 1.0
352 0 : DO i = 1, 3
353 0 : IF(i.EQ.ind) CYCLE
354 0 : denom = denom*(e(i)-e(ind))
355 : ENDDO
356 : !First term
357 0 : resWeightFilm = (z-e(ind))/denom
358 0 : DO j = 1, 3
359 0 : IF(j.EQ.ind) CYCLE
360 : denom = 1.0
361 0 : DO i = 1, 3
362 0 : IF(i.EQ.j) CYCLE
363 0 : denom = denom*(e(i)-e(j))
364 : ENDDO
365 0 : resWeightFilm = resWeightFilm + (z-e(j))**2/denom*LOG(ABS((z-e(j))/(z-e(ind))))/(e(ind)-e(j))
366 : ENDDO
367 0 : ELSE IF(ndeg.EQ.1) THEN
368 :
369 0 : IF(ALL(ideg(:,:).NE.3)) THEN
370 0 : IF(ind.NE.3) THEN
371 : resWeightFilm = (z-e(3))**2/((e(1)-e(3))**3)*LOG(ABS((z-e(3))/(z-e(1)))) &
372 0 : -(z-e(1))/(e(3)-e(1))**2+1.5/(e(3)-e(1))
373 : ELSE
374 : resWeightFilm = (z-e(3))/((e(1)-e(3))*(e(2)-e(3))) &
375 : - 2*(z-e(1))*(z-e(3))/(e(3)-e(1))**3 * LOG(ABS((z-e(1))/(z-e(3)))) &
376 0 : +(z-e(1))/(e(3)-e(1))**2
377 : ENDIF
378 : ELSE
379 0 : IF(ind.EQ.1) THEN
380 : resWeightFilm = (z-e(1))/((e(2)-e(1))*(e(3)-e(1))) &
381 : - 2*(z-e(1))*(z-e(3))/(e(1)-e(3))**3 * LOG(ABS((z-e(3))/(z-e(1)))) &
382 0 : +(z-e(3))/(e(1)-e(3))**2
383 : ELSE
384 : resWeightFilm = (z-e(1))**2/((e(3)-e(1))**3)*LOG(ABS((z-e(1))/(z-e(3)))) &
385 0 : -(z-e(1))/(e(1)-e(3))**2+0.5/(e(1)-e(3))
386 : ENDIF
387 : ENDIF
388 :
389 : ELSE
390 0 : resWeightFilm = 1.0/3.0 * 1/(z-e(1))
391 : ENDIF
392 : ENDIF
393 :
394 0 : END FUNCTION resWeightFilm
395 :
396 : END MODULE m_resWeight
|
