pwdft  0.1
PW-DFT code for education
libtetrabz_polstat_mod.F90
Go to the documentation of this file.
1 !
2 ! Copyright (c) 2014 Mitsuaki Kawamura
3 !
4 ! Permission is hereby granted, free of charge, to any person obtaining a
5 ! copy of this software and associated documentation files (the
6 ! "Software"), to deal in the Software without restriction, including
7 ! without limitation the rights to use, copy, modify, merge, publish,
8 ! distribute, sublicense, and/or sell copies of the Software, and to
9 ! permit persons to whom the Software is furnished to do so, subject to
10 ! the following conditions:
11 !
12 ! The above copyright notice and this permission notice shall be included
13 ! in all copies or substantial portions of the Software.
14 !
15 ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
16 ! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17 ! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
18 ! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
19 ! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
20 ! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
21 ! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
22 !
23 MODULE libtetrabz_polstat_mod
24  !
25  IMPLICIT NONE
26  !
27  PRIVATE
28  PUBLIC libtetrabz_polstat
29  !
30 CONTAINS
31 !
32 ! Compute Static polalization function
33 !
34 SUBROUTINE libtetrabz_polstat(ltetra,bvec,nb,nge,eig1,eig2,ngw,wght,comm) BIND(C)
35  !
36  USE iso_c_binding
37  USE libtetrabz_common, ONLY : libtetrabz_initialize, libtetrabz_interpol_indx, libtetrabz_mpisum_dv
38  IMPLICIT NONE
39  !
40  INTEGER(C_INT),INTENT(IN) :: ltetra, nb, nge(3), ngw(3)
41  REAL(c_double),INTENT(IN) :: bvec(9), eig1(nb,product(nge(1:3))), eig2(nb,product(nge(1:3)))
42  REAL(c_double),INTENT(OUT) :: wght(nb*nb,product(ngw(1:3)))
43  INTEGER(C_INT),INTENT(IN),OPTIONAL :: comm
44  !
45  LOGICAL :: linterpol
46  INTEGER :: nt_local, nk_local, nkbz, ik, kintp(8)
47  INTEGER,ALLOCATABLE :: ik_global(:,:), ik_local(:,:)
48  REAL(8) :: wlsm(4,20), wintp(1,8)
49  REAL(8),ALLOCATABLE :: wghtd(:,:,:), kvec(:,:)
50  !
51  IF(PRESENT(comm)) THEN
52  CALL libtetrabz_initialize(ltetra,nge,ngw,bvec,linterpol,wlsm,nk_local,&
53  & nt_local,nkbz,ik_global,ik_local,kvec,comm)
54  ELSE
55  CALL libtetrabz_initialize(ltetra,nge,ngw,bvec,linterpol,wlsm,nk_local,&
56  & nt_local,nkbz,ik_global,ik_local,kvec)
57  END IF
58  !
59  IF(linterpol) THEN
60  !
61  ALLOCATE(wghtd(nb*nb,1,nk_local))
62  CALL libtetrabz_polstat_main(wlsm,nt_local,ik_global,ik_local,nb,nkbz,eig1,eig2,nk_local,wghtd)
63  !
64  ! Interpolation
65  !
66  wght(1:nb*nb,1:product(ngw(1:3))) = 0d0
67  DO ik = 1, nk_local
68  CALL libtetrabz_interpol_indx(ngw,kvec(1:3,ik),kintp,wintp)
69  wght(1:nb*nb,kintp(1:8)) = wght(1:nb*nb, kintp(1:8)) &
70  & + matmul(wghtd(1:nb*nb,1:1,ik), wintp(1:1,1:8))
71  END DO ! ik = 1, nk_local
72  DEALLOCATE(wghtd, kvec)
73  !
74  IF(PRESENT(comm)) CALL libtetrabz_mpisum_dv(comm, nb*nb*product(ngw(1:3)), wght)
75  !
76  ELSE
77  CALL libtetrabz_polstat_main(wlsm,nt_local,ik_global,ik_local,nb,nkbz,eig1,eig2,nk_local,wght)
78  END IF
79  !
80  DEALLOCATE(ik_global, ik_local)
81  !
82 END SUBROUTINE libtetrabz_polstat
83 !
84 ! Main SUBROUTINE for polalization function : Theta(- E1) * Theta(E2) / (E2 - E1)
85 !
86 SUBROUTINE libtetrabz_polstat_main(wlsm,nt_local,ik_global,ik_local,nb,nkBZ,eig1,eig2,nk_local,polstat)
87  !
88  USE libtetrabz_common, ONLY : libtetrabz_sort, &
89  & libtetrabz_tsmall_a1, libtetrabz_tsmall_b1, &
90  & libtetrabz_tsmall_b2, libtetrabz_tsmall_b3, &
91  & libtetrabz_tsmall_c1, libtetrabz_tsmall_c2, &
92  & libtetrabz_tsmall_c3
93  IMPLICIT NONE
94  !
95  INTEGER,INTENT(IN) :: nt_local, nb, nkBZ, nk_local, &
96  & ik_global(20,nt_local), ik_local(20,nt_local)
97  REAL(8),INTENT(IN) :: wlsm(4,20), eig1(nb,nkBZ), eig2(nb,nkBZ)
98  REAL(8),INTENT(OUT) :: polstat(nb,nb,nk_local)
99  !
100  INTEGER :: ib, it, indx(4)
101  REAL(8) :: e(4), ei1(4,nb), ei2(4), ej1(4,nb), ej2(4,nb), thr = 1d-10, tsmall(4,4), v, w1(nb,4), w2(nb,4)
102  !
103  polstat(1:nb,1:nb,1:nk_local) = 0d0
104  !
105  !$OMP PARALLEL DEFAULT(NONE) &
106  !$OMP & SHARED(eig1,eig2,ik_global,ik_local,nb,nt_local,polstat,thr,wlsm) &
107  !$OMP & PRIVATE(e,ei1,ei2,ej1,ej2,ib,indx,it,tsmall,V,w1,w2)
108  !
109  DO it = 1, nt_local
110  !
111  DO ib = 1, nb
112  ei1(1:4, ib) = matmul(wlsm(1:4,1:20), eig1(ib, ik_global(1:20,it)))
113  ej1(1:4, ib) = matmul(wlsm(1:4,1:20), eig2(ib, ik_global(1:20,it)))
114  END DO
115  !
116  !$OMP DO
117  DO ib = 1, nb
118  !
119  w1(1:nb,1:4) = 0d0
120  e(1:4) = ei1(1:4, ib)
121  CALL libtetrabz_sort(4,e,indx)
122  !
123  IF(e(1) <= 0d0 .AND. 0d0 < e(2)) THEN
124  !
125  CALL libtetrabz_tsmall_a1(e,v,tsmall)
126  !
127  IF(v > thr) THEN
128  !
129  ei2(1:4 ) = matmul(tsmall(1:4,1:4), ei1(indx(1:4), ib))
130  ej2(1:4,1:nb) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),1:nb))
131  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
132  w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
133  & + v * matmul(w2(1:nb,1:4), tsmall(1:4,1:4))
134  !
135  END IF
136  !
137  ELSE IF(e(2) <= 0d0 .AND. 0d0 < e(3)) THEN
138  !
139  CALL libtetrabz_tsmall_b1(e,v,tsmall)
140  !
141  IF(v > thr) THEN
142  !
143  ei2(1:4 ) = matmul(tsmall(1:4,1:4), ei1(indx(1:4), ib))
144  ej2(1:4,1:nb) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),1:nb))
145  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
146  w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
147  & + v * matmul(w2(1:nb,1:4), tsmall(1:4,1:4))
148  !
149  END IF
150  !
151  CALL libtetrabz_tsmall_b2(e,v,tsmall)
152  !
153  IF(v > thr) THEN
154  !
155  ei2(1:4 ) = matmul(tsmall(1:4,1:4), ei1(indx(1:4), ib))
156  ej2(1:4,1:nb) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),1:nb))
157  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
158  w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
159  & + v * matmul(w2(1:nb,1:4), tsmall(1:4,1:4))
160  !
161  END IF
162  !
163  CALL libtetrabz_tsmall_b3(e,v,tsmall)
164  !
165  IF(v > thr) THEN
166  !
167  ei2(1:4 ) = matmul(tsmall(1:4,1:4), ei1(indx(1:4), ib))
168  ej2(1:4,1:nb) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),1:nb))
169  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
170  w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
171  & + v * matmul(w2(1:nb,1:4), tsmall(1:4,1:4))
172  !
173  END IF
174  !
175  ELSE IF(e(3) <= 0d0 .AND. 0d0 < e(4)) THEN
176  !
177  CALL libtetrabz_tsmall_c1(e,v,tsmall)
178  !
179  IF(v > thr) THEN
180  !
181  ei2(1:4 ) = matmul(tsmall(1:4,1:4), ei1(indx(1:4), ib))
182  ej2(1:4,1:nb) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),1:nb))
183  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
184  w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
185  & + v * matmul(w2(1:nb,1:4), tsmall(1:4,1:4))
186  !
187  END IF
188  !
189  CALL libtetrabz_tsmall_c2(e,v,tsmall)
190  !
191  IF(v > thr) THEN
192  !
193  ei2(1:4 ) = matmul(tsmall(1:4,1:4), ei1(indx(1:4), ib))
194  ej2(1:4,1:nb) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),1:nb))
195  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
196  w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
197  & + v * matmul(w2(1:nb,1:4), tsmall(1:4,1:4))
198  !
199  END IF
200  !
201  CALL libtetrabz_tsmall_c3(e,v,tsmall)
202  !
203  IF(v > thr) THEN
204  !
205  ei2(1:4 ) = matmul(tsmall(1:4,1:4), ei1(indx(1:4), ib))
206  ej2(1:4,1:nb) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),1:nb))
207  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
208  w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
209  & + v * matmul(w2(1:nb,1:4), tsmall(1:4,1:4))
210  !
211  END IF
212  !
213  ELSE IF(e(4) <= 0d0) THEN
214  !
215  ei2(1:4 ) = ei1(1:4, ib)
216  ej2(1:4,1:nb) = ej1(1:4,1:nb)
217  CALL libtetrabz_polstat2(nb,ei2,ej2,w2)
218  w1(1:nb,1:4) = w1(1:nb,1:4) + w2(1:nb,1:4)
219  !
220  ELSE
221  !
222  cycle
223  !
224  END IF
225  !
226  polstat(1:nb,ib,ik_local(1:20,it)) = polstat(1:nb,ib, ik_local(1:20,it)) &
227  & + matmul(w1(1:nb,1:4), wlsm(1:4,1:20))
228  !
229  END DO ! ib
230  !$OMP END DO NOWAIT
231  !
232  END DO ! it
233  !
234  !$OMP END PARALLEL
235  !
236  polstat(1:nb,1:nb,1:nk_local) = polstat(1:nb,1:nb,1:nk_local) / dble(6 * nkbz)
237  !
238 END SUBROUTINE libtetrabz_polstat_main
239 !
240 ! Tetrahedra method for theta( - E2)
241 !
242 SUBROUTINE libtetrabz_polstat2(nb,ei1,ej1,w1)
243  !
244  USE libtetrabz_common, ONLY : libtetrabz_sort, &
245  & libtetrabz_tsmall_a1, libtetrabz_tsmall_b1, &
246  & libtetrabz_tsmall_b2, libtetrabz_tsmall_b3, &
247  & libtetrabz_tsmall_c1, libtetrabz_tsmall_c2, &
248  & libtetrabz_tsmall_c3
249  IMPLICIT NONE
250  !
251  INTEGER,INTENT(IN) :: nb
252  REAL(8),INTENT(IN) :: ei1(4), ej1(4,nb)
253  REAL(8),INTENT(INOUT) :: w1(nb,4)
254  !
255  INTEGER :: ib, indx(4)
256  REAL(8) :: de(4), e(4), thr = 1d-8, tsmall(4,4), v, w2(4)
257  !
258  DO ib = 1, nb
259  !
260  w1(ib,1:4) = 0d0
261  e(1:4) = - ej1(1:4, ib)
262  CALL libtetrabz_sort(4,e,indx)
263  !
264  IF((e(1) <= 0d0 .AND. 0d0 < e(2)) .OR. (e(1) < 0d0 .AND. 0d0 <= e(2))) THEN
265  !
266  CALL libtetrabz_tsmall_a1(e,v,tsmall)
267  !
268  IF(v > thr) THEN
269  !
270  de(1:4) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),ib) - ei1(indx(1:4)))
271  CALL libtetrabz_polstat3(de,w2)
272  w1(ib,indx(1:4)) = w1(ib, indx(1:4)) &
273  & + v * matmul(w2(1:4), tsmall(1:4,1:4))
274  !
275  END IF
276  !
277  ELSE IF((e(2) <= 0d0 .AND. 0d0 < e(3)) .OR. (e(2) < 0d0 .AND. 0d0 <= e(3))) THEN
278  !
279  CALL libtetrabz_tsmall_b1(e,v,tsmall)
280  !
281  IF(v > thr) THEN
282  !
283  de(1:4) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),ib) - ei1(indx(1:4)))
284  CALL libtetrabz_polstat3(de,w2)
285  w1(ib,indx(1:4)) = w1(ib,indx(1:4)) &
286  & + v * matmul(w2( 1:4 ), tsmall(1:4,1:4))
287  !
288  END IF
289  !
290  CALL libtetrabz_tsmall_b2(e,v,tsmall)
291  !
292  IF(v > thr) THEN
293  !
294  de(1:4) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),ib) - ei1(indx(1:4)))
295  CALL libtetrabz_polstat3(de,w2)
296  w1(ib,indx(1:4)) = w1(ib, indx(1:4)) &
297  & + v * matmul(w2(1:4), tsmall(1:4,1:4))
298  !
299  END IF
300  !
301  CALL libtetrabz_tsmall_b3(e,v,tsmall)
302  !
303  IF(v > thr) THEN
304  !
305  de(1:4) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),ib) - ei1(indx(1:4)))
306  CALL libtetrabz_polstat3(de,w2)
307  w1(ib,indx(1:4)) = w1(ib, indx(1:4)) &
308  & + v * matmul(w2(1:4), tsmall(1:4,1:4))
309  !
310  END IF
311  !
312  ELSE IF((e(3) <= 0d0 .AND. 0d0 < e(4)) .OR. (e(3) < 0d0 .AND. 0d0 <= e(4))) THEN
313  !
314  CALL libtetrabz_tsmall_c1(e,v,tsmall)
315  !
316  IF(v > thr) THEN
317  !
318  de(1:4) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),ib) - ei1(indx(1:4)))
319  CALL libtetrabz_polstat3(de,w2)
320  w1(ib,indx(1:4)) = w1(ib, indx(1:4)) &
321  & + v * matmul(w2(1:4), tsmall(1:4,1:4))
322  !
323  END IF
324  !
325  CALL libtetrabz_tsmall_c2(e,v,tsmall)
326  !
327  IF(v > thr) THEN
328  !
329  de(1:4) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),ib) - ei1(indx(1:4)))
330  CALL libtetrabz_polstat3(de,w2)
331  w1(ib,indx(1:4)) = w1(ib, indx(1:4)) &
332  & + v * matmul(w2(1:4), tsmall(1:4,1:4))
333  !
334  END IF
335  !
336  CALL libtetrabz_tsmall_c3(e,v,tsmall)
337  !
338  IF(v > thr) THEN
339  !
340  de(1:4) = matmul(tsmall(1:4,1:4), ej1(indx(1:4),ib) - ei1(indx(1:4)))
341  CALL libtetrabz_polstat3(de,w2)
342  w1(ib,indx(1:4)) = w1(ib, indx(1:4)) &
343  & + v * matmul(w2(1:4), tsmall(1:4,1:4))
344  !
345  END IF
346  !
347  ELSE IF( e(4) <= 0d0 ) THEN
348  !
349  de(1:4) = ej1(1:4,ib) - ei1(1:4)
350  CALL libtetrabz_polstat3(de,w2)
351  w1(ib,1:4) = w1(ib,1:4) + w2(1:4)
352  !
353  END IF
354  !
355  END DO ! ib = 1, nb
356  !
357 END SUBROUTINE libtetrabz_polstat2
358 !
359 ! Tetarahedra method for delta(om - ep + e)
360 !
361 SUBROUTINE libtetrabz_polstat3(de,w1)
362  !
363  USE libtetrabz_common, ONLY : libtetrabz_sort
364  IMPLICIT NONE
365  !
366  REAL(8),INTENT(IN) :: de(4)
367  REAL(8),INTENT(INOUT) :: w1(4)
368  !
369  INTEGER :: ii, indx(4)
370  REAL(8) :: e(4), ln(4), thr, thr2
371  !
372  e(1:4) = de(1:4)
373  CALL libtetrabz_sort(4,e,indx)
374  !
375  thr = maxval(e(1:4)) * 1d-3
376  thr2 = 1d-8
377  !
378  DO ii = 1, 4
379  IF(e(ii) < thr2) THEN
380  IF(ii == 3) THEN
381  stop " Nesting ! "
382  END IF
383  ln(ii) = 0d0
384  e(ii) = 0d0
385  ELSE
386  ln(ii) = log(e(ii))
387  END IF
388  END DO
389  !
390  IF(abs(e(4) - e(3)) < thr ) THEN
391  IF(abs(e(4) - e(2)) < thr ) THEN
392  IF(abs(e(4) - e(1)) < thr ) THEN
393  !
394  ! e(4) = e(3) = e(2) = e(1)
395  !
396  w1(indx(4)) = 0.25d0 / e(4)
397  w1(indx(3)) = w1(indx(4))
398  w1(indx(2)) = w1(indx(4))
399  w1(indx(1)) = w1(indx(4))
400  !
401  ELSE
402  !
403  ! e(4) = e(3) = e(2)
404  !
405  w1(indx(4)) = libtetrabz_polstat_1211(e(4),e(1),ln(4),ln(1))
406  w1(indx(3)) = w1(indx(4))
407  w1(indx(2)) = w1(indx(4))
408  w1(indx(1)) = libtetrabz_polstat_1222(e(1),e(4),ln(1),ln(4))
409  !
410  IF(any(w1(1:4) < 0d0)) THEN
411  WRITE(*,'(100e15.5)') e(1:4)
412  WRITE(*,'(100e15.5)') w1(indx(1:4))
413  stop "weighting 4=3=2"
414  END IF
415  !
416  END IF
417  ELSE IF(abs(e(2) - e(1)) < thr) THEN
418  !
419  ! e(4) = e(3), e(2) = e(1)
420  !
421  w1(indx(4)) = libtetrabz_polstat_1221(e(4),e(2), ln(4),ln(2))
422  w1(indx(3)) = w1(indx(4))
423  w1(indx(2)) = libtetrabz_polstat_1221(e(2),e(4), ln(2),ln(4))
424  w1(indx(1)) = w1(indx(2))
425  !
426  IF(any(w1(1:4) < 0d0)) THEN
427  WRITE(*,'(100e15.5)') e(1:4)
428  WRITE(*,'(100e15.5)') w1(indx(1:4))
429  stop "weighting 4=3 2=1"
430  END IF
431  !
432  ELSE
433  !
434  ! e(4) = e(3)
435  !
436  w1(indx(4)) = libtetrabz_polstat_1231(e(4),e(1),e(2),ln(4),ln(1),ln(2))
437  w1(indx(3)) = w1(indx(4))
438  w1(indx(2)) = libtetrabz_polstat_1233(e(2),e(1),e(4),ln(2),ln(1),ln(4))
439  w1(indx(1)) = libtetrabz_polstat_1233(e(1),e(2),e(4),ln(1),ln(2),ln(4))
440  !
441  IF(any(w1(1:4) < 0d0)) THEN
442  WRITE(*,'(100e15.5)') e(1:4)
443  WRITE(*,'(100e15.5)') w1(indx(1:4))
444  stop "weighting 4=3"
445  END IF
446  !
447  END IF
448  ELSE IF(abs(e(3) - e(2)) < thr) THEN
449  IF(abs(e(3) - e(1)) < thr) THEN
450  !
451  ! e(3) = e(2) = e(1)
452  !
453  w1(indx(4)) = libtetrabz_polstat_1222(e(4),e(3), ln(4),ln(3))
454  w1(indx(3)) = libtetrabz_polstat_1211(e(3),e(4), ln(3),ln(4))
455  w1(indx(2)) = w1(indx(3))
456  w1(indx(1)) = w1(indx(3))
457  !
458  IF(any(w1(1:4) < 0d0)) THEN
459  WRITE(*,'(100e15.5)') e(1:4)
460  WRITE(*,'(100e15.5)') w1(indx(1:4))
461  stop "weighting 3=2=1"
462  END IF
463  !
464  ELSE
465  !
466  ! e(3) = e(2)
467  !
468  w1(indx(4)) = libtetrabz_polstat_1233(e(4),e(1),e(3),ln(4),ln(1),ln(3))
469  w1(indx(3)) = libtetrabz_polstat_1231(e(3),e(1),e(4),ln(3),ln(1),ln(4))
470  w1(indx(2)) = w1(indx(3))
471  w1(indx(1)) = libtetrabz_polstat_1233(e(1),e(4),e(3),ln(1),ln(4),ln(3))
472  !
473  IF(any(w1(1:4) < 0d0)) THEN
474  WRITE(*,'(100e15.5)') e(1:4)
475  WRITE(*,'(100e15.5)') w1(indx(1:4))
476  stop "weighting 3=2"
477  END IF
478  !
479  END IF
480  ELSE IF(abs(e(2) - e(1)) < thr) THEN
481  !
482  ! e(2) = e(1)
483  !
484  w1(indx(4)) = libtetrabz_polstat_1233(e(4),e(3),e(2),ln(4),ln(3),ln(2))
485  w1(indx(3)) = libtetrabz_polstat_1233(e(3),e(4),e(2),ln(3),ln(4),ln(2))
486  w1(indx(2)) = libtetrabz_polstat_1231(e(2),e(3),e(4),ln(2),ln(3),ln(4))
487  w1(indx(1)) = w1(indx(2))
488  !
489  IF(any(w1(1:4) < 0d0)) THEN
490  WRITE(*,'(100e15.5)') e(1:4)
491  WRITE(*,'(100e15.5)') w1(indx(1:4))
492  stop "weighting 2=1"
493  END IF
494  !
495  ELSE
496  !
497  ! Different each other.
498  !
499  w1(indx(4)) = libtetrabz_polstat_1234(e(4),e(1),e(2),e(3),ln(4),ln(1),ln(2),ln(3))
500  w1(indx(3)) = libtetrabz_polstat_1234(e(3),e(1),e(2),e(4),ln(3),ln(1),ln(2),ln(4))
501  w1(indx(2)) = libtetrabz_polstat_1234(e(2),e(1),e(3),e(4),ln(2),ln(1),ln(3),ln(4))
502  w1(indx(1)) = libtetrabz_polstat_1234(e(1),e(2),e(3),e(4),ln(1),ln(2),ln(3),ln(4))
503  !
504  IF(any(w1(1:4) < 0d0)) THEN
505  WRITE(*,'(100e15.5)') e(1:4)
506  WRITE(*,'(100e15.5)') w1(indx(1:4))
507  stop "weighting"
508  END IF
509  !
510  END IF
511  !
512 END SUBROUTINE libtetrabz_polstat3
513 !
514 ! Results of Integration (1-x-y-z)/(g0+(g1-g0)x+(g2-g0)y+(g3-g0))
515 ! for 0<x<1, 0<y<1-x, 0<z<1-x-y
516 !
517 ! 1, Different each other
518 !
519 FUNCTION libtetrabz_polstat_1234(g1,g2,g3,g4,lng1,lng2,lng3,lng4) RESULT(w)
520  !
521  IMPLICIT NONE
522  !
523  REAL(8),INTENT(IN) :: g1,g2,g3,g4,lng1,lng2,lng3,lng4
524  REAL(8) :: w
525  !
526  REAL(8) :: w2, w3, w4
527  !
528  w2 = ((lng2 - lng1)/(g2 - g1)*g2 - 1d0)*g2/(g2 - g1)
529  w3 = ((lng3 - lng1)/(g3 - g1)*g3 - 1d0)*g3/(g3 - g1)
530  w4 = ((lng4 - lng1)/(g4 - g1)*g4 - 1d0)*g4/(g4 - g1)
531  w2 = ((w2 - w3)*g2)/(g2 - g3)
532  w4 = ((w4 - w3)*g4)/(g4 - g3)
533  w = (w4 - w2)/(g4 - g2)
534  !
535 END FUNCTION libtetrabz_polstat_1234
536 !
537 ! 2, g4 = g1
538 !
539 FUNCTION libtetrabz_polstat_1231(g1,g2,g3,lng1,lng2,lng3) RESULT(w)
540  !
541  IMPLICIT NONE
542  !
543  REAL(8),INTENT(IN) :: g1,g2,g3,lng1,lng2,lng3
544  REAL(8) :: w
545  !
546  REAL(8) :: w2, w3
547  !
548  w2 = ((lng2 - lng1)/(g2 - g1)*g2 - 1d0)*g2**2/(g2 - g1) - g1/( &
549  & 2d0)
550  w2 = w2/(g2 - g1)
551  w3 = ((lng3 - lng1)/(g3 - g1)*g3 - 1d0)*g3**2/(g3 - g1) - g1/( &
552  & 2d0)
553  w3 = w3/(g3 - g1)
554  w = (w3 - w2)/(g3 - g2)
555  !
556 END FUNCTION libtetrabz_polstat_1231
557 !
558 ! 3, g4 = g3
559 !
560 FUNCTION libtetrabz_polstat_1233(g1,g2,g3,lng1,lng2,lng3) RESULT(w)
561  !
562  IMPLICIT NONE
563  !
564  REAL(8),INTENT(IN) :: g1,g2,g3,lng1,lng2,lng3
565  REAL(8) :: w
566  !
567  REAL(8) :: w2, w3
568  !
569  w2 = (lng2 - lng1)/(g2 - g1)*g2 - 1d0
570  w2 = (g2*w2)/(g2 - g1)
571  w3 = (lng3 - lng1)/(g3 - g1)*g3 - 1d0
572  w3 = (g3*w3)/(g3 - g1)
573  w2 = (w3 - w2)/(g3 - g2)
574  w3 = (lng3 - lng1)/(g3 - g1)*g3 - 1d0
575  w3 = 1d0 - (2d0*w3*g1)/(g3 - g1)
576  w3 = w3/(g3 - g1)
577  w = (g3*w3 - g2*w2)/(g3 - g2)
578  !
579 END FUNCTION libtetrabz_polstat_1233
580 !
581 ! 4, g4 = g1 and g3 = g2
582 !
583 FUNCTION libtetrabz_polstat_1221(g1,g2,lng1,lng2) RESULT(w)
584  !
585  IMPLICIT NONE
586  !
587  REAL(8),INTENT(IN) :: g1, g2, lng1, lng2
588  REAL(8) :: w
589  !
590  w = 1d0 - (lng2 - lng1)/(g2 - g1)*g1
591  w = -1d0 + (2d0*g2*w)/(g2 - g1)
592  w = -1d0 + (3d0*g2*w)/(g2 - g1)
593  w = w/(2d0*(g2 - g1))
594  !
595 END FUNCTION libtetrabz_polstat_1221
596 !
597 ! 5, g4 = g3 = g2
598 !
599 FUNCTION libtetrabz_polstat_1222(g1,g2,lng1,lng2) RESULT(w)
600  !
601  IMPLICIT NONE
602  !
603  REAL(8),INTENT(IN) :: g1, g2, lng1, lng2
604  REAL(8) :: w
605  !
606  w = (lng2 - lng1)/(g2 - g1)*g2 - 1d0
607  w = (2d0*g1*w)/(g2 - g1) - 1d0
608  w = (3d0*g1*w)/(g2 - g1) + 1d0
609  w = w/(2d0*(g2 - g1))
610  !
611 END FUNCTION libtetrabz_polstat_1222
612 !
613 ! 6, g4 = g3 = g1
614 !
615 FUNCTION libtetrabz_polstat_1211(g1,g2,lng1,lng2) RESULT(w)
616  !
617  IMPLICIT NONE
618  !
619  REAL(8),INTENT(IN) :: g1,g2,lng1,lng2
620  REAL(8) :: w
621  !
622  w = -1d0 + (lng2 - lng1)/(g2 - g1)*g2
623  w = -1d0 + (2d0*g2*w)/(g2 - g1)
624  w = -1d0 + (3d0*g2*w)/(2d0*(g2 - g1))
625  w = w/(3d0*(g2 - g1))
626  !
627 END FUNCTION libtetrabz_polstat_1211
628 !
629 END MODULE libtetrabz_polstat_mod
libtetrabz_polstat_mod::libtetrabz_polstat
subroutine, public libtetrabz_polstat(ltetra, bvec, nb, nge, eig1, eig2, ngw, wght, comm)
Definition: libtetrabz_polstat_mod.F90:35
libtetrabz_common::libtetrabz_tsmall_b2
subroutine, public libtetrabz_tsmall_b2(e, V, tsmall)
Definition: libtetrabz_common.F90:432
libtetrabz_common::libtetrabz_tsmall_c1
subroutine, public libtetrabz_tsmall_c1(e, V, tsmall)
Definition: libtetrabz_common.F90:484
libtetrabz_common::libtetrabz_mpisum_dv
subroutine, public libtetrabz_mpisum_dv(comm, ndim, vector)
Definition: libtetrabz_common.F90:687
libtetrabz_common::libtetrabz_tsmall_a1
subroutine, public libtetrabz_tsmall_a1(e, V, tsmall)
Definition: libtetrabz_common.F90:380
libtetrabz_common::libtetrabz_tsmall_b1
subroutine, public libtetrabz_tsmall_b1(e, V, tsmall)
Definition: libtetrabz_common.F90:406
libtetrabz_common::libtetrabz_initialize
subroutine, public libtetrabz_initialize(ltetra, nge, ngw, bvec, linterpol, wlsm, nk_local, nt_local, nkBZ, ik_global, ik_local, kvec, comm)
Definition: libtetrabz_common.F90:40
libtetrabz_polstat_mod
Definition: libtetrabz_polstat_mod.F90:23
libtetrabz_polstat_mod::libtetrabz_polstat_main
subroutine libtetrabz_polstat_main(wlsm, nt_local, ik_global, ik_local, nb, nkBZ, eig1, eig2, nk_local, polstat)
Definition: libtetrabz_polstat_mod.F90:87
libtetrabz_common::libtetrabz_interpol_indx
subroutine, public libtetrabz_interpol_indx(ng, kvec, kintp, wintp)
Definition: libtetrabz_common.F90:340
libtetrabz_common::libtetrabz_tsmall_c3
subroutine, public libtetrabz_tsmall_c3(e, V, tsmall)
Definition: libtetrabz_common.F90:536
libtetrabz_common::libtetrabz_tsmall_c2
subroutine, public libtetrabz_tsmall_c2(e, V, tsmall)
Definition: libtetrabz_common.F90:510
libtetrabz_common
Definition: libtetrabz_common.F90:23
libtetrabz_polstat_mod::libtetrabz_polstat3
subroutine libtetrabz_polstat3(de, w1)
Definition: libtetrabz_polstat_mod.F90:362
libtetrabz_common::libtetrabz_tsmall_b3
subroutine, public libtetrabz_tsmall_b3(e, V, tsmall)
Definition: libtetrabz_common.F90:458
libtetrabz_polstat_mod::libtetrabz_polstat2
subroutine libtetrabz_polstat2(nb, ei1, ej1, w1)
Definition: libtetrabz_polstat_mod.F90:243
libtetrabz_common::libtetrabz_sort
subroutine, public libtetrabz_sort(n, key, indx)
Definition: libtetrabz_common.F90:308
libtetrabz_polstat_mod::libtetrabz_polstat_1234
real(8) function libtetrabz_polstat_1234(g1, g2, g3, g4, lng1, lng2, lng3, lng4)
Definition: libtetrabz_polstat_mod.F90:520