TABLE OF CONTENTS


ABINIT/d2c_weights [ Functions ]

[ Top ] [ Functions ]

NAME

  d2c_weights

FUNCTION

  This routine calculates the integration weights on a coarse k-grid
  using the integration weights from a denser k-gird. The weights of
  the extra k points that being shared are evenly distributed. It also
  condenses the velocity*wtk and velcity^2*wtk.

INPUTS

  elph_ds%k_fine%nkpt = number of fine FS k-points
  elph_ds%k_fine%wtk = integration weights of the fine FS k-grid
  elph_ds%k_phon%nkpt = number of coarse FS k-points
  elph_tr_ds%el_veloc = electronic velocities from the fine k-grid

OUTPUT

  elph_ds%k_phon%wtk = integration weights of the coarse FS k-grid
  elph_ds%k_phon%velocwtk = velocity time integration weights of the coarse FS k-grid
  elph_ds%k_phon%vvelocwtk = velocity^2 time integration weights of the coarse FS k-grid

PARENTS

      elphon,get_nv_fs_en

CHILDREN

SOURCE

 81 subroutine d2c_weights(elph_ds,elph_tr_ds)
 82 
 83 
 84 !This section has been created automatically by the script Abilint (TD).
 85 !Do not modify the following lines by hand.
 86 #undef ABI_FUNC
 87 #define ABI_FUNC 'd2c_weights'
 88 !End of the abilint section
 89 
 90  implicit none
 91 
 92 !Arguments ------------------------------------
 93  type(elph_type),intent(inout) :: elph_ds
 94  type(elph_tr_type),intent(inout),optional :: elph_tr_ds
 95 
 96 !Local variables-------------------------------
 97  integer :: ii, jj, kk
 98  integer :: ikpt, jkpt, kkpt
 99  integer :: iikpt, jjkpt, kkkpt
