1C $Id$ 2************************************************************************ 3* * 4 subroutine ecp_int1 (mem_max,DryRun, 5 & R_AC,X_AC,Y_AC,Z_AC,l_a,n_prim_a,n_cont_a,coef_a,zeta_a, 6 & R_BC,X_BC,Y_BC,Z_BC,l_b,n_prim_b,n_cont_b,coef_b,zeta_b, 7 & l_ecp,n_prim_c,n_coef_c, 8 & ind_z,ind_c,n_z_tot,n_c_tot,l_ecp_max,n_cent, 9 & zeta_c,coef_c,p_min,p_max,tol,sphcart,tmp,l_tmp, 10 & csco,lcsco, 11 & ecp_ints,n_int,n_blk,transpose,ibug) 12* * 13* Calculate Type 1 radial integrals for a given ECP centre * 14* * 15* Argument (status) - description * 16* * 17* mem_max (out) - maximum scratch memory required * 18* DryRun (inp) - logical to only return memory if true * 19* R_AC (inp) - distance between centres A and C * 20* X_AC,Y_AC,Z_AC (inp) - cartesian coordinates of centre C relative * 21* to centre A, X_AC = X_C - X_A, etc. * 22* l_a (inp) - (maximum) angular momentum of functions on centre A * 23* n_prim_a (inp) - number of primitive functions on centre A * 24* n_cont_a (inp) - number of contracted functions on centre A * 25* coef_a (inp) - centre A contraction coefficients * 26* zeta_a (inp) - centre A exponents * 27* R_BC (inp) - distance between centres B and C * 28* X_BC,Y_BC,Z_BC (inp) - cartesian coordinates of centre C relative * 29* to centre B, X_BC = X_C - X_B, etc. * 30* l_b (inp) - (maximum) angular momentum of functions on centre B * 31* n_prim_b (inp) - number of primitive functions on centre B * 32* n_cont_b (inp) - number of contracted functions on centre B * 33* coef_b (inp) - centre B contraction coefficients * 34* zeta_b (inp) - centre B exponents * 35* l_ecp_max (inp) - maximum angular momentum of ECP projector * 36* n_prim_c (inp) - number of primitive functions for each power of r * 37* in ECP expansion * 38* n_coef_c (inp) - array of number of coefficients/exponents for * 39* each projector * 40* ind_z (inp) - array of addresses of first exponent for each l * 41* ind_c (inp) - array of addresses of first coefficient for each l * 42* n_z_tot (inp) - total number of ECP exponents * 43* n_c_tot (inp) - total number of ECP coefficients * 44* n_cent (inp) - number of ECP centres * 45* zeta_c (inp) - ECP exponents * 46* coef_c (inp) - ECP contraction coefficients * 47* p_min (inp) - minimum power of r in ECP expansion * 48* p_max (inp) - maximum power of r in ECP expansion * 49* tol (inp) - maximum relative error in bessel functions * 50* sphcart (inp) - 0 for cartesian integrals, 1 for spherical * 51* tmp (scr) - work array * 52* l_tmp (inp) - length of tmp * 53* ecp_ints (out) - integrals over ECP * 54* n_int (inp) - number of ECP integrals * 55* n_blk (inp) - 1 for scalar only, 3 for s-o only, 4 for both * 56* transpose (inp) - true if centres A and B are to be transposed. * 57* ibug (inp) - debug flag. 0 for no debug, 1 for address printing, * 58* 2 for array printing, 3 for both. * 59* * 60* Notes: * 61* ----- * 62* * 63* The ECP centre is centre C. Centre B is assumed to coincide with * 64* centre C. * 65* The integrals come out in the order cont_a, cont_b, cmpt_a, cmpt_b * 66* where cont = contracted functions, cmpt = cartesian components * 67* The integrals are added to the array ecp_ints, i.e. ecp_ints is * 68* incremented by the integrals from this routine. * 69* * 70* Written by K. G. Dyall * 71* * 72************************************************************************ 73 implicit none 74#include "stdio.fh" 75#include "ecp_consts.fh" 76#include "util.fh" 77#include "errquit.fh" 78 integer l_a,n_prim_a,n_cont_a,l_b,n_prim_b,n_cont_b,l_ecp, 79 & l_ecp_max,n_z_tot,n_c_tot,n_cent,l_tmp,n_int,p_min,p_max, 80 & sphcart,mem_max,ibug 81 integer n_coef_c(-1:l_ecp_max,n_cent,2), 82 & n_prim_c(p_min:p_max,-1:l_ecp_max,n_cent,2), 83 & ind_z(-1:l_ecp_max,n_cent,2),ind_c(-1:l_ecp_max,n_cent,2) 84 integer i,i_a,i_b,i_c,i_ai,i_ang,i_ang_a,i_ca,i_Ga,i_free,i_ind, 85 & i_cont_c,i_co,i_off,i_int_a,i_k,i_t,i_x,i_y,i_z,i_Q_int,i_Qa, 86 & i_Qabc,i_gam,i_alp,i_pre,i_tmp,i_wrk,i_coef,k,l,l_c,l_c_min, 87 & l_c_max,l_sa,l_max,l_min,ll,l_lo,m_count,m_skip,n,n_na,n_nsa, 88 & n_ta,n_all_a,n_nb,n_all_b,n_abc,n_ab,nc_ab,ncab,n_rad,n_l, 89 & n_l_c,n_ang_a,n_ang_t,n_Qa,n_blk,n_co,n_co_tot,n_co_max, 90 & n_cont_c,n_x,n_pass 91 integer lcsco 92 logical DryRun,transpose,debug_gen,debug_addresses,debug_arrays, 93 & skip 94 double precision zeta_c(n_z_tot),coef_c(n_c_tot), 95 & coef_a(n_prim_a,n_cont_a),coef_b(n_prim_b,n_cont_b), 96 & zeta_a(n_prim_a),zeta_b(n_prim_b), 97 & tmp(l_tmp),ecp_ints(n_int,n_blk), 98 & R_AC,X_AC,Y_AC,Z_AC,R_BC,X_BC,Y_BC,Z_BC, 99 & tol,fac,log_prefactor 100 double precision csco(lcsco) 101* 102 logical ecp_skipint 103 external ecp_skipint 104* 105 debug_gen = ibug .gt. 0 106 debug_addresses = mod(ibug,2) .eq. 1 107 debug_arrays = (mod(ibug,10)/2 .eq. 1) .and. .not.DryRun 108* 109 if (DryRun) mem_max = 0 110 if ((sphcart .ne. 0) .and. (l_b .gt. l_ecp)) return 111* 112 if (debug_gen) write (LuOut,'(//A,/)') 'Entering ecp_int1 ...' 113 if (debug_gen) write (LuOut,*) 'ibug =',ibug 114 if (debug_addresses) then 115 write (LuOut,*) p_min,p_max,l_ecp 116 write (LuOut,*) ((n_prim_c(i,k,1,1),i = p_min,p_max), 117 & k = -1,l_ecp) 118 end if 119* 120* Set up angular momentum ranges and check whether to do block 121* 122 if (sphcart .eq. 0) then 123 l_c_min = mod(l_b,2) 124 l_c_max = l_b-2*(max(l_b-l_ecp+1,0)/2) 125 else 126 l_c_min = l_b 127 l_c_max = l_b 128 end if 129 if (debug_addresses) write (LuOut,*) 'l_c_min,l_c_max', 130 & l_c_min,l_c_max 131 l_max = l_a+l_b 132 l_min = l_b 133* 134* Allocate memory for ecp-independent quantities 135* 136 n_na = (l_a+1)*(l_a+2)/2 137 n_all_a = n_na*(l_a+3)/3 138 l_sa = l_ecp+l_a 139 n_nsa = (l_sa+1)**2 140 n_ta = (l_sa+1)*(l_sa+2)/2 141 n_nb = (l_b+1)*(l_b+2)/2 142 if (debug_addresses) 143 & write (LuOut,*) 'n_na,n_all_a,l_sa,n_nsa,n_ta,n_nb', 144 & n_na,n_all_a,l_sa,n_nsa,n_ta,n_nb 145* 146 n_ab = n_prim_a*n_prim_b 147 nc_ab = n_prim_a*n_cont_b 148 ncab = n_cont_a*n_cont_b 149 if (debug_addresses) 150 & write (LuOut,*) 'n_ab,nc_ab,ncab',n_ab,nc_ab,ncab 151* 152 i_ca = 1 153 i_Ga = i_ca+n_na*n_all_a 154 i_x = i_Ga+n_nsa 155 i_y = i_x+l_sa+1 156 i_z = i_y+l_sa+1 157 i_t = i_z+l_sa+1 158 i_free = i_t+n_ta 159 if (debug_addresses) then 160 write (LuOut,*) 'i_ca,i_Ga,i_x,i_y,i_z,i_t,i_free', 161 & i_ca,i_Ga,i_x,i_y,i_z,i_t,i_free 162 end if 163 if (DryRun) then 164 mem_max = max(mem_max,i_free-1) 165 if (debug_addresses) write (LuOut,*) 'mem_max',mem_max 166 else 167 if (i_free-1 .gt. l_tmp) call errquit( 168 & ' Insufficient memory in ecp_int1',99, MEM_ERR) 169* 170* Expand cartesian basis about ECP centre in spherical tensors 171* 172 if (debug_arrays) 173 & write (LuOut,*) 'X_AC,Y_AC,Z_AC',X_AC,Y_AC,Z_AC 174 call ecp_cart_xpd (l_a,n_na,n_all_a,X_AC,Y_AC,Z_AC, 175 & tmp(i_x),tmp(i_y),tmp(i_z),tmp(i_t),tmp(i_ca),1, 176 & csco,lcsco) 177 if (debug_arrays) call ecp_matpr(tmp(i_ca),1,n_na,1,n_all_a, 178 & 1,n_na,1,n_all_a,'Cartesian expansion','E',78,4) 179* 180* Set up spherical tensors which multiply bessel functions 181* 182 call ecp_sph_tens (l_sa,n_nsa,n_ta,R_AC,X_AC,Y_AC,Z_AC, 183 & tmp(i_x),tmp(i_y),tmp(i_z),tmp(i_t),tmp(i_Ga), 184 & csco,lcsco) 185 if (debug_arrays) call ecp_matpr(tmp(i_Ga),1,n_nsa,1,1, 186 & 1,n_nsa,1,1,'Spherical tensors','E',78,4) 187 end if 188* 189* Loop over ECP angular momentum projectors 190* 191 do l_c = l_c_max,l_c_min,-2 192* 193* Define parameters for the different integral class cases 194* 195 call ecp_ind_var (l_c,n_blk,n_coef_c(l_c,1,1),n_coef_c(l_c,1,2), 196 & ind_c(l_c,1,1),ind_c(l_c,1,2),ind_z(l_c,1,1),ind_z(l_c,1,2), 197 & n_x,n_co_tot,n_co_max,i_off,n_pass,i_cont_c,n_cont_c,skip) 198 if (skip) go to 10 199* 200* Check magnitude of integrals 201* 202 if (.not.DryRun) then 203 if (n_blk .eq. 1) then 204 if (ecp_skipint ( 205 & l_a,n_prim_a,n_cont_a,coef_a,zeta_a,R_AC, 206 & l_b,n_prim_b,n_cont_b,coef_b,zeta_b,R_BC, 207 & n_coef_c(l_c,1,1),zeta_c(ind_z(l_c,1,1)), 208 & coef_c(ind_c(l_c,1,1)))) go to 10 209 else if (n_blk .eq. 3) then 210 if (ecp_skipint ( 211 & l_a,n_prim_a,n_cont_a,coef_a,zeta_a,R_AC, 212 & l_b,n_prim_b,n_cont_b,coef_b,zeta_b,R_BC, 213 & n_coef_c(l_c,1,2),zeta_c(ind_z(l_c,1,2)), 214 & coef_c(ind_c(l_c,1,2)))) go to 10 215 else 216 if (ecp_skipint ( 217 & l_a,n_prim_a,n_cont_a,coef_a,zeta_a,R_AC, 218 & l_b,n_prim_b,n_cont_b,coef_b,zeta_b,R_BC, 219 & n_coef_c(l_c,1,1),zeta_c(ind_z(l_c,1,1)), 220 & coef_c(ind_c(l_c,1,1))) .and. 221 & ecp_skipint ( 222 & l_a,n_prim_a,n_cont_a,coef_a,zeta_a,R_AC, 223 & l_b,n_prim_b,n_cont_b,coef_b,zeta_b,R_BC, 224 & n_coef_c(l_c,1,2),zeta_c(ind_z(l_c,1,2)), 225 & coef_c(ind_c(l_c,1,2)))) go to 10 226 end if 227 end if 228* 229* Set up array dimensions for integrals 230* 231 ll = min(l_a,l_c) 232 n_abc = n_ab*n_co_max 233 n_rad = (ll+1)*(l_a+1)-ll*(ll+1)/2 234 if (debug_addresses) write (LuOut,*) 235 & 'l_a,l_c,ll,n_rad',l_a,l_c,ll,n_rad 236 if (sphcart .eq. 0) then 237 m_count = (l_a+l_b)/2 238 ll = max(l_a-l_c,0) 239 l = ll/2 240 n_rad = n_rad+l*(ll-l) 241 else 242 m_count = ll 243 end if 244 m_skip = (l_b-l_c)/2 245 l_lo = l_min-m_skip 246C write (LuOut,*) 'm_count',m_count 247 n_Qa = l_a+m_skip+1 248* 249* Allocate scratch memory for integrals 250* 251 i_Q_int = i_free 252 i_Qa = i_Q_int+ncab*n_rad*n_cont_c 253 i_Qabc = i_Qa+n_ab*n_Qa 254 i_ai = i_Qabc+n_abc 255 i_gam = i_ai+n_abc 256 i_alp = i_gam+n_abc 257 i_pre = i_alp+n_abc 258 i_tmp = i_pre+n_abc 259 i_ind = i_tmp+n_abc*6 260 i_wrk = i_ind+n_abc 261 i_coef = i_wrk+n_ab 262 i_free = i_coef+n_co_tot 263 if (debug_addresses) then 264 write (LuOut,*) 'i_Q_int,i_Qa,i_Qabc,i_ai,i_gam,i_alp,i_pre,', 265 & 'i_tmp,i_ind,i_wrk,i_free' 266 write (LuOut,*) i_Q_int,i_Qa,i_Qabc,i_ai,i_gam,i_alp,i_pre, 267 & i_tmp,i_ind,i_wrk,i_free 268 end if 269 if (DryRun) then 270 mem_max = max(mem_max,i_free-1) 271 if (debug_addresses) write (LuOut,*) 'mem_max',mem_max 272 else 273 if (i_free-1 .gt. l_tmp) call errquit( 274 & ' Insufficient memory in ecp_int1',99, MEM_ERR) 275 i_t = i_cont_c 276 i_x = i_Q_int 277 do i_y = 1,n_pass 278 i_z = ind_z(l_c,1,i_t) 279 i_co = ind_c(l_c,1,i_t) 280 n_co = n_coef_c(l_c,1,i_t) 281* 282* Set up argument values etc 283* 284 i = 0 285 do i_c = 0,n_co-1 286 do i_b = 1,n_prim_b 287 do i_a = 1,n_prim_a 288 tmp(i_gam+i) = one/sqrt(zeta_c(i_z+i_c)+zeta_b(i_b) 289 & +zeta_a(i_a)) 290 tmp(i_alp+i) = R_ac*zeta_a(i_a)*tmp(i_gam+i) 291 tmp(i_ai+i) = one/(two*R_ac*zeta_a(i_a)) 292 log_prefactor = tmp(i_alp+i)**2-zeta_a(i_a)*R_ac**2 293 tmp(i_pre+i) = exp(log_prefactor) 294 if (debug_arrays) write (LuOut,*) i,tmp(i_gam+i), 295 & tmp(i_alp+i),tmp(i_ai+i),tmp(i_pre+i) 296 i = i+1 297 end do 298 end do 299 end do 300 call dcopy (n_co,coef_c(i_co),1,tmp(i_coef),1) 301 if (n_cont_c .eq. 2) call dcopy 302 & (n_co,coef_c(ind_c(l_c,1,2)),1,tmp(i_coef+n_co),1) 303 n_abc = n_ab*n_co 304* 305* Calculate radial integrals 306* 307 call ecp_radint1 (p_min,p_max, 308 & l_lo,l_min,l_max,m_count,m_skip, 309 & n_prim_a,n_cont_a,coef_a,n_prim_b,n_cont_b,coef_b, 310 & n_co,n_prim_c(p_min,l_c,1,i_cont_c),n_cont_c, 311 & tmp(i_coef),tmp(i_ai),tmp(i_gam),tmp(i_alp),tmp(i_pre), 312 & tol,sphcart, 313 & n_ab,nc_ab,n_abc,n_rad,tmp(i_tmp),tmp(i_ind),tmp(i_wrk), 314 & tmp(i_Qabc),tmp(i_Qa),tmp(i_x),transpose,ibug/10) 315 if (debug_arrays) call ecp_matpr(tmp(i_x),1,ncab,1,n_rad, 316 & 1,ncab,1,n_rad,'Radial integrals','E',78,4) 317 i_t = i_t+1 318 i_x = i_x+ncab*n_rad 319 end do 320 end if 321* 322* Allocate memory for the contraction of radial and angular parts 323* 324 n_l_c = 2*l_c+1 325 ll = min(l_a,l_c) 326 n_ang_a = (ll+1)*(l_a+1)**2-ll*(ll+1)*(2*ll+1)/6 327 n_ang_t = n_ang_a*n_l_c 328 if (debug_addresses) write (LuOut,*) 'n_ang_a',n_ang_a 329 i_ang = i_Qa 330 i_ang_a = i_ang+max(n_ang_a*n_nb*n_x,n_ang_t) 331 i_z = i_ang_a 332 i_free = i_ang_a+max(n_ang_t*n_x,n_na*n_nb) 333 if (DryRun) then 334 mem_max = max(mem_max,i_free-1) 335 if (debug_addresses) write (LuOut,*) 'mem_max',mem_max 336 else 337 if (i_free-1 .gt. l_tmp) call errquit( 338 & ' Insufficient memory in ecp_int1',99, MEM_ERR) 339* 340* Set up angular coefficients for centre A and contract over 341* components of spherical tensors (sum over q). 342* 343 i = i_ang_a 344 do l = 0,l_a 345 n_l = 2*l+1 346 do k = l+l_c,abs(l-l_c),-2 347 i_k = k**2 348C write (LuOut,*) i,loc(tmp(i)) 349 call ecp_angint (tmp(i),l,k,l_c,tmp(i_Ga+i_k)) 350 i = i+n_l_c*n_l 351 end do 352 end do 353 if (debug_arrays) call ecp_matpr (tmp(i_ang_a),-l_c,l_c,1, 354 & n_ang_a,-l_c,l_c,1,n_ang_a,'Angular integrals','E',78,4) 355* 356* Perform sum over m. This involves transformation of the 357* cartesian functions to a spherical basis on centre B, followed 358* by evaluation of an overlap integral, which gives a factor of 359* 4\pi/(2\ell+1). This factor cancels with the factor from the 360* projector, and the cartesian to spherical transformation 361* transfers to the other angular integral, leaving a factor of 362* 2\pi(1+\delta_{m,0}). For the spin-orbit integrals, the matrix 363* element must also be inserted, but this is done in ecp_so_ang. 364* 365 if (n_x .gt. 1) then 366 call dcopy (n_ang_t,tmp(i_ang_a),1,tmp(i_ang),1) 367 i = i_ang_a+(n_x-3)*n_ang_t 368 call ecp_so_ang (l_c,tmp(i_ang),tmp(i),n_ang_a) 369 if (debug_arrays) then 370 do n = 1,3 371 call ecp_matpr (tmp(i),-l_c,l_c,1,n_ang_a,-l_c,l_c, 372 & 1,n_ang_a,'Spin-orbit angular integrals','E',78,4) 373 i = i+n_ang_t 374 end do 375 end if 376 end if 377 if (n_x .ne. 3) call dscal (n_ang_a,two,tmp(i_ang_a+l_c), 378 & n_l_c) 379 fac = pi+pi 380 call dscal (n_ang_t*n_x,fac,tmp(i_ang_a),1) 381 if (debug_arrays) call ecp_matpr (tmp(i_ang_a),-l_c,l_c,1, 382 & n_ang_a,-l_c,l_c,1,n_ang_a,'Scaled angular integrals', 383 & 'E',78,4) 384 call ecp_cstrans (l_b,n_nb,n_ang_a*n_x,l_c,l_c,l,tmp(i_ang), 385 & n_nb,tmp(i_ang_a),n_l_c,csco,lcsco,csco,-1,-1,1) 386 n_ang_t = n_ang_a*n_nb 387 if (debug_arrays) then 388 i = i_ang 389 do n = 1,n_x 390 call ecp_matpr (tmp(i),1,n_nb,1,n_ang_a, 391 & 1,n_nb,1,n_ang_a,'Ang ints summed over m','E',78,4) 392 i = i+n_ang_t 393 end do 394 end if 395* 396* Now loop over angular momenta of expanded function a and 397* perform contraction of angular intgrals with radial integrals 398* and expansion coefficients 399* 400 do i = 1,n_x 401 i_x = i+i_off 402 i_a = i_ca 403 i_c = i_Q_int 404 do n = 0,l_a 405 i_Qa = i_c 406 if (sphcart .eq. 0) then 407 l_lo = mod(n,2) 408 else 409 l_lo = n 410 end if 411 do l = n,l_lo,-2 412 if (debug_gen) write (LuOut,*) 'n,l',n,l 413 ll = min(l-1,l_c) 414 i_ang_a = i_ang+n_nb*((ll+1)*l**2-ll*(ll+1)*(2*ll+1)/6) 415 n_l = 2*l+1 416 i_int_a = i_Qa 417 do k = l+l_c,abs(l-l_c),-2 418 if (debug_addresses) then 419 write (LuOut,*) 'ang coef',i_ang_a-i_ang+1 420 write (LuOut,*) 'cart exp',i_a-i_ca+1 421 write (LuOut,*) 'integral',i_int_a-i_Q_int+1 422 end if 423 fac = 2*k+1 424 if (transpose) then 425 if (i .gt. abs(n_x-3)) fac = -fac 426 call dgemm ('N','T',n_nb,n_na,n_l,fac, 427 & tmp(i_ang_a),n_nb,tmp(i_a),n_na,zero, 428 & tmp(i_z),n_nb) 429 if (debug_arrays) call ecp_matpr (tmp(i_z),1,n_nb, 430 & 1,n_na,1,n_nb,1,n_na,'Completed angular ints', 431 & 'F',78,4) 432 call ecp_angrad (n_nb,n_cont_b,n_na,n_cont_a, 433 & tmp(i_z),tmp(i_int_a),ecp_ints(1,i_x)) 434 else 435 call dgemm ('N','T',n_na,n_nb,n_l,fac, 436 & tmp(i_a),n_na,tmp(i_ang_a),n_nb,zero, 437 & tmp(i_z),n_na) 438 if (debug_arrays) call ecp_matpr (tmp(i_z),1,n_na, 439 & 1,n_nb,1,n_na,1,n_nb,'Completed angular ints', 440 & 'F',78,4) 441 call ecp_angrad (n_na,n_cont_a,n_nb,n_cont_b, 442 & tmp(i_z),tmp(i_int_a),ecp_ints(1,i_x)) 443 end if 444 i_ang_a = i_ang_a+n_l*n_nb 445 i_int_a = i_int_a+ncab 446 end do 447 i_a = i_a+n_l*n_na 448 i_Qa = i_Qa+ncab 449 end do 450 if (sphcart .eq. 0) then 451 i_c = i_c+(min((n+l_c)/2,n)+1)*ncab 452 else 453 i_c = i_c+(min(n,l_c)+1)*ncab 454 end if 455 end do 456 i_ang = i_ang+n_ang_t 457 if ((n_x .eq. 4) .and. (i .eq. 1)) 458 & i_Q_int = i_Q_int+ncab*n_rad 459 end do 460* 461 end if 462 10 continue 463 end do 464 if (debug_arrays) then 465 n_all_a = n_na*n_cont_a 466 n_all_b = n_nb*n_cont_b 467 do i = 1,n_blk 468 call ecp_matpr (ecp_ints(1,i),1,n_all_b,1,n_all_a, 469 & 1,n_all_b,1,n_all_a,'ECP integrals','E',78,4) 470 end do 471 end if 472 if (debug_gen) write (LuOut,*) 'Exiting ecp_int1' 473* 474 return 475 end 476 477