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

[skip-CI] Reorganization of rho core code (again)

This commit is contained in:
Paolo Giannozzi 2023-10-07 10:18:17 +02:00
parent c125310efa
commit 28d73e2205
10 changed files with 251 additions and 204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
PW/src/force_cc.f90
View File

@ -27,6 +27,7 @@ SUBROUTINE force_cc( forcecc )
USE noncollin_module, ONLY : noncolin
USE mp_bands, ONLY : intra_bgrp_comm
USE mp, ONLY : mp_sum
USE rhoc_mod, ONLY : interp_rhc
!
IMPLICIT NONE
!

2
PW/src/init_vloc.f90
View File

@ -45,7 +45,7 @@ SUBROUTINE init_vloc()
qmax = tpiba2 * MAXVAL ( gl )
CALL mp_max (qmax, intra_bgrp_comm)
!! this is the actual maximum |G|^2 needed in the interpolation table
!! for variable-cell calculations this may exceed ecutrho, so we use
!! for variable-cell calculations. It may exceed ecutrho, so we use
!! "cell_factor" (1.2 or so) as below, in order to avoid too frequent
!! re-allocations of the interpolation table
!

10
PW/src/set_rhoc.f90
View File

@ -22,17 +22,19 @@ SUBROUTINE set_rhoc
USE cell_base, ONLY : omega, tpiba2
USE fft_base, ONLY : dfftp
USE fft_rho, ONLY : rho_g2r
USE gvect, ONLY : ngm, ngl, gl, igtongl
USE gvect, ONLY : ngm, ngl, gl, igtongl, ecutrho
USE vlocal, ONLY : strf
USE mp_bands, ONLY : intra_bgrp_comm
USE mp, ONLY : mp_sum
USE scf, ONLY : rho_core, rhog_core
USE cellmd, ONLY : cell_factor
USE rhoc_mod, ONLY : init_tab_rhc, interp_rhc
!
IMPLICIT NONE
!
REAL(DP) , ALLOCATABLE :: rhocg(:)
! the radial fourier transform
REAL(DP) :: rhoneg
REAL(DP) :: qmax, rhoneg
! used to check the core charge
INTEGER :: ir, nt, ng
! counter on mesh points
@ -43,6 +45,10 @@ SUBROUTINE set_rhoc
rho_core(:) = 0.0_DP
IF ( ANY( upf(1:ntyp)%nlcc ) ) THEN
!
qmax = sqrt(ecutrho)*cell_factor
CALL init_tab_rhc ( qmax, omega, intra_bgrp_comm, ir )
!
ALLOCATE (rhocg( ngl))
!$acc data create(rhocg) copyin(gl, igtongl, strf, rhog_core, rho_core)
!

1
PW/src/stres_cc.f90
View File

@ -20,6 +20,7 @@ SUBROUTINE stres_cc( sigmaxcc )
USE gvect, ONLY : ngm, gstart, ngl, gl, igtongl, g, gg
USE ener, ONLY : etxc, vtxc
USE lsda_mod, ONLY : nspin
USE rhoc_mod, ONLY : interp_rhc, interp_drhc
USE scf, ONLY : rho, rho_core, rhog_core
USE vlocal, ONLY : strf
USE control_flags, ONLY : gamma_only

4
upflib/Makefile
View File

@ -11,6 +11,7 @@ MODFLAGS= $(MOD_FLAG)../UtilXlib $(MOD_FLAG)../external/devxlib/src
OBJS_DEP= \
vloc_mod.o \
rhoc_mod.o \
gen_us_dj.o \
gen_us_dy.o \
init_us_0.o \
@ -20,7 +21,6 @@ init_us_2_base.o \
init_tab_atwfc.o \
init_tab_beta.o \
init_tab_rho.o \
init_tab_rhc.o \
init_tab_qrad.o
OBJS_NODEP= \
@ -30,8 +30,6 @@ dqvan2.o \
dylmr2.o \
gth.o \
interp_atwfc.o \
interp_rhc.o \
interp_drhc.o \
paw_variables.o \
pseudo_types.o \
qvan2.o \

71
upflib/init_tab_rhc.f90
View File

