1* $Id$
2c=======================================================================
3c     For the set of NBLS1 quartets of primitive shells I,J,K,L
4c  (where NBLS1 is already reduced after some integrals have
5c  been neglected) two-electron integrals of the type (i+j,s|k+l,s)
6c  are calculated according to the OBara-SAIka (OBASAI) method.
7c  The algorythm uses also the method proposed by Tracy Hamilton
8c  (TRACY) to shift angular momenta from electron 1 (centers 1, 2 ) to
9c  electron 2 (centers 3,4).
10c  This set of routines is called when the total angular momentum (+1)
11c  mmax = i+j+k+l +1 is GT.2 (for mmax.le.2 a special code is used)
12c
13c     As a first step the (s,s|s,s)(m) , m=1,mmax integrals
14c  are calculated in the SSSSM routine using the FM routine.
15c  All necessary data for FM are sent by the common block RYS.
16c  These integrals are stored in the wt1(nbls1,mmax,l12) matrix
17c  and then used as a input for OBASAI routines which are called
18c  from here. Part of (ss|ss)(m) itegrals, namely these with
19c  m=1, are stored in the wt0(nbls1,l01,l02) matrix where all
20c  final primitive integrals are kept for futher contraction.
21c
22c     In a second step integrals of the form (i+j+k+l,s|s,s)(m)
23c  or (s,s|i+j+k+l,s)(m) are calculated in the OBASAI routines.
24c  The first type of integrals is calculacted if the total angular
25c  momentum of the first pair ij is greater than the second pair kl
26c  ( nsij >= nskl ). Otherwise, integrals in the second form are
27c  constructed. Moreover, the presence of L-shells in the first and
28c  second pair positions is taken into account. Overall, there are
29c  4 cases for which the OBASAI routines are called with different set
30c  of parameters. These cases are determined in the HOW_2_SHIFT routine.
31c
32c     As a last step the final (i+j,s|k+l,s) integrals are calculated
33c  from the previous ones by shifting angular momenta from position
34c  1 to 3 or 3 to 1. This is performed in the TRAC12 or TRAC34 routines.
35c
36c     Final primitive integrals (i+j,s|k+l,s) return in the wt0
37c  matrix and then they are contracted in the ASSEMBLX routine.
38c
39c  Everything is essentially doubled because of two different blocking
40c  strategies (iroute=1 (texas93) & iroute=2 (texas95) ).
41c=======================================================================
42c For IROUTE=2 :
43c
44c According to the new blocking criterions six arrays
45c  abnia(nbls1,mmax), cdnia(nbls1,mmax),
46c             and
47c  rhoapb(nbls1), rhocpd(nbls1)
48c             and
49c    abcd(nbls1), habcd(nbls1,3,*)
50c
51c are now different. The first dimension is ALWAYS one.
52c Thus, now we have :
53c  abnia(mmax) and cdnia(mmax)
54c  rhoapb , rhocpd
55c  abcd   , habcd(3,*)
56c=======================================================================
57      subroutine trobsa_2(bl,nbls1,l11,l12,mem2)
58C-----------------------------------------------------------------------
59c trobsa_2 is used for iroute=2 only
60C-----------------------------------------------------------------------
61      implicit real*8 (a-h,o-z)
62      logical stable
63      common /tracy_stability/ stable(5000) ! this is limit limblks
64c
65      common/obarai/
66     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
67     * nqi,nqj,nqk,nql,nsij,nskl,
68     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg
69c
70c     common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
71c
72      common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2,
73     * ibfij1,ibfij2,ibfkl1,ibfkl2,
74     * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3,
75     * ibf3l,issss,
76     * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4,
77     * ixij,iyij,izij, iwij,ivij,iuij,isij
78c
79      common /memor5b/ irppq,
80     * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234,
81     * idx1,idx2,indx
82c
83      common /memor5c/ itxab,itxcd,iabcd,ihabcd
84      common /memor5d/ iabnix,icdnix,ixpnx,ixqnx,ihabcdx
85c
86      dimension bl(*)
87C-----------------------------------------------------------------------
88c test : counter
89      common /lobsa_times/ xlobsa2,xlobsa4
90C-----------------------------------------------------------------------
91c PARAMETERS
92c-----------
93c Input
94c
95c 1. bl(*) - storage for everything (common big in leit)
96c 2. nbls1 - reduced block-size
97c 3. l11   _  mmax total ang.mom. +1  (i+j+k+l+1)
98c 4. l12   _ lensm(mmax) total number of function up to mmax (see Iobara)
99c 5. mem2  - dimension for wt2(nbls1,mem2) used in Trac12, Trac34
100c
101c l11,l12,mem2 are setup in ithe Memo4a routine and give dimensions for
102c wt1(nbls1,l11,l12), wt2(nbls1,mem2) ;
103c
104c Output
105c
106c wt0(nbls1,lnij,lnkl) - containing (i+j,s|k+l,s) integrals
107c                        This is located in bl(*) from bl(iwt0)
108C-----------------------------------------------------------------------
109c decide how to shift angular momentum : from ij to kl or vice versa
110c (it depends on stability and l-shells) :
111c
112      call how_2_shift(stable(1),nsij,nskl,mmax,mmax1,immax1,
113     $     kmmax1,lobsa)
114c
115c output : mmax1,immax1,kmmax1 and lobsa showing how to shift
116C-----------------------------------------------------------------------
117c calculate (ss|ss)(m) integrals :
118c
119         call ssssm(nbls1,bl(irys),bl(iconst),mmax1,
120     *              bl(iwt1),l11,l12,bl(iwt0),lnij,lnkl)
121C-----------------------------------------------------------------------
122c now going by 10/20 or 30/40 depends upon numerical stability
123c                and presence of l-shells
124ctest: counter
125c
126      if(lobsa.le.2) then
127         xlobsa2=xlobsa2+nbls1
128      else
129         xlobsa4=xlobsa4+nbls1
130      endif
131C-----------------------------------------------------------------------
132c
133#if defined(XLF11) || defined(XLFLINUX)
134      if ( lobsa .eq. 1 ) goto 10
135      if ( lobsa .eq. 2 ) goto 20
136      if ( lobsa .eq. 3 ) goto 30
137      if ( lobsa .eq. 4 ) goto 40
138#else
139      go to (10,20,30,40) lobsa
140#endif
141c
142   10 continue
143c........(s,s|k+l,s)
144         call obasai_2(bl(irhocpd),bl(icdnix),bl(ixqnx),bl(ixwq),
145     *                  mmax1,kmmax1,bl(iwt1),l11,l12,nbls1)
146         call wt0wt1(bl(iwt0),lnij,lnkl,nbls1,bl(iwt1),l11,l12,
147     *               nsij,nskl,2)
148c
149   20 continue
150c........(i+j+k+l,s|s,s)
151         call obasai_2(bl(irhoapb),bl(iabnix),bl(ixpnx),bl(ixwp),
152     *                  mmax1,immax1,bl(iwt1),l11,l12,nbls1)
153         call wt0wt1(bl(iwt0),lnij,lnkl,nbls1,bl(iwt1),l11,l12,
154     *               nsij,nskl,1)
155       if(nskl.gt.1) then
156         call wt2wt1(bl(iwt2),mem2,nbls1,bl(iwt1),l11,l12,mmax1)
157            call trac12_2(bl(iwt0),lnij,lnkl,nbls1,bl(iwt2),mem2,
158     *                    bl(ip1234),bl(iabcd),bl(ihabcdx) )
159       endif
160ctry98
161cccccc   call rescale_wt0(nbls1,bl(iwt0),lnij*lnkl,bl(iconst))
162ctry98
163       return
164c
165   30 continue
166c........(i+j,s|s,s)
167         call obasai_2(bl(irhoapb),bl(iabnix),bl(ixpnx),bl(ixwp),
168     *                  mmax1,immax1,bl(iwt1),l11,l12,nbls1)
169         call wt0wt1(bl(iwt0),lnij,lnkl,nbls1,bl(iwt1),l11,l12,
170     *               nsij,nskl,1)
171   40 continue
172c........(s,s|i+j+k+l,s)
173         call obasai_2(bl(irhocpd),bl(icdnix),bl(ixqnx),bl(ixwq),
174     *                  mmax1,kmmax1,bl(iwt1),l11,l12,nbls1)
175         call wt0wt1(bl(iwt0),lnij,lnkl,nbls1,bl(iwt1),l11,l12,
176     *               nsij,nskl,2)
177       if(nsij.gt.1) then
178         call wt2wt1(bl(iwt2),mem2,nbls1,bl(iwt1),l11,l12,mmax1)
179            call trac34_2(bl(iwt0),lnij,lnkl,nbls1,bl(iwt2),mem2,
180     *                    bl(ip1234),bl(iabcd),bl(ihabcdx) )
181       endif
182c
183c---------------------------------------------------------------
184      end
185c=======================================================================
186      subroutine trobsa_1(bl,nbls1,l11,l12,mem2)
187C-----------------------------------------------------------------------
188c trobsa_1 is used for iroute=1 only
189C-----------------------------------------------------------------------
190      implicit real*8 (a-h,o-z)
191      logical stable
192      common /tracy_stability/ stable(5000) ! this is limit limblks
193c
194      common/obarai/
195     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,mmax,
196     * nqi,nqj,nqk,nql,nsij,nskl,
197     * nqij,nqij1,nsij1,nqkl,nqkl1,nskl1,ijbeg,klbeg
198c
199c     common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
200c
201      common /memor4/ iwt0,iwt1,iwt2,ibuf,ibuf2,
202     * ibfij1,ibfij2,ibfkl1,ibfkl2,
203     * ibf2l1,ibf2l2,ibf2l3,ibf2l4,ibfij3,ibfkl3,
204     * ibf3l,issss,
205     * ix2l1,ix2l2,ix2l3,ix2l4,ix3l1,ix3l2,ix3l3,ix3l4,
206     * ixij,iyij,izij, iwij,ivij,iuij,isij
207c
208      common /memor5b/ irppq,
209     * irho,irr1,irys,irhoapb,irhocpd,iconst,ixwp,ixwq,ip1234,
210     * idx1,idx2,indx
211c
212      common /memor5c/ itxab,itxcd,iabcd,ihabcd
213      common /memor5d/ iabnix,icdnix,ixpnx,ixqnx,ihabcdx
214c
215      dimension bl(*)
216C-----------------------------------------------------------------------
217c test : counter
218      common /lobsa_times/ xlobsa2,xlobsa4
219C-----------------------------------------------------------------------
220C-----------------------------------------------------------------------
221c PARAMETERS
222c-----------
223c Input
224c
225c 1. bl(*) - storage for everything (common big in leit)
226c 2. nbls1 - reduced block-size
227c 3. l11   _  mmax total ang.mom. +1  (i+j+k+l+1)
228c 4. l12   _ lensm(mmax) total number of function up to mmax (see Iobara)
229c 5. mem2  - dimension for wt2(nbls1,mem2) used in Trac12, Trac34
230c
231c l11,l12,mem2 are setup in ithe Memo4a routine and give dimensions for
232c wt1(nbls1,l11,l12), wt2(nbls1,mem2) ;
233c
234c Output
235c
236c wt0(nbls1,lnij,lnkl) - containing (i+j,s|k+l,s) integrals
237c                        This is located in bl(*) from bl(iwt0)
238C-----------------------------------------------------------------------
239c The block of NBLS1 primitive quartets has to be split into two parts
240c (1) contaning stable quartets, (2) contaning unstable quartets
241c These two could be calculated seperatly. To do so some data must be
242c copied and some additional memory must be allocated.
243c However, the number of unstable quartets is very small (1%) which makes
244c cheaper to calculated the whole block as stable and then recalculate
245c only unstable quartets. This way we can avoid copying data and diminish
246c an additional memory requiment.
247C-----------------------------------------------------------------------
248c Thus, calculate first all quartets in a block as stable ones.
249c it will go by lobsa=1,2 or 2
250c
251         call how_2_shift(.true.,nsij,nskl,mmax,mmax1,immax1,
252     $     kmmax1,lobsa)
253c
254c output : mmax1,immax1,kmmax1 and lobsa showing how to shift
255C-----------------------------------------------------------------------
256c calculate (ss|ss)(m) integrals :
257c
258         call ssssm(nbls1,bl(irys),bl(iconst),mmax1,
259     *              bl(iwt1),l11,l12,bl(iwt0),lnij,lnkl)
260c
261#if defined(XLF11) || defined(XLFLINUX)
262      if ( lobsa .eq. 1 ) goto 10
263      if ( lobsa .eq. 2 ) goto 20
264      if ( lobsa .eq. 3 ) goto 30
265      if ( lobsa .eq. 4 ) goto 40
266#else
267      go to (10,20,30,40) lobsa
268#endif
269c
270   10    continue
271c........(s,s|k+l,s)
272         call obasai_1(bl(irhocpd),bl(icdnix),bl(ixqnx),bl(ixwq),
273     *                  mmax1,kmmax1,bl(iwt1),l11,l12,nbls1)
274         call wt0wt1(bl(iwt0),lnij,lnkl,nbls1,bl(iwt1),l11,l12,
275     *               nsij,nskl,2)
276c
277   20    continue
278c........(i+j+k+l,s|s,s)
279         call obasai_1(bl(irhoapb),bl(iabnix),bl(ixpnx),bl(ixwp),
280     *                  mmax1,immax1,bl(iwt1),l11,l12,nbls1)
281         call wt0wt1(bl(iwt0),lnij,lnkl,nbls1,bl(iwt1),l11,l12,
282     *               nsij,nskl,1)
283       if(nskl.gt.1) then
284         call wt2wt1(bl(iwt2),mem2,nbls1,bl(iwt1),l11,l12,mmax1)
285            call trac12_1(bl(iwt0),lnij,lnkl,nbls1,bl(iwt2),mem2,
286     *                    bl(ip1234),bl(iabcd),bl(ihabcdx) )
287       endif
288C-----------------------------------------------------------------------
289c take care of unstable quartets :
290c
291      call memo1_int(nbls1,iq_u)
292      call find_unstable(stable,nbls1, bl(iq_u),nbls1_u)
293C-----------------------------------------------------------------------
294ctest: counter
295c
296         xlobsa2=xlobsa2+(nbls1-nbls1_u)
297         xlobsa4=xlobsa4+nbls1_u
298C-----------------------------------------------------------------------
299c
300c output: iq_u(nbls1_u) - list of unstable quartets in a block
301c
302      IF(nbls1_u.EQ.0) THEN
303         call retmem(1)
304         RETURN
305      ENDIF
306C-----------------------------------------------------------------------
307c to copy corresponding data for calculations of unstable quartets
308c memory is needed :
309c
310c rys,const & habcd,abcd,p1234 can be overwritten but
311c rhoapb, abnix, xpnx, xwp   and   rhocpd, cdnix, xqnx, xwq
312c and final integrals wt0() must have new allocations :
313c
314       mmax_1=mmax-1
315       call getmem(nbls1_u*lnij*lnkl, iwt0_u)
316       call getmem(nbls1_u, irhoapb_u )              ! rhoapb & rhocpd
317       call getmem(nbls1_u*(mmax_1), iabnix_u )      ! abnix  & cdnix
318       call getmem(nbls1_u*3, ixpnx_u )              ! xpnx   & xqnx
319       call getmem(nbls1_u*3, ixwp_u )               ! xwp    & xwq
320c-----------------------------------------------------------------------
321c
322         call how_2_shift(.false.,nsij,nskl,mmax,mmax1,immax1,
323     $      kmmax1,lobsa)
324c
325c it will go by lobsa=3,4 or 4
326c-----------------------------------------------------------------------
327c calculate (ss|ss)(m) integrals :
328c
329       call copy_unstable_1(nbls1,bl(iq_u),nbls1_u,bl(irys),bl(iconst))
330c
331       call ssssm(nbls1_u,bl(irys),bl(iconst),mmax1,
332     *              bl(iwt1),l11,l12,bl(iwt0_u),lnij,lnkl)
333c
334#if defined(XLF11) || defined(XLFLINUX)
335      if ( lobsa .eq. 1 ) goto 10
336      if ( lobsa .eq. 2 ) goto 20
337      if ( lobsa .eq. 3 ) goto 30
338      if ( lobsa .eq. 4 ) goto 40
339#else
340      go to (10,20,30,40) lobsa
341#endif
342c
343   30  continue
344c........(i+j,s|s,s)
345c
346       call copy_unstable_2(nbls1, bl(iq_u),nbls1_u, mmax_1,
347     *              bl(irhoapb  ),bl(iabnix  ),bl(ixpnx  ),bl(ixwp  ),
348     *              bl(irhoapb_u),bl(iabnix_u),bl(ixpnx_u),bl(ixwp_u))
349c
350       call obasai_1(bl(irhoapb_u),bl(iabnix_u),bl(ixpnx_u),bl(ixwp_u),
351     *               mmax1,immax1,bl(iwt1),l11,l12,nbls1_u)
352       call wt0wt1(bl(iwt0_u),lnij,lnkl,nbls1_u,bl(iwt1),l11,l12,
353     *               nsij,nskl,1)
354   40  continue
355c........(s,s|i+j+k+l,s)
356c
357       irhocpd_u= irhoapb_u
358       icdnix_u = iabnix_u
359       ixqnx_u  = ixpnx_u
360       ixwq_u   = ixwp_u
361c
362       call copy_unstable_2(nbls1, bl(iq_u),nbls1_u, mmax_1,
363     *              bl(irhocpd  ),bl(icdnix  ),bl(ixqnx  ),bl(ixwq  ),
364     *              bl(irhocpd_u),bl(icdnix_u),bl(ixqnx_u),bl(ixwq_u))
365c
366       call obasai_1(bl(irhocpd_u),bl(icdnix_u),bl(ixqnx_u),bl(ixwq_u),
367     *               mmax1,kmmax1,bl(iwt1),l11,l12,nbls1_u)
368       call wt0wt1(bl(iwt0_u),lnij,lnkl,nbls1_u,bl(iwt1),l11,l12,
369     *             nsij,nskl,2)
370      if(nsij.gt.1) then
371         call copy_unstable_3(nbls1, bl(iq_u),nbls1_u,mmax,
372     *                        bl(ip1234),bl(iabcd),bl(ihabcdx) )
373         call wt2wt1(bl(iwt2),mem2,nbls1_u,bl(iwt1),l11,l12,mmax1)
374         call trac34_1(bl(iwt0_u),lnij,lnkl,nbls1_u,bl(iwt2),mem2,
375     *                 bl(ip1234),bl(iabcd),bl(ihabcdx) )
376      endif
377C-----------------------------------------------------------------------
378c Replace unstable (incorrect) integrals in wt0(nbls,lnij,lnkl) by
379c correct ones from wt0_u(nbls1_u,lnij,lnkl) :
380c
381      call copy_unstable_4(nbls1, bl(iq_u),nbls1_u,nqij,nqkl,
382     *                     bl(iwt0),bl(iwt0_u),lnij,lnkl)
383C-----------------------------------------------------------------------
384c release memory allocated for unstable quartets
385c
386      call retmem(6)
387C-----------------------------------------------------------------------
388      end
389c=======================================================================
390      subroutine find_unstable(stable,nbls1, iq_u,nbls1_u)
391      logical stable(nbls1)
392      dimension iq_u(*)
393c
394      nbls1_u=0
395      do i=1,nbls1
396         if(.not.stable(i)) then
397            nbls1_u=nbls1_u+1
398            iq_u(nbls1_u)=i
399         endif
400      enddo
401c
402      end
403c=======================================================================
404      subroutine copy_unstable_1(nbls1,iq_u,nbls1_u,rys,const)
405      implicit real*8 (a-h,o-z)
406      dimension iq_u(*)
407      dimension rys(*),const(*)
408c
409      do i=1,nbls1_u
410         ijkl=iq_u(i)
411         rys(i)  =rys(ijkl)
412         const(i)=const(ijkl)
413      enddo
414c
415      end
416c=======================================================================
417      subroutine copy_unstable_2(nbls1, iq_u,nbls1_u, mmax_1,
418     *                            rhoapb  ,abnix  ,xpnx  ,xwp  ,
419     *                            rhoapb_u,abnix_u,xpnx_u,xwp_u )
420      implicit real*8 (a-h,o-z)
421      dimension iq_u(*)
422      dimension rhoapb(nbls1)        ,rhoapb_u(nbls1_u)
423      dimension  abnix(nbls1,mmax_1) , abnix_u(nbls1_u,mmax_1)
424      dimension   xpnx(nbls1,3)      ,  xpnx_u(nbls1_u,3)
425      dimension    xwp(nbls1,3)      ,   xwp_u(nbls1_u,3)
426c
427      do i=1,nbls1_u
428         ijkl=iq_u(i)
429         rhoapb_u(i)=rhoapb(ijkl)
430      enddo
431      do m=1,mmax_1
432         do i=1,nbls1_u
433            ijkl=iq_u(i)
434            abnix_u(i,m)=abnix(ijkl,m)
435         enddo
436      enddo
437      do icart=1,3
438         do i=1,nbls1_u
439            ijkl=iq_u(i)
440            xpnx_u(i,icart)= xpnx(ijkl,icart)
441             xwp_u(i,icart)=  xwp(ijkl,icart)
442         enddo
443      enddo
444c
445      end
446c=======================================================================
447      subroutine copy_unstable_3(nbls1, iq_u,nbls1_u,mmax,
448     *                           p1234,abcd,habcdx )
449      implicit real*8 (a-h,o-z)
450#include "texas_lpar.fh"
451      dimension iq_u(*)
452c     dimension p1234(nbls1,3),abcd(nbls1),habcdx(nbls1,3,nfu(mmax))
453      dimension p1234(*)      ,abcd(*)    ,habcdx(*)
454c
455      do i=1,nbls1_u
456         ijkl=iq_u(i)
457         abcd(i)=abcd(ijkl)
458      enddo
459c
460      do icart=1,3
461         icart_f=(icart-1)*nbls1
462         icart_u=(icart-1)*nbls1_u
463         do i=1,nbls1_u
464            ijkl=iq_u(i)
465            ij_f=icart_f +ijkl
466            ij_u=icart_u +i
467            p1234(ij_u)=p1234(ij_f)
468         enddo
469      enddo
470c
471      mdim=nfu(mmax)
472      do m=1,mdim
473         m1=(m-1)*3
474         do icart=1,3
475            jk_f=( m1+(icart-1) )*nbls1
476            jk_u=( m1+(icart-1) )*nbls1_u
477            do i=1,nbls1_u
478               ijkl=iq_u(i)
479               ijk_f=jk_f +ijkl
480               ijk_u=jk_u +i
481               habcdx(ijk_u)=habcdx(ijk_f)
482            enddo
483         enddo
484      enddo
485c
486      end
487c=======================================================================
488      subroutine copy_unstable_4(nbls1,iq_u,nbls1_u,nqij,nqkl,
489     *                           wt0,wt0_u, lnij,lnkl)
490      implicit real*8 (a-h,o-z)
491#include "texas_lpar.fh"
492c
493      dimension iq_u(*)
494      dimension wt0(nbls1,lnij,lnkl), wt0_u(nbls1_u,lnij,lnkl)
495c
496      do kl=nfu(nqkl)+1,lnkl
497         do ij=nfu(nqij)+1,lnij
498            do i=1,nbls1_u
499               ijkl=iq_u(i)
500               wt0(ijkl,ij,kl)=wt0_u(i,ij,kl)
501            enddo
502         enddo
503      enddo
504c
505      end
506c=======================================================================
507      subroutine ssssm(nbls,rysx,const,mmax,wt1,l11,l12,wt0,l01,l02)
508      implicit real*8 (a-h,o-z)
509      dimension const(*),rysx(*)
510      dimension wt0(nbls,l01,l02),wt1(nbls,l11,l12)
511      dimension f0m(0:30)
512c--------------------------------------------------------------------
513c This subroutine calculates (s,s|s,s)(m) integrals with m=1,MMAX
514c where MMAX is the total angular momentum (+1). These integrals are
515c needed in the Obara-Saika method.
516c
517c INPUT:
518c
519c NBLS     - reduced block-size /number of quartets of primitive shells/
520c RYSX(nbls) - (P-Q)**2*(a+b)*(c+d)/(a+b+c+d) an parameter for Fm
521c CONST(nbls)- see PREC4NEG and PRECAL2A
522c              CONST=PI3*SABCD/(PQ*SQRT(PPQ)) for all int.
523c MMAX     - total ang. mom. +1
524c l11,l12 - dimensions for wt1
525c l01,l02 - dimensions for wt0
526c
527c OUTPUT:
528c
529c wt1(ijkl,m,1) - integrals (s,s|s,s)(m=1,mmax)
530c wt0(ijkl,1,1) - integrals (s,s|s,s)(m=1)
531c
532c--------------------------------------------------------------------
533C  CALCULATE (SS,SS)(M)  M=1,MMAX
534C    S ORBITALS ARE NORMALIZED
535C
536         do 411 i=1,nbls
537         xrys=rysx(i)
538         call fm(xrys,mmax-1,f0m)
539         do 411 m=1,mmax
540         wt1(i,m,1)=const(i)*f0m(m-1)
541ccccc    write(6,*)'   m=',m,' ssssm=',wt1(i,m,1)
542 411     CONTINUE
543c
544         do 420 i=1,nbls
545         wt0(i,1,1)=wt1(i,1,1)
546  420    continue
547cxx   call dcopy(nbls,wt1(1,1,1),1,wt0(1,1,1),1)
548c
549      end
550c====================================================================
551      subroutine wt0wt1(wt0,l01,l02,nbls,wt1,l11,l12,nsij,nskl,lab)
552      implicit real*8 (a-h,o-z)
553#include "texas_lpar.fh"
554      dimension wt0(nbls,l01,l02),wt1(nbls,l11,l12)
555c
556#ifdef XLF11
557      if ( lab .eq. 1 ) goto 10
558      if ( lab .eq. 2 ) goto 20
559#else
560      go to (10,20) lab
561#endif
562c
563   10 continue
564      do 100 inp=2,nfu(nsij +1)
565         do 100 i=1,nbls
566         wt0(i,inp,1)=wt1(i,1,inp)
567  100 continue
568c
569      return
570c
571   20 continue
572      do 200 knp=2,nfu(nskl +1)
573         do 200 i=1,nbls
574         wt0(i,1,knp)=wt1(i,1,knp)
575  200 continue
576c
577      end
578c=================================================================
579c obasai subroutines :
580c----------------------------------------------------------------
581c PARAMETERS
582c-----------
583c
584c Input :
585c 1. RHOAPB(nbls1) -  (c+d)/(a+b+c+d) - exponents
586c 2. ABNIA(nbls1,*) -   L*( 0.5/(a+b) )  with L=1,2,...MMAX-1
587c 3. XPA(nbls1,3)   -  (P-A) coordinates
588c 4. XWP(nbls1,3)   -  (W-P) coordinates
589c 5. MMAX         - total ang. mom. +1
590c 6. IMMAX        - see Trobsa
591c 7. XT(nbls1,l11,l12) - integrals (s,s|s,s)(m)
592c 8. NBLS - reduced block-size (nbls1)
593c
594c Output
595c
596c 1. XT(nbls1,l11,l12) - integrals (i+j+k+l,s|s,s)(m)
597c                        only integrals with m=1 are used later
598c----------------------------------------------------------------
599c     This is the OBARA-SAIKA recursive method to generate integrals
600c  with all angular momenta placed in the position 1 /(i+j+k+l,s|s,s)/
601c  or in the position 2 /(s,s|i+j+k+l,s)/. This subroutine is called
602c  for both cases with different set of parameters from OBSA1-OBSA4
603c  routines. It is accesable from any place.
604c     As a input, integrals (s,s|s,s)(m) from SSSSM routine are used
605c  /in xt(nbls,mmax,1)/. Integrals (i+j+k+l,s|s,s)(m) return also in
606c  the array xt(nbls,mmax,l12) with l12=nfu(mmax+1)
607c----------------------------------------------------------------
608c========================
609      subroutine obasai_1(rhoapb,abnia,xpa,xwp,mmax,immax,xt,
610     *                    l11,l12,nbls)
611      implicit real*8 (a-h,o-z)
612#include "texas_lpar.fh"
613      dimension xt(nbls,l11,l12)
614      dimension abnia(nbls,*),rhoapb(*)
615      dimension xpa(nbls,3),xwp(nbls,3)
616c---------------------------------------
617       MMM=MMAX-1
618       do 100 inp=2,4
619       in0=ifrst(inp)
620       icr=icoor(inp)
621          call recur1(nbls,l11,mmm,xt(1,1,inp),xt(1,1,in0),
622     *                xpa(1,icr),xwp(1,icr))
623  100  continue
624c
625       mmm=mmax-2
626c
627       do 105 im=1,immax
628          do 110 inm=nfu(im)+1,nfu(im+1)
629          icrm=icoor(inm)
630             do 115 ixyz=icrm,3
631             in0=npxyz(ixyz,inm)
632                 do 120 jxyz=ixyz,3
633                 inp=npxyz(jxyz,in0)
634                 call recur1(nbls,l11,mmm,xt(1,1,inp),xt(1,1,in0),
635     *                       xpa(1,jxyz),xwp(1,jxyz))
636  120            continue
637             inpa=npxyz(ixyz,in0)
638             nia0= nia(ixyz,in0)
639             call recur2_1(nbls,l11,mmm,
640     *                   xt(1,1,inpa),xt(1,1,inm),abnia(1,nia0),rhoapb)
641  115        continue
642  110     continue
643       mmm=mmm-1
644  105  continue
645c
646      end
647c========================
648      subroutine obasai_2(rhoapb,abnia,xpa,xwp,mmax,immax,xt,
649     *                    l11,l12,nbls)
650      implicit real*8 (a-h,o-z)
651#include "texas_lpar.fh"
652c
653      dimension xt(nbls,l11,l12)
654      dimension xpa(nbls,3),xwp(nbls,3)
655ccccc dimension abnia(nbls,*),rhoapb(*)
656      dimension abnia(*)
657c--------------------------------------
658       MMM=MMAX-1
659       do 100 inp=2,4
660       in0=ifrst(inp)
661       icr=icoor(inp)
662          call recur1(nbls,l11,mmm,xt(1,1,inp),xt(1,1,in0),
663     *                xpa(1,icr),xwp(1,icr))
664  100  continue
665c
666       mmm=mmax-2
667c
668       do 105 im=1,immax
669          do 110 inm=nfu(im)+1,nfu(im+1)
670          icrm=icoor(inm)
671             do 115 ixyz=icrm,3
672             in0=npxyz(ixyz,inm)
673                 do 120 jxyz=ixyz,3
674                 inp=npxyz(jxyz,in0)
675                 call recur1(nbls,l11,mmm,xt(1,1,inp),xt(1,1,in0),
676     *                       xpa(1,jxyz),xwp(1,jxyz))
677  120            continue
678             inpa=npxyz(ixyz,in0)
679             nia0= nia(ixyz,in0)
680             call recur2_2(nbls,l11,mmm,
681     *                   xt(1,1,inpa),xt(1,1,inm),abnia(nia0),rhoapb)
682  115        continue
683  110     continue
684       mmm=mmm-1
685  105  continue
686c
687      end
688c=================================================================
689      subroutine recur1(nbls,l11,mmm,xtp,xt0,xpa,xwp)
690      implicit real*8 (a-h,o-z)
691      dimension xtp(nbls,l11),xt0(nbls,l11),xpa(nbls),xwp(nbls)
692c
693  150 continue
694          do 1501 m=1,mmm
695          m1=m+1
696          do 1501 i=1,nbls
697          xtp(i,m)=xpa(i)*xt0(i,m)+xwp(i)*xt0(i,m1)
698 1501     continue
699c
700      end
701c=================================================================
702c recur2 routines :
703c
704c========================
705      subroutine recur2_1(nbls,l11,mmm,xtp,xtm,abnia,rhoapb)
706      implicit real*8 (a-h,o-z)
707      dimension xtp(nbls,l11),xtm(nbls,l11),abnia(nbls)
708      dimension rhoapb(*)
709c
710           do 1501 m=1,mmm
711           m1=m+1
712           do 1501 i=1,nbls
713           xtp(i,m)=xtp(i,m)+abnia(i)*(xtm(i,m) - xtm(i,m1)*rhoapb(i))
714 1501      continue
715c
716      end
717c========================
718      subroutine recur2_2(nbls,l11,mmm,xtp,xtm,abnia,rhoapb)
719      implicit real*8 (a-h,o-z)
720      dimension xtp(nbls,l11),xtm(nbls,l11)
721c
722           do 1501 m=1,mmm
723           m1=m+1
724           do 1501 i=1,nbls
725           xtp(i,m)=xtp(i,m)+abnia   *(xtm(i,m) - xtm(i,m1)*rhoapb   )
726 1501      continue
727      end
728c=================================================================
729      subroutine wt2wt1(wt2,mem2,nbls,wt1,l11,l12,mmax)
730      implicit real*8 (a-h,o-z)
731#include "texas_lpar.fh"
732      dimension wt1(nbls,l11,l12),wt2(nbls,mem2)
733c
734       do 145 inp=1,nfu(mmax+1)
735          do 145 i=1,nbls
736          wt2(i,inp)=wt1(i,1,inp)
737  145  continue
738c
739      end
740c=================================================================
741C*****************************************************************
742c trac12 routines :
743c
744C  This subroutine performs calculations for a block (nbls)
745C  of quartets of primitive shells. The integrals calculated
746C  here are of the type :
747C               (i+j,s | k+l,s)
748C  This subroutine is called from TROBSA.
749C
750C     Calculations are performed according to the Tracy's
751C  recursive formula. It is made in the loop over an angular
752C  momentum increasing on the the center no. 3. For each such
753C  a recursive step the tracij routine is called with the wt0 and
754C  the wt2 matrix at three different locations. The wt2 matrix is
755C  two-domensional here but it is three-dim. in tracij. This
756C  to execute Tracy's recursives. At the very begining wt2
757C  contains wt1(1,nbls,l12)= (i+j+k+l,s|s,s) (m=1) integrals.
758C  Final intgrals (i+j,s|k+l,s) return from tracij in the wt0 matrix
759C  and nothing else is done with them here. They go back through
760C  the routine TROBSA from which OBSAIJ (OBSAKL) was called to
761C  the routine TWOE (where Trobsa is called) and then are
762C  contracted in the routine ASSEMBLE.
763C
764C  INPUT
765C  -------
766C
767C  information from the obarai common block
768C
769C   and precalculated quantities :
770C
771C  p1234 (nbls,3) - geometry stuff
772C
773C  OUTPUT
774C  -------
775C  wt0(nbls,l01,l02) - contains final (i+j,s|k+l,s) integrals
776C
777C  Locally in use - the wt2(nbls,mem2) matrix for tracij
778C
779C*****************************************************************
780c========================
781      subroutine trac12_1(wt0,l01,l02,nbls,wt2,mem2, p1234, abcd,habcd)
782      IMPLICIT REAL*8 (A-H,O-Z)
783cnmr
784      character*11 scftype
785      character*8 where
786      common /runtype/ scftype,where
787c
788      common /tracy/ kbeg,kend,i0b,i0e,kp
789      common/obarai/
790     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
791     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
792     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
793#include "texas_lpar.fh"
794C
795      dimension wt0(nbls,l01,l02),wt2(nbls,mem2)
796      dimension p1234(nbls,3)
797      dimension abcd(nbls),habcd(nbls,3,*)
798c---------------------------------------------------------------
799c   calculate (i+j,s|k+l,s) integrals from (i+j+k+l,s|s,s) ones
800c            according to the Tracy's recursive formula
801c   All target classes needed for further shifts (A->B) are
802c   constructed.
803c     target classes :  from (i,s|k,s) to (i+j,s|k+l,s)
804C
805c   The target classes appear in last nskl-nqkl+1 recursive steps
806c   total number of recursive steps is nrs=nskl-2+1
807c---------------------------------------------------------------
808cderivatives:
809c
810      mmax1=mmax
811      if(where.eq.'shif'.or. where.eq.'forc') then
812          mmax1=mmax-1
813      endif
814      if(where.eq.'hess') then
815          mmax1=mmax-2
816      endif
817c
818c---------------------------------------------------------------
819c
820c  addressing in the wt2 matrix for recurcive in Tracy
821c
822      ia3=0
823      k31=nfu(mmax+1)
824c-?-- k31=nfu(mmax1+1)
825      k32=1
826      ia2=ia3
827      k21=k31
828      k22=k32
829c
830      do 2000 kp=2,nskl
831      kbeg=nfu(kp)+1
832      kend=nfu(kp+1)
833c
834      i0b=mmax1+1-kp
835      i0e=nqij-nskl+kp
836      if(i0e.le.0) i0e=1
837c
838      ia1=ia2+k21*k22
839      k11=nfu(i0b+1)
840      k12=nfu(kp+1)
841c
842      i11=ia1+1
843      i21=ia2+1
844      i31=ia3+1
845      call tracij_1(wt2(1,i11),k11,k12,
846     *              wt2(1,i21),k21,k22,
847     *              wt2(1,i31),k31,k32,
848     *              p1234,wt0,l01,l02,nbls,
849     *              abcd,habcd)
850      ia3=ia2
851      ia2=ia1
852      k31=k21
853      k32=k22
854      k21=k11
855      k22=k12
856 2000 continue
857c
858      END
859c========================
860      subroutine trac12_2(wt0,l01,l02,nbls,wt2,mem2, p1234, abcd,habcd)
861      IMPLICIT REAL*8 (A-H,O-Z)
862cnmr
863      character*11 scftype
864      character*8 where
865      common /runtype/ scftype,where
866cnmr
867      common /tracy/ kbeg,kend,i0b,i0e,kp
868c
869      common/obarai/
870     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
871     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
872     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
873#include "texas_lpar.fh"
874      dimension wt0(nbls,l01,l02),wt2(nbls,mem2)
875      dimension p1234(nbls,3)
876ccccc dimension abcd(nbls),habcd(nbls,3,*)
877      dimension habcd(3,*)
878c-------------------------------------------------------
879c   calculate (i+j,s|k+l,s) integrals from (i+j+k+l,s|s,s) ones
880c            according to the Tracy's recursive formula
881c   All target classes needed for further shifts (A->B) are
882c   constructed.
883c     target classes :  from (i,s|k,s) to (i+j,s|k+l,s)
884C
885c   The target classes appear in last nskl-nqkl+1 recursive steps
886c   total number of recursive steps is nrs=nskl-2+1
887c-------------------------------------------------------
888cderivatives
889c
890      mmax1=mmax
891      if(where.eq.'shif'.or. where.eq.'forc') then
892          mmax1=mmax-1
893      endif
894      if(where.eq.'hess') then
895          mmax1=mmax-2
896      endif
897c
898c  addressing in the wt2 matrix for recurcive in Tracy
899c
900      ia3=0
901      k31=nfu(mmax+1)
902c-?-- k31=nfu(mmax1+1)
903      k32=1
904      ia2=ia3
905      k21=k31
906      k22=k32
907c
908      do 2000 kp=2,nskl
909      kbeg=nfu(kp)+1
910      kend=nfu(kp+1)
911c
912c--->  i0b=mmax+1-kp
913c--->  i0e=nqij-nskl+kp
914       i0b=mmax1+1-kp
915       i0e=nqij-nskl+kp
916       if(i0e.le.0) i0e=1
917c
918      ia1=ia2+k21*k22
919      k11=nfu(i0b+1)
920      k12=nfu(kp+1)
921c
922      i11=ia1+1
923      i21=ia2+1
924      i31=ia3+1
925c
926      call tracij_2(wt2(1,i11),k11,k12,
927     *              wt2(1,i21),k21,k22,
928     *              wt2(1,i31),k31,k32,
929     *              p1234,wt0,l01,l02,nbls,
930     *              abcd,habcd)
931      ia3=ia2
932      ia2=ia1
933      k31=k21
934      k32=k22
935      k21=k11
936      k22=k12
937 2000 continue
938c
939      END
940c=================================================================
941c trac34 routines :
942C*****************************************************************
943C  This subroutine performs calculations for a block (nbls)
944C  of quartets of primitive shells. The integrals calculated
945C  here are of the type :
946C               (i+j,s | k+l,s)
947C  in the case when the angular momentum i+j .LT. k+l
948C  This subroutine is called from TROBSA.
949C  see description in TRAC12.
950C*****************************************************************
951c========================
952      subroutine trac34_1(wt0,l01,l02,nbls,wt2,mem2, p1234, abcd,habcd)
953      IMPLICIT REAL*8 (A-H,O-Z)
954      character*11 scftype
955      character*8 where
956      common /runtype/ scftype,where
957c
958      common /tracy/ ibeg,iend,k0b,k0e,ip
959      common/obarai/
960     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
961     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
962     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
963#include "texas_lpar.fh"
964C
965      dimension wt0(nbls,l01,l02),wt2(nbls,mem2)
966      dimension p1234(nbls,3)
967      dimension abcd(nbls),habcd(nbls,3,*)
968c-----------------------------------------------------
969cderivatives:
970c
971      mmax1=mmax
972      if(where.eq.'shif'.or. where.eq.'forc') then
973          mmax1=mmax-1
974      endif
975      if(where.eq.'hess') then
976          mmax1=mmax-2
977      endif
978c
979c-----------------------------------------------------
980c  addressing in the wt2 matrix for recursive in Tracy
981c
982      ia3=0
983      k31=1
984      k32=nfu(mmax+1)
985      ia2=ia3
986      k21=k31
987      k22=k32
988c
989      do 2000 ip=2,nsij
990      ibeg=nfu(ip)+1
991      iend=nfu(ip+1)
992cccc
993c98    k0b=mmax+1-ip
994c98    k0e=nqkl-nsij+ip
995c98
996       k0b=mmax1+1-ip
997       k0e=nqkl-nsij+ip
998       if(k0e.le.0) k0e=1
999c
1000       ia1=ia2+k21*k22
1001       k11=nfu(ip+1)
1002       k12=nfu(k0b+1)
1003c
1004      i11=ia1+1
1005      i21=ia2+1
1006      i31=ia3+1
1007      call trackl_1(wt2(1,i11),k11,k12,
1008     *              wt2(1,i21),k21,k22,
1009     *              wt2(1,i31),k31,k32,
1010     *              p1234,wt0,l01,l02,nbls,
1011     *              abcd,habcd)
1012      ia3=ia2
1013      ia2=ia1
1014      k31=k21
1015      k32=k22
1016      k21=k11
1017      k22=k12
1018 2000 continue
1019c
1020      END
1021c========================
1022      subroutine trac34_2(wt0,l01,l02,nbls,wt2,mem2, p1234, abcd,habcd)
1023      IMPLICIT REAL*8 (A-H,O-Z)
1024      character*11 scftype
1025      character*8 where
1026      common /runtype/ scftype,where
1027c
1028      common /tracy/ ibeg,iend,k0b,k0e,ip
1029      common/obarai/
1030     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
1031     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
1032     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
1033#include "texas_lpar.fh"
1034C
1035      dimension wt0(nbls,l01,l02),wt2(nbls,mem2)
1036      dimension p1234(nbls,3)
1037ccccc dimension abcd(nbls),habcd(nbls,3,*)
1038      dimension habcd(3,*)
1039c-----------------------------------------------------
1040cderivatives :
1041c
1042      mmax1=mmax
1043      if(where.eq.'shif'.or. where.eq.'forc') then
1044          mmax1=mmax-1
1045      endif
1046      if(where.eq.'hess') then
1047          mmax1=mmax-2
1048      endif
1049c
1050c-----------------------------------------------------
1051c  addressing in the wt2 matrix for recursive in Tracy
1052c
1053      ia3=0
1054      k31=1
1055      k32=nfu(mmax+1)
1056      ia2=ia3
1057      k21=k31
1058      k22=k32
1059c
1060      do 2000 ip=2,nsij
1061      ibeg=nfu(ip)+1
1062      iend=nfu(ip+1)
1063c
1064c98    k0b=mmax+1-ip
1065c98    k0e=nqkl-nsij+ip
1066c98
1067       k0b=mmax1+1-ip
1068       k0e=nqkl-nsij+ip
1069       if(k0e.le.0) k0e=1
1070c
1071       ia1=ia2+k21*k22
1072       k11=nfu(ip+1)
1073       k12=nfu(k0b+1)
1074c
1075      i11=ia1+1
1076      i21=ia2+1
1077      i31=ia3+1
1078      call trackl_2(wt2(1,i11),k11,k12,
1079     *              wt2(1,i21),k21,k22,
1080     *              wt2(1,i31),k31,k32,
1081     *              p1234,wt0,l01,l02,nbls,
1082     *              abcd,habcd)
1083      ia3=ia2
1084      ia2=ia1
1085      k31=k21
1086      k32=k22
1087      k21=k11
1088      k22=k12
1089 2000 continue
1090c
1091      END
1092c=================================================================
1093c tracij routines :
1094c
1095C*****************************************************************
1096C  This subroutine shifts BY ONE an angular momentum from
1097C  position 1 to 3 according to the recursive formula of Tracy.
1098C  This is called ones from OBSAIJ for every new type of
1099C  orbitals constructed on the center 3 i.e.
1100C  (i+j+k+l-1,s|p,s),then (= -2|d,s),(= -3,s|fs),..(i+j,s|k+l,s)
1101C  The result is stored in the xt1 matrix and then used in next
1102C  step as a xt2, xt2 as xt3 and new xt1 is calculated. At the
1103C  end the final integrals are stored in the xt0 matrix. It is
1104C  done since the current recursive step (nrec) approaches nqkl
1105C  which is the first target class.
1106C
1107C
1108C  INPUT
1109C  ------
1110C  xt2, xt3 matrices - contain informations from previous
1111C                      recursive step
1112C
1113C  OUTPUT
1114C  -------
1115C  xt1 and xt0 matrices - xt0 contains at the end integrals
1116C                         of the type (i+j,s|k+l,s)
1117C  Comments
1118C ----------
1119C  p1234(nbls,3) - geometry parameters
1120C
1121C
1122C  Other important informations needed here are sent by
1123C         common /tracy/ kbeg,kend,i0b,i0e,nrec
1124C  They are set up in the calling place in trac12_1
1125C*****************************************************************
1126c========================
1127      subroutine tracij_1(xt1,l1b,l1e,xt2,l2b,l2e,xt3,l3b,l3e,
1128     *                    p1234,xt0,l01,l02,nbls,abcd,habcd)
1129      IMPLICIT REAL*8 (A-H,O-Z)
1130      COMMON /NUMBER/ ZERO,HALF,ONE,TWO,THREE,FOUR,FIVE,TEN,TEN6,TENM8,P
1131     1I,ACC
1132      common /flops/ iflop(20)
1133      common/obarai/
1134     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
1135     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
1136     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
1137#include "texas_lpar.fh"
1138      common /tracy/ kbeg,kend,i0b,i0e,nrec
1139c
1140      dimension xt0(nbls,l01,l02),xt1(nbls,l1b,l1e),
1141     *                            xt2(nbls,l2b,l2e),
1142     *                            xt3(nbls,l3b,l3e)
1143      dimension p1234(nbls,3)
1144      dimension abcd(nbls),habcd(nbls,3,*)
1145c------------------------------------------------------
1146          do 2005 knp=kbeg,kend
1147          kn0=ifrst(knp)
1148          kcr=icoor(knp)
1149          knm=nmxyz(kcr,kn0)
1150c
1151          if(knm.gt.0) then
1152             do in0=nfu(i0e)+1,nfu(i0b+1)
1153                inp=npxyz(kcr,in0)
1154                inm=nmxyz(kcr,in0)
1155      if(inm.gt.0) then
1156                do  i=1,nbls
1157       xt1(i,in0,knp)=p1234(i,kcr)*xt2(i,in0,kn0)-abcd(i)*xt2(i,inp,kn0)
1158     +  +habcd(i,kcr,kn0)*xt3(i,in0,knm)
1159     +  +habcd(i,kcr,in0)*xt2(i,inm,kn0)
1160                enddo
1161             else
1162                do  i=1,nbls
1163       xt1(i,in0,knp)=p1234(i,kcr)*xt2(i,in0,kn0)-abcd(i)*xt2(i,inp,kn0)
1164     +  +habcd(i,kcr,kn0)*xt3(i,in0,knm)
1165                   enddo
1166                endif
1167c
1168      enddo
1169          else
1170             do 1000 in0=nfu(i0e)+1,nfu(i0b+1)
1171             inp=npxyz(kcr,in0)
1172             inm=nmxyz(kcr,in0)
1173      if(inm.gt.0) then
1174                do 1002 i=1,nbls
1175       xt1(i,in0,knp)=p1234(i,kcr)*xt2(i,in0,kn0)-abcd(i)*xt2(i,inp,kn0)
1176     +                  +habcd(i,kcr,in0)*xt2(i,inm,kn0)
1177 1002           continue
1178             else
1179                do 1001 i=1,nbls
1180       xt1(i,in0,knp)=p1234(i,kcr)*xt2(i,in0,kn0)-abcd(i)*xt2(i,inp,kn0)
1181 1001           continue
1182             endif
1183c
1184 1000        continue
1185          endif
1186 2005     continue
1187c
1188      if(nrec.ge.nqkl) then
1189        do 150 knp=kbeg,kend
1190        do 150 in0=nfu(nqij)+1,nfu(nsij+1)
1191           do 150 i=1,nbls
1192           xt0(i,in0,knp)=xt1(i,in0,knp)
1193  150   continue
1194      endif
1195c
1196      end
1197c========================
1198      subroutine tracij_2(xt1,l1b,l1e,xt2,l2b,l2e,xt3,l3b,l3e,
1199     *                    p1234,xt0,l01,l02,nbls,abcd,habcd)
1200      implicit real*8 (a-h,o-z)
1201      common/obarai/
1202     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
1203     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
1204     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
1205#include "texas_lpar.fh"
1206      common /tracy/ kbeg,kend,i0b,i0e,nrec
1207c
1208      dimension xt0(nbls,l01,l02),xt1(nbls,l1b,l1e),
1209     *                            xt2(nbls,l2b,l2e),
1210     *                            xt3(nbls,l3b,l3e)
1211      dimension p1234(nbls,3)
1212cccc  dimension abcd(nbls),habcd(nbls,3,*)
1213      dimension habcd(3,*)
1214c------------------------------------------------------
1215c establish beginning & end for the loops over in0 & inm:
1216c
1217      if(i0e.gt.1) then
1218         inm_beg=nfu(i0e-1)+1
1219      else
1220         inm_beg= 1
1221      endif
1222c
1223      inm_end=nfu(i0b  )
1224c
1225      in0_beg=nfu(i0e)+1
1226      in0_end=nfu(i0b+1)
1227c------------------------------------------------------
1228      do knp=kbeg,kend
1229         kn0=ifrst(knp)
1230         kcr=icoor(knp)
1231         knm=nmxyz(kcr,kn0)
1232         habcdkn0=habcd(kcr,kn0)
1233c
1234c do in0 & inp, do not do inm :
1235c
1236         if(knm.gt.0) then
1237c.ok..     do in0=nfu(i0b+1),nfu(i0e)+1,-1    ! this is ok too
1238c.ok..      do in0=nfu(i0e)+1,nfu(i0b+1)       ! ok
1239            do in0=in0_beg,in0_end
1240               inp=npxyz(kcr,in0)
1241               do i=1,nbls
1242                  xt1(i,in0,knp)=p1234(i,kcr)*xt2(i,in0,kn0)
1243     *                               -abcd   *xt2(i,inp,kn0)
1244     *                              +habcdkn0*xt3(i,in0,knm)
1245               enddo
1246            enddo
1247         else
1248c.ok.       do in0=nfu(i0e)+1,nfu(i0b+1)
1249            do in0=in0_beg,in0_end
1250               inp=npxyz(kcr,in0)
1251               do i=1,nbls
1252                  xt1(i,in0,knp)=p1234(i,kcr)*xt2(i,in0,kn0)
1253     *                               -abcd   *xt2(i,inp,kn0)
1254               enddo
1255            enddo
1256         endif      !      if(knm.gt.0) then
1257c
1258c do inm only :
1259c
1260c.ok.    do inm=inm_end,inm_beg,-1    ! this is ok too
1261         do inm=inm_beg,inm_end
1262            in0=npxyz(kcr,inm)
1263            habcdin0=habcd(kcr,in0)
1264            do i=1,nbls
1265               xt1(i,in0,knp)=xt1(i,in0,knp)+habcdin0*xt2(i,inm,kn0)
1266            enddo
1267         enddo
1268      enddo
1269c------------------------------------------------------
1270c------------ saving  target classes ------------------
1271      if(nrec.ge.nqkl) then
1272         do knp=kbeg,kend
1273            do in0=nfu(nqij)+1,nfu(nsij+1)
1274               do i=1,nbls
1275                  xt0(i,in0,knp)=xt1(i,in0,knp)
1276               enddo
1277            enddo
1278         enddo
1279      endif
1280c---------------- target classes ----------------------
1281c
1282      end
1283c=================================================================
1284c trackl routines :
1285c
1286C
1287C  This subroutine shifts BY ONE an angular momentum from
1288C  position 3 to 1 according to the recursive formula of Tracy.
1289C  This is called ones from OBSAKL for every new type of
1290C  orbitals constructed on the center 1 i.e.
1291C  (p,s|i+j+k+l-1,s),then (d,s|=-2),...(i+j,s|k+l,s)
1292C  see description in the tracij subroutine.
1293C
1294C 1998 : the trackl_1 _2 routines are NOT analogues to tracij_1 &_2
1295C        Now both tracij_ & trackl_ have the same P1234() factors
1296C        regardless of nsij & nskl relation (it used to be different)
1297C        These p1234() factors are calculated in xwpq_ routines
1298C        Because of having the same factors in trackl_ there is
1299C        one more step involving multiplication of final
1300C        xt1(i,inp,kn0) integrals by abcd() .
1301c=================================================================
1302c
1303      subroutine trackl_1(xt1,l1b,l1e,xt2,l2b,l2e,xt3,l3b,l3e,
1304     *                    p1234,xt0,l01,l02,nbls, abcd,habcd)
1305      IMPLICIT REAL*8 (A-H,O-Z)
1306      COMMON /NUMBER/ ZERO,HALF,ONE,TWO,THREE,FOUR,FIVE,TEN,TEN6,TENM8,P
1307     1I,ACC
1308      common /flops/ iflop(20)
1309      common/obarai/
1310     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
1311     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
1312     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
1313#include "texas_lpar.fh"
1314      common /tracy/ ibeg,iend,k0b,k0e,nrec
1315c
1316      dimension xt0(nbls,l01,l02),xt1(nbls,l1b,l1e),
1317     *                            xt2(nbls,l2b,l2e),
1318     *                            xt3(nbls,l3b,l3e)
1319      dimension p1234(nbls,3)
1320      dimension abcd(nbls),habcd(nbls,3,*)
1321c------------------------------------------------------
1322c
1323          do 2005 inp=ibeg,iend
1324          in0=ifrst(inp)
1325          icr=icoor(inp)
1326          inm=nmxyz(icr,in0)
1327c
1328      if(inm.gt.0) then
1329        do  kn0=nfu(k0e)+1,nfu(k0b+1)
1330          knp=npxyz(icr,kn0)
1331          knm=nmxyz(icr,kn0)
1332          if(knm.gt.0) then
1333            do  i=1,nbls
1334      xt1(i,inp,kn0)=(p1234(i,icr)*xt2(i,in0,kn0)-xt2(i,in0,knp)
1335     + +habcd(i,icr,in0)*xt3(i,inm,kn0)
1336     + +habcd(i,icr,kn0)*xt2(i,in0,knm))*abcd(i)
1337            enddo
1338          else
1339            do  i=1,nbls
1340      xt1(i,inp,kn0)=(p1234(i,icr)*xt2(i,in0,kn0)-xt2(i,in0,knp)
1341     + +habcd(i,icr,in0)*xt3(i,inm,kn0))*abcd(i)
1342            enddo
1343          endif
1344c
1345c
1346        enddo
1347      else
1348              do 1000 kn0=nfu(k0e)+1,nfu(k0b+1)
1349              knp=npxyz(icr,kn0)
1350              knm=nmxyz(icr,kn0)
1351      if(knm.gt.0) then
1352                 do 1001 i=1,nbls
1353      xt1(i,inp,kn0)=(p1234(i,icr)*xt2(i,in0,kn0)-xt2(i,in0,knp)+
1354     + habcd(i,icr,kn0)*xt2(i,in0,knm))*abcd(i)
1355 1001            continue
1356              else
1357                 do i=1,nbls
1358      xt1(i,inp,kn0)=(p1234(i,icr)*xt2(i,in0,kn0)-xt2(i,in0,knp))*
1359     * abcd(i)
1360                 enddo
1361      endif
1362c
1363 1000        continue
1364      endif
1365 2005     continue
1366c
1367      if(nrec.ge.nqij) then
1368        do 150 kn0=nfu(nqkl)+1,nfu(nskl+1)
1369        do 150 inp=ibeg,iend
1370           do 150 i=1,nbls
1371           xt0(i,inp,kn0)=xt1(i,inp,kn0)
1372  150   continue
1373      endif
1374c
1375      end
1376c=========================
1377      subroutine trackl_2(xt1,l1b,l1e,xt2,l2b,l2e,xt3,l3b,l3e,
1378     *                    p1234,xt0,l01,l02,nbls, abcd,habcd)
1379      IMPLICIT REAL*8 (A-H,O-Z)
1380      common/obarai/
1381     * lni,lnj,lnk,lnl,lnij,lnkl,lnijkl,MMAX,
1382     * NQI,NQJ,NQK,NQL,NSIJ,NSKL,
1383     * NQIJ,NQIJ1,NSIJ1,NQKL,NQKL1,NSKL1,ijbex,klbex
1384#include "texas_lpar.fh"
1385      common /tracy/ ibeg,iend,k0b,k0e,nrec
1386c
1387      dimension xt0(nbls,l01,l02),xt1(nbls,l1b,l1e),
1388     *                            xt2(nbls,l2b,l2e),
1389     *                            xt3(nbls,l3b,l3e)
1390      dimension p1234(nbls,3)
1391cccc  dimension abcd(nbls),habcd(nbls,3,*)
1392      dimension habcd(3,*)
1393c------------------------------------------------------
1394c Note : shifting from 1 to 3 goes like this:
1395c
1396c    x1(i0,kp)=      [-b*AB-d*CD]/(c+d)*x2(i0,k0)
1397c              -             (a+b)/c+d)*x2(ip,k0)
1398c              + 0.5*nia(coor,i0)/(c+d)*x2(im,k0)
1399c              + 0.5*nia(coor,k0)/(c+d)*x2(i0,km)
1400c
1401c and corresponding shifting from 3 to 1 :
1402c
1403c    x1(ip,k0)=      [-b*AB-d*CD]/(a+b)*x2(i0,k0)
1404c              -             (c+d)/a+b)*x2(i0,kp)
1405c              + 0.5*nia(coor,i0)/(a+b)*x2(im,k0)
1406c              + 0.5*nia(coor,k0)/(a+b)*x2(i0,km)
1407c
1408c Taking factor (c+d)/(a+b) at the front yields
1409c the expression for x1(ip,kn0)
1410c
1411c    x1(ip,k0)=(c+d)/(a+b)
1412c                   *{[-b*AB-d*CD]/(c+d)*x2(i0,k0)
1413c               -                        x2(i0,kp)
1414c               + 0.5*nia(coor,i0)/(c+d)*x2(im,k0)
1415c               + 0.5*nia(coor,k0)/(c+d)*x2(i0,km) }
1416c
1417c which has THE SAME internal multiplicative factors
1418c as the one for x1(i0,kp) .
1419c------------------------------------------------------
1420c establish beginning & end for the loops over kn0 & knm:
1421c
1422      if(k0e.gt.1) then
1423         knm_beg=nfu(k0e-1)+1
1424      else
1425         knm_beg= 1
1426      endif
1427c
1428      knm_end=nfu(k0b  )
1429c
1430      kn0_beg=nfu(k0e)+1
1431      kn0_end=nfu(k0b+1)
1432c------------------------------------------------------
1433      do inp=ibeg,iend
1434         in0=ifrst(inp)
1435         icr=icoor(inp)
1436         inm=nmxyz(icr,in0)
1437         habcdin0=habcd(icr,in0)
1438c
1439c do kn0 & knp, do not do knm :
1440c
1441         if(inm.gt.0) then
1442c.ok.       do kn0=nfu(k0e)+1,nfu(k0b+1)
1443            do kn0=kn0_beg,kn0_end
1444               knp=npxyz(icr,kn0)
1445               do i=1,nbls
1446                  xt1(i,inp,kn0)=p1234(i,icr)*xt2(i,in0,kn0)
1447     *                                       -xt2(i,in0,knp)
1448     *                              +habcdin0*xt3(i,inm,kn0)
1449               enddo
1450            enddo
1451         else
1452c.ok.       do kn0=nfu(k0e)+1,nfu(k0b+1)
1453            do kn0=kn0_beg,kn0_end
1454               knp=npxyz(icr,kn0)
1455               do i=1,nbls
1456                  xt1(i,inp,kn0)=p1234(i,icr)*xt2(i,in0,kn0)
1457     *                                       -xt2(i,in0,knp)
1458               enddo
1459            enddo
1460         endif            !  if(inm.gt.0) then
1461c
1462c do knm only :
1463c
1464         do knm=knm_beg,knm_end
1465            kn0=npxyz(icr,knm)
1466            habcdkn0=habcd(icr,kn0)
1467            do i=1,nbls
1468               xt1(i,inp,kn0)=xt1(i,inp,kn0)
1469     *              +habcdkn0*xt2(i,in0,knm)
1470            enddo
1471         enddo
1472c
1473c because of different formulation for trackl_ than for tracij_
1474c recursive multipy everything by abcd(=(c+d)/(a+b) )
1475c
1476         do kn0=kn0_beg,kn0_end
1477            do i=1,nbls
1478               xt1(i,inp,kn0)=xt1(i,inp,kn0)*abcd
1479            enddo
1480         enddo
1481c
1482      enddo               !  do inp=ibeg,iend
1483c
1484c------------ saving  target classes ------------------
1485      if(nrec.ge.nqij) then
1486         do kn0=nfu(nqkl)+1,nfu(nskl+1)
1487            do inp=ibeg,iend
1488               do i=1,nbls
1489                  xt0(i,inp,kn0)=xt1(i,inp,kn0)
1490               enddo
1491            enddo
1492         enddo
1493      endif
1494c---------------- target classes ----------------------
1495c
1496      end
1497c=================================================================
1498c     subroutine rescale_wt0(nbls1,wt0,ninteg,const)
1499c     implicit real*8 (A-H,O-Z)
1500c     common /resc_ssssm/ resc(1000)
1501c     dimension const(nbls1)
1502c     dimension wt0(nbls1,ninteg)
1503c
1504c     do ijkl=1,nbls1
1505c        xscale=const(ijkl)*resc(ijkl)
1506c        do integ=1,ninteg
1507c           wt0(ijkl,integ)=wt0(ijkl,integ)*xscale
1508c        enddo
1509c     enddo
1510c
1511c     end
1512c=================================================================
1513      subroutine how_2_shift(stable,nsij,nskl,mmax,mmax1,immax1,
1514     $     kmmax1,lobsa)
1515      character*11 scftype
1516      character*8 where
1517      common /runtype/ scftype,where
1518      logical stable
1519      common/shell/lshellt,lshelij,lshelkl,lhelp,lcas2(4),lcas3(4)
1520c-----------------------------------------------------------------
1521c Since Tracy's recursive might suffer from numerical instability
1522c if an exponent on the second center is much bigger than others,
1523c one has to shift an ang.mom. from the THIRD to the FIRST position
1524c in such cases (the usual way is just opposite).
1525c This decision is made here using logical variable stable() which
1526c is determined in the XWPQ_1 (_2) subroutine.
1527c-----------------------------------------------------------------
1528c
1529cccc  if(nsij.ge.nskl) then
1530      if(stable) then
1531         immax=mmax-2
1532         kmmax=nskl-2
1533         lobsa=2
1534         if(lshelij.gt.0) lobsa=1
1535      else
1536         immax=nsij-2
1537         kmmax=mmax-2
1538         lobsa=4
1539         if(lshelkl.gt.0) lobsa=3
1540      endif
1541C-----------------------------------------------------------------------
1542         mmax1 =mmax
1543         immax1=immax
1544         kmmax1=kmmax
1545c
1546c giao derivatives:
1547         if(where.eq.'shif') then
1548            mmax1 = mmax-1
1549            immax1=immax-1
1550c           kmmax1=kmmax-1
1551            if(immax1.le.0) immax1=1
1552            if(kmmax1.le.0) kmmax1=1
1553         endif
1554c grad.derivatives:
1555         if(where.eq.'forc') then
1556            mmax1 = mmax-1
1557            immax1=immax-1    ! that's ok
1558            kmmax1=kmmax-1    ! try this ?
1559            if(immax1.le.0) immax1=1
1560            if(kmmax1.le.0) kmmax1=1
1561         endif
1562c hess.derivatives:
1563         if(where.eq.'hess') then
1564            mmax1 = mmax-2
1565c ?         immax1=immax-1    ! ?
1566c ?         kmmax1=kmmax-1
1567            if(immax1.le.0) immax1=1
1568            if(kmmax1.le.0) kmmax1=1
1569         endif
1570c---------------------------------------------------------------
1571      end
1572c=================================================================
1573