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

Removal of duplicated PP code 

*** PLEASE TEST MACROSCOPIC ELECTRIC FIELDS AND BERRY PHASES WITH US-PP ***
This commit is contained in:
Paolo Giannozzi 2021-04-05 22:21:48 +02:00
parent 074fe831a3
commit 437f630e80
11 changed files with 33 additions and 319 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PW/src/Makefile
View File

@ -33,8 +33,6 @@ newd.o \
beef.o \
bp_mod.o \
bp_c_phase.o \
bp_calc_btq.o \
bp_qvan3.o \
bp_strings.o \
buffers.o \
c_bands.o \

19
PW/src/bp_c_phase.f90
View File

@ -9,6 +9,8 @@
! from c_phase_field.f90
! May 2012, A. Dal Corso: Noncollinear/spin-orbit case allowed (experimental).
! January 2019 Ronald E Cohen, fixed mods and averages on same branch
! April 2021 Paolo Giannozzi replaced calc_btq and qvan3 with "standard"
! QE way to compute Q(|G|) via interpolation table qrad + qvan2
!
!##############################################################################!
!# #!
@ -31,10 +33,10 @@
!# LIST OF FILES #!
!# ~~~~~~~~~~~~~ #!
!# The complete list of files added to the PWSCF distribution is: #!
!# * ../PW/bp_calc_btq.f90 #!
!# * ../PW/bp_c_phase.f90 #!
!# * ../PW/bp_qvan3.f90 #!
!# * ../PW/bp_strings.f90 #!
!# * ../PW/bp_calc_btq.f90 REMOVED - replaced by interpolation array qrad #!
!# * ../PW/bp_qvan3.f90 REMOVED - replaced by interpolation in qvan2 #!
!# #!
!# The PWSCF files that needed (minor) modifications were: #!
!# * ../PW/electrons.f90 #!
@ -241,6 +243,7 @@ SUBROUTINE c_phase
LOGICAL, ALLOCATABLE :: l_cal(:) ! flag for occupied/empty states
REAL(DP) :: t1,t !!REC
REAL(DP) :: dk(3)
REAL(DP) :: dk2
REAL(DP) :: dkmod
REAL(DP) :: el_loc
REAL(DP) :: eps
@ -251,7 +254,6 @@ SUBROUTINE c_phase
REAL(DP), ALLOCATABLE :: pdl_elec(:)
REAL(DP), ALLOCATABLE :: phik(:)
REAL(DP) :: phik_ave
REAL(DP) :: qrad_dk(nbetam,nbetam,lmaxq,ntyp)
REAL(DP) :: weight
REAL(DP) :: upol(3)
REAL(DP) :: pdl_elec_dw
@ -416,18 +418,19 @@ SUBROUTINE c_phase
! electronic polarization: form factor !
! ------------------------------------------------------------------------- !
if(okvan) then
! --- Calculate Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] ---
CALL calc_btq(dkmod,qrad_dk,0)
! --- Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] is in array qrad---
! CALL calc_btq(dkmod,qrad_dk,0) is no longer needed
! --- Calculate the q-space real spherical harmonics at dk [Y_LM] ---
dkmod=dk(1)**2+dk(2)**2+dk(3)**2
CALL ylmr2(lmaxq*lmaxq, 1, dk, dkmod, ylm_dk)
dk2=dk(1)**2+dk(2)**2+dk(3)**2
CALL ylmr2(lmaxq*lmaxq, 1, dk, dk2, ylm_dk)
! --- Form factor: 4 pi sum_LM c_ij^LM Y_LM(Omega) Q_ij^L(|r|) ---
q_dk = (0.d0, 0.d0)
DO np =1, ntyp
if( upf(np)%tvanp ) then
DO iv = 1, nh(np)
DO jv = iv, nh(np)
call qvan3(iv,jv,np,pref,ylm_dk,qrad_dk)
! call to qvan3 no longer needed
CALL qvan2(1,iv,jv,np,dkmod,pref,ylm_dk)
q_dk(iv,jv,np) = omega*pref
q_dk(jv,iv,np) = omega*pref
ENDDO

76
PW/src/bp_calc_btq.f90
View File

