From 2f9ecea6c279b0955941baa2472316b45189d1f8 Mon Sep 17 00:00:00 2001 From: Paolo Giannozzi Date: Thu, 14 Sep 2023 15:29:01 +0200 Subject: [PATCH] Cleanup of vloc_of_g --- PW/src/init_vloc.f90 | 27 ++------ upflib/vloc_mod.f90 | 147 ++++++++++++++++++------------------------- 2 files changed, 65 insertions(+), 109 deletions(-) diff --git a/PW/src/init_vloc.f90 b/PW/src/init_vloc.f90 index a771d5e9c..056f6926a 100644 --- a/PW/src/init_vloc.f90 +++ b/PW/src/init_vloc.f90 @@ -13,8 +13,7 @@ SUBROUTINE init_vloc() !! potential vloc(ig,it) for each type of atom. ! USE kinds, ONLY : dp - USE vloc_mod, ONLY : vloc_of_g, vloc_coul - USE m_gth, ONLY : vloc_gth + USE vloc_mod, ONLY : vloc_of_g USE uspp_param, ONLY : upf USE ions_base, ONLY : ntyp => nsp USE cell_base, ONLY : omega, tpiba2 @@ -37,32 +36,14 @@ SUBROUTINE init_vloc() modified_coulomb = do_cutoff_2D .OR. (do_comp_esm .and. ( esm_bc .ne. 'pbc' )) ! DO nt = 1, ntyp - IF (do_cutoff_2D .AND. rgrid(nt)%r(msh(nt)) > lz) THEN + IF (rgrid(nt)%r(msh(nt)) > lz) THEN call errore('init_vloc','2D cutoff smaller than pseudo cutoff radius: & - & increase interlayer distance (or see Modules/read_pseudo.f90)',1) + & increase interlayer distance (or see Modules/read_pseudo.f90)',nt) END IF ! ! compute V_loc(G) for a given type of atom ! - IF ( upf(nt)%is_gth ) THEN - ! - ! special case: GTH pseudopotential - ! - CALL vloc_gth( nt, upf(nt)%zp, tpiba2, ngl, gl, omega, vloc(1,nt) ) - ! - ELSE IF ( upf(nt)%tcoulombp ) THEN - ! - ! special case: pseudopotential is coulomb 1/r potential - ! - CALL vloc_coul( upf(nt)%zp, tpiba2, ngl, gl, omega, vloc(1,nt) ) - ! - ELSE - ! - ! normal case - ! - CALL vloc_of_g( nt, ngl, gl, tpiba2, modified_coulomb, omega, vloc(1,nt) ) - ! - ENDIF + CALL vloc_of_g( nt, ngl, gl, tpiba2, modified_coulomb, omega, vloc(1,nt) ) ! ENDDO ! diff --git a/upflib/vloc_mod.f90 b/upflib/vloc_mod.f90 index d9ab4afd0..a797e0a7f 100644 --- a/upflib/vloc_mod.f90 +++ b/upflib/vloc_mod.f90 @@ -17,7 +17,7 @@ MODULE vloc_mod IMPLICIT NONE ! PRIVATE - PUBLIC :: vloc_of_g, vloc_coul + PUBLIC :: vloc_of_g PUBLIC ::dvloc_of_g,dvloc_coul ! CONTAINS @@ -35,7 +35,7 @@ SUBROUTINE vloc_of_g( nt, ngl, gl, tpiba2, modified_coulomb, omega, vloc ) ! USE atom, ONLY : rgrid, msh USE uspp_param, ONLY : upf - ! USE m_gth, ONLY : vloc_gth + USE m_gth, ONLY : vloc_gth ! IMPLICIT NONE ! @@ -65,101 +65,76 @@ SUBROUTINE vloc_of_g( nt, ngl, gl, tpiba2, modified_coulomb, omega, vloc ) ! igl0: position of first nonzero G ! ir : counter on mesh points ! - ALLOCATE (aux(msh(nt)), aux1(msh(nt))) - ! IF (gl (1) < eps8) THEN + igl0 = 2 + ELSE + igl0 = 1 + END IF + IF ( upf(nt)%is_gth ) THEN + ! special case: GTH pseudopotential + CALL vloc_gth( nt, upf(nt)%zp, tpiba2, ngl, gl, omega, vloc ) + ELSE IF ( upf(nt)%tcoulombp ) THEN + ! special case: pure Coulomb pseudopotential + IF ( igl0 > 1 ) vloc(1) = 0.0_dp + vloc (igl0:ngl) = - fpi * upf(nt)%zp*e2 / omega / tpiba2 / gl (igl0:ngl) + ELSE + ! normal case + ALLOCATE (aux(msh(nt)), aux1(msh(nt))) ! - IF ( modified_coulomb ) THEN - ! G=0 limit for the modified Coulomb potential (ESM, 2D) + IF ( igl0 == 2 ) THEN + ! + IF ( modified_coulomb ) THEN + ! G=0 limit for the modified Coulomb potential (ESM, 2D) + DO ir = 1, msh(nt) + aux (ir) = rgrid(nt)%r(ir) * (rgrid(nt)%r(ir)*upf(nt)%vloc (ir) +& + upf(nt)%zp * e2 * erf (rgrid(nt)%r(ir)) ) + END DO + ELSE + ! G=0 limit for Coulomb potential + DO ir = 1, msh(nt) + aux (ir) = rgrid(nt)%r(ir) * (rgrid(nt)%r(ir)*upf(nt)%vloc (ir) +& + upf(nt)%zp * e2 ) + END DO + END IF + ! + CALL simpson (msh(nt), aux, upf(nt)%rab, vloc(1)) + ! + END IF + ! here the G<>0 terms with long range removed + ! + ! G-independent term + DO ir = 1, msh(nt) + aux (ir) = rgrid(nt)%r(ir)*upf(nt)%vloc (ir) + & + upf(nt)%zp * e2 * erf (rgrid(nt)%r(ir)) + END DO + ! + ! and here we perform the integral, after multiplying for the |G| + ! dependent part + ! + DO igl = igl0, ngl + gx = sqrt ( gl(igl)*tpiba2 ) DO ir = 1, msh(nt) - aux (ir) = rgrid(nt)%r(ir) * ( rgrid(nt)%r(ir)*upf(nt)%vloc (ir) + & - upf(nt)%zp * e2 * erf (rgrid(nt)%r(ir)) ) + aux1(ir) = aux (ir) * sin ( gx*upf(nt)%r(ir) ) / gx END DO - ELSE - ! G=0 limit for Coulomb potential - DO ir = 1, msh(nt) - aux (ir) = rgrid(nt)%r(ir) * ( rgrid(nt)%r(ir)*upf(nt)%vloc (ir) + & - upf(nt)%zp * e2 ) + CALL simpson (msh(nt), aux1, upf(nt)%rab, vloc(igl)) + END DO + DEALLOCATE (aux, aux1) + ! + IF ( .not. modified_coulomb ) THEN + fac = upf(nt)%zp * e2 / tpiba2 + DO igl = igl0, ngl + ! + ! here we re-add the analytic fourier transform of the erf function + ! + vloc(igl) = vloc(igl) - fac * exp (-gl (igl)*tpiba2*0.25d0) / gl (igl) END DO END IF - ! - CALL simpson (msh(nt), aux, upf(nt)%rab, vloc(1)) - ! - igl0 = 2 - ! - ELSE - ! - igl0 = 1 - ! + vloc (:) = vloc(:) * fpi / omega END IF - ! here the G<>0 terms with long range removed - ! - ! G-independent term - DO ir = 1, msh(nt) - aux (ir) = rgrid(nt)%r(ir)*upf(nt)%vloc (ir) + & - upf(nt)%zp * e2 * erf (rgrid(nt)%r(ir)) - END DO - ! - ! and here we perform the integral, after multiplying for the |G| - ! dependent part - ! - DO igl = igl0, ngl - gx = sqrt ( gl(igl)*tpiba2 ) - DO ir = 1, msh(nt) - aux1(ir) = aux (ir) * sin ( gx*upf(nt)%r(ir) ) / gx - END DO - CALL simpson (msh(nt), aux1, upf(nt)%rab, vloc(igl)) - END DO - ! - IF ( .not. modified_coulomb ) THEN - fac = upf(nt)%zp * e2 / tpiba2 - DO igl = igl0, ngl - ! - ! here we re-add the analytic fourier transform of the erf function - ! - vloc(igl) = vloc(igl) - fac * exp (-gl (igl)*tpiba2*0.25d0) / gl (igl) - END DO - END IF - vloc (:) = vloc(:) * fpi / omega - DEALLOCATE (aux, aux1) END SUBROUTINE vloc_of_g ! !---------------------------------------------------------------------- -SUBROUTINE vloc_coul( zp, tpiba2, ngl, gl, omega, vloc ) - !---------------------------------------------------------------------- - !! Fourier transform of the Coulomb potential - For all-electron - !! calculations, in specific cases only, for testing purposes. - ! - INTEGER, INTENT(IN) :: ngl - !! the number of shells of G vectors - REAL(DP), INTENT(IN) :: zp - !! valence pseudocharge - REAL(DP), INTENT(IN) :: tpiba2 - !! 2 pi / alat - REAL(DP), INTENT(IN) :: omega - !! the volume of the unit cell - REAL(DP), INTENT(IN) :: gl(ngl) - !! the moduli of g vectors for each shell - ! - ! ... local variables - ! - REAL(DP), INTENT(OUT) :: vloc(ngl) - ! the fourier transform of the potential - INTEGER :: igl0 - ! - IF (gl (1) < eps8) THEN - igl0 = 2 - vloc(1) = 0.0_dp - ELSE - igl0 = 1 - END IF - - vloc (igl0:ngl) = - fpi * zp *e2 / omega / tpiba2 / gl (igl0:ngl) - -END SUBROUTINE vloc_coul -! -!---------------------------------------------------------------------- SUBROUTINE dvloc_of_g( mesh, msh, rab, r, vloc_at, zp, tpiba2, ngl, gl, & modified_coulomb, omega, dvloc ) !----------------------------------------------------------------------