@ -1,71 +0,0 @@
!
! Copyright (C) 2023 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------
SUBROUTINE init_tab_rhc (omega, comm)
!----------------------------------------------------------------------
!
!! Compute interpolation table for atomic core (pseudo-)charge density
!! tab_rhc(i,nt) = \rhoc_nt(q_i) for atom of type nt, on grid q_i
!! Output in tab_rhc, module uspp_data
!
USE upf_kinds, ONLY : dp
USE upf_const, ONLY : fpi, eps8
USE atom, ONLY : rgrid, msh
USE uspp_param, ONLY : upf, nsp
USE uspp_data, ONLY : nqxq, dq, tab_rhc
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: comm
!! MPI communicator, to split the workload
REAL(dp), INTENT(IN) :: omega
!! Unit-cell volume
!
INTEGER :: ndm, startq, lastq, nt, iq, ir
!! Various indices
REAL(dp) :: q
REAL(dp), ALLOCATABLE :: aux (:)
!! Work space
!
IF ( .NOT. ANY(upf(1:nsp)%nlcc) ) RETURN
ndm = MAXVAL( msh(1:nsp) )
ALLOCATE (aux(ndm))
!
CALL divide (comm, nqxq, startq, lastq)
!
DO nt = 1, nsp
!
tab_rhc(:,nt)= 0.d0
IF ( .NOT. upf(nt)%nlcc ) CYCLE
!
DO iq = startq, lastq
!
q = (iq - 1) * dq
DO ir = 1, msh(nt)
IF ( iq > 1 .AND. rgrid(nt)%r(ir) > eps8 ) then
aux (ir) = upf(nt)%rho_atc(ir) * rgrid(nt)%r2(ir) * &
sin(q*rgrid(nt)%r(ir)) / (q*rgrid(nt)%r(ir))
ELSE
aux (ir) = upf(nt)%rho_atc(ir) * rgrid(nt)%r2(ir)
ENDIF
ENDDO
!
CALL simpson ( msh(nt), aux, rgrid(nt)%rab, tab_rhc(iq,nt) )
tab_rhc (iq,nt) = fpi * tab_rhc (iq,nt) / omega
!
ENDDO
!
END DO
!
CALL mp_sum ( tab_rhc (:,1:nsp), comm )
!$acc update device(tab_rhc)
!
DEALLOCATE (aux)
!
END SUBROUTINE init_tab_rhc

11
upflib/init_us_1.f90
View File

@ -309,7 +309,6 @@ subroutine init_us_1( nat, ityp, omega, ngm, g, gg, intra_bgrp_comm )
!
CALL init_tab_beta ( omega, intra_bgrp_comm )
CALL init_tab_rho ( omega, intra_bgrp_comm )
CALL init_tab_rhc ( omega, intra_bgrp_comm )
!
#if defined __CUDA
!
@ -322,9 +321,11 @@ subroutine init_us_1( nat, ityp, omega, ngm, g, gg, intra_bgrp_comm )
nhtoj_d=nhtoj
ijtoh_d=ijtoh
qq_nt_d=qq_nt
!$acc update device(qq_at)
if (is_spinorbit) then
dvan_so_d=dvan_so
fcoef_d=fcoef
!$acc update device(qq_so)
else
dvan_d=dvan
endif
@ -332,14 +333,6 @@ subroutine init_us_1( nat, ityp, omega, ngm, g, gg, intra_bgrp_comm )
ofsbeta_d=ofsbeta
!
#endif
!
if (nhm>0) then
!$acc update device(qq_at)
if (is_spinorbit) then
!$acc update device(qq_so)
endif
endif
!
call stop_clock ('init_us_1')
return
!

60
upflib/interp_drhc.f90
View File

