1      Subroutine nbf_to_mbf(basis_hand,  mbf, mcenters,
2     &     ctr_quad, rad_quad, xyz,
3     &                  ibf,iniz, ifin, zprim,  acc_gauss)
4c
5C$Id$
6c
7      implicit none
8#include "errquit.fh"
9c
10      integer basis_hand
11      integer mcenters ! [input]
12      double precision acc_gauss ! [input]
13      double precision xyz(3,*) ! [input]
14      integer ctr_quad ! grid center  [input]
15      integer mbf  ! [input/output]
16      integer ibf(*) ! [output]
17      double precision rad_quad ! [input]
18c
19#include "bas.fh"
20c
21c     Distance Squared between Sampling Points and Centers
22c
23      integer  iniz(mcenters),ifin(mcenters)
24      double precision zprim(*)
25      integer ncontrset, n1, icset, ictr,  nprim, ncontr,
26     &        isphere,  l, iprimo,npt,nshbf
27      double precision zmin,acc_loc,r2,x,y,z
28      integer n,ic1,ic2,m
29      double precision alpha,logeps
30      double precision gaussian_range,r_arg
31      gaussian_range(n,logeps,alpha) =
32     $     (n*log(-logeps) - n*log(alpha) - 4.0d0*logeps) /
33     $     sqrt(-16.0d0*alpha*logeps)
34c
35      if (.not.bas_numcont(basis_hand, ncontrset))
36     &   call errquit('Exiting in nbf_to_mbf.',1, BASIS_ERR)
37c
38      call ifill(mcenters,0,iniz,1)
39      call ifill(mcenters,0,ifin,1)
40      acc_loc=-acc_gauss
41c
42      n1 = 0
43      npt=0
44      do 400 ictr=1,mcenters
45        if (.not.bas_ce2cnr(basis_hand,ictr,ic1,ic2))
46     &       call errquit('Exiting in xc_signf.',11, BASIS_ERR)
47        r2=0d0
48        if (ictr.ne.ctr_quad) then
49          x = xyz(1,ctr_quad) - xyz(1,ictr)
50          y = xyz(2,ctr_quad) - xyz(2,ictr)
51          z = xyz(3,ctr_quad) - xyz(3,ictr)
52          r2 = sqrt(x*x + y*y + z*z)
53        endif
54        r_arg=0d0
55        if (rad_quad.lt.r2) r_arg = (r2-rad_quad)
56        do icset = ic1,ic2
57c
58c       get info about current contraction set
59c
60          if (.not. bas_continfo(basis_hand, icset,  l ,nprim,
61     &         ncontr, isphere))
62     &         call errquit('Exiting in xc_signf.',4, BASIS_ERR)
63c
64c
65c       get exponents and contraction coefficients for this contraction set
66c
67        if (.not.bas_get_exponent(basis_hand, icset, zprim))
68     &       call errquit('Exiting in xc_signf.',7, BASIS_ERR)
69c
70c       Determine the minimum Gaussian exponent.
71c
72
73        zmin = 1.D+06
74        do iprimo = 1,nprim
75          zmin = min(zprim(iprimo),zmin)
76        enddo
77c
78c       Only include those basis functions that are "non-zero" for at least
79c       one  point in the sampling set.
80c
81          if(isphere.eq.1) then
82            nshbf=ncontr*(2*l+1)
83          else
84            nshbf=ncontr*(((l+1)*(l+2))/2)
85          endif
86        if (r_arg.lt.gaussian_range(l,acc_loc,zmin)) then
87          do m=1,nshbf
88            ibf(n1+m ) = npt+m
89          enddo
90          if (iniz(ictr).eq.0) iniz(ictr)=n1+1
91          n1=n1+nshbf
92        endif
93
94        npt=npt+nshbf
95        enddo
96c
97        ifin(ictr)= n1
98c
99  400 continue
100c
101      mbf = n1
102c
103      return
104      end
105