1  IMPLICIT NONE
2  INTEGER, PARAMETER :: dp=KIND(0.0D0)
3  REAL(KIND=dp)      :: res
4
5  res=exp_radius_very_extended(  0    ,      1   ,       0      ,    1, &
6                               (/0.0D0,0.0D0,0.0D0/),&
7                               (/1.0D0,0.0D0,0.0D0/),&
8                               (/1.0D0,0.0D0,0.0D0/),&
9                                 1.0D0,1.0D0,1.0D0,1.0D0)
10  if (res.ne.1.0d0) call abort()
11
12CONTAINS
13
14 FUNCTION exp_radius_very_extended(la_min,la_max,lb_min,lb_max,ra,rb,rp,&
15                          zetp,eps,prefactor,cutoff) RESULT(radius)
16
17    INTEGER, INTENT(IN)                      :: la_min, la_max, lb_min, lb_max
18    REAL(KIND=dp), INTENT(IN)                :: ra(3), rb(3), rp(3), zetp, &
19                                                eps, prefactor, cutoff
20    REAL(KIND=dp)                            :: radius
21
22    INTEGER                                  :: i, ico, j, jco, la(3), lb(3), &
23                                                lxa, lxb, lya, lyb, lza, lzb
24    REAL(KIND=dp)                            :: bini, binj, coef(0:20), &
25                                                epsin_local, polycoef(0:60), &
26                                                prefactor_local, rad_a, &
27                                                rad_b, s1, s2
28
29    epsin_local=1.0E-2_dp
30
31    prefactor_local=prefactor*MAX(1.0_dp,cutoff)
32    rad_a=SQRT(SUM((ra-rp)**2))
33    rad_b=SQRT(SUM((rb-rp)**2))
34
35    polycoef(0:la_max+lb_max)=0.0_dp
36    DO lxa=0,la_max
37    DO lxb=0,lb_max
38       coef(0:la_max+lb_max)=0.0_dp
39       bini=1.0_dp
40       s1=1.0_dp
41       DO i=0,lxa
42          binj=1.0_dp
43          s2=1.0_dp
44          DO j=0,lxb
45             coef(lxa+lxb-i-j)=coef(lxa+lxb-i-j) + bini*binj*s1*s2
46             binj=(binj*(lxb-j))/(j+1)
47             s2=s2*(rad_b)
48          ENDDO
49          bini=(bini*(lxa-i))/(i+1)
50          s1=s1*(rad_a)
51       ENDDO
52       DO i=0,lxa+lxb
53          polycoef(i)=MAX(polycoef(i),coef(i))
54       ENDDO
55    ENDDO
56    ENDDO
57
58    polycoef(0:la_max+lb_max)=polycoef(0:la_max+lb_max)*prefactor_local
59    radius=0.0_dp
60    DO i=0,la_max+lb_max
61          radius=MAX(radius,polycoef(i)**(i+1))
62    ENDDO
63
64  END FUNCTION exp_radius_very_extended
65
66END
67