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