More pseudopotential reorganization

Calculation of interpolation tables moved to separate routines init_tab_*.f90
This commit is contained in:
Paolo Giannozzi 2021-04-17 08:25:38 +02:00 committed by Sasha Fonari
parent c56463b5bf
commit 9419b5e8a8
11 changed files with 222 additions and 195 deletions

View File

@ -42,10 +42,10 @@ Fixed in development version:
Incompatible changes in development version:
* lfcpopt & lfcpdyn are replaced by lfcp. Only static optimization of the
Fermi energy works. Currently, molecular dynamics with FCP does not work.
* Exchange-correlation code moved to XClib/
* More pseudopotential-related code moved to upflib/
* Pseudopotential-related duplicated variable indv_ijkb0 merged into
ofsbeta ("offset of beta functions")
* Exchange-correlation code (vdW excepted) moved to XClib/
* Much more pseudopotential-related code re-organized and moved to upflib/
Variables ofsbeta ("offset of beta functions") and indv_ijkb0 merged
Now only ofsbeta is used
* Code computing [H,x] commutator moved from LR_Modules/ to PW/src/ and
disentangled from linear-response variables - now is used in PP/ as well

View File

@ -244,7 +244,7 @@
!
qrad(:, :, :, :) = zero
! RM - need to call init_us_1 to re-calculate qrad
! PG - maybe it would be preferable to call compute_qrad?
! PG - maybe it would be sufficient to call init_tab_qrad?
CALL init_us_1(nat, ityp, omega, ngm, g, gg, intra_bgrp_comm)
ENDIF
ENDIF

View File

@ -20,7 +20,7 @@ SUBROUTINE scale_h
USE constants, ONLY : eps8
USE gvect, ONLY : g, gg, ngm
USE klist, ONLY : xk, wk, nkstot
USE uspp_data, ONLY : nqxq, qrad, tab, tab_at, dq, scale_uspp_data
USE uspp_data, ONLY : nqxq, dq, scale_uspp_data
USE control_flags, ONLY : iverbosity
USE start_k, ONLY : nks_start, xk_start, nk1,nk2,nk3
USE exx_base, ONLY : exx_grid_init, exx_mp_init

View File