@ -1,60 +0,0 @@
!
! Copyright (C) 2023 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------------
SUBROUTINE interp_drhc( nt, ngl, gl, tpiba2, drhocg )
!--------------------------------------------------------------------------
!! Calculates the Fourier transform of \(d\text{Rho}_c/dG\).
!
USE upf_kinds, ONLY : dp
USE uspp_data, ONLY : tab_rhc, dq
!
IMPLICIT NONE
!
INTEGER :: nt
!! input: atomic type
INTEGER :: ngl
!! input: the number of g shell
REAL(DP), INTENT(IN) :: gl(ngl)
!! input: the number of G shells
REAL(DP), INTENT(IN) :: tpiba2
!! input: 2 times pi / alat
REAL(DP), INTENT(OUT) :: drhocg(ngl)
!! Fourier transform of d Rho_c/dG
!
! ... local variables
!
REAL(DP) :: gx, px, ux, vx, wx
! the modulus of g for a given shell
! variables used for interpolation
INTEGER :: igl, i0, i1, i2, i3
! counters
!
!$acc data present_or_copyin(gl) present_or_copyout(drhocg)
!$acc parallel loop
DO igl = 1, ngl
gx = SQRT( gl(igl) * tpiba2 )
px = gx / dq - int (gx/dq)
ux = 1.d0 - px
vx = 2.d0 - px
wx = 3.d0 - px
i0 = INT(gx/dq) + 1
i1 = i0 + 1
i2 = i0 + 2
i3 = i0 + 3
drhocg (igl) = (- tab_rhc(i0, nt) * (ux*vx + vx*wx + ux*wx) / 6.0_dp &
+ tab_rhc(i1, nt) * (wx*vx - px*wx - px*vx) / 2.0_dp &
- tab_rhc(i2, nt) * (wx*ux - px*wx - px*ux) / 2.0_dp &
+ tab_rhc(i3, nt) * (ux*vx - px*ux - px*vx) / 6.0_dp ) / dq
ENDDO
!
!$acc end data
!
RETURN
!
END SUBROUTINE interp_drhc

58
upflib/interp_rhc.f90
View File

@ -1,58 +0,0 @@
!
! Copyright (C) 2023 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!-----------------------------------------------------------------------
SUBROUTINE interp_rhc( nt, ngl, gl, tpiba2, rhocg )
!-----------------------------------------------------------------------
!! Calculates the radial Fourier transform of the core charge.
!
USE upf_kinds, ONLY : dp
USE uspp_data, ONLY : tab_rhc, dq
!
IMPLICIT NONE
!
INTEGER :: nt
!! input: atomic type
INTEGER :: ngl
!! input: the number of g shell
REAL(DP) :: gl(ngl)
!! input: the number of G shells
REAL(DP) :: tpiba2
!! input: 2 times pi / alat
REAL(DP) :: rhocg(ngl)
!! output: the Fourier transform of the core charge
!
! ... local variables
!
REAL(DP) :: gx, px, ux, vx, wx
! the modulus of g for a given shell
! variables used for interpolation
INTEGER :: igl, i0, i1, i2, i3
! counters
!
!$acc data present_or_copyin(gl) present_or_copyout(rhocg) present(tab_rhc)
!$acc parallel loop
DO igl = 1, ngl
gx = SQRT(gl(igl) * tpiba2)
px = gx / dq - int (gx/dq)
ux = 1.d0 - px
vx = 2.d0 - px
wx = 3.d0 - px
i0 = INT(gx/dq) + 1
i1 = i0 + 1
i2 = i0 + 2
i3 = i0 + 3
rhocg (igl) = tab_rhc(i0, nt) * ux * vx * wx / 6.d0 + &
tab_rhc(i1, nt) * px * vx * wx / 2.d0 - &
tab_rhc(i2, nt) * px * ux * wx / 2.d0 + &
tab_rhc(i3, nt) * px * ux * vx / 6.d0
ENDDO
!$acc end data
!
END SUBROUTINE interp_rhc

237
upflib/rhoc_mod.f90 Normal file
View File

