pwdft  0.1
PW-DFT code for education
libtetrabz_common.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_common
24  !
25  IMPLICIT NONE
26  !
27  PRIVATE
28  PUBLIC libtetrabz_initialize, libtetrabz_sort, libtetrabz_interpol_indx, &
29  & libtetrabz_tsmall_a1, libtetrabz_tsmall_b1, libtetrabz_tsmall_b2, libtetrabz_tsmall_b3, &
30  & libtetrabz_tsmall_c1, libtetrabz_tsmall_c2, libtetrabz_tsmall_c3, &
31  & libtetrabz_triangle_a1, libtetrabz_triangle_b1, &
32  & libtetrabz_triangle_b2, libtetrabz_triangle_c1, &
33  & libtetrabz_mpisum_d, libtetrabz_mpisum_dv, libtetrabz_mpisum_zv
34  !
35 CONTAINS
36 !
37 ! define shortest diagonal line & define type of tetragonal
38 !
39 SUBROUTINE libtetrabz_initialize(ltetra,nge,ngw,bvec,linterpol,wlsm,nk_local,nt_local,nkBZ,ik_global,ik_local,kvec,comm)
40  !
41  IMPLICIT NONE
42  !
43  INTEGER,INTENT(IN) :: ltetra, nge(3), ngw(3)
44  REAL(8),INTENT(IN) :: bvec(3,3)
45  LOGICAL,INTENT(OUT) :: linterpol
46  REAL(8),INTENT(OUT) :: wlsm(4,20)
47  INTEGER,INTENT(OUT) :: nk_local, nt_local, nkbz
48  INTEGER,INTENT(OUT),ALLOCATABLE :: ik_global(:,:), ik_local(:,:)
49  REAL(8),INTENT(OUT),ALLOCATABLE :: kvec(:,:)
50  INTEGER,INTENT(IN),OPTIONAL :: comm
51  !
52  INTEGER :: itype, i1, i2, i3, it, divvec(4,4), ivvec0(4), ivvec(3,20,6)
53  REAL(8) :: l(4), bvec2(3,3), bvec3(3,4)
54  !
55  nkbz = product(nge(1:3))
56  linterpol = .NOT. all(nge(1:3) == ngw(1:3))
57  !
58  DO i1 = 1, 3
59  bvec2(1:3,i1) = bvec(1:3,i1) / dble(nge(i1))
60  END DO
61  !
62  bvec3(1:3,1) = -bvec2(1:3,1) + bvec2(1:3,2) + bvec2(1:3,3)
63  bvec3(1:3,2) = bvec2(1:3,1) - bvec2(1:3,2) + bvec2(1:3,3)
64  bvec3(1:3,3) = bvec2(1:3,1) + bvec2(1:3,2) - bvec2(1:3,3)
65  bvec3(1:3,4) = bvec2(1:3,1) + bvec2(1:3,2) + bvec2(1:3,3)
66  !
67  ! length of delta bvec
68  !
69  DO i1 = 1, 4
70  l(i1) = dot_product(bvec3(1:3,i1),bvec3(1:3,i1))
71  END DO
72  !
73  itype = minloc(l(1:4),1)
74  !
75  ! start & last
76  !
77  ivvec0(1:4) = (/ 0, 0, 0, 0 /)
78  !
79  divvec(1:4,1) = (/ 1, 0, 0, 0 /)
80  divvec(1:4,2) = (/ 0, 1, 0, 0 /)
81  divvec(1:4,3) = (/ 0, 0, 1, 0 /)
82  divvec(1:4,4) = (/ 0, 0, 0, 1 /)
83  !
84  ivvec0(itype) = 1
85  divvec(itype, itype) = - 1
86  !
87  ! Corners of tetrahedra
88  !
89  it = 0
90  DO i1 = 1, 3
91  DO i2 = 1, 3
92  IF(i2 == i1) cycle
93  DO i3 = 1, 3
94  IF(i3 == i1 .OR. i3 == i2) cycle
95  !
96  it = it + 1
97  !
98  ivvec(1:3,1,it) = ivvec0(1:3)
99  ivvec(1:3,2,it) = ivvec(1:3,1,it) + divvec(1:3,i1)
100  ivvec(1:3,3,it) = ivvec(1:3,2,it) + divvec(1:3,i2)
101  ivvec(1:3,4,it) = ivvec(1:3,3,it) + divvec(1:3,i3)
102  !
103  END DO
104  END DO
105  END DO
106  !
107  ! Additional points
108  !
109  ivvec(1:3, 5,1:6) = 2 * ivvec(1:3,1,1:6) - ivvec(1:3,2,1:6)
110  ivvec(1:3, 6,1:6) = 2 * ivvec(1:3,2,1:6) - ivvec(1:3,3,1:6)
111  ivvec(1:3, 7,1:6) = 2 * ivvec(1:3,3,1:6) - ivvec(1:3,4,1:6)
112  ivvec(1:3, 8,1:6) = 2 * ivvec(1:3,4,1:6) - ivvec(1:3,1,1:6)
113  !
114  ivvec(1:3, 9,1:6) = 2 * ivvec(1:3,1,1:6) - ivvec(1:3,3,1:6)
115  ivvec(1:3,10,1:6) = 2 * ivvec(1:3,2,1:6) - ivvec(1:3,4,1:6)
116  ivvec(1:3,11,1:6) = 2 * ivvec(1:3,3,1:6) - ivvec(1:3,1,1:6)
117  ivvec(1:3,12,1:6) = 2 * ivvec(1:3,4,1:6) - ivvec(1:3,2,1:6)
118  !
119  ivvec(1:3,13,1:6) = 2 * ivvec(1:3,1,1:6) - ivvec(1:3,4,1:6)
120  ivvec(1:3,14,1:6) = 2 * ivvec(1:3,2,1:6) - ivvec(1:3,1,1:6)
121  ivvec(1:3,15,1:6) = 2 * ivvec(1:3,3,1:6) - ivvec(1:3,2,1:6)
122  ivvec(1:3,16,1:6) = 2 * ivvec(1:3,4,1:6) - ivvec(1:3,3,1:6)
123  !
124  ivvec(1:3,17,1:6) = ivvec(1:3,4,1:6) - ivvec(1:3,1,1:6) + ivvec(1:3,2,1:6)
125  ivvec(1:3,18,1:6) = ivvec(1:3,1,1:6) - ivvec(1:3,2,1:6) + ivvec(1:3,3,1:6)
126  ivvec(1:3,19,1:6) = ivvec(1:3,2,1:6) - ivvec(1:3,3,1:6) + ivvec(1:3,4,1:6)
127  ivvec(1:3,20,1:6) = ivvec(1:3,3,1:6) - ivvec(1:3,4,1:6) + ivvec(1:3,1,1:6)
128  !
129  IF(ltetra == 1) THEN
130  !
131  !WRITE(*,*) "[libtetrabz] Linear tetrahedron method is used."
132  !
133  wlsm(1:4,1:20) = 0.0d0
134  wlsm(1,1) = 1.0d0
135  wlsm(2,2) = 1.0d0
136  wlsm(3,3) = 1.0d0
137  wlsm(4,4) = 1.0d0
138  !
139  ELSE IF(ltetra == 2) THEN
140  !
141  !WRITE(*,*) "[libtetrabz] Improved tetrahedron method is used."
142  !
143  wlsm(1, 1: 4) = dble((/1440, 0, 30, 0/))
144  wlsm(2, 1: 4) = dble((/ 0, 1440, 0, 30/))
145  wlsm(3, 1: 4) = dble((/ 30, 0, 1440, 0/))
146  wlsm(4, 1: 4) = dble((/ 0, 30, 0, 1440/))
147  !
148  wlsm(1, 5: 8) = dble((/ -38, 7, 17, -28/))
149  wlsm(2, 5: 8) = dble((/ -28, -38, 7, 17/))
150  wlsm(3, 5: 8) = dble((/ 17, -28, -38, 7/))
151  wlsm(4, 5: 8) = dble((/ 7, 17, -28, -38/))
152  !
153  wlsm(1, 9:12) = dble((/ -56, 9, -46, 9/))
154  wlsm(2, 9:12) = dble((/ 9, -56, 9, -46/))
155  wlsm(3, 9:12) = dble((/ -46, 9, -56, 9/))
156  wlsm(4, 9:12) = dble((/ 9, -46, 9, -56/))
157  !
158  wlsm(1,13:16) = dble((/ -38, -28, 17, 7/))
159  wlsm(2,13:16) = dble((/ 7, -38, -28, 17/))
160  wlsm(3,13:16) = dble((/ 17, 7, -38, -28/))
161  wlsm(4,13:16) = dble((/ -28, 17, 7, -38/))
162  !
163  wlsm(1,17:20) = dble((/ -18, -18, 12, -18/))
164  wlsm(2,17:20) = dble((/ -18, -18, -18, 12/))
165  wlsm(3,17:20) = dble((/ 12, -18, -18, -18/))
166  wlsm(4,17:20) = dble((/ -18, 12, -18, -18/))
167  !
168  wlsm(1:4,1:20) = wlsm(1:4,1:20) / 1260d0
169  !
170  ELSE
171  !
172  WRITE(*,*) "[libtetrabz] STOP! ltetrta is invalid."
173  stop
174  !
175  END IF
176  !
177  IF (PRESENT(comm)) THEN
178  CALL libtetrabz_kgrid(linterpol,ivvec,nge,nkbz,nk_local,nt_local,ik_global,ik_local,kvec,comm)
179  ELSE
180  CALL libtetrabz_kgrid(linterpol,ivvec,nge,nkbz,nk_local,nt_local,ik_global,ik_local,kvec)
181  END IF
182  !
183 END SUBROUTINE libtetrabz_initialize
184 !
185 ! Initialize grid
186 !
187 SUBROUTINE libtetrabz_kgrid(linterpol,ivvec,ng,nkBZ,nk_local,nt_local,ik_global,ik_local,kvec,comm)
188  !
189  IMPLICIT NONE
190  !
191  LOGICAL,INTENT(INOUT) :: linterpol
192  INTEGER,INTENT(IN) :: ivvec(3,20,6), ng(3), nkbz
193  INTEGER,INTENT(OUT) :: nk_local, nt_local
194  INTEGER,INTENT(OUT),ALLOCATABLE :: ik_global(:,:), ik_local(:,:)
195  REAL(8),INTENT(OUT),ALLOCATABLE :: kvec(:,:)
196  INTEGER,INTENT(IN),OPTIONAL :: comm
197  !
198  INTEGER :: it, i1, i2, i3, ii, ikv(3), nt, ik, nt_front, loc2glob(nkbz), glob2loc(nkbz)
199  !
200  IF(PRESENT(comm)) THEN
201  CALL libtetrabz_dividempi(comm,6 * nkbz,nt_front,nt_local)
202  linterpol = linterpol .OR. (6*nkbz /= nt_local)
203  ELSE
204  nt_front = 0
205  nt_local = 6 * nkbz
206  END IF
207  ALLOCATE(ik_global(20, nt_local), ik_local(20, nt_local))
208  !
209  ! k-index for energy (Global index)
210  !
211  nt = 0
212  DO i3 = 1, ng(3)
213  DO i2 = 1, ng(2)
214  DO i1 = 1, ng(1)
215  !
216  DO it = 1, 6
217  !
218  nt = nt + 1
219  IF(nt <= nt_front .OR. nt_front + nt_local < nt) cycle
220  !
221  DO ii = 1, 20
222  !
223  ikv(1:3) = (/i1, i2, i3/) + ivvec(1:3,ii,it) - 1
224  ikv(1:3) = modulo(ikv(1:3), ng(1:3))
225  !
226  ik_global(ii,nt - nt_front) = 1 + ikv(1) + ng(1) * ikv(2) + ng(1) * ng(2) * ikv(3)
227  !
228  END DO
229  !
230  END DO
231  !
232  END DO
233  END DO
234  END DO
235  !
236  ! k-index for weight (Local index)
237  !
238  IF(.NOT. linterpol) THEN
239  nk_local = nkbz
240  ik_local(1:20,1:nt_local) = ik_global(1:20,1:nt_local)
241  RETURN
242  END IF
243  !
244  glob2loc(1:nkbz) = 0
245  nk_local = 0
246  DO nt = 1, nt_local
247  DO ii = 1, 20
248  !
249  IF(glob2loc(ik_global(ii,nt)) /= 0) THEN
250  ik_local(ii,nt) = glob2loc(ik_global(ii,nt))
251  ELSE
252  !
253  nk_local = nk_local + 1
254  loc2glob(nk_local) = ik_global(ii,nt)
255  glob2loc(ik_global(ii,nt)) = nk_local
256  ik_local(ii,nt) = nk_local
257  !
258  END IF
259  !
260  END DO
261  END DO
262  !
263  ! k-vector in the fractional coordinate
264  !
265  ALLOCATE(kvec(3,nk_local))
266  DO ik = 1, nk_local
267  ! loc2glob(ik) - 1 = i1 + ng(1) * i2 + ng(1) * ng(2) * i3
268  i1 = mod(loc2glob(ik) - 1, ng(1))
269  i2 = mod((loc2glob(ik) - 1) / ng(1), ng(2))
270  i3 = (loc2glob(ik) - 1) / (ng(1) * ng(2))
271  kvec(1:3,ik) = dble((/i1, i2, i3/)) / dble(ng(1:3))
272  END DO
273  !
274 END SUBROUTINE libtetrabz_kgrid
275 !
276 ! Compute cnt and dsp
277 !
278 SUBROUTINE libtetrabz_dividempi(comm,nt,nt_front,nt_local)
279  !
280 #if defined(__MPI)
281  USE mpi, ONLY : mpi_comm_size, mpi_comm_rank
282 #endif
283  IMPLICIT NONE
284  !
285  INTEGER,INTENT(IN) :: comm, nt
286  INTEGER,INTENT(OUT) :: nt_front, nt_local
287  !
288  INTEGER :: petot = 1, my_rank = 0
289 #if defined(__MPI)
290  INTEGER :: ierr
291  CALL mpi_comm_size(comm, petot, ierr)
292  CALL mpi_comm_rank(comm, my_rank, ierr)
293 #endif
294  !
295  IF(my_rank < mod(nt, petot)) THEN
296  nt_local = nt / petot + 1
297  nt_front = my_rank * nt_local
298  ELSE
299  nt_local = nt / petot
300  nt_front = my_rank * nt_local + mod(nt, petot)
301  END IF
302  !
303 END SUBROUTINE libtetrabz_dividempi
304 !
305 ! Simple sort
306 !
307 SUBROUTINE libtetrabz_sort(n,key,indx)
308  !
309  IMPLICIT NONE
310  !
311  integer,INTENT(IN) :: n
312  REAL(8),INTENT(inout) :: key(n)
313  INTEGER,INTENT(OUT) :: indx(n)
314  !
315  INTEGER :: i, i0, indx0
316  REAL(8) :: key0
317  !
318  DO i = 1, n
319  indx(i) = i
320  END DO
321  !
322  DO i = 1, n - 1
323  key0 = minval(key(i+1:n))
324  i0 = minloc(key(i+1:n),1) + i
325  IF(key(i) > key0) THEN
326  key(i0) = key(i)
327  key(i) = key0
328  !
329  indx0 = indx(i0)
330  indx(i0) = indx(i)
331  indx(i) = indx0
332  END IF
333  END DO
334  !
335 END SUBROUTINE libtetrabz_sort
336 !
337 ! Linear interpolation
338 !
339 SUBROUTINE libtetrabz_interpol_indx(ng,kvec,kintp,wintp)
340  !
341  IMPLICIT NONE
342  !
343  INTEGER,INTENT(in) :: ng(3)
344  REAL(8),INTENT(in) :: kvec(3)
345  INTEGER,INTENT(out) :: kintp(8)
346  REAL(8),INTENT(out) :: wintp(8)
347  !
348  INTEGER :: ii
349  REAL(8) :: xv(3)
350  integer :: ikv0(3), ikv1(3), i1, i2, i3
351  !
352  xv(1:3) = kvec(1:3) * dble(ng(1:3))
353  ikv0(1:3) = floor(xv(1:3))
354  xv(1:3) = xv(1:3) - dble(ikv0(1:3))
355  !
356  ii = 0
357  DO i1 = 0, 1
358  DO i2 = 0, 1
359  DO i3 = 0, 1
360  !
361  ii = ii + 1
362  !
363  ikv1(1:3) = ikv0(1:3) + (/i1, i2, i3/)
364  ikv1(1:3) = modulo(ikv1(1:3), ng(1:3))
365  kintp(ii) = 1 + ikv1(1) + ng(1) * ikv1(2) + ng(1) * ng(2) * ikv1(3)
366  !
367  wintp(ii) = xv(1)**i1 * (1.0d0 - xv(1))**(1 - i1) &
368  & * xv(2)**i2 * (1.0d0 - xv(2))**(1 - i2) &
369  & * xv(3)**i3 * (1.0d0 - xv(3))**(1 - i3)
370  !
371  END DO
372  END DO
373  END DO
374  !
375 END SUBROUTINE libtetrabz_interpol_indx
376 !
377 ! Cut small tetrahedron A1
378 !
379 SUBROUTINE libtetrabz_tsmall_a1(e,V,tsmall)
380  !
381  IMPLICIT NONE
382  !
383  REAL(8),INTENT(IN) :: e(4)
384  REAL(8),INTENT(OUT) :: v
385  REAL(8),INTENT(OUT) :: tsmall(4,4)
386  !
387  INTEGER :: ii
388  REAL(8) :: a(4,4)
389  !
390  DO ii = 1, 4
391  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
392  END DO
393  !
394  v = a(2,1) * a(3,1) * a(4,1)
395  !
396  tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
397  tsmall(2, 1:4) = (/a(1,2), a(2,1), 0d0, 0d0/)
398  tsmall(3, 1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
399  tsmall(4, 1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
400  !
401 END SUBROUTINE libtetrabz_tsmall_a1
402 !
403 ! Cut small tetrahedron B1
404 !
405 SUBROUTINE libtetrabz_tsmall_b1(e,V,tsmall)
406  !
407  IMPLICIT NONE
408  !
409  REAL(8),INTENT(IN) :: e(4)
410  REAL(8),INTENT(OUT) :: v
411  REAL(8),INTENT(OUT) :: tsmall(4,4)
412  !
413  INTEGER :: ii
414  REAL(8) :: a(4,4)
415  !
416  DO ii = 1, 4
417  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
418  END DO
419  !
420  v = a(3,1) * a(4,1) * a(2,4)
421  !
422  tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
423  tsmall(2, 1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
424  tsmall(3, 1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
425  tsmall(4, 1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
426  !
427 END SUBROUTINE libtetrabz_tsmall_b1
428 !
429 ! Cut small tetrahedron B2
430 !
431 SUBROUTINE libtetrabz_tsmall_b2(e,V,tsmall)
432  !
433  IMPLICIT NONE
434  !
435  REAL(8),INTENT(IN) :: e(4)
436  REAL(8),INTENT(OUT) :: v
437  REAL(8),INTENT(OUT) :: tsmall(4,4)
438  !
439  INTEGER :: ii
440  REAL(8) :: a(4,4)
441  !
442  DO ii = 1, 4
443  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
444  END DO
445  !
446  v = a(3,2) * a(4,2)
447  !
448  tsmall(1, 1:4) = (/1d0, 0d0, 0d0, 0d0/)
449  tsmall(2, 1:4) = (/0d0, 1d0, 0d0, 0d0/)
450  tsmall(3, 1:4) = (/0d0, a(2,3), a(3,2), 0d0/)
451  tsmall(4, 1:4) = (/0d0, a(2,4), 0d0, a(4,2)/)
452  !
453 END SUBROUTINE libtetrabz_tsmall_b2
454 !
455 ! Cut small tetrahedron B3
456 !
457 SUBROUTINE libtetrabz_tsmall_b3(e,V,tsmall)
458  !
459  IMPLICIT NONE
460  !
461  REAL(8),INTENT(IN) :: e(4)
462  REAL(8),INTENT(OUT) :: v
463  REAL(8),INTENT(OUT) :: tsmall(4,4)
464  !
465  INTEGER :: ii
466  REAL(8) :: a(4,4)
467  !
468  DO ii = 1, 4
469  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
470  END DO
471  !
472  v = a(2,3) * a(3,1) * a(4,2)
473  !
474  tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
475  tsmall(2, 1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
476  tsmall(3, 1:4) = (/ 0d0, a(2,3), a(3,2), 0d0/)
477  tsmall(4, 1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
478  !
479 END SUBROUTINE libtetrabz_tsmall_b3
480 !
481 ! Cut small tetrahedron C1
482 !
483 SUBROUTINE libtetrabz_tsmall_c1(e,V,tsmall)
484  !
485  IMPLICIT NONE
486  !
487  REAL(8),INTENT(IN) :: e(4)
488  REAL(8),INTENT(OUT) :: v
489  REAL(8),INTENT(OUT) :: tsmall(4,4)
490  !
491  INTEGER :: ii
492  REAL(8) :: a(4,4)
493  !
494  DO ii = 1, 4
495  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
496  END DO
497  !
498  v = a(4,3)
499  !
500  tsmall(1, 1:4) = (/1d0, 0d0, 0d0, 0d0/)
501  tsmall(2, 1:4) = (/0d0, 1d0, 0d0, 0d0/)
502  tsmall(3, 1:4) = (/0d0, 0d0, 1d0, 0d0/)
503  tsmall(4, 1:4) = (/0d0, 0d0, a(3,4), a(4,3)/)
504  !
505 END SUBROUTINE libtetrabz_tsmall_c1
506 !
507 ! Cut small tetrahedron C2
508 !
509 SUBROUTINE libtetrabz_tsmall_c2(e,V,tsmall)
510  !
511  IMPLICIT NONE
512  !
513  REAL(8),INTENT(IN) :: e(4)
514  REAL(8),INTENT(OUT) :: v
515  REAL(8),INTENT(OUT) :: tsmall(4,4)
516  !
517  INTEGER :: ii
518  REAL(8) :: a(4,4)
519  !
520  DO ii = 1, 4
521  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
522  END DO
523  !
524  v = a(3,4) * a(4,2)
525  !
526  tsmall(1, 1:4) = (/1d0, 0d0, 0d0, 0d0/)
527  tsmall(2, 1:4) = (/0d0, 1d0, 0d0, 0d0/)
528  tsmall(3, 1:4) = (/0d0, a(2,4), 0d0, a(4,2)/)
529  tsmall(4, 1:4) = (/0d0, 0d0, a(3,4), a(4,3)/)
530  !
531 END SUBROUTINE libtetrabz_tsmall_c2
532 !
533 ! Cut small tetrahedron C3
534 !
535 SUBROUTINE libtetrabz_tsmall_c3(e,V,tsmall)
536  !
537  IMPLICIT NONE
538  !
539  REAL(8),INTENT(IN) :: e(4)
540  REAL(8),INTENT(OUT) :: v
541  REAL(8),INTENT(OUT) :: tsmall(4,4)
542  !
543  INTEGER :: ii
544  REAL(8) :: a(4,4)
545  !
546  DO ii = 1, 4
547  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
548  END DO
549  !
550  v = a(3,4) * a(2,4) * a(4,1)
551  !
552  tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
553  tsmall(2, 1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
554  tsmall(3, 1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
555  tsmall(4, 1:4) = (/ 0d0, 0d0, a(3,4), a(4,3)/)
556  !
557 END SUBROUTINE libtetrabz_tsmall_c3
558 !
559 ! Cut triangle A1
560 !
561 SUBROUTINE libtetrabz_triangle_a1(e,V,tsmall)
562  !
563  IMPLICIT NONE
564  !
565  REAL(8),INTENT(IN) :: e(4)
566  REAL(8),INTENT(OUT) :: v
567  REAL(8),INTENT(OUT) :: tsmall(3,4)
568  !
569  INTEGER :: ii
570  REAL(8) :: a(4,4)
571  !
572  DO ii = 1, 4
573  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
574  END DO
575  !
576  !V = 3d0 * a(2,1) * a(3,1) * a(4,1) / (0d0 - e(1))
577  v = 3d0 * a(2,1) * a(3,1) / (e(4) - e(1))
578  !
579  tsmall(1,1:4) = (/a(1,2), a(2,1), 0d0, 0d0/)
580  tsmall(2,1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
581  tsmall(3,1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
582  !
583 END SUBROUTINE libtetrabz_triangle_a1
584 !
585 ! Cut triangle B1
586 !
587 SUBROUTINE libtetrabz_triangle_b1(e,V,tsmall)
588  !
589  IMPLICIT NONE
590  !
591  REAL(8),INTENT(IN) :: e(4)
592  REAL(8),INTENT(OUT) :: v
593  REAL(8),INTENT(OUT) :: tsmall(3,4)
594  !
595  INTEGER :: ii
596  REAL(8) :: a(4,4)
597  !
598  DO ii = 1, 4
599  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
600  END DO
601  !
602  !V = 3d0 * a(3,1) * a(4,1) * a(2,4) / (0d0 - e(1))
603  v = 3d0 * a(4,1) * a(2,4) / (e(3) - e(1))
604  !
605  tsmall(1,1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
606  tsmall(2,1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
607  tsmall(3,1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
608  !
609 END SUBROUTINE libtetrabz_triangle_b1
610 !
611 ! Cut triangle B2
612 !
613 SUBROUTINE libtetrabz_triangle_b2(e,V,tsmall)
614  !
615  IMPLICIT NONE
616  !
617  REAL(8),INTENT(IN) :: e(4)
618  REAL(8),INTENT(OUT) :: v
619  REAL(8),INTENT(OUT) :: tsmall(3,4)
620  !
621  INTEGER :: ii
622  REAL(8) :: a(4,4)
623  !
624  DO ii = 1, 4
625  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
626  END DO
627  !
628  !V = 3d0 * a(2,3) * a(3,1) * a(4,2) / (0d0 - e(1))
629  v = 3d0 * a(2,3) * a(4,2) / (e(3) - e(1))
630  !
631  tsmall(1,1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
632  tsmall(2,1:4) = (/ 0d0, a(2,3), a(3,2), 0d0/)
633  tsmall(3,1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
634  !
635 END SUBROUTINE libtetrabz_triangle_b2
636 !
637 ! Cut triangle C1
638 !
639 SUBROUTINE libtetrabz_triangle_c1(e,V,tsmall)
640  !
641  IMPLICIT NONE
642  !
643  REAL(8),INTENT(IN) :: e(4)
644  REAL(8),INTENT(OUT) :: v
645  REAL(8),INTENT(OUT) :: tsmall(3,4)
646  !
647  INTEGER :: ii
648  REAL(8) :: a(4,4)
649  !
650  DO ii = 1, 4
651  a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
652  END DO
653  !
654  !V = 3d0 * a(1,4) * a(2,4) * a(3,4) / (e(4) - 0d0)
655  v = 3d0 * a(1,4) * a(2,4) / (e(4) - e(3))
656  !
657  tsmall(1,1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
658  tsmall(2,1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
659  tsmall(3,1:4) = (/ 0d0, 0d0, a(3,4), a(4,3)/)
660  !
661 END SUBROUTINE libtetrabz_triangle_c1
662 !
663 ! MPI_Allreduce for double scaler
664 !
665 SUBROUTINE libtetrabz_mpisum_d(comm,scaler)
666  !
667 #if defined(__MPI)
668  USE mpi, ONLY : mpi_double_precision, mpi_in_place, mpi_sum
669 #endif
670  IMPLICIT NONE
671  !
672  INTEGER :: comm
673  REAL(8) :: scaler
674  !
675 #if defined(__MPI)
676  INTEGER :: ierr
677  !
678  CALL mpi_allreduce(mpi_in_place, scaler, 1, &
679  & mpi_double_precision, mpi_sum, comm, ierr)
680 #endif
681  !
682 END SUBROUTINE libtetrabz_mpisum_d
683 !
684 ! MPI_Allreduce for double vector
685 !
686 SUBROUTINE libtetrabz_mpisum_dv(comm,ndim,vector)
687  !
688 #if defined(__MPI)
689  USE mpi, ONLY : mpi_double_precision, mpi_in_place, mpi_sum
690 #endif
691  IMPLICIT NONE
692  !
693  INTEGER :: comm, ndim
694  REAL(8) :: vector(ndim)
695  !
696 #if defined(__MPI)
697  INTEGER :: ierr
698  !
699  CALL mpi_allreduce(mpi_in_place, vector, ndim, &
700  & mpi_double_precision, mpi_sum, comm, ierr)
701 #endif
702  !
703 END SUBROUTINE libtetrabz_mpisum_dv
704 !
705 ! MPI_Allreduce for double complex vector
706 !
707 SUBROUTINE libtetrabz_mpisum_zv(comm,ndim,vector)
708  !
709 #if defined(__MPI)
710  USE mpi, ONLY : mpi_double_complex, mpi_in_place, mpi_sum
711 #endif
712  IMPLICIT NONE
713  !
714  INTEGER :: comm, ndim
715  COMPLEX(8) :: vector(ndim)
716  !
717 #if defined(__MPI)
718  INTEGER :: ierr
719  !
720  CALL mpi_allreduce(mpi_in_place, vector, ndim, &
721  & mpi_double_complex, mpi_sum, comm, ierr)
722 #endif
723  !
724 END SUBROUTINE libtetrabz_mpisum_zv
725 !
726 END MODULE libtetrabz_common
libtetrabz_common::libtetrabz_triangle_b2
subroutine, public libtetrabz_triangle_b2(e, V, tsmall)
Definition: libtetrabz_common.F90:614
libtetrabz_common::libtetrabz_tsmall_b2
subroutine, public libtetrabz_tsmall_b2(e, V, tsmall)
Definition: libtetrabz_common.F90:432
libtetrabz_common::libtetrabz_mpisum_d
subroutine, public libtetrabz_mpisum_d(comm, scaler)
Definition: libtetrabz_common.F90:666
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_triangle_c1
subroutine, public libtetrabz_triangle_c1(e, V, tsmall)
Definition: libtetrabz_common.F90:640
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_common::libtetrabz_interpol_indx
subroutine, public libtetrabz_interpol_indx(ng, kvec, kintp, wintp)
Definition: libtetrabz_common.F90:340
libtetrabz_common::libtetrabz_mpisum_zv
subroutine, public libtetrabz_mpisum_zv(comm, ndim, vector)
Definition: libtetrabz_common.F90:708
libtetrabz_common::libtetrabz_triangle_a1
subroutine, public libtetrabz_triangle_a1(e, V, tsmall)
Definition: libtetrabz_common.F90:562
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_common::libtetrabz_tsmall_b3
subroutine, public libtetrabz_tsmall_b3(e, V, tsmall)
Definition: libtetrabz_common.F90:458
libtetrabz_common::libtetrabz_sort
subroutine, public libtetrabz_sort(n, key, indx)
Definition: libtetrabz_common.F90:308
libtetrabz_common::libtetrabz_triangle_b1
subroutine, public libtetrabz_triangle_b1(e, V, tsmall)
Definition: libtetrabz_common.F90:588