100  integer :: icomp, jcomp
101  integer :: iFSband, iband
102  integer :: ikpt_fine, ikpt_phon
103  integer :: nkpt_fine1, nkpt_phon1
104  integer :: nkpt_fine2, nkpt_phon2
105  integer :: nkpt_fine3, nkpt_phon3
106  integer :: nscale1, nscale2, nscale3
107 
108 ! *************************************************************************
109  nkpt_phon1 = elph_ds%kptrlatt(1,1)
110  nkpt_phon2 = elph_ds%kptrlatt(2,2)
111  nkpt_phon3 = elph_ds%kptrlatt(3,3)
112  nkpt_fine1 = elph_ds%kptrlatt_fine(1,1)
113  nkpt_fine2 = elph_ds%kptrlatt_fine(2,2)
114  nkpt_fine3 = elph_ds%kptrlatt_fine(3,3)
115  nscale1 = dble(nkpt_fine1/nkpt_phon1)
116  nscale2 = dble(nkpt_fine2/nkpt_phon2)
117  nscale3 = dble(nkpt_fine3/nkpt_phon3)
118  if (abs(INT(nscale1)-nscale1) > 0.01) then
119    MSG_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
120  end if
121  if (abs(INT(nscale2)-nscale2) > 0.01) then
122    MSG_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
123  end if
124  if (abs(INT(nscale3)-nscale3) > 0.01) then
125    MSG_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
126  end if
127  nscale1 = INT(nscale1)
128  nscale2 = INT(nscale2)
129  nscale3 = INT(nscale3)
130 
131 !bxu, get wtk of coarse grid from fine grid
132  elph_ds%k_phon%wtk = zero
133  if (present(elph_tr_ds)) then
134    elph_ds%k_phon%velocwtk = zero
135    elph_ds%k_phon%vvelocwtk = zero
136  end if
137 
138  do ikpt = 1, nkpt_phon1
139    do jkpt = 1, nkpt_phon2
140      do kkpt = 1, nkpt_phon3
141        ikpt_phon = kkpt + (jkpt-1)*nkpt_phon3 + (ikpt-1)*nkpt_phon2*nkpt_phon3
142 !      inside the paralellepipe
143        do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
144          do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
145            do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
146              iikpt = 1 + (ikpt-1)*nscale1 + ii
147              jjkpt = 1 + (jkpt-1)*nscale2 + jj
148              kkkpt = 1 + (kkpt-1)*nscale3 + kk
149              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
150              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
151              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
152              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
153              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
154              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
155              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
156              elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
157 &             elph_ds%k_fine%wtk(:,ikpt_fine,:)
158              if (present(elph_tr_ds)) then
159                do iFSband=1,elph_ds%ngkkband !FS bands
160                  iband=iFSband+elph_ds%minFSband-1 ! full bands
161                  do icomp = 1, 3
162                    elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
163 &                   elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
164                    do jcomp = 1, 3
165                      elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
166 &                     elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
167 &                     elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
168 &                     elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
169                    end do
170                  end do
171                end do
172              end if
173            end do
174          end do
175        end do
176 !      on the 6 faces
177        if (MOD(nscale3,2) == 0) then ! when nscale3 is an even number
178          do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
179            do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
180              iikpt = 1 + (ikpt-1)*nscale1 + ii
181              jjkpt = 1 + (jkpt-1)*nscale2 + jj
182              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
183              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
184              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
185              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
186 
187              kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
188              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
189              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
190              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
191              elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
192 &             0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
193              if (present(elph_tr_ds)) then
194                do iFSband=1,elph_ds%ngkkband !FS bands
195                  iband=iFSband+elph_ds%minFSband-1 ! full bands
196                  do icomp = 1, 3
197                    elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
198 &                   0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
199                    do jcomp = 1, 3
200                      elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
201 &                     elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
202 &                     0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
203 &                     elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
204                    end do
205                  end do
206                end do
207              end if
208 
209              kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
210              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
211              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
212              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
213              elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
214 &             0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
215              if (present(elph_tr_ds)) then
216                do iFSband=1,elph_ds%ngkkband !FS bands
217                  iband=iFSband+elph_ds%minFSband-1 ! full bands
218                  do icomp = 1, 3
219                    elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
220 &                   0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
221                    do jcomp = 1, 3
222                      elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
223 &                     elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
224 &                     0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
225 &                     elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
226                    end do
227                  end do
228                end do
229              end if
230            end do
231          end do
232        end if
233        if (MOD(nscale2,2) == 0) then ! when nscale2 is an even number
234          do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
235            do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
236              iikpt = 1 + (ikpt-1)*nscale1 + ii
237              kkkpt = 1 + (kkpt-1)*nscale3 + kk
238              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
239              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
240              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
241              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
242 
243              jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
244              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
245              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
246              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
247              elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
248 &             0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
249              if (present(elph_tr_ds)) then
250                do iFSband=1,elph_ds%ngkkband !FS bands
251                  iband=iFSband+elph_ds%minFSband-1 ! full bands
252                  do icomp = 1, 3
253                    elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
254 &                   0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
255                    do jcomp = 1, 3
256                      elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
257 &                     elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
258 &                     0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
259 &                     elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
260                    end do
261                  end do
262                end do
263              end if
264 
265              jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
266              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
267              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
268              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
269              elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
270 &             0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
271              if (present(elph_tr_ds)) then
272                do iFSband=1,elph_ds%ngkkband !FS bands
273                  iband=iFSband+elph_ds%minFSband-1 ! full bands
274                  do icomp = 1, 3
275                    elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
276 &                   0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
277                    do jcomp = 1, 3
278                      elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
279 &                     elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
280 &                     0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
281 &                     elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
282                    end do
283                  end do
284                end do
285              end if
286            end do
287          end do
288        end if
289        if (MOD(nscale1,2) == 0) then ! when nscale1 is an even number
290          do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
291            do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
292              kkkpt = 1 + (kkpt-1)*nscale3 + kk
293              jjkpt = 1 + (jkpt-1)*nscale2 + jj
294              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
295              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
296              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
297              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
298 
299              iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
300              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
301              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
302              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
303              elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
304 &             0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
305              if (present(elph_tr_ds)) then
306                do iFSband=1,elph_ds%ngkkband !FS bands
307                  iband=iFSband+elph_ds%minFSband-1 ! full bands
308                  do icomp = 1, 3
309                    elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
310 &                   0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
311                    do jcomp = 1, 3
312                      elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
313 &                     elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
314 &                     0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
315 &                     elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
316                    end do
317                  end do
318                end do
319              end if
320 
321              iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
322              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
323              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
324              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
325              elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
326 &             0.5_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
327              if (present(elph_tr_ds)) then
328                do iFSband=1,elph_ds%ngkkband !FS bands
329                  iband=iFSband+elph_ds%minFSband-1 ! full bands
330                  do icomp = 1, 3
331                    elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
332 &                   0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
333                    do jcomp = 1, 3
334                      elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
335 &                     elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
336 &                     0.5_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
337 &                     elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
338                    end do
339                  end do
340                end do
341              end if
342            end do
343          end do
344 !        on the 12 sides
345        end if
346        if (MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
347          do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
348            iikpt = 1 + (ikpt-1)*nscale1 + ii
349            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
350            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
351 
352            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
353            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
354            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
355            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
356            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
357            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
358            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
359            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
360 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
361            if (present(elph_tr_ds)) then
362              do iFSband=1,elph_ds%ngkkband !FS bands
363                iband=iFSband+elph_ds%minFSband-1 ! full bands
364                do icomp = 1, 3
365                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
366 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
367                  do jcomp = 1, 3
368                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
369 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
370 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
371 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
372                  end do
373                end do
374              end do
375            end if
376 
377            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
378            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
379            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
380            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
381            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
382            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
383            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
384            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
385 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
386            if (present(elph_tr_ds)) then
387              do iFSband=1,elph_ds%ngkkband !FS bands
388                iband=iFSband+elph_ds%minFSband-1 ! full bands
389                do icomp = 1, 3
390                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
391 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
392                  do jcomp = 1, 3
393                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
394 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
395 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
396 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
397                  end do
398                end do
399              end do
400            end if
401 
402            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
403            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
404            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
405            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
406            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
407            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
408            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
409            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
410 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
411            if (present(elph_tr_ds)) then
412              do iFSband=1,elph_ds%ngkkband !FS bands
413                iband=iFSband+elph_ds%minFSband-1 ! full bands
414                do icomp = 1, 3
415                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
416 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
417                  do jcomp = 1, 3
418                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
419 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
420 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
421 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
422                  end do
423                end do
424              end do
425            end if
426 
427            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
428            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
429            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
430            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
431            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
432            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
433            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
434            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
435 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
436            if (present(elph_tr_ds)) then
437              do iFSband=1,elph_ds%ngkkband !FS bands
438                iband=iFSband+elph_ds%minFSband-1 ! full bands
439                do icomp = 1, 3
440                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
441 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
442                  do jcomp = 1, 3
443                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
444 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
445 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
446 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
447                  end do
448                end do
449              end do
450            end if
451          end do
452        end if
453        if (MOD(nscale1,2) == 0 .and. MOD(nscale3,2) == 0) then
454          do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
455            jjkpt = 1 + (jkpt-1)*nscale2 + jj
456            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
457            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
458 
459            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
460            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
461            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
462            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
463            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
464            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
465            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
466            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
467 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
468            if (present(elph_tr_ds)) then
469              do iFSband=1,elph_ds%ngkkband !FS bands
470                iband=iFSband+elph_ds%minFSband-1 ! full bands
471                do icomp = 1, 3
472                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
473 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
474                  do jcomp = 1, 3
475                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
476 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
477 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
478 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
479                  end do
480                end do
481              end do
482            end if
483 
484            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
485            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
486            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
487            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
488            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
489            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
490            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
491            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
492 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
493            if (present(elph_tr_ds)) then
494              do iFSband=1,elph_ds%ngkkband !FS bands
495                iband=iFSband+elph_ds%minFSband-1 ! full bands
496                do icomp = 1, 3
497                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
498 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
499                  do jcomp = 1, 3
500                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
501 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
502 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
503 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
504                  end do
505                end do
506              end do
507            end if
508 
509            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
510            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
511            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
512            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
513            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
514            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
515            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
516            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
517 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
518            if (present(elph_tr_ds)) then
519              do iFSband=1,elph_ds%ngkkband !FS bands
520                iband=iFSband+elph_ds%minFSband-1 ! full bands
521                do icomp = 1, 3
522                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
523 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
524                  do jcomp = 1, 3
525                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
526 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
527 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
528 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
529                  end do
530                end do
531              end do
532            end if
533 
534            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
535            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
536            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
537            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
538            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
539            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
540            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
541            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
542 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
543            if (present(elph_tr_ds)) then
544              do iFSband=1,elph_ds%ngkkband !FS bands
545                iband=iFSband+elph_ds%minFSband-1 ! full bands
546                do icomp = 1, 3
547                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
548 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
549                  do jcomp = 1, 3
550                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
551 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
552 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
553 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
554                  end do
555                end do
556              end do
557            end if
558          end do
559        end if
560        if (MOD(nscale2,2) == 0 .and. MOD(nscale1,2) == 0) then
561          do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
562            kkkpt = 1 + (kkpt-1)*nscale3 + kk
563            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
564            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
565 
566            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
567            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
568            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
569            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
570            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
571            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
572            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
573            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
574 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
575            if (present(elph_tr_ds)) then
576              do iFSband=1,elph_ds%ngkkband !FS bands
577                iband=iFSband+elph_ds%minFSband-1 ! full bands
578                do icomp = 1, 3
579                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
580 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
581                  do jcomp = 1, 3
582                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
583 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
584 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
585 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
586                  end do
587                end do
588              end do
589            end if
590 
591            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
592            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
593            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
594            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
595            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
596            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
597            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
598            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
599 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
600            if (present(elph_tr_ds)) then
601              do iFSband=1,elph_ds%ngkkband !FS bands
602                iband=iFSband+elph_ds%minFSband-1 ! full bands
603                do icomp = 1, 3
604                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
605 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
606                  do jcomp = 1, 3
607                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
608 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
609 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
610 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
611                  end do
612                end do
613              end do
614            end if
615 
616            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
617            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
618            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
619            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
620            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
621            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
622            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
623            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
624 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
625            if (present(elph_tr_ds)) then
626              do iFSband=1,elph_ds%ngkkband !FS bands
627                iband=iFSband+elph_ds%minFSband-1 ! full bands
628                do icomp = 1, 3
629                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
630 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
631                  do jcomp = 1, 3
632                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
633 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
634 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
635 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
636                  end do
637                end do
638              end do
639            end if
640 
641            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
642            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
643            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
644            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
645            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
646            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
647            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
648            elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
649 &           0.25_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
650            if (present(elph_tr_ds)) then
651              do iFSband=1,elph_ds%ngkkband !FS bands
652                iband=iFSband+elph_ds%minFSband-1 ! full bands
653                do icomp = 1, 3
654                  elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
655 &                 0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
656                  do jcomp = 1, 3
657                    elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
658 &                   elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
659 &                   0.25_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
660 &                   elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
661                  end do
662                end do
663              end do
664            end if
665          end do
666 !        on the 8 corners
667        end if
668        if (MOD(nscale1,2) == 0 .and. MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
669          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
670          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
671          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
672          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
673          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
674          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
675          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
676          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
677          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
678          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
679          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
680 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
681          if (present(elph_tr_ds)) then
682            do iFSband=1,elph_ds%ngkkband !FS bands
683              iband=iFSband+elph_ds%minFSband-1 ! full bands
684              do icomp = 1, 3
685                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
686 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
687                do jcomp = 1, 3
688                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
689 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
690 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
691 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
692                end do
693              end do
694            end do
695          end if
696 
697          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
698          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
699          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
700          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
701          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
702          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
703          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
704          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
705          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
706          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
707          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
708 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
709          if (present(elph_tr_ds)) then
710            do iFSband=1,elph_ds%ngkkband !FS bands
711              iband=iFSband+elph_ds%minFSband-1 ! full bands
712              do icomp = 1, 3
713                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
714 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
715                do jcomp = 1, 3
716                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
717 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
718 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
719 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
720                end do
721              end do
722            end do
723          end if
724 
725          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
726          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
727          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
728          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
729          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
730          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
731          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
732          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
733          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
734          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
735          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
736 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
737          if (present(elph_tr_ds)) then
738            do iFSband=1,elph_ds%ngkkband !FS bands
739              iband=iFSband+elph_ds%minFSband-1 ! full bands
740              do icomp = 1, 3
741                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
742 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
743                do jcomp = 1, 3
744                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
745 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
746 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
747 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
748                end do
749              end do
750            end do
751          end if
752 
753          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
754          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
755          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
756          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
757          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
758          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
759          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
760          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
761          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
762          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
763          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
764 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
765          if (present(elph_tr_ds)) then
766            do iFSband=1,elph_ds%ngkkband !FS bands
767              iband=iFSband+elph_ds%minFSband-1 ! full bands
768              do icomp = 1, 3
769                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
770 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
771                do jcomp = 1, 3
772                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
773 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
774 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
775 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
776                end do
777              end do
778            end do
779          end if
780 
781          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
782          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
783          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
784          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
785          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
786          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
787          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
788          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
789          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
790          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
791          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
792 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
793          if (present(elph_tr_ds)) then
794            do iFSband=1,elph_ds%ngkkband !FS bands
795              iband=iFSband+elph_ds%minFSband-1 ! full bands
796              do icomp = 1, 3
797                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
798 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
799                do jcomp = 1, 3
800                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
801 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
802 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
803 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
804                end do
805              end do
806            end do
807          end if
808 
809          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
810          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
811          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
812          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
813          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
814          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
815          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
816          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
817          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
818          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
819          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
820 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
821          if (present(elph_tr_ds)) then
822            do iFSband=1,elph_ds%ngkkband !FS bands
823              iband=iFSband+elph_ds%minFSband-1 ! full bands
824              do icomp = 1, 3
825                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
826 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
827                do jcomp = 1, 3
828                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
829 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
830 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
831 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
832                end do
833              end do
834            end do
835          end if
836 
837          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
838          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
839          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
840          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
841          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
842          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
843          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
844          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
845          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
846          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
847          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
848 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
849          if (present(elph_tr_ds)) then
850            do iFSband=1,elph_ds%ngkkband !FS bands
851              iband=iFSband+elph_ds%minFSband-1 ! full bands
852              do icomp = 1, 3
853                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
854 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
855                do jcomp = 1, 3
856                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
857 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
858 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
859 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
860                end do
861              end do
862            end do
863          end if
864 
865          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
866          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
867          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
868          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
869          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
870          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
871          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
872          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
873          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
874          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
875          elph_ds%k_phon%wtk(:,ikpt_phon,:)=elph_ds%k_phon%wtk(:,ikpt_phon,:)+&
876 &         0.125_dp*elph_ds%k_fine%wtk(:,ikpt_fine,:)
877          if (present(elph_tr_ds)) then
878            do iFSband=1,elph_ds%ngkkband !FS bands
879              iband=iFSband+elph_ds%minFSband-1 ! full bands
880              do icomp = 1, 3
881                elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)=elph_ds%k_phon%velocwtk(iFSband,ikpt_phon,icomp,:)+&
882 &               0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)
883                do jcomp = 1, 3
884                  elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:) = &
885 &                 elph_ds%k_phon%vvelocwtk(iFSband,ikpt_phon,icomp,jcomp,:)+&
886 &                 0.125_dp*elph_ds%k_fine%wtk(iFSband,ikpt_fine,:)* &
887 &                 elph_tr_ds%el_veloc(ikpt_fine,iband,icomp,:)*elph_tr_ds%el_veloc(ikpt_fine,iband,jcomp,:)
888                end do
889              end do
890            end do
891          end if
892        end if
893      end do
894    end do
895  end do
896 
897 !bxu, divide by nscale^3 to be consistent with the normalization of kpt_phon
898  elph_ds%k_phon%wtk = elph_ds%k_phon%wtk/nscale1/nscale2/nscale3
899  if (present(elph_tr_ds)) then
900    elph_ds%k_phon%velocwtk = elph_ds%k_phon%velocwtk/nscale1/nscale2/nscale3
901    elph_ds%k_phon%vvelocwtk = elph_ds%k_phon%vvelocwtk/nscale1/nscale2/nscale3
902  end if
903 
904 end subroutine d2c_weights