@ -1,76 +0,0 @@
!
! Copyright (C) 2004 Vanderbilt's group at Rutgers University, NJ
! 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 calc_btq( ql, qr_k, idbes )
!----------------------------------------------------------------------
!! Calculates the Bessel-transform (or its derivative if idbes=1)
!! of the augmented qrad charges at a given ql point.
!! Rydberg atomic units are used.
!
USE kinds, ONLY: DP
USE atom, ONLY: rgrid
USE ions_base, ONLY: ntyp => nsp
USE cell_base, ONLY: omega
USE constants, ONLY: fpi
USE uspp_param, ONLY: upf, nbetam, lmaxq
!
IMPLICIT NONE
!
REAL(DP) :: ql, qr_k(nbetam,nbetam,lmaxq,ntyp)
INTEGER :: idbes
!
INTEGER :: i, np, l, ilmin, ilmax, iv, jv, ijv, ilast
REAL(DP) :: qrk
REAL(DP), ALLOCATABLE :: jl(:), aux(:)
!
DO np = 1, ntyp
!
IF ( upf(np)%tvanp ) THEN
!
ALLOCATE( jl(upf(np)%kkbeta), aux(upf(np)%kkbeta) )
DO iv = 1, upf(np)%nbeta
DO jv = iv, upf(np)%nbeta
ijv = jv * (jv-1) / 2 + iv
ilmin = ABS( upf(np)%lll(iv) - upf(np)%lll(jv) )
ilmax = upf(np)%lll(iv) + upf(np)%lll(jv)
! only need to calculate for l=lmin,lmin+2 ...lmax-2,lmax
DO l = ilmin,ilmax,2
aux(:) = 0.0_DP
aux(1:upf(np)%kkbeta) = upf(np)%qfuncl(1:upf(np)%kkbeta,ijv,l)
IF (idbes == 1) THEN
!
CALL sph_dbes( upf(np)%kkbeta, rgrid(np)%r, ql, l, jl )
!
ELSE
!
CALL sph_bes( upf(np)%kkbeta, rgrid(np)%r, ql, l, jl )
!
ENDIF
!
! jl is the Bessel function (or its derivative) calculated at ql
! now integrate qfunc*jl*r^2 = Bessel transform of qfunc
!
DO i = 1, upf(np)%kkbeta
aux(i) = jl(i)*aux(i)
ENDDO
! if (tlog(np)) then
CALL simpson( upf(np)%kkbeta, aux, rgrid(np)%rab, qrk )
!
qr_k(iv,jv,l+1,np) = qrk * fpi / omega
qr_k(jv,iv,l+1,np) = qr_k(iv,jv,l+1,np)
!
ENDDO
ENDDO
ENDDO
DEALLOCATE( aux, jl )
ENDIF
ENDDO
!
RETURN
!
END SUBROUTINE calc_btq

98
PW/src/bp_qvan3.f90
View File

