mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 22:25:48 -04:00
initial caching of AO integrals for triples; SO integrals; performance measurement for the triples
This commit is contained in:
Robert Harrison
parent
820cc9430e
commit
375e612fd2
7 changed files with 1124 additions and 51 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
9
src/uccsdt/cflopsP.fh
Normal file
9
src/uccsdt/cflopsP.fh
Normal file
|
|
@ -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
|
||||
|
||||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
||||
|
|
|
|||
782
src/uccsdt/soints.F
Normal file
782
src/uccsdt/soints.F
Normal file
|
|
@ -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) = <ij|kl>
|
||||
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 <ij|kl>
|
||||
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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue