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

General Cleanup. Some problems related to the use of local pseudopotentials have been fixed. 

C.S.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@502 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sbraccia 2004-01-22 12:01:13 +00:00
parent 9c66ff989f
commit 0b32850e1f
6 changed files with 534 additions and 410 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
PW/c_bands.f90
View File

@ -149,7 +149,8 @@ SUBROUTINE c_bands( iter, ik_, dr2 )
!
! ... various initializations
!
CALL init_us_2( npw, igk, xk(1,ik), vkb )
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
!
! ... read in wavefunctions from the previous iteration
!
@ -355,7 +356,8 @@ SUBROUTINE c_bands( iter, ik_, dr2 )
!
! ... various initializations
!
CALL init_us_2( npw, igk, xk(1,ik), vkb )
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
!
! ... read in wavefunctions from the previous iteration
!

27
PW/force_us.f90
View File

@ -60,7 +60,7 @@ SUBROUTINE force_us( forcenl )
!
REAL(KIND=DP), ALLOCATABLE :: becp(:,:), dbecp (:,:,:)
! auxiliary variables contain <beta|psi> and <dbeta|psi>
COMPLEX(KIND=DP), ALLOCATABLE :: vkb1 (:,:)
COMPLEX(KIND=DP), ALLOCATABLE :: vkb1(:,:)
! auxiliary variable contains g*|beta>
REAL(KIND=DP) :: ps
INTEGER :: ik, ipol, ibnd, ig, ih, jh, na, nt, ikb, jkb, ijkb0
@ -82,20 +82,23 @@ SUBROUTINE force_us( forcenl )
IF ( nks > 1 ) THEN
READ( iunigk ) npw, igk
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
CALL init_us_2( npw, igk, xk(1,ik), vkb )
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
END IF
!
CALL pw_gemm( 'Y', nkb, nbnd, npw, vkb, npwx, evc, npwx, becp, nkb )
IF ( nkb > 0 ) &
CALL pw_gemm( 'Y', nkb, nbnd, npw, vkb, npwx, evc, npwx, becp, nkb )
!
DO ipol = 1, 3
DO jkb = 1,nkb
DO jkb = 1, nkb
DO ig = 1, npw
vkb1(ig,jkb) = vkb(ig,jkb) * (0.D0,-1.D0) * g(ipol,igk(ig))
END DO
END DO
!
CALL pw_gemm( 'Y', nkb, nbnd, npw, vkb1, npwx, evc, npwx, &
dbecp(1,1,ipol), nkb )
IF ( nkb > 0 ) &
CALL pw_gemm( 'Y', nkb, nbnd, npw, vkb1, npwx, evc, npwx, &
dbecp(1,1,ipol), nkb )
!
END DO
!
@ -213,20 +216,23 @@ SUBROUTINE force_us( forcenl )
IF ( nks > 1 ) THEN
READ( iunigk ) npw, igk
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
CALL init_us_2( npw, igk, xk(1,ik), vkb )
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
END IF
!
CALL ccalbec( nkb, npwx, npw, nbnd, becp, vkb, evc )
!
DO ipol = 1, 3
DO jkb = 1,nkb
DO jkb = 1, nkb
DO ig = 1, npw
vkb1(ig,jkb) = vkb(ig,jkb) * (0.D0,-1.D0) * g(ipol,igk(ig))
END DO
END DO
!
CALL ZGEMM( 'C', 'N', nkb, nbnd, npw, (1.D0, 0.D0), vkb1, npwx, &
evc, npwx, (0.D0, 0.D0), dbecp(1,1,ipol), nkb )
IF ( nkb > 0 ) &
CALL ZGEMM( 'C', 'N', nkb, nbnd, npw, ( 1.D0, 0.D0 ), &
vkb1, npwx, evc, npwx, ( 0.D0, 0.D0 ), &
dbecp(1,1,ipol), nkb )
END DO
!
ijkb0 = 0
@ -319,4 +325,3 @@ SUBROUTINE force_us( forcenl )
END SUBROUTINE force_us_k
!
END SUBROUTINE force_us