@ -6,6 +6,8 @@ set(src_upflib
init_us_b0.f90
init_us_1.f90
init_tab_atwfc.f90
init_tab_beta.f90
init_tab_qrad.f90
interp_atwfc.f90
init_us_2_base.f90
gth.f90

View File

@ -16,6 +16,8 @@ init_us_0.o \
init_us_b0.o \
init_us_1.o \
init_tab_atwfc.o \
init_tab_beta.o \
init_tab_qrad.o \
init_us_2_base.o \
interp_atwfc.o \
paw_variables.o \

View File

@ -35,9 +35,9 @@
- Names of interpolation tables and related routines are random:
CP PW new name? contains computed in
betagx tab tab_beta beta(G) functions compute_betagx,
tab_d2y tab_beta_d2y splines for beta(G) compute_beta
tab_d2y tab_beta_d2y splines for beta(G) init_tab_beta
dbetagx dbeta(G)/dG compute_betagx
qradx qrad tab_q Q(G) for USPP/PAW compute_qrad,
qradx qrad tab_q Q(G) for USPP/PAW init_tab_qrad
dqradx dQ(G)/dG compute_qradx
tab_atwfc atomic R_nl(G) init_tab_atwfc
(tab_atrho atomic rho(G) to be done)

93
upflib/init_tab_beta.f90 Normal file
View File

@ -0,0 +1,93 @@
!
! Copyright (C) 2021 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_beta ( omega, intra_bgrp_comm )
!----------------------------------------------------------------------
!
! Compute interpolation table for beta(G) radial functions
!
USE upf_kinds, ONLY : dp
USE upf_const, ONLY : fpi
USE atom, ONLY : rgrid
USE uspp_param, ONLY : upf, lmaxq, nbetam, nsp
USE uspp_data, ONLY : nqx, dq, tab, tab_d2y, spline_ps, tab_d, tab_d2y_d
USE mp, ONLY : mp_sum
USE splinelib, ONLY : spline
!
IMPLICIT NONE
!
real(DP), intent(in) :: omega
integer, intent(in) :: intra_bgrp_comm
!
INTEGER :: ndm, startq, lastq, nt, l, nb, iq, ir
REAL(dp) :: qi
! q-point grid for interpolation
REAL(dp) :: pref
! the prefactor of the Q functions
real(DP) :: vqint, d1
!
real(DP), allocatable :: xdata(:)
! work space for spline
REAL(dp), allocatable :: aux (:)
! work space
REAL(dp), allocatable :: besr(:)
! work space
!
ndm = MAXVAL ( upf(:)%kkbeta )
allocate( aux (ndm) )
allocate (besr( ndm))
pref = fpi / sqrt (omega)
call divide (intra_bgrp_comm, nqx, startq, lastq)
tab (:,:,:) = 0.d0
do nt = 1, nsp
if ( upf(nt)%is_gth ) cycle
do nb = 1, upf(nt)%nbeta
l = upf(nt)%lll (nb)
do iq = startq, lastq
qi = (iq - 1) * dq
call sph_bes (upf(nt)%kkbeta, rgrid(nt)%r, qi, l, besr)
do ir = 1, upf(nt)%kkbeta
aux (ir) = upf(nt)%beta (ir, nb) * besr (ir) * rgrid(nt)%r(ir)
enddo
call simpson (upf(nt)%kkbeta, aux, rgrid(nt)%rab, vqint)
tab (iq, nb, nt) = vqint * pref
enddo
enddo
enddo
deallocate (besr)
deallocate (aux)
!
call mp_sum( tab, intra_bgrp_comm )
!
! initialize spline interpolation
!
if (spline_ps) then
allocate( xdata(nqx) )
do iq = 1, nqx
xdata(iq) = (iq - 1) * dq
enddo
do nt = 1, nsp
do nb = 1, upf(nt)%nbeta
d1 = (tab(2,nb,nt) - tab(1,nb,nt)) / dq
call spline(xdata, tab(:,nb,nt), 0.d0, d1, tab_d2y(:,nb,nt))
enddo
enddo
deallocate(xdata)
!
endif
!
! update GPU memory (taking care of zero-dim allocations)
!
#if defined __CUDA
if ( nbetam > 0 ) then
tab_d=tab
if (spline_ps) tab_d2y_d=tab_d2y
endif
#endif
!
END SUBROUTINE init_tab_beta

100
upflib/init_tab_qrad.f90 Normal file
View File

@ -0,0 +1,100 @@
!
! Copyright (C) 2021 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_qrad (omega, intra_bgrp_comm)
!----------------------------------------------------------------------
!
! Compute interpolation table qrad(i,nm,l+1,nt) = Q^{(L)}_{nm,nt}(q_i)
! of angular momentum L, for atom of type nt, on grid q_i, where
! nm = combined (n,m) index; n,m = 1,...,nbeta (number of beta functions)
!
USE upf_kinds, ONLY : dp
USE upf_const, ONLY : fpi
USE atom, ONLY : rgrid
USE uspp_param, ONLY : upf, lmaxq, nbetam, nsp
USE uspp_data, ONLY : nqxq, dq, qrad, qrad_d
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
real(DP), intent(in) :: omega
integer, intent(in) :: intra_bgrp_comm
!
INTEGER :: ndm, startq, lastq, nt, l, nb, mb, ijv, iq, ir
! various indices
REAL(dp) :: prefr
! the prefactor of the Q functions
REAL(dp) :: q
REAL(dp), ALLOCATABLE :: aux (:), besr (:)
! various work space
!
prefr = fpi / omega
ndm = MAXVAL ( upf(:)%kkbeta )
ALLOCATE (aux ( ndm))
ALLOCATE (besr( ndm))
!
CALL divide (intra_bgrp_comm, nqxq, startq, lastq)
!
qrad(:,:,:,:)= 0.d0
DO nt = 1, nsp
if ( upf(nt)%tvanp ) then
DO l = 0, upf(nt)%nqlc -1
!
! l is the true (combined) angular momentum
! Note that the index of array qfuncl runs from 0 to l,
! while the same index for qrad runs from 1 to l+1
! FIXME: qrad has "holes" if USPP/PAW do not precede NCPP
!
DO iq = startq, lastq
!
q = (iq - 1) * dq
!
! here we compute the spherical bessel function for each q_i
!
CALL sph_bes ( upf(nt)%kkbeta, rgrid(nt)%r, q, l, besr)
!
DO nb = 1, upf(nt)%nbeta
!
! the Q are symmetric with respect to nb,nm indices
!
DO mb = nb, upf(nt)%nbeta
ijv = mb * (mb - 1) / 2 + nb
IF ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .AND. &
( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .AND. &
(mod(l+upf(nt)%lll(nb)+upf(nt)%lll(mb),2)==0) ) THEN
DO ir = 1, upf(nt)%kkbeta
aux (ir) = besr (ir) * upf(nt)%qfuncl(ir,ijv,l)
ENDDO
!
! and then we integrate with all the Q functions
!
CALL simpson ( upf(nt)%kkbeta, aux, rgrid(nt)%rab, &
qrad(iq,ijv,l+1, nt) )
ENDIF
ENDDO
ENDDO
! igl
ENDDO
! l
ENDDO
qrad (:, :, :, nt) = qrad (:, :, :, nt)*prefr
CALL mp_sum ( qrad (:, :, :, nt), intra_bgrp_comm )
ENDIF
! nsp
ENDDO
!
DEALLOCATE (besr)
DEALLOCATE (aux)
!
! update GPU memory (taking care of zero-dim allocations)
!
#if defined __CUDA
if ( nbetam > 0 .and. lmaxq > 0 ) qrad_d=qrad
#endif
!
END SUBROUTINE init_tab_qrad

View File

@ -180,7 +180,7 @@ SUBROUTINE init_us_0(ecutrho,intra_bgrp_comm)
ENDDO
!
! 2) compute the fourier transform of the Qs and their integrated power spectum in
! reciprocal space - FIXME: use routine compute_qrad in init_us_1
! reciprocal space - FIXME: use routine init_tab_qrad
!
DO iq = startq, lastq
!

