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