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