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_dbldelta_mod
30 : : !
31 : :
32 : : use m_abicore
33 : : use m_errors
34 : :
35 : : !IMPLICIT NONE
36 : : !
37 : : PRIVATE
38 : : PUBLIC :: libtetrabz_dbldelta
39 : : !
40 : : CONTAINS
41 : : !
42 : : ! Compute doubledelta
43 : : !
44 : 35 : SUBROUTINE libtetrabz_dbldelta(ltetra,bvec,nb,nge,eig1,eig2,ngw,wght,comm)
45 : : !
46 : : use, intrinsic :: iso_c_binding
47 : : USE libtetrabz_common, ONLY : libtetrabz_initialize, libtetrabz_interpol_indx, libtetrabz_mpisum_dv
48 : : !IMPLICIT NONE
49 : : !
50 : : INTEGER(C_INT),INTENT(IN) :: ltetra, nb, nge(3), ngw(3)
51 : : REAL(C_DOUBLE),INTENT(IN) :: bvec(9), eig1(nb,PRODUCT(nge(1:3))), eig2(nb,PRODUCT(nge(1:3)))
52 : : REAL(C_DOUBLE),INTENT(OUT) :: wght(nb*nb,PRODUCT(ngw(1:3)))
53 : : INTEGER(C_INT),INTENT(IN),OPTIONAL :: comm
54 : : !
55 : : LOGICAL :: linterpol
56 : : INTEGER :: nt_local, nk_local, nkBZ, ik, kintp(20), nintp
57 : 35 : INTEGER,ALLOCATABLE :: ik_global(:,:), ik_local(:,:)
58 : : REAL(8) :: wlsm(4,20), wintp(1,20)
59 : 35 : REAL(8),ALLOCATABLE :: wghtd(:,:,:), kvec(:,:)
60 : : !
61 : 35 : nintp = 16 * ltetra - 12
62 : : !
63 [ + - ]: 35 : IF(PRESENT(comm)) THEN
64 : : CALL libtetrabz_initialize(ltetra,nge,ngw,bvec,linterpol,wlsm,nk_local,&
65 : 35 : & nt_local,nkBZ,ik_global,ik_local,kvec,comm)
66 : : ELSE
67 : : CALL libtetrabz_initialize(ltetra,nge,ngw,bvec,linterpol,wlsm,nk_local,&
68 : 0 : & nt_local,nkBZ,ik_global,ik_local,kvec)
69 : : END IF
70 : : !
71 [ - + ]: 35 : IF(linterpol) THEN
72 : : !
73 [ # # # # : 0 : ABI_MALLOC(wghtd, (nb*nb,1,nk_local))
# # # # #
# # # # #
# # ]
74 : 0 : CALL libtetrabz_dbldelta_main(wlsm,nt_local,ik_global,ik_local,nb,nkBZ,eig1,eig2,nk_local,wghtd)
75 : : !
76 : : ! Interpolation
77 : : !
78 [ # # # # : 0 : wght(1:nb*nb,1:PRODUCT(ngw(1:3))) = 0d0
# # ]
79 [ # # ]: 0 : DO ik = 1, nk_local
80 : 0 : CALL libtetrabz_interpol_indx(nintp,ngw,kvec(1:3,ik),kintp,wintp)
81 : : wght(1:nb*nb,kintp(1:nintp)) = wght(1:nb*nb, kintp(1:nintp)) &
82 [ # # # # : 0 : & + MATMUL(wghtd(1:nb*nb,1:1,ik), wintp(1:1,1:nintp))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
83 : : END DO ! ik = 1, nk_local
84 [ # # ]: 0 : ABI_FREE(wghtd)
85 [ # # ]: 0 : ABI_FREE(kvec)
86 : : !
87 [ # # # # ]: 0 : IF(PRESENT(comm)) CALL libtetrabz_mpisum_dv(comm, nb*nb*PRODUCT(ngw(1:3)), wght)
88 : : !
89 : : ELSE
90 : 35 : CALL libtetrabz_dbldelta_main(wlsm,nt_local,ik_global,ik_local,nb,nkBZ,eig1,eig2,nk_local,wght)
91 : : END IF
92 : : !
93 [ - + ]: 35 : ABI_FREE(ik_global)
94 [ - + ]: 35 : ABI_FREE(ik_local)
95 : : !
96 [ - + ]: 35 : END SUBROUTINE libtetrabz_dbldelta
97 : : !
98 : : ! Main SUBROUTINE for Delta(E1) * Delta(E2)
99 : : !
100 : 35 : SUBROUTINE libtetrabz_dbldelta_main(wlsm,nt_local,ik_global,ik_local,nb,nkBZ,eig1,eig2,nk_local,dbldelta)
101 : : !
102 : : USE libtetrabz_common, ONLY : libtetrabz_sort, &
103 : : & libtetrabz_triangle_a1, libtetrabz_triangle_b1, &
104 : : & libtetrabz_triangle_b2, libtetrabz_triangle_c1
105 : : !IMPLICIT NONE
106 : : !
107 : : INTEGER,INTENT(IN) :: nt_local, nb, nkBZ, nk_local, &
108 : : & ik_global(20,nt_local), ik_local(20,nt_local)
109 : : REAL(8),INTENT(IN) :: wlsm(4,20), eig1(nb,nkBZ), eig2(nb,nkBZ)
110 : : REAL(8),INTENT(OUT) :: dbldelta(nb,nb,nk_local)
111 : : !
112 : : INTEGER :: ib, indx(4), it
113 : 70 : REAL(8) :: e(4), ei1(4,nb), ej1(4,nb), ej2(3,nb), V, thr = 1d-10, &
114 : 35 : & tsmall(3,4), w1(nb,4), w2(nb,3)
115 : : !
116 [ + + + + : 125475 : dbldelta(1:nb,1:nb,1:nk_local) = 0d0
+ + ]
117 : : !
118 : : !$OMP PARALLEL DEFAULT(NONE) &
119 : : !$OMP & SHARED(dbldelta,eig1,eig2,ik_global,ik_local,nb,nt_local,thr,wlsm) &
120 : : !$OMP & PRIVATE(e,ei1,ej1,ej2,ib,indx,it,tsmall,V,w1,w2)
121 : : !
122 [ + + ]: 107555 : DO it = 1, nt_local
123 : : !
124 [ + + ]: 322560 : DO ib = 1, nb
125 [ + + + + : 27095040 : ei1(1:4,ib) = MATMUL(wlsm(1:4,1:20), eig1(ib, ik_global(1:20,it)))
+ + + + ]
126 [ + + + + : 27417600 : ej1(1:4,ib) = MATMUL(wlsm(1:4,1:20), eig2(ib, ik_global(1:20,it)))
+ + + + ]
127 : : END DO
128 : : !
129 : : !$OMP DO
130 [ + + ]: 322595 : DO ib = 1, nb
131 : : !
132 [ + + + + ]: 2795520 : w1(1:nb,1:4) = 0d0
133 [ + + ]: 1075200 : e(1:4) = ei1(1:4, ib)
134 : 215040 : CALL libtetrabz_sort(4,e,indx)
135 : : !
136 [ + + + + ]: 215040 : IF(e(1) < 0d0 .AND. 0d0 <= e(2)) THEN
137 : : !
138 : 28560 : CALL libtetrabz_triangle_a1(e,V,tsmall)
139 : : !
140 [ + - ]: 28560 : IF(V > thr) THEN
141 : : !
142 [ + + + + : 1570800 : ej2(1:3,1:nb) = MATMUL(tsmall(1:3,1:4), ej1(indx(1:4),1:nb))
+ - + + +
+ + + + +
+ + + - ]
143 : 28560 : CALL libtetrabz_dbldelta2(nb,ej2,w2)
144 : : w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
145 [ - + + - : 1942080 : & + V * MATMUL(w2(1:nb,1:3), tsmall(1:3,1:4))
+ - - + -
+ + + + +
+ + + + +
+ + + + +
+ - ]
146 : : !
147 : : END IF
148 : : !
149 [ + + + + ]: 186480 : ELSE IF( e(2) < 0d0 .AND. 0d0 <= e(3)) THEN
150 : : !
151 : 17640 : CALL libtetrabz_triangle_b1(e,V,tsmall)
152 : : !
153 [ + - ]: 17640 : IF(V > thr) THEN
154 : : !
155 [ + + + + : 970200 : ej2(1:3,1:nb) = MATMUL(tsmall(1:3,1:4), ej1(indx(1:4),1:nb))
+ - + + +
+ + + + +
+ + + - ]
156 : 17640 : CALL libtetrabz_dbldelta2(nb,ej2,w2)
157 : : w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
158 [ - + + - : 1199520 : & + V * MATMUL(w2(1:nb,1:3), tsmall(1:3,1:4))
+ - - + -
+ + + + +
+ + + + +
+ + + + +
+ - ]
159 : : !
160 : : END IF
161 : : !
162 : 17640 : CALL libtetrabz_triangle_b2(e,V,tsmall)
163 : : !
164 [ + - ]: 17640 : IF(V > thr) THEN
165 : : !
166 [ + + + + : 970200 : ej2(1:3,1:nb) = MATMUL(tsmall(1:3,1:4), ej1(indx(1:4),1:nb))
+ - + + +
+ + + + +
+ + + - ]
167 : 17640 : CALL libtetrabz_dbldelta2(nb,ej2,w2)
168 : : w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
169 [ - + + - : 1199520 : & + V * MATMUL(w2(1:nb,1:3), tsmall(1:3,1:4))
+ - - + -
+ + + + +
+ + + + +
+ + + + +
+ - ]
170 : : !
171 : : END IF
172 : : !
173 [ + + + + ]: 168840 : ELSE IF(e(3) < 0d0 .AND. 0d0 < e(4)) THEN
174 : : !
175 : 10920 : CALL libtetrabz_triangle_c1(e,V,tsmall)
176 : : !
177 [ + - ]: 10920 : IF(V > thr) THEN
178 : : !
179 [ + + + + : 600600 : ej2(1:3,1:nb) = MATMUL(tsmall(1:3,1:4), ej1(indx(1:4),1:nb))
+ - + + +
+ + + + +
+ + + - ]
180 : 10920 : CALL libtetrabz_dbldelta2(nb,ej2,w2)
181 : : w1(1:nb,indx(1:4)) = w1(1:nb, indx(1:4)) &
182 [ - + + - : 742560 : & + V * MATMUL(w2(1:nb,1:3), tsmall(1:3,1:4))
+ - - + -
+ + + + +
+ + + + +
+ + + + +
+ - ]
183 : : !
184 : : END IF
185 : : !
186 : : END IF
187 : : !
188 : : dbldelta(1:nb,ib,ik_local(1:20,it)) = dbldelta(1:nb,ib, ik_local(1:20,it)) &
189 [ - + + - : 95800320 : & + MATMUL(w1(1:nb,1:4), wlsm(1:4,1:20))
+ - - + -
+ + + + +
+ + + + +
+ + + + +
+ + + + +
- ]
190 : : !
191 : : END DO ! ib
192 : : !$OMP END DO NOWAIT
193 : : !
194 : : END DO ! it
195 : : !
196 : : !$OMP END PARALLEL
197 : : !
198 [ + + + + : 125475 : dbldelta(1:nb,1:nb,1:nk_local) = dbldelta(1:nb,1:nb,1:nk_local) / DBLE(6 * nkBZ)
+ + ]
199 : : !
200 : 35 : END SUBROUTINE libtetrabz_dbldelta_main
201 : : !
202 : : ! 2nd step of tetrahedra method.
203 : : !
204 : 74760 : SUBROUTINE libtetrabz_dbldelta2(nb,ej,w)
205 : : !
206 : : USE libtetrabz_common, ONLY : libtetrabz_sort
207 : : !IMPLICIT NONE
208 : : !
209 : : INTEGER,INTENT(IN) :: nb
210 : : REAL(8),INTENT(IN) :: ej(3,nb)
211 : : REAL(8),INTENT(INOUT) :: w(nb,3)
212 : : !
213 : : INTEGER :: ib, ii, indx(3)
214 : : REAL(8) :: a(3,3), e(3), V
215 : : REAL(8) :: ediff(3)
216 : : !character(len=500) :: msg
217 : : !
218 [ + + ]: 224280 : DO ib = 1, nb
219 : : !
220 [ - + + - : 747600 : IF(maxval(ABS(ej(1:3,ib))) < 1d-10) then
+ + - + ]
221 : : ! MG reduce tolerance wrt original version.
222 : : !IF(maxval(ABS(ej(1:3,ib))) < 1d-22) then
223 : : !write(std_out, ej(1:3,ib))
224 : : !write(msg, *)"Nesting for band index:", ib, "ej:", ej(1:3,ib)
225 : : !ABI_WARNING(msg)
226 : : !ABI_ERROR(msg)
227 : : !ABI_ERROR("STOP Nesting !!")
228 [ # # ]: 0 : w(ib,1:3) = 0d0
229 : : cycle
230 : : end if
231 : : !
232 [ + + ]: 598080 : w(ib,1:3) = 0d0
233 [ + + ]: 598080 : e(1:3) = ej(1:3, ib)
234 : 149520 : CALL libtetrabz_sort(3,e,indx)
235 : : !
236 [ + + ]: 598080 : DO ii = 1, 3
237 [ + + ]: 1794240 : ediff = e(1:3) - e(ii)
238 [ + + + + ]: 1794240 : where (abs(ediff) < 1.e-10)
239 : : ediff = 1.e-10
240 : : end where
241 [ + + ]: 1943760 : a(1:3,ii) = (0d0 - e(ii)) / ediff
242 : : !a(1:3,ii) = (0d0 - e(ii)) / (e(1:3) - e(ii))
243 : : END DO
244 : : !
245 [ + + + + : 224280 : IF((e(1) < 0d0 .AND. 0d0 <= e(2)) .OR. (e(1) <= 0d0 .AND. 0d0 < e(2))) THEN
+ + - + ]
246 : : !
247 : : !V = a(2,1) * a(3,1) / (0d0 - e(1))
248 : 7355 : V = a(2,1) / (e(3) - e(1))
249 : : !
250 : 7355 : w(ib,indx(1)) = V * (a(1,2) + a(1,3))
251 : 7355 : w(ib,indx(2)) = V * a(2,1)
252 : 7355 : w(ib,indx(3)) = V * a(3,1)
253 : : !
254 [ + + + + : 142165 : ELSE IF((e(2) <= 0d0 .AND. 0d0 < e(3)) .OR. (e(2) < 0d0 .AND. 0d0 <= e(3))) THEN
+ + - + ]
255 : : !
256 : : !V = a(1,3) * a(2,3) / (e(3) - 0d0)
257 : 6848 : V = a(2,3) / (e(3) - e(1))
258 : : !
259 : 6848 : w(ib,indx(1)) = V * a(1,3)
260 : 6848 : w(ib,indx(2)) = V * a(2,3)
261 : 6848 : w(ib,indx(3)) = V * (a(3,1) + a(3,2))
262 : : !
263 : : END IF
264 : : !
265 : : END DO ! ib
266 : : !
267 : 74760 : END SUBROUTINE libtetrabz_dbldelta2
268 : : !
269 [ - - + - : 289800 : END MODULE libtetrabz_dbldelta_mod
+ - + - +
- + - ]
|
