1c $Id$
2*
3C> \ingroup nwint
4C> @{
5C>
6C> \brief Compute projected multipole integrals up to the specified level
7C>
8C> This routine computes projected multipole integrals up to level
9C> \f$l_{max}\f$ (\f$0 \rightarrow l_{max}\f$).
10C> The general form is \f$\langle pole | shell \rangle\f$.
11C>
12C> Integrals are returned in \f$\langle pole | shell \rangle\f$ blocks
13C> one block for each \f$L\f$ value \f$0 \rightarrow l_{max}\f$.
14C>
15C> For example, a multipole, L = 1 and a D shell
16C> would yield (1+3)*6 = 24 integrals.
17C>
18C> The order would be:
19C> <table>
20C> <tr><td> <0|xx> </td><td> <0|xy> </td><td> <0|xz> </td><td> <0|yy> </td><td> \f$\ldots\f$ </td><td> <0|zz> </td><td> first six </td></tr>
21C> <tr><td> <x|xx> </td><td> <x|xy> </td><td> <x|xz> </td><td> <x|yy> </td><td> \f$\ldots\f$ </td><td> <x|zz> </td><td> second six </td></tr>
22C> <tr><td> <y|xx> </td><td> <y|xy> </td><td> <y|xz> </td><td> <y|yy> </td><td> \f$\ldots\f$ </td><td> <y|zz> </td><td> third six </td></tr>
23C> <tr><td> <z|xx> </td><td> <z|xy> </td><td> <z|xz> </td><td> <z|yy> </td><td> \f$\ldots\f$ </td><td> <z|zz> </td><td> third six </td></tr>
24C> </table>
25C>
26c:tex-% This is part of the standard API routines
27c:tex-\subsection{int\_projpole}
28c:tex-This routine computes projected multipole integrals up to level
29c:tex-lmax (0$\rightarrow$lmax):\\
30c:tex-The general form is $<$pole$|$shell$>$
31c:tex-
32c:tex-Integrals are returned in $<$pole$|$shell$>$ blocks one block for each
33c:tex-L value 0$\rightarrow$lmax.
34c:tex-
35c:tex-For example, a multipole, L = 1 and a d shell
36c:tex-would yield (1+3)*6 = 24 integrals.
37c:tex-
38c:tex-The order would be:\\
39c:tex-\begin{tabular}{ccccccc}
40c:tex-$<0|xx>$ & $<0|xy>$ & $<0|xz>$ & $<0|yy>$ & $\ldots$ & $<0|zz>$ & first six\\
41c:tex-$<x|xx>$ & $<x|xy>$ & $<x|xz>$ & $<x|yy>$ & $\ldots$ & $<x|zz>$ & second six\\
42c:tex-$<y|xx>$ & $<y|xy>$ & $<y|xz>$ & $<y|yy>$ & $\ldots$ & $<y|zz>$ & third six \\
43c:tex-$<z|xx>$ & $<z|xy>$ & $<z|xz>$ & $<z|yy>$ & $\ldots$ & $<z|zz>$ & fourth six
44c:tex-\end{tabular}
45c:tex-
46c:tex-{\it Syntax:}
47c:tex-\begin{verbatim}
48      subroutine int_projpole(i_basis, ish, centerl, lmax,
49     &    lscr, scr, lmpint, MP)
50c:tex-\end{verbatim}
51      implicit none
52#include "apiP.fh"
53#include "errquit.fh"
54#include "nwc_const.fh"
55#include "basP.fh"
56#include "basdeclsP.fh"
57#include "geobasmapP.fh"
58#include "geomP.fh"
59#include "mafdecls.fh"
60#include "bas_exndcf_dec.fh"
61#include "bas_ibs_dec.fh"
62#include "stdio.fh"
63#include "int_nbf.fh"
64c
65c::functions
66      logical int_chk_init
67      logical cando_nw_1e, cando_nw
68      integer int_nint_cart, int_nint
69      external int_chk_init
70      external cando_nw_1e, cando_nw
71      external int_nint_cart, int_nint
72c::passed
73c:tex-\begin{verbatim}
74      integer i_basis             !< [Input] basis set handle for ish
75      integer ish                 !< [Input] i shell/contraction
76      integer lmax                !< [Input] maximum lvalue for
77                                  !< multipole ints in this batch
78      double precision centerl(3) !< [Input] coordinates of multipole
79      integer lscr                !< [Input] length of scratch array
80      double precision scr(lscr)  !< [Input] scratch array
81      integer lmpint              !< [Input/Output] length of
82                                  !< multipole integrals array
83      double precision MP(lmpint) !< [Output] multipole integrals
84c:tex-\end{verbatim}
85c::local
86      logical shells_ok
87      integer ibas, Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom
88      integer i_gen_save
89      integer ucont, lval
90      integer l_int, i_int, i_int_x,num_int
91      integer int_point, int_size
92      integer lpole
93      logical inline_chk_sh
94c
95#include "bas_exndcf_sfn.fh"
96#include "bas_ibs_sfn.fh"
97c... statement function for int_chk_sh
98      inline_chk_sh(ibas,ish) =
99     $     ((ish.gt.0) .and. (ish.le.ncont_tot_gb(ibas)))
100c
101c check initialization
102c
103      if (.not.int_chk_init('int_projpole'))
104     &       call errquit('int_projpole: int_init was not called' ,0,
105     &          INT_ERR)
106c
107      ibas = i_basis + BASIS_HANDLE_OFFSET
108c
109      shells_ok = inline_chk_sh(ibas,ish)
110      if (.not. shells_ok)
111     &       call errquit('int_projpole: invalid contraction/shell',0,
112     &            BASIS_ERR)
113c
114***   set defNxyz such that it can handle the maximum multi-pole
115c
116      lpole = lmax/4 + 1
117      call defNxyz(lpole)
118c
119      ucont   = (sf_ibs_cn2ucn(ish,ibas))
120      Li      = infbs_cont(CONT_TYPE ,ucont,ibas)
121      i_prim  = infbs_cont(CONT_NPRIM,ucont,ibas)
122      i_gen   = infbs_cont(CONT_NGEN ,ucont,ibas)
123      i_iexp  = infbs_cont(CONT_IEXP ,ucont,ibas)
124      i_icfp  = infbs_cont(CONT_ICFP ,ucont,ibas)
125      i_cent  = (sf_ibs_cn2ce(ish,ibas))
126      i_geom  = ibs_geom(ibas)
127c
128      num_int = 0
129      do 00100 lval = 0,lmax
130        l_int   = (lval+1)*(lval+2)/2
131        i_int_x = int_nint_cart(i_basis, ish, 0, 0, 0,0, 0,0)
132        i_int   = int_nint(i_basis, ish, 0, 0, 0,0, 0,0)
133        int_point = num_int + 1
134        int_size  = lmpint - num_int
135        if ((num_int+l_int*i_int_x).gt.lmpint) then
136          write(luout,*)
137     &        ' int_projpole: lmpint  = ',lmpint
138          write(luout,*)
139     &        ' int_projpole: (num_int+l_int*i_int_x) = ',
140     &        (num_int+l_int*i_int_x)
141          call errquit('int_projpole: lmpint too small ',911, INT_ERR)
142        endif
143        num_int = num_int + l_int*i_int
144        if (cando_nw_1e(i_basis,ish,0)) then
145          call int_hf1sp(
146     &        centerl,DCexp,
147*                                   ! no ECPs allowed when projpole?
148     &        DCcoeff, 1, 1, lval, -1,
149     &        coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)),
150     &        dbl_mb(mb_exndcf(i_icfp,ibas)),
151     &        i_prim, i_gen, Li, i_cent,
152     &        coords(1,1,i_geom),charge(1,i_geom),
153     &        geom_invnucexp(1,i_geom),ncenter(i_geom),
154     &        MP(int_point),scr,scr,int_size,
155c............. doS    doT     doV     canonical
156     &        .true.,.false.,.false.,.false.,
157c.............. dryrun
158     &        .false.,scr,lscr,'int_projpole')
159          if (bas_spherical(ibas)) then    ! only one basis set
160            if (Li.eq.-1) then
161              i_gen_save = i_gen
162              i_gen = 1
163            endif
164            call spcart_2ctran(MP(int_point),scr,lscr,
165     &          int_nbf_x(lval),int_nbf_x(lval),lval,1,.false.,
166     &          int_nbf_x(Li),int_nbf_s(Li),Li,i_gen,
167     &          bas_spherical(ibas),.false.)
168            if (Li.eq.-1) i_gen = i_gen_save
169          endif
170        else
171          call errquit
172     &        ('int_projpole: could not do sp or nw integrals',0,
173     &           INT_ERR)
174        endif
175c
17600100 continue
177      end
178C> @}
179