Obscure k+G ordering problem, presents since v.6.0! see comments in routines

This commit is contained in:
Paolo Giannozzi 2022-02-18 15:55:35 +01:00
parent 2d01c25adc
commit f2dfc26039
3 changed files with 37 additions and 58 deletions

View File

@ -16,6 +16,10 @@ Fixed in development version:
* Bugfix in DFPT+U with PAW pseudos and when fildvscf/=''
* Makov-Payne correction wasn't working with ibrav=0 and cubic cell (A. Fonari)
* Card ADDITIONAL_KPOINTS wasn't working as expected (Prasenjit Ghosh)
* Subtle bug in G-vector ordering, usually triggered by almost-but-not-quite
symmetric primitive lattice vectors, was affecting k=0 calculations (not CP)
since v.6.0. The ultimate fix would require to change routine hpsort_eps;
the current workaround consists in not recomputing k+G indices if k=0.
New in 7.0 version:
* Unscreened and screened (range-separated) vdW-DF hybrids

View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2022 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -10,7 +10,15 @@ subroutine hpsort_eps (n, ra, ind, eps)
!---------------------------------------------------------------------
!! Sort an array ra(1:n) into ascending order using heapsort algorithm,
!! and considering two elements being equal if their values differ
!! for less than "eps".
!! for less than "eps". IMPORTANT NOTICE (PG February 2022):
!! Assume you have in input a,b,c, in this order, with |a-b| < eps,
!! |b-c| < eps, but |a-c| > eps. The resulting output order may be a,b,c,
!! not c,a,b as it should. I think this is a bug. I don't know how to fix
!! it, I am not sure it does any harm, but re-ordering k+G with this same
!! routine may yield a different ordering for k+G and G vectors even if k=0.
!! This was definitely a bug that has been around for years and that was
!! worked around by avoiding to recompute the indices for k+G if k=0.
!!
!! \(\text{n}\) is input, \(\text{ra}\) is replaced on output by its
!! sorted rearrangement.
!! Create an index table (ind) by making an exchange in the index array
@ -124,33 +132,12 @@ subroutine hpsort_eps (n, ra, ind, eps)
end do sorting
!
end subroutine hpsort_eps
!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!---------------------------------------------------------------------
subroutine hpsort (n, ra, ind)
!---------------------------------------------------------------------
!! Sort an array ra(1:n) into ascending order using heapsort algorithm.
!! \(\text{n}\) is input, \(\text{ra}\) is replaced on output by its
!! sorted rearrangement.
!! Create an index table (ind) by making an exchange in the index array
!! whenever an exchange is made on the sorted data array (ra).
!! In case of equal values in the data array (ra) the values in the
!! index array (ind) are used to order the entries.
!! If on input ind(1) = 0 then indices are initialized in the routine,
!! If on input ind(1) != 0 then indices are assumed to have been
!! initialized before entering the routine and these indices are carried
!! around during the sorting process.
!
! no work space needed !
! free us from machine-dependent sorting-routines !
!
! adapted from Numerical Recipes pg. 329 (new edition)
!! As hpsort_eps, without the "eps" trick (or equivalently, eps=0)
!
use kinds, only : DP
implicit none
@ -241,34 +228,11 @@ subroutine hpsort (n, ra, ind)
goto 10
!
end subroutine hpsort
!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!---------------------------------------------------------------------
subroutine ihpsort (n, ia, ind)
!---------------------------------------------------------------------
!! Sort an integer array ia(1:n) into ascending order using heapsort
!! algorithm. \(\text{n}\) is input, \(\text{ia}\) is replaced on output
!! by its sorted rearrangement.
!! Create an index table (ind) by making an exchange in the index array
!! whenever an exchange is made on the sorted data array (ia).
!! in case of equal values in the data array (ia) the values in the
!! index array (ind) are used to order the entries.
!! If on input ind(1) = 0 then indices are initialized in the routine,
!! If on input ind(1) != 0 then indices are assumed to have been
!! initialized before entering the routine and these indices are carried
!! around during the sorting process.
!
! no work space needed !
! free us from machine-dependent sorting-routines !
!
! adapted from Numerical Recipes pg. 329 (new edition)
!! As "hpsort", for integer array ia(1:n)
!
implicit none
!-input/output variables

View File

@ -10,9 +10,17 @@ SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
!----------------------------------------------------------------------------
!! Sorts k+g in order of increasing magnitude, up to ecut.
!
!! NB: this version should yield the same ordering for different ecut
!! and the same ordering in all machines AS LONG AS INPUT DATA
!! IS EXACTLY THE SAME
!! NB1: this version should yield the same ordering for different ecut
!! and the same ordering in all machines AS LONG AS INPUT DATA
!! IS EXACTLY THE SAME
!! NB2: this version assumes that input G-vectors are ordered by the same
!! routine, "hpsort_eps", used here. In principle for k=0 this should
!! guarantee that the ordering of k+G is the same as for G-vectors.
!! In practice, in some special cases (primitive lattice vectors that
!! are close to but not exactly equal to a symmetric Bravais lattice)
!! this does not hold, presumably due to a limitation of "hpsort_eps".
!! This is ia source of troouble for Gamma-only calculations, so one
!! explicitly sets igk(i)=i for k=0
!
USE kinds, ONLY: DP
USE constants, ONLY: eps8
@ -75,13 +83,16 @@ SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
CALL infomsg( 'gk_sort', 'unexpected exit from do-loop' )
!
! ... order vector gk keeping initial position in index
! ... see comments above about the k=0 case
!
CALL hpsort_eps( ngk, gk, igk, eps8 )
!
! ... now order true |k+G|
!
DO nk = 1, ngk
gk(nk) = SUM( (k(:) + g(:,igk(nk)) )**2 )
ENDDO
IF ( k(1)**2 + k(2)**2 + k(3)**2 > eps8 ) THEN
CALL hpsort_eps( ngk, gk, igk, eps8 )
!
! ... now order true |k+G|
!
DO nk = 1, ngk
gk(nk) = SUM( (k(:) + g(:,igk(nk)) )**2 )
ENDDO
END IF
!
END SUBROUTINE gk_sort