ABINIT/d2c_wtq [ Functions ]

[ Top ] [ Functions ]

NAME

  d2c_wtq

FUNCTION

  This routine calculates the integration weights on a coarse k-grid
  using the integration weights from a denser k-gird. The weights of
  the extra k points that being shared are evenly distributed.

INPUTS

  elph_ds%k_fine%nkpt = number of fine q-points
  elph_ds%k_fine%wtq = integration weights of the fine q-grid
  elph_ds%k_phon%nkpt = number of coarse q-points

OUTPUT

  elph_ds%k_phon%wtq = integration weights of the coarse k-grid

PARENTS

      mka2f,mka2f_tr

CHILDREN

SOURCE

 931 subroutine d2c_wtq(elph_ds)
 932 
 933 
 934 !This section has been created automatically by the script Abilint (TD).
 935 !Do not modify the following lines by hand.
 936 #undef ABI_FUNC
 937 #define ABI_FUNC 'd2c_wtq'
 938 !End of the abilint section
 939 
 940  implicit none
 941 
 942 !Arguments ------------------------------------
 943  type(elph_type),intent(inout) :: elph_ds
 944 
 945 !Local variables-------------------------------
 946  integer :: ii, jj, kk
 947  integer :: ikpt, jkpt, kkpt
 948  integer :: iikpt, jjkpt, kkkpt
 949  integer :: ikpt_fine, ikpt_phon
 950  integer :: nkpt_fine1, nkpt_phon1
 951  integer :: nkpt_fine2, nkpt_phon2
 952  integer :: nkpt_fine3, nkpt_phon3
 953  integer :: nscale1, nscale2, nscale3
 954 
 955 ! *************************************************************************
 956  nkpt_phon1 = elph_ds%kptrlatt(1,1)
 957  nkpt_phon2 = elph_ds%kptrlatt(2,2)
 958  nkpt_phon3 = elph_ds%kptrlatt(3,3)
 959  nkpt_fine1 = elph_ds%kptrlatt_fine(1,1)
 960  nkpt_fine2 = elph_ds%kptrlatt_fine(2,2)
 961  nkpt_fine3 = elph_ds%kptrlatt_fine(3,3)
 962  nscale1 = dble(nkpt_fine1/nkpt_phon1)
 963  nscale2 = dble(nkpt_fine2/nkpt_phon2)
 964  nscale3 = dble(nkpt_fine3/nkpt_phon3)
 965  if (abs(INT(nscale1)-nscale1) > 0.01) then
 966    MSG_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
 967  end if
 968  if (abs(INT(nscale2)-nscale2) > 0.01) then
 969    MSG_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
 970  end if
 971  if (abs(INT(nscale3)-nscale3) > 0.01) then
 972    MSG_ERROR('The denser k-gird MUST be multiples of the phon k-grid')
 973  end if
 974  nscale1 = INT(nscale1)
 975  nscale2 = INT(nscale2)
 976  nscale3 = INT(nscale3)
 977 
 978 !bxu, get wtq of coarse grid from fine grid
 979  elph_ds%k_phon%wtq = zero
 980 
 981  do ikpt = 1, nkpt_phon1
 982    do jkpt = 1, nkpt_phon2
 983      do kkpt = 1, nkpt_phon3
 984        ikpt_phon = kkpt + (jkpt-1)*nkpt_phon3 + (ikpt-1)*nkpt_phon2*nkpt_phon3
 985 !      inside the paralellepipe
 986        do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
 987          do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
 988            do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
 989              iikpt = 1 + (ikpt-1)*nscale1 + ii
 990              jjkpt = 1 + (jkpt-1)*nscale2 + jj
 991              kkkpt = 1 + (kkpt-1)*nscale3 + kk
 992              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
 993              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
 994              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
 995              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
 996              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
 997              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
 998              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
 999              elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1000 &             elph_ds%k_fine%wtq(:,ikpt_fine,:)
