1C $Id$
2************************************************************************
3C> \ingroup nwint
4C> @{
5C>
6C> \brief Generate the SSLL-type two-electron integrals
7C>
8C> This routine generates the SSLL-type two-electron integrals for
9C> a relativistic basis set, within the modified Dirac formalism.
10C> The integrals returned are \f$\nabla_A\cdot\nabla_B (a^Sb^S|cd)\f$,
11C> \f$\nabla_A\times\nabla_B (a^Sb^S|cd)\f$, or the 9 derivatives.
12C>
13C> Author: K. G. Dyall
14C>
15c:tex-\subsection{rel\_SSLL}
16c:tex-This routine generates the SSLL-type two-electron integrals for
17c:tex-a relativistic basis set, within the modified Dirac formalism.
18c:tex-The integrals returned are $\nabla_A\cdot\nabla_B (a^Sb^S|cd)$,
19c:tex-$\nabla_A\times\nabla_B (a^Sb^S|cd)$, or the 9 derivatives.
20c:tex-
21c:tex-\noindent Author: K. G. Dyall
22c:tex-
23c:tex-{\it Syntax:}
24c:tex-\begin{verbatim}
25      subroutine rel_SSLL (
26     &    Axyz,zeta_A,coef_A,NPA,NCA,l_A,
27     &    Bxyz,zeta_B,coef_B,NPB,NCB,l_B,
28     &    Cxyz,zeta_C,coef_C,NPC,NCC,l_C,
29     &    Dxyz,zeta_D,coef_D,NPD,NCD,l_D,
30     &    rel_ints,n_ints,canAB,canCD,canPQ,
31     &    DryRun,scr,lscr,ibug,ntyp)
32c:tex-\end{verbatim}
33      implicit none
34#include "stdio.fh"
35#include "rel_consts.fh"
36#include "errquit.fh"
37*
38c:tex-{\it Argument list:}
39c:tex-\begin{verbatim}
40      integer NPA     !< [Input] num. prims in shell A
41      integer NCA     !< [Input] num general conts in shell A
42      integer l_A     !< [Input] angular momentum of shell A
43      integer NPB     !< [Input] num. prims in shell B
44      integer NCB     !< [Input] num general conts in shell B
45      integer l_B     !< [Input] angular momentum of shell B
46      integer NPC     !< [Input] num. prims in shell C
47      integer NCC     !< [Input] num general conts in shell C
48      integer l_C     !< [Input] angular momentum of shell C
49      integer NPD     !< [Input] num. prims in shell D
50      integer NCD     !< [Input] num general conts in shell D
51      integer l_D     !< [Input] angular momentum of shell D
52      integer lscr    !< [Input] size of scratch array
53      integer n_ints  !< [Input] size of any integral buffer
54      integer ibug    !< [Input] debug variable
55      integer ntyp    !< [Input] potential energy integral type
56      double precision Axyz(3)         !< [Input] position of center A
57      double precision zeta_A(NPA)     !< [Input] exponents of shell A
58      double precision coef_A(NPA,NCA) !< [Input] A small coeffs
59      double precision Bxyz(3)         !< [Input] position of center B
60      double precision zeta_B(NPB)     !< [Input] exponents of shell B
61      double precision coef_B(NPB,NCB) !< [Input] B small coeffs
62      double precision Cxyz(3)         !< [Input] position of center C
63      double precision zeta_C(NPC)     !< [Input] exponents of shell C
64      double precision coef_C(NPC,NCC) !< [Input] C small coeffs
65      double precision Dxyz(3)         !< [Input] position of center D
66      double precision zeta_D(NPD)     !< [Input] exponents of shell D
67      double precision coef_D(NPD,NCD) !< [Input] D small coeffs
68      double precision scr(lscr)       !< [Scratch] scratch buffers
69      double precision rel_ints(n_ints,ntyp)  !< [Output] LLSS integrals
70      logical canAB   !< [Input] compute only canonical ints (false only)
71      logical canCD   !< [Input] compute only canonical ints (false only)
72      logical canPQ   !< [Input] compute only canonical ints (false only)
73      logical DryRun  !< [Input] true means only compute required memory
74c:tex-\end{verbatim}
75c:tex-See rel_pot for a description of the allowed values of ibug and ntyp
76c:tex-Note that in the current version of this routine, the call to rel_pot
77c:tex-uses a dummy ntyp=1. It is kept in the input so that in future, the
78c:tex-spin-orbit integrals can also be obtained with a call to this routine.
79c:tex-
80c:tex-{\it Subroutines called:} hf2, rel\_pot, daxpy
81*
82      integer n_cart_a   ! cartesian components of shell A
83      integer n_cart_b   ! cartesian components of shell B
84      integer n_cart_c   ! cartesian components of shell C
85      integer n_cart_d   ! cartesian components of shell D
86      integer n_cart_ab  ! n_cart_a*n_cart_b
87      integer n_cont_ab  ! NCA*NCB
88      integer n_all_a    ! n_cart_a*NCA
89      integer n_all_b    ! n_cart_b*NCB
90      integer n_all_c    ! n_cart_c*NCC
91      integer n_all_d    ! n_cart_d*NCD
92      integer n_ab       ! number of ab densities
93      integer n_cd       ! number of cd densities
94      integer n_abcd     ! number of integrals
95      integer n_cartp_a  ! cartesian components for l_A+1
96      integer n_cartp_b  ! cartesian components for l_B+1
97      integer n_cartm_a  ! cartesian components for l_A-1
98      integer n_cartm_b  ! cartesian components for l_B-1
99      integer n_intpp    ! number of integrals for l_A+1,l_B+1
100      integer n_intpm    ! number of integrals for l_A-1,l_B+1
101      integer n_intmp    ! number of integrals for l_A+1,l_B-1
102      integer n_intmm    ! number of integrals for l_A-1,l_B-1
103      integer i_xca      ! address in scr of exp*coef for shell A
104      integer i_xcb      ! address in scr of exp*coef for shell B
105      integer i_pp       ! address in scr of integrals for l_A+1,l_B+1
106      integer i_pm       ! address in scr of integrals for l_A-1,l_B+1
107      integer i_mp       ! address in scr of integrals for l_A+1,l_B-1
108      integer i_mm       ! address in scr of integrals for l_A-1,l_B-1
109      integer i_scr      ! address of free space in scr
110      integer memscr     ! free space in scr
111      integer max_mem    ! maximum memory used
112      integer i,j,k      ! loop indices etc.
113      double precision one ! Obvious!
114      parameter (one = 1.0D0)
115*
116      logical debug_gen       ! do general debug printing
117      logical debug_addresses ! do address debug printing
118      logical debug_arrays    ! do array debug printing
119*
120      debug_gen = ibug .gt. 0
121      debug_addresses = mod(ibug,2) .eq. 1
122      debug_arrays = mod(ibug,10)/2 .eq. 1
123      max_mem = 0
124      if (debug_gen) write (LuOut,*) 'Entering rel_SSLL'
125*
126      if ((ntyp .ne. 1) .and. (ntyp .ne. 3) .and. (ntyp .ne. 4) .and.
127     &    (ntyp .ne. 9))
128     &    call errquit('Invalid value of ntyp in rel_SSLL',99,
129     &       UNKNOWN_ERR)
130*
131      n_cart_a = (l_a+1)*(l_a+2)/2
132      n_cart_b = (l_b+1)*(l_b+2)/2
133      n_cart_c = (l_c+1)*(l_c+2)/2
134      n_cart_d = (l_d+1)*(l_d+2)/2
135      n_cart_ab = n_cart_a*n_cart_b
136      n_cont_ab = NCA*NCB
137      n_all_a = n_cart_a*NCA
138      n_all_b = n_cart_b*NCB
139      n_all_c = n_cart_c*NCC
140      n_all_d = n_cart_d*NCD
141      n_ab = n_all_a*n_all_b
142      n_cd = n_all_c*n_all_d
143      n_abcd = n_ab*n_cd
144      if ((n_ints .lt. n_abcd) .and. (.not.DryRun)) call errquit (
145     &    'Integral buffer n_ints too small in rel_SSLL',99, MEM_ERR)
146      if (debug_addresses) then
147        write (LuOut,*) 'n_cart_a',n_cart_a
148        write (LuOut,*) 'n_cart_b',n_cart_b
149        write (LuOut,*) 'n_cart_c',n_cart_c
150        write (LuOut,*) 'n_cart_d',n_cart_d
151        write (LuOut,*) 'n_cart_ab',n_cart_ab
152        write (LuOut,*) 'n_cont_ab',n_cont_ab
153        write (LuOut,*) 'n_all_a',n_all_a
154        write (LuOut,*) 'n_all_b',n_all_b
155        write (LuOut,*) 'n_all_c',n_all_c
156        write (LuOut,*) 'n_all_d',n_all_d
157        write (LuOut,*) 'n_ab',n_ab
158        write (LuOut,*) 'n_cd',n_cd
159      end if
160*
161*   Set up pointers to scratch space for coefficients multiplied by
162*   exponents and for integrals with shifted l values
163*
164      n_cartp_a = n_cart_a+l_A+2
165      n_cartp_b = n_cart_b+l_B+2
166      n_cartm_a = n_cart_a-l_A-1
167      n_cartm_b = n_cart_b-l_B-1
168      n_intpp = n_cartp_a*n_cartp_b*n_cont_ab*n_cd
169      n_intpm = n_cartm_a*n_cartp_b*n_cont_ab*n_cd
170      n_intmp = n_cartp_a*n_cartm_b*n_cont_ab*n_cd
171      n_intmm = n_cartm_a*n_cartm_b*n_cont_ab*n_cd
172      i_xca = 1
173      i_xcb = i_xca+NPA*NCA
174      i_pp = i_xcb+NPB*NCB
175      i_pm = i_pp+n_intpp
176      i_mp = i_pm+n_intpm
177      i_mm = i_mp+n_intmp
178      i_scr = i_mm+n_intmm
179      memscr = lscr-i_scr+1
180*
181      if (debug_addresses) then
182        write (LuOut,*) 'n_cartp_a',n_cartp_a
183        write (LuOut,*) 'n_cartp_b',n_cartp_b
184        write (LuOut,*) 'n_cartm_a',n_cartm_a
185        write (LuOut,*) 'n_cartm_b',n_cartm_b
186        write (LuOut,*) 'n_intpp',n_intpp
187        write (LuOut,*) 'n_intpm',n_intpm
188        write (LuOut,*) 'n_intmp',n_intmp
189        write (LuOut,*) 'n_intmm',n_intmm
190        write (LuOut,*) 'i_xca,i_xcb',i_xca,i_xcb
191        write (LuOut,*) 'i_pp,i_pm,i_mp,i_mm',i_pp,i_pm,i_mp,i_mm
192        write (LuOut,*) 'i_scr',i_scr
193      end if
194*
195*   Set up coefficients multiplied by exponents
196*
197      if (.not.DryRun) then
198        if (memscr .lt. 0) call errquit (
199     &      'Insufficient scratch memory in rel_SSLL',99, MEM_ERR)
200        k = i_xca-1
201        do j = 1,NCA
202          do i = 1,NPA
203            scr(k+i) = zeta_A(i)*coef_A(i,j)
204          end do
205          k = k+NPA
206        end do
207        k = i_xcb-1
208        do j = 1,NCB
209          do i = 1,NPB
210            scr(k+i) = zeta_B(i)*coef_B(i,j)
211          end do
212          k = k+NPB
213        end do
214      end if
215*
216*         Calculate integrals for l_A+1, l_B+1
217*
218      call hf2(
219     &    Axyz,zeta_A,scr(i_xca),NPA,NCA,l_A+1,
220     &    Bxyz,zeta_B,scr(i_xcb),NPB,NCB,l_B+1,
221     &    Cxyz,zeta_C,coef_C,NPC,NCC,l_C,
222     &    Dxyz,zeta_D,coef_D,NPD,NCD,l_D,
223     &    scr(i_pp),n_intpp,canAB,canCD,canPQ,
224     &    DryRun,scr(i_scr),memscr)
225      if (DryRun) then
226        max_mem = max(max_mem,i_scr+memscr-1)
227      else if (debug_arrays) then
228        call ecp_matpr(scr(i_pp),1,n_intpp,1,1,
229     &      1,n_intpp,1,1,'++ ints','E',120,6)
230      end if
231*
232*         Calculate integrals for l_A-1, l_B+1
233*
234      if (l_A .gt. 0) then
235        memscr = lscr-i_scr+1
236        call hf2(
237     &      Axyz,zeta_A,coef_A,NPA,NCA,l_A-1,
238     &      Bxyz,zeta_B,scr(i_xcb),NPB,NCB,l_B+1,
239     &      Cxyz,zeta_C,coef_C,NPC,NCC,l_C,
240     &      Dxyz,zeta_D,coef_D,NPD,NCD,l_D,
241     &      scr(i_pm),n_intpm,canAB,canCD,canPQ,
242     &      DryRun,scr(i_scr),memscr)
243        if (DryRun) then
244          max_mem = max(max_mem,i_scr+memscr-1)
245        else if (debug_arrays) then
246          call ecp_matpr(scr(i_pm),1,n_intpm,1,1,
247     &        1,n_intpm,1,1,'+- ints','E',120,6)
248        end if
249      end if
250*
251*         Calculate integrals for l_A+1, l_B-1
252*
253      if (l_B .gt. 0) then
254        memscr = lscr-i_scr+1
255        call hf2(
256     &      Axyz,zeta_A,scr(i_xca),NPA,NCA,l_A+1,
257     &      Bxyz,zeta_B,coef_B,NPB,NCB,l_B-1,
258     &      Cxyz,zeta_C,coef_C,NPC,NCC,l_C,
259     &      Dxyz,zeta_D,coef_D,NPD,NCD,l_D,
260     &      scr(i_mp),n_intmp,canAB,canCD,canPQ,
261     &      DryRun,scr(i_scr),memscr)
262        if (DryRun) then
263          max_mem = max(max_mem,i_scr+memscr-1)
264        else if (debug_arrays) then
265          call ecp_matpr(scr(i_mp),1,n_intmp,1,1,
266     &        1,n_intmp,1,1,'-+ ints','E',120,6)
267        end if
268*
269*         Calculate integrals for l_A-1, l_B-1
270*
271        if (l_A .gt. 0) then
272          memscr = lscr-i_scr+1
273          call hf2(
274     &        Axyz,zeta_A,coef_A,NPA,NCA,l_A-1,
275     &        Bxyz,zeta_B,coef_B,NPB,NCB,l_B-1,
276     &        Cxyz,zeta_C,coef_C,NPC,NCC,l_C,
277     &        Dxyz,zeta_D,coef_D,NPD,NCD,l_D,
278     &        scr(i_mm),n_intmm,canAB,canCD,canPQ,
279     &        DryRun,scr(i_scr),memscr)
280          if (DryRun) then
281            max_mem = max(max_mem,i_scr+memscr-1)
282          else if (debug_arrays) then
283            call ecp_matpr(scr(i_mm),1,n_intmm,1,1,
284     &          1,n_intmm,1,1,'-- ints','E',120,6)
285          end if
286        end if
287      end if
288*
289*     Compute the relativistic integrals
290*
291      memscr = lscr-i_scr+1
292      call rel_pot2 (scr(i_pp),scr(i_pm),scr(i_mp),scr(i_mm),
293     &    rel_ints,n_ints,ntyp,n_cd,
294     &    l_A,n_cartp_a,n_cart_a,n_cartm_a,NCA,
295     &    l_B,n_cartp_b,n_cart_b,n_cartm_b,NCB,
296     &    DryRun,scr(i_scr),memscr,ibug)
297      if (DryRun) then
298        max_mem = max(max_mem,i_scr+memscr-1)
299        lscr = max_mem
300      end if
301*
302      return
303      end
304C>
305C> @}
306