6
PW/h_psi.f90
View File

@ -99,7 +99,8 @@ SUBROUTINE h_psi( lda, n, m, psi, hpsi )
!
! ... Here the product with the non local potential V_NL psi
!
CALL pw_gemm( 'Y', nkb, m, n, vkb, lda, psi, lda, becp, nkb )
IF ( nkb > 0 ) &
CALL pw_gemm( 'Y', nkb, m, n, vkb, lda, psi, lda, becp, nkb )
!
IF ( nkb > 0 ) CALL add_vuspsi( lda, n, m, psi, hpsi )
!
@ -125,7 +126,8 @@ SUBROUTINE h_psi( lda, n, m, psi, hpsi )
!
CALL start_clock( 'init' )
!
CALL ccalbec( nkb, lda, n, m, becp, vkb, psi )
IF ( nkb > 0 ) &
CALL ccalbec( nkb, lda, n, m, becp, vkb, psi )
!
! ... Here we apply the kinetic energy (k+G)^2 psi
!

6
PW/pw_gemm.f90
View File

@ -39,13 +39,13 @@ SUBROUTINE pw_gemm( sum_over_nodes, na, nb, n, a, lda, b, ldb, c, ldc )
CALL start_clock( 'pw_gemm' )
!
CALL DGEMM( 'C', 'N', na, nb, 2 * n, 2.D0, a, 2 * lda, b, &
2 * ldb, 0.d0, c, ldc )
2 * ldb, 0.D0, c, ldc )
!
IF ( gstart == 2 ) &
CALL DGER( na, nb, -1.D0, a, 2 * lda, b, 2 * ldb, c, ldc )
CALL DGER( na, nb, -1.D0, a, ( 2 * lda ), b, ( 2 * ldb ), c, ldc )
!
#ifdef __PARA
if ( sum_over_nodes == 'y' .OR. sum_over_nodes == 'Y' ) &
IF ( sum_over_nodes == 'y' .OR. sum_over_nodes == 'Y' ) &
CALL reduce( ldc * nb, c )
#endif
!

893
PW/sum_band.f90
View File

