1!
2! Copyright (C) 2001-2010 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!----------------------------------------------------------------------------
9SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
10   !----------------------------------------------------------------------------
11   !! Sorts k+g in order of increasing magnitude, up to ecut.
12   !
13   !! NB: this version should yield the same ordering for different ecut
14   !!     and the same ordering in all machines AS LONG AS INPUT DATA
15   !!     IS EXACTLY THE SAME
16   !
17   USE kinds,      ONLY: DP
18   USE constants,  ONLY: eps8
19   USE wvfct,      ONLY: npwx
20   !
21   IMPLICIT NONE
22   !
23   REAL(DP), INTENT(IN)  :: k(3)
24   !! the k point
25   INTEGER,  INTENT(IN)  :: ngm
26   !! the number of g vectors
27   REAL(DP), INTENT(IN)  :: g(3,ngm)
28   !! the coordinates of G vectors
29   REAL(DP), INTENT(IN)  :: ecut
30   !! the cut-off energy
31   INTEGER,  INTENT(OUT) :: ngk
32   !! the number of k+G vectors inside the "ecut sphere"
33   INTEGER,  INTENT(OUT) :: igk(npwx)
34   !! the correspondence k+G <-> G
35   REAL(DP), INTENT(OUT) :: gk(npwx)
36   !! the moduli of k+G
37   !
38   !  ... local variables
39   !
40   INTEGER :: ng   ! counter on   G vectors
41   INTEGER :: nk   ! counter on k+G vectors
42   REAL(DP) :: q   ! |k+G|^2
43   REAL(DP) :: q2x ! upper bound for |G|
44   !
45   ! ... first we count the number of k+G vectors inside the cut-off sphere
46   !
47   q2x = ( SQRT( SUM(k(:)**2) ) + SQRT( ecut ) )**2
48   !
49   ngk = 0
50   igk(:) = 0
51   gk (:) = 0.0_DP
52   !
53   DO ng = 1, ngm
54      q = SUM( ( k(:) + g(:,ng) )**2 )
55      IF ( q <= eps8 ) q = 0.0_DP
56      !
57      ! ... here if |k+G|^2 <= Ecut
58      !
59      IF ( q <= ecut ) THEN
60         ngk = ngk + 1
61         IF ( ngk > npwx ) &
62            CALL errore( 'gk_sort', 'array gk out-of-bounds', 1 )
63         !
64         gk(ngk) = q
65         !
66         ! set the initial value of index array
67         igk(ngk) = ng
68      ELSE
69         ! if |G| > |k| + SQRT( Ecut )  stop search and order vectors
70         IF ( SUM( g(:,ng)**2 ) > ( q2x + eps8 ) ) EXIT
71      ENDIF
72   ENDDO
73   !
74   IF ( ng > ngm ) &
75      CALL infomsg( 'gk_sort', 'unexpected exit from do-loop' )
76   !
77   ! ... order vector gk keeping initial position in index
78   !
79   CALL hpsort_eps( ngk, gk, igk, eps8 )
80   !
81   ! ... now order true |k+G|
82   !
83   DO nk = 1, ngk
84      gk(nk) = SUM( (k(:) + g(:,igk(nk)) )**2 )
85   ENDDO
86   !
87END SUBROUTINE gk_sort
88