1C $Id$
2************************************************************************
3*                                                                      *
4      logical function ecp_skipint(
5     &    l_a,n_prim_a,n_cont_a,coef_a,zeta_a,R_AC,
6     &    l_b,n_prim_b,n_cont_b,coef_b,zeta_b,R_BC,
7     &    n_coef_c,zeta_c,coef_c)
8*                                                                      *
9*   Do screening check for ECP integrals                               *
10*                                                                      *
11*   Argument (status) - description                                    *
12*                                                                      *
13*   l_a (inp) - (maximum) angular momentum of functions on centre A    *
14*   n_prim_a (inp) - number of primitive functions on centre A         *
15*   n_cont_a (inp) - number of contracted functions on centre A        *
16*   coef_a (inp) - centre A contraction coefficients                   *
17*   zeta_a (inp) - centre A exponents                                  *
18*   R_AC (inp) - distance between centres A and C                      *
19*   l_b (inp) - (maximum) angular momentum of functions on centre B    *
20*   n_prim_b (inp) - number of primitive functions on centre B         *
21*   n_cont_b (inp) - number of contracted functions on centre B        *
22*   coef_b (inp) - centre B contraction coefficients                   *
23*   zeta_b (inp) - centre B exponents                                  *
24*   R_BC (inp) - distance between centres B and C                      *
25*   n_coef_c (inp) - number of coefficients/exponents for ECP          *
26*   zeta_c (inp) - ECP exponents                                       *
27*   coef_c (inp) - ECP contraction coefficients                        *
28*                                                                      *
29*                                                                      *
30*   Written by K. G. Dyall                                             *
31*                                                                      *
32************************************************************************
33      implicit none
34#include "ecp_consts.fh"
35      integer l_a,n_prim_a,n_cont_a
36      integer l_b,n_prim_b,n_cont_b
37      Integer n_coef_c
38      integer i_a,i_b,i_c
39      double precision zeta_c(n_coef_c)
40      double precision coef_c(n_coef_c)
41      double precision coef_a(n_prim_a,n_cont_a)
42      double precision coef_b(n_prim_b,n_cont_b)
43      double precision zeta_a(n_prim_a)
44      double precision zeta_b(n_prim_b)
45      double precision R_AB
46      double precision R_AC
47      double precision R_BC
48      double precision log_prefac,pr,pre,pref,sum
49      double precision a,b,c,ca,cb,aA,bB,ab,aAbB,abc,P,wa,wb
50      logical first
51*
52      integer idamax
53      external idamax
54*
55      wa = l_a
56      wb = l_b
57      first = .true.
58      log_prefac = zero
59*
60      R_AB = abs(R_AC-R_BC)
61      do i_b = 1,n_prim_b
62        b = zeta_b(i_b)
63        i_c = idamax(n_cont_b,coef_b(i_b,1),n_prim_b)
64        if (abs(coef_b(i_b,i_c)) .ne. zero) then
65          cb = log(abs(coef_b(i_b,i_c)))
66        else
67          cb = zero
68        end if
69        bB = b*R_BC
70        do i_a = 1,n_prim_a
71          a = zeta_a(i_a)
72          i_c = idamax(n_cont_a,coef_a(i_a,1),n_prim_a)
73          if (abs(coef_a(i_a,i_c)) .ne. zero) then
74            ca = log(abs(coef_a(i_a,i_c)))
75          else
76            ca = zero
77          end if
78          aA = a*R_AC
79          aAbB = aA+bB
80          ab = zeta_a(i_a)+zeta_b(i_b)
81          pref = zero
82          sum = zero
83          do i_c = 1,n_coef_c
84            c = zeta_c(i_c)
85            abc = ab+c
86            P = aAbB/abc
87            pre = 0.5d0*(log(pi)-log(abc))
88     &          - (a*b*R_AB**2 + b*c*R_BC**2 + a*c*R_AC**2)/abc
89            if (abs(R_AC-P) .gt. 1.0d-10) pre = pre+wa*log(abs(R_AC-P))
90            if (abs(R_BC-P) .gt. 1.0d-10) pre = pre+wb*log(abs(R_BC-P))
91            if (pref .eq. zero) then
92              pref = pre
93              sum = coef_c(i_c)
94            else
95              pr = pre - pref
96              if (pr .lt. zero) then
97                if (pr .gt. ln_thr_ecp)
98     &              sum = sum+coef_c(i_c)*exp(pr)
99              else
100                if (-pr .gt. ln_thr_ecp) then
101                  sum = sum*exp(-pr)
102                else
103                  sum = zero
104                end if
105                sum = sum+coef_c(i_c)
106                pref = pre
107              end if
108            end if
109          end do
110          if (sum .ne. zero) pref = pref+log(abs(sum))
111          pref = pref+ca+cb
112          if (first) then
113            log_prefac = pref
114            first = .false.
115          else
116            log_prefac = max(log_prefac, pref)
117          end if
118        end do
119      end do
120*
121      log_prefac = log_prefac + log(dble(n_prim_a*n_prim_b))
122      ecp_skipint = log_prefac .lt. ln_thr_ecp
123*
124      return
125      end
126
127
128