Commit 797172ac authored by Matthias Redies's avatar Matthias Redies

split off intgrf routines

parent 62d823df
......@@ -26,6 +26,7 @@ hybrid/wavefproducts.f90
hybrid/wavefproducts_aux.f90
hybrid/wavefproducts_inv.f90
hybrid/wavefproducts_noinv.f90
hybrid/util.f90
hybrid/rorder.f90
hybrid/util/util.f90
hybrid/util/rorder.f90
hybrid/util/intgrf.f90
)
......@@ -8,6 +8,7 @@ CONTAINS
USE m_types
USE m_hybrid_core
USE m_util
use m_intgrf
USE m_io_hybrid
IMPLICIT NONE
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
......
......@@ -45,7 +45,7 @@ CONTAINS
SUBROUTINE add_vnonlocal(nk, lapw, atoms, hybrid, dimension, kpts, jsp, results, xcpot, noco, hmat)
USE m_symm_hf, ONLY: symm_hf
USE m_util, ONLY: intgrf, intgrf_init
USE m_intgrf, ONLY: intgrf, intgrf_init
USE m_exchange_valence_hf
USE m_exchange_core
USE m_symmetrizeh
......
......@@ -8,7 +8,8 @@
& dimension, mpi, &
& input, sym, noco,&
& cell, lapw, jsp)
USE m_util, ONLY: intgrf, intgrf_init, chr, sphbessel, harmonicsr
USE m_util, ONLY: chr, sphbessel, harmonicsr
use m_intgrf, only: intgrf, intgrf_init
USE m_constants
USE m_types
USE m_io_hybrid
......
......@@ -43,7 +43,8 @@ CONTAINS
USE m_olap, ONLY: olap_pw
use m_types_hybrid, only: gptnorm
USE m_trafo, ONLY: symmetrize, bramat_trafo
USE m_util, ONLY: intgrf, intgrf_init, primitivef
USE m_intgrf, ONLY: intgrf, intgrf_init
use m_util, only: primitivef
USE m_hsefunctional, ONLY: change_coulombmatrix
USE m_wrapper
USE m_io_hybrid
......@@ -1532,6 +1533,7 @@ CONTAINS
USE m_wrapper
USE m_trafo
USE m_util
use m_intgrf
USE m_olap
IMPLICIT NONE
......
......@@ -24,7 +24,8 @@ CONTAINS
ex_vv_r, ex_vv_c, l_real)
USE m_constants
USE m_util
use m_util
USE m_intgrf
USE m_wrapper
USE m_types
USE m_io_hybrid
......@@ -239,6 +240,7 @@ CONTAINS
USE m_constants
USE m_util
use m_intgrf
USE m_wrapper
USE m_types
USE m_io_hybrid
......@@ -413,6 +415,7 @@ CONTAINS
USE m_constants
USE m_util
use m_intgrf
USE m_wrapper
USE m_gaunt
USE m_trafo
......@@ -555,6 +558,7 @@ CONTAINS
USE m_constants
USE m_util
use m_intgrf
USE m_wrapper
USE m_gaunt
USE m_trafo
......
......@@ -13,6 +13,7 @@ CONTAINS
USE m_types
USE m_eig66_io
USE m_util
USE m_intgrf
USE m_checkolap
USE m_hybrid_core
USE m_gen_wavf
......
......@@ -972,8 +972,8 @@ CONTAINS
USE m_constants
USE m_types_hybrid, ONLY: gptnorm
USE m_util, ONLY: sphbessel, pure_intgrf, intgrf_init, intgrf_out, NEGATIVE_EXPONENT_WARNING, NEGATIVE_EXPONENT_ERROR
USE m_util, ONLY: sphbessel
use m_intgrf, only: pure_intgrf, intgrf_init, intgrf_out,NEGATIVE_EXPONENT_WARNING, NEGATIVE_EXPONENT_ERROR
IMPLICIT NONE
! scalar input
......@@ -1304,7 +1304,8 @@ CONTAINS
potential, fourier_trafo)
USE m_constants
USE m_util, ONLY: sphbessel, pure_intgrf, intgrf_init, intgrf_out, NEGATIVE_EXPONENT_WARNING, NEGATIVE_EXPONENT_ERROR
USE m_util, ONLY: sphbessel
use m_intgrf, only: pure_intgrf, intgrf_init, intgrf_out, NEGATIVE_EXPONENT_WARNING, NEGATIVE_EXPONENT_ERROR
USE m_trafo, ONLY: symmetrize
IMPLICIT NONE
......@@ -2132,6 +2133,7 @@ CONTAINS
USE m_constants
USE m_util
USE m_intgrf
USE m_wrapper
IMPLICIT NONE
......@@ -2382,6 +2384,7 @@ CONTAINS
USE m_constants
USE m_util
use m_intgrf
USE m_wrapper
USE m_gaunt
USE m_trafo
......
......@@ -49,7 +49,7 @@ CONTAINS
USE m_types
USE m_symm_hf
USE m_util, ONLY: intgrf, intgrf_init
USE m_intgrf, ONLY: intgrf, intgrf_init
USE m_exchange_valence_hf
USE m_exchange_core
USE m_symmetrizeh
......
......@@ -17,6 +17,7 @@ CONTAINS
USE m_ylm
USE m_gaunt
USE m_util
use m_intgrf
USE m_types
USE m_io_hybrid
IMPLICIT NONE
......@@ -786,7 +787,8 @@ CONTAINS
USE m_olap
USE m_wrapper
USE m_util, only: derivative, intgrf_init, intgrf
USE m_util, only: derivative
use m_intgrf, only: intgrf_init, intgrf
USE m_dr2fdr
USE m_constants
USE m_types
......
......@@ -42,7 +42,7 @@ CONTAINS
USE m_judft
USE m_loddop, ONLY: loddop
USE m_util, ONLY: intgrf_init, intgrf
USE m_intgrf, ONLY: intgrf_init, intgrf
use m_rorder, only: rorderpf
USE m_hybrid_core
USE m_wrapper
......@@ -588,7 +588,7 @@ CONTAINS
use m_types
USE m_radfun, ONLY: radfun
USE m_radflo, ONLY: radflo
USE m_util, ONLY: intgrf
USE m_intgrf, ONLY: intgrf
implicit NONE
type(t_atoms), intent(in) :: atoms
type(t_enpara), intent(in) :: enpara
......@@ -667,7 +667,7 @@ CONTAINS
subroutine gen_gvec(cell, kpts, mpbasis)
use m_types
USE m_util, ONLY: intgrf_init, intgrf
USE m_intgrf, ONLY: intgrf_init, intgrf
use m_rorder, only: rorderpf
implicit NONE
type(t_cell), intent(in) :: cell
......
......@@ -157,7 +157,7 @@ CONTAINS
& atoms, hybrid, cell,&
& bas1, bas2)
USE m_util, ONLY: intgrf, intgrf_init
USE m_intgrf, ONLY: intgrf, intgrf_init
USE m_types
IMPLICIT NONE
TYPE(t_hybrid), INTENT(IN) :: hybrid
......
......@@ -19,7 +19,6 @@ CONTAINS
USE m_types
USE m_util, ONLY: modulo1
IMPLICIT NONE
TYPE(t_sym), INTENT(IN) :: sym
......@@ -73,7 +72,8 @@ CONTAINS
USE m_constants
USE m_types
USE m_util, ONLY: modulo1, intgrf, intgrf_init
USE m_util, ONLY: modulo1
use m_intgrf, only: intgrf, intgrf_init
USE m_olap, ONLY: wfolap_inv, wfolap_noinv, wfolap1, wfolap_init
USE m_trafo, ONLY: waveftrafo_symm
USE m_io_hybrid
......
module m_intgrf
implicit none
TYPE :: intgrf_out
REAL :: value ! value of the integration
INTEGER :: ierror ! error code
END TYPE intgrf_out
! error and warning codes for intgrf function
INTEGER, PARAMETER :: NO_ERROR = 0
INTEGER, PARAMETER :: NEGATIVE_EXPONENT_WARNING = 1
INTEGER, PARAMETER :: NEGATIVE_EXPONENT_ERROR = 2
CONTAINS
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Integrates function f numerically (Lagrange and Simpson integration)
! on grid(itype) and is much faster than intgr.
! (Only normal outward integration.)
! Before first use of this function it has to be initialized with
! intgrf_init.
FUNCTION intgrf(f, jri, jmtd, rmsh, dx, ntype, itype, gridf)
use m_juDFT
IMPLICIT NONE
REAL :: intgrf
INTEGER, INTENT(IN) :: itype, ntype, jmtd
INTEGER, INTENT(IN) :: jri(ntype)
REAL, INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
REAL, INTENT(IN) :: gridf(jmtd, ntype)
REAL, INTENT(IN) :: f(*)
! - local -
TYPE(intgrf_out) :: fct_res
fct_res = pure_intgrf(f, jri, jmtd, rmsh, dx, ntype, itype, gridf)
IF (fct_res%ierror == NEGATIVE_EXPONENT_WARNING) THEN
write (6, *) 'intgrf: Warning!'// &
'Negative exponent x in extrapolation a+c*r**x'
ELSEIF (fct_res%ierror == NEGATIVE_EXPONENT_ERROR) THEN
write (6, *) &
'intgrf: Negative exponent x in extrapolation a+c*r**x'
CALL juDFT_error( &
'intgrf: Negative exponent x in extrapolation a+c*r**x')
END IF
intgrf = fct_res%value
END FUNCTION intgrf
! pure wrapper for intgrf with same functionality
! can be used within forall loops
PURE FUNCTION pure_intgrf(f, jri, jmtd, rmsh, dx, ntype, itype, gridf)
IMPLICIT NONE
TYPE(intgrf_out) :: pure_intgrf
INTEGER, INTENT(IN) :: itype, ntype, jmtd
INTEGER, INTENT(IN) :: jri(ntype)
REAL, INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
REAL, INTENT(IN) :: gridf(jmtd, ntype)
REAL, INTENT(IN) :: f(*)
! - local -
INTEGER :: n
REAL :: r1, h, a, x
n = jri(itype)
r1 = rmsh(1, itype)
h = dx(itype)
pure_intgrf%ierror = NO_ERROR
! integral from 0 to r1 approximated by leading term in power series expansion
IF (f(1)*f(2) > 1e-10 .AND. h > 0) THEN
IF (f(2) == f(1)) THEN
pure_intgrf%value = r1*f(1)
ELSE
x = (f(3) - f(2))/(f(2) - f(1))
a = (f(2) - x*f(1))/(1 - x)
x = log(x)/h
IF (x < 0) THEN
IF (x > -1) THEN
pure_intgrf%ierror = NEGATIVE_EXPONENT_WARNING
ELSE IF (x <= -1) THEN
pure_intgrf%ierror = NEGATIVE_EXPONENT_ERROR
RETURN
END IF
END IF
pure_intgrf%value = r1*(f(1) + x*a)/(x + 1)
! x = f(2) / f(1)
! x = log(x)/h
! IF(x.lt.0) THEN
! IF(x.gt.-1) write(6,'(A,ES9.1)') 'intgrf: Warning!&
! & Negative exponent x in&
! & extrapolation c*r**x:',x
! IF(x.le.-1) write(6,'(A,ES9.1)') 'intgrf: Negative exponent&
! & x in extrapolation&
! & c*r**x:',x
! IF(x.le.-1) call juDFT_error('intgrf: Negative exponent&
! & x in extrapolation &
! & c*r**x')
! END IF
! intgrf = (r1*f(1))/(x+1)
END IF
ELSE
pure_intgrf%value = 0
END IF
! integrate from r(1) to r(n) by multiplying with gridf
pure_intgrf%value = pure_intgrf%value &
+ dot_product(gridf(:n, itype), f(:n))
END FUNCTION pure_intgrf
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Initializes fast numerical integration intgrf (see below).
SUBROUTINE intgrf_init(ntype, jmtd, jri, dx, rmsh, gridf)
IMPLICIT NONE
INTEGER, INTENT(IN) :: ntype, jmtd
INTEGER, INTENT(IN) :: jri(ntype)
REAL, INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
REAL, INTENT(OUT), ALLOCATABLE :: gridf(:, :)
! - local -
INTEGER :: i, j, itype
INTEGER :: n, nstep, n0 = 6
INTEGER, PARAMETER :: simpson(7) = (/41, 216, 27, 272, 27, 216, 41/)
REAL :: r1, h, dr
REAL :: r(7)
REAL, PARAMETER :: lagrange(7, 6) = reshape( &
(/19087., 65112., -46461., 37504., -20211., 6312., -863., &
-863., 25128., 46989., -16256., 7299., -2088., 271., &
& 271., -2760., 30819., 37504., -6771., 1608., -191., &
-191., 1608., -6771., 37504., 30819., -2760., 271., &
& 271., -2088., 7299., -16256., 46989., 25128., -863., &
-863., 6312., -20211., 37504., -46461., 65112., 19087./), & ! The last row is actually never used.
(/7, 6/))
n = jmtd
ALLOCATE (gridf(n, ntype))
gridf = 0
DO itype = 1, ntype
n = jri(itype)
r1 = rmsh(1, itype)
h = dx(itype)
nstep = (n - 1)/6
n0 = n - 6*nstep
dr = exp(h)
! Calculate Lagrange-integration coefficients from r(1) to r(n0)
r(1) = r1
IF (n0 > 1) THEN
DO i = 2, 7
r(i) = r(i - 1)*dr
END DO
DO i = 1, 7
gridf(i, itype) = h/60480*r(i)*sum(lagrange(i, 1:n0 - 1))
END DO
r(1) = r(n0)
END IF
! Calculate Simpson-integration coefficients from r(n0) to r(n)
DO i = 1, nstep
DO j = 2, 7
r(j) = r(j - 1)*dr
END DO
DO j = n0, n0 + 6
gridf(j, itype) = gridf(j, itype) + h/140*r(j - n0 + 1)* &
simpson(j - n0 + 1)
END DO
n0 = n0 + 6
r(1) = r(7)
END DO
END DO
END SUBROUTINE intgrf_init
end module m_intgrf
......@@ -2,13 +2,8 @@ MODULE m_util
USE m_juDFT
! error and warning codes for intgrf function
INTEGER, PARAMETER :: NO_ERROR = 0
INTEGER, PARAMETER :: NEGATIVE_EXPONENT_WARNING = 1
INTEGER, PARAMETER :: NEGATIVE_EXPONENT_ERROR = 2
! return type of the pure intgrf function
TYPE :: intgrf_out
REAL :: value ! value of the integration
INTEGER :: ierror ! error code
END TYPE intgrf_out
INTERFACE derivative
MODULE PROCEDURE derivative_t, derivative_nt
......@@ -98,181 +93,7 @@ CONTAINS
END IF
END FUNCTION wigner3j
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Integrates function f numerically (Lagrange and Simpson integration)
! on grid(itype) and is much faster than intgr.
! (Only normal outward integration.)
! Before first use of this function it has to be initialized with
! intgrf_init.
FUNCTION intgrf(f, jri, jmtd, rmsh, dx, ntype, itype, gridf)
IMPLICIT NONE
REAL :: intgrf
INTEGER, INTENT(IN) :: itype, ntype, jmtd
INTEGER, INTENT(IN) :: jri(ntype)
REAL , INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
REAL , INTENT(IN) :: gridf(jmtd, ntype)
REAL , INTENT(IN) :: f(*)
! - local -
TYPE(intgrf_out) :: fct_res
fct_res = pure_intgrf(f, jri, jmtd, rmsh, dx, ntype, itype, gridf)
IF (fct_res%ierror == NEGATIVE_EXPONENT_WARNING) THEN
write (6, *) 'intgrf: Warning!'// &
'Negative exponent x in extrapolation a+c*r**x'
ELSEIF (fct_res%ierror == NEGATIVE_EXPONENT_ERROR) THEN
write (6, *) &
'intgrf: Negative exponent x in extrapolation a+c*r**x'
CALL juDFT_error( &
'intgrf: Negative exponent x in extrapolation a+c*r**x')
END IF
intgrf = fct_res%value
END FUNCTION intgrf
! pure wrapper for intgrf with same functionality
! can be used within forall loops
PURE FUNCTION pure_intgrf(f, jri, jmtd, rmsh, dx, ntype, itype, gridf)
IMPLICIT NONE
TYPE(intgrf_out) :: pure_intgrf
INTEGER, INTENT(IN) :: itype, ntype, jmtd
INTEGER, INTENT(IN) :: jri(ntype)
REAL , INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
REAL , INTENT(IN) :: gridf(jmtd, ntype)
REAL , INTENT(IN) :: f(*)
! - local -
INTEGER :: n
REAL :: r1, h, a, x
n = jri(itype)
r1 = rmsh(1, itype)
h = dx(itype)
pure_intgrf%ierror = NO_ERROR
! integral from 0 to r1 approximated by leading term in power series expansion
IF (f(1)*f(2) > 1e-10 .AND. h > 0) THEN
IF (f(2) == f(1)) THEN
pure_intgrf%value = r1*f(1)
ELSE
x = (f(3) - f(2))/(f(2) - f(1))
a = (f(2) - x*f(1))/(1 - x)
x = log(x)/h
IF (x < 0) THEN
IF (x > -1) THEN
pure_intgrf%ierror = NEGATIVE_EXPONENT_WARNING
ELSE IF (x <= -1) THEN
pure_intgrf%ierror = NEGATIVE_EXPONENT_ERROR
RETURN
END IF
END IF
pure_intgrf%value = r1*(f(1) + x*a)/(x + 1)
! x = f(2) / f(1)
! x = log(x)/h
! IF(x.lt.0) THEN
! IF(x.gt.-1) write(6,'(A,ES9.1)') 'intgrf: Warning!&
! & Negative exponent x in&
! & extrapolation c*r**x:',x
! IF(x.le.-1) write(6,'(A,ES9.1)') 'intgrf: Negative exponent&
! & x in extrapolation&
! & c*r**x:',x
! IF(x.le.-1) call juDFT_error('intgrf: Negative exponent&
! & x in extrapolation &
! & c*r**x')
! END IF
! intgrf = (r1*f(1))/(x+1)
END IF
ELSE
pure_intgrf%value = 0
END IF
! integrate from r(1) to r(n) by multiplying with gridf
pure_intgrf%value = pure_intgrf%value &
+ dot_product(gridf(:n, itype), f(:n))
END FUNCTION pure_intgrf
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Initializes fast numerical integration intgrf (see below).
SUBROUTINE intgrf_init(ntype, jmtd, jri, dx, rmsh, gridf)
IMPLICIT NONE
INTEGER, INTENT(IN) :: ntype, jmtd
INTEGER, INTENT(IN) :: jri(ntype)
REAL, INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
REAL, INTENT(OUT), ALLOCATABLE :: gridf(:, :)
! - local -
INTEGER :: i, j, itype
INTEGER :: n, nstep, n0 = 6
INTEGER, PARAMETER :: simpson(7) = (/41, 216, 27, 272, 27, 216, 41/)
REAL :: r1, h, dr
REAL :: r(7)
REAL, PARAMETER :: lagrange(7, 6) = reshape( &
(/19087., 65112., -46461., 37504., -20211., 6312., -863., &
-863., 25128., 46989., -16256., 7299., -2088., 271., &
& 271., -2760., 30819., 37504., -6771., 1608., -191., &
-191., 1608., -6771., 37504., 30819., -2760., 271., &
& 271., -2088., 7299., -16256., 46989., 25128., -863., &
-863., 6312., -20211., 37504., -46461., 65112., 19087./), & ! The last row is actually never used.
(/7, 6/))
n = jmtd
ALLOCATE (gridf(n, ntype))
gridf = 0
DO itype = 1, ntype
n = jri(itype)
r1 = rmsh(1, itype)
h = dx(itype)
nstep = (n - 1)/6
n0 = n - 6*nstep
dr = exp(h)
! Calculate Lagrange-integration coefficients from r(1) to r(n0)
r(1) = r1
IF (n0 > 1) THEN
DO i = 2, 7
r(i) = r(i - 1)*dr
END DO
DO i = 1, 7
gridf(i, itype) = h/60480*r(i)*sum(lagrange(i, 1:n0 - 1))
END DO
r(1) = r(n0)
END IF
! Calculate Simpson-integration coefficients from r(n0) to r(n)
DO i = 1, nstep
DO j = 2, 7
r(j) = r(j - 1)*dr
END DO
DO j = n0, n0 + 6
gridf(j, itype) = gridf(j, itype) + h/140*r(j - n0 + 1)* &
simpson(j - n0 + 1)
END DO
n0 = n0 + 6
r(1) = r(7)
END DO
END DO
END SUBROUTINE intgrf_init
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment