1c
2c $Id$
3c
4c Modified from HONDO
5c
6      subroutine hnd_hlf_ij2(xyzi,expi,coefi,i_nprim,i_ngen, Li,
7     1  xyzj,expj,coefj, j_nprim, j_ngen, Lj, kat, dvij, ddvij,
8     2  zan, xyz, nder, nint, scr,lscr)
9c
10c This is a wrapper routine, setting up scratch blocks used in actual
11c integral routine
12c
13      implicit none
14      integer i_nprim  ! [input] num. prims on function i
15      integer i_ngen   ! [input] num general conts on func. i
16      integer Li       ! [input] angular momentum of func. i
17      integer j_nprim  ! [input] num. prims on function j
18      integer j_ngen   ! [input] num general conts on func. j
19      integer Lj       ! [input] angular momentum of func. j
20      integer nder     ! [input] 1=1rst der; 2=2nd der
21      integer nint     ! [input] number of base integrals
22      integer kat      ! [input] lexical number of an atom
23      integer lscr     ! [input] length of scratch space
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 xyz(3,*)  ! [input] all atom positions
31      double precision dvij(nint,3) ! [output] 1rst. der integrals
32      double precision ddvij(nint,3,3) ! [output] 2nd der integrals
33      double precision zan(*)   ! [input] nuclear charges
34      double precision scr(lscr) ! [input] scratch buffer
35c
36# include "hnd_pointers.fh"
37c
38c Use scratch for temps in integral calculation
39c Scratch needs are dependent on nder:
40c nder=1 (6 * (Li+1)*(Lj+1)*((Li+Lj+2)/2+1) + 3 0.0d0 length blocks
41c nder=2 (9 * (Li+1)*(Lj+1)*((Li+Lj+2)/2+1)
42c
43c The pointers are defined in hnd_pointers.fh
44c
45      call hnd_hlf_ij21(xyzi,expi,coefi,i_nprim,i_ngen, Li,
46     2     xyzj,expj,coefj, j_nprim, j_ngen, Lj, kat, dvij, ddvij,
47     3     zan, xyz, nder, nint,
48     4     scr(tvp(22)),scr(tvp(23)),scr(tvp(24)),
49     5     scr(tvp(25)),scr(tvp(26)),scr(tvp(27)),
50     6     scr(tvp(46)),scr(tvp(47)),scr(tvp(48)))
51c
52      return
53      end
54c
55      subroutine hnd_hlf_ij21(xyzi,expi,coefi,i_nprim,i_ngen, Li,
56     1  xyzj,expj,coefj, j_nprim, j_ngen, Lj, kat, dvij, ddvij,
57     2  zan, xyz, nder, nint,xv,yv,zv,dxv,dyv,dzv,ddxv,ddyv,ddzv)
58c
59c  This is the routine that actually computes the 1rst and 2nd
60c  derivatives for the Helman-Feynman term.  It is assumed that the buffers
61c  for the integrals have been 0.0d0ed before entering this routine.
62c
63      implicit none
64#include "stdio.fh"
65#include "hnd_rys.fh"
66#include "hnd_tol.fh"
67#include "errquit.fh"
68      integer i_nprim  ! [input] num. prims on function i
69      integer i_ngen   ! [input] num general conts on func. i
70      integer Li       ! [input] angular momentum of func. i
71      integer j_nprim  ! [input] num. prims on function j
72      integer j_ngen   ! [input] num general conts on func. j
73      integer Lj       ! [input] angular momentum of func. j
74      integer nder     ! [input] 1=1rst der; 2=2nd der
75      integer nint     ! [input] number of base integrals
76      integer kat      ! [input] lexical number of an atom
77      double precision xyzi(3)  ! [input] position of center i
78      double precision expi(i_nprim) ! [input] exponents on i
79      double precision coefi(i_nprim,i_ngen) ! [input] i coeffs
80      double precision xyzj(3)  ! [input] position of center j
81      double precision expj(j_nprim)  ! [input] exponents on j
82      double precision coefj(j_nprim,j_ngen)  ! [input] j coeffs
83      double precision xyz(3,*)  ! [input] all atom positions
84      double precision dvij(nint,3) ! [output] 1rst. der integrals
85      double precision ddvij(nint,3,3) ! [output] 2nd der integrals
86      double precision zan(*)   ! [input] nuclear charges
87c
88      character*8 errmsg
89      common/hnd_xyzder/xint,yint,zint,tx,x0,y0,z0,xi,yi,zi,xj,yj,zj,
90     1                             ni,nj,cx,cy,cz
91      double precision xint, yint, zint, tx, x0, y0, z0, xi, yi, zi
92      double precision xj, yj, zj, cx, cy, cz
93      double precision ijx, ijy, ijz
94      double precision rln10, tol, rr, ai, aj, arri
95      double precision axi, ayi, azi, csi, cpi, cdi, cfi, cgi
96      double precision aa, aa1, dum, fac, csj, cpj, cdj, cfj, cgj
97      double precision ax, ay, az, dum1, dum2, pij
98      double precision dumx, dumy, dumz, dumxx, dumyy, dumzz
99      double precision dumxy, dumxz, dumyz
100      double precision pi212, aax, aay, aaz, znuc
101      double precision uu, u2, u4, ww, w2, w4, tt, xv, yv, zv
102      double precision dxv, dyv, dzv, ddxv, ddyv, ddzv
103      integer ni, nj, iroot, Nxyz
104      integer lit,maxi, ljt, maxj
105      integer ig, jg, ij, i, j, ix, iy, iz, jx, jy, jz
106      dimension   xv(Li+1,Lj+1,*),  yv(Li+1,Lj+1,*),  zv(Li+1,Lj+1,*)
107      dimension  dxv(Li+1,Lj+1,*), dyv(Li+1,Lj+1,*), dzv(Li+1,Lj+1,*)
108      dimension ddxv(Li+1,Lj+1,*),ddyv(Li+1,Lj+1,*),ddzv(Li+1,Lj+1,*)
109      dimension Nxyz(3)
110      dimension w2(maxrys),w4(maxrys)
111      dimension errmsg(3)
112      data errmsg /'program ','stop in ','-hlfspd-'/
113      data rln10  /2.30258d+00/
114      data pi212  /1.1283791670955d+00/
115c
116      tol =rln10*itol
117c
118c     ----- calculate -helfey- term -----
119c
120c     ----- ishell -----
121c
122      xi=xyzi(1)
123      yi=xyzi(2)
124      zi=xyzi(3)
125      lit = Li + 1
126      maxi=lit*(lit+1)/2
127c
128c     ----- jshell -----
129c
130      xj=xyzj(1)
131      yj=xyzj(2)
132      zj=xyzj(3)
133      ljt = Lj + 1
134      maxj=ljt*(ljt+1)/2
135c
136      rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2
137      nroots=(lit+ljt+nder-2)/2+1
138      if(nroots.gt.maxrys) then
139         write(luout,9997) maxrys,lit,ljt,nroots
140         call hnd_hnderr(3,errmsg)
141      endif
142c
143c     ----- i primitive -----
144c
145      do 7000 ig=1,i_nprim
146      ai=expi(ig)
147      arri=ai*rr
148      axi=ai*xi
149      ayi=ai*yi
150      azi=ai*zi
151      csi=coefi(ig,i_ngen)
152c
153c     ----- j primitive -----
154c
155      do 6000 jg=1,j_nprim
156      aj=expj(jg)
157      aa=ai+aj
158      aa1=1.0d0/aa
159      dum=aj*arri*aa1
160      if(dum.gt.tol) go to 6000
161      fac= exp(-dum)
162      csj=coefj(jg,j_ngen)
163      ax=(axi+aj*xj)*aa1
164      ay=(ayi+aj*yj)*aa1
165      az=(azi+aj*zj)*aa1
166c
167c     ----- density factor -----
168c
169      pij=fac*csi*csj*pi212*aa1
170c
171c     ----- hellmann-feynman term -----
172c
173      aax=aa*ax
174      aay=aa*ay
175      aaz=aa*az
176c
177c     ----- kat -----
178c
179      znuc=-zan(kat)
180      cx=xyz(1,kat)
181      cy=xyz(2,kat)
182      cz=xyz(3,kat)
183      yy=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2)
184      call hnd_droot
185      do 420 iroot=1,nroots
186      uu=u9(iroot)*aa
187      u2=uu
188      u4=uu*uu
189      ww=w9(iroot)*znuc
190      w2(iroot)=ww*u2*2.0d0
191      w4(iroot)=ww*u4*4.0d0
192      tt=1.0d0/(aa+uu)
193      tx= sqrt(tt)
194      x0=(aax+uu*cx)*tt
195      y0=(aay+uu*cy)*tt
196      z0=(aaz+uu*cz)*tt
197      do 410 j=1,ljt
198      nj=j
199      do 410 i=1,lit
200      ni=i
201
202      call hnd_sxyz
203
204      xv(i,j,iroot)=xint
205      yv(i,j,iroot)=yint
206      zv(i,j,iroot)=zint
207
208      call hnd_dervxyz(1)
209
210      dxv(i,j,iroot)=xint
211      dyv(i,j,iroot)=yint
212      dzv(i,j,iroot)=zint
213
214      if (nder .gt. 1) then
215         call hnd_dervxyz(2)
216
217         ddxv(i,j,iroot)=xint
218         ddyv(i,j,iroot)=yint
219         ddzv(i,j,iroot)=zint
220      endif
221  410 continue
222  420 continue
223c
224      ij=0
225      do 450 i=1,maxi
226      call getNxyz(Li,i,Nxyz)
227      ix = Nxyz(1) + 1
228      iy = Nxyz(2) + 1
229      iz = Nxyz(3) + 1
230      do 440 j=1,maxj
231      call getNxyz(Lj,j,Nxyz)
232      jx = Nxyz(1) + 1
233      jy = Nxyz(2) + 1
234      jz = Nxyz(3) + 1
235      ij=ij+1
236      if (nder.eq.1) then
237        dumx=0.0d0
238        dumy=0.0d0
239        dumz=0.0d0
240        do 430 iroot=1,nroots
241           dumx = dumx + dxv(ix,jx,iroot)*  yv(iy,jy,iroot)*
242     1             zv(iz,jz,iroot)*w2(iroot)
243           dumy = dumy + xv(ix,jx,iroot)*  dyv(iy,jy,iroot)*
244     1             zv(iz,jz,iroot)*w2(iroot)
245           dumz = dumz + xv(ix,jx,iroot)*   yv(iy,jy,iroot)*
246     1            dzv(iz,jz,iroot)*w2(iroot)
247  430   continue
248        dumx =dumx *pij
249        dumy =dumy *pij
250        dumz =dumz *pij
251        dvij(ij,1)=dvij(ij,1)+dumx
252        dvij(ij,2)=dvij(ij,2)+dumy
253        dvij(ij,3)=dvij(ij,3)+dumz
254      elseif (nder.eq.2) then
255        dumxx=0.0d0
256        dumyy=0.0d0
257        dumzz=0.0d0
258        dumxy=0.0d0
259        dumxz=0.0d0
260        dumyz=0.0d0
261        do 431 iroot=1,nroots
262           dum=xv(ix,jx,iroot)*yv(iy,jy,iroot)*zv(iz,jz,iroot)*w2(iroot)
263           dumxx=dumxx-dum+ ddxv(ix,jx,iroot)*  yv(iy,jy,iroot)*
264     1             zv(iz,jz,iroot)*w4(iroot)
265           dumyy=dumyy-dum+   xv(ix,jx,iroot)*ddyv(iy,jy,iroot)*
266     1             zv(iz,jz,iroot)*w4(iroot)
267           dumzz=dumzz-dum+   xv(ix,jx,iroot)*  yv(iy,jy,iroot)*
268     1           ddzv(iz,jz,iroot)*w4(iroot)
269           dumxy=dumxy+ dxv(ix,jx,iroot)* dyv(iy,jy,iroot)*
270     1             zv(iz,jz,iroot)*w4(iroot)
271           dumxz=dumxz+ dxv(ix,jx,iroot)*  yv(iy,jy,iroot)*
272     1            dzv(iz,jz,iroot)*w4(iroot)
273           dumyz=dumyz+  xv(ix,jx,iroot)* dyv(iy,jy,iroot)*
274     1            dzv(iz,jz,iroot)*w4(iroot)
275  431   continue
276        dumxx=dumxx*pij
277        dumyy=dumyy*pij
278        dumzz=dumzz*pij
279        dumxy=dumxy*pij
280        dumxz=dumxz*pij
281        dumyz=dumyz*pij
282        ddvij(ij,1,1)=ddvij(ij,1,1)+dumxx
283        ddvij(ij,1,2)=ddvij(ij,1,2)+dumxy
284        ddvij(ij,1,3)=ddvij(ij,1,3)+dumxz
285        ddvij(ij,2,1)=ddvij(ij,2,1)+dumxy
286        ddvij(ij,2,2)=ddvij(ij,2,2)+dumyy
287        ddvij(ij,2,3)=ddvij(ij,2,3)+dumyz
288        ddvij(ij,3,1)=ddvij(ij,3,1)+dumxz
289        ddvij(ij,3,2)=ddvij(ij,3,2)+dumyz
290        ddvij(ij,3,3)=ddvij(ij,3,3)+dumzz
291      endif
292c
293  440 continue
294  450 continue
295c
296 6000 continue
297 7000 continue
298c
299      return
300 9997 format(' in -hlf- , the rys quadrature is not implemented',
301     1       ' beyond -nroots- = ',i3,/,
302     2       ' lit,ljt,nroots = ',3i3)
303      end
304