1      subroutine dim_grid_quadv0b(ictr_buf, qxyz, qwght, nqpts, rad,
2     $                            ncoef_max, natoms, iniz, ifin, rho,
3     $                            xyz, expo, Bmat, Fmat, Pmat, rq,
4     $                            cetobfr, ibf_ao, bas_cent_info,
5     $                            bas_cset_info, docset, rchi_atom,
6     $                            g_vdim, nmat, do_scf, do_cpks_l, l3d,
7     $                            ipm, imag, dimxyz, muind)
8c
9c$Id: dim_grid_quadv0b.F 19821 2010-12-14 07:46:49Z d3y133 $
10c
11      implicit none
12#include "errquit.fh"
13#include "cgridfile.fh"
14#include "stdio.fh"
15#include "cdft.fh"
16#include "bas.fh"
17#include "mafdecls.fh"
18#include "global.fh"
19#include "util.fh"
20#include "grid_sic.fh"
21#include "dftps.fh"
22#include "geom.fh"
23#include "dimqm.fh"
24c
25c   Input Variables
26      integer ictr_buf ! ctr of grid
27      double precision qxyz(3, n_per_rec) ! Quadrature point coordinates
28      double precision qwght(n_per_rec)   ! Quadrature point weights
29      integer nqpts                       ! Number of quadrature points
30      double precision rad
31      integer ncoef_max
32      integer natoms    ! Number of atoms
33      integer iniz(natoms)
34      integer ifin(natoms)
35      double precision rho(*)
36      double precision xyz(3, natoms)   ! Atom coordinates
37      double precision expo(nbf_ao_mxprim)
38      double precision Bmat(nbf_ao_mxnbf_ce*n_per_rec)
39      double precision Fmat(*)
40      double precision Pmat(*)
41      double precision rq(n_per_rec, natoms)
42      integer cetobfr(2, natoms)
43      integer ibf_ao(nbf_ao)
44      integer bas_cent_info(3,natoms)
45      integer bas_cset_info(6,*)
46      logical docset(*)
47      double precision rchi_atom(natoms)
48      integer g_vdim ! DIM potential for this block of the quadrature
49      integer nmat ! Number of density matrices (number of perturbing directions)
50      logical do_scf
51      logical do_cpks_l
52      logical l3d
53      integer ipm
54      integer imag
55      double precision dimxyz(3, nDIM) ! DIM Coordinates
56      double precision muind(3, nDIM, nmat) ! DIM Coordinates
57c
58c   Local Variables
59      integer nbf_ao_mxnbf_ce2, mbf_ao, npol, ncontrset, maxdbas
60      integer lchi_ao, ichi_ao, ldelchi_ao, idelchi_ao
61      integer ldmat,idmat,i, dim_grid_nbfm
62      double precision dabsmax
63      double precision acc_ao_gauss
64      external dabsmax
65      external dim_grid_nbfm
66      integer iscf_rho,iscf_delrho
67      integer iscf_tau,iscf_ttau
68      logical grid_sicinit,grid_sicend
69      external grid_sicinit,grid_sicend
70      double precision acc_xc_gauss
71      integer lemat,iemat,lfmat,ifmat,k_scr,l_scr
72      logical do_2nd
73      logical stat
74      logical dprint
75c
76      npol = (ipol*(ipol+1))/2
77      acc_ao_gauss = dble(iaoacc)
78      nbf_ao_mxnbf_ce2 = nbf_ao_mxnbf_ce*nbf_ao_mxnbf_ce
79      if(nqpts.eq.0) return
80      dprint = ldebug
81c      if(dprint) then
82c        write(luout,*) "dim_grid_quadv0b start"
83c      end if
84       maxdbas = 0
85       idelchi_ao = 0
86c
87c     Evaluate the AO basis set at each of the quad. points.
88c     allocate arrays for exponents and contraction coefficients
89c     within int_eval_gbsets
90c     Evaluate the distances (squared) between each center and the points
91c
92      mbf_ao=nbf_ao
93      if(natoms.gt.1) then
94         call icopy(mbf_ao, 0,0, ibf_ao, 1)
95
96         mbf_ao=dim_grid_nbfm(AO_bas_han,  natoms,
97     &        ictr_buf,rad,xyz,qxyz,nqpts,
98     &        ibf_ao, docset, iniz, ifin, expo,
99     &        minexp,ldiff,acc_ao_gauss,iatype_pt_chg)
100         if (mbf_ao.eq.0) return
101      else
102        if (.not.bas_numcont(AO_bas_han, ncontrset))
103     &     call errquit('Exiting in grid_quadv0b',0, BASIS_ERR)
104        iniz(1)=1
105        ifin(1)=nbf_ao
106        do i=1,ncontrset
107          docset(i)=.true.
108        enddo
109        do i=1,nbf_ao
110          ibf_ao(i)=i
111        enddo
112      endif
113
114      if (.not.MA_Push_Get(mt_dbl, nqpts*mbf_ao, 'chi_ao',
115     &     lchi_ao, ichi_ao))
116     &     call errquit('grid_quadv0b: cannot allocate chi_ao',0,
117     &       MA_ERR)
118
119      call qdist(rchi_atom, rq, qxyz, xyz, nqpts, natoms)
120
121      call xc_eval_basis(ao_bas_han, maxdbas, dbl_mb(ichi_ao),
122     &     0d0, 0d0, 0d0, rq, qxyz, xyz, nqpts, natoms,
123     &     iniz, docset, bas_cent_info, bas_cset_info)
124c
125c     get reduced Xi(r) and dXi(r) over atoms
126c
127      call util_rmsatbf(nqpts, natoms,iniz,ifin,
128     ,     dbl_mb(ichi_ao),rchi_atom)
129
130      if(nqpts.eq.0) goto 2010
131
132c
133c     Evaluate the DIM potential
134c     Set up pointer to the SCF density for the CPKS LHS case.
135c
136      call dim_eval_fnl(rho, nqpts, qwght,
137     &     dimxyz, qxyz,
138     $     nmat, ipm, imag, muind)
139c       rho(1:nqpts*nmat) = 0.0d0
140c
141c   Transform the potential from over the quadrature points to the shape of the Fock Matrix
142      if(do_scf) then
143c
144        call dim_tabcd(0, l3d, Fmat, Pmat, rho, Bmat,
145     $                   dbl_mb(ichi_ao), 1, nqpts, mbf_ao,
146     $                   nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2,
147     $                   AO_bas_han, natoms, iniz, ifin, g_vdim,
148     $                   ibf_ao, cetobfr)
149c
150      else if (do_cpks_l) then
151c
152        call dim_tabcd(1, l3d, Fmat, Pmat, rho, Bmat,
153     $                   dbl_mb(ichi_ao), nmat, nqpts, mbf_ao,
154     $                   nbf_ao_mxnbf_ce, nbf_ao_mxnbf_ce2,
155     $                   AO_bas_han, natoms, iniz, ifin, g_vdim,
156     $                   ibf_ao, cetobfr)
157c
158      endif
159c         write(luout,*) "g_dim AFTER dim_tabcd"
160c         call ga_print(g_vxc)
161c
162c      endif
163
164 2010 continue
165
166      if (sic_orb_index.eq.1) then
167         if(.not.grid_sicend(l_vect1,ldelrho_sig))
168     ,        call errquit(' grid_quadv0b: grid_sicend failed',0,
169     &       CALC_ERR)
170      endif
171
172      if (.not.ma_pop_stack(lchi_ao))
173     &     call errquit('grid_quadv0b: cannot pop stack', 3, MA_ERR)
174
175c      if(dprint) then
176c        write(luout,*) "end dim_grid_quadv0b"
177c      end if
178      return
179      end subroutine dim_grid_quadv0b
180c
181c     function dim_grid_nbfm
182c
183      integer function dim_grid_nbfm(basis_hand,   mcenters,
184     &     ctr_quad,  rad_quad, xyz, qxyz,nq,
185     &                  ibf, docset, iniz, ifin, zprim,
186     .     minexp,ldiff,acc_gauss,iatype_pt_chg)
187c
188C$Id: grid_quadv0b.F 19821 2010-12-14 07:46:49Z d3y133 $
189c
190      implicit none
191#include "errquit.fh"
192c
193      integer basis_hand
194      integer mcenters ! [input]
195      double precision acc_gauss ! [input]
196      double precision xyz(3,*) ! [input]
197      double precision rad_quad ! something about radius of this integration shell [input]
198      integer nq ! [in] # grid pts
199      double precision qxyz(3,nq) ! [input] coord grid pts
200      integer ctr_quad ! grid center  [input]
201      integer ibf(*) ! [output]
202      logical docset(*) ! [output]
203      integer ldiff(*) ! [in]
204      double precision minexp(*) ! [in]
205      logical iatype_pt_chg(*) ! [in]
206c
207#include "bas.fh"
208c
209c     Distance Squared between Sampling Points and Centers
210c
211      double precision r_q0 ! min dist
212      integer  iniz(mcenters),ifin(mcenters)
213      double precision zprim(*)
214      integer n1, icset, ictr,  nprim, ncontr,
215     &        isphere,  l, iprimo,npt,nshbf
216      double precision zmin,acc_loc,x,y,z
217      integer n,ic1,ic2,m,ibfhi,ibflo
218      double precision alpha,logeps
219      double precision gaussian_range,r2,r_arg
220      logical veryfar
221c
222      gaussian_range(n,logeps,alpha) =
223     $    (n*log(-logeps) - n*log(alpha) - 4.0d0*logeps) /
224     $    sqrt(-16.0d0*alpha*logeps)
225c
226      call ifill(mcenters,0,iniz,1)
227      call ifill(mcenters,0,ifin,1)
228      if(acc_gauss.gt.25d0.or.acc_gauss.lt.3d0) call errquit(
229     ' ' grid_nbfm: silly accgauss ',nint(acc_gauss), UNKNOWN_ERR)
230      acc_loc=-acc_gauss
231c
232      n1 = 0
233      npt=0
234      do 400 ictr=1,mcenters ! Loop over all atoms
235        if(iatype_pt_chg(ictr)) goto 400 ! Cycle if this atom is a point charge
236        if (.not.bas_ce2cnr(basis_hand,ictr,ic1,ic2))
237     &      call errquit('Exiting in xc_signf.',11, BASIS_ERR)
238
239        r2=0d0
240        if (ictr.ne.ctr_quad) then ! Calculate r-squared if this atom isn't the center of the quadrature
241          x = xyz(1,ctr_quad) - xyz(1,ictr)
242          y = xyz(2,ctr_quad) - xyz(2,ictr)
243          z = xyz(3,ctr_quad) - xyz(3,ictr)
244          r2 = sqrt(x*x + y*y + z*z)
245        endif
246        r_arg=0d0
247        if (rad_quad.lt.r2) r_arg = (r2-rad_quad)
248c
249c     check on most diffuse fn
250c
251        veryfar=r_arg.gt.gaussian_range(ldiff(ictr),acc_loc,
252     .                                  minexp(ictr))
253        if(veryfar) then
254          if(.not.bas_ce2bfr(basis_hand, ictr, ibflo, ibfhi))
255     &        call errquit('Exiting in grid_nbfm',4, BASIS_ERR)
256          nshbf=ibfhi-ibflo+1
257          npt=npt+nshbf
258        else
259          r_q0=1d10
260          do n=1,nq
261            x = qxyz(1,n) - xyz(1,ictr)
262            y = qxyz(2,n) - xyz(2,ictr)
263            z = qxyz(3,n) - xyz(3,ictr)
264            r_q0 = min(r_q0,sqrt(x*x + y*y + z*z))
265          enddo
266          do icset = ic1,ic2
267            docset(icset) = .false.
268c
269c       get info about current contraction set
270            if (.not. bas_continfo(basis_hand, icset,  l ,nprim,
271     &                             ncontr, isphere))
272     &          call errquit('Exiting in grid_nbfm.',5, BASIS_ERR)
273c
274c       get exponents and contraction coefficients for this contraction set
275            if (.not.bas_get_exponent(basis_hand, icset, zprim))
276     &          call errquit('Exiting in grid_nbfm.',7, BASIS_ERR)
277c
278c       Determine the minimum Gaussian exponent.
279            zmin = 1.D+06
280            do iprimo = 1,nprim
281              zmin = min(zprim(iprimo),zmin)
282            enddo
283c
284c       Only include those basis functions that are "non-zero" for at least
285c       one  point in the sampling set.
286            nshbf=ncontr*(((l+1)*(l+2))/2)
287            if(isphere.eq.1) then
288              nshbf=ncontr*(2*l+1)
289            endif
290c
291c     pre-screening for worst case (max radius)
292c
293            if (r_q0.lt.gaussian_range(l,acc_loc,zmin)) then
294              do m=1,nshbf
295                ibf(n1+m) = npt+m
296              enddo
297              docset(icset) = .true.
298              if (iniz(ictr).eq.0) iniz(ictr)=n1+1
299              n1=n1+nshbf
300            endif
301            npt=npt+nshbf
302          enddo ! icset = ic1, ic2
303          ifin(ictr)= n1
304        endif ! If veryfar
305  400 continue ! End loop over atoms
306c
307      dim_grid_nbfm = n1
308c
309      return
310      end
311c
312c     Precompute relevant basis info for XC
313c
314c     BGJ - 9/00
315c
316      Subroutine dim_make_basis_info(basis_hand, bas_cent_info,
317     &     bas_cset_info, mcenters)
318c
319C$Id: grid_quadv0b.F 19821 2010-12-14 07:46:49Z d3y133 $
320c
321      implicit none
322#include "errquit.fh"
323c
324      integer basis_hand        ! [input]
325      integer mcenters          ! [input]
326      integer bas_cent_info(3,mcenters) ! [output]
327      integer bas_cset_info(6,*) ! [output]
328c
329#include "bas.fh"
330c
331      integer ictr, icset, ifirst, ilast, nprim, ncontr, l, isphere
332      integer ic1, ic2
333c
334      do ictr = 1,mcenters
335         bas_cent_info(1,ictr) = 0 ! max angular momentum
336         if (.not.bas_ce2cnr(basis_hand,ictr,ic1,ic2))
337     &        call errquit('Exiting in xc_make_basis_info',1, BASIS_ERR)
338         bas_cent_info(2,ictr) = ic1
339         bas_cent_info(3,ictr) = ic2
340c
341         do icset = ic1,ic2
342c
343c     Get info about current contraction set - first and last basis function,
344c     angular momentum, number of primitives, number of contractions and
345c     whether spherical harmomic
346c
347            if (.not. bas_cn2bfr(basis_hand, icset, ifirst, ilast))
348     &           call errquit('Exiting in xc_make_basis_info',1,
349     &       BASIS_ERR)
350            if (.not. bas_continfo(basis_hand, icset, l, nprim,
351     &           ncontr, isphere))
352     &           call errquit('Exiting in xc_make_basis_info',2,
353     &       BASIS_ERR)
354c
355            if (l.lt.0) then
356               call errquit('L code < 0 in xc_make_basis_info',1,
357     &       BASIS_ERR)
358            endif
359c
360            bas_cent_info(1,ictr) = max(bas_cent_info(1,ictr),l)
361c
362            bas_cset_info(1,icset) = ifirst
363            bas_cset_info(2,icset) = ilast
364            bas_cset_info(3,icset) = l
365            bas_cset_info(4,icset) = nprim
366            bas_cset_info(5,icset) = ncontr
367            bas_cset_info(6,icset) = isphere
368         enddo
369      enddo
370c
371      return
372      end
373