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