1      Subroutine grid_gen_pts(rtdb)
2c
3c$Id$
4c
5      implicit none
6#include "errquit.fh"
7c
8#include "inp.fh"
9#include "rtdb.fh"
10#include "mafdecls.fh"
11#include "cdft.fh"
12#include "geom.fh"
13#include "cgridfile.fh"
14c
15      integer rtdb
16      integer lcoord, icoord, lcharge, icharge, ltags, itags
17      integer mxnrad, nqlen, lqs, iqs
18      integer ang_quad_buf_size, nq_tot
19      logical grid_file_balance,grid_getnsubb,oldgrid
20      external  grid_file_balance,grid_getnsubb
21c
22c     Need coordinates of all atoms.
23c
24c     allocate space for atomic coordinates and charges
25c
26      if (.not. Ma_Push_Get(MT_Dbl,ncenters*3,'coordinates',lcoord,
27     &   icoord))call errquit(
28     (     'grid_gen_pts: failed to alloc coordinates',0, MA_ERR)
29      if (.not. Ma_Push_Get(MT_Dbl,ncenters,'charges',lcharge,
30     &   icharge))call errquit(
31     (     'grid_gen_pts: failed to alloc charges',0, MA_ERR)
32      if (.not. Ma_Push_Get(MT_Byte, ncenters*16, 'center tags',
33     &   ltags, itags))call errquit(
34     (     'grid_gen_pts: failed to alloc center tags',0, MA_ERR)
35c
36c     Get ncenter tags, coordinates, and charges from the geometry object.
37c
38      if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
39     &                        Dbl_MB(icoord), Dbl_MB(icharge)))
40     &     call errquit('grid_gen_pts: geom_cart_get failed',74,
41     &       GEOM_ERR)
42      oldgrid=.true.
43c
44c     A list of pruned (or unpruned) r, theta and phi is generated here.
45c
46      if (.not. rtdb_get(rtdb, 'dft:mxnrad', mt_int, 1, mxnrad))
47     &   call errquit('grid_gen: rtdb_get failed', 119, RTDB_ERR)
48c
49      if(oldgrid) then
50      nqlen = ncenters*mxnrad
51      if (.not.MA_Push_get(MT_int,4*nqlen,'shell list',
52     &                     lqs,iqs))
53     &   call errquit('grid_gen: cannot allocate shell list',0, MA_ERR)
54      call grid_list(rtdb, int_mb(iqs), nqlen, nq_tot)
55      endif
56c
57c     The points (rotationally invariant if desired) are generated and written to disk
58c     with fixed number of points (buffer size).
59c
60c     buffer size needed to hold largest angular quadrature is nqmax
61c
62      ang_quad_buf_size = nqmax
63c
64      if (.not. MA_Pop_Stack(lqs))
65     &   call errquit('grid_gen_pts: pop stack failed.',0, MA_ERR)
66      if (.not. MA_Pop_Stack(ltags))
67     &   call errquit('grid_gen_pts: pop stack failed.',0, MA_ERR)
68      if (.not. MA_Pop_Stack(lcharge))
69     &   call errquit('grid_gen_pts: pop stack failed.',0, MA_ERR)
70      if (.not. MA_Pop_Stack(lcoord))
71     &   call errquit('grid_gen_pts: pop stack failed.',0, MA_ERR)
72      return
73      end
74