1      subroutine hnd_giaotv10(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,
2     & expj,coefj,j_nprim,j_ngen,Lj,coord,zan,nat,nint,tv10,scr,lscr)
3c
4c $Id$
5c
6c     ----- Wrapper routine that sets the sizes of scratch blocks -----
7c
8      implicit double precision (a-h,o-z)
9#include "hnd_pointers.fh"
10      dimension scr(lscr)
11      dimension xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim)
12      dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen)
13      dimension tv10(nint,3),coord(3,nat),zan(nat)
14c
15c Use scratch for temps in integral calculation
16c Scratch needs are
17c 6*(Li+2)*(Lj+3) + 6*(Li+2)*(Lj+1)*((Li+Lj+1)/2+1)
18c
19c The pointers are defined in hdn_pointers.fh
20c
21      call hnd_giaotv101(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,expj,
22     1 coefj,j_nprim,j_ngen,Lj,coord,zan,nat,nint,tv10,scr(gh01(1)),
23     2 scr(gh01(12)),scr(gh01(2)),scr(gh01(3)),scr(gh01(4)),
24     3 scr(gh01(5)),scr(gh01(6)),scr(gh01(7)),scr(gh01(8)),scr(gh01(9)),
25     4 scr(gh01(10)),scr(gh01(11)))
26c
27      return
28      end
29c
30      subroutine hnd_giaotv101(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,
31     1                 expj,coefj,j_nprim,j_ngen,Lj,coord,zan,nat,nint,
32     2                 tv10,xs,ys,zs,xxs,yys,zzs,xv,yv,zv,xxv,yyv,zzv)
33      implicit double precision (a-h,o-z)
34#include "nwc_const.fh"
35#include "hnd_rys.fh"
36#include "hnd_tol.fh"
37#include "stdio.fh"
38      integer NoKinetic             ! FA-11-08-10
39      common /skipKinetic/NoKinetic ! FA-11-08-10
40
41      common/hnd_xyzder/xint,yint,zint,t,x0,y0,z0,xi,yi,zi,xj,yj,zj,
42     1                  ni,nj,ccx,ccy,ccz
43      dimension Nxyz(3),xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim)
44      dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen)
45      dimension tv10(nint,3),coord(3,nat),zan(nat)
46      dimension  xs(Li+2,Lj+3),ys(Li+2,Lj+3),zs(Li+2,Lj+3)
47      dimension  xv(Li+2,Lj+1,*),yv(Li+2,Lj+1,*),zv(Li+2,Lj+1,*)
48      dimension xxs(Li+2,Lj+3),yys(Li+2,Lj+3),zzs(Li+2,Lj+3)
49      dimension xxv(Li+2,Lj+1,*),yyv(Li+2,Lj+1,*),zzv(Li+2,Lj+1,*)
50      data rln10  /2.30258d+00/
51      data pi212  /1.1283791670955d+00/
52c
53      tol=rln10*itol
54c
55c     Zero integral array
56c
57      call dcopy(nint*3,0.0d0,0,tv10,1)
58c
59c     ----- ishell -----
60c
61      xi=xyzi(1)
62      yi=xyzi(2)
63      zi=xyzi(3)
64      lit = Li + 1
65      maxi = lit*(lit+1)/2
66      litmod=lit+1
67c
68c     ----- jshell -----
69c
70      xj=xyzj(1)
71      yj=xyzj(2)
72      zj=xyzj(3)
73      ljt = Lj + 1
74      maxj = ljt*(ljt+1)/2
75      ljtmod=ljt+2
76c
77      qijx=yi*zj-zi*yj
78      qijy=zi*xj-xi*zj
79      qijz=xi*yj-yi*xj
80      tijx=xi-xj
81      tijy=yi-yj
82      tijz=zi-zj
83c
84      rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2
85c
86      nroots=(1+lit+ljt-2)/2+1
87      if(nroots.gt.maxrys) then
88         write(luout,9997) maxrys,lit,ljt,nroots
89         call errquit('hnd_giaotv10: need higher Rys rt',nroots,INT_ERR)
90      endif
91c
92c     ----- i primitive -----
93c
94      do ig=1, i_nprim
95      ai=expi(ig)
96      arri=ai*rr
97      axi=ai*xi
98      ayi=ai*yi
99      azi=ai*zi
100      csi=coefi(ig,i_ngen)
101c
102c     ----- j primitive -----
103c
104      do jg=1,j_nprim
105      aj=expj(jg)
106      aa=ai+aj
107      aa1=1.0d0/aa
108      dum=aj*arri*aa1
109      if(dum.gt.tol) goto 1000
110      fac= exp(-dum)
111      csj=coefj(jg,j_ngen)
112      ax=(axi+aj*xj)*aa1
113      ay=(ayi+aj*yj)*aa1
114      az=(azi+aj*zj)*aa1
115c
116c     ----- density factor -----
117c
118      cij=csi*csj*fac
119c
120c     ----- kinetic energy -----
121c
122      if (NoKinetic .eq. 1) goto 20 ! SKIP KINETIC ENERGY-FA-11-08-10
123      t = sqrt(aa1)
124      x0=ax
125      y0=ay
126      z0=az
127      do j=1,ljtmod
128         nj=j
129         do i=1,litmod
130            ni=i
131            call hnd_sxyz
132            xs(i,j)=xint*t
133            ys(i,j)=yint*t
134            zs(i,j)=zint*t
135         enddo
136         do i=1,lit
137            xxs(i,j)=xs(i+1,j)
138            yys(i,j)=ys(i+1,j)
139            zzs(i,j)=zs(i+1,j)
140         enddo
141      enddo
142c
143      do i=1,lit
144         xv(i,1,1)=(xs(i,1)-xs(i,3)*(aj+aj))*aj
145         yv(i,1,1)=(ys(i,1)-ys(i,3)*(aj+aj))*aj
146         zv(i,1,1)=(zs(i,1)-zs(i,3)*(aj+aj))*aj
147         xxv(i,1,1)=(xxs(i,1)-xxs(i,3)*(aj+aj))*aj
148         yyv(i,1,1)=(yys(i,1)-yys(i,3)*(aj+aj))*aj
149         zzv(i,1,1)=(zzs(i,1)-zzs(i,3)*(aj+aj))*aj
150         if (ljt.gt.1) then
151            xv(i,2,1)=(xs(i,2)*dble(2+2-1)-xs(i,4)*(aj+aj))*aj
152            yv(i,2,1)=(ys(i,2)*dble(2+2-1)-ys(i,4)*(aj+aj))*aj
153            zv(i,2,1)=(zs(i,2)*dble(2+2-1)-zs(i,4)*(aj+aj))*aj
154            xxv(i,2,1)=(xxs(i,2)*dble(2+2-1)-xxs(i,4)*(aj+aj))*aj
155            yyv(i,2,1)=(yys(i,2)*dble(2+2-1)-yys(i,4)*(aj+aj))*aj
156            zzv(i,2,1)=(zzs(i,2)*dble(2+2-1)-zzs(i,4)*(aj+aj))*aj
157         endif
158         do j=3,ljt
159            xv(i,j,1)=(xs(i,j)*dble(j+j-1)-xs(i,j+2)*(aj+aj))*aj
160     1                 -xs(i,j-2)*dble(((j-1)*(j-2))/2)
161            yv(i,j,1)=(ys(i,j)*dble(j+j-1)-ys(i,j+2)*(aj+aj))*aj
162     1                 -ys(i,j-2)*dble(((j-1)*(j-2))/2)
163            zv(i,j,1)=(zs(i,j)*dble(j+j-1)-zs(i,j+2)*(aj+aj))*aj
164     1                 -zs(i,j-2)*dble(((j-1)*(j-2))/2)
165            xxv(i,j,1)=(xxs(i,j)*dble(j+j-1)-xxs(i,j+2)*(aj+aj))*aj
166     1                 -xxs(i,j-2)*dble(((j-1)*(j-2))/2)
167            yyv(i,j,1)=(yys(i,j)*dble(j+j-1)-yys(i,j+2)*(aj+aj))*aj
168     1                 -yys(i,j-2)*dble(((j-1)*(j-2))/2)
169            zzv(i,j,1)=(zzs(i,j)*dble(j+j-1)-zzs(i,j+2)*(aj+aj))*aj
170     1                 -zzs(i,j-2)*dble(((j-1)*(j-2))/2)
171         enddo
172      enddo
173c
174      ij=0
175      do j=1,maxj
176         call getNxyz(Lj,j,Nxyz)
177         jx = Nxyz(1) + 1
178         jy = Nxyz(2) + 1
179         jz = Nxyz(3) + 1
180         do i=1,maxi
181            call getNxyz(Li,i,Nxyz)
182            ix = Nxyz(1) + 1
183            iy = Nxyz(2) + 1
184            iz = Nxyz(3) + 1
185            dum = xv(ix,jx,1)* ys(iy,jy)  * zs(iz,jz)
186     1          + xs(ix,jx)  * yv(iy,jy,1)* zs(iz,jz)
187     2          + xs(ix,jx)  * ys(iy,jy)  * zv(iz,jz,1)
188            dumx=xxv(ix,jx,1)* ys(iy,jy)  * zs(iz,jz)
189     1          +xxs(ix,jx)  * yv(iy,jy,1)* zs(iz,jz)
190     2          +xxs(ix,jx)  * ys(iy,jy)  * zv(iz,jz,1)
191            dumy= xv(ix,jx,1)*yys(iy,jy)  * zs(iz,jz)
192     1          + xs(ix,jx)  *yyv(iy,jy,1)* zs(iz,jz)
193     2          + xs(ix,jx)  *yys(iy,jy)  * zv(iz,jz,1)
194            dumz= xv(ix,jx,1)* ys(iy,jy)  *zzs(iz,jz)
195     1          + xs(ix,jx)  * yv(iy,jy,1)*zzs(iz,jz)
196     2          + xs(ix,jx)  * ys(iy,jy)  *zzv(iz,jz,1)
197            ij=ij+1
198            tv10(ij,1)=tv10(ij,1)+cij*(qijx*dum+tijy*dumz-tijz*dumy)
199            tv10(ij,2)=tv10(ij,2)+cij*(qijy*dum+tijz*dumx-tijx*dumz)
200            tv10(ij,3)=tv10(ij,3)+cij*(qijz*dum+tijx*dumy-tijy*dumx)
201         enddo
202      enddo
203 20   continue ! SKIP KINETIC ENERGY-FA-11-08-10
204c
205c     ----- nuclear attraction -----
206c
207      aax=aa*ax
208      aay=aa*ay
209      aaz=aa*az
210      do ic=1,nat
211      znuc=-zan(ic)
212      cx=coord(1,ic)
213      cy=coord(2,ic)
214      cz=coord(3,ic)
215      yy=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2)
216      call hnd_droot
217      do iroot=1,nroots
218         uu=u9(iroot)*aa
219         ww=w9(iroot)*znuc
220         tt=1.0d0/(aa+uu)
221         t = sqrt(tt)
222         x0=(aax+uu*cx)*tt
223         y0=(aay+uu*cy)*tt
224         z0=(aaz+uu*cz)*tt
225         do j=1,ljt
226            nj=j
227            do i=1,litmod
228               ni=i
229               call hnd_sxyz
230               xv(i,j,iroot)=xint
231               yv(i,j,iroot)=yint
232               zv(i,j,iroot)=zint*ww
233            enddo
234            do i=1,lit
235               xxv(i,j,iroot)=xv(i+1,j,iroot)
236               yyv(i,j,iroot)=yv(i+1,j,iroot)
237               zzv(i,j,iroot)=zv(i+1,j,iroot)
238            enddo
239         enddo
240      enddo
241c
242      ij=0
243      do j=1,maxj
244         call getNxyz(Lj,j,Nxyz)
245         jx = Nxyz(1) + 1
246         jy = Nxyz(2) + 1
247         jz = Nxyz(3) + 1
248         do i=1,maxi
249            call getNxyz(Li,i,Nxyz)
250            ix = Nxyz(1) + 1
251            iy = Nxyz(2) + 1
252            iz = Nxyz(3) + 1
253            dum =0.0d0
254            dumx=0.0d0
255            dumy=0.0d0
256            dumz=0.0d0
257            do iroot=1,nroots
258            dum =dum + xv(ix,jx,iroot)* yv(iy,jy,iroot)* zv(iz,jz,iroot)
259            dumx=dumx+xxv(ix,jx,iroot)* yv(iy,jy,iroot)* zv(iz,jz,iroot)
260            dumy=dumy+ xv(ix,jx,iroot)*yyv(iy,jy,iroot)* zv(iz,jz,iroot)
261            dumz=dumz+ xv(ix,jx,iroot)* yv(iy,jy,iroot)*zzv(iz,jz,iroot)
262            enddo
263            dum =dum *(aa1*pi212)
264            dumx=dumx*(aa1*pi212)
265            dumy=dumy*(aa1*pi212)
266            dumz=dumz*(aa1*pi212)
267            ij=ij+1
268            tv10(ij,1)=tv10(ij,1)+cij*(qijx*dum+tijy*dumz-tijz*dumy)
269            tv10(ij,2)=tv10(ij,2)+cij*(qijy*dum+tijz*dumx-tijx*dumz)
270            tv10(ij,3)=tv10(ij,3)+cij*(qijz*dum+tijx*dumy-tijy*dumx)
271         enddo
272      enddo
273c
274c     ----- end loop over atoms -----
275c
276      enddo
277c
278c     ----- end loop over primitives -----
279c
280 1000 continue
281      enddo
282      enddo
283c
284      return
285 9997 format(' in -gitv10- the rys quadrature is not implemented',
286     1       ' beyond -nroots- = ',i2,/,' lit,ljt,nroots= ',3i3)
287      end
288