scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Normalize whitespace and keywords (N. Nemec) 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6443 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
nn245 2010-02-22 08:13:22 +00:00
parent bcd81e4e90
commit e8b49314d3
16 changed files with 2080 additions and 2080 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
EE/data_structure_coarse.f90
View File

@ -11,7 +11,7 @@
!
!
!-----------------------------------------------------------------------
subroutine data_structure_coarse( lgamma, nr1, nr2, nr3, ecutwfc )
SUBROUTINE data_structure_coarse( lgamma, nr1, nr2, nr3, ecutwfc )
!-----------------------------------------------------------------------
! this routine sets the data structure for the fft arrays.
! In the parallel case distributes columns to processes, too
@ -29,26 +29,26 @@ subroutine data_structure_coarse( lgamma, nr1, nr2, nr3, ecutwfc )
!
IMPLICIT NONE
!
LOGICAL, INTENT(IN) :: lgamma
INTEGER, INTENT(IN) :: nr1, nr2, nr3
REAL(DP), INTENT(IN) :: ecutwfc
LOGICAL, INTENT(in) :: lgamma
INTEGER, INTENT(in) :: nr1, nr2, nr3
REAL(DP), INTENT(in) :: ecutwfc
!
integer :: n1, n2, n3, i1, i2, i3
INTEGER :: n1, n2, n3, i1, i2, i3
! counters on G space
!
real(DP) :: amod
! modulus of G vectors
integer, allocatable :: st(:,:), stw(:,:), sts(:,:)
INTEGER, ALLOCATABLE :: st(:,:), stw(:,:), sts(:,:)
! sticks maps
integer :: ub(3), lb(3)
INTEGER :: ub(3), lb(3)
! upper and lower bounds for maps
! cut-off for the wavefunctions
integer :: np, nps1, nq, nqs, max1, min1, max2, min2, kpoint, m1, &
INTEGER :: np, nps1, nq, nqs, max1, min1, max2, min2, kpoint, m1, &
m2, i, mc, nct_, ic, ics
!
! ... Sets the dimensions of the coarse grdi
@ -60,8 +60,8 @@ subroutine data_structure_coarse( lgamma, nr1, nr2, nr3, ecutwfc )
!
! ...
!
CALL fft_dlay_allocate( dfftc, nproc_pool, MAX( nrx1c, nrx3c ), nrx2c )
!
CALL fft_dlay_allocate( dfftc, nproc_pool, max( nrx1c, nrx3c ), nrx2c )
!
! ... Computes the number of g necessary to the calculation
!
n1 = nr1 + 1
@ -84,15 +84,15 @@ subroutine data_structure_coarse( lgamma, nr1, nr2, nr3, ecutwfc )
(i1 * bg (3, 1) + i2 * bg (3, 2) + i3 * bg (3, 3) ) **2
IF( amod <= gcutmc ) ngmc = ngmc + 1
IF( amod <= gcutmc ) stw( i2, i3 ) = 1
END DO
END DO
END DO
ENDDO
ENDDO
ENDDO
!
CALL fft_dlay_scalar( dfftc, ub, lb, nr1c, nr2c, nr3c, nrx1c, nrx2c, nrx3c, stw )
!
DEALLOCATE( stw )
!
! ...
! ...
!
ngmc_l = ngmc
ngmc_g = ngmc
@ -100,5 +100,5 @@ subroutine data_structure_coarse( lgamma, nr1, nr2, nr3, ecutwfc )
RETURN
!
end subroutine data_structure_coarse
END SUBROUTINE data_structure_coarse

380
Modules/fft_base.f90
View File

@ -9,7 +9,7 @@
!
!----------------------------------------------------------------------
! FFT base Module.
! Written by Carlo Cavazzoni
! Written by Carlo Cavazzoni
!----------------------------------------------------------------------
!
!=----------------------------------------------------------------------=!
@ -24,9 +24,9 @@
IMPLICIT NONE
! ... data structure containing all information
! ... about fft data distribution for a given
! ... about fft data distribution for a given
! ... potential grid, and its wave functions sub-grid.
TYPE ( fft_dlay_descriptor ) :: dfftp ! descriptor for dense grid
TYPE ( fft_dlay_descriptor ) :: dffts ! descriptor for smooth grid
TYPE ( fft_dlay_descriptor ) :: dfftb ! descriptor for box grids
@ -49,11 +49,11 @@
!
#if defined __NONBLOCKING_FFT
!
! NON BLOCKING SCATTER, should be better on switched network
! NON BLOCKING SCATTER, should be better on switched network
! like infiniband, ethernet, myrinet
!
!-----------------------------------------------------------------------
subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
SUBROUTINE fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!-----------------------------------------------------------------------
!
! transpose the fft grid across nodes
@ -86,15 +86,15 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
#ifdef __PARA
USE parallel_include
#endif
use mp_global, ONLY : nproc_pool, me_pool, intra_pool_comm, nproc, &
USE mp_global, ONLY : nproc_pool, me_pool, intra_pool_comm, nproc, &
my_image_id, nogrp, pgrp_comm, nplist, me_pgrp, npgrp
USE kinds, ONLY : DP
implicit none
IMPLICIT NONE
integer, intent(in) :: nrx3, nxx_, sign, ncp_ (:), npp_ (:)
complex (DP) :: f_in (nxx_), f_aux (nxx_)
logical, optional, intent(in) :: use_tg
INTEGER, INTENT(in) :: nrx3, nxx_, sign, ncp_ (:), npp_ (:)
COMPLEX (DP) :: f_in (nxx_), f_aux (nxx_)
LOGICAL, OPTIONAL, INTENT(in) :: use_tg
#ifdef __PARA
@ -114,26 +114,26 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!
! Task Groups
use_tg_ = .FALSE.
use_tg_ = .false.
IF( present( use_tg ) ) use_tg_ = use_tg
IF( PRESENT( use_tg ) ) use_tg_ = use_tg
me = me_pool + 1
!
IF( use_tg_ ) THEN
! This is the number of procs. in the plane-wave group
nprocp = npgrp
nprocp = npgrp
ipoffset = me_pgrp
gcomm = pgrp_comm
ELSE
nprocp = nproc_pool
ipoffset = me_pool
gcomm = intra_pool_comm
END IF
ENDIF
!
if ( nprocp == 1 ) return
IF ( nprocp == 1 ) RETURN
!
call start_clock ('fft_scatter')
CALL start_clock ('fft_scatter')
!
! sendcount(proc): amount of data processor "me" must send to processor
! recvcount(proc): amount of data processor "me" must receive from
@ -143,45 +143,45 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
! rdispls+1 is the beginning of data that must be received from pr
!
IF( use_tg_ ) THEN
do proc = 1, nprocp
DO proc = 1, nprocp
gproc = nplist( proc ) + 1
sendcount (proc) = npp_ ( gproc ) * ncp_ (me)
recvcount (proc) = npp_ (me) * ncp_ ( gproc )
end do
ENDDO
offset(1) = 0
do proc = 2, nprocp
DO proc = 2, nprocp
gproc = nplist( proc - 1 ) + 1
offset(proc) = offset(proc - 1) + npp_ ( gproc )
end do
ENDDO
ELSE
do proc = 1, nprocp
DO proc = 1, nprocp
sendcount (proc) = npp_ (proc) * ncp_ (me)
recvcount (proc) = npp_ (me) * ncp_ (proc)
end do
ENDDO
offset(1) = 0
do proc = 2, nprocp
DO proc = 2, nprocp
offset(proc) = offset(proc - 1) + npp_ (proc - 1)
end do
END IF
ENDDO
ENDIF
!
sdispls (1) = 0
rdispls (1) = 0
do proc = 2, nprocp
DO proc = 2, nprocp
sdispls (proc) = sdispls (proc - 1) + sendcount (proc - 1)
rdispls (proc) = rdispls (proc - 1) + recvcount (proc - 1)
enddo
ENDDO
!
ierr = 0
!
if ( sign > 0 ) then
IF ( sign > 0 ) THEN
!
! "forward" scatter from columns to planes
!
! step one: store contiguously the slices and send
!
do ip = 1, nprocp
DO ip = 1, nprocp
! the following two lines make the loop iterations different on each
! the following two lines make the loop iterations different on each
! proc in order to avoid that all procs send a msg at the same proc
! at the same time.
!
@ -198,112 +198,112 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!
SELECT CASE ( npp_ ( gproc ) )
CASE ( 1 )
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
f_aux (dest + (k - 1) ) = f_in (from + (k - 1) * nrx3 )
enddo
ENDDO
CASE ( 2 )
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
f_aux ( dest + (k - 1) * 2 - 1 + 1 ) = f_in ( from + (k - 1) * nrx3 - 1 + 1 )
f_aux ( dest + (k - 1) * 2 - 1 + 2 ) = f_in ( from + (k - 1) * nrx3 - 1 + 2 )
enddo
ENDDO
CASE ( 3 )
!$omp parallel do
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
f_aux ( dest + (k - 1) * 3 - 1 + 1 ) = f_in ( from + (k - 1) * nrx3 - 1 + 1 )
f_aux ( dest + (k - 1) * 3 - 1 + 2 ) = f_in ( from + (k - 1) * nrx3 - 1 + 2 )
f_aux ( dest + (k - 1) * 3 - 1 + 3 ) = f_in ( from + (k - 1) * nrx3 - 1 + 3 )
enddo
ENDDO
CASE ( 4 )
!$omp parallel do
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
f_aux ( dest + (k - 1) * 4 - 1 + 1 ) = f_in ( from + (k - 1) * nrx3 - 1 + 1 )
f_aux ( dest + (k - 1) * 4 - 1 + 2 ) = f_in ( from + (k - 1) * nrx3 - 1 + 2 )
f_aux ( dest + (k - 1) * 4 - 1 + 3 ) = f_in ( from + (k - 1) * nrx3 - 1 + 3 )
f_aux ( dest + (k - 1) * 4 - 1 + 4 ) = f_in ( from + (k - 1) * nrx3 - 1 + 4 )
enddo
ENDDO
CASE DEFAULT
!$omp parallel do default(shared), private(i, kdest, kfrom)
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
kdest = dest + (k - 1) * npp_ ( gproc ) - 1
kfrom = from + (k - 1) * nrx3 - 1
do i = 1, npp_ ( gproc )
DO i = 1, npp_ ( gproc )
f_aux ( kdest + i ) = f_in ( kfrom + i )
enddo
enddo
ENDDO
ENDDO
END SELECT
!
! post the non-blocking send, f_aux can't be overwritten until operation has completed
!
call mpi_isend( f_aux( sdispls( proc ) + 1 ), sendcount( proc ), MPI_DOUBLE_COMPLEX, &
CALL mpi_isend( f_aux( sdispls( proc ) + 1 ), sendcount( proc ), MPI_DOUBLE_COMPLEX, &
proc-1, me, gcomm, sh( proc ), ierr )
!
if( ABS(ierr) /= 0 ) call errore ('fft_scatter', ' forward send info<>0', ABS(ierr) )
IF( abs(ierr) /= 0 ) CALL errore ('fft_scatter', ' forward send info<>0', abs(ierr) )
!
!
end do
ENDDO
!
! step two: receive
!
do ip = 1, nprocp
DO ip = 1, nprocp
!
proc = ipoffset + 1 - ip
IF( proc < 1 ) proc = proc + nprocp
!
! now post the receive
! now post the receive
!
CALL mpi_irecv( f_in( rdispls( proc ) + 1 ), recvcount( proc ), MPI_DOUBLE_COMPLEX, &
proc-1, MPI_ANY_TAG, gcomm, rh( proc ), ierr )
!
if( ABS(ierr) /= 0 ) call errore ('fft_scatter', ' forward receive info<>0', ABS(ierr) )
IF( abs(ierr) /= 0 ) CALL errore ('fft_scatter', ' forward receive info<>0', abs(ierr) )
!
tstr( proc ) = .false.
tsts( proc ) = .false.
!
end do
ENDDO
!
! maybe useless; ensures that no garbage is present in the output
!
f_in( rdispls( nprocp ) + recvcount( nprocp ) + 1 : SIZE( f_in ) ) = 0.0_DP
f_in( rdispls( nprocp ) + recvcount( nprocp ) + 1 : size( f_in ) ) = 0.0_DP
!
lrcv = .false.
lsnd = .false.
!
! exit only when all test are true: message operation have completed
!
do while ( .not. lrcv .or. .not. lsnd )
DO WHILE ( .not. lrcv .or. .not. lsnd )
lrcv = .true.
lsnd = .true.
do proc = 1, nprocp
DO proc = 1, nprocp
!
IF( .not. tstr( proc ) ) THEN
call mpi_test( rh( proc ), tstr( proc ), istat, ierr )
END IF
CALL mpi_test( rh( proc ), tstr( proc ), istat, ierr )
ENDIF
!
IF( .not. tsts( proc ) ) THEN
call mpi_test( sh( proc ), tsts( proc ), istat, ierr )
END IF
CALL mpi_test( sh( proc ), tsts( proc ), istat, ierr )
ENDIF
!
lrcv = lrcv .and. tstr( proc )
lsnd = lsnd .and. tsts( proc )
!
end do
!
end do
ENDDO
!
ENDDO
!
else
ELSE
!
! "backward" scatter from planes to columns
!
do ip = 1, nprocp
DO ip = 1, nprocp
! post the non blocking send
proc = ipoffset + 1 + ip
IF( proc > nprocp ) proc = proc - nprocp
call mpi_isend( f_in( rdispls( proc ) + 1 ), recvcount( proc ), MPI_DOUBLE_COMPLEX, &
CALL mpi_isend( f_in( rdispls( proc ) + 1 ), recvcount( proc ), MPI_DOUBLE_COMPLEX, &
proc-1, me, gcomm, sh( proc ), ierr )
if( ABS(ierr) /= 0 ) call errore ('fft_scatter', ' backward send info<>0', ABS(ierr) )
IF( abs(ierr) /= 0 ) CALL errore ('fft_scatter', ' backward send info<>0', abs(ierr) )
! post the non blocking receive
@ -312,105 +312,105 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
CALL mpi_irecv( f_aux( sdispls( proc ) + 1 ), sendcount( proc ), MPI_DOUBLE_COMPLEX, &
proc-1, MPI_ANY_TAG, gcomm, rh(proc), ierr )
if( ABS(ierr) /= 0 ) call errore ('fft_scatter', ' backward receive info<>0', ABS(ierr) )
IF( abs(ierr) /= 0 ) CALL errore ('fft_scatter', ' backward receive info<>0', abs(ierr) )
tstr( ip ) = .false.
tsts( ip ) = .false.
end do
ENDDO
!
lrcv = .false.
lsnd = .false.
!
! exit only when all test are true: message hsve been sent and received
!
do while ( .not. lsnd )
DO WHILE ( .not. lsnd )
!
lsnd = .true.
!
do proc = 1, nprocp
DO proc = 1, nprocp
!
IF( .not. tsts( proc ) ) THEN
call mpi_test( sh( proc ), tsts( proc ), istat, ierr )
END IF
CALL mpi_test( sh( proc ), tsts( proc ), istat, ierr )
ENDIF
lsnd = lsnd .and. tsts( proc )
end do
ENDDO
end do
ENDDO
!
lrcv = .false.
!
do while ( .not. lrcv )
DO WHILE ( .not. lrcv )
!
lrcv = .true.
!
do proc = 1, nprocp
DO proc = 1, nprocp
gproc = proc
IF( use_tg_ ) gproc = nplist(proc)+1
IF( .not. tstr( proc ) ) THEN
call mpi_test( rh( proc ), tstr( proc ), istat, ierr )
CALL mpi_test( rh( proc ), tstr( proc ), istat, ierr )
IF( tstr( proc ) ) THEN
from = 1 + sdispls( proc )
dest = 1 + offset( proc )
!
!
! optimize for large parallel execution, where npp_ ( gproc ) ~ 1
!
SELECT CASE ( npp_ ( gproc ) )
CASE ( 1 )
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
f_in ( dest + (k - 1) * nrx3 ) = f_aux ( from + k - 1 )
end do
ENDDO
CASE ( 2 )
do k = 1, ncp_ ( me )
DO k = 1, ncp_ ( me )
f_in ( dest + (k - 1) * nrx3 - 1 + 1 ) = f_aux( from + (k - 1) * 2 - 1 + 1 )
f_in ( dest + (k - 1) * nrx3 - 1 + 2 ) = f_aux( from + (k - 1) * 2 - 1 + 2 )
enddo
ENDDO
CASE ( 3 )
!$omp parallel do
do k = 1, ncp_ ( me )
DO k = 1, ncp_ ( me )
f_in ( dest + (k - 1) * nrx3 - 1 + 1 ) = f_aux( from + (k - 1) * 3 - 1 + 1 )
f_in ( dest + (k - 1) * nrx3 - 1 + 2 ) = f_aux( from + (k - 1) * 3 - 1 + 2 )
f_in ( dest + (k - 1) * nrx3 - 1 + 3 ) = f_aux( from + (k - 1) * 3 - 1 + 3 )
enddo
ENDDO
CASE ( 4 )
!$omp parallel do
do k = 1, ncp_ ( me )
DO k = 1, ncp_ ( me )
f_in ( dest + (k - 1) * nrx3 - 1 + 1 ) = f_aux( from + (k - 1) * 4 - 1 + 1 )
f_in ( dest + (k - 1) * nrx3 - 1 + 2 ) = f_aux( from + (k - 1) * 4 - 1 + 2 )
f_in ( dest + (k - 1) * nrx3 - 1 + 3 ) = f_aux( from + (k - 1) * 4 - 1 + 3 )
f_in ( dest + (k - 1) * nrx3 - 1 + 4 ) = f_aux( from + (k - 1) * 4 - 1 + 4 )
enddo
ENDDO
CASE DEFAULT
!$omp parallel do default(shared), private(i, kdest, kfrom)
do k = 1, ncp_ ( me )
DO k = 1, ncp_ ( me )
kdest = dest + (k - 1) * nrx3 - 1
kfrom = from + (k - 1) * npp_ ( gproc ) - 1
do i = 1, npp_ ( gproc )
DO i = 1, npp_ ( gproc )
f_in ( kdest + i ) = f_aux( kfrom + i )
enddo
enddo
ENDDO
ENDDO
END SELECT
END IF
ENDIF
END IF
ENDIF
lrcv = lrcv .and. tstr( proc )
end do
ENDDO
end do
ENDDO
endif
ENDIF
call stop_clock ('fft_scatter')
CALL stop_clock ('fft_scatter')
#endif
@ -418,19 +418,19 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
! CALL f_hpmstop( 10 )
#endif
return
RETURN
end subroutine fft_scatter
END SUBROUTINE fft_scatter
!
!
!
#else
!
! ALLTOALL based SCATTER, should be better on network
! ALLTOALL based SCATTER, should be better on network
! with a defined topology, like on bluegene and cray machine
!
!-----------------------------------------------------------------------
subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
SUBROUTINE fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!-----------------------------------------------------------------------
!
! transpose the fft grid across nodes
@ -463,20 +463,20 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
#ifdef __PARA
USE parallel_include
#endif
use mp_global, ONLY : nproc_pool, me_pool, intra_pool_comm, nproc, &
USE mp_global, ONLY : nproc_pool, me_pool, intra_pool_comm, nproc, &
my_image_id, nogrp, pgrp_comm, nplist
USE kinds, ONLY : DP
implicit none
IMPLICIT NONE
integer, intent(in) :: nrx3, nxx_, sign, ncp_ (:), npp_ (:)
complex (DP) :: f_in (nxx_), f_aux (nxx_)
logical, optional, intent(in) :: use_tg
INTEGER, INTENT(in) :: nrx3, nxx_, sign, ncp_ (:), npp_ (:)
COMPLEX (DP) :: f_in (nxx_), f_aux (nxx_)
LOGICAL, OPTIONAL, INTENT(in) :: use_tg
#ifdef __PARA
integer :: dest, from, k, offset, proc, ierr, me, nprocp, gproc, gcomm, i, kdest, kfrom
integer :: sendcount (nproc_pool), sdispls (nproc_pool), recvcount (nproc_pool), rdispls (nproc_pool)
INTEGER :: dest, from, k, offset, proc, ierr, me, nprocp, gproc, gcomm, i, kdest, kfrom
INTEGER :: sendcount (nproc_pool), sdispls (nproc_pool), recvcount (nproc_pool), rdispls (nproc_pool)
!
LOGICAL :: use_tg_
@ -487,22 +487,22 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!
! Task Groups
use_tg_ = .FALSE.
use_tg_ = .false.
IF( present( use_tg ) ) use_tg_ = use_tg
IF( PRESENT( use_tg ) ) use_tg_ = use_tg
me = me_pool + 1
!
IF( use_tg_ ) THEN
! This is the number of procs. in the plane-wave group
nprocp = nproc_pool / nogrp
nprocp = nproc_pool / nogrp
ELSE
nprocp = nproc_pool
END IF
ENDIF
!
if (nprocp.eq.1) return
IF (nprocp.eq.1) RETURN
!
call start_clock ('fft_scatter')
CALL start_clock ('fft_scatter')
!
! sendcount(proc): amount of data processor "me" must send to processor
! recvcount(proc): amount of data processor "me" must receive from
@ -512,28 +512,28 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!
!
IF( use_tg_ ) THEN
do proc = 1, nprocp
DO proc = 1, nprocp
gproc = nplist( proc ) + 1
sendcount (proc) = npp_ ( gproc ) * ncp_ (me)
recvcount (proc) = npp_ (me) * ncp_ ( gproc )
end do
ENDDO
ELSE
do proc = 1, nprocp
DO proc = 1, nprocp
sendcount (proc) = npp_ (proc) * ncp_ (me)
recvcount (proc) = npp_ (me) * ncp_ (proc)
end do
END IF
ENDDO
ENDIF
!
sdispls (1) = 0
rdispls (1) = 0
do proc = 2, nprocp
DO proc = 2, nprocp
sdispls (proc) = sdispls (proc - 1) + sendcount (proc - 1)
rdispls (proc) = rdispls (proc - 1) + recvcount (proc - 1)
enddo
ENDDO
!
ierr = 0
if (sign.gt.0) then
IF (sign.gt.0) THEN
!
! "forward" scatter from columns to planes
!
@ -541,37 +541,37 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!
offset = 1
do proc = 1, nprocp
DO proc = 1, nprocp
from = offset
dest = 1 + sdispls (proc)
IF( use_tg_ ) THEN
gproc = nplist(proc)+1
ELSE
gproc = proc
END IF
ENDIF
!
! optimize for large parallel execution, where npp_ ( gproc ) ~ 1
!
IF( npp_ ( gproc ) > 128 ) THEN
do k = 1, ncp_ (me)
call DCOPY (2 * npp_ ( gproc ), f_in (from + (k - 1) * nrx3), &
DO k = 1, ncp_ (me)
CALL DCOPY (2 * npp_ ( gproc ), f_in (from + (k - 1) * nrx3), &
1, f_aux (dest + (k - 1) * npp_ ( gproc ) ), 1)
enddo
ELSE IF( npp_ ( gproc ) == 1 ) THEN
do k = 1, ncp_ (me)
ENDDO
ELSEIF( npp_ ( gproc ) == 1 ) THEN
DO k = 1, ncp_ (me)
f_aux (dest + (k - 1) ) = f_in (from + (k - 1) * nrx3 )
enddo
ENDDO
ELSE
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
kdest = dest + (k - 1) * npp_ ( gproc ) - 1
kfrom = from + (k - 1) * nrx3 - 1
do i = 1, npp_ ( gproc )
DO i = 1, npp_ ( gproc )
f_aux ( kdest + i ) = f_in ( kfrom + i )
enddo
enddo
END IF
ENDDO
ENDDO
ENDIF
offset = offset + npp_ ( gproc )
enddo
ENDDO
!
! maybe useless; ensures that no garbage is present in the output
!
@ -583,16 +583,16 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
gcomm = pgrp_comm
ELSE
gcomm = intra_pool_comm
END IF
ENDIF
call mpi_barrier (gcomm, ierr) ! why barrier? for buggy openmpi over ib
CALL mpi_barrier (gcomm, ierr) ! why barrier? for buggy openmpi over ib
call mpi_alltoallv (f_aux(1), sendcount, sdispls, MPI_DOUBLE_COMPLEX, f_in(1), &
CALL mpi_alltoallv (f_aux(1), sendcount, sdispls, MPI_DOUBLE_COMPLEX, f_in(1), &
recvcount, rdispls, MPI_DOUBLE_COMPLEX, gcomm, ierr)
if( ABS(ierr) /= 0 ) call errore ('fft_scatter', 'info<>0', ABS(ierr) )
IF( abs(ierr) /= 0 ) CALL errore ('fft_scatter', 'info<>0', abs(ierr) )
!
else
ELSE
!
! "backward" scatter from planes to columns
!
@ -602,14 +602,14 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
gcomm = pgrp_comm
ELSE
gcomm = intra_pool_comm
END IF
ENDIF
call mpi_barrier (gcomm, ierr) ! why barrier? for buggy openmpi over ib
CALL mpi_barrier (gcomm, ierr) ! why barrier? for buggy openmpi over ib
call mpi_alltoallv (f_in(1), recvcount, rdispls, MPI_DOUBLE_COMPLEX, f_aux(1), &
CALL mpi_alltoallv (f_in(1), recvcount, rdispls, MPI_DOUBLE_COMPLEX, f_aux(1), &
sendcount, sdispls, MPI_DOUBLE_COMPLEX, gcomm, ierr)
if( ABS(ierr) /= 0 ) call errore ('fft_scatter', 'info<>0', ABS(ierr) )
IF( abs(ierr) /= 0 ) CALL errore ('fft_scatter', 'info<>0', abs(ierr) )
!
! step one: store contiguously the columns
!
@ -617,43 +617,43 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
!
offset = 1
!
do proc = 1, nprocp
DO proc = 1, nprocp
from = 1 + sdispls (proc)
dest = offset
IF( use_tg_ ) THEN
gproc = nplist(proc)+1
ELSE
gproc = proc
END IF
!
ENDIF
!
! optimize for large parallel execution, where npp_ ( gproc ) ~ 1
!
IF( npp_ ( gproc ) > 128 ) THEN
do k = 1, ncp_ (me)
call DCOPY ( 2 * npp_ ( gproc ), f_aux (from + (k - 1) * npp_ ( gproc ) ), 1, &
DO k = 1, ncp_ (me)
CALL DCOPY ( 2 * npp_ ( gproc ), f_aux (from + (k - 1) * npp_ ( gproc ) ), 1, &
f_in (dest + (k - 1) * nrx3 ), 1 )
enddo
ELSE IF ( npp_ ( gproc ) == 1 ) THEN
do k = 1, ncp_ (me)
ENDDO
ELSEIF ( npp_ ( gproc ) == 1 ) THEN
DO k = 1, ncp_ (me)
f_in ( dest + (k - 1) * nrx3 ) = f_aux ( from + (k - 1) )
end do
ENDDO
ELSE
do k = 1, ncp_ (me)
DO k = 1, ncp_ (me)
kdest = dest + (k - 1) * nrx3 - 1
kfrom = from + (k - 1) * npp_ ( gproc ) - 1
do i = 1, npp_ ( gproc )
DO i = 1, npp_ ( gproc )
f_in ( kdest + i ) = f_aux( kfrom + i )
enddo
enddo
END IF
ENDDO
ENDDO
ENDIF
!
offset = offset + npp_ ( gproc )
!
enddo
ENDDO
endif
ENDIF
call stop_clock ('fft_scatter')
CALL stop_clock ('fft_scatter')
#endif
@ -661,9 +661,9 @@ subroutine fft_scatter ( f_in, nrx3, nxx_, f_aux, ncp_, npp_, sign, use_tg )
! CALL f_hpmstop( 10 )
#endif
return
RETURN
end subroutine fft_scatter
END SUBROUTINE fft_scatter
#endif
@ -689,8 +689,8 @@ SUBROUTINE grid_gather( f_in, f_out )
INTEGER :: proc, info
INTEGER :: displs(0:nproc_pool-1), recvcount(0:nproc_pool-1)
!
IF( SIZE( f_in ) < dfftp%nnr ) &
CALL errore( ' grid_gather ', ' f_in too small ', dfftp%nnr - SIZE( f_in ) )
IF( size( f_in ) < dfftp%nnr ) &
CALL errore( ' grid_gather ', ' f_in too small ', dfftp%nnr - size( f_in ) )
!
CALL start_clock( 'gather' )
!
@ -706,11 +706,11 @@ SUBROUTINE grid_gather( f_in, f_out )
!
displs(proc) = displs(proc-1) + recvcount(proc-1)
!
END IF
ENDIF
!
END DO
ENDDO
!
info = SIZE( f_out ) - displs( nproc_pool - 1 ) - recvcount( nproc_pool - 1 )
info = size( f_out ) - displs( nproc_pool - 1 ) - recvcount( nproc_pool - 1 )
!
IF( info < 0 ) &
CALL errore( ' grid_gather ', ' f_out too small ', -info )
@ -755,8 +755,8 @@ SUBROUTINE grid_scatter( f_in, f_out )
INTEGER :: proc, info
INTEGER :: displs(0:nproc_pool-1), sendcount(0:nproc_pool-1)
!
IF( SIZE( f_out ) < dfftp%nnr ) &
CALL errore( ' grid_scatter ', ' f_out too small ', dfftp%nnr - SIZE( f_in ) )
IF( size( f_out ) < dfftp%nnr ) &
CALL errore( ' grid_scatter ', ' f_out too small ', dfftp%nnr - size( f_in ) )
!
CALL start_clock( 'scatter' )
!
@ -772,11 +772,11 @@ SUBROUTINE grid_scatter( f_in, f_out )
!
displs(proc) = displs(proc-1) + sendcount(proc-1)
!
END IF
ENDIF
!
END DO
ENDDO
!
info = SIZE( f_in ) - displs( nproc_pool - 1 ) - sendcount( nproc_pool - 1 )
info = size( f_in ) - displs( nproc_pool - 1 ) - sendcount( nproc_pool - 1 )
!
IF( info < 0 ) &
CALL errore( ' grid_scatter ', ' f_in too small ', -info )
@ -812,13 +812,13 @@ SUBROUTINE cgather_sym( f_in, f_out )
USE mp_global, ONLY : intra_pool_comm, intra_image_comm, &
nproc_pool, me_pool
USE mp, ONLY : mp_barrier
USE parallel_include
USE parallel_include
!
IMPLICIT NONE
!
COMPLEX(DP) :: f_in( : ), f_out(:)
!
#if defined (__PARA)
#if defined (__PARA)
!
INTEGER :: proc, info
INTEGER :: displs(0:nproc_pool-1), recvcount(0:nproc_pool-1)
@ -838,9 +838,9 @@ SUBROUTINE cgather_sym( f_in, f_out )
!
displs(proc) = displs(proc-1) + recvcount(proc-1)
!
END IF
ENDIF
!
END DO
ENDDO
!
CALL mp_barrier( intra_pool_comm )
!
@ -874,13 +874,13 @@ SUBROUTINE cgather_smooth ( f_in, f_out )
USE mp_global, ONLY : intra_pool_comm, nproc_pool, me_pool, root_pool
USE mp, ONLY : mp_barrier
USE kinds, ONLY : DP
USE parallel_include
USE parallel_include
!
IMPLICIT NONE
!
COMPLEX(DP) :: f_in(:), f_out(:)
!
#if defined (__PARA)
#if defined (__PARA)
!
INTEGER :: proc, info
INTEGER :: displs(0:nproc_pool-1), recvcount(0:nproc_pool-1)
@ -900,9 +900,9 @@ SUBROUTINE cgather_smooth ( f_in, f_out )
!
displs(proc) = displs(proc-1) + recvcount(proc-1)
!
END IF
ENDIF
!
END DO
ENDDO
!
CALL mp_barrier( intra_pool_comm )
!
@ -935,13 +935,13 @@ SUBROUTINE cscatter_sym( f_in, f_out )
me_pool, root_pool
USE mp, ONLY : mp_barrier
USE kinds, ONLY : DP
USE parallel_include
USE parallel_include
!
IMPLICIT NONE
!
COMPLEX(DP) :: f_in(:), f_out(:)
!
#if defined (__PARA)
#if defined (__PARA)
!
INTEGER :: proc, info
INTEGER :: displs(0:nproc_pool-1), sendcount(0:nproc_pool-1)
@ -961,12 +961,12 @@ SUBROUTINE cscatter_sym( f_in, f_out )
!
displs(proc) = displs(proc-1) + sendcount(proc-1)
!
END IF
ENDIF
!
END DO
ENDDO
!
CALL mp_barrier( intra_pool_comm )
!
!
CALL MPI_SCATTERV( f_in, sendcount, displs, MPI_DOUBLE_PRECISION, &
f_out, sendcount(me_pool), MPI_DOUBLE_PRECISION, &
root_pool, intra_pool_comm, info )
@ -997,13 +997,13 @@ SUBROUTINE cscatter_smooth( f_in, f_out )
me_pool, root_pool
USE mp, ONLY : mp_barrier
USE kinds, ONLY : DP
USE parallel_include
USE parallel_include
!
IMPLICIT NONE
!
COMPLEX(DP) :: f_in(:), f_out(:)
!
#if defined (__PARA)
#if defined (__PARA)
!
INTEGER :: proc, info
INTEGER :: displs(0:nproc_pool-1), sendcount(0:nproc_pool-1)
@ -1023,12 +1023,12 @@ SUBROUTINE cscatter_smooth( f_in, f_out )
!
displs(proc) = displs(proc-1) + sendcount(proc-1)
!
END IF
ENDIF
!
END DO
ENDDO
!
CALL mp_barrier( intra_pool_comm )
!
!
CALL MPI_SCATTERV( f_in, sendcount, displs, MPI_DOUBLE_PRECISION, &
f_out, sendcount(me_pool), MPI_DOUBLE_PRECISION, &
root_pool, intra_pool_comm, info )

174
Modules/fft_parallel.f90
View File

@ -60,7 +60,7 @@ SUBROUTINE tg_cft3s( f, dfft, isgn, use_task_groups )
USE fft_scalar, ONLY : cft_1z, cft_2xy
USE fft_base, ONLY : fft_scatter
USE kinds, ONLY : DP
USE mp_global, only : me_pool, nproc_pool, ogrp_comm, npgrp, nogrp, &
USE mp_global, ONLY : me_pool, nproc_pool, ogrp_comm, npgrp, nogrp, &
intra_pool_comm, nolist, nplist
USE fft_types, ONLY : fft_dlay_descriptor
USE parallel_include
@ -68,11 +68,11 @@ SUBROUTINE tg_cft3s( f, dfft, isgn, use_task_groups )
!
IMPLICIT NONE
!
COMPLEX(DP), INTENT(INOUT) :: f( : ) ! array containing data to be transformed
type (fft_dlay_descriptor), intent(in) :: dfft
COMPLEX(DP), INTENT(inout) :: f( : ) ! array containing data to be transformed
TYPE (fft_dlay_descriptor), INTENT(in) :: dfft
! descriptor of fft data layout
INTEGER, INTENT(IN) :: isgn ! fft direction
LOGICAL, OPTIONAL, INTENT(IN) :: use_task_groups
INTEGER, INTENT(in) :: isgn ! fft direction
LOGICAL, OPTIONAL, INTENT(in) :: use_task_groups
! specify if you want to use task groups parallelization
!
INTEGER :: me_p
@ -84,13 +84,13 @@ SUBROUTINE tg_cft3s( f, dfft, isgn, use_task_groups )
!
CALL start_clock( 'cft3s' )
!
IF( PRESENT( use_task_groups ) ) THEN
IF( present( use_task_groups ) ) THEN
use_tg = use_task_groups
ELSE
use_tg = .FALSE.
END IF
use_tg = .false.
ENDIF
!
IF( use_tg .AND. .NOT. dfft%have_task_groups ) &
IF( use_tg .and. .not. dfft%have_task_groups ) &
CALL errore( ' tg_cft3s ', ' call requiring task groups for a descriptor without task groups ', 1 )
!
n1 = dfft%nr1
@ -105,7 +105,7 @@ SUBROUTINE tg_cft3s( f, dfft, isgn, use_task_groups )
ALLOCATE( YF ( nogrp * dfft%nnrx ) )
ELSE
ALLOCATE( aux( dfft%nnrx ) )
END IF
ENDIF
!
me_p = me_pool + 1
!
@ -116,7 +116,7 @@ SUBROUTINE tg_cft3s( f, dfft, isgn, use_task_groups )
IF( use_tg ) &
CALL errore( ' tg_cfft ', ' task groups on large mesh not implemented ', 1 )
!
call cft_1z( f, dfft%nsp( me_p ), n3, nx3, isgn, aux )
CALL cft_1z( f, dfft%nsp( me_p ), n3, nx3, isgn, aux )
!
planes = dfft%iplp
!
@ -127,12 +127,12 @@ SUBROUTINE tg_cft3s( f, dfft, isgn, use_task_groups )
IF( use_tg ) THEN
CALL cft_1z( yf, dfft%tg_nsw( me_p ), n3, nx3, isgn, aux )
ELSE
call cft_1z( f, dfft%nsw( me_p ), n3, nx3, isgn, aux )
END IF
CALL cft_1z( f, dfft%nsw( me_p ), n3, nx3, isgn, aux )
ENDIF
!
planes = dfft%iplw
!
END IF
ENDIF
!
CALL fw_scatter( isgn ) ! forwart scatter from stick to planes
!
@ -140,54 +140,54 @@ SUBROUTINE tg_cft3s( f, dfft, isgn, use_task_groups )
CALL cft_2xy( f, dfft%tg_npp( me_p ), n1, n2, nx1, nx2, isgn, planes )
ELSE
CALL cft_2xy( f, dfft%npp( me_p ), n1, n2, nx1, nx2, isgn, planes )
END IF
ENDIF
!
ELSE
!
if ( isgn /= -2 ) then
IF ( isgn /= -2 ) THEN
!
IF( use_tg ) &
CALL errore( ' tg_cfft ', ' task groups on large mesh not implemented ', 1 )
!
planes = dfft%iplp
!
else
ELSE
!
planes = dfft%iplw
!
endif
ENDIF
IF( use_tg ) THEN
CALL cft_2xy( f, dfft%tg_npp( me_p ), n1, n2, nx1, nx2, isgn, planes )
ELSE
call cft_2xy( f, dfft%npp( me_p ), n1, n2, nx1, nx2, isgn, planes)
END IF
CALL cft_2xy( f, dfft%npp( me_p ), n1, n2, nx1, nx2, isgn, planes)
ENDIF
!
CALL bw_scatter( isgn )
!
IF ( isgn /= -2 ) THEN
!
call cft_1z( aux, dfft%nsp( me_p ), n3, nx3, isgn, f )
CALL cft_1z( aux, dfft%nsp( me_p ), n3, nx3, isgn, f )
!
ELSE
!
IF( use_tg ) THEN
CALL cft_1z( aux, dfft%tg_nsw( me_p ), n3, nx3, isgn, yf )
ELSE
call cft_1z( aux, dfft%nsw( me_p ), n3, nx3, isgn, f )
END IF
CALL cft_1z( aux, dfft%nsw( me_p ), n3, nx3, isgn, f )
ENDIF
!
CALL unpack_group_sticks()
!
END IF
ENDIF
!
END IF
ENDIF
!
DEALLOCATE( aux )
!
IF( use_tg ) THEN
DEALLOCATE( yf )
END IF
ENDIF
!
CALL stop_clock( 'cft3s' )
!
@ -200,14 +200,14 @@ CONTAINS
INTEGER :: ierr
!
if( .NOT. use_tg ) return
IF( .not. use_tg ) RETURN
!
IF( dfft%tg_rdsp(nogrp) + dfft%tg_rcv(nogrp) > SIZE( yf ) ) THEN
IF( dfft%tg_rdsp(nogrp) + dfft%tg_rcv(nogrp) > size( yf ) ) THEN
CALL errore( ' tg_cfft ', ' inconsistent size ', 1 )
END IF
IF( dfft%tg_psdsp(nogrp) + dfft%tg_snd(nogrp) > SIZE( f ) ) THEN
ENDIF
IF( dfft%tg_psdsp(nogrp) + dfft%tg_snd(nogrp) > size( f ) ) THEN
CALL errore( ' tg_cfft ', ' inconsistent size ', 2 )
END IF
ENDIF
CALL start_clock( 'ALLTOALL' )
!
@ -219,8 +219,8 @@ CONTAINS
CALL MPI_ALLTOALLV( f(1), dfft%tg_snd, dfft%tg_psdsp, MPI_DOUBLE_COMPLEX, yf(1), dfft%tg_rcv, &
& dfft%tg_rdsp, MPI_DOUBLE_COMPLEX, ogrp_comm, IERR)
IF( ierr /= 0 ) THEN
CALL errore( ' tg_cfft ', ' alltoall error 1 ', ABS(ierr) )
END IF
CALL errore( ' tg_cfft ', ' alltoall error 1 ', abs(ierr) )
ENDIF
#endif
@ -239,14 +239,14 @@ CONTAINS
!
INTEGER :: ierr
!
if( .NOT. use_tg ) return
IF( .not. use_tg ) RETURN
!
IF( dfft%tg_usdsp(nogrp) + dfft%tg_snd(nogrp) > SIZE( f ) ) THEN
IF( dfft%tg_usdsp(nogrp) + dfft%tg_snd(nogrp) > size( f ) ) THEN
CALL errore( ' tg_cfft ', ' inconsistent size ', 3 )
END IF
IF( dfft%tg_rdsp(nogrp) + dfft%tg_rcv(nogrp) > SIZE( yf ) ) THEN
ENDIF
IF( dfft%tg_rdsp(nogrp) + dfft%tg_rcv(nogrp) > size( yf ) ) THEN
CALL errore( ' tg_cfft ', ' inconsistent size ', 4 )
END IF
ENDIF
CALL start_clock( 'ALLTOALL' )
@ -255,8 +255,8 @@ CONTAINS
dfft%tg_rcv, dfft%tg_rdsp, MPI_DOUBLE_COMPLEX, f(1), &
dfft%tg_snd, dfft%tg_usdsp, MPI_DOUBLE_COMPLEX, ogrp_comm, IERR)
IF( ierr /= 0 ) THEN
CALL errore( ' tg_cfft ', ' alltoall error 2 ', ABS(ierr) )
END IF
CALL errore( ' tg_cfft ', ' alltoall error 2 ', abs(ierr) )
ENDIF
#endif
CALL stop_clock( 'ALLTOALL' )
@ -279,9 +279,9 @@ CONTAINS
!dfft%npp: number of planes per processor
!
!
use fft_base, only: fft_scatter
USE fft_base, ONLY: fft_scatter
!
INTEGER, INTENT(IN) :: iopt
INTEGER, INTENT(in) :: iopt
INTEGER :: nppx, ip, nnp, npp, ii, i, mc, j, ioff, it
!
IF( iopt == 2 ) THEN
@ -301,67 +301,67 @@ CONTAINS
npp = dfft%npp( me_p )
nnp = dfft%nnp
!
call fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsw, dfft%npp, iopt )
CALL fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsw, dfft%npp, iopt )
!
END IF
ENDIF
!
!
!$omp parallel default(shared), private( ii, mc, j, i, ioff, ip, it )
!$omp do
do i = 1, SIZE( f )
DO i = 1, size( f )
f(i) = (0.d0, 0.d0)
end do
ENDDO
!
ii = 0
!
do ip = 1, nproc_pool
DO ip = 1, nproc_pool
!
ioff = dfft%iss( ip )
!
!$omp do
do i = 1, dfft%nsw( ip )
DO i = 1, dfft%nsw( ip )
!
mc = dfft%ismap( i + ioff )
!
it = ( ii + i - 1 ) * nppx
!
do j = 1, npp
DO j = 1, npp
f( mc + ( j - 1 ) * nnp ) = aux( j + it )
end do
ENDDO
!
end do
ENDDO
!
ii = ii + dfft%nsw( ip )
!
end do
ENDDO
!$omp end parallel
!
ELSE IF( iopt == 1 ) THEN
ELSEIF( iopt == 1 ) THEN
!
if ( nproc_pool == 1 ) then
IF ( nproc_pool == 1 ) THEN
nppx = dfft%nr3x
else
ELSE
nppx = dfft%npp( me_p )
end if
ENDIF
!
call fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsp, dfft%npp, iopt )
CALL fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsp, dfft%npp, iopt )
!
!$omp parallel default(shared)
!$omp do
do i = 1, SIZE(f)
!$omp do
DO i = 1, size(f)
f(i) = (0.d0, 0.d0)
end do
ENDDO
!
!$omp do private(mc,j)
do i = 1, dfft%nst
DO i = 1, dfft%nst
mc = dfft%ismap( i )
do j = 1, dfft%npp( me_p )
DO j = 1, dfft%npp( me_p )
f( mc + ( j - 1 ) * dfft%nnp ) = aux( j + ( i - 1 ) * nppx )
end do
end do
ENDDO
ENDDO
!$omp end parallel
!
END IF
ENDIF
!
RETURN
END SUBROUTINE fw_scatter
@ -370,9 +370,9 @@ CONTAINS
SUBROUTINE bw_scatter( iopt )
!
use fft_base, only: fft_scatter
USE fft_base, ONLY: fft_scatter
!
INTEGER, INTENT(IN) :: iopt
INTEGER, INTENT(in) :: iopt
INTEGER :: nppx, ip, nnp, npp, ii, i, mc, j, it
!
!
@ -391,22 +391,22 @@ CONTAINS
npp = dfft%npp( me_p )
nnp = dfft%nnp
!
END IF
ENDIF
!$omp parallel default(shared), private( mc, j, i, ii, ip, it )
ii = 0
do ip = 1, nproc_pool
DO ip = 1, nproc_pool
!$omp do
do i = 1, dfft%nsw( ip )
DO i = 1, dfft%nsw( ip )
mc = dfft%ismap( i + dfft%iss( ip ) )
it = (ii + i - 1)*nppx
do j = 1, npp
DO j = 1, npp
aux( j + it ) = f( mc + ( j - 1 ) * nnp )
end do
end do
ENDDO
ENDDO
ii = ii + dfft%nsw( ip )
end do
ENDDO
!$omp end parallel
!
IF( use_tg ) THEN
@ -415,30 +415,30 @@ CONTAINS
!
ELSE
!
call fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsw, dfft%npp, iopt )
CALL fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsw, dfft%npp, iopt )
!
END IF
ENDIF
!
ELSE IF( iopt == -1 ) THEN
ELSEIF( iopt == -1 ) THEN
!
if ( nproc_pool == 1 ) then
IF ( nproc_pool == 1 ) THEN
nppx = dfft%nr3x
else
ELSE
nppx = dfft%npp( me_p )
end if
ENDIF
!$omp parallel default(shared), private( mc, j, i )
!$omp do
do i = 1, dfft%nst
!$omp do
DO i = 1, dfft%nst
mc = dfft%ismap( i )
do j = 1, dfft%npp( me_p )
DO j = 1, dfft%npp( me_p )
aux( j + ( i - 1 ) * nppx ) = f( mc + ( j - 1 ) * dfft%nnp )
end do
end do
ENDDO
ENDDO
!$omp end parallel
!
call fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsp, dfft%npp, iopt )
CALL fft_scatter( aux, nx3, dfft%nnr, f, dfft%nsp, dfft%npp, iopt )
!
END IF
ENDIF
!
RETURN
END SUBROUTINE bw_scatter

