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

Splines are disabled by default, but now they are 

compiled in. Set spline_ps = .true. in the &system namelist.
(D.C.)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3717 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
ceresoli 2007-01-25 11:42:52 +00:00
parent eae6a402ca
commit 21f7d5ebca
16 changed files with 196 additions and 268 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
GIPAW/gipaw_module.f90
View File

@ -103,13 +103,14 @@ CONTAINS
SUBROUTINE gipaw_readin() SUBROUTINE gipaw_readin()
USE io_files, ONLY : nd_nmbr, prefix, tmp_dir USE io_files, ONLY : nd_nmbr, prefix, tmp_dir
USE io_global, ONLY : ionode USE io_global, ONLY : ionode
USE us, ONLY : spline_ps
IMPLICIT NONE IMPLICIT NONE
INTEGER :: ios INTEGER :: ios
NAMELIST /inputgipaw/ job, prefix, tmp_dir, conv_threshold, & NAMELIST /inputgipaw/ job, prefix, tmp_dir, conv_threshold, &
q_gipaw, iverbosity, filcurr, filfield, & q_gipaw, iverbosity, filcurr, filfield, &
read_recon_in_paratec_fmt, & read_recon_in_paratec_fmt, &
file_reconstruction, use_nmr_macroscopic_shape, & file_reconstruction, use_nmr_macroscopic_shape, &
nmr_macroscopic_shape nmr_macroscopic_shape, spline_ps
if ( .not. ionode ) goto 400 if ( .not. ionode ) goto 400
@ -124,6 +125,7 @@ CONTAINS
read_recon_in_paratec_fmt = .FALSE. read_recon_in_paratec_fmt = .FALSE.
file_reconstruction ( : ) = " " file_reconstruction ( : ) = " "
nmr_macroscopic_shape = 2.0_dp / 3.0_dp nmr_macroscopic_shape = 2.0_dp / 3.0_dp
spline_ps = .true. ! TRUE in this case!!!!!
read( 5, inputgipaw, err = 200, iostat = ios ) read( 5, inputgipaw, err = 200, iostat = ios )

60
GIPAW/init_paw_2_no_phase.f90
View File

@ -21,15 +21,10 @@ subroutine init_paw_2_no_phase (npw_, igk_, q_, vkb_)
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau
USE gvect , ONLY : eigts1, eigts2, eigts3, g, ig1, ig2, ig3 USE gvect , ONLY : eigts1, eigts2, eigts3, g, ig1, ig2, ig3
USE us, ONLY : dq USE us, ONLY : dq
#ifdef USE_SPLINES
USE paw, ONLY : paw_nkb, paw_lmaxkb, paw_nhm, paw_nh, paw_nhtol, & USE paw, ONLY : paw_nkb, paw_lmaxkb, paw_nhm, paw_nh, paw_nhtol, &
paw_nhtom, paw_indv, paw_nbeta, paw_tab, paw_tab_d2y paw_nhtom, paw_indv, paw_nbeta, paw_tab, paw_tab_d2y
USE us, ONLY : nqxq, dq USE us, ONLY : nqxq, dq, spline_ps
USE splinelib USE splinelib
#else
USE paw, ONLY : paw_nkb, paw_lmaxkb, paw_nhm, paw_nh, paw_nhtol, &
paw_nhtom, paw_indv, paw_tab, paw_nbeta
#endif
! !
implicit none implicit none
! !
@ -51,10 +46,8 @@ subroutine init_paw_2_no_phase (npw_, igk_, q_, vkb_)
complex(DP) :: phase, pref complex(DP) :: phase, pref
complex(DP), allocatable :: sk(:) complex(DP), allocatable :: sk(:)
#ifdef USE_SPLINES
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
integer :: startq, lastq, iq integer :: iq
#endif
! !
! !
if (paw_lmaxkb.lt.0) return if (paw_lmaxkb.lt.0) return
@ -82,38 +75,35 @@ subroutine init_paw_2_no_phase (npw_, igk_, q_, vkb_)
qg(ig) = sqrt(qg(ig))*tpiba qg(ig) = sqrt(qg(ig))*tpiba
enddo enddo
#ifdef USE_SPLINES if (spline_ps) then
startq = 1 allocate(xdata(nqxq))
lastq = nqxq do iq = 1, nqxq
!call divide (nqxq, startq, lastq) xdata(iq) = (iq - 1) * dq
allocate(xdata(lastq-startq+1)) enddo
do iq = startq, lastq endif
xdata(iq) = (iq - 1) * dq
enddo
#endif
jkb = 0 jkb = 0
do nt = 1, ntyp do nt = 1, ntyp
! calculate beta in G-space using an interpolation table ! calculate beta in G-space using an interpolation table
do nb = 1, paw_nbeta (nt) do nb = 1, paw_nbeta (nt)
do ig = 1, npw_ do ig = 1, npw_
#ifdef USE_SPLINES if (spline_ps) then
vq(ig) = splint(xdata, paw_tab(:,nb,nt), paw_tab_d2y(:,nb,nt), & vq(ig) = splint(xdata, paw_tab(:,nb,nt), &
qg(ig)) paw_tab_d2y(:,nb,nt), qg(ig))
#else else
px = qg (ig) / dq - int (qg (ig) / dq) px = qg (ig) / dq - int (qg (ig) / dq)
ux = 1.d0 - px ux = 1.d0 - px
vx = 2.d0 - px vx = 2.d0 - px
wx = 3.d0 - px wx = 3.d0 - px
i0 = qg (ig) / dq + 1 i0 = qg (ig) / dq + 1
i1 = i0 + 1 i1 = i0 + 1
i2 = i0 + 2 i2 = i0 + 2
i3 = i0 + 3 i3 = i0 + 3
vq (ig) = paw_tab (i0, nb, nt) * ux * vx * wx / 6.d0 + & vq (ig) = paw_tab (i0, nb, nt) * ux * vx * wx / 6.d0 + &
paw_tab (i1, nb, nt) * px * vx * wx / 2.d0 - & paw_tab (i1, nb, nt) * px * vx * wx / 2.d0 - &
paw_tab (i2, nb, nt) * px * ux * wx / 2.d0 + & paw_tab (i2, nb, nt) * px * ux * wx / 2.d0 + &
paw_tab (i3, nb, nt) * px * ux * vx / 6.d0 paw_tab (i3, nb, nt) * px * ux * vx / 6.d0
#endif endif
enddo enddo
! add spherical harmonic part ! add spherical harmonic part
do ih = 1, paw_nh (nt) do ih = 1, paw_nh (nt)