1001            end do
1002          end do
1003        end do
1004 !      on the 6 faces
1005        if (MOD(nscale3,2) == 0) then ! when nscale3 is an even number
1006          do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
1007            do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
1008              iikpt = 1 + (ikpt-1)*nscale1 + ii
1009              jjkpt = 1 + (jkpt-1)*nscale2 + jj
1010              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1011              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1012              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1013              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1014 
1015              kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1016              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1017              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1018              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1019              elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1020 &             0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1021 
1022              kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1023              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1024              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1025              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1026              elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1027 &             0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1028            end do
1029          end do
1030        end if
1031        if (MOD(nscale2,2) == 0) then ! when nscale2 is an even number
1032          do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
1033            do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
1034              iikpt = 1 + (ikpt-1)*nscale1 + ii
1035              kkkpt = 1 + (kkpt-1)*nscale3 + kk
1036              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1037              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1038              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1039              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1040 
1041              jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1042              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1043              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1044              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1045              elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1046 &             0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1047 
1048              jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1049              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1050              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1051              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1052              elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1053 &             0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1054            end do
1055          end do
1056        end if
1057        if (MOD(nscale1,2) == 0) then ! when nscale1 is an even number
1058          do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
1059            do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
1060              kkkpt = 1 + (kkpt-1)*nscale3 + kk
1061              jjkpt = 1 + (jkpt-1)*nscale2 + jj
1062              if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1063              if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1064              if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1065              if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1066 
1067              iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1068              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1069              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1070              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1071              elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1072 &             0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1073 
1074              iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1075              if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1076              if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1077              ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1078              elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1079 &             0.5_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1080            end do
1081          end do
1082 !        on the 12 sides
1083        end if
1084        if (MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
1085          do ii = -((nscale1+1)/2-1), ((nscale1+1)/2-1)
1086            iikpt = 1 + (ikpt-1)*nscale1 + ii
1087            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1088            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1089 
1090            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1091            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1092            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1093            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1094            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1095            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1096            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1097            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1098 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1099 
1100            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1101            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1102            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1103            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1104            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1105            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1106            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1107            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1108 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1109 
1110            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1111            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1112            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1113            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1114            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1115            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1116            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1117            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1118 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1119 
1120            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1121            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1122            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1123            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1124            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1125            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1126            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1127            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1128 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1129          end do
1130        end if
1131        if (MOD(nscale1,2) == 0 .and. MOD(nscale3,2) == 0) then
1132          do jj = -((nscale2+1)/2-1), ((nscale2+1)/2-1)
1133            jjkpt = 1 + (jkpt-1)*nscale2 + jj
1134            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1135            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1136 
1137            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1138            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1139            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1140            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1141            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1142            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1143            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1144            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1145 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1146 
1147            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1148            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1149            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1150            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1151            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1152            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1153            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1154            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1155 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1156 
1157            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1158            kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1159            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1160            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1161            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1162            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1163            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1164            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1165 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1166 
1167            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1168            kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1169            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1170            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1171            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1172            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1173            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1174            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1175 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1176          end do
1177        end if
1178        if (MOD(nscale2,2) == 0 .and. MOD(nscale1,2) == 0) then
1179          do kk = -((nscale3+1)/2-1), ((nscale3+1)/2-1)
1180            kkkpt = 1 + (kkpt-1)*nscale3 + kk
1181            if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1182            if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1183 
1184            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1185            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1186            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1187            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1188            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1189            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1190            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1191            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1192 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1193 
1194            jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1195            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1196            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1197            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1198            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1199            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1200            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1201            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1202 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1203 
1204            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1205            iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1206            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1207            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1208            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1209            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1210            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1211            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1212 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1213 
1214            jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1215            iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1216            if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1217            if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1218            if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1219            if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1220            ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1221            elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1222 &           0.25_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1223          end do
1224 !        on the 8 corners
1225        end if
1226        if (MOD(nscale1,2) == 0 .and. MOD(nscale2,2) == 0 .and. MOD(nscale3,2) == 0) then
1227          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1228          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1229          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1230          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1231          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1232          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1233          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1234          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1235          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1236          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1237          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1238 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1239 
1240          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1241          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1242          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1243          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1244          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1245          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1246          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1247          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1248          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1249          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1250          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1251 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1252 
1253          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1254          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1255          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1256          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1257          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1258          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1259          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1260          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1261          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1262          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1263          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1264 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1265 
1266          iikpt = 1 + (ikpt-1)*nscale1 + nscale1/2
1267          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1268          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1269          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1270          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1271          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1272          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1273          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1274          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1275          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1276          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1277 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1278 
1279          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1280          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1281          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1282          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1283          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1284          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1285          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1286          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1287          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1288          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1289          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1290 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1291 
1292          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1293          jjkpt = 1 + (jkpt-1)*nscale2 + nscale2/2
1294          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1295          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1296          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1297          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1298          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1299          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1300          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1301          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1302          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1303 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1304 
1305          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1306          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1307          kkkpt = 1 + (kkpt-1)*nscale3 + nscale3/2
1308          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1309          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1310          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1311          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1312          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1313          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1314          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1315          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1316 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1317 
1318          iikpt = 1 + (ikpt-1)*nscale1 - nscale1/2
1319          jjkpt = 1 + (jkpt-1)*nscale2 - nscale2/2
1320          kkkpt = 1 + (kkpt-1)*nscale3 - nscale3/2
1321          if (iikpt .le. 0) iikpt = iikpt + nkpt_fine1
1322          if (jjkpt .le. 0) jjkpt = jjkpt + nkpt_fine2
1323          if (kkkpt .le. 0) kkkpt = kkkpt + nkpt_fine3
1324          if (iikpt .gt. nkpt_fine1) iikpt = iikpt - nkpt_fine1
1325          if (jjkpt .gt. nkpt_fine2) jjkpt = jjkpt - nkpt_fine2
1326          if (kkkpt .gt. nkpt_fine3) kkkpt = kkkpt - nkpt_fine3
1327          ikpt_fine = kkkpt + (jjkpt-1)*nkpt_fine3 + (iikpt-1)*nkpt_fine2*nkpt_fine3
1328          elph_ds%k_phon%wtq(:,ikpt_phon,:)=elph_ds%k_phon%wtq(:,ikpt_phon,:)+&
1329 &         0.125_dp*elph_ds%k_fine%wtq(:,ikpt_fine,:)
1330        end if
1331      end do
1332    end do
1333  end do
1334 
1335 !bxu, divide by nscale^3 to be consistent with the normalization of kpt_phon
1336  elph_ds%k_phon%wtq = elph_ds%k_phon%wtq/nscale1/nscale2/nscale3
1337 
1338 end subroutine d2c_wtq

ABINIT/ep_el_weights [ Functions ]

[ Top ] [ Functions ]

NAME

 ep_el_weights

FUNCTION

 This routine calculates the Fermi Surface integration weights
 for the electron phonon routines, by different methods

    1) Gaussian smearing
    2) Tetrahedron method
    3) Window in bands for all k-points
    4) Fermi Dirac smearing, follows gaussian with a different smearing function

