1      subroutine hnd_giaoh11(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,
2     & expj,coefj,j_nprim,j_ngen,Lj,catms,nat,nint,e11,scr,lscr,
3     & para,dia)
4c
5c $Id$
6c
7c     ----- Wrapper routine that sets the sizes of scratch blocks -----
8c
9      implicit double precision (a-h,o-z)
10#include "hnd_pointers.fh"
11      dimension scr(lscr)
12      logical para,dia
13      dimension xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim)
14      dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen)
15      dimension e11(nint,3,3,*)
16      dimension catms(3,nat)
17c
18c Use scratch for temps in integral calculation
19c Scratch needs are
20c 11*3*(Li+2)*(Lj+2)*((Li+Lj+3)/2+1)
21c
22c The pointers are defined in hdn_pointers.fh
23c
24      call hnd_giaoh111(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,expj,
25     1 coefj,j_nprim,j_ngen,Lj,catms,nat,nint,e11,para,dia,
26     2 scr(gh01(1)),scr(gh01(2)),scr(gh01(3)),scr(gh01(4)),
27     3 scr(gh01(5)),scr(gh01(6)),scr(gh01(7)),scr(gh01(8)),
28     4 scr(gh01(9)),scr(gh01(10)),scr(gh01(11)))
29c
30      return
31      end
32c
33      subroutine hnd_giaoh111(xyzi,expi,coefi,i_nprim,i_ngen,Li,xyzj,
34     1 expj,coefj,j_nprim,j_ngen,Lj,catms,nat,nint,e11,para,dia,
35     4 v0,dv0,v0d,v1,dv1,v11,dv11,v1d,v12,dv12,v11d)
36c
37      implicit double precision (a-h,o-z)
38#include "nwc_const.fh"
39#include "hnd_rys.fh"
40#include "hnd_tol.fh"
41#include "stdio.fh"
42      common/hnd_xyzder/xint,yint,zint,t,x0,y0,z0,xi,yi,zi,xj,yj,zj,
43     1                  ni,nj,cx,cy,cz
44      dimension Nxyz(3),xyzi(3),xyzj(3),expi(i_nprim),expj(j_nprim)
45      dimension coefi(i_nprim,i_ngen),coefj(j_nprim,j_ngen)
46      dimension e11(nint,3,3,*)
47      dimension tm(21),catms(3,nat)
48      dimension   v0(3,Li+2,Lj+2,*)    ! blocks for vx0,vy0,vz0
49      dimension  dv0(3,Li+2,Lj+2,*)    ! blocks for dvx0,dvy0,dvz0
50      dimension  v0d(3,Li+2,Lj+2,*)    ! blocks for vx0d,vy0d,vz0d
51      dimension   v1(3,Li+2,Lj+2,*)    ! blocks for vx1,vy1,vz1
52      dimension  dv1(3,Li+2,Lj+2,*)    ! blocks for dvx1,dvy1,dvz1
53      dimension  v11(3,Li+2,Lj+2,*)    ! blocks for vx1x,vy1y,vz1z
54      dimension dv11(3,Li+2,Lj+2,*)    ! blocks for dvx1x,dvy1y,dvz1z
55      dimension  v1d(3,Li+2,Lj+2,*)    ! blocks for vx1d,vy1d,vz1d
56      dimension  v12(3,Li+2,Lj+2,*)    ! blocks for xvx1,yvy1,zvz1
57      dimension dv12(3,Li+2,Lj+2,*)    ! blocks for xdvx1,ydvy1,zdvz1
58      dimension v11d(3,Li+2,Lj+2,*)    ! blocks for xvx1d,yvy1d,zvz1d
59      logical para,dia
60      data rln10  /2.30258d+00/
61      data pi212  /1.1283791670955d+00/
62c
63      tol=rln10*itol
64c
65c     Zero integral array
66c
67      call dcopy(nint*9*nat,0.0d0,0,e11,1)
68c
69c     ----- ishell -----
70c
71      xi=xyzi(1)
72      yi=xyzi(2)
73      zi=xyzi(3)
74      lit = Li + 1
75      maxi = lit*(lit+1)/2
76      litmod=lit+1
77c
78c     ----- jshell -----
79c
80      xj=xyzj(1)
81      yj=xyzj(2)
82      zj=xyzj(3)
83      ljt = Lj + 1
84      maxj = ljt*(ljt+1)/2
85      ljtmod=ljt+1
86c
87      qijx=yi*zj-zi*yj
88      qijy=zi*xj-xi*zj
89      qijz=xi*yj-yi*xj
90      tijx=xi   -   xj
91      tijy=yi   -   yj
92      tijz=zi   -   zj
93c
94      rr=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2
95c
96      nroots=(lit+ljt+3-2)/2+1
97      if(nroots.gt.maxrys) then
98         write(luout,9999) maxrys,lit,ljt,nroots
99         call errquit('hnd_giaoh11: need higher Rys rt',nroots,INT_ERR)
100      endif
101c
102c     ----- i primitive -----
103c
104      do ig=1, i_nprim
105      ai=expi(ig)
106      arri=ai*rr
107      axi=ai*xi
108      ayi=ai*yi
109      azi=ai*zi
110      csi=coefi(ig,i_ngen)
111c
112c     ----- j primitive -----
113c
114      do jg=1,j_nprim
115      aj=expj(jg)
116      aa=ai+aj
117      aa1=1.0d0/aa
118      dum=aj*arri*aa1
119      if(dum.gt.tol) goto 1000
120      fac= exp(-dum)
121      csj=coefj(jg,j_ngen)
122      ax=(axi+aj*xj)*aa1
123      ay=(ayi+aj*yj)*aa1
124      az=(azi+aj*zj)*aa1
125c
126c     ----- density factor -----
127c
128      cij=csi*csj*fac*pi212*aa1
129c
130c     ----- -h11- dia + paramagnetic terms -----
131c
132      aax=aa*ax
133      aay=aa*ay
134      aaz=aa*az
135      do kat=1,nat
136         cx=catms(1,kat)
137         cy=catms(2,kat)
138         cz=catms(3,kat)
139         xx=aa*((ax-cx)**2+(ay-cy)**2+(az-cz)**2)
140         yy=xx
141         call hnd_droot
142         do ir=1,nroots
143            uu=u9(ir)*aa
144            ww=w9(ir)
145            vv=ww
146            ww=ww*(uu+uu)
147            tt=1.0d0/(aa+uu)
148            t = sqrt(tt)
149            x0=(aax+uu*cx)*tt
150            y0=(aay+uu*cy)*tt
151            z0=(aaz+uu*cz)*tt
152            do j=1,ljtmod
153               nj=j
154               do i=1,litmod
155                  ni=i
156                  call hnd_sxyz
157c
158c     ----- for 1/r -----
159c
160                  if(kat.eq.1) then
161                      v0(1,i,j,ir)=xint     ! vx0
162                      v0(2,i,j,ir)=yint     ! vy0
163                      v0(3,i,j,ir)=zint*vv  ! vz0
164                  endif
165c
166c     ----- for x/r**3 -----
167c
168                  v1(1,i,j,ir)=xint         ! vx1
169                  v1(2,i,j,ir)=yint         ! vy1
170                  v1(3,i,j,ir)=zint*ww      ! vz1
171                  call hnd_dervxyz(1)
172                  dv1(1,i,j,ir)=xint        ! dvx1
173                  dv1(2,i,j,ir)=yint        ! dvy1
174                  dv1(3,i,j,ir)=zint*ww     ! dvz1
175               enddo ! i-loop
176            enddo    ! j-loop
177         enddo       ! ir-loop
178c
179         do ir=1,nroots
180            do j=1,ljt
181               do i=1,lit
182                    v11(1,i,j,ir)= v1(1,i,j+1,ir)   ! vx1x
183                    v11(2,i,j,ir)= v1(2,i,j+1,ir)   ! vy1y
184                    v11(3,i,j,ir)= v1(3,i,j+1,ir)   ! vz1z
185                   dv11(1,i,j,ir)=dv1(1,i,j+1,ir)   ! dvx1x
186                   dv11(2,i,j,ir)=dv1(2,i,j+1,ir)   ! dvy1y
187                   dv11(3,i,j,ir)=dv1(3,i,j+1,ir)   ! dvz1z
188               enddo ! i-loop
189            enddo    ! j-loop
190            do j=1,ljtmod
191               do i=1,lit
192                   v12(1,i,j,ir)= v1(1,i+1,j,ir)    ! xvx1
193                   v12(2,i,j,ir)= v1(2,i+1,j,ir)    ! yvy1
194                   v12(3,i,j,ir)= v1(3,i+1,j,ir)    ! zvz1
195                  dv12(1,i,j,ir)=dv1(1,i+1,j,ir)    ! xdvx1
196                  dv12(2,i,j,ir)=dv1(2,i+1,j,ir)    ! ydvy1
197                  dv12(3,i,j,ir)=dv1(3,i+1,j,ir)    ! zdvz1
198               enddo ! i-loop
199            enddo    ! j-loop
200         enddo       ! ir-loop
201c
202         do ir=1,nroots
203c
204c     ----- for 1/r derivatives wrt. centers -i- and -j- -----
205c
206         if (kat.eq.1) then
207c
208c     ----- derivatives with respect to xj ... -----
209c
210            do i=1,lit
211               v0d(1,i,1,ir)=-(-(aj+aj)*v0(1,i,2,ir))          ! vx0d and vx0
212               v0d(2,i,1,ir)=-(-(aj+aj)*v0(2,i,2,ir))          ! vy0d and vy0
213               v0d(3,i,1,ir)=-(-(aj+aj)*v0(3,i,2,ir))          ! vz0d and vz0
214               do j=2,ljt
215                  v0d(1,i,j,ir)=-(dble(j-1)*v0(1,i,j-1,ir)-    ! vx0d and vx0
216     &                              (aj+aj)*v0(1,i,j+1,ir))
217                  v0d(2,i,j,ir)=-(dble(j-1)*v0(2,i,j-1,ir)-    ! vy0d and vy0
218     &                              (aj+aj)*v0(2,i,j+1,ir))
219                  v0d(3,i,j,ir)=-(dble(j-1)*v0(3,i,j-1,ir)-    ! vz0d and vz0
220     &                              (aj+aj)*v0(3,i,j+1,ir))
221               enddo ! j-loop
222            enddo    ! i-loop
223c
224c     ----- derivatives with respect to xi ... -----
225c
226            do j=1,ljt
227               dv0(1,1,j,ir)=-(-(ai+ai)*v0(1,2,j,ir))          ! dvx0 and vx0
228               dv0(2,1,j,ir)=-(-(ai+ai)*v0(2,2,j,ir))          ! dvy0 and vy0
229               dv0(3,1,j,ir)=-(-(ai+ai)*v0(3,2,j,ir))          ! dvz0 and vz0
230               do i=2,lit
231                  dv0(1,i,j,ir)=-(dble(i-1)*v0(1,i-1,j,ir)-    ! dvx0 and vx0
232     &                              (ai+ai)*v0(1,i+1,j,ir))
233                  dv0(2,i,j,ir)=-(dble(i-1)*v0(2,i-1,j,ir)-    ! dvy0 and vy0
234     &                              (ai+ai)*v0(2,i+1,j,ir))
235                  dv0(3,i,j,ir)=-(dble(i-1)*v0(3,i-1,j,ir)-    ! vz0d and vz0
236     &                              (ai+ai)*v0(3,i+1,j,ir))
237               enddo ! i-loop
238            enddo    ! j-loop
239         endif
240c
241c     ----- d/dx ... operators -----
242c
243         do i=1,lit
244            v1d(1,i,1,ir)= (-(aj+aj)*v1(1,i,2,ir))             ! vx1d and vx1
245            v1d(2,i,1,ir)= (-(aj+aj)*v1(2,i,2,ir))             ! vy1d and vy1
246            v1d(3,i,1,ir)= (-(aj+aj)*v1(3,i,2,ir))             ! vz1d and vz1
247            do j=2,ljt
248               v1d(1,i,j,ir)= (dble(j-1)*v1(1,i,j-1,ir)-       ! vx1d and vx1
249     &                           (aj+aj)*v1(1,i,j+1,ir))
250               v1d(2,i,j,ir)= (dble(j-1)*v1(2,i,j-1,ir)-       ! vy1d and vy1
251     &                           (aj+aj)*v1(2,i,j+1,ir))
252               v1d(3,i,j,ir)= (dble(j-1)*v1(3,i,j-1,ir)-       ! vz1d and vz1
253     &                           (aj+aj)*v1(3,i,j+1,ir))
254            enddo ! j-loop
255         enddo    ! i-loop
256c
257         do i=1,lit
258            v11d(1,i,1,ir)= (-(aj+aj)*v12(1,i,2,ir))           ! xvx1d and xvx1
259            v11d(2,i,1,ir)= (-(aj+aj)*v12(2,i,2,ir))           ! yvy1d and yvy1
260            v11d(3,i,1,ir)= (-(aj+aj)*v12(3,i,2,ir))           ! zvz1d and zvz1
261            do j=2,ljt
262               v11d(1,i,j,ir)= (dble(j-1)*v12(1,i,j-1,ir)-     ! xvx1d and xvx1
263     &                            (aj+aj)*v12(1,i,j+1,ir))
264               v11d(2,i,j,ir)= (dble(j-1)*v12(2,i,j-1,ir)-     ! yvy1d and yvy1
265     &                            (aj+aj)*v12(2,i,j+1,ir))
266               v11d(3,i,j,ir)= (dble(j-1)*v12(3,i,j-1,ir)-     ! zvz1d and zvz1
267     &                            (aj+aj)*v12(3,i,j+1,ir))
268            enddo ! j-loop
269         enddo    ! i-loop
270c
271         enddo    ! ir-loop
272c
273         ij=0
274         do j=1,maxj
275            call getNxyz(Lj,j,Nxyz)
276            jx = Nxyz(1) + 1
277            jy = Nxyz(2) + 1
278            jz = Nxyz(3) + 1
279            do i=1,maxi
280               call getNxyz(Li,i,Nxyz)
281               ix = Nxyz(1) + 1
282               iy = Nxyz(2) + 1
283               iz = Nxyz(3) + 1
284               transx=0.0d0
285               transy=0.0d0
286               transz=0.0d0
287               call dcopy(21,0.0d0,0,tm,1)
288               do ir=1,nroots
289c
290                if(kat.eq.1) then
291c
292c     ----- translation invariance for nuclear derivatives of 1/r -----
293c
294                  transx=transx
295     1                   +dv0(1,ix,jx,ir)*v0(2,iy,jy,ir)*v0(3,iz,jz,ir)
296     2                   +dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1(3,iz,jz,ir)
297     3                   +v0d(1,ix,jx,ir)*v0(2,iy,jy,ir)*v0(3,iz,jz,ir)
298                  transy=transy
299     1                   +v0(1,ix,jx,ir)*dv0(2,iy,jy,ir)*v0(3,iz,jz,ir)
300     2                   +v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir)
301     3                   +v0(1,ix,jx,ir)*v0d(2,iy,jy,ir)*v0(3,iz,jz,ir)
302                  transz=transz
303     1                   +v0(1,ix,jx,ir)*v0(2,iy,jy,ir)*dv0(3,iz,jz,ir)
304     2                   +v1(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir)
305     3                   +v0(1,ix,jx,ir)*v0(2,iy,jy,ir)*v0d(3,iz,jz,ir)
306                endif
307c
308c     ----- for h(1,1) diamagnetic -----
309c
310c       xx = xx + dvx1x *   vy1 *   vz1
311c       xy = xy +  vx1x *  dvy1 *   vz1
312c       xz = xz +  vx1x *   vy1 *  dvz1
313c       yx = yx +  dvx1 *  vy1y *   vz1
314c       yy = yy +   vx1 * dvy1y *   vz1
315c       yz = yz +   vx1 *  vy1y *   vz1
316c       zx = zx +  dvx1 *   vy1 *  dvz1
317c       zy = zy +   vx1 *  dvy1 *  vz1z
318c       zz = zz +   vx1 *   vy1 * dvz1z
319c
320                if (dia) then
321              tm(1)=tm(1)+dv11(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1(3,iz,jz,ir)
322              tm(2)=tm(2)+v11(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir)
323              tm(3)=tm(3)+v11(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir)
324              tm(4)=tm(4)+dv1(1,ix,jx,ir)*v11(2,iy,jy,ir)*v1(3,iz,jz,ir)
325              tm(5)=tm(5)+v1(1,ix,jx,ir)*dv11(2,iy,jy,ir)*v1(3,iz,jz,ir)
326              tm(6)=tm(6)+v1(1,ix,jx,ir)*v11(2,iy,jy,ir)*dv1(3,iz,jz,ir)
327              tm(7)=tm(7)+dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v11(3,iz,jz,ir)
328              tm(8)=tm(8)+v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v11(3,iz,jz,ir)
329              tm(9)=tm(9)+v1(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv11(3,iz,jz,ir)
330                endif
331c
332c     ----- for h(1,1) paramagnetic -----
333c
334c               t10 = t10 +  vx1 * dvy1 * vz1d -  vx1 * vy1d * dvz1
335c               t11 = t11 + vx1d *  vy1 * vz1d - dvx1 *  vy1 * vz1d
336c               t12 = t12 + dvx1 * vy1d *  vz1 - vx1d * dvy1 *  vz1
337c
338                if (para) then
339                tm(10)=tm(10)
340     1                +v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1d(3,iz,jz,ir)
341     2                -v1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*dv1(3,iz,jz,ir)
342                tm(11)=tm(11)
343     1                +v1d(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir)
344     2                -dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1d(3,iz,jz,ir)
345                tm(12)=tm(12)
346     1                +dv1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*v1(3,iz,jz,ir)
347     2                -v1d(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir)
348c
349c       t13 = t10 +  xvx1 *  dvy1 *  vz1d -  xvx1 *  vy1d *  dvz1
350c       t14 = t11 + xvx1d *   vy1 *  dvz1 - xdvx1 *   vy1 *  vz1d
351c       t15 = t12 + xdvx1 *  vy1d *   vz1 - xvx1d *  dvy1 *   vz1
352c       t16 = t12 +   vx1 * ydvy1 *  vz1d -   vx1 * yvy1d *  dvz1
353c       t17 = t12 +  vx1d *  yvy1 *  dvz1 -  dvx1 *  yvy1 *  vz1d
354c       t18 = t12 +  dvx1 * yvy1d *   vz1 -  vx1d * ydvy1 *   vz1
355c       t19 = t12 +   vx1 *  dvy1 * zvz1d -   vx1 *  vy1d * zdvz1
356c       t20 = t12 +  vx1d *   vy1 * zdvz1 -  dvx1 *   vy1 * zvz1d
357c       t21 = t12 +  dvx1 *  vy1d *  zvz1 -  vx1d *  dvy1 *  zvz1
358c
359                tm(13)=tm(13)
360     1                +v12(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1d(3,iz,jz,ir)
361     2                -v12(1,ix,jx,ir)*v1d(2,iy,jy,ir)*dv1(3,iz,jz,ir)
362                tm(14)=tm(14)
363     1                +v11d(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv1(3,iz,jz,ir)
364     2                -dv12(1,ix,jx,ir)*v1(2,iy,jy,ir)*v1d(3,iz,jz,ir)
365                tm(15)=tm(15)
366     1                +dv12(1,ix,jx,ir)*v1d(2,iy,jy,ir)*v1(3,iz,jz,ir)
367     2                -v11d(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v1(3,iz,jz,ir)
368                tm(16)=tm(16)
369     1                +v1(1,ix,jx,ir)*dv12(2,iy,jy,ir)*v1d(3,iz,jz,ir)
370     2                -v1(1,ix,jx,ir)*v11d(2,iy,jy,ir)*dv1(3,iz,jz,ir)
371                tm(17)=tm(17)
372     1                +v1d(1,ix,jx,ir)*v12(2,iy,jy,ir)*dv1(3,iz,jz,ir)
373     2                -dv1(1,ix,jx,ir)*v12(2,iy,jy,ir)*v1d(3,iz,jz,ir)
374                tm(18)=tm(18)
375     1                +dv1(1,ix,jx,ir)*v11d(2,iy,jy,ir)*v1(3,iz,jz,ir)
376     2                -v1d(1,ix,jx,ir)*dv12(2,iy,jy,ir)*v1(3,iz,jz,ir)
377                tm(19)=tm(19)
378     1                +v1(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v11d(3,iz,jz,ir)
379     2                -v1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*dv12(3,iz,jz,ir)
380                tm(20)=tm(20)
381     1                +v1d(1,ix,jx,ir)*v1(2,iy,jy,ir)*dv12(3,iz,jz,ir)
382     2                -dv1(1,ix,jx,ir)*v1(2,iy,jy,ir)*v11d(3,iz,jz,ir)
383                tm(21)=tm(21)
384     1                +dv1(1,ix,jx,ir)*v1d(2,iy,jy,ir)*v12(3,iz,jz,ir)
385     2                -v1d(1,ix,jx,ir)*dv1(2,iy,jy,ir)*v12(3,iz,jz,ir)
386                endif
387c
388               enddo
389               ij=ij+1
390c
391c     ----- h(1,1) diamagnetic -----
392c
393               if (dia) then
394                  e11(ij,1,1,kat)=e11(ij,1,1,kat)+(tm(5)+tm(9))*cij
395                  e11(ij,1,2,kat)=e11(ij,1,2,kat)-tm(2)*cij
396                  e11(ij,1,3,kat)=e11(ij,1,3,kat)-tm(3)*cij
397                  e11(ij,2,1,kat)=e11(ij,2,1,kat)-tm(4)*cij
398                  e11(ij,2,2,kat)=e11(ij,2,2,kat)+(tm(1)+tm(9))*cij
399                  e11(ij,2,3,kat)=e11(ij,2,3,kat)-tm(6)*cij
400                  e11(ij,3,1,kat)=e11(ij,3,1,kat)-tm(7)*cij
401                  e11(ij,3,2,kat)=e11(ij,3,2,kat)-tm(8)*cij
402                  e11(ij,3,3,kat)=e11(ij,3,3,kat)+(tm(1)+tm(5))*cij
403               endif
404c
405c     ----- h(1,1) paramagnetic -----
406c
407c     ----- a bit unclear ... yx ... instead of ... xy ... -----
408c
409               if (para) then
410                  tm(1)=(qijx*tm(10)+tijy*tm(19)-tijz*tm(16))*cij
411                  tm(2)=(qijx*tm(11)+tijy*tm(20)-tijz*tm(17))*cij
412                  tm(3)=(qijx*tm(12)+tijy*tm(21)-tijz*tm(18))*cij
413                  tm(4)=(qijy*tm(10)+tijz*tm(13)-tijx*tm(19))*cij
414                  tm(5)=(qijy*tm(11)+tijz*tm(14)-tijx*tm(20))*cij
415                  tm(6)=(qijy*tm(12)+tijz*tm(15)-tijx*tm(21))*cij
416                  tm(7)=(qijz*tm(10)+tijx*tm(16)-tijy*tm(13))*cij
417                  tm(8)=(qijz*tm(11)+tijx*tm(17)-tijy*tm(14))*cij
418                  tm(9)=(qijz*tm(12)+tijx*tm(18)-tijy*tm(15))*cij
419c
420                  e11(ij,1,1,kat)=e11(ij,1,1,kat)+tm(1)
421                  e11(ij,1,2,kat)=e11(ij,1,2,kat)+tm(4)
422                  e11(ij,1,3,kat)=e11(ij,1,3,kat)+tm(7)
423                  e11(ij,2,1,kat)=e11(ij,2,1,kat)+tm(2)
424                  e11(ij,2,2,kat)=e11(ij,2,2,kat)+tm(5)
425                  e11(ij,2,3,kat)=e11(ij,2,3,kat)+tm(8)
426                  e11(ij,3,1,kat)=e11(ij,3,1,kat)+tm(3)
427                  e11(ij,3,2,kat)=e11(ij,3,2,kat)+tm(6)
428                  e11(ij,3,3,kat)=e11(ij,3,3,kat)+tm(9)
429               endif
430c
431c              if(kat.eq.1) then
432c                if(abs(transx*cij).ge.1.0d-08.or.abs(transy*cij)
433c    &              .ge.1.0d-08.or.abs(transz*cij).ge.1.0d-08) then
434c                  write(luout,9993)
435c                  call errquit('hnd_giaoh11: no transl inv',0,INT_ERR)
436c                endif
437c              endif
438c
439            enddo  ! j-loop final summation
440         enddo     ! i-loop final summation
441c
442      enddo        ! kat-loop
443c
444 1000 continue
445      enddo        ! jprim loop
446      enddo        ! iprim loop
447c
448      return
449 9999 format(' in -giaoh11- , the rys quadrature is not implemented',
450     1       ' beyond -nroots- = ',i3,/,' lit,ljt,nroots= ',3i3)
451 9993 format(' something wrong with translational',
452     1       ' invariance in -giaoh11- . stop. ')
453      end
454