63
GIPAW/init_us_2_no_phase.f90
View File

@ -20,12 +20,8 @@ subroutine init_us_2_no_phase (npw_, igk_, q_, vkb_)
USE constants, ONLY : tpi USE constants, ONLY : tpi
USE gvect, ONLY : eigts1, eigts2, eigts3, ig1, ig2, ig3, g USE gvect, ONLY : eigts1, eigts2, eigts3, ig1, ig2, ig3, g
USE wvfct, ONLY : npw, npwx, igk USE wvfct, ONLY : npw, npwx, igk
#ifdef USE_SPLINES USE us, ONLY : nqxq, dq, tab, tab_d2y, spline_ps
USE us, ONLY : nqxq, dq, tab, tab_d2y
USE splinelib USE splinelib
#else
USE us, ONLY : dq, tab
#endif
USE uspp, ONLY : nkb, vkb, nhtol, nhtolm, indv USE uspp, ONLY : nkb, vkb, nhtol, nhtolm, indv
USE uspp_param, ONLY : lmaxkb, nbeta, nhm, nh USE uspp_param, ONLY : lmaxkb, nbeta, nhm, nh
! !
@ -49,10 +45,8 @@ subroutine init_us_2_no_phase (npw_, igk_, q_, vkb_)
complex(DP) :: phase, pref complex(DP) :: phase, pref
complex(DP), allocatable :: sk(:) complex(DP), allocatable :: sk(:)
#ifdef USE_SPLINES
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
integer :: startq, lastq, iq integer :: iq
#endif
! !
! !
@ -80,42 +74,39 @@ subroutine init_us_2_no_phase (npw_, igk_, q_, vkb_)
qg(ig) = sqrt(qg(ig))*tpiba qg(ig) = sqrt(qg(ig))*tpiba
enddo enddo
#ifdef USE_SPLINES if (spline_ps) then
startq = 1 allocate(xdata(nqxq))
lastq = nqxq do iq = 1, nqxq
!call divide (nqxq, startq, lastq) xdata(iq) = (iq - 1) * dq
allocate(xdata(lastq-startq+1)) enddo
do iq = startq, lastq endif
xdata(iq) = (iq - 1) * dq
enddo
#endif
jkb = 0 jkb = 0
do nt = 1, ntyp do nt = 1, ntyp
! calculate beta in G-space using an interpolation table ! calculate beta in G-space using an interpolation table
do nb = 1, nbeta (nt) do nb = 1, nbeta (nt)
do ig = 1, npw_ do ig = 1, npw_
#ifdef USE_SPLINES if (spline_ps) then
vq(ig) = splint(xdata, tab(:,nb,nt), tab_d2y(:,nb,nt), qg(ig)) vq(ig) = splint(xdata, tab(:,nb,nt), tab_d2y(:,nb,nt), qg(ig))
#else
px = qg (ig) / dq - int (qg (ig) / dq)
ux = 1.d0 - px
vx = 2.d0 - px
wx = 3.d0 - px
i0 = INT( qg (ig) / dq ) + 1
i1 = i0 + 1
i2 = i0 + 2
i3 = i0 + 3
!*apsi TMPTMPTMP
if ( i3 > size(tab,1) ) then
vq(ig) = 0.0_dp
else else
vq (ig) = tab (i0, nb, nt) * ux * vx * wx / 6.d0 + & px = qg (ig) / dq - int (qg (ig) / dq)
tab (i1, nb, nt) * px * vx * wx / 2.d0 - & ux = 1.d0 - px
tab (i2, nb, nt) * px * ux * wx / 2.d0 + & vx = 2.d0 - px
tab (i3, nb, nt) * px * ux * vx / 6.d0 wx = 3.d0 - px
i0 = INT( qg (ig) / dq ) + 1
i1 = i0 + 1
i2 = i0 + 2
i3 = i0 + 3
!*apsi TMPTMPTMP
if ( i3 > size(tab,1) ) then
vq(ig) = 0.0_dp
else
vq (ig) = tab (i0, nb, nt) * ux * vx * wx / 6.d0 + &
tab (i1, nb, nt) * px * vx * wx / 2.d0 - &
tab (i2, nb, nt) * px * ux * wx / 2.d0 + &
tab (i3, nb, nt) * px * ux * vx / 6.d0
endif
endif endif
#endif
enddo enddo
! add spherical harmonic part ! add spherical harmonic part
do ih = 1, nh (nt) do ih = 1, nh (nt)