INPUTS

   ep_b_min = minimal band to include in FS window integration
   ep_b_max = maximal band to include in FS window integration
   eigenGS = Ground State eigenvalues
   elphsmear = smearing width for Gaussian method
   fermie = Fermi level
   gprimd = Reciprocal lattice vectors (dimensionful)
   irredtoGS = mapping of elph k-points to ground state grid
   kptrlatt = k-point grid vectors (if divided by determinant of present matrix)
   max_occ = maximal occupancy for a band
   minFSband = minimal band included for Fermi Surface integration in Gaussian and Tetrahedron cases
   nFSband = number of bands in FS integration
   nsppol = number of spin polarizations
   telphint = option for FS integration:
      0 Tetrahedron method
      1 Gaussian smearing
      2 Window in bands for all k-points
      3 Fermi Dirac smearing
   k_obj%nkpt = number of FS k-points
   k_obj%kpt = FS k-points
   k_obj%full2irr = mapping of FS k-points from full grid to irred points
   k_obj%full2full = mapping of FS k-points in full grid under symops

OUTPUT

TODO

   weights should be recalculated on-the-fly! The present implementation is not flexible!

PARENTS

      get_nv_fs_en,get_tau_k

CHILDREN

      destroy_tetra,get_tetra_weight,init_tetra,matr3inv

SOURCE