292
Modules/fft_types.f90
View File

@ -17,16 +17,16 @@ MODULE fft_types
TYPE fft_dlay_descriptor
INTEGER :: nst ! total number of sticks
INTEGER, POINTER :: nsp(:) ! number of sticks per processor ( potential )
! using proc index starting from 1 !!
! using proc index starting from 1 !!
! on proc mpime -> nsp( mpime + 1 )
INTEGER, POINTER :: nsw(:) ! number of sticks per processor ( wave func )
! using proc index as above
INTEGER :: nr1 !
INTEGER :: nr2 ! effective FFT dimensions
INTEGER :: nr3 !
INTEGER :: nr1x !
INTEGER :: nr3 !
INTEGER :: nr1x !
INTEGER :: nr2x ! FFT grids leading dimensions
INTEGER :: nr3x !
INTEGER :: nr3x !
INTEGER :: npl ! number of "Z" planes for this processor = npp( mpime + 1 )
INTEGER :: nnp ! number of 0 and non 0 sticks in a plane ( ~nr1*nr2/nproc )
INTEGER :: nnr ! local number of FFT grid elements ( ~nr1*nr2*nr3/proc )
@ -75,7 +75,7 @@ CONTAINS
SUBROUTINE fft_dlay_allocate( desc, nproc, nx, ny )
TYPE (fft_dlay_descriptor) :: desc
INTEGER, INTENT(IN) :: nproc, nx, ny
INTEGER, INTENT(in) :: nproc, nx, ny
ALLOCATE( desc%nsp( nproc ) )
ALLOCATE( desc%nsw( nproc ) )
ALLOCATE( desc%ngl( nproc ) )
@ -102,7 +102,7 @@ CONTAINS
desc%id = 0
desc%have_task_groups = .FALSE.
desc%have_task_groups = .false.
NULLIFY( desc%tg_nsw )
NULLIFY( desc%tg_npp )
NULLIFY( desc%tg_snd )
@ -116,35 +116,35 @@ CONTAINS
SUBROUTINE fft_dlay_deallocate( desc )
TYPE (fft_dlay_descriptor) :: desc
IF ( ASSOCIATED( desc%nsp ) ) DEALLOCATE( desc%nsp )
IF ( ASSOCIATED( desc%nsw ) ) DEALLOCATE( desc%nsw )
IF ( ASSOCIATED( desc%ngl ) ) DEALLOCATE( desc%ngl )
IF ( ASSOCIATED( desc%nwl ) ) DEALLOCATE( desc%nwl )
IF ( ASSOCIATED( desc%npp ) ) DEALLOCATE( desc%npp )
IF ( ASSOCIATED( desc%ipp ) ) DEALLOCATE( desc%ipp )
IF ( ASSOCIATED( desc%iss ) ) DEALLOCATE( desc%iss )
IF ( ASSOCIATED( desc%isind ) ) DEALLOCATE( desc%isind )
IF ( ASSOCIATED( desc%ismap ) ) DEALLOCATE( desc%ismap )
IF ( ASSOCIATED( desc%iplp ) ) DEALLOCATE( desc%iplp )
IF ( ASSOCIATED( desc%iplw ) ) DEALLOCATE( desc%iplw )
IF ( associated( desc%nsp ) ) DEALLOCATE( desc%nsp )
IF ( associated( desc%nsw ) ) DEALLOCATE( desc%nsw )
IF ( associated( desc%ngl ) ) DEALLOCATE( desc%ngl )
IF ( associated( desc%nwl ) ) DEALLOCATE( desc%nwl )
IF ( associated( desc%npp ) ) DEALLOCATE( desc%npp )
IF ( associated( desc%ipp ) ) DEALLOCATE( desc%ipp )
IF ( associated( desc%iss ) ) DEALLOCATE( desc%iss )
IF ( associated( desc%isind ) ) DEALLOCATE( desc%isind )
IF ( associated( desc%ismap ) ) DEALLOCATE( desc%ismap )
IF ( associated( desc%iplp ) ) DEALLOCATE( desc%iplp )
IF ( associated( desc%iplw ) ) DEALLOCATE( desc%iplw )
desc%id = 0
IF( desc%have_task_groups ) THEN
IF ( ASSOCIATED( desc%tg_nsw ) ) DEALLOCATE( desc%tg_nsw )
IF ( ASSOCIATED( desc%tg_npp ) ) DEALLOCATE( desc%tg_npp )
IF ( ASSOCIATED( desc%tg_snd ) ) DEALLOCATE( desc%tg_snd )
IF ( ASSOCIATED( desc%tg_rcv ) ) DEALLOCATE( desc%tg_rcv )
IF ( ASSOCIATED( desc%tg_psdsp ) ) DEALLOCATE( desc%tg_psdsp )
IF ( ASSOCIATED( desc%tg_usdsp ) ) DEALLOCATE( desc%tg_usdsp )
IF ( ASSOCIATED( desc%tg_rdsp ) ) DEALLOCATE( desc%tg_rdsp )
END IF
desc%have_task_groups = .FALSE.
IF ( associated( desc%tg_nsw ) ) DEALLOCATE( desc%tg_nsw )
IF ( associated( desc%tg_npp ) ) DEALLOCATE( desc%tg_npp )
IF ( associated( desc%tg_snd ) ) DEALLOCATE( desc%tg_snd )
IF ( associated( desc%tg_rcv ) ) DEALLOCATE( desc%tg_rcv )
IF ( associated( desc%tg_psdsp ) ) DEALLOCATE( desc%tg_psdsp )
IF ( associated( desc%tg_usdsp ) ) DEALLOCATE( desc%tg_usdsp )
IF ( associated( desc%tg_rdsp ) ) DEALLOCATE( desc%tg_rdsp )
ENDIF
desc%have_task_groups = .false.
END SUBROUTINE fft_dlay_deallocate
!=----------------------------------------------------------------------------=!
SUBROUTINE fft_box_allocate( desc, nproc, nat )
TYPE (fft_dlay_descriptor) :: desc
INTEGER, INTENT(IN) :: nat, nproc
INTEGER, INTENT(in) :: nat, nproc
ALLOCATE( desc%irb( 3, nat ) )
ALLOCATE( desc%imin3( nat ) )
ALLOCATE( desc%imax3( nat ) )
@ -157,18 +157,18 @@ CONTAINS
desc%npp = 0
desc%ipp = 0
desc%np3 = 0
desc%have_task_groups = .FALSE.
desc%have_task_groups = .false.
END SUBROUTINE fft_box_allocate
SUBROUTINE fft_box_deallocate( desc )
TYPE (fft_dlay_descriptor) :: desc
IF( ASSOCIATED( desc%irb ) ) DEALLOCATE( desc%irb )
IF( ASSOCIATED( desc%imin3 ) ) DEALLOCATE( desc%imin3 )
IF( ASSOCIATED( desc%imax3 ) ) DEALLOCATE( desc%imax3 )
IF( ASSOCIATED( desc%npp ) ) DEALLOCATE( desc%npp )
IF( ASSOCIATED( desc%ipp ) ) DEALLOCATE( desc%ipp )
IF( ASSOCIATED( desc%np3 ) ) DEALLOCATE( desc%np3 )
desc%have_task_groups = .FALSE.
IF( associated( desc%irb ) ) DEALLOCATE( desc%irb )
IF( associated( desc%imin3 ) ) DEALLOCATE( desc%imin3 )
IF( associated( desc%imax3 ) ) DEALLOCATE( desc%imax3 )
IF( associated( desc%npp ) ) DEALLOCATE( desc%npp )
IF( associated( desc%ipp ) ) DEALLOCATE( desc%ipp )
IF( associated( desc%np3 ) ) DEALLOCATE( desc%np3 )
desc%have_task_groups = .false.
END SUBROUTINE fft_box_deallocate
@ -179,22 +179,22 @@ CONTAINS
TYPE (fft_dlay_descriptor) :: desc
LOGICAL, INTENT(IN) :: tk
INTEGER, INTENT(IN) :: nst
INTEGER, INTENT(IN) :: nr1, nr2, nr3, nr1x, nr2x, nr3x
INTEGER, INTENT(IN) :: me ! processor index Starting from 1
INTEGER, INTENT(IN) :: nproc ! number of processor
INTEGER, INTENT(IN) :: nogrp ! number of groups for task-grouping
INTEGER, INTENT(IN) :: idx(:)
INTEGER, INTENT(IN) :: in1(:)
INTEGER, INTENT(IN) :: in2(:)
INTEGER, INTENT(IN) :: ncp(:)
INTEGER, INTENT(IN) :: ncpw(:)
INTEGER, INTENT(IN) :: ngp(:)
INTEGER, INTENT(IN) :: ngpw(:)
INTEGER, INTENT(IN) :: lb(:), ub(:)
INTEGER, INTENT(IN) :: st( lb(1) : ub(1), lb(2) : ub(2) )
INTEGER, INTENT(IN) :: stw( lb(1) : ub(1), lb(2) : ub(2) )
LOGICAL, INTENT(in) :: tk
INTEGER, INTENT(in) :: nst
INTEGER, INTENT(in) :: nr1, nr2, nr3, nr1x, nr2x, nr3x
INTEGER, INTENT(in) :: me ! processor index Starting from 1
INTEGER, INTENT(in) :: nproc ! number of processor
INTEGER, INTENT(in) :: nogrp ! number of groups for task-grouping
INTEGER, INTENT(in) :: idx(:)
INTEGER, INTENT(in) :: in1(:)
INTEGER, INTENT(in) :: in2(:)
INTEGER, INTENT(in) :: ncp(:)
INTEGER, INTENT(in) :: ncpw(:)
INTEGER, INTENT(in) :: ngp(:)
INTEGER, INTENT(in) :: ngpw(:)
INTEGER, INTENT(in) :: lb(:), ub(:)
INTEGER, INTENT(in) :: st( lb(1) : ub(1), lb(2) : ub(2) )
INTEGER, INTENT(in) :: stw( lb(1) : ub(1), lb(2) : ub(2) )
INTEGER :: npp( nproc ), n3( nproc ), nsp( nproc )
INTEGER :: np, nq, i, is, iss, i1, i2, m1, m2, n1, n2, ip
@ -203,20 +203,20 @@ CONTAINS
!
INTEGER :: sm
IF( ( SIZE( desc%ngl ) < nproc ) .OR. ( SIZE( desc%npp ) < nproc ) .OR. &
( SIZE( desc%ipp ) < nproc ) .OR. ( SIZE( desc%iss ) < nproc ) ) &
IF( ( size( desc%ngl ) < nproc ) .or. ( size( desc%npp ) < nproc ) .or. &
( size( desc%ipp ) < nproc ) .or. ( size( desc%iss ) < nproc ) ) &
CALL errore( ' fft_dlay_set ', ' wrong descriptor dimensions ', 1 )
IF( ( nr1 > nr1x ) .OR. ( nr2 > nr2x ) .OR. ( nr3 > nr3x ) ) &
IF( ( nr1 > nr1x ) .or. ( nr2 > nr2x ) .or. ( nr3 > nr3x ) ) &
CALL errore( ' fft_dlay_set ', ' wrong fft dimensions ', 2 )
IF( ( SIZE( idx ) < nst ) .OR. ( SIZE( in1 ) < nst ) .OR. ( SIZE( in2 ) < nst ) ) &
IF( ( size( idx ) < nst ) .or. ( size( in1 ) < nst ) .or. ( size( in2 ) < nst ) ) &
CALL errore( ' fft_dlay_set ', ' wrong number of stick dimensions ', 3 )
IF( ( SIZE( ncp ) < nproc ) .OR. ( SIZE( ngp ) < nproc ) ) &
IF( ( size( ncp ) < nproc ) .or. ( size( ngp ) < nproc ) ) &
CALL errore( ' fft_dlay_set ', ' wrong stick dimensions ', 4 )
desc%have_task_groups = .FALSE.
desc%have_task_groups = .false.
! Set the number of "xy" planes for each processor
! in other word do a slab partition along the z axis
@ -225,23 +225,23 @@ CONTAINS
npp = 0
IF ( nproc == 1 ) THEN
npp(1) = nr3
ELSE IF( nproc <= nr3 ) THEN
ELSEIF( nproc <= nr3 ) THEN
np = nr3 / nproc
nq = nr3 - np * nproc
DO i = 1, nproc
npp(i) = np
IF ( i <= nq ) npp(i) = np + 1
END DO
ENDDO
ELSE
DO ip = 1, nr3 ! some compiler complains for empty DO loops
DO i = 1, nproc, nogrp
npp(i) = npp(i) + 1
sm = sm + 1
IF ( sm == nr3 ) EXIT
END DO
IF ( sm == nr3 ) EXIT
END DO
END IF
IF ( sm == nr3 ) exit
ENDDO
IF ( sm == nr3 ) exit
ENDDO
ENDIF
desc%npp( 1:nproc ) = npp
desc%npl = npp( me )
@ -251,17 +251,17 @@ CONTAINS
n3 = 0
DO i = 2, nproc
n3(i) = n3(i-1) + npp(i-1)
END DO
ENDDO
desc%ipp( 1:nproc ) = n3
! Set the proper number of sticks
IF( .NOT. tk ) THEN
IF( .not. tk ) THEN
desc%nst = 2*nst - 1
ELSE
desc%nst = nst
END IF
ENDIF
! Set fft actual and leading dimensions
@ -277,33 +277,33 @@ CONTAINS
IF ( nproc == 1 ) THEN
desc%nnr = nr1x * nr2x * nr3x
desc%nnrx = desc%nnr
desc%nnrx = desc%nnr
ELSE
desc%nnr = MAX( nr3x * ncp(me), nr1x * nr2x * npp(me) )
desc%nnr = MAX( 1, desc%nnr ) ! ensure that desc%nrr > 0 ( for extreme parallelism )
desc%nnr = max( nr3x * ncp(me), nr1x * nr2x * npp(me) )
desc%nnr = max( 1, desc%nnr ) ! ensure that desc%nrr > 0 ( for extreme parallelism )
desc%nnrx = desc%nnr
DO i = 1, nproc
desc%nnrx = MAX( desc%nnrx, nr3x * ncp( i ) )
desc%nnrx = MAX( desc%nnrx, nr1x * nr2x * npp( i ) )
END DO
desc%nnrx = MAX( 1, desc%nnrx ) ! ensure that desc%nrr > 0 ( for extreme parallelism )
END IF
desc%nnrx = max( desc%nnrx, nr3x * ncp( i ) )
desc%nnrx = max( desc%nnrx, nr1x * nr2x * npp( i ) )
ENDDO
desc%nnrx = max( 1, desc%nnrx ) ! ensure that desc%nrr > 0 ( for extreme parallelism )
ENDIF
desc%ngl( 1:nproc ) = ngp( 1:nproc )
desc%nwl( 1:nproc ) = ngpw( 1:nproc )
IF( SIZE( desc%isind ) < ( nr1x * nr2x ) ) &
IF( size( desc%isind ) < ( nr1x * nr2x ) ) &
CALL errore( ' fft_dlay_set ', ' wrong descriptor dimensions, isind ', 5 )
IF( SIZE( desc%iplp ) < ( nr1x ) .OR. SIZE( desc%iplw ) < ( nr1x ) ) &
IF( size( desc%iplp ) < ( nr1x ) .or. size( desc%iplw ) < ( nr1x ) ) &
CALL errore( ' fft_dlay_set ', ' wrong descriptor dimensions, ipl ', 5 )
!
! 1. Temporarily store in the array "desc%isind" the index of the processor
! that own the corresponding stick (index of proc starting from 1)
! 2. Set the array elements of "desc%iplw" and "desc%iplp" to one
! 2. Set the array elements of "desc%iplw" and "desc%iplp" to one
! for that index corresponding to YZ planes containing at least one stick
! this are used in the FFT transform along Y
!
@ -317,28 +317,28 @@ CONTAINS
i1 = in1( is )
i2 = in2( is )
IF( st( i1, i2 ) > 0 ) THEN
m1 = i1 + 1; if ( m1 < 1 ) m1 = m1 + nr1
m2 = i2 + 1; if ( m2 < 1 ) m2 = m2 + nr2
m1 = i1 + 1; IF ( m1 < 1 ) m1 = m1 + nr1
m2 = i2 + 1; IF ( m2 < 1 ) m2 = m2 + nr2
IF( stw( i1, i2 ) > 0 ) THEN
desc%isind( m1 + ( m2 - 1 ) * nr1x ) = st( i1, i2 )
desc%iplw( m1 ) = 1
ELSE
ELSE
desc%isind( m1 + ( m2 - 1 ) * nr1x ) = -st( i1, i2 )
END IF
ENDIF
desc%iplp( m1 ) = 1
IF( .NOT. tk ) THEN
n1 = -i1 + 1; if ( n1 < 1 ) n1 = n1 + nr1
n2 = -i2 + 1; if ( n2 < 1 ) n2 = n2 + nr2
IF( .not. tk ) THEN
n1 = -i1 + 1; IF ( n1 < 1 ) n1 = n1 + nr1
n2 = -i2 + 1; IF ( n2 < 1 ) n2 = n2 + nr2
IF( stw( -i1, -i2 ) > 0 ) THEN
desc%isind( n1 + ( n2 - 1 ) * nr1x ) = st( -i1, -i2 )
desc%iplw( n1 ) = 1
ELSE
ELSE
desc%isind( n1 + ( n2 - 1 ) * nr1x ) = -st( -i1, -i2 )
END IF
ENDIF
desc%iplp( n1 ) = 1
END IF
END IF
END DO
ENDIF
ENDIF
ENDDO
!
! Compute for each proc the global index ( starting from 0 ) of the first
@ -350,16 +350,16 @@ CONTAINS
desc%iss( i ) = 0
ELSE
desc%iss( i ) = desc%iss( i - 1 ) + ncp( i - 1 )
END IF
END DO
ENDIF
ENDDO
IF( SIZE( desc%ismap ) < ( nst ) ) &
IF( size( desc%ismap ) < ( nst ) ) &
CALL errore( ' fft_dlay_set ', ' wrong descriptor dimensions ', 6 )
!
! 1. Set the array desc%ismap which maps stick indexes to
! 1. Set the array desc%ismap which maps stick indexes to
! position in the palne ( iss )
! 2. Re-set the array "desc%isind", that maps position
! 2. Re-set the array "desc%isind", that maps position
! in the plane with stick indexes (it is the inverse of desc%ismap )
!
@ -367,7 +367,7 @@ CONTAINS
desc%ismap = 0
nsp = 0
DO iss = 1, SIZE( desc%isind )
DO iss = 1, size( desc%isind )
ip = desc%isind( iss )
IF( ip > 0 ) THEN
nsp( ip ) = nsp( ip ) + 1
@ -376,24 +376,24 @@ CONTAINS
desc%isind( iss ) = nsp( ip )
ELSE
desc%isind( iss ) = 0
END IF
END IF
END DO
ENDIF
ENDIF
ENDDO
! chack number of stick against the input value
IF( ANY( nsp( 1:nproc ) /= ncpw( 1:nproc ) ) ) THEN
IF( any( nsp( 1:nproc ) /= ncpw( 1:nproc ) ) ) THEN
DO ip = 1, nproc
WRITE( stdout,*) ' * ', ip, ' * ', nsp( ip ), ' /= ', ncpw( ip )
END DO
ENDDO
CALL errore( ' fft_dlay_set ', ' inconsistent number of sticks ', 7 )
END IF
ENDIF
desc%nsw( 1:nproc ) = nsp
! then add pseudopotential stick
! then add pseudopotential stick
DO iss = 1, SIZE( desc%isind )
DO iss = 1, size( desc%isind )
ip = desc%isind( iss )
IF( ip < 0 ) THEN
nsp( -ip ) = nsp( -ip ) + 1
@ -402,18 +402,18 @@ CONTAINS
desc%isind( iss ) = nsp( -ip )
ELSE
desc%isind( iss ) = 0
END IF
END IF
END DO
ENDIF
ENDIF
ENDDO
! chack number of stick against the input value
IF( ANY( nsp( 1:nproc ) /= ncp( 1:nproc ) ) ) THEN
IF( any( nsp( 1:nproc ) /= ncp( 1:nproc ) ) ) THEN
DO ip = 1, nproc
WRITE( stdout,*) ' * ', ip, ' * ', nsp( ip ), ' /= ', ncp( ip )
END DO
ENDDO
CALL errore( ' fft_dlay_set ', ' inconsistent number of sticks ', 8 )
END IF
ENDIF
desc%nsp( 1:nproc ) = nsp
@ -421,7 +421,7 @@ CONTAINS
desc%id = icount
! Initialize the pointer to the fft tables
desc%tptr = icount
RETURN
@ -436,20 +436,20 @@ CONTAINS
TYPE (fft_dlay_descriptor) :: desc
INTEGER, INTENT(IN) :: nat, me, nproc
INTEGER, INTENT(IN) :: irb( :, : )
INTEGER, INTENT(IN) :: npp( : )
INTEGER, INTENT(IN) :: ipp( : )
INTEGER, INTENT(IN) :: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx
INTEGER, INTENT(in) :: nat, me, nproc
INTEGER, INTENT(in) :: irb( :, : )
INTEGER, INTENT(in) :: npp( : )
INTEGER, INTENT(in) :: ipp( : )
INTEGER, INTENT(in) :: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx
INTEGER :: ir3, ibig3, irb3, imin3, imax3, nr3, isa
IF( nat > SIZE( desc%irb, 2 ) ) THEN
WRITE( stdout, fmt="( ///,'NAT, SIZE = ',2I10)" ) nat, SIZE( desc%irb, 2 )
IF( nat > size( desc%irb, 2 ) ) THEN
WRITE( stdout, fmt="( ///,'NAT, SIZE = ',2I10)" ) nat, size( desc%irb, 2 )
CALL errore(" fft_box_set ", " inconsistent dimensions ", 1 )
END IF
ENDIF
IF( nproc > SIZE( desc%npp ) ) &
IF( nproc > size( desc%npp ) ) &
CALL errore(" fft_box_set ", " inconsistent dimensions ", 2 )
desc%nr1 = nr1b
@ -459,11 +459,11 @@ CONTAINS
desc%nr2x = nr2bx
desc%nr3x = nr3bx
desc%irb( 1:3, 1:nat ) = irb( 1:3, 1:nat )
desc%npp( 1:nproc ) = npp( 1:nproc )
desc%ipp( 1:nproc ) = ipp( 1:nproc )
desc%irb( 1:3, 1:nat ) = irb( 1:3, 1:nat )
desc%npp( 1:nproc ) = npp( 1:nproc )
desc%ipp( 1:nproc ) = ipp( 1:nproc )
nr3 = SUM( npp( 1:nproc ) )
nr3 = sum( npp( 1:nproc ) )
DO isa = 1, nat
@ -471,24 +471,24 @@ CONTAINS
imax3 = 1
irb3 = irb( 3, isa )
do ir3 = 1, nr3b
ibig3 = 1 + MOD( irb3 + ir3 - 2, nr3 )
if( ibig3 < 1 .or. ibig3 > nr3 ) &
& call errore(' fft_box_set ',' ibig3 wrong ', ibig3 )
DO ir3 = 1, nr3b
ibig3 = 1 + mod( irb3 + ir3 - 2, nr3 )
IF( ibig3 < 1 .or. ibig3 > nr3 ) &
& CALL errore(' fft_box_set ',' ibig3 wrong ', ibig3 )
ibig3 = ibig3 - ipp( me )
if ( ibig3 > 0 .and. ibig3 <= npp(me) ) then
IF ( ibig3 > 0 .and. ibig3 <= npp(me) ) THEN
imin3 = min( imin3, ir3 )
imax3 = max( imax3, ir3 )
end if
end do
ENDIF
ENDDO
desc%imin3( isa ) = imin3
desc%imax3( isa ) = imax3
desc%np3( isa ) = imax3 - imin3 + 1
END DO
ENDDO
desc%have_task_groups = .FALSE.
desc%have_task_groups = .false.
END SUBROUTINE fft_box_set
@ -497,16 +497,16 @@ CONTAINS
SUBROUTINE fft_dlay_scalar( desc, ub, lb, nr1, nr2, nr3, nr1x, nr2x, nr3x, stw )
implicit none
IMPLICIT NONE
TYPE (fft_dlay_descriptor) :: desc
INTEGER, INTENT(IN) :: lb(:), ub(:)
INTEGER, INTENT(IN) :: stw( lb(2) : ub(2), lb(3) : ub(3) )
INTEGER, INTENT(in) :: lb(:), ub(:)
INTEGER, INTENT(in) :: stw( lb(2) : ub(2), lb(3) : ub(3) )
integer :: nr1, nr2, nr3, nr1x, nr2x, nr3x
integer :: m1, m2, i2, i3
INTEGER :: nr1, nr2, nr3, nr1x, nr2x, nr3x
INTEGER :: m1, m2, i2, i3
IF( SIZE( desc%iplw ) < nr3x .OR. SIZE( desc%isind ) < nr2x * nr3x ) &
IF( size( desc%iplw ) < nr3x .or. size( desc%isind ) < nr2x * nr3x ) &
CALL errore(' fft_dlay_scalar ', ' wrong dimensions ', 1 )
desc%isind = 0
@ -525,21 +525,21 @@ CONTAINS
DO i2 = lb( 2 ), ub( 2 )
DO i3 = lb( 3 ), ub( 3 )
m1 = i2 + 1; if ( m1 < 1 ) m1 = m1 + nr2
m2 = i3 + 1; if ( m2 < 1 ) m2 = m2 + nr3
m1 = i2 + 1; IF ( m1 < 1 ) m1 = m1 + nr2
m2 = i3 + 1; IF ( m2 < 1 ) m2 = m2 + nr3
IF( stw( i2, i3 ) > 0 ) THEN
desc%isind( m1 + ( m2 - 1 ) * nr2x ) = 1 ! st( i1, i2 )
desc%iplw( m2 ) = 1
END IF
END DO
END DO
ENDIF
ENDDO
ENDDO
desc%nnr = nr1x * nr2x * nr3x
desc%npl = nr3
desc%nnp = nr1x * nr2x
desc%npp = nr3
desc%ipp = 0
desc%have_task_groups = .FALSE.
desc%have_task_groups = .false.
desc%nnrx = desc%nnr
RETURN

