LCOV - code coverage report
Current view: top level - math - util.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 333 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 19 0.0 %

          Line data    Source code
       1             :       MODULE m_util
       2             :       USE m_juDFT
       3             : c     error and warning codes for intgrf function
       4             :       INTEGER, PARAMETER :: NO_ERROR                  = 0
       5             :       INTEGER, PARAMETER :: NEGATIVE_EXPONENT_WARNING = 1
       6             :       INTEGER, PARAMETER :: NEGATIVE_EXPONENT_ERROR   = 2
       7             : c     return type of the pure intgrf function
       8             :       TYPE :: intgrf_out
       9             :         REAL    :: value    ! value of the integration
      10             :         INTEGER :: ierror   ! error code
      11             :       END TYPE intgrf_out
      12             : 
      13             :       INTERFACE derivative
      14             :         MODULE PROCEDURE derivative_t,derivative_nt
      15             :       END INTERFACE
      16             :       
      17             :       CONTAINS
      18             : 
      19             : 
      20             : c     Calculates Gaunt coefficients, i.e. the integrals of three spherical harmonics
      21             : c     integral ( conjg(Y(l1,m1)) * Y(l2,m2) * conjg(Y(l3,m3)) )
      22             : c     They are also the coefficients C(l1,l2,l3,m1,m2,m3) in
      23             : c     conjg(Y(l1,m1)) * Y(l2,m2) = sum(l3,m3) C(l1,l2,l3,m1,m2,m3) Y(l3,m3)
      24             : c     fac contains factorial up to maxfac, i.e. fac(i)= i! 
      25             : C     sfac contains square root of fac, i.e. sfac(i)= sqrt(i!) 
      26             : 
      27           0 :       FUNCTION gaunt(l1,l2,l3,m1,m2,m3,maxfac,fac,sfac)
      28             :       
      29             :       USE m_constants ,ONLY: pimach
      30             :       
      31             :       IMPLICIT NONE
      32             :       
      33             :       REAL                :: gaunt
      34             :       INTEGER, INTENT(IN) :: l1,l2,l3,m1,m2,m3,maxfac
      35             :       REAL   , INTENT(IN) :: fac(0:maxfac)
      36             :       REAL   , INTENT(IN) :: sfac(0:maxfac)
      37             : 
      38           0 :       gaunt = 0
      39           0 :       IF(m3.ne.m2-m1)   RETURN
      40           0 :       IF(abs(m1).gt.l1) RETURN
      41           0 :       IF(abs(m2).gt.l2) RETURN
      42           0 :       IF(abs(m3).gt.l3) RETURN
      43           0 :       IF(l3.lt.abs(l1-l2).or.l3.gt.l1+l2) RETURN
      44             :       gaunt = (-1)**(m1+m3) *
      45             :      &        sqrt((2*l1+1)*(2*l2+1)*(2*l3+1)/pimach()/4)*
      46             :      &        wigner3j(l1,l2,l3,-m1,m2,-m3,maxfac,fac,sfac)*
      47           0 :      &        wigner3j(l1,l2,l3, 0,  0,  0,maxfac,fac,sfac)
      48           0 :       END FUNCTION gaunt
      49             :  
      50             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      51             : 
      52             : c     Calculates the Wigner 3j symbols using Racah's formula
      53             : 
      54           0 :       FUNCTION wigner3j(l1,l2,l3,m1,m2,m3,maxfac,fac,sfac)
      55             : 
      56             :       IMPLICIT NONE
      57             : 
      58             :       REAL                :: wigner3j
      59             :       INTEGER, INTENT(IN) :: l1,l2,l3,m1,m2,m3,maxfac
      60             :       REAL   , INTENT(IN) :: fac(0:maxfac)
      61             :       REAL   , INTENT(IN) :: sfac(0:maxfac)
      62             : c     - local - 
      63             :       INTEGER             :: tmin,tmax,t,f1,f2,f3,f4,f5
      64             : 
      65           0 :       wigner3j = 0
      66             : 
      67             : c     The following IF clauses should be in the calling routine and commented here.
      68             : c      if(-m3.ne.m1+m2)  return
      69             : c      if(abs(m1).gt.l1) return
      70             : c      if(abs(m2).gt.l2) return
      71             : c      if(abs(m3).gt.l3) return
      72             : c      if(l3.lt.abs(l1-l2).or.l3.gt.l1+l2) return
      73             : 
      74           0 :       f1   = l3-l2+m1
      75           0 :       f2   = l3-l1-m2
      76           0 :       f3   = l1+l2-l3
      77           0 :       f4   = l1-m1
      78           0 :       f5   = l2+m2
      79           0 :       tmin = max(0,-f1,-f2) ! The arguments to fac (see below)
      80           0 :       tmax = min(f3,f4,f5)  ! must not be negative.
      81             : 
      82             :       ! The following line is only for testing and should be removed at a later time.      
      83           0 :       IF(tmax-tmin .ne. min(l1+m1,l1-m1,l2+m2,l2-m2,l3+m3,l3-m3,
      84             :      &                      l1+l2-l3,l1-l2+l3,-l1+l2+l3))
      85           0 :      &  STOP 'wigner3j: Number of terms incorrect.'
      86             :       
      87           0 :       IF(tmin.le.tmax) THEN
      88           0 :         DO t = tmin,tmax
      89             :           wigner3j = wigner3j + (-1)**t /
      90             :      &                          (  fac(t)    * fac(f1+t) * fac(f2+t) 
      91           0 :      &                            *fac(f3-t) * fac(f4-t) * fac(f5-t) )
      92             :         END DO
      93             :         wigner3j = wigner3j * (-1)**(l1-l2-m3) * sfac(l1+l2-l3)
      94             :      &                      * sfac(l1-l2+l3)   * sfac(-l1+l2+l3) 
      95             :      &                      / sfac(l1+l2+l3+1) *
      96             :      &                        sfac(l1+m1)      * sfac(l1-m1) *
      97             :      &                        sfac(l2+m2)      * sfac(l2-m2) *
      98           0 :      &                        sfac(l3+m3)      * sfac(l3-m3)
      99             :       END IF
     100           0 :       END FUNCTION wigner3j
     101             : 
     102             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     103             : 
     104             : c     Integrates function f numerically (Lagrange and Simpson integration)
     105             : c     on grid(itype) and is much faster than intgr. 
     106             : c     (Only normal outward integration.)
     107             : c     Before first use of this function it has to be initialized with 
     108             : c     intgrf_init.
     109             : 
     110           0 :       FUNCTION intgrf(f,jri,jmtd,rmsh,dx,ntype,itype,gridf)
     111             : 
     112             : 
     113             :       IMPLICIT NONE
     114             : 
     115             :       REAL                 :: intgrf
     116             :       INTEGER, INTENT(IN)  :: itype,ntype,jmtd
     117             :       INTEGER, INTENT(IN)  :: jri(ntype)
     118             :       REAL   , INTENT(IN)  :: dx(ntype),rmsh(jmtd,ntype)
     119             :       REAL   , INTENT(IN)  :: gridf(jmtd,ntype)
     120             :       REAL   , INTENT(IN)  :: f(*)
     121             : c     - local -
     122             :       TYPE(intgrf_out)     :: fct_res
     123             : 
     124             : 
     125           0 :       fct_res = pure_intgrf(f,jri,jmtd,rmsh,dx,ntype,itype,gridf)
     126           0 :       IF (fct_res%ierror == NEGATIVE_EXPONENT_WARNING) THEN
     127             :         write(6,*) 'intgrf: Warning!'//
     128           0 :      +               'Negative exponent x in extrapolation a+c*r**x'
     129           0 :       ELSEIF (fct_res%ierror == NEGATIVE_EXPONENT_ERROR) THEN
     130             :         write(6,*)
     131           0 :      +       'intgrf: Negative exponent x in extrapolation a+c*r**x'
     132             :         CALL juDFT_error(
     133           0 :      +          'intgrf: Negative exponent x in extrapolation a+c*r**x')
     134             :       END IF
     135           0 :       intgrf = fct_res%value
     136             : 
     137           0 :       END FUNCTION intgrf
     138             : 
     139             : c     pure wrapper for intgrf with same functionality
     140             : c     can be used within forall loops
     141             : 
     142           0 :       PURE FUNCTION pure_intgrf(f,jri,jmtd,rmsh,dx,ntype,itype,gridf)
     143             : 
     144             :       IMPLICIT NONE
     145             : 
     146             :       TYPE(intgrf_out)     :: pure_intgrf
     147             : 
     148             :       INTEGER, INTENT(IN)  :: itype,ntype,jmtd
     149             :       INTEGER, INTENT(IN)  :: jri(ntype)
     150             :       REAL   , INTENT(IN)  :: dx(ntype),rmsh(jmtd,ntype)
     151             :       REAL   , INTENT(IN)  :: gridf(jmtd,ntype)
     152             :       REAL   , INTENT(IN)  :: f(*)
     153             : c     - local -
     154             :       INTEGER              :: n
     155             :       REAL                 :: r1,h,a,x
     156             : 
     157           0 :       n  = jri(itype)
     158           0 :       r1 = rmsh(1,itype)
     159           0 :       h  = dx(itype)
     160             : 
     161           0 :       pure_intgrf%ierror = NO_ERROR
     162             : 
     163             :       ! integral from 0 to r1 approximated by leading term in power series expansion
     164           0 :       IF (f(1)*f(2).gt.1e-10.and.h.gt.0) THEN
     165           0 :         IF(f(2).eq.f(1)) THEN
     166           0 :           pure_intgrf%value = r1*f(1)
     167             :         ELSE
     168           0 :           x      = (f(3)-f(2)) / (f(2)-f(1))
     169           0 :           a      = (f(2)-x*f(1)) / (1-x)
     170           0 :           x      = log(x)/h
     171           0 :           IF(x.lt.0) THEN
     172           0 :             IF(x.gt.-1) THEN
     173             :               pure_intgrf%ierror = NEGATIVE_EXPONENT_WARNING
     174           0 :             ELSE IF(x.le.-1) THEN
     175           0 :               pure_intgrf%ierror = NEGATIVE_EXPONENT_ERROR
     176           0 :               RETURN
     177             :             END IF
     178             :           END IF
     179             : 
     180           0 :           pure_intgrf%value = r1*(f(1)+x*a) / (x+1)
     181             : 
     182             : !           x      = f(2) / f(1)
     183             : !           x      = log(x)/h
     184             : !           IF(x.lt.0) THEN
     185             : !             IF(x.gt.-1) write(6,'(A,ES9.1)') 'intgrf: Warning!&
     186             : !      &                                        Negative exponent x in& 
     187             : !      &                                        extrapolation c*r**x:',x
     188             : !             IF(x.le.-1) write(6,'(A,ES9.1)') 'intgrf: Negative exponent&
     189             : !      &                                        x in extrapolation&
     190             : !      &                                        c*r**x:',x
     191             : !             IF(x.le.-1) STOP                 'intgrf: Negative exponent& 
     192             : !      &                                        x in extrapolation &
     193             : !      &                                        c*r**x'
     194             : !           END IF
     195             : !           intgrf = (r1*f(1))/(x+1)
     196             : 
     197             :         END IF
     198             :       ELSE
     199             :         pure_intgrf%value = 0
     200             :       END IF 
     201             :       
     202             :       ! integrate from r(1) to r(n) by multiplying with gridf
     203             :       pure_intgrf%value = pure_intgrf%value 
     204           0 :      +                  + dot_product(gridf(:n,itype),f(:n))
     205             :       
     206           0 :       END FUNCTION pure_intgrf
     207             : 
     208             : 
     209             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     210             : 
     211             : c     Initializes fast numerical integration intgrf (see below).
     212             :       
     213           0 :       SUBROUTINE intgrf_init(ntype,jmtd,jri,dx,rmsh,gridf)
     214             : 
     215             :       IMPLICIT NONE
     216             : 
     217             :       INTEGER, INTENT(IN)   :: ntype,jmtd
     218             :       INTEGER, INTENT(IN)   :: jri(ntype)
     219             :       REAL,    INTENT(IN)   :: dx(ntype),rmsh(jmtd,ntype)
     220             :       REAL, INTENT(OUT), ALLOCATABLE  :: gridf(:,:)
     221             : 
     222             : c     - local -
     223             :       INTEGER              :: i,j,itype
     224             :       INTEGER              :: n,nstep,n0 = 6
     225             :       INTEGER, PARAMETER   :: simpson(7) = (/41,216,27,272,27,216,41/)
     226             :       REAL                 :: r1,h,dr
     227             :       REAL                 :: r(7)     
     228             :       REAL, PARAMETER      :: lagrange(7,6)= reshape(
     229             :      &          (/19087.,65112.,-46461., 37504.,-20211., 6312., -863.,
     230             :      &             -863.,25128., 46989.,-16256.,  7299.,-2088.,  271.,
     231             :      &              271.,-2760., 30819., 37504., -6771., 1608., -191.,
     232             :      &             -191., 1608., -6771., 37504., 30819.,-2760.,  271.,
     233             :      &              271.,-2088.,  7299.,-16256., 46989.,25128., -863.,
     234             :      &             -863., 6312.,-20211., 37504.,-46461.,65112.,19087./), ! The last row is actually never used.
     235             :      &                                      (/7,6/) )
     236             : 
     237           0 :       n = jmtd
     238           0 :       ALLOCATE ( gridf(n,ntype) )
     239             : 
     240           0 :       gridf = 0
     241             : 
     242           0 :       DO itype = 1,ntype
     243             : 
     244           0 :         n  = jri(itype)
     245           0 :         r1 = rmsh(1,itype)
     246           0 :         h  = dx(itype)
     247             : 
     248           0 :         nstep = (n-1)/6
     249           0 :         n0    = n-6*nstep
     250           0 :         dr    = exp(h)
     251             : 
     252             :         ! Calculate Lagrange-integration coefficients from r(1) to r(n0)
     253           0 :         r(1)=r1
     254           0 :         IF(n0.gt.1) THEN
     255           0 :           DO i=2,7
     256           0 :             r(i) = r(i-1)*dr
     257             :           END DO
     258           0 :           DO i=1,7
     259           0 :             gridf(i,itype) = h/60480 * r(i) * sum(lagrange(i,1:n0-1))
     260             :           END DO
     261           0 :           r(1)  = r(n0)
     262             :         END IF
     263             : 
     264             :         ! Calculate Simpson-integration coefficients from r(n0) to r(n)
     265           0 :         DO i = 1,nstep
     266           0 :           DO j = 2,7
     267           0 :             r(j) = r(j-1)*dr
     268             :           END DO
     269           0 :           DO j = n0,n0+6
     270             :             gridf(j,itype) = gridf(j,itype) + h/140 * r(j-n0+1) *
     271           0 :      &                       simpson(j-n0+1)
     272             :           END DO
     273           0 :           n0   = n0 + 6
     274           0 :           r(1) = r(7)
     275             :         END DO
     276             : 
     277             :       END DO
     278             :         
     279           0 :       END SUBROUTINE intgrf_init
     280             : 
     281             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     282             : 
     283             : c     Calculates the primitive of f, on grid(itypein):
     284             : c
     285             : c                   r
     286             : c     primf(r) = integral f(r') dr'   ( r = grid point )
     287             : c                   0
     288             : c
     289             : c     If itypein is negative, the primitive
     290             : c     
     291             : c                   R
     292             : c     primf(r) = integral f(r') dr'   ( R = MT sphere radius )
     293             : c                   r
     294             : c
     295             : c     is calculated instead.
     296             : c
     297             : 
     298             : c     -----------------------------
     299             : 
     300             : c     Fast calculation of primitive.
     301             : c     Only Lagrange integration is used
     302             : 
     303           0 :       SUBROUTINE primitivef(primf,fin,rmsh,dx,jri,jmtd,itypein,ntype)
     304             : 
     305             :       IMPLICIT NONE
     306             : 
     307             : c     - scalars -
     308             :       INTEGER, INTENT(IN)   :: itypein,jmtd,ntype
     309             : 
     310             : c     - arrays -
     311             :       INTEGER, INTENT(IN)   :: jri(ntype)
     312             :       REAL,    INTENT(OUT)  :: primf( jri( abs(itypein) ) )
     313             :       REAL,    INTENT(IN)   :: fin( jri( abs(itypein) ) )
     314             :       REAL,    INTENT(IN)   :: rmsh(jmtd,ntype),dx(ntype)
     315             : 
     316             : c     - local scalars -
     317             :       INTEGER               :: itype,n,i,n0
     318             :       REAL                  :: h,x,h1
     319             :       REAL                  :: intgr,r1,a,dr
     320             : 
     321             : c     - local arrays -
     322             :       REAL                  :: fr(7)
     323           0 :       REAL                  :: f( jri( abs(itypein) ))
     324           0 :       REAL                  :: r( jri( abs(itypein) ) )
     325             :       REAL, PARAMETER       :: lagrange(7,6)= reshape(
     326             :      &          (/19087.,65112.,-46461., 37504.,-20211., 6312., -863.,
     327             :      &             -863.,25128., 46989.,-16256.,  7299.,-2088.,  271.,
     328             :      &              271.,-2760., 30819., 37504., -6771., 1608., -191.,
     329             :      &             -191., 1608., -6771., 37504., 30819.,-2760.,  271.,
     330             :      &              271.,-2088.,  7299.,-16256., 46989.,25128., -863.,
     331             :      &             -863., 6312.,-20211., 37504.,-46461.,65112.,19087./),
     332             :      &                                       (/7,6/) )
     333             :     
     334             : 
     335           0 :       itype = abs(itypein)
     336             : 
     337           0 :       primf = 0
     338             : 
     339           0 :       n = jri(itype)
     340           0 :       h = dx(itype)
     341             : 
     342           0 :       IF(itypein.gt.0) THEN
     343           0 :         r1 = rmsh(1,itype)      ! perform outward integration
     344           0 :         f  = fin                ! (from 0 to r)
     345             :       ELSE
     346           0 :         r1 = rmsh(jri(itype),itype)         ! perform inward integration
     347           0 :         h  = -h                             ! (from MT sphere radius to r)
     348           0 :         f  = fin(n:1:-1)                    !
     349             :       END IF
     350             :       
     351             :       ! integral from 0 to r1 approximated by leading term in power series expansion (only if h>0)
     352           0 :       IF(h.gt.0.and.f(1)*f(2).gt.1e-10) THEN
     353           0 :         IF(f(2).eq.f(1)) THEN
     354           0 :           intgr = r1*f(1)
     355             :         ELSE
     356           0 :           x     = (f(3)-f(2))/(f(2)-f(1))
     357           0 :           a     = (f(2)-x*f(1)) / (1-x)
     358           0 :           x     = log(x)/h
     359           0 :           IF(x.lt.0) THEN
     360           0 :              IF(x>-1) WRITE(6,'(A,ES9.1)')
     361             :      +            '+intgr: Warning! Negative &exponent x in'//
     362           0 :      +        'extrapolation a+c*r**x:',x
     363           0 :              IF(x<=-1) WRITE(6,'(A,ES9.1)')'intgr: Negative exponent,'//
     364           0 :      +             'x in extrapolation a+c*r**x:',x
     365           0 :              IF(x<=-1) STOP 'intgr:Negative exponent x in extrapolation'
     366             :           END IF
     367           0 :           intgr = r1*(f(1)+x*a) / (x+1)
     368             :         END IF
     369             :       ELSE
     370             :         intgr = 0
     371             :       END IF
     372             : 
     373           0 :       primf(1) = intgr
     374           0 :       dr       = exp(h)
     375           0 :       r(1)     = r1
     376           0 :       n0       = 1
     377           0 :       h1       = h/60480
     378             :       
     379             :       ! Lagrange integration from r(n0) to r(n0+5)
     380           0 :  1    DO i=2,7
     381           0 :         r(i) = r(i-1)*dr
     382             :       END DO
     383           0 :       fr = f(n0:n0+6) * r(:7)
     384           0 :       DO i=1,6
     385           0 :         intgr = intgr + h1 * dot_product(lagrange(:,i),fr)
     386           0 :         IF(primf(n0+i).eq.0) primf(n0+i) = intgr ! avoid double-definition
     387             :       END DO
     388           0 :       IF(n0+12.le.n) THEN
     389           0 :         r(1) = r(7)
     390           0 :         n0   = n0 + 6
     391           0 :         GOTO 1
     392           0 :       ELSE IF(n0+6.lt.n) THEN
     393           0 :         r(1)  = r(n-5-n0)
     394           0 :         n0    = n-6
     395           0 :         intgr = primf(n-6)
     396           0 :         GOTO 1
     397             :       END IF
     398             : 
     399           0 :       IF(itypein.lt.0) THEN    !
     400           0 :         primf = -primf(n:1:-1) ! Inward integration
     401             :       END IF                    !
     402             : 
     403           0 :       END SUBROUTINE primitivef
     404             : 
     405             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     406             : 
     407             :       !function modulo1 maps kpoint into first BZ
     408           0 :       FUNCTION modulo1(kpoint,nkpt3)
     409             : 
     410             :       IMPLICIT NONE
     411             : 
     412             :       INTEGER,INTENT(IN)  :: nkpt3(3)
     413             :       REAL   ,INTENT(IN)  :: kpoint(3)
     414             :       REAL                :: modulo1(3)
     415             :       INTEGER             :: help(3)
     416             : 
     417           0 :       modulo1 = kpoint*nkpt3
     418           0 :       help    = nint(modulo1)
     419           0 :       IF(any(abs(help-modulo1).gt.1e-10)) THEN
     420           0 :         WRITE(6,'(A,F5.3,2('','',F5.3),A)')'modulo1: argument (',kpoint,
     421           0 :      +                         ') is not an element of the k-point set.'
     422             :         CALL juDFT_error(
     423             :      +     'modulo1: argument not an element of k-point set.', 
     424           0 :      +     calledby = 'util:modulo1')
     425             :       END IF
     426           0 :       modulo1 = modulo(help,nkpt3)*1.0/nkpt3
     427             : 
     428             :       END FUNCTION modulo1
     429             :       
     430             : 
     431             : 
     432             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
     433             : 
     434             : c     Returns derivative of f in df.
     435             : 
     436           0 :       SUBROUTINE derivative_t(df,f,atoms,itype)
     437             :       USE m_types
     438             :       IMPLICIT NONE
     439             :       REAL,   INTENT(IN)   ::   f(:)
     440             :       REAL,   INTENT(OUT)  ::   df(:)
     441             :       TYPE(t_atoms),INTENT(IN)::atoms
     442             :       INTEGER,INTENT(IN)   ::   itype
     443             : 
     444             :       call derivative_nt(df,f,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,
     445           0 :      +     atoms%ntype,itype)
     446           0 :       END SUBROUTINE
     447             :       
     448           0 :       SUBROUTINE derivative_nt(df,f,jmtd,jri,dx,rmsh,ntype,itype)
     449             : 
     450             :       IMPLICIT NONE
     451             : 
     452             :       INTEGER,INTENT(IN)   :: ntype,itype,jmtd
     453             :       INTEGER,INTENT(IN)   :: jri(ntype)
     454             :       REAL,   INTENT(IN)   :: dx(ntype),rmsh(jmtd,ntype)
     455             :       REAL,   INTENT(IN)   :: f(jri(itype))
     456             :       REAL,   INTENT(OUT)  :: df(jri(itype))
     457             :       REAL                 :: x,h,r,d21,d32,d43,d31,d42,d41,df1,df2
     458             :       INTEGER              :: i,n
     459             :       
     460           0 :       n = jri(itype)
     461           0 :       h = dx(itype)
     462           0 :       r = rmsh(1,itype)
     463             :       ! use power series expansion a+c**x for first point
     464           0 :       IF(f(2).eq.f(1)) THEN
     465           0 :         df(1) = 0.0
     466             :       ELSE
     467           0 :         x     = (f(3)-f(2))/(f(2)-f(1))
     468           0 :         df(1) = (f(2)-f(1)) / (x-1) * log(x)/h / r
     469             :       END IF
     470             :       ! use Lagrange interpolation of 3rd order for all other points (and averaging)
     471           0 :       d21 = r * (exp(h)-1) ; d32 = d21 * exp(h) ; d43 = d32 * exp(h)
     472           0 :       d31 = d21 + d32        ; d42 = d32 + d43
     473           0 :       d41 = d31 + d43
     474             :       df(2) = - d32*d42 / (d21*d31*d41)         * f(1) 
     475             :      &        + ( 1.0/d21 - 1.0/d32 - 1.0/d42)  * f(2)
     476             :      &        + d21*d42 / (d31*d32*d43)         * f(3)
     477           0 :      &        - d21*d32 / (d41*d42*d43)         * f(4)
     478             :       df1   =   d32*d43 / (d21*d31*d41)         * f(1) 
     479             :      &        - d31*d43 / (d21*d32*d42)         * f(2)
     480             :      &        + ( 1.0/d31 + 1.0/d32 - 1.0/d43 ) * f(3)
     481           0 :      &        + d31*d32 / (d41*d42*d43)         * f(4)
     482           0 :       DO i = 3,n-2
     483           0 :         d21 = d32 ; d32 = d43 ; d43 = d43 * exp(h)
     484           0 :         d31 = d42 ; d42 = d42 * exp(h)
     485           0 :         d41 = d41 * exp(h)
     486             :         df2   = - d32*d42 / (d21*d31*d41)         * f(i-1) 
     487             :      &          + ( 1.0/d21 - 1.0/d32 - 1.0/d42)  * f(i)
     488             :      &          + d21*d42 / (d31*d32*d43)         * f(i+1)
     489           0 :      &          - d21*d32 / (d41*d42*d43)         * f(i+2)
     490           0 :         df(i) = ( df1 + df2 ) / 2
     491             :         df1   =   d32*d43 / (d21*d31*d41)         * f(i-1) 
     492             :      &          - d31*d43 / (d21*d32*d42)         * f(i)
     493             :      &          + ( 1.0/d31 + 1.0/d32 - 1.0/d43 ) * f(i+1)
     494           0 :      &          + d31*d32 / (d41*d42*d43)         * f(i+2)
     495             :       END DO
     496           0 :       df(n-1) = df1
     497             :       df(n)   = - d42*d43 / (d21*d31*d41)         * f(n-3) 
     498             :      &          + d41*d43 / (d21*d32*d42)         * f(n-2) 
     499             :      &          - d41*d42 / (d31*d32*d43)         * f(n-1) 
     500           0 :      &          + ( 1.0/d41 + 1.0/d42 + 1.0/d43 ) * f(n)
     501             : 
     502           0 :       END SUBROUTINE derivative_nt
     503             : 
     504             : 
     505             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     506             : c     Orders iarr(1:n) according to size and returns a correspondingly defined pointer in pnt
     507             : 
     508           0 :       SUBROUTINE iorderp(pnt,iarr,n)
     509             : 
     510             :       IMPLICIT NONE
     511             : 
     512             :       INTEGER, INTENT(IN)  :: n
     513             :       INTEGER, INTENT(OUT) :: pnt(1:n)
     514             :       INTEGER, INTENT(IN)  :: iarr(1:n)
     515             :       INTEGER              :: i,j,k
     516             : 
     517           0 :       DO i=1,n
     518           0 :         pnt(i) = i
     519           0 :         DO j=1,i-1
     520           0 :           IF(iarr(pnt(j)).gt.iarr(i)) THEN
     521           0 :             DO k=i,j+1,-1
     522           0 :               pnt(k) = pnt(k-1)
     523             :             END DO
     524           0 :             pnt(j) = i
     525           0 :             EXIT
     526             :           END IF
     527             :         END DO
     528             :       END DO
     529             : 
     530           0 :       END SUBROUTINE iorderp
     531             : 
     532             : c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     533             : c     Orders rarr(1:n) according to size and returns a correspondingly defined pointer in pnt
     534             : 
     535           0 :       SUBROUTINE rorderp(pnt,rarr,n)
     536             : 
     537             :       IMPLICIT NONE
     538             : 
     539             :       INTEGER, INTENT(IN)  :: n
     540             :       INTEGER, INTENT(OUT) :: pnt(1:n)
     541             :       REAL,    INTENT(IN)  :: rarr(1:n)
     542             :       INTEGER              :: i,j,k
     543             : 
     544           0 :       DO i=1,n
     545           0 :         pnt(i) = i
     546           0 :         DO j=1,i-1
     547           0 :           IF(rarr(pnt(j)).gt.rarr(i)) THEN
     548           0 :             DO k=i,j+1,-1
     549           0 :               pnt(k) = pnt(k-1)
     550             :             END DO
     551           0 :             pnt(j) = i
     552           0 :             EXIT
     553             :           END IF
     554             :         END DO
     555             :       END DO
     556             : 
     557           0 :       END SUBROUTINE rorderp
     558             : 
     559             : 
     560             : c     Same as rorderp but divides the problem in halves np times (leading to 2**np intervals) and is
     561             : c     much faster than rorderp (devide and conquer algorithm).
     562             : c     There is an optimal np, while for larger np the overhead (also memory-wise) outweights the speed-up.
     563             : c     np = max(0,int(log(n*0.001)/log(2.0))) should be a safe choice.
     564           0 :       RECURSIVE SUBROUTINE rorderpf(pnt,rarr,n,np)
     565             : 
     566             :       IMPLICIT NONE
     567             : 
     568             :       INTEGER, INTENT(IN)  :: n,np
     569             :       INTEGER, INTENT(OUT) :: pnt(n)
     570             :       REAL   , INTENT(IN)  :: rarr(n)
     571           0 :       REAL                 :: rarr1(n)
     572             :       REAL                 :: ravg
     573           0 :       INTEGER              :: pnt1(n),pnt2(n)
     574             :       INTEGER              :: n1,n2,i
     575             :       
     576           0 :       IF(np.eq.0) THEN
     577           0 :         CALL rorderp(pnt,rarr,n)
     578             :         RETURN
     579           0 :       ELSE IF(np.lt.0) THEN
     580           0 :         STOP 'rorderpf: fourth argument must be non-negative (bug?).'
     581             :       END IF
     582           0 :       ravg = sum(rarr)/n
     583             :       ! first half
     584           0 :       n1 = 0
     585           0 :       DO i = 1,n
     586           0 :         IF(rarr(i).le.ravg) THEN
     587           0 :           n1        = n1 + 1
     588           0 :           rarr1(n1) = rarr(i)
     589           0 :           pnt1(n1)  = i
     590             :         END IF
     591             :       END DO
     592           0 :       CALL rorderpf(pnt2,rarr1,n1,np-1)
     593           0 :       pnt(:n1) = pnt1(pnt2(:n1))
     594             :       ! second half
     595           0 :       n2 = 0
     596           0 :       DO i = 1,n
     597           0 :         IF(rarr(i).gt.ravg) THEN
     598           0 :           n2        = n2 + 1
     599           0 :           rarr1(n2) = rarr(i)
     600           0 :           pnt1(n2)  = i
     601             :         END IF
     602             :       END DO
     603           0 :       CALL rorderpf(pnt2,rarr1,n2,np-1)
     604           0 :       pnt(n1+1:) = pnt1(pnt2(:n2))
     605             :       END  SUBROUTINE rorderpf
     606             : 
     607             : 
     608             : c     Calculates the spherical Bessel functions of orders 0 to l at x
     609             : c     by backward recurrence using j_l(x) = (2l+3)/x j_l+1(x) - j_l+2(x) .
     610             : c     (Starting points are calculated according to Zhang, Min, 
     611             : c      "Computation of Special Functions".)
     612             : 
     613           0 :       pure SUBROUTINE sphbessel(sphbes,x,l)
     614             : 
     615             :       IMPLICIT NONE
     616             : 
     617             :       INTEGER, INTENT(IN)    :: l
     618             :       REAL   , INTENT(IN)    :: x
     619             :       REAL   , INTENT(INOUT) :: sphbes(0:l)
     620             :       REAL                   :: s0,s1,f,f0,f1,cs
     621             :       INTEGER                :: ll,lsta,lmax,msta2
     622             : 
     623             : !      IF( x.lt.0 ) THEN
     624             : !        STOP 'sphbes: negative argument (bug?).' 
     625             : !      ELSE 
     626           0 :       IF( x.eq.0 ) THEN
     627           0 :         sphbes(0)  = 1.0
     628           0 :         DO ll = 1,l
     629           0 :           sphbes(ll) = 0.0
     630             :         END DO
     631             :         RETURN
     632             :       ENDIF
     633           0 :       sphbes(0) = sin(x) / x
     634           0 :       IF( l .eq. 0 ) RETURN
     635           0 :       sphbes(1) = ( sphbes(0) - cos(x) ) / x
     636             : !       IF(l.le.1) RETURN
     637           0 :       s0   = sphbes(0)
     638           0 :       s1   = sphbes(1)
     639           0 :       lsta = lsta1(x,200)      !
     640           0 :       lmax = l                 !
     641           0 :       IF(lsta.lt.l) THEN       !
     642           0 :         lmax            = lsta ! determine starting point lsta
     643           0 :         sphbes(lmax+1:) = 0.0  ! for backward recurrence
     644             :       ELSE                     !
     645           0 :         lsta = lsta2(x,l,15)   ! 
     646             :       END IF                    !
     647           0 :       f0 = 0.0                                                      !
     648           0 :       f1 = 1e-100                                                   !
     649           0 :       DO ll = lsta,0,-1                                             ! backward recurrence
     650           0 :         f  = f1 / x * (2*ll+3) - f0 ; IF(ll.le.lmax) sphbes(ll) = f ! with arbitrary start values
     651           0 :         f0 = f1                                                     !
     652           0 :         f1 = f                                                      !
     653             :       END DO                                                        !
     654           0 :       IF(abs(s0).gt.abs(s1)) THEN ; cs = s0 / f   !
     655           0 :       ELSE                        ; cs = s1 / f0  ! scale to correct values
     656             :       END IF                                      ! 
     657           0 :       sphbes = cs * sphbes                        !
     658             : 
     659             :       CONTAINS
     660             : 
     661             : c     Test starting point
     662             :       PURE FUNCTION lsta0(x,mp)
     663             : 
     664             :       IMPLICIT NONE
     665             : 
     666             :       INTEGER             :: lsta0
     667             :       INTEGER, INTENT(IN) :: mp
     668             :       REAL   , INTENT(IN) :: x
     669             :       REAL                :: f,lgx
     670             : 
     671             :       lgx   = log10(x)
     672             :       lsta0 = 0
     673             :       f     = lgx
     674             :       DO WHILE(f.gt.-mp)
     675             :         lsta0 = lsta0 + 1
     676             :         f     = f + lgx - log10(2.0*lsta0+1)
     677             :       END DO
     678             :       END FUNCTION lsta0
     679             : 
     680             : c     Returns starting point lsta1 for backward recurrence such that sphbes(lsta1) approx. 10^(-mp).
     681           0 :       PURE FUNCTION lsta1(x,mp)
     682             : 
     683             :       IMPLICIT NONE
     684             : 
     685             :       INTEGER             :: lsta1
     686             :       INTEGER, INTENT(IN) :: mp
     687             :       REAL   , INTENT(IN) :: x
     688             :       REAL                :: f0,f1,f
     689             :       INTEGER             :: n0,n1,nn,it
     690             : 
     691           0 :       n0 = int(1.1*x) + 1
     692           0 :       f0 = envj(n0,x) - mp
     693           0 :       n1 = n0 + 5
     694           0 :       f1 = envj(n1,x) - mp
     695           0 :       DO it = 1,20
     696           0 :         nn = n1 - (n1-n0) / (1.0-f0/f1)
     697           0 :         f  = envj(nn,x) - mp
     698           0 :         IF(abs(nn-n1).lt.1) EXIT
     699           0 :         n0 = n1
     700           0 :         f0 = f1
     701           0 :         n1 = nn
     702           0 :         f1 = f
     703             :       END DO
     704           0 :       lsta1 = nn
     705           0 :       END FUNCTION lsta1
     706             : 
     707             : c     Returns the starting point lsta2 for backward recurrence such that all sphbes(l) have mp significant digits. 
     708           0 :       PURE FUNCTION lsta2(x,l,mp)
     709             : 
     710             :       IMPLICIT NONE
     711             : 
     712             :       INTEGER             :: lsta2
     713             :       INTEGER, INTENT(IN) :: l,mp
     714             :       REAL   , INTENT(IN) :: x
     715             :       REAL                :: f0,f1,f,hmp,ejn,obj
     716             :       INTEGER             :: n0,n1,nn,it
     717             : 
     718           0 :       hmp = 0.5 * mp
     719           0 :       ejn = envj(l,x)
     720           0 :       IF( ejn.le.hmp ) THEN
     721           0 :         obj = mp
     722           0 :         n0  = int(1.1*x) + 1
     723             :       ELSE
     724           0 :         obj = hmp + ejn
     725           0 :         n0  = l
     726             :       END IF
     727           0 :       f0 = envj(n0,x) - obj
     728           0 :       n1 = n0 + 5
     729           0 :       f1 = envj(n1,x) - obj
     730           0 :       DO it = 1,20
     731           0 :         nn = n1 - (n1-n0) / (10-f0/f1)
     732           0 :         f  = envj(nn,x) - obj
     733           0 :         IF(abs(nn-n1).lt.1) EXIT
     734           0 :         n0 = n1
     735           0 :         f0 = f1
     736           0 :         n1 = nn
     737           0 :         f1 = f
     738             :       END DO
     739           0 :       lsta2 = nn + 10
     740           0 :       END FUNCTION lsta2
     741             : 
     742             : 
     743           0 :       PURE FUNCTION envj(l,x)
     744             : 
     745             :       IMPLICIT NONE
     746             : 
     747             :       REAL                :: envj
     748             :       REAL   , INTENT(IN) :: x
     749             :       INTEGER, INTENT(IN) :: l
     750             : 
     751           0 :       envj = 0.5 * log10(6.28*l) - l*log10(1.36*x/l)
     752             : 
     753           0 :       END FUNCTION envj
     754             : 
     755             :       END SUBROUTINE sphbessel
     756             : 
     757             : 
     758             : 
     759             : c     Returns the spherical harmonics Y_lm(^rvec)
     760             : c     for l = 0,...,ll in Y(1,...,(ll+1)**2).
     761             : 
     762           0 :       PURE SUBROUTINE harmonicsr(Y,rvec,ll)
     763             : 
     764             :       IMPLICIT NONE
     765             : 
     766             :       integer , intent(in)  :: ll
     767             :       real    , intent(in)  :: rvec(3)
     768             :       complex , intent(out) :: Y((ll+1)**2)
     769             :       complex               :: c
     770             :       real                  :: stheta,ctheta,sphi,cphi,r,rvec1(3)
     771             :       integer               :: l,m,lm
     772             :       complex , parameter   :: img = (0.0,1.0)
     773             :  
     774           0 :       Y(1) = 0.282094791773878
     775           0 :       IF(ll.eq.0) RETURN
     776             : 
     777             :       stheta = 0
     778             :       ctheta = 0
     779             :       sphi   = 0
     780             :       cphi   = 0
     781           0 :       r      = sqrt(sum(rvec**2))
     782           0 :       IF(r.gt.1e-16) THEN
     783           0 :         rvec1  = rvec / r
     784           0 :         ctheta = rvec1(3)
     785           0 :         stheta = sqrt(rvec1(1)**2+rvec1(2)**2)
     786           0 :         IF(stheta.gt.1e-16) THEN
     787           0 :           cphi = rvec1(1) / stheta
     788           0 :           sphi = rvec1(2) / stheta
     789             :         END IF
     790             :       ELSE
     791           0 :         Y(2:) = 0.0
     792             :         RETURN
     793             :       END IF
     794             : 
     795             : c     define Y,l,-l and Y,l,l
     796           0 :       r = Y(1)
     797           0 :       c = 1
     798           0 :       DO l=1,ll
     799           0 :         r           = r*stheta*sqrt(1.0+1.0/(2*l))
     800           0 :         c           = c * (cphi + img*sphi)
     801           0 :         Y(l**2+1)   = r*conjg(c)  ! l,-l
     802           0 :         Y((l+1)**2) = r*c*(-1)**l ! l,l
     803             :       END DO
     804             : 
     805             : c     define Y,l,-l+1 and Y,l,l-1
     806           0 :       Y(3) = 0.48860251190292*ctheta
     807           0 :       DO l=2,ll
     808           0 :         r          = sqrt(2.0*l+1) * ctheta
     809           0 :         Y(l**2+2)  = r*Y((l-1)**2+1) ! l,-l+1
     810           0 :         Y(l*(l+2)) = r*Y(l**2)       ! l,l-1
     811             :       END DO
     812             : 
     813             : c     define Y,l,m, |m|<l-1
     814           0 :       DO l=2,ll
     815           0 :         lm = l**2 + 2
     816           0 :         DO m=-l+2,l-2
     817           0 :           lm = lm + 1
     818             :           Y(lm) = sqrt((2.0*l+1)/(l+m)/(l-m)) * (
     819             :      &            sqrt(2.0*l-1)*ctheta           *Y(lm-2*l)-
     820           0 :      &            sqrt((l+m-1.0)*(l-m-1)/(2*l-3))*Y(lm-4*l+2) )
     821             :         END DO
     822             :       END DO
     823             : 
     824             :       END SUBROUTINE harmonicsr
     825             : 
     826             : c     Returns the complex error function.
     827           0 :       FUNCTION cerf(z)
     828             :       
     829             :       USE m_constants , ONLY: pimach
     830             :       
     831             :       IMPLICIT NONE
     832             :       COMPLEX, INTENT(IN) ::  z
     833             :       COMPLEX             ::  cerf
     834             :       COMPLEX             ::  z1,z2,c,d,delta
     835             :       REAL                ::  pi
     836             :       INTEGER             ::  i
     837             :       
     838             :       
     839           0 :       pi = pimach()
     840           0 :       z1 = z
     841           0 :       IF(real(z).lt.0) z1 = -z1
     842           0 :       IF(real(z1).lt.2.0) THEN ! McLaurin series
     843           0 :         z2   = z1**2
     844           0 :         i    = 0
     845           0 :         c    = z1
     846           0 :         cerf = z1
     847             :         DO
     848           0 :           i    = i + 1
     849           0 :           c    = -c * z2 / i
     850           0 :           cerf = cerf + c/(2*i+1)
     851           0 :           IF(abs(c/(2*i+1)).lt.1e-20) EXIT
     852             :         END DO
     853           0 :         cerf = cerf * 2/sqrt(pi)
     854             :       ELSE ! continued fraction using Lentz's method
     855             :         d    = 0.0
     856             :         c    = z1
     857             :         cerf = z1
     858             :         i    = 0
     859             :         DO
     860           0 :           i     = i + 1
     861           0 :           c     =   2*z1 + i/c
     862           0 :           d     = ( 2*z1 + i*d )**(-1)
     863           0 :           delta = c*d
     864           0 :           cerf  = cerf * delta
     865           0 :           IF(abs(1-delta).lt.1e-15) EXIT
     866           0 :           i     = i + 1
     867           0 :           c     =   z1 + i/c
     868           0 :           d     = ( z1 + i*d )**(-1)
     869           0 :           delta = c*d
     870           0 :           cerf  = cerf * delta
     871           0 :           IF(abs(1-delta).lt.1e-15) EXIT
     872           0 :           IF(i.eq.10000) 
     873           0 :      &        STOP 'cerf: Lentz method not converged after 10000 steps.'
     874             :         END DO
     875           0 :         cerf = 1 - exp(-z1**2) / cerf / sqrt(pi)        
     876             :       END IF
     877           0 :       IF(real(z).lt.0) cerf = -cerf
     878           0 :       END FUNCTION cerf
     879             : 
     880           0 :       FUNCTION chr(int)
     881             :       
     882             :       IMPLICIT NONE
     883             :       
     884             :       CHARACTER(5) :: chr
     885             :       INTEGER      :: int
     886             :       
     887           0 :       WRITE(chr,'(I5)') int
     888           0 :       END FUNCTION chr
     889             : 
     890             :       END MODULE m_util

Generated by: LCOV version 1.13