1391 subroutine ep_el_weights(ep_b_min, ep_b_max, eigenGS, elphsmear, enemin, enemax, nene, gprimd, &
1392 &    irredtoGS, kptrlatt, max_occ, minFSband, nband, nFSband, nsppol, telphint, k_obj, tmp_wtk)
1393 
1394 
1395 !This section has been created automatically by the script Abilint (TD).
1396 !Do not modify the following lines by hand.
1397 #undef ABI_FUNC
1398 #define ABI_FUNC 'ep_el_weights'
1399 !End of the abilint section
1400 
1401  implicit none
1402 
1403 !Arguments ------------------------------------
1404 !scalars
1405  type(elph_kgrid_type), intent(in) :: k_obj
1406  integer, intent(in) :: ep_b_min
1407  integer, intent(in) :: ep_b_max
1408  integer,intent(in) :: nband,nene
1409  real(dp), intent(in) :: elphsmear
1410  real(dp), intent(in) :: enemin,enemax
1411  real(dp), intent(in) :: gprimd(3,3)
1412  integer, intent(in) :: kptrlatt(3,3)
1413  real(dp), intent(in) :: max_occ
1414  integer, intent(in) :: minFSband
1415  integer, intent(in) :: nFSband
1416  integer, intent(in) :: nsppol
1417  integer, intent(in) :: telphint
1418 
1419 ! arrays
1420  real(dp), intent(in) :: eigenGS(nband,k_obj%nkptirr,nsppol)
1421  real(dp), intent(out) :: tmp_wtk(nFSband,k_obj%nkpt,nsppol,nene)
1422  integer, intent(in) :: irredtoGS(k_obj%nkptirr)
1423 
1424 !Local variables-------------------------------
1425 !scalars
1426  integer,parameter :: bcorr0=0
1427  integer :: ikpt, ikptgs, ib1, iband
1428  integer :: ierr, ie, isppol
1429  real(dp) :: deltaene, rcvol, fermie
1430  real(dp) :: smdeltaprefactor, smdeltafactor, xx
1431 
1432 ! arrays
1433  real(dp) :: rlatt(3,3), klatt(3,3)
1434  real(dp), allocatable :: tmp_eigen(:), tweight(:,:), dtweightde(:,:)
1435  character (len=500) :: message
1436  character (len=80) :: errstr
1437  type(t_tetrahedron) :: tetrahedra
1438 
1439 ! *************************************************************************
1440 
1441  ! Initialize tmp_wtk with zeros
1442  tmp_wtk = zero
1443 
1444  !write(std_out,*) 'ep_el : nkpt ', k_obj%nkpt
1445 !===================================
1446 !Set up integration weights for FS
1447 !===================================
1448  deltaene = (enemax-enemin)/dble(nene-1)
1449 
1450  if (telphint == 0) then
1451 
1452 !  =========================
1453 !  Tetrahedron integration
1454 !  =========================
1455 
1456    rlatt(:,:) = kptrlatt(:,:)
1457    call matr3inv(rlatt,klatt)
1458 
1459    call init_tetra(k_obj%full2full(1,1,:), gprimd,klatt,k_obj%kpt, k_obj%nkpt,&
1460 &   tetrahedra, ierr, errstr)
1461    ABI_CHECK(ierr==0,errstr)
1462 
1463    rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
1464 &   -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
1465 &   +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
1466 
1467 !  fix small window around fermie for tetrahedron weight calculation
1468    deltaene = (enemax-enemin)/dble(nene-1)
1469 
1470    ABI_ALLOCATE(tmp_eigen,(k_obj%nkpt))
1471    ABI_ALLOCATE(tweight,(k_obj%nkpt,nene))
1472    ABI_ALLOCATE(dtweightde,(k_obj%nkpt,nene))
1473 
1474    do iband = 1,nFSband
1475 !    for each spin pol
1476      do isppol=1,nsppol
1477 !    For this band get its contribution
1478        tmp_eigen(:) = zero
1479        do ikpt=1,k_obj%nkpt
1480          ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1481          tmp_eigen(ikpt) = eigenGS(minFSband+iband-1,ikptgs,isppol)
1482        end do
1483 !      calculate general integration weights at each irred kpoint 
1484 !      as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
1485        call get_tetra_weight(tmp_eigen,enemin,enemax,&
1486 &       max_occ,nene,k_obj%nkpt,tetrahedra,bcorr0,&
1487 &       tweight,dtweightde,xmpi_comm_self)
1488 
1489        tmp_wtk(iband,:,isppol,:) = dtweightde(:,:)*k_obj%nkpt
1490      end do
1491    end do
1492    ABI_DEALLOCATE(tmp_eigen)
1493    ABI_DEALLOCATE(tweight)
1494    ABI_DEALLOCATE(dtweightde)
1495 
1496    call destroy_tetra(tetrahedra)
1497 
1498  else if (telphint == 1) then
1499 
1500 !  ==============================================================
1501 !  Gaussian or integration:
1502 !  Each kpt contributes a gaussian of integrated weight 1
1503 !  for each band. The gaussian being centered at the Fermi level.
1504 !  ===============================================================
1505 
1506 !  took out factor 1/k_obj%nkpt which intervenes only at integration time
1507 
1508 !  MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1509 
1510 !  gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1511    smdeltaprefactor = max_occ*sqrt(piinv)/elphsmear
1512    smdeltafactor = one/elphsmear
1513 
1514 !  SPPOL loop on isppol as well to get 2 sets of weights
1515    do isppol=1,nsppol
1516      fermie = enemin
1517      do ie = 1, nene
1518        fermie = fermie + deltaene
1519 !      fine grid
1520        do ikpt=1, k_obj%nkpt
1521          ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1522          do ib1=1,nFSband
1523            xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1524            if (abs(xx) < 40._dp) then
1525              tmp_wtk(ib1,ikpt,isppol,ie) = exp(-xx*xx)*smdeltaprefactor
1526            end if
1527          end do
1528        end do
1529      end do
1530    end do
1531 
1532 
1533  else if (telphint == 2) then ! range of bands occupied
1534 
1535 !  SPPOL eventually be able to specify bands for up and down separately
1536    fermie = enemin
1537    do ie = 1, nene
1538      fermie = fermie + deltaene
1539      do ikpt=1,k_obj%nkpt
1540        do ib1=ep_b_min, ep_b_max
1541 !        for the moment both spin channels same
1542          tmp_wtk(ib1,ikpt,:,ie) = max_occ
1543        end do
1544      end do
1545    end do
1546 
1547    write(std_out,*) ' ep_el_weights : DOS is calculated from states in bands ',ep_b_min,' to ',ep_b_max
1548 
1549  else if (telphint == 3) then
1550 
1551 !  ==============================================================
1552 !  Fermi Dirac integration:
1553 !  Each kpt contributes a Fermi Dirac smearing function of integrated weight 1
1554 !  for each band. The function being centered at the Fermi level.
1555 !  ===============================================================
1556 
1557 !  took out factor 1/k_obj%nkpt which intervenes only at integration time
1558 
1559 !  MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1560 
1561 !  gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1562    smdeltaprefactor = half*max_occ/elphsmear
1563    smdeltafactor = one/elphsmear
1564 
1565 !  SPPOL loop on isppol as well to get 2 sets of weights
1566    do isppol=1,nsppol
1567      fermie = enemin
1568      do ie = 1, nene
1569        fermie = fermie + deltaene
1570 !      fine grid
1571        do ikpt=1, k_obj%nkpt
1572          ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1573          do ib1=1,nFSband
1574            xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1575            tmp_wtk(ib1,ikpt,isppol,ie) = smdeltaprefactor / (one + cosh(xx))
1576          end do
1577        end do
1578      end do
1579    end do
1580 
1581 
1582  else
1583    write (message,'(a,i0)')" telphint should be between 0 and 3, found: ",telphint
1584    MSG_BUG(message)
1585  end if ! if telphint
1586 
1587 end subroutine ep_el_weights

ABINIT/ep_fs_weights [ Functions ]

[ Top ] [ Functions ]

NAME

 ep_fs_weights

FUNCTION

 This routine calculates the Fermi Surface integration weights
  for the electron phonon routines, by different methods
    1) Gaussian smearing
    2) Tetrahedron method
    3) Window in bands for all k-points
    4) Fermi Dirac smearing, follows gaussian with a different smearing function

INPUTS

   ep_b_min = minimal band to include in FS window integration
   ep_b_max = maximal band to include in FS window integration
   eigenGS = Ground State eigenvalues
   elphsmear = smearing width for Gaussian method
   fermie = Fermi level
   gprimd = Reciprocal lattice vectors (dimensionful)
   irredtoGS = mapping of elph k-points to ground state grid
   kptrlatt = k-point grid vectors (if divided by determinant of present matrix)
   max_occ = maximal occupancy for a band
   minFSband = minimal band included for Fermi Surface integration in Gaussian and Tetrahedron cases
   nFSband = number of bands in FS integration
   nsppol = number of spin polarizations
   telphint = option for FS integration:
      0 Tetrahedron method
      1 Gaussian smearing
      2 Window in bands for all k-points
      3 Fermi Dirac smearing
   k_obj%nkpt = number of FS k-points
   k_obj%kpt = FS k-points
   k_obj%full2irr = mapping of FS k-points from full grid to irred points
   k_obj%full2full = mapping of FS k-points in full grid under symops

OUTPUT

   k_obj%wtk = integration weights

TODO

   weights should be recalculated on-the-fly! The present implementation is not flexible!

PARENTS

      elphon,get_nv_fs_en,get_nv_fs_temp

CHILDREN

      destroy_tetra,get_tetra_weight,init_tetra,matr3inv,wrtout

SOURCE

