1C $Id$ 2************************************************************************ 3* * 4 subroutine ecp_radint2 (p_min,p_max, 5 & l_a,n_prim_a,n_cont_a,n_int_a,coeff_a,ai,alpha, 6 & l_b,n_prim_b,n_cont_b,n_int_b,coeff_b,bi,beta, 7 & l_c,n_prim_c,n_cont_c,n_c,coeff_c,gamma,prefactor,tol,sphcart, 8 & n_ab,n_abc,n_Q,m_count,temp,Qabc,Qab,Q,Q_int,ibug) 9* * 10* Calculate Type 2 radial integrals for a given ECP centre, angular * 11* projector and exponent p * 12* * 13* Argument (status) - description * 14* * 15* p_min - minimum power of r in ECP expansion * 16* p_max - maximum power of r in ECP expansion * 17* l_a (inp) - (maximum) angular momentum of functions on centre a * 18* n_prim_a (inp) - number of primitive functions on centre a * 19* n_cont_a (inp) - number of contracted functions on centre a * 20* n_int_a - total number of m values for centre a * 21* coeff_a - contraction coefficients of basis functions on centre a * 22* ai (inp) - values of 1/2*zeta_a*R_ac * 23* alpha (inp) - values of a/2*sqrt(c) * 24* l_b (inp) - (maximum) angular momentum of functions on centre b * 25* n_prim_b (inp) - number of primitive functions on centre b * 26* n_cont_b (inp) - number of contracted functions on centre b * 27* n_int_b - total number of m values for centre b * 28* coeff_b - contraction coefficients of basis functions on centre b * 29* bi (inp) - values of 1/2*zeta_b*R_bc * 30* beta (inp) - values of b/2*sqrt(c) * 31* l_c (inp) - angular momentum of ECP projector on centre c * 32* n_c (inp) - total number of ECP primitive functions * 33* n_prim_c (inp) - number of primitive functions for each power of r * 34* in ECP expansion * 35* coeff_c - contraction coefficients of potential on centre c * 36* gamma (inp) - values of 1/sqrt(c) * 37* prefactor (inp) - exp[(alpha+beta)^2-zeta_a*R_ac^2-zeta_b*R_bc^2] * 38* tol (inp) - maximum relative error in bessel functions * 39* sphcart - 1 for spherical basis, 0 for cartesian basis. * 40* n_ab (inp) - n_prim_a*n_prim_b * 41* n_abc (inp) - n_prim_a*n_prim_b*n_c * 42* n_Q - dimension of intermediate array Q * 43* m_count (inp) - number of Q^2_{mm} to be generated * 44* temp - work array * 45* Qabc - uncontracted Q integrals used in upward recursion in k * 46* Qab - Q integrals contracted over core potential with maximum m, n * 47* values * 48* Q - work array used in downward recursion * 49* Q_int - final fully contracted Q integrals * 50* ibug - debug flag. 0 for no debug, 1 for address printing, 2 for * 51* array printing, 3 for both. * 52* * 53* Written by K. G. Dyall * 54* * 55************************************************************************ 56 implicit none 57#include "stdio.fh" 58#include "ecp_consts.fh" 59#include "errquit.fh" 60 integer h,i,j,ibug,ind_a,ind_b,ind_p,ind_even,ind_odd,i_c, 61 & ind_max,ind_mid,ind_min,ind_hi,ind_mi,ind_lo, 62 & k,k_max,k_p, 63 & l_a,l_b,l_c,l_max,l_mid,l_min, 64 & l,lambda_a,lambda_b,la_min,lb_min, 65 & m,m_aeqb,m_count,m_max,m_min,m_min1, 66 & mc_a,mc_aeqb,mc_b,mc_even,mc_odd, 67 & mt_max,mt_mid,mt_min, 68 & n,n_c,n_cont_c,n_ab,n_abc,n_abp,n_Q,n_max,n_abm, 69 & n_prim_a,n_prim_b,n_cont_a,n_cont_b,n_int_a,n_int_b, 70 & p,p_max,p_min,sphcart 71 integer n_prim_c(p_min:p_max) 72 logical debug_gen,debug_addresses,debug_arrays 73 double precision ai(n_abc),bi(n_abc),coeff_c(n_c), 74 & alpha(n_abc),beta(n_abc),gamma(n_abc),prefactor(n_abc), 75 & coeff_a(n_prim_a,n_cont_a),coeff_b(n_prim_b,n_cont_b), 76 & Qabc(n_abc,m_count,3),Qab(n_ab,0:l_a,0:l_b,n_cont_c), 77 & Q(n_ab,n_Q,n_cont_c),temp(n_abc,20), 78 & Q_int(n_cont_a*n_cont_b,n_int_a,n_int_b,n_cont_c), 79 & tol 80* 81 debug_gen = ibug .gt. 0 82 debug_addresses = mod(ibug,2) .eq. 1 83 debug_arrays = mod(ibug,10)/2 .eq. 1 84* 85 if (debug_gen) then 86 write (LuOut,'(//A,/)') 'Entering ecp_radint2 ...' 87 write (LuOut,*) 'ibug =',ibug 88 end if 89 if (debug_addresses) then 90 write (LuOut,*) 'l_a,l_b,l_c',l_a,l_b,l_c 91 write (LuOut,*) 'n_ab,n_abc,n_Q,m_count', 92 & n_ab,n_abc,n_Q,m_count 93 end if 94 k_max = l_a+l_b 95 l_max = max(l_a,l_b) 96 l_min = min(l_a,l_b) 97 if (debug_addresses) write (LuOut,*) 'k_max,l_max,l_min', 98 & k_max,l_max,l_min 99 mt_max = l_max*(l_max+1)/2 100 mt_min = l_min*(l_min-1)/2 101 mt_mid = mt_max 102 if (l_a .eq. l_b) then 103 m_aeqb = 0 104 else 105 mt_max = mt_max-l_max 106 m_aeqb = 1 107 end if 108 if (debug_addresses) write (LuOut,*) 'mt_max,mt_mid,mt_min', 109 & mt_max,mt_mid,mt_min 110* 111 ind_max = 1 112 ind_mid = ind_max+mt_max 113 ind_min = ind_mid+mt_mid 114 if (debug_addresses) write (LuOut,*) 'ind_max,ind_mid,ind_min', 115 & ind_max,ind_mid,ind_min 116 if (l_a .ge. l_b) then 117 j = 1 118 h = 0 119 else 120 j = 0 121 h = 1 122 end if 123 if (debug_addresses) write (LuOut,*) 'j,h',j,h 124* 125* Zero out arrays for accumulation 126* 127 n_abm = n_ab*(mt_min+mt_mid+mt_max) 128 call dcopy(n_abm,zero,0,Q,1) 129 n_abm = (l_a+1)*(l_b+1)*n_ab*n_cont_c 130 call dcopy(n_abm,zero,0,Qab,1) 131* 132* Set up initial values for K = 0 and 1 133* 134 m_min = l_c 135 m_max = m_min+m_count-1 136 if (debug_addresses) write (LuOut,*) 'm_count,m_min,m_max', 137 & m_count,m_min,m_max 138* 139* Loop over p values to generate primitive integrals 140* 141 ind_p = 1 142 if (debug_addresses) write (LuOut,*) 'n_abc',n_abc 143 do p = p_min,p_max 144 n_abp = n_ab*n_prim_c(p) 145 if (debug_addresses) write (LuOut,*) 'p =',p 146 if (debug_addresses) write (LuOut,*) n_abp,n_ab,n_prim_c(p) 147 if (p .eq. 4) then 148 if (n_abp .gt. 0) call ecp_t2_init4 (n_abp,n_abc,m_min,m_max, 149 & j,h,tol,ai,bi,alpha(ind_p),beta(ind_p),gamma(ind_p), 150 & prefactor(ind_p),temp(1,2),temp(1,1),Qabc(ind_p,1,1)) 151 else if (p .eq. 3) then 152 if (n_abp .gt. 0) call ecp_t2_init3 (n_abp,n_abc,m_min,m_max, 153 & j,h,tol,ai,bi,alpha(ind_p),beta(ind_p),gamma(ind_p), 154 & prefactor(ind_p),temp(1,2),temp(1,1),Qabc(ind_p,1,1)) 155 else if (p .eq. 2) then 156 if (n_abp .gt. 0) call ecp_t2_init2 (n_abp,n_abc,m_min,m_max, 157 & h,tol,alpha(ind_p),beta(ind_p),gamma(ind_p), 158 & prefactor(ind_p),temp,temp(1,6),temp(1,7),Qabc(ind_p,1,1)) 159 else if (p .eq. 1) then 160 if (n_abp .gt. 0) call ecp_t2_init1 (n_abp,n_abc,m_min,m_max, 161 & j,h,tol,ai,bi,alpha(ind_p),beta(ind_p),gamma(ind_p), 162 & prefactor(ind_p),temp(1,2),temp(1,1),Qabc(ind_p,1,1)) 163 else if (p .eq. 0) then 164 if (n_abp .gt. 0) call ecp_t2_init0 (n_abp,n_abc,m_min,m_max, 165 & j,h,tol,ai,bi,alpha(ind_p),beta(ind_p),gamma(ind_p), 166 & prefactor(ind_p),temp(1,2),temp(1,1),Qabc(ind_p,1,1)) 167 else 168 call errquit( 169 & 'Illegal p value in routine ecp_radint2', 170 & p, BASIS_ERR) 171 end if 172 if (debug_addresses) write (LuOut,*) 'end p loop' 173 ind_p = ind_p+n_abp 174 end do 175 if (debug_arrays) then 176 call ecp_matpr (Qabc(1,1,1),1,n_abc,m_min,m_max, 177 & 1,n_abc,m_min,m_max,'Qabc(i,m,1)','E',78,4) 178 call ecp_matpr (Qabc(1,1,2),1,n_abc,m_min+1,m_max, 179 & 1,n_abc,m_min+1,m_max,'Qabc(i,m,2)','E',78,4) 180 end if 181* 182* Contract over ECP functions for lambda_a, lambda_b = 0,0 183* 184 do i_c = 1,n_cont_c 185 call ecp_contract (n_ab,n_c,1,Qabc,coeff_c,Qab(1,0,0,i_c)) 186 if (debug_arrays) call ecp_matpr (Qab(1,0,0,i_c),1,n_prim_a, 187 & 1,n_prim_b,1,n_prim_a,1,n_prim_b,'Qab(i,0,0)','E',78,4) 188 if (k_max .ge. 2) call ecp_contract (n_ab,n_c,l_max,Qabc(1,2,1), 189 & coeff_c,Q(1,ind_mid,i_c)) 190* 191* Contract over ECP functions for lambda_a, lambda_b = 0,1 or 1,0 192* 193 if (k_max .gt. 0) then 194 call ecp_contract (n_ab,n_c,1,Qabc(1,1,2),coeff_c, 195 & Qab(1,j,h,i_c)) 196 if (debug_arrays) then 197 if (j .gt. h) then 198 call ecp_matpr (Qab(1,j,h,i_c),1,n_prim_a,1,n_prim_b, 199 & 1,n_prim_a,1,n_prim_b,'Qab(i,1,0)','E',78,4) 200 else 201 call ecp_matpr (Qab(1,j,h,i_c),1,n_prim_a,1,n_prim_b, 202 & 1,n_prim_a,1,n_prim_b,'Qab(i,0,1)','E',78,4) 203 end if 204 end if 205 if (k_max .ge. 2) call ecp_contract (n_ab,n_c,l_max-m_aeqb, 206 & Qabc(1,2,2),coeff_c,Q(1,ind_max,i_c)) 207 if (debug_arrays) then 208 if (mt_max .gt. 0) write (LuOut,'(1p5e15.7)') 209 & (Q(1,ind_max+i,i_c),i=0,mt_max-1) 210 if (mt_mid .gt. 0) write (LuOut,'(1p5e15.7)') 211 & (Q(1,ind_mid+i,i_c),i=0,mt_mid-1) 212 if (mt_min .gt. 0) write (LuOut,'(1p5e15.7)') 213 & (Q(1,ind_min+i,i_c),i=0,mt_min-1) 214 end if 215 end if 216 end do 217* 218* Upward recursion in k to maximum. Recursion is performed for 219* m if l_a > l_b otherwise for n. 220* 221 if (debug_gen) write (LuOut,'(/A,/)') 'Starting upward recursion' 222 ind_lo = 1 223 ind_mi = 2 224 ind_hi = 3 225 mc_odd = l_max-m_aeqb 226 mc_even = l_max-1 227 ind_odd = ind_max+mc_odd 228 mc_odd = mc_odd-1 229 ind_even = ind_mid+l_max 230* 231 do k_p = 2,k_max 232 l_mid = k_p/2 233 if (debug_gen) write (LuOut,*) 'k_p,l_mid',k_p,l_mid 234 if (mod(k_p,2) .eq. 0) then 235* 236* Even k recursion to obtain Q^k_{mm} 237* 238 m_max = m_max-1 239 m_min = m_min+1 240 ind_p = 1 241 if (debug_addresses) write (LuOut,*) 'k,m_min,m_max', 242 & k_p,m_min,m_max 243 do p = p_min,p_max 244 n_abp = n_ab*n_prim_c(p) 245 k = k_p+p-2 246 if (debug_arrays) then 247 write(LuOut,'(1p5e15.7)') (Qabc(1,i+2,ind_lo),i=0, 248 & m_max-m_min) 249 write(LuOut,'(1p5e15.7)') (Qabc(1,i+1,ind_mi),i=0, 250 & m_max-m_min) 251 end if 252 call ecp_up_k (m_min,m_max,k,0,j,n_abp,n_abc, 253 & alpha(ind_p),beta(ind_p),gamma(ind_p), 254 & Qabc(ind_p,2,ind_lo),Qabc(ind_p,1,ind_mi), 255 & Qabc(ind_p,1,ind_hi)) 256 if (debug_arrays) write (LuOut,'(1p5e15.7)') 257 & (Qabc(1,i+1,ind_hi),i=0,m_max-m_min) 258 ind_p = ind_p+n_abp 259 end do 260* 261* Contract over ECP functions 262* 263 if (debug_addresses) then 264 write (LuOut,*) 'ind_hi',ind_hi 265 write (LuOut,*) 'ind_even,mc_even',ind_even,mc_even 266 end if 267 do i_c = 1,n_cont_c 268 if (mc_even .gt. 0) call ecp_contract (n_ab,n_c,mc_even, 269 & Qabc(1,2,ind_hi),coeff_c,Q(1,ind_even,i_c)) 270 if ((l_mid .le. l_a) .and. (l_mid .le. l_b)) then 271 call ecp_contract (n_ab,n_c,1,Qabc(1,1,ind_hi), 272 & coeff_c,Qab(1,l_mid,l_mid,i_c)) 273 if (debug_arrays) call ecp_matpr (Qab(1,l_mid,l_mid,i_c), 274 & 1,n_prim_a,1,n_prim_b,1,n_prim_a,1,n_prim_b, 275 & 'Qab(*,l_mid,l_mid)','E',78,4) 276 end if 277 end do 278 ind_even = ind_even+mc_even 279 mc_even = mc_even-1 280 else 281* 282* Odd k recursion to obtain Q^k_{m+1m} or Q^k_{mm+1} 283* 284 ind_p = 1 285 do p = p_min,p_max 286 n_abp = n_ab*n_prim_c(p) 287 k = k_p+p-2 288 call ecp_up_k (m_min,m_max,k,1,h,n_abp,n_abc, 289 & alpha(ind_p),beta(ind_p),gamma(ind_p), 290 & Qabc(ind_p,2,ind_lo),Qabc(ind_p,1,ind_mi), 291 & Qabc(ind_p,1,ind_hi)) 292 ind_p = ind_p+n_abp 293 end do 294* 295* Contract over ECP functions 296* 297 do i_c = 1,n_cont_c 298 if (mc_odd .gt. 0) call ecp_contract (n_ab,n_c,mc_odd, 299 & Qabc(1,2,ind_hi),coeff_c,Q(1,ind_odd,i_c)) 300 if ((l_mid+j .le. l_a) .and. (l_mid+h .le. l_b)) then 301 call ecp_contract (n_ab,n_c,1,Qabc(1,1,ind_hi), 302 & coeff_c,Qab(1,l_mid+j,l_mid+h,i_c)) 303 if (debug_arrays) call ecp_matpr (Qab(1,l_mid+j,l_mid+h, 304 & i_c),1,n_prim_a,1,n_prim_b,1,n_prim_a,1,n_prim_b, 305 & 'Qab(*,l_mid+j,l_mid+h)','E',78,4) 306 end if 307 end do 308 ind_odd = ind_odd+mc_odd 309 mc_odd = mc_odd-1 310 end if 311 i = ind_lo 312 ind_lo = ind_mi 313 ind_mi = ind_hi 314 ind_hi = i 315 end do 316* 317* Downward recursions in m,n 318* 319 m_min = l_c 320 ind_odd = ind_max 321 ind_even = ind_mid 322 if (l_a .ge. l_b) then 323 ind_a = ind_max 324 ind_b = ind_min 325 else 326 ind_a = ind_min 327 ind_b = ind_max 328 end if 329 mc_a = l_a-1 330 mc_aeqb = l_a-m_aeqb 331 mc_b = l_b-1 332 mc_even = l_max 333 if (debug_gen) write(LuOut,'(/A,/)') 'Starting downward recursion' 334 if (debug_addresses) then 335 write (LuOut,*) 'ind_max,ind_mid,ind_min', 336 & ind_max,ind_mid,ind_min 337 write (LuOut,*) 'mt_max,mt_mid,mt_min',mt_max,mt_mid,mt_min 338 end if 339C write (LuOut,'(1p5e15.7)') (Q(1,ind_max+i),i=0,mt_max-1) 340C write (LuOut,'(1p5e15.7)') (Q(1,ind_mid+i),i=0,mt_mid-1) 341C write (LuOut,'(1p5e15.7)') (Q(1,ind_min+i),i=0,mt_min-1) 342 do k_p = 1,k_max 343 l_mid = k_p/2 344 if (debug_gen) write (LuOut,*) 'k_p,l_mid',k_p,l_mid 345* 346* Recursion from Q^k_{mm} for even k 347* 348 if (mod(k_p,2) .eq. 0) then 349 m_min = m_min+1 350* 351* Shift angular momentum from centre b to centre a, by 352* downward recursion on m_b 353* 354 m_max = m_min+mc_a-1 355 n_max = min(l_mid,l_a-l_mid)-1 356 ind_hi = ind_even 357 ind_lo = ind_a 358 if (debug_addresses) 359 & write (LuOut,*) 'k,m_min,m_max',k_p,m_min,m_max 360C write (LuOut,'(1p5e15.7)') (Q(1,ind_max+i),i=0,mt_max-1) 361C write (LuOut,'(1p5e15.7)') (Q(1,ind_mid+i),i=0,mt_mid-1) 362C write (LuOut,'(1p5e15.7)') (Q(1,ind_min+i),i=0,mt_min-1) 363 do n = 0,n_max 364 lambda_a = l_mid+n+1 365 lambda_b = l_mid-n-1 366 do i_c = 1,n_cont_c 367 call ecp_down_m (m_min-n,m_max-2*n,n_ab,bi, 368 & Q(1,ind_hi,i_c),Q(1,ind_lo,i_c),Q(1,ind_lo,i_c)) 369 if ((lambda_a .le. l_a) .and. (lambda_b .le. l_b)) 370 & call dcopy (n_ab,Q(1,ind_lo,i_c),1, 371 & Qab(1,lambda_a,lambda_b,i_c),1) 372 end do 373 ind_hi = ind_lo+1 374 ind_lo = ind_hi-mc_aeqb-n 375 end do 376* 377* .. and shift angular momentum from centre a to centre b, by 378* downward recursion on m_a 379* 380 m_max = m_min+mc_b-1 381 n_max = min(l_mid,l_b-l_mid)-1 382 ind_hi = ind_even 383 ind_lo = ind_b 384 do n = 0,n_max 385 lambda_a = l_mid-n-1 386 lambda_b = l_mid+n+1 387 do i_c = 1,n_cont_c 388 call ecp_down_m (m_min-n,m_max-2*n,n_ab,ai, 389 & Q(1,ind_hi,i_c),Q(1,ind_lo,i_c),Q(1,ind_lo,i_c)) 390 if ((lambda_a .le. l_a) .and. (lambda_b .le. l_b)) 391 & call dcopy (n_ab,Q(1,ind_lo,i_c),1, 392 & Qab(1,lambda_a,lambda_b,i_c),1) 393 end do 394 ind_hi = ind_lo+1 395 ind_lo = ind_hi-mc_b-n 396 end do 397 ind_a = ind_a+mc_aeqb 398 ind_b = ind_b+mc_b 399 mc_a = mc_a-1 400 mc_aeqb = mc_aeqb-1 401 mc_b = mc_b-1 402 else 403* 404* Odd k : generate Q_mm+1 from Qm+1m or vice-versa, depending on 405* whether l_a > l_b or v.v. 406* 407 m_min1 = m_min+1 408 if (l_a .ge. l_b) then 409 m_max = m_min1+mc_b 410 if (debug_addresses) then 411 write (LuOut,*) 'ind_a+1,ind_even+1,ind_b', 412 & ind_a+1,ind_even+1,ind_b 413 write (LuOut,*) 'm_min1,m_max',m_min1,m_max 414 end if 415 do i_c = 1,n_cont_c 416 if (m_min1 .lt. m_max) call ecp_down_m (m_min1+1,m_max, 417 & n_ab,ai,Q(1,ind_a+1,i_c),Q(1,ind_even+1,i_c), 418 & Q(1,ind_b,i_c)) 419 if ((l_mid .le. l_a) .and. (l_mid+1 .le. l_b)) 420 & call ecp_down_m (m_min1,m_min1,n_ab,ai, 421 & Q(1,ind_a,i_c),Q(1,ind_even,i_c), 422 & Qab(1,l_mid,l_mid+1,i_c)) 423 end do 424 else 425 m_max = m_min1+mc_a 426 if (debug_addresses) then 427 write (LuOut,*) 'ind_b+1,ind_even+1,ind_a', 428 & ind_b+1,ind_even+1,ind_a 429 write (LuOut,*) 'm_min1,m_max',m_min1,m_max 430 end if 431 do i_c = 1,n_cont_c 432 if (m_min1 .lt. m_max) call ecp_down_m (m_min1+1,m_max, 433 & n_ab,bi,Q(1,ind_b+1,i_c),Q(1,ind_even+1,i_c), 434 & Q(1,ind_a,i_c)) 435 if ((l_mid+1 .le. l_a) .and. (l_mid .le. l_b)) 436 & call ecp_down_m (m_min1,m_min1,n_ab,bi, 437 & Q(1,ind_b,i_c),Q(1,ind_even,i_c), 438 & Qab(1,l_mid+1,l_mid,i_c)) 439 end do 440 end if 441 ind_even = ind_even+mc_even 442 mc_even = mc_even-1 443* 444* Shift angular momentum from centre b to centre a, by 445* downward recursion on m_b 446* 447 m_max = m_min+mc_a-1 448 n_max = min(l_mid,l_a-l_mid-1)-1 449 ind_hi = ind_a 450 ind_lo = ind_a-mc_aeqb 451 if (debug_addresses) 452 & write (LuOut,*) 'k,m_min,m_max',k_p,m_min,m_max 453 do n = 0,n_max 454 lambda_a = l_mid+n+2 455 lambda_b = l_mid-n-1 456 do i_c = 1,n_cont_c 457 call ecp_down_m (m_min-n,m_max-2*n,n_ab,bi, 458 & Q(1,ind_hi,i_c),Q(1,ind_lo,i_c),Q(1,ind_lo,i_c)) 459 if ((lambda_a .le. l_a) .and. (lambda_b .le. l_b)) 460 & call dcopy (n_ab,Q(1,ind_lo,i_c),1, 461 & Qab(1,lambda_a,lambda_b,i_c),1) 462 end do 463 ind_hi = ind_lo+1 464 ind_lo = ind_lo-mc_aeqb-n 465 end do 466* 467* .. and shift angular momentum from centre a to centre b, by 468* downward recursion on m_a 469* 470 m_max = m_min+mc_b-1 471 n_max = min(l_mid,l_b-l_mid-1)-1 472 ind_hi = ind_b 473 ind_lo = ind_b-mc_b 474 do n = 0,n_max 475 lambda_a = l_mid-n-1 476 lambda_b = l_mid+n+2 477 do i_c = 1,n_cont_c 478 call ecp_down_m (m_min-n,m_max-2*n,n_ab,ai, 479 & Q(1,ind_hi,i_c),Q(1,ind_lo,i_c),Q(1,ind_lo,i_c)) 480 if ((lambda_a .le. l_a) .and. (lambda_b .le. l_b)) 481 & call dcopy (n_ab,Q(1,ind_lo,i_c),1, 482 & Qab(1,lambda_a,lambda_b,i_c),1) 483 end do 484 ind_hi = ind_lo+1 485 ind_lo = ind_lo-mc_b-n 486 end do 487 end if 488 end do 489* 490* At this point, the highest m,n values in each block have been 491* generated. Now use downward recursion again to get all members of 492* the block. This has to be done in a different order because the 493* final result is indexed by m and n, not by (m-n). 494* 495 do i_c = 1,n_cont_c 496 ind_b = 0 497 do lambda_b = 0,l_b 498 if (sphcart .eq. 0) then 499 lb_min = max(0,lambda_b-(lambda_b+l_c)/2) 500 else 501 lb_min = max(0,lambda_b-l_c) 502 end if 503* 504* First do recursion in b to generate top row for all a. 505* 506 do lambda_a = 0,l_a 507 m = lambda_b+l_c+1 508 do l = lambda_b-1,lb_min,-1 509 m = m-2 510 call ecp_down_m (m,m,n_ab,bi, 511 & Qab(1,lambda_a,l+1,i_c),Qab(1,lambda_a,l,i_c), 512 & Qab(1,lambda_a,l,i_c)) 513 end do 514 end do 515* 516* Now contract over b and do recursion over a to generate final 517* integrals, contracting over a. 518* 519 do l = lambda_b,lb_min,-1 520 ind_b = ind_b+1 521 ind_a = 0 522 do lambda_a = 0,l_a 523 if (sphcart .eq. 0) then 524 la_min = max(0,lambda_a-(lambda_a+l_c)/2) 525 else 526 la_min = max(0,lambda_a-l_c) 527 end if 528 if (debug_addresses) write (LuOut,*) 'lambda_a,l', 529 & lambda_a,l 530 call dgemm ('N','N',n_prim_a,n_cont_b,n_prim_b,one, 531 & Qab(1,lambda_a,l,i_c),n_prim_a,coeff_b,n_prim_b,zero, 532 & Q(1,lambda_a+1,i_c),n_prim_a) 533 ind_a = ind_a+1 534 call dgemm ('T','N',n_cont_a,n_cont_b,n_prim_a,one, 535 & coeff_a,n_prim_a,Q(1,lambda_a+1,i_c),n_prim_a,zero, 536 & Q_int(1,ind_a,ind_b,i_c),n_cont_a) 537 if (debug_gen) 538 & write (LuOut,'(A,2I5)') 'ind_a,ind_b',ind_a,ind_b 539 if (debug_arrays) call ecp_matpr(Q_int(1,ind_a,ind_b,i_c), 540 & 1,n_cont_a,1,n_cont_b,1,n_cont_a,1,n_cont_b, 541 & 'Q_int','E',78,4) 542 m = lambda_a+l_c+1 543 do k = lambda_a,la_min+1,-1 544 m = m-2 545 call ecp_down_m (m,m,n_prim_a*n_cont_b, 546 & ai,Q(1,k+1,i_c),Q(1,k,i_c),Q(1,k,i_c)) 547 ind_a = ind_a+1 548 call dgemm ('T','N',n_cont_a,n_cont_b,n_prim_a,one, 549 & coeff_a,n_prim_a,Q(1,k,i_c),n_prim_a,zero, 550 & Q_int(1,ind_a,ind_b,i_c),n_cont_a) 551 if (debug_gen) 552 & write (LuOut,'(A,2I5)') 'ind_a,ind_b',ind_a,ind_b 553 if (debug_arrays) call ecp_matpr (Q_int(1,ind_a,ind_b, 554 & i_c),1,n_cont_a,1,n_cont_b,1,n_cont_a,1,n_cont_b, 555 & 'Q_int','E',78,4) 556 end do 557 end do 558 end do 559 end do 560 end do 561 if (debug_gen) write (LuOut,'(//A,/)') 'Exiting ecp_radint2' 562* 563 return 564 end 565