Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : #include "abi_common.h"
7 : !
8 : ! Copyright (C) 2014 Mitsuaki Kawamura
9 : !
10 : ! Permission is hereby granted, free of charge, to any person obtaining a
11 : ! copy of this software and associated documentation files (the
12 : ! "Software"), to deal in the Software without restriction, including
13 : ! without limitation the rights to use, copy, modify, merge, publish,
14 : ! distribute, sublicense, and/or sell copies of the Software, and to
15 : ! permit persons to whom the Software is furnished to do so, subject to
16 : ! the following conditions:
17 : !
18 : ! The above copyright notice and this permission notice shall be included
19 : ! in all copies or substantial portions of the Software.
20 : !
21 : ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : ! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
23 : ! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
24 : ! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
25 : ! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
26 : ! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
27 : ! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
28 : !
29 : MODULE libtetrabz_common
30 : !
31 : use m_abicore
32 : use m_errors
33 : USE_MPI
34 :
35 : IMPLICIT NONE
36 :
37 : #if defined HAVE_MPI1
38 : include 'mpif.h'
39 : #endif
40 : !
41 : PRIVATE
42 :
43 : PUBLIC :: libtetrabz_initialize, libtetrabz_sort, libtetrabz_interpol_indx, &
44 : & libtetrabz_tsmall_a1, libtetrabz_tsmall_b1, libtetrabz_tsmall_b2, libtetrabz_tsmall_b3, &
45 : & libtetrabz_tsmall_c1, libtetrabz_tsmall_c2, libtetrabz_tsmall_c3, &
46 : & libtetrabz_triangle_a1, libtetrabz_triangle_b1, &
47 : & libtetrabz_triangle_b2, libtetrabz_triangle_c1, &
48 : & libtetrabz_mpisum_d, libtetrabz_mpisum_dv, libtetrabz_mpisum_zv
49 : !
50 :
51 : CONTAINS
52 : !
53 : ! define shortest diagonal line & define type of tetragonal
54 : !
55 35 : SUBROUTINE libtetrabz_initialize(ltetra,nge,ngw,bvec,linterpol,wlsm,nk_local,nt_local,nkBZ,ik_global,ik_local,kvec,comm)
56 : !
57 : !IMPLICIT NONE
58 : !
59 : INTEGER,INTENT(IN) :: ltetra, nge(3), ngw(3)
60 : REAL(8),INTENT(IN) :: bvec(3,3)
61 : LOGICAL,INTENT(OUT) :: linterpol
62 : REAL(8),INTENT(OUT) :: wlsm(4,20)
63 : INTEGER,INTENT(OUT) :: nk_local, nt_local, nkBZ
64 : INTEGER,INTENT(OUT),ALLOCATABLE :: ik_global(:,:), ik_local(:,:)
65 : REAL(8),INTENT(OUT),ALLOCATABLE :: kvec(:,:)
66 : INTEGER,INTENT(IN),OPTIONAL :: comm
67 : !
68 : INTEGER :: itype, i1, i2, i3, it, divvec(4,4), ivvec0(4), ivvec(3,20,6)
69 : REAL(8) :: l(4), bvec2(3,3), bvec3(3,4)
70 : !
71 140 : nkBZ = PRODUCT(nge(1:3))
72 140 : linterpol = .NOT. ALL(nge(1:3) == ngw(1:3))
73 : !
74 140 : DO i1 = 1, 3
75 455 : bvec2(1:3,i1) = bvec(1:3,i1) / DBLE(nge(i1))
76 : END DO
77 : !
78 140 : bvec3(1:3,1) = -bvec2(1:3,1) + bvec2(1:3,2) + bvec2(1:3,3)
79 140 : bvec3(1:3,2) = bvec2(1:3,1) - bvec2(1:3,2) + bvec2(1:3,3)
80 140 : bvec3(1:3,3) = bvec2(1:3,1) + bvec2(1:3,2) - bvec2(1:3,3)
81 140 : bvec3(1:3,4) = bvec2(1:3,1) + bvec2(1:3,2) + bvec2(1:3,3)
82 : !
83 : ! length of delta bvec
84 : !
85 175 : DO i1 = 1, 4
86 595 : l(i1) = DOT_PRODUCT(bvec3(1:3,i1),bvec3(1:3,i1))
87 : END DO
88 : !
89 210 : itype = MINLOC(l(1:4),1)
90 : !
91 : ! start & last
92 : !
93 35 : ivvec0(1:4) = (/ 0, 0, 0, 0 /)
94 : !
95 175 : divvec(1:4,1) = (/ 1, 0, 0, 0 /)
96 175 : divvec(1:4,2) = (/ 0, 1, 0, 0 /)
97 175 : divvec(1:4,3) = (/ 0, 0, 1, 0 /)
98 : divvec(1:4,4) = (/ 0, 0, 0, 1 /)
99 : !
100 35 : ivvec0(itype) = 1
101 35 : divvec(itype, itype) = - 1
102 : !
103 : ! Corners of tetrahedra
104 : !
105 35 : it = 0
106 140 : DO i1 = 1, 3
107 455 : DO i2 = 1, 3
108 315 : IF(i2 == i1) CYCLE
109 945 : DO i3 = 1, 3
110 630 : IF(i3 == i1 .OR. i3 == i2) CYCLE
111 : !
112 210 : it = it + 1
113 : !
114 840 : ivvec(1:3,1,it) = ivvec0(1:3)
115 840 : ivvec(1:3,2,it) = ivvec(1:3,1,it) + divvec(1:3,i1)
116 840 : ivvec(1:3,3,it) = ivvec(1:3,2,it) + divvec(1:3,i2)
117 1155 : ivvec(1:3,4,it) = ivvec(1:3,3,it) + divvec(1:3,i3)
118 : !
119 : END DO
120 : END DO
121 : END DO
122 : !
123 : ! Additional points
124 : !
125 875 : ivvec(1:3, 5,1:6) = 2 * ivvec(1:3,1,1:6) - ivvec(1:3,2,1:6)
126 875 : ivvec(1:3, 6,1:6) = 2 * ivvec(1:3,2,1:6) - ivvec(1:3,3,1:6)
127 875 : ivvec(1:3, 7,1:6) = 2 * ivvec(1:3,3,1:6) - ivvec(1:3,4,1:6)
128 875 : ivvec(1:3, 8,1:6) = 2 * ivvec(1:3,4,1:6) - ivvec(1:3,1,1:6)
129 : !
130 875 : ivvec(1:3, 9,1:6) = 2 * ivvec(1:3,1,1:6) - ivvec(1:3,3,1:6)
131 875 : ivvec(1:3,10,1:6) = 2 * ivvec(1:3,2,1:6) - ivvec(1:3,4,1:6)
132 875 : ivvec(1:3,11,1:6) = 2 * ivvec(1:3,3,1:6) - ivvec(1:3,1,1:6)
133 875 : ivvec(1:3,12,1:6) = 2 * ivvec(1:3,4,1:6) - ivvec(1:3,2,1:6)
134 : !
135 875 : ivvec(1:3,13,1:6) = 2 * ivvec(1:3,1,1:6) - ivvec(1:3,4,1:6)
136 875 : ivvec(1:3,14,1:6) = 2 * ivvec(1:3,2,1:6) - ivvec(1:3,1,1:6)
137 875 : ivvec(1:3,15,1:6) = 2 * ivvec(1:3,3,1:6) - ivvec(1:3,2,1:6)
138 875 : ivvec(1:3,16,1:6) = 2 * ivvec(1:3,4,1:6) - ivvec(1:3,3,1:6)
139 : !
140 875 : ivvec(1:3,17,1:6) = ivvec(1:3,4,1:6) - ivvec(1:3,1,1:6) + ivvec(1:3,2,1:6)
141 875 : ivvec(1:3,18,1:6) = ivvec(1:3,1,1:6) - ivvec(1:3,2,1:6) + ivvec(1:3,3,1:6)
142 875 : ivvec(1:3,19,1:6) = ivvec(1:3,2,1:6) - ivvec(1:3,3,1:6) + ivvec(1:3,4,1:6)
143 875 : ivvec(1:3,20,1:6) = ivvec(1:3,3,1:6) - ivvec(1:3,4,1:6) + ivvec(1:3,1,1:6)
144 : !
145 35 : IF(ltetra == 1) THEN
146 : !
147 : !WRITE(*,*) "[libtetrabz] Linear tetrahedron method is used."
148 : !
149 0 : wlsm(1:4,1:20) = 0.0d0
150 0 : wlsm(1,1) = 1.0d0
151 0 : wlsm(2,2) = 1.0d0
152 0 : wlsm(3,3) = 1.0d0
153 0 : wlsm(4,4) = 1.0d0
154 : !
155 35 : ELSE IF(ltetra == 2) THEN
156 : !
157 : !WRITE(*,*) "[libtetrabz] Improved tetrahedron method is used."
158 : !
159 175 : wlsm(1, 1: 4) = DBLE((/1440, 0, 30, 0/))
160 175 : wlsm(2, 1: 4) = DBLE((/ 0, 1440, 0, 30/))
161 175 : wlsm(3, 1: 4) = DBLE((/ 30, 0, 1440, 0/))
162 175 : wlsm(4, 1: 4) = DBLE((/ 0, 30, 0, 1440/))
163 : !
164 175 : wlsm(1, 5: 8) = DBLE((/ -38, 7, 17, -28/))
165 175 : wlsm(2, 5: 8) = DBLE((/ -28, -38, 7, 17/))
166 175 : wlsm(3, 5: 8) = DBLE((/ 17, -28, -38, 7/))
167 175 : wlsm(4, 5: 8) = DBLE((/ 7, 17, -28, -38/))
168 : !
169 175 : wlsm(1, 9:12) = DBLE((/ -56, 9, -46, 9/))
170 175 : wlsm(2, 9:12) = DBLE((/ 9, -56, 9, -46/))
171 175 : wlsm(3, 9:12) = DBLE((/ -46, 9, -56, 9/))
172 175 : wlsm(4, 9:12) = DBLE((/ 9, -46, 9, -56/))
173 : !
174 175 : wlsm(1,13:16) = DBLE((/ -38, -28, 17, 7/))
175 175 : wlsm(2,13:16) = DBLE((/ 7, -38, -28, 17/))
176 175 : wlsm(3,13:16) = DBLE((/ 17, 7, -38, -28/))
177 175 : wlsm(4,13:16) = DBLE((/ -28, 17, 7, -38/))
178 : !
179 175 : wlsm(1,17:20) = DBLE((/ -18, -18, 12, -18/))
180 175 : wlsm(2,17:20) = DBLE((/ -18, -18, -18, 12/))
181 175 : wlsm(3,17:20) = DBLE((/ 12, -18, -18, -18/))
182 175 : wlsm(4,17:20) = DBLE((/ -18, 12, -18, -18/))
183 : !
184 3535 : wlsm(1:4,1:20) = wlsm(1:4,1:20) / 1260d0
185 : !
186 : ELSE
187 : !
188 0 : ABI_ERROR("[libtetrabz] STOP! ltetrta is invalid.")
189 : !
190 : END IF
191 : !
192 35 : IF (PRESENT(comm)) THEN
193 35 : CALL libtetrabz_kgrid(linterpol,ivvec,nge,nkBZ,nk_local,nt_local,ik_global,ik_local,kvec,comm)
194 : ELSE
195 0 : CALL libtetrabz_kgrid(linterpol,ivvec,nge,nkBZ,nk_local,nt_local,ik_global,ik_local,kvec)
196 : END IF
197 : !
198 35 : END SUBROUTINE libtetrabz_initialize
199 : !
200 : ! Initialize grid
201 : !
202 35 : SUBROUTINE libtetrabz_kgrid(linterpol,ivvec,ng,nkBZ,nk_local,nt_local,ik_global,ik_local,kvec,comm)
203 : !
204 : !IMPLICIT NONE
205 : !
206 : LOGICAL,INTENT(INOUT) :: linterpol
207 : INTEGER,INTENT(IN) :: ivvec(3,20,6), ng(3), nkBZ
208 : INTEGER,INTENT(OUT) :: nk_local, nt_local
209 : INTEGER,INTENT(OUT),ALLOCATABLE :: ik_global(:,:), ik_local(:,:)
210 : REAL(8),INTENT(OUT),ALLOCATABLE :: kvec(:,:)
211 : INTEGER,INTENT(IN),OPTIONAL :: comm
212 : !
213 35 : INTEGER :: it, i1, i2, i3, ii, ikv(3), nt, ik, nt_front, loc2glob(nkBZ), glob2loc(nkBZ)
214 : !
215 35 : IF(PRESENT(comm)) THEN
216 35 : CALL libtetrabz_divideMPI(comm,6 * nkBZ,nt_front,nt_local)
217 35 : linterpol = linterpol .OR. (6*nkBZ /= nt_local)
218 : ELSE
219 0 : nt_front = 0
220 0 : nt_local = 6 * nkBZ
221 : END IF
222 105 : ABI_MALLOC(ik_global, (20, nt_local))
223 70 : ABI_MALLOC(ik_local, (20, nt_local))
224 : !
225 : ! k-index for energy (Global index)
226 : !
227 35 : nt = 0
228 315 : DO i3 = 1, ng(3)
229 2555 : DO i2 = 1, ng(2)
230 20440 : DO i1 = 1, ng(1)
231 : !
232 127680 : DO it = 1, 6
233 : !
234 107520 : nt = nt + 1
235 107520 : IF(nt <= nt_front .OR. nt_front + nt_local < nt) CYCLE
236 : !
237 2275840 : DO ii = 1, 20
238 : !
239 8601600 : ikv(1:3) = (/i1, i2, i3/) + ivvec(1:3,ii,it) - 1
240 8601600 : ikv(1:3) = MODULO(ikv(1:3), ng(1:3))
241 : !
242 2257920 : ik_global(ii,nt - nt_front) = 1 + ikv(1) + ng(1) * ikv(2) + ng(1) * ng(2) * ikv(3)
243 : !
244 : END DO
245 : !
246 : END DO
247 : !
248 : END DO
249 : END DO
250 : END DO
251 : !
252 : ! k-index for weight (Local index)
253 : !
254 35 : IF(.NOT. linterpol) THEN
255 35 : nk_local = nkBZ
256 2257955 : ik_local(1:20,1:nt_local) = ik_global(1:20,1:nt_local)
257 : RETURN
258 : END IF
259 : !
260 0 : glob2loc(1:nkBZ) = 0
261 0 : nk_local = 0
262 0 : DO nt = 1, nt_local
263 0 : DO ii = 1, 20
264 : !
265 0 : IF(glob2loc(ik_global(ii,nt)) /= 0) THEN
266 0 : ik_local(ii,nt) = glob2loc(ik_global(ii,nt))
267 : ELSE
268 : !
269 0 : nk_local = nk_local + 1
270 0 : loc2glob(nk_local) = ik_global(ii,nt)
271 0 : glob2loc(ik_global(ii,nt)) = nk_local
272 0 : ik_local(ii,nt) = nk_local
273 : !
274 : END IF
275 : !
276 : END DO
277 : END DO
278 : !
279 : ! k-vector in the fractional coordinate
280 : !
281 0 : ABI_MALLOC(kvec, (3,nk_local))
282 0 : DO ik = 1, nk_local
283 : ! loc2glob(ik) - 1 = i1 + ng(1) * i2 + ng(1) * ng(2) * i3
284 0 : i1 = MOD(loc2glob(ik) - 1, ng(1))
285 0 : i2 = MOD((loc2glob(ik) - 1) / ng(1), ng(2))
286 0 : i3 = (loc2glob(ik) - 1) / (ng(1) * ng(2))
287 0 : kvec(1:3,ik) = DBLE((/i1, i2, i3/)) / DBLE(ng(1:3))
288 : END DO
289 : !
290 : END SUBROUTINE libtetrabz_kgrid
291 : !
292 : ! Compute cnt and dsp
293 : !
294 35 : SUBROUTINE libtetrabz_divideMPI(comm,nt,nt_front,nt_local)
295 : !
296 : !IMPLICIT NONE
297 : !
298 : INTEGER,INTENT(IN) :: comm, nt
299 : INTEGER,INTENT(OUT) :: nt_front, nt_local
300 : !
301 : INTEGER :: petot = 1, my_rank = 0
302 : #if defined(HAVE_MPI)
303 : INTEGER :: ierr
304 35 : CALL MPI_COMM_SIZE(comm, petot, ierr)
305 35 : CALL MPI_COMM_RANK(comm, my_rank, ierr)
306 : #endif
307 : !
308 35 : IF(my_rank < MOD(nt, petot)) THEN
309 0 : nt_local = nt / petot + 1
310 0 : nt_front = my_rank * nt_local
311 : ELSE
312 35 : nt_local = nt / petot
313 35 : nt_front = my_rank * nt_local + MOD(nt, petot)
314 : END IF
315 : !
316 35 : END SUBROUTINE libtetrabz_divideMPI
317 : !
318 : ! Simple sort
319 : !
320 364560 : pure SUBROUTINE libtetrabz_sort(n,key,indx)
321 : !
322 : !IMPLICIT NONE
323 : !
324 : integer,INTENT(IN) :: n
325 : REAL(8),INTENT(inout) :: key(n)
326 : INTEGER,INTENT(OUT) :: indx(n)
327 : !
328 : INTEGER :: i, i0, indx0
329 : REAL(8) :: key0
330 : !
331 1673280 : DO i = 1, n
332 1673280 : indx(i) = i
333 : END DO
334 : !
335 1308720 : DO i = 1, n - 1
336 3627120 : key0 = MINVAL(key(i+1:n))
337 2682960 : i0 = MINLOC(key(i+1:n),1) + i
338 1308720 : IF(key(i) > key0) THEN
339 552154 : key(i0) = key(i)
340 552154 : key(i) = key0
341 : !
342 552154 : indx0 = indx(i0)
343 552154 : indx(i0) = indx(i)
344 552154 : indx(i) = indx0
345 : END IF
346 : END DO
347 : !
348 364560 : END SUBROUTINE libtetrabz_sort
349 : !
350 : ! Linear interpolation
351 : !
352 0 : pure SUBROUTINE libtetrabz_interpol_indx(nintp,ng,kvec,kintp,wintp)
353 : !
354 : !IMPLICIT NONE
355 : !
356 : INTEGER,INTENT(in) :: nintp, ng(3)
357 : REAL(8),INTENT(in) :: kvec(3)
358 : INTEGER,INTENT(out) :: kintp(nintp)
359 : REAL(8),INTENT(out) :: wintp(nintp)
360 : !
361 : INTEGER :: ikv(3,20), dikv(3,3), ii
362 : REAL(8) :: x, y, z, xv(3)
363 : !
364 : ! Search nearest neighbor grid points.
365 : !
366 0 : xv(1:3) = kvec(1:3) * DBLE(ng(1:3))
367 0 : ikv(1:3,1) = NINT(xv(1:3))
368 0 : dikv(1:3,1:3) = 0
369 0 : DO ii = 1, 3
370 0 : dikv(ii,ii) = ikv(ii,1) - FLOOR(xv(ii))
371 0 : dikv(ii,ii) = 1 - 2 * dikv(ii,ii)
372 : END DO
373 0 : xv(1:3) = ABS(xv(1:3) - DBLE(ikv(1:3,1)))
374 0 : x = xv(1)
375 0 : y = xv(2)
376 0 : z = xv(3)
377 : !
378 0 : ikv(1:3, 2) = ikv(1:3,1) + dikv(1:3,1)
379 0 : ikv(1:3, 3) = ikv(1:3,1) + dikv(1:3,2)
380 0 : ikv(1:3, 4) = ikv(1:3,1) + dikv(1:3,3)
381 : !
382 0 : IF(nintp == 4) THEN
383 : !
384 0 : wintp(1) = 1d0 - x - y - z
385 0 : wintp(2) = x
386 0 : wintp(3) = y
387 0 : wintp(4) = z
388 : !
389 : ELSE
390 : !
391 0 : ikv(1:3, 5) = ikv(1:3,1) + SUM(dikv(1:3,1:3), 2)
392 : !
393 0 : ikv(1:3, 6) = ikv(1:3,1) - dikv(1:3,1)
394 0 : ikv(1:3, 7) = ikv(1:3,1) - dikv(1:3,2)
395 0 : ikv(1:3, 8) = ikv(1:3,1) - dikv(1:3,3)
396 : !
397 0 : ikv(1:3, 9) = ikv(1:3,1) + 2*dikv(1:3,1)
398 0 : ikv(1:3,10) = ikv(1:3,1) + 2*dikv(1:3,2)
399 0 : ikv(1:3,11) = ikv(1:3,1) + 2*dikv(1:3,3)
400 : !
401 0 : ikv(1:3,12) = ikv(1:3,1) + dikv(1:3,2) + dikv(1:3,3)
402 0 : ikv(1:3,13) = ikv(1:3,1) + dikv(1:3,3) + dikv(1:3,1)
403 0 : ikv(1:3,14) = ikv(1:3,1) + dikv(1:3,1) + dikv(1:3,2)
404 : !
405 0 : ikv(1:3,15) = ikv(1:3,1) - dikv(1:3,1) + dikv(1:3,3)
406 0 : ikv(1:3,16) = ikv(1:3,1) - dikv(1:3,2) + dikv(1:3,1)
407 0 : ikv(1:3,17) = ikv(1:3,1) - dikv(1:3,3) + dikv(1:3,2)
408 : !
409 0 : ikv(1:3,18) = ikv(1:3,1) + dikv(1:3,1) - dikv(1:3,3)
410 0 : ikv(1:3,19) = ikv(1:3,1) + dikv(1:3,2) - dikv(1:3,1)
411 0 : ikv(1:3,20) = ikv(1:3,1) + dikv(1:3,3) - dikv(1:3,2)
412 : !
413 : wintp( 1) = ( (x - 2d0)*(x - 1d0)*(1d0 + x) &
414 : & + (y - 2d0)*(y - 1d0)*(1d0 + y) &
415 : & + (z - 2d0)*(z - 1d0)*(1d0 + z) &
416 : & + 2d0*(x*y + y*z + z*x)*(x + y + z - 1d0) &
417 0 : & - 8d0*x*y*z - 4d0) * 0.5d0
418 : wintp( 2) = x * ( 2d0 + x*(1d0 - x - y - z) &
419 : & + y*(1d0 - 2d0*y + z) &
420 0 : & + z*(1d0 - 2d0*z + y)) * 0.5d0
421 : wintp( 3) = y * ( 2d0 + y*(1d0 - x - y - z) &
422 : & + x*(1d0 - 2d0*x + z) &
423 0 : & + z*(1d0 - 2d0*z + x)) * 0.5d0
424 : wintp( 4) = z * ( 2d0 + z*(1d0 - x - y - z) &
425 : & + y*(1d0 - 2d0*y + x) &
426 0 : & + x*(1d0 - 2d0*x + y)) * 0.5d0
427 0 : wintp( 5) = x * y * z
428 0 : wintp( 6) = x * (1d0 - x) * ( x + 3d0*y + 3d0*z - 2d0) / 6d0
429 0 : wintp( 7) = y * (1d0 - y) * (3d0*x + y + 3d0*z - 2d0) / 6d0
430 0 : wintp( 8) = z * (1d0 - z) * (3d0*x + 3d0*y + z - 2d0) / 6d0
431 0 : wintp( 9) = x * (x - 1d0) * (x + 1d0) / 6d0
432 0 : wintp(10) = y * (y - 1d0) * (y + 1d0) / 6d0
433 0 : wintp(11) = z * (z - 1d0) * (z + 1d0) / 6d0
434 0 : wintp(12) = y * z * (y + z - 2d0 * x) * 0.5d0
435 0 : wintp(13) = z * x * (z + x - 2d0 * y) * 0.5d0
436 0 : wintp(14) = x * y * (x + y - 2d0 * z) * 0.5d0
437 0 : wintp(15) = x * z * (x - 1d0) * 0.5d0
438 0 : wintp(16) = x * y * (y - 1d0) * 0.5d0
439 0 : wintp(17) = y * z * (z - 1d0) * 0.5d0
440 0 : wintp(18) = x * z * (z - 1d0) * 0.5d0
441 0 : wintp(19) = x * y * (x - 1d0) * 0.5d0
442 0 : wintp(20) = y * z * (y - 1d0) * 0.5d0
443 : !
444 : END IF
445 : !
446 0 : DO ii = 1, nintp
447 0 : ikv(1:3,ii) = MODULO(ikv(1:3,ii), ng(1:3))
448 0 : kintp(ii) = 1 + ikv(1,ii) + ng(1) * ikv(2,ii) + ng(1) * ng(2) * ikv(3,ii)
449 : END DO
450 : !
451 0 : END SUBROUTINE libtetrabz_interpol_indx
452 :
453 74760 : pure function a_from_e(e) result(a)
454 :
455 : !IMPLICIT NONE
456 : REAL(8),INTENT(IN) :: e(4)
457 : REAL(8) :: a(4,4)
458 :
459 : INTEGER :: ii
460 : REAL(8) :: ediff(4)
461 :
462 373800 : DO ii = 1, 4
463 1495200 : ediff = e(1:4) - e(ii)
464 1495200 : where (abs(ediff) < 1.e-10)
465 : ediff = 1.e-10
466 : end where
467 1569960 : a(1:4,ii) = (0d0 - e(ii)) / ediff
468 : END DO
469 :
470 74760 : end function a_from_e
471 :
472 : !
473 : ! Cut small tetrahedron A1
474 : !
475 0 : pure SUBROUTINE libtetrabz_tsmall_a1(e,V,tsmall)
476 : !
477 : !IMPLICIT NONE
478 : !
479 : REAL(8),INTENT(IN) :: e(4)
480 : REAL(8),INTENT(OUT) :: V
481 : REAL(8),INTENT(OUT) :: tsmall(4,4)
482 : !
483 : !INTEGER :: ii
484 : REAL(8) :: a(4,4)
485 : !
486 0 : a = a_from_e(e)
487 : !DO ii = 1, 4
488 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
489 : !END DO
490 : !
491 0 : V = a(2,1) * a(3,1) * a(4,1)
492 : !
493 0 : tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
494 0 : tsmall(2, 1:4) = (/a(1,2), a(2,1), 0d0, 0d0/)
495 0 : tsmall(3, 1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
496 0 : tsmall(4, 1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
497 : !
498 0 : END SUBROUTINE libtetrabz_tsmall_a1
499 : !
500 : ! Cut small tetrahedron B1
501 : !
502 0 : pure SUBROUTINE libtetrabz_tsmall_b1(e,V,tsmall)
503 : !
504 : !IMPLICIT NONE
505 : !
506 : REAL(8),INTENT(IN) :: e(4)
507 : REAL(8),INTENT(OUT) :: V
508 : REAL(8),INTENT(OUT) :: tsmall(4,4)
509 : !
510 : !INTEGER :: ii
511 : REAL(8) :: a(4,4)
512 : !
513 0 : a = a_from_e(e)
514 : !DO ii = 1, 4
515 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
516 : !END DO
517 : !
518 0 : V = a(3,1) * a(4,1) * a(2,4)
519 : !
520 0 : tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
521 0 : tsmall(2, 1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
522 0 : tsmall(3, 1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
523 0 : tsmall(4, 1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
524 : !
525 0 : END SUBROUTINE libtetrabz_tsmall_b1
526 : !
527 : ! Cut small tetrahedron B2
528 : !
529 0 : pure SUBROUTINE libtetrabz_tsmall_b2(e,V,tsmall)
530 : !
531 : !IMPLICIT NONE
532 : !
533 : REAL(8),INTENT(IN) :: e(4)
534 : REAL(8),INTENT(OUT) :: V
535 : REAL(8),INTENT(OUT) :: tsmall(4,4)
536 : !
537 : !INTEGER :: ii
538 : REAL(8) :: a(4,4)
539 : !
540 0 : a = a_from_e(e)
541 : !DO ii = 1, 4
542 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
543 : !END DO
544 : !
545 0 : V = a(3,2) * a(4,2)
546 : !
547 0 : tsmall(1, 1:4) = (/1d0, 0d0, 0d0, 0d0/)
548 0 : tsmall(2, 1:4) = (/0d0, 1d0, 0d0, 0d0/)
549 0 : tsmall(3, 1:4) = (/0d0, a(2,3), a(3,2), 0d0/)
550 0 : tsmall(4, 1:4) = (/0d0, a(2,4), 0d0, a(4,2)/)
551 : !
552 0 : END SUBROUTINE libtetrabz_tsmall_b2
553 : !
554 : ! Cut small tetrahedron B3
555 : !
556 0 : pure SUBROUTINE libtetrabz_tsmall_b3(e,V,tsmall)
557 : !
558 : !IMPLICIT NONE
559 : !
560 : REAL(8),INTENT(IN) :: e(4)
561 : REAL(8),INTENT(OUT) :: V
562 : REAL(8),INTENT(OUT) :: tsmall(4,4)
563 : !
564 : !INTEGER :: ii
565 : REAL(8) :: a(4,4)
566 : !
567 0 : a = a_from_e(e)
568 : !DO ii = 1, 4
569 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
570 : !END DO
571 : !
572 0 : V = a(2,3) * a(3,1) * a(4,2)
573 : !
574 0 : tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
575 0 : tsmall(2, 1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
576 0 : tsmall(3, 1:4) = (/ 0d0, a(2,3), a(3,2), 0d0/)
577 0 : tsmall(4, 1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
578 : !
579 0 : END SUBROUTINE libtetrabz_tsmall_b3
580 : !
581 : ! Cut small tetrahedron C1
582 : !
583 0 : pure SUBROUTINE libtetrabz_tsmall_c1(e,V,tsmall)
584 : !
585 : !IMPLICIT NONE
586 : !
587 : REAL(8),INTENT(IN) :: e(4)
588 : REAL(8),INTENT(OUT) :: V
589 : REAL(8),INTENT(OUT) :: tsmall(4,4)
590 : !
591 : !INTEGER :: ii
592 : REAL(8) :: a(4,4)
593 : !
594 0 : a = a_from_e(e)
595 : !DO ii = 1, 4
596 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
597 : !END DO
598 : !
599 0 : V = a(4,3)
600 : !
601 0 : tsmall(1, 1:4) = (/1d0, 0d0, 0d0, 0d0/)
602 0 : tsmall(2, 1:4) = (/0d0, 1d0, 0d0, 0d0/)
603 0 : tsmall(3, 1:4) = (/0d0, 0d0, 1d0, 0d0/)
604 0 : tsmall(4, 1:4) = (/0d0, 0d0, a(3,4), a(4,3)/)
605 : !
606 0 : END SUBROUTINE libtetrabz_tsmall_c1
607 : !
608 : ! Cut small tetrahedron C2
609 : !
610 0 : pure SUBROUTINE libtetrabz_tsmall_c2(e,V,tsmall)
611 : !
612 : !IMPLICIT NONE
613 : !
614 : REAL(8),INTENT(IN) :: e(4)
615 : REAL(8),INTENT(OUT) :: V
616 : REAL(8),INTENT(OUT) :: tsmall(4,4)
617 : !
618 : !INTEGER :: ii
619 : REAL(8) :: a(4,4)
620 : !
621 0 : a = a_from_e(e)
622 : !DO ii = 1, 4
623 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
624 : !END DO
625 : !
626 0 : V = a(3,4) * a(4,2)
627 : !
628 0 : tsmall(1, 1:4) = (/1d0, 0d0, 0d0, 0d0/)
629 0 : tsmall(2, 1:4) = (/0d0, 1d0, 0d0, 0d0/)
630 0 : tsmall(3, 1:4) = (/0d0, a(2,4), 0d0, a(4,2)/)
631 0 : tsmall(4, 1:4) = (/0d0, 0d0, a(3,4), a(4,3)/)
632 : !
633 0 : END SUBROUTINE libtetrabz_tsmall_c2
634 : !
635 : ! Cut small tetrahedron C3
636 : !
637 0 : pure SUBROUTINE libtetrabz_tsmall_c3(e,V,tsmall)
638 : !
639 : !IMPLICIT NONE
640 : !
641 : REAL(8),INTENT(IN) :: e(4)
642 : REAL(8),INTENT(OUT) :: V
643 : REAL(8),INTENT(OUT) :: tsmall(4,4)
644 : !
645 : !INTEGER :: ii
646 : REAL(8) :: a(4,4)
647 : !
648 0 : a = a_from_e(e)
649 : !DO ii = 1, 4
650 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
651 : !END DO
652 : !
653 0 : V = a(3,4) * a(2,4) * a(4,1)
654 : !
655 0 : tsmall(1, 1:4) = (/ 1d0, 0d0, 0d0, 0d0/)
656 0 : tsmall(2, 1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
657 0 : tsmall(3, 1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
658 0 : tsmall(4, 1:4) = (/ 0d0, 0d0, a(3,4), a(4,3)/)
659 : !
660 0 : END SUBROUTINE libtetrabz_tsmall_c3
661 : !
662 : ! Cut triangle A1
663 : !
664 28560 : pure SUBROUTINE libtetrabz_triangle_a1(e,V,tsmall)
665 : !
666 : !IMPLICIT NONE
667 : !
668 : REAL(8),INTENT(IN) :: e(4)
669 : REAL(8),INTENT(OUT) :: V
670 : REAL(8),INTENT(OUT) :: tsmall(3,4)
671 : !
672 : !INTEGER :: ii
673 : REAL(8) :: a(4,4)
674 : !
675 28560 : a = a_from_e(e)
676 : !DO ii = 1, 4
677 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
678 : !END DO
679 : !
680 : !V = 3d0 * a(2,1) * a(3,1) * a(4,1) / (0d0 - e(1))
681 28560 : V = 3d0 * a(2,1) * a(3,1) / (e(4) - e(1))
682 : !
683 142800 : tsmall(1,1:4) = (/a(1,2), a(2,1), 0d0, 0d0/)
684 142800 : tsmall(2,1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
685 142800 : tsmall(3,1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
686 : !
687 28560 : END SUBROUTINE libtetrabz_triangle_a1
688 : !
689 : ! Cut triangle B1
690 : !
691 17640 : pure SUBROUTINE libtetrabz_triangle_b1(e,V,tsmall)
692 : !
693 : !IMPLICIT NONE
694 : !
695 : REAL(8),INTENT(IN) :: e(4)
696 : REAL(8),INTENT(OUT) :: V
697 : REAL(8),INTENT(OUT) :: tsmall(3,4)
698 : !
699 : !INTEGER :: ii
700 : REAL(8) :: a(4,4)
701 : !
702 17640 : a = a_from_e(e)
703 : !DO ii = 1, 4
704 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
705 : !END DO
706 : !
707 : !V = 3d0 * a(3,1) * a(4,1) * a(2,4) / (0d0 - e(1))
708 17640 : V = 3d0 * a(4,1) * a(2,4) / (e(3) - e(1))
709 : !
710 88200 : tsmall(1,1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
711 88200 : tsmall(2,1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
712 88200 : tsmall(3,1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
713 : !
714 17640 : END SUBROUTINE libtetrabz_triangle_b1
715 : !
716 : ! Cut triangle B2
717 : !
718 17640 : pure SUBROUTINE libtetrabz_triangle_b2(e,V,tsmall)
719 : !
720 : !IMPLICIT NONE
721 : !
722 : REAL(8),INTENT(IN) :: e(4)
723 : REAL(8),INTENT(OUT) :: V
724 : REAL(8),INTENT(OUT) :: tsmall(3,4)
725 : !
726 : !INTEGER :: ii
727 : REAL(8) :: a(4,4)
728 : !
729 17640 : a = a_from_e(e)
730 : !DO ii = 1, 4
731 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
732 : !END DO
733 : !
734 : !V = 3d0 * a(2,3) * a(3,1) * a(4,2) / (0d0 - e(1))
735 17640 : V = 3d0 * a(2,3) * a(4,2) / (e(3) - e(1))
736 : !
737 88200 : tsmall(1,1:4) = (/a(1,3), 0d0, a(3,1), 0d0/)
738 88200 : tsmall(2,1:4) = (/ 0d0, a(2,3), a(3,2), 0d0/)
739 88200 : tsmall(3,1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
740 : !
741 17640 : END SUBROUTINE libtetrabz_triangle_b2
742 : !
743 : ! Cut triangle C1
744 : !
745 10920 : pure SUBROUTINE libtetrabz_triangle_c1(e,V,tsmall)
746 : !
747 : !IMPLICIT NONE
748 : !
749 : REAL(8),INTENT(IN) :: e(4)
750 : REAL(8),INTENT(OUT) :: V
751 : REAL(8),INTENT(OUT) :: tsmall(3,4)
752 : !
753 : !INTEGER :: ii
754 : REAL(8) :: a(4,4)
755 : !
756 10920 : a = a_from_e(e)
757 : !DO ii = 1, 4
758 : ! a(1:4,ii) = (0d0 - e(ii)) / (e(1:4) - e(ii))
759 : !END DO
760 : !
761 : !V = 3d0 * a(1,4) * a(2,4) * a(3,4) / (e(4) - 0d0)
762 10920 : V = 3d0 * a(1,4) * a(2,4) / (e(4) - e(3))
763 : !
764 54600 : tsmall(1,1:4) = (/a(1,4), 0d0, 0d0, a(4,1)/)
765 54600 : tsmall(2,1:4) = (/ 0d0, a(2,4), 0d0, a(4,2)/)
766 54600 : tsmall(3,1:4) = (/ 0d0, 0d0, a(3,4), a(4,3)/)
767 : !
768 10920 : END SUBROUTINE libtetrabz_triangle_c1
769 : !
770 : ! MPI_Allreduce for double scaler
771 : !
772 0 : SUBROUTINE libtetrabz_mpisum_d(comm,scaler)
773 : !
774 : !IMPLICIT NONE
775 : !
776 : INTEGER :: comm
777 : REAL(8) :: scaler
778 : !
779 : #if defined(HAVE_MPI)
780 : INTEGER :: ierr
781 : REAL(8) :: arr_scaler(1)
782 : !
783 : CALL MPI_allREDUCE([scaler], arr_scaler, 1, &
784 0 : & MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr)
785 0 : scaler=arr_scaler(1)
786 : #endif
787 : !
788 0 : END SUBROUTINE libtetrabz_mpisum_d
789 : !
790 : ! MPI_Allreduce for double vector
791 : !
792 0 : SUBROUTINE libtetrabz_mpisum_dv(comm,ndim,vector)
793 : !
794 : !IMPLICIT NONE
795 : !
796 : INTEGER :: comm, ndim
797 : REAL(8) :: vector(ndim)
798 : !
799 : #if defined(HAVE_MPI)
800 : INTEGER :: ierr
801 : #ifndef HAVE_MPI2_INPLACE
802 : REAL(8) :: vector_out(ndim)
803 : #endif
804 : !
805 : #ifdef HAVE_MPI2_INPLACE
806 : CALL MPI_allREDUCE(MPI_IN_PLACE, vector, ndim, &
807 0 : & MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr)
808 : #else
809 : CALL MPI_allREDUCE(vector, vector_out, ndim, &
810 : & MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr)
811 : vector(1:ndim)=vector_out(1:ndim)
812 : #endif
813 : #endif
814 : !
815 0 : END SUBROUTINE libtetrabz_mpisum_dv
816 : !
817 : ! MPI_Allreduce for double complex vector
818 : !
819 0 : SUBROUTINE libtetrabz_mpisum_zv(comm,ndim,vector)
820 : !
821 : !IMPLICIT NONE
822 : !
823 : INTEGER :: comm, ndim
824 : COMPLEX(8) :: vector(ndim)
825 : !
826 : #if defined(HAVE_MPI)
827 : INTEGER :: ierr
828 : #ifndef HAVE_MPI2_INPLACE
829 : COMPLEX(8) :: vector_out(ndim)
830 : #endif
831 : !
832 : #ifdef HAVE_MPI2_INPLACE
833 : CALL MPI_allREDUCE(MPI_IN_PLACE, vector, ndim, &
834 0 : & MPI_DOUBLE_COMPLEX, MPI_SUM, comm, ierr)
835 : #else
836 : CALL MPI_allREDUCE(vector, vector_out, ndim, &
837 : & MPI_DOUBLE_COMPLEX, MPI_SUM, comm, ierr)
838 : vector(1:ndim)=vector_out(1:ndim)
839 : #endif
840 : #endif
841 : !
842 0 : END SUBROUTINE libtetrabz_mpisum_zv
843 : !
844 : END MODULE libtetrabz_common
|
