1C $Id$
2************************************************************************
3C> \ingroup nwint
4C> @{
5C>
6C> \brief Generate LLSS-type two-electron integrals
7C>
8C> This routine generates the LLSS-type two-electron integrals for
9C> a relativistic basis set, within the modified Dirac formalism.
10C> The integrals returned are \f$\nabla_C\cdot\nabla_D (ab|c^Sd^S)\f$,
11C> \f$\nabla_C\times\nabla_D (ab|c^Sd^S)\f$, or the 9 derivatives.
12C>
13C> Author: K. G. Dyall
14C>
15c:tex-\subsection{rel\_LLSS}
16c:tex-This routine generates the LLSS-type two-electron integrals for
17c:tex-a relativistic basis set, within the modified Dirac formalism.
18c:tex-The integrals returned are $\nabla_C\cdot\nabla_D (ab|c^Sd^S)$,
19c:tex-$\nabla_C\times\nabla_D (ab|c^Sd^S)$, 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_LLSS (
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 large 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 large 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_cd  ! n_cart_c*n_cart_d
87      integer n_cont_cd  ! NCC*NCD
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_c  ! cartesian components for l_C+1
96      integer n_cartp_d  ! cartesian components for l_D+1
97      integer n_cartm_c  ! cartesian components for l_C-1
98      integer n_cartm_d  ! cartesian components for l_D-1
99      integer n_intpp    ! number of integrals for l_C+1,l_D+1
100      integer n_intpm    ! number of integrals for l_C-1,l_D+1
101      integer n_intmp    ! number of integrals for l_C+1,l_D-1
102      integer n_intmm    ! number of integrals for l_C-1,l_D-1
103      integer i_xcc      ! address in scr of exp*coef for shell C
104      integer i_xcd      ! address in scr of exp*coef for shell D
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_LLSS'
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_LLSS',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_cd = n_cart_c*n_cart_d
136      n_cont_cd = NCC*NCD
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_LLSS',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_cd',n_cart_cd
152        write (LuOut,*) 'n_cont_cd',n_cont_cd
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_c = n_cart_c+l_C+2
165      n_cartp_d = n_cart_d+l_D+2
166      n_cartm_c = n_cart_c-l_C-1
167      n_cartm_d = n_cart_d-l_D-1
168      n_intpp = n_cartp_c*n_cartp_d*n_cont_cd*n_ab
169      n_intpm = n_cartm_c*n_cartp_d*n_cont_cd*n_ab
170      n_intmp = n_cartp_c*n_cartm_d*n_cont_cd*n_ab
171      n_intmm = n_cartm_c*n_cartm_d*n_cont_cd*n_ab
172      i_xcc = 1
173      i_xcd = i_xcc+NPC*NCC
174      i_pp = i_xcd+NPD*NCD
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_c',n_cartp_c
183        write (LuOut,*) 'n_cartp_d',n_cartp_d
184        write (LuOut,*) 'n_cartm_c',n_cartm_c
185        write (LuOut,*) 'n_cartm_d',n_cartm_d
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_xcc,i_xcd',i_xcc,i_xcd
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        write (LuOut,*) 'memscr,lscr',memscr,lscr
194      end if
195*
196*   Set up coefficients multiplied by exponents
197*
198      if (.not.DryRun) then
199        if (memscr .lt. 0) call errquit (
200     &      'Insufficient scratch memory in rel_LLSS',99, MEM_ERR)
201        k = i_xcc-1
202        do j = 1,NCC
203          do i = 1,NPC
204            scr(k+i) = zeta_C(i)*coef_C(i,j)
205          end do
206          k = k+NPC
207        end do
208        k = i_xcd-1
209        do j = 1,NCD
210          do i = 1,NPD
211            scr(k+i) = zeta_D(i)*coef_D(i,j)
212          end do
213          k = k+NPD
214        end do
215      end if
216*
217*         Calculate integrals for l_C+1, l_D+1
218*
219      if (debug_gen) write (LuOut,*) 'calling hf2 ++'
220      call hf2(
221     &    Axyz,zeta_A,coef_A,NPA,NCA,l_A,
222     &    Bxyz,zeta_B,coef_B,NPB,NCB,l_B,
223     &    Cxyz,zeta_C,scr(i_xcc),NPC,NCC,l_C+1,
224     &    Dxyz,zeta_D,scr(i_xcd),NPD,NCD,l_D+1,
225     &    scr(i_pp),n_intpp,canAB,canCD,canPQ,
226     &    DryRun,scr(i_scr),memscr)
227      if (DryRun) then
228        max_mem = max(max_mem,i_scr+memscr-1)
229      else if (debug_arrays) then
230        call ecp_matpr(scr(i_pp),1,n_intpp,1,1,
231     &      1,n_intpp,1,1,'++ ints','E',120,6)
232      end if
233*
234*         Calculate integrals for l_C-1, l_D+1
235*
236      if (l_C .gt. 0) then
237        memscr = lscr-i_scr+1
238        if (debug_gen) write (LuOut,*) 'calling hf2 +-'
239        call hf2(
240     &      Axyz,zeta_A,coef_A,NPA,NCA,l_A,
241     &      Bxyz,zeta_B,coef_B,NPB,NCB,l_B,
242     &      Cxyz,zeta_C,coef_C,NPC,NCC,l_C-1,
243     &      Dxyz,zeta_D,scr(i_xcd),NPD,NCD,l_D+1,
244     &      scr(i_pm),n_intpm,canAB,canCD,canPQ,
245     &      DryRun,scr(i_scr),memscr)
246        if (DryRun) then
247          max_mem = max(max_mem,i_scr+memscr-1)
248        else if (debug_arrays) then
249          call ecp_matpr(scr(i_pm),1,n_intpm,1,1,
250     &        1,n_intpm,1,1,'+- ints','E',120,6)
251        end if
252      end if
253*
254*         Calculate integrals for l_C+1, l_D-1
255*
256      if (l_D .gt. 0) then
257        memscr = lscr-i_scr+1
258        if (debug_gen) write (LuOut,*) 'calling hf2 -+'
259        call hf2(
260     &      Axyz,zeta_A,coef_A,NPA,NCA,l_A,
261     &      Bxyz,zeta_B,coef_B,NPB,NCB,l_B,
262     &      Cxyz,zeta_C,scr(i_xcc),NPC,NCC,l_C+1,
263     &      Dxyz,zeta_D,coef_D,NPD,NCD,l_D-1,
264     &      scr(i_mp),n_intmp,canAB,canCD,canPQ,
265     &      DryRun,scr(i_scr),memscr)
266        if (DryRun) then
267          max_mem = max(max_mem,i_scr+memscr-1)
268        else if (debug_arrays) then
269          call ecp_matpr(scr(i_mp),1,n_intmp,1,1,
270     &        1,n_intmp,1,1,'-+ ints','E',120,6)
271        end if
272*
273*         Calculate integrals for l_C-1, l_D-1
274*
275        if (l_C .gt. 0) then
276          memscr = lscr-i_scr+1
277          if (debug_gen) write (LuOut,*) 'calling hf2 --'
278          call hf2(
279     &        Axyz,zeta_A,coef_A,NPA,NCA,l_A,
280     &        Bxyz,zeta_B,coef_B,NPB,NCB,l_B,
281     &        Cxyz,zeta_C,coef_C,NPC,NCC,l_C-1,
282     &        Dxyz,zeta_D,coef_D,NPD,NCD,l_D-1,
283     &        scr(i_mm),n_intmm,canAB,canCD,canPQ,
284     &        DryRun,scr(i_scr),memscr)
285          if (DryRun) then
286            max_mem = max(max_mem,i_scr+memscr-1)
287          else if (debug_arrays) then
288            call ecp_matpr(scr(i_mm),1,n_intmm,1,1,
289     &          1,n_intmm,1,1,'-- ints','E',120,6)
290          end if
291        end if
292      end if
293*
294*     Compute the relativistic integrals
295*
296      memscr = lscr-i_scr+1
297      if (debug_gen) write (LuOut,*) 'calling rel_pot'
298      call rel_pot (
299     &    scr(i_pp),scr(i_pm),scr(i_mp),scr(i_mm),rel_ints,n_ints,ntyp,
300     &    l_C,n_cartp_c,n_cart_c,n_cartm_c,NCC*n_ab,
301     &    l_D,n_cartp_d,n_cart_d,n_cartm_d,NCD,
302     &    DryRun,scr(i_scr),memscr,ibug)
303      if (DryRun) then
304        max_mem = max(max_mem,i_scr+memscr-1)
305        lscr = max_mem
306      end if
307      if (debug_gen) write (LuOut,*) 'Exiting rel_LLSS'
308*
309      return
310      end
311C>
312C> @}
313