@ -1,98 +0,0 @@
!
! Copyright (C) 2004 Vanderbilt's group at Rutgers University, NJ
! 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 .
!
! Modified by PG - Oct.2007: removed obsolete comments
!--------------------------------------------------------------------------
SUBROUTINE qvan3( iv, jv, is, qg, ylm_k, qr )
!---------------------------------------------------------------------
!! It calculates:
!! \[ \text{qg} = \sum_{LM} (-I)^L \text{AP}(LM,\text{iv},
!! \text{jv}) \text{YR}_{LM} \text{QRAD}(\text{iv},
!! \text{jv},L,\text{is}) \]
!
! It calculates: qg = SUM_LM (-I)^L AP(LM,iv,jv) YR_LM QRAD(iv,jv,L,is)
!
USE kinds, ONLY: DP
USE ions_base, ONLY: ntyp => nsp
USE uspp_param, ONLY: lmaxq, nbetam
USE uspp, ONLY: nlx, lpl, lpx, ap, indv, nhtol, nhtolm
!
IMPLICIT NONE
!
INTEGER :: iv
!! beta function index
INTEGER :: jv
!! beta function index
INTEGER :: is
!! atomic type
COMPLEX(DP) :: qg
!! output: see routine comments
REAL(DP) :: ylm_k(lmaxq*lmaxq)
!! q-space real spherical harmonics at dk [\(Y_{LM}\)]
REAL(DP) :: qr(nbetam,nbetam,lmaxq,ntyp)
!! Bessel transform of \(Q_{ij}(|r|)\) at dk [\(Q_{ij}^L(|r|)\)]
!
! ... local variables
!
COMPLEX(DP) :: sig
INTEGER :: ivs, jvs, ivl, jvl, lp, l, i
!
!
ivs = indv(iv,is)
jvs = indv(jv,is)
ivl = nhtolm(iv,is)
jvl = nhtolm(jv,is)
!
IF (ivs > nbetam .OR. jvs > nbetam) &
CALL errore( ' qvan3 ', ' wrong dimensions (1)', MAX(ivs,jvs) )
IF (ivl > nlx .OR. jvl > nlx) &
CALL errore( ' qvan3 ', ' wrong dimensions (2)', MAX(ivl,jvl) )
!
qg = (0.0_DP,0.0_DP)
!
!odl Write(*,*) 'QVAN3 -- ivs jvs = ',ivs,jvs
!odl Write(*,*) 'QVAN3 -- ivl jvl = ',ivl,jvl
DO i = 1, lpx(ivl,jvl)
!odl Write(*,*) 'QVAN3 -- i = ',i
lp = lpl(ivl,jvl,i)
!odl Write(*,*) 'QVAN3 -- lp = ',lp
!
! ... EXTRACTION OF ANGULAR MOMENT L FROM LP:
!
IF (lp == 1) THEN
l = 1
ELSEIF ((lp >= 2).AND.(lp <= 4)) THEN
l = 2
ELSEIF ((lp >= 5).AND.(lp <= 9)) THEN
l = 3
ELSEIF ((lp >= 10).AND.(lp <= 16)) THEN
l = 4
ELSEIF ((lp >= 17).AND.(lp <= 25)) THEN
l = 5
ELSEIF ((lp >= 26).AND.(lp <= 36)) THEN
l = 6
ELSEIF ((lp >= 37).AND.(lp <= 49)) THEN
l = 7
ELSEIF (lp > 49) THEN
CALL errore( ' qvan3 ',' l not programmed ', lp )
ENDIF
!
sig = (0.0_DP,-1.0_DP)**(l-1)
sig = sig * ap(lp,ivl,jvl)
!
!odl Write(*,*) 'QVAN3 -- sig = ',sig
!
! WRITE( stdout,*) 'qvan3',ng1,LP,L,ivs,jvs
!
qg = qg + sig * ylm_k(lp) * qr(ivs,jvs,l,is)
!
ENDDO
!
!
RETURN
!
END SUBROUTINE qvan3

12
PW/src/c_phase_field.f90
View File

@ -97,6 +97,7 @@ SUBROUTINE c_phase_field( el_pola, ion_pola, fact_pola, pdir )
INTEGER :: nt
INTEGER :: nspinnc
REAL(DP) :: dk(3)
REAL(DP) :: dk2
REAL(DP) :: dkmod
REAL(DP) :: el_loc
REAL(DP) :: eps
@ -107,7 +108,6 @@ SUBROUTINE c_phase_field( el_pola, ion_pola, fact_pola, pdir )
REAL(DP), ALLOCATABLE :: loc_k(:)
REAL(DP), ALLOCATABLE :: pdl_elec(:)
REAL(DP), ALLOCATABLE :: phik(:)
REAL(DP) :: qrad_dk(nbetam,nbetam,lmaxq,ntyp)
REAL(DP) :: weight
REAL(DP) :: pola, pola_ion
REAL(DP), ALLOCATABLE :: wstring(:)
@ -306,12 +306,12 @@ SUBROUTINE c_phase_field( el_pola, ion_pola, fact_pola, pdir )
! electronic polarization: form factor !
! ------------------------------------------------------------------------- !
IF (okvan) THEN
! --- Calculate Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] ---
CALL calc_btq(dkmod,qrad_dk,0)
! --- Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] in array qrad ---
! CALL calc_btq(dkmod,qrad_dk,0) is no longer needed, see bp_c_phase
!
! --- Calculate the q-space real spherical harmonics at dk [Y_LM] ---
dkmod = dk(1)**2+dk(2)**2+dk(3)**2
CALL ylmr2(lmaxq*lmaxq, 1, dk, dkmod, ylm_dk)
dk2 = dk(1)**2+dk(2)**2+dk(3)**2
CALL ylmr2(lmaxq*lmaxq, 1, dk, dk2, ylm_dk)
!
! --- Form factor: 4 pi sum_LM c_ij^LM Y_LM(Omega) Q_ij^L(|r|) ---
q_dk=(0.d0,0.d0)
@ -319,7 +319,7 @@ SUBROUTINE c_phase_field( el_pola, ion_pola, fact_pola, pdir )
IF ( upf(np)%tvanp ) THEN
DO iv = 1, nh(np)
DO jv = iv, nh(np)
CALL qvan3( iv,jv,np,pref,ylm_dk,qrad_dk )
CALL qvan2(1,iv,jv,np,dkmod,pref,ylm_dk)
q_dk(iv,jv,np) = omega*pref
q_dk(jv,iv,np) = omega*pref
ENDDO

