1c
2c $Id$
3c
4c  Modified from the HONDO dipole integral code
5c
6      subroutine hnd_dipder(xyzi, expi, coefi,i_nprim,i_ngen, Li,
7     *  xyzj, expj, coefj, j_nprim, j_ngen, Lj,
8     *  center, buf, lbuf, nint, scr, lscr)
9c
10c This routine calculates the dipole derivative integrals.  This is a
11c wrapper to the work routine hnd_ddipint.
12c
13      implicit none
14c
15      integer i_nprim  ! [input] num. prims on function i
16      integer i_ngen   ! [input] num general conts on func. i
17      integer Li       ! [input] angular momentum of func. i
18      integer j_nprim  ! [input] num. prims on function j
19      integer j_ngen   ! [input] num general conts on func. j
20      integer Lj       ! [input] angular momentum of func. j
21      integer lscr     ! [input] size of scratch array
22      integer lbuf     ! [input] size of any integral
23      integer nint     ! [input] size of any integral buffer
24      double precision xyzi(3)  ! [input] position of center i
25      double precision expi(i_nprim) ! [input] exponents on i
26      double precision coefi(i_nprim,i_ngen) ! [input] i coeffs
27      double precision xyzj(3)  ! [input] position of center j
28      double precision expj(j_nprim)  ! [input] exponents on j
29      double precision coefj(j_nprim,j_ngen)  ! [input] j coeffs
30      double precision center(3)  ! [input] center for expansion
31      double precision scr(lscr)  ! [scratch] scratch buffer
32      double precision buf(lbuf)  ! [output] dipole derivative integrals
33c
34      double precision zero
35      data zero   /0.0d+00/
36c
37# include "hnd_pointers.fh"
38c
39c Zero out the integral
40c
41      call dcopy(nint*18,zero,0,buf,1)
42c
43c Use scratch for temps in integral calculation
44c Scratch needs are 6*(Li+2)*(Lj+2)+12*(Li+1)*(Lj+1)
45c
46c The pointers are defined in hdn_pointers.fh
47c
48      call hnd_ddipint(xyzi, expi, coefi,i_nprim,i_ngen, Li,
49     *  xyzj, expj, coefj, j_nprim, j_ngen, Lj,
50     *  center, buf, buf((nint*9)+1),nint,
51     *  scr(tvp(1)),scr(tvp(2)),scr(tvp(3)),scr(tvp(4)),scr(tvp(5)),
52     *  scr(tvp(6)),scr(tvp(7)),scr(tvp(8)),scr(tvp(9)),
53     *  scr(tvp(10)),scr(tvp(11)),scr(tvp(12)),scr(tvp(13)),
54     *  scr(tvp(14)),scr(tvp(15)),scr(tvp(16)),scr(tvp(17)),
55     *  scr(tvp(18)))
56c
57      return
58      end
59c
60      subroutine hnd_ddipint(xyzi, expi, coefi,i_nprim,i_ngen, Li,
61     *  xyzj, expj, coefj, j_nprim, j_ngen, Lj, center, didij, djdij,
62     *  nint, xin, yin,zin, xxin, yyin, zzin, dxsdi, dysdi, dzsdi,
63     *  dxxdi, dyydi, dzzdi, dxsdj, dysdj, dzsdj, dxxdj, dyydj, dzzdj)
64c
65c This is the main work routine for the dipole derivative integrals.
66c It is assumed the didij and djdij have been zeroed before entering.
67c
68      implicit none
69#include "hnd_tol.fh"
70#include "stdio.fh"
71#include "errquit.fh"
72      integer i_nprim  ! [input] num. prims on function i
73      integer i_ngen   ! [input] num general conts on func. i
74      integer Li       ! [input] angular momentum of func. i
75      integer j_nprim  ! [input] num. prims on function j
76      integer j_ngen   ! [input] num general conts on func. j
77      integer Lj       ! [input] angular momentum of func. j
78      integer nint     ! [input] number of base integrals
79      double precision xyzi(3)  ! [input] position of center i
80      double precision expi(i_nprim) ! [input] exponents on i
81      double precision coefi(i_nprim,i_ngen) ! [input] i coeffs
82      double precision xyzj(3)  ! [input] position of center j
83      double precision expj(j_nprim)  ! [input] exponents on j
84      double precision coefj(j_nprim,j_ngen)  ! [input] j coeffs
85      double precision center(3)  ! [input] center for expansion
86      double precision didij(nint,9) ! [output] dip.der. integrals wrt i
87      double precision djdij(nint,9) ! [ouptut] dip.der. integrals wrt j
88c
89      integer ni, nj
90      integer ijx, ijy, ijz, ix, iy, iz, jx, jy, jz
91      integer nder, lit, ljt, litder, ljtder
92      integer maxi, maxj
93      integer ig, jg, ij, i, j
94      double precision zero, one
95      double precision rln10, tol
96      double precision rr, ai, arri, axi, ayi, azi
97      double precision aj, aa, aa1, dum, fac
98      double precision csi, cpi, cdi, cfi, cgi
99      double precision csj, cpj, cdj, cfj, cgj
100      double precision ax, ay, az, dum1, dum2
101      double precision dij
102      double precision xin, yin,zin, xxin, yyin, zzin
103      double precision dxsdi, dysdi, dzsdi, dxxdi, dyydi, dzzdi
104      double precision dxsdj, dysdj, dzsdj, dxxdj, dyydj, dzzdj
105      double precision dumxdx, dumydx, dumzdx
106      double precision dumxdy, dumydy, dumzdy
107      double precision dumxdz, dumydz, dumzdz
108      double precision xint, yint, zint, xintx, yinty, zintz
109      double precision t, xc, yc, zc, x0, y0, z0
110      double precision xi, yi, zi, xj, yj, zj
111      double precision t1
112c
113      common/hnd_xyzdip/xint,yint,zint,xintx,yinty,zintz,t,xc,yc,zc,
114     1              x0,y0,z0,xi,yi,zi,xj,yj,zj,ni,nj
115      dimension   xin(Li+2,Lj+2),  yin(Li+2,Lj+2),  zin(Li+2,Lj+2)
116      dimension  xxin(Li+2,Lj+2), yyin(Li+2,Lj+2), zzin(Li+2,Lj+2)
117      dimension dxsdi(Li+1,Lj+1),dysdi(Li+1,Lj+1),dzsdi(Li+1,Lj+1)
118      dimension dxsdj(Li+1,Lj+1),dysdj(Li+1,Lj+1),dzsdj(Li+1,Lj+1)
119      dimension dxxdi(Li+1,Lj+1),dyydi(Li+1,Lj+1),dzzdi(Li+1,Lj+1)
120      dimension dxxdj(Li+1,Lj+1),dyydj(Li+1,Lj+1),dzzdj(Li+1,Lj+1)
121      dimension Nxyz(3)
122      data zero  /0.0d+00/
123      data one   /1.0e+00/
124      data rln10 /2.30258d+00/
125c
126      tol=rln10*itol
127      nder = 1
128c
129      xc=center(1)
130      yc=center(2)
131      zc=center(3)
132c
133c     ----- ishell -----
134c
135      xi=xyzi(1)
136      yi=xyzi(2)
137      zi=xyzi(3)
138      lit = Li + 1
139      maxi=lit*(lit+1)/2
140c
141      litder=lit+nder
142c
143c     ----- jshell -----
144c
145      xj=xyzj(1)
146      yj=xyzj(2)
147      zj=xyzj(3)
148      ljt = Lj + 1
149      maxj=ljt*(ljt+1)/2
150c
151      ljtder=ljt+nder
152c
153      rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2
154      nroots=(lit+ljt-2)/2+1
155      if(nroots.gt.maxrys) then
156         write(luout,9997) maxrys,lit,ljt,nroots
157         call errquit('hnd_dipder: angular momentum too high',555,
158     &       BASIS_ERR)
159      endif
160c
161c     ----- i primitive -----
162c
163      do 7000 ig=1,i_nprim
164      ai=expi(ig)
165      arri=ai*rr
166      axi=ai*xi
167      ayi=ai*yi
168      azi=ai*zi
169      csi=coefi(ig,i_ngen)
170c
171c     ----- j primitive -----
172c
173      do 6000 jg=1,j_nprim
174      aj=expj(jg)
175      aa=ai+aj
176      aa1=one/aa
177      dum=aj*arri*aa1
178      if(dum.gt.tol) go to 6000  ! the integral is zero
179      fac= exp(-dum)
180      csj=coefj(jg,j_ngen)
181      ax=(axi+aj*xj)*aa1
182      ay=(ayi+aj*yj)*aa1
183      az=(azi+aj*zj)*aa1
184c
185c     ----- density factor -----
186c
187      dij=fac*csi*csj
188c
189c     ----- dipole moment integrals -----
190c
191      t = sqrt(aa)
192      t1=one/t
193      x0=ax
194      y0=ay
195      z0=az
196      do 370 j=1,ljtder
197      nj=j
198      do 370 i=1,litder
199      ni=i
200      call hnd_dipxyz
201       xin(i,j)=xint*t1
202       yin(i,j)=yint*t1
203       zin(i,j)=zint*t1
204      xxin(i,j)=xintx*t1
205      yyin(i,j)=yinty*t1
206      zzin(i,j)=zintz*t1
207  370 continue
208c
209      call hnd_deriaj(dxsdi,dysdi,dzsdi,dxsdj,dysdj,dzsdj,
210     *                xin,yin,zin,
211     *                lit,ljt,ai,aj)
212      call hnd_deriaj(dxxdi,dyydi,dzzdi,dxxdj,dyydj,dzzdj,
213     *                xxin,yyin,zzin,
214     *                lit,ljt,ai,aj)
215      ij=0
216      do 390 i=1,maxi
217      call getNxyz(Li,i,Nxyz)
218      ix = Nxyz(1) + 1
219      iy = Nxyz(2) + 1
220      iz = Nxyz(3) + 1
221      do 390 j=1,maxj
222      call getNxyz(Lj,j,Nxyz)
223      jx = Nxyz(1) + 1
224      jy = Nxyz(2) + 1
225      jz = Nxyz(3) + 1
226      ij=ij+1
227c
228c  First do derivative wrt the first atom
229c
230      dumxdx = dxxdi(ix,jx)*  yin(iy,jy)*  zin(iz,jz)
231      dumydx = dxsdi(ix,jx)* yyin(iy,jy)*  zin(iz,jz)
232      dumzdx = dxsdi(ix,jx)*  yin(iy,jy)* zzin(iz,jz)
233      dumxdy =  xxin(ix,jx)*dysdi(iy,jy)*  zin(iz,jz)
234      dumydy =   xin(ix,jx)*dyydi(iy,jy)*  zin(iz,jz)
235      dumzdy =   xin(ix,jx)*dysdi(iy,jy)* zzin(iz,jz)
236      dumxdz =  xxin(ix,jx)*  yin(iy,jy)*dzsdi(iz,jz)
237      dumydz =   xin(ix,jx)* yyin(iy,jy)*dzsdi(iz,jz)
238      dumzdz =   xin(ix,jx)*  yin(iy,jy)*dzzdi(iz,jz)
239c
240      didij(ij,1) = didij(ij,1) + dij*dumxdx
241      didij(ij,2) = didij(ij,2) + dij*dumydx
242      didij(ij,3) = didij(ij,3) + dij*dumzdx
243      didij(ij,4) = didij(ij,4) + dij*dumxdy
244      didij(ij,5) = didij(ij,5) + dij*dumydy
245      didij(ij,6) = didij(ij,6) + dij*dumzdy
246      didij(ij,7) = didij(ij,7) + dij*dumxdz
247      didij(ij,8) = didij(ij,8) + dij*dumydz
248      didij(ij,9) = didij(ij,9) + dij*dumzdz
249c
250c  Next do derivative wrt the second atom
251c
252      dumxdx = dxxdj(ix,jx)*  yin(iy,jy)*  zin(iz,jz)
253      dumydx = dxsdj(ix,jx)* yyin(iy,jy)*  zin(iz,jz)
254      dumzdx = dxsdj(ix,jx)*  yin(iy,jy)* zzin(iz,jz)
255      dumxdy =  xxin(ix,jx)*dysdj(iy,jy)*  zin(iz,jz)
256      dumydy =   xin(ix,jx)*dyydj(iy,jy)*  zin(iz,jz)
257      dumzdy =   xin(ix,jx)*dysdj(iy,jy)* zzin(iz,jz)
258      dumxdz =  xxin(ix,jx)*  yin(iy,jy)*dzsdj(iz,jz)
259      dumydz =   xin(ix,jx)* yyin(iy,jy)*dzsdj(iz,jz)
260      dumzdz =   xin(ix,jx)*  yin(iy,jy)*dzzdj(iz,jz)
261c
262      djdij(ij,1) = djdij(ij,1) + dij*dumxdx
263      djdij(ij,2) = djdij(ij,2) + dij*dumydx
264      djdij(ij,3) = djdij(ij,3) + dij*dumzdx
265      djdij(ij,4) = djdij(ij,4) + dij*dumxdy
266      djdij(ij,5) = djdij(ij,5) + dij*dumydy
267      djdij(ij,6) = djdij(ij,6) + dij*dumzdy
268      djdij(ij,7) = djdij(ij,7) + dij*dumxdz
269      djdij(ij,8) = djdij(ij,8) + dij*dumydz
270      djdij(ij,9) = djdij(ij,9) + dij*dumzdz
271  390 continue
272c
273 6000 continue
274 7000 continue
275      return
276c
277 9997 format(' in -hnd_stvint- the rys quadrature is not implemented',
278     1       ' beyond -nroots- = ',i2,/,
279     2       ' lit,ljt,nroots= ',3i3)
280      end
281