Branch data 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 : : !!***
|