1640 subroutine ep_fs_weights(ep_b_min, ep_b_max, eigenGS, elphsmear, fermie, gprimd, &
1641 &    irredtoGS, kptrlatt, max_occ, minFSband, nband, nFSband, nsppol, telphint, k_obj)
1642 
1643 
1644 !This section has been created automatically by the script Abilint (TD).
1645 !Do not modify the following lines by hand.
1646 #undef ABI_FUNC
1647 #define ABI_FUNC 'ep_fs_weights'
1648 !End of the abilint section
1649 
1650  implicit none
1651 
1652 !Arguments ------------------------------------
1653 !scalars
1654  type(elph_kgrid_type), intent(inout) :: k_obj
1655  integer, intent(in) :: ep_b_min
1656  integer, intent(in) :: ep_b_max
1657  integer,intent(in) :: nband
1658  real(dp), intent(in) :: elphsmear
1659  real(dp), intent(in) :: fermie
1660  real(dp), intent(in) :: gprimd(3,3)
1661  integer, intent(in) :: kptrlatt(3,3)
1662  real(dp), intent(in) :: max_occ
1663  integer, intent(in) :: minFSband
1664  integer, intent(in) :: nFSband
1665  integer, intent(in) :: nsppol
1666  integer, intent(in) :: telphint
1667 
1668 ! arrays
1669  real(dp), intent(in) :: eigenGS(nband,k_obj%nkptirr,nsppol)
1670  integer, intent(in) :: irredtoGS(k_obj%nkptirr)
1671 
1672 !Local variables-------------------------------
1673 !scalars
1674  integer,parameter :: bcorr0=0
1675  integer :: ikpt, ikptgs, ib1, isppol, iband
1676  integer :: nene, ifermi
1677  integer :: ierr
1678 
1679  real(dp) :: enemin, enemax, deltaene, rcvol
1680  real(dp) :: smdeltaprefactor, smdeltafactor, xx
1681 
1682 ! arrays
1683  real(dp) :: rlatt(3,3), klatt(3,3)
1684  real(dp), allocatable :: tmp_eigen(:), tweight(:,:), dtweightde(:,:)
1685 
1686  character (len=500) :: message
1687  character (len=80) :: errstr
1688 
1689  type(t_tetrahedron) :: tetrahedra
1690 
1691 ! *************************************************************************
1692 
1693  write(std_out,*) 'ep_fs : nkpt ', k_obj%nkpt
1694  write(message, '(a)' ) '- ep_fs_weights  1  = '
1695  call wrtout(std_out,message,'PERS')
1696 
1697 !===================================
1698 !Set up integration weights for FS
1699 !===================================
1700 
1701  if (telphint == 0) then
1702 
1703 !  =========================
1704 !  Tetrahedron integration
1705 !  =========================
1706 
1707    rlatt(:,:) = kptrlatt(:,:)
1708    call matr3inv(rlatt,klatt)
1709 
1710    call init_tetra(k_obj%full2full(1,1,:), gprimd,klatt,k_obj%kpt, k_obj%nkpt,&
1711 &   tetrahedra, ierr, errstr)
1712    ABI_CHECK(ierr==0,errstr)
1713 
1714    rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
1715 &   -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
1716 &   +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
1717 
1718 !  just do weights at FS
1719    nene = 100
1720 
1721 !  fix small window around fermie for tetrahedron weight calculation
1722    deltaene = 2*elphsmear/dble(nene-1)
1723    ifermi = int(nene/2)
1724    enemin = fermie - dble(ifermi-1)*deltaene
1725    enemax = enemin + dble(nene-1)*deltaene
1726 
1727    ABI_ALLOCATE(tmp_eigen,(k_obj%nkpt))
1728    ABI_ALLOCATE(tweight,(k_obj%nkpt,nene))
1729    ABI_ALLOCATE(dtweightde,(k_obj%nkpt,nene))
1730 
1731    do iband = 1,nFSband
1732 !    for each spin pol
1733      do isppol=1,nsppol
1734 !      For this band get its contribution
1735        tmp_eigen(:) = zero
1736        do ikpt=1,k_obj%nkpt
1737          ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1738          tmp_eigen(ikpt) = eigenGS(minFSband+iband-1,ikptgs,isppol)
1739        end do
1740 !      calculate general integration weights at each irred kpoint 
1741 !      as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
1742        call get_tetra_weight(tmp_eigen,enemin,enemax,&
1743 &       max_occ,nene,k_obj%nkpt,tetrahedra,bcorr0,&
1744 &       tweight,dtweightde,xmpi_comm_self)
1745 
1746        k_obj%wtk(iband,:,isppol) = dtweightde(:,ifermi)*k_obj%nkpt
1747      end do
1748 
1749    end do
1750    ABI_DEALLOCATE(tmp_eigen)
1751    ABI_DEALLOCATE(tweight)
1752    ABI_DEALLOCATE(dtweightde)
1753 
1754    call destroy_tetra(tetrahedra)
1755 
1756  else if (telphint == 1) then
1757 
1758 !  ==============================================================
1759 !  Gaussian or integration:
1760 !  Each kpt contributes a gaussian of integrated weight 1
1761 !  for each band. The gaussian being centered at the Fermi level.
1762 !  ===============================================================
1763 
1764 !  took out factor 1/k_obj%nkpt which intervenes only at integration time
1765 
1766 !  MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1767 
1768 !  gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1769    smdeltaprefactor = max_occ*sqrt(piinv)/elphsmear
1770    smdeltafactor = one/elphsmear
1771 
1772    k_obj%wtk = zero
1773 !  SPPOL loop on isppol as well to get 2 sets of weights
1774    do isppol=1,nsppol
1775 !    fine grid
1776      do ikpt=1, k_obj%nkpt
1777        ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1778        do ib1=1,nFSband
1779          xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1780          if (abs(xx) < 40._dp) then
1781            k_obj%wtk(ib1,ikpt,isppol) = exp(-xx*xx)*smdeltaprefactor
1782          end if
1783        end do
1784      end do
1785    end do
1786 
1787 
1788  else if (telphint == 2) then ! range of bands occupied
1789 
1790 !  SPPOL eventually be able to specify bands for up and down separately
1791    k_obj%wtk = zero
1792    do ikpt=1,k_obj%nkpt
1793      do ib1=ep_b_min, ep_b_max
1794 !      for the moment both spin channels same
1795        k_obj%wtk(ib1,ikpt,:) = max_occ
1796      end do
1797    end do
1798 
1799    write(std_out,*) ' ep_fs_weights : DOS is calculated from states in bands ',ep_b_min,' to ',ep_b_max
1800 
1801  else if (telphint == 3) then
1802 
1803 !  ==============================================================
1804 !  Fermi Dirac integration:
1805 !  Each kpt contributes a Fermi Dirac smearing function of integrated weight 1
1806 !  for each band. The function being centered at the Fermi level.
1807 !  ===============================================================
1808 
1809 !  took out factor 1/k_obj%nkpt which intervenes only at integration time
1810 
1811 !  MJV 18/5/2008 does smdeltaprefactor need to contain max_occ?
1812 
1813 !  gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1814    smdeltaprefactor = half*max_occ/elphsmear
1815    smdeltafactor = one/elphsmear
1816 
1817    k_obj%wtk = zero
1818 !  SPPOL loop on isppol as well to get 2 sets of weights
1819    do isppol=1,nsppol
1820 !    fine grid
1821      do ikpt=1, k_obj%nkpt
1822        ikptgs = irredtoGS(k_obj%full2irr(1,ikpt))
1823        do ib1=1,nFSband
1824          xx = smdeltafactor*(eigenGS(minFSband-1+ib1,ikptgs,isppol) - fermie)
1825          k_obj%wtk(ib1,ikpt,isppol) = smdeltaprefactor / (one + cosh(xx))
1826        end do
1827      end do
1828    end do
1829 
1830 
1831  else
1832    write (message,'(a,i0)')" telphint should be between 0 and 3, found: ",telphint
1833    MSG_BUG(message)
1834  end if ! if telphint
1835 
1836 end subroutine ep_fs_weights

