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
|
