1c $Id$
2c----------------------------------------------------
3C*
4C*  THESE ROUTINES SHIFT THE ANGULAR MOMENTUM
5C*
6C*         FROM POSITION 1 TO POSITION 2
7C*
8C*                   AND
9C*
10C*         FROM POSITION 3 TO POSITION 4
11c----------------------------------------------------
12c    for re-ordered basis set
13c   nqi.ge.nqj  and  nqk.ge.nql
14c
15c   other cases are not included here !
16c----------------------------------------------------
17      subroutine amtfer (a,b,n)
18      implicit none
19      double precision a(*),b(*)
20      integer n, i
21c
22c     local copy of tfer for automatic inlining by compiler
23c
24      do i = 1, n
25         b(i) = a(i)
26      enddo
27c
28      end
29      subroutine amshift(bl,nbls,l01,l02,npij,npkl,ngcd)
30      implicit real*8 (a-h,o-z)
31      character*11 scftype
32      character*8 where
33#include "texas_lpar.fh"
34      common /runtype/ scftype,where
35c
36      common /cpu/ intsize,iacc,icache,memreal
37c
38      COMMON/SHELL/LSHELLT,LSHELIJ,LSHELKL,LHELP,LCAS2(4),LCAS3(4)
39      common/obarai/
40     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
41     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
42     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbeg,klbeg
43c
44ccccc common /big/ bl(1)
45      common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2,
46     * ibfij1,ibfij2,ibfkl1,ibfkl2,
47     * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3,
48     * ibf3l,issss,
49     * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4,
50     * ixij,iyij,izij, iwij,ivij,iuij,isij
51c
52      common /memor4a/ ibf3l1,ibf3l2,ibf3l3,ibf3l4
53c only for first & second derivatives (for use in amshift):
54cccc  common /memor4b/ibuf0
55      common /memor4b/ider0,ider1,ider2
56c
57      common /memor5a/ iaa,ibb,icc,idd,icis,icjs,icks,icls,
58     * ixab,ixp,ixpn,ixpp,iabnia,iapb,i1apb,ifij,icij,isab,
59     * ixcd,ixq,ixqn,ixqq,icdnia,icpd,i1cpd,ifkl,ickl,iscd
60      common /memor5b/ irppq,
61     * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234,
62     * idx1,idx2,indx
63c
64      dimension bl(*)
65c
66C**********************************
67c* dimensions for "shifts"
68c
69      lsmx=max(lnij,lnkl)
70      lsjl=max(nfu(nqj+1),nfu(nql+1))
71c
72      lqij=nfu(nqij+1)
73      lqkl=nfu(nqkl+1)
74      lqmx=max(lqij,lqkl)
75c
76c* memory for array INDXX :
77c
78      mindxx=lnijkl
79      if(intsize.ne.1) mindxx=lnijkl/intsize+1
80c
81      call getmem(mindxx,indxx)
82      call make_indxx(bl(indxx),lni,lnj,lnk,lnl)
83c
84c------------------------------------------------
85c for odinary integrals :
86c
87          nbuf2=ibuf2
88          nbuf=ibuf
89          m=1
90c
91c for giao integral derivatives :
92          if(where.eq.'shif') then
93            nmr =ngcd*nbls*lnijkl
94            nbuf=ibuf+nmr
95            m=6
96          endif
97c
98c for gradient integral derivatives :
99          if(where.eq.'forc') then
100            m=9
101          endif
102c
103c for hessian  integral derivatives :
104          if(where.eq.'hess') then
105            m=45
106          endif
107c-
108          mnbls=m*nbls
109c-
110c------------------------------------------------
111c------------------------------------------------
112C************************************************
113C**   SPECIAL CASES WHERE THE SHIFTING OF ANGULAR
114C**   MOMENTUM IS NOT NEEDED AT ALL / (DS|SS)..(SS|SD),
115C**   (XS|YS),(XS|SY),(SX|YS),(SX|SY) /
116C
117      IF(NQIJ.EQ.NSIJ .AND. NQKL.EQ.NSKL) THEN
118c
119c         find apropriate matrices with l-shells
120c
121          ibfijx=ibfij1
122          if(lshelij.eq.2) ibfijx=ibfij2
123          ibfklx=ibfkl1
124          if(lshelkl.eq.2) ibfklx=ibfkl2
125          ibf2lx=ibf2l1
126          if(lcas2(2).eq.1) ibf2lx=ibf2l2
127          if(lcas2(3).eq.1) ibf2lx=ibf2l3
128          if(lcas2(4).eq.1) ibf2lx=ibf2l4
129c
130c-
131          incre =mnbls*lnijkl
132          incre2=mnbls*l01*l02
133c-
134            do 100 iqu=1,ngcd
135              jbuf =nbuf +(iqu-1)*incre
136              jbuf2=nbuf2+(iqu-1)*incre2
137              call noshift(l01,l02,mnbls,  bl(jbuf), bl(jbuf2),
138     *          bl(ibfijx),           bl(ibfklx),
139     *          bl(ibf2lx),
140     *          lqij,lqkl,
141     *          bl(indxx),lni*lnj,lnk*lnl)
142  100       continue
143         call retmem(1)
144         return
145      ENDIF
146CC
147ccccccccccccccccccccc
148c
149       call convr3(bl,m,nbls,npij,npkl,bl(idx1),bl(idx2),
150     *             bl(ixab),bl(ixcd),ixabn,ixcdn)
151c
152c****
153c
154       if (lshellt.eq.0) then
155c-
156        incre =mnbls*lnijkl
157        incre2=mnbls*l01*l02
158c-
159         do 200 iqu=1,ngcd
160           jbuf =nbuf +(iqu-1)*incre
161           jbuf2=nbuf2+(iqu-1)*incre2
162ccccc      if(where.ne.'forc') then
163           if(where.eq.'buff' .or. where.eq.'shif') then
164              call shift0l(bl(jbuf),bl(jbuf2),
165     *                     l01,l02,bl(iwij),lsmx,lsjl,
166     *                     bl(ixij),nfu(nqi+1),lqij,lnkl,mnbls,
167     *                     bl(ixabn),bl(ixcdn),
168     *                     bl(indxx),lni,lnj,lnk,lnl)
169           endif
170           if(where.eq.'forc') then
171cccccc        jbuf0=ibuf0+(iqu-1)*nbls*l01*l02
172              jbuf0=ider0+(iqu-1)*nbls*l01*l02
173              iwij0=iwij +mnbls*lsmx*lsjl
174              ixij0=ixij +mnbls*nfu(nqi+1)*lqij*lnkl
175              call shift0l_der1(bl(jbuf),bl(jbuf2),bl(jbuf0),
176     *                         l01,l02,bl(iwij),bl(iwij0),lsmx,lsjl,
177     *                         bl(ixij),bl(ixij0),nfu(nqi+1),lqij,lnkl,
178     *                         mnbls,nbls,
179     *                         bl(ixabn),bl(ixcdn),
180     *                         bl(indxx),lni,lnj,lnk,lnl)
181           endif
182           if(where.eq.'hess') then
183              jbuf0=ider0+(iqu-1)*nbls*l01*l02
184              jbuf1=ider1+(iqu-1)*9*nbls*l01*l02
185              iwij1=iwij +mnbls*lsmx*lsjl
186              ixij1=ixij +mnbls*nfu(nqi+1)*lqij*lnkl
187              iwij0=iwij1 +    9*nbls*lsmx*lsjl
188              ixij0=ixij1 +    9*nbls*nfu(nqi+1)*lqij*lnkl
189              call shift0l_der2(bl(jbuf),bl(jbuf2),bl(jbuf1),bl(jbuf0),
190     *                         l01,l02,bl(iwij),bl(iwij1),bl(iwij0),
191     *                         lsmx,lsjl,
192     *                         bl(ixij),bl(ixij1),bl(ixij0),
193     *                         nfu(nqi+1),lqij,lnkl, mnbls,nbls,
194     *                         bl(ixabn),bl(ixcdn),
195     *                         bl(indxx),lni,lnj,lnk,lnl)
196           endif
197  200    continue
198         call retmem(3)
199         return
200       endif
201c
202c- 1 l-shell
203       if (lshellt.eq.1) then
204c---
205          jbuf  = nbuf
206          jbuf2 =nbuf2
207          jbfijx=ibfij1
208          if(lshelij.eq.2) jbfijx=ibfij2
209          jbfklx=ibfkl1
210          if(lshelkl.eq.2) jbfklx=ibfkl2
211c---
212          call shift1l(bl(jbuf),bl(jbuf2),l01,l02,
213     *                 bl(jbfijx),bl(jbfklx),
214     *                 lqij,lqkl,
215     *                 bl(iwij),lsmx, bl(ivij),lsjl,
216     *                 bl(ixij),nfu(nqi+1),lnkl,bl(iyij),nfu(nqj+1),
217     *                 mnbls,
218     *                 bl(ixabn),bl(ixcdn),
219     *                 bl(indxx),lni,lnj,lnk,lnl)
220c
221          call retmem(3)
222          return
223       endif
224ccc
225c- 2 l-shell
226       if (lshellt.eq.2) then
227c-
228          jbuf  = nbuf
229          jbuf2 =nbuf2
230          jbfij1=ibfij1
231          jbfij2=ibfij2
232          jbfkl1=ibfkl1
233          jbfkl2=ibfkl2
234c-
235          jbfij3=ibfij3
236          jbfkl3=ibfkl3
237          jbf2l12=ibf2l1
238          if(lcas2(2).eq.1) jbf2l12=ibf2l2
239          jbf2l34=ibf2l3
240          if(lcas2(4).eq.1) jbf2l34=ibf2l4
241c---
242          call shift2l(bl(jbuf),bl(jbuf2),l01,l02,
243     *                 bl(jbfij1),bl(jbfij2),bl(jbfkl1),bl(jbfkl2),
244     *                 lqij,lqkl,
245     *                 bl(jbfij3),bl(jbfkl3),
246     *                 bl(jbf2l12),bl(jbf2l34),
247     *                 bl(iwij),lsmx, bl(ivij),bl(iuij),bl(isij),lsjl,
248     *          bl(ixij),nfu(nqi+1),lnkl,bl(iyij),bl(izij),nfu(nqj+1),
249     *                 bl(ix2l1),mnbls,
250     *                 bl(ixabn),bl(ixcdn),
251     *                 bl(indxx),lni,lnj,lnk,lnl)
252c-
253          call retmem(3)
254          return
255       endif
256ccc
257c- 3 l-shell
258c
259       if (lshellt.eq.3) then
260c---
261          jbuf  = nbuf
262          jbuf2 =nbuf2
263          jbfij1=ibfij1
264          jbfij2=ibfij2
265          jbfkl1=ibfkl1
266          jbfkl2=ibfkl2
267c-
268          jbfij3=ibfij3
269          jbfkl3=ibfkl3
270          jbf2l1=ibf2l1
271          jbf2l2=ibf2l2
272          jbf2l3=ibf2l3
273          jbf2l4=ibf2l4
274c-
275          jbf3l12=ibf3l1
276          ix3l12=ix3l1
277          if(lcas3(2).eq.1) then
278              jbf3l12=ibf3l2
279              ix3l12=ix3l2
280          endif
281          jbf3l34=ibf3l3
282          ix3l34=ix3l3
283          if(lcas3(4).eq.1) then
284              jbf3l34=ibf3l4
285              ix3l34=ix3l4
286          endif
287c---
288          call shift3l(bl(jbuf),bl(jbuf2),l01,l02,
289     *                 bl(jbfij1),bl(jbfij2),bl(jbfkl1),bl(jbfkl2),
290     *                 lqij,lqkl,
291     *                 bl(jbfij3),bl(jbfkl3),
292     *                 bl(jbf2l1),bl(jbf2l2),bl(jbf2l3),bl(jbf2l4),
293     *                 bl(jbf3l12),bl(jbf3l34),lqmx,
294     *                 bl(iwij),lsmx, bl(ivij),bl(iuij),bl(isij),lsjl,
295     *           bl(ixij),nfu(nqi+1),lnkl,bl(iyij),bl(izij),nfu(nqj+1),
296     *                 bl(ix2l1),bl(ix2l2),bl(ix2l3),bl(ix2l4),
297     *                 bl(ix3l12),bl(ix3l34),mnbls,
298     *                 bl(ixabn),bl(ixcdn),
299     *                 bl(indxx),lni,lnj,lnk,lnl)
300c-
301          call retmem(3)
302          return
303       endif
304cc
305c- 4 l-shell
306       if (lshellt.eq.4) then
307c-
308          jbuf  = nbuf
309          jbuf2 =nbuf2
310          jbfij1=ibfij1
311          jbfij2=ibfij2
312          jbfkl1=ibfkl1
313          jbfkl2=ibfkl2
314c-
315          jbfij3=ibfij3
316          jbfkl3=ibfkl3
317          jbf2l1=ibf2l1
318          jbf2l2=ibf2l2
319          jbf2l3=ibf2l3
320          jbf2l4=ibf2l4
321c
322          jbf3l1=ibf3l1
323          jbf3l2=ibf3l2
324          jbf3l3=ibf3l3
325          jbf3l4=ibf3l4
326c-
327          jssss =issss
328c---
329          call shift4l(bl(jbuf),bl(jbuf2),l01,l02,
330     *                 bl(jbfij1),bl(jbfij2),bl(jbfkl1),bl(jbfkl2),
331     *                 lqij,lqkl,
332     *                 bl(jbfij3),bl(jbfkl3),
333     *                 bl(jbf2l1),bl(jbf2l2),bl(jbf2l3),bl(jbf2l4),
334     *                 bl(jbf3l1),bl(jbf3l2),bl(jbf3l3),bl(jbf3l4),lqmx,
335     *                 bl(jssss),
336     *                 bl(iwij),lsmx, bl(ivij),bl(iuij),bl(isij),lsjl,
337     *           bl(ixij),nfu(nqi+1),lnkl,bl(iyij),bl(izij),nfu(nqj+1),
338     *                 bl(ix2l1),bl(ix2l2),bl(ix2l3),bl(ix2l4),
339     *                 bl(ix3l1),bl(ix3l2),bl(ix3l3),bl(ix3l4),mnbls,
340     *                 bl(ixabn),bl(ixcdn),
341     *                 bl(indxx),lni,lnj,lnk,lnl)
342c-
343          call retmem(3)
344          return
345       endif
346c----------------------
347      return
348      end
349c=======================================================
350c moved into the convert.f file
351c     subroutine convr3(bl,m,nbls,npij,npkl,idx1,idx2,
352c    *                   xab,xcd, ixabn,ixcdn)
353c=======================================================
354      subroutine daxpy3(n,a,z1,z2,z3,y1,y2,y3,x)
355c------------------------------------------------
356c* performs the vector operations with a stride=1
357c*
358c*     Z = Y + A*X
359c*
360c*   for three values of a and three matrices Z and Y and this same X
361c*
362c------------------------------------------------
363      implicit real*8 (a-h,o-z)
364      dimension z1(n),z2(n),z3(n),y1(n),y2(n),y3(n),x(n),a(n,3)
365c
366      do 10 i=1,n
367      z1(i)=y1(i) + a(i,1)*x(i)
368      z2(i)=y2(i) + a(i,2)*x(i)
369      z3(i)=y3(i) + a(i,3)*x(i)
370   10 continue
371c*
372      end
373c=====================================================
374      subroutine noshift(lt1,lt2,mnbls,
375     *                   buf,buf2,
376     *                   bfijx,      bfklx,
377     *                   bf2lx,
378     *                   lt3,lt4,
379     *                   indxx,ln12,ln34)
380c**
381c**   special cases where the shifting of angular
382c**   momentum is not needed at all / (ds|ss)..(ss|sd),
383c**   (xs|ys),(xs|sy),(sx|ys),(sx|sy) /
384c**
385      implicit real*8 (a-h,o-z)
386      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
387cc
388      common/obarai/
389     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
390     * nqi,nqj,nqk,nql,nsij,nskl,
391     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg
392      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
393cc
394      dimension buf2(mnbls,lt1,lt2),
395c    * bfij1(mnbls,lt3,lt2),bfij2(mnbls,lt3,lt2),
396c    * bfkl1(mnbls,lt1,lt4),bfkl2(mnbls,lt1,lt4),
397c    * bf2l1(mnbls,lt3,lt4),bf2l2(mnbls,lt3,lt4),
398c    * bf2l3(mnbls,lt3,lt4),bf2l4(mnbls,lt3,lt4)
399c---
400     * bfijx(mnbls,lt3,lt2),
401     * bfklx(mnbls,lt1,lt4),
402     * bf2lx(mnbls,lt3,lt4)
403       dimension buf(mnbls,*)
404       dimension indxx(ln12,ln34)
405c-----------------------------------
406c
407       ijb1=ijbeg-1
408       klb1=klbeg-1
409c2002
410c        ijkl=0
411c        do 5031 i=1,lni
412c        ii=(i-1)*lnj
413c        do 5031 j=1,lnj
414c        ij=ii+j
415c        do 5031 k=1,lnk
416c        kk=(k-1)*lnl
417c        do 5031 l=1,lnl
418c        kl=kk+l
419c        ijkl=ijkl+1
420c        indxx(ij     ,kl     )=ijkl
421c5031    continue
422c2002
423         do 5034 ij=ijbeg,lnij
424         do 5034 kl=klbeg,lnkl
425         ijkl=indxx(ij-ijb1,kl-klb1)
426            do 5034 i=1,mnbls
427         buf(i,ijkl)=buf2(i,ij,kl)
428 5034    continue
429c
430           if(lshelkl.eq.1 .or. lshelkl.eq.2) then
431              do 5035 ij=ijbeg,lnij
432              ijkl=indxx(ij-ijb1,1)
433                 do 5035 i=1,mnbls
434c--->         buf(i,ijkl)=bfkl1(i,ij,1)
435              buf(i,ijkl)=bfklx(i,ij,1)
436 5035         continue
437           endif
438c-------
439           if(lshelij.eq.1 .or. lshelij.eq.2) then
440              do 6035 kl=klbeg,lnkl
441              ijkl=indxx(1,kl-klb1)
442                 do 6035 i=1,mnbls
443c----->       buf(i,ijkl)=bfij1(i,1,kl)
444              buf(i,ijkl)=bfijx(i,1,kl)
445 6035         continue
446           endif
447c---------
448           if(lshellt.eq.2) then
449                ijkl=indxx(1,1)
450                   do 1010 i=1,mnbls
451c------>           buf(i,ijkl)=bf2l1(i,1,1)
452                   buf(i,ijkl)=bf2lx(i,1,1)
453 1010              continue
454cc
455           endif
456      return
457      end
458c=====================================================
459      subroutine shift0l(buf,buf2,lt1,lt2,
460     *                   wij,lt3,lsjl,xij,lt4,lt5,lt6,mnbls,xab,xcd,
461     *                   indxx,lni1,lnj1,lnk1,lnl1)
462c------------------------------------
463c  when l-shells are not present
464c------------------------------------
465      implicit real*8 (a-h,o-z)
466#include "texas_lpar.fh"
467      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
468      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
469      common/obarai/
470     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
471     * nqi,nqj,nqk,nql,nsij,nskl,
472     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbex,klbex
473c
474       dimension buf2(mnbls,lt1,lt2),wij(mnbls,lt3,lsjl)
475       dimension xij(mnbls,lt4,lt5,lt6)
476       dimension buf(mnbls,*)
477       dimension xab(mnbls,3),xcd(mnbls,3)
478       dimension indxx(lni1,lnj1,lnk1,lnl1)
479c------------------------------------
480c
481       nibeg=nfu(nqi)+1
482       niend=nfu(nqi+1)
483       njbeg=nfu(nqj)+1
484       njend=nfu(nqj+1)
485ccc
486       nkbeg=nfu(nqk)+1
487       nkend=nfu(nqk+1)
488       nlbeg=nfu(nql)+1
489       nlend=nfu(nql+1)
490c
491       nqix=nqi
492ccccc  nqjx=nqj
493c
494       nqkx=nqk
495ccccc  nqlx=nql
496       nqklx=nqkl
497c
498      do 100 nkl=nfu(nqklx)+1,nfu(nskl+1)
499c
500           do 102 nij=nfu(nqix)+1,nfu(nsij+1)
501           call amtfer(buf2(1,nij,nkl),wij(1,nij,1),mnbls)
502  102      continue
503c
504       call horiz12(wij,lt3,lsjl,xab,mnbls,nqi,nqj,nsij1)
505c
506       do 107 nj=njbeg,njend
507       do 107 ni=nibeg,niend
508          call amtfer(wij(1,ni,nj),xij(1,ni,nj,nkl),mnbls)
509  107  continue
510  100 continue
511c
512c------------------------------------
513c this part shifts angular momentum
514c   from position 3 to position 4
515c------------------------------------
516c
517      ixyz=0
518      do 300 ni=nibeg,niend
519      ixyz=ixyz+1
520      jxyz=0
521      do 300 nj=njbeg,njend
522      jxyz=jxyz+1
523c
524         do 301 nkl=nfu(nqkx)+1,nfu(nskl+1)
525         call amtfer(xij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
526  301    continue
527c
528       call horiz12(wij,lt3,lsjl,xcd,mnbls,nqk,nql,nskl1)
529c
530      kxyz=0
531      do 305 nk=nkbeg,nkend
532      kxyz=kxyz+1
533      lxyz=0
534      do 305 nl=nlbeg,nlend
535      lxyz=lxyz+1
536      indx=indxx(ixyz,jxyz,kxyz,lxyz)
537      call amtfer(wij(1,nk,nl),buf(1,indx),mnbls)
538  305 continue
539  300 continue
540c
541      end
542c=============================================================
543      subroutine shift1l(buf,
544     *                buf2,lt1,lt2,bijx,bklx,lt3,lt4,
545     *                wij,lt5,vij,lt6,
546     *                xij,lt7,lt8,yij,lt9,mnbls,xab,xcd,
547     *                indxx,lni1,lnj1,lnk1,lnl1)
548c***********************************************************
549c*
550c*  when 1 l-shell is present somewhere
551c*          in position 1,2 or 3,4
552c***********************************************************
553      implicit real*8 (a-h,o-z)
554#include "texas_lpar.fh"
555cxxx
556      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
557c
558      common/obarai/
559     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
560     * nqi,nqj,nqk,nql,nsij,nskl,
561     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbex,klbex
562c
563      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
564c
565c
566      dimension buf2(mnbls,lt1,lt2),
567     * bijx(mnbls,lt3,lt2),
568     * bklx(mnbls,lt1,lt4)
569      dimension wij(mnbls,lt5,lt6),vij(mnbls,lt5,lt6)
570      dimension xij(mnbls,lt7,lt3,lt8),yij(mnbls,lt7,lt9,lt4)
571      dimension buf(mnbls,*)
572      dimension xab(mnbls,3),xcd(mnbls,3)
573c
574      dimension indxx(lni1,lnj1,lnk1,lnl1)
575c***********************************************************
576c
577       nibeg=nfu(nqi)+1
578       niend=nfu(nqi+1)
579       njbeg=nfu(nqj)+1
580       njend=nfu(nqj+1)
581c
582       nkbeg=nfu(nqk)+1
583       nkend=nfu(nqk+1)
584       nlbeg=nfu(nql)+1
585       nlend=nfu(nql+1)
586c
587       nqix=nqi
588ccccc  nqjx=nqj
589       if(ityp.eq.3) then
590          nibeg=1
591          if(jtyp.le.3) nqix=1
592       endif
593       if(jtyp.eq.3) then
594          njbeg=1
595ccccc     if(ityp.le.3) nqjx=1
596       endif
597c
598       nqkx=nqk
599ccccc  nqlx=nql
600       nqklx=nqkl
601       if(ktyp.eq.3) then
602          nkbeg=1
603          if(ltyp.le.3) then
604             nqklx=1
605             nqkx=1
606          endif
607       endif
608       if(ltyp.eq.3) then
609          nlbeg=1
610          if(ktyp.le.3) then
611             nqklx=1
612ccccc        nqlx=1
613          endif
614       endif
615c
616c*****
617c
618      do 10 nkl=nfu(nqkl)+1,nfu(nskl+1)
619c
620       do 11 ij=nqi,nsij
621       ijbeg=nfu(ij)+1
622       ijend=nfu(ij+1)
623           do 12 nij=ijbeg,ijend
624           call amtfer(buf2(1,nij,nkl),wij(1,nij,1),mnbls)
625   12      continue
626c
627   11  continue
628c
629cccccccc
630       call horiz12(wij,lt5,lt6,xab,mnbls,nqi,nqj,nsij1)
631cccccccc
632c
633       do 16 nj=nfu(nqj)+1,nfu(nqj+1)
634       do 16 ni=nfu(nqi)+1,nfu(nqi+1)
635       call amtfer(wij(1,ni,nj),xij(1,ni,nj,nkl),mnbls)
636   16  continue
637c
638ccc  here lshelij can be eq. 0, 1 or 2 only  ccc
639c
640      if(lshelij.gt.0) then
641          if(lshelij.eq.1) then
642              if(jtyp.eq.1) then
643                call amtfer(bijx(1,1,nkl),xij(1,1,1,nkl),mnbls)
644              else
645                call daxpy3(mnbls,xab,
646     *          xij(1,1,2,nkl),xij(1,1,3,nkl),xij(1,1,4,nkl),
647     *          bijx(1,2,nkl),bijx(1,3,nkl),bijx(1,4,nkl),bijx(1,1,nkl))
648              endif
649          else
650               do 17 nij=nfu(nqi )+1,nfu(nqi +1)
651               call amtfer(bijx(1,nij,nkl),xij(1,nij,1,nkl),mnbls)
652   17          continue
653          endif
654      endif
655c
656   10 continue
657c******
658ccc  here lshelkl can be eq. 0, 1 or 2 only  ccc
659      if(lshelkl.gt.0) then
660ccc
661      do 100 nkl=nfu(nqklx)+1,nfu(nqkl+1)
662c
663       do 101 ij=nqix,nsij
664       ijbeg=nfu(ij)+1
665       ijend=nfu(ij+1)
666c
667             do 103 nij=ijbeg,ijend
668             call amtfer(bklx(1,nij,nkl),vij(1,nij,1),mnbls)
669  103        continue
670  101  continue
671c
672cccccc
673        call horiz12(vij,lt5,lt6,xab,mnbls,nqi,nqj,nsij1)
674cccccc
675            do 1071 nj=njbeg,njend
676            do 1071 ni=nibeg,niend
677            call amtfer(vij(1,ni,nj),yij(1,ni,nj,nkl),mnbls)
678 1071       continue
679c
680  100 continue
681c
682      endif
683c
684c***********************************************************
685c*    this part    shifts the angular momentum
686c*         from position 3 to position 4
687c***********************************************************
688c
689      ixyz=0
690      do 300 ni=nibeg,niend
691      ixyz=ixyz+1
692      jxyz=0
693      do 300 nj=njbeg,njend
694      jxyz=jxyz+1
695c
696         do 301 nkl=nfu(nqkx)+1,nfu(nskl+1)
697         call amtfer(xij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
698  301    continue
699c
700        call horiz12(wij,lt5,lt6,xcd,mnbls,nqk,nql,nskl1)
701c
702      if(lshelkl.gt.0)  then
703c
704         if(lshelkl.eq.1) then
705            if(ltyp.eq.1) then
706             call amtfer(yij(1,ni,nj,1),wij(1,1,1),mnbls)
707            else
708      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
709     *    yij(1,ni,nj,2),yij(1,ni,nj,3),yij(1,ni,nj,4),yij(1,ni,nj,1))
710            endif
711         else
712             do 312 nkl=nfu(nqkx)+1,nfu(nqkl+1)
713             call amtfer(yij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
714  312        continue
715         endif
716      endif
717c
718c
719ccccccccccccccccccccccccccccccccccc
720      kxyz=0
721      do 305 nk=nkbeg,nkend
722      kxyz=kxyz+1
723      lxyz=0
724      do 305 nl=nlbeg,nlend
725      lxyz=lxyz+1
726      indx=indxx(ixyz,jxyz,kxyz,lxyz)
727      call amtfer(wij(1,nk,nl),buf(1,indx),mnbls)
728  305 continue
729  300 continue
730c****
731c
732      return
733      end
734c=============================================================
735      subroutine shift2l(buf,
736     *                buf2,lt1,lt2,bij1,bij2,bkl1,bkl2,lt3,lt4,
737     *                bij3,bkl3,b2l12,b2l34,
738     *                wij,lt5,vij,uij,sij,lt6,
739     *                xij,lt7,lt8,yij,zij,lt9,
740     *                x2l,mnbls,xab,xcd,
741     *                indxx,lni1,lnj1,lnk1,lnl1)
742c***********************************************************
743c*
744c*  when 2 l-shells are present somewhere
745c*          in position 1,2 or 3,4
746c***********************************************************
747      implicit real*8 (a-h,o-z)
748#include "texas_lpar.fh"
749cxxxx
750      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
751c
752      common/obarai/
753     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
754     * nqi,nqj,nqk,nql,nsij,nskl,
755     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbex,klbex
756c
757      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
758c
759c
760      dimension x2l(mnbls,lt3,lt4)
761      dimension buf2(mnbls,lt1,lt2),
762     * bij1(mnbls,lt3,lt2),bij2(mnbls,lt3,lt2),
763     * bkl1(mnbls,lt1,lt4),bkl2(mnbls,lt1,lt4),
764     * bij3(mnbls,lt2),bkl3(mnbls,lt1),
765     * b2l12(mnbls,lt3,lt4),
766     * b2l34(mnbls,lt3,lt4)
767       dimension wij(mnbls,lt5,lt6),
768     * vij(mnbls,lt5,lt6),uij(mnbls,lt5,lt6),sij(mnbls,lt5,lt6)
769      dimension xij(mnbls,lt7,lt3,lt8),yij(mnbls,lt7,lt9,lt4),
770     *                             zij(mnbls,lt7,lt9,lt4)
771      dimension buf(mnbls,*)
772      dimension xab(mnbls,3),xcd(mnbls,3)
773      dimension indxx(lni1,lnj1,lnk1,lnl1)
774c***********************************************************
775c
776       nibeg=nfu(nqi)+1
777       niend=nfu(nqi+1)
778       njbeg=nfu(nqj)+1
779       njend=nfu(nqj+1)
780c
781       nkbeg=nfu(nqk)+1
782       nkend=nfu(nqk+1)
783       nlbeg=nfu(nql)+1
784       nlend=nfu(nql+1)
785c
786       nqix=nqi
787ccccc  nqjx=nqj
788       if(ityp.eq.3) then
789          nibeg=1
790          if(jtyp.le.3) nqix=1
791       endif
792       if(jtyp.eq.3) then
793          njbeg=1
794ccccc     if(ityp.le.3) nqjx=1
795       endif
796c
797       nqkx=nqk
798ccccc  nqlx=nql
799       nqklx=nqkl
800       if(ktyp.eq.3) then
801          nkbeg=1
802          if(ltyp.le.3) then
803             nqklx=1
804             nqkx=1
805          endif
806       endif
807       if(ltyp.eq.3) then
808          nlbeg=1
809          if(ktyp.le.3) then
810             nqklx=1
811ccccc        nqlx=1
812          endif
813       endif
814c
815c*****
816      do 100 nkl=nfu(nqklx)+1,nfu(nskl+1)
817c
818       do 101 ij=nqix,nsij
819       ijbeg=nfu(ij)+1
820       ijend=nfu(ij+1)
821           do 102 nij=ijbeg,ijend
822           call amtfer(buf2(1,nij,nkl),wij(1,nij,1),mnbls)
823  102      continue
824c
825       if( nkl.le.nfu(nqkl+1)) then
826         if(lshelkl.eq.1.or.lshelkl.eq.3) then
827           do 103 nij=ijbeg,ijend
828           call amtfer(bkl1(1,nij,nkl),vij(1,nij,1),mnbls)
829  103      continue
830         endif
831         if(lshelkl.eq.2.or.lshelkl.eq.3) then
832           do 104 nij=ijbeg,ijend
833           call amtfer(bkl2(1,nij,nkl),uij(1,nij,1),mnbls)
834  104      continue
835         endif
836         if(lshelkl.eq.3.and.nkl.eq.1) then
837           do 105 nij=ijbeg,ijend
838           call amtfer(bkl3(1,nij),sij(1,nij,1),mnbls)
839  105      continue
840         endif
841       endif
842c
843  101  continue
844c
845ccccccccccccccccccccccc
846c
847c
848       call horiz12(wij,lt5,lt6,xab,mnbls,nqi,nqj,nsij1)
849c
850      if( nkl.le.nfu(nqkl+1)) then
851          if(lshelkl.eq.1.or.lshelkl.eq.3) then
852            call horiz12(vij,lt5,lt6,xab,mnbls,nqi,nqj,nsij1)
853          endif
854          if(lshelkl.eq.2.or.lshelkl.eq.3) then
855            call horiz12(uij,lt5,lt6,xab,mnbls,nqi,nqj,nsij1)
856          endif
857      endif
858      if(nkl.eq.1) then
859          if(lshelkl.eq.3) then
860            call horiz12(sij,lt5,lt6,xab,mnbls,nqi,nqj,nsij1)
861          endif
862      endif
863cccc
864c
865       do 107 nj=njbeg,njend
866       do 107 ni=nibeg,niend
867       call amtfer(wij(1,ni,nj),xij(1,ni,nj,nkl),mnbls)
868  107  continue
869       if( nkl.le.nfu(nqkl+1)) then
870           if(lshelkl.eq.1.or.lshelkl.eq.3) then
871                do 1071 nj=njbeg,njend
872                do 1071 ni=nibeg,niend
873                call amtfer(vij(1,ni,nj),yij(1,ni,nj,nkl),mnbls)
874 1071           continue
875           endif
876           if(lshelkl.eq.2.or.lshelkl.eq.3) then
877                do 1072 nj=njbeg,njend
878                do 1072 ni=nibeg,niend
879                call amtfer(uij(1,ni,nj),zij(1,ni,nj,nkl),mnbls)
880 1072           continue
881           endif
882       endif
883c
884       if(lshelij.eq.1.or.lshelij.eq.3) then
885          if(jtyp.eq.1) then
886             call amtfer(bij1(1,1,nkl),xij(1,1,1,nkl),mnbls)
887          else
888      call daxpy3(mnbls,xab,
889     *          xij(1,1,2,nkl),xij(1,1,3,nkl),xij(1,1,4,nkl),
890     *          bij1(1,2,nkl),bij1(1,3,nkl),bij1(1,4,nkl),bij1(1,1,nkl))
891          endif
892       endif
893       if(lshelij.eq.2.or.lshelij.eq.3) then
894           do 108 nij=nfu(nqi )+1,nfu(nqi +1)
895           call amtfer(bij2(1,nij,nkl),xij(1,nij,1,nkl),mnbls)
896  108      continue
897       endif
898       if(lshelij.eq.3) then
899           call amtfer(bij3(1,nkl),xij(1,1,1,nkl),mnbls)
900       endif
901c*****
902      if( nkl.le.nfu(nqkl+1)) then
903c****  2 l-shells ****
904c
905           if(lshelij.eq.1) then
906              if(jtyp.eq.1) then
907                if(lcas2(1).eq.1 .or. lcas2(2).eq.1) then
908                   call amtfer(b2l12(1,1,nkl),x2l(1,1,nkl),mnbls)
909                endif
910              else
911                if(lcas2(1).eq.1 .or. lcas2(2).eq.1) then
912                   call daxpy3(mnbls,xab,
913     *             x2l(1,2,nkl),x2l(1,3,nkl),x2l(1,4,nkl),
914     *             b2l12(1,2,nkl),b2l12(1,3,nkl),b2l12(1,4,nkl),
915     *             b2l12(1,1,nkl))
916                endif
917              endif
918           endif
919           if(lshelij.eq.2) then
920                if(lcas2(3).eq.1 .or. lcas2(4).eq.1) then
921                  do 109 nij=nfu(nqix)+1,nfu(nqij+1)
922                  call amtfer(b2l34(1,nij,nkl),x2l(1,nij,nkl),mnbls)
923  109             continue
924                endif
925           endif
926c
927      endif
928c
929  100 continue
930c
931c***********************************************************
932c*    this part    shifts the angular momentum
933c*         from position 3 to position 4
934c***********************************************************
935c
936      ixyz=0
937      do 300 ni=nibeg,niend
938      ixyz=ixyz+1
939      jxyz=0
940      do 300 nj=njbeg,njend
941      jxyz=jxyz+1
942c
943         do 301 nkl=nfu(nqkx)+1,nfu(nskl+1)
944         call amtfer(xij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
945  301    continue
946c
947cccccccccc
948c
949       call horiz12(wij,lt5,lt6,xcd,mnbls,nqk,nql,nskl1)
950cccccccccc
951c
952       if(lshelkl.eq.1.or.lshelkl.eq.3) then
953         if(ltyp.eq.1) then
954           call amtfer(yij(1,ni,nj,1),wij(1,1,1),mnbls)
955         else
956      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
957     * yij(1,ni,nj,2),yij(1,ni,nj,3),yij(1,ni,nj,4),yij(1,ni,nj,1))
958         endif
959       endif
960       if(lshelkl.eq.2.or.lshelkl.eq.3) then
961         do 312 nkl=nfu(nqkx)+1,nfu(nqkl+1)
962         call amtfer(zij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
963  312    continue
964       endif
965c
966       if(lshelkl.eq.3) then
967         call amtfer(sij(1,ni,nj),wij(1,1,1),mnbls)
968       endif
969c
970       if( ni.le.nfu(nqi +1).and.nj.le.nfu(nqj +1) ) then
971c
972c****  2 l-shells ****
973         if(lshelkl.eq.1) then
974            if(ltyp.eq.1) then
975              if(lcas2(1).eq.1) then
976               if(ni.eq.1.and.nj.ge.nfu(nqj)+1) then
977                   call amtfer(x2l(1,nj,1),wij(1,1,1),mnbls)
978               endif
979              endif
980              if(lcas2(3).eq.1) then
981               if(nj.eq.1.and.ni.ge.nfu(nqi)+1) then
982                   call amtfer(x2l(1,ni,1),wij(1,1,1),mnbls)
983               endif
984              endif
985            else
986              if(lcas2(1).eq.1.and.ni.eq.1.and.nj.ge.nfu(nqj)+1) then
987      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
988     *            x2l(1,nj,2),x2l(1,nj,3),x2l(1,nj,4),x2l(1,nj,1))
989              endif
990              if(lcas2(3).eq.1.and.nj.eq.1.and.ni.ge.nfu(nqi)+1) then
991      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
992     *            x2l(1,ni,2),x2l(1,ni,3),x2l(1,ni,4),x2l(1,ni,1))
993              endif
994            endif
995         endif
996         if(lshelkl.eq.2) then
997              do 313 nkl=nfu(nqkx)+1,nfu(nqkl+1)
998              if(lcas2(2).eq.1.and.ni.eq.1) then
999                  call amtfer(x2l(1,nj,nkl),wij(1,nkl,1),mnbls)
1000              endif
1001              if(lcas2(4).eq.1.and.nj.eq.1) then
1002                  call amtfer(x2l(1,ni,nkl),wij(1,nkl,1),mnbls)
1003              endif
1004  313         continue
1005         endif
1006       endif
1007c****
1008c
1009ccccccccccccccccccccccccccccccccccc
1010      kxyz=0
1011      do 305 nk=nkbeg,nkend
1012      kxyz=kxyz+1
1013      lxyz=0
1014      do 305 nl=nlbeg,nlend
1015      lxyz=lxyz+1
1016      indx=indxx(ixyz,jxyz,kxyz,lxyz)
1017      call amtfer(wij(1,nk,nl),buf(1,indx),mnbls)
1018  305 continue
1019  300 continue
1020c
1021      return
1022      end
1023c=============================================================
1024      subroutine shift3l(buf,
1025     *                buf2,lt1,lt2,bij1,bij2,bkl1,bkl2,lt3,lt4,
1026     *                bij3,bkl3,b2l1,b2l2,b2l3,b2l4,
1027     *                b3l12,b3l34,lt5,
1028     *                wij,lt6,vij,uij,sij,lt7,
1029     *                xij,lt8,lt9,yij,zij,lt10,
1030     *                x2l1,x2l2,x2l3,x2l4,x3l12,x3l34,mnbls,
1031     *                xab,xcd,
1032     *                indxx,lni1,lnj1,lnk1,lnl1)
1033c***********************************************************
1034c*
1035c*  when 3 l-shells are present somewhere
1036c*          in positions 1,2 or 3,4
1037c***********************************************************
1038      implicit real*8 (a-h,o-z)
1039#include "texas_lpar.fh"
1040cxxx
1041      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
1042c
1043      common/obarai/
1044     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
1045     * nqi,nqj,nqk,nql,nsij,nskl,
1046     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbex,klbex
1047c
1048      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
1049c
1050c
1051      dimension x2l1(mnbls,lt3,lt4),x2l2(mnbls,lt3,lt4),
1052     *          x2l3(mnbls,lt3,lt4),x2l4(mnbls,lt3,lt4)
1053      dimension x3l12(mnbls,lt4),x3l34(mnbls,lt3)
1054c
1055      dimension buf2(mnbls,lt1,lt2),
1056     * bij1(mnbls,lt3,lt2),bij2(mnbls,lt3,lt2),
1057     * bkl1(mnbls,lt1,lt4),bkl2(mnbls,lt1,lt4),
1058     * bij3(mnbls,lt2),bkl3(mnbls,lt1),
1059     * b2l1(mnbls,lt3,lt4),b2l2(mnbls,lt3,lt4),
1060     * b2l3(mnbls,lt3,lt4),b2l4(mnbls,lt3,lt4),
1061     * b3l12(mnbls,lt5),b3l34(mnbls,lt5)
1062      dimension wij(mnbls,lt6,lt7),
1063     * vij(mnbls,lt6,lt7),uij(mnbls,lt6,lt7),sij(mnbls,lt6,lt7)
1064      dimension xij(mnbls,lt8,lt3,lt9),yij(mnbls,lt8,lt10,lt4)
1065     *                            ,zij(mnbls,lt8,lt10,lt4)
1066      dimension buf(mnbls,*)
1067      dimension xab(mnbls,3),xcd(mnbls,3)
1068      dimension indxx(lni1,lnj1,lnk1,lnl1)
1069c***********************************************************
1070c
1071       nibeg=nfu(nqi)+1
1072       niend=nfu(nqi+1)
1073       njbeg=nfu(nqj)+1
1074       njend=nfu(nqj+1)
1075c
1076       nkbeg=nfu(nqk)+1
1077       nkend=nfu(nqk+1)
1078       nlbeg=nfu(nql)+1
1079       nlend=nfu(nql+1)
1080c
1081       nqix=nqi
1082ccccc  nqjx=nqj
1083       if(ityp.eq.3) then
1084          nibeg=1
1085          if(jtyp.le.3) nqix=1
1086       endif
1087       if(jtyp.eq.3) then
1088          njbeg=1
1089ccccc     if(ityp.le.3) nqjx=1
1090       endif
1091c
1092       nqkx=nqk
1093ccccc  nqlx=nql
1094       nqklx=nqkl
1095       if(ktyp.eq.3) then
1096          nkbeg=1
1097          if(ltyp.le.3) then
1098             nqklx=1
1099             nqkx=1
1100          endif
1101       endif
1102       if(ltyp.eq.3) then
1103          nlbeg=1
1104          if(ktyp.le.3) then
1105             nqklx=1
1106ccccc        nqlx=1
1107          endif
1108       endif
1109c
1110c*****
1111c
1112      do 100 nkl=nfu(nqklx)+1,nfu(nskl+1)
1113c
1114       do 101 ij=nqix,nsij
1115       ijbeg=nfu(ij)+1
1116       ijend=nfu(ij+1)
1117           do 102 nij=ijbeg,ijend
1118           call amtfer(buf2(1,nij,nkl),wij(1,nij,1),mnbls)
1119  102      continue
1120c
1121       if( nkl.le.nfu(nqkl+1)) then
1122         if(lshelkl.eq.1.or.lshelkl.eq.3) then
1123           do 103 nij=ijbeg,ijend
1124           call amtfer(bkl1(1,nij,nkl),vij(1,nij,1),mnbls)
1125  103      continue
1126         endif
1127         if(lshelkl.eq.2.or.lshelkl.eq.3) then
1128           do 104 nij=ijbeg,ijend
1129           call amtfer(bkl2(1,nij,nkl),uij(1,nij,1),mnbls)
1130  104      continue
1131         endif
1132         if(lshelkl.eq.3.and.nkl.eq.1) then
1133           do 105 nij=ijbeg,ijend
1134           call amtfer(bkl3(1,nij),sij(1,nij,1),mnbls)
1135  105      continue
1136         endif
1137       endif
1138c
1139  101  continue
1140c
1141cccccccccccc
1142c
1143       call horiz12(wij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1144c
1145      if( nkl.le.nfu(nqkl+1)) then
1146          if(lshelkl.eq.1.or.lshelkl.eq.3) then
1147            call horiz12(vij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1148          endif
1149          if(lshelkl.eq.2.or.lshelkl.eq.3) then
1150            call horiz12(uij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1151          endif
1152      endif
1153      if(nkl.eq.1) then
1154          if(lshelkl.eq.3) then
1155            call horiz12(sij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1156          endif
1157      endif
1158cccccccccccc
1159c
1160       do 107 nj=njbeg,njend
1161       do 107 ni=nibeg,niend
1162       call amtfer(wij(1,ni,nj),xij(1,ni,nj,nkl),mnbls)
1163  107  continue
1164       if( nkl.le.nfu(nqkl+1)) then
1165           if(lshelkl.eq.1.or.lshelkl.eq.3) then
1166                do 1071 nj=njbeg,njend
1167                do 1071 ni=nibeg,niend
1168                call amtfer(vij(1,ni,nj),yij(1,ni,nj,nkl),mnbls)
1169 1071           continue
1170           endif
1171           if(lshelkl.eq.2.or.lshelkl.eq.3) then
1172                do 1072 nj=njbeg,njend
1173                do 1072 ni=nibeg,niend
1174                call amtfer(uij(1,ni,nj),zij(1,ni,nj,nkl),mnbls)
1175 1072           continue
1176           endif
1177       endif
1178c
1179       if(lshelij.eq.1.or.lshelij.eq.3) then
1180          if(jtyp.eq.1) then
1181             call amtfer(bij1(1,1,nkl),xij(1,1,1,nkl),mnbls)
1182          else
1183      call daxpy3(mnbls,xab,
1184     *    xij(1,1,2,nkl),xij(1,1,3,nkl),xij(1,1,4,nkl),
1185     *    bij1(1,2,nkl),bij1(1,3,nkl),bij1(1,4,nkl),bij1(1,1,nkl))
1186          endif
1187       endif
1188       if(lshelij.eq.2.or.lshelij.eq.3) then
1189           do 108 nij=nfu(nqi )+1,nfu(nqi +1)
1190           call amtfer(bij2(1,nij,nkl),xij(1,nij,1,nkl),mnbls)
1191  108      continue
1192       endif
1193       if(lshelij.eq.3) then
1194           call amtfer(bij3(1,nkl),xij(1,1,1,nkl),mnbls)
1195       endif
1196c*****
1197      if( nkl.le.nfu(nqkl+1)) then
1198c****  2 or 3 l-shells ****
1199c
1200           if(lshelij.eq.1) then
1201              if(jtyp.eq.1) then
1202                   call amtfer(b2l1(1,1,nkl),x2l1(1,1,nkl),mnbls)
1203                   call amtfer(b2l2(1,1,nkl),x2l2(1,1,nkl),mnbls)
1204c-- test ?      if(nkl.eq.1) then
1205                if(nkl.eq.1 .and. lcas3(3).eq.1) then
1206                   call amtfer(b3l34(1,1),x3l34(1,1),mnbls)
1207                endif
1208              else
1209                call daxpy3(mnbls,xab,
1210     *          x2l1(1,2,nkl),x2l1(1,3,nkl),x2l1(1,4,nkl),
1211     *          b2l1(1,2,nkl),b2l1(1,3,nkl),b2l1(1,4,nkl),b2l1(1,1,nkl))
1212cc
1213                call daxpy3(mnbls,xab,
1214     *          x2l2(1,2,nkl),x2l2(1,3,nkl),x2l2(1,4,nkl),
1215     *          b2l2(1,2,nkl),b2l2(1,3,nkl),b2l2(1,4,nkl),b2l2(1,1,nkl))
1216c-- test ?      if(nkl.eq.1) then
1217                if(nkl.eq.1 .and. lcas3(3).eq.1) then
1218                  call daxpy3(mnbls,xab,
1219     *            x3l34(1,2),x3l34(1,3),x3l34(1,4),
1220     *            b3l34(1,2),b3l34(1,3),b3l34(1,4),b3l34(1,1))
1221                endif
1222              endif
1223           endif
1224           if(lshelij.eq.2) then
1225                do 109 nij=nfu(nqix)+1,nfu(nqij+1)
1226                   call amtfer(b2l3(1,nij,nkl),x2l3(1,nij,nkl),mnbls)
1227                   call amtfer(b2l4(1,nij,nkl),x2l4(1,nij,nkl),mnbls)
1228  109           continue
1229c-- test ?      if(nkl.eq.1) then
1230                if(nkl.eq.1 .and. lcas3(4).eq.1) then
1231                   do 110 nij=nfu(nqi )+1,nfu(nqi +1)
1232                   call amtfer(b3l34(1,nij),x3l34(1,nij),mnbls)
1233  110              continue
1234                endif
1235           endif
1236           if(lshelij.eq.3) then
1237                if(lcas2(1).eq.1) then
1238                  call daxpy3(mnbls,xab,
1239     *            x2l1(1,2,nkl),x2l1(1,3,nkl),x2l1(1,4,nkl),
1240     *            b2l1(1,2,nkl),b2l1(1,3,nkl),b2l1(1,4,nkl),
1241     *            b2l1(1,1,nkl))
1242cc
1243                  call amtfer(b2l3(1,2,nkl),x2l3(1,2,nkl),mnbls)
1244                  call amtfer(b2l3(1,3,nkl),x2l3(1,3,nkl),mnbls)
1245                  call amtfer(b2l3(1,4,nkl),x2l3(1,4,nkl),mnbls)
1246                endif
1247                if(lcas2(2).eq.1) then
1248                  call daxpy3(mnbls,xab,
1249     *            x2l2(1,2,nkl),x2l2(1,3,nkl),x2l2(1,4,nkl),
1250     *            b2l2(1,2,nkl),b2l2(1,3,nkl),b2l2(1,4,nkl),
1251     *            b2l2(1,1,nkl))
1252cc
1253                  call amtfer(b2l4(1,2,nkl),x2l4(1,2,nkl),mnbls)
1254                  call amtfer(b2l4(1,3,nkl),x2l4(1,3,nkl),mnbls)
1255                  call amtfer(b2l4(1,4,nkl),x2l4(1,4,nkl),mnbls)
1256                endif
1257                if(lcas3(1).eq.1) then
1258                  call amtfer(b3l12(1,nkl),x3l12(1,nkl),mnbls)
1259                endif
1260                if(lcas3(2).eq.1) then
1261                  call amtfer(b3l12(1,nkl),x3l12(1,nkl),mnbls)
1262                endif
1263           endif
1264      endif
1265c
1266  100 continue
1267c
1268c***********************************************************
1269c*    this part    shifts the angular momentum
1270c*         from position 3 to position 4
1271c***********************************************************
1272c
1273      ixyz=0
1274      do 300 ni=nibeg,niend
1275      ixyz=ixyz+1
1276      jxyz=0
1277      do 300 nj=njbeg,njend
1278      jxyz=jxyz+1
1279c
1280         do 301 nkl=nfu(nqkx)+1,nfu(nskl+1)
1281         call amtfer(xij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
1282  301    continue
1283c
1284ccccccc
1285        call horiz12(wij,lt6,lt7,xcd,mnbls,nqk,nql,nskl1)
1286ccccccc
1287cc
1288       if(lshelkl.eq.1.or.lshelkl.eq.3) then
1289         if(ltyp.eq.1) then
1290           call amtfer(yij(1,ni,nj,1),wij(1,1,1),mnbls)
1291         else
1292      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1293     *   yij(1,ni,nj,2),yij(1,ni,nj,3),yij(1,ni,nj,4),yij(1,ni,nj,1))
1294         endif
1295       endif
1296       if(lshelkl.eq.2.or.lshelkl.eq.3) then
1297         do 312 nkl=nfu(nqkx)+1,nfu(nqkl+1)
1298         call amtfer(zij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
1299  312    continue
1300       endif
1301c
1302       if(lshelkl.eq.3) then
1303         call amtfer(sij(1,ni,nj),wij(1,1,1),mnbls)
1304       endif
1305c
1306       if( ni.le.nfu(nqi +1).and.nj.le.nfu(nqj +1) ) then
1307c
1308c****  2,3  l-shells ****
1309         if(lshelkl.eq.1) then
1310            if(ltyp.eq.1) then
1311               if(ni.eq.1.and.nj.ge.nfu(nqj)+1) then
1312                  call amtfer(x2l1(1,nj,1),wij(1,1,1),mnbls)
1313               endif
1314               if(nj.eq.1.and.ni.ge.nfu(nqi)+1) then
1315                  call amtfer(x2l3(1,ni,1),wij(1,1,1),mnbls)
1316               endif
1317c--test ?      if(ni.eq.1.and.nj.eq.1) then
1318               if(ni.eq.1.and.nj.eq.1 .and. lcas3(1).eq.1) then
1319                  call amtfer(x3l12(1,1),wij(1,1,1),mnbls)
1320               endif
1321            else
1322              if(ni.eq.1.and.nj.ge.nfu(nqj)+1) then
1323      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1324     *  x2l1(1,nj,2),x2l1(1,nj,3),x2l1(1,nj,4),x2l1(1,nj,1))
1325              endif
1326              if(nj.eq.1.and.ni.ge.nfu(nqi)+1) then
1327      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1328     *  x2l3(1,ni,2),x2l3(1,ni,3),x2l3(1,ni,4),x2l3(1,ni,1))
1329              endif
1330c--test ?     if(ni.eq.1.and.nj.eq.1) then
1331              if(ni.eq.1.and.nj.eq.1 .and. lcas3(1).eq.1) then
1332      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1333     *               x3l12(1,2),x3l12(1,3),x3l12(1,4),x3l12(1,1))
1334              endif
1335            endif
1336c
1337         endif
1338         if(lshelkl.eq.2) then
1339              do 313 nkl=nfu(nqkx)+1,nfu(nqkl+1)
1340              if(ni.eq.1)  then
1341                  call amtfer(x2l2(1,nj,nkl),wij(1,nkl,1),mnbls)
1342              endif
1343              if(nj.eq.1)  then
1344                  call amtfer(x2l4(1,ni,nkl),wij(1,nkl,1),mnbls)
1345              endif
1346c--test ?     if(ni.eq.1.and.nj.eq.1) then
1347              if(ni.eq.1.and.nj.eq.1 .and. lcas3(2).eq.1) then
1348                  call amtfer(x3l12(1,nkl),wij(1,nkl,1),mnbls)
1349              endif
1350  313         continue
1351         endif
1352c
1353c****  3 l-shells ****
1354         if(lshelkl.eq.3) then
1355              if(lcas2(1).eq.1.and.ni.eq.1) then
1356      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1357     *  x2l1(1,nj,2),x2l1(1,nj,3),x2l1(1,nj,4),x2l1(1,nj,1))
1358cc
1359              call amtfer(x2l2(1,nj,2),wij(1,2,1),mnbls)
1360              call amtfer(x2l2(1,nj,3),wij(1,3,1),mnbls)
1361              call amtfer(x2l2(1,nj,4),wij(1,4,1),mnbls)
1362              endif
1363              if(lcas2(3).eq.1.and.nj.eq.1) then
1364      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1365     *  x2l3(1,ni,2),x2l3(1,ni,3),x2l3(1,ni,4),x2l3(1,ni,1))
1366cc
1367              call amtfer(x2l4(1,ni,2),wij(1,2,1),mnbls)
1368              call amtfer(x2l4(1,ni,3),wij(1,3,1),mnbls)
1369              call amtfer(x2l4(1,ni,4),wij(1,4,1),mnbls)
1370              endif
1371c
1372              if(lcas3(3).eq.1.and.ni.eq.1) then
1373                  call amtfer(x3l34(1,nj),wij(1,1,1),mnbls)
1374              endif
1375              if(lcas3(4).eq.1.and.nj.eq.1) then
1376                  call amtfer(x3l34(1,ni),wij(1,1,1),mnbls)
1377              endif
1378         endif
1379       endif
1380c
1381ccccccccccccccccccccccccccccccccccc
1382      kxyz=0
1383      do 305 nk=nkbeg,nkend
1384      kxyz=kxyz+1
1385      lxyz=0
1386      do 305 nl=nlbeg,nlend
1387      lxyz=lxyz+1
1388      indx=indxx(ixyz,jxyz,kxyz,lxyz)
1389      call amtfer(wij(1,nk,nl),buf(1,indx),mnbls)
1390  305 continue
1391  300 continue
1392c
1393      return
1394      end
1395c=============================================================
1396      subroutine shift4l(buf,
1397     *                buf2,lt1,lt2,bij1,bij2,bkl1,bkl2,lt3,lt4,
1398     *                bij3,bkl3,b2l1,b2l2,b2l3,b2l4,
1399     *                b3l1,b3l2,b3l3,b3l4,lt5,ssss,
1400     *                wij,lt6,vij,uij,sij,lt7,
1401     *                xij,lt8,lt9,yij,zij,lt10,
1402     *                x2l1,x2l2,x2l3,x2l4,x3l1,x3l2,x3l3,x3l4,mnbls,
1403     *                xab,xcd,
1404     *                indxx,lni1,lnj1,lnk1,lnl1)
1405c***********************************************************
1406c*        when 4 l-shells are present
1407c***********************************************************
1408      implicit real*8 (a-h,o-z)
1409#include "texas_lpar.fh"
1410ctest
1411      character*11 scftype
1412      character*8 where
1413      common /runtype/ scftype,where
1414ctest
1415      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
1416c
1417      common/obarai/
1418     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
1419     * nqi,nqj,nqk,nql,nsij,nskl,
1420     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbex,klbex
1421      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
1422c
1423      dimension x2l1(mnbls,lt3,lt4),x2l2(mnbls,lt3,lt4),
1424     *          x2l3(mnbls,lt3,lt4),x2l4(mnbls,lt3,lt4)
1425      dimension x3l1(mnbls,lt4),x3l2(mnbls,lt4),
1426     *          x3l3(mnbls,lt3),x3l4(mnbls,lt3)
1427c
1428      dimension buf2(mnbls,lt1,lt2),
1429     * bij1(mnbls,lt3,lt2),bij2(mnbls,lt3,lt2),
1430     * bkl1(mnbls,lt1,lt4),bkl2(mnbls,lt1,lt4),
1431     * bij3(mnbls,lt2),bkl3(mnbls,lt1),
1432     * b2l1(mnbls,lt3,lt4),b2l2(mnbls,lt3,lt4),
1433     * b2l3(mnbls,lt3,lt4),b2l4(mnbls,lt3,lt4),
1434cc-> * b3l(mnbls,lt5,4)
1435     * b3l1(mnbls,lt5),b3l2(mnbls,lt5),b3l3(mnbls,lt5),b3l4(mnbls,lt5)
1436      dimension wij(mnbls,lt6,lt7),
1437     * vij(mnbls,lt6,lt7),uij(mnbls,lt6,lt7),sij(mnbls,lt6,lt7)
1438      dimension xij(mnbls,lt8,lt3,lt9),yij(mnbls,lt8,lt10,lt4)
1439     *                            ,zij(mnbls,lt8,lt10,lt4)
1440c   ? dimension ssss(nbls)
1441      dimension ssss(mnbls)
1442      dimension buf(mnbls,*)
1443      dimension xab(mnbls,3),xcd(mnbls,3)
1444      dimension indxx(lni1,lnj1,lnk1,lnl1)
1445c***********************************************************
1446c
1447       niend=nfu(nqi+1)
1448       njend=nfu(nqj+1)
1449c
1450       nkend=nfu(nqk+1)
1451       nlend=nfu(nql+1)
1452c
1453       nqix=1
1454ccccc  nqjx=1
1455       nibeg=1
1456       njbeg=1
1457c
1458       nqkx=1
1459ccccc  nqlx=1
1460       nqklx=1
1461       nkbeg=1
1462       nlbeg=1
1463c
1464c*****
1465c     do 100 nkl=nfu(nqklx)+1,nfu(nskl+1)
1466      do 100 nkl=1,10
1467c
1468c      do 101 ij=nqix,nsij
1469       do 101 ij=1,3
1470       ijbeg=nfu(ij)+1
1471       ijend=nfu(ij+1)
1472           do 102 nij=ijbeg,ijend
1473           call amtfer(buf2(1,nij,nkl),wij(1,nij,1),mnbls)
1474  102      continue
1475c
1476       if( nkl.le.nfu(nqkl+1)) then
1477           do 103 nij=ijbeg,ijend
1478           call amtfer(bkl1(1,nij,nkl),vij(1,nij,1),mnbls)
1479  103      continue
1480           do 104 nij=ijbeg,ijend
1481           call amtfer(bkl2(1,nij,nkl),uij(1,nij,1),mnbls)
1482  104      continue
1483         if(                 nkl.eq.1) then
1484           do 105 nij=ijbeg,ijend
1485           call amtfer(bkl3(1,nij),sij(1,nij,1),mnbls)
1486  105      continue
1487         endif
1488       endif
1489c
1490  101  continue
1491c
1492ccccccccccc
1493c
1494        call horiz12(wij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1495c
1496      if( nkl.le.nfu(nqkl+1)) then
1497            call horiz12(vij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1498            call horiz12(uij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1499      endif
1500      if( nkl.eq.1) then
1501            call horiz12(sij,lt6,lt7,xab,mnbls,nqi,nqj,nsij1)
1502      endif
1503ccccccccccc
1504c
1505        do 107 nj=njbeg,njend
1506        do 107 ni=nibeg,niend
1507        call amtfer(wij(1,ni,nj),xij(1,ni,nj,nkl),mnbls)
1508  107   continue
1509c
1510      if( nkl.le.nfu(nqkl+1)) then
1511        do 1071 nj=njbeg,njend
1512        do 1071 ni=nibeg,niend
1513        call amtfer(vij(1,ni,nj),yij(1,ni,nj,nkl),mnbls)
1514        call amtfer(uij(1,ni,nj),zij(1,ni,nj,nkl),mnbls)
1515 1071   continue
1516      endif
1517c
1518      call daxpy3(mnbls,xab,
1519     *   xij(1,1,2,nkl),xij(1,1,3,nkl),xij(1,1,4,nkl),
1520     *   bij1(1,2,nkl),bij1(1,3,nkl),bij1(1,4,nkl),bij1(1,1,nkl))
1521c
1522           do 108 nij=nfu(nqi )+1,nfu(nqi +1)
1523           call amtfer(bij2(1,nij,nkl),xij(1,nij,1,nkl),mnbls)
1524  108      continue
1525c
1526           call amtfer(bij3(1,nkl),xij(1,1,1,nkl),mnbls)
1527c
1528c*****
1529      if( nkl.le.nfu(nqkl+1)) then
1530      call daxpy3(mnbls,xab,x2l1(1,2,nkl),x2l1(1,3,nkl),x2l1(1,4,nkl),
1531     *    b2l1(1,2,nkl),b2l1(1,3,nkl),b2l1(1,4,nkl),b2l1(1,1,nkl))
1532                call amtfer(b2l3(1,2,nkl),x2l3(1,2,nkl),mnbls)
1533                call amtfer(b2l3(1,3,nkl),x2l3(1,3,nkl),mnbls)
1534                call amtfer(b2l3(1,4,nkl),x2l3(1,4,nkl),mnbls)
1535cc
1536      call daxpy3(mnbls,xab,x2l2(1,2,nkl),x2l2(1,3,nkl),x2l2(1,4,nkl),
1537     *  b2l2(1,2,nkl),b2l2(1,3,nkl),b2l2(1,4,nkl),b2l2(1,1,nkl))
1538cc
1539                call amtfer(b2l4(1,2,nkl),x2l4(1,2,nkl),mnbls)
1540                call amtfer(b2l4(1,3,nkl),x2l4(1,3,nkl),mnbls)
1541                call amtfer(b2l4(1,4,nkl),x2l4(1,4,nkl),mnbls)
1542c
1543                call amtfer(b3l1(1,nkl),x3l1(1,nkl),mnbls)
1544                call amtfer(b3l2(1,nkl),x3l2(1,nkl),mnbls)
1545c
1546      endif
1547      if(nkl.eq.1) then
1548c
1549      call daxpy3(mnbls,xab,x3l3(1,2),x3l3(1,3),x3l3(1,4),
1550     *  b3l3(1,2),b3l3(1,3),b3l3(1,4),b3l3(1,1))
1551c
1552                do 1101 nij=nfu(nqix)+1,nfu(nqij+1)
1553                call amtfer(b3l4(1,nij ),x3l4(1,nij),mnbls)
1554 1101           continue
1555      endif
1556c
1557  100 continue
1558c
1559c***********************************************************
1560c*    this part    shifts the angular momentum
1561c*         from position 3 to position 4
1562c*                         or
1563c*         from position 4 to position 3
1564c***********************************************************
1565c
1566      ixyz=0
1567      do 300 ni=nibeg,niend
1568      ixyz=ixyz+1
1569      jxyz=0
1570      do 300 nj=njbeg,njend
1571      jxyz=jxyz+1
1572c
1573         do 301 nkl=nfu(nqkx)+1,nfu(nskl+1)
1574         call amtfer(xij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
1575  301    continue
1576c
1577ccccccccc
1578c
1579           call horiz12(wij,lt6,lt7,xcd,mnbls,nqk,nql,nskl1)
1580c
1581cccccc
1582      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1583     * yij(1,ni,nj,2),yij(1,ni,nj,3),yij(1,ni,nj,4),yij(1,ni,nj,1))
1584cc
1585         do 312 nkl=nfu(nqkx)+1,nfu(nqkl+1)
1586         call amtfer(zij(1,ni,nj,nkl),wij(1,nkl,1),mnbls)
1587  312    continue
1588c
1589         call amtfer(sij(1,ni,nj),wij(1,1,1),mnbls)
1590c
1591c****
1592       if( ni.le.nfu(nqi +1).and.nj.le.nfu(nqj +1) ) then
1593              if(ni.eq.1 .and. nj.ge.2) then
1594      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1595     *  x2l1(1,nj,2),x2l1(1,nj,3),x2l1(1,nj,4),x2l1(1,nj,1))
1596cc
1597              call amtfer(x2l2(1,nj,2),wij(1,2,1),mnbls)
1598              call amtfer(x2l2(1,nj,3),wij(1,3,1),mnbls)
1599              call amtfer(x2l2(1,nj,4),wij(1,4,1),mnbls)
1600cc
1601              endif
1602              if(nj.eq.1 .and. ni.ge.2) then
1603      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1604     * x2l3(1,ni,2),x2l3(1,ni,3),x2l3(1,ni,4),x2l3(1,ni,1))
1605cc
1606              call amtfer(x2l4(1,ni,2),wij(1,2,1),mnbls)
1607              call amtfer(x2l4(1,ni,3),wij(1,3,1),mnbls)
1608              call amtfer(x2l4(1,ni,4),wij(1,4,1),mnbls)
1609              endif
1610c
1611              if(ni.eq.1) call amtfer(x3l3(1,nj),wij(1,1,1),mnbls)
1612              if(nj.eq.1) call amtfer(x3l4(1,ni),wij(1,1,1),mnbls)
1613       endif
1614       if( ni.eq.1.and.nj.eq.1 ) then
1615      call daxpy3(mnbls,xcd,wij(1,1,2),wij(1,1,3),wij(1,1,4),
1616     *  x3l1(1,2),x3l1(1,3),x3l1(1,4),x3l1(1,1))
1617cc
1618          call amtfer(x3l2(1,2),wij(1,2,1),mnbls)
1619          call amtfer(x3l2(1,3),wij(1,3,1),mnbls)
1620          call amtfer(x3l2(1,4),wij(1,4,1),mnbls)
1621cc
1622          call amtfer(ssss(1),wij(1,1,1),mnbls)
1623       endif
1624c
1625ccccccccccccccccccccccccccccccccccc
1626      kxyz=0
1627      do 305 nk=nkbeg,nkend
1628      kxyz=kxyz+1
1629      lxyz=0
1630      do 305 nl=nlbeg,nlend
1631      lxyz=lxyz+1
1632      indx=indxx(ixyz,jxyz,kxyz,lxyz)
1633      call amtfer(wij(1,nk,nl),buf(1,indx),mnbls)
1634  305 continue
1635  300 continue
1636c
1637      return
1638      end
1639c=====================================================
1640      subroutine horiz12(wij,lw1,lw2,xab,mnbls,nqi,nqj,nsij1)
1641      implicit real*8 (a-h,o-z)
1642c
1643#include "texas_lpar.fh"
1644c
1645      dimension wij(mnbls,lw1,lw2),xab(mnbls,3)
1646c---------------------------------------------------
1647          do 110 j=2,nqj
1648          jbeg=nfu(j)+1
1649          jend=nfu(j+1)
1650             do 115 i=nsij1-j,nqi,-1
1651             ibeg=nfu(i)+1
1652             iend=nfu(i+1)
1653                do 120 nj=jbeg,jend
1654                njm=ilast(nj)
1655                kcr=icool(nj)
1656                   do 125 ni=ibeg,iend
1657                   nij=npxyz(kcr,ni)
1658                      do 130 n=1,mnbls
1659        wij(n,ni,nj)=wij(n,nij,njm)+ xab(n,kcr) *wij(n,ni,njm)
1660  130                 continue
1661  125              continue
1662  120          continue
1663  115        continue
1664  110     continue
1665c
1666      end
1667c=====================================================
1668      subroutine shift0l_der1(buf,buf2,buf0,lt1,lt2,
1669     *                   wij,wij0,lt3,lsjl,
1670     *                   xij,xij0,lt4,lt5,lt6,mnbls,nbls,
1671     *                   xab,xcd,
1672     *                   indxx,lni1,lnj1,lnk1,lnl1)
1673c------------------------------------
1674c  when l-shells are not present
1675c------------------------------------
1676      implicit real*8 (a-h,o-z)
1677#include "texas_lpar.fh"
1678      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
1679      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
1680      common/obarai/
1681     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
1682     * nqi,nqj,nqk,nql,nsij,nskl,
1683     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbex,klbex
1684c
1685      dimension buf0(nbls,lt1,lt2),wij0(nbls,lt3,lsjl)
1686      dimension buf2(mnbls,lt1,lt2),wij(mnbls,lt3,lsjl)
1687      dimension xij(mnbls,lt4,lt5,lt6), xij0(nbls,lt4,lt5,lt6)
1688      dimension buf(mnbls,*)
1689      dimension xab(mnbls,3),xcd(mnbls,3)
1690      dimension indxx(lni1,lnj1,lnk1,lnl1)
1691c------------------------------------
1692c
1693      nibeg=nfu(nqi)+1
1694      niend=nfu(nqi+1)
1695      njbeg=nfu(nqj)+1
1696      njend=nfu(nqj+1)
1697cc
1698      nkbeg=nfu(nqk)+1
1699      nkend=nfu(nqk+1)
1700      nlbeg=nfu(nql)+1
1701      nlend=nfu(nql+1)
1702c
1703      nqix=nqi
1704ccccc nqjx=nqj
1705c
1706      nqkx=nqk
1707ccccc nqlx=nql
1708      nqklx=nqkl
1709c
1710      do 100 nkl=nfu(nqklx)+1,nfu(nskl+1)
1711c
1712           do 102 nij=nfu(nqix)+1,nfu(nsij+1)
1713           call amtfer(buf2(1,nij,nkl),wij(1,nij,1),mnbls)
1714           call amtfer(buf0(1,nij,nkl),wij0(1,nij,1),nbls)
1715  102      continue
1716c
1717           call horiz12_der1(wij0,wij,lt3,lsjl,xab, nbls,nqi,nqj,nsij1)
1718c
1719           do 107 nj=njbeg,njend
1720           do 107 ni=nibeg,niend
1721           call amtfer(wij(1,ni,nj),xij(1,ni,nj,nkl),mnbls)
1722           call amtfer(wij0(1,ni,nj),xij0(1,ni,nj,nkl), nbls)
1723  107      continue
1724  100 continue
1725c
1726c------------------------------------
1727c this part shifts angular momentum
1728c   from position 3 to position 4
1729c------------------------------------
1730c
1731      ixyz=0
1732      do 300 ni=nibeg,niend
1733      ixyz=ixyz+1
1734      jxyz=0
1735      do 300 nj=njbeg,njend
1736      jxyz=jxyz+1
1737c
1738         do 301 nkl=nfu(nqkx)+1,nfu(nskl+1)
1739         call amtfer( xij(1,ni,nj,nkl), wij(1,nkl,1),mnbls)
1740         call amtfer(xij0(1,ni,nj,nkl),wij0(1,nkl,1), nbls)
1741  301    continue
1742c
1743         call horiz34_der1(wij0,wij,lt3,lsjl,xcd, nbls,nqk,nql,nskl1)
1744c
1745      kxyz=0
1746      do 305 nk=nkbeg,nkend
1747      kxyz=kxyz+1
1748      lxyz=0
1749      do 305 nl=nlbeg,nlend
1750      lxyz=lxyz+1
1751      indx=indxx(ixyz,jxyz,kxyz,lxyz)
1752      call amtfer(wij(1,nk,nl),buf(1,indx),mnbls)
1753  305 continue
1754  300 continue
1755c
1756      end
1757c=============================================================
1758      subroutine horiz12_der1(wij0,wij,lw1,lw2,xab,nbls,nqi,nqj,nsij1)
1759      implicit real*8 (a-h,o-z)
1760c
1761#include "texas_lpar.fh"
1762c
1763      dimension wij0(nbls,lw1,lw2)
1764      dimension wij(9,nbls,lw1,lw2),xab(9,nbls,3)
1765c---------------------------------------------------
1766          do 110 j=2,nqj
1767          jbeg=nfu(j)+1
1768          jend=nfu(j+1)
1769             do 115 i=nsij1-j,nqi,-1
1770             ibeg=nfu(i)+1
1771             iend=nfu(i+1)
1772                do 120 nj=jbeg,jend
1773                njm=ilast(nj)
1774                kcr=icool(nj)
1775                   do 125 ni=ibeg,iend
1776                   nij=npxyz(kcr,ni)
1777                      do 130 n=1, nbls
1778c Ax
1779                      wij(1,n,ni,nj)=wij(1,n,nij,njm)
1780     *                              + xab(1,n,kcr)*wij(1,n,ni,njm)
1781c Bx
1782                      wij(2,n,ni,nj)=wij(2,n,nij,njm)
1783     *                              + xab(2,n,kcr)*wij(2,n,ni,njm)
1784c Cx
1785                      wij(3,n,ni,nj)=wij(3,n,nij,njm)
1786     *                              + xab(3,n,kcr)*wij(3,n,ni,njm)
1787c Ay
1788                      wij(4,n,ni,nj)=wij(4,n,nij,njm)
1789     *                              + xab(4,n,kcr)*wij(4,n,ni,njm)
1790c By
1791                      wij(5,n,ni,nj)=wij(5,n,nij,njm)
1792     *                              + xab(5,n,kcr)*wij(5,n,ni,njm)
1793c Cy
1794                      wij(6,n,ni,nj)=wij(6,n,nij,njm)
1795     *                              + xab(6,n,kcr)*wij(6,n,ni,njm)
1796c Az
1797                      wij(7,n,ni,nj)=wij(7,n,nij,njm)
1798     *                              + xab(7,n,kcr)*wij(7,n,ni,njm)
1799c Bz
1800                      wij(8,n,ni,nj)=wij(8,n,nij,njm)
1801     *                              + xab(8,n,kcr)*wij(8,n,ni,njm)
1802c Cz
1803                      wij(9,n,ni,nj)=wij(9,n,nij,njm)
1804     *                              + xab(9,n,kcr)*wij(9,n,ni,njm)
1805c
1806c add an extra term arising from differentiation of the shifting formula
1807c ONLY for derivatives over center A and B , not C :
1808c
1809                      wij0(n,ni,nj)=wij0(n,nij,njm)
1810     *                              + xab(1,n,kcr)*wij0(n,ni,njm)
1811                      if(kcr.eq.1) then
1812                         wij(1,n,ni,nj)=wij(1,n,ni,nj) + wij0(n,ni,njm)
1813                         wij(2,n,ni,nj)=wij(2,n,ni,nj) - wij0(n,ni,njm)
1814                      endif
1815                      if(kcr.eq.2) then
1816                         wij(4,n,ni,nj)=wij(4,n,ni,nj) + wij0(n,ni,njm)
1817                         wij(5,n,ni,nj)=wij(5,n,ni,nj) - wij0(n,ni,njm)
1818                      endif
1819                      if(kcr.eq.3) then
1820                         wij(7,n,ni,nj)=wij(7,n,ni,nj) + wij0(n,ni,njm)
1821                         wij(8,n,ni,nj)=wij(8,n,ni,nj) - wij0(n,ni,njm)
1822                      endif
1823  130                 continue
1824  125              continue
1825  120          continue
1826  115        continue
1827  110     continue
1828c
1829      end
1830c=====================================================
1831      subroutine horiz34_der1(wij0,wij,lw1,lw2,xab,nbls,nqi,nqj,nsij1)
1832      implicit real*8 (a-h,o-z)
1833c
1834#include "texas_lpar.fh"
1835c
1836      dimension wij0(nbls,lw1,lw2)
1837      dimension wij(9,nbls,lw1,lw2),xab(9,nbls,3)
1838c---------------------------------------------------
1839          do 110 j=2,nqj
1840          jbeg=nfu(j)+1
1841          jend=nfu(j+1)
1842             do 115 i=nsij1-j,nqi,-1
1843             ibeg=nfu(i)+1
1844             iend=nfu(i+1)
1845                do 120 nj=jbeg,jend
1846                njm=ilast(nj)
1847                kcr=icool(nj)
1848                   do 125 ni=ibeg,iend
1849                   nij=npxyz(kcr,ni)
1850                      do 130 n=1, nbls
1851c Ax
1852                      wij(1,n,ni,nj)=wij(1,n,nij,njm)
1853     *                              + xab(1,n,kcr)*wij(1,n,ni,njm)
1854c Bx
1855                      wij(2,n,ni,nj)=wij(2,n,nij,njm)
1856     *                              + xab(2,n,kcr)*wij(2,n,ni,njm)
1857c Cx
1858                      wij(3,n,ni,nj)=wij(3,n,nij,njm)
1859     *                              + xab(3,n,kcr)*wij(3,n,ni,njm)
1860c Ay
1861                      wij(4,n,ni,nj)=wij(4,n,nij,njm)
1862     *                              + xab(4,n,kcr)*wij(4,n,ni,njm)
1863c By
1864                      wij(5,n,ni,nj)=wij(5,n,nij,njm)
1865     *                              + xab(5,n,kcr)*wij(5,n,ni,njm)
1866c Cy
1867                      wij(6,n,ni,nj)=wij(6,n,nij,njm)
1868     *                              + xab(6,n,kcr)*wij(6,n,ni,njm)
1869c Az
1870                      wij(7,n,ni,nj)=wij(7,n,nij,njm)
1871     *                              + xab(7,n,kcr)*wij(7,n,ni,njm)
1872c Bz
1873                      wij(8,n,ni,nj)=wij(8,n,nij,njm)
1874     *                              + xab(8,n,kcr)*wij(8,n,ni,njm)
1875c Cz
1876                      wij(9,n,ni,nj)=wij(9,n,nij,njm)
1877     *                              + xab(9,n,kcr)*wij(9,n,ni,njm)
1878c
1879c add an extra term arising from differentiation of the shifting formula
1880c ONLY for derivatives over center C, not over A and B:
1881c
1882                      wij0(n,ni,nj)=wij0(n,nij,njm)
1883     *                              + xab(1,n,kcr)*wij0(n,ni,njm)
1884                      if(kcr.eq.1) then
1885                         wij(3,n,ni,nj)=wij(3,n,ni,nj) + wij0(n,ni,njm)
1886                      endif
1887                      if(kcr.eq.2) then
1888                         wij(6,n,ni,nj)=wij(6,n,ni,nj) + wij0(n,ni,njm)
1889                      endif
1890                      if(kcr.eq.3) then
1891                         wij(9,n,ni,nj)=wij(9,n,ni,nj) + wij0(n,ni,njm)
1892                      endif
1893  130                 continue
1894  125              continue
1895  120          continue
1896  115        continue
1897  110     continue
1898c
1899      end
1900c=====================================================
1901      subroutine shift0l_der2(buf,buf2,buf1,buf0,lt1,lt2,
1902     *                   wij,wij1,wij0,lt3,lsjl,
1903     *                   xij,xij1,xij0,lt4,lt5,lt6,mnbls,nbls,
1904     *                   xab,xcd,
1905     *                   indxx,lni1,lnj1,lnk1,lnl1)
1906c------------------------------------
1907c  when l-shells are not present
1908c------------------------------------
1909      implicit real*8 (a-h,o-z)
1910      common /types/iityp,jjtyp,kktyp,lltyp, ityp,jtyp,ktyp,ltyp
1911      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
1912      common/obarai/
1913     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
1914     * nqi,nqj,nqk,nql,nsij,nskl,
1915     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbex,klbex
1916#include "texas_lpar.fh"
1917c
1918      dimension buf0(nbls,lt1,lt2),wij0(nbls,lt3,lsjl)
1919      dimension buf1(9*nbls,lt1,lt2),wij1(9*nbls,lt3,lsjl)
1920      dimension buf2(mnbls,lt1,lt2),wij(mnbls,lt3,lsjl)
1921      dimension xij(mnbls,lt4,lt5,lt6),xij1(9*nbls,lt4,lt5,lt6),
1922     *                                 xij0(  nbls,lt4,lt5,lt6)
1923      dimension buf(mnbls,*)
1924      dimension xab(mnbls,3),xcd(mnbls,3)
1925      dimension indxx(lni1,lnj1,lnk1,lnl1)
1926c------------------------------------
1927      nbls9=9*nbls
1928c
1929      nibeg=nfu(nqi)+1
1930      niend=nfu(nqi+1)
1931      njbeg=nfu(nqj)+1
1932      njend=nfu(nqj+1)
1933cc
1934      nkbeg=nfu(nqk)+1
1935      nkend=nfu(nqk+1)
1936      nlbeg=nfu(nql)+1
1937      nlend=nfu(nql+1)
1938c
1939      nqix=nqi
1940ccccc nqjx=nqj
1941c
1942      nqkx=nqk
1943ccccc nqlx=nql
1944      nqklx=nqkl
1945c
1946      do 100 nkl=nfu(nqklx)+1,nfu(nskl+1)
1947c
1948           do 102 nij=nfu(nqix)+1,nfu(nsij+1)
1949           call amtfer(buf2(1,nij,nkl),wij(1,nij,1),mnbls)
1950           call amtfer(buf1(1,nij,nkl),wij1(1,nij,1),nbls9)
1951           call amtfer(buf0(1,nij,nkl),wij0(1,nij,1),nbls)
1952  102      continue
1953c
1954           call horiz12_der2(wij0,wij1,wij,
1955     *                       lt3,lsjl,xab, nbls,nqi,nqj,nsij1)
1956c
1957           do 107 nj=njbeg,njend
1958           do 107 ni=nibeg,niend
1959           call amtfer(wij(1,ni,nj),xij(1,ni,nj,nkl),mnbls)
1960           call amtfer(wij1(1,ni,nj),xij1(1,ni,nj,nkl),nbls9)
1961           call amtfer(wij0(1,ni,nj),xij0(1,ni,nj,nkl),nbls)
1962  107      continue
1963  100 continue
1964c
1965c------------------------------------
1966c this part shifts angular momentum
1967c   from position 3 to position 4
1968c------------------------------------
1969c
1970      ixyz=0
1971      do 300 ni=nibeg,niend
1972      ixyz=ixyz+1
1973      jxyz=0
1974      do 300 nj=njbeg,njend
1975      jxyz=jxyz+1
1976c
1977         do 301 nkl=nfu(nqkx)+1,nfu(nskl+1)
1978         call amtfer( xij(1,ni,nj,nkl), wij(1,nkl,1),mnbls)
1979         call amtfer(xij1(1,ni,nj,nkl),wij1(1,nkl,1),nbls9)
1980         call amtfer(xij0(1,ni,nj,nkl),wij0(1,nkl,1), nbls)
1981  301    continue
1982c
1983         call horiz34_der2(wij0,wij1,wij,
1984     *                     lt3,lsjl,xcd, nbls,nqk,nql,nskl1)
1985c
1986      kxyz=0
1987      do 305 nk=nkbeg,nkend
1988      kxyz=kxyz+1
1989      lxyz=0
1990      do 305 nl=nlbeg,nlend
1991      lxyz=lxyz+1
1992      indx=indxx(ixyz,jxyz,kxyz,lxyz)
1993      call amtfer(wij(1,nk,nl),buf(1,indx),mnbls)
1994  305 continue
1995  300 continue
1996c
1997      end
1998c=============================================================
1999      subroutine horiz12_der2(wij0,wij1,wij2,
2000     *                        lw1,lw2,xab,nbls,nqi,nqj,nsij1)
2001      implicit real*8 (a-h,o-z)
2002c
2003#include "texas_lpar.fh"
2004c
2005      dimension wij0(nbls,lw1,lw2)
2006      dimension wij1(9,nbls,lw1,lw2)
2007      dimension wij2(45,nbls,lw1,lw2),xab(45,nbls,3)
2008c---------------------------------------------------
2009          do 110 j=2,nqj
2010          jbeg=nfu(j)+1
2011          jend=nfu(j+1)
2012             do 115 i=nsij1-j,nqi,-1
2013             ibeg=nfu(i)+1
2014             iend=nfu(i+1)
2015                do 120 nj=jbeg,jend
2016                njm=ilast(nj)
2017                kcr=icool(nj)
2018                   do 125 ni=ibeg,iend
2019                   nij=npxyz(kcr,ni)
2020                      do 130 n=1, nbls
2021      wij0(n,ni,nj)=wij0(n,nij,njm) + xab(1,n,kcr)*wij0(n,ni,njm)
2022                      do m=1,9
2023      wij1(m,n,ni,nj)=wij1(m,n,nij,njm) + xab(m,n,kcr)*wij1(m,n,ni,njm)
2024                      enddo
2025                      do m=1,45
2026      wij2(m,n,ni,nj)=wij2(m,n,nij,njm) + xab(m,n,kcr)*wij2(m,n,ni,njm)
2027                      enddo
2028c
2029c add an extra term arising from differentiation of the shifting formula
2030c ONLY for derivatives over center A and B , not C :
2031c
2032c  first derivatives:
2033c
2034          if(kcr.eq.1) then
2035             wij1(1,n,ni,nj)=wij1(1,n,ni,nj) + wij0(n,ni,njm)
2036             wij1(2,n,ni,nj)=wij1(2,n,ni,nj) - wij0(n,ni,njm)
2037c
2038          endif
2039          if(kcr.eq.2) then
2040             wij1(4,n,ni,nj)=wij1(4,n,ni,nj) + wij0(n,ni,njm)
2041             wij1(5,n,ni,nj)=wij1(5,n,ni,nj) - wij0(n,ni,njm)
2042          endif
2043          if(kcr.eq.3) then
2044             wij1(7,n,ni,nj)=wij1(7,n,ni,nj) + wij0(n,ni,njm)
2045             wij1(8,n,ni,nj)=wij1(8,n,ni,nj) - wij0(n,ni,njm)
2046          endif
2047c
2048c second derivatives:
2049c
2050          if(kcr.eq.1) then
2051c block aa:
2052             wij2(1,n,ni,nj)=wij2(1,n,ni,nj) + wij1(1,n,ni,njm)*2.d0
2053             wij2(2,n,ni,nj)=wij2(2,n,ni,nj) + wij1(4,n,ni,njm)
2054             wij2(3,n,ni,nj)=wij2(3,n,ni,nj) + wij1(7,n,ni,njm)
2055c block ab:
2056             wij2(7,n,ni,nj) =wij2(7,n,ni,nj) - wij1(1,n,ni,njm)
2057     *                                        + wij1(2,n,ni,njm)
2058             wij2(8,n,ni,nj) =wij2(8,n,ni,nj) + wij1(5,n,ni,njm)
2059             wij2(9,n,ni,nj) =wij2(9,n,ni,nj) + wij1(8,n,ni,njm)
2060             wij2(10,n,ni,nj)=wij2(10,n,ni,nj)- wij1(4,n,ni,njm)
2061             wij2(13,n,ni,nj)=wij2(13,n,ni,nj)- wij1(7,n,ni,njm)
2062c block ac:
2063             wij2(16,n,ni,nj)=wij2(16,n,ni,nj)+ wij1(3,n,ni,njm)
2064             wij2(17,n,ni,nj)=wij2(17,n,ni,nj)+ wij1(6,n,ni,njm)
2065             wij2(18,n,ni,nj)=wij2(18,n,ni,nj)+ wij1(9,n,ni,njm)
2066c block bb:
2067             wij2(25,n,ni,nj)=wij2(25,n,ni,nj)- wij1(2,n,ni,njm)*2.d0
2068             wij2(26,n,ni,nj)=wij2(26,n,ni,nj)- wij1(5,n,ni,njm)
2069             wij2(27,n,ni,nj)=wij2(27,n,ni,nj)- wij1(8,n,ni,njm)
2070c block bc:
2071             wij2(31,n,ni,nj)=wij2(31,n,ni,nj)- wij1(3,n,ni,njm)
2072             wij2(32,n,ni,nj)=wij2(32,n,ni,nj)- wij1(6,n,ni,njm)
2073             wij2(33,n,ni,nj)=wij2(33,n,ni,nj)- wij1(9,n,ni,njm)
2074c block cc: no contributions of this kind
2075          endif
2076          if(kcr.eq.2) then
2077c block aa:
2078             wij2(2,n,ni,nj)=wij2(2,n,ni,nj) + wij1(1,n,ni,njm)
2079             wij2(4,n,ni,nj)=wij2(4,n,ni,nj) + wij1(4,n,ni,njm)*2.d0
2080             wij2(5,n,ni,nj)=wij2(5,n,ni,nj) + wij1(7,n,ni,njm)
2081c block ab:
2082             wij2(8,n,ni,nj) =wij2(8,n,ni,nj) - wij1(1,n,ni,njm)
2083             wij2(10,n,ni,nj)=wij2(10,n,ni,nj)+ wij1(2,n,ni,njm)
2084             wij2(11,n,ni,nj)=wij2(11,n,ni,nj)- wij1(4,n,ni,njm)
2085     *                                        + wij1(5,n,ni,njm)
2086             wij2(12,n,ni,nj)=wij2(12,n,ni,nj)+ wij1(8,n,ni,njm)
2087             wij2(14,n,ni,nj)=wij2(14,n,ni,nj)- wij1(7,n,ni,njm)
2088c block ac:
2089             wij2(19,n,ni,nj)=wij2(19,n,ni,nj)+ wij1(3,n,ni,njm)
2090             wij2(20,n,ni,nj)=wij2(20,n,ni,nj)+ wij1(6,n,ni,njm)
2091             wij2(21,n,ni,nj)=wij2(21,n,ni,nj)+ wij1(9,n,ni,njm)
2092c block bb:
2093             wij2(26,n,ni,nj)=wij2(26,n,ni,nj)- wij1(2,n,ni,njm)
2094             wij2(28,n,ni,nj)=wij2(28,n,ni,nj)- wij1(5,n,ni,njm)*2.d0
2095             wij2(29,n,ni,nj)=wij2(29,n,ni,nj)- wij1(8,n,ni,njm)
2096c block bc:
2097             wij2(34,n,ni,nj)=wij2(34,n,ni,nj)- wij1(3,n,ni,njm)
2098             wij2(35,n,ni,nj)=wij2(35,n,ni,nj)- wij1(6,n,ni,njm)
2099             wij2(36,n,ni,nj)=wij2(36,n,ni,nj)- wij1(9,n,ni,njm)
2100c block cc: no contributions of this kind
2101          endif
2102          if(kcr.eq.3) then
2103c block aa:
2104             wij2(3,n,ni,nj)=wij2(3,n,ni,nj) + wij1(1,n,ni,njm)
2105             wij2(5,n,ni,nj)=wij2(5,n,ni,nj) + wij1(4,n,ni,njm)
2106             wij2(6,n,ni,nj)=wij2(6,n,ni,nj) + wij1(7,n,ni,njm)*2.d0
2107c block ab:
2108             wij2(9,n,ni,nj) =wij2(9,n,ni,nj) - wij1(1,n,ni,njm)
2109             wij2(12,n,ni,nj)=wij2(12,n,ni,nj)- wij1(4,n,ni,njm)
2110             wij2(13,n,ni,nj)=wij2(13,n,ni,nj)+ wij1(2,n,ni,njm)
2111             wij2(14,n,ni,nj)=wij2(14,n,ni,nj)+ wij1(5,n,ni,njm)
2112             wij2(15,n,ni,nj)=wij2(15,n,ni,nj)- wij1(7,n,ni,njm)
2113     *                                        + wij1(8,n,ni,njm)
2114c block ac:
2115             wij2(22,n,ni,nj)=wij2(22,n,ni,nj)+ wij1(3,n,ni,njm)
2116             wij2(23,n,ni,nj)=wij2(23,n,ni,nj)+ wij1(6,n,ni,njm)
2117             wij2(24,n,ni,nj)=wij2(24,n,ni,nj)+ wij1(9,n,ni,njm)
2118c block bb:
2119             wij2(27,n,ni,nj)=wij2(27,n,ni,nj)- wij1(2,n,ni,njm)
2120             wij2(29,n,ni,nj)=wij2(29,n,ni,nj)- wij1(5,n,ni,njm)
2121             wij2(30,n,ni,nj)=wij2(30,n,ni,nj)- wij1(8,n,ni,njm)*2.d0
2122c block bc:
2123             wij2(37,n,ni,nj)=wij2(37,n,ni,nj)- wij1(3,n,ni,njm)
2124             wij2(38,n,ni,nj)=wij2(38,n,ni,nj)- wij1(6,n,ni,njm)
2125             wij2(39,n,ni,nj)=wij2(39,n,ni,nj)- wij1(9,n,ni,njm)
2126c block cc: no contributions of this kind
2127          endif
2128  130                 continue
2129  125              continue
2130  120          continue
2131  115        continue
2132  110     continue
2133c
2134      end
2135c=====================================================
2136      subroutine horiz34_der2(wij0,wij1,wij2,
2137     *                        lw1,lw2,xab,nbls,nqi,nqj,nsij1)
2138      implicit real*8 (a-h,o-z)
2139c
2140#include "texas_lpar.fh"
2141c
2142      dimension wij0(nbls,lw1,lw2)
2143      dimension wij1(9,nbls,lw1,lw2)
2144      dimension wij2(45,nbls,lw1,lw2),xab(45,nbls,3)
2145c---------------------------------------------------
2146          do 110 j=2,nqj
2147          jbeg=nfu(j)+1
2148          jend=nfu(j+1)
2149             do 115 i=nsij1-j,nqi,-1
2150             ibeg=nfu(i)+1
2151             iend=nfu(i+1)
2152                do 120 nj=jbeg,jend
2153                njm=ilast(nj)
2154                kcr=icool(nj)
2155                   do 125 ni=ibeg,iend
2156                   nij=npxyz(kcr,ni)
2157                      do 130 n=1, nbls
2158      wij0(n,ni,nj)=wij0(n,nij,njm) + xab(1,n,kcr)*wij0(n,ni,njm)
2159                         do m=1,9
2160      wij1(m,n,ni,nj)=wij1(m,n,nij,njm) + xab(m,n,kcr)*wij1(m,n,ni,njm)
2161                         enddo
2162                         do m=1,45
2163      wij2(m,n,ni,nj)=wij2(m,n,nij,njm) + xab(m,n,kcr)*wij2(m,n,ni,njm)
2164                         enddo
2165c
2166c add an extra term arising from differentiation of the shifting formula
2167c ONLY for derivatives over center C, not over A and B:
2168c
2169c first derivatives:
2170c
2171      if(kcr.eq.1) then
2172         wij1(3,n,ni,nj)=wij1(3,n,ni,nj) + wij0(n,ni,njm)
2173      endif
2174      if(kcr.eq.2) then
2175         wij1(6,n,ni,nj)=wij1(6,n,ni,nj) + wij0(n,ni,njm)
2176      endif
2177      if(kcr.eq.3) then
2178         wij1(9,n,ni,nj)=wij1(9,n,ni,nj) + wij0(n,ni,njm)
2179      endif
2180c
2181c second derivatives:
2182c
2183          if(kcr.eq.1) then
2184c block ac:
2185             wij2(16,n,ni,nj)=wij2(16,n,ni,nj)+ wij1(1,n,ni,njm)
2186             wij2(19,n,ni,nj)=wij2(19,n,ni,nj)+ wij1(4,n,ni,njm)
2187             wij2(22,n,ni,nj)=wij2(22,n,ni,nj)+ wij1(7,n,ni,njm)
2188c block bc:
2189             wij2(31,n,ni,nj)=wij2(31,n,ni,nj)+ wij1(2,n,ni,njm)
2190             wij2(34,n,ni,nj)=wij2(34,n,ni,nj)+ wij1(5,n,ni,njm)
2191             wij2(37,n,ni,nj)=wij2(37,n,ni,nj)+ wij1(8,n,ni,njm)
2192c block cc:
2193             wij2(40,n,ni,nj)=wij2(40,n,ni,nj)+ wij1(3,n,ni,njm)*2.d0
2194             wij2(41,n,ni,nj)=wij2(41,n,ni,nj)+ wij1(6,n,ni,njm)
2195             wij2(42,n,ni,nj)=wij2(42,n,ni,nj)+ wij1(9,n,ni,njm)
2196          endif
2197          if(kcr.eq.2) then
2198c block ac:
2199             wij2(17,n,ni,nj)=wij2(17,n,ni,nj)+ wij1(1,n,ni,njm)
2200             wij2(20,n,ni,nj)=wij2(20,n,ni,nj)+ wij1(4,n,ni,njm)
2201             wij2(23,n,ni,nj)=wij2(23,n,ni,nj)+ wij1(7,n,ni,njm)
2202c block bc:
2203             wij2(32,n,ni,nj)=wij2(32,n,ni,nj)+ wij1(2,n,ni,njm)
2204             wij2(35,n,ni,nj)=wij2(35,n,ni,nj)+ wij1(5,n,ni,njm)
2205             wij2(38,n,ni,nj)=wij2(38,n,ni,nj)+ wij1(8,n,ni,njm)
2206c block cc:
2207             wij2(41,n,ni,nj)=wij2(41,n,ni,nj)+ wij1(3,n,ni,njm)
2208             wij2(43,n,ni,nj)=wij2(43,n,ni,nj)+ wij1(6,n,ni,njm)*2.d0
2209             wij2(44,n,ni,nj)=wij2(44,n,ni,nj)+ wij1(9,n,ni,njm)
2210          endif
2211          if(kcr.eq.3) then
2212c block ac:
2213             wij2(18,n,ni,nj)=wij2(18,n,ni,nj)+ wij1(1,n,ni,njm)
2214             wij2(21,n,ni,nj)=wij2(21,n,ni,nj)+ wij1(4,n,ni,njm)
2215             wij2(24,n,ni,nj)=wij2(24,n,ni,nj)+ wij1(7,n,ni,njm)
2216c block bc:
2217             wij2(33,n,ni,nj)=wij2(33,n,ni,nj)+ wij1(2,n,ni,njm)
2218             wij2(36,n,ni,nj)=wij2(36,n,ni,nj)+ wij1(5,n,ni,njm)
2219             wij2(39,n,ni,nj)=wij2(39,n,ni,nj)+ wij1(8,n,ni,njm)
2220c block cc:
2221             wij2(42,n,ni,nj)=wij2(42,n,ni,nj)+ wij1(3,n,ni,njm)
2222             wij2(44,n,ni,nj)=wij2(44,n,ni,nj)+ wij1(6,n,ni,njm)
2223             wij2(45,n,ni,nj)=wij2(45,n,ni,nj)+ wij1(9,n,ni,njm)*2.d0
2224          endif
2225  130                 continue
2226  125              continue
2227  120          continue
2228  115        continue
2229  110     continue
2230c
2231      end
2232c=======================================================================
2233      subroutine make_indxx(indxx,lni,lnj,lnk,lnl)
2234      dimension indxx(lni,lnj,lnk,lnl)
2235c
2236      ncount=0
2237      do 20 i=1,lni
2238      do 20 j=1,lnj
2239      do 20 k=1,lnk
2240      do 20 l=1,lnl
2241         ncount=ncount+1
2242      indxx(i,j,k,l)=ncount
2243   20 continue
2244c
2245      end
2246c=======================================================================
2247