@ -58,440 +58,553 @@ SUBROUTINE sum_band()
!
RETURN
!
CONTAINS
!
!--------------------------------------------------------------------------
SUBROUTINE sum_band_gamma()
!--------------------------------------------------------------------------
!
! ... gamma version
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: ikb, jkb, ijkb0, ih, jh, ijh, na, np
! counters on beta functions, atoms, pseudopotentials
INTEGER :: ir, is, ig, ibnd, ik
! counter on 3D r points
! counter on spin polarizations
! counter on g vectors
! counter on bands
! counter on k points
REAL(KIND=DP) :: w1, w2
! weight
REAL(KIND=DP), ALLOCATABLE :: becp(:,:)
! contain <beta|psi>
!
!
ALLOCATE( becp( nkb, nbnd ) )
!
becsum(:,:,:) = 0.D0
rho(:,:) = 0.D0
eband = 0.D0
demet = 0.D0
!
! ... calculate weights for the insulator case
!
IF ( .NOT. lgauss .AND. .NOT. ltetra .AND. .NOT. tfixed_occ ) THEN
CONTAINS
!
!--------------------------------------------------------------------------
SUBROUTINE sum_band_gamma()
!--------------------------------------------------------------------------
!
CALL iweights( nks, wk, nbnd, nelec, et, ef, wg )
! ... gamma version
!
! ... calculate weights for the metallic case
IMPLICIT NONE
!
ELSE IF ( ltetra ) THEN
! ... local variables
!
INTEGER :: ikb, jkb, ijkb0, ih, jh, ijh, na, np
! counters on beta functions, atoms, pseudopotentials
INTEGER :: ir, is, ig, ibnd, ik
! counter on 3D r points
! counter on spin polarizations
! counter on g vectors
! counter on bands
! counter on k points
REAL(KIND=DP) :: w1, w2
! weight
REAL(KIND=DP), ALLOCATABLE :: becp(:,:)
! contains <beta|psi>
!
!
ALLOCATE( becp( nkb, nbnd ) )
!
becsum(:,:,:) = 0.D0
rho(:,:) = 0.D0
eband = 0.D0
demet = 0.D0
!
! ... calculate weights for the insulator case
!
IF ( .NOT. lgauss .AND. .NOT. ltetra .AND. .NOT. tfixed_occ ) THEN
!
CALL iweights( nks, wk, nbnd, nelec, et, ef, wg )
!
! ... calculate weights for the metallic case
!
ELSE IF ( ltetra ) THEN
!
#ifdef __PARA
CALL poolrecover( et, nbnd, nkstot, nks )
!
IF ( me == 1 .AND. mypool == 1 ) THEN
CALL poolrecover( et, nbnd, nkstot, nks )
!
IF ( me == 1 .AND. mypool == 1 ) THEN
!
#endif
CALL tweights( nkstot, nspin, nbnd, nelec, ntetra, &
tetra, et, ef, wg )
#ifdef __PARA
!
END IF
!
CALL poolscatter( nbnd, nkstot, wg, nks, wg )
!
IF ( me == 1 ) CALL poolbcast( 1, ef )
!
CALL broadcast( 1, ef )
!
#endif
CALL tweights( nkstot, nspin, nbnd, nelec, ntetra, &
tetra, et, ef, wg )
#ifdef __PARA
ELSE IF ( lgauss ) THEN
!
CALL gweights( nks, wk, nbnd, nelec, &
degauss, ngauss, et, ef, demet, wg )
!
ELSE IF ( tfixed_occ ) THEN
!
ef = - 1.0D+20
!
DO is = 1, nspin
!
DO ibnd = 1, nbnd
!
wg(ibnd,is) = f_inp(ibnd,is)
!
IF ( wg(ibnd,is) > 0.D0 ) ef = MAX( ef, et(ibnd,is) )
!
END DO
!
END DO
!
END IF
!
CALL poolscatter( nbnd, nkstot, wg, nks, wg )
! ... Needed for LDA+U
!
IF ( me == 1 ) CALL poolbcast( 1, ef )
IF ( lda_plus_u ) CALL new_ns()
!
CALL broadcast( 1, ef )
! ... here we sum for each k point the contribution
! ... of the wavefunctions to the charge
!
#endif
ELSE IF ( lgauss ) THEN
CALL gweights( nks, wk, nbnd, nelec, degauss, ngauss, &
et, ef, demet, wg )
ELSE IF ( tfixed_occ ) THEN
ef = - 1.0e+20
DO is = 1, nspin
IF ( nks > 1 ) REWIND( iunigk )
!
k_loop: DO ik = 1, nks
!
IF ( lsda ) current_spin = isk(ik)
!
IF ( nks > 1 ) THEN
!
READ( iunigk ) npw, igk
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
!
END IF
!
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
!
! ... here we compute the band energy: the sum of the eigenvalues
!
DO ibnd = 1, nbnd
wg(ibnd,is) = f_inp(ibnd,is)
if (wg(ibnd,is) > 0.d0) ef = MAX (ef, et (ibnd, is) )
!
eband = eband + et(ibnd,ik) * wg(ibnd,ik)
!
! ... the sum of eband and demet is the integral for
! ... e < ef of e n(e) which reduces for degauss=0 to the sum of
! ... the eigenvalues.
!
END DO
END DO
END IF
!
! ... Needed for LDA+U
!
IF ( lda_plus_u ) CALL new_ns()
!
! ... here we sum for each k point the contribution
! ... of the wavefunctions to the charge
!
IF ( nks > 1 ) REWIND( iunigk )
!
DO ik = 1, nks
IF ( lsda ) current_spin = isk(ik)
IF ( nks > 1 ) THEN
READ( iunigk ) npw, igk
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
END IF
CALL init_us_2( npw, igk, xk(1,ik), vkb )
!
! ... here we compute the band energy: the sum of the eigenvalues
!
DO ibnd = 1, nbnd
!
eband = eband + et(ibnd,ik) * wg(ibnd,ik)
!
! ... the sum of eband and demet is the integral for e < ef of e n(e)
! ... which reduces for degauss=0 to the sum of the eigenvalues.
!
END DO
DO ibnd = 1, nbnd, 2
psic(:) = (0.D0,0.D0)
IF ( ibnd < nbnd ) THEN
DO ibnd = 1, nbnd, 2
!
! ... two ffts at the same time
psic(:) = ( 0.D0, 0.D0 )
!
DO ig = 1, npw
psic(nls(igk(ig))) = &
evc(ig,ibnd) + (0.D0,1.D0) * evc(ig,ibnd+1)
psic(nlsm(igk(ig))) = &
CONJG( evc(ig,ibnd) - (0.D0,1.D0) * evc(ig,ibnd+1) )
IF ( ibnd < nbnd ) THEN
!
! ... two ffts at the same time
!
psic(nls(igk(1:npw))) = evc(1:npw,ibnd) + &
( 0.D0, 1.D0 ) * evc(1:npw,ibnd+1)
psic(nlsm(igk(1:npw))) = CONJG( evc(1:npw,ibnd) - &
( 0.D0, 1.D0 ) * evc(1:npw,ibnd+1) )
!
ELSE
!
psic(nls(igk(1:npw))) = evc(1:npw,ibnd)
psic(nlsm(igk(1:npw))) = CONJG( evc(1:npw,ibnd) )
!
END IF
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
!
w1 = wg(ibnd,ik) / omega
!
! ... increment the charge density ...
!
IF ( ibnd < nbnd ) THEN
!
! ... two ffts at the same time
!
w2 = wg(ibnd+1,ik) / omega
!
ELSE
!
w2 = w1
!
END IF
!
DO ir = 1, nrxxs
!
rho(ir,current_spin) = rho(ir,current_spin) + &
w1 * REAL( psic(ir) )**2 + &
w2 * IMAG( psic(ir) )**2
!
END DO
ELSE
DO ig = 1, npw
psic(nls(igk(ig))) = evc(ig,ibnd)
psic(nlsm(igk(ig))) = CONJG( evc(ig,ibnd) )
END DO
END IF
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
!
w1 = wg(ibnd,ik) / omega
!
! ... increment the charge density ...
!
IF ( ibnd < nbnd ) THEN
!
! ... two ffts at the same time
!
w2 = wg(ibnd+1,ik) / omega
ELSE
w2 = w1
END IF
DO ir = 1, nrxxs
rho(ir,current_spin) = rho(ir,current_spin) + &
w1 * REAL( psic(ir) )**2 + w2 * IMAG( psic(ir) )**2
END DO
END DO
!
! ... If we have a US pseudopotential we compute here the sumbec term
!
IF ( .NOT. okvan ) GO TO 10
!
CALL pw_gemm( 'Y', nkb, nbnd, npw, vkb, npwx, evc, npwx, becp, nkb )
!
CALL start_clock( 'sumbec' )
!
DO ibnd = 1, nbnd
w1 = wg(ibnd,ik)
ijkb0 = 0
DO np = 1, ntyp
IF ( tvanp(np) ) THEN
DO na = 1, nat
IF ( ityp(na) == np ) THEN
ijh = 1
DO ih = 1, nh(np)
ikb = ijkb0 + ih
becsum(ijh,na,current_spin) = &
!
! ... If we have a US pseudopotential we compute here the sumbec term
!
IF ( .NOT. okvan ) CYCLE k_loop
!
IF ( nkb > 0 ) &
CALL pw_gemm( 'Y', nkb, nbnd, npw, &
vkb, npwx, evc, npwx, becp, nkb )
!
CALL start_clock( 'sumbec' )
!
DO ibnd = 1, nbnd
!
w1 = wg(ibnd,ik)
ijkb0 = 0
!
DO np = 1, ntyp
!
IF ( tvanp(np) ) THEN
!
DO na = 1, nat
!
IF ( ityp(na) == np ) THEN
!
ijh = 1
!
DO ih = 1, nh(np)
!
ikb = ijkb0 + ih
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
w1 * becp(ikb,ibnd) * becp(ikb,ibnd)
ijh = ijh + 1
DO jh = ( ih + 1 ), nh(np)
jkb = ijkb0 + jh
becsum(ijh,na,current_spin) = &
!
ijh = ijh + 1
!
DO jh = ( ih + 1 ), nh(np)
!
jkb = ijkb0 + jh
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
w1 * 2.D0 * becp(ikb,ibnd) * becp(jkb,ibnd)
ijh = ijh + 1
!
ijh = ijh + 1
!
END DO
!
END DO
END DO
ijkb0 = ijkb0 + nh(np)
END IF
END DO
ELSE
DO na = 1, nat
IF ( ityp(na) == np ) ijkb0 = ijkb0 + nh(np)
END DO
END IF
!
ijkb0 = ijkb0 + nh(np)
!
END IF
!
END DO
!
ELSE
!
DO na = 1, nat
!
IF ( ityp(na) == np ) ijkb0 = ijkb0 + nh(np)
!
END DO
!
END IF
!
END DO
!
END DO
END DO
!
CALL stop_clock( 'sumbec' )
!
10 CONTINUE
!
END DO
!
DEALLOCATE( becp )
!
! ... If a double grid is used, interpolate onto the fine grid
!
IF ( doublegrid ) THEN
DO is = 1, nspin
CALL interpolate( rho(1,is), rho(1,is), 1 )
END DO
END IF
!
! ... Here we add the Ultrasoft contribution to the charge
!
IF ( okvan ) CALL addusdens()
!
#ifdef __PARA
CALL poolreduce( 1, eband )
CALL poolreduce( 1, demet )
#endif
!
! ... symmetrization of the charge density (and local magnetization)
!
#ifdef __PARA
!
! ... reduce charge density across pools
!
CALL poolreduce( nspin * nrxx, rho )
!
DO is = 1, nspin
CALL psymrho( rho(1,is), nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, &
s, ftau )
END DO
#else
DO is = 1, nspin
CALL symrho( rho(1,is), nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, &
s, ftau )
END DO
#endif
!
RETURN
!
END SUBROUTINE sum_band_gamma
!
!
!--------------------------------------------------------------------
subroutine sum_band_k()
!--------------------------------------------------------------------
!
! ... k-points version
!
USE becmod, ONLY : becp
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: ikb, jkb, ijkb0, ih, jh, ijh, na, np
! counters on beta functions, atoms, pseudopotentials
INTEGER :: ir, is, ig, ibnd, ik
! counter on 3D r points
! counter on spin polarizations
! counter on g vectors
! counter on bands
! counter on k points
REAL(KIND=DP) :: w1
! weight
!
!
becsum(:,:,:) = 0.D0
rho(:,:) = 0.D0
eband = 0.D0
demet = 0.D0
!
! ... calculate weights for the insulator case
!
IF ( .NOT. lgauss .AND. .NOT. ltetra .AND. .NOT. tfixed_occ ) THEN
!
CALL iweights( nks, wk, nbnd, nelec, et, ef, wg )
!
! ... calculate weights for the metallic case
!
ELSE IF ( ltetra ) THEN
#ifdef __PARA
!
CALL poolrecover( et, nbnd, nkstot, nks )
!
IF ( me == 1 .AND. mypool == 1 ) THEN
!
#endif
CALL tweights( nkstot, nspin, nbnd, nelec, ntetra, tetra, et, ef, wg )
#ifdef __PARA
CALL stop_clock( 'sumbec' )
!
ENDIF
END DO k_loop
!
CALL poolscatter( nbnd, nkstot, wg, nks, wg )
DEALLOCATE( becp )
!
IF ( me == 1 ) CALL poolbcast( 1, ef )
! ... If a double grid is used, interpolate onto the fine grid
!
CALL broadcast( 1, ef )
#endif
!
ELSE IF ( lgauss ) THEN
!
CALL gweights( nks, wk, nbnd, nelec, degauss, ngauss, et, ef, demet, wg )
!
ELSE IF ( tfixed_occ ) THEN
!
ef = - 1.0e+20
DO is = 1, nspin
DO ibnd = 1, nbnd
wg(ibnd,is) = f_inp(ibnd,is)
if (wg(ibnd,is) > 0.d0) ef = MAX (ef, et (ibnd, is) )
IF ( doublegrid ) THEN
!
DO is = 1, nspin
!
CALL interpolate( rho(1,is), rho(1,is), 1 )
!
END DO
END DO
END IF
!
! ... Needed for LDA+U
!
IF ( lda_plus_u ) call new_ns()
!
! ... here we sum for each k point the contribution
! ... of the wavefunctions to the charge
!
IF ( nks > 1 ) REWIND( iunigk )
!
DO ik = 1, nks
IF ( lsda ) current_spin = isk(ik)
IF ( nks > 1 ) THEN
READ( iunigk ) npw, igk
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
!
END IF
!
CALL init_us_2( npw, igk, xk(1,ik), vkb )
! ... Here we add the Ultrasoft contribution to the charge
!
! ... here we compute the band energy: the sum of the eigenvalues
IF ( okvan ) CALL addusdens()
!
DO ibnd = 1, nbnd
!
eband = eband + et(ibnd,ik) * wg(ibnd,ik)
!
! ... the sum of eband and demet is the integral for e < ef of e n(e)
! ... which reduces for degauss=0 to the sum of the eigenvalues.
! ... the factors two is for spin degeneracy
!
psic(:) = (0.D0,0.D0)
!
DO ig = 1, npw
psic(nls(igk(ig))) = evc(ig,ibnd)
END DO
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
!
w1 = wg(ibnd,ik) / omega
!
! ... increment the charge density ...
!
DO ir = 1, nrxxs
rho(ir,current_spin) = rho(ir,current_spin) + &
w1 * ( REAL( psic(ir) )**2 + IMAG( psic(ir) )**2 )
END DO
END DO
#ifdef __PARA
CALL poolreduce( 1, eband )
CALL poolreduce( 1, demet )
#endif
!
! ... If we have a US pseudopotential we compute here the sumbec term
! ... symmetrization of the charge density (and local magnetization)
!
IF ( .NOT. okvan ) GO TO 10
#ifdef __PARA
!
CALL ccalbec( nkb, npwx, npw, nbnd, becp, vkb, evc )
! ... reduce charge density across pools
!
CALL start_clock( 'sumbec' )
CALL poolreduce( nspin * nrxx, rho )
!
DO ibnd = 1, nbnd
w1 = wg(ibnd, ik)
ijkb0 = 0
DO np = 1, ntyp
IF ( tvanp(np) ) THEN
DO na = 1, nat
IF ( ityp(na) == np ) THEN
ijh = 1
DO ih = 1, nh(np)
ikb = ijkb0 + ih
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
w1 * REAL( CONJG( becp(ikb,ibnd) ) * becp(ikb,ibnd))
ijh = ijh + 1
DO jh = ( ih + 1 ), nh(np)
jkb = ijkb0 + jh
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + w1 * 2.D0 * &
REAL( CONJG( becp(ikb,ibnd) ) * becp(jkb,ibnd) )
ijh = ijh + 1
END DO
END DO
ijkb0 = ijkb0 + nh(np)
END IF
END DO
ELSE
DO na = 1, nat
IF ( ityp(na) == np ) ijkb0 = ijkb0 + nh(np)
END DO
END IF
END DO
END DO
!
CALL stop_clock( 'sumbec' )
!
10 CONTINUE
!
END DO
!
! ... If a double grid is used, interpolate onto the fine grid
!
IF ( doublegrid ) THEN
DO is = 1, nspin
CALL interpolate( rho(1,is), rho(1,is), 1 )
!
CALL psymrho( rho(1,is), nrx1, nrx2, nrx3, &
nr1, nr2, nr3, nsym, s, ftau )
!
END DO
END IF
!
! ... Here we add the Ultrasoft contribution to the charge
!
IF ( okvan ) CALL addusdens()
!
#ifdef __PARA
CALL poolreduce( 1, eband )
CALL poolreduce( 1, demet )
#endif
!
! ... symmetrization of the charge density (and local magnetization)
!
#ifdef __PARA
!
! ... reduce charge density across pools
!
CALL poolreduce( nspin * nrxx, rho )
!
DO is = 1, nspin
CALL psymrho( rho(1,is), nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, &
s, ftau )
END DO
#else
DO is = 1, nspin
CALL symrho( rho(1,is), nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, &
s, ftau )
END DO
DO is = 1, nspin
!
CALL symrho( rho(1,is), nrx1, nrx2, nrx3, &
nr1, nr2, nr3, nsym, s, ftau )
!
END DO
#endif
!
RETURN
!
END SUBROUTINE sum_band_k
!
!
RETURN
!
END SUBROUTINE sum_band_gamma
!
!
!--------------------------------------------------------------------
SUBROUTINE sum_band_k()
!--------------------------------------------------------------------
!
! ... k-points version
!
USE becmod, ONLY : becp
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: ikb, jkb, ijkb0, ih, jh, ijh, na, np
! counters on beta functions, atoms, pseudopotentials
INTEGER :: ir, is, ig, ibnd, ik
! counter on 3D r points
! counter on spin polarizations
! counter on g vectors
! counter on bands
! counter on k points
REAL(KIND=DP) :: w1
! weight
!
!
becsum(:,:,:) = 0.D0
rho(:,:) = 0.D0
eband = 0.D0
demet = 0.D0
!
! ... calculate weights for the insulator case
!
IF ( .NOT. lgauss .AND. .NOT. ltetra .AND. .NOT. tfixed_occ ) THEN
!
CALL iweights( nks, wk, nbnd, nelec, et, ef, wg )
!
! ... calculate weights for the metallic case
!
ELSE IF ( ltetra ) THEN
#ifdef __PARA
!
CALL poolrecover( et, nbnd, nkstot, nks )
!
IF ( me == 1 .AND. mypool == 1 ) THEN
!
#endif
CALL tweights( nkstot, nspin, nbnd, &
nelec, ntetra, tetra, et, ef, wg )
#ifdef __PARA
!
ENDIF
!
CALL poolscatter( nbnd, nkstot, wg, nks, wg )
!
IF ( me == 1 ) CALL poolbcast( 1, ef )
!
CALL broadcast( 1, ef )
#endif
!
ELSE IF ( lgauss ) THEN
!
CALL gweights( nks, wk, nbnd, nelec, &
degauss, ngauss, et, ef, demet, wg )
!
ELSE IF ( tfixed_occ ) THEN
!
ef = - 1.0D+20
!
DO is = 1, nspin
!
DO ibnd = 1, nbnd
!
wg(ibnd,is) = f_inp(ibnd,is)
!
IF ( wg(ibnd,is) > 0.D0 ) ef = MAX( ef, et(ibnd,is) )
!
END DO
!
END DO
!
END IF
!
! ... Needed for LDA+U
!
IF ( lda_plus_u ) CALL new_ns()
!
! ... here we sum for each k point the contribution
! ... of the wavefunctions to the charge
!
IF ( nks > 1 ) REWIND( iunigk )
!
k_loop: DO ik = 1, nks
!
IF ( lsda ) current_spin = isk(ik)
!
IF ( nks > 1 ) THEN
!
READ( iunigk ) npw, igk
CALL davcio( evc, nwordwfc, iunwfc, ik, -1 )
!
END IF
!
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
!
! ... here we compute the band energy: the sum of the eigenvalues
!
DO ibnd = 1, nbnd
!
eband = eband + et(ibnd,ik) * wg(ibnd,ik)
!
! ... the sum of eband and demet is the integral for e < ef of
! ... e n(e) which reduces for degauss=0 to the sum of the
! ... eigenvalues
! ... the factor two is for spin degeneracy
!
psic(:) = ( 0.D0, 0.D0 )
!
psic(nls(igk(1:npw))) = evc(1:npw,ibnd)
!
CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
!
w1 = wg(ibnd,ik) / omega
!
! ... increment the charge density ...
!
DO ir = 1, nrxxs
!
rho(ir,current_spin) = rho(ir,current_spin) + &
w1 * ( REAL( psic(ir) )**2 + &
IMAG( psic(ir) )**2 )
!
END DO
!
END DO
!
! ... If we have a US pseudopotential we compute here the sumbec term
!
IF ( .NOT. okvan ) CYCLE k_loop
!
IF ( nkb > 0 ) &
CALL ccalbec( nkb, npwx, npw, nbnd, becp, vkb, evc )
!
CALL start_clock( 'sumbec' )
!
DO ibnd = 1, nbnd
!
w1 = wg(ibnd,ik)
ijkb0 = 0
!
DO np = 1, ntyp
!
IF ( tvanp(np) ) THEN
!
DO na = 1, nat
!
IF ( ityp(na) == np ) THEN
!
ijh = 1
!
DO ih = 1, nh(np)
!
ikb = ijkb0 + ih
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + &
w1 * REAL( CONJG( becp(ikb,ibnd) ) * &
becp(ikb,ibnd) )
!
ijh = ijh + 1
!
DO jh = ( ih + 1 ), nh(np)
!
jkb = ijkb0 + jh
!
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + w1 * 2.D0 * &
REAL( CONJG( becp(ikb,ibnd) ) * &
becp(jkb,ibnd) )
!
ijh = ijh + 1
!
END DO
!
END DO
!
ijkb0 = ijkb0 + nh(np)
!
END IF
!
END DO
!
ELSE
!
DO na = 1, nat
!
IF ( ityp(na) == np ) ijkb0 = ijkb0 + nh(np)
!
END DO
!
END IF
!
END DO
!
END DO
!
CALL stop_clock( 'sumbec' )
!
END DO k_loop
!
! ... If a double grid is used, interpolate onto the fine grid
!
IF ( doublegrid ) THEN
!
DO is = 1, nspin
!
CALL interpolate( rho(1,is), rho(1,is), 1 )
!
END DO
!
END IF
!
! ... Here we add the Ultrasoft contribution to the charge
!
IF ( okvan ) CALL addusdens()
!
#ifdef __PARA
CALL poolreduce( 1, eband )
CALL poolreduce( 1, demet )
#endif
!
! ... symmetrization of the charge density (and local magnetization)
!
#ifdef __PARA
!
! ... reduce charge density across pools
!
CALL poolreduce( nspin * nrxx, rho )
!
DO is = 1, nspin
!
CALL psymrho( rho(1,is), nrx1, nrx2, nrx3, &
nr1, nr2, nr3, nsym, s, ftau )
!
END DO
#else
DO is = 1, nspin
!
CALL symrho( rho(1,is), nrx1, nrx2, nrx3, &
nr1, nr2, nr3, nsym, s, ftau )
!
END DO
#endif
!
RETURN
!
END SUBROUTINE sum_band_k
!
END SUBROUTINE sum_band

6
PW/wfcinit.f90
View File

@ -175,7 +175,8 @@ SUBROUTINE wfcinit()
END DO
END IF
!
CALL init_us_2( npw, igk, xk(1,ik), vkb )
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
!
! ... Diagonalize the Hamiltonian on the basis of atomic wfcs
!
@ -336,7 +337,8 @@ SUBROUTINE wfcinit()
!
END IF
!
CALL init_us_2( npw, igk, xk(1,ik), vkb )
IF ( nkb > 0 ) &
CALL init_us_2( npw, igk, xk(1,ik), vkb )
!
! ... Diagonalize the Hamiltonian on the basis of atomic wfcs
!