From 0b32850e1f027d3ad6bb6389e62ba8812ecc57dc Mon Sep 17 00:00:00 2001 From: sbraccia Date: Thu, 22 Jan 2004 12:01:13 +0000 Subject: [PATCH] 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 --- PW/c_bands.f90 | 6 +- PW/force_us.f90 | 27 +- PW/h_psi.f90 | 6 +- PW/pw_gemm.f90 | 6 +- PW/sum_band.f90 | 893 +++++++++++++++++++++++++++--------------------- PW/wfcinit.f90 | 6 +- 6 files changed, 534 insertions(+), 410 deletions(-) diff --git a/PW/c_bands.f90 b/PW/c_bands.f90 index 7b2e06f27..3d15d61b2 100644 --- a/PW/c_bands.f90 +++ b/PW/c_bands.f90 @@ -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 ! diff --git a/PW/force_us.f90 b/PW/force_us.f90 index 3da7a3c4d..34e10f8d6 100644 --- a/PW/force_us.f90 +++ b/PW/force_us.f90 @@ -60,7 +60,7 @@ SUBROUTINE force_us( forcenl ) ! REAL(KIND=DP), ALLOCATABLE :: becp(:,:), dbecp (:,:,:) ! auxiliary variables contain and - 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 - diff --git a/PW/h_psi.f90 b/PW/h_psi.f90 index 82f715b14..e69acda85 100644 --- a/PW/h_psi.f90 +++ b/PW/h_psi.f90 @@ -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 ! diff --git a/PW/pw_gemm.f90 b/PW/pw_gemm.f90 index 1b67712ca..b356140a3 100644 --- a/PW/pw_gemm.f90 +++ b/PW/pw_gemm.f90 @@ -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 ! diff --git a/PW/sum_band.f90 b/PW/sum_band.f90 index 89d71d202..51976a8db 100644 --- a/PW/sum_band.f90 +++ b/PW/sum_band.f90 @@ -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 - ! - ! - 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 + ! + ! + 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 diff --git a/PW/wfcinit.f90 b/PW/wfcinit.f90 index 2f1781c02..1cffe66ea 100644 --- a/PW/wfcinit.f90 +++ b/PW/wfcinit.f90 @@ -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 !