1C $Id$
2************************************************************************
3C> \ingroup nwint
4C> @{
5C>
6C> \brief Generate the integral for various relativistic Hamiltonians
7C>
8C> This routine genrates the \f$pV\cdot p\f$ integrals needed for
9C> various relativistic Hamiltonians.
10C> \f{eqnarray*}{
11C>   \langle a | pV\cdot p| b \rangle_{ab} &=& \nabla_A\cdot\nabla_B V_{ab}^{SS}
12C> \f}
13C>
14C> Author: K. G. Dyall
15C>
16c:tex-\subsection{rel\_pvp}
17c:tex-This routine generates the pV.p integrals needed for various
18c:tex-relativistic Hamiltonians,
19c:tex-\begin{equation}
20c:tex- {\langle a | pV.p | b \rangle}_{ab} =
21c:tex-   \nabla_A\cdot\nabla_B V_{ab}^{SS}
22c:tex- \nonumber
23c:tex-\end{equation}
24c:tex-
25c:tex-\noindent Author: K. G. Dyall
26c:tex-
27c:tex-{\it Syntax:}
28c:tex-\begin{verbatim}
29      subroutine rel_pvp (
30     &    Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A,ictr_A,
31     &    Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B,ictr_B,
32     &    Cxyz,zan,exinv,nat,pVp,lpvp,canAB,DryRun,scr,lscr,
33     &    msg,ibug,ntyp)
34c:tex-\end{verbatim}
35      implicit none
36#include "stdio.fh"
37#include "rel_consts.fh"
38#include "errquit.fh"
39*
40c:tex-{\it Argument list:}
41c:tex-\begin{verbatim}
42      integer n_prim_A !< [Input] num. prims in shell A
43      integer n_cont_A !< [Input] num general conts in shell A
44      integer l_A      !< [Input] angular momentum of shell A
45      integer ictr_A   !< [Input] lexical atom index for shell A
46      integer n_prim_B !< [Input] num. prims in shell B
47      integer n_cont_B !< [Input] num general conts in shell B
48      integer l_B      !< [Input] angular momentum of shell B
49      integer ictr_B   !< [Input] lexical atom index for shell B
50      integer nat      !< [Input] number of atoms
51      integer lscr     !< [Input] size of scratch array
52      integer lpvp     !< [Input] size of integral buffer
53      integer ibug     !< [Input] debug variable
54      integer ntyp     !< [Input] potential energy integral type
55      double precision Axyz(3)          !< [Input] position of center A
56      double precision zeta_A(n_prim_A) !< [Input] exponents of shell A
57      double precision coefS_A(n_prim_A,n_cont_A) !< [Input] A small coeffs
58      double precision Bxyz(3)          !< [Input] position of center B
59      double precision zeta_B(n_prim_B) !< [Input] exponents of shell B
60      double precision coefS_B(n_prim_B,n_cont_B) !< [Input] B small coeffs
61      double precision Cxyz(3,nat)  !< [Input] all atom positions
62      double precision zan(nat)     !< [Input] charges on all atoms
63      double precision exinv(nat)   !< [Input] inverse nuclear exponents
64      double precision scr(lscr)    !< [Scratch] scratch buffers
65      double precision pVp(lpvp,ntyp) !< [Output] pV.p integrals
66      logical canAB   !< [Input] compute only canonical ints (false only)
67      logical DryRun  !< [Input] true means only compute required memory
68      character*(*) msg !< [Input] calling func. identification message
69c:tex-\end{verbatim}
70c:tex-See rel_pot for a description of the allowed values of ibug and ntyp
71c:tex-Note that in the current version of this routine, the call to rel_pot
72c:tex-uses a dummy ntyp=1. It is kept in the input so that in future, the
73c:tex-spin-orbit integrals can also be obtained with a call to this routine.
74c:tex-
75c:tex-{\it Subroutines called:} int\_hf1sp, rel\_pot, daxpy
76*
77      integer n_cart_a  ! cartesian components of shell A
78      integer n_cart_b  ! cartesian components of shell B
79      integer n_cart_ab ! n_cart_a*n_cart_b
80      integer n_cont_ab ! n_cont_a*n_cont_b
81      integer n_all_b   ! n_cart_b*n_cont_b
82      integer n_all_a   ! n_cart_a*n_cont_a
83      integer n_ab      ! number of integrals
84      integer n_cartp_a ! cartesian components for l_A+1
85      integer n_cartp_b ! cartesian components for l_B+1
86      integer n_cartm_a ! cartesian components for l_A-1
87      integer n_cartm_b ! cartesian components for l_B-1
88      integer n_intpp   ! number of integrals for l_A+1,l_B+1
89      integer n_intpm   ! number of integrals for l_A-1,l_B+1
90      integer n_intmp   ! number of integrals for l_A+1,l_B-1
91      integer n_intmm   ! number of integrals for l_A-1,l_B-1
92      integer i_xca     ! address in scr of exp*coef for shell A
93      integer i_xcb     ! address in scr of exp*coef for shell B
94      integer i_pp      ! address in scr of integrals for l_A+1,l_B+1
95      integer i_pm      ! address in scr of integrals for l_A-1,l_B+1
96      integer i_mp      ! address in scr of integrals for l_A+1,l_B-1
97      integer i_mm      ! address in scr of integrals for l_A-1,l_B-1
98      integer i_scr     ! address of free space in scr
99      integer memscr    ! free space in scr
100      integer max_mem   ! maximum memory used
101      integer i,j,k     ! loop indices etc.
102      double precision one ! Obvious!
103      parameter (one = 1.0D0)
104*
105      logical debug_gen       ! do general debug printing
106      logical debug_addresses ! do address debug printing
107      logical debug_arrays    ! do array debug printing
108      logical doS     ! compute overlap (True/False)
109      logical doT     ! compute kinetic (True/False)
110      logical doV     ! compute potential (True/False)
111*
112      debug_gen = ibug .gt. 0
113      debug_addresses = mod(ibug,2) .eq. 1
114      debug_arrays = (mod(ibug,10)/2 .eq. 1) .and. .not.DryRun
115      max_mem = 0
116*
117      if (debug_gen) then
118        write (LuOut,*) 'l_A',l_A
119        write (LuOut,*) 'n_prim_A',n_prim_A
120        write (LuOut,*) 'n_cont_A',n_cont_A
121        write (LuOut,*) 'ictr_A',ictr_A
122        write (LuOut,*) 'l_B',l_B
123        write (LuOut,*) 'n_prim_B',n_prim_B
124        write (LuOut,*) 'n_cont_B',n_cont_B
125        write (LuOut,*) 'ictr_B',ictr_B
126        write (LuOut,*) 'msg ',msg
127      end if
128*
129      n_cart_a = (l_a+1)*(l_a+2)/2
130      n_cart_b = (l_b+1)*(l_b+2)/2
131      n_cart_ab = n_cart_a*n_cart_b
132      n_cont_ab = n_cont_a*n_cont_b
133      n_all_a = n_cart_a*n_cont_a
134      n_all_b = n_cart_b*n_cont_b
135      n_ab = n_cart_ab*n_cont_ab
136      if (lpvp .lt. n_ab .and. .not.DryRun) call errquit (
137     &      'Integral buffer length too small in rel_pvp',99, INT_ERR)
138      if (debug_addresses) then
139        write (LuOut,*) 'n_cart_a',n_cart_a
140        write (LuOut,*) 'n_cart_b',n_cart_b
141        write (LuOut,*) 'n_cart_ab',n_cart_ab
142        write (LuOut,*) 'n_cont_ab',n_cont_ab
143        write (LuOut,*) 'n_all_a',n_all_a
144        write (LuOut,*) 'n_all_b',n_all_b
145        write (LuOut,*) 'n_ab',n_ab
146      end if
147      if (debug_arrays) then
148        call ecp_matpr (coefS_A,1,n_prim_a,1,n_cont_a,
149     &      1,n_prim_a,1,n_cont_a,'S coef A','E',120,6)
150        call ecp_matpr (coefS_B,1,n_prim_b,1,n_cont_b,
151     &      1,n_prim_b,1,n_cont_b,'S coef B','E',120,6)
152      end if
153*
154*   Generate small component potential arrays
155*
156*
157*   Set up pointers to scratch space for coefficients multiplied by
158*   exponents and for integrals with shifted l values
159*
160      n_cartp_a = n_cart_a+l_A+2
161      n_cartp_b = n_cart_b+l_B+2
162      n_cartm_a = n_cart_a-l_A-1
163      n_cartm_b = n_cart_b-l_B-1
164      n_intpp = n_cartp_a*n_cartp_b*n_cont_ab
165      n_intpm = n_cartm_a*n_cartp_b*n_cont_ab
166      n_intmp = n_cartp_a*n_cartm_b*n_cont_ab
167      n_intmm = n_cartm_a*n_cartm_b*n_cont_ab
168      i_xca = 1
169      i_xcb = i_xca+n_prim_A*n_cont_A
170      i_pp = i_xcb+n_prim_B*n_cont_B
171      i_pm = i_pp+n_intpp
172      i_mp = i_pm+n_intpm
173      i_mm = i_mp+n_intmp
174      i_scr = i_mm+n_intmm
175*
176      if (debug_addresses) then
177        write (LuOut,*) 'n_cartp_a',n_cartp_a
178        write (LuOut,*) 'n_cartp_b',n_cartp_b
179        write (LuOut,*) 'n_cartm_a',n_cartm_a
180        write (LuOut,*) 'n_cartm_b',n_cartm_b
181        write (LuOut,*) 'n_intpp',n_intpp
182        write (LuOut,*) 'n_intpm',n_intpm
183        write (LuOut,*) 'n_intmp',n_intmp
184        write (LuOut,*) 'n_intmm',n_intmm
185        write (LuOut,*) 'i_xca,i_xcb',i_xca,i_xcb
186        write (LuOut,*) 'i_pp,i_pm,i_mp,i_mm',i_pp,i_pm,i_mp,i_mm
187        write (LuOut,*) 'i_scr',i_scr
188      end if
189*
190*   Set up coefficients multiplied by exponents
191*
192      memscr = lscr-i_scr+1
193      if (.not.DryRun) then
194        if (memscr .lt. 0) call errquit (
195     &      'Insufficient scratch memory in rel_pvp',99, MEM_ERR)
196        k = i_xca-1
197        do j = 1,n_cont_a
198          do i = 1,n_prim_A
199            scr(k+i) = zeta_A(i)*coefS_A(i,j)
200          end do
201          k = k+n_prim_A
202        end do
203        k = i_xcb-1
204        do j = 1,n_cont_B
205          do i = 1,n_prim_B
206            scr(k+i) = zeta_B(i)*coefS_B(i,j)
207          end do
208          k = k+n_prim_A
209        end do
210      end if
211      doS = .false.
212      doT = .false.
213      doV = .true.
214*
215*         Calculate integrals for l_A+1, l_B+1
216*
217      call int_hf1sp(
218     &    Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1,ictr_A,
219     &    Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1,ictr_B,
220     &    Cxyz,zan,exinv,nat,scr,scr,scr(i_pp),n_intpp,doS,doT,doV,
221     &    canAB,DryRun,scr(i_scr),memscr,msg)
222      if (DryRun) then
223        max_mem = max(max_mem,i_scr+memscr-1)
224        memscr = lscr-i_scr+1
225      end if
226*
227*         Calculate integrals for l_A-1, l_B+1
228*
229      if (l_A .gt. 0) then
230        call int_hf1sp(
231     &      Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1,ictr_A,
232     &      Bxyz,zeta_B,scr(i_xcb),n_prim_B,n_cont_B,l_B+1,ictr_B,
233     &      Cxyz,zan,exinv,nat,scr,scr,scr(i_pm),n_intpm,doS,doT,doV,
234     &      canAB,DryRun,scr(i_scr),memscr,msg)
235        if (DryRun) then
236          max_mem = max(max_mem,i_scr+memscr-1)
237          memscr = lscr-i_scr+1
238        end if
239      end if
240*
241*         Calculate integrals for l_A+1, l_B-1
242*
243      if (l_B .gt. 0) then
244        call int_hf1sp(
245     &      Axyz,zeta_A,scr(i_xca),n_prim_A,n_cont_A,l_A+1,ictr_A,
246     &      Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1,ictr_B,
247     &      Cxyz,zan,exinv,nat,scr,scr,scr(i_mp),n_intmp,doS,doT,doV,
248     &      canAB,DryRun,scr(i_scr),memscr,msg)
249        if (DryRun) then
250          max_mem = max(max_mem,i_scr+memscr-1)
251          memscr = lscr-i_scr+1
252        end if
253*
254*         Calculate integrals for l_A-1, l_B-1
255*
256        if (l_A .gt. 0) then
257          call int_hf1sp(
258     &        Axyz,zeta_A,coefS_A,n_prim_A,n_cont_A,l_A-1,ictr_A,
259     &        Bxyz,zeta_B,coefS_B,n_prim_B,n_cont_B,l_B-1,ictr_B,
260     &        Cxyz,zan,exinv,nat,scr,scr,scr(i_mm),n_intmm,doS,doT,doV,
261     &        canAB,DryRun,scr(i_scr),memscr,msg)
262          if (DryRun) then
263            max_mem = max(max_mem,i_scr+memscr-1)
264            memscr = lscr-i_scr+1
265          end if
266        end if
267      end if
268*
269*     Compute the relativistic potential energy integrals
270*
271      call rel_pot (scr(i_pp),scr(i_pm),scr(i_mp),scr(i_mm),
272     &    pVp,lpvp,ntyp,
273     &    l_A,n_cartp_a,n_cart_a,n_cartm_a,n_cont_A,
274     &    l_B,n_cartp_b,n_cart_b,n_cartm_b,n_cont_B,
275     &    DryRun,scr(i_scr),memscr,ibug/10)
276      if (DryRun) then
277        max_mem = max(max_mem,i_scr+memscr-1)
278        lscr = max_mem
279      else if (debug_arrays) then
280        call ecp_matpr (pVp,1,n_all_b,1,n_all_a,
281     &      1,n_all_b,1,n_all_a,'pV.p integrals','E',120,6)
282      end if
283*
284      return
285      end
286C>
287C> @}
288