@ -0,0 +1,237 @@
!
! Copyright (C) 2023 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
MODULE rhoc_mod
!
!! Variables and routines for nonlinear core correction
!! Contains generation of interpolation tables in reciprocal space,
!! interpolation routines and other utility routines
!! Code moved to upflib and restructured by Paolo Giannozzi, 2023
!
USE upf_kinds, ONLY : dp
USE upf_const, ONLY : fpi, e2, eps8
!
IMPLICIT NONE
!
PRIVATE
PUBLIC :: init_tab_rhc
PUBLIC :: deallocate_tab_rhc
PUBLIC :: scale_tab_rhc
PUBLIC :: interp_rhc
PUBLIC :: interp_drhc
!
SAVE
!
INTEGER :: nqx = 0
!! size of interpolation table
REAL(DP), PARAMETER:: dq = 0.01_dp
!! grid step for interpolation table
REAL(DP) :: qmax = 0.0_dp
!! max q covered by the interpolation table
REAL(DP), ALLOCATABLE :: tab_rhc(:,:)
!! interpolation table for atomic pseudo-core charge density
!
CONTAINS
!
!----------------------------------------------------------------------
SUBROUTINE init_tab_rhc (qmax_, omega, comm, ierr)
!----------------------------------------------------------------------
!
!! Compute interpolation table for atomic core (pseudo-)charge density
!! tab_rhc(i,nt) = \rhoc_nt(q_i) for atom of type nt, on grid q_i
!! Output in tab_rhc
!
USE atom, ONLY : rgrid, msh
USE uspp_param, ONLY : upf, nsp
USE mp, ONLY : mp_sum
!
INTEGER, INTENT(IN) :: comm
!! MPI communicator, to split the workload
INTEGER, INTENT(OUT) :: ierr
!! return code: ierr = 0 if interpolation table (IT) was allocated
!! ierr =-1 if IT had insufficent dimension and was re-allocated
!! ierr =-2 if IT was already present and nothing is done
REAL(dp), INTENT(IN) :: qmax_
!! Interpolate q up to qmax_ (sqrt(Ry), q^2 is an energy)
REAL(dp), INTENT(IN) :: omega
!! Unit-cell volume
!
INTEGER :: ndm, startq, lastq, nt, iq, ir
!! Various indices
REAL(dp) :: q
REAL(dp), ALLOCATABLE :: aux (:)
!! Work space
!
ierr = 0
IF ( .NOT. ANY(upf(1:nsp)%nlcc) ) RETURN
!
IF ( .NOT. ALLOCATED(tab_rhc) ) THEN
!! table not yet allocated
qmax = qmax_
nqx = INT( qmax/dq + 4)
ALLOCATE ( tab_rhc(0:nqx,nsp) )
!$acc enter data create(tab_rhc)
ELSE IF ( qmax_ > qmax ) THEN
DEALLOCATE ( tab_rhc )
!! table ìs allocated but dimension insufficient: re-allocate
!! (with some margin so that this does not happen too often)
qmax = qmax_ + MAX(dq,qmax_-qmax) * 10
nqx = INT( qmax/dq + 4)
ALLOCATE ( tab_rhc(0:nqx,nsp) )
!$acc enter data create(tab_rhc)
ierr =-1
ELSE
ierr =-2
RETURN
END IF
!
ndm = MAXVAL( msh(1:nsp) )
ALLOCATE (aux(ndm))
!
CALL divide (comm, nqx, startq, lastq)
!
DO nt = 1, nsp
!
tab_rhc(:,nt)= 0.d0
IF ( .NOT. upf(nt)%nlcc ) CYCLE
!
DO iq = startq, lastq
!
q = (iq - 1) * dq
DO ir = 1, msh(nt)
IF ( iq > 1 .AND. rgrid(nt)%r(ir) > eps8 ) then
!! check above prevents divide by zero
aux (ir) = upf(nt)%rho_atc(ir) * rgrid(nt)%r2(ir) * &
sin(q*rgrid(nt)%r(ir)) / (q*rgrid(nt)%r(ir))
ELSE
aux (ir) = upf(nt)%rho_atc(ir) * rgrid(nt)%r2(ir)
ENDIF
ENDDO
!
CALL simpson ( msh(nt), aux, rgrid(nt)%rab, tab_rhc(iq,nt) )
tab_rhc (iq,nt) = fpi * tab_rhc (iq,nt) / omega
!
ENDDO
!
END DO
!
CALL mp_sum ( tab_rhc (:,1:nsp), comm )
!$acc update device(tab_rhc)
!
DEALLOCATE (aux)
RETURN
!
END SUBROUTINE init_tab_rhc
!
!-----------------------------------------------------------------------
SUBROUTINE interp_rhc( nt, ngl, gl, tpiba2, rhocg )
!-----------------------------------------------------------------------
!! Calculates the radial Fourier transform of the core charge.
!
INTEGER :: nt
!! input: atomic type
INTEGER :: ngl
!! input: the number of g shell
REAL(DP) :: gl(ngl)
!! input: the number of G shells
REAL(DP) :: tpiba2
!! input: 2 times pi / alat
REAL(DP) :: rhocg(ngl)
!! output: the Fourier transform of the core charge
!
! ... local variables
!
REAL(DP) :: gx, px, ux, vx, wx
! the modulus of g for a given shell
! variables used for interpolation
INTEGER :: igl, i0, i1, i2, i3
! counters
!
!$acc data present_or_copyin(gl) present_or_copyout(rhocg) present(tab_rhc)
!$acc parallel loop
DO igl = 1, ngl
gx = SQRT(gl(igl) * tpiba2)
px = gx / dq - int (gx/dq)
ux = 1.d0 - px
vx = 2.d0 - px
wx = 3.d0 - px
i0 = INT(gx/dq) + 1
i1 = i0 + 1
i2 = i0 + 2
i3 = i0 + 3
rhocg (igl) = tab_rhc(i0, nt) * ux * vx * wx / 6.d0 + &
tab_rhc(i1, nt) * px * vx * wx / 2.d0 - &
tab_rhc(i2, nt) * px * ux * wx / 2.d0 + &
tab_rhc(i3, nt) * px * ux * vx / 6.d0
ENDDO
!$acc end data
!
END SUBROUTINE interp_rhc
!
!----------------------------------------------------------------------------
SUBROUTINE interp_drhc( nt, ngl, gl, tpiba2, drhocg )
!--------------------------------------------------------------------------
!! Calculates the Fourier transform of \(d\text{Rho}_c/dG\).
!
INTEGER :: nt
!! input: atomic type
INTEGER :: ngl
!! input: the number of g shell
REAL(DP), INTENT(IN) :: gl(ngl)
!! input: the number of G shells
REAL(DP), INTENT(IN) :: tpiba2
!! input: 2 times pi / alat
REAL(DP), INTENT(OUT) :: drhocg(ngl)
!! Fourier transform of d Rho_c/dG
!
! ... local variables
!
REAL(DP) :: gx, px, ux, vx, wx
! the modulus of g for a given shell
! variables used for interpolation
INTEGER :: igl, i0, i1, i2, i3
! counters
!
!$acc data present_or_copyin(gl) present_or_copyout(drhocg)
!$acc parallel loop
DO igl = 1, ngl
gx = SQRT( gl(igl) * tpiba2 )
px = gx / dq - int (gx/dq)
ux = 1.d0 - px
vx = 2.d0 - px
wx = 3.d0 - px
i0 = INT(gx/dq) + 1
i1 = i0 + 1
i2 = i0 + 2
i3 = i0 + 3
drhocg (igl) = (- tab_rhc(i0, nt) * (ux*vx + vx*wx + ux*wx) / 6.0_dp &
+ tab_rhc(i1, nt) * (wx*vx - px*wx - px*vx) / 2.0_dp &
- tab_rhc(i2, nt) * (wx*ux - px*wx - px*ux) / 2.0_dp &
+ tab_rhc(i3, nt) * (ux*vx - px*ux - px*vx) / 6.0_dp ) / dq
ENDDO
!$acc end data
!
END SUBROUTINE interp_drhc
!
subroutine deallocate_tab_rhc
!$acc exit data delete(tab_rhc)
if( allocated( tab_rhc) ) deallocate( tab_rhc)
end subroutine deallocate_tab_rhc
!
subroutine scale_tab_rhc( vol_ratio_m1 )
! vol_ratio_m1 = omega_old / omega
real(DP), intent(in) :: vol_ratio_m1
!
tab_rhc(:,:) = tab_rhc(:,:) * vol_ratio_m1
!$acc update device (tab_rhc)
end subroutine scale_tab_rhc
!
END MODULE rhoc_mod