View File

@ -238,7 +238,7 @@ subroutine init_us_1( nat, ityp, omega, ngm, g, gg, intra_bgrp_comm )
! here for the US types we compute the Fourier transform of the
! Q functions.
!
IF ( lmaxq > 0 ) CALL compute_qrad(omega, intra_bgrp_comm)
IF ( lmaxq > 0 ) CALL init_tab_qrad(omega, intra_bgrp_comm)
!
! and finally we compute the qq coefficients by integrating the Q.
! The qq are the g=0 components of Q
@ -308,11 +308,11 @@ subroutine init_us_1( nat, ityp, omega, ngm, g, gg, intra_bgrp_comm )
!
! fill interpolation table tab (and tab_d2y for spline interpolation)
!
CALL compute_beta ( omega, intra_bgrp_comm )
CALL init_tab_beta ( omega, intra_bgrp_comm )
!
#if defined __CUDA
!
! update GPU memory (taking care of zero-dim allocations)ù
! update GPU memory (taking care of zero-dim allocations)
!
if (nhm>0) then
indv_d=indv
@ -338,184 +338,3 @@ subroutine init_us_1( nat, ityp, omega, ngm, g, gg, intra_bgrp_comm )
return
!
end subroutine init_us_1
!
!----------------------------------------------------------------------
SUBROUTINE compute_qrad (omega, intra_bgrp_comm)
!----------------------------------------------------------------------
!
! Compute interpolation table qrad(i,nm,l+1,nt) = Q^{(L)}_{nm,nt}(q_i)
! of angular momentum L, for atom of type nt, on grid q_i, where
! nm = combined (n,m) index; n,m = 1,...,nbeta (number of beta functions)
!
USE upf_kinds, ONLY : dp
USE upf_const, ONLY : fpi
USE atom, ONLY : rgrid
USE uspp_param, ONLY : upf, lmaxq, nbetam, nsp
USE uspp_data, ONLY : nqxq, dq, qrad, qrad_d
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
real(DP), intent(in) :: omega
integer, intent(in) :: intra_bgrp_comm
!
INTEGER :: ndm, startq, lastq, nt, l, nb, mb, ijv, iq, ir
! various indices
REAL(dp) :: prefr
! the prefactor of the Q functions
REAL(dp) :: q
REAL(dp), ALLOCATABLE :: aux (:), besr (:)
! various work space
!
prefr = fpi / omega
ndm = MAXVAL ( upf(:)%kkbeta )
ALLOCATE (aux ( ndm))
ALLOCATE (besr( ndm))
!
CALL divide (intra_bgrp_comm, nqxq, startq, lastq)
!
qrad(:,:,:,:)= 0.d0
DO nt = 1, nsp
if ( upf(nt)%tvanp ) then
DO l = 0, upf(nt)%nqlc -1
!
! l is the true (combined) angular momentum
! Note that the index of array qfuncl runs from 0 to l,
! while the same index for qrad runs from 1 to l+1
! FIXME: qrad has "holes" if USPP/PAW do not precede NCPP
!
DO iq = startq, lastq
!
q = (iq - 1) * dq
!
! here we compute the spherical bessel function for each q_i
!
CALL sph_bes ( upf(nt)%kkbeta, rgrid(nt)%r, q, l, besr)
!
DO nb = 1, upf(nt)%nbeta
!
! the Q are symmetric with respect to nb,nm indices
!
DO mb = nb, upf(nt)%nbeta
ijv = mb * (mb - 1) / 2 + nb
IF ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .AND. &
( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .AND. &
(mod(l+upf(nt)%lll(nb)+upf(nt)%lll(mb),2)==0) ) THEN
DO ir = 1, upf(nt)%kkbeta
aux (ir) = besr (ir) * upf(nt)%qfuncl(ir,ijv,l)
ENDDO
!
! and then we integrate with all the Q functions
!
CALL simpson ( upf(nt)%kkbeta, aux, rgrid(nt)%rab, &
qrad(iq,ijv,l+1, nt) )
ENDIF
ENDDO
ENDDO
! igl
ENDDO
! l
ENDDO
qrad (:, :, :, nt) = qrad (:, :, :, nt)*prefr
CALL mp_sum ( qrad (:, :, :, nt), intra_bgrp_comm )
ENDIF
! nsp
ENDDO
!
DEALLOCATE (besr)
DEALLOCATE (aux)
!
! update GPU memory (taking care of zero-dim allocations)
!
#if defined __CUDA
if ( nbetam > 0 .and. lmaxq > 0 ) qrad_d=qrad
#endif
!
END SUBROUTINE compute_qrad
!----------------------------------------------------------------------
SUBROUTINE compute_beta ( omega, intra_bgrp_comm )
!----------------------------------------------------------------------
!
! Compute interpolation table for beta(G) radial functions
!
USE upf_kinds, ONLY : dp
USE upf_const, ONLY : fpi
USE atom, ONLY : rgrid
USE uspp_param, ONLY : upf, lmaxq, nbetam, nsp
USE uspp_data, ONLY : nqx, dq, tab, tab_d2y, spline_ps, tab_d, tab_d2y_d
USE mp, ONLY : mp_sum
USE splinelib, ONLY : spline
!
IMPLICIT NONE
!
real(DP), intent(in) :: omega
integer, intent(in) :: intra_bgrp_comm
!
INTEGER :: ndm, startq, lastq, nt, l, nb, iq, ir
REAL(dp) :: qi
! q-point grid for interpolation
REAL(dp) :: pref
! the prefactor of the Q functions
real(DP) :: vqint, d1
!
real(DP), allocatable :: xdata(:)
! work space for spline
REAL(dp), allocatable :: aux (:)
! work space
REAL(dp), allocatable :: besr(:)
! work space
!
ndm = MAXVAL ( upf(:)%kkbeta )
allocate( aux (ndm) )
allocate (besr( ndm))
pref = fpi / sqrt (omega)
call divide (intra_bgrp_comm, nqx, startq, lastq)
tab (:,:,:) = 0.d0
do nt = 1, nsp
if ( upf(nt)%is_gth ) cycle
do nb = 1, upf(nt)%nbeta
l = upf(nt)%lll (nb)
do iq = startq, lastq
qi = (iq - 1) * dq
call sph_bes (upf(nt)%kkbeta, rgrid(nt)%r, qi, l, besr)
do ir = 1, upf(nt)%kkbeta
aux (ir) = upf(nt)%beta (ir, nb) * besr (ir) * rgrid(nt)%r(ir)
enddo
call simpson (upf(nt)%kkbeta, aux, rgrid(nt)%rab, vqint)
tab (iq, nb, nt) = vqint * pref
enddo
enddo
enddo
deallocate (besr)
deallocate (aux)
!
call mp_sum( tab, intra_bgrp_comm )
!
! initialize spline interpolation
!
if (spline_ps) then
allocate( xdata(nqx) )
do iq = 1, nqx
xdata(iq) = (iq - 1) * dq
enddo
do nt = 1, nsp
do nb = 1, upf(nt)%nbeta
d1 = (tab(2,nb,nt) - tab(1,nb,nt)) / dq
call spline(xdata, tab(:,nb,nt), 0.d0, d1, tab_d2y(:,nb,nt))
enddo
enddo
deallocate(xdata)
!
endif
!
! update GPU memory (taking care of zero-dim allocations)
!
#if defined __CUDA
if ( nbetam > 0 ) then
tab_d=tab
if (spline_ps) tab_d2y_d=tab_d2y
endif
#endif
!
END SUBROUTINE compute_beta

