Commit 1730364b authored by Matthias Redies's avatar Matthias Redies
Browse files

format all the hybrid code

parent 567d369d
MODULE m_hf_init MODULE m_hf_init
! !
! preparations for HF and hybrid functional calculation ! preparations for HF and hybrid functional calculation
! !
CONTAINS CONTAINS
SUBROUTINE hf_init(hybrid,kpts,atoms,input,DIMENSION,hybdat,l_real) SUBROUTINE hf_init(hybrid, kpts, atoms, input, DIMENSION, hybdat, l_real)
USE m_types USE m_types
USE m_read_core USE m_read_core
USE m_util USE m_util
USE m_io_hybrid USE m_io_hybrid
IMPLICIT NONE IMPLICIT NONE
TYPE(t_hybrid),INTENT(INOUT) :: hybrid TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_kpts),INTENT(IN) :: kpts TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_input),INTENT(IN) :: input TYPE(t_input), INTENT(IN) :: input
TYPE(t_dimension),INTENT(IN) :: DIMENSION TYPE(t_dimension), INTENT(IN) :: DIMENSION
TYPE(t_hybdat),INTENT(OUT) :: hybdat TYPE(t_hybdat), INTENT(OUT) :: hybdat
LOGICAL,INTENT(IN) :: l_real LOGICAL, INTENT(IN) :: l_real
INTEGER:: itype,ieq,l,m,i,nk,l1,l2,m1,m2,ok
INTEGER:: itype, ieq, l, m, i, nk, l1, l2, m1, m2, ok
!initialize hybdat%gridf for radial integration
CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
!Alloc variables
ALLOCATE(hybdat%lmaxc(atoms%ntype))
ALLOCATE(hybdat%bas1(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
ALLOCATE(hybdat%bas2(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
ALLOCATE(hybdat%bas1_MT(hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
ALLOCATE(hybdat%drbas1_MT(hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
!sym%tau = oneD%ods%tau
! preparations for core states !initialize hybdat%gridf for radial integration
CALL core_init(dimension,input,atoms, hybdat%lmaxcd,hybdat%maxindxc) CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, hybdat%gridf)
ALLOCATE( hybdat%nindxc(0:hybdat%lmaxcd,atoms%ntype), stat = ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation hybdat%nindxc'
ALLOCATE( hybdat%core1(atoms%jmtd,hybdat%maxindxc,0:hybdat%lmaxcd,atoms%ntype), stat=ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation core1'
ALLOCATE( hybdat%core2(atoms%jmtd,hybdat%maxindxc,0:hybdat%lmaxcd,atoms%ntype), stat=ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation core2'
ALLOCATE( hybdat%eig_c(hybdat%maxindxc,0:hybdat%lmaxcd,atoms%ntype), stat=ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation hybdat%eig_c'
hybdat%nindxc = 0 ; hybdat%core1 = 0 ; hybdat%core2 = 0 ; hybdat%eig_c = 0
! pre-calculate gaunt coefficients !Alloc variables
ALLOCATE (hybdat%lmaxc(atoms%ntype))
ALLOCATE (hybdat%bas1(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
ALLOCATE (hybdat%bas2(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
ALLOCATE (hybdat%bas1_MT(hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
ALLOCATE (hybdat%drbas1_MT(hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
hybdat%maxfac = max(2*atoms%lmaxd+hybrid%maxlcutm1+1,2*hybdat%lmaxcd+2*atoms%lmaxd+1) !sym%tau = oneD%ods%tau
ALLOCATE ( hybdat%fac( 0:hybdat%maxfac),hybdat%sfac( 0:hybdat%maxfac),stat=ok)
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation fac,hybdat%sfac'
hybdat%fac(0) = 1
hybdat%sfac(0) = 1
DO i=1,hybdat%maxfac
hybdat%fac(i) = hybdat%fac(i-1)*i ! hybdat%fac(i) = i!
hybdat%sfac(i) = hybdat%sfac(i-1)*sqrt(i*1.0) ! hybdat%sfac(i) = sqrt(i!)
END DO
! preparations for core states
CALL core_init(dimension, input, atoms, hybdat%lmaxcd, hybdat%maxindxc)
ALLOCATE (hybdat%nindxc(0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation hybdat%nindxc'
ALLOCATE (hybdat%core1(atoms%jmtd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation core1'
ALLOCATE (hybdat%core2(atoms%jmtd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation core2'
ALLOCATE (hybdat%eig_c(hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation hybdat%eig_c'
hybdat%nindxc = 0; hybdat%core1 = 0; hybdat%core2 = 0; hybdat%eig_c = 0
ALLOCATE ( hybdat%gauntarr( 2, 0:atoms%lmaxd, 0:atoms%lmaxd, 0:hybrid%maxlcutm1, -atoms%lmaxd:atoms%lmaxd ,-hybrid%maxlcutm1:hybrid%maxlcutm1 ),stat=ok) ! pre-calculate gaunt coefficients
IF( ok .ne. 0 ) STOP 'eigen: failure allocation hybdat%gauntarr'
hybdat%gauntarr = 0
DO l2 = 0,atoms%lmaxd
DO l1 = 0,atoms%lmaxd
DO l = abs(l1-l2),min(l1+l2,hybrid%maxlcutm1)
DO m = -l,l
DO m1 = -l1,l1
m2 = m1 + m ! Gaunt condition -m1+m2-m = 0
IF(abs(m2).le.l2) hybdat%gauntarr(1,l1,l2,l,m1,m) = gaunt(l1,l2,l,m1,m2,m,hybdat%maxfac,hybdat%fac,hybdat%sfac)
m2 = m1 - m ! switch role of l2-index
IF(abs(m2).le.l2) hybdat%gauntarr(2,l1,l2,l,m1,m) = gaunt(l2,l1,l,m2,m1,m,hybdat%maxfac,hybdat%fac,hybdat%sfac)
END DO
END DO
END DO
END DO
END DO
!skip_kpt = .false. hybdat%maxfac = max(2*atoms%lmaxd + hybrid%maxlcutm1 + 1, 2*hybdat%lmaxcd + 2*atoms%lmaxd + 1)
ALLOCATE (hybdat%fac(0:hybdat%maxfac), hybdat%sfac(0:hybdat%maxfac), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation fac,hybdat%sfac'
hybdat%fac(0) = 1
hybdat%sfac(0) = 1
DO i = 1, hybdat%maxfac
hybdat%fac(i) = hybdat%fac(i - 1)*i ! hybdat%fac(i) = i!
hybdat%sfac(i) = hybdat%sfac(i - 1)*sqrt(i*1.0) ! hybdat%sfac(i) = sqrt(i!)
END DO
END SUBROUTINE hf_init ALLOCATE (hybdat%gauntarr(2, 0:atoms%lmaxd, 0:atoms%lmaxd, 0:hybrid%maxlcutm1, -atoms%lmaxd:atoms%lmaxd, -hybrid%maxlcutm1:hybrid%maxlcutm1), stat=ok)
IF (ok /= 0) STOP 'eigen: failure allocation hybdat%gauntarr'
hybdat%gauntarr = 0
DO l2 = 0, atoms%lmaxd
DO l1 = 0, atoms%lmaxd
DO l = abs(l1 - l2), min(l1 + l2, hybrid%maxlcutm1)
DO m = -l, l
DO m1 = -l1, l1
m2 = m1 + m ! Gaunt condition -m1+m2-m = 0
IF (abs(m2) <= l2) hybdat%gauntarr(1, l1, l2, l, m1, m) = gaunt(l1, l2, l, m1, m2, m, hybdat%maxfac, hybdat%fac, hybdat%sfac)
m2 = m1 - m ! switch role of l2-index
IF (abs(m2) <= l2) hybdat%gauntarr(2, l1, l2, l, m1, m) = gaunt(l2, l1, l, m2, m1, m, hybdat%maxfac, hybdat%fac, hybdat%sfac)
END DO
END DO
END DO
END DO
END DO
!skip_kpt = .false.
END SUBROUTINE hf_init
END MODULE m_hf_init END MODULE m_hf_init
...@@ -40,95 +40,95 @@ MODULE m_add_vnonlocal ...@@ -40,95 +40,95 @@ MODULE m_add_vnonlocal
! c ! c
! M.Betzinger (09/07) c ! M.Betzinger (09/07) c
! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c ! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c
CONTAINS CONTAINS
SUBROUTINE add_vnonlocal(nk,lapw,atoms,hybrid,dimension,kpts,jsp,results,xcpot,noco,hmat) SUBROUTINE add_vnonlocal(nk, lapw, atoms, hybrid, dimension, kpts, jsp, results, xcpot, noco, hmat)
USE m_symm_hf, ONLY: symm_hf USE m_symm_hf, ONLY: symm_hf
USE m_util, ONLY: intgrf,intgrf_init USE m_util, ONLY: intgrf, intgrf_init
USE m_exchange_valence_hf USE m_exchange_valence_hf
USE m_exchange_core USE m_exchange_core
USE m_symmetrizeh USE m_symmetrizeh
USE m_wrapper USE m_wrapper
USE m_hsefunctional, ONLY: exchange_vccvHSE,exchange_ccccHSE USE m_hsefunctional, ONLY: exchange_vccvHSE, exchange_ccccHSE
USE m_types USE m_types
USE m_io_hybrid USE m_io_hybrid
IMPLICIT NONE IMPLICIT NONE
TYPE(t_results), INTENT(INOUT) :: results TYPE(t_results), INTENT(INOUT) :: results
CLASS(t_xcpot), INTENT(IN) :: xcpot CLASS(t_xcpot), INTENT(IN) :: xcpot
TYPE(t_dimension), INTENT(IN) :: dimension TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_hybrid), INTENT(INOUT) :: hybrid TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_kpts), INTENT(IN) :: kpts TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_lapw), INTENT(IN) :: lapw TYPE(t_lapw), INTENT(IN) :: lapw
TYPE(t_atoms), INTENT(IN) :: atoms TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_noco), INTENT(IN) :: noco TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_mat), INTENT(INOUT) :: hmat TYPE(t_mat), INTENT(INOUT) :: hmat
INTEGER, INTENT(IN) :: jsp INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: nk INTEGER, INTENT(IN) :: nk
! local scalars ! local scalars
INTEGER :: n,nn,iband,nbasfcn INTEGER :: n, nn, iband, nbasfcn
REAL :: a_ex REAL :: a_ex
TYPE(t_mat) :: olap,tmp,v_x,z TYPE(t_mat) :: olap, tmp, v_x, z
COMPLEX :: exch(dimension%neigd,dimension%neigd) COMPLEX :: exch(dimension%neigd, dimension%neigd)
! initialize weighting factor for HF exchange part ! initialize weighting factor for HF exchange part
a_ex=xcpot%get_exchange_weight() a_ex = xcpot%get_exchange_weight()
nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco) nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
CALL v_x%init(hmat%l_real,nbasfcn,nbasfcn) CALL v_x%init(hmat%l_real, nbasfcn, nbasfcn)
CALL read_v_x(v_x,kpts%nkpt*(jsp-1)+nk) CALL read_v_x(v_x, kpts%nkpt*(jsp - 1) + nk)
! add non-local x-potential to the hamiltonian hmat ! add non-local x-potential to the hamiltonian hmat
DO n = 1, v_x%matsize1 DO n = 1, v_x%matsize1
DO nn = 1, n DO nn = 1, n
IF (hmat%l_real) THEN IF (hmat%l_real) THEN
hmat%data_r(nn,n) = hmat%data_r(nn,n) - a_ex*v_x%data_r(nn,n) hmat%data_r(nn, n) = hmat%data_r(nn, n) - a_ex*v_x%data_r(nn, n)
ELSE ELSE
hmat%data_c(nn,n) = hmat%data_c(nn,n) - a_ex*v_x%data_c(nn,n) hmat%data_c(nn, n) = hmat%data_c(nn, n) - a_ex*v_x%data_c(nn, n)
ENDIF ENDIF
END DO END DO
END DO END DO
! calculate HF energy ! calculate HF energy
IF(hybrid%l_calhf) THEN IF (hybrid%l_calhf) THEN
WRITE(6,'(A)') new_line('n')//new_line('n')//' ### '// ' diagonal HF exchange elements (eV) ###' WRITE (6, '(A)') new_line('n')//new_line('n')//' ### '//' diagonal HF exchange elements (eV) ###'
WRITE(6,'(A)') new_line('n') // ' k-point '// 'band tail pole total(valence+core)' WRITE (6, '(A)') new_line('n')//' k-point '//'band tail pole total(valence+core)'
END IF END IF
! read in lower triangle part of overlap matrix from direct acces file olap ! read in lower triangle part of overlap matrix from direct acces file olap
CALL olap%init(hmat%l_real,nbasfcn,nbasfcn) CALL olap%init(hmat%l_real, nbasfcn, nbasfcn)
CALL read_olap(olap,kpts%nkpt*(jsp-1)+nk) CALL read_olap(olap, kpts%nkpt*(jsp - 1) + nk)
IF (.NOT.olap%l_real) olap%data_c=conjg(olap%data_c) IF (.NOT. olap%l_real) olap%data_c = conjg(olap%data_c)
CALL z%init(olap%l_real, nbasfcn, dimension%neigd)
CALL z%init(olap%l_real,nbasfcn,dimension%neigd) CALL read_z(z, kpts%nkpt*(jsp - 1) + nk)
CALL read_z(z,kpts%nkpt*(jsp-1)+nk)
! calculate exchange contribution of current k-point nk to total energy (te_hfex) ! calculate exchange contribution of current k-point nk to total energy (te_hfex)
! in the case of a spin-unpolarized calculation the factor 2 is added in eigen.F90 ! in the case of a spin-unpolarized calculation the factor 2 is added in eigen.F90
IF (.NOT.v_x%l_real) v_x%data_c=conjg(v_x%data_c) IF (.NOT. v_x%l_real) v_x%data_c = conjg(v_x%data_c)
exch = 0 exch = 0
z%matsize1=MIN(z%matsize1,v_x%matsize2) z%matsize1 = MIN(z%matsize1, v_x%matsize2)
CALL v_x%multiply(z,tmp) CALL v_x%multiply(z, tmp)
DO iband = 1, hybrid%nbands(nk) DO iband = 1, hybrid%nbands(nk)
IF (z%l_real) THEN IF (z%l_real) THEN
exch(iband,iband) = dot_product(z%data_r(:z%matsize1,iband),tmp%data_r(:,iband)) exch(iband, iband) = dot_product(z%data_r(:z%matsize1, iband), tmp%data_r(:, iband))
ELSE ELSE
exch(iband,iband) = dot_product(z%data_c(:z%matsize1,iband),tmp%data_c(:,iband)) exch(iband, iband) = dot_product(z%data_c(:z%matsize1, iband), tmp%data_c(:, iband))
END IF END IF
IF(iband.LE.hybrid%nobd(nk)) THEN IF (iband <= hybrid%nobd(nk)) THEN
results%te_hfex%valence = results%te_hfex%valence -a_ex*results%w_iks(iband,nk,jsp)*exch(iband,iband) results%te_hfex%valence = results%te_hfex%valence - a_ex*results%w_iks(iband, nk, jsp)*exch(iband, iband)
END IF END IF
IF(hybrid%l_calhf) THEN IF (hybrid%l_calhf) THEN
WRITE(6, '( '' ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,4X,3F15.5)')& WRITE (6, '( '' ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,4X,3F15.5)') &
kpts%bkf(:,nk),iband, (REAL(exch(iband,iband))-hybrid%div_vv(iband,nk,jsp))*(-27.211608),& kpts%bkf(:, nk), iband, (REAL(exch(iband, iband)) - hybrid%div_vv(iband, nk, jsp))*(-27.211608), &
hybrid%div_vv(iband,nk,jsp)*(-27.211608),REAL(exch(iband,iband))*(-27.211608) hybrid%div_vv(iband, nk, jsp)*(-27.211608), REAL(exch(iband, iband))*(-27.211608)
END IF END IF
END DO END DO
END SUBROUTINE add_vnonlocal END SUBROUTINE add_vnonlocal
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
...@@ -2,120 +2,120 @@ ...@@ -2,120 +2,120 @@
! [1] Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51 ! [1] Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
module m_exponential_integral module m_exponential_integral
implicit none implicit none
real, parameter :: series_laguerre = 4.0 real, parameter :: series_laguerre = 4.0
contains contains
! Calculate the exponential integral E_1(x): ! Calculate the exponential integral E_1(x):
! !
! inf ! inf
! / -t ! / -t
! | e ! | e
! E (x) = | dt ----- ! E (x) = | dt -----
! 1 | t ! 1 | t
! / ! /
! x ! x
! !
! Input: arg - position at which exponential integral is evaluated (arg > 0) ! Input: arg - position at which exponential integral is evaluated (arg > 0)
! Output: res - E_1(arg) ! Output: res - E_1(arg)
pure subroutine calculateExponentialIntegral(arg, res) pure subroutine calculateExponentialIntegral(arg, res)
implicit none implicit none
real, intent(in) :: arg real, intent(in) :: arg
real, intent(out) :: res real, intent(out) :: res
! For arguments smaller than 4 the series expansion is used ! For arguments smaller than 4 the series expansion is used
if (arg < series_laguerre) then if (arg < series_laguerre) then
res = seriesExpansion(arg) res = seriesExpansion(arg)
! otherwise a Gauss-Laguerre expansion is better ! otherwise a Gauss-Laguerre expansion is better
else else
res = exp(-arg) * gauss_laguerre(arg) res = exp(-arg)*gauss_laguerre(arg)
endif endif
end subroutine calculateExponentialIntegral end subroutine calculateExponentialIntegral
! Series expansion of the exponential integral ! Series expansion of the exponential integral
! !
! n_cut ! n_cut
! ----- n n ! ----- n n
! \ (-1) x ! \ (-1) x
! E (x) = -gamma - ln(x) - ) -------- ! E (x) = -gamma - ln(x) - ) --------
! 1 / n * n! ! 1 / n * n!
! ----- ! -----
! n = 1 ! n = 1
! !
! where gamma is the Euler constant. ! where gamma is the Euler constant.
! n_cut is set to 25 ! n_cut is set to 25
! Input: arg - argument for which the exponential integral is approximated ! Input: arg - argument for which the exponential integral is approximated
! Return: approximation by series expansion for E_1(arg) ! Return: approximation by series expansion for E_1(arg)
pure real function seriesExpansion(arg) pure real function seriesExpansion(arg)
implicit none implicit none
real, intent(in) :: arg real, intent(in) :: arg
real :: res, fact ! result of the summation, 1 / n real :: res, fact ! result of the summation, 1 / n
integer :: i ! counter variable integer :: i ! counter variable
real, parameter :: EULER_GAMMA = 0.57721566490153286060651209008241 ! Euler constant real, parameter :: EULER_GAMMA = 0.57721566490153286060651209008241 ! Euler constant
integer, parameter :: ITERATION = 25 ! Cutoff for series expansion integer, parameter :: ITERATION = 25 ! Cutoff for series expansion
! initialize summation result ! initialize summation result
res = 0.0 res = 0.0
! perform the summation ! perform the summation
do i = ITERATION, 2, -1 do i = ITERATION, 2, -1
! calculate 1/n ! calculate 1/n
fact = 1.0 / i fact = 1.0/i
! add next term of summation ! add next term of summation
res = arg * fact * (fact - res) res = arg*fact*(fact - res)
end do end do
! calculate the final result ! calculate the final result
seriesExpansion = -EULER_GAMMA - log(arg) + arg * (1.0 - res) seriesExpansion = -EULER_GAMMA - log(arg) + arg*(1.0 - res)
end function seriesExpansion end function seriesExpansion
! The Gauss Laguerre expansion of the exponential integral can be written as ! The Gauss Laguerre expansion of the exponential integral can be written as
! !
! N ! N
! E (arg) ----- a ! E (arg) ----- a
! 1 \ n ! 1 \ n
! ------- = ) -------- ! ------- = ) --------
! -arg / x + arg ! -arg / x + arg
! e ----- n ! e ----- n
! n=1 ! n=1
! !
! where the a_n and x_n are determined by least quadrature and are given in [1] ! where the a_n and x_n are determined by least quadrature and are given in [1]
! Input: arg - point at which Gaussian Laguerre quadrature is calculated ! Input: arg - point at which Gaussian Laguerre quadrature is calculated
! Return: E_1(arg) in this approximation ! Return: E_1(arg) in this approximation
pure real function gauss_laguerre(arg) pure real function gauss_laguerre(arg)
implicit none implicit none
real, intent(in) :: arg real, intent(in) :: arg
! the quadrature constants a_n and x_n from [1] ! the quadrature constants a_n and x_n from [1]
real, parameter :: a(1:15) = (/& real, parameter :: a(1:15) = (/ &
0.2182348859400869e+00, 0.3422101779228833e+00, 0.2630275779416801e+00, & 0.2182348859400869e+00, 0.3422101779228833e+00, 0.2630275779416801e+00, &
0.1264258181059305e+00, 0.4020686492100091e-01, 0.8563877803611838e-02, & 0.1264258181059305e+00, 0.4020686492100091e-01, 0.8563877803611838e-02, &
0.1212436147214252e-02, 0.1116743923442519e-03, 0.6459926762022901e-05, & 0.1212436147214252e-02, 0.1116743923442519e-03, 0.6459926762022901e-05, &
0.2226316907096273e-06, 0.4227430384979365e-08, 0.3921897267041089e-10, & 0.2226316907096273e-06, 0.4227430384979365e-08, 0.3921897267041089e-10, &
0.1456515264073126e-12, 0.1483027051113301e-15, 0.1600594906211133e-19/) 0.1456515264073126e-12, 0.1483027051113301e-15, 0.1600594906211133e-19/)
real, parameter :: x(1:15) = (/& real, parameter :: x(1:15) = (/ &
0.9330781201728180e-01, 0.4926917403018839e+00, 0.1215595412070949e+01, & 0.9330781201728180e-01, 0.4926917403018839e+00, 0.1215595412070949e+01, &
0.2269949526203743e+01, 0.3667622721751437e+01, 0.5425336627413553e+01, & 0.2269949526203743e+01, 0.3667622721751437e+01, 0.5425336627413553e+01, &
0.7565916226613068e+01, 0.1012022856801911e+02, 0.1313028248217572e+02, & 0.7565916226613068e+01, 0.1012022856801911e+02, 0.1313028248217572e+02, &
0.1665440770832996e+02, 0.2077647889944877e+02, 0.2562389422672878e+02, & 0.1665440770832996e+02, 0.2077647889944877e+02, 0.2562389422672878e+02, &
0.3140751916975394e+02, 0.3853068330648601e+02, 0.4802608557268579e+02/) 0.3140751916975394e+02, 0.3853068330648601e+02, 0.4802608557268579e+02/)
! Calculate the summation ! Calculate the summation
gauss_laguerre = sum( a / (x + arg) ) gauss_laguerre = sum(a/(x + arg))
end function gauss_laguerre end function gauss_laguerre
end module m_exponential_integral end module m_exponential_integral
...@@ -10,22 +10,22 @@ ...@@ -10,22 +10,22 @@
! and writes them out in cmt and z, respectively. ! ! and writes them out in cmt and z, respectively. !