Line data    Source code

       1             : module m_sphbessel_integral
       2             : contains
       3     1893426 :    FUNCTION sphbessel_integral(atoms, itype, qnrm, nqnrm, iqnrm1, iqnrm2, l, hybinp, &
       4     1893426 :                                sphbes0, l_warnin, l_warnout)
       5             : 
       6             :       USE m_types
       7             :       USE m_constants
       8             : 
       9             :       IMPLICIT NONE
      10             : 
      11             :       TYPE(t_hybinp), INTENT(IN)   :: hybinp
      12             :       TYPE(t_atoms), INTENT(IN)   :: atoms
      13             : 
      14             :       INTEGER, INTENT(IN)  :: itype, nqnrm, iqnrm1, iqnrm2, l
      15             :       REAL, INTENT(IN)  :: qnrm(nqnrm), sphbes0(-1:hybinp%lexp + 2, atoms%ntype, nqnrm)
      16             :       LOGICAL, INTENT(IN), OPTIONAL  ::  l_warnin
      17             :       LOGICAL, INTENT(INOUT), OPTIONAL  ::  l_warnout
      18             :       REAL                  :: sphbessel_integral
      19             :       REAL                  :: q1, q2, dq, s, sb01, sb11, sb21, sb31, sb02, sb12
      20             :       REAL                  :: sb22, sb32, a1, a2, da, b1, b2, db, c1, c2, dc, r1, r2
      21             :       LOGICAL               :: l_warn, l_warned
      22             : 
      23     1893426 :       IF (PRESENT(l_warnin)) THEN
      24     1893426 :          l_warn = l_warnin
      25             :       ELSE
      26             :          l_warn = .TRUE.
      27             :       END IF
      28     1893426 :       l_warned = .FALSE.
      29             : 
      30     1893426 :       q1 = qnrm(iqnrm1)
      31     1893426 :       q2 = qnrm(iqnrm2)
      32     1893426 :       s = atoms%rmt(itype)
      33     1893426 :       IF (abs(q1) < 1e-12 .AND. abs(q2) < 1e-12) THEN
      34         170 :       IF (l > 0) THEN
      35             :          sphbessel_integral = 0
      36             :       ELSE
      37          10 :          sphbessel_integral = 2*s**5/15
      38             :       ENDIF
      39     1893256 :       ELSE IF (abs(q1) < 1e-12 .OR. abs(q2) < 1e-12) THEN
      40        2584 :       IF (l > 0) THEN
      41             :          sphbessel_integral = 0
      42         152 :       ELSE IF (abs(q1) < 1e-12) THEN
      43             :          sphbessel_integral = s**3/(3*q2**2)*(q2*s*sphbes0(1, itype, iqnrm2) &
      44         152 :                                               + sphbes0(2, itype, iqnrm2))
      45             :       ELSE
      46             :          sphbessel_integral = s**3/(3*q1**2)*(q1*s*sphbes0(1, itype, iqnrm1) &
      47           0 :                                               + sphbes0(2, itype, iqnrm1))
      48             :       ENDIF
      49     1890672 :       ELSE IF (abs(q1 - q2) < 1e-12) THEN
      50             :       sphbessel_integral = s**3/(2*q1**2)*((2*l + 3)*sphbes0(l + 1, itype, iqnrm1)**2 - &
      51      226780 :                                            (2*l + 1)*sphbes0(l, itype, iqnrm1)*sphbes0(l + 2, itype, iqnrm1))
      52             :       ELSE ! We use either if two fromulas that are stable for high and small q1/q2 respectively
      53     1663892 :       sb01 = sphbes0(l - 1, itype, iqnrm1)
      54     1663892 :       sb11 = sphbes0(l, itype, iqnrm1)
      55     1663892 :       sb21 = sphbes0(l + 1, itype, iqnrm1)
      56     1663892 :       sb31 = sphbes0(l + 2, itype, iqnrm1)
      57     1663892 :       sb02 = sphbes0(l - 1, itype, iqnrm2)
      58     1663892 :       sb12 = sphbes0(l, itype, iqnrm2)
      59     1663892 :       sb22 = sphbes0(l + 1, itype, iqnrm2)
      60     1663892 :       sb32 = sphbes0(l + 2, itype, iqnrm2)
      61     1663892 :       dq = q1**2 - q2**2
      62     1663892 :       a1 = q2/q1*sb21*sb02
      63     1663892 :       a2 = q1/q2*sb22*sb01
      64     1663892 :       da = a1 - a2
      65     1663892 :       b1 = sb31*sb12
      66     1663892 :       b2 = sb32*sb11
      67     1663892 :       db = b1 - b2
      68     1663892 :       c1 = sb21*sb22/(q1*q2)
      69     1663892 :       c2 = db/dq*(2*l + 1)/(2*l + 3)
      70     1663892 :       dc = c1 + c2
      71     1663892 :       r1 = ABS(da/a1)
      72     1663892 :       r2 = MIN(ABS(db/b1), ABS(dc/c1))
      73             : ! Ensure numerical stability. If both formulas are not sufficiently stable, the program stops.
      74     1663892 :       IF (r1 > r2) THEN
      75      228124 :       IF (r1 < 1e-6 .AND. l_warn) THEN
      76           0 :          WRITE (oUnit, '(A,E12.5,A,E12.5,A)') 'sphbessel_integral: Warning! Formula One possibly unstable. Ratios:', &
      77           0 :             r1, '(', r2, ')'
      78           0 :          WRITE (oUnit, '(A,2F15.10,I4)') '                    Current qnorms and atom type:', q1, q2, itype
      79           0 :          l_warned = .TRUE.
      80             :       END IF
      81      228124 :       sphbessel_integral = s**3/dq*da
      82             :       ELSE
      83     1435768 :       IF (r2 < 1e-6 .AND. l_warn) THEN
      84           0 :          WRITE (oUnit, '(A,E13.5,A,E13.5,A)') 'sphbessel_integral: Warning! Formula Two possibly unstable. Ratios:', &
      85           0 :             r2, '(', r1, ')'
      86           0 :          WRITE (oUnit, '(A,2F15.10,I4)') '                    Current qnorms and atom type:', &
      87           0 :             q1, q2, itype
      88           0 :          l_warned = .TRUE.
      89             :       END IF
      90     1435768 :       sphbessel_integral = s**3*dc
      91             :       END IF
      92             :       END IF
      93             : 
      94     1893426 :       IF (PRESENT(l_warnout)) l_warnout = l_warned
      95             : 
      96     1893426 :    END FUNCTION sphbessel_integral
      97             : end module m_sphbessel_integral