diff --git a/src/uccsdt/GNUmakefile b/src/uccsdt/GNUmakefile index 2e86a8fb77..4c9b9dac90 100644 --- a/src/uccsdt/GNUmakefile +++ b/src/uccsdt/GNUmakefile @@ -1,7 +1,7 @@ LIBRARY = libccsd.a -OBJ = uccsdt.o moints_shpair.o +OBJ_OPTIMIZE = uccsdt.o moints_shpair.o soints.o include ../config/makefile.h include ../config/makelib.h diff --git a/src/uccsdt/cbassym.fh b/src/uccsdt/cbassym.fh index 21d31af040..eeaf4275a5 100644 --- a/src/uccsdt/cbassym.fh +++ b/src/uccsdt/cbassym.fh @@ -7,6 +7,8 @@ c integer sh_uq_op(8,nw_max_shells), bf_uq_op(8,nw_max_nbf) integer sh_n_uq, bf_n_uq integer sh_uq(nw_max_shells), bf_uq(nw_max_nbf) + integer sh_nat(2,nw_max_shells), bf_nat(2,nw_max_nbf) + integer sh_uq_bf(4,nw_max_shells) double precision char_tab(8,0:7) integer bf_per_ir(0:7) integer bf_per_ir_cum(0:7) @@ -14,9 +16,10 @@ c character*8 zir(0:7) integer so_uq_to_lab(2,0:7,nw_max_nbf) integer so_lab1(3,nw_max_nbf) - integer so_lab2(0:7,nw_max_nbf) + integer so_lab2(nw_max_nbf) common /cbassym/ sh_op_map, bf_op_map, bf_phase, sh_n_uq_op, $ bf_n_uq_op, sh_uq_op, bf_uq_op, sh_n_uq, bf_n_uq, sh_uq, + $ sh_nat, bf_nat, sh_uq_bf, $ bf_uq, char_tab, bf_per_ir, bf_per_ir_cum, bf_so_ir, $ so_uq_to_lab, so_lab1, so_lab2 common /ccbassym/zir diff --git a/src/uccsdt/cflopsP.fh b/src/uccsdt/cflopsP.fh new file mode 100644 index 0000000000..dc7449f5fa --- /dev/null +++ b/src/uccsdt/cflopsP.fh @@ -0,0 +1,9 @@ +c +c 19th entry for measuring timing overhead +c + double precision flops(19) + double precision times(19) + double precision calls(19) + double precision enter + common /cflopstrp/flops, times, calls, enter + diff --git a/src/uccsdt/cuccsdtP.fh b/src/uccsdt/cuccsdtP.fh index 0a1755fb76..da8a210fc8 100644 --- a/src/uccsdt/cuccsdtP.fh +++ b/src/uccsdt/cuccsdtP.fh @@ -68,18 +68,23 @@ c logical ocheck ! If true check the integral lists c double precision eps(nw_max_nbf,1:2) ! Eigenvalues for each spin +c + integer aoints_file ! EAF file handle for AO integrals + logical aoints_file_read, aoints_file_write + double precision aoints_file_offset c common /cuccsdt/ c c doubles c - $ eps, + $ eps, aoints_file_offset, c c integers c $ basis, geom, g_movecs, k_movecs, l_movecs, $ nsh, nbf, nmo, nact, nir, nc, no, nv, nc_sym, no_sym, nv_sym, - $ c_sym, o_sym, v_sym, map, ocheck, + $ c_sym, o_sym, v_sym, map, ocheck, + $ aoints_file, aoints_file_read, aoints_file_write, c $ asuper, bsuper, csuper, $ asuper_spin, csuper_spin, diff --git a/src/uccsdt/moints_shpair.F b/src/uccsdt/moints_shpair.F index 47ed18121f..90fd8d0e14 100644 --- a/src/uccsdt/moints_shpair.F +++ b/src/uccsdt/moints_shpair.F @@ -603,6 +603,51 @@ c$$$ else if (list .eq. 23) then ! t(m,i,c,b) c$$$ else c$$$ call errquit('uccsdt_triples_get_list: bad list', list) c$$$ end if +c + end + subroutine uccsdt_triples_moints_init() + implicit none +#include "cuccsdtP.fh" +#include "util.fh" +#include "eaf.fh" +c + integer availkw + integer util_scratch_dir_avail_for_me + character*255 filename +c + availkw = util_scratch_dir_avail_for_me()/8 + write(6,*) ' AVAILKW ', availkw , 1.1d0*dble(nbf)**4/(4000d0*nir) + if (availkw .gt. 1.1d0*dble(nbf)**4/(4000d0*nir)) then + call util_file_name('.trpao', .true., .true., filename) + if (eaf_open(filename, eaf_rw, aoints_file) .ne. 0) call + $ errquit('uccsdt_triples: failed to open ao file',0) + aoints_file_read = .false. + aoints_file_write = .true. + else + aoints_file = -1 + aoints_file_read = .false. + aoints_file_write = .false. + end if +c + end + subroutine uccsdt_triples_moints_tidy() + implicit none +#include "cuccsdtP.fh" +#include "util.fh" +#include "eaf.fh" +c + character*255 filename +c + if (aoints_file .ne. -1) then + if (eaf_close(aoints_file) .ne. 0) call errquit + $ ('uccsdt_triples: failed to close aoints file',0) + call util_file_name('.trpao', .true., .true., filename) + if (eaf_delete(filename) .ne. 0) call errquit + $ ('uccsdt_triples: failed to delete aoints file',0) + end if + aoints_file = -1 + aoints_file_read = .false. + aoints_file_write = .false. c end subroutine uccsdt_triples_moints() @@ -670,6 +715,18 @@ c c ublockdim = max. no. of basis functions in a shell block c maxsh = max. no. of shells in a shell block c shmaxbf= max. no of functions in a shell +c +c SYMMETRY IS USED TO REDUCE THE FLOP COUNT BUT IT SHOULD +c ALSO BE USED TO REDUCE THE MEMORY REQUIRED. NEED TO ADJUST +c THE GLOBAL ARRAY SO THAT IT STORES SYMMETRY UNIQUE USH ONLY. +c +c !!! IF BUFFERING ON FILE NEED TO USE SAME TASK DISTRIBUTION +c !!! AFTER FIRST PASS WRITING THE FILE AND ALSO THE SAME +c !!! BLOCKING FACTORS FOR THE AO INTEGRALS ... THIS IS NOT +c !!! YET BEING DONE. +c + if (aoints_file_write .and. nproc.gt.1) call errquit + $ ('moints_triples: see comments above here!!', 0) c if (.not. bas_nbf_cn_max(basis, shmaxbf)) call errquit $ ('uccsdt: bas_nbf_cn_max?', 0) @@ -693,6 +750,8 @@ c end do call errquit('moints_shpair_one_index: insufficient memory',0) 10 continue +c + call moints_aoblock_rewind() c do ublk = 1, nshblock ushlo = shlist(1,ublk) @@ -1279,6 +1338,11 @@ c endif c * call ma_summarize_allocated_blocks() +c + if (aoints_file_write) then + aoints_file_write = .false. + aoints_file_read = .true. + end if c end subroutine moints_triples_final( @@ -1557,7 +1621,6 @@ c double precision g1(nbf,ilo1:ihi1,vlo:vhi,ulo:uhi) double precision g2(nbf,ilo2:ihi2,vlo:vhi,ulo:uhi) c . ! [output] g(w,i,v,u,iset) = (wi|vu) -c c (wi|vu) = sum(l=1,nbf) (wl|vu) Cli c c For given AO shells ush & vsh return in g all integrals with one @@ -1577,7 +1640,7 @@ c Could add the use of symmetry but in initial application the c symmetry is being handled by using unique ush. c integer idim1, idim2, udim, vdim - integer nsh, nbf + integer nsh integer wshlo, wshhi, lshlo, lshhi integer wlo, whi, llo, lhi, junk, l_work, k_work, lenwork integer avail_ma, need, wblk, lblk @@ -1702,25 +1765,31 @@ c . ! [output] g(w,i,v,u,iset) = (wi|vu) double precision ints(llo:lhi,wlo:whi,vlo:vhi,ulo:uhi) c integer u, v, idim1, idim2, wdim, ldim -c -c Form in ints(l,w,v,u) the integrals (lw|vu) with l<=w -c - call moints_aoblock( - $ basis, oskel, - $ lshlo, lshhi, llo, lhi, - $ wshlo, wshhi, wlo, whi, - $ vsh, vsh, vlo, vhi, - $ ush, ush, ulo, uhi, - $ tol, - $ ints) -c -c Now have all of the AO integrals ... transform as necessary. -c Since we are restricting w>=l we have to do l->i and w->i. + integer aodim + logical moints_aoblock_read c idim1 = ihi1 - ilo1 + 1 if (nset .eq. 2) idim2 = ihi2 - ilo2 + 1 wdim = whi - wlo + 1 ldim = lhi - llo + 1 + aodim = ldim*wdim*(vhi-vlo+1)*(uhi-ulo+1) +c +c Form in ints(l,w,v,u) the integrals (lw|vu) with l<=w +c + if (.not. moints_aoblock_read(ush,vsh,aodim,ints)) then + call moints_aoblock( + $ basis, oskel, + $ lshlo, lshhi, llo, lhi, + $ wshlo, wshhi, wlo, whi, + $ vsh, vsh, vlo, vhi, + $ ush, ush, ulo, uhi, + $ tol, + $ ints) + call moints_aoblock_write(ush,vsh,aodim,ints) + endif +c +c Now have all of the AO integrals ... transform as necessary. +c Since we are restricting w>=l we have to do l->i and w->i. c do u = ulo, uhi do v = vlo, vhi @@ -2097,3 +2166,73 @@ c ok = iok .eq. 0 c end + subroutine moints_aoblock_rewind() + implicit none +#include "cuccsdtP.fh" +c + aoints_file_offset = 0.0d0 +c + end + logical function moints_aoblock_read(ush,vsh,aodim,ints) + implicit none +#include "cuccsdtP.fh" +#include "eaf.fh" + integer ush, vsh, aodim + double precision ints(aodim) +c + integer err + character*80 msg +c +c Needs logic adding to +c . check that the correct data is being read +c . compress the data +c . use a GA instead of a file if there is room +c + if (aoints_file_read) then +c +* write(6,*) ' READING ', ush, vsh, aoints_file_offset, aodim +c + err = eaf_read(aoints_file,aoints_file_offset,ints,8*aodim) + if (err .ne. 0) then + call eaf_errmsg(err, msg) + write(6,*) ' !!!! EAF error ', msg + call errquit('moints_aoblock_read: eaf_read?', 0) + end if + aoints_file_offset = aoints_file_offset + 8*aodim + moints_aoblock_read = .true. + else + moints_aoblock_read = .false. + end if +c + end + subroutine moints_aoblock_write(ush,vsh,aodim,ints) + implicit none +#include "cuccsdtP.fh" +#include "eaf.fh" + integer ush, vsh, aodim + double precision ints(aodim) +c + integer err + character*80 msg +c +c Needs logic adding to +c . enable a check that the correct data is being read +c . compress the data +c . use a GA instead of a file if there is room +c + if (aoints_file_write) then +c +* write(6,*) ' WRITING ', ush, vsh, aoints_file_offset, aodim +c + err = eaf_write(aoints_file,aoints_file_offset,ints,8*aodim) + if (err .ne. 0) then + call eaf_errmsg(err, msg) + write(6,*) ' !!!! EAF error ', msg + call errquit('moints_aoblock_write: eaf_read?', 0) + end if + aoints_file_offset = aoints_file_offset + 8*aodim + end if +c + end + + diff --git a/src/uccsdt/soints.F b/src/uccsdt/soints.F new file mode 100644 index 0000000000..e8da62a932 --- /dev/null +++ b/src/uccsdt/soints.F @@ -0,0 +1,782 @@ + subroutine soints_block(basis, nbf, + $ ishuqlo, ishuqhi, ioff, irange, + $ jshuqlo, jshuqhi, joff, jrange, + $ kshuqlo, kshuqhi, koff, krange, + $ lshuqlo, lshuqhi, loff, lrange, + $ odirac, + $ tol, g) + implicit none +#include "mafdecls.fh" +#include "bas.fh" +#include "schwarz.fh" +#include "nwc_const.fh" +#include "cbassym.fh" + integer basis + integer nbf + integer ishuqlo, ishuqhi, irange(2,0:7), ioff(nbf) + integer jshuqlo, jshuqhi, jrange(2,0:7), joff(nbf,0:7) + integer kshuqlo, kshuqhi, krange(2,0:7), koff(nbf,0:7) + integer lshuqlo, lshuqhi, lrange(2,0:7), loff(nbf,0:7) + logical odirac + double precision tol + double precision g(*) +c +c Assumes that int_init(), schwarz_init(), sym_abelian_init(), +c have all been called and that soints_block_info() was used +c to generate the addressing arrays *off() and *range(). +c +c For SOs arising from unique AO shells in the given ranges, +c compute the 2-electron integrals in the SO basis using +c the symmetry-blocked, or second, labelling for the SOs. +c In this scheme SOs are blocked by irrep and within the +c irrep are in the natural order of symmetry unique basis +c functions. See sym_abelian() for more info on the labelling. +c +c if (odirac) then +c . g(i,j,k,l) = +c else +c . g(i,j,k,l) = (ij|kl) +c endif +c +c The integrals non-zero by symmetry are stored densely as if +c g() were declared as g(i,j,symj,k,l,syml,symkl) corresponding +c to the below loops. The array can also be thought as +c being in the order g(ij,kl,symkl) where ij and kl run +c over the pairs with symmetry kl. The offset to the first +c element of a given symkl is symkloff(symkl). +c +c ind = 0 +c do symkl = 0, nir-1 +c . symij = symkl +c . do syml = 0, nir-1 +c . symk = ieor(syml,symkl) +c . do l = SOs of syml from l-shell-block +c . do k = SOs of symk from k-shell-block +c . do symj = 0, nir-1 +c . symi = ieor(symij,symj) +c . do j = SOs of symj from j-shell-block +c . do i = SOs of symi from i-shell-block +c . ind = ind + 1 +c . g(ind) = (ij|kl) or +c +c The block of integrals may be addressed using arrays generated +c by soints_block_info() +c (this will only work for indices non-zero by symmetry and +c with the same loop restrictions) +c +c ijkl = ioff(i) + joff(j,symij) + koff(k,symkl) +c . + loff(l,symkl) +c integral = g(ijkl) +c +c Also useful, +c +c irange(1,symi) = first orbital i of symi in i shell block +c irange(2,symi) = last orbital i of symi in i shell block +c ditto for j, k, l +c lenij(symij) = no. of ij pairs in the block with symmetry ij +c lenkl(symkl) = ditto kl +c symkloff(symkl) = first integral in the symkl block. +c + integer xishuq, xjshuq, xkshuq, xlshuq + integer ishuq, jshuq, kshuq, lshuq, lopuq, kopuq, jopuq + integer ish, jsh, ksh, lsh, jdim, kdim, ldim + integer i, j, k, l, i2, j2, k2, l2, ijkl, ind, iuq +* integer ngot + integer symi, symj, symk, syml, symij, symkl + integer juq, kuq, luq, jnw, knw, lnw, jshnw, kshnw, lshnw, nsh + integer ilo, ihi, jlo, jhi, klo, khi, llo, lhi + integer nopuq + integer l_q4, l_scr, l_eri, l_lab, l_q, l_ao + integer k_q4, k_scr, k_eri, k_lab, k_q, k_ao + integer k_iq, k_jq, k_kq, k_lq + integer k_ilab, k_jlab, k_klab, k_llab + integer leneri, lenscr, lenao + integer maxq, nq, idimuq + double precision smax, skl, gijkl, c(8,8) + integer nops, geom + integer sym_number_ops + external sym_number_ops + integer mapi(nw_max_nbf), mapj(nw_max_nbf), mapk(nw_max_nbf), + $ mapl(nw_max_nbf) +c + if (.not. bas_geom(basis,geom)) call errquit + $ ('soints_block: basis?', 0) + if (.not. bas_numcont(basis,nsh)) call errquit + $ ('soints_block: basis?', 0) +c +c AO integrals are stored with symmetry unique i, and all symmetry +c related j, k, l in natural order. map{ijkl} maps from NWChem AO +c labels into this order offset and modified to address the array. +c + call ifill(nbf,99999999,mapi,1) + call ifill(nbf,99999999,mapj,1) + call ifill(nbf,99999999,mapk,1) + call ifill(nbf,99999999,mapl,1) + jdim = sh_uq_bf(2,jshuqhi) - sh_uq_bf(1,jshuqlo) + 1 + kdim = sh_uq_bf(2,kshuqhi) - sh_uq_bf(1,kshuqlo) + 1 + ldim = sh_uq_bf(2,lshuqhi) - sh_uq_bf(1,lshuqlo) + 1 + idimuq = 0 + do xishuq = ishuqlo, ishuqhi + ish = sh_uq(xishuq) + if (.not. bas_cn2bfr(basis, ish, ilo, ihi)) + $ call errquit('soints_block: cn2bfr?', ish) + do i = ilo,ihi + mapi(i) = idimuq + i - ilo + end do + idimuq = idimuq + ihi - ilo + 1 + end do + lenao = idimuq*jdim*kdim*ldim +c + jlo = sh_uq_bf(1,jshuqlo) + jhi = sh_uq_bf(2,jshuqhi) + do j = jlo, jhi + mapj(bf_nat(2,j)) = (j-jlo)*idimuq + end do + klo = sh_uq_bf(1,kshuqlo) + khi = sh_uq_bf(2,kshuqhi) + do k = klo, khi + mapk(bf_nat(2,k)) = (k-klo)*(idimuq*jdim) + end do + llo = sh_uq_bf(1,lshuqlo) + lhi = sh_uq_bf(2,lshuqhi) + do l = llo, lhi + mapl(bf_nat(2,l)) = (l-llo)*(idimuq*jdim*kdim) + end do +c +c Allocate memory for the integrals. Cannot do them all at once since +c Texas internal memory allocation fails if we give it more than +c 10000 quartets at once. +c + nops = sym_number_ops(geom) + 1 + maxq = 5000 +* maxq = (ishuqhi-ishuqlo+1)*(jshuqhi-jshuqlo+1)* +* $ (kshuqhi-kshuqlo+1)*(lshuqhi-lshuqlo+1)*(nops**3) +c + call intb_mem_2e4c(leneri, lenscr) + leneri = max(leneri,100*(6**4)) ! Room for 100 (dd|dd) + lenscr = lenscr + lenscr/2 +C + if (.not. ma_push_get(mt_dbl, maxq, 'q4', l_q4, k_q4)) + $ call errquit('soints_block: ma q4?',maxq) + if (.not. ma_push_get(mt_int, 4*maxq, 'quartets', l_q, k_q)) + $ call errquit('soints_block: ma q?', 4*maxq) + if (.not. ma_push_get(mt_int, 4*leneri, 'labels', l_lab, k_lab)) + $ call errquit('soints_block: ma lab?', 4*leneri) + if (.not. ma_push_get(mt_dbl, leneri, 'eri', l_eri, k_eri)) + $ call errquit('soints_block: ma eri?', leneri) + if (.not. ma_push_get(mt_dbl, lenscr, 'scr', l_scr, k_scr)) + $ call errquit('soints_block: ma scr?', lenscr) + if (.not. ma_push_get(mt_dbl, lenao, 'ao', l_ao, k_ao)) + $ call errquit('soints_block: ma ao?',lenao) + k_iq = k_q + k_jq = k_q + maxq + k_kq = k_q + maxq*2 + k_lq = k_q + maxq*3 + k_ilab = k_lab + k_jlab = k_lab + leneri + k_klab = k_lab + leneri*2 + k_llab = k_lab + leneri*3 + call dfill(lenao, 0d0, dbl_mb(k_ao), 1) +c +c Loop thru symmetry unique shells, then symmetry related shells jkl +c + nq = 0 + smax = schwarz_max() + do xlshuq = lshuqlo, lshuqhi + lshuq = sh_uq(xlshuq) + do xkshuq = kshuqlo, kshuqhi + kshuq = sh_uq(xkshuq) + do xjshuq = jshuqlo, jshuqhi + jshuq = sh_uq(xjshuq) + do xishuq = ishuqlo, ishuqhi + ishuq = sh_uq(xishuq) + ish = ishuq ! Don't use related ish + do lopuq = 1, sh_n_uq_op(lshuq) + lsh = sh_op_map(sh_uq_op(lopuq,lshuq),lshuq) + do kopuq = 1, sh_n_uq_op(kshuq) + ksh = sh_op_map(sh_uq_op(kopuq,kshuq),kshuq) + skl = schwarz_shell(ksh,lsh) + if (skl*smax.gt.tol) then + do jopuq = 1, sh_n_uq_op(jshuq) + jsh=sh_op_map(sh_uq_op(jopuq,jshuq),jshuq) + if (skl*schwarz_shell(ish,jsh).gt.tol)then + + if ( ish.le.0.or.ish.gt.nsh.or. + $ jsh.le.0.or.jsh.gt.nsh.or. + $ ksh.le.0.or.ksh.gt.nsh.or. + $ lsh.le.0.or.lsh.gt.nsh) then + write(6,*) ish,jsh,ksh,lsh + call errquit('bad sh',0) + end if + int_mb(k_iq + nq) = ish + int_mb(k_jq + nq) = jsh + int_mb(k_kq + nq) = ksh + int_mb(k_lq + nq) = lsh + dbl_mb(k_q4 + nq) = 1.0d0 + nq = nq + 1 +c + if (nq .eq. maxq) then + call soints_aoints(basis, nq, + $ int_mb(k_iq), int_mb(k_jq), + $ int_mb(k_kq), int_mb(k_lq), + $ dbl_mb(k_q4), tol, + $ int_mb(k_ilab), int_mb(k_jlab), + $ int_mb(k_klab), int_mb(k_llab), + $ dbl_mb(k_eri), leneri, + $ dbl_mb(k_scr), lenscr, + $ mapi, mapj, mapk, mapl, + $ dbl_mb(k_ao), lenao) + nq = 0 + end if + end if + end do + end if + end do + end do + end do + end do + end do + end do + call soints_aoints(basis, nq, + $ int_mb(k_iq), int_mb(k_jq), + $ int_mb(k_kq), int_mb(k_lq), + $ dbl_mb(k_q4), tol, + $ int_mb(k_ilab), int_mb(k_jlab), + $ int_mb(k_klab), int_mb(k_llab), + $ dbl_mb(k_eri), leneri, + $ dbl_mb(k_scr), lenscr, + $ mapi, mapj, mapk, mapl, + $ dbl_mb(k_ao), lenao) +c +c$$$ 789 format(i6,4x,4i5,f12.6) +c$$$ write(6,*) ' The AO integrals' +c$$$ ind = 0 +c$$$ ngot = 0 +c$$$ do l = sh_uq_bf(1,lshuqlo), sh_uq_bf(2,lshuqhi) +c$$$ do k = sh_uq_bf(1,kshuqlo), sh_uq_bf(2,kshuqhi) +c$$$ do j = sh_uq_bf(1,jshuqlo), sh_uq_bf(2,jshuqhi) +c$$$ do iuq = sh_uq_bf(3,ishuqlo), sh_uq_bf(4,ishuqhi) +c$$$ gijkl = dbl_mb(k_ao+ind) ! The integral +c$$$ ind = ind + 1 +c$$$ if (abs(gijkl).gt.1e-6) then +c$$$ write(6,789) ind, iuq,j,k,l,gijkl +c$$$ ngot = ngot + 1 +c$$$ end if +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ write(6,*) ngot +c +c Now have all of the AO integrals. +c +c For each of lsh, ksh, and jsh loop thru the generators +c (symmetry unique/lexically highest functions) in the block of +c unique shells and transform them to the SO basis. +c + do lshuq = lshuqlo, lshuqhi ! Loop thru unique shells + lshnw = sh_uq(lshuq) ! NWChem order shell + nopuq = sh_n_uq_op(lshnw) ! No. of related functions + do luq = sh_uq_bf(3,lshuq), sh_uq_bf(4,lshuq) + lnw = bf_uq(luq) + call soints_uq_to_so_matrix(luq, c) +c$$$ write(6,*) ' luq ', luq, ' nopuq ', nopuq, ' lnw ', lnw, +c$$$ $ ' mapl ', mapl(lnw) +c$$$ call output(c,1,nopuq,1,nopuq,8,8,1) + call soints_inplace_mxm(idimuq*jdim*kdim, nopuq, 1, + $ dbl_mb(k_ao+mapl(lnw)), ldim, c) + end do + end do +c$$$ write(6,*) ' After L' +c$$$ ind = 0 +c$$$ ngot = 0 +c$$$ do l = sh_uq_bf(1,lshuqlo), sh_uq_bf(2,lshuqhi) +c$$$ do k = sh_uq_bf(1,kshuqlo), sh_uq_bf(2,kshuqhi) +c$$$ do j = sh_uq_bf(1,jshuqlo), sh_uq_bf(2,jshuqhi) +c$$$ do iuq = sh_uq_bf(3,ishuqlo), sh_uq_bf(4,ishuqhi) +c$$$ gijkl = dbl_mb(k_ao+ind) ! The integral +c$$$ ind = ind + 1 +c$$$ if (abs(gijkl).gt.1e-6) then +c$$$ write(6,789) ind, iuq,j,k,l,gijkl +c$$$ ngot = ngot + 1 +c$$$ end if +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ write(6,*) ngot + do kshuq = kshuqlo, kshuqhi ! Loop thru unique shells + kshnw = sh_uq(kshuq) ! NWChem order shell + nopuq = sh_n_uq_op(kshnw) ! No. of related functions + do kuq = sh_uq_bf(3,kshuq), sh_uq_bf(4,kshuq) + knw = bf_uq(kuq) + call soints_uq_to_so_matrix(kuq, c) + call soints_inplace_mxm(idimuq*jdim, nopuq, ldim, + $ dbl_mb(k_ao+mapk(knw)), kdim, c) + end do + end do +c$$$ write(6,*) ' After K' +c$$$ ind = 0 +c$$$ ngot = 0 +c$$$ do l = sh_uq_bf(1,lshuqlo), sh_uq_bf(2,lshuqhi) +c$$$ do k = sh_uq_bf(1,kshuqlo), sh_uq_bf(2,kshuqhi) +c$$$ do j = sh_uq_bf(1,jshuqlo), sh_uq_bf(2,jshuqhi) +c$$$ do iuq = sh_uq_bf(3,ishuqlo), sh_uq_bf(4,ishuqhi) +c$$$ gijkl = dbl_mb(k_ao+ind) ! The integral +c$$$ ind = ind + 1 +c$$$ if (abs(gijkl).gt.1e-6) then +c$$$ write(6,789) ind, iuq,j,k,l,gijkl +c$$$ ngot = ngot + 1 +c$$$ end if +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ write(6,*) ngot + do jshuq = jshuqlo, jshuqhi ! Loop thru unique shells + jshnw = sh_uq(jshuq) ! NWChem order shell + nopuq = sh_n_uq_op(jshnw) ! No. of related functions + do juq = sh_uq_bf(3,jshuq), sh_uq_bf(4,jshuq) + jnw = bf_uq(juq) + call soints_uq_to_so_matrix(juq, c) + call soints_inplace_mxm(idimuq, nopuq, kdim*ldim, + $ dbl_mb(k_ao+mapj(jnw)), jdim, c) + end do + end do +c$$$ write(6,*) ' After K' +c$$$ ind = 0 +c$$$ ngot = 0 +c$$$ do l = sh_uq_bf(1,lshuqlo), sh_uq_bf(2,lshuqhi) +c$$$ do k = sh_uq_bf(1,kshuqlo), sh_uq_bf(2,kshuqhi) +c$$$ do j = sh_uq_bf(1,jshuqlo), sh_uq_bf(2,jshuqhi) +c$$$ do iuq = sh_uq_bf(3,ishuqlo), sh_uq_bf(4,ishuqhi) +c$$$ gijkl = dbl_mb(k_ao+ind) ! The integral +c$$$ ind = ind + 1 +c$$$ if (abs(gijkl).gt.1e-6) then +c$$$ write(6,789) ind, iuq,j,k,l,gijkl +c$$$ ngot = ngot + 1 +c$$$ end if +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ write(6,*) ngot +c +c The transformation is complete. Only the unique ish were +c computed and they were not transformed, but the transformation +c of the other three indices act as a projector for their +c symmetry. +c +c A range of unique shells will give rise to the same ranges +c for both basis functions and SOs in natural order. +c + ind = 0 + do l = sh_uq_bf(1,lshuqlo), sh_uq_bf(2,lshuqhi) + syml = so_lab1(1,l) + l2 = so_lab1(2,l) + do k = sh_uq_bf(1,kshuqlo), sh_uq_bf(2,kshuqhi) + symk = so_lab1(1,k) + k2 = so_lab1(2,k) + symkl= ieor(symk,syml) + symij= symkl + do j = sh_uq_bf(1,jshuqlo), sh_uq_bf(2,jshuqhi) + symj = so_lab1(1,j) + j2 = so_lab1(2,j) + symi = ieor(symij,symj) + do iuq = sh_uq_bf(3,ishuqlo), sh_uq_bf(4,ishuqhi) + i2 = so_uq_to_lab(2,symi,iuq) +c + gijkl = dbl_mb(k_ao+ind) ! The integral + ind = ind + 1 + if (i2 .lt. 0) then + if (abs(gijkl) .gt. 1d-6) then + write(6,*) i2,j2,k2,l2,gijkl + call errquit('should be zero?',ijkl) + end if + else + ijkl = ioff(i2) + joff(j2,symij) + + $ koff(k2,symkl) + loff(l2,symkl) + g(ijkl) = gijkl + end if + end do + end do + end do + end do +c + if (.not. ma_chop_stack(l_q4)) call errquit + $ ('soints_aoblock: ma corrupted?',0) +c + end + subroutine soints_block_info(basis, nbf, + $ ishuqlo, ishuqhi, ioff, irange, + $ jshuqlo, jshuqhi, joff, jrange, + $ kshuqlo, kshuqhi, koff, krange, + $ lshuqlo, lshuqhi, loff, lrange, + $ lenij, lenkl, symkloff, leng) + implicit none + integer basis + integer nbf + integer ishuqlo, ishuqhi, irange(2,0:7), ioff(nbf) + integer jshuqlo, jshuqhi, jrange(2,0:7), joff(nbf,0:7) + integer kshuqlo, kshuqhi, krange(2,0:7), koff(nbf,0:7) + integer lshuqlo, lshuqhi, lrange(2,0:7), loff(nbf,0:7) + integer lenij(0:7), lenkl(0:7), symkloff(0:7), leng +c +c See the comments in soints_block() for detailed info. +c +c lenij(symij) is the number of ij pairs for a given symmetry +c lenkl(symkl) ................ kl .......................... +c symkloff(symkl) = index of the first element in the symkl block. +c leng returns the no. of integrals non-zero by symmetry. +c + integer symi, symj, symk, syml, symij, symkl + integer i, j, k, l, ijkl, count + integer leni, lenk +c +c Generate the ranges of SOs in each block of shells +c + call soints_shblock_info(basis, ishuqlo, ishuqhi, irange) + call soints_shblock_info(basis, jshuqlo, jshuqhi, jrange) + call soints_shblock_info(basis, kshuqlo, kshuqhi, krange) + call soints_shblock_info(basis, lshuqlo, lshuqhi, lrange) +c +c Generate the addressing arrays for the integrals. Only entries +c within the given blocks will be valid. +c + call ifill( nbf, 99999999, ioff, 1) + call ifill(8*nbf, 99999999, joff, 1) + call ifill(8*nbf, 99999999, koff, 1) + call ifill(8*nbf, 99999999, loff, 1) +c +c ioff(i) is just the index within its symmetry block +c + do symi = 0, 7 + do i = irange(1,symi), irange(2,symi) + ioff(i) = 1 + i - irange(1,symi) + end do + end do +c +c For given (symij, symj) the matrix (i,j) is dense +c ... joff(j,symij) is the offset to the top of the +c j'th column within the given symmetry, plus the size of +c all the preceeding i symmetries. +c + do symij = 0, 7 + lenij(symij) = 0 + do symj = 0, 7 + symi = ieor(symij,symj) + do j = jrange(1,symj),jrange(2,symj) + joff(j,symij) = lenij(symij) + leni = max(0, irange(2,symi)-irange(1,symi)+1) + lenij(symij) = lenij(symij) + leni + end do + end do + end do +c +c koff is offset within its symmetry block * lenij +c + do symkl = 0,7 + symij = symkl + do syml = 0, 7 + symk = ieor(symkl,syml) + do k = krange(1,symk), krange(2,symk) + koff(k,symkl) = (k - krange(1,symk))*lenij(symij) + end do + end do + end do +c +c loff as for joff but * lenij and accumulate all previous kl symmetries +c + count = 0 + do symkl = 0, 7 + symij = symkl + lenkl(symkl) = 0 + symkloff(symkl) = count + 1 + do syml = 0, 7 + symk = ieor(symkl,syml) + do l = lrange(1,syml),lrange(2,syml) + loff(l,symkl) = count + lenk = max(0, krange(2,symk)-krange(1,symk)+1) + lenkl(symkl) = lenkl(symkl) + lenk + count = count + lenk*lenij(symij) + end do + end do + end do +c +c Verify the indexing ... comment this out when we are sure +c it is working since it is O(N^4). +c + count = 0 + do symkl = 0, 7 + symij = symkl + do syml = 0, 7 + symk = ieor(symkl,syml) + do l = lrange(1,syml),lrange(2,syml) + do k = krange(1,symk),krange(2,symk) + do symj = 0, 7 + symi = ieor(symij,symj) + do j = jrange(1,symj),jrange(2,symj) + do i = irange(1,symi),irange(2,symi) + count = count + 1 + ijkl = ioff(i) + joff(j,symij) + + $ koff(k,symkl) + loff(l,symkl) + if (ijkl .ne. count) then + write(6,*) symi, symj, symk, syml, + $ symij, symkl + write(6,*) i,j,k,l,ioff(i),joff(j,symij), + $ koff(k,symkl),loff(l,symkl), + $ ijkl, count + call errquit(' soints_block_info: bad',0) + endif + end do + end do + end do + end do + end do + end do + end do +c + leng = count +c + end + subroutine soints_shblock_info(basis, ishuqlo, ishuqhi, irange) + implicit none +#include "bas.fh" +#include "nwc_const.fh" +#include "cbassym.fh" + integer basis, ishuqlo, ishuqhi ! [input] + integer irange(2,0:7) ! [output] +c +c Given a range of symmetry unique AO shells determine +c the range of SO functions within each irrep (in the +c second labelling scheme). +c +c Loop thru all of the unique basis functions. Find the +c first/last function in natural in the first/last shell. +c If there are no functions in an irrep, the range is set +c to [0,-1]. +c + integer ir, lab1, lab1lo, lab1hi +c + do ir = 0, 7 + irange(1,ir) = 0 + irange(2,ir) = -1 + end do +c +c Determine the range of basis functions in natural order +c ... the range of SOs in natual order (lab1) will be the same. +c + lab1lo = sh_uq_bf(1,ishuqlo) + lab1hi = sh_uq_bf(2,ishuqhi) +c +c Loop forward thru the SOs to determine the first in each irrep +c then backward to get the last, mapping into the symmetry +c blocked order (lab2) at the same time. +c + do ir = 0,7 + do lab1 = lab1lo,lab1hi + if (so_lab1(1,lab1).eq.ir .and. irange(1,ir).eq.0) then + irange(1,ir) = so_lab1(2,lab1) + goto 10 + endif + end do + 10 continue + end do + do ir = 0,7 + do lab1 = lab1hi,lab1lo,-1 + if (so_lab1(1,lab1).eq.ir .and. irange(2,ir).eq.-1) then + irange(2,ir) = so_lab1(2,lab1) + goto 20 + end if + end do + 20 continue + end do +c +c$$$ write(6,*) ' SO UQSHRANGE: ', ishuqlo, ishuqhi +c$$$ write(6,1) (irange(1,ir),irange(2,ir),ir=0,7) +c$$$ 1 format(8(2i4,4x)) +c + end + subroutine soints_uq_to_so_matrix(luq, c) + implicit none +#include "nwc_const.fh" +#include "cbassym.fh" + integer luq ! [input] + double precision c(8,8) ! [output] +c +c For symmetry unique basis function number luq return +c in C the matrix that transforms the symmetry related +c AO functions into SO basis +c +c |iruq> = sum(k=1..nopuq) C(iruq,luq) |k> +c +c where iruq runs over the irreducible representations arising +c from this generator (the actual irrep is ir=bf_so_ir(iruq,luq)) +c +c and k runs over the functions related by symmetry to the +c unique basis function (their NWChem numbering is given +c as bf_op_map(bf_uq_op(opuq),bf_uq(luq)), opuq=1,nopuq. +c + integer nopuq, opuq, op, iruq, ir, lnw + double precision scale, phase +c + call dfill(8*8, 1d300, c, 1) ! Make unitialized numbers bad news +c + lnw = bf_uq(luq) ! NWChem label for unique function + nopuq = bf_n_uq_op(lnw) ! No. of related functions + scale = 1d0/dble(nopuq) + do opuq = 1,nopuq + op = bf_uq_op(opuq,lnw) + phase = bf_phase(op,lnw)*scale + do iruq = 1, nopuq + ir = bf_so_ir(iruq,luq) + c(iruq,opuq) = char_tab(op,ir)*phase + end do + end do +c + end + subroutine soints_inplace_mxm(n, nopuq, m, a, dim2, c) + implicit none + integer n, nopuq, m, dim2 + double precision a(n,dim2,m), c(8,8) +c +c c() is a small matrix that transforms a set of symmetry +c related functions from the AO to SO bases and a() is +c a set of integrals that are to be transformed. +c +c a(i,iruq,j) <= sum(opuq=1,nopuq) a(i,opuq,j)*c(iruq,opuq) +c +c This routine does the transformation inplace via a +c small buffer. +c + integer maxbuf + parameter (maxbuf = 64) + double precision buf(maxbuf,8) + integer ilo, ihi, idim, i, j, ii, opuq, iruq +c + do j = 1, m + do ilo = 1, n, maxbuf + ihi = min(n,ilo+maxbuf-1) + idim = ihi-ilo+1 + do iruq = 1, nopuq + do i = 1, idim + buf(i,iruq) = 0d0 + end do + end do + do opuq = 1, nopuq + do iruq = 1, nopuq + do i = ilo,ihi + ii = i-ilo+1 + buf(ii,iruq) = buf(ii,iruq) + + $ a(i,opuq,j)*c(iruq,opuq) + end do + end do + end do + do iruq = 1, nopuq + do i = ilo,ihi + ii = i-ilo+1 + a(i,iruq,j) = buf(ii,iruq) + end do + end do + end do + end do +c + end + subroutine soints_aoints(basis, nq, iq, jq, kq, lq, + $ q4, tol, ilab, jlab, klab, llab, eri, leneri, + $ scr, lenscr, mapi, mapj, mapk, mapl, ao, lenao) + implicit none + integer basis + integer nq + integer iq(*), jq(*), kq(*), lq(*) + integer ilab(*), jlab(*), klab(*), llab(*) + integer mapi(*), mapj(*), mapk(*), mapl(*) + integer leneri, lenscr, lenao + double precision q4(*), tol, eri(*), scr(*), ao(*) +c + integer neri, ieri, ijkl + logical more + double precision block_eff + logical intb_init4c, intb_2e4c + external intb_init4c, intb_2e4c +c + if (nq .le. 0) return +c + if (.not. intb_init4c( + $ basis, iq, jq, basis, kq, lq, + $ nq, q4, .true., + $ lenscr, scr, leneri, block_eff)) + $ call errquit('intb_init?',nq) +c + 10 more = intb_2e4c( + $ basis, iq, jq, basis, kq, lq, + $ nq, q4, .true., tol, .false., + $ ilab, jlab, klab, llab, + $ eri, leneri, neri, lenscr, scr) +c +c Map the labels from NWChem order into symmetry natural +c order and then into the SO blocked ordering and then put the +c AO integrals temporaily into where the SO integrals will end up. +c We are exploiting here that there are as many symmetry related +c functions as there are symmetry adapted orbitals. +c + do ieri = 1, neri + ijkl = 1 + + $ mapi(ilab(ieri)) + + $ mapj(jlab(ieri)) + + $ mapk(klab(ieri)) + + $ mapl(llab(ieri)) + if (ijkl.lt.1 .or. ijkl.gt.lenao) call errquit('ijkl',ijkl) + ao(ijkl) = eri(ieri) + end do +c + if (more) goto 10 +c + end + subroutine soints_sovecs(nbf,nmo,c) + implicit none +#include "nwc_const.fh" +#include "cbassym.fh" + integer nbf, nmo + double precision c(nbf,nmo) +c +c On input C() contains a set of molecular orbital coefficients +c in the atomic orbital basis ... transform the AO index into the +c symmetry-adapted orbital basis. The resulting SOs are blocked +c by symmetry (second labelling scheme). +c +c Assumes cbassym has been initialized. +c + integer i, j, juq, nopuq, iruq, opuq + double precision t(nw_max_nbf), s(8,8), sum +c + write(6,*) ' INPUT' + call output(c,1,nbf,1,nmo,nbf,nmo,1) +c +c First reorder the AO functions into symmetry natural order, +c and transform into SOs in natural order, then copy back +c in blocked order. +c + do i = 1, nmo + do j = 1, nbf + t(bf_nat(1,j)) = c(j,i) + end do + j = 0 + do juq = 1, bf_n_uq + nopuq = bf_n_uq_op(bf_uq(juq)) + call soints_uq_to_so_matrix(juq, s) + do iruq = 1, nopuq + sum = 0d0 + do opuq = 1, nopuq + sum = sum + s(iruq,opuq)*t(j+opuq) + end do + c(so_lab1(2,j+iruq),i) = sum + end do + j = j + nopuq + end do + end do +c + write(6,*) ' OUTPUT' + call output(c,1,nbf,1,nmo,nbf,nmo,1) +c + end diff --git a/src/uccsdt/uccsdt.F b/src/uccsdt/uccsdt.F index cc4dbd3bd6..1feb218deb 100644 --- a/src/uccsdt/uccsdt.F +++ b/src/uccsdt/uccsdt.F @@ -60,6 +60,24 @@ c $ l_irs_b, k_irs_b, k_irs(2), l_work, k_work integer i, ispin, lo, hi, ir character*4 uline(8) +c$$$c +c$$$c For debug +c$$$c +c$$$ integer +c$$$ $ ishuqlo, ishuqhi, ioff(nw_max_nbf), +c$$$ $ irange(2,0:7), +c$$$ $ jshuqlo, jshuqhi, joff(nw_max_nbf,0:7), +c$$$ $ jrange(2,0:7), +c$$$ $ kshuqlo, kshuqhi, koff(nw_max_nbf,0:7), +c$$$ $ krange(2,0:7), +c$$$ $ lshuqlo, lshuqhi, loff(nw_max_nbf,0:7), +c$$$ $ lrange(2,0:7), +c$$$ $ lenij(0:7), lenkl(0:7), leng, symkloff(0:7) +c$$$ double precision g(100000) +c$$$c +c$$$c end debug +c$$$c + c logical int_normalize external int_normalize @@ -152,6 +170,7 @@ c $ sh_uq_op, bf_uq_op, $ sh_n_uq, bf_n_uq, $ sh_uq, bf_uq, + $ sh_nat, bf_nat, sh_uq_bf, $ char_tab, $ bf_per_ir, bf_per_ir_cum, $ bf_so_ir, @@ -421,6 +440,46 @@ c call output(dbl_mb(k_movecs(2)),1,nbf,1,nmo,nbf,nmo,1) end if end if +c +c$$$c$$$ call soints_sovecs(nbf,nmo,dbl_mb(k_movecs(1))) +c$$$c$$$ stop +c$$$c +c$$$c Test the SO integral code +c$$$c +c$$$ do lshuqlo = 1, sh_n_uq, 4 +c$$$ lshuqhi = min(sh_n_uq,lshuqlo+4-1) +c$$$ do kshuqlo = 1, lshuqlo, 4 +c$$$ kshuqhi = min(sh_n_uq,kshuqlo+4-1) +c$$$ do jshuqlo = 1, sh_n_uq, 4 +c$$$ jshuqhi = min(sh_n_uq,jshuqlo+4-1) +c$$$ do ishuqlo = 1, jshuqlo, 4 +c$$$ ishuqhi = min(sh_n_uq,ishuqlo+4-1) +c$$$ +c$$$ write(6,*) ishuqlo, jshuqlo, kshuqlo, lshuqlo +c$$$ +c$$$ call soints_block_info(basis, nbf, +c$$$ $ ishuqlo, ishuqhi, ioff, irange, +c$$$ $ jshuqlo, jshuqhi, joff, jrange, +c$$$ $ kshuqlo, kshuqhi, koff, krange, +c$$$ $ lshuqlo, lshuqhi, loff, lrange, +c$$$ $ lenij, lenkl, symkloff, leng) +c$$$ if (leng .gt. 100000) call errquit('leng',0) +c$$$c +c$$$ call soints_block(basis, nbf, +c$$$ $ ishuqlo, ishuqhi, ioff, irange, +c$$$ $ jshuqlo, jshuqhi, joff, jrange, +c$$$ $ kshuqlo, kshuqhi, koff, krange, +c$$$ $ lshuqlo, lshuqhi, loff, lrange, +c$$$ $ .false., +c$$$ $ 1d-12, g) +c$$$ +c$$$c +c$$$ end do +c$$$ end do +c$$$ end do +c$$$ end do +c$$$c +c$$$ write(6,*) ' SOINTS seem OK' c end subroutine uccsdt_tidy(rtdb) @@ -485,6 +544,7 @@ c $ int_bcik(5000), c $ int_acik(5000) c c double precision int_buf(5000), t_buf(5000), w(100000) + double precision tmpr(20*20*20*20) integer btop, ctop, ktop, vlo, vhi, olo, ohi, $ vlo_ck, vhi_ck, olo_ck, ohi_ck integer lentbuf, l_tbuf, k_tbuf @@ -624,8 +684,8 @@ c $ sym3(syma,symb,symc)) if (test_sym .eq. 0) then if (next .eq. count) then - write(6,*) ga_nodeid(), ' doing ', - $ iblk, jblk, kblk, ablk, bblk, cblk +* write(6,*) ga_nodeid(), ' doing ', +* $ iblk, jblk, kblk, ablk, bblk, cblk call util_flush(6) if (oprint) write(6,*) $ 'TSB:kblk, klo, khi, symk ', @@ -773,7 +833,7 @@ c $ clo, chi, elo, ehi, $ ilo, ihi, jlo, jhi, $ klo, khi, - $ spin, spinck) + $ spin, spinck, tmpr) end if c spine = spinck @@ -1281,7 +1341,7 @@ c $ int_acik(a,c,i,k) c if (spin .eq. spinck) then - v = v + + v = v $ + t1(c,i)* $ int_abjk(a,b,j,k) $ - t1(c,j)* @@ -1363,7 +1423,7 @@ c logical oprint, oprint_more integer maxsh, wblockdim, ublockdim integer a1, a2, c1, c2, numasuper, numcsuper, - $ a_block_max, c_block_max, e_block_max + $ a_block_max, c_block_max integer numasuper_best, numcsuper_best, npass_best double precision avail_ma, avail_ga, ga_max, ma_max, $ ma_moints, ga_moints, ma_list, ga_list, ma_wv, @@ -1658,8 +1718,8 @@ c c For debug set them to something sensible and then c recompute their actual values. c - oblock_max = 15 - vblock_max = 15 + oblock_max = 20 + vblock_max = 20 c c Sub-block the 'a' super-block ... sub-blocks will be restricted c to orbitals of the same symmetry whereas super-blocks can contain @@ -2355,8 +2415,6 @@ c integer nasuper, ncsuper ! No. of super blocks for current spin case integer a_super_block(2,max_nsuper_block) ! Super block info integer c_super_block(2,max_nsuper_block) - integer l_t2aa, l_t2ab, l_t2bb, l_t1a, l_t1b - integer k_t2aa, k_t2ab, k_t2bb, k_t1a, k_t1b integer l_ttemp, k_ttemp integer noa, nob, nva, nvb c @@ -2366,18 +2424,18 @@ c c For debug purposes generate the full list of Dirac integrals c in a common block c - if (.not. ma_push_get(mt_dbl, nbf*nbf, 'most', l_most, k_most)) - $ call errquit('memory for most?',nbf*nbf) - if (.not. ma_push_get(mt_dbl, nbf**4, 'moint', l_moint, k_moint)) - $ call errquit('memory for moint?',nbf**4) - call util_transpose(dbl_mb(k_movecs(1)),nbf, - $ dbl_mb(k_most),nbf,nbf,nbf) - call jan_full_transform_noinit( - $ rtdb, basis, - $ nbf, nbf, nbf, nbf, - $ nbf, nbf, nbf, nbf, - $ dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most), - $ dbl_mb(k_moint), 'Dirac') +c$$$ if (.not. ma_push_get(mt_dbl, nbf*nbf, 'most', l_most, k_most)) +c$$$ $ call errquit('memory for most?',nbf*nbf) +c$$$ if (.not. ma_push_get(mt_dbl, nbf**4, 'moint', l_moint, k_moint)) +c$$$ $ call errquit('memory for moint?',nbf**4) +c$$$ call util_transpose(dbl_mb(k_movecs(1)),nbf, +c$$$ $ dbl_mb(k_most),nbf,nbf,nbf) +c$$$ call jan_full_transform_noinit( +c$$$ $ rtdb, basis, +c$$$ $ nbf, nbf, nbf, nbf, +c$$$ $ nbf, nbf, nbf, nbf, +c$$$ $ dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most), +c$$$ $ dbl_mb(k_moint), 'Dirac') c c For debug purposes generate the full list of t1 and t2 amplitudes c @@ -2432,6 +2490,10 @@ c debug: temp zero t1a c c call dfill (nva*noa, 0.0d0, dbl_mb(k_t1a), 1) c +c + call uccsdt_triples_moints_init() + call uccsdt_flops_init() +c c Loop over the four spin cases c total = 0.0d0 @@ -2495,19 +2557,22 @@ c total = total + energy(spina, spinc) end do end do -c +c if (ga_nodeid() .eq. 0) then write(6,77) total 77 format(//' The total (T) energy is ', f20.12,//) end if +c + call uccsdt_triples_moints_tidy() + call uccsdt_flops_print() c if (.not. ma_pop_stack(l_t2ab)) call errquit('t2ab?',0) if (.not. ma_pop_stack(l_t2bb)) call errquit('t2bb?',0) if (.not. ma_pop_stack(l_t2aa)) call errquit('t2aa?',0) if (.not. ma_pop_stack(l_t1b)) call errquit('t1b?',0) if (.not. ma_pop_stack(l_t1a)) call errquit('t1a?',0) - if (.not. ma_pop_stack(l_moint)) call errquit('moint?',0) - if (.not. ma_pop_stack(l_most)) call errquit('most?',0) +* if (.not. ma_pop_stack(l_moint)) call errquit('moint?',0) +* if (.not. ma_pop_stack(l_most)) call errquit('most?',0) c end subroutine sym_map_amps_t1( @@ -2551,22 +2616,24 @@ c $ call errquit('x',1) double precision t2(nbf, nbf, nbf, nbf) double precision t_buf(nbf, nbf, nbf, nbf) integer v1, v2, o1, o2 + integer vv2, oo1, oo2 integer spin1, spin2 call dfill(nbf**4, 0.0d0, t_buf, 1) call dcopy(nbf**4, t2, 1, t_buf, 1) call dfill(nbf**4, 0.0d0, t2, 1) - do v1 = 1, nbf - do v2 = 1, nbf - do o1 = 1, nbf - do o2 = 1, nbf + do o2 = 1, nbf + oo2 = map(o2,spin2) + do o1 = 1, nbf + oo1 = map(o1,spin1) + do v2 = 1, nbf + vv2 = map(v2,spin2) + do v1 = 1, nbf c c write(6,*)' v1, v2, map(v1,spin1), map(v2,spin2): ', c & v1, v2, map(v1,spin1), map(v2,spin2) c write(6,*)' o1, o2, map(o1,spin1), map(o2,spin2): ', c & o1, o2, map(o1,spin1), map(o2,spin2) - t2(map(v1,spin1),map(v2,spin2), - & map(o1,spin1),map(o2,spin2)) = - & t_buf(v1,v2,o1,o2) + t2(map(v1,spin1),vv2,oo1,oo2) = t_buf(v1,v2,o1,o2) c end do end do @@ -2576,3 +2643,71 @@ c if (.not. ma_verify_allocator_stuff()) c $ call errquit('x',2) return end + subroutine uccsdt_flops_init() + implicit none +#include "cflopsP.fh" +c + call dfill(19, 0d0, flops, 1) + call dfill(19, 0d0, times, 1) + call dfill(19, 0d0, calls, 1) +c + end + subroutine uccsdt_flops_print() + implicit none +#include "cflopsP.fh" +#include "util.fh" +c + double precision overhead, start, mflops, used + integer ntime, i, j +c +c Estimate the timing overhead +c + start = util_cpusec() + ntime = 0 + do i = 1, 100 + do j = 1, 1000 + call uccsdt_flops_enter() + call uccsdt_flops_exit(19,0d0) + end do + ntime = ntime + 100 + overhead = util_cpusec() - start + if (overhead .gt. 1.0d0) goto 10 + enddo + 10 overhead = overhead / dble(ntime) + write(6,1) overhead + 1 format(//,'Triples kernel performance information ',//, + $ ' Timing overhead = ', 1p,d9.1,'s per call',/) +c + write(6,4) + 4 format(' Routine Calls FLOPs Time MFLOP/s',/, + $ ' ------- ------- -------- -------- --------') + do i = 1, 18 + used = times(i) - calls(i)*overhead + if (used .gt. 0d0) then + mflops = 1d-6*flops(i)/used + write(6,2) i, calls(i), flops(i), used, mflops + 2 format(1x,i5,2x,f9.0, 1x, 1p, d9.1, 0p, 1x, f9.2, 1x, f9.2) + end if + end do +c + end + subroutine uccsdt_flops_enter() + implicit none +#include "cflopsP.fh" +#include "util.fh" +c + enter = util_cpusec() +c + end + subroutine uccsdt_flops_exit(i, flopcount) + implicit none +#include "cflopsP.fh" +#include "util.fh" + integer i + double precision flopcount +c + times(i) = times(i) + util_cpusec() - enter + flops(i) = flops(i) + flopcount + calls(i) = calls(i) + 1 +c + end