View File

@ -18,6 +18,19 @@ init_tab_atwfc.o : upf_const.o
init_tab_atwfc.o : upf_kinds.o
init_tab_atwfc.o : uspp.o
init_tab_atwfc.o : uspp_data.o
init_tab_beta.o : ../UtilXlib/mp.o
init_tab_beta.o : atom.o
init_tab_beta.o : splinelib.o
init_tab_beta.o : upf_const.o
init_tab_beta.o : upf_kinds.o
init_tab_beta.o : uspp.o
init_tab_beta.o : uspp_data.o
init_tab_qrad.o : ../UtilXlib/mp.o
init_tab_qrad.o : atom.o
init_tab_qrad.o : upf_const.o
init_tab_qrad.o : upf_kinds.o
init_tab_qrad.o : uspp.o
init_tab_qrad.o : uspp_data.o
init_us_0.o : ../UtilXlib/mp.o
init_us_0.o : atom.o
init_us_0.o : upf_const.o
@ -28,12 +41,10 @@ init_us_0.o : uspp_data.o
init_us_1.o : ../UtilXlib/mp.o
init_us_1.o : atom.o
init_us_1.o : paw_variables.o
init_us_1.o : splinelib.o
init_us_1.o : upf_const.o
init_us_1.o : upf_kinds.o
init_us_1.o : upf_spinorb.o
init_us_1.o : uspp.o
init_us_1.o : uspp_data.o
init_us_2_base.o : gth.o
init_us_2_base.o : splinelib.o
init_us_2_base.o : upf_const.o