Line data Source code
1 : !!****m* ABINIT/m_epweights
2 : !! NAME
3 : !! m_epweights
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2012-2026 ABINIT group (BXU, MVer)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_epweights
23 :
24 : use defs_basis
25 : use defs_elphon
26 : use m_abicore
27 : use m_errors
28 : use m_tetrahedron
29 : !use m_htetra
30 : use m_xmpi
31 :
32 : use m_matrix, only : matr3inv
33 :
34 : implicit none
35 :
36 : private
37 : !!***
38 :
39 : public :: d2c_weights
40 : public :: d2c_wtq!
41 : public :: ep_el_weights
42 : public :: ep_fs_weights
43 : public :: ep_ph_weights
44 : !!***
45 :
46 : contains
47 : !!***
48 :
49 : !!****f* ABINIT/d2c_weights
50 : !! NAME
51 : !! d2c_weights
52 : !!
53 : !! FUNCTION
54 : !! This routine calculates the integration weights on a coarse k-grid
55 : !! using the integration weights from a denser k-gird. The weights of
56 : !! the extra k points that being shared are evenly distributed. It also
57 : !! condenses the velocity*wtk and velcity^2*wtk.
58 : !!
59 : !! INPUTS
60 : !! elph_ds%k_fine%nkpt = number of fine FS k-points
61 : !! elph_ds%k_fine%wtk = integration weights of the fine FS k-grid
62 : !! elph_ds%k_phon%nkpt = number of coarse FS k-points
63 : !! elph_tr_ds%el_veloc = electronic velocities from the fine k-grid
64 : !!
65 : !! OUTPUT
66 : !! elph_ds%k_phon%wtk = integration weights of the coarse FS k-grid
67 : !! elph_ds%k_phon%velocwtk = velocity time integration weights of the coarse FS k-grid
68 : !! elph_ds%k_phon%vvelocwtk = velocity^2 time integration weights of the coarse FS k-grid
69 : !!
70 : !! SOURCE
71 :
72 0 : subroutine d2c_weights(elph_ds,elph_tr_ds)
73 :
74 : !Arguments ------------------------------------
75 : type(elph_type),intent(inout) :: elph_ds
76 : type(elph_tr_type),intent(inout),optional :: elph_tr_ds
77 :
78 : !Local variables-------------------------------
79 : integer :: ii, jj, kk
80 : integer :: ikpt, jkpt, kkpt
81 : integer :: iikpt, jjkpt, kkkpt
82 : integer :: icomp, jcomp
83 : integer :: iFSband, iband
84 : integer :: ikpt_fine, ikpt_phon
85 : integer :: nkpt_fine1, nkpt_phon1
86 : integer :: nkpt_fine2, nkpt_phon2
87 : integer :: nkpt_fine3, nkpt_phon3
88 : integer :: nscale1, nscale2, nscale3
89 :
90 : ! *************************************************************************
91 0 : nkpt_phon1 = elph_ds%kptrlatt(1,1)
92 0 : nkpt_phon2 = elph_ds%kptrlatt(2,2)
93 0 : nkpt_phon3 = elph_ds%kptrlatt(3,3)
94 0 : nkpt_fine1 = elph_ds%kptrlatt_fine(1,1)
95 0 : nkpt_fine2 = elph_ds%kptrlatt_fine(2,2)
96 0 : nkpt_fine3 = elph_ds%kptrlatt_fine(3,3)
97 0 : nscale1 = dble(nkpt_fine1/nkpt_phon1)
98 0 : nscale2 = dble(nkpt_fine2/nkpt_phon2)
99 0 : nscale3 = dble(nkpt_fine3/nkpt_phon3)
100 : if (abs(INT(nscale1)-nscale1) > 0.01) then
101 : ABI_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
102 : end if
103 : if (abs(INT(nscale2)-nscale2) > 0.01) then
104 : ABI_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
105 : end if
106 : if (abs(INT(nscale3)-nscale3) > 0.01) then
107 : ABI_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
108 : end if
109 0 : nscale1 = INT(nscale1)
110 0 : nscale2 = INT(nscale2)
111 0 : nscale3 = INT(nscale3)
112 :
113 : !bxu, get wtk of coarse grid from fine grid
114 0 : elph_ds%k_phon%wtk = zero
115 0 : if (present(elph_tr_ds)) then
116 0 : elph_ds%k_phon%velocwtk = zero
117 0 : elph_ds%k_phon%vvelocwtk = zero
118 : end if
119 :
120 0 : do ikpt = 1, nkpt_phon1
121 0 : do jkpt = 1, nkpt_phon2
122 0 : do kkpt = 1, nkpt_phon3
123 0 : ikpt_phon = kkpt + (jkpt-1)*nkpt_phon3 + (ikpt-1)*nkpt_phon2*nkpt_phon3
124 : ! inside the paralellepipe
125 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
126 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
127 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
128 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
129 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
130 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
131 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
132 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
133 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
134 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
135 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
136 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
137 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
138 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
139 0 : & elph_ds%k_fine%wtk(:,ikpt_fine,:)
140 0 : if (present(elph_tr_ds)) then
141 0 : do iFSband=1,elph_ds%ngkkband !FS bands
142 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
143 0 : do icomp = 1, 3
144 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
145 0 : & elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
146 0 : do jcomp = 1, 3
147 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
148 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
149 : & elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
150 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
151 : end do
152 : end do
153 : end do
154 : end if
155 : end do
156 : end do
157 : end do
158 : ! on the 6 faces
159 0 : if (MOD(nscale3,2) == 0) then ! when nscale3 is an even number
160 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
161 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
162 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
163 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
164 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
165 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
166 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
167 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
168 :
169 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
170 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
171 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
172 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
173 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
174 0 : & 0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
175 0 : if (present(elph_tr_ds)) then
176 0 : do iFSband=1,elph_ds%ngkkband !FS bands
177 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
178 0 : do icomp = 1, 3
179 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
180 0 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
181 0 : do jcomp = 1, 3
182 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
183 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
184 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
185 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
186 : end do
187 : end do
188 : end do
189 : end if
190 :
191 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
192 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
193 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
194 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
195 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
196 0 : & 0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
197 0 : if (present(elph_tr_ds)) then
198 0 : do iFSband=1,elph_ds%ngkkband !FS bands
199 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
200 0 : do icomp = 1, 3
201 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
202 0 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
203 0 : do jcomp = 1, 3
204 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
205 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
206 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
207 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
208 : end do
209 : end do
210 : end do
211 : end if
212 : end do
213 : end do
214 : end if
215 0 : if (MOD(nscale2,2) == 0) then ! when nscale2 is an even number
216 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
217 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
218 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
219 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
220 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
221 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
222 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
223 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
224 :
225 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
226 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
227 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
228 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
229 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
230 0 : & 0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
231 0 : if (present(elph_tr_ds)) then
232 0 : do iFSband=1,elph_ds%ngkkband !FS bands
233 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
234 0 : do icomp = 1, 3
235 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
236 0 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
237 0 : do jcomp = 1, 3
238 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
239 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
240 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
241 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
242 : end do
243 : end do
244 : end do
245 : end if
246 :
247 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
248 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
249 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
250 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
251 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
252 0 : & 0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
253 0 : if (present(elph_tr_ds)) then
254 0 : do iFSband=1,elph_ds%ngkkband !FS bands
255 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
256 0 : do icomp = 1, 3
257 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
258 0 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
259 0 : do jcomp = 1, 3
260 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
261 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
262 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
263 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
264 : end do
265 : end do
266 : end do
267 : end if
268 : end do
269 : end do
270 : end if
271 0 : if (MOD(nscale1,2) == 0) then ! when nscale1 is an even number
272 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
273 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
274 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
275 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
276 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
277 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
278 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
279 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
280 :
281 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
282 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
283 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
284 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
285 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
286 0 : & 0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
287 0 : if (present(elph_tr_ds)) then
288 0 : do iFSband=1,elph_ds%ngkkband !FS bands
289 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
290 0 : do icomp = 1, 3
291 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
292 0 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
293 0 : do jcomp = 1, 3
294 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
295 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
296 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
297 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
298 : end do
299 : end do
300 : end do
301 : end if
302 :
303 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
304 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
305 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
306 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
307 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
308 0 : & 0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
309 0 : if (present(elph_tr_ds)) then
310 0 : do iFSband=1,elph_ds%ngkkband !FS bands
311 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
312 0 : do icomp = 1, 3
313 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
314 0 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
315 0 : do jcomp = 1, 3
316 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
317 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
318 : & 0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
319 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
320 : end do
321 : end do
322 : end do
323 : end if
324 : end do
325 : end do
326 : ! on the 12 sides
327 : end if
328 0 : if (MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
329 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
330 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
331 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
332 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
333 :
334 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
335 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
336 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
337 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
338 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
339 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
340 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
341 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
342 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
343 0 : if (present(elph_tr_ds)) then
344 0 : do iFSband=1,elph_ds%ngkkband !FS bands
345 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
346 0 : do icomp = 1, 3
347 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
348 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
349 0 : do jcomp = 1, 3
350 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
351 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
352 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
353 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
354 : end do
355 : end do
356 : end do
357 : end if
358 :
359 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
360 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
361 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
362 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
363 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
364 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
365 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
366 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
367 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
368 0 : if (present(elph_tr_ds)) then
369 0 : do iFSband=1,elph_ds%ngkkband !FS bands
370 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
371 0 : do icomp = 1, 3
372 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
373 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
374 0 : do jcomp = 1, 3
375 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
376 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
377 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
378 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
379 : end do
380 : end do
381 : end do
382 : end if
383 :
384 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
385 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
386 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
387 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
388 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
389 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
390 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
391 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
392 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
393 0 : if (present(elph_tr_ds)) then
394 0 : do iFSband=1,elph_ds%ngkkband !FS bands
395 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
396 0 : do icomp = 1, 3
397 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
398 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
399 0 : do jcomp = 1, 3
400 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
401 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
402 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
403 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
404 : end do
405 : end do
406 : end do
407 : end if
408 :
409 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
410 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
411 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
412 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
413 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
414 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
415 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
416 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
417 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
418 0 : if (present(elph_tr_ds)) then
419 0 : do iFSband=1,elph_ds%ngkkband !FS bands
420 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
421 0 : do icomp = 1, 3
422 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
423 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
424 0 : do jcomp = 1, 3
425 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
426 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
427 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
428 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
429 : end do
430 : end do
431 : end do
432 : end if
433 : end do
434 : end if
435 0 : if (MOD(nscale1,2) == 0 .and. MOD(nscale3,2) == 0) then
436 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
437 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
438 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
439 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
440 :
441 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
442 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
443 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
444 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
445 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
446 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
447 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
448 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
449 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
450 0 : if (present(elph_tr_ds)) then
451 0 : do iFSband=1,elph_ds%ngkkband !FS bands
452 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
453 0 : do icomp = 1, 3
454 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
455 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
456 0 : do jcomp = 1, 3
457 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
458 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
459 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
460 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
461 : end do
462 : end do
463 : end do
464 : end if
465 :
466 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
467 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
468 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
469 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
470 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
471 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
472 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
473 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
474 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
475 0 : if (present(elph_tr_ds)) then
476 0 : do iFSband=1,elph_ds%ngkkband !FS bands
477 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
478 0 : do icomp = 1, 3
479 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
480 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
481 0 : do jcomp = 1, 3
482 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
483 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
484 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
485 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
486 : end do
487 : end do
488 : end do
489 : end if
490 :
491 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
492 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
493 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
494 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
495 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
496 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
497 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
498 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
499 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
500 0 : if (present(elph_tr_ds)) then
501 0 : do iFSband=1,elph_ds%ngkkband !FS bands
502 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
503 0 : do icomp = 1, 3
504 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
505 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
506 0 : do jcomp = 1, 3
507 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
508 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
509 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
510 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
511 : end do
512 : end do
513 : end do
514 : end if
515 :
516 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
517 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
518 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
519 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
520 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
521 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
522 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
523 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
524 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
525 0 : if (present(elph_tr_ds)) then
526 0 : do iFSband=1,elph_ds%ngkkband !FS bands
527 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
528 0 : do icomp = 1, 3
529 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
530 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
531 0 : do jcomp = 1, 3
532 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
533 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
534 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
535 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
536 : end do
537 : end do
538 : end do
539 : end if
540 : end do
541 : end if
542 0 : if (MOD(nscale2,2) == 0 .and. MOD(nscale1,2) == 0) then
543 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
544 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
545 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
546 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
547 :
548 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
549 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
550 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
551 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
552 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
553 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
554 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
555 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
556 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
557 0 : if (present(elph_tr_ds)) then
558 0 : do iFSband=1,elph_ds%ngkkband !FS bands
559 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
560 0 : do icomp = 1, 3
561 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
562 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
563 0 : do jcomp = 1, 3
564 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
565 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
566 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
567 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
568 : end do
569 : end do
570 : end do
571 : end if
572 :
573 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
574 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
575 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
576 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
577 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
578 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
579 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
580 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
581 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
582 0 : if (present(elph_tr_ds)) then
583 0 : do iFSband=1,elph_ds%ngkkband !FS bands
584 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
585 0 : do icomp = 1, 3
586 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
587 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
588 0 : do jcomp = 1, 3
589 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
590 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
591 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
592 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
593 : end do
594 : end do
595 : end do
596 : end if
597 :
598 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
599 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
600 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
601 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
602 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
603 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
604 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
605 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
606 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
607 0 : if (present(elph_tr_ds)) then
608 0 : do iFSband=1,elph_ds%ngkkband !FS bands
609 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
610 0 : do icomp = 1, 3
611 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
612 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
613 0 : do jcomp = 1, 3
614 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
615 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
616 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
617 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
618 : end do
619 : end do
620 : end do
621 : end if
622 :
623 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
624 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
625 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
626 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
627 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
628 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
629 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
630 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
631 0 : & 0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
632 0 : if (present(elph_tr_ds)) then
633 0 : do iFSband=1,elph_ds%ngkkband !FS bands
634 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
635 0 : do icomp = 1, 3
636 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
637 0 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
638 0 : do jcomp = 1, 3
639 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
640 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
641 : & 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
642 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
643 : end do
644 : end do
645 : end do
646 : end if
647 : end do
648 : ! on the 8 corners
649 : end if
650 0 : if (MOD(nscale1,2) == 0 .and. MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
651 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
652 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
653 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
654 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
655 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
656 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
657 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
658 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
659 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
660 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
661 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
662 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
663 0 : if (present(elph_tr_ds)) then
664 0 : do iFSband=1,elph_ds%ngkkband !FS bands
665 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
666 0 : do icomp = 1, 3
667 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
668 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
669 0 : do jcomp = 1, 3
670 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
671 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
672 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
673 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
674 : end do
675 : end do
676 : end do
677 : end if
678 :
679 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
680 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
681 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
682 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
683 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
684 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
685 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
686 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
687 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
688 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
689 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
690 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
691 0 : if (present(elph_tr_ds)) then
692 0 : do iFSband=1,elph_ds%ngkkband !FS bands
693 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
694 0 : do icomp = 1, 3
695 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
696 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
697 0 : do jcomp = 1, 3
698 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
699 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
700 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
701 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
702 : end do
703 : end do
704 : end do
705 : end if
706 :
707 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
708 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
709 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
710 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
711 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
712 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
713 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
714 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
715 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
716 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
717 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
718 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
719 0 : if (present(elph_tr_ds)) then
720 0 : do iFSband=1,elph_ds%ngkkband !FS bands
721 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
722 0 : do icomp = 1, 3
723 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
724 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
725 0 : do jcomp = 1, 3
726 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
727 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
728 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
729 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
730 : end do
731 : end do
732 : end do
733 : end if
734 :
735 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
736 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
737 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
738 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
739 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
740 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
741 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
742 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
743 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
744 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
745 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
746 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
747 0 : if (present(elph_tr_ds)) then
748 0 : do iFSband=1,elph_ds%ngkkband !FS bands
749 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
750 0 : do icomp = 1, 3
751 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
752 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
753 0 : do jcomp = 1, 3
754 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
755 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
756 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
757 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
758 : end do
759 : end do
760 : end do
761 : end if
762 :
763 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
764 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
765 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
766 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
767 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
768 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
769 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
770 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
771 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
772 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
773 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
774 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
775 0 : if (present(elph_tr_ds)) then
776 0 : do iFSband=1,elph_ds%ngkkband !FS bands
777 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
778 0 : do icomp = 1, 3
779 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
780 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
781 0 : do jcomp = 1, 3
782 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
783 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
784 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
785 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
786 : end do
787 : end do
788 : end do
789 : end if
790 :
791 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
792 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
793 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
794 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
795 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
796 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
797 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
798 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
799 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
800 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
801 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
802 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
803 0 : if (present(elph_tr_ds)) then
804 0 : do iFSband=1,elph_ds%ngkkband !FS bands
805 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
806 0 : do icomp = 1, 3
807 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
808 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
809 0 : do jcomp = 1, 3
810 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
811 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
812 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
813 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
814 : end do
815 : end do
816 : end do
817 : end if
818 :
819 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
820 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
821 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
822 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
823 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
824 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
825 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
826 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
827 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
828 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
829 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
830 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
831 0 : if (present(elph_tr_ds)) then
832 0 : do iFSband=1,elph_ds%ngkkband !FS bands
833 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
834 0 : do icomp = 1, 3
835 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
836 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
837 0 : do jcomp = 1, 3
838 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
839 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
840 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
841 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
842 : end do
843 : end do
844 : end do
845 : end if
846 :
847 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
848 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
849 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
850 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
851 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
852 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
853 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
854 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
855 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
856 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
857 : elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
858 0 : & 0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
859 0 : if (present(elph_tr_ds)) then
860 0 : do iFSband=1,elph_ds%ngkkband !FS bands
861 0 : iband=iFSband+elph_ds%minFSband-1 ! full bands
862 0 : do icomp = 1, 3
863 : elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
864 0 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
865 0 : do jcomp = 1, 3
866 : elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
867 : & elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
868 : & 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
869 0 : & elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
870 : end do
871 : end do
872 : end do
873 : end if
874 : end if
875 : end do
876 : end do
877 : end do
878 :
879 : !bxu, divide by nscale^3 to be consistent with the normalization of kpt_phon
880 0 : elph_ds%k_phon%wtk = elph_ds%k_phon%wtk/nscale1/nscale2/nscale3
881 0 : if (present(elph_tr_ds)) then
882 0 : elph_ds%k_phon%velocwtk = elph_ds%k_phon%velocwtk/nscale1/nscale2/nscale3
883 0 : elph_ds%k_phon%vvelocwtk = elph_ds%k_phon%vvelocwtk/nscale1/nscale2/nscale3
884 : end if
885 :
886 0 : end subroutine d2c_weights
887 : !!***
888 :
889 : !!****f* ABINIT/d2c_wtq
890 : !! NAME
891 : !! d2c_wtq
892 : !!
893 : !! FUNCTION
894 : !! This routine calculates the integration weights on a coarse k-grid
895 : !! using the integration weights from a denser k-gird. The weights of
896 : !! the extra k points that being shared are evenly distributed.
897 : !!
898 : !! INPUTS
899 : !! elph_ds%k_fine%nkpt = number of fine q-points
900 : !! elph_ds%k_fine%wtq = integration weights of the fine q-grid
901 : !! elph_ds%k_phon%nkpt = number of coarse q-points
902 : !!
903 : !! OUTPUT
904 : !! elph_ds%k_phon%wtq = integration weights of the coarse k-grid
905 : !!
906 : !! SOURCE
907 :
908 0 : subroutine d2c_wtq(elph_ds)
909 :
910 : !Arguments ------------------------------------
911 : type(elph_type),intent(inout) :: elph_ds
912 :
913 : !Local variables-------------------------------
914 : integer :: ii, jj, kk
915 : integer :: ikpt, jkpt, kkpt
916 : integer :: iikpt, jjkpt, kkkpt
917 : integer :: ikpt_fine, ikpt_phon
918 : integer :: nkpt_fine1, nkpt_phon1
919 : integer :: nkpt_fine2, nkpt_phon2
920 : integer :: nkpt_fine3, nkpt_phon3
921 : integer :: nscale1, nscale2, nscale3
922 :
923 : ! *************************************************************************
924 0 : nkpt_phon1 = elph_ds%kptrlatt(1,1)
925 0 : nkpt_phon2 = elph_ds%kptrlatt(2,2)
926 0 : nkpt_phon3 = elph_ds%kptrlatt(3,3)
927 0 : nkpt_fine1 = elph_ds%kptrlatt_fine(1,1)
928 0 : nkpt_fine2 = elph_ds%kptrlatt_fine(2,2)
929 0 : nkpt_fine3 = elph_ds%kptrlatt_fine(3,3)
930 0 : nscale1 = dble(nkpt_fine1/nkpt_phon1)
931 0 : nscale2 = dble(nkpt_fine2/nkpt_phon2)
932 0 : nscale3 = dble(nkpt_fine3/nkpt_phon3)
933 : if (abs(INT(nscale1)-nscale1) > 0.01) then
934 : ABI_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
935 : end if
936 : if (abs(INT(nscale2)-nscale2) > 0.01) then
937 : ABI_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
938 : end if
939 : if (abs(INT(nscale3)-nscale3) > 0.01) then
940 : ABI_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
941 : end if
942 0 : nscale1 = INT(nscale1)
943 0 : nscale2 = INT(nscale2)
944 0 : nscale3 = INT(nscale3)
945 :
946 : !bxu, get wtq of coarse grid from fine grid
947 0 : elph_ds%k_phon%wtq = zero
948 :
949 0 : do ikpt = 1, nkpt_phon1
950 0 : do jkpt = 1, nkpt_phon2
951 0 : do kkpt = 1, nkpt_phon3
952 0 : ikpt_phon = kkpt + (jkpt-1)*nkpt_phon3 + (ikpt-1)*nkpt_phon2*nkpt_phon3
953 : ! inside the paralellepipe
954 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
955 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
956 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
957 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
958 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
959 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
960 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
961 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
962 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
963 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
964 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
965 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
966 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
967 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
968 0 : & elph_ds%k_fine%wtq(:,ikpt_fine,:)
969 : end do
970 : end do
971 : end do
972 : ! on the 6 faces
973 0 : if (MOD(nscale3,2) == 0) then ! when nscale3 is an even number
974 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
975 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
976 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
977 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
978 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
979 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
980 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
981 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
982 :
983 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
984 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
985 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
986 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
987 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
988 0 : & 0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
989 :
990 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
991 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
992 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
993 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
994 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
995 0 : & 0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
996 : end do
997 : end do
998 : end if
999 0 : if (MOD(nscale2,2) == 0) then ! when nscale2 is an even number
1000 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
1001 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
1002 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
1003 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
1004 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1005 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1006 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1007 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1008 :
1009 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1010 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1011 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1012 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1013 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1014 0 : & 0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1015 :
1016 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1017 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1018 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1019 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1020 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1021 0 : & 0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1022 : end do
1023 : end do
1024 : end if
1025 0 : if (MOD(nscale1,2) == 0) then ! when nscale1 is an even number
1026 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
1027 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
1028 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
1029 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
1030 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1031 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1032 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1033 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1034 :
1035 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1036 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1037 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1038 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1039 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1040 0 : & 0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1041 :
1042 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1043 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1044 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1045 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1046 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1047 0 : & 0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1048 : end do
1049 : end do
1050 : ! on the 12 sides
1051 : end if
1052 0 : if (MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
1053 0 : do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
1054 0 : iikpt = 1 + (ikpt-1)*nscale1 + ii
1055 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1056 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1057 :
1058 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1059 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1060 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1061 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1062 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1063 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1064 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1065 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1066 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1067 :
1068 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1069 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1070 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1071 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1072 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1073 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1074 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1075 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1076 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1077 :
1078 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1079 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1080 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1081 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1082 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1083 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1084 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1085 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1086 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1087 :
1088 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1089 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1090 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1091 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1092 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1093 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1094 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1095 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1096 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1097 : end do
1098 : end if
1099 0 : if (MOD(nscale1,2) == 0 .and. MOD(nscale3,2) == 0) then
1100 0 : do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
1101 0 : jjkpt = 1 + (jkpt-1)*nscale2 + jj
1102 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1103 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1104 :
1105 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1106 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1107 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1108 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1109 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1110 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1111 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1112 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1113 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1114 :
1115 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1116 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1117 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1118 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1119 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1120 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1121 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1122 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1123 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1124 :
1125 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1126 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1127 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1128 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1129 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1130 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1131 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1132 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1133 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1134 :
1135 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1136 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1137 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1138 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1139 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1140 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1141 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1142 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1143 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1144 : end do
1145 : end if
1146 0 : if (MOD(nscale2,2) == 0 .and. MOD(nscale1,2) == 0) then
1147 0 : do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
1148 0 : kkkpt = 1 + (kkpt-1)*nscale3 + kk
1149 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1150 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1151 :
1152 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1153 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1154 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1155 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1156 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1157 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1158 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1159 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1160 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1161 :
1162 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1163 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1164 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1165 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1166 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1167 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1168 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1169 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1170 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1171 :
1172 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1173 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1174 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1175 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1176 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1177 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1178 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1179 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1180 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1181 :
1182 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1183 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1184 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1185 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1186 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1187 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1188 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1189 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1190 0 : & 0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1191 : end do
1192 : ! on the 8 corners
1193 : end if
1194 0 : if (MOD(nscale1,2) == 0 .and. MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
1195 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1196 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1197 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1198 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1199 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1200 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1201 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1202 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1203 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1204 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1205 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1206 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1207 :
1208 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1209 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1210 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1211 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1212 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1213 0 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1214 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1215 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1216 0 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1217 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1218 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1219 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1220 :
1221 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1222 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1223 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1224 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1225 0 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1226 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1227 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1228 0 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1229 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1230 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1231 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1232 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1233 :
1234 0 : iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1235 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1236 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1237 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1238 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1239 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1240 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1241 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1242 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1243 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1244 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1245 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1246 :
1247 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1248 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1249 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1250 0 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1251 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1252 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1253 0 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1254 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1255 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1256 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1257 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1258 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1259 :
1260 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1261 0 : jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1262 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1263 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1264 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1265 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1266 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1267 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1268 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1269 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1270 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1271 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1272 :
1273 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1274 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1275 0 : kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1276 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1277 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1278 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1279 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1280 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1281 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1282 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1283 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1284 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1285 :
1286 0 : iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1287 0 : jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1288 0 : kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1289 : if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1290 : if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1291 : if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1292 : if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1293 : if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1294 : if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1295 0 : ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1296 : elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1297 0 : & 0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1298 : end if
1299 : end do
1300 : end do
1301 : end do
1302 :
1303 : !bxu, divide by nscale^3 to be consistent with the normalization of kpt_phon
1304 0 : elph_ds%k_phon%wtq = elph_ds%k_phon%wtq/nscale1/nscale2/nscale3
1305 :
1306 0 : end subroutine d2c_wtq
1307 : !!***
1308 :
1309 : !!****f* ABINIT/ep_el_weights
1310 : !!
1311 : !! NAME
1312 : !! ep_el_weights
1313 : !!
1314 : !! FUNCTION
1315 : !! This routine calculates the Fermi Surface integration weights
1316 : !! for the electron phonon routines, by different methods
1317 : !!
1318 : !! 1) Gaussian smearing
1319 : !! 2) Tetrahedron method
1320 : !! 3) Window in bands for all k-points
1321 : !! 4) Fermi Dirac smearing, follows gaussian with a different smearing function
1322 : !!
1323 : !! INPUTS
1324 : !! ep_b_min = minimal band to include in FS window integration
1325 : !! ep_b_max = maximal band to include in FS window integration
1326 : !! eigenGS = Ground State eigenvalues
1327 : !! elphsmear = smearing width for Gaussian method
1328 : !! fermie = Fermi level
1329 : !! gprimd = Reciprocal lattice vectors (dimensionful)
1330 : !! irredtoGS = mapping of elph k-points to ground state grid
1331 : !! kptrlatt = k-point grid vectors (if divided by determinant of present matrix)
1332 : !! max_occ = maximal occupancy for a band
1333 : !! minFSband = minimal band included for Fermi Surface integration in Gaussian and Tetrahedron cases
1334 : !! nFSband = number of bands in FS integration
1335 : !! nsppol = number of spin polarizations
1336 : !! telphint = option for FS integration:
1337 : !! 0 Tetrahedron method
1338 : !! 1 Gaussian smearing
1339 : !! 2 Window in bands for all k-points
1340 : !! 3 Fermi Dirac smearing
1341 : !! k_obj%nkpt = number of FS k-points
1342 : !! k_obj%kpt = FS k-points
1343 : !! k_obj%full2irr = mapping of FS k-points from full grid to irred points
1344 : !! k_obj%full2full = mapping of FS k-points in full grid under symops
1345 : !!
1346 : !! OUTPUT
1347 : !!
1348 : !! TODO
1349 : !! weights should be recalculated on-the-fly! The present implementation is not flexible!
1350 : !!
1351 : !! SOURCE
1352 :
1353 0 : subroutine ep_el_weights(ep_b_min, ep_b_max, eigenGS, elphsmear, enemin, enemax, nene, gprimd, &
1354 0 : & irredtoGS, kptrlatt, max_occ, minFSband, nband, nFSband, nsppol, telphint, k_obj, tmp_wtk)
1355 :
1356 : !Arguments ------------------------------------
1357 : !scalars
1358 : type(elph_kgrid_type), intent(in) :: k_obj
1359 : integer, intent(in) :: ep_b_min
1360 : integer, intent(in) :: ep_b_max
1361 : integer,intent(in) :: nband,nene
1362 : real(dp), intent(in) :: elphsmear
1363 : real(dp), intent(in) :: enemin,enemax
1364 : real(dp), intent(in) :: gprimd(3,3)
1365 : integer, intent(in) :: kptrlatt(3,3)
1366 : real(dp), intent(in) :: max_occ
1367 : integer, intent(in) :: minFSband
1368 : integer, intent(in) :: nFSband
1369 : integer, intent(in) :: nsppol
1370 : integer, intent(in) :: telphint
1371 :
1372 : ! arrays
1373 : real(dp), intent(in) :: eigenGS(nband,k_obj%nkptirr,nsppol)
1374 : real(dp), intent(out) :: tmp_wtk(nFSband,k_obj%nkpt,nsppol,nene)
1375 : integer, intent(in) :: irredtoGS(k_obj%nkptirr)
1376 :
1377 : !Local variables-------------------------------
1378 : !scalars
1379 : integer,parameter :: bcorr0=0
1380 : integer :: ikpt, ikptgs, ib1, iband
1381 : integer :: ierr, ie, isppol
1382 : real(dp) :: deltaene, rcvol, fermie
1383 : real(dp) :: smdeltaprefactor, smdeltafactor, xx
1384 :
1385 : ! arrays
1386 : real(dp) :: rlatt(3,3), klatt(3,3)
1387 0 : real(dp), allocatable :: tmp_eigen(:), tweight(:,:), dtweightde(:,:)
1388 : character (len=500) :: message
1389 : character (len=80) :: errstr
1390 0 : type(t_tetrahedron) :: tetrahedra
1391 : !type(htetra_t) :: tetrahedra
1392 :
1393 : ! *************************************************************************
1394 :
1395 : ! Initialize tmp_wtk with zeros
1396 0 : tmp_wtk = zero
1397 :
1398 : !write(std_out,*) 'ep_el : nkpt ', k_obj%nkpt
1399 : !===================================
1400 : !Set up integration weights for FS
1401 : !===================================
1402 0 : deltaene = (enemax-enemin)/dble(nene-1)
1403 :
1404 0 : if (telphint == 0) then
1405 :
1406 : ! =========================
1407 : ! Tetrahedron integration
1408 : ! =========================
1409 :
1410 0 : rlatt(:,:) = kptrlatt(:,:)
1411 0 : call matr3inv(rlatt,klatt)
1412 :
1413 : call init_tetra(k_obj%full2full(1,1,:), gprimd,klatt,k_obj%kpt, k_obj%nkpt,&
1414 0 : tetrahedra, ierr, errstr, xmpi_comm_self)
1415 : !call htetra_init(tetra, k_obj%full2full(1,1,:), gprimd, klatt, k_obj%kpt, k_obj%nkpt, &
1416 : ! k_obk%nkptirr, nkpt_ibz, ierr, errstr, xmpi_comm_self)
1417 :
1418 0 : ABI_CHECK(ierr==0,errstr)
1419 :
1420 : rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
1421 : & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
1422 0 : & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
1423 :
1424 : ! fix small window around fermie for tetrahedron weight calculation
1425 0 : deltaene = (enemax-enemin)/dble(nene-1)
1426 :
1427 0 : ABI_MALLOC(tmp_eigen,(k_obj%nkpt))
1428 0 : ABI_MALLOC(tweight,(k_obj%nkpt,nene))
1429 0 : ABI_MALLOC(dtweightde,(k_obj%nkpt,nene))
1430 :
1431 0 : do iband = 1,nFSband
1432 : ! for each spin pol
1433 0 : do isppol=1,nsppol
1434 : ! For this band get its contribution
1435 0 : tmp_eigen(:) = zero
1436 0 : do ikpt=1,k_obj%nkpt
1437 0 : ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1438 0 : tmp_eigen(ikpt) = eigenGS(minFSband+iband-1,ikptgs,isppol)
1439 : end do
1440 : ! calculate general integration weights at each irred kpoint
1441 : ! as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
1442 : call get_tetra_weight(tmp_eigen,enemin,enemax,&
1443 : & max_occ,nene,k_obj%nkpt,tetrahedra,bcorr0,&
1444 0 : & tweight,dtweightde,xmpi_comm_self)
1445 :
1446 0 : tmp_wtk(iband,:,isppol,:) = dtweightde(:,:)*k_obj%nkpt
1447 : end do
1448 : end do
1449 0 : ABI_FREE(tmp_eigen)
1450 0 : ABI_FREE(tweight)
1451 0 : ABI_FREE(dtweightde)
1452 :
1453 0 : call destroy_tetra(tetrahedra)
1454 : !call tetrahedra%free()
1455 :
1456 0 : else if (telphint == 1) then
1457 :
1458 : ! ==============================================================
1459 : ! Gaussian or integration:
1460 : ! Each kpt contributes a gaussian of integrated weight 1
1461 : ! for each band. The gaussian being centered at the Fermi level.
1462 : ! ===============================================================
1463 :
1464 : ! took out factor 1/k_obj%nkpt which intervenes only at integration time
1465 :
1466 : ! MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1467 :
1468 : ! gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1469 0 : smdeltaprefactor = max_occ*sqrt(piinv)/elphsmear
1470 0 : smdeltafactor = one/elphsmear
1471 :
1472 : ! SPPOL loop on isppol as well to get 2 sets of weights
1473 0 : do isppol=1,nsppol
1474 : fermie = enemin
1475 0 : do ie = 1, nene
1476 0 : fermie = fermie + deltaene
1477 : ! fine grid
1478 0 : do ikpt=1, k_obj%nkpt
1479 0 : ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1480 0 : do ib1=1,nFSband
1481 0 : xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1482 0 : if (abs(xx) < 40._dp) then
1483 0 : tmp_wtk(ib1,ikpt,isppol,ie) = exp(-xx*xx)*smdeltaprefactor
1484 : end if
1485 : end do
1486 : end do
1487 : end do
1488 : end do
1489 :
1490 :
1491 0 : else if (telphint == 2) then ! range of bands occupied
1492 :
1493 : ! SPPOL eventually be able to specify bands for up and down separately
1494 : fermie = enemin
1495 0 : do ie = 1, nene
1496 0 : fermie = fermie + deltaene
1497 0 : do ikpt=1,k_obj%nkpt
1498 0 : do ib1=ep_b_min, ep_b_max
1499 : ! for the moment both spin channels same
1500 0 : tmp_wtk(ib1,ikpt,:,ie) = max_occ
1501 : end do
1502 : end do
1503 : end do
1504 :
1505 0 : write(std_out,*) ' ep_el_weights : DOS is calculated from states in bands ',ep_b_min,' to ',ep_b_max
1506 :
1507 0 : else if (telphint == 3) then
1508 :
1509 : ! ==============================================================
1510 : ! Fermi Dirac integration:
1511 : ! Each kpt contributes a Fermi Dirac smearing function of integrated weight 1
1512 : ! for each band. The function being centered at the Fermi level.
1513 : ! ===============================================================
1514 :
1515 : ! took out factor 1/k_obj%nkpt which intervenes only at integration time
1516 :
1517 : ! MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1518 :
1519 : ! gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1520 0 : smdeltaprefactor = half*max_occ/elphsmear
1521 0 : smdeltafactor = one/elphsmear
1522 :
1523 : ! SPPOL loop on isppol as well to get 2 sets of weights
1524 0 : do isppol=1,nsppol
1525 : fermie = enemin
1526 0 : do ie = 1, nene
1527 0 : fermie = fermie + deltaene
1528 : ! fine grid
1529 0 : do ikpt=1, k_obj%nkpt
1530 0 : ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1531 0 : do ib1=1,nFSband
1532 0 : xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1533 0 : tmp_wtk(ib1,ikpt,isppol,ie) = smdeltaprefactor / (one + cosh(xx))
1534 : end do
1535 : end do
1536 : end do
1537 : end do
1538 :
1539 :
1540 : else
1541 0 : write (message,'(a,i0)')" telphint should be between 0 and 3, found: ",telphint
1542 0 : ABI_BUG(message)
1543 : end if ! if telphint
1544 :
1545 0 : end subroutine ep_el_weights
1546 : !!***
1547 :
1548 : !!****f* ABINIT/ep_fs_weights
1549 : !!
1550 : !! NAME
1551 : !! ep_fs_weights
1552 : !!
1553 : !! FUNCTION
1554 : !! This routine calculates the Fermi Surface integration weights
1555 : !! for the electron phonon routines, by different methods
1556 : !! 1) Gaussian smearing
1557 : !! 2) Tetrahedron method
1558 : !! 3) Window in bands for all k-points
1559 : !! 4) Fermi Dirac smearing, follows gaussian with a different smearing function
1560 : !!
1561 : !! INPUTS
1562 : !! ep_b_min = minimal band to include in FS window integration
1563 : !! ep_b_max = maximal band to include in FS window integration
1564 : !! eigenGS = Ground State eigenvalues
1565 : !! elphsmear = smearing width for Gaussian method
1566 : !! fermie = Fermi level
1567 : !! gprimd = Reciprocal lattice vectors (dimensionful)
1568 : !! irredtoGS = mapping of elph k-points to ground state grid
1569 : !! kptrlatt = k-point grid vectors (if divided by determinant of present matrix)
1570 : !! max_occ = maximal occupancy for a band
1571 : !! minFSband = minimal band included for Fermi Surface integration in Gaussian and Tetrahedron cases
1572 : !! nFSband = number of bands in FS integration
1573 : !! nsppol = number of spin polarizations
1574 : !! telphint = option for FS integration:
1575 : !! 0 Tetrahedron method
1576 : !! 1 Gaussian smearing
1577 : !! 2 Window in bands for all k-points
1578 : !! 3 Fermi Dirac smearing
1579 : !! k_obj%nkpt = number of FS k-points
1580 : !! k_obj%kpt = FS k-points
1581 : !! k_obj%full2irr = mapping of FS k-points from full grid to irred points
1582 : !! k_obj%full2full = mapping of FS k-points in full grid under symops
1583 : !!
1584 : !! OUTPUT
1585 : !! k_obj%wtk = integration weights
1586 : !!
1587 : !! TODO
1588 : !! weights should be recalculated on-the-fly! The present implementation is not flexible!
1589 : !!
1590 : !! SOURCE
1591 :
1592 15 : subroutine ep_fs_weights(ep_b_min, ep_b_max, eigenGS, elphsmear, fermie, gprimd, &
1593 15 : & irredtoGS, kptrlatt, max_occ, minFSband, nband, nFSband, nsppol, telphint, k_obj)
1594 :
1595 : !Arguments ------------------------------------
1596 : !scalars
1597 : type(elph_kgrid_type), intent(inout) :: k_obj
1598 : integer, intent(in) :: ep_b_min
1599 : integer, intent(in) :: ep_b_max
1600 : integer,intent(in) :: nband
1601 : real(dp), intent(in) :: elphsmear
1602 : real(dp), intent(in) :: fermie
1603 : real(dp), intent(in) :: gprimd(3,3)
1604 : integer, intent(in) :: kptrlatt(3,3)
1605 : real(dp), intent(in) :: max_occ
1606 : integer, intent(in) :: minFSband
1607 : integer, intent(in) :: nFSband
1608 : integer, intent(in) :: nsppol
1609 : integer, intent(in) :: telphint
1610 :
1611 : ! arrays
1612 : real(dp), intent(in) :: eigenGS(nband,k_obj%nkptirr,nsppol)
1613 : integer, intent(in) :: irredtoGS(k_obj%nkptirr)
1614 :
1615 : !Local variables-------------------------------
1616 : !scalars
1617 : integer,parameter :: bcorr0=0
1618 : integer :: ikpt, ikptgs, ib1, isppol, iband
1619 : integer :: nene, ifermi
1620 : integer :: ierr
1621 :
1622 : real(dp) :: enemin, enemax, deltaene, rcvol
1623 : real(dp) :: smdeltaprefactor, smdeltafactor, xx
1624 :
1625 : ! arrays
1626 : real(dp) :: rlatt(3,3), klatt(3,3)
1627 15 : real(dp), allocatable :: tmp_eigen(:), tweight(:,:), dtweightde(:,:)
1628 :
1629 : character (len=500) :: message
1630 : character (len=80) :: errstr
1631 :
1632 15 : type(t_tetrahedron) :: tetrahedra
1633 :
1634 : ! *************************************************************************
1635 :
1636 15 : write(std_out,*) 'ep_fs : nkpt ', k_obj%nkpt
1637 15 : write(message, '(a)' ) '- ep_fs_weights 1 = '
1638 15 : call wrtout(std_out,message,'PERS')
1639 :
1640 : !===================================
1641 : !Set up integration weights for FS
1642 : !===================================
1643 :
1644 15 : if (telphint == 0) then
1645 :
1646 : ! =========================
1647 : ! Tetrahedron integration
1648 : ! =========================
1649 :
1650 26 : rlatt(:,:) = kptrlatt(:,:)
1651 2 : call matr3inv(rlatt,klatt)
1652 :
1653 : call init_tetra(k_obj%full2full(1,1,:), gprimd,klatt,k_obj%kpt, k_obj%nkpt,&
1654 434 : & tetrahedra, ierr, errstr, xmpi_comm_self)
1655 2 : ABI_CHECK(ierr==0,errstr)
1656 :
1657 : rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
1658 : & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
1659 2 : & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
1660 :
1661 : ! just do weights at FS
1662 2 : nene = 100
1663 :
1664 : ! fix small window around fermie for tetrahedron weight calculation
1665 2 : deltaene = 2*elphsmear/dble(nene-1)
1666 2 : ifermi = int(nene/2)
1667 2 : enemin = fermie - dble(ifermi-1)*deltaene
1668 2 : enemax = enemin + dble(nene-1)*deltaene
1669 :
1670 6 : ABI_MALLOC(tmp_eigen,(k_obj%nkpt))
1671 8 : ABI_MALLOC(tweight,(k_obj%nkpt,nene))
1672 6 : ABI_MALLOC(dtweightde,(k_obj%nkpt,nene))
1673 :
1674 14 : do iband = 1,nFSband
1675 : ! for each spin pol
1676 26 : do isppol=1,nsppol
1677 : ! For this band get its contribution
1678 2604 : tmp_eigen(:) = zero
1679 2604 : do ikpt=1,k_obj%nkpt
1680 2592 : ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1681 2604 : tmp_eigen(ikpt) = eigenGS(minFSband+iband-1,ikptgs,isppol)
1682 : end do
1683 : ! calculate general integration weights at each irred kpoint
1684 : ! as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
1685 : call get_tetra_weight(tmp_eigen,enemin,enemax,&
1686 : & max_occ,nene,k_obj%nkpt,tetrahedra,bcorr0,&
1687 12 : & tweight,dtweightde,xmpi_comm_self)
1688 :
1689 2616 : k_obj%wtk(iband,:,isppol) = dtweightde(:,ifermi)*k_obj%nkpt
1690 : end do
1691 :
1692 : end do
1693 2 : ABI_FREE(tmp_eigen)
1694 2 : ABI_FREE(tweight)
1695 2 : ABI_FREE(dtweightde)
1696 :
1697 2 : call destroy_tetra(tetrahedra)
1698 :
1699 13 : else if (telphint == 1) then
1700 :
1701 : ! ==============================================================
1702 : ! Gaussian or integration:
1703 : ! Each kpt contributes a gaussian of integrated weight 1
1704 : ! for each band. The gaussian being centered at the Fermi level.
1705 : ! ===============================================================
1706 :
1707 : ! took out factor 1/k_obj%nkpt which intervenes only at integration time
1708 :
1709 : ! MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1710 :
1711 : ! gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1712 12 : smdeltaprefactor = max_occ*sqrt(piinv)/elphsmear
1713 12 : smdeltafactor = one/elphsmear
1714 :
1715 6241 : k_obj%wtk = zero
1716 : ! SPPOL loop on isppol as well to get 2 sets of weights
1717 25 : do isppol=1,nsppol
1718 : ! fine grid
1719 849 : do ikpt=1, k_obj%nkpt
1720 824 : ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1721 6229 : do ib1=1,nFSband
1722 5392 : xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1723 6216 : if (abs(xx) < 40._dp) then
1724 3921 : k_obj%wtk(ib1,ikpt,isppol) = exp(-xx*xx)*smdeltaprefactor
1725 : end if
1726 : end do
1727 : end do
1728 : end do
1729 :
1730 :
1731 1 : else if (telphint == 2) then ! range of bands occupied
1732 :
1733 : ! SPPOL eventually be able to specify bands for up and down separately
1734 450 : k_obj%wtk = zero
1735 65 : do ikpt=1,k_obj%nkpt
1736 449 : do ib1=ep_b_min, ep_b_max
1737 : ! for the moment both spin channels same
1738 832 : k_obj%wtk(ib1,ikpt,:) = max_occ
1739 : end do
1740 : end do
1741 :
1742 1 : write(std_out,*) ' ep_fs_weights : DOS is calculated from states in bands ',ep_b_min,' to ',ep_b_max
1743 :
1744 0 : else if (telphint == 3) then
1745 :
1746 : ! ==============================================================
1747 : ! Fermi Dirac integration:
1748 : ! Each kpt contributes a Fermi Dirac smearing function of integrated weight 1
1749 : ! for each band. The function being centered at the Fermi level.
1750 : ! ===============================================================
1751 :
1752 : ! took out factor 1/k_obj%nkpt which intervenes only at integration time
1753 :
1754 : ! MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1755 :
1756 : ! gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1757 0 : smdeltaprefactor = half*max_occ/elphsmear
1758 0 : smdeltafactor = one/elphsmear
1759 :
1760 0 : k_obj%wtk = zero
1761 : ! SPPOL loop on isppol as well to get 2 sets of weights
1762 0 : do isppol=1,nsppol
1763 : ! fine grid
1764 0 : do ikpt=1, k_obj%nkpt
1765 0 : ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1766 0 : do ib1=1,nFSband
1767 0 : xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1768 0 : k_obj%wtk(ib1,ikpt,isppol) = smdeltaprefactor / (one + cosh(xx))
1769 : end do
1770 : end do
1771 : end do
1772 :
1773 :
1774 : else
1775 0 : write (message,'(a,i0)')" telphint should be between 0 and 3, found: ",telphint
1776 0 : ABI_BUG(message)
1777 : end if ! if telphint
1778 :
1779 15 : end subroutine ep_fs_weights
1780 : !!***
1781 :
1782 : !!****f* ABINIT/ep_ph_weights
1783 : !!
1784 : !! NAME
1785 : !! ep_ph_weights
1786 : !!
1787 : !! FUNCTION
1788 : !! This routine calculates the phonon integration weights
1789 : !! for the electron phonon routines, by different methods
1790 : !! 1) Gaussian smearing
1791 : !! 0) Tetrahedron method
1792 : !!
1793 : !! INPUTS
1794 : !! phfrq = phonon energies
1795 : !! elphsmear = smearing width for Gaussian method
1796 : !! omega = input phonon energy
1797 : !! gprimd = Reciprocal lattice vectors (dimensionful)
1798 : !! kptrlatt = k-point grid vectors (if divided by determinant of present matrix)
1799 : !! telphint = option for FS integration:
1800 : !! 0 Tetrahedron method
1801 : !! 1 Gaussian smearing
1802 : !! k_obj%nkpt = number of FS k-points
1803 : !! k_obj%kpt = FS k-points
1804 : !! k_obj%full2full = mapping of FS k-points in full grid under symops
1805 : !!
1806 : !! OUTPUT
1807 : !! tmp_wtq = integration weights
1808 : !!
1809 : !! TODO
1810 : !! weights should be recalculated on-the-fly! The present implementation is not flexible!
1811 : !!
1812 : !! SOURCE
1813 :
1814 15 : subroutine ep_ph_weights(phfrq,elphsmear,omega_min,omega_max,nomega,gprimd,kptrlatt,nbranch,telphint,k_obj,tmp_wtq)
1815 :
1816 : !Arguments ------------------------------------
1817 : !scalars
1818 : type(elph_kgrid_type), intent(inout) :: k_obj
1819 : integer,intent(in) :: nbranch
1820 : real(dp), intent(in) :: elphsmear
1821 : real(dp), intent(in) :: omega_min,omega_max
1822 : real(dp), intent(in) :: gprimd(3,3)
1823 : integer, intent(in) :: kptrlatt(3,3)
1824 : integer, intent(in) :: nomega
1825 : integer, intent(in) :: telphint
1826 :
1827 : ! arrays
1828 : real(dp), intent(in) :: phfrq(nbranch,k_obj%nkpt)
1829 : real(dp), intent(out) :: tmp_wtq(nbranch,k_obj%nkpt,nomega)
1830 :
1831 : !Local variables-------------------------------
1832 : !scalars
1833 : integer,parameter :: bcorr0=0
1834 : integer :: ikpt, ib1, ibranch
1835 : integer :: ierr, iomega
1836 : real(dp) :: rcvol, max_occ
1837 : real(dp) :: smdeltaprefactor, smdeltafactor, xx, gaussmaxarg
1838 : real(dp) :: domega,omega
1839 :
1840 : ! arrays
1841 : real(dp) :: rlatt(3,3), klatt(3,3)
1842 15 : real(dp), allocatable :: tweight(:,:), dtweightde(:,:)
1843 : character (len=80) :: errstr
1844 15 : type(t_tetrahedron) :: tetrahedra
1845 :
1846 : ! *************************************************************************
1847 :
1848 : !write(std_out,*) 'ep_ph : nqpt ', k_obj%nkpt
1849 : !===================================
1850 : !Set up integration weights for FS
1851 : !===================================
1852 15 : max_occ = one
1853 15 : gaussmaxarg = sqrt(-log(1.d-100))
1854 15 : domega = (omega_max - omega_min)/(nomega - 1)
1855 :
1856 15 : if (telphint == 0) then
1857 :
1858 : ! =========================
1859 : ! Tetrahedron integration
1860 : ! =========================
1861 :
1862 26 : rlatt(:,:) = kptrlatt(:,:)
1863 2 : call matr3inv(rlatt,klatt)
1864 :
1865 : call init_tetra(k_obj%full2full(1,1,:), gprimd,klatt,k_obj%kpt, k_obj%nkpt,&
1866 434 : & tetrahedra, ierr, errstr, xmpi_comm_self)
1867 2 : ABI_CHECK(ierr==0,errstr)
1868 :
1869 : rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
1870 : & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
1871 2 : & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
1872 :
1873 : ! do all the omega points for tetrahedron weight calculation
1874 :
1875 8 : ABI_MALLOC(tweight,(k_obj%nkpt,nomega))
1876 6 : ABI_MALLOC(dtweightde,(k_obj%nkpt,nomega))
1877 :
1878 8 : do ibranch = 1,nbranch
1879 : call get_tetra_weight(phfrq(ibranch,:),omega_min,omega_max,&
1880 : & max_occ,nomega,k_obj%nkpt,tetrahedra,bcorr0,&
1881 1302 : & tweight,dtweightde,xmpi_comm_self)
1882 :
1883 522110 : tmp_wtq(ibranch,:,:) = dtweightde(:,:)*k_obj%nkpt
1884 : end do
1885 2 : ABI_FREE(tweight)
1886 2 : ABI_FREE(dtweightde)
1887 :
1888 2 : call destroy_tetra(tetrahedra)
1889 :
1890 13 : else if (telphint == 1) then
1891 :
1892 : ! ==============================================================
1893 : ! Gaussian or integration:
1894 : ! Each kpt contributes a gaussian of integrated weight 1
1895 : ! for each branch. The gaussian being centered at the input energy
1896 : ! ===============================================================
1897 :
1898 : ! took out factor 1/k_obj%nkpt which intervenes only at integration time
1899 :
1900 : ! gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1901 12 : smdeltaprefactor = max_occ*sqrt(piinv)/elphsmear
1902 12 : smdeltafactor = one/elphsmear
1903 :
1904 1486920 : tmp_wtq = zero
1905 : omega = omega_min
1906 4824 : do iomega = 1, nomega
1907 4812 : omega = omega + domega
1908 332040 : do ikpt=1, k_obj%nkpt
1909 1486908 : do ib1=1,nbranch
1910 1154880 : xx = smdeltafactor*(phfrq(ib1,ikpt)-omega)
1911 1482096 : if (abs(xx) < gaussmaxarg) then
1912 302668 : tmp_wtq(ib1,ikpt,iomega) = exp(-xx*xx)*smdeltaprefactor
1913 : end if
1914 : end do
1915 : end do
1916 : end do
1917 : end if ! if telphint
1918 :
1919 15 : end subroutine ep_ph_weights
1920 : !!***
1921 :
1922 : end module m_epweights
1923 : !!***
|
