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