6
Modules/input_parameters.f90
View File

@ -567,6 +567,9 @@ MODULE input_parameters
LOGICAL :: assume_molsys = .FALSE. LOGICAL :: assume_molsys = .FALSE.
! use spline interpolation for pseudopotential
LOGICAL :: spline_ps = .false.
NAMELIST / system / ibrav, celldm, a, b, c, cosab, cosac, cosbc, nat, & NAMELIST / system / ibrav, celldm, a, b, c, cosab, cosac, cosbc, nat, &
ntyp, nbnd, nelec, ecutwfc, ecutrho, nr1, nr2, nr3, nr1s, nr2s, & ntyp, nbnd, nelec, ecutwfc, ecutrho, nr1, nr2, nr3, nr1s, nr2s, &
nr3s, nr1b, nr2b, nr3b, nosym, starting_magnetization, & nr3s, nr1b, nr2b, nr3b, nosym, starting_magnetization, &
@ -580,7 +583,8 @@ MODULE input_parameters
noncolin, lspinorb, lambda, angle1, angle2, report, & noncolin, lspinorb, lambda, angle1, angle2, report, &
constrained_magnetization, B_field, fixed_magnetization, & constrained_magnetization, B_field, fixed_magnetization, &
sic, sic_epsilon, force_pairing, sic_alpha, & sic, sic_epsilon, force_pairing, sic_alpha, &
tot_charge, multiplicity, tot_magnetization tot_charge, multiplicity, tot_magnetization, &
spline_ps
!=----------------------------------------------------------------------------=! !=----------------------------------------------------------------------------=!
! ELECTRONS Namelist Input Parameters ! ELECTRONS Namelist Input Parameters

2
Modules/read_namelists.f90
View File

@ -201,6 +201,7 @@ MODULE read_namelists_module
report = 1 report = 1
! !
assume_molsys = .FALSE. assume_molsys = .FALSE.
spline_ps = .false.
! !
RETURN RETURN
! !
@ -739,6 +740,7 @@ MODULE read_namelists_module
CALL mp_bcast( lambda, ionode_id ) CALL mp_bcast( lambda, ionode_id )
! !
CALL mp_bcast( assume_molsys, ionode_id ) CALL mp_bcast( assume_molsys, ionode_id )
CALL mp_bcast( spline_ps, ionode_id )
! !
RETURN RETURN
! !

11
PW/allocate_nlpot.f90
View File

@ -35,11 +35,8 @@ subroutine allocate_nlpot
USE ldaU, ONLY : Hubbard_lmax, ns, nsnew USE ldaU, ONLY : Hubbard_lmax, ns, nsnew
USE noncollin_module, ONLY : noncolin USE noncollin_module, ONLY : noncolin
USE wvfct, ONLY : npwx, npw, igk, g2kin USE wvfct, ONLY : npwx, npw, igk, g2kin
#ifdef USE_SPLINES USE us, ONLY : qrad, tab, tab_d2y, tab_at, dq, nqx, &
USE us, ONLY : qrad, tab, tab_d2y, tab_at, dq, nqx, nqxq nqxq, spline_ps
#else
USE us, ONLY : qrad, tab, tab_at, dq, nqx, nqxq
#endif
USE uspp, ONLY : indv, nhtol, nhtolm, qq, dvan, deeq, vkb, nkb, & USE uspp, ONLY : indv, nhtol, nhtolm, qq, dvan, deeq, vkb, nkb, &
nkbus, nhtoj, becsum, qq_so, dvan_so, deeq_nc nkbus, nhtoj, becsum, qq_so, dvan_so, deeq_nc
USE uspp_param, ONLY : lmaxq, lmaxkb, lll, nbeta, nh, nhm, nbetam,tvanp USE uspp_param, ONLY : lmaxq, lmaxkb, lll, nbeta, nh, nhm, nbetam,tvanp
@ -134,10 +131,8 @@ subroutine allocate_nlpot
allocate (tab( nqx , nbetam , nsp)) allocate (tab( nqx , nbetam , nsp))
#ifdef USE_SPLINES
! d2y is for the cubic splines ! d2y is for the cubic splines
allocate (tab_d2y( nqx , nbetam , nsp)) if (spline_ps) allocate (tab_d2y( nqx , nbetam , nsp))
#endif
allocate (tab_at( nqx , nchix , nsp)) allocate (tab_at( nqx , nchix , nsp))

12
PW/clean_pw.f90
View File

@ -26,11 +26,7 @@ SUBROUTINE clean_pw( lflag )
USE scf, ONLY : rho, rhog, vr, vltot, & USE scf, ONLY : rho, rhog, vr, vltot, &
rho_core, rhog_core,vrs rho_core, rhog_core,vrs
USE wavefunctions_module, ONLY : evc, psic, psic_nc USE wavefunctions_module, ONLY : evc, psic, psic_nc
#ifdef USE_SPLINES USE us, ONLY : qrad, tab, tab_at, tab_d2y, spline_ps
USE us, ONLY : qrad, tab, tab_at, tab_d2y
#else
USE us, ONLY : qrad, tab, tab_at
#endif
USE uspp, ONLY : deallocate_uspp USE uspp, ONLY : deallocate_uspp
USE ldaU, ONLY : ns, nsnew, swfcatom USE ldaU, ONLY : ns, nsnew, swfcatom
USE extfield, ONLY : forcefield USE extfield, ONLY : forcefield
@ -88,9 +84,9 @@ SUBROUTINE clean_pw( lflag )
IF ( ALLOCATED( psic ) ) DEALLOCATE( psic ) IF ( ALLOCATED( psic ) ) DEALLOCATE( psic )
IF ( ALLOCATED( psic_nc ) ) DEALLOCATE( psic_nc ) IF ( ALLOCATED( psic_nc ) ) DEALLOCATE( psic_nc )
IF ( ALLOCATED( vrs ) ) DEALLOCATE( vrs ) IF ( ALLOCATED( vrs ) ) DEALLOCATE( vrs )
#ifdef USE_SPLINES if (spline_ps) then
IF ( ALLOCATED( tab_d2y) ) DEALLOCATE( tab_d2y ) IF ( ALLOCATED( tab_d2y) ) DEALLOCATE( tab_d2y )
#endif endif
IF ( doublegrid ) THEN IF ( doublegrid ) THEN
IF ( ASSOCIATED( nls ) ) DEALLOCATE( nls ) IF ( ASSOCIATED( nls ) ) DEALLOCATE( nls )
END IF END IF

