1C> \ingroup nwint
2C> @{
3C>
4C> \brief Compute 1-electron Gaussian periodic multipole
5C>  integral derivatives for \f$l\f$ from \f$0\f$ to \f$l_{max}\f$
6C>
7C> See [1] for details.
8C>
9C> [1] JE Jaffe, AC Hess,
10C>     <i>"Gaussian basis density functional theory for systems
11C>     periodic in two or three dimensions: Energy and forces"</i>,
12C>    J.Chem.Phys. <b>105</b>, 10983-10998 (1996), DOI:
13C>    <a href="https://doi.org/10.1063/1.472866">
14C>    10.1063/1.472866</a>
15C>
16      subroutine intpd_mpolel(i_basis, ish, j_basis, jsh, R,
17     &    lval, centerl,
18     &    lscr, scr, lmpint, MP, num_mpint,
19     &    idatom)
20*
21* $Id$
22*
23c
24c routine to compute multipole integral derivatives at a given lvalue
25c with the jsh translated by the fractional coordinate vector R.
26c The general form is <shell|pole|shell|3|3>
27c
28c the returned buffer is logically:
29c           (mpole range, jlo:jhi, ilo:ihi, 3(xyz), 3(atom))
30c where mpole range is 1:((lval+1)*(lval+2)/2)
31c       3(xyz) is the x,y,z coordinate derivative for the atom index
32c       3(atom) is at most three centers one is where ish exists
33c                                        two is where jsh exists
34c                                        three is the center of the multipole
35c
36c  Integral derivatives are returned in shell blocks of <L|ish|jsh> L=lval
37c  one block for the given L value.  EACH block repeated for by 9 (xyz,atoms)
38c  for ish = d and Lval = 1 and jsh = p you would get:
39c      (6*3*3)*3*3=486 integral derivatives
40c  order would be
41c   <xx|x|x> <xx|y|x> <xx|z|x> ( 1- 3)(x,atom1)
42c   <xx|x|y> <xx|y|y> <xx|z|y> ( 4- 6)(x,atom1)
43c   <xx|x|z> <xx|y|z> <xx|z|z> ( 7- 9)(x,atom1)
44c   <xy|x|x> <xy|y|x> <xy|z|x> (10-12)(x,atom1)
45c   <xy|x|y> <xy|y|y> <xy|z|y> (13-15)(x,atom1)
46c   <xy|x|z> <xy|y|z> <xy|z|z> (16-18)(x,atom1)
47c   <xz|x|x> <xz|y|x> <xz|z|x> (19-21)(x,atom1)
48c   <xz|x|y> <xz|y|y> <xz|z|y> (22-24)(x,atom1)
49c   <xz|x|z> <xz|y|z> <xz|z|z> (25-27)(x,atom1)
50c   <yy|x|x> <yy|y|x> <yy|z|x> (28-30)(x,atom1)
51c   <yy|x|y> <yy|y|y> <yy|z|y> (31-33)(x,atom1)
52c   <yy|x|z> <yy|y|z> <yy|z|z> (34-36)(x,atom1)
53c   <yz|x|x> <yz|y|x> <yz|z|x> (37-39)(x,atom1)
54c   <yz|x|y> <yz|y|y> <yz|z|y> (40-42)(x,atom1)
55c   <yz|x|z> <yz|y|z> <yz|z|z> (43-45)(x,atom1)
56c   <zz|x|x> <zz|y|x> <zz|z|x> (46-48)(x,atom1)
57c   <zz|x|y> <zz|y|y> <zz|z|y> (49-51)(x,atom1)
58c   <zz|x|z> <zz|y|z> <zz|z|z> (52-54)(x,atom1)
59c   repeat above for (y,atom1), (z,atom1)
60c   repeat above for atom2 and multipole center
61c
62c
63c
64c  for ish = p and Lval = 2 and jsh = p you would get:
65c      (3*6*3)*3*3 = 486 integral derivatives
66c  order would be
67c   <x|xx|x> <x|xy|x> <x|xz|x> <x|yy|x> <x|yz|x> <x|zz|x>  ( 1- 6)(x,atom1)
68c   <x|xx|y> <x|xy|y> <x|xz|y> <x|yy|y> <x|yz|y> <x|zz|y>  ( 7-12)(x,atom1)
69c   <x|xx|z> <x|xy|z> <x|xz|z> <x|yy|z> <x|yz|z> <x|zz|z>  (13-18)(x,atom1)
70c   <y|xx|x> <y|xy|x> <y|xz|x> <y|yy|x> <y|yz|x> <y|zz|x>  (19-24)(x,atom1)
71c   <y|xx|y> <y|xy|y> <y|xz|y> <y|yy|y> <y|yz|y> <y|zz|y>  (25-30)(x,atom1)
72c   <y|xx|z> <y|xy|z> <y|xz|z> <y|yy|z> <y|yz|z> <y|zz|z>  (31-36)(x,atom1)
73c   <z|xx|x> <z|xy|x> <z|xz|x> <z|yy|x> <z|yz|x> <z|zz|x>  (37-42)(x,atom1)
74c   <z|xx|y> <z|xy|y> <z|xz|y> <z|yy|y> <z|yz|y> <z|zz|y>  (43-48)(x,atom1)
75c   <z|xx|z> <z|xy|z> <z|xz|z> <z|yy|z> <z|yz|z> <z|zz|z>  (49-54)(x,atom1)
76c   repeat above for (y,atom1), (z,atom1)
77c   repeat above for atom2 and multipole center
78c
79c  for ish = s and lval = 4 and jsh = p you would get:
80c     (1*15*3)*3*3 = 405 integral derivatives
81c   <s|xxxx|x> <s|xxxy|x> <s|xxxz|x> <s|xxyy|x> <s|xxyz|x> <s|xxzz|x> ( 1- 6)(x,atom1)
82c   <s|xyyy|x> <s|xyyz|x> <s|xyzz|x> <s|xzzz|x> <s|yyyy|x> <s|yyyz|x> ( 7-12)(x,atom1)
83c   <s|yyzz|x> <s|yzzz|x> <s|zzzz|x>                                  (13-15)(x,atom1)
84c   <s|xxxx|y> <s|xxxy|y> <s|xxxz|y> <s|xxyy|y> <s|xxyz|y> <s|xxzz|y> (16-21)(x,atom1)
85c   <s|xyyy|y> <s|xyyz|y> <s|xyzz|y> <s|xzzz|y> <s|yyyy|y> <s|yyyz|y> (22-27)(x,atom1)
86c   <s|yyzz|y> <s|yzzz|y> <s|zzzz|y>                                  (28-30)(x,atom1)
87c   <s|xxxx|z> <s|xxxy|z> <s|xxxz|z> <s|xxyy|z> <s|xxyz|z> <s|xxzz|z> (31-36)(x,atom1)
88c   <s|xyyy|z> <s|xyyz|z> <s|xyzz|z> <s|xzzz|z> <s|yyyy|z> <s|yyyz|z> (37-42)(x,atom1)
89c   <s|yyzz|z> <s|yzzz|z> <s|zzzz|z>                                  (43-45)(x,atom1)
90c   repeat above for (y,atom1), (z,atom1)
91c   repeat above for atom2 and multipole center
92c
93c
94      implicit none
95#include "apiP.fh"
96#include "errquit.fh"
97#include "nwc_const.fh"
98#include "basP.fh"
99#include "basdeclsP.fh"
100#include "geobasmapP.fh"
101#include "geomP.fh"
102#include "stdio.fh"
103#include "mafdecls.fh"
104#include "bas_exndcf_dec.fh"
105#include "bas_ibs_dec.fh"
106#include "int_nbf.fh"
107c
108c::functions
109      logical int_chk_init
110      integer int_nint_cart
111      external int_chk_init
112      external int_nint_cart
113c::passed
114      integer i_basis             !< [Input] basis set handle for ish
115      integer ish                 !< [Input] i shell/contraction
116      integer j_basis             !< [Input] basis set handle for jsh
117      integer jsh                 !< [Input] j shell/contraction
118c...     translation vectors are in fractional coordinates
119      double precision R(3)       !< [Input] translation vec on j cont.
120      integer lval                !< [Input] maximum lvalue for
121                                  !< multipole integrals in this batch
122      double precision centerl(3) !< [Input] coordinates of multipole
123      integer lscr                !< [Input] length of scratch array
124      double precision scr(lscr)  !< [Input] scratch array
125      integer lmpint              !< [Input] length of multipole
126                                  !< integrals array
127      double precision MP(lmpint) !< [Output] multipole integrals
128      integer num_mpint           !< [Output] number of multipole integrals
129      integer idatom(3) !< [Output] array identifying centers for derivatives
130c                       ! e.g., the first nint*3  derivatives go to center idatom(1)
131c                       !       the second nint*3 derivatives go to center idatom(2)
132c::local
133      logical shells_ok
134      integer ibas, Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom
135      integer jbas, Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom
136      integer ucont
137      integer l_int, ij_int, num_int, num_intd
138      integer lpole
139      integer ioffj
140      double precision xyz_new(3)
141      logical any_spherical
142      logical inline_chk_sh
143c
144      integer WarnP
145      save WarnP
146      data WarnP /0/
147c
148#include "bas_exndcf_sfn.fh"
149#include "bas_ibs_sfn.fh"
150c
151c... statement function for int_chk_sh
152      inline_chk_sh(ibas,ish) =
153     $     ((ish.gt.0) .and. (ish.le.ncont_tot_gb(ibas)))
154c
155c check initialization
156c
157      if (.not.int_chk_init('intpd_mpolel'))
158     &       call errquit('intpd_mpolel: int_init was not called' ,0,
159     &          INT_ERR)
160c
161c  check if gencon/sp shells
162c
163      call int_nogencont_check(i_basis,'intpd_mpolel:i_basis')
164      call int_nogencont_check(j_basis,'intpd_mpolel:j_basis')
165      call int_nospshell_check(i_basis,'intpd_mpolel:i_basis')
166      call int_nospshell_check(j_basis,'intpd_mpolel:j_basis')
167c
168      ibas = i_basis + BASIS_HANDLE_OFFSET
169      jbas = j_basis + BASIS_HANDLE_OFFSET
170c
171      shells_ok = inline_chk_sh(ibas,ish)
172      shells_ok = shells_ok .and. inline_chk_sh(jbas,jsh)
173      if (.not. shells_ok)
174     &       call errquit('intpd_mpolel: invalid contraction/shell',0,
175     &           BASIS_ERR)
176c
177***   set defNxyz such that it can handle the maximum multi-pole
178c
179      lpole = lval/4 + 1
180      lpole = lpole + 1  ! for derivative
181      call defNxyz(lpole)
182c
183      ucont   = (sf_ibs_cn2ucn(ish,ibas))
184      Li      = infbs_cont(CONT_TYPE ,ucont,ibas)
185      i_prim  = infbs_cont(CONT_NPRIM,ucont,ibas)
186      i_gen   = infbs_cont(CONT_NGEN ,ucont,ibas)
187      i_iexp  = infbs_cont(CONT_IEXP ,ucont,ibas)
188      i_icfp  = infbs_cont(CONT_ICFP ,ucont,ibas)
189      i_cent  = (sf_ibs_cn2ce(ish,ibas))
190      i_geom  = ibs_geom(ibas)
191c
192      ucont   = (sf_ibs_cn2ucn(jsh,jbas))
193      Lj      = infbs_cont(CONT_TYPE ,ucont,jbas)
194      j_prim  = infbs_cont(CONT_NPRIM,ucont,jbas)
195      j_gen   = infbs_cont(CONT_NGEN ,ucont,jbas)
196      j_iexp  = infbs_cont(CONT_IEXP ,ucont,jbas)
197      j_icfp  = infbs_cont(CONT_ICFP ,ucont,jbas)
198      j_cent  = (sf_ibs_cn2ce(jsh,jbas))
199      j_geom  = ibs_geom(jbas)
200c
201      any_spherical = bas_spherical(ibas).or.bas_spherical(jbas)
202c
203      if (i_geom.ne.j_geom.and.WarnP.eq.0) then
204        write(luout,*)'intpd_mpolel: WARNING: possible geometry',
205     &      ' inconsistency'
206        write(luout,*)'i_basis geometry handle:',i_geom
207        write(luout,*)'j_basis geometry handle:',j_geom
208        WarnP = 1
209      endif
210c
211      call intp_txyz(j_cent,j_geom,R,xyz_new)
212c
213      if (i_gen.gt.1 .or. j_gen.gt.1) then
214        write(luout,*)
215     &      ' hf3ois does not handle general contractions yet'
216        call errquit('intpd_mpolel: general contraction error ',911,
217     &            BASIS_ERR)
218      endif
219c
220      l_int    = (lval+1)*(lval+2)/2
221      ij_int   = int_nint_cart(i_basis, ish, j_basis, jsh, 0,0, 0,0)
222      num_int  =  l_int*ij_int
223      num_intd = num_int*3*3
224      if (num_intd.gt.lmpint) then
225        write(luout,*)' intpd_mpolel: lmpint   = ',lmpint
226        write(luout,*)' intpd_mpolel: num_intd = ',num_intd
227        call errquit('intpd_mpolel: lmpint too small ',911,
228     &          BASIS_ERR)
229      endif
230      if ((lval.eq.0).and.(i_cent.eq.j_cent)) then
231        call dcopy(num_intd,0.0d00,0,MP,1)
232        call ifill(3,0,idatom,1)
233        return
234      endif
235c
236      call hfd3ois(
237     &    coords(1,i_cent,i_geom),
238     &    dbl_mb(mb_exndcf(i_iexp,ibas)),
239     &    dbl_mb(mb_exndcf(i_icfp,ibas)),
240     &    i_prim, 1, Li,
241     &    xyz_new,
242     &    dbl_mb(mb_exndcf(j_iexp,jbas)),
243     &    dbl_mb(mb_exndcf(j_icfp,jbas)),
244     &    j_prim, 1, Lj,
245     &    centerl,
246     &    0.0d00,
247     &    1.0d00,
248     &    1, 1, lval,
249c.....................DryRun
250     &    MP,num_int,.false.,scr,lscr)
251      num_mpint = num_intd
252      if (i_cent.eq.j_cent) then
253        idatom(1) = i_cent
254        idatom(2) = 0
255        idatom(3) = -3
256        ioffj = num_int*3 + 1
257        call daxpy(num_int*3,1.0d00,MP(ioffj),1,MP,1)
258      else
259        idatom(1) = i_cent
260        idatom(2) = j_cent
261        idatom(3) = -3
262      endif
263      if (any_spherical) then
264        if (Li.eq.-1) i_gen = 1
265        if (Lj.eq.-1) j_gen = 1
266        call spcart_3cBtran(MP,scr,lscr,
267     &      int_nbf_x(Li),int_nbf_s(Li),Li,i_gen,bas_spherical(ibas),
268     &      int_nbf_x(Lj),int_nbf_s(Lj),Lj,j_gen,bas_spherical(jbas),
269     &      int_nbf_x(lval),int_nbf_x(lval),lval,1,.false.,
270     &      9,.false.)
271      endif
272      end
273C> @}
274