450
Modules/stick_base.f90
View File

@ -19,11 +19,11 @@
PUBLIC :: sticks_maps, sticks_sort, sticks_countg, sticks_dist, sticks_pairup
PUBLIC :: sticks_owner, sticks_deallocate, pstickset, sticks_maps_scalar
! ... sticks_owner : stick owner, sticks_owner( i, j ) is the index of the processor
! ... (starting from 1) owning the stick whose x and y coordinate are i and j.
INTEGER, ALLOCATABLE, TARGET :: sticks_owner( : , : )
INTEGER, ALLOCATABLE, TARGET :: sticks_owner( : , : )
INTERFACE sticks_dist
MODULE PROCEDURE sticks_dist1
@ -38,16 +38,16 @@
USE mp, ONLY: mp_sum
USE mp_global, ONLY: me_pool, nproc_pool, intra_pool_comm
LOGICAL, INTENT(IN) :: tk ! if true use the full space grid
INTEGER, INTENT(IN) :: ub(:) ! upper bounds for i-th grid dimension
INTEGER, INTENT(IN) :: lb(:) ! lower bounds for i-th grid dimension
REAL(DP) , INTENT(IN) :: b1(:), b2(:), b3(:) ! reciprocal space base vectors
REAL(DP) , INTENT(IN) :: gcut ! cut-off for potentials
REAL(DP) , INTENT(IN) :: gcutw ! cut-off for plane waves
REAL(DP) , INTENT(IN) :: gcuts ! cut-off for smooth mesh
INTEGER, INTENT(OUT) :: st( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(OUT) :: stw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(OUT) :: sts(lb(1): ub(1), lb(2):ub(2) ) ! stick map for smooth mesh
LOGICAL, INTENT(in) :: tk ! if true use the full space grid
INTEGER, INTENT(in) :: ub(:) ! upper bounds for i-th grid dimension
INTEGER, INTENT(in) :: lb(:) ! lower bounds for i-th grid dimension
REAL(DP) , INTENT(in) :: b1(:), b2(:), b3(:) ! reciprocal space base vectors
REAL(DP) , INTENT(in) :: gcut ! cut-off for potentials
REAL(DP) , INTENT(in) :: gcutw ! cut-off for plane waves
REAL(DP) , INTENT(in) :: gcuts ! cut-off for smooth mesh
INTEGER, INTENT(out) :: st( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(out) :: stw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(out) :: sts(lb(1): ub(1), lb(2):ub(2) ) ! stick map for smooth mesh
INTEGER :: i, j, k, kip
REAL(DP) :: gsq
@ -60,111 +60,111 @@
! ... cut-off gcut, wavefunction cut-off gcutw, and smooth mesh cut-off gcuts
! ... st(i,j) will contain the number of G vectors of the stick whose
! ... indices are (i,j).
#if defined (__EKO)
write(*,*) ! Workaround for EKOPath compiler bug
#endif
IF( .NOT. tk ) THEN
! ... indices are (i,j).
kip = 0 + ABS(lb(3)) + 1
IF( MOD( kip, nproc_pool ) == me_pool ) THEN
#if defined (__EKO)
WRITE(*,*) ! Workaround for EKOPath compiler bug
#endif
IF( .not. tk ) THEN
kip = 0 + abs(lb(3)) + 1
IF( mod( kip, nproc_pool ) == me_pool ) THEN
st (0,0) = st (0,0) + 1
stw(0,0) = stw(0,0) + 1
sts(0,0) = sts(0,0) + 1
END IF
ENDIF
DO i= 0, 0
DO j= 0, 0
DO i= 0, 0
DO j= 0, 0
DO k= 1, ub(3)
kip = k + ABS(lb(3)) + 1
IF( MOD( kip, nproc_pool ) == me_pool ) THEN
gsq= (DBLE(i)*b1(1)+DBLE(j)*b2(1)+DBLE(k)*b3(1) )**2
gsq=gsq+(DBLE(i)*b1(2)+DBLE(j)*b2(2)+DBLE(k)*b3(2) )**2
gsq=gsq+(DBLE(i)*b1(3)+DBLE(j)*b2(3)+DBLE(k)*b3(3) )**2
IF(gsq.LE.gcut ) THEN
kip = k + abs(lb(3)) + 1
IF( mod( kip, nproc_pool ) == me_pool ) THEN
gsq= (dble(i)*b1(1)+dble(j)*b2(1)+dble(k)*b3(1) )**2
gsq=gsq+(dble(i)*b1(2)+dble(j)*b2(2)+dble(k)*b3(2) )**2
gsq=gsq+(dble(i)*b1(3)+dble(j)*b2(3)+dble(k)*b3(3) )**2
IF(gsq.le.gcut ) THEN
st(i,j) = st(i,j) + 1
IF(gsq.LE.gcutw) THEN
IF(gsq.le.gcutw) THEN
stw(i,j) = stw(i,j) + 1
END IF
IF(gsq.LE.gcuts) THEN
ENDIF
IF(gsq.le.gcuts) THEN
sts(i,j) = sts(i,j) + 1
END IF
END IF
END IF
END DO
END DO
END DO
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
DO i = 0, 0
DO j = 1, ub(2)
DO k = lb(3), ub(3)
kip = k + ABS(lb(3)) + 1
IF( MOD( kip, nproc_pool) == me_pool ) THEN
gsq= (DBLE(i)*b1(1)+DBLE(j)*b2(1)+DBLE(k)*b3(1) )**2
gsq=gsq+(DBLE(i)*b1(2)+DBLE(j)*b2(2)+DBLE(k)*b3(2) )**2
gsq=gsq+(DBLE(i)*b1(3)+DBLE(j)*b2(3)+DBLE(k)*b3(3) )**2
IF(gsq.LE.gcut ) THEN
kip = k + abs(lb(3)) + 1
IF( mod( kip, nproc_pool) == me_pool ) THEN
gsq= (dble(i)*b1(1)+dble(j)*b2(1)+dble(k)*b3(1) )**2
gsq=gsq+(dble(i)*b1(2)+dble(j)*b2(2)+dble(k)*b3(2) )**2
gsq=gsq+(dble(i)*b1(3)+dble(j)*b2(3)+dble(k)*b3(3) )**2
IF(gsq.le.gcut ) THEN
st(i,j) = st(i,j) + 1
IF(gsq.LE.gcutw) THEN
IF(gsq.le.gcutw) THEN
stw(i,j) = stw(i,j) + 1
END IF
IF(gsq.LE.gcuts) THEN
ENDIF
IF(gsq.le.gcuts) THEN
sts(i,j) = sts(i,j) + 1
END IF
END IF
END IF
END DO
END DO
END DO
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
DO i = 1, ub(1)
DO j = lb(2), ub(2)
DO k = lb(3), ub(3)
kip = k + ABS(lb(3)) + 1
IF( MOD( kip, nproc_pool) == me_pool ) THEN
gsq= (DBLE(i)*b1(1)+DBLE(j)*b2(1)+DBLE(k)*b3(1) )**2
gsq=gsq+(DBLE(i)*b1(2)+DBLE(j)*b2(2)+DBLE(k)*b3(2) )**2
gsq=gsq+(DBLE(i)*b1(3)+DBLE(j)*b2(3)+DBLE(k)*b3(3) )**2
IF(gsq.LE.gcut ) THEN
kip = k + abs(lb(3)) + 1
IF( mod( kip, nproc_pool) == me_pool ) THEN
gsq= (dble(i)*b1(1)+dble(j)*b2(1)+dble(k)*b3(1) )**2
gsq=gsq+(dble(i)*b1(2)+dble(j)*b2(2)+dble(k)*b3(2) )**2
gsq=gsq+(dble(i)*b1(3)+dble(j)*b2(3)+dble(k)*b3(3) )**2
IF(gsq.le.gcut ) THEN
st(i,j) = st(i,j) + 1
IF(gsq.LE.gcutw) THEN
IF(gsq.le.gcutw) THEN
stw(i,j) = stw(i,j) + 1
END IF
IF(gsq.LE.gcuts) THEN
ENDIF
IF(gsq.le.gcuts) THEN
sts(i,j) = sts(i,j) + 1
END IF
END IF
END IF
END DO
END DO
END DO
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
ELSE
DO i= lb(1), ub(1)
DO j= lb(2), ub(2)
DO k= lb(3), ub(3)
kip = k + ABS(lb(3)) + 1
IF( MOD( kip, nproc_pool ) == me_pool ) THEN
gsq= (DBLE(i)*b1(1)+DBLE(j)*b2(1)+DBLE(k)*b3(1) )**2
gsq=gsq+(DBLE(i)*b1(2)+DBLE(j)*b2(2)+DBLE(k)*b3(2) )**2
gsq=gsq+(DBLE(i)*b1(3)+DBLE(j)*b2(3)+DBLE(k)*b3(3) )**2
IF(gsq.LE.gcut ) THEN
kip = k + abs(lb(3)) + 1
IF( mod( kip, nproc_pool ) == me_pool ) THEN
gsq= (dble(i)*b1(1)+dble(j)*b2(1)+dble(k)*b3(1) )**2
gsq=gsq+(dble(i)*b1(2)+dble(j)*b2(2)+dble(k)*b3(2) )**2
gsq=gsq+(dble(i)*b1(3)+dble(j)*b2(3)+dble(k)*b3(3) )**2
IF(gsq.le.gcut ) THEN
st(i,j) = st(i,j) + 1
END IF
IF(gsq.LE.gcutw) THEN
ENDIF
IF(gsq.le.gcutw) THEN
stw(i,j) = stw(i,j) + 1
END IF
IF(gsq.LE.gcuts) THEN
ENDIF
IF(gsq.le.gcuts) THEN
sts(i,j) = sts(i,j) + 1
END IF
END IF
END DO
END DO
END DO
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
END IF
ENDIF
CALL mp_sum(st ,intra_pool_comm )
CALL mp_sum(stw ,intra_pool_comm )
@ -175,7 +175,7 @@
WRITE( 6,*) 'testtesttesttesttesttesttesttesttesttest'
WRITE( 6,*) 'lb = ', lb(1), lb(2)
WRITE( 6,*) 'ub = ', ub(1), ub(2)
WRITE( 6,*) 'counts = ', COUNT( st > 0 ), COUNT( stw > 0 ), COUNT( sts > 0 )
WRITE( 6,*) 'counts = ', count( st > 0 ), count( stw > 0 ), count( sts > 0 )
WRITE( 6,*) 'cut-offs = ', gcut, gcutw, gcuts
WRITE( 6,*) 'b1 = ', b1(1:3)
WRITE( 6,*) 'b2 = ', b2(1:3)
@ -183,8 +183,8 @@
DO i = lb(1), ub(1)
DO j = lb(2), ub(2)
WRITE( 6,'(2I4,3I6)') i,j,st(i,j),stw(i,j),sts(i,j)
END DO
END DO
ENDDO
ENDDO
WRITE( 6,*) 'testtesttesttesttesttesttesttesttesttest'
! Test sticks
#endif
@ -196,19 +196,19 @@
SUBROUTINE sticks_maps_scalar( lgamma, ub, lb, b1, b2, b3, gcutm, gkcut, gcutms, stw, ngm, ngms )
LOGICAL, INTENT(IN) :: lgamma ! if true use gamma point simmetry
INTEGER, INTENT(IN) :: ub(:) ! upper bounds for i-th grid dimension
INTEGER, INTENT(IN) :: lb(:) ! lower bounds for i-th grid dimension
REAL(DP) , INTENT(IN) :: b1(:), b2(:), b3(:) ! reciprocal space base vectors
REAL(DP) , INTENT(IN) :: gcutm ! cut-off for potentials
REAL(DP) , INTENT(IN) :: gkcut ! cut-off for plane waves
REAL(DP) , INTENT(IN) :: gcutms ! cut-off for smooth mesh
LOGICAL, INTENT(in) :: lgamma ! if true use gamma point simmetry
INTEGER, INTENT(in) :: ub(:) ! upper bounds for i-th grid dimension
INTEGER, INTENT(in) :: lb(:) ! lower bounds for i-th grid dimension
REAL(DP) , INTENT(in) :: b1(:), b2(:), b3(:) ! reciprocal space base vectors
REAL(DP) , INTENT(in) :: gcutm ! cut-off for potentials
REAL(DP) , INTENT(in) :: gkcut ! cut-off for plane waves
REAL(DP) , INTENT(in) :: gcutms ! cut-off for smooth mesh
!
INTEGER, INTENT(OUT) :: ngm, ngms
INTEGER, INTENT(out) :: ngm, ngms
!
! stick map for wave functions, note that map is taken in YZ plane
!
INTEGER, INTENT(OUT) :: stw( lb(2) : ub(2), lb(3) : ub(3) )
INTEGER, INTENT(out) :: stw( lb(2) : ub(2), lb(3) : ub(3) )
INTEGER :: i1, i2, i3, n1, n2, n3
REAL(DP) :: amod
@ -216,40 +216,40 @@
ngm = 0
ngms = 0
n1 = MAX( ABS( lb(1) ), ABS( ub(1) ) )
n2 = MAX( ABS( lb(2) ), ABS( ub(2) ) )
n3 = MAX( ABS( lb(3) ), ABS( ub(3) ) )
n1 = max( abs( lb(1) ), abs( ub(1) ) )
n2 = max( abs( lb(2) ), abs( ub(2) ) )
n3 = max( abs( lb(3) ), abs( ub(3) ) )
loop1: do i1 = - n1, n1
loop1: DO i1 = - n1, n1
!
! Gamma-only: exclude space with x<0
!
if (lgamma .and. i1 < 0) cycle loop1
IF (lgamma .and. i1 < 0) CYCLE loop1
!
loop2: do i2 = - n2, n2
loop2: DO i2 = - n2, n2
!
! Gamma-only: exclude plane with x=0, y<0
!
if(lgamma .and. i1 == 0.and. i2 < 0) cycle loop2
IF(lgamma .and. i1 == 0.and. i2 < 0) CYCLE loop2
!
loop3: do i3 = - n3, n3
loop3: DO i3 = - n3, n3
!
! Gamma-only: exclude line with x=0, y=0, z<0
!
if(lgamma .and. i1 == 0 .and. i2 == 0 .and. i3 < 0) cycle loop3
IF(lgamma .and. i1 == 0 .and. i2 == 0 .and. i3 < 0) CYCLE loop3
!
amod = (i1 * b1 (1) + i2 * b2 (1) + i3 * b3 (1) ) **2 + &
(i1 * b1 (2) + i2 * b2 (2) + i3 * b3 (2) ) **2 + &
(i1 * b1 (3) + i2 * b2 (3) + i3 * b3 (3) ) **2
if (amod <= gcutm) ngm = ngm + 1
if (amod <= gcutms) ngms = ngms + 1
if (amod <= gkcut ) then
IF (amod <= gcutm) ngm = ngm + 1
IF (amod <= gcutms) ngms = ngms + 1
IF (amod <= gkcut ) THEN
stw( i2, i3 ) = 1
if (lgamma) stw( -i2, -i3 ) = 1
end if
enddo loop3
enddo loop2
enddo loop1
IF (lgamma) stw( -i2, -i3 ) = 1
ENDIF
ENDDO loop3
ENDDO loop2
ENDDO loop1
RETURN
END SUBROUTINE sticks_maps_scalar
@ -259,7 +259,7 @@
SUBROUTINE sticks_sort( ngc, ngcw, ngcs, nct, idx )
! ... This subroutine sorts the sticks indexes, according to
! ... This subroutine sorts the sticks indexes, according to
! ... the length and type of the sticks, wave functions sticks
! ... first, then smooth mesh sticks, and finally potential
! ... sticks
@ -269,21 +269,21 @@
! lengths of sticks, ngc for potential mesh, ngcw for wave functions mesh
! and ngcs for smooth mesh
INTEGER, INTENT(IN) :: ngc(:), ngcw(:), ngcs(:)
INTEGER, INTENT(in) :: ngc(:), ngcw(:), ngcs(:)
! nct, total number of sticks
INTEGER, INTENT(IN) :: nct
INTEGER, INTENT(in) :: nct
! index, on output, new sticks indexes
INTEGER, INTENT(OUT) :: idx(:)
INTEGER, INTENT(out) :: idx(:)
INTEGER :: mc, nr3x, ic
REAL(DP) :: dn3
REAL(DP), ALLOCATABLE :: aux(:)
nr3x = MAXVAL( ngc(1:nct) ) + 1
nr3x = maxval( ngc(1:nct) ) + 1
dn3 = REAL( nr3x )
IF( nproc_pool > 1 ) THEN
@ -294,37 +294,37 @@
aux(mc) = dn3 * aux(mc) + ngc(mc)
aux(mc) = -aux(mc)
idx(mc) = 0
END DO
ENDDO
CALL hpsort( nct, aux(1), idx(1))
DEALLOCATE( aux )
ELSE
ic = 0
do mc = 1, nct
if( ngcw(mc) > 0 ) then
DO mc = 1, nct
IF( ngcw(mc) > 0 ) THEN
ic = ic + 1
idx(ic) = mc
endif
end do
do mc = 1, nct
if( ngcs(mc) > 0 .AND. ngcw(mc) == 0 ) then
ENDIF
ENDDO
DO mc = 1, nct
IF( ngcs(mc) > 0 .and. ngcw(mc) == 0 ) THEN
ic = ic + 1
idx(ic) = mc
endif
end do
do mc = 1, nct
if( ngc(mc) > 0 .AND. ngcs(mc) == 0 .AND. ngcw(mc) == 0 ) then
ENDIF
ENDDO
DO mc = 1, nct
IF( ngc(mc) > 0 .and. ngcs(mc) == 0 .and. ngcw(mc) == 0 ) THEN
ic = ic + 1
idx(ic) = mc
endif
end do
END IF
ENDIF
ENDDO
ENDIF
#if defined __STICKS_DEBUG
WRITE( 6,*) '-----------------'
WRITE( 6,*) 'STICKS_SORT DEBUG'
DO mc = 1, nct
WRITE( 6, fmt="(4I10)" ) idx(mc), ngcw( idx(mc) ), ngcs( idx(mc) ), ngc( idx(mc) )
END DO
ENDDO
WRITE( 6,*) '-----------------'
#endif
@ -335,14 +335,14 @@
SUBROUTINE sticks_countg( tk, ub, lb, st, stw, sts, in1, in2, ngc, ngcw, ngcs )
INTEGER, INTENT(IN) :: ub(:), lb(:)
INTEGER, INTENT(IN) :: st( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(IN) :: stw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(IN) :: sts(lb(1): ub(1), lb(2):ub(2) ) ! stick map for smooth mesh
LOGICAL, INTENT(IN) :: tk
INTEGER, INTENT(in) :: ub(:), lb(:)
INTEGER, INTENT(in) :: st( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(in) :: stw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(in) :: sts(lb(1): ub(1), lb(2):ub(2) ) ! stick map for smooth mesh
LOGICAL, INTENT(in) :: tk
INTEGER, INTENT(OUT) :: in1(:), in2(:)
INTEGER, INTENT(OUT) :: ngc(:), ngcw(:), ngcs(:)
INTEGER, INTENT(out) :: in1(:), in2(:)
INTEGER, INTENT(out) :: ngc(:), ngcw(:), ngcs(:)
INTEGER :: j1, j2, i1, i2, nct, min_size
@ -357,16 +357,16 @@
ngcs = 0
ngcw = 0
min_size = MIN( SIZE( in1 ), SIZE( in2 ), SIZE( ngc ), SIZE( ngcw ), SIZE( ngcs ) )
min_size = min( size( in1 ), size( in2 ), size( ngc ), size( ngcw ), size( ngcs ) )
DO j2 = 0, ( ub(2) - lb(2) )
DO j1 = 0, ( ub(1) - lb(1) )
i1 = j1
if( i1 > ub(1) ) i1 = lb(1) + ( i1 - ub(1) ) - 1
IF( i1 > ub(1) ) i1 = lb(1) + ( i1 - ub(1) ) - 1
i2 = j2
if( i2 > ub(2) ) i2 = lb(2) + ( i2 - ub(2) ) - 1
IF( i2 > ub(2) ) i2 = lb(2) + ( i2 - ub(2) ) - 1
IF( st( i1, i2 ) > 0 ) THEN
@ -380,15 +380,15 @@
in2(nct) = i2
ngc(nct) = st( i1 , i2)
IF( stw( i1, i2 ) .GT. 0 ) ngcw(nct) = stw( i1 , i2)
IF( sts( i1, i2 ) .GT. 0 ) ngcs(nct) = sts( i1 , i2)
IF( stw( i1, i2 ) .gt. 0 ) ngcw(nct) = stw( i1 , i2)
IF( sts( i1, i2 ) .gt. 0 ) ngcs(nct) = sts( i1 , i2)
END IF
ENDIF
! WRITE(7,fmt="(5I5)") i1, i2, nct, ngc(nct), ngcw( nct )
END DO
END DO
ENDDO
ENDDO
RETURN
END SUBROUTINE sticks_countg
@ -400,18 +400,18 @@
USE mp_global, ONLY: nproc_pool
LOGICAL, INTENT(IN) :: tk
LOGICAL, INTENT(in) :: tk
INTEGER, INTENT(IN) :: ub(:), lb(:), idx(:)
INTEGER, INTENT(OUT) :: stown( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(OUT) :: stownw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(OUT) :: stowns(lb(1): ub(1), lb(2):ub(2) ) ! stick map for smooth mesh
INTEGER, INTENT(in) :: ub(:), lb(:), idx(:)
INTEGER, INTENT(out) :: stown( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(out) :: stownw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(out) :: stowns(lb(1): ub(1), lb(2):ub(2) ) ! stick map for smooth mesh
INTEGER, INTENT(IN) :: in1(:), in2(:)
INTEGER, INTENT(IN) :: ngc(:), ngcw(:), ngcs(:)
INTEGER, INTENT(IN) :: nct
INTEGER, INTENT(OUT) :: ncp(:), ncpw(:), ncps(:)
INTEGER, INTENT(OUT) :: ngp(:), ngpw(:), ngps(:)
INTEGER, INTENT(in) :: in1(:), in2(:)
INTEGER, INTENT(in) :: ngc(:), ngcw(:), ngcs(:)
INTEGER, INTENT(in) :: nct
INTEGER, INTENT(out) :: ncp(:), ncpw(:), ncps(:)
INTEGER, INTENT(out) :: ngp(:), ngpw(:), ngps(:)
INTEGER :: mc, i1, i2, i, j, jj
@ -435,33 +435,33 @@
i1 = in1( i )
i2 = in2( i )
!
if ( ( .NOT. tk ) .AND. ( (i1 < 0) .or. ( (i1 == 0) .and. (i2 < 0) ) ) ) go to 30
IF ( ( .not. tk ) .and. ( (i1 < 0) .or. ( (i1 == 0) .and. (i2 < 0) ) ) ) GOTO 30
!
jj = 1
if ( ngcw(i) > 0 ) then
IF ( ngcw(i) > 0 ) THEN
!
! this is an active sticks: find which processor has currently
! the smallest number of plane waves
!
do j = 1, nproc_pool
if ( ngpw(j) < ngpw(jj) ) then
DO j = 1, nproc_pool
IF ( ngpw(j) < ngpw(jj) ) THEN
jj = j
else if ( ( ngpw(j) == ngpw(jj) ) .AND. ( ncpw(j) < ncpw(jj) ) ) then
ELSEIF ( ( ngpw(j) == ngpw(jj) ) .and. ( ncpw(j) < ncpw(jj) ) ) THEN
jj = j
end if
end do
ENDIF
ENDDO
else
ELSE
!
! this is an inactive sticks: find which processor has currently
! the smallest number of G-vectors
!
do j = 1, nproc_pool
if ( ngp(j) < ngp(jj) ) jj = j
end do
DO j = 1, nproc_pool
IF ( ngp(j) < ngp(jj) ) jj = j
ENDDO
end if
ENDIF
!
! potential mesh
@ -471,50 +471,50 @@
! smooth mesh
if ( ngcs(i) > 0 ) then
IF ( ngcs(i) > 0 ) THEN
ncps(jj) = ncps(jj) + 1
ngps(jj) = ngps(jj) + ngcs(i)
stowns(i1,i2) = jj
endif
ENDIF
! wave functions mesh
if ( ngcw(i) > 0 ) then
IF ( ngcw(i) > 0 ) THEN
ncpw(jj) = ncpw(jj) + 1
ngpw(jj) = ngpw(jj) + ngcw(i)
stownw(i1,i2) = jj
endif
ENDIF
30 continue
30 CONTINUE
END DO
ENDDO
RETURN
END SUBROUTINE sticks_dist1
!=----------------------------------------------------------------------=
SUBROUTINE sticks_pairup( tk, ub, lb, idx, in1, in2, ngc, ngcw, ngcs, nct, &
ncp, ncpw, ncps, ngp, ngpw, ngps, stown, stownw, stowns )
USE mp_global, ONLY: nproc_pool
LOGICAL, INTENT(IN) :: tk
LOGICAL, INTENT(in) :: tk
INTEGER, INTENT(IN) :: ub(:), lb(:), idx(:)
INTEGER, INTENT(INOUT) :: stown( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(INOUT) :: stownw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(INOUT) :: stowns(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(in) :: ub(:), lb(:), idx(:)
INTEGER, INTENT(inout) :: stown( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential
INTEGER, INTENT(inout) :: stownw(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(inout) :: stowns(lb(1): ub(1), lb(2):ub(2) ) ! stick map for wave functions
INTEGER, INTENT(IN) :: in1(:), in2(:)
INTEGER, INTENT(IN) :: ngc(:), ngcw(:), ngcs(:)
INTEGER, INTENT(IN) :: nct
INTEGER, INTENT(OUT) :: ncp(:), ncpw(:), ncps(:)
INTEGER, INTENT(OUT) :: ngp(:), ngpw(:), ngps(:)
INTEGER, INTENT(in) :: in1(:), in2(:)
INTEGER, INTENT(in) :: ngc(:), ngcw(:), ngcs(:)
INTEGER, INTENT(in) :: nct
INTEGER, INTENT(out) :: ncp(:), ncpw(:), ncps(:)
INTEGER, INTENT(out) :: ngp(:), ngpw(:), ngps(:)
INTEGER :: mc, i1, i2, i, jj
IF ( .NOT. tk ) THEN
IF ( .not. tk ) THEN
! when gamma symmetry is used only the sticks of half reciprocal space
! are generated, then here we pair-up the sticks with those of the other
@ -527,39 +527,39 @@
i2 = in2(i)
IF( i1 == 0 .and. i2 == 0 ) THEN
jj = stown( i1, i2 )
if( jj > 0 ) ngp( jj ) = ngp( jj ) + ngc( i ) - 1
IF( jj > 0 ) ngp( jj ) = ngp( jj ) + ngc( i ) - 1
jj = stowns( i1, i2 )
if( jj > 0 ) ngps( jj ) = ngps( jj ) + ngcs( i ) - 1
IF( jj > 0 ) ngps( jj ) = ngps( jj ) + ngcs( i ) - 1
jj = stownw( i1, i2 )
if( jj > 0 ) ngpw( jj ) = ngpw( jj ) + ngcw( i ) - 1
IF( jj > 0 ) ngpw( jj ) = ngpw( jj ) + ngcw( i ) - 1
ELSE
jj = stown( i1, i2 )
if( jj > 0 ) then
IF( jj > 0 ) THEN
stown( -i1, -i2 ) = jj
ncp( jj ) = ncp( jj ) + 1
ngp( jj ) = ngp( jj ) + ngc( i )
end if
ENDIF
jj = stowns( i1, i2 )
if( jj > 0 ) then
IF( jj > 0 ) THEN
stowns( -i1, -i2 ) = jj
ncps( jj ) = ncps( jj ) + 1
ngps( jj ) = ngps( jj ) + ngcs( i )
end if
ENDIF
jj = stownw( i1, i2 )
if( jj > 0 ) then
IF( jj > 0 ) THEN
stownw( -i1, -i2 ) = jj
ncpw( jj ) = ncpw( jj ) + 1
ngpw( jj ) = ngpw( jj ) + ngcw( i )
end if
END IF
END DO
ENDIF
ENDIF
ENDDO
END IF
ENDIF
IF( ALLOCATED( sticks_owner ) ) DEALLOCATE( sticks_owner )
IF( allocated( sticks_owner ) ) DEALLOCATE( sticks_owner )
ALLOCATE( sticks_owner( lb(1): ub(1), lb(2):ub(2) ) )
sticks_owner( :, : ) = ABS( stown( :, :) )
sticks_owner( :, : ) = abs( stown( :, :) )
RETURN
END SUBROUTINE sticks_pairup
@ -580,12 +580,12 @@
fft_dlay_scalar
TYPE(fft_dlay_descriptor), INTENT(INOUT) :: dfftp, dffts
REAL(DP), INTENT(IN) :: a1(3), a2(3), a3(3), alat
REAL(DP), INTENT(IN) :: gcut, gkcut, gcuts
INTEGER, INTENT(IN) :: nr1, nr2, nr3, nr1x, nr2x, nr3x
INTEGER, INTENT(IN) :: nr1s, nr2s, nr3s, nr1sx, nr2sx, nr3sx
INTEGER, INTENT(OUT) :: ngw, ngm, ngs
TYPE(fft_dlay_descriptor), INTENT(inout) :: dfftp, dffts
REAL(DP), INTENT(in) :: a1(3), a2(3), a3(3), alat
REAL(DP), INTENT(in) :: gcut, gkcut, gcuts
INTEGER, INTENT(in) :: nr1, nr2, nr3, nr1x, nr2x, nr3x
INTEGER, INTENT(in) :: nr1s, nr2s, nr3s, nr1sx, nr2sx, nr3sx
INTEGER, INTENT(out) :: ngw, ngm, ngs
LOGICAL :: tk
! ... tk logical flag, TRUE if the symulation does not have the
@ -667,7 +667,7 @@
INTEGER :: ip, ngm_ , ngs_
INTEGER, ALLOCATABLE :: idx(:)
tk = .NOT. gamma_only
tk = .not. gamma_only
ub(1) = ( nr1 - 1 ) / 2
ub(2) = ( nr2 - 1 ) / 2
ub(3) = ( nr3 - 1 ) / 2
@ -696,9 +696,9 @@
! ... Now count the number of stick nst and nstw
nst = COUNT( st > 0 )
nstw = COUNT( stw > 0 )
nsts = COUNT( sts > 0 )
nst = count( st > 0 )
nstw = count( stw > 0 )
nsts = count( sts > 0 )
IF (ionode) THEN
WRITE( stdout,*)
@ -707,7 +707,7 @@
3X,'----------')
WRITE( stdout,15) nst, nstw, nsts
15 FORMAT( 3X, 'nst =', I6, ', nstw =', I6, ', nsts =', I6 )
END IF
ENDIF
ALLOCATE(ist(nst,5))
@ -769,11 +769,11 @@
CALL sticks_maps_scalar( (.not.tk), ub, lb, b1, b2, b3, gcut, gkcut, gcuts, stw, ngm_ , ngs_ )
IF( ngm_ /= ngm ) CALL errore( ' pstickset ', ' inconsistent ngm ', ABS( ngm - ngm_ ) )
IF( ngs_ /= ngs ) CALL errore( ' pstickset ', ' inconsistent ngs ', ABS( ngs - ngs_ ) )
IF( ngm_ /= ngm ) CALL errore( ' pstickset ', ' inconsistent ngm ', abs( ngm - ngm_ ) )
IF( ngs_ /= ngs ) CALL errore( ' pstickset ', ' inconsistent ngs ', abs( ngs - ngs_ ) )
CALL fft_dlay_allocate( dfftp, nproc_pool, MAX(nr1x, nr3x), nr2x )
CALL fft_dlay_allocate( dffts, nproc_pool, MAX(nr1sx, nr3sx), nr2sx )
CALL fft_dlay_allocate( dfftp, nproc_pool, max(nr1x, nr3x), nr2x )
CALL fft_dlay_allocate( dffts, nproc_pool, max(nr1sx, nr3sx), nr2sx )
CALL fft_dlay_scalar( dfftp, ub, lb, nr1, nr2, nr3, nr1x, nr2x, nr3x, stw )
CALL fft_dlay_scalar( dffts, ub, lb, nr1s, nr2s, nr3s, nr1sx, nr2sx, nr3sx, stw )
@ -781,23 +781,23 @@
#endif
! ... Maximum number of sticks (potentials)
nstpx = MAXVAL( nstp )
nstpx = maxval( nstp )
! ... Maximum number of sticks (wave func.)
nstpwx = MAXVAL( nstpw )
nstpwx = maxval( nstpw )
IF (ionode) WRITE( stdout,118)
118 FORMAT(3X,' n.st n.stw n.sts n.g n.gw n.gs')
119 FORMAT(3X,' PEs n.st n.stw n.sts n.g n.gw n.gs')
WRITE( stdout,121) MINVAL(nstp), MINVAL(nstpw), MINVAL(nstps), MINVAL(sstp), MINVAL(sstpw), MINVAL(sstps)
WRITE( stdout,122) MAXVAL(nstp), MAXVAL(nstpw), MAXVAL(nstps), MAXVAL(sstp), MAXVAL(sstpw), MAXVAL(sstps)
WRITE( stdout,121) minval(nstp), minval(nstpw), minval(nstps), minval(sstp), minval(sstpw), minval(sstps)
WRITE( stdout,122) maxval(nstp), maxval(nstpw), maxval(nstps), maxval(sstp), maxval(sstpw), maxval(sstps)
! DO ip = 1, nproc_pool
! IF (ionode) THEN
! WRITE( stdout,120) ip, nstp(ip), nstpw(ip), nstps(ip), sstp(ip), sstpw(ip), sstps(ip)
! END IF
! END DO
IF (ionode) THEN
WRITE( stdout,120) SUM(nstp), SUM(nstpw), SUM(nstps), SUM(sstp), SUM(sstpw), SUM(sstps)
END IF
WRITE( stdout,120) sum(nstp), sum(nstpw), sum(nstps), sum(sstp), sum(sstpw), sum(sstps)
ENDIF
120 FORMAT(3X,7I8)
121 FORMAT(3X,'min ',6I8)
122 FORMAT(3X,'max ',6I8)
@ -824,10 +824,10 @@
!=----------------------------------------------------------------------=
SUBROUTINE sticks_deallocate
IF( ALLOCATED( sticks_owner ) ) DEALLOCATE( sticks_owner )
IF( allocated( sticks_owner ) ) DEALLOCATE( sticks_owner )
RETURN
END SUBROUTINE sticks_deallocate
!=----------------------------------------------------------------------=
END MODULE stick_base

68
Modules/task_groups.f90
View File

@ -50,7 +50,7 @@ SUBROUTINE task_groups_init( dffts )
IMPLICIT NONE
TYPE(fft_dlay_descriptor), INTENT(INOUT) :: dffts
TYPE(fft_dlay_descriptor), INTENT(inout) :: dffts
!----------------------------------
!Local Variables declaration
@ -64,20 +64,20 @@ SUBROUTINE task_groups_init( dffts )
INTEGER :: strd
!
write( stdout, 100 ) nogrp, npgrp
WRITE( stdout, 100 ) nogrp, npgrp
100 format( /,3X,'Task Groups are in USE',/,3X,'groups and procs/group : ',I5,I5 )
100 FORMAT( /,3X,'Task Groups are in USE',/,3X,'groups and procs/group : ',I5,I5 )
!Find maximum chunk of local data concerning coefficients of eigenfunctions in g-space
#if defined __MPI
call MPI_Allgather( dffts%nnr, 1, MPI_INTEGER, nnrsx_vec, 1, MPI_INTEGER, intra_pool_comm, IERR)
CALL MPI_Allgather( dffts%nnr, 1, MPI_INTEGER, nnrsx_vec, 1, MPI_INTEGER, intra_pool_comm, IERR)
strd = maxval( nnrsx_vec( 1:nproc_pool ) )
#else
strd = dffts%nnr
#endif
if( strd /= dffts%nnrx ) call errore( ' task_groups_init ', ' inconsistent nnrx ', 1 )
IF( strd /= dffts%nnrx ) CALL errore( ' task_groups_init ', ' inconsistent nnrx ', 1 )
!-------------------------------------------------------------------------------------
!C. Bekas...TASK GROUP RELATED. FFT DATA STRUCTURES ARE ALREADY DEFINED ABOVE
@ -89,45 +89,45 @@ SUBROUTINE task_groups_init( dffts )
!we choose to do the latter one.
!-------------------------------------------------------------------------------------
!
allocate( dffts%tg_nsw(nproc_pool))
allocate( dffts%tg_npp(nproc_pool))
ALLOCATE( dffts%tg_nsw(nproc_pool))
ALLOCATE( dffts%tg_npp(nproc_pool))
num_sticks = 0
num_planes = 0
do i = 1, nogrp
DO i = 1, nogrp
num_sticks = num_sticks + dffts%nsw( nolist(i) + 1 )
num_planes = num_planes + dffts%npp( nolist(i) + 1 )
enddo
ENDDO
#if defined __MPI
call MPI_ALLGATHER(num_sticks, 1, MPI_INTEGER, dffts%tg_nsw(1), 1, MPI_INTEGER, intra_pool_comm, IERR)
call MPI_ALLGATHER(num_planes, 1, MPI_INTEGER, dffts%tg_npp(1), 1, MPI_INTEGER, intra_pool_comm, IERR)
CALL MPI_ALLGATHER(num_sticks, 1, MPI_INTEGER, dffts%tg_nsw(1), 1, MPI_INTEGER, intra_pool_comm, IERR)
CALL MPI_ALLGATHER(num_planes, 1, MPI_INTEGER, dffts%tg_npp(1), 1, MPI_INTEGER, intra_pool_comm, IERR)
#else
dffts%tg_nsw(1) = num_sticks
dffts%tg_npp(1) = num_planes
#endif
allocate( dffts%tg_snd( nogrp ) )
allocate( dffts%tg_rcv( nogrp ) )
allocate( dffts%tg_psdsp( nogrp ) )
allocate( dffts%tg_usdsp( nogrp ) )
allocate( dffts%tg_rdsp( nogrp ) )
ALLOCATE( dffts%tg_snd( nogrp ) )
ALLOCATE( dffts%tg_rcv( nogrp ) )
ALLOCATE( dffts%tg_psdsp( nogrp ) )
ALLOCATE( dffts%tg_usdsp( nogrp ) )
ALLOCATE( dffts%tg_rdsp( nogrp ) )
dffts%tg_snd(1) = dffts%nr3x * dffts%nsw( me_pool + 1 )
if( dffts%nr3x * dffts%nsw( me_pool + 1 ) > dffts%nnrx ) then
call errore( ' task_groups_init ', ' inconsistent dffts%nnrx ', 1 )
endif
IF( dffts%nr3x * dffts%nsw( me_pool + 1 ) > dffts%nnrx ) THEN
CALL errore( ' task_groups_init ', ' inconsistent dffts%nnrx ', 1 )
ENDIF
dffts%tg_psdsp(1) = 0
dffts%tg_usdsp(1) = 0
dffts%tg_rcv(1) = dffts%nr3x * dffts%nsw( nolist(1) + 1 )
dffts%tg_rdsp(1) = 0
do i = 2, nogrp
DO i = 2, nogrp
dffts%tg_snd(i) = dffts%nr3x * dffts%nsw( me_pool + 1 )
dffts%tg_psdsp(i) = dffts%tg_psdsp(i-1) + dffts%nnrx
dffts%tg_usdsp(i) = dffts%tg_usdsp(i-1) + dffts%tg_snd(i-1)
dffts%tg_rcv(i) = dffts%nr3x * dffts%nsw( nolist(i) + 1 )
dffts%tg_rdsp(i) = dffts%tg_rdsp(i-1) + dffts%tg_rcv(i-1)
enddo
ENDDO
! ALLOCATE( dffts%tg_sca_snd( nproc_pool / nogrp ) )
! ALLOCATE( dffts%tg_sca_rcv( nproc_pool / nogrp ) )
@ -152,7 +152,7 @@ SUBROUTINE task_groups_init( dffts )
dffts%have_task_groups = .true.
return
RETURN
END SUBROUTINE task_groups_init
@ -170,7 +170,7 @@ SUBROUTINE tg_gather( dffts, v, tg_v )
IMPLICIT NONE
TYPE(fft_dlay_descriptor), INTENT(IN) :: dffts
TYPE(fft_dlay_descriptor), INTENT(in) :: dffts
REAL(DP) :: v(:)
REAL(DP) :: tg_v(:)
@ -180,13 +180,13 @@ SUBROUTINE tg_gather( dffts, v, tg_v )
nsiz_tg = dffts%nnrx * nogrp
if( size( tg_v ) < nsiz_tg ) &
call errore( ' tg_gather ', ' tg_v too small ', ( nsiz_tg - size( tg_v ) ) )
IF( size( tg_v ) < nsiz_tg ) &
CALL errore( ' tg_gather ', ' tg_v too small ', ( nsiz_tg - size( tg_v ) ) )
nsiz = dffts%npp( me_pool+1 ) * dffts%nr1x * dffts%nr2x
if( size( v ) < nsiz ) &
call errore( ' tg_gather ', ' v too small ', ( nsiz - size( v ) ) )
IF( size( v ) < nsiz ) &
CALL errore( ' tg_gather ', ' v too small ', ( nsiz - size( v ) ) )
!
! The potential in v is distributed accros all processors
@ -195,24 +195,24 @@ SUBROUTINE tg_gather( dffts, v, tg_v )
!
recv_cnt(1) = dffts%npp( nolist(1) + 1 ) * dffts%nr1x * dffts%nr2x
recv_displ(1) = 0
do i = 2, nogrp
DO i = 2, nogrp
recv_cnt(i) = dffts%npp( nolist(i) + 1 ) * dffts%nr1x * dffts%nr2x
recv_displ(i) = recv_displ(i-1) + recv_cnt(i-1)
enddo
ENDDO
! clean only elements that will not be overwritten
!
do i = recv_displ(nogrp) + recv_cnt( nogrp ) + 1, size( tg_v )
DO i = recv_displ(nogrp) + recv_cnt( nogrp ) + 1, size( tg_v )
tg_v( i ) = 0.0d0
enddo
ENDDO
#if defined (__PARA) && defined (__MPI)
call MPI_Allgatherv( v(1), nsiz, MPI_DOUBLE_PRECISION, &
CALL MPI_Allgatherv( v(1), nsiz, MPI_DOUBLE_PRECISION, &
tg_v(1), recv_cnt, recv_displ, MPI_DOUBLE_PRECISION, ogrp_comm, IERR)
if( ierr /= 0 ) &
call errore( ' tg_gather ', ' MPI_Allgatherv ', abs( ierr ) )
IF( ierr /= 0 ) &
CALL errore( ' tg_gather ', ' MPI_Allgatherv ', abs( ierr ) )
#endif

6
PP/pw2casino.f90
View File

@ -61,7 +61,7 @@ PROGRAM pw2casino
!
prefix = 'pwscf'
CALL get_env( 'ESPRESSO_TMPDIR', outdir )
IF ( TRIM( outdir ) == ' ' ) outdir = './'
IF ( trim( outdir ) == ' ' ) outdir = './'
casino_gather = .false.
ios = 0
IF ( ionode ) THEN
@ -69,10 +69,10 @@ PROGRAM pw2casino
READ (5, inputpp, iostat=ios)
tmp_dir = trimcheck (outdir)
!
END IF
ENDIF
!
CALL mp_bcast( ios, ionode_id )
IF ( ios/=0 ) CALL errore('pw2casino', 'reading inputpp namelist', ABS(ios))
IF ( ios/=0 ) CALL errore('pw2casino', 'reading inputpp namelist', abs(ios))
!
! ... Broadcast variables
!

482
PP/pw2casino_write.f90
View File

@ -48,127 +48,127 @@ SUBROUTINE write_casino_pwfn(gather)
REAL(DP), ALLOCATABLE :: g_l(:,:), g_g(:,:), g2(:)
COMPLEX(DP), ALLOCATABLE :: evc_l(:), evc_g(:)
call init_us_1
call newd
CALL init_us_1
CALL newd
! four times npwx should be enough
allocate (idx (4*npwx) )
allocate (igtog (4*npwx) )
ALLOCATE (idx (4*npwx) )
ALLOCATE (igtog (4*npwx) )
idx(:) = 0
igtog(:) = 0
if( lsda )then
IF( lsda )THEN
nbndup = nbnd
nbnddown = nbnd
nk = nks/2
! nspin = 2
else
ELSE
nbndup = nbnd
nbnddown = 0
nk = nks
! nspin = 1
endif
ENDIF
call calc_energies
CALL calc_energies
do ispin = 1, nspin
do ik = 1, nk
DO ispin = 1, nspin
DO ik = 1, nk
ikk = ik + nk*(ispin-1)
call gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
do ig =1, npw
if( igk(ig) > 4*npwx ) &
call errore ('pw2casino','increase allocation of index', ig)
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
DO ig =1, npw
IF( igk(ig) > 4*npwx ) &
CALL errore ('pw2casino','increase allocation of index', ig)
idx( igk(ig) ) = 1
enddo
enddo
enddo
ENDDO
ENDDO
ENDDO
ngtot = 0
do ig = 1, 4*npwx
if( idx(ig) >= 1 )then
DO ig = 1, 4*npwx
IF( idx(ig) >= 1 )THEN
ngtot = ngtot + 1
igtog(ngtot) = ig
endif
enddo
ENDIF
ENDDO
if(ionode.or..not.gather)then
IF(ionode.or..not.gather)THEN
io = 77
write (6,'(/,5x,''Writing file pwfn.data for program CASINO'')')
call seqopn( 77, 'pwfn.data', 'formatted',exst)
call write_header
endif
WRITE (6,'(/,5x,''Writing file pwfn.data for program CASINO'')')
CALL seqopn( 77, 'pwfn.data', 'formatted',exst)
CALL write_header
ENDIF
allocate ( g_l(3,ngtot), evc_l(ngtot) )
do ig = 1, ngtot
ALLOCATE ( g_l(3,ngtot), evc_l(ngtot) )
DO ig = 1, ngtot
g_l(:,ig) = g(:,igtog(ig))
enddo
ENDDO
if(gather)then
allocate ( ngtot_d(nproc_pool), ngtot_cumsum(nproc_pool) )
call mp_gather( ngtot, ngtot_d, ionode_id, intra_pool_comm )
call mp_bcast( ngtot_d, ionode_id, intra_pool_comm )
IF(gather)THEN
ALLOCATE ( ngtot_d(nproc_pool), ngtot_cumsum(nproc_pool) )
CALL mp_gather( ngtot, ngtot_d, ionode_id, intra_pool_comm )
CALL mp_bcast( ngtot_d, ionode_id, intra_pool_comm )
id = 0
do ip = 1,nproc_pool
DO ip = 1,nproc_pool
ngtot_cumsum(ip) = id
id = id + ngtot_d(ip)
enddo
ENDDO
ngtot_g = id
allocate ( g_g(3,ngtot_g), evc_g(ngtot_g) )
call mp_gather( g_l, g_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
ALLOCATE ( g_g(3,ngtot_g), evc_g(ngtot_g) )
CALL mp_gather( g_l, g_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
if(ionode)then
allocate ( indx(ngtot_g) )
call create_index2(g_g,indx)
call write_gvecs(g_g,indx)
endif
else
allocate ( indx(ngtot) )
call create_index2(g_l,indx)
call write_gvecs(g_l,indx)
endif
IF(ionode)THEN
ALLOCATE ( indx(ngtot_g) )
CALL create_index2(g_g,indx)
CALL write_gvecs(g_g,indx)
ENDIF
ELSE
ALLOCATE ( indx(ngtot) )
CALL create_index2(g_l,indx)
CALL write_gvecs(g_l,indx)
ENDIF
if(ionode.or..not.gather)call write_wfn_head
IF(ionode.or..not.gather)CALL write_wfn_head
do ik = 1, nk
if(ionode.or..not.gather)call write_kpt_head
DO ik = 1, nk
IF(ionode.or..not.gather)CALL write_kpt_head
do ispin = 1, nspin
DO ispin = 1, nspin
ikk = ik + nk*(ispin-1)
if( nks > 1 )then
call gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
call davcio(evc,nwordwfc,iunwfc,ikk,-1)
endif
do ibnd = 1, nbnd
if(ionode.or..not.gather)call write_bnd_head
IF( nks > 1 )THEN
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio(evc,nwordwfc,iunwfc,ikk,-1)
ENDIF
DO ibnd = 1, nbnd
IF(ionode.or..not.gather)CALL write_bnd_head
evc_l(:) = (0.d0, 0d0)
do ig=1, ngtot
DO ig=1, ngtot
! now for all G vectors find the PW coefficient for this k-point
find_ig: do ig7 = 1, npw
if( igk(ig7) == igtog(ig) )then
find_ig: DO ig7 = 1, npw
IF( igk(ig7) == igtog(ig) )THEN
evc_l(ig) = evc(ig7,ibnd)
exit find_ig
endif
enddo find_ig
enddo
if(gather)then
call mp_gather( evc_l, evc_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
ENDIF
ENDDO find_ig
ENDDO
IF(gather)THEN
CALL mp_gather( evc_l, evc_g, ngtot_d, ngtot_cumsum, ionode_id, intra_pool_comm)
if(ionode)call write_wfn_data(evc_g,indx)
else
call write_wfn_data(evc_l,indx)
endif
enddo
enddo
enddo
if(ionode.or..not.gather)close(io)
IF(ionode)CALL write_wfn_data(evc_g,indx)
ELSE
CALL write_wfn_data(evc_l,indx)
ENDIF
ENDDO
ENDDO
ENDDO
IF(ionode.or..not.gather)CLOSE(io)
deallocate (igtog)
deallocate (idx)
DEALLOCATE (igtog)
DEALLOCATE (idx)
deallocate ( g_l, evc_l )
if(gather) deallocate ( ngtot_d, ngtot_cumsum, g_g, evc_g )
if(ionode.or..not.gather) deallocate (indx)
DEALLOCATE ( g_l, evc_l )
IF(gather) DEALLOCATE ( ngtot_d, ngtot_cumsum, g_g, evc_g )
IF(ionode.or..not.gather) DEALLOCATE (indx)
CONTAINS
@ -180,89 +180,89 @@ CONTAINS
REAL(DP) :: charge, etotefield
allocate (aux(nrxx))
call allocate_bec_type ( nkb, nbnd, becp )
ALLOCATE (aux(nrxx))
CALL allocate_bec_type ( nkb, nbnd, becp )
ek = 0.d0
eloc= 0.d0
enl = 0.d0
demet=0.d0
!
do ispin = 1, nspin
DO ispin = 1, nspin
!
! calculate the local contribution to the total energy
!
! bring rho to G-space
!
aux(:) = cmplx( rho%of_r(:,ispin), 0.d0,kind=DP)
call cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
CALL cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
!
do nt=1,ntyp
do ig = 1, ngm
DO nt=1,ntyp
DO ig = 1, ngm
eloc = eloc + vloc(igtongl(ig),nt) * strf(ig,nt) &
* conjg(aux(nl(ig)))
enddo
enddo
ENDDO
ENDDO
do ik = 1, nk
DO ik = 1, nk
ikk = ik + nk*(ispin-1)
call gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
call davcio (evc, nwordwfc, iunwfc, ikk, - 1)
call init_us_2 (npw, igk, xk (1, ikk), vkb)
call calbec ( npw, vkb, evc, becp )
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio (evc, nwordwfc, iunwfc, ikk, - 1)
CALL init_us_2 (npw, igk, xk (1, ikk), vkb)
CALL calbec ( npw, vkb, evc, becp )
!
! -TS term for metals (ifany)
!
if( degauss > 0.0_dp)then
do ibnd = 1, nbnd
IF( degauss > 0.0_dp)THEN
DO ibnd = 1, nbnd
demet = demet + wk (ik) * &
degauss * w1gauss ( (ef-et(ibnd,ik)) / degauss, ngauss)
enddo
endif
ENDDO
ENDIF
!
! calculate the kinetic energy
!
do ibnd = 1, nbnd
do j = 1, npw
DO ibnd = 1, nbnd
DO j = 1, npw
ek = ek + conjg(evc(j,ibnd)) * evc(j,ibnd) * &
g2kin(j) * wg(ibnd,ikk)
enddo
ENDDO
!
! Calculate Non-local energy
!
ijkb0 = 0
do nt = 1, ntyp
do na = 1, nat
if(ityp (na) .eq. nt)then
do ih = 1, nh (nt)
DO nt = 1, ntyp
DO na = 1, nat
IF(ityp (na) .eq. nt)THEN
DO ih = 1, nh (nt)
ikb = ijkb0 + ih
enl=enl+conjg(becp%k(ikb,ibnd))*becp%k(ikb,ibnd) &
*wg(ibnd,ikk)* dvan(ih,ih,nt)
do jh = ( ih + 1 ), nh(nt)
DO jh = ( ih + 1 ), nh(nt)
jkb = ijkb0 + jh
enl=enl + &
(conjg(becp%k(ikb,ibnd))*becp%k(jkb,ibnd)+&
conjg(becp%k(jkb,ibnd))*becp%k(ikb,ibnd))&
* wg(ibnd,ikk) * dvan(ih,jh,nt)
enddo
ENDDO
enddo
ENDDO
ijkb0 = ijkb0 + nh (nt)
endif
enddo
enddo
enddo
enddo
enddo
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
#ifdef __PARA
call mp_sum( eloc, intra_pool_comm )
call mp_sum( ek, intra_pool_comm )
call mp_sum( ek, inter_pool_comm )
call mp_sum( enl, inter_pool_comm )
call mp_sum( demet, inter_pool_comm )
CALL mp_sum( eloc, intra_pool_comm )
CALL mp_sum( ek, intra_pool_comm )
CALL mp_sum( ek, inter_pool_comm )
CALL mp_sum( enl, inter_pool_comm )
CALL mp_sum( demet, inter_pool_comm )
#endif
eloc = eloc * omega
ek = ek * tpiba2
@ -275,23 +275,23 @@ CONTAINS
!
! compute hartree and xc contribution
!
call v_of_rho( rho, rho_core, rhog_core, &
CALL v_of_rho( rho, rho_core, rhog_core, &
ehart, etxc, vtxc, eth, etotefield, charge, vnew )
!
etot=(ek + (etxc-etxcc)+ehart+eloc+enl+ewld)+demet
call deallocate_bec_type (becp)
deallocate (aux)
CALL deallocate_bec_type (becp)
DEALLOCATE (aux)
write (stdout,*) 'Kinetic energy ', ek/e2
write (stdout,*) 'Local energy ', eloc/e2
write (stdout,*) 'Non-Local energy ', enl/e2
write (stdout,*) 'Ewald energy ', ewld/e2
write (stdout,*) 'xc contribution ',(etxc-etxcc)/e2
write (stdout,*) 'hartree energy ', ehart/e2
if( degauss > 0.0_dp ) &
write (stdout,*) 'Smearing (-TS) ', demet/e2
write (stdout,*) 'Total energy ', etot/e2
WRITE (stdout,*) 'Kinetic energy ', ek/e2
WRITE (stdout,*) 'Local energy ', eloc/e2
WRITE (stdout,*) 'Non-Local energy ', enl/e2
WRITE (stdout,*) 'Ewald energy ', ewld/e2
WRITE (stdout,*) 'xc contribution ',(etxc-etxcc)/e2
WRITE (stdout,*) 'hartree energy ', ehart/e2
IF( degauss > 0.0_dp ) &
WRITE (stdout,*) 'Smearing (-TS) ', demet/e2
WRITE (stdout,*) 'Total energy ', etot/e2
END SUBROUTINE calc_energies
@ -299,62 +299,62 @@ CONTAINS
SUBROUTINE write_header
INTEGER j, na, nt, at_num
write(io,'(a)') title
write(io,'(a)')
write(io,'(a)') ' BASIC INFO'
write(io,'(a)') ' ----------'
write(io,'(a)') ' Generated by:'
write(io,'(a)') ' PWSCF'
write(io,'(a)') ' Method:'
write(io,'(a)') ' DFT'
write(io,'(a)') ' DFT Functional:'
write(io,'(a)') ' unknown'
write(io,'(a)') ' Pseudopotential'
write(io,'(a)') ' unknown'
write(io,'(a)') ' Plane wave cutoff (au)'
write(io,*) ecutwfc/2
write(io,'(a)') ' Spin polarized:'
write(io,*)lsda
if( degauss > 0.0_dp )then
write(io,'(a)') ' Total energy (au per primitive cell; includes -TS term)'
write(io,*)etot/e2, demet/e2
else
write(io,'(a)') ' Total energy (au per primitive cell)'
write(io,*)etot/e2
endif
write(io,'(a)') ' Kinetic energy (au per primitive cell)'
write(io,*)ek/e2
write(io,'(a)') ' Local potential energy (au per primitive cell)'
write(io,*)eloc/e2
write(io,'(a)') ' Non local potential energy(au per primitive cel)'
write(io,*)enl/e2
write(io,'(a)') ' Electron electron energy (au per primitive cell)'
write(io,*)ehart/e2
write(io,'(a)') ' Ion ion energy (au per primitive cell)'
write(io,*)ewld/e2
write(io,'(a)') ' Number of electrons per primitive cell'
write(io,*)nint(nelec)
WRITE(io,'(a)') title
WRITE(io,'(a)')
WRITE(io,'(a)') ' BASIC INFO'
WRITE(io,'(a)') ' ----------'
WRITE(io,'(a)') ' Generated by:'
WRITE(io,'(a)') ' PWSCF'
WRITE(io,'(a)') ' Method:'
WRITE(io,'(a)') ' DFT'
WRITE(io,'(a)') ' DFT Functional:'
WRITE(io,'(a)') ' unknown'
WRITE(io,'(a)') ' Pseudopotential'
WRITE(io,'(a)') ' unknown'
WRITE(io,'(a)') ' Plane wave cutoff (au)'
WRITE(io,*) ecutwfc/2
WRITE(io,'(a)') ' Spin polarized:'
WRITE(io,*)lsda
IF( degauss > 0.0_dp )THEN
WRITE(io,'(a)') ' Total energy (au per primitive cell; includes -TS term)'
WRITE(io,*)etot/e2, demet/e2
ELSE
WRITE(io,'(a)') ' Total energy (au per primitive cell)'
WRITE(io,*)etot/e2
ENDIF
WRITE(io,'(a)') ' Kinetic energy (au per primitive cell)'
WRITE(io,*)ek/e2
WRITE(io,'(a)') ' Local potential energy (au per primitive cell)'
WRITE(io,*)eloc/e2
WRITE(io,'(a)') ' Non local potential energy(au per primitive cel)'
WRITE(io,*)enl/e2
WRITE(io,'(a)') ' Electron electron energy (au per primitive cell)'
WRITE(io,*)ehart/e2
WRITE(io,'(a)') ' Ion ion energy (au per primitive cell)'
WRITE(io,*)ewld/e2
WRITE(io,'(a)') ' Number of electrons per primitive cell'
WRITE(io,*)nint(nelec)
! uncomment the following ifyou want the Fermi energy - KN 2/4/09
! write(io,'(a)') ' Fermi energy (au)'
! write(io,*) ef/e2
write(io,'(a)') ' '
write(io,'(a)') ' GEOMETRY'
write(io,'(a)') ' -------- '
write(io,'(a)') ' Number of atoms per primitive cell '
write(io,*) nat
write(io,'(a)')' Atomic number and position of the atoms(au) '
do na = 1, nat
WRITE(io,'(a)') ' '
WRITE(io,'(a)') ' GEOMETRY'
WRITE(io,'(a)') ' -------- '
WRITE(io,'(a)') ' Number of atoms per primitive cell '
WRITE(io,*) nat
WRITE(io,'(a)')' Atomic number and position of the atoms(au) '
DO na = 1, nat
nt = ityp(na)
at_num = atomic_number(trim(atm(nt)))
write(io,'(i6,3f20.14)') at_num, (alat*tau(j,na),j=1,3)
enddo
write(io,'(a)') ' Primitive lattice vectors (au) '
write(io,100) alat*at(1,1), alat*at(2,1), alat*at(3,1)
write(io,100) alat*at(1,2), alat*at(2,2), alat*at(3,2)
write(io,100) alat*at(1,3), alat*at(2,3), alat*at(3,3)
write(io,'(a)') ' '
WRITE(io,'(i6,3f20.14)') at_num, (alat*tau(j,na),j=1,3)
ENDDO
WRITE(io,'(a)') ' Primitive lattice vectors (au) '
WRITE(io,100) alat*at(1,1), alat*at(2,1), alat*at(3,1)
WRITE(io,100) alat*at(1,2), alat*at(2,2), alat*at(3,2)
WRITE(io,100) alat*at(1,3), alat*at(2,3), alat*at(3,3)
WRITE(io,'(a)') ' '
100 format (3(1x,f20.15))
100 FORMAT (3(1x,f20.15))
END SUBROUTINE write_header
@ -364,48 +364,48 @@ CONTAINS
INTEGER,INTENT(in) :: indx(:)
INTEGER ig
write(io,'(a)') ' G VECTORS'
write(io,'(a)') ' ---------'
write(io,'(a)') ' Number of G-vectors'
write(io,*) size(g,2)
write(io,'(a)') ' Gx Gy Gz (au)'
do ig = 1, size(g,2)
write(io,100) tpi/alat*g(1,indx(ig)), tpi/alat*g(2,indx(ig)), &
WRITE(io,'(a)') ' G VECTORS'
WRITE(io,'(a)') ' ---------'
WRITE(io,'(a)') ' Number of G-vectors'
WRITE(io,*) size(g,2)
WRITE(io,'(a)') ' Gx Gy Gz (au)'
DO ig = 1, size(g,2)
WRITE(io,100) tpi/alat*g(1,indx(ig)), tpi/alat*g(2,indx(ig)), &
tpi/alat*g(3,indx(ig))
enddo
ENDDO
100 format (3(1x,f20.15))
100 FORMAT (3(1x,f20.15))
write(io,'(a)') ' '
WRITE(io,'(a)') ' '
END SUBROUTINE write_gvecs
SUBROUTINE write_wfn_head
write(io,'(a)') ' WAVE FUNCTION'
write(io,'(a)') ' -------------'
write(io,'(a)') ' Number of k-points'
write(io,*) nk
WRITE(io,'(a)') ' WAVE FUNCTION'
WRITE(io,'(a)') ' -------------'
WRITE(io,'(a)') ' Number of k-points'
WRITE(io,*) nk
END SUBROUTINE write_wfn_head
SUBROUTINE write_kpt_head
INTEGER j
write(io,'(a)') ' k-point # ; # of bands (up spin/down spin); &
WRITE(io,'(a)') ' k-point # ; # of bands (up spin/down spin); &
& k-point coords (au)'
write(io,'(3i4,3f20.16)') ik, nbndup, nbnddown, &
WRITE(io,'(3i4,3f20.16)') ik, nbndup, nbnddown, &
(tpi/alat*xk(j,ik),j=1,3)
END SUBROUTINE write_kpt_head
SUBROUTINE write_bnd_head
! KN: if you want to print occupancies, replace these two lines ...
write(io,'(a)') ' Band, spin, eigenvalue (au)'
write(io,*) ibnd, ispin, et(ibnd,ikk)/e2
WRITE(io,'(a)') ' Band, spin, eigenvalue (au)'
WRITE(io,*) ibnd, ispin, et(ibnd,ikk)/e2
! ...with the following two - KN 2/4/09
! write(io,'(a)') ' Band, spin, eigenvalue (au), occupation number'
! write(io,*) ibnd, ispin, et(ibnd,ikk)/e2, wg(ibnd,ikk)/wk(ikk)
write(io,'(a)') ' Eigenvectors coefficients'
WRITE(io,'(a)') ' Eigenvectors coefficients'
END SUBROUTINE write_bnd_head
@ -414,9 +414,9 @@ CONTAINS
INTEGER,INTENT(in) :: indx(:)
INTEGER ig
do ig=1, size(evc,1)
write(io,*)evc(indx(ig))
enddo
DO ig=1, size(evc,1)
WRITE(io,*)evc(indx(ig))
ENDDO
END SUBROUTINE write_wfn_data
@ -425,10 +425,10 @@ CONTAINS
INTEGER,INTENT(out) :: x_index(size(y,2))
DOUBLE PRECISION y2(size(y,2))
INTEGER i
do i = 1,size(y,2)
DO i = 1,size(y,2)
y2(i) = sum(y(:,i)**2)
enddo
call create_index(y2,x_index)
ENDDO
CALL create_index(y2,x_index)
END SUBROUTINE create_index2
@ -445,77 +445,77 @@ CONTAINS
INTEGER n,i,x_indexj,ir,itemp,j,jstack,k,l,lp1,istack(stacksize)
DOUBLE PRECISION yj
n=size(x_index)
do j=1,n
DO j=1,n
x_index(j)=j
enddo ! j
if(n<=1)return
ENDDO ! j
IF(n<=1)RETURN
jstack=0
l=1
ir=n
do
if(ir-l<ins_sort_thresh)then
jloop : do j=l+1,ir
DO
IF(ir-l<ins_sort_thresh)THEN
jloop : DO j=l+1,ir
x_indexj=x_index(j) ; yj=y(x_indexj)
do i=j-1,l,-1
if(y(x_index(i))<=yj)then
DO i=j-1,l,-1
IF(y(x_index(i))<=yj)THEN
x_index(i+1)=x_indexj
cycle jloop
endif! y(x_index(i))<=yj
CYCLE jloop
ENDIF! y(x_index(i))<=yj
x_index(i+1)=x_index(i)
enddo ! i
ENDDO ! i
x_index(l)=x_indexj
enddo jloop ! j
if(jstack==0)return
ENDDO jloop ! j
IF(jstack==0)RETURN
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
ELSE
k=(l+ir)/2
lp1=l+1
itemp=x_index(k) ; x_index(k)=x_index(lp1) ; x_index(lp1)=itemp
if(y(x_index(l))>y(x_index(ir)))then
IF(y(x_index(l))>y(x_index(ir)))THEN
itemp=x_index(l) ; x_index(l)=x_index(ir) ; x_index(ir)=itemp
endif
if(y(x_index(lp1))>y(x_index(ir)))then
ENDIF
IF(y(x_index(lp1))>y(x_index(ir)))THEN
itemp=x_index(lp1) ; x_index(lp1)=x_index(ir) ; x_index(ir)=itemp
endif
if(y(x_index(l))>y(x_index(lp1)))then
ENDIF
IF(y(x_index(l))>y(x_index(lp1)))THEN
itemp=x_index(l) ; x_index(l)=x_index(lp1) ; x_index(lp1)=itemp
endif
ENDIF
i=lp1
j=ir
x_indexj=x_index(lp1)
yj=y(x_indexj)
do
do
DO
DO
i=i+1
if(y(x_index(i))>=yj)exit
enddo ! i
do
IF(y(x_index(i))>=yj)exit
ENDDO ! i
DO
j=j-1
if(y(x_index(j))<=yj)exit
enddo ! j
if(j<i)exit
IF(y(x_index(j))<=yj)exit
ENDDO ! j
IF(j<i)exit
itemp=x_index(i) ; x_index(i)=x_index(j) ; x_index(j)=itemp
enddo
ENDDO
x_index(lp1)=x_index(j)
x_index(j)=x_indexj
jstack=jstack+2
if(jstack>stacksize)then
write(6,*)'stacksize is too small.'
stop
endif! jstack>stacksize
if(ir-i+1>=j-l)then
IF(jstack>stacksize)THEN
WRITE(6,*)'stacksize is too small.'
STOP
ENDIF! jstack>stacksize
IF(ir-i+1>=j-l)THEN
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
ELSE
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif! ir-i+1>=j-l
endif! ir-l<ins_sort_thresh
enddo
ENDIF! ir-i+1>=j-l
ENDIF! ir-l<ins_sort_thresh
ENDDO
END SUBROUTINE create_index
END SUBROUTINE write_casino_pwfn

108
PW/allocate_fft.f90
View File

@ -34,80 +34,80 @@ SUBROUTINE allocate_fft
!
! determines the data structure for fft arrays
!
call data_structure( gamma_only )
CALL data_structure( gamma_only )
!
! DCC
if( do_coarse ) call data_structure_coarse( gamma_only, nr1,nr2,nr3, ecutwfc )
IF( do_coarse ) CALL data_structure_coarse( gamma_only, nr1,nr2,nr3, ecutwfc )
!
if (nrxx.lt.ngm) then
write( stdout, '(/,4x," nr1=",i4," nr2= ", i4, " nr3=",i4, &
IF (nrxx.lt.ngm) THEN
WRITE( stdout, '(/,4x," nr1=",i4," nr2= ", i4, " nr3=",i4, &
&" nrxx = ",i8," ngm=",i8)') nr1, nr2, nr3, nrxx, ngm
call errore ('allocate_fft', 'the nr"s are too small!', 1)
CALL errore ('allocate_fft', 'the nr"s are too small!', 1)
endif
if (nrxxs.lt.ngms) then
write( stdout, '(/,4x," nr1s=",i4," nr2s= ", i4, " nr3s=",i4, &
ENDIF
IF (nrxxs.lt.ngms) THEN
WRITE( stdout, '(/,4x," nr1s=",i4," nr2s= ", i4, " nr3s=",i4, &
&" nrxxs = ",i8," ngms=",i8)') nr1s, nr2s, nr3s, nrxxs, ngms
call errore ('allocate_fft', 'the nrs"s are too small!', 1)
CALL errore ('allocate_fft', 'the nrs"s are too small!', 1)
endif
if (ngm <= 0) call errore ('allocate_fft', 'wrong ngm', 1)
if (ngms <= 0) call errore ('allocate_fft', 'wrong ngms', 1)
if (nrxx <= 0) call errore ('allocate_fft', 'wrong nrxx', 1)
if (nrxxs<= 0) call errore ('allocate_fft', 'wrong nrxxs', 1)
if (nspin<= 0) call errore ('allocate_fft', 'wrong nspin', 1)
ENDIF
IF (ngm <= 0) CALL errore ('allocate_fft', 'wrong ngm', 1)
IF (ngms <= 0) CALL errore ('allocate_fft', 'wrong ngms', 1)
IF (nrxx <= 0) CALL errore ('allocate_fft', 'wrong nrxx', 1)
IF (nrxxs<= 0) CALL errore ('allocate_fft', 'wrong nrxxs', 1)
IF (nspin<= 0) CALL errore ('allocate_fft', 'wrong nspin', 1)
!
! Allocate memory for all kind of stuff.
!
allocate (g( 3, ngm))
allocate (gg( ngm))
allocate (nl( ngm))
if (gamma_only) allocate (nlm(ngm))
allocate (igtongl( ngm))
allocate (ig1( ngm))
allocate (ig2( ngm))
allocate (ig3( ngm))
ALLOCATE (g( 3, ngm))
ALLOCATE (gg( ngm))
ALLOCATE (nl( ngm))
IF (gamma_only) ALLOCATE (nlm(ngm))
ALLOCATE (igtongl( ngm))
ALLOCATE (ig1( ngm))
ALLOCATE (ig2( ngm))
ALLOCATE (ig3( ngm))
call create_scf_type(rho)
call create_scf_type(v, do_not_allocate_becsum = .true.)
call create_scf_type(vnew, do_not_allocate_becsum = .true.)
allocate (vltot( nrxx))
allocate (rho_core( nrxx))
if (dft_is_meta() ) then
allocate ( kedtau(nrxxs,nspin) )
else
allocate ( kedtau(1,nspin) )
endif
allocate( rhog_core( ngm ) )
allocate (psic( nrxx))
allocate (vrs( nrxx, nspin))
if (doublegrid) then
allocate (nls( ngms))
if (gamma_only) allocate (nlsm(ngm))
else
CALL create_scf_type(rho)
CALL create_scf_type(v, do_not_allocate_becsum = .true.)
CALL create_scf_type(vnew, do_not_allocate_becsum = .true.)
ALLOCATE (vltot( nrxx))
ALLOCATE (rho_core( nrxx))
IF (dft_is_meta() ) THEN
ALLOCATE ( kedtau(nrxxs,nspin) )
ELSE
ALLOCATE ( kedtau(1,nspin) )
ENDIF
ALLOCATE( rhog_core( ngm ) )
ALLOCATE (psic( nrxx))
ALLOCATE (vrs( nrxx, nspin))
IF (doublegrid) THEN
ALLOCATE (nls( ngms))
IF (gamma_only) ALLOCATE (nlsm(ngm))
ELSE
nls => nl
if (gamma_only) nlsm=> nlm
endif
IF (gamma_only) nlsm=> nlm
ENDIF
! DCC
if( do_coarse ) then
allocate (nlc( ngmc))
if (gamma_only) allocate (nlcm(ngmc))
endif
IF( do_coarse ) THEN
ALLOCATE (nlc( ngmc))
IF (gamma_only) ALLOCATE (nlcm(ngmc))
ENDIF
if (noncolin) allocate (psic_nc( nrxx, npol))
IF (noncolin) ALLOCATE (psic_nc( nrxx, npol))
if ( ((report.ne.0).or.(i_cons.ne.0)) .and. (noncolin.and.domag) .or. (i_cons.eq.1) ) then
IF ( ((report.ne.0).or.(i_cons.ne.0)) .and. (noncolin.and.domag) .or. (i_cons.eq.1) ) THEN
!
! In order to print out local quantities, integrated around the atoms,
! we need the following variables
!
allocate(pointlist(nrxx))
allocate(factlist(nrxx))
allocate(r_loc(nat))
call make_pointlists
endif
ALLOCATE(pointlist(nrxx))
ALLOCATE(factlist(nrxx))
ALLOCATE(r_loc(nat))
CALL make_pointlists
ENDIF
return
RETURN
END SUBROUTINE allocate_fft

134
PW/becmod.f90
View File

@ -22,13 +22,13 @@ MODULE becmod
SAVE
!
#ifdef __STD_F95
#define __ALLOCATABLE pointer
#define __ALLOCATABLE POINTER
#define __ALLOCATED associated
#else
#define __ALLOCATABLE allocatable
#define __ALLOCATABLE ALLOCATABLE
#define __ALLOCATED allocated
#endif
TYPE bec_type
TYPE bec_type
REAL(DP), __ALLOCATABLE :: r(:,:) ! appropriate for gammaonly
COMPLEX(DP),__ALLOCATABLE :: k(:,:) ! appropriate for generic k
COMPLEX(DP),__ALLOCATABLE :: nc(:,:,:) ! appropriate for noncolin
@ -39,7 +39,7 @@ MODULE becmod
PRIVATE
REAL(DP), ALLOCATABLE :: &
becp_r(:,:) ! <beta|psi> for real (at Gamma) wavefunctions
becp_r(:,:) ! <beta|psi> for real (at Gamma) wavefunctions
COMPLEX(DP), ALLOCATABLE :: &
becp_k (:,:), & ! as above for complex wavefunctions
becp_nc(:,:,:) ! as above for spinors
@ -52,7 +52,7 @@ MODULE becmod
!
PUBLIC :: bec_type, becp, allocate_bec_type, deallocate_bec_type, calbec, &
beccopy
! becp_, becp_k, becp_nc,
! becp_, becp_k, becp_nc,
!
CONTAINS
!-----------------------------------------------------------------------
@ -60,24 +60,24 @@ CONTAINS
!-----------------------------------------------------------------------
!
IMPLICIT NONE
COMPLEX (DP), INTENT (IN) :: beta(:,:), psi(:,:)
TYPE (bec_type), INTENT (INOUT) :: betapsi ! NB: must be INOUT otherwise
COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:)
TYPE (bec_type), INTENT (inout) :: betapsi ! NB: must be INOUT otherwise
! the allocatd array is lost
INTEGER, INTENT (IN) :: npw
INTEGER, INTENT (in) :: npw
INTEGER, OPTIONAL :: nbnd
!
INTEGER :: local_nbnd
!
IF ( PRESENT (nbnd) ) THEN
IF ( present (nbnd) ) THEN
local_nbnd = nbnd
ELSE
local_nbnd = SIZE ( psi, 2)
END IF
local_nbnd = size ( psi, 2)
ENDIF
IF ( gamma_only ) THEN
IF ( gamma_only ) THEN
CALL calbec_gamma ( npw, beta, psi, betapsi%r, local_nbnd )
!
ELSE IF ( noncolin) THEN
ELSEIF ( noncolin) THEN
!
CALL calbec_nc ( npw, beta, psi, betapsi%nc, local_nbnd )
!
@ -85,7 +85,7 @@ CONTAINS
!
CALL calbec_k ( npw, beta, psi, betapsi%k, local_nbnd )
!
END IF
ENDIF
!
RETURN
!
@ -94,7 +94,7 @@ CONTAINS
SUBROUTINE calbec_gamma ( npw, beta, psi, betapsi, nbnd )
!-----------------------------------------------------------------------
!
! ... matrix times matrix with summation index (k=1,npw) running on
! ... matrix times matrix with summation index (k=1,npw) running on
! ... half of the G-vectors or PWs - assuming k=0 is the G=0 component:
! ... betapsi(i,j) = 2Re(\sum_k beta^*(i,k)psi(k,j)) + beta^*(i,0)psi(0,j)
!
@ -102,30 +102,30 @@ CONTAINS
USE mp, ONLY : mp_sum
IMPLICIT NONE
COMPLEX (DP), INTENT (IN) :: beta(:,:), psi(:,:)
REAL (DP), INTENT (OUT) :: betapsi(:,:)
INTEGER, INTENT (IN) :: npw
COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:)
REAL (DP), INTENT (out) :: betapsi(:,:)
INTEGER, INTENT (in) :: npw
INTEGER, OPTIONAL :: nbnd
!
INTEGER :: nkb, npwx, m
!
nkb = SIZE (beta, 2)
nkb = size (beta, 2)
IF ( nkb == 0 ) RETURN
!
CALL start_clock( 'calbec' )
npwx= SIZE (beta, 1)
IF ( npwx /= SIZE (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1)
npwx= size (beta, 1)
IF ( npwx /= size (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1)
IF ( npwx < npw ) CALL errore ('calbec', 'size mismatch', 2)
IF ( PRESENT (nbnd) ) THEN
IF ( present (nbnd) ) THEN
m = nbnd
ELSE
m = SIZE ( psi, 2)
END IF
m = size ( psi, 2)
ENDIF
#ifdef DEBUG
write (*,*) 'calbec gamma'
write (*,*) nkb, SIZE (betapsi,1) , m , SIZE (betapsi, 2)
WRITE (*,*) 'calbec gamma'
WRITE (*,*) nkb, size (betapsi,1) , m , size (betapsi, 2)
#endif
IF ( nkb /= SIZE (betapsi,1) .OR. m > SIZE (betapsi, 2) ) &
IF ( nkb /= size (betapsi,1) .or. m > size (betapsi, 2) ) &
CALL errore ('calbec', 'size mismatch', 3)
!
IF ( m == 1 ) THEN
@ -141,7 +141,7 @@ CONTAINS
IF ( gstart == 2 ) &
CALL DGER( nkb, m, -1.0_DP, beta, 2*npwx, psi, 2*npwx, betapsi, nkb )
!
END IF
ENDIF
!
CALL mp_sum( betapsi( :, 1:m ), intra_pool_comm )
!
@ -155,37 +155,37 @@ CONTAINS
SUBROUTINE calbec_k ( npw, beta, psi, betapsi, nbnd )
!-----------------------------------------------------------------------
!
! ... matrix times matrix with summation index (k=1,npw) running on
! ... matrix times matrix with summation index (k=1,npw) running on
! ... G-vectors or PWs : betapsi(i,j) = \sum_k beta^*(i,k) psi(k,j)
!
USE mp_global, ONLY : intra_pool_comm
USE mp, ONLY : mp_sum
IMPLICIT NONE
COMPLEX (DP), INTENT (IN) :: beta(:,:), psi(:,:)
COMPLEX (DP), INTENT (OUT) :: betapsi(:,:)
INTEGER, INTENT (IN) :: npw
COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:)
COMPLEX (DP), INTENT (out) :: betapsi(:,:)
INTEGER, INTENT (in) :: npw
INTEGER, OPTIONAL :: nbnd
!
INTEGER :: nkb, npwx, m
!
nkb = SIZE (beta, 2)
nkb = size (beta, 2)
IF ( nkb == 0 ) RETURN
!
CALL start_clock( 'calbec' )
npwx= SIZE (beta, 1)
IF ( npwx /= SIZE (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1)
npwx= size (beta, 1)
IF ( npwx /= size (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1)
IF ( npwx < npw ) CALL errore ('calbec', 'size mismatch', 2)
IF ( PRESENT (nbnd) ) THEN
IF ( present (nbnd) ) THEN
m = nbnd
ELSE
m = SIZE ( psi, 2)
END IF
m = size ( psi, 2)
ENDIF
#ifdef DEBUG
write (*,*) 'calbec k'
write (*,*) nkb, SIZE (betapsi,1) , m , SIZE (betapsi, 2)
WRITE (*,*) 'calbec k'
WRITE (*,*) nkb, size (betapsi,1) , m , size (betapsi, 2)
#endif
IF ( nkb /= SIZE (betapsi,1) .OR. m > SIZE (betapsi, 2) ) &
IF ( nkb /= size (betapsi,1) .or. m > size (betapsi, 2) ) &
CALL errore ('calbec', 'size mismatch', 3)
!
IF ( m == 1 ) THEN
@ -198,7 +198,7 @@ CONTAINS
CALL ZGEMM( 'C', 'N', nkb, m, npw, (1.0_DP,0.0_DP), &
beta, npwx, psi, npwx, (0.0_DP,0.0_DP), betapsi, nkb )
!
END IF
ENDIF
!
CALL mp_sum( betapsi( :, 1:m ), intra_pool_comm )
!
@ -212,7 +212,7 @@ CONTAINS
SUBROUTINE calbec_nc ( npw, beta, psi, betapsi, nbnd )
!-----------------------------------------------------------------------
!
! ... matrix times matrix with summation index (k below) running on
! ... matrix times matrix with summation index (k below) running on
! ... G-vectors or PWs corresponding to two different polarizations:
! ... betapsi(i,1,j) = \sum_k=1,npw beta^*(i,k) psi(k,j)
! ... betapsi(i,2,j) = \sum_k=1,npw beta^*(i,k) psi(k+npwx,j)
@ -221,34 +221,34 @@ CONTAINS
USE mp, ONLY : mp_sum
IMPLICIT NONE
COMPLEX (DP), INTENT (IN) :: beta(:,:), psi(:,:)
COMPLEX (DP), INTENT (OUT) :: betapsi(:,:,:)
INTEGER, INTENT (IN) :: npw
COMPLEX (DP), INTENT (in) :: beta(:,:), psi(:,:)
COMPLEX (DP), INTENT (out) :: betapsi(:,:,:)
INTEGER, INTENT (in) :: npw
INTEGER, OPTIONAL :: nbnd
!
INTEGER :: nkb, npwx, npol, m
!
nkb = SIZE (beta, 2)
nkb = size (beta, 2)
IF ( nkb == 0 ) RETURN
!
call start_clock ('calbec')
npwx= SIZE (beta, 1)
IF ( 2*npwx /= SIZE (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1)
CALL start_clock ('calbec')
npwx= size (beta, 1)
IF ( 2*npwx /= size (psi, 1) ) CALL errore ('calbec', 'size mismatch', 1)
IF ( npwx < npw ) CALL errore ('calbec', 'size mismatch', 2)
IF ( PRESENT (nbnd) ) THEN
IF ( present (nbnd) ) THEN
m = nbnd
ELSE
m = SIZE ( psi, 2)
END IF
npol= SIZE (betapsi, 2)
m = size ( psi, 2)
ENDIF
npol= size (betapsi, 2)
#ifdef DEBUG
write (*,*) 'calbec nc'
write (*,*) nkb, SIZE (betapsi,1) , m , SIZE (betapsi, 3)
WRITE (*,*) 'calbec nc'
WRITE (*,*) nkb, size (betapsi,1) , m , size (betapsi, 3)
#endif
IF ( nkb /= SIZE (betapsi,1) .OR. m > SIZE (betapsi, 3) ) &
IF ( nkb /= size (betapsi,1) .or. m > size (betapsi, 3) ) &
CALL errore ('calbec', 'size mismatch', 3)
!
call ZGEMM ('C', 'N', nkb, m*npol, npw, (1.0_DP, 0.0_DP), beta, &
CALL ZGEMM ('C', 'N', nkb, m*npol, npw, (1.0_DP, 0.0_DP), beta, &
npwx, psi, npwx, (0.0_DP, 0.0_DP), betapsi, nkb)
!
CALL mp_sum( betapsi( :, :, 1:m ), intra_pool_comm )
@ -264,18 +264,18 @@ CONTAINS
!-----------------------------------------------------------------------
IMPLICIT NONE
TYPE (bec_type) :: bec
INTEGER, INTENT (IN) :: nkb, nbnd
INTEGER, INTENT (in) :: nkb, nbnd
!
#ifdef __STD_F95
! otherwise they might still be defined
! nullify (bec%r, bec%k, bec%nc)
#endif
IF ( gamma_only ) THEN
IF ( gamma_only ) THEN
!
ALLOCATE( bec%r( nkb, nbnd ) )
bec%r(:,:)=0.0D0
!
ELSE IF ( noncolin) THEN
ELSEIF ( noncolin) THEN
!
ALLOCATE( bec%nc( nkb, npol, nbnd ) )
bec%nc(:,:,:)=(0.0D0,0.0D0)
@ -285,7 +285,7 @@ CONTAINS
ALLOCATE( bec%k( nkb, nbnd ) )
bec%k(:,:)=(0.0D0,0.0D0)
!
END IF
ENDIF
!
RETURN
!
@ -298,9 +298,9 @@ CONTAINS
IMPLICIT NONE
TYPE (bec_type) :: bec
!
if (__ALLOCATED(bec%r)) deallocate(bec%r)
if (__ALLOCATED(bec%nc)) deallocate(bec%nc)
if (__ALLOCATED(bec%k)) deallocate(bec%k)
IF (__ALLOCATED(bec%r)) DEALLOCATE(bec%r)
IF (__ALLOCATED(bec%nc)) DEALLOCATE(bec%nc)
IF (__ALLOCATED(bec%k)) DEALLOCATE(bec%k)
!
RETURN
!
@ -308,9 +308,9 @@ CONTAINS
SUBROUTINE beccopy(bec, bec1, nkb, nbnd)
IMPLICIT NONE
TYPE(bec_type), INTENT(IN) :: bec
TYPE(bec_type), INTENT(in) :: bec
TYPE(bec_type) :: bec1
INTEGER, INTENT(IN) :: nkb, nbnd
INTEGER, INTENT(in) :: nkb, nbnd
IF (gamma_only) THEN
CALL dcopy(nkb*nbnd, bec1%r, 1, bec%r, 1)

178
PW/data_structure.f90
View File

@ -7,7 +7,7 @@
!
!
!-----------------------------------------------------------------------
subroutine data_structure( lgamma )
SUBROUTINE data_structure( lgamma )
!-----------------------------------------------------------------------
! this routine sets the data structure for the fft arrays
! (both the smooth and the hard mesh)
@ -31,34 +31,34 @@ subroutine data_structure( lgamma )
!
USE task_groups, ONLY : task_groups_init
!
implicit none
logical, intent(in) :: lgamma
integer :: n1, n2, n3, i1, i2, i3
IMPLICIT NONE
LOGICAL, INTENT(in) :: lgamma
INTEGER :: n1, n2, n3, i1, i2, i3
! counters on G space
!
real(DP) :: amod
! modulus of G vectors
integer, allocatable :: stw(:,:)
INTEGER, ALLOCATABLE :: stw(:,:)
! sticks maps
integer :: ub(3), lb(3)
INTEGER :: ub(3), lb(3)
! upper and lower bounds for maps
real(DP) :: gkcut
! cut-off for the wavefunctions
integer :: ncplane, nxx
integer :: ncplanes, nxxs
INTEGER :: ncplane, nxx
INTEGER :: ncplanes, nxxs
#ifdef __PARA
integer, allocatable :: st(:,:), sts(:,:)
INTEGER, ALLOCATABLE :: st(:,:), sts(:,:)
! sticks maps
integer, allocatable :: ngc (:), ngcs (:), ngkc (:)
integer :: ncp (nproc), nct, nkcp (nproc), ncts, ncps(nproc)
integer :: ngp (nproc), ngps(nproc), ngkp (nproc), ncp_(nproc),&
INTEGER, ALLOCATABLE :: ngc (:), ngcs (:), ngkc (:)
INTEGER :: ncp (nproc), nct, nkcp (nproc), ncts, ncps(nproc)
INTEGER :: ngp (nproc), ngps(nproc), ngkp (nproc), ncp_(nproc),&
i, j, jj, idum
! nxx ! local fft data dim
@ -67,16 +67,16 @@ subroutine data_structure( lgamma )
! ncp(nproc), &! number of (density) columns per proc
logical :: tk = .TRUE.
LOGICAL :: tk = .true.
! map type: true for full space sticks map, false for half space sticks map
integer, allocatable :: in1(:), in2(:), idx(:)
INTEGER, ALLOCATABLE :: in1(:), in2(:), idx(:)
! sticks coordinates
!
! Subroutine body
!
tk = .NOT. lgamma
tk = .not. lgamma
!
! set the values of fft arrays
@ -100,17 +100,17 @@ subroutine data_structure( lgamma )
!
! check the number of plane per process
!
if ( nr3 < nproc_pool ) &
call infomsg ('data_structure', 'some processors have no planes ')
IF ( nr3 < nproc_pool ) &
CALL infomsg ('data_structure', 'some processors have no planes ')
if ( nr3s < nproc_pool ) &
call infomsg ('data_structure', 'some processors have no smooth planes ')
IF ( nr3s < nproc_pool ) &
CALL infomsg ('data_structure', 'some processors have no smooth planes ')
!
! compute gkcut calling an internal procedure
!
call calculate_gkcut()
CALL calculate_gkcut()
#ifdef DEBUG
WRITE( stdout, '(5x,"ecutrho & ecutwfc",2f12.2)') tpiba2 * gcutm, &
@ -141,20 +141,20 @@ subroutine data_structure( lgamma )
CALL sticks_maps( tk, ub, lb, bg(:,1), bg(:,2), bg(:,3), gcutm, gkcut, gcutms, st, stw, sts )
nct = COUNT( st > 0 )
ncts = COUNT( sts > 0 )
nct = count( st > 0 )
ncts = count( sts > 0 )
if ( nct > ncplane ) &
& call errore('data_structure','too many sticks',1)
IF ( nct > ncplane ) &
& CALL errore('data_structure','too many sticks',1)
if ( ncts > ncplanes ) &
& call errore('data_structure','too many sticks',2)
IF ( ncts > ncplanes ) &
& CALL errore('data_structure','too many sticks',2)
if ( nct == 0 ) &
& call errore('data_structure','number of sticks 0', 1)
IF ( nct == 0 ) &
& CALL errore('data_structure','number of sticks 0', 1)
if ( ncts == 0 ) &
& call errore('data_structure','number smooth sticks 0', 1)
IF ( ncts == 0 ) &
& CALL errore('data_structure','number smooth sticks 0', 1)
!
! local pointers deallocated at the end
@ -191,8 +191,8 @@ subroutine data_structure( lgamma )
ELSE
ngm = ngp ( me_pool + 1 ) / 2
ngms = ngps( me_pool + 1 ) / 2
END IF
END IF
ENDIF
ENDIF
CALL fft_dlay_allocate( dfftp, nproc_pool, nrx1, nrx2 )
CALL fft_dlay_allocate( dffts, nproc_pool, nrx1s, nrx2s )
@ -210,10 +210,10 @@ subroutine data_structure( lgamma )
! if tk = .FALSE. only half reciprocal space is considered, then we
! need to correct the number of sticks
IF( .NOT. tk ) THEN
IF( .not. tk ) THEN
nct = nct*2 - 1
ncts = ncts*2 - 1
END IF
ENDIF
!
! set the number of plane per process
@ -222,32 +222,32 @@ subroutine data_structure( lgamma )
! npp ( 1 : nproc_pool ) = dfftp%npp ( 1 : nproc_pool )
! npps( 1 : nproc_pool ) = dffts%npp ( 1 : nproc_pool )
if ( dfftp%nnp /= ncplane ) &
& call errore('data_structure','inconsistent plane dimension on dense grid', ABS(dfftp%nnp-ncplane) )
IF ( dfftp%nnp /= ncplane ) &
& CALL errore('data_structure','inconsistent plane dimension on dense grid', abs(dfftp%nnp-ncplane) )
if ( dffts%nnp /= ncplanes ) &
& call errore('data_structure','inconsistent plane dimension on smooth grid', ABS(dffts%nnp-ncplanes) )
IF ( dffts%nnp /= ncplanes ) &
& CALL errore('data_structure','inconsistent plane dimension on smooth grid', abs(dffts%nnp-ncplanes) )
WRITE( stdout, '(/5x,"Planes per process (thick) : nr3 =", &
& i4," npp = ",i4," ncplane =",i6)') nr3, dfftp%npp (me_pool + 1) , ncplane
if ( nr3s /= nr3 ) WRITE( stdout, '(5x,"Planes per process (smooth): nr3s=",&
IF ( nr3s /= nr3 ) WRITE( stdout, '(5x,"Planes per process (smooth): nr3s=",&
&i4," npps= ",i4," ncplanes=",i6)') nr3s, dffts%npp (me_pool + 1) , ncplanes
WRITE( stdout,*)
WRITE( stdout,'(5X, &
& "Proc/ planes cols G planes cols G columns G"/5X, &
& "Pool (dense grid) (smooth grid) (wavefct grid)")')
do i=1,nproc_pool
DO i=1,nproc_pool
WRITE( stdout,'(5x,i4,1x,2(i5,i7,i9),i7,i9)') i, dfftp%npp(i), ncp(i), ngp(i), &
& dffts%npp(i), ncps(i), ngps(i), nkcp(i), ngkp(i)
end do
ENDDO
IF ( nproc_pool > 1 ) THEN
WRITE( stdout,'(5x,"tot",2x,2(i5,i7,i9),i7,i9)') &
SUM(dfftp%npp(1:nproc_pool)), SUM(ncp(1:nproc_pool)), SUM(ngp(1:nproc_pool)), &
SUM(dffts%npp(1:nproc_pool)),SUM(ncps(1:nproc_pool)),SUM(ngps(1:nproc_pool)),&
SUM(nkcp(1:nproc_pool)), SUM(ngkp(1:nproc_pool))
END IF
sum(dfftp%npp(1:nproc_pool)), sum(ncp(1:nproc_pool)), sum(ngp(1:nproc_pool)), &
sum(dffts%npp(1:nproc_pool)),sum(ncps(1:nproc_pool)),sum(ngps(1:nproc_pool)),&
sum(nkcp(1:nproc_pool)), sum(ngkp(1:nproc_pool))
ENDIF
WRITE( stdout,*)
DEALLOCATE( stw, st, sts, in1, in2, idx, ngc, ngcs, ngkc )
@ -275,9 +275,9 @@ subroutine data_structure( lgamma )
nrxx = dfftp%nnr
nrxxs = dffts%nnr
nrxx = MAX( nrxx, nrxxs )
nrxx = max( nrxx, nrxxs )
!
! nxx is just a copy
! nxx is just a copy
!
nxx = nrxx
nxxs = nrxxs
@ -298,13 +298,13 @@ subroutine data_structure( lgamma )
nrx3s = nr3s
nrxxs = nrx1s * nrx2s * nrx3s
! nxx is just a copy
! nxx is just a copy
!
nxx = nrxx
nxxs = nrxxs
CALL fft_dlay_allocate( dfftp, nproc_pool, MAX(nrx1, nrx3), nrx2 )
CALL fft_dlay_allocate( dffts, nproc_pool, MAX(nrx1s, nrx3s), nrx2s )
CALL fft_dlay_allocate( dfftp, nproc_pool, max(nrx1, nrx3), nrx2 )
CALL fft_dlay_allocate( dffts, nproc_pool, max(nrx1s, nrx3s), nrx2s )
CALL calculate_gkcut()
@ -324,51 +324,51 @@ subroutine data_structure( lgamma )
ALLOCATE( stw ( lb(2):ub(2), lb(3):ub(3) ) )
stw = 0
do i1 = - n1, n1
DO i1 = - n1, n1
!
! Gamma-only: exclude space with x<0
!
if (lgamma .and. i1 < 0) go to 10
IF (lgamma .and. i1 < 0) GOTO 10
!
do i2 = - n2, n2
DO i2 = - n2, n2
!
! Gamma-only: exclude plane with x=0, y<0
!
if(lgamma .and. i1 == 0.and. i2 < 0) go to 20
IF(lgamma .and. i1 == 0.and. i2 < 0) GOTO 20
!
do i3 = - n3, n3
DO i3 = - n3, n3
!
! Gamma-only: exclude line with x=0, y=0, z<0
!
if(lgamma .and. i1 == 0 .and. i2 == 0 .and. i3 < 0) go to 30
IF(lgamma .and. i1 == 0 .and. i2 == 0 .and. i3 < 0) GOTO 30
!
amod = (i1 * bg (1, 1) + i2 * bg (1, 2) + i3 * bg (1, 3) ) **2 + &
(i1 * bg (2, 1) + i2 * bg (2, 2) + i3 * bg (2, 3) ) **2 + &
(i1 * bg (3, 1) + i2 * bg (3, 2) + i3 * bg (3, 3) ) **2
if (amod <= gcutm) ngm = ngm + 1
if (amod <= gcutms) ngms = ngms + 1
if (amod <= gkcut ) then
IF (amod <= gcutm) ngm = ngm + 1
IF (amod <= gcutms) ngms = ngms + 1
IF (amod <= gkcut ) THEN
stw( i2, i3 ) = 1
if (lgamma) stw( -i2, -i3 ) = 1
end if
30 continue
enddo
20 continue
enddo
10 continue
enddo
IF (lgamma) stw( -i2, -i3 ) = 1
ENDIF
30 CONTINUE
ENDDO
20 CONTINUE
ENDDO
10 CONTINUE
ENDDO
call fft_dlay_scalar( dfftp, ub, lb, nr1, nr2, nr3, nrx1, nrx2, nrx3, stw )
call fft_dlay_scalar( dffts, ub, lb, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, stw )
CALL fft_dlay_scalar( dfftp, ub, lb, nr1, nr2, nr3, nrx1, nrx2, nrx3, stw )
CALL fft_dlay_scalar( dffts, ub, lb, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, stw )
deallocate( stw )
DEALLOCATE( stw )
#endif
IF( nxx < dfftp%nnr ) &
CALL errore( ' data_structure ', ' inconsistent value for nxx ', ABS( nxx - dfftp%nnr ) )
CALL errore( ' data_structure ', ' inconsistent value for nxx ', abs( nxx - dfftp%nnr ) )
IF( nxxs /= dffts%nnr ) &
CALL errore( ' data_structure ', ' inconsistent value for nxxs ', ABS( nxxs - dffts%nnr ) )
CALL errore( ' data_structure ', ' inconsistent value for nxxs ', abs( nxxs - dffts%nnr ) )
!
! compute the global number of g, i.e. the sum over all processors
! within a pool
@ -377,8 +377,8 @@ subroutine data_structure( lgamma )
ngms_l = ngms
ngm_g = ngm
ngms_g = ngms
call mp_sum( ngm_g , intra_pool_comm )
call mp_sum( ngms_g, intra_pool_comm )
CALL mp_sum( ngm_g , intra_pool_comm )
CALL mp_sum( ngms_g, intra_pool_comm )
IF( use_task_groups ) THEN
!
@ -387,18 +387,18 @@ subroutine data_structure( lgamma )
!
CALL task_groups_init( dffts )
!
END IF
ENDIF
return
RETURN
contains
CONTAINS
subroutine calculate_gkcut()
SUBROUTINE calculate_gkcut()
!
integer :: kpoint
INTEGER :: kpoint
!
if (nks == 0) then
IF (nks == 0) THEN
!
! if k-points are automatically generated (which happens later)
! use max(bg)/2 as an estimate of the largest k-point
@ -407,22 +407,22 @@ contains
sqrt (bg (1, 1) **2 + bg (2, 1) **2 + bg (3, 1) **2), &
sqrt (bg (1, 2) **2 + bg (2, 2) **2 + bg (3, 2) **2), &
sqrt (bg (1, 3) **2 + bg (2, 3) **2 + bg (3, 3) **2) )
else
ELSE
gkcut = 0.0d0
do kpoint = 1, nks
DO kpoint = 1, nks
gkcut = max (gkcut, sqrt (ecutwfc) / tpiba + sqrt ( &
xk (1, kpoint) **2 + xk (2, kpoint) **2 + xk (3, kpoint) **2) )
enddo
endif
gkcut = gkcut**2
ENDDO
ENDIF
gkcut = gkcut**2
!
! ... find maximum value among all the processors
!
call mp_max (gkcut, inter_pool_comm )
CALL mp_max (gkcut, inter_pool_comm )
!
return
end subroutine calculate_gkcut
RETURN
END SUBROUTINE calculate_gkcut
end subroutine data_structure
END SUBROUTINE data_structure

78
PW/ggen.f90
View File

@ -82,17 +82,17 @@ SUBROUTINE ggen()
!
! Gamma-only: exclude space with x < 0
!
IF ( gamma_only .AND. i < 0) go to 10
IF ( gamma_only .and. i < 0) GOTO 10
DO j = - n2, n2
!
! exclude plane with x = 0, y < 0
!
IF ( gamma_only .AND. i == 0 .AND. j < 0) go to 11
IF ( gamma_only .and. i == 0 .and. j < 0) GOTO 11
DO k = - n3, n3
!
! exclude line with x = 0, y = 0, z < 0
!
IF ( gamma_only .AND. i == 0 .AND. j == 0 .AND. k < 0) go to 12
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)
@ -110,7 +110,7 @@ SUBROUTINE ggen()
ELSE
g2sort_g(ngm) = 0.d0
ENDIF
END IF
ENDIF
12 CONTINUE
ENDDO
11 CONTINUE
@ -119,9 +119,9 @@ SUBROUTINE ggen()
ENDDO
IF (ngm /= ngm_g ) &
CALL errore ('ggen', 'g-vectors missing !', ABS(ngm - ngm_g))
CALL errore ('ggen', 'g-vectors missing !', abs(ngm - ngm_g))
IF (ngms /= ngms_g) &
CALL errore ('ggen', 'smooth g-vectors missing !', ABS(ngms - ngms_g))
CALL errore ('ggen', 'smooth g-vectors missing !', abs(ngms - ngms_g))
igsrt(1) = 0
CALL hpsort_eps( ngm_g, g2sort_g, igsrt, eps8 )
@ -134,16 +134,16 @@ SUBROUTINE ggen()
iswap = mill_g(i,indsw)
mill_g(i,indsw) = mill_g(i,igsrt(indsw))
mill_g(i,igsrt(indsw)) = iswap
END DO
ENDDO
! .. swap indices
iswap = indsw; indsw = igsrt(indsw); igsrt(iswap) = iswap
IF(igsrt(indsw) == ng) THEN
igsrt(indsw)=indsw
ELSE
GOTO 7
END IF
END IF
END DO
ENDIF
ENDIF
ENDDO
DEALLOCATE( igsrt )
@ -159,12 +159,12 @@ SUBROUTINE ggen()
k = mill_g(3, ng)
#ifdef __PARA
m1 = MOD (i, nr1) + 1
IF (m1.LT.1) m1 = m1 + nr1
m2 = MOD (j, nr2) + 1
IF (m2.LT.1) m2 = m2 + nr2
m1 = mod (i, nr1) + 1
IF (m1.lt.1) m1 = m1 + nr1
m2 = mod (j, nr2) + 1
IF (m2.lt.1) m2 = m2 + nr2
mc = m1 + (m2 - 1) * nrx1
IF ( dfftp%isind ( mc ) .EQ.0) GOTO 1
IF ( dfftp%isind ( mc ) .eq.0) GOTO 1
#endif
tt = 0.d0
@ -185,7 +185,7 @@ SUBROUTINE ggen()
gg (ngm) = tt
IF (tt > eps8) THEN
esort (ngm) = tt
esort (ngm) = tt
ELSE
esort (ngm) = 0.d0
ENDIF
@ -193,8 +193,8 @@ SUBROUTINE ggen()
1 CONTINUE
ENDDO
IF (ngm.NE.ngmx) &
CALL errore ('ggen', 'g-vectors missing !', ABS(ngm - ngmx))
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.
@ -210,7 +210,7 @@ SUBROUTINE ggen()
!
DO ng = 1, ngm - 1
20 indsw = nl (ng)
IF (indsw.NE.ng) THEN
IF (indsw.ne.ng) THEN
DO ipol = 1, 3
swap = g (ipol, indsw)
g (ipol, indsw) = g (ipol, nl (indsw) )
@ -243,41 +243,41 @@ SUBROUTINE ggen()
!
! determine first nonzero g vector
!
IF (gg(1).LE.eps8) THEN
IF (gg(1).le.eps8) THEN
gstart=2
ELSE
gstart=1
END IF
ENDIF
!
! 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, &
n1 = nint (g (1, ng) * at (1, 1) + g (2, ng) * at (2, 1) + g (3, &
ng) * at (3, 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, &
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
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, &
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
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 (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) &
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) &
IF (ng.le.ngms) nls (ng) = n1s + (n2s - 1) * nrx1s + (n3s - 1) &
* nrx1s * nr2s
#endif
ELSE
@ -311,7 +311,7 @@ SUBROUTINE ggen()
igtongl (ng) = ngl
ENDDO
ALLOCATE (gl( ngl))
ALLOCATE (gl( ngl))
gl (1) = gg (1)
igl = 1
DO ng = 2, ngm
@ -321,7 +321,7 @@ SUBROUTINE ggen()
ENDIF
ENDDO
IF (igl.NE.ngl) CALL errore ('setup', 'igl <> ngl', ngl)
IF (igl.ne.ngl) CALL errore ('setup', 'igl <> ngl', ngl)
ENDIF
@ -360,21 +360,21 @@ SUBROUTINE index_minusg()
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 (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) &
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) &
IF (ng.le.ngms) nlsm(ng) = n1s + (n2s - 1) * nrx1s + (n3s - 1) &
* nrx1s * nr2s
#endif
ELSE
CALL errore('index_minusg','Mesh too small?',ng)
ENDIF
ENDDO
RETURN
END SUBROUTINE index_minusg

28
PW/gk_sort.f90
View File

@ -21,18 +21,18 @@ SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
!
! ... Here the dummy variables
!
INTEGER, INTENT(IN) :: ngm
INTEGER, INTENT(in) :: ngm
! input : the number of g vectors
INTEGER, INTENT(INOUT) :: ngk
INTEGER, INTENT(inout) :: ngk
! input/output : the number of k+G vectors inside the "ecut sphere"
INTEGER, INTENT(OUT) :: igk(npwx)
INTEGER, INTENT(out) :: igk(npwx)
! output : the correspondence k+G <-> G
REAL(DP), INTENT(IN) :: k(3), g(3,ngm), ecut
REAL(DP), INTENT(in) :: k(3), g(3,ngm), ecut
! input : the k point
! input : the coordinates of G vectors
! input : the cut-off energy
REAL(DP), INTENT(OUT) :: gk(npwx)
REAL(DP), INTENT(out) :: gk(npwx)
! output : the moduli of k+G
!
INTEGER :: ng, nk
@ -45,7 +45,7 @@ SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
!
! ... first we count the number of k+G vectors inside the cut-off sphere
!
q2x = ( SQRT( k(1)**2 + k(2)**2 + k(3)**2 ) + SQRT( ecut ) )**2
q2x = ( sqrt( k(1)**2 + k(2)**2 + k(3)**2 ) + sqrt( ecut ) )**2
!
ngk = 0
!
@ -64,15 +64,15 @@ SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
IF ( ngk > npwx ) &
CALL errore( 'gk_sort', 'array gk out-of-bounds', 1 )
!
IF ( q > eps8 ) THEN
IF ( q > eps8 ) THEN
!
gk(ngk) = q
gk(ngk) = q
!
ELSE
!
gk(ngk) = 0.D0
!
END IF
ENDIF
!
! ... set the initial value of index array
!
@ -82,11 +82,11 @@ SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
!
! ... if |G| > |k| + SQRT( Ecut ) stop search and order vectors
!
IF ( ( g(1,ng)**2 + g(2,ng)**2 + g(3,ng)**2 ) > ( q2x + eps8 ) ) EXIT
IF ( ( g(1,ng)**2 + g(2,ng)**2 + g(3,ng)**2 ) > ( q2x + eps8 ) ) exit
!
END IF
ENDIF
!
END DO
ENDDO
!
IF ( ng > ngm ) &
CALL infomsg( 'gk_sort', 'unexpected exit from do-loop')
@ -102,8 +102,8 @@ SUBROUTINE gk_sort( k, ngm, g, ecut, ngk, igk, gk )
gk(nk) = ( k(1) + g(1,igk(nk) ) )**2 + &
( k(2) + g(2,igk(nk) ) )**2 + &
( k(3) + g(3,igk(nk) ) )**2
!
END DO
!
ENDDO
!
RETURN
!

92
PW/pwcom.f90
View File

@ -15,7 +15,7 @@ MODULE basis
!
INTEGER :: &
natomwfc ! number of starting wavefunctions
CHARACTER(LEN=30) :: &!
CHARACTER(len=30) :: &!
starting_wfc, &! 'random' or 'atomic' or 'atomic+randm' or 'file'
starting_pot, &! 'atomic' or 'file'
startingconfig ! 'input' or 'file'
@ -67,7 +67,7 @@ MODULE gvect
ecfixed, &!
qcutz = 0.0_DP,&! For the modified Ekin functional
q2sigma !
complex(DP), ALLOCATABLE :: &
COMPLEX(DP), ALLOCATABLE :: &
eigts1(:,:), &!
eigts2(:,:), &! the phases e^{-iG*tau_s}
eigts3(:,:) !
@ -90,13 +90,13 @@ MODULE gsmooth
!
INTEGER :: &
ngms, &! the number of smooth G vectors
ngms_g, &! the global number of smooth G vectors
ngms_g, &! the global number of smooth G vectors
! (sum over all processors)
ngms_l, &! the local number of smooth G vectors
ngms_l, &! the local number of smooth G vectors
! (only present processor)
nr1s, &!
nr2s, &! the dimension of the smooth grid
nr3s, &!
nr3s, &!
nrx1s, &! maximum dimension of the smooth grid
nrx2s, &! maximum dimension of the smooth grid
nrx3s, &! maximum dimension of the smooth grid
@ -115,7 +115,7 @@ END MODULE gsmooth
MODULE klist
!
! ... The variables for the k-points
!
!
USE kinds, ONLY : DP
USE parameters, ONLY : npk
!
@ -130,8 +130,8 @@ MODULE klist
neldw, &! number of spin-dw electrons (if two_fermi_energies=t)
tot_magnetization, &! nelup-neldw >= 0 (negative value means unspecified)
tot_charge
REAL(DP) :: &
qnorm= 0.0_dp ! |q|, used in phonon+US calculations only
REAL(DP) :: &
qnorm= 0.0_dp ! |q|, used in phonon+US calculations only
INTEGER, ALLOCATABLE :: &
ngk(:) ! number of plane waves for each k point
INTEGER :: &
@ -140,8 +140,8 @@ MODULE klist
ngauss ! type of smearing technique
LOGICAL :: &
lgauss, &! if .TRUE.: use gaussian broadening
lxkcry=.FALSE., &! if .TRUE.:k-pnts in cryst. basis accepted in input
two_fermi_energies ! if .TRUE.: nelup and neldw set ef_up and ef_dw
lxkcry=.false., &! if .TRUE.:k-pnts in cryst. basis accepted in input
two_fermi_energies ! if .TRUE.: nelup and neldw set ef_up and ef_dw
! separately
!
END MODULE klist
@ -150,7 +150,7 @@ END MODULE klist
MODULE lsda_mod
!
! ... The variables needed for the lsda calculation
!
!
USE kinds, ONLY : DP
USE parameters, ONLY : ntypx, npk
!
@ -160,20 +160,20 @@ MODULE lsda_mod
lsda
REAL(DP) :: &
magtot, &! total magnetization
absmag, &! total absolute magnetization
absmag, &! total absolute magnetization
starting_magnetization(ntypx) ! the magnetization used to start with
INTEGER :: &
nspin, &! number of spin polarization: 2 if lsda, 1 other
current_spin, &! spin of the current kpoint
isk(npk) ! for each k-point: 1=spin up, 2=spin down
!
!
END MODULE lsda_mod
!
!
MODULE ktetra
!
! ... The variables for the tetrahedron method
!
!
SAVE
!
INTEGER :: &
@ -198,16 +198,16 @@ MODULE rap_point_group
nclass, & ! The number of classes of the point group
nelem(12), & ! The elements of each class
elem(8,12), & ! Which elements in the smat list for each class
which_irr(12) ! For each class gives its position in the
which_irr(12) ! For each class gives its position in the
! character table.
!
COMPLEX(DP) :: char_mat(12,12) ! the character tables
CHARACTER(LEN=15) :: name_rap(12) ! the name of the representation
CHARACTER(LEN=3) :: ir_ram(12) ! a string I, R or I+R for infrared,
CHARACTER(len=15) :: name_rap(12) ! the name of the representation
CHARACTER(len=3) :: ir_ram(12) ! a string I, R or I+R for infrared,
! Raman, or infrared+raman modes.
CHARACTER(LEN=11) :: gname ! the name of the group
CHARACTER(LEN=5) :: name_class(12) ! the name of the class
CHARACTER(len=11) :: gname ! the name of the group
CHARACTER(len=5) :: name_class(12) ! the name of the class
!
END MODULE rap_point_group
@ -220,14 +220,14 @@ MODULE rap_point_group_so
nelem_so(24), &! The elements of each class
elem_so(12,24), &! Which elements in the smat list for each class
has_e(12,24), & ! if -1 the smat is multiplied by -E
which_irr_so(24) ! For each class gives its position in the
which_irr_so(24) ! For each class gives its position in the
! character table.
!
COMPLEX(DP) :: char_mat_so(12,24), & ! the character tables
d_spin(2,2,48) ! the rotation in spin space
CHARACTER(LEN=15) :: name_rap_so(12) ! the name of the representation
CHARACTER(LEN=5) :: name_class_so(24), & ! the name of the class
CHARACTER(len=15) :: name_rap_so(12) ! the name of the representation
CHARACTER(len=5) :: name_class_so(24), & ! the name of the class
name_class_so1(24) ! the name of the class
!
END MODULE rap_point_group_so
@ -247,15 +247,15 @@ MODULE rap_point_group_is
COMPLEX(DP) :: &
d_spin_is(2,2,48) ! the rotation in spin space
CHARACTER(LEN=45) :: sname_is(48) ! name of the symmetries
CHARACTER(LEN=11) :: gname_is ! the name of the invariant group
CHARACTER(len=45) :: sname_is(48) ! name of the symmetries
CHARACTER(len=11) :: gname_is ! the name of the invariant group
!
END MODULE rap_point_group_is
!
MODULE vlocal
!
! ... The variables needed for the local potential in reciprocal space
!
!
USE kinds, ONLY : DP
!
SAVE
@ -271,7 +271,7 @@ END MODULE vlocal
MODULE wvfct
!
! ... The variables needed to compute the band structure
!
!
USE kinds, ONLY : DP
!
SAVE
@ -298,7 +298,7 @@ END MODULE wvfct
MODULE ener
!
! ... The variables needed to compute the energies
!
!
USE kinds, ONLY : DP
!
SAVE
@ -324,7 +324,7 @@ END MODULE ener
MODULE force_mod
!
! ... The variables for the first derivative of the energy
!
!
USE kinds, ONLY : DP
!
SAVE
@ -342,7 +342,7 @@ END MODULE force_mod
MODULE relax
!
! ... The variables used to control ionic relaxations
!
!
USE kinds, ONLY : DP
!
SAVE
@ -359,7 +359,7 @@ END MODULE relax
MODULE cellmd
!
! ... The variables used to control cell relaxation
!
!
USE kinds, ONLY : DP
!
SAVE
@ -376,7 +376,7 @@ MODULE cellmd
ntcheck ! # of steps between thermalizations
LOGICAL :: lmovecell ! used in cell relaxation
!
CHARACTER(LEN=2) :: &
CHARACTER(len=2) :: &
calc=' ' ! main switch for variable cell shape MD
! see readin, vcsmd and/or INPUT files
!
@ -386,7 +386,7 @@ END MODULE cellmd
MODULE us
!
! ... These parameters are needed with the US pseudopotentials
!
!
USE kinds, ONLY : DP
!
SAVE
@ -400,7 +400,7 @@ MODULE us
qrad(:,:,:,:), &! radial FT of Q functions
tab(:,:,:), &! interpolation table for PPs
tab_at(:,:,:) ! interpolation table for atomic wfc
LOGICAL :: spline_ps = .FALSE.
LOGICAL :: spline_ps = .false.
REAL(DP), ALLOCATABLE :: &
tab_d2y(:,:,:) ! for cubic splines
!
@ -410,7 +410,7 @@ END MODULE us
MODULE ldaU
!
! ... The quantities needed in lda+U calculations
!
!
USE kinds, ONLY : DP
USE parameters, ONLY : lqmax, ntypx
!
@ -421,19 +421,19 @@ MODULE ldaU
swfcatom(:,:) ! orthogonalized atomic wfcs
! REAL(DP), ALLOCATABLE :: &
! v_hub(:,:,:,:) ! the hubbard contribution to the potential
REAL(DP) :: &
REAL(DP) :: &
eth, &! the Hubbard contribution to the energy
Hubbard_U(ntypx), &! the Hubbard U
Hubbard_alpha(ntypx), &! the Hubbard alpha (used to calculate U)
starting_ns(lqmax,nspinx,ntypx) !
INTEGER :: &
INTEGER :: &
niter_with_fixed_ns, &! no. of iterations with fixed ns
Hubbard_l(ntypx), &! the angular momentum of Hubbard states
Hubbard_lmax = 0 ! maximum angular momentum of Hubbard states
LOGICAL :: &
LOGICAL :: &
lda_plus_u, &! .TRUE. if lda+u calculation is performed
conv_ns ! .TRUE. if ns are converged
CHARACTER(LEN=30) :: & ! 'atomic', 'ortho-atomic', 'file'
CHARACTER(len=30) :: & ! 'atomic', 'ortho-atomic', 'file'
U_projection ! specifies how input coordinates are given
INTEGER, ALLOCATABLE :: &
oatwfc(:) ! offset of atomic wfcs used for projections
@ -444,7 +444,7 @@ END MODULE ldaU
MODULE extfield
!
! ... The quantities needed in calculations with external field
!
!
USE kinds, ONLY : DP
!
SAVE
@ -479,8 +479,8 @@ MODULE fixed_occ
f_inp(:,:) ! the occupations for each spin
LOGICAL :: &
tfixed_occ, & ! if .TRUE. the occupations are fixed.
one_atom_occupations ! if .TRUE. the occupations are decided
! according to the projections of the
one_atom_occupations ! if .TRUE. the occupations are decided
! according to the projections of the
! wavefunctions on the initial atomic
! wavefunctions (to be used only
! for an isolated atom)
@ -488,10 +488,10 @@ MODULE fixed_occ
END MODULE fixed_occ
MODULE spin_orb
USE kinds, ONLY: DP
USE parameters, ONLY : lmaxx
SAVE
LOGICAL :: &
@ -514,8 +514,8 @@ MODULE bp
SAVE
!
LOGICAL :: &
lberry =.FALSE., & ! if .TRUE. calculate polarization using Berry phase
lelfield=.FALSE. ! if .TRUE. finite electric field using Berry phase
lberry =.false., & ! if .TRUE. calculate polarization using Berry phase
lelfield=.false. ! if .TRUE. finite electric field using Berry phase
INTEGER :: &
gdir, &! G-vector for polarization calculation
nppstr, &! number of k-points (parallel vector)
@ -525,7 +525,7 @@ MODULE bp
COMPLEX(DP), ALLOCATABLE , TARGET :: evcelm(:,:,:) ! wave function for storing projectors for electric field operator
COMPLEX(DP), ALLOCATABLE , TARGET :: evcelp(:,:,:) ! wave function for storing projectors for electric field operator
COMPLEX(DP), ALLOCATABLE, TARGET :: fact_hepsi(:,:)!factors for hermitean electric field operators
COMPLEX(DP), ALLOCATABLE, TARGET :: bec_evcel(:,:)!for storing bec's factors with evcel
COMPLEX(DP), ALLOCATABLE, TARGET :: bec_evcel(:,:)!for storing bec's factors with evcel
INTEGER, ALLOCATABLE, TARGET :: mapgp_global(:,:)! map for G'= G+1 correspondence
INTEGER, ALLOCATABLE, TARGET :: mapgm_global(:,:)! map for G'= G-1 correspondence
REAL(DP), ALLOCATABLE, TARGET :: forces_bp_efield(:,:)!ionic and US contributions to the atomic forces due to el. fields

1330
PW/realus.f90
View File

File diff suppressed because it is too large Load Diff

330
PW/vloc_psi.f90
View File

@ -24,7 +24,7 @@ SUBROUTINE vloc_psi_gamma(lda, n, m, psi, v, hpsi)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: lda, n, m
INTEGER, INTENT(in) :: lda, n, m
COMPLEX(DP), INTENT(in) :: psi (lda, m)
COMPLEX(DP), INTENT(inout):: hpsi (lda, m)
REAL(DP), INTENT(in) :: v(nrxxs)
@ -43,142 +43,142 @@ SUBROUTINE vloc_psi_gamma(lda, n, m, psi, v, hpsi)
!
use_tg = ( use_task_groups ) .and. ( m >= nogrp )
!
if( use_tg ) then
IF( use_tg ) THEN
!
v_siz = dffts%nnrx * nogrp
!
allocate( tg_v ( v_siz ) )
allocate( tg_psic( v_siz ) )
ALLOCATE( tg_v ( v_siz ) )
ALLOCATE( tg_psic( v_siz ) )
!
call tg_gather( dffts, v, tg_v )
CALL tg_gather( dffts, v, tg_v )
!
incr = 2 * nogrp
!
endif
ENDIF
!
! the local potential V_Loc psi. First bring psi to real space
!
do ibnd = 1, m, incr
DO ibnd = 1, m, incr
!
if( use_tg ) then
IF( use_tg ) THEN
!
tg_psic = (0.d0, 0.d0)
ioff = 0
!
do idx = 1, 2*nogrp, 2
if( idx + ibnd - 1 < m ) then
do j = 1, n
DO idx = 1, 2*nogrp, 2
IF( idx + ibnd - 1 < m ) THEN
DO j = 1, n
tg_psic(nls (igk(j))+ioff) = psi(j,idx+ibnd-1) + &
(0.0d0,1.d0) * psi(j,idx+ibnd)
tg_psic(nlsm(igk(j))+ioff) = conjg( psi(j,idx+ibnd-1) - &
(0.0d0,1.d0) * psi(j,idx+ibnd) )
enddo
elseif( idx + ibnd - 1 == m ) then
do j = 1, n
ENDDO
ELSEIF( idx + ibnd - 1 == m ) THEN
DO j = 1, n
tg_psic(nls (igk(j))+ioff) = psi(j,idx+ibnd-1)
tg_psic(nlsm(igk(j))+ioff) = conjg( psi(j,idx+ibnd-1) )
enddo
endif
ENDDO
ENDIF
ioff = ioff + dffts%nnrx
enddo
ENDDO
!
else
ELSE
!
psic(:) = (0.d0, 0.d0)
if (ibnd < m) then
IF (ibnd < m) THEN
! two ffts at the same time
do j = 1, n
DO j = 1, n
psic (nls (igk(j))) = psi(j, ibnd) + (0.0d0,1.d0)*psi(j, ibnd+1)
psic (nlsm(igk(j))) = conjg(psi(j, ibnd) - (0.0d0,1.d0)*psi(j, ibnd+1))
enddo
else
do j = 1, n
ENDDO
ELSE
DO j = 1, n
psic (nls (igk(j))) = psi(j, ibnd)
psic (nlsm(igk(j))) = conjg(psi(j, ibnd))
enddo
endif
ENDDO
ENDIF
!
endif
ENDIF
!
! fft to real space
! product with the potential v on the smooth grid
! back to reciprocal space
!
if( use_tg ) then
IF( use_tg ) THEN
!
call tg_cft3s ( tg_psic, dffts, 2, use_tg )
CALL tg_cft3s ( tg_psic, dffts, 2, use_tg )
!
do j = 1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
DO j = 1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
tg_psic (j) = tg_psic (j) * tg_v(j)
enddo
ENDDO
!
call tg_cft3s ( tg_psic, dffts, -2, use_tg )
CALL tg_cft3s ( tg_psic, dffts, -2, use_tg )
!
else
ELSE
!
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
CALL cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
!
do j = 1, nrxxs
DO j = 1, nrxxs
psic (j) = psic (j) * v(j)
enddo
ENDDO
!
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)
CALL cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)
!
endif
ENDIF
!
! addition to the total product
!
if( use_tg ) then
IF( use_tg ) THEN
!
ioff = 0
!
do idx = 1, 2*nogrp, 2
DO idx = 1, 2*nogrp, 2
!
if( idx + ibnd - 1 < m ) then
do j = 1, n
IF( idx + ibnd - 1 < m ) THEN
DO j = 1, n
fp= ( tg_psic( nls(igk(j)) + ioff ) + tg_psic( nlsm(igk(j)) + ioff ) ) * 0.5d0
fm= ( tg_psic( nls(igk(j)) + ioff ) - tg_psic( nlsm(igk(j)) + ioff ) ) * 0.5d0
hpsi (j, ibnd+idx-1) = hpsi (j, ibnd+idx-1) + cmplx( dble(fp), aimag(fm),kind=DP)
hpsi (j, ibnd+idx ) = hpsi (j, ibnd+idx ) + cmplx(aimag(fp),- dble(fm),kind=DP)
enddo
elseif( idx + ibnd - 1 == m ) then
do j = 1, n
ENDDO
ELSEIF( idx + ibnd - 1 == m ) THEN
DO j = 1, n
hpsi (j, ibnd+idx-1) = hpsi (j, ibnd+idx-1) + tg_psic( nls(igk(j)) + ioff )
enddo
endif
ENDDO
ENDIF
!
ioff = ioff + dffts%nr3x * dffts%nsw( me_pool + 1 )
!
enddo
ENDDO
!
else
if (ibnd < m) then
ELSE
IF (ibnd < m) THEN
! two ffts at the same time
do j = 1, n
DO j = 1, n
fp = (psic (nls(igk(j))) + psic (nlsm(igk(j))))*0.5d0
fm = (psic (nls(igk(j))) - psic (nlsm(igk(j))))*0.5d0
hpsi (j, ibnd) = hpsi (j, ibnd) + cmplx( dble(fp), aimag(fm),kind=DP)
hpsi (j, ibnd+1) = hpsi (j, ibnd+1) + cmplx(aimag(fp),- dble(fm),kind=DP)
enddo
else
do j = 1, n
ENDDO
ELSE
DO j = 1, n
hpsi (j, ibnd) = hpsi (j, ibnd) + psic (nls(igk(j)))
enddo
endif
endif
ENDDO
ENDIF
ENDIF
!
enddo
ENDDO
!
if( use_tg ) then
IF( use_tg ) THEN
!
deallocate( tg_psic )
deallocate( tg_v )
DEALLOCATE( tg_psic )
DEALLOCATE( tg_v )
!
endif
ENDIF
!
return
RETURN
END SUBROUTINE vloc_psi_gamma
!
!-----------------------------------------------------------------------
@ -199,7 +199,7 @@ SUBROUTINE vloc_psi_k(lda, n, m, psi, v, hpsi)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: lda, n, m
INTEGER, INTENT(in) :: lda, n, m
COMPLEX(DP), INTENT(in) :: psi (lda, m)
COMPLEX(DP), INTENT(inout):: hpsi (lda, m)
REAL(DP), INTENT(in) :: v(nrxxs)
@ -217,116 +217,116 @@ SUBROUTINE vloc_psi_k(lda, n, m, psi, v, hpsi)
!
incr = 1
!
if( use_tg ) then
IF( use_tg ) THEN
!
v_siz = dffts%nnrx * nogrp
!
allocate( tg_v ( v_siz ) )
allocate( tg_psic( v_siz ) )
ALLOCATE( tg_v ( v_siz ) )
ALLOCATE( tg_psic( v_siz ) )
!
call tg_gather( dffts, v, tg_v )
CALL tg_gather( dffts, v, tg_v )
incr = nogrp
!
endif
ENDIF
!
! the local potential V_Loc psi. First bring psi to real space
!
do ibnd = 1, m, incr
DO ibnd = 1, m, incr
!
if( use_tg ) then
IF( use_tg ) THEN
!
tg_psic = (0.d0, 0.d0)
ioff = 0
!
do idx = 1, nogrp
DO idx = 1, nogrp
if( idx + ibnd - 1 <= m ) then
IF( idx + ibnd - 1 <= m ) THEN
!$omp parallel do
do j = 1, n
DO j = 1, n
tg_psic(nls (igk(j))+ioff) = psi(j,idx+ibnd-1)
enddo
ENDDO
!$omp end parallel do
endif
ENDIF
ioff = ioff + dffts%nnrx
enddo
ENDDO
!
call tg_cft3s ( tg_psic, dffts, 2, use_tg )
CALL tg_cft3s ( tg_psic, dffts, 2, use_tg )
!
else
ELSE
!
psic(:) = (0.d0, 0.d0)
psic (nls (igk(1:n))) = psi(1:n, ibnd)
!
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
CALL cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
!
endif
ENDIF
!
! fft to real space
! product with the potential v on the smooth grid
! back to reciprocal space
!
if( use_tg ) then
IF( use_tg ) THEN
!
!$omp parallel do
do j = 1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
DO j = 1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
tg_psic (j) = tg_psic (j) * tg_v(j)
enddo
ENDDO
!$omp end parallel do
!
call tg_cft3s ( tg_psic, dffts, -2, use_tg )
CALL tg_cft3s ( tg_psic, dffts, -2, use_tg )
!
else
ELSE
!
!$omp parallel do
do j = 1, nrxxs
DO j = 1, nrxxs
psic (j) = psic (j) * v(j)
enddo
ENDDO
!$omp end parallel do
!
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)
CALL cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)
!
endif
ENDIF
!
! addition to the total product
!
if( use_tg ) then
IF( use_tg ) THEN
!
ioff = 0
!
do idx = 1, nogrp
DO idx = 1, nogrp
!
if( idx + ibnd - 1 <= m ) then
IF( idx + ibnd - 1 <= m ) THEN
!$omp parallel do
do j = 1, n
DO j = 1, n
hpsi (j, ibnd+idx-1) = hpsi (j, ibnd+idx-1) + tg_psic( nls(igk(j)) + ioff )
enddo
ENDDO
!$omp end parallel do
endif
ENDIF
!
ioff = ioff + dffts%nr3x * dffts%nsw( me_pool + 1 )
!
enddo
ENDDO
!
else
ELSE
!$omp parallel do
do j = 1, n
DO j = 1, n
hpsi (j, ibnd) = hpsi (j, ibnd) + psic (nls(igk(j)))
enddo
ENDDO
!$omp end parallel do
endif
ENDIF
!
enddo
ENDDO
!
if( use_tg ) then
IF( use_tg ) THEN
!
deallocate( tg_psic )
deallocate( tg_v )
DEALLOCATE( tg_psic )
DEALLOCATE( tg_v )
!
endif
ENDIF
!
return
RETURN
END SUBROUTINE vloc_psi_k
!
!-----------------------------------------------------------------------
@ -349,7 +349,7 @@ SUBROUTINE vloc_psi_nc (lda, n, m, psi, v, hpsi)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: lda, n, m
INTEGER, INTENT(in) :: lda, n, m
REAL(DP), INTENT(in) :: v(nrxx,4)
COMPLEX(DP), INTENT(in) :: psi (lda*npol, m)
COMPLEX(DP), INTENT(inout):: hpsi (lda,npol,m)
@ -368,125 +368,125 @@ SUBROUTINE vloc_psi_nc (lda, n, m, psi, v, hpsi)
!
use_tg = ( use_task_groups ) .and. ( m >= nogrp )
!
if( use_tg ) then
IF( use_tg ) THEN
v_siz = dffts%nnrx * nogrp
allocate( tg_v( v_siz, 4 ) )
call tg_gather( dffts, v(:,1), tg_v(:,1) )
call tg_gather( dffts, v(:,2), tg_v(:,2) )
call tg_gather( dffts, v(:,3), tg_v(:,3) )
call tg_gather( dffts, v(:,4), tg_v(:,4) )
allocate( tg_psic( v_siz, npol ) )
ALLOCATE( tg_v( v_siz, 4 ) )
CALL tg_gather( dffts, v(:,1), tg_v(:,1) )
CALL tg_gather( dffts, v(:,2), tg_v(:,2) )
CALL tg_gather( dffts, v(:,3), tg_v(:,3) )
CALL tg_gather( dffts, v(:,4), tg_v(:,4) )
ALLOCATE( tg_psic( v_siz, npol ) )
incr = nogrp
endif
ENDIF
!
! the local potential V_Loc psi. First the psi in real space
!
do ibnd = 1, m, incr
DO ibnd = 1, m, incr
if( use_tg ) then
IF( use_tg ) THEN
!
do ipol = 1, npol
DO ipol = 1, npol
!
tg_psic(:,ipol) = ( 0.D0, 0.D0 )
ioff = 0
!
do idx = 1, nogrp
DO idx = 1, nogrp
!
if( idx + ibnd - 1 <= m ) then
do j = 1, n
IF( idx + ibnd - 1 <= m ) THEN
DO j = 1, n
tg_psic( nls( igk(j) ) + ioff, ipol ) = psi( j +(ipol-1)*lda, idx+ibnd-1 )
enddo
endif
ENDDO
ENDIF
ioff = ioff + dffts%nnrx
enddo
ENDDO
!
call tg_cft3s ( tg_psic(:,ipol), dffts, 2, use_tg )
CALL tg_cft3s ( tg_psic(:,ipol), dffts, 2, use_tg )
!
enddo
ENDDO
!
else
ELSE
psic_nc = (0.d0,0.d0)
do ipol=1,npol
do j = 1, n
DO ipol=1,npol
DO j = 1, n
psic_nc(nls(igk(j)),ipol) = psi(j+(ipol-1)*lda,ibnd)
enddo
call cft3s (psic_nc(1,ipol), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
enddo
endif
ENDDO
CALL cft3s (psic_nc(1,ipol), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
ENDDO
ENDIF
!
! product with the potential v = (vltot+vr) on the smooth grid
!
if( use_tg ) then
do j=1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
IF( use_tg ) THEN
DO j=1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
sup = tg_psic(j,1) * (tg_v(j,1)+tg_v(j,4)) + &
tg_psic(j,2) * (tg_v(j,2)-(0.d0,1.d0)*tg_v(j,3))
sdwn = tg_psic(j,2) * (tg_v(j,1)-tg_v(j,4)) + &
tg_psic(j,1) * (tg_v(j,2)+(0.d0,1.d0)*tg_v(j,3))
tg_psic(j,1)=sup
tg_psic(j,2)=sdwn
enddo
else
do j=1, nrxxs
ENDDO
ELSE
DO j=1, nrxxs
sup = psic_nc(j,1) * (v(j,1)+v(j,4)) + &
psic_nc(j,2) * (v(j,2)-(0.d0,1.d0)*v(j,3))
sdwn = psic_nc(j,2) * (v(j,1)-v(j,4)) + &
psic_nc(j,1) * (v(j,2)+(0.d0,1.d0)*v(j,3))
psic_nc(j,1)=sup
psic_nc(j,2)=sdwn
enddo
endif
ENDDO
ENDIF
!
! back to reciprocal space
!
if( use_tg ) then
IF( use_tg ) THEN
!
do ipol = 1, npol
DO ipol = 1, npol
call tg_cft3s ( tg_psic(:,ipol), dffts, -2, use_tg )
CALL tg_cft3s ( tg_psic(:,ipol), dffts, -2, use_tg )
!
ioff = 0
!
do idx = 1, nogrp
DO idx = 1, nogrp
!
if( idx + ibnd - 1 <= m ) then
do j = 1, n
IF( idx + ibnd - 1 <= m ) THEN
DO j = 1, n
hpsi (j, ipol, ibnd+idx-1) = hpsi (j, ipol, ibnd+idx-1) + tg_psic( nls(igk(j)) + ioff, ipol )
enddo
endif
ENDDO
ENDIF
!
ioff = ioff + dffts%nr3x * dffts%nsw( me_pool + 1 )
!
enddo
ENDDO
enddo
ENDDO
!
else
ELSE
do ipol=1,npol
call cft3s (psic_nc(1,ipol), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2)
enddo
DO ipol=1,npol
CALL cft3s (psic_nc(1,ipol), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2)
ENDDO
!
! addition to the total product
!
do ipol=1,npol
do j = 1, n
DO ipol=1,npol
DO j = 1, n
hpsi(j,ipol,ibnd) = hpsi(j,ipol,ibnd) + psic_nc(nls(igk(j)),ipol)
enddo
enddo
ENDDO
ENDDO
endif
ENDIF
enddo
ENDDO
if( use_tg ) then
IF( use_tg ) THEN
!
deallocate( tg_v )
deallocate( tg_psic )
DEALLOCATE( tg_v )
DEALLOCATE( tg_psic )
!
endif
ENDIF
!
return
RETURN
END SUBROUTINE vloc_psi_nc