Harmonization of the three h*psi routines (two were contain-ed, one

was not: now none is contain-ed)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@5409 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2009-02-13 16:11:32 +00:00
parent 77396eb935
commit 017abc0fd8
1 changed files with 266 additions and 273 deletions

View File

@ -22,9 +22,64 @@ SUBROUTINE h_psi( lda, n, m, psi, hpsi )
! ... hpsi H*psi
!
USE kinds, ONLY : DP
USE bp, ONLY : lelfield,l3dstring,gdir, efield, efield_cry
USE control_flags, ONLY : gamma_only
USE noncollin_module, ONLY: npol, noncolin
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: lda, n, m
COMPLEX(DP), INTENT(INOUT) :: psi(lda*npol,m)
COMPLEX(DP), INTENT(OUT) :: hpsi(lda*npol,m)
!
INTEGER :: ipol
!
CALL start_clock( 'h_psi' )
!
IF ( gamma_only ) THEN
!
CALL h_psi_gamma( lda, n, m, psi, hpsi )
!
ELSE IF ( noncolin ) THEN
!
CALL h_psi_nc( lda, n, m, psi, hpsi )
!
ELSE
!
CALL h_psi_k ( lda, n, m, psi, hpsi )
!
END IF
!
! ... electric enthalpy if required
!
IF ( lelfield ) THEN
!
IF ( .NOT.l3dstring ) THEN
CALL h_epsi_her_apply( lda, n, m, psi, hpsi,gdir, efield )
ELSE
DO ipol=1,3
CALL h_epsi_her_apply( lda, n, m, psi, hpsi,ipol,efield_cry(ipol) )
END DO
END IF
!
END IF
!
CALL stop_clock( 'h_psi' )
!
RETURN
!
END SUBROUTINE h_psi
!
!-----------------------------------------------------------------------
SUBROUTINE h_psi_gamma ( lda, n, m, psi, hpsi )
!-----------------------------------------------------------------------
!
! ... gamma version
!
USE kinds, ONLY : DP
USE becmod, ONLY : rbecp, calbec
USE uspp, ONLY : vkb, nkb
USE wvfct, ONLY : igk, g2kin
USE gsmooth, ONLY : nls, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, nrxxs
USE ldaU, ONLY : lda_plus_u
USE lsda_mod, ONLY : current_spin
USE scf, ONLY : vrs
@ -34,299 +89,237 @@ SUBROUTINE h_psi( lda, n, m, psi, hpsi )
USE funct, ONLY : exx_is_active
#endif
USE funct, ONLY : dft_is_meta
USE bp, ONLY : lelfield,l3dstring,gdir, efield, efield_cry
USE control_flags, ONLY : gamma_only
USE noncollin_module, ONLY: npol, noncolin
!
IMPLICIT NONE
!
! ... input/output arguments
INTEGER, INTENT(IN) :: lda, n, m
COMPLEX(DP), INTENT(INOUT) :: psi(lda,m)
COMPLEX(DP), INTENT(OUT) :: hpsi(lda,m)
!
INTEGER :: lda, n, m
COMPLEX(DP) :: psi(lda*npol,m)
COMPLEX(DP) :: hpsi(lda*npol,m)
INTEGER :: ipol
INTEGER :: ibnd, j
!
!
CALL start_clock( 'h_psi' )
!
IF ( gamma_only ) THEN
!
#if defined __PGI
CALL h_psi_gamma( psi )
#else
CALL h_psi_gamma( )
#endif
!
ELSE IF ( noncolin ) THEN
!
! ... Unlike h_psi_gamma and h_psi_k, this is an external routine
!
CALL h_psi_nc( lda, n, m, psi, hpsi )
!
ELSE
!
#if defined __PGI
CALL h_psi_k( psi )
#else
CALL h_psi_k( )
#endif
!
END IF
CALL start_clock( 'h_psi:init' )
!
! ... electric enthalpy if required
! ... Here we apply the kinetic energy (k+G)^2 psi
!
DO ibnd = 1, m
!
! ... set to zero the imaginary part of psi at G=0
! ... absolutely needed for numerical stability
!
IF ( gstart == 2 ) psi(1,ibnd) = CMPLX( DBLE( psi(1,ibnd) ), 0.D0 )
!
hpsi(1:n,ibnd) = g2kin(1:n) * psi(1:n,ibnd)
!
END DO
!
CALL stop_clock( 'h_psi:init' )
!
if (dft_is_meta()) call h_psi_meta (lda, n, m, psi, hpsi)
!
! ... Here we add the Hubbard potential times psi
!
IF ( lda_plus_u ) CALL vhpsi( lda, n, m, psi, hpsi )
!
! ... the local potential V_Loc psi
!
CALL vloc_psi( lda, n, m, psi, vrs(1,current_spin), hpsi )
!
! ... Here the product with the non local potential V_NL psi
!
IF ( nkb > 0 ) THEN
!
CALL calbec ( n, vkb, psi, rbecp, m )
!
CALL add_vuspsi( lda, n, m, psi, hpsi )
!
IF ( lelfield ) THEN
if(.not.l3dstring) then
CALL h_epsi_her_apply( lda, n, m, psi, hpsi,gdir, efield )
else
do ipol=1,3
CALL h_epsi_her_apply( lda, n, m, psi, hpsi,ipol,efield_cry(ipol) )
enddo
endif
END IF
!
CALL stop_clock( 'h_psi' )
#ifdef EXX
IF ( exx_is_active() ) CALL vexx( lda, n, m, psi, hpsi )
#endif
!
RETURN
!
CONTAINS
!
!-----------------------------------------------------------------------
#if defined __PGI
SUBROUTINE h_psi_gamma( psi )
#else
SUBROUTINE h_psi_gamma( )
#endif
!-----------------------------------------------------------------------
!
! ... gamma version
!
USE becmod, ONLY : rbecp, calbec
!
IMPLICIT NONE
!
#if defined __PGI
COMPLEX(DP) :: psi(:,:)
#endif
INTEGER :: ibnd, j
!
!
CALL start_clock( 'h_psi:init' )
!
! ... Here we apply the kinetic energy (k+G)^2 psi
!
DO ibnd = 1, m
!
! ... set to zero the imaginary part of psi at G=0
! ... absolutely needed for numerical stability
!
IF ( gstart == 2 ) psi(1,ibnd) = CMPLX( DBLE( psi(1,ibnd) ), 0.D0 )
!
hpsi(1:n,ibnd) = g2kin(1:n) * psi(1:n,ibnd)
!
END DO
!
CALL stop_clock( 'h_psi:init' )
if (dft_is_meta()) call h_psi_meta (lda, n, m, psi, hpsi)
!
! ... Here we add the Hubbard potential times psi
!
IF ( lda_plus_u ) CALL vhpsi( lda, n, m, psi, hpsi )
!
! ... the local potential V_Loc psi
!
CALL vloc_psi( lda, n, m, psi, vrs(1,current_spin), hpsi )
!
! ... Here the product with the non local potential V_NL psi
!
IF ( nkb > 0 ) THEN
!
CALL calbec ( n, vkb, psi, rbecp, m )
!
CALL add_vuspsi( lda, n, m, psi, hpsi )
!
END IF
!
END SUBROUTINE h_psi_gamma
!
!-----------------------------------------------------------------------
SUBROUTINE h_psi_k ( lda, n, m, psi, hpsi )
!-----------------------------------------------------------------------
!
! ... k-points version
!
USE wavefunctions_module, ONLY : psic
USE becmod, ONLY : becp, calbec
USE control_flags, ONLY : use_task_groups
USE task_groups, ONLY : tg_gather
USE fft_parallel, ONLY : tg_cft3s
USE fft_base, ONLY : dffts
USE mp_global, ONLY : nogrp, me_pool
USE kinds, ONLY : DP
USE uspp, ONLY : vkb, nkb
USE wvfct, ONLY : igk, g2kin
USE gsmooth, ONLY : nls, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, nrxxs
USE ldaU, ONLY : lda_plus_u
USE lsda_mod, ONLY : current_spin
USE scf, ONLY : vrs
#ifdef EXX
IF ( exx_is_active() ) CALL vexx( lda, n, m, psi, hpsi )
USE exx, ONLY : vexx
USE funct, ONLY : exx_is_active
#endif
!
RETURN
!
END SUBROUTINE h_psi_gamma
USE funct, ONLY : dft_is_meta
USE bp, ONLY : lelfield,l3dstring,gdir, efield, efield_cry
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: lda, n, m
COMPLEX(DP), INTENT(IN) :: psi(lda,m)
COMPLEX(DP), INTENT(OUT):: hpsi(lda,m)
!
INTEGER :: ibnd, j, incr, v_siz, ioff, idx
!
COMPLEX(DP), ALLOCATABLE :: tg_psi(:)
REAL(DP), ALLOCATABLE :: tg_vrs(:)
LOGICAL :: use_tg
!
CALL start_clock( 'h_psi:init' )
!
use_tg = ( use_task_groups ) .AND. ( m >= nogrp )
!
incr = 1
!
IF( use_tg ) THEN
v_siz = dffts%nnrx * nogrp
ALLOCATE( tg_vrs( v_siz ) )
ALLOCATE( tg_psi( v_siz ) )
CALL tg_gather( dffts, vrs(:,current_spin), tg_vrs )
incr = nogrp
END IF
!
! ... Here we apply the kinetic energy (k+G)^2 psi
!
DO ibnd = 1, m
!
hpsi(1:n,ibnd) = g2kin(1:n) * psi(1:n,ibnd)
!
!-----------------------------------------------------------------------
#if defined __PGI
SUBROUTINE h_psi_k( psi )
#else
SUBROUTINE h_psi_k( )
#endif
!-----------------------------------------------------------------------
!
! ... k-points version
!
USE wavefunctions_module, ONLY : psic
USE becmod, ONLY : becp, calbec
USE control_flags, ONLY : use_task_groups
USE task_groups, ONLY : tg_gather
USE fft_parallel, ONLY : tg_cft3s
USE fft_base, ONLY : dffts
USE mp_global, ONLY : nogrp, me_pool
END DO
!
CALL stop_clock( 'h_psi:init' )
!
if (dft_is_meta()) call h_psi_meta (lda, n, m, psi, hpsi)
!
! ... Here we add the Hubbard potential times psi
!
IF ( lda_plus_u ) CALL vhpsi( lda, n, m, psi, hpsi )
!
! ... the local potential V_Loc psi. First the psi in real space
!
DO ibnd = 1, m, incr
!
CALL start_clock( 'h_psi:firstfft' )
!
IF( use_tg ) THEN
!
tg_psi = ( 0.D0, 0.D0 )
ioff = 0
!
DO idx = 1, nogrp
!
IF( idx + ibnd - 1 <= m ) THEN
DO j = 1, n
tg_psi( nls( igk(j) ) + ioff ) = psi(j,idx+ibnd-1)
END DO
END IF
!
IMPLICIT NONE
!
#if defined __PGI
COMPLEX(DP) :: psi(:,:)
#endif
INTEGER :: ibnd, j, incr, v_siz, ioff, idx
! counters
COMPLEX(DP), ALLOCATABLE :: tg_psi(:)
REAL(DP), ALLOCATABLE :: tg_vrs(:)
LOGICAL :: use_tg
!
CALL start_clock( 'h_psi:init' )
!
use_tg = ( use_task_groups ) .AND. ( m >= nogrp )
!
incr = 1
!
IF( use_tg ) THEN
v_siz = dffts%nnrx * nogrp
ALLOCATE( tg_vrs( v_siz ) )
ALLOCATE( tg_psi( v_siz ) )
CALL tg_gather( dffts, vrs(:,current_spin), tg_vrs )
incr = nogrp
END IF
!
! ... Here we apply the kinetic energy (k+G)^2 psi
!
DO ibnd = 1, m
!
hpsi(1:n,ibnd) = g2kin(1:n) * psi(1:n,ibnd)
!
END DO
!
CALL stop_clock( 'h_psi:init' )
if (dft_is_meta()) call h_psi_meta (lda, n, m, psi, hpsi)
!
! ... Here we add the Hubbard potential times psi
!
IF ( lda_plus_u ) CALL vhpsi( lda, n, m, psi, hpsi )
!
! ... the local potential V_Loc psi. First the psi in real space
!
DO ibnd = 1, m, incr
!
CALL start_clock( 'h_psi:firstfft' )
!
IF( use_tg ) THEN
!
tg_psi = ( 0.D0, 0.D0 )
ioff = 0
!
DO idx = 1, nogrp
!
IF( idx + ibnd - 1 <= m ) THEN
DO j = 1, n
tg_psi( nls( igk(j) ) + ioff ) = psi(j,idx+ibnd-1)
END DO
END IF
ioff = ioff + dffts%nnrx
ioff = ioff + dffts%nnrx
END DO
!
call tg_cft3s ( tg_psi, dffts, 2, use_tg )
!
ELSE
!
psic(1:nrxxs) = ( 0.D0, 0.D0 )
!
psic(nls(igk(1:n))) = psi(1:n,ibnd)
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
!
ENDIF
!
CALL stop_clock( 'h_psi:firstfft' )
!
! ... product with the potential vrs = (vltot+vr) on the smooth grid
!
IF( use_tg ) THEN
do j = 1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
tg_psi (j) = tg_psi (j) * tg_vrs(j)
enddo
ELSE
psic(1:nrxxs) = psic(1:nrxxs) * vrs(1:nrxxs,current_spin)
END IF
!
! ... back to reciprocal space
!
CALL start_clock( 'h_psi:secondfft' )
!
IF( use_tg ) THEN
!
call tg_cft3s ( tg_psi, dffts, -2, use_tg )
!
ioff = 0
!
DO idx = 1, nogrp
!
IF( idx + ibnd - 1 <= m ) THEN
DO j = 1, n
hpsi (j, ibnd+idx-1) = hpsi (j, ibnd+idx-1) + tg_psi( nls(igk(j)) + ioff )
END DO
!
call tg_cft3s ( tg_psi, dffts, 2, use_tg )
!
ELSE
!
psic(1:nrxxs) = ( 0.D0, 0.D0 )
!
psic(nls(igk(1:n))) = psi(1:n,ibnd)
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
!
ENDIF
!
CALL stop_clock( 'h_psi:firstfft' )
!
! ... product with the potential vrs = (vltot+vr) on the smooth grid
!
IF( use_tg ) THEN
do j = 1, nrx1s * nrx2s * dffts%tg_npp( me_pool + 1 )
tg_psi (j) = tg_psi (j) * tg_vrs(j)
enddo
ELSE
psic(1:nrxxs) = psic(1:nrxxs) * vrs(1:nrxxs,current_spin)
END IF
!
! ... back to reciprocal space
!
CALL start_clock( 'h_psi:secondfft' )
!
IF( use_tg ) THEN
!
call tg_cft3s ( tg_psi, dffts, -2, use_tg )
!
ioff = 0
!
DO idx = 1, nogrp
!
IF( idx + ibnd - 1 <= m ) THEN
DO j = 1, n
hpsi (j, ibnd+idx-1) = hpsi (j, ibnd+idx-1) + tg_psi( nls(igk(j)) + ioff )
END DO
END IF
!
ioff = ioff + dffts%nr3x * dffts%nsw( me_pool + 1 )
!
END DO
!
ELSE
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2 )
!
! ... addition to the total product
!
hpsi(1:n,ibnd) = hpsi(1:n,ibnd) + psic(nls(igk(1:n)))
!
END IF
!
CALL stop_clock( 'h_psi:secondfft' )
!
END DO
!
IF( use_tg ) THEN
!
DEALLOCATE( tg_vrs )
DEALLOCATE( tg_psi )
!
END IF
!
! ... Here the product with the non local potential V_NL psi
!
IF ( nkb > 0 ) THEN
!
CALL calbec( n, vkb, psi, becp, m )
!
CALL add_vuspsi( lda, n, m, psi, hpsi )
!
END IF
!
#ifdef EXX
IF ( exx_is_active() ) CALL vexx( lda, n, m, psi, hpsi )
#endif
!
RETURN
!
END SUBROUTINE h_psi_k
END IF
!
ioff = ioff + dffts%nr3x * dffts%nsw( me_pool + 1 )
!
END DO
!
ELSE
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -2 )
!
! ... addition to the total product
!
hpsi(1:n,ibnd) = hpsi(1:n,ibnd) + psic(nls(igk(1:n)))
!
END IF
!
END SUBROUTINE h_psi
CALL stop_clock( 'h_psi:secondfft' )
!
END DO
!
IF( use_tg ) THEN
!
DEALLOCATE( tg_vrs )
DEALLOCATE( tg_psi )
!
END IF
!
! ... Here the product with the non local potential V_NL psi
!
IF ( nkb > 0 ) THEN
!
CALL calbec( n, vkb, psi, becp, m )
!
CALL add_vuspsi( lda, n, m, psi, hpsi )
!
END IF
!
#ifdef EXX
IF ( exx_is_active() ) CALL vexx( lda, n, m, psi, hpsi )
#endif
!
RETURN
!
END SUBROUTINE h_psi_k
!
!-----------------------------------------------------------------------
subroutine h_psi_nc (lda, n, m, psi, hpsi)