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