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