1 subroutine vband_ibz ( nk1,nk2,nk3, nbnd, nksfit, etk, eqk, at, vk, dfk ) 2 ! 3 ! Written by Burak Himmetoglu (burakhmmtgl@gmail.com) 4 ! Uses some parts of the PW distribution 5 ! 6 ! Computes band velocities and forward derivatives using finite 7 ! differences. 8 ! 9 10!$ use omp_lib 11 12 implicit none 13 14 integer, intent(in) :: nk1,nk2,nk3, nbnd, nksfit, eqk(nk1*nk2*nk3) 15 16 double precision, intent(in) :: etk(nbnd,nksfit), at(3,3) 17 18 double precision, intent(out) :: vk(nbnd,nk1*nk2*nk3,3), dfk(nbnd,nk1*nk2*nk3,3) 19 ! vk is band velocity 20 ! dfk is forward derivative (used in adaptive smearing) 21 22 integer :: i,j,k,n,ik,ibnd,nktot,nm,np 23 24 double precision :: temp, temp2, vaux(3,nk1*nk2*nk3) 25 26 ! Create the k-points, then compute velocities 27 28 nktot = nk1*nk2*nk3 29 30 !$omp parallel do default(none) & 31 !$omp collapse(4) & 32 !$omp private(ibnd,i,j,k,n,np,nm) & 33 !$omp shared(nbnd,nk1,nk2,nk3,vk,dfk,etk,eqk) 34 do ibnd=1,nbnd 35 do i=1,nk1 36 do j=1,nk2 37 do k=1,nk3 38 n = (k-1) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 39 ! 40 ! velocity along k1 41 ! 42 if ( i .eq. 1 ) then 43 np = (k-1) + (j-1)*nk3 + (i)*nk2*nk3 + 1 44 nm = (k-1) + (j-1)*nk3 + (i-2+nk1)*nk2*nk3 + 1 45 else if (i .eq. nk1) then 46 np = (k-1) + (j-1)*nk3 + (i-nk1)*nk2*nk3 + 1 47 nm = (k-1) + (j-1)*nk3 + (i-2)*nk2*nk3 + 1 48 else 49 np = (k-1) + (j-1)*nk3 + (i)*nk2*nk3 + 1 50 nm = (k-1) + (j-1)*nk3 + (i-2)*nk2*nk3 + 1 51 end if 52 vk(ibnd,n,1) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(nm)) )/(2.0/nk1) 53 dfk(ibnd,n,1) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(n)) )/(2.0/nk1) 54 ! 55 ! velocity along k2 56 ! 57 if ( j .eq. 1 ) then 58 np = (k-1) + (j)*nk3 + (i-1)*nk2*nk3 + 1 59 nm = (k-1) + (j-2+nk2)*nk3 + (i-1)*nk2*nk3 + 1 60 else if (j .eq. nk2) then 61 np = (k-1) + (j-nk2)*nk3 + (i-1)*nk2*nk3 + 1 62 nm = (k-1) + (j-2)*nk3 + (i-1)*nk2*nk3 + 1 63 else 64 np = (k-1) + (j)*nk3 + (i-1)*nk2*nk3 + 1 65 nm = (k-1) + (j-2)*nk3 + (i-1)*nk2*nk3 + 1 66 end if 67 vk(ibnd,n,2) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(nm)) )/(2.0/nk2) 68 dfk(ibnd,n,2) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(n)) )/(2.0/nk2) 69 ! 70 ! velocity along k3 71 ! 72 if ( k .eq. 1 ) then 73 np = (k) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 74 nm = (k-2+nk3) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 75 else if (k .eq. nk3) then 76 np = (k-nk3) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 77 nm = (k-2) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 78 else 79 np = (k) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 80 nm = (k-2) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 81 end if 82 vk(ibnd,n,3) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(nm)) )/(2.0/nk3) 83 dfk(ibnd,n,3) = ( etk(ibnd,eqk(np))-etk(ibnd,eqk(n)) )/(2.0/nk3) 84 ! 85 end do ! k 86 end do ! j 87 end do ! i 88 end do ! ibnd 89 !$omp end parallel do 90 ! 91 ! Convert the derivatives in crystal coordinates to cartesian coordinates 92 ! 93 do ibnd=1,nbnd 94 ! 95 do j=1,3 96 vaux(j,:) = vk(ibnd,:,j) ! temporary vector 97 end do 98 ! 99 ! vaux is updated after this call 100 ! 101 call cryst_to_cart (nktot,vaux,at,1) 102 ! 103 do j=1,3 104 vk(ibnd,:,j) = vaux(j,:) 105 end do 106 ! 107 end do 108 ! 109 end subroutine vband_ibz 110