Branch data 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
|
