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.

COPYRIGHT

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

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

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