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