1      subroutine xc_signf(basis, nbf, iandex,mcenters,
2     &                      xyz,xyzm,ctr_quad, rad_quad,
3     &                       zprim,  acc_gauss,iatype_pt_chg)
4c
5C$Id$
6c
7      implicit none
8#include "errquit.fh"
9c
10      integer basis
11      integer mcenters ! [input]
12      integer iandex(*) ! [output]
13      double precision xyz(3,*) ! [input]
14      double precision xyzm(3,*) ! [input/output]
15      double precision acc_gauss ! [input]
16      double precision rad_quad ! [input]
17      logical iatype_pt_chg(*) ! [input]
18      integer nbf  ! [input/output]
19
20c
21#include "bas.fh"
22c
23c     Distance Squared between Sampling Points and Centers
24c
25      double precision zprim(*),acc_sigf
26      integer ctr_quad ! atom where the grid is centered [input]
27      integer icset, ictr,
28     &        ifirst, ilast, nprim,  iprimo
29      double precision zmin,r2,x,y,z,r_arg
30      integer mcenters_scr
31      integer ic1,ic2,dum
32      integer n,l
33      double precision alpha,logeps
34      double precision gaussian_range
35      gaussian_range(n,logeps,alpha) =
36     $     (n*log(-logeps) - n*log(alpha) - 4.0d0*logeps) /
37     $     sqrt(-16.0d0*alpha*logeps)
38
39c
40c      acc_sigf=max(acc_gauss,20d0)
41      acc_sigf=log(1d-10)
42      mcenters_scr = 0
43c
44      do ictr=1,mcenters
45
46        if(iatype_pt_chg(ictr)) goto 2001
47        if (ictr.eq.ctr_quad) then
48          mcenters_scr=mcenters_scr+1
49          iandex(mcenters_scr)=ictr
50          goto 2001
51        endif
52        if (.not.bas_ce2cnr(basis,ictr,ic1,ic2))
53     &       call errquit('Exiting in xc_signf.',11, BASIS_ERR)
54c
55c     no bf on center ie but we keep grid
56c
57        if(ic1.eq.0) then
58           mcenters_scr=mcenters_scr+1
59           iandex(mcenters_scr)=ictr
60           goto 2001
61        endif
62
63        x = xyz(1,ctr_quad) - xyz(1,ictr)
64        y = xyz(2,ctr_quad) - xyz(2,ictr)
65        z = xyz(3,ctr_quad) - xyz(3,ictr)
66        r2 = sqrt(x*x + y*y + z*z)
67        r_arg=0d0
68c        if (rad_quad.lt.r2) r_arg = (r2-rad_quad)**2
69        if (rad_quad.lt.r2) r_arg = (r2-rad_quad)
70        do icset = ic1,ic2
71c
72c       get info about current contraction set
73c
74        if (.not. bas_cn2bfr(basis, icset, ifirst, ilast))
75     &       call errquit('Exiting in xc_signf.',4, BASIS_ERR)
76            if (.not. bas_continfo(basis, icset,  l ,nprim,
77     &                             dum, dum))
78     &       call errquit('Exiting in xc_signf.',4, BASIS_ERR)
79c
80c
81c       get exponents and contraction coefficients for this contraction set
82c
83        if (.not.bas_get_exponent(basis, icset, zprim))
84     &       call errquit('Exiting in xc_signf.',7, BASIS_ERR)
85c
86c       Determine the minimum Gaussian exponent.
87c
88
89        zmin = 1.D+06
90        do iprimo = 1,nprim
91          zmin = min(zprim(iprimo),zmin)
92        enddo
93c
94c       Only include those basis functions that are "non-zero" for at least
95c       one  point in the sampling set.
96c
97c        if (zmin*r_arg.lt.acc_sigf) then
98        if (r_arg.lt.gaussian_range(l,acc_sigf,zmin)) then
99          mcenters_scr=mcenters_scr+1
100          iandex(mcenters_scr)=ictr
101          goto 2001
102        endif
103      enddo
104 2001 continue
105      enddo
106      do ictr = 1, mcenters_scr
107        ic1=iandex(ictr)
108        xyzm(1,ictr) = xyzm(1,ic1)
109        xyzm(2,ictr) = xyzm(2,ic1)
110        xyzm(3,ictr) = xyzm(3,ic1)
111      enddo
112      mcenters=mcenters_scr
113      return
114      end
115