12
PW/src/forces_bp_efield.f90
View File

@ -125,6 +125,7 @@ SUBROUTINE forces_us_efield( forces_bp, pdir, e_field )
INTEGER :: nstring
INTEGER :: nt
REAL(dp) :: dk(3)
REAL(dp) :: dk2
REAL(dp) :: dkmod
REAL(dp) :: el_loc
REAL(dp) :: eps
@ -135,7 +136,6 @@ SUBROUTINE forces_us_efield( forces_bp, pdir, e_field )
REAL(dp), ALLOCATABLE :: loc_k(:)
REAL(dp), ALLOCATABLE :: pdl_elec(:)
REAL(dp), ALLOCATABLE :: phik(:)
REAL(dp) :: qrad_dk(nbetam,nbetam,lmaxq,ntyp)
REAL(dp) :: weight
REAL(dp) :: pola, pola_ion
REAL(dp), ALLOCATABLE :: wstring(:)
@ -340,12 +340,12 @@ SUBROUTINE forces_us_efield( forces_bp, pdir, e_field )
! electronic polarization: form factor !
! ------------------------------------------------------------------------- !
IF (okvan) THEN
! --- Calculate Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] ---
CALL calc_btq( dkmod, qrad_dk, 0 )
! --- Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] in array qrad ---
! CALL calc_btq( dkmod, qrad_dk, 0 ) no longer needed
!
! --- Calculate the q-space real spherical harmonics at dk [Y_LM] ---
dkmod = dk(1)**2+dk(2)**2+dk(3)**2
CALL ylmr2( lmaxq*lmaxq, 1, dk, dkmod, ylm_dk )
dk2 = dk(1)**2+dk(2)**2+dk(3)**2
CALL ylmr2( lmaxq*lmaxq, 1, dk, dk2, ylm_dk )
!
! --- Form factor: 4 pi sum_LM c_ij^LM Y_LM(Omega) Q_ij^L(|r|) ---
q_dk = (0.d0,0.d0)
@ -353,7 +353,7 @@ SUBROUTINE forces_us_efield( forces_bp, pdir, e_field )
IF ( upf(np)%tvanp ) THEN
DO iv = 1, nh(np)
DO jv = iv, nh(np)
CALL qvan3( iv, jv, np, pref, ylm_dk, qrad_dk )
CALL qvan2( 1, iv, jv, np, pref, ylm_dk )
q_dk(iv,jv,np) = omega*pref
q_dk(jv,iv,np) = omega*pref
ENDDO

18
PW/src/h_epsi_her_set.f90
View File

