1! ---
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt .
6! See Docs/Contributors.txt for a list of contributors.
7! ---
8subroutine compute_pw_matrix( nncount, bvectorsfrac )
9!
10! In this subroutine we compute the matrix elements of the plane wave,
11! for all the wave vectors that connect a given k-point to its nearest
12! neighbours
13!
14! ------------------------------ INPUT -----------------------------------------
15! integer nncount            : The number of nearest neighbours belonging to
16!                                each k-point of the Monkhorst-Pack mesh
17! real(dp) bvectorsfrac(:,:) : Vectors that connect each mesh k-point
18!                              to its nearest neighbours.
19! ------------------------------ OUTPUT ----------------------------------------
20! Matrix elements of the plane wave for all the required wave vectors
21! connecting a given k-point to its nearest neighbours
22! This array (delkmatgen) is defined in the module m_planewavematrixvar,
23! file m_planewavematrixvar.F90
24! ------------------------------ BEHAVIOUR -------------------------------------
25! First of all, we allocate the arrays where the matrix elements of the
26! plane waves will be stored
27!
28! Loop on all the neighbour k-points
29!
30!   For each neighbour, compute the vector that connects each mesh point
31!   to its nearest neighbour in bohr^-1.
32!   Remember that bvectorsfrac are read in reduced units, so we have
33!   to multiply then by the reciprocal lattice vector.
34!   This is done in the function getkvector.
35!
36!   Call to the subroutine that computes the matrix elements of the plane
37!   wave in the grid (planewavematrix)
38!
39!   Store the matrix elements of the plane wave for a given wave vector
40!   (array delkmat)
41!   in a permanent array,
42!   (array delkmatgen)
43!   So delkmat can be rewritten for the next wave vector.
44!
45! ------------------------------ UNITS -----------------------------------------
46!   bvectorsfrac are given in reduced units.
47! ------------------------------------------------------------------------------
48
49! Used module procedures
50  use alloc,             only: re_alloc        ! Re-allocation routine
51  use m_planewavematrix, only: planewavematrix ! To compute the plane wave
52                                               !    matrix elements in the grid
53
54! Used module variables
55  use precision,            only: dp         ! Real douple precision type
56  use m_planewavematrixvar, only: delkmat    ! Matrix elements of a plane wave
57  use m_planewavematrixvar, only: delkmatgen ! Matrix elements of a plane wave
58                                             !   (saved)
59  use sparse_matrices,      only: maxnh      ! Maximum number of orbitals
60                                             !   overlapping
61  use parallel,             only: IOnode     ! Input/Output node
62
63!! For debugging
64!  use m_writedelk,          only: write_delk
65!  use atomlist,             only: no_u, no_s, iaorb, iphorb
66!  use atomlist,             only: indxuo, qtot
67!  use sparse_matrices,      only: maxnh, listh, listhptr, numh, xijo, S
68!  use m_spin,               only: nspin
69!  use siesta_options,       only: temp
70!  use sys,                  only: die
71!! End debugging
72
73  implicit none
74
75! The number of nearest neighbours belonging to
76!   each k-point of the Monkhorst-Pack mesh
77  integer,  intent(in)          :: nncount
78! The vectors b that connect each mesh-point k to its nearest neighbours
79  real(dp), intent(in)          :: bvectorsfrac(3,nncount)
80
81!
82! Internal variables
83!
84  integer  :: inn
85  real(dp) :: bvector(3)
86
87!! For debugging
88!  logical  :: onlygamma
89!! End debugging
90
91  if( IOnode )              &
92 &        write(6,'(/,a)')  &
93 &       'compute_pw_matrix: Computing the matrix elements of a plane wave'
94
95! First of all, we allocate the arrays where the matrix elements of the
96! plane waves will be stored
97  if ( .not. associated(delkmatgen) ) then
98    nullify(delkmatgen)
99    call re_alloc(delkmatgen, 1, nncount, 1, maxnh,                  &
100 &                name='delkmatgen', routine='compute_pw_matrix')
101  endif
102  if ( .not. associated(delkmat) ) then
103    nullify(delkmat)
104    call re_alloc(delkmat, 1, maxnh,  &
105 &                name='delkmat',routine='compute_pw_matrix')
106  endif
107
108! Loop on neighbour k-points
109  do inn = 1, nncount
110
111!   For each neighbour, compute the vector that connects each mesh point
112!   to its nearest neighbour in bohr^-1 (done in the subroutine getkvector).
113!   Remember that bvectorsfrac are read in reduced units, so we have
114!   to multiply then by the reciprocal lattice vector.
115    call getkvector( bvectorsfrac(:,inn), bvector )
116
117!!$    if( IOnode ) then
118!!$      write(6,'(a)')         &
119!!$ &      'compute_pw_matrix: Vector connecting a k-point with its neighbours'
120!!$      write(6,'(a,3f12.5)')  &
121!!$ &      'compute_pw_matrix:', bvector(:)
122!!$    endif
123
124!   Call to the subroutine that computes the matrix elements of the plane
125!   wave in the grid (planewavematrix)
126    call planewavematrix( -1, bvector )
127
128!   Store the matrix elements of the plane wave for a given wave vector
129!   (array delkmat)
130!   in a permanent array,
131!   (array delkmatgen)
132!   So delkmat can be rewritten for the next wave vector.
133    delkmatgen(inn,:) = delkmat(:)
134
135!! For debugging
136!      onlygamma = .false.
137!      call write_delk( onlygamma, no_u, no_s, nspin, indxuo, maxnh, &
138! &                     numh, listhptr, listh, delkmat, S, qtot, temp, xijo)
139!      call die()
140!! End debugging
141
142  enddo
143
144end subroutine compute_pw_matrix
145
146subroutine getkvector(frac, vector)
147!
148! Converts from fractional to Bohr^-1 coordinates in reciprocal space.
149!
150
151use m_siesta2wannier90, only : reclatvec ! Reciprocal lattice vectors
152use precision,          only : dp        ! Real double precision type
153
154integer              :: i,j
155real(dp),intent(in)  :: frac(3)
156real(dp),intent(out) :: vector(3)
157
158
159  do i = 1, 3
160    vector(i) = 0.0_dp
161    do j = 1, 3
162      vector(i) = frac(j)*reclatvec(i,j) + vector(i)
163    enddo
164  enddo
165
166end subroutine getkvector
167
168