ABINIT/ep_ph_weights [ Functions ]

[ Top ] [ Functions ]

NAME

 ep_ph_weights

FUNCTION

 This routine calculates the phonon integration weights
  for the electron phonon routines, by different methods
    1) Gaussian smearing
    0) Tetrahedron method

INPUTS

   phfrq = phonon energies
   elphsmear = smearing width for Gaussian method
   omega = input phonon energy
   gprimd = Reciprocal lattice vectors (dimensionful)
   kptrlatt = k-point grid vectors (if divided by determinant of present matrix)
   telphint = option for FS integration:
      0 Tetrahedron method
      1 Gaussian smearing
   k_obj%nkpt = number of FS k-points
   k_obj%kpt = FS k-points
   k_obj%full2full = mapping of FS k-points in full grid under symops

OUTPUT

   tmp_wtq = integration weights

TODO

   weights should be recalculated on-the-fly! The present implementation is not flexible!

PARENTS

      get_tau_k,mka2f,mka2f_tr

CHILDREN

      destroy_tetra,get_tetra_weight,init_tetra,matr3inv

SOURCE

1877 subroutine ep_ph_weights(phfrq,elphsmear,omega_min,omega_max,nomega,gprimd,kptrlatt,nbranch,telphint,k_obj,tmp_wtq)
1878 
1879 
1880 !This section has been created automatically by the script Abilint (TD).
1881 !Do not modify the following lines by hand.
1882 #undef ABI_FUNC
1883 #define ABI_FUNC 'ep_ph_weights'
1884 !End of the abilint section
1885 
1886  implicit none
1887 
1888 !Arguments ------------------------------------
1889 !scalars
1890  type(elph_kgrid_type), intent(inout) :: k_obj
1891  integer,intent(in) :: nbranch
1892  real(dp), intent(in) :: elphsmear
1893  real(dp), intent(in) :: omega_min,omega_max
1894  real(dp), intent(in) :: gprimd(3,3)
1895  integer, intent(in) :: kptrlatt(3,3)
1896  integer, intent(in) :: nomega
1897  integer, intent(in) :: telphint
1898 
1899 ! arrays
1900  real(dp), intent(in) :: phfrq(nbranch,k_obj%nkpt)
1901  real(dp), intent(out) :: tmp_wtq(nbranch,k_obj%nkpt,nomega)
1902 
1903 !Local variables-------------------------------
1904 !scalars
1905  integer,parameter :: bcorr0=0
1906  integer :: ikpt, ib1, ibranch
1907  integer :: ierr, iomega
1908  real(dp) :: rcvol, max_occ
1909  real(dp) :: smdeltaprefactor, smdeltafactor, xx, gaussmaxarg
1910  real(dp) :: domega,omega
1911 
1912 ! arrays
1913  real(dp) :: rlatt(3,3), klatt(3,3)
1914  real(dp), allocatable :: tweight(:,:), dtweightde(:,:)
1915  character (len=80) :: errstr
1916  type(t_tetrahedron) :: tetrahedra
1917 
1918 ! *************************************************************************
1919 
1920  !write(std_out,*) 'ep_ph : nqpt ', k_obj%nkpt
1921 !===================================
1922 !Set up integration weights for FS
1923 !===================================
1924  max_occ = one
1925  gaussmaxarg = sqrt(-log(1.d-100))
1926  domega = (omega_max - omega_min)/(nomega - 1)
1927 
1928  if (telphint == 0) then
1929 
1930 !  =========================
1931 !  Tetrahedron integration
1932 !  =========================
1933 
1934    rlatt(:,:) = kptrlatt(:,:)
1935    call matr3inv(rlatt,klatt)
1936 
1937    call init_tetra(k_obj%full2full(1,1,:), gprimd,klatt,k_obj%kpt, k_obj%nkpt,&
1938 &   tetrahedra, ierr, errstr)
1939    ABI_CHECK(ierr==0,errstr)
1940 
1941    rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
1942 &   -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
1943 &   +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
1944 
1945 !  do all the omega points for tetrahedron weight calculation
1946 
1947    ABI_ALLOCATE(tweight,(k_obj%nkpt,nomega))
1948    ABI_ALLOCATE(dtweightde,(k_obj%nkpt,nomega))
1949 
1950    do ibranch = 1,nbranch
1951      call get_tetra_weight(phfrq(ibranch,:),omega_min,omega_max,&
1952 &     max_occ,nomega,k_obj%nkpt,tetrahedra,bcorr0,&
1953 &     tweight,dtweightde,xmpi_comm_self)
1954 
1955      tmp_wtq(ibranch,:,:) = dtweightde(:,:)*k_obj%nkpt
1956    end do
1957    ABI_DEALLOCATE(tweight)
1958    ABI_DEALLOCATE(dtweightde)
1959 
1960    call destroy_tetra(tetrahedra)
1961 
1962  else if (telphint == 1) then
1963 
1964 !  ==============================================================
1965 !  Gaussian or integration:
1966 !  Each kpt contributes a gaussian of integrated weight 1
1967 !  for each branch. The gaussian being centered at the input energy
1968 !  ===============================================================
1969 
1970 !  took out factor 1/k_obj%nkpt which intervenes only at integration time
1971 
1972 !  gaussian smdeltaprefactor = sqrt(piinv)/elphsmear/k_obj%nkpt
1973    smdeltaprefactor = max_occ*sqrt(piinv)/elphsmear
1974    smdeltafactor = one/elphsmear
1975 
1976    tmp_wtq = zero
1977    omega = omega_min
1978    do iomega = 1, nomega
1979      omega = omega + domega
1980      do ikpt=1, k_obj%nkpt
1981        do ib1=1,nbranch
1982          xx = smdeltafactor*(phfrq(ib1,ikpt)-omega)
1983          if (abs(xx) < gaussmaxarg) then
1984            tmp_wtq(ib1,ikpt,iomega) = exp(-xx*xx)*smdeltaprefactor
1985          end if
1986        end do
1987      end do
1988    end do
1989  end if ! if telphint
1990 
1991 end subroutine ep_ph_weights

ABINIT/m_epweights [ Modules ]

[ Top ] [ Modules ]

NAME

  m_epweights

FUNCTION

COPYRIGHT

  Copyright (C) 2012-2018 ABINIT group (BXU, MVer)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_epweights
28 
29  use defs_basis
30  use defs_elphon
31  use m_abicore
32  use m_errors
33  use m_tetrahedron
34  use m_xmpi
35 
36  use m_symtk,           only : matr3inv
37 
38  implicit none
39 
40  private