@ -85,6 +85,7 @@ SUBROUTINE h_epsi_her_set( pdir, e_field )
INTEGER :: ik_stringa!k-point index inside string
REAL(DP) :: dk(3)
REAL(DP) :: dkm(3)! -dk
REAL(DP) :: dk2
REAL(DP) :: dkmod
REAL(DP) :: eps
REAL(DP) :: fac
@ -93,7 +94,6 @@ SUBROUTINE h_epsi_her_set( pdir, e_field )
!REAL(DP) :: gvec
REAL(DP), ALLOCATABLE :: ln(:,:,:)
REAL(DP), ALLOCATABLE :: ln0(:,:,:)!map g-space global to g-space k-point dependent
REAL(DP) :: qrad_dk(nbetam,nbetam,lmaxq,ntyp)
REAL(DP) :: ylm_dk(lmaxq*lmaxq)
COMPLEX(DP), ALLOCATABLE :: aux(:)
COMPLEX(DP), ALLOCATABLE :: aux0(:)
@ -278,12 +278,12 @@ SUBROUTINE h_epsi_her_set( pdir, e_field )
! electronic polarization: form factor !
!-------------------------------------------------------------------------!
! --- Calculate Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] ---
CALL calc_btq( dkmod, qrad_dk, 0 )
! --- Bessel transform of Q_ij(|r|) at dk [Q_ij^L(|r|)] in array qrad ---
! CALL calc_btq( dkmod, qrad_dk, 0 ) no longer needed
!
! --- Calculate the q-space real spherical harmonics at dk [Y_LM] ---
dkmod = dk(1)**2 + dk(2)**2 + dk(3)**2
CALL ylmr2( lmaxq*lmaxq, 1, dk, dkmod, ylm_dk )
dk2 = dk(1)**2 + dk(2)**2 + dk(3)**2
CALL ylmr2( lmaxq*lmaxq, 1, dk, dk2, ylm_dk )
!
! --- Form factor: 4 pi sum_LM c_ij^LM Y_LM(Omega) Q_ij^L(|r|) ---
q_dk = (0.d0,0.d0)
@ -291,7 +291,7 @@ SUBROUTINE h_epsi_her_set( pdir, e_field )
IF ( upf(np)%tvanp ) THEN
DO iv = 1, nh(np)
DO jv = iv, nh(np)
CALL qvan3( iv, jv, np, pref, ylm_dk, qrad_dk )
CALL qvan2( 1, iv, jv, np, dkmod, pref, ylm_dk )
q_dk(iv,jv,np) = omega*pref
q_dk(jv,iv,np) = omega*pref
ENDDO
@ -302,8 +302,8 @@ SUBROUTINE h_epsi_her_set( pdir, e_field )
!
! --- Calculate the q-space real spherical harmonics at -dk [Y_LM] ---
!
dkmod = dkm(1)**2 + dkm(2)**2 + dkm(3)**2
CALL ylmr2( lmaxq*lmaxq, 1, dkm, dkmod, ylm_dk )
dk2 = dkm(1)**2 + dkm(2)**2 + dkm(3)**2
CALL ylmr2( lmaxq*lmaxq, 1, dkm, dk2, ylm_dk )
!
! --- Form factor: 4 pi sum_LM c_ij^LM Y_LM(Omega) Q_ij^L(|r|) ---
!
@ -312,7 +312,7 @@ SUBROUTINE h_epsi_her_set( pdir, e_field )
IF ( upf(np)%tvanp ) THEN
DO iv = 1, nh(np)
DO jv = iv, nh(np)
CALL qvan3( iv, jv, np, pref, ylm_dk, qrad_dk )
CALL qvan2( 1, iv, jv, np, dkmod, pref, ylm_dk )
q_dkp(iv,jv,np) = omega*pref
q_dkp(jv,iv,np) = omega*pref
ENDDO

9
PW/src/make.depend
View File

@ -270,12 +270,6 @@ bp_c_phase.o : ../../upflib/uspp.o
bp_c_phase.o : bp_mod.o
bp_c_phase.o : buffers.o
bp_c_phase.o : pwcom.o
bp_calc_btq.o : ../../Modules/cell_base.o
bp_calc_btq.o : ../../Modules/constants.o
bp_calc_btq.o : ../../Modules/ions_base.o
bp_calc_btq.o : ../../Modules/kind.o
bp_calc_btq.o : ../../upflib/atom.o
bp_calc_btq.o : ../../upflib/uspp.o
bp_mod.o : ../../Modules/becmod.o
bp_mod.o : ../../Modules/cell_base.o
bp_mod.o : ../../Modules/fft_base.o
@ -283,9 +277,6 @@ bp_mod.o : ../../Modules/kind.o
bp_mod.o : ../../Modules/mp_images.o
bp_mod.o : ../../Modules/recvec.o
bp_mod.o : ../../UtilXlib/mp.o
bp_qvan3.o : ../../Modules/ions_base.o
bp_qvan3.o : ../../Modules/kind.o
bp_qvan3.o : ../../upflib/uspp.o
bp_strings.o : ../../Modules/kind.o
buffers.o : ../../Modules/io_files.o
buffers.o : ../../Modules/kind.o

103
upflib/sph_bes.f90
View File

