mirror of https://gitlab.com/QEF/q-e.git
New ggen by Norbert verified and re-commited; useless re-ordering
removed from CP as well; redundant sorting routine kb07ad_cp90 removed; routines sort_gvec moved into cp_fpmd.f90, together with the other G-vector related routines taht shlould one day merged with those of PW git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6474 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
parent
affdef4e48
commit
102b7d2c8c
|
@ -113,7 +113,8 @@ subroutine ggenb (b1b, b2b, b3b, nr1b ,nr2b, nr3b, nr1bx ,nr2bx, nr3bx, gcutb )
|
|||
170 format(' ggenb: # of gb vectors < gcutb ngb = ',i6)
|
||||
END IF
|
||||
|
||||
call kb07ad_cp90 (gb,ngb,idx)
|
||||
idx(1)=0
|
||||
call hpsort (ngb,gb,idx)
|
||||
|
||||
do ig=1,ngb-1
|
||||
icurr=ig
|
||||
|
@ -734,36 +735,6 @@ SUBROUTINE glocal( ng, g, ig_l2g, mill_l, ng_g, g2_g, mill_g, nr1, nr2, nr3, isi
|
|||
|
||||
if( ng /= ng_l ) call errore( ' glocal ', ' inconsistent number of G vectors ', ng_l )
|
||||
|
||||
allocate(idx(ng))
|
||||
!
|
||||
! reorder the local g's in order of increasing magnitude.
|
||||
!
|
||||
call kb07ad_cp90(g,ng,idx)
|
||||
!
|
||||
do ig=1,ng-1
|
||||
icurr=ig
|
||||
30 if(idx(icurr).ne.ig) then
|
||||
|
||||
it=ig_l2g(icurr)
|
||||
ig_l2g(icurr)=ig_l2g(idx(icurr))
|
||||
ig_l2g(idx(icurr))=it
|
||||
|
||||
mill=mill_l(:,icurr)
|
||||
mill_l(:,icurr)=mill_l(:,idx(icurr))
|
||||
mill_l(:,idx(icurr))=mill
|
||||
!
|
||||
it=icurr
|
||||
icurr=idx(icurr)
|
||||
idx(it)=it
|
||||
if(idx(icurr).eq.ig) then
|
||||
idx(icurr)=icurr
|
||||
goto 35
|
||||
endif
|
||||
goto 30
|
||||
endif
|
||||
35 continue
|
||||
end do
|
||||
|
||||
! ... Uncomment to make tests and comparisons with other codes
|
||||
! IF ( ionode ) THEN
|
||||
! DO ig=1,ng
|
||||
|
@ -774,7 +745,7 @@ SUBROUTINE glocal( ng, g, ig_l2g, mill_l, ng_g, g2_g, mill_g, nr1, nr2, nr3, isi
|
|||
! END IF
|
||||
|
||||
|
||||
deallocate( idx )
|
||||
! deallocate( idx )
|
||||
|
||||
RETURN
|
||||
END SUBROUTINE glocal
|
||||
|
@ -813,8 +784,8 @@ SUBROUTINE gchkrefold( ng, mill_l, nr1, nr2, nr3 )
|
|||
nr3m1=(nr3-1)/2
|
||||
end if
|
||||
do ig=1,ng
|
||||
if ( mill_l(1,ig).lt.-nr1m1.or.mill_l(1,ig).gt.nr1m1 .or. &
|
||||
& mill_l(2,ig).lt.-nr2m1.or.mill_l(2,ig).gt.nr2m1 .or. &
|
||||
if ( mill_l(1,ig).lt.-nr1m1.or.mill_l(1,ig).gt.nr1m1 .or. &
|
||||
& mill_l(2,ig).lt.-nr2m1.or.mill_l(2,ig).gt.nr2m1 .or. &
|
||||
& mill_l(3,ig).lt.-nr3m1.or.mill_l(3,ig).gt.nr3m1 ) &
|
||||
& nrefold=nrefold+1
|
||||
end do
|
||||
|
@ -1007,6 +978,65 @@ END SUBROUTINE gshcount
|
|||
!
|
||||
return
|
||||
end subroutine gcal
|
||||
!
|
||||
!---------------------------------------------------------------------
|
||||
|
||||
subroutine sort_gvec( ng, g2, mill )
|
||||
!---------------------------------------------------------------------
|
||||
!
|
||||
! first the input variables
|
||||
!
|
||||
use kinds, ONLY: DP
|
||||
use constants, ONLY: eps8
|
||||
implicit none
|
||||
INTEGER, INTENT(IN) :: ng
|
||||
REAL(DP) :: g2( * )
|
||||
INTEGER :: mill( 3, * )
|
||||
|
||||
REAL(DP), ALLOCATABLE :: gsort( : )
|
||||
INTEGER, ALLOCATABLE :: idx( : )
|
||||
INTEGER :: ig, icurr, it, im
|
||||
REAL(DP) :: gsq
|
||||
|
||||
ALLOCATE( gsort( ng ) )
|
||||
ALLOCATE( idx( ng ) )
|
||||
|
||||
DO ig = 1, ng
|
||||
IF ( g2(ig) > eps8 ) THEN
|
||||
gsort(ig) = g2(ig)
|
||||
ELSE
|
||||
gsort(ig) = 0.d0
|
||||
END IF
|
||||
END DO
|
||||
|
||||
idx(1) = 0
|
||||
CALL hpsort_eps( ng, gsort( 1 ), idx( 1 ), eps8 )
|
||||
|
||||
! ... sort indices accordingly
|
||||
DO ig = 1, ng-1
|
||||
icurr = ig
|
||||
30 IF( idx(icurr) /= ig ) THEN
|
||||
! ... swap g-vec indices
|
||||
im = mill(1,icurr); mill(1,icurr) = mill(1,idx(icurr)); mill(1,idx(icurr)) = im
|
||||
im = mill(2,icurr); mill(2,icurr) = mill(2,idx(icurr)); mill(2,idx(icurr)) = im
|
||||
im = mill(3,icurr); mill(3,icurr) = mill(3,idx(icurr)); mill(3,idx(icurr)) = im
|
||||
! ... swap modules
|
||||
gsq = g2( icurr ); g2( icurr ) = g2( idx(icurr) ); g2( idx(icurr) ) = gsq
|
||||
! ... swap indices
|
||||
it = icurr; icurr = idx(icurr); idx(it) = it
|
||||
IF( idx(icurr) == ig ) THEN
|
||||
idx(icurr) = icurr
|
||||
ELSE
|
||||
GOTO 30
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
|
||||
DEALLOCATE( gsort )
|
||||
DEALLOCATE( idx )
|
||||
|
||||
return
|
||||
end subroutine sort_gvec
|
||||
|
||||
|
||||
!=----------------------------------------------------------------------------=!
|
||||
|
|
242
PW/ggen.f90
242
PW/ggen.f90
|
@ -33,13 +33,12 @@ SUBROUTINE ggen()
|
|||
! here a few local variables
|
||||
!
|
||||
REAL(DP) :: t (3), tt, swap
|
||||
REAL(DP), ALLOCATABLE :: esort (:)
|
||||
!
|
||||
INTEGER :: ngmx, n1, n2, n3, n1s, n2s, n3s
|
||||
!
|
||||
REAL(DP), ALLOCATABLE :: g2sort_g(:)
|
||||
! array containing all g vectors, on all processors: replicated data
|
||||
INTEGER, ALLOCATABLE :: mill_g(:,:)
|
||||
INTEGER, ALLOCATABLE :: mill_g(:,:), mill_unsorted(:,:)
|
||||
! array containing all g vectors generators, on all processors:
|
||||
! replicated data
|
||||
INTEGER, ALLOCATABLE :: igsrt(:)
|
||||
|
@ -63,60 +62,48 @@ SUBROUTINE ggen()
|
|||
! and computes all the g vectors inside a sphere
|
||||
!
|
||||
ALLOCATE( ig_l2g( ngm_l ) )
|
||||
ALLOCATE( mill_g( 3, ngm_g ) )
|
||||
ALLOCATE( mill_g( 3, ngm_g ),mill_unsorted( 3, ngm_g ) )
|
||||
ALLOCATE( igsrt( ngm_g ) )
|
||||
ALLOCATE( g2sort_g( ngm_g ) )
|
||||
g2sort_g(:) = 1.0d20
|
||||
!
|
||||
n1 = nr1 + 1
|
||||
n2 = nr2 + 1
|
||||
n3 = nr3 + 1
|
||||
!
|
||||
! save present value of ngm in ngmx variable
|
||||
!
|
||||
ngmx = ngm
|
||||
!
|
||||
ngm = 0
|
||||
ngms = 0
|
||||
DO i = - n1, n1
|
||||
iloop: DO i = -nr1-1, nr1+1
|
||||
!
|
||||
! Gamma-only: exclude space with x < 0
|
||||
! gamma-only: exclude space with x < 0
|
||||
!
|
||||
IF ( gamma_only .and. i < 0) GOTO 10
|
||||
DO j = - n2, n2
|
||||
IF ( gamma_only .and. i < 0) CYCLE iloop
|
||||
jloop: DO j = -nr2-1, nr2+1
|
||||
!
|
||||
! exclude plane with x = 0, y < 0
|
||||
! gamma-only: exclude plane with x = 0, y < 0
|
||||
!
|
||||
IF ( gamma_only .and. i == 0 .and. j < 0) GOTO 11
|
||||
DO k = - n3, n3
|
||||
IF ( gamma_only .and. i == 0 .and. j < 0) CYCLE jloop
|
||||
kloop: DO k = -nr3-1, nr3+1
|
||||
!
|
||||
! exclude line with x = 0, y = 0, z < 0
|
||||
! gamma-only: exclude line with x = 0, y = 0, z < 0
|
||||
!
|
||||
IF ( gamma_only .and. i == 0 .and. j == 0 .and. k < 0) GOTO 12
|
||||
tt = 0.d0
|
||||
DO ipol = 1, 3
|
||||
t (ipol) = i * bg (ipol, 1) + j * bg (ipol, 2) + k * bg (ipol, 3)
|
||||
tt = tt + t (ipol) * t (ipol)
|
||||
ENDDO
|
||||
IF ( gamma_only .and. i == 0 .and. j == 0 .and. k < 0) CYCLE kloop
|
||||
t(:) = i * bg (:,1) + j * bg (:,2) + k * bg (:,3)
|
||||
tt = sum(t(:)**2)
|
||||
IF (tt <= gcutm) THEN
|
||||
ngm = ngm + 1
|
||||
IF (tt <= gcutms) ngms = ngms + 1
|
||||
IF (ngm > ngm_g) CALL errore ('ggen', 'too many g-vectors', ngm)
|
||||
mill_g( 1, ngm ) = i
|
||||
mill_g( 2, ngm ) = j
|
||||
mill_g( 3, ngm ) = k
|
||||
mill_unsorted( :, ngm ) = (/ i,j,k /)
|
||||
IF ( tt > eps8 ) THEN
|
||||
g2sort_g(ngm) = tt
|
||||
ELSE
|
||||
g2sort_g(ngm) = 0.d0
|
||||
ENDIF
|
||||
ENDIF
|
||||
12 CONTINUE
|
||||
ENDDO
|
||||
11 CONTINUE
|
||||
ENDDO
|
||||
10 CONTINUE
|
||||
ENDDO
|
||||
ENDDO kloop
|
||||
ENDDO jloop
|
||||
ENDDO iloop
|
||||
|
||||
IF (ngm /= ngm_g ) &
|
||||
CALL errore ('ggen', 'g-vectors missing !', abs(ngm - ngm_g))
|
||||
|
@ -125,118 +112,42 @@ SUBROUTINE ggen()
|
|||
|
||||
igsrt(1) = 0
|
||||
CALL hpsort_eps( ngm_g, g2sort_g, igsrt, eps8 )
|
||||
DEALLOCATE( g2sort_g )
|
||||
DO ng = 1, ngm_g-1
|
||||
indsw = ng
|
||||
7 IF(igsrt(indsw) /= ng) THEN
|
||||
! .. swap indices
|
||||
DO i = 1, 3
|
||||
iswap = mill_g(i,indsw)
|
||||
mill_g(i,indsw) = mill_g(i,igsrt(indsw))
|
||||
mill_g(i,igsrt(indsw)) = iswap
|
||||
ENDDO
|
||||
! .. swap indices
|
||||
iswap = indsw; indsw = igsrt(indsw); igsrt(iswap) = iswap
|
||||
IF(igsrt(indsw) == ng) THEN
|
||||
igsrt(indsw)=indsw
|
||||
ELSE
|
||||
GOTO 7
|
||||
ENDIF
|
||||
ENDIF
|
||||
ENDDO
|
||||
mill_g(1,:) = mill_unsorted(1,igsrt(:))
|
||||
mill_g(2,:) = mill_unsorted(2,igsrt(:))
|
||||
mill_g(3,:) = mill_unsorted(3,igsrt(:))
|
||||
DEALLOCATE( g2sort_g, igsrt, mill_unsorted )
|
||||
|
||||
DEALLOCATE( igsrt )
|
||||
|
||||
! WRITE( stdout, fmt="(//,' --- Executing new GGEN Loop ---',//)" )
|
||||
|
||||
ALLOCATE(esort(ngm) )
|
||||
esort(:) = 1.0d20
|
||||
ngm = 0
|
||||
ngms = 0
|
||||
DO ng = 1, ngm_g
|
||||
ngloop: DO ng = 1, ngm_g
|
||||
i = mill_g(1, ng)
|
||||
j = mill_g(2, ng)
|
||||
k = mill_g(3, ng)
|
||||
|
||||
#ifdef __PARA
|
||||
m1 = mod (i, nr1) + 1
|
||||
IF (m1.lt.1) m1 = m1 + nr1
|
||||
IF (m1 < 1) m1 = m1 + nr1
|
||||
m2 = mod (j, nr2) + 1
|
||||
IF (m2.lt.1) m2 = m2 + nr2
|
||||
IF (m2 < 1) m2 = m2 + nr2
|
||||
mc = m1 + (m2 - 1) * nrx1
|
||||
IF ( dfftp%isind ( mc ) .eq.0) GOTO 1
|
||||
IF ( dfftp%isind ( mc ) == 0) CYCLE ngloop
|
||||
#endif
|
||||
|
||||
tt = 0.d0
|
||||
DO ipol = 1, 3
|
||||
t (ipol) = i * bg (ipol, 1) + j * bg (ipol, 2) + k * bg (ipol, 3)
|
||||
tt = tt + t (ipol) * t (ipol)
|
||||
ENDDO
|
||||
|
||||
ngm = ngm + 1
|
||||
IF (tt <= gcutms) ngms = ngms + 1
|
||||
IF (ngm > ngmx) CALL errore ('ggen', 'too many g-vectors', ngm)
|
||||
!
|
||||
|
||||
! Here map local and global g index !!!
|
||||
!
|
||||
ig_l2g( ngm ) = ng
|
||||
!
|
||||
g (1:3, ngm) = t (1:3)
|
||||
gg (ngm) = tt
|
||||
|
||||
IF (tt > eps8) THEN
|
||||
esort (ngm) = tt
|
||||
ELSE
|
||||
esort (ngm) = 0.d0
|
||||
ENDIF
|
||||
g (1:3, ngm) = i * bg (:, 1) + j * bg (:, 2) + k * bg (:, 3)
|
||||
gg (ngm) = sum(g (1:3, ngm)**2)
|
||||
|
||||
1 CONTINUE
|
||||
ENDDO
|
||||
IF (gg (ngm) <= gcutms) ngms = ngms + 1
|
||||
IF (ngm > ngmx) CALL errore ('ggen', 'too many g-vectors', ngm)
|
||||
ENDDO ngloop
|
||||
|
||||
IF (ngm.ne.ngmx) &
|
||||
CALL errore ('ggen', 'g-vectors missing !', abs(ngm - ngmx))
|
||||
!
|
||||
! reorder the g's in order of increasing magnitude. On exit
|
||||
! from hpsort esort is ordered, and nl contains the new order.
|
||||
!
|
||||
! initialize the index inside sorting routine
|
||||
IF (ngm /= ngmx) &
|
||||
CALL errore ('ggen', 'g-vectors missing !', abs(ngm - ngmx))
|
||||
|
||||
nl (1) = 0
|
||||
CALL hpsort_eps ( ngm, esort, nl, eps8 )
|
||||
!
|
||||
DEALLOCATE( esort )
|
||||
!
|
||||
! reorder also the g vectors, and nl
|
||||
!
|
||||
DO ng = 1, ngm - 1
|
||||
20 indsw = nl (ng)
|
||||
IF (indsw.ne.ng) THEN
|
||||
DO ipol = 1, 3
|
||||
swap = g (ipol, indsw)
|
||||
g (ipol, indsw) = g (ipol, nl (indsw) )
|
||||
g (ipol, nl (indsw) ) = swap
|
||||
ENDDO
|
||||
swap = gg (indsw)
|
||||
gg (indsw) = gg (nl (indsw) )
|
||||
gg (nl (indsw) ) = swap
|
||||
|
||||
!
|
||||
! Remember: ig_l2g is the index of a given G vectors in the
|
||||
! sorted global array containing all G vectors, it is used to
|
||||
! collect all wave function components
|
||||
!
|
||||
iswap = ig_l2g( indsw )
|
||||
ig_l2g( indsw ) = ig_l2g( nl(indsw) )
|
||||
ig_l2g( nl(indsw) ) = iswap
|
||||
|
||||
iswap = nl (ng)
|
||||
nl (ng) = nl (indsw)
|
||||
nl (indsw) = iswap
|
||||
|
||||
GOTO 20
|
||||
ENDIF
|
||||
|
||||
ENDDO
|
||||
!
|
||||
! here to initialize berry_phase
|
||||
! CALL berry_setup(ngm, ngm_g, nr1, nr2, nr3, mill_g)
|
||||
|
@ -252,40 +163,39 @@ SUBROUTINE ggen()
|
|||
! Now set nl and nls with the correct fft correspondence
|
||||
!
|
||||
DO ng = 1, ngm
|
||||
n1 = nint (g (1, ng) * at (1, 1) + g (2, ng) * at (2, 1) + g (3, &
|
||||
ng) * at (3, 1) ) + 1
|
||||
n1 = nint (sum(g (:, ng) * at (:, 1))) + 1
|
||||
ig1 (ng) = n1 - 1
|
||||
n1s = n1
|
||||
IF (n1.lt.1) n1 = n1 + nr1
|
||||
IF (n1s.lt.1) n1s = n1s + nr1s
|
||||
n2 = nint (g (1, ng) * at (1, 2) + g (2, ng) * at (2, 2) + g (3, &
|
||||
ng) * at (3, 2) ) + 1
|
||||
IF (n1<1) n1 = n1 + nr1
|
||||
IF (n1s<1) n1s = n1s + nr1s
|
||||
|
||||
n2 = nint (sum(g (:, ng) * at (:, 2))) + 1
|
||||
ig2 (ng) = n2 - 1
|
||||
n2s = n2
|
||||
IF (n2.lt.1) n2 = n2 + nr2
|
||||
IF (n2s.lt.1) n2s = n2s + nr2s
|
||||
n3 = nint (g (1, ng) * at (1, 3) + g (2, ng) * at (2, 3) + g (3, &
|
||||
ng) * at (3, 3) ) + 1
|
||||
IF (n2<1) n2 = n2 + nr2
|
||||
IF (n2s<1) n2s = n2s + nr2s
|
||||
|
||||
n3 = nint (sum(g (:, ng) * at (:, 3))) + 1
|
||||
ig3 (ng) = n3 - 1
|
||||
n3s = n3
|
||||
IF (n3.lt.1) n3 = n3 + nr3
|
||||
IF (n3s.lt.1) n3s = n3s + nr3s
|
||||
IF (n1.le.nr1.and.n2.le.nr2.and.n3.le.nr3) THEN
|
||||
#if defined (__PARA) && !defined (__USE_3D_FFT)
|
||||
nl (ng) = n3 + ( dfftp%isind (n1 + (n2 - 1) * nrx1) - 1) * nrx3
|
||||
IF (ng.le.ngms) nls (ng) = n3s + ( dffts%isind (n1s + (n2s - 1) &
|
||||
* nrx1s) - 1) * nrx3s
|
||||
#else
|
||||
nl (ng) = n1 + (n2 - 1) * nrx1 + (n3 - 1) * nrx1 * nrx2
|
||||
IF (ng.le.ngms) nls (ng) = n1s + (n2s - 1) * nrx1s + (n3s - 1) &
|
||||
* nrx1s * nr2s
|
||||
#endif
|
||||
ELSE
|
||||
IF (n3<1) n3 = n3 + nr3
|
||||
IF (n3s<1) n3s = n3s + nr3s
|
||||
|
||||
IF (n1>nr1 .or. n2>nr2 .or. n3>nr3) &
|
||||
CALL errore('ggen','Mesh too small?',ng)
|
||||
ENDIF
|
||||
|
||||
#if defined (__PARA) && !defined (__USE_3D_FFT)
|
||||
nl (ng) = n3 + ( dfftp%isind (n1 + (n2 - 1) * nrx1) - 1) * nrx3
|
||||
IF (ng <= ngms) &
|
||||
nls (ng) = n3s + ( dffts%isind (n1s + (n2s - 1) * nrx1s) - 1) * nrx3s
|
||||
#else
|
||||
nl (ng) = n1 + (n2 - 1) * nrx1 + (n3 - 1) * nrx1 * nrx2
|
||||
IF (ng <= ngms) &
|
||||
nls (ng) = n1s + (n2s - 1) * nrx1s + (n3s - 1) * nrx1s * nr2s
|
||||
#endif
|
||||
ENDDO
|
||||
!
|
||||
DEALLOCATE( mill_g )
|
||||
DEALLOCATE( mill_g )
|
||||
!
|
||||
! calculate number of G shells: ngl
|
||||
!
|
||||
|
@ -321,7 +231,7 @@ SUBROUTINE ggen()
|
|||
ENDIF
|
||||
ENDDO
|
||||
|
||||
IF (igl.ne.ngl) CALL errore ('setup', 'igl <> ngl', ngl)
|
||||
IF (igl /= ngl) CALL errore ('setup', 'igl <> ngl', ngl)
|
||||
|
||||
ENDIF
|
||||
|
||||
|
@ -332,11 +242,11 @@ END SUBROUTINE ggen
|
|||
!
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE index_minusg()
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
! compute indices nlm and nlms giving the correspondence
|
||||
! between the fft mesh points and -G (for gamma-only calculations)
|
||||
!
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
! compute indices nlm and nlms giving the correspondence
|
||||
! between the fft mesh points and -G (for gamma-only calculations)
|
||||
!
|
||||
USE gvect, ONLY : ngm, nr1, nr2, nr3, &
|
||||
nrx1, nrx2, nrx3, nlM, ig1, ig2, ig3
|
||||
USE gsmooth, ONLY : nr1s, nr2s, nr3s, nrx1s, nrx3s, nlsm, ngms
|
||||
|
@ -351,28 +261,30 @@ SUBROUTINE index_minusg()
|
|||
n1s = n1
|
||||
IF (n1 < 1) n1 = n1 + nr1
|
||||
IF (n1s < 1) n1s = n1s + nr1s
|
||||
|
||||
n2 = -ig2 (ng) + 1
|
||||
n2s = n2
|
||||
IF (n2 < 1) n2 = n2 + nr2
|
||||
IF (n2s < 1) n2s = n2s + nr2s
|
||||
|
||||
n3 = -ig3 (ng) + 1
|
||||
n3s = n3
|
||||
IF (n3 < 1) n3 = n3 + nr3
|
||||
IF (n3s < 1) n3s = n3s + nr3s
|
||||
IF (n1.le.nr1 .and. n2.le.nr2 .and. n3.le.nr3) THEN
|
||||
#if defined (__PARA) && !defined (__USE_3D_FFT)
|
||||
nlm(ng) = n3 + (dfftp%isind (n1 + (n2 - 1) * nrx1) - 1) * nrx3
|
||||
IF (ng.le.ngms) nlsm(ng) = n3s + (dffts%isind (n1s + (n2s - 1) &
|
||||
* nrx1s) - 1) * nrx3s
|
||||
#else
|
||||
nlm(ng) = n1 + (n2 - 1) * nrx1 + (n3 - 1) * nrx1 * nrx2
|
||||
IF (ng.le.ngms) nlsm(ng) = n1s + (n2s - 1) * nrx1s + (n3s - 1) &
|
||||
* nrx1s * nr2s
|
||||
#endif
|
||||
ELSE
|
||||
|
||||
IF (n1>nr1 .or. n2>nr2 .or. n3>nr3) THEN
|
||||
CALL errore('index_minusg','Mesh too small?',ng)
|
||||
ENDIF
|
||||
|
||||
#if defined (__PARA) && !defined (__USE_3D_FFT)
|
||||
nlm(ng) = n3 + (dfftp%isind (n1 + (n2 - 1) * nrx1) - 1) * nrx3
|
||||
IF (ng<=ngms) &
|
||||
nlsm(ng) = n3s + (dffts%isind (n1s + (n2s - 1) * nrx1s) - 1) * nrx3s
|
||||
#else
|
||||
nlm(ng) = n1 + (n2 - 1) * nrx1 + (n3 - 1) * nrx1 * nrx2
|
||||
IF (ng<=ngms) &
|
||||
nlsm(ng) = n1s + (n2s - 1) * nrx1s + (n3s - 1) * nrx1s * nr2s
|
||||
#endif
|
||||
ENDDO
|
||||
|
||||
END SUBROUTINE index_minusg
|
||||
|
||||
|
|
1
TODO
1
TODO
|
@ -216,7 +216,6 @@ TODO LIST - February 2010
|
|||
|
||||
4.5 Trouble-makers. inconsistencies, etc
|
||||
4.5.1 Negative Charge problems (see qe-forge, H on graphene)
|
||||
4.5.2 Miller indices: a potential memory bottleneck, used only for GWW!!!
|
||||
4.5.3 G-vector shells, especially in the variable-cell case, and the
|
||||
various tricks to reduce cpu by not re-calculating things that
|
||||
depend on |G| only (see e.g. qvan2). Maybe we should move to
|
||||
|
|
|
@ -36,7 +36,6 @@ sph_bes.o \
|
|||
sph_dbes.o \
|
||||
transto.o \
|
||||
date_and_tim.o \
|
||||
sort_gvec.o \
|
||||
volume.o \
|
||||
dylmr2.o \
|
||||
ylmr2.o
|
||||
|
|
127
flib/sort.f90
127
flib/sort.f90
|
@ -355,130 +355,3 @@ subroutine ihpsort (n, ia, ind)
|
|||
goto 10
|
||||
!
|
||||
end subroutine ihpsort
|
||||
!-------------------------------------------------------------------------
|
||||
subroutine kb07ad_cp90(count,n,idx)
|
||||
!-------------------------------------------------------------------------
|
||||
!
|
||||
! kb07ad handles double precision variables
|
||||
! standard fortran 66 (a verified pfort subroutine)
|
||||
! the work-space 'mark' of length 50 permits up to 2**(50/2) numbers
|
||||
! to be sorted.
|
||||
implicit none
|
||||
integer :: n, idx(*)
|
||||
real(8) :: count(*)
|
||||
real(8) :: av, x
|
||||
integer :: k1, ifk, lngth, ip, k, it, ifka, intest, iy
|
||||
integer :: i, m, la, is, idf, mloop, is1, j, mark(50)
|
||||
! set index array to original order .
|
||||
do i=1,n
|
||||
idx(i)=i
|
||||
end do
|
||||
! check that a trivial case has not been entered .
|
||||
if(n.eq.1) go to 10
|
||||
if(n.gt.1) go to 30
|
||||
write(6,20)
|
||||
20 format(///20x,'***kb07ad***no numbers to be sorted ** return to', &
|
||||
& ' calling program' )
|
||||
goto 10
|
||||
! 'm' is the length of segment which is short enough to enter
|
||||
! the final sorting routine. it may be easily changed.
|
||||
30 m=12
|
||||
! set up initial values.
|
||||
la=2
|
||||
is=1
|
||||
idf=n
|
||||
do 190 mloop=1,n
|
||||
! if segment is short enough sort with final sorting routine .
|
||||
ifka=idf-is
|
||||
if((ifka+1).gt.m)goto 70
|
||||
!********* final sorting ***
|
||||
! ( a simple bubble sort )
|
||||
is1=is+1
|
||||
do 60 j=is1,idf
|
||||
i=j
|
||||
40 if(count(i-1).lt.count(i))goto 60
|
||||
if(count(i-1).gt.count(i))goto 50
|
||||
if(idx(i-1).lt.idx(i))goto 60
|
||||
50 av=count(i-1)
|
||||
count(i-1)=count(i)
|
||||
count(i)=av
|
||||
it=idx(i-1)
|
||||
idx(i-1)=idx(i)
|
||||
idx(i)=it
|
||||
i=i-1
|
||||
if(i.gt.is)goto 40
|
||||
60 continue
|
||||
la=la-2
|
||||
goto 170
|
||||
! ******* quicksort ********
|
||||
! select the number in the central position in the segment as
|
||||
! the test number.replace it with the number from the segment's
|
||||
! highest address.
|
||||
70 iy=(is+idf)/2
|
||||
x=count(iy)
|
||||
intest=idx(iy)
|
||||
count(iy)=count(idf)
|
||||
idx(iy)=idx(idf)
|
||||
! the markers 'i' and 'ifk' are used for the beginning and end
|
||||
! of the section not so far tested against the present value
|
||||
! of x .
|
||||
k=1
|
||||
ifk=idf
|
||||
! we alternate between the outer loop that increases i and the
|
||||
! inner loop that reduces ifk, moving numbers and indices as
|
||||
! necessary, until they meet .
|
||||
do 110 i=is,idf
|
||||
if(x.gt.count(i))goto 110
|
||||
if(x.lt.count(i))goto 80
|
||||
if(intest.gt.idx(i))goto 110
|
||||
80 if(i.ge.ifk)goto 120
|
||||
count(ifk)=count(i)
|
||||
idx(ifk)=idx(i)
|
||||
k1=k
|
||||
do 100 k=k1,ifka
|
||||
ifk=idf-k
|
||||
if(count(ifk).gt.x)goto 100
|
||||
if(count(ifk).lt.x)goto 90
|
||||
if(intest.le.idx(ifk))goto 100
|
||||
90 if(i.ge.ifk)goto 130
|
||||
count(i)=count(ifk)
|
||||
idx(i)=idx(ifk)
|
||||
go to 110
|
||||
100 continue
|
||||
goto 120
|
||||
110 continue
|
||||
! return the test number to the position marked by the marker
|
||||
! which did not move last. it divides the initial segment into
|
||||
! 2 parts. any element in the first part is less than or equal
|
||||
! to any element in the second part, and they may now be sorted
|
||||
! independently .
|
||||
120 count(ifk)=x
|
||||
idx(ifk)=intest
|
||||
ip=ifk
|
||||
goto 140
|
||||
130 count(i)=x
|
||||
idx(i)=intest
|
||||
ip=i
|
||||
! store the longer subdivision in workspace.
|
||||
140 if((ip-is).gt.(idf-ip))goto 150
|
||||
mark(la)=idf
|
||||
mark(la-1)=ip+1
|
||||
idf=ip-1
|
||||
goto 160
|
||||
150 mark(la)=ip-1
|
||||
mark(la-1)=is
|
||||
is=ip+1
|
||||
! find the length of the shorter subdivision.
|
||||
160 lngth=idf-is
|
||||
if(lngth.le.0)goto 180
|
||||
! if it contains more than one element supply it with workspace .
|
||||
la=la+2
|
||||
goto 190
|
||||
170 if(la.le.0)goto 10
|
||||
! obtain the address of the shortest segment awaiting quicksort
|
||||
180 idf=mark(la)
|
||||
is=mark(la-1)
|
||||
190 continue
|
||||
10 return
|
||||
end subroutine kb07ad_cp90
|
||||
|
||||
|
|
|
@ -1,66 +0,0 @@
|
|||
!
|
||||
! 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 sort_gvec( ng, g2, mill )
|
||||
!---------------------------------------------------------------------
|
||||
!
|
||||
! first the input variables
|
||||
!
|
||||
use kinds, ONLY: DP
|
||||
use constants, ONLY: eps8
|
||||
implicit none
|
||||
INTEGER, INTENT(IN) :: ng
|
||||
REAL(DP) :: g2( * )
|
||||
INTEGER :: mill( 3, * )
|
||||
|
||||
REAL(DP), ALLOCATABLE :: gsort( : )
|
||||
INTEGER, ALLOCATABLE :: idx( : )
|
||||
INTEGER :: ig, icurr, it, im
|
||||
REAL(DP) :: gsq
|
||||
|
||||
ALLOCATE( gsort( ng ) )
|
||||
ALLOCATE( idx( ng ) )
|
||||
|
||||
DO ig = 1, ng
|
||||
IF ( g2(ig) > eps8 ) THEN
|
||||
gsort(ig) = g2(ig)
|
||||
ELSE
|
||||
gsort(ig) = 0.d0
|
||||
END IF
|
||||
END DO
|
||||
|
||||
idx(1) = 0
|
||||
CALL hpsort_eps( ng, gsort( 1 ), idx( 1 ), eps8 )
|
||||
|
||||
! ... sort indices accordingly
|
||||
DO ig = 1, ng-1
|
||||
icurr = ig
|
||||
30 IF( idx(icurr) /= ig ) THEN
|
||||
! ... swap g-vec indices
|
||||
im = mill(1,icurr); mill(1,icurr) = mill(1,idx(icurr)); mill(1,idx(icurr)) = im
|
||||
im = mill(2,icurr); mill(2,icurr) = mill(2,idx(icurr)); mill(2,idx(icurr)) = im
|
||||
im = mill(3,icurr); mill(3,icurr) = mill(3,idx(icurr)); mill(3,idx(icurr)) = im
|
||||
! ... swap modules
|
||||
gsq = g2( icurr ); g2( icurr ) = g2( idx(icurr) ); g2( idx(icurr) ) = gsq
|
||||
! ... swap indices
|
||||
it = icurr; icurr = idx(icurr); idx(it) = it
|
||||
IF( idx(icurr) == ig ) THEN
|
||||
idx(icurr) = icurr
|
||||
ELSE
|
||||
GOTO 30
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
|
||||
DEALLOCATE( gsort )
|
||||
DEALLOCATE( idx )
|
||||
|
||||
return
|
||||
end subroutine sort_gvec
|
Loading…
Reference in New Issue