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