59
PW/gen_us_dj.f90
View File

@ -22,12 +22,8 @@ subroutine gen_us_dj (ik, dvkb)
USE gvect, ONLY : ig1, ig2, ig3, eigts1, eigts2, eigts3, g USE gvect, ONLY : ig1, ig2, ig3, eigts1, eigts2, eigts3, g
USE wvfct, ONLY : npw, npwx, igk USE wvfct, ONLY : npw, npwx, igk
USE uspp, ONLY : nkb, indv, nhtol, nhtolm USE uspp, ONLY : nkb, indv, nhtol, nhtolm
#ifdef USE_SPLINES USE us, ONLY : nqxq, tab, tab_d2y, dq, spline_ps
USE us, ONLY : nqxq, tab, tab_d2y, dq
USE splinelib USE splinelib
#else
USE us, ONLY : tab, dq
#endif
USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh
! !
implicit none implicit none
@ -59,10 +55,8 @@ subroutine gen_us_dj (ik, dvkb)
complex(DP), allocatable :: sk (:) complex(DP), allocatable :: sk (:)
#ifdef USE_SPLINES integer :: iq
integer :: startq, lastq, iq
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
#endif
if (nkb.eq.0) return if (nkb.eq.0) return
@ -85,36 +79,35 @@ subroutine gen_us_dj (ik, dvkb)
call stop_clock('stres_us32') call stop_clock('stres_us32')
call start_clock('stres_us33') call start_clock('stres_us33')
#ifdef USE_SPLINES
startq = 1 if (spline_ps) then
lastq = nqxq allocate(xdata(nqxq))
!call divide (nqxq, startq, lastq) do iq = 1, nqxq
allocate(xdata(lastq-startq+1)) xdata(iq) = (iq - 1) * dq
do iq = startq, lastq enddo
xdata(iq) = (iq - 1) * dq endif
enddo
#endif
do nt = 1, ntyp do nt = 1, ntyp
do nb = 1, nbeta (nt) do nb = 1, nbeta (nt)
do ig = 1, npw do ig = 1, npw
qt = sqrt(q (ig)) * tpiba qt = sqrt(q (ig)) * tpiba
#ifdef USE_SPLINES if (spline_ps) then
djl(ig,nb,nt) = splint_deriv(xdata, tab(:,nb,nt), tab_d2y(:,nb,nt), qt) djl(ig,nb,nt) = splint_deriv(xdata, tab(:,nb,nt), &
#else tab_d2y(:,nb,nt), qt)
px = qt / dq - int (qt / dq) else
ux = 1.d0 - px px = qt / dq - int (qt / dq)
vx = 2.d0 - px ux = 1.d0 - px
wx = 3.d0 - px vx = 2.d0 - px
i0 = qt / dq + 1 wx = 3.d0 - px
i1 = i0 + 1 i0 = qt / dq + 1
i2 = i0 + 2 i1 = i0 + 1
i3 = i0 + 3 i2 = i0 + 2
djl(ig,nb,nt) = ( tab (i0, nb, nt) * (-vx*wx-ux*wx-ux*vx)/6.d0 + & i3 = i0 + 3
tab (i1, nb, nt) * (+vx*wx-px*wx-px*vx)/2.d0 - & djl(ig,nb,nt) = ( tab (i0, nb, nt) * (-vx*wx-ux*wx-ux*vx)/6.d0 + &
tab (i2, nb, nt) * (+ux*wx-px*wx-px*ux)/2.d0 + & tab (i1, nb, nt) * (+vx*wx-px*wx-px*vx)/2.d0 - &
tab (i3, nb, nt) * (+ux*vx-px*vx-px*ux)/6.d0 )/dq tab (i2, nb, nt) * (+ux*wx-px*wx-px*ux)/2.d0 + &
#endif tab (i3, nb, nt) * (+ux*vx-px*vx-px*ux)/6.d0 )/dq
endif
enddo enddo
enddo enddo
enddo enddo

58
PW/gen_us_dy.f90
View File

@ -23,12 +23,8 @@ subroutine gen_us_dy (ik, u, dvkb)
USE gvect, ONLY : ig1, ig2, ig3, eigts1, eigts2, eigts3, g USE gvect, ONLY : ig1, ig2, ig3, eigts1, eigts2, eigts3, g
USE wvfct, ONLY : npw, npwx, igk USE wvfct, ONLY : npw, npwx, igk
USE uspp, ONLY : nkb, indv, nhtol, nhtolm USE uspp, ONLY : nkb, indv, nhtol, nhtolm
#ifdef USE_SPLINES USE us, ONLY : nqxq, tab, tab_d2y, dq, spline_ps
USE us, ONLY : nqxq, tab, tab_d2y, dq
USE splinelib USE splinelib
#else
USE us, ONLY : tab, dq
#endif
USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh
! !
implicit none implicit none
@ -49,10 +45,8 @@ subroutine gen_us_dy (ik, u, dvkb)
complex(DP), allocatable :: sk (:) complex(DP), allocatable :: sk (:)
complex(DP) :: phase, pref complex(DP) :: phase, pref
#ifdef USE_SPLINES integer :: iq
integer :: startq, lastq, iq
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
#endif
dvkb(:,:) = (0.d0, 0.d0) dvkb(:,:) = (0.d0, 0.d0)
if (lmaxkb.le.0) return if (lmaxkb.le.0) return
@ -79,36 +73,34 @@ subroutine gen_us_dy (ik, u, dvkb)
q (ig) = sqrt ( q(ig) ) * tpiba q (ig) = sqrt ( q(ig) ) * tpiba
end do end do
#ifdef USE_SPLINES if (spline_ps) then
startq = 1 allocate(xdata(nqxq))
lastq = nqxq do iq = 1, nqxq
!call divide (nqxq, startq, lastq) xdata(iq) = (iq - 1) * dq
allocate(xdata(lastq-startq+1)) enddo
do iq = startq, lastq endif
xdata(iq) = (iq - 1) * dq
enddo
#endif
do nt = 1, ntyp do nt = 1, ntyp
! calculate beta in G-space using an interpolation table ! calculate beta in G-space using an interpolation table
do nb = 1, nbeta (nt) do nb = 1, nbeta (nt)
do ig = 1, npw do ig = 1, npw
#ifdef USE_SPLINES if (spline_ps) then
vkb0(ig,nb,nt) = splint(xdata, tab(:,nb,nt), tab_d2y(:,nb,nt), q(ig)) vkb0(ig,nb,nt) = splint(xdata, tab(:,nb,nt), &
#else tab_d2y(:,nb,nt), q(ig))
px = q (ig) / dq - int (q (ig) / dq) else
ux = 1.d0 - px px = q (ig) / dq - int (q (ig) / dq)
vx = 2.d0 - px ux = 1.d0 - px
wx = 3.d0 - px vx = 2.d0 - px
i0 = q (ig) / dq + 1 wx = 3.d0 - px
i1 = i0 + 1 i0 = q (ig) / dq + 1
i2 = i0 + 2 i1 = i0 + 1
i3 = i0 + 3 i2 = i0 + 2
vkb0 (ig, nb, nt) = tab (i0, nb, nt) * ux * vx * wx / 6.d0 + & i3 = i0 + 3
tab (i1, nb, nt) * px * vx * wx / 2.d0 - & vkb0 (ig, nb, nt) = tab (i0, nb, nt) * ux * vx * wx / 6.d0 + &
tab (i2, nb, nt) * px * ux * wx / 2.d0 + & tab (i1, nb, nt) * px * vx * wx / 2.d0 - &
tab (i3, nb, nt) * px * ux * vx / 6.d0 tab (i2, nb, nt) * px * ux * wx / 2.d0 + &
#endif tab (i3, nb, nt) * px * ux * vx / 6.d0
endif
enddo enddo
enddo enddo
enddo enddo

44
PW/init_paw_1.f90
View File

@ -19,19 +19,12 @@ subroutine init_paw_1
USE cell_base , ONLY : omega USE cell_base , ONLY : omega
USE ions_base, ONLY : nat, ntyp => nsp, ityp USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE constants , ONLY : fpi USE constants , ONLY : fpi
#ifdef USE_SPLINES USE us, ONLY : nqxq, dq, nqx, tab, tab_d2y, qrad, spline_ps
USE us, ONLY : nqxq, dq, nqx, tab, tab_d2y, qrad
USE paw , ONLY : paw_nhm, paw_nh, paw_lmaxkb, paw_nkb, paw_nl, & USE paw , ONLY : paw_nhm, paw_nh, paw_lmaxkb, paw_nkb, paw_nl, &
paw_iltonh, paw_tab, aephi, paw_betar, psphi, & paw_iltonh, paw_tab, aephi, paw_betar, psphi, &
paw_indv, paw_nhtom, paw_nhtol, paw_nbeta, & paw_indv, paw_nhtom, paw_nhtol, paw_nbeta, &
paw_tab_d2y paw_tab_d2y
USE splinelib USE splinelib
#else
USE us, ONLY : nqx, dq
USE paw , ONLY : paw_nhm, paw_nh, paw_lmaxkb, paw_nkb, paw_nl, &
paw_iltonh, paw_tab, aephi, paw_betar, psphi, &
paw_indv, paw_nhtom, paw_nhtol, paw_nbeta
#endif
USE uspp, ONLY : ap, aainit USE uspp, ONLY : ap, aainit
USE atom , ONLY : r, rab, msh USE atom , ONLY : r, rab, msh
USE io_global, ONLY : stdout USE io_global, ONLY : stdout
@ -60,11 +53,9 @@ subroutine init_paw_1
integer :: n_overlap_warnings integer :: n_overlap_warnings
#ifdef USE_SPLINES
integer :: paw_nbeta_max integer :: paw_nbeta_max
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
real(DP) :: d1 real(DP) :: d1
#endif
call start_clock ('init_paw_1') call start_clock ('init_paw_1')
! !
@ -270,29 +261,24 @@ subroutine init_paw_1
call reduce (nqx * nbrx * ntyp, paw_tab) call reduce (nqx * nbrx * ntyp, paw_tab)
#endif #endif
#ifdef USE_SPLINES
!<ceres>
! initialize spline interpolation ! initialize spline interpolation
paw_nbeta_max = maxval ( paw_nbeta (:) ) if (spline_ps) then
allocate ( paw_tab_d2y ( nqx, paw_nbeta_max, ntyp ) ) paw_nbeta_max = maxval ( paw_nbeta (:) )
allocate ( paw_tab_d2y ( nqx, paw_nbeta_max, ntyp ) )
paw_tab_d2y = 0.0 paw_tab_d2y = 0.0
startq = 1 allocate(xdata(nqxq))
lastq = nqxq do iq = 1, nqxq
allocate(xdata(lastq-startq+1)) xdata(iq) = (iq - 1) * dq
do iq = startq, lastq enddo
xdata(iq) = (iq - 1) * dq do nt = 1, ntyp
enddo do nb = 1, paw_nbeta (nt)
do nt = 1, ntyp
do nb = 1, paw_nbeta (nt)
l = aephi(nt, nb)%label%l l = aephi(nt, nb)%label%l
d1 = (paw_tab(2,nb,nt) - paw_tab(1,nb,nt)) / dq d1 = (paw_tab(2,nb,nt) - paw_tab(1,nb,nt)) / dq
call spline(xdata, paw_tab(:,nb,nt), 0.d0, d1, paw_tab_d2y(:,nb,nt)) call spline(xdata, paw_tab(:,nb,nt), 0.d0, d1, paw_tab_d2y(:,nb,nt))
enddo enddo
enddo enddo
deallocate(xdata) deallocate(xdata)
!</ceres> endif
#endif
deallocate (ylmk0) deallocate (ylmk0)
deallocate (besr) deallocate (besr)

58
PW/init_paw_2.f90
View File

@ -21,15 +21,10 @@ subroutine init_paw_2 (npw_, igk_, q_, vkb_)
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau
USE gvect , ONLY : eigts1, eigts2, eigts3, g, ig1, ig2, ig3 USE gvect , ONLY : eigts1, eigts2, eigts3, g, ig1, ig2, ig3
USE us, ONLY : dq USE us, ONLY : dq
#ifdef USE_SPLINES
USE paw, ONLY : paw_nkb, paw_lmaxkb, paw_nhm, paw_nh, paw_nhtol, & USE paw, ONLY : paw_nkb, paw_lmaxkb, paw_nhm, paw_nh, paw_nhtol, &
paw_nhtom, paw_indv, paw_nbeta, paw_tab, paw_tab_d2y paw_nhtom, paw_indv, paw_nbeta, paw_tab, paw_tab_d2y
USE us, ONLY : nqxq, dq USE us, ONLY : nqxq, dq, spline_ps
USE splinelib USE splinelib
#else
USE paw, ONLY : paw_nkb, paw_lmaxkb, paw_nhm, paw_nh, paw_nhtol, &
paw_nhtom, paw_indv, paw_tab, paw_nbeta
#endif
! !
implicit none implicit none
! !
@ -51,10 +46,8 @@ subroutine init_paw_2 (npw_, igk_, q_, vkb_)
complex(DP) :: phase, pref complex(DP) :: phase, pref
complex(DP), allocatable :: sk(:) complex(DP), allocatable :: sk(:)
#ifdef USE_SPLINES
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
integer :: startq, lastq, iq integer :: iq
#endif
! !
! !
if (paw_lmaxkb.lt.0) return if (paw_lmaxkb.lt.0) return
@ -82,36 +75,35 @@ subroutine init_paw_2 (npw_, igk_, q_, vkb_)
qg(ig) = sqrt(qg(ig))*tpiba qg(ig) = sqrt(qg(ig))*tpiba
enddo enddo
#ifdef USE_SPLINES if (spline_ps) then
call divide (nqxq, startq, lastq) allocate(xdata(nqxq))
allocate(xdata(lastq-startq+1)) do iq = 1, nqxq
do iq = startq, lastq xdata(iq) = (iq - 1) * dq
xdata(iq) = (iq - 1) * dq enddo
enddo endif
#endif
jkb = 0 jkb = 0
do nt = 1, ntyp do nt = 1, ntyp
! calculate beta in G-space using an interpolation table ! calculate beta in G-space using an interpolation table
do nb = 1, paw_nbeta (nt) do nb = 1, paw_nbeta (nt)
do ig = 1, npw_ do ig = 1, npw_
#ifdef USE_SPLINES if (spline_ps) then
vq(ig) = splint(xdata, paw_tab(:,nb,nt), paw_tab_d2y(:,nb,nt), & vq(ig) = splint(xdata, paw_tab(:,nb,nt), paw_tab_d2y(:,nb,nt), &
qg(ig)) qg(ig))
#else else
px = qg (ig) / dq - int (qg (ig) / dq) px = qg (ig) / dq - int (qg (ig) / dq)
ux = 1.d0 - px ux = 1.d0 - px
vx = 2.d0 - px vx = 2.d0 - px
wx = 3.d0 - px wx = 3.d0 - px
i0 = qg (ig) / dq + 1 i0 = qg (ig) / dq + 1
i1 = i0 + 1 i1 = i0 + 1
i2 = i0 + 2 i2 = i0 + 2
i3 = i0 + 3 i3 = i0 + 3
vq (ig) = paw_tab (i0, nb, nt) * ux * vx * wx / 6.d0 + & vq (ig) = paw_tab (i0, nb, nt) * ux * vx * wx / 6.d0 + &
paw_tab (i1, nb, nt) * px * vx * wx / 2.d0 - & paw_tab (i1, nb, nt) * px * vx * wx / 2.d0 - &
paw_tab (i2, nb, nt) * px * ux * wx / 2.d0 + & paw_tab (i2, nb, nt) * px * ux * wx / 2.d0 + &
paw_tab (i3, nb, nt) * px * ux * vx / 6.d0 paw_tab (i3, nb, nt) * px * ux * vx / 6.d0
#endif endif
enddo enddo
! add spherical harmonic part ! add spherical harmonic part
do ih = 1, paw_nh (nt) do ih = 1, paw_nh (nt)

21
PW/init_us_1.f90
View File

@ -35,12 +35,8 @@ subroutine init_us_1
USE cell_base, ONLY : omega, tpiba USE cell_base, ONLY : omega, tpiba
USE gvect, ONLY : g, gg USE gvect, ONLY : g, gg
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
#ifdef USE_SPLINES USE us, ONLY : nqxq, dq, nqx, tab, tab_d2y, qrad, spline_ps
USE us, ONLY : nqxq, dq, nqx, tab, tab_d2y, qrad
USE splinelib USE splinelib
#else
USE us, ONLY : nqxq, dq, nqx, tab, qrad
#endif
USE uspp, ONLY : nhtol, nhtoj, nhtolm, dvan, qq, indv, ap, aainit, & USE uspp, ONLY : nhtol, nhtoj, nhtolm, dvan, qq, indv, ap, aainit, &
qq_so, dvan_so, okvan qq_so, dvan_so, okvan
USE uspp_param, ONLY : lmaxq, dion, betar, qfunc, qfcoef, rinner, nbeta, & USE uspp_param, ONLY : lmaxq, dion, betar, qfunc, qfcoef, rinner, nbeta, &
@ -73,10 +69,10 @@ subroutine init_us_1
integer, external :: sph_ind integer, external :: sph_ind
complex(DP) :: coeff, qgm(1) complex(DP) :: coeff, qgm(1)
real(DP) :: spinor, ji, jk real(DP) :: spinor, ji, jk
#ifdef USE_SPLINES
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
real(DP) :: d1 real(DP) :: d1
#endif
call start_clock ('init_us_1') call start_clock ('init_us_1')
! !
@ -373,12 +369,11 @@ subroutine init_us_1
call reduce (nqx * nbetam * ntyp, tab) call reduce (nqx * nbetam * ntyp, tab)
#endif #endif
#ifdef USE_SPLINES
! initialize spline interpolation ! initialize spline interpolation
startq = 1 if (spline_ps) then
lastq = nqxq allocate( xdata(nqxq) )
allocate(xdata(lastq-startq+1)) do iq = 1, nqxq
do iq = startq, lastq
xdata(iq) = (iq - 1) * dq xdata(iq) = (iq - 1) * dq
enddo enddo
do nt = 1, ntyp do nt = 1, ntyp
@ -389,7 +384,7 @@ subroutine init_us_1
enddo enddo
enddo enddo
deallocate(xdata) deallocate(xdata)
#endif endif
deallocate (ylmk0) deallocate (ylmk0)
deallocate (qtot) deallocate (qtot)

57
PW/init_us_2.f90
View File

@ -20,12 +20,8 @@ subroutine init_us_2 (npw_, igk_, q_, vkb_)
USE constants, ONLY : tpi USE constants, ONLY : tpi
USE gvect, ONLY : eigts1, eigts2, eigts3, ig1, ig2, ig3, g USE gvect, ONLY : eigts1, eigts2, eigts3, ig1, ig2, ig3, g
USE wvfct, ONLY : npw, npwx, igk USE wvfct, ONLY : npw, npwx, igk
#ifdef USE_SPLINES USE us, ONLY : nqxq, dq, tab, tab_d2y, spline_ps
USE us, ONLY : nqxq, dq, tab, tab_d2y
USE splinelib USE splinelib
#else
USE us, ONLY : dq, tab
#endif
USE uspp, ONLY : nkb, vkb, nhtol, nhtolm, indv USE uspp, ONLY : nkb, vkb, nhtol, nhtolm, indv
USE uspp_param, ONLY : lmaxkb, nbeta, nhm, nh USE uspp_param, ONLY : lmaxkb, nbeta, nhm, nh
! !
@ -49,10 +45,8 @@ subroutine init_us_2 (npw_, igk_, q_, vkb_)
complex(DP) :: phase, pref complex(DP) :: phase, pref
complex(DP), allocatable :: sk(:) complex(DP), allocatable :: sk(:)
#ifdef USE_SPLINES
real(DP), allocatable :: xdata(:) real(DP), allocatable :: xdata(:)
integer :: startq, lastq, iq integer :: iq
#endif
! !
! !
@ -80,37 +74,34 @@ subroutine init_us_2 (npw_, igk_, q_, vkb_)
qg(ig) = sqrt(qg(ig))*tpiba qg(ig) = sqrt(qg(ig))*tpiba
enddo enddo
#ifdef USE_SPLINES if (spline_ps) then
startq = 1 allocate(xdata(nqxq))
lastq = nqxq do iq = 1, nqxq
!call divide (nqxq, startq, lastq) xdata(iq) = (iq - 1) * dq
allocate(xdata(lastq-startq+1)) enddo
do iq = startq, lastq endif
xdata(iq) = (iq - 1) * dq
enddo
#endif
jkb = 0 jkb = 0
do nt = 1, ntyp do nt = 1, ntyp
! calculate beta in G-space using an interpolation table ! calculate beta in G-space using an interpolation table
do nb = 1, nbeta (nt) do nb = 1, nbeta (nt)
do ig = 1, npw_ do ig = 1, npw_
#ifdef USE_SPLINES if (spline_ps) then
vq(ig) = splint(xdata, tab(:,nb,nt), tab_d2y(:,nb,nt), qg(ig)) vq(ig) = splint(xdata, tab(:,nb,nt), tab_d2y(:,nb,nt), qg(ig))
#else else
px = qg (ig) / dq - int (qg (ig) / dq) px = qg (ig) / dq - int (qg (ig) / dq)
ux = 1.d0 - px ux = 1.d0 - px
vx = 2.d0 - px vx = 2.d0 - px
wx = 3.d0 - px wx = 3.d0 - px
i0 = INT( qg (ig) / dq ) + 1 i0 = INT( qg (ig) / dq ) + 1
i1 = i0 + 1 i1 = i0 + 1
i2 = i0 + 2 i2 = i0 + 2
i3 = i0 + 3 i3 = i0 + 3
vq (ig) = tab (i0, nb, nt) * ux * vx * wx / 6.d0 + & vq (ig) = tab (i0, nb, nt) * ux * vx * wx / 6.d0 + &
tab (i1, nb, nt) * px * vx * wx / 2.d0 - & tab (i1, nb, nt) * px * vx * wx / 2.d0 - &
tab (i2, nb, nt) * px * ux * wx / 2.d0 + & tab (i2, nb, nt) * px * ux * wx / 2.d0 + &
tab (i3, nb, nt) * px * ux * vx / 6.d0 tab (i3, nb, nt) * px * ux * vx / 6.d0
#endif endif
enddo enddo
! add spherical harmonic part ! add spherical harmonic part
do ih = 1, nh (nt) do ih = 1, nh (nt)

4
PW/input.f90
View File

@ -201,7 +201,7 @@ SUBROUTINE iosys()
edir, emaxpos, eopreg, eamp, noncolin, lambda, & edir, emaxpos, eopreg, eamp, noncolin, lambda, &
angle1, angle2, constrained_magnetization, & angle1, angle2, constrained_magnetization, &
B_field, fixed_magnetization, report, lspinorb,& B_field, fixed_magnetization, report, lspinorb,&
assume_molsys assume_molsys, spline_ps
! !
! ... ELECTRONS namelist ! ... ELECTRONS namelist
! !
@ -243,6 +243,7 @@ SUBROUTINE iosys()
USE constraints_module, ONLY : init_constraint USE constraints_module, ONLY : init_constraint
USE metadyn_vars, ONLY : init_metadyn_vars USE metadyn_vars, ONLY : init_metadyn_vars
USE read_namelists_module, ONLY : read_namelists, sm_not_set USE read_namelists_module, ONLY : read_namelists, sm_not_set
USE us, ONLY : spline_ps_ => spline_ps
! !
IMPLICIT NONE IMPLICIT NONE
! !
@ -1125,6 +1126,7 @@ SUBROUTINE iosys()
lambda_ = lambda lambda_ = lambda
! !
assume_molsys_ = assume_molsys assume_molsys_ = assume_molsys
spline_ps_ = spline_ps
! !
Hubbard_U_(1:ntyp) = hubbard_u(1:ntyp) Hubbard_U_(1:ntyp) = hubbard_u(1:ntyp)
Hubbard_alpha_(1:ntyp) = hubbard_alpha(1:ntyp) Hubbard_alpha_(1:ntyp) = hubbard_alpha(1:ntyp)

2
PW/paw.f90
View File

@ -39,12 +39,10 @@ MODULE paw
paw_becp(:,:) ! products of wavefunctions and proj paw_becp(:,:) ! products of wavefunctions and proj
REAL(DP), ALLOCATABLE :: & REAL(DP), ALLOCATABLE :: &
paw_tab(:,:,:) ! interpolation table for PPs paw_tab(:,:,:) ! interpolation table for PPs
#ifdef USE_SPLINES
!<ceres> !<ceres>
REAL(DP), ALLOCATABLE :: & REAL(DP), ALLOCATABLE :: &
paw_tab_d2y(:,:,:) ! for cubic splines paw_tab_d2y(:,:,:) ! for cubic splines
!</ceres> !</ceres>
#endif
! !
type wfc_label type wfc_label
integer :: na , & ! Atom number integer :: na , & ! Atom number

3
PW/pwcom.f90
View File

@ -489,10 +489,9 @@ MODULE us
qrad(:,:,:,:), &! radial FT of Q functions qrad(:,:,:,:), &! radial FT of Q functions
tab(:,:,:), &! interpolation table for PPs tab(:,:,:), &! interpolation table for PPs
tab_at(:,:,:) ! interpolation table for atomic wfc tab_at(:,:,:) ! interpolation table for atomic wfc
#ifdef USE_SPLINES LOGICAL :: spline_ps
REAL(DP), ALLOCATABLE :: & REAL(DP), ALLOCATABLE :: &
tab_d2y(:,:,:) ! for cubic splines tab_d2y(:,:,:) ! for cubic splines
#endif
! !
END MODULE us END MODULE us
! !