@ -296,106 +296,3 @@ SUBROUTINE sph_dbes1 ( nr, r, xg, l, jl, djl )
end if
!
end SUBROUTINE sph_dbes1
!
!----------------------------------------------------------------------------
SUBROUTINE sph_dbes( MMAX, R, XG, L, DJL )
!----------------------------------------------------------------------------
!! Calculates derivatives of spherical Bessel functions \(j_l(Gr)\)
!! \(-x D(j_l(x))/dx \) - Called in bp_calc_btq but not actually used
!! FIXME: MERGE WITH sph_dbes1
!
USE upf_kinds, only: DP
USE upf_const, ONLY : eps8
!
IMPLICIT NONE
!
INTEGER :: MMAX, L
REAL(DP) :: XG
REAL(DP) :: DJL(MMAX), R(MMAX)
!
INTEGER :: IR
REAL(DP) :: XRG, XRG2
!
!
IF ( L == 1 ) THEN ! S PART
IF( XG < eps8 ) THEN
DO IR=1,MMAX
DJL(IR) = 0.D0
END DO
ELSE
DJL(1) = 0.D0
DO IR=2,MMAX
XRG=R(IR)*XG
DJL(IR) = SIN(XRG)/XRG-COS(XRG)
END DO
ENDIF
ENDIF
!
IF ( L == 2 ) THEN ! P PART
IF( XG < eps8 ) THEN
DO IR=1,MMAX
DJL(IR) = 0.D0
END DO
ELSE
DJL(1) = 0.D0
DO IR=2,MMAX
XRG=R(IR)*XG
DJL(IR) = 2.D0*(SIN(XRG)/XRG-COS(XRG))/XRG - SIN(XRG)
END DO
ENDIF
ENDIF
!
IF ( L == 3 ) THEN ! D PART
IF ( XG < eps8 ) THEN
DO IR=1,MMAX
DJL(IR) = 0.D0
END DO
ELSE
DJL(1) = 0.D0
DO IR=2,MMAX
XRG=R(IR)*XG
DJL(IR) = ( SIN(XRG)*(9.D0/(XRG*XRG)-4.D0) - &
9.D0*COS(XRG)/XRG ) /XRG + COS(XRG)
END DO
END IF
END IF
!
IF ( L == 4 ) THEN ! F PART
IF ( XG < eps8 ) THEN
DO IR=1,MMAX
DJL(IR) = 0.D0
END DO
ELSE
DJL(1) = 0.D0
DO IR=2,MMAX
XRG=R(IR)*XG
XRG2=XRG*XRG
DJL(IR) = SIN(XRG)*(60.D0/(XRG2*XRG2)-27.D0/XRG2+1.d0) - &
COS(XRG)*(60.D0/XRG2-7.D0)/XRG
END DO
END IF
END IF
!
IF ( L == 5 ) THEN ! G PART
IF ( XG < eps8 ) THEN
DO IR=1,MMAX
DJL(IR) = 0.D0
END DO
ELSE
DJL(1) = 0.D0
DO IR=2,MMAX
XRG=R(IR)*XG
XRG2=XRG*XRG
DJL(IR) = SIN(XRG)*(525.D0/(XRG2*XRG2)-240.D0/XRG2+11.D0)/XRG - &
COS(XRG)*(525.D0/(XRG2*XRG2)-65.D0/XRG2+1.D0)
END DO
END IF
END IF
!
IF ( L <= 0 .OR. L >= 6 ) &
CALL errore( 'sph_dbes', ' L NOT PROGRAMMED, L= ',L )
!
RETURN
!
END SUBROUTINE sph_dbes

1
upflib/upf_to_internal.f90
View File

@ -68,7 +68,6 @@ SUBROUTINE set_upf_q (upf)
! cases. This is already the case when upf%tpawp or upf%q_with_l are .true.
! For vanderbilt US pseudos, where nqf and rinner are non zero, we do here
! what otherwise would be done multiple times in many parts of the code
! (such as in init_us_1, addusforce_r, bp_calc_btq, compute_qdipol)
! whenever the q_l(r) were to be constructed.
! For simple rrkj3 pseudos we duplicate the information contained in q(r)
! for all q_l(r).

2
upflib/ylmr2.f90
View File

@ -24,7 +24,7 @@ subroutine ylmr2 (lmax2, ng, g, gg, ylm)
real(DP), intent(in) :: g (3, ng), gg (ng)
!
! BEWARE: gg = g(1)^2 + g(2)^2 +g(3)^2 is not checked on input
! incorrect results will ensue if the above does not hold
! Incorrect results will ensue if gg != g(1)^2 + g(2)^2 +g(3)^2
!
real(DP), intent(out) :: ylm (ng,lmax2)
!