1c 2c 3c $Id$ 4c 5c 6c This routine returns the Coulomb and exchange integral 7c operator matrices for the range of MO-indices as mo_indx_hi, mo_indx_lo 8c The g_coul, g_exch global arrays are ordered as 9c 10c 11c 12c ij 13c (ij|ab) = ( J ) = g_coul[ ij : (a-1)*N2 + b ] = g_coul [ ij : (b-1)*N2 + a ] 14c ab 15c 16c ij 17c (ia|jb) = ( K ) = g_exch[ ij : (a-1)*N2 + b ] 18c ab 19c 20c ----------------- 21c Note: This routine differs from the standard in that 22c it computes the integrals 6 times as opposed to 2 times 23c in the standard. However, it has lower communication scaling 24c for K, O(N^4), compared to O(N^5) in the standard. Therefore 25c this routine is more appropriate for very high parallelism 26c ----------------- 27c 28c 29c 30 subroutine moints_build_6x(basis, ousesym, 31 $ occ_start, mo1_lo, mo1_hi, 32 $ mo2_lo, mo2_hi, 33 $ g_movecs, 34 $ g_coul, ocoul, 35 $ g_exch, oexch, 36 $ blen, oblk ) 37 implicit none 38#include "errquit.fh" 39#include "tcgmsg.fh" 40#include "global.fh" 41#include "mafdecls.fh" 42#include "bas.fh" 43#include "util.fh" 44#include "sym.fh" 45#include "schwarz.fh" 46#include "msgids.fh" 47c 48c Arguments 49c 50 integer basis ! Basis handle 51 logical ousesym ! Symmetry toggle 52 integer occ_start ! Offset from frozen core 53 integer mo1_lo, mo1_hi ! 1st Pair Index range 54 integer mo2_lo, mo2_hi ! 2nd Pair Index range 55 integer g_movecs ! MO coefficients 56 integer g_coul ! Coulomb operator 57 integer g_exch ! Exchange operator 58 logical ocoul,oexch ! Type selection 59 integer blen ! Blocking length 60 logical oblk ! Toggle AO blocking 61c 62c Local variables 63c 64 integer geom, nmo, nmo1, nmo2, nbf, nsh, maxbfsh 65 integer qlo, qhi 66 integer nblk, bsize, ngrp 67 integer ii, jj, iz, kz, jz 68 integer l_gmap, k_gmap 69 integer ish, jsh, ilen, jlen 70 integer ibflo,ibfhi,jbflo,jbfhi,kbflo,kbfhi,lbflo,lbfhi 71 integer kshlo, kshhi, lshlo, lshhi 72 integer kblen, lblen 73 integer kgr, lgr 74 integer i, nmixed 75 integer l_ssbb, k_ssbb, l_ssbbt, k_ssbbt 76 integer l_hlp, k_hlp, l_ssni,k_ssni 77 integer l_eri, k_eri, l_iscr,k_iscr 78 integer l_mo, k_mo, l_xmo, k_xmo 79 integer l_shmap, k_shmap, l_bfmap, k_bfmap, l_rbfmap, k_rbfmap 80 integer l_glo, k_glo, l_ghi, k_ghi 81 integer l_sym, k_sym 82 integer n_ssbb, n_ssbb1, n_ssni, n_hlp, n_ijni 83 integer mem2, max2e 84 integer num_nodes, ploop, next 85 double precision tol2e, scale, scale0, schw_ij, integ_acc 86 double precision tz, thalf, tint, t1jidx, t2jidx, t1kidx, t2kidx 87 double precision t34kidx, t34jidx, thalf0, ttotal, tsynch 88 double precision ttask, ttaskmax, ttaskmin, ttaskagg 89 double precision flop1, flop2k, flop2j, q2 90 integer tottask 91 logical status, st, osym, odoit 92#include "moints_stats.fh" 93#include "moints.fh" 94c 95c 96c 97c integer moints_numgr, gr_len, moints_nxttask 98c external moints_numgr, gr_len, moints_nxttask 99 integer gr_len 100 external gr_len 101c 102c 103c 104 data tol2e/1.d-12/ 105c 106c 107c 108 integ_acc = min(1d-9, max(0.01d0*tol2e,1d-17)) 109 call int_acc_set(integ_acc) 110c 111c General basis info 112c 113 ttotal = util_cpusec() 114 num_nodes = ga_nnodes() 115 if (.not. bas_geom(basis, geom)) call errquit 116 $ ('moints: cannot get geometry', 0, GEOM_ERR) 117 status = bas_numbf(basis,nbf) 118 status = status.and.bas_numcont(basis,nsh) 119 status = status.and.bas_nbf_cn_max(basis,maxbfsh) 120 if (.not.status) call errquit('moints: cannot get basis info',0, 121 & BASIS_ERR) 122 qlo = min(occ_start,mo2_lo,mo1_lo) 123 qhi = max(mo1_hi,mo2_hi) 124 nmo = qhi - qlo + 1 125 nmo1 = mo1_hi - mo1_lo + 1 126 nmo2 = mo2_hi - mo2_lo + 1 127c 128c Symmetry adapt the MOs and renumber irreps to start at zero 129c 130 osym = (ousesym.and.(sym_number_ops(geom).gt.0)) 131 if (osym) then 132 if (.not. ma_push_get(MT_INT, nbf, 'movec syms',l_sym, k_sym)) 133 $ call errquit('moints_6x: no memory for syms?',0, 134 & MA_ERR) 135 call sym_movecs_adapt(basis, 1d-8, g_movecs, int_mb(k_sym), 136 $ nmixed) 137 if (nmixed .ne. 0) 138 $ call errquit('moints_6x: symmetry contamination', nmixed, 139 & INPUT_ERR) 140 do i =0, nbf-1 141 int_mb(k_sym+i) = int_mb(k_sym+i) - 1 142 enddo 143 if (util_print('orbital symmetry',print_debug)) then 144 write(6,887) 145 887 format('Symmetry of MOs') 146 write(6,888) (int_mb(k_sym+i),i=0,nbf-1) 147 888 format(16i3) 148 endif 149 endif 150c 151c Initialize AO disk caching 152c 153 status = status .and. moints_aodisk_prep() 154c 155c Local MO coefficients 156c 157 status = status.and.ma_push_get(MT_DBL,(nbf*nmo), 158 $ 'movecs cols',l_mo,k_mo) 159 call ga_get(g_movecs, 1, nbf, qlo, qhi, dbl_mb(k_mo), nbf) 160c 161c Integrals allocation 162c 163 if (oblk) then 164 call intb_mem_2e4c( max2e, mem2 ) 165 max2e = max( max2e, maxbfsh*maxbfsh*blen*blen ) 166 else 167 call int_mem_2e4c( max2e, mem2) 168 endif 169 status = ma_push_get(MT_DBL, max2e,'moints: buf', l_eri, k_eri) 170 status = ma_push_get(MT_DBL, mem2, 'moints: scr', l_iscr, k_iscr) 171c 172c Shell group mapping 173c 174 nblk = moints_numgr( basis, blen ) 175 status = ma_push_get(MT_INT,(nblk*4),'shell group map', 176 $ l_gmap, k_gmap) 177 call moints_grmap( basis, blen, nblk, int_mb(k_gmap)) 178c 179c Reorder and group shells by type 180c 181 status = ma_push_get(MT_INT,nsh,'shell order map', 182 $ l_shmap, k_shmap) 183 status = ma_push_get(MT_INT,nsh,'group lo', l_glo, k_glo ) 184 status = ma_push_get(MT_INT,nsh,'group hi', l_ghi, k_ghi) 185 status = ma_push_get(MT_INT,nbf,'basis map', 186 $ l_bfmap, k_bfmap) 187 status = ma_push_get(MT_INT,nbf,'rev basis map', 188 $ l_rbfmap, k_rbfmap) 189 call moints_shorder( basis, nsh, nbf, blen, ngrp, 190 $ int_mb(k_glo), int_mb(k_ghi), 191 $ int_mb(k_shmap), 192 $ int_mb(k_bfmap), int_mb(k_rbfmap) ) 193c 194c Copy of MO coefficients with reordered rows 195c 196 status = ma_push_get(MT_DBL,(nbf*nmo),'reorder mos', 197 $ l_xmo, k_xmo) 198 call row_exch( nbf, nmo, int_mb(k_rbfmap), dbl_mb(k_mo), 199 $ dbl_mb(k_xmo) ) 200c 201c Temporary partially-transformed arrays 202c 203 bsize = max(blen,maxbfsh) 204 n_ssbb = maxbfsh*maxbfsh*bsize*bsize 205 n_ssbb1 = max((nmo1*nmo1),n_ssbb) 206 n_hlp = max((bsize*maxbfsh*maxbfsh*nmo1),(maxbfsh*nbf)) 207 n_ssni = maxbfsh*maxbfsh*nbf*nmo1 208 status = ma_push_get(MT_DBL,n_ssbb1,'ssbb block',l_ssbb,k_ssbb) 209 status = ma_push_get(MT_DBL,n_ssbb,'ssbbt block',l_ssbbt,k_ssbbt) 210 status = ma_push_get(MT_DBL,n_hlp,'hlp block',l_hlp,k_hlp) 211 status = ma_push_get(MT_DBL,n_ssni,'ssni block',l_ssni,k_ssni) 212 if (.not.(status)) call errquit 213 $ ('moints_6x: cannot allocate local memory',0, MA_ERR) 214c 215c Initialize 216c 217 flop1 = 0.d0 218 flop2k = 0.d0 219 flop2j = 0.d0 220 tint = 0.d0 221 t1jidx = 0.d0 222 t2jidx = 0.d0 223 t1kidx = 0.d0 224 t2kidx = 0.d0 225 t34jidx = 0.d0 226 t34kidx = 0.d0 227 tottask = 0 228 ttaskmax = 0.d0 229 ttaskmin = 1.d24 230 ttaskagg = 0.d0 231 thalf0 = util_cpusec() 232 ploop = 0 233 next = -1 234 if (ocoul) call ga_zero(g_coul) 235 if (oexch) call ga_zero(g_exch) 236c 237c Exchange part 238c Four fold shell loop 239c 240c 241 if (oexch) then 242 if (next.lt.0) next = moints_nxttask(num_nodes) 243 244 do ii=1,nsh 245 ish = int_mb(k_shmap+ii-1) 246 status = bas_cn2bfr(basis,ish,ibflo,ibfhi) 247 ilen = ibfhi - ibflo + 1 248 do jj=1,ii 249 jsh = int_mb(k_shmap+jj-1) 250 status = bas_cn2bfr(basis,jsh,jbflo,jbfhi) 251 jlen = jbfhi - jbflo + 1 252 scale0 = 1.d0 253 schw_ij = schwarz_shell(ish,jsh) 254 if (ish.eq.jsh) scale0 = 0.5d0 255 if (next.eq.ploop) then 256 ttask = util_cpusec() 257 tottask = tottask + 1 258 call dfill(n_ssni,0.d0,dbl_mb(k_ssni),1) 259 do kgr=1,ngrp 260 kshlo = int_mb(k_glo+kgr-1) 261 kshhi = int_mb(k_ghi+kgr-1) 262 st = bas_cn2bfr(basis,int_mb(k_shmap+kshlo-1),iz,kz) 263 st = bas_cn2bfr(basis,int_mb(k_shmap+kshhi-1),kz,jz) 264 kbflo = int_mb(k_rbfmap+iz-1) 265 kbfhi = int_mb(k_rbfmap+jz-1) 266 kblen = kbfhi - kbflo + 1 267 do lgr=1,kgr 268 lshlo = int_mb(k_glo+lgr-1) 269 lshhi = int_mb(k_ghi+lgr-1) 270 st = bas_cn2bfr(basis,int_mb(k_shmap+lshlo-1),iz,kz) 271 st = bas_cn2bfr(basis,int_mb(k_shmap+lshhi-1),kz,jz) 272 lbflo = int_mb(k_rbfmap+iz-1) 273 lbfhi = int_mb(k_rbfmap+jz-1) 274 lblen = lbfhi - lbflo + 1 275 tz = util_cpusec() 276 scale = 1.d0 277 if (kgr.eq.lgr) scale = 0.5d0 278 call moints6x_gblkk( basis, ish, jsh, kshlo, kshhi, 279 $ lshlo, lshhi, int_mb(k_shmap), 280 $ int_mb(k_rbfmap), 281 $ schw_ij, tol2e, osym, oblk, 282 $ max2e, dbl_mb(k_eri), 283 $ mem2, dbl_mb(k_iscr), 284 $ ibflo, ibfhi, jbflo, jbfhi, 285 $ kbflo, kbfhi, lbflo, lbfhi, 286 $ dbl_mb(k_ssbb), 287 $ dbl_mb(k_ssbbt) ) 288 tint = tint + util_cpusec() - tz 289 flop1 = flop1 + 4*lblen*ilen*jlen*kblen*nmo1 290 tz = util_cpusec() 291 call moints6x_trf1K(nbf, qlo, qhi, mo1_lo, mo1_hi, 292 $ ilen, jlen, kbflo, kbfhi, lbflo, 293 $ lbfhi, scale, dbl_mb(k_ssbb), 294 $ dbl_mb(k_ssbbt), dbl_mb(k_xmo), 295 $ dbl_mb(k_ssni), dbl_mb(k_hlp)) 296 t1kidx = t1kidx + util_cpusec() - tz 297 enddo 298 enddo 299 300 tz = util_cpusec() 301 call moints6x_trf2K( nbf, qlo, qhi, 302 $ occ_start, mo1_lo, mo1_hi, 303 $ ibflo, ibfhi, jbflo, jbfhi, 304 $ dbl_mb(k_ssni), 305 $ dbl_mb(k_hlp), dbl_mb(k_ssbb), 306 $ dbl_mb(k_hlp), dbl_mb(k_ssbb), 307 $ dbl_mb(k_xmo), g_exch ) 308 309 310 t2kidx = t2kidx + util_cpusec() - tz 311 next = moints_nxttask(num_nodes) 312 ttask = util_cpusec() - ttask 313 ttaskmax = max(ttaskmax,ttask) 314 ttaskmin = min(ttaskmin,ttask) 315 ttaskagg = ttaskagg + ttask 316 endif 317 ploop = ploop + 1 318 enddo 319 enddo 320 endif 321c 322 if (oexch) then 323 if (util_print('exchange half integral',print_debug)) then 324 call moints_op_print(occ_start,mo1_lo,mo1_hi,nbf,g_exch) 325 endif 326 endif 327c 328c Coulomb part 329c 4-fold shell loop 330c 331 if (ocoul) then 332 if (next.lt.0) next = moints_nxttask(num_nodes) 333 334 do ii=1,nsh 335 do jj=1,ii 336 ish = max(int_mb(k_shmap+ii-1),int_mb(k_shmap+jj-1)) 337 jsh = min(int_mb(k_shmap+ii-1),int_mb(k_shmap+jj-1)) 338 status = bas_cn2bfr(basis,ish,ibflo,ibfhi) 339 status = bas_cn2bfr(basis,jsh,jbflo,jbfhi) 340 ilen = ibfhi - ibflo + 1 341 jlen = jbfhi - jbflo + 1 342 schw_ij = schwarz_shell(ish,jsh) 343 scale = 1.d0 344 if (ish.eq.jsh) scale = scale*0.5d0 345 odoit = schw_ij*schwarz_max().ge.tol2e 346 if (odoit .and. osym) then 347 odoit = sym_shell_pair(basis, ish, jsh, q2) 348 endif 349 350 if (odoit) then 351 if (next.eq.ploop) then 352 ttask = util_cpusec() 353 tottask = tottask + 1 354 n_ijni = ilen*jlen*nbf*nmo1 355 call dfill(n_ijni,0.d0,dbl_mb(k_ssni),1) 356 357 do kgr=1,ngrp 358 kshlo = int_mb(k_glo+kgr-1) 359 kshhi = int_mb(k_ghi+kgr-1) 360 st = bas_cn2bfr(basis,int_mb(k_shmap+kshlo-1),iz,kz) 361 st = bas_cn2bfr(basis,int_mb(k_shmap+kshhi-1),kz,jz) 362 kbflo = int_mb(k_rbfmap+iz-1) 363 kbfhi = int_mb(k_rbfmap+jz-1) 364 kblen = kbfhi - kbflo + 1 365 do lgr=1,kgr 366 lshlo = int_mb(k_glo+lgr-1) 367 lshhi = int_mb(k_ghi+lgr-1) 368 st = bas_cn2bfr(basis,int_mb(k_shmap+lshlo-1),iz,kz) 369 st = bas_cn2bfr(basis,int_mb(k_shmap+lshhi-1),kz,jz) 370 lbflo = int_mb(k_rbfmap+iz-1) 371 lbfhi = int_mb(k_rbfmap+jz-1) 372 lblen = lbfhi - lbflo + 1 373 tz = util_cpusec() 374 call moints_gblk( basis, ish, jsh, 375 $ kshlo, kshhi, lshlo, lshhi, 376 $ int_mb(k_shmap),int_mb(k_rbfmap), 377 $ schw_ij, tol2e, osym, oblk, 378 $ max2e, dbl_mb(k_eri), 379 $ mem2, dbl_mb(k_iscr), 380 $ ibflo, ibfhi, jbflo, jbfhi, 381 $ kbflo, kbfhi, lbflo, lbfhi, 382 $ dbl_mb(k_ssbb) ) 383 tint = tint + util_cpusec() - tz 384 flop1 = flop1 + 4*lblen*ilen*jlen*kblen*nmo1 385 if (lgr.ne.kgr) then 386 tz = util_cpusec() 387 call moints_blktr(ilen, jlen, kblen, lblen, 388 $ dbl_mb(k_ssbb), 389 $ dbl_mb(k_ssbbt)) 390 call moints_trf1(nbf, qlo, qhi, mo1_lo, mo1_hi, 391 $ ilen, jlen, kbflo, kbfhi, 392 $ lbflo, lbfhi, 1.d0, 393 $ dbl_mb(k_ssbb), 394 $ dbl_mb(k_ssbbt), 395 $ dbl_mb(k_xmo), 396 $ dbl_mb(k_ssni), 397 $ dbl_mb(k_hlp)) 398 else 399 tz = util_cpusec() 400 call moints_trf1(nbf, qlo, qhi, mo1_lo, mo1_hi, 401 $ ilen, jlen, kbflo, kbfhi, 402 $ lbflo, lbfhi, 0.5d0, 403 $ dbl_mb(k_ssbb), 404 $ dbl_mb(k_ssbb), 405 $ dbl_mb(k_xmo), 406 $ dbl_mb(k_ssni), 407 $ dbl_mb(k_hlp)) 408 endif 409 t1jidx = t1jidx + util_cpusec() - tz 410 enddo 411 enddo 412 flop2j = flop2j+2*nbf*ilen*jlen*((nmo1*(nmo1+1))/2) 413 tz = util_cpusec() 414 call moints_trf2J(nbf, qlo, qhi, 415 $ occ_start, mo1_lo, mo1_hi, 416 $ ibflo, ibfhi, jbflo, jbfhi, 417 $ dbl_mb(k_ssni), dbl_mb(k_hlp), 418 $ dbl_mb(k_ssbb), dbl_mb(k_xmo), 419 $ g_coul) 420 t2jidx = t2jidx + util_cpusec() - tz 421 next = moints_nxttask(num_nodes) 422 ttask = util_cpusec() - ttask 423 ttaskmax = max(ttaskmax,ttask) 424 ttaskmin = min(ttaskmin,ttask) 425 ttaskagg = ttaskagg + ttask 426 endif 427 ploop = ploop + 1 428 endif 429 enddo 430 enddo 431 endif 432c 433 if (ocoul) then 434 if (util_print('coulomb half integral',print_debug)) then 435 call moints_op_print(occ_start,mo1_lo,mo1_hi,nbf,g_coul) 436 endif 437 endif 438c 439c Synchronize 440c 441 tz = util_cpusec() 442 next = moints_nxttask(-num_nodes) 443 tsynch = util_cpusec() - tz 444 thalf = util_cpusec() - thalf0 445 call moints_aodisk_tidy() 446c 447c Clean-up 448c 449 if (.not. ma_pop_stack(l_ssni)) 450 $ call errquit('moints: failed to pop', l_ssni, MA_ERR) 451 if (.not. ma_pop_stack(l_hlp)) 452 $ call errquit('moints: failed to pop', l_hlp, MA_ERR) 453 if (.not. ma_pop_stack(l_ssbbt)) 454 $ call errquit('moints: failed to pop', l_ssbbt, MA_ERR) 455 if (.not. ma_pop_stack(l_ssbb)) 456 $ call errquit('moints: failed to pop', l_ssbb, MA_ERR) 457 if (.not. ma_pop_stack(l_xmo)) 458 $ call errquit('moints: failed to pop', l_xmo, MA_ERR) 459 if (.not. ma_pop_stack(l_rbfmap)) 460 $ call errquit('moints: failed to pop', l_rbfmap, MA_ERR) 461 if (.not. ma_pop_stack(l_bfmap)) 462 $ call errquit('moints: failed to pop', l_bfmap, MA_ERR) 463 if (.not. ma_pop_stack(l_ghi)) 464 $ call errquit('moints: failed to pop', l_ghi, MA_ERR) 465 if (.not. ma_pop_stack(l_glo)) 466 $ call errquit('moints: failed to pop', l_glo, MA_ERR) 467 if (.not. ma_pop_stack(l_shmap)) 468 $ call errquit('moints: failed to pop', l_shmap, MA_ERR) 469 if (.not. ma_pop_stack(l_gmap)) 470 $ call errquit('moints: failed to pop', l_gmap, MA_ERR) 471 if (.not. ma_pop_stack(l_iscr)) 472 $ call errquit('moints: failed to pop', l_iscr, MA_ERR) 473 if (.not. ma_pop_stack(l_eri)) 474 $ call errquit('moints: failed to pop', l_eri, MA_ERR) 475c 476c 477c 478 status = ma_push_get(MT_DBL,(nbf*nbf),'hlp',l_hlp,k_hlp) 479c 480c Exchange 3 & 4 qtr transformation 481c 482 if (oexch) then 483 tz = util_cpusec() 484 call moints_Ktrf34(nbf, qlo, qhi, occ_start, mo1_lo, mo1_hi, 485 $ mo2_lo, mo2_hi, .false., dbl_mb(k_mo), 486 $ dbl_mb(k_hlp), osym, int_mb(k_sym), 487 $ g_exch) 488 t34kidx = util_cpusec() - tz 489 if (util_print('exchange integral',print_debug)) then 490 call moints_op_print(occ_start,mo1_lo,mo1_hi,nmo2,g_exch) 491 endif 492C call MOINTS_OP_PRINT(OCC_START,MO1_LO,MO1_HI,NMO2,G_EXCH) 493 endif 494c 495c Coulomb 3 & 4 qtr transformation 496c 497 498 if (ocoul) then 499 tz = util_cpusec() 500 call moints_Jtrf34(nbf, qlo, qhi, occ_start, mo1_lo, mo1_hi, 501 $ mo2_lo, mo2_hi, dbl_mb(k_mo), 502 $ dbl_mb(k_hlp), osym, int_mb(k_sym), 503 $ g_coul) 504 if (util_print('coulomb integral',print_debug)) then 505 call moints_op_print(occ_start,mo1_lo,mo1_hi,nmo2,g_coul) 506 endif 507 t34jidx = util_cpusec() - tz 508 endif 509 call ga_sync() 510c 511c Clean-up 512c 513 if (.not. ma_pop_stack(l_hlp)) 514 $ call errquit('moints: failed to pop', l_hlp, MA_ERR) 515 if (.not. ma_pop_stack(l_mo)) 516 $ call errquit('moints: failed to pop', l_mo, MA_ERR) 517 if (osym) then 518 if (.not. ma_pop_stack(l_sym)) 519 $ call errquit('moints_2x: memory corrupt',0, MA_ERR) 520 endif 521c 522c 523c 524#ifdef BLOCK_TRANSF 525 if (ga_nodeid().eq.0) write(6,965) 526 965 format(/,10x,'**** BLOCK TRANSFER ENABLED ****') 527#endif 528#ifdef NOCOMMS 529 if (ga_nodeid().eq.0) write(6,334) 530 334 format(/,10x,'**** COMMUNICATION DISABLED ****') 531#endif 532 533 534c 535c Timings and Statistics bookkeeping 536c 537 flop1 = flop1*1.d-6 538 flop2k = flop2k*1.d-6 539 flop2j = flop2j*1.d-6 540 ttotal = util_cpusec() - ttotal 541 mi_npass = mi_npass + 1.d0 542 mi_ttotal = mi_ttotal + ttotal 543 mi_thalf = mi_thalf + thalf 544 mi_tint = mi_tint + tint 545 mi_t1j = mi_t1j + t1jidx 546 mi_t2j = mi_t2j + t2jidx 547 mi_t1k = mi_t1k + t1kidx 548 mi_t2k = mi_t2k + t2kidx 549 mi_t34k = mi_t34k + t34kidx 550 mi_t34j = mi_t34j + t34jidx 551 mi_flop1 = mi_flop1 + flop1 552 mi_synch = mi_synch + tsynch 553 mi_maxsynch = max(mi_maxsynch,tsynch) 554 mi_minsynch = min(mi_minsynch,tsynch) 555 mi_aggsynch = mi_aggsynch + tsynch 556 mi_nsynchs = mi_nsynchs + 1 557 mi_maxtask = max(mi_maxtask,ttaskmax) 558 mi_mintask = min(mi_mintask,ttaskmin) 559 mi_aggtask = mi_aggtask + ttaskagg 560 mi_ntasks = mi_ntasks + tottask 561 562 return 563 end 564 565 566 567 568 569 570 571 572 573 574c 575c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 576c 577c OLD VERSIONS WITH NO REORDERING 578c 579c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 580c 581c 582c 583c Generate AO integrals in suitable form 584c exchange transformation: 585c Branch according to block or non-block 586c scheme 587c 588c 589 590#ifdef NOCOMPILE 591 subroutine moints6x_gblkk_old( basis, ish, jsh, 592 $ kshlo, kshhi, lshlo, lshhi, 593 $ schw_ij, tol2e, erilen, eri, 594 $ scrlen, iscr, 595 $ ibflo, ibfhi, jbflo, jbfhi, 596 $ kblo, kbhi, lblo, lbhi, 597 $ sbsb, sbsbt, oblk ) 598 implicit none 599#include "errquit.fh" 600#include "mafdecls.fh" 601#include "bas.fh" 602#include "schwarz.fh" 603 integer basis, ish, jsh, kshlo, kshhi, lshlo, lshhi 604 integer erilen, scrlen 605 double precision schw_ij, tol2e, eri(erilen), iscr(scrlen) 606 integer ibflo, ibfhi, jbflo, jbfhi 607 integer kblo, kbhi, lblo, lbhi 608 double precision sbsb(lblo:lbhi,kblo:kbhi,jbflo:jbfhi,ibflo:ibfhi) 609 double precision sbsbt(kblo:kbhi,lblo:lbhi, 610 $ jbflo:jbfhi,ibflo:ibfhi) 611 logical oblk 612c 613 integer mxshq, imem 614 integer k_tmp, l_tmp 615 integer k_iv, k_jv, k_kv, k_lv 616 integer k_il, k_jl, k_kl, k_ll 617c 618 619 620 if (.not.(oblk)) then 621 call moints6x_gblkk_sq_old( basis, ish, jsh, 622 $ kshlo, kshhi, lshlo, lshhi, 623 $ schw_ij, tol2e, erilen, eri, 624 $ scrlen, iscr, 625 $ ibflo, ibfhi, jbflo, jbfhi, 626 $ kblo, kbhi, lblo, lbhi, 627 $ sbsb, sbsbt ) 628 else 629 mxshq = (kshhi - kshlo + 1)*(lshhi - lshlo + 1) 630 imem = 4*mxshq + 4*erilen 631 if (.not. ma_push_get(MT_INT, imem, 'gblk tmp',l_tmp, k_tmp)) 632 $ call errquit('moints6x_gblkk: no memory for labels',imem, 633 & MA_ERR) 634 k_iv = k_tmp 635 k_jv = k_tmp + mxshq 636 k_kv = k_tmp + 2*mxshq 637 k_lv = k_tmp + 3*mxshq 638 k_il = k_tmp + 4*mxshq 639 k_jl = k_tmp + 4*mxshq + erilen 640 k_kl = k_tmp + 4*mxshq + 2*erilen 641 k_ll = k_tmp + 4*mxshq + 3*erilen 642 call moints6x_gblkk_mq_old( basis, ish, jsh, 643 $ kshlo, kshhi, lshlo, lshhi, 644 $ schw_ij, tol2e, erilen, eri, 645 $ scrlen, iscr, 646 $ ibflo, ibfhi, jbflo, jbfhi, 647 $ kblo, kbhi, lblo, lbhi, 648 $ sbsb, sbsbt, mxshq, 649 $ int_mb(k_iv), int_mb(k_jv), 650 $ int_mb(k_kv), int_mb(k_lv), 651 $ int_mb(k_il), int_mb(k_jl), 652 $ int_mb(k_kl), int_mb(k_ll) ) 653 if (.not.ma_pop_stack(l_tmp)) 654 $ call errquit('moints6x_gblkk: failed to pop',0, MA_ERR) 655 endif 656 return 657 end 658 659 660 661 662 663c 664c This version uses single shell 665c quartet AO integrals. 666c 667 subroutine moints6x_gblkk_sq_old( basis, ish, jsh, 668 $ kshlo, kshhi, lshlo, lshhi, 669 $ schw_ij, tol2e, erilen, eri, 670 $ scrlen, iscr, 671 $ ibflo, ibfhi, jbflo, jbfhi, 672 $ kblo, kbhi, lblo, lbhi, 673 $ sbsb, sbsbt ) 674 implicit none 675#include "bas.fh" 676#include "schwarz.fh" 677 integer basis, ish, jsh, kshlo, kshhi, lshlo, lshhi 678 integer erilen, scrlen 679 double precision schw_ij, tol2e, eri(erilen), iscr(scrlen) 680 integer ibflo, ibfhi, jbflo, jbfhi 681 integer kblo, kbhi, lblo, lbhi 682 double precision sbsb(lblo:lbhi,kblo:kbhi,jbflo:jbfhi,ibflo:ibfhi) 683 double precision sbsbt(kblo:kbhi,lblo:lbhi, 684 $ jbflo:jbfhi,ibflo:ibfhi) 685c 686 integer ilen, jlen 687 integer ksh, lsh, kbflo, kbfhi, lbflo, lbfhi, ltop 688 integer klen, llen, kblen, lblen, bsize, k1, k2, j, i 689 logical status 690c 691c 692c 693 ilen = ibfhi - ibflo + 1 694 jlen = jbfhi - jbflo + 1 695 kblen = kbhi - kblo + 1 696 lblen = lbhi - lblo + 1 697 bsize = kblen*lblen*ilen*jlen 698 call dfill(bsize,0.d0,sbsb,1) 699 call dfill(bsize,0.d0,sbsbt,1) 700 do ksh=kshlo,kshhi 701 status = bas_cn2bfr(basis,ksh,kbflo,kbfhi) 702 klen = kbfhi - kbflo + 1 703 ltop = lshhi 704 if (kshlo.eq.lshlo) ltop = ksh 705 do lsh=lshlo,ltop 706C do lsh=lshlo,lshhi 707 status = bas_cn2bfr(basis,lsh,lbflo,lbfhi) 708 llen = lbfhi - lbflo + 1 709 if (schwarz_shell(ish,ksh)* 710 $ schwarz_shell(jsh,lsh).ge.tol2e) then 711 call int_2e4c(basis, ish, ksh, basis, jsh, lsh, 712 $ scrlen, iscr, erilen, eri ) 713 call moints6x_eri2blkk(ilen, jlen, klen, llen, 714 $ eri, sbsb(lbflo,kbflo,jbflo,ibflo), 715 $ lblen, kblen ) 716 endif 717 if (schwarz_shell(ish,lsh)* 718 $ schwarz_shell(jsh,ksh).ge.tol2e) then 719 call int_2e4c(basis, ish, lsh, basis, jsh, ksh, 720 $ scrlen, iscr, erilen, eri ) 721 call moints6x_eri2blkk(ilen, jlen, llen, klen, 722 $ eri, sbsbt(kbflo,lbflo,jbflo,ibflo), 723 $ kblen, lblen ) 724 endif 725 enddo 726 enddo 727 if (kshlo.eq.lshlo) then 728 do i=ibflo,ibfhi 729 do j=jbflo,jbfhi 730 do k1=kblo,kbhi 731 do k2=kblo,k1-1 732 sbsb(k1,k2,j,i) = sbsbt(k1,k2,j,i) 733 sbsbt(k2,k1,j,i) = sbsb(k2,k1,j,i) 734 enddo 735 enddo 736 enddo 737 enddo 738 endif 739c 740c 741c 742 return 743 end 744 745 746 747 748 749 750c 751c This uses multiple quartets at a time 752c AO integrals with labels 753c 754c 755c 756 subroutine moints6x_gblkk_mq_old( basis, ish, jsh, 757 $ kshlo, kshhi, lshlo, lshhi, 758 $ schw_ij, tol2e, erilen, eri, 759 $ scrlen, iscr, 760 $ ibflo, ibfhi, jbflo, jbfhi, 761 $ kblo, kbhi, lblo, lbhi, 762 $ sbsb, sbsbt, mxshq, 763 $ ishv, jshv, kshv, lshv, 764 $ ilab, jlab, klab, llab ) 765 implicit none 766#include "errquit.fh" 767#include "bas.fh" 768#include "schwarz.fh" 769 integer basis, ish, jsh, kshlo, kshhi, lshlo, lshhi 770 integer erilen, scrlen 771 double precision schw_ij, tol2e, eri(erilen), iscr(scrlen) 772 integer ibflo, ibfhi, jbflo, jbfhi 773 integer kblo, kbhi, lblo, lbhi 774 double precision sbsb(lblo:lbhi,kblo:kbhi,jbflo:jbfhi,ibflo:ibfhi) 775 double precision sbsbt(kblo:kbhi,lblo:lbhi, 776 $ jbflo:jbfhi,ibflo:ibfhi) 777 integer mxshq 778 integer ishv(mxshq), jshv(mxshq), kshv(mxshq), lshv(mxshq) 779 integer ilab(erilen), jlab(erilen), klab(erilen), llab(erilen) 780c 781 integer bsize, ieri, neri, neritotal, nshq 782 integer ii, jj, kk, ll 783 double precision beffect, q4 784 logical use_q4, more 785c 786 logical intb_2e4c, intb_init4c 787 external intb_2e4c, intb_init4c 788c 789 data use_q4/.false./ 790 data q4/1.d0/ 791c 792c 793 neritotal = 0 794 bsize = (kbhi - kblo + 1)*(lbhi - lblo + 1)* 795 $ (ibfhi - ibflo + 1)*(jbfhi - jbflo + 1) 796 call dfill(bsize,0.d0,sbsb,1) 797 call dfill(bsize,0.d0,sbsbt,1) 798c 799c 800c First time around ( I K | J L ) type 801c 802c Fill out shell arrays with interacting quartet labels 803c 804 nshq = 0 805 do kk=kshlo,kshhi 806 do ll=lshlo,lshhi 807 if ((schwarz_shell(ish,kk)*schwarz_shell(jsh,ll)) 808 $ .ge.tol2e) then 809 if (nshq.eq.mxshq) 810 $ call errquit('moints6x_gblkk: internal error',0, 811 & UNKNOWN_ERR) 812 nshq = nshq + 1 813 kshv(nshq) = kk 814 lshv(nshq) = ll 815 endif 816 enddo 817 enddo 818 call ifill( nshq, ish, ishv, 1 ) 819 call ifill( nshq, jsh, jshv, 1 ) 820c 821c Prepare texas for shell block 822c 823 if (.not. intb_init4c( basis, ishv, kshv, basis, jshv, lshv, 824 $ nshq, q4, use_q4, 825 $ scrlen, iscr, erilen, beffect )) 826 $ call errquit('moints6x_gblkk: intb_init4c failed',0, 827 & INT_ERR) 828c 829c Get some batch of integrals 830c 831 100 more = intb_2e4c( basis, ishv, kshv, basis, jshv, lshv, 832 $ nshq, q4, use_q4, tol2e, .false., 833 $ ilab, klab, jlab, llab, 834 $ eri, erilen, neri, scrlen, iscr ) 835c 836c Unpack labels and assign to integrals 837c 838 if (neri.gt.0) then 839 neritotal = neritotal + neri 840 do ieri=1,neri 841 ii = ilab(ieri) 842 kk = klab(ieri) 843 jj = jlab(ieri) 844 ll = llab(ieri) 845 sbsb(ll,kk,jj,ii) = eri(ieri) 846 enddo 847 endif 848 if (more) goto 100 849c 850c 851c Second time around ( I L | J K ) type 852c 853c Fill out shell arrays with interacting quartet labels 854c 855 nshq = 0 856 do kk=kshlo,kshhi 857 do ll=lshlo,lshhi 858 if ((schwarz_shell(ish,ll)*schwarz_shell(jsh,kk)) 859 $ .ge.tol2e) then 860 if (nshq.eq.mxshq) 861 $ call errquit('moints6x_gblkk: internal error',0, INT_ERR) 862 nshq = nshq + 1 863 kshv(nshq) = kk 864 lshv(nshq) = ll 865 endif 866 enddo 867 enddo 868 call ifill( nshq, ish, ishv, 1 ) 869 call ifill( nshq, jsh, jshv, 1 ) 870c 871c Prepare texas for shell block 872c 873 if (.not. intb_init4c( basis, ishv, lshv, basis, jshv, kshv, 874 $ nshq, q4, use_q4, 875 $ scrlen, iscr, erilen, beffect )) 876 $ call errquit('moints6x_gblkk: intb_init4c failed',0, INT_ERR) 877c 878c Get some batch of integrals 879c 880 200 more = intb_2e4c( basis, ishv, lshv, basis, jshv, kshv, 881 $ nshq, q4, use_q4, tol2e, .false., 882 $ ilab, llab, jlab, klab, 883 $ eri, erilen, neri, scrlen, iscr ) 884c 885c Unpack labels and assign to integrals 886c 887 if (neri.gt.0) then 888 neritotal = neritotal + neri 889 do ieri=1,neri 890 ii = ilab(ieri) 891 ll = llab(ieri) 892 jj = jlab(ieri) 893 kk = klab(ieri) 894 sbsbt(kk,ll,jj,ii) = eri(ieri) 895 enddo 896 endif 897 if (more) goto 200 898c 899c 900c 901 return 902 end 903 904 905 906 907 908 909 subroutine moints6x_eri2blkk( ilen, jlen, klen, llen, 910 $ eri, blk, lblen, kblen ) 911 implicit none 912 integer ilen, jlen, klen, llen, lblen, kblen 913 double precision blk(lblen, kblen, jlen, ilen ) 914 double precision eri(llen, jlen, klen, ilen ) 915 integer k,l,i,j 916 917 do i=1,ilen 918 do j=1,jlen 919 do k=1,klen 920 do l=1,llen 921 blk(l,k,j,i) = eri(l,j,k,i) 922 enddo 923 enddo 924 enddo 925 enddo 926 927 return 928 end 929 930#endif 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946c 947c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 948c 949c NEW VERSIONS WITH SHELL REORDERING 950c 951c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 952c 953c 954c 955c Generate AO integrals in suitable form 956c exchange transformation: 957c Branch according to block or non-block 958c scheme 959c 960c 961 subroutine moints6x_gblkk( basis, ish, jsh, kshlo, kshhi, 962 $ lshlo, lshhi, shmap, rbfmap, 963 $ schw_ij, tol2e, osym, oblk, 964 $ erilen, eri, scrlen, iscr, 965 $ ibflo, ibfhi, jbflo, jbfhi, 966 $ kblo, kbhi, lblo, lbhi, 967 $ sbsb, sbst ) 968 implicit none 969#include "errquit.fh" 970#include "mafdecls.fh" 971#include "bas.fh" 972#include "schwarz.fh" 973 integer basis, ish, jsh, kshlo, kshhi, lshlo, lshhi 974 integer shmap(*), rbfmap(*) 975 integer erilen, scrlen 976 double precision schw_ij, tol2e, eri(erilen), iscr(scrlen) 977 logical osym, oblk 978 integer ibflo, ibfhi, jbflo, jbfhi 979 integer kblo, kbhi, lblo, lbhi 980 double precision sbsb(lblo:lbhi,kblo:kbhi,jbflo:jbfhi,ibflo:ibfhi) 981 double precision sbst(kblo:kbhi,lblo:lbhi,jbflo:jbfhi,ibflo:ibfhi) 982c 983c 984 integer mxshq, imem 985 integer k_tmp, l_tmp 986 integer k_iv, k_jv, k_kv, k_lv 987 integer k_il, k_jl, k_kl, k_ll 988 integer blklen 989 logical status 990 integer blkid, blkidt 991c 992 logical moints_gblk_fromdisk, moints_gblk_todisk 993 external moints_gblk_fromdisk, moints_gblk_todisk 994c 995 blklen = (ibfhi-ibflo+1)*(jbfhi-jbflo+1)* 996 $ (kbhi-kblo+1)*(lbhi-lblo+1) 997c 998c Check to see if we can retrieve from disk 999c For diagonal blocks, we just need to copy since (ik|jl) = (jl|ik) 1000c 1001 blkid = kshlo*1000 + lshlo 1002 blkidt = (kshlo*1000 + lshlo)*1000 1003 if (moints_gblk_fromdisk( blkid, ish, jsh, kshlo, lshlo, 1004 $ ibflo, ibfhi, jbflo, jbfhi, 1005 $ kblo, kbhi, lblo, lbhi, 1006 $ sbsb )) then 1007 if (kshlo.ne.lshlo) then 1008 status = moints_gblk_fromdisk( blkidt, ish, jsh, 1009 $ kshlo, lshlo, 1010 $ ibflo, ibfhi, jbflo, jbfhi, 1011 $ lblo, lbhi, kblo, kbhi, 1012 $ sbst ) 1013 if (.not.(status)) call errquit( 1014 $ 'moints_6x_gblkk: inconsistent AO disk records',0, 1015 & INT_ERR) 1016 else 1017 call dcopy( blklen, sbsb, 1, sbst, 1 ) 1018 endif 1019 return 1020 endif 1021c 1022c Otherwise compute directly 1023c 1024 if (oblk) then 1025 mxshq = (kshhi - kshlo + 1)*(lshhi - lshlo + 1) 1026 imem = 4*mxshq + 4*erilen 1027 if (.not. ma_push_get(MT_INT, imem, 'gblk tmp',l_tmp, k_tmp)) 1028 $ call errquit('moints6x_gblkk: no memory for labels',imem, 1029 & INT_ERR) 1030 k_iv = k_tmp 1031 k_jv = k_tmp + mxshq 1032 k_kv = k_tmp + 2*mxshq 1033 k_lv = k_tmp + 3*mxshq 1034 k_il = k_tmp + 4*mxshq 1035 k_jl = k_tmp + 4*mxshq + erilen 1036 k_kl = k_tmp + 4*mxshq + 2*erilen 1037 k_ll = k_tmp + 4*mxshq + 3*erilen 1038 call moints6x_gblkk_mq( basis, ish, jsh, 1039 $ kshlo, kshhi, lshlo, lshhi, 1040 $ shmap, rbfmap, 1041 $ schw_ij, tol2e, erilen, eri, 1042 $ scrlen, iscr, 1043 $ ibflo, ibfhi, jbflo, jbfhi, 1044 $ kblo, kbhi, lblo, lbhi, 1045 $ sbsb, sbst, mxshq, 1046 $ int_mb(k_iv), int_mb(k_jv), 1047 $ int_mb(k_kv), int_mb(k_lv), 1048 $ int_mb(k_il), int_mb(k_jl), 1049 $ int_mb(k_kl), int_mb(k_ll) ) 1050 if (.not.ma_pop_stack(l_tmp)) 1051 $ call errquit('moints6x_gblkk: failed to pop',0, MA_ERR) 1052 else 1053 call moints6x_gblkk_sq( basis, ish, jsh, 1054 $ kshlo, kshhi, lshlo, lshhi, 1055 $ shmap, rbfmap, 1056 $ schw_ij, tol2e, erilen, eri, 1057 $ scrlen, iscr, 1058 $ ibflo, ibfhi, jbflo, jbfhi, 1059 $ kblo, kbhi, lblo, lbhi, 1060 $ sbsb, sbst ) 1061 endif 1062c 1063c Cache this block to disk if required 1064c Only one copy for diagonal blocks 1065c 1066 status = moints_gblk_todisk( blkid, ish, jsh, kshlo, lshlo, 1067 $ ibflo, ibfhi, jbflo, jbfhi, 1068 $ kblo, kbhi, lblo, lbhi, 1069 $ sbsb ) 1070 if (kshlo.ne.lshlo) then 1071 status = moints_gblk_todisk( blkidt, ish, jsh, kshlo, lshlo, 1072 $ ibflo, ibfhi, jbflo, jbfhi, 1073 $ lblo, lbhi, kblo, kbhi, 1074 $ sbst ) 1075 endif 1076 1077 1078 return 1079 end 1080 1081 1082 1083 1084 1085c 1086c This version uses single shell 1087c quartet AO integrals. 1088c 1089 subroutine moints6x_gblkk_sq( basis, ish, jsh, 1090 $ kshlo, kshhi, lshlo, lshhi, 1091 $ shmap, rbfmap, 1092 $ schw_ij, tol2e, erilen, eri, 1093 $ scrlen, iscr, 1094 $ ibflo, ibfhi, jbflo, jbfhi, 1095 $ kblo, kbhi, lblo, lbhi, 1096 $ sbsb, sbst ) 1097 implicit none 1098#include "bas.fh" 1099#include "schwarz.fh" 1100 integer basis, ish, jsh, kshlo, kshhi, lshlo, lshhi 1101 integer shmap(*), rbfmap(*) 1102 integer erilen, scrlen 1103 double precision schw_ij, tol2e, eri(erilen), iscr(scrlen) 1104 integer ibflo, ibfhi, jbflo, jbfhi 1105 integer kblo, kbhi, lblo, lbhi 1106 double precision sbsb(lblo:lbhi,kblo:kbhi,jbflo:jbfhi,ibflo:ibfhi) 1107 double precision sbst(kblo:kbhi,lblo:lbhi,jbflo:jbfhi,ibflo:ibfhi) 1108c 1109 integer ilen, jlen 1110 integer ksh, lsh, kbflo, kbfhi, lbflo, lbfhi, ltop, kk, ll 1111 integer klen, llen, kblen, lblen, bsize, k1, k2, j, i 1112 logical status 1113c 1114c 1115c 1116 ilen = ibfhi - ibflo + 1 1117 jlen = jbfhi - jbflo + 1 1118 kblen = kbhi - kblo + 1 1119 lblen = lbhi - lblo + 1 1120 bsize = kblen*lblen*ilen*jlen 1121 call dfill(bsize,0.d0,sbsb,1) 1122 call dfill(bsize,0.d0,sbst,1) 1123c 1124 do kk=kshlo,kshhi 1125 ksh = shmap(kk) 1126 status = bas_cn2bfr(basis,ksh,kbflo,kbfhi) 1127 klen = kbfhi - kbflo + 1 1128 ltop = lshhi 1129 if (kshlo.eq.lshlo) ltop = kk 1130 do ll=lshlo,ltop 1131 lsh = shmap(ll) 1132 status = bas_cn2bfr(basis,lsh,lbflo,lbfhi) 1133 llen = lbfhi - lbflo + 1 1134 if (schwarz_shell(ish,ksh)* 1135 $ schwarz_shell(jsh,lsh).ge.tol2e) then 1136 call int_2e4c(basis, ish, ksh, basis, jsh, lsh, 1137 $ scrlen, iscr, erilen, eri ) 1138 call moints6x_eri2blkkx( rbfmap, ilen, jlen, 1139 $ kbflo, kbfhi, lbflo, lbfhi, 1140 $ kblo, kbhi, lblo, lbhi, 1141 $ eri, sbsb ) 1142 endif 1143 if (schwarz_shell(ish,lsh)* 1144 $ schwarz_shell(jsh,ksh).ge.tol2e) then 1145 call int_2e4c(basis, ish, lsh, basis, jsh, ksh, 1146 $ scrlen, iscr, erilen, eri ) 1147 call moints6x_eri2blkkx( rbfmap, ilen, jlen, 1148 $ lbflo, lbfhi, kbflo, kbfhi, 1149 $ lblo, lbhi, kblo, kbhi, 1150 $ eri, sbst ) 1151 endif 1152 enddo 1153 enddo 1154c 1155 if (kshlo.eq.lshlo) then 1156 do i=ibflo,ibfhi 1157 do j=jbflo,jbfhi 1158 do k1=kblo,kbhi 1159 do k2=kblo,k1-1 1160 sbsb(k1,k2,j,i) = sbst(k1,k2,j,i) 1161 sbst(k2,k1,j,i) = sbsb(k2,k1,j,i) 1162 enddo 1163 enddo 1164 enddo 1165 enddo 1166 endif 1167c 1168 return 1169 end 1170 1171 1172 1173 1174 1175 1176c 1177c This uses multiple quartets at a time 1178c AO integrals with labels 1179c 1180 subroutine moints6x_gblkk_mq( basis, ish, jsh, 1181 $ kshlo, kshhi, lshlo, lshhi, 1182 $ shmap, rbfmap, 1183 $ schw_ij, tol2e, erilen, eri, 1184 $ scrlen, iscr, 1185 $ ibflo, ibfhi, jbflo, jbfhi, 1186 $ kblo, kbhi, lblo, lbhi, 1187 $ sbsb, sbst, mxshq, 1188 $ ishv, jshv, kshv, lshv, 1189 $ ilab, jlab, klab, llab ) 1190 implicit none 1191#include "errquit.fh" 1192#include "bas.fh" 1193#include "schwarz.fh" 1194 integer basis, ish, jsh, kshlo, kshhi, lshlo, lshhi 1195 integer shmap(*), rbfmap(*) 1196 integer erilen, scrlen 1197 double precision schw_ij, tol2e, eri(erilen), iscr(scrlen) 1198 integer ibflo, ibfhi, jbflo, jbfhi 1199 integer kblo, kbhi, lblo, lbhi 1200 double precision sbsb(lblo:lbhi,kblo:kbhi,jbflo:jbfhi,ibflo:ibfhi) 1201 double precision sbst(kblo:kbhi,lblo:lbhi,jbflo:jbfhi,ibflo:ibfhi) 1202 integer mxshq 1203 integer ishv(mxshq), jshv(mxshq), kshv(mxshq), lshv(mxshq) 1204 integer ilab(erilen), jlab(erilen), klab(erilen), llab(erilen) 1205c 1206 integer bsize, ieri, neri, neritotal, nshq 1207 integer ii, jj, kk, ll, ksh, lsh 1208 double precision beffect, q4 1209 logical use_q4, more 1210c 1211 logical intb_2e4c, intb_init4c 1212 external intb_2e4c, intb_init4c 1213c 1214 data use_q4/.false./ 1215 data q4/1.d0/ 1216c 1217c 1218 neritotal = 0 1219 bsize = (kbhi - kblo + 1)*(lbhi - lblo + 1)* 1220 $ (ibfhi - ibflo + 1)*(jbfhi - jbflo + 1) 1221 call dfill(bsize,0.d0,sbsb,1) 1222 call dfill(bsize,0.d0,sbst,1) 1223c 1224c 1225c First time around ( I K | J L ) type 1226c 1227c Fill out shell arrays with interacting quartet labels 1228c 1229 nshq = 0 1230 do kk=kshlo,kshhi 1231 ksh = shmap(kk) 1232 do ll=lshlo,lshhi 1233 lsh = shmap(ll) 1234 if ((schwarz_shell(ish,ksh)*schwarz_shell(jsh,lsh)) 1235 $ .ge.tol2e) then 1236 if (nshq.eq.mxshq) 1237 $ call errquit('moints6x_gblkk: internal error',0, INT_ERR) 1238 nshq = nshq + 1 1239 kshv(nshq) = ksh 1240 lshv(nshq) = lsh 1241 endif 1242 enddo 1243 enddo 1244 call ifill( nshq, ish, ishv, 1 ) 1245 call ifill( nshq, jsh, jshv, 1 ) 1246c 1247c Prepare texas for shell block 1248c 1249 if (.not. intb_init4c( basis, ishv, kshv, basis, jshv, lshv, 1250 $ nshq, q4, use_q4, 1251 $ scrlen, iscr, erilen, beffect )) 1252 $ call errquit('moints6x_gblkk: intb_init4c failed',0, INT_ERR) 1253c 1254c Get some batch of integrals 1255c 1256 100 more = intb_2e4c( basis, ishv, kshv, basis, jshv, lshv, 1257 $ nshq, q4, use_q4, tol2e, .false., 1258 $ ilab, klab, jlab, llab, 1259 $ eri, erilen, neri, scrlen, iscr ) 1260c 1261c Unpack labels and assign to integrals 1262c 1263 if (neri.gt.0) then 1264 neritotal = neritotal + neri 1265 do ieri=1,neri 1266 ii = ilab(ieri) 1267 jj = jlab(ieri) 1268 kk = rbfmap(klab(ieri)) 1269 ll = rbfmap(llab(ieri)) 1270 sbsb(ll,kk,jj,ii) = eri(ieri) 1271 enddo 1272 endif 1273 if (more) goto 100 1274c 1275c 1276c Second time around ( I L | J K ) type 1277c 1278c Fill out shell arrays with interacting quartet labels 1279c 1280 nshq = 0 1281 do kk=kshlo,kshhi 1282 ksh = shmap(kk) 1283 do ll=lshlo,lshhi 1284 lsh = shmap(ll) 1285 if ((schwarz_shell(ish,lsh)*schwarz_shell(jsh,ksh)) 1286 $ .ge.tol2e) then 1287 if (nshq.eq.mxshq) 1288 $ call errquit('moints6x_gblkk: internal error',0, INT_ERR) 1289 nshq = nshq + 1 1290 kshv(nshq) = ksh 1291 lshv(nshq) = lsh 1292 endif 1293 enddo 1294 enddo 1295 call ifill( nshq, ish, ishv, 1 ) 1296 call ifill( nshq, jsh, jshv, 1 ) 1297c 1298c Prepare texas for shell block 1299c 1300 if (.not. intb_init4c( basis, ishv, lshv, basis, jshv, kshv, 1301 $ nshq, q4, use_q4, 1302 $ scrlen, iscr, erilen, beffect )) 1303 $ call errquit('moints6x_gblkk: intb_init4c failed',0, INT_ERR) 1304c 1305c Get some batch of integrals 1306c 1307 200 more = intb_2e4c( basis, ishv, lshv, basis, jshv, kshv, 1308 $ nshq, q4, use_q4, tol2e, .false., 1309 $ ilab, llab, jlab, klab, 1310 $ eri, erilen, neri, scrlen, iscr ) 1311c 1312c Unpack labels and assign to integrals 1313c 1314 if (neri.gt.0) then 1315 neritotal = neritotal + neri 1316 do ieri=1,neri 1317 ii = ilab(ieri) 1318 jj = jlab(ieri) 1319 ll = rbfmap(llab(ieri)) 1320 kk = rbfmap(klab(ieri)) 1321 sbst(kk,ll,jj,ii) = eri(ieri) 1322 enddo 1323 endif 1324 if (more) goto 200 1325c 1326c 1327c 1328 return 1329 end 1330 1331 1332 1333 1334 1335 1336 subroutine moints6x_eri2blkkx( rbfmap, ilen, jlen, 1337 $ kbflo, kbfhi, lbflo, lbfhi, 1338 $ kblo, kbhi, lblo, lbhi, 1339 $ eri, blk ) 1340 implicit none 1341 integer rbfmap(*) 1342 integer ilen, jlen 1343 integer kbflo, kbfhi, lbflo, lbfhi 1344 integer kblo, kbhi, lblo, lbhi 1345 double precision blk(lblo:lbhi,kblo:kbhi,jlen,ilen) 1346 double precision eri(lbflo:lbfhi,jlen,kbflo:kbfhi,ilen) 1347 integer kk,ll,k,l,i,j 1348 1349 do i=1,ilen 1350 do j=1,jlen 1351 do k=kbflo,kbfhi 1352 kk = rbfmap(k) 1353 do l=lbflo,lbfhi 1354 ll = rbfmap(l) 1355 blk(ll,kk,j,i) = eri(l,j,k,i) 1356 enddo 1357 enddo 1358 enddo 1359 enddo 1360 1361 return 1362 end 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 subroutine moints6x_pushK( nbf, ilen, jlen, klo, khi, 1378 $ nmo, x, ssni ) 1379 implicit none 1380 integer nbf, ilen, jlen, klo, khi, nmo 1381 double precision x(klo:khi,jlen,ilen,nmo) 1382 double precision ssni(nbf,jlen,ilen,nmo) 1383 integer a,i,j,k 1384 1385 do a=1,nmo 1386 do i=1,ilen 1387 do j=1,jlen 1388 do k=klo,khi 1389 ssni(k,j,i,a) = ssni(k,j,i,a) + x(k,j,i,a) 1390 enddo 1391 enddo 1392 enddo 1393 enddo 1394 return 1395 end 1396 1397 1398 1399 1400 1401 1402 1403 subroutine moints6x_trf1K(nbf, qlo, qhi, molo, mohi, 1404 $ ilen, jlen, klo, khi, llo, lhi, 1405 $ scale, ssbb, ssbbt, c, ssni, hlp ) 1406 implicit none 1407 integer nbf, qlo, qhi 1408 integer molo, mohi, ilen, jlen, klo, khi, llo, lhi 1409 double precision scale 1410 double precision ssbb(llo:lhi,klo:khi,jlen,ilen) 1411 double precision ssbbt(klo:khi,llo:lhi,jlen,ilen) 1412 double precision c(nbf,qlo:qhi) 1413 double precision ssni(nbf,jlen,ilen,molo:mohi) 1414 double precision hlp(*) 1415c 1416c 1417 integer nmo, llen, klen, kjilen, ljilen 1418c 1419c 1420 nmo = mohi - molo + 1 1421 llen = lhi - llo + 1 1422 klen = khi - klo + 1 1423 kjilen = klen*jlen*ilen 1424 ljilen = llen*jlen*ilen 1425 call dgemm( 't', 'n', kjilen, nmo, llen, scale, ssbb, llen, 1426 $ c(llo,molo), nbf, 0.d0, hlp, kjilen ) 1427 call moints6x_pushK( nbf, ilen, jlen, klo, khi, 1428 $ nmo, hlp, ssni ) 1429 call dgemm( 't', 'n', ljilen, nmo, klen, scale, ssbbt, klen, 1430 $ c(klo,molo), nbf, 0.d0, hlp, ljilen ) 1431 call moints6x_pushK( nbf, ilen, jlen, llo, lhi, 1432 $ nmo, hlp, ssni ) 1433 1434 return 1435 end 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 subroutine moints6x_trf2K(nbf, qlo, qhi, ostart, olo, ohi, 1446 $ ilo, ihi, jlo, jhi, ssni, h1, h2, 1447 $ g1, g2, c, g_exch ) 1448 implicit none 1449 integer nbf, qlo, qhi 1450 integer ostart, olo, ohi, ilo, ihi, jlo, jhi 1451 double precision ssni(nbf,jlo:jhi,ilo:ihi,olo:ohi) 1452#ifdef SCATTER_TRANSF 1453 double precision h1(*) 1454 integer iv(100),jv(100) 1455#endif 1456 double precision h1(nbf,ilo:ihi) 1457 double precision h2(jlo:jhi,ilo:ihi) 1458 double precision g1(nbf,jlo:jhi) 1459 double precision g2(ilo:ihi,jlo:jhi) 1460 double precision c(nbf,qlo:qhi) 1461 integer g_exch 1462c 1463 integer ilen, jlen, ijlen 1464 integer nni, nnj, ijlo, ijhi, jilo, jihi 1465 integer aoff, ofroz, aa, bb, a, b, ab, i, j 1466c 1467 ofroz = ostart - 1 1468 aoff = ((olo-ofroz)*(olo-ofroz-1))/2 1469 ilen = ihi - ilo + 1 1470 jlen = jhi - jlo + 1 1471 ijlen = ilen*jlen 1472 nni = ilen*nbf 1473 nnj = jlen*nbf 1474 1475#ifdef BLOCK_TRANSF 1476 ijlo = (ilo-1)*nbf + 1 1477 ijhi = ihi*nbf 1478#endif 1479 do a=olo,ohi 1480 do b=ostart,a 1481 call dgemm('t','n',ijlen,1,nbf,1.d0,ssni(1,jlo,ilo,b), 1482 $ nbf,c(1,a),nbf,0.d0,h2,ijlen) 1483#ifndef NOCOMMS 1484 aa = a - ofroz 1485 bb = b - ofroz 1486 ab = (aa*(aa-1))/2 + bb - aoff 1487#ifdef BLOCK_TRANSF 1488 call dfill(nni,0.d0,h1,1) 1489 do i=ilo,ihi 1490 do j=jlo,jhi 1491 h1(j,i) = h2(j,i) 1492 enddo 1493 enddo 1494 call ga_acc(g_exch,ijlo,ijhi,ab,ab,h1,nni,1.d0) 1495#else 1496#if SCATTER_TRANSF 1497 ij = 0 1498 do i=ilo,ihi 1499 ii = (i-1)*nbf 1500 do j=jlo,jhi 1501 ij = ij + 1 1502 h1(ij) = h2(j,i) 1503 iv(ij) = ii + j 1504 jv(ij) = ab 1505 enddo 1506 enddo 1507 call ga_dscatter(g_exch,h1,iv,jv,ijlen) 1508#else 1509 do i=ilo,ihi 1510 ijlo = (i-1)*nbf + jlo 1511 ijhi = (i-1)*nbf + jhi 1512 call ga_acc(g_exch,ijlo,ijhi,ab,ab,h2(jlo,i),ilen,1.d0) 1513 enddo 1514#endif 1515#endif 1516#endif 1517 enddo 1518 enddo 1519 1520 1521 1522 if (ilo.eq.jlo) goto 120 1523#ifdef BLOCK_TRANSF 1524 jilo = (jlo-1)*nbf + 1 1525 jihi = jhi*nbf 1526#endif 1527 do a=olo,ohi 1528 do b=ostart,a 1529 call dgemm('t','n',ijlen,1,nbf,1.d0,ssni(1,jlo,ilo,a), 1530 $ nbf,c(1,b),nbf,0.d0,h2,ijlen) 1531#ifndef NOCOMMS 1532 aa = a - ofroz 1533 bb = b - ofroz 1534 ab = (aa*(aa-1))/2 + bb - aoff 1535#ifdef BLOCK_TRANSF 1536 call dfill(nnj,0.d0,g1,1) 1537 do j=jlo,jhi 1538 do i=ilo,ihi 1539 g1(i,j) = h2(j,i) 1540 enddo 1541 enddo 1542 call ga_acc(g_exch,jilo,jihi,ab,ab,g1,nnj,1.d0) 1543#else 1544#if SCATTER_TRANSF 1545 ij = 0 1546 do j=jlo,jhi 1547 jj = (j-1)*nbf 1548 do i=ilo,ihi 1549 ij = ij + 1 1550 h1(ij) = h2(j,i) 1551 iv(ij) = jj + i 1552 jv(ij) = ab 1553 enddo 1554 enddo 1555 call ga_dscatter(g_exch,h1,iv,jv,ijlen) 1556#else 1557 do j=jlo,jhi 1558 do i=ilo,ihi 1559 g1(i,j) = h2(j,i) 1560 enddo 1561 enddo 1562 do j=jlo,jhi 1563 jilo = (j-1)*nbf + ilo 1564 jihi = (j-1)*nbf + ihi 1565 call ga_acc(g_exch,jilo,jihi,ab,ab,g1(ilo,j),ilen,1.d0) 1566 enddo 1567#endif 1568#endif 1569#endif 1570 enddo 1571 enddo 1572 1573 120 continue 1574 return 1575 end 1576