comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-29 06:35:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

clean push of changes necessary for COLUMBUS integration

Browse source
This commit is contained in:
Lachlan Thomas Belcher 2023-02-23 15:04:46 +00:00
parent 38c055a2a4
commit 09ca931164
38 changed files with 7317 additions and 29 deletions
Download patch file Download diff file Expand all files Collapse all files

24
src/NWints/aoints/GNUmakefile Normal file
View file

@ -0,0 +1,24 @@
# $Id
include ../../config/makefile.h
OBJ_OPTIMIZE = dft_aoints.o wrt_dft_aoints.o int_1e_sifs.o \
int_so_sifs.o int_2e_sifs.o int_2e_sifs_a.o \
sifs_2e_task.o int_2e_sifs_b.o nadct_trace.o \
rdhcid.o egrad_trace.o rd1mat.o sif2ga.o \
sif2da_2e_sort.o sif2arr.o rd_d2bl.o \
d2geri_trace.o dint_block_trc.o sequ.o \
sopgrdtrc.o asif2ga.o \
print_dint_block.o print_soblock.o trc_soblock.o \
int_mom_sifs.o hdoverlap.o
LIBRARY = libnwints.a
LIB_INCLUDES = -I../../nwdft/grid -I../../nwdft/include \
-I../../ddscf # I../../nwdft/util
include ../../config/makelib.h
# FOPTIONS += -Minform=inform -Mdclchk

92
src/NWints/aoints/asif2ga.F Normal file
View file

@ -0,0 +1,92 @@
subroutine asif2ga(g,locar,basis,nbft)
#include "global.fh"
#include "bas.fh"
#include "nwc_const.fh"
#include "errquit.fh"
#include "mafdecls.fh"
integer g ! atom blocked global array
double precision locar(*) ! local array from sifs with all values
integer basis
integer nbft
integer nshell
integer ishell, jshell
integer ijshell, ilo, ihi, jlo, jhi, blkdim
integer max1e, mem1
integer l_buf
integer k_buf
integer iloc, i, j, ijmap, il, jl
logical shells_ok
logical int_chk_sh
external int_chk_sh
cgk debug
* write(*,*)'gk: entered sif2ga'
cgk end
if (.not. bas_numcont(basis, nshell)) call errquit
& ('sif2ga: bas_numcont failed for basis', basis,
& BASIS_ERR)
call int_mem_1e(max1e, mem1)
cgk debug
* write(*,*)'gk: max1e = ', max1e
cgk end
if (.not. MA_push_get(MT_DBL,max1e,'sif2ga:buf',l_buf,k_buf))
$ call errquit('sif2ga: ma failed', max1e, MA_ERR)
call ga_zero(g)
cgk debug
* write(*,*)'gk: nshell = ', nshell
cgk end
ijmap = 0
do jshell = 1, nshell
* write(*,*)'gk: jshell = ', jshell
if (.not. bas_cn2bfr(basis, jshell, jlo, jhi))
$ call errquit('sif2ga: bas_cn2bfr jshell', jshell,
& BASIS_ERR)
* write(*,*)'gk: jlo = ', jlo, ' jhi = ', jhi
do ishell = 1, nshell
* write(*,*)'gk: ishell = ', ishell
if (.not. bas_cn2bfr(basis, ishell, ilo, ihi))
$ call errquit('sif2ga: bas_cn2bfr ishell', ishell,
& BASIS_ERR)
* write(*,*)'gk: ilo = ', ilo, ' ihi = ', ihi
call dfill(max1e, 0.0d0, dbl_mb(k_buf), 1)
blkdim=ihi-ilo+1
iloc=0
do j = jlo, jhi
do i = ilo, ihi
jl=max(i,j)
il=min(i,j)
ijmap=(jl-1)*jl/2+il
* write(*,2000) 'gk: j=',j,' i=',i,' ijmap=',ijmap,
* & ' iloc=',iloc
2000 format(4(a,i4))
* write(*,2010) 'gk: ',locar(ijmap)
2010 format(a,f20.12)
if(i.gt.j) then
* original phase
* dbl_mb(k_buf+iloc)=locar(ijmap)
dbl_mb(k_buf+iloc)=-locar(ijmap)
else
* original phase
* dbl_mb(k_buf+iloc)=-locar(ijmap)
dbl_mb(k_buf+iloc)=locar(ijmap)
endif
iloc=iloc+1
enddo
enddo
call ga_put(g,ilo,ihi,jlo,jhi,dbl_mb(k_buf),blkdim)
enddo
enddo
if (.not. MA_pop_stack(l_buf)) call errquit
$ ('sif2ga: pop failed', 0, GA_ERR)
return
end

104
src/NWints/aoints/d2geri_trace.F Normal file
View file

@ -0,0 +1,104 @@
subroutine d2geri_trace(mx_d2, d2, lscr, scr, lbuf, buf, nbft,
& b2s, nshell, sminb, smaxb, nsp, shoff, nat, frc_2el, ibas,
& dbg)
integer mx_d2
double precision d2(mx_d2)
integer lscr
double precision scr(lscr)
integer lbuf
double precision buf(lbuf)
integer nbft
double precision b2s(nbft)
integer nshell
integer sminb(nshell)
integer smaxb(nshell)
integer nsp
integer shoff(nsp,nsp)
integer nat
double precision frc_2el(3,nat)
integer ibas
integer dbg(mx_d2)
integer ish, jsh, ksh, lsh, ijsh, klsh
integer ilo, ihi, jlo, jhi, klo, khi, llo, lhi
integer slen, tlen, ulen, vlen
integer i, j, k, l
integer shf, blklen
integer idatom(4)
double precision fact
cgk debug
* write(*,*)'gk: in d2geri_trace'
* open (unit=27,form='formatted')
cgk end
do ksh=1,nshell
klo=sminb(ksh)
khi=smaxb(ksh)
ulen=khi-klo+1
do lsh=1,ksh
* klsh=nshell*(lsh-1)-((lsh-1)*lsh)/2+ksh
klsh=(ksh*(ksh-1))/2+lsh
llo=sminb(lsh)
lhi=smaxb(lsh)
vlen=lhi-llo+1
do ish=ksh,nshell
ilo=sminb(ish)
ihi=smaxb(ish)
slen=ihi-ilo+1
do jsh=1,ish
* ijsh=nshell*(jsh-1)-((jsh-1)*jsh)/2+ish
ijsh=(ish*(ish-1))/2+jsh
jlo=sminb(jsh)
jhi=smaxb(jsh)
tlen=jhi-jlo+1
if(ijsh.ge.klsh) then
cgk debug
* write(*,*)'#####################################'
* write(6,2312) 'ijsh=', ijsh, ' klsh=',klsh
* write(*,2321) 'd2geri s=',ish,' t=',jsh,' u=',ksh,' v=', lsh
* write(*,2321) 'ilo=',ilo,' ihi=',ihi,' jlo=',jlo,' jhi=',jhi
* write(*,2321) 'klo=',klo,' khi=',khi,' llo=',llo,' lhi=',lhi
cgk end
blklen=slen*tlen*ulen*vlen*4*3
if(blklen.gt.lbuf) write(*,*)'ERROR in lbuf'
* call dfill (lbuf, 0d0, buf, 1)
* call dfill (lscr, 0d0, scr, 1)
call intd_2e4c(ibas, ish, jsh, ibas, ksh, lsh, lscr, scr,
& lbuf, buf, idatom)
* call print_dint_block(ilo, ihi, jlo, jhi, klo, khi,
* & llo, lhi, buf, idatom)
shf=shoff(ijsh,klsh)
* write(6,'(a,i5)') 'shoff=',shf
2312 format(2(a,i2))
fact=5d-1
* if(ijsh.lt.klsh) fact=fact*2d0
if(ish.ne.jsh) fact=fact*2d0
if(ksh.ne.lsh) fact=fact*2d0
if(ish.ne.ksh.or.jsh.ne.lsh) fact=fact*2d0
* write(*,*)'gk: calling dint_block_trc'
call dint_block_trc(blklen, d2, ilo, ihi, jlo,
& jhi, klo, khi, llo, lhi, nshell, buf, idatom, nat,
& frc_2el, fact, shf, ish, jsh, ksh, lsh, nbft, b2s, dbg)
2321 format(4(a,i2))
* write(*,*)'gk: back from dint_block_trc'
endif
enddo
enddo
enddo
enddo
* write(*,*)'gk: leaving d2geri_trace'
* close(27)
return
end

380
src/NWints/aoints/dft_aoints.F Normal file
View file

@ -0,0 +1,380 @@
* logical function sodft_aoints(rtdb)
* implicit none
*#include "rtdb.fh"
* integer rtdb
*c
* logical nwdft
* external nwdft
* logical status
*c
*cgk deub
* write(*,*)'gk: in sodft_aoints in (nwdft.F) calling nwdft'
*cgk end
* status = rtdb_cput(rtdb,'dft:theory', 1, 'sodft')
* sodft_energy = nwdft(rtdb)
*cgk deub
* write(*,*)'gk: in sodft_aoints back from nwdft'
*cgk end
* call grid_cleanup(.true.)
*c
* end
logical function sodft_aoints(rtdb)
c
c>>> driver
c
C$Id: nwdft.F,v 1.76 2008/11/11 00:40:27 niri Exp $
implicit none
#include "errquit.fh"
c****
c**** nwchem handles
c****
integer rtdb
c
c--> Include standard input and memory checking common structures.
c
c
c>>> original NWdft includes ...
c
#include "mafdecls.fh"
#include "global.fh"
#include "tcgmsg.fh"
#include "stdio.fh"
c
c>>> includes added on addition of NWdft to nwchem
c
#include "bas.fh"
#include "geom.fh"
#include "rtdb.fh"
#include "inp.fh"
#include "cdft.fh"
#include "util.fh"
#include "cgridfile.fh"
#include "cosmo.fh"
c
c local declarations
c
double precision energy
integer iproc,itype
integer bases(3), i3
logical LResult, oprint, converged
logical dyall_mod_dir ! do modified Dirac calculation
c
logical wrt_dft_aoints
external wrt_dft_aoints
logical nadct_trace
external nadct_trace
logical dft_main0d, movecs_converged,grid_reopen,xc_gotxc
external dft_main0d, movecs_converged,grid_reopen,xc_gotxc
logical grid_ok,l1ecache
integer igok
character*80 theory
character*30 operation
c
cgk debug
* write(*,*)'gk: entered sodft_aoints...'
* write(*,*)'gk: do some hand shaking so we can transfer some of'
* & ,' this to aoints'
cgk end
sodft_aoints = .false.
iproc = ga_nodeid()
c
c Make sure database is parallel, push context down to DFT.
c
LResult = rtdb_parallel(.true.)
call util_print_push
call util_print_rtdb_load(rtdb, 'dft')
if (.not. rtdb_cget(rtdb, 'task:operation', 1, operation))
$ operation = ' '
if(.not.rtdb_cput(rtdb,'dft:theory', 1, 'sodft'))
& call errquit('dft_aoints: can not put theory on rtdb',
& 0,RTDB_ERR)
c
c init pstat
c
cgk debug
* write(*,*)'gk: calling dft_pstat_init from sodft_aoints'
cgk end
call dft_pstat_init(rtdb)
cgk debug
* write(*,*)'gk: back from dft_pstat_init from sodft_aoints'
cgk end
c
c init fdist
c
cgk debug
* write(*,*)'gk: calling dft_fdist_init from sodft_aoints'
cgk end
call dft_fdist_init(rtdb)
cgk debug
* write(*,*)'gk: back from dft_fdist_init from sodft_aoints'
cgk end
c
c Check for aoints...probably not neccessary now.
c
* if (.not. rtdb_cget(rtdb, 'task:operation',1,operation))
* $ call errquit('sodft_aoints: task:operation is not in rtdb',
* & 0,RTDB_ERR)
c
oprint = util_print('information', print_low)
if (oprint) then
if (.not. rtdb_cget(rtdb, 'title', 1, title))
& title = ' '
if(iproc.eq.0) then
call util_print_centered(
& LuOut, 'NWChem DFT Module', 40, .true.)
write(LuOut,*)
write(LuOut,*)
if (title .ne. ' ') then
call util_print_centered(LuOut, title, 40, .false.)
write(LuOut,*)
write(LuOut,*)
endif
endif
endif
c
c Gather input from all available resources e.g. input unit,
c old checkpoint files, etc.
c
cgk debug
* write(*,*)'gk: calling dft_rdinput from sodft_aoints'
cgk end
call dft_rdinput(rtdb)
cgk debug
* write(*,*)'gk: back from dft_rdinput from sodft_aoints'
cgk end
c
c ----- cosmo initialization ----
c
cosmo_last = .false.
if ( rtdb_get(rtdb,'slv:cosmo',mt_log,1,cosmo_on)) then
if(cosmo_on) then
call cosmo_initialize(rtdb,geom,ao_bas_han,oprint)
c
c Turn cosmo on, we want to run the calculation
c Start with gas_phase run
c
cosmo_last = .true.
cosmo_on = .true.
if(.not.rtdb_get(rtdb,'cosmo_phase',mt_int,1,cosmo_phase))
> cosmo_phase = 1
endif
endif
c
c Can we check here to see if system if already converged?
c
* if (movecs_in .eq. movecs_out) then
* if (.not. rtdb_cget(rtdb, 'dft:theory', 1, theory))
* $ call errquit('dft: theory not specified',0, RTDB_ERR)
* if(theory .eq. 'dft'.or.theory.eq.'sodft')then
* if(.not.xc_gotxc()) then
* grid_ok=.true.
* else
* grid_ok=grid_reopen(geom)
* endif
* igok=0
* if(grid_ok) igok=1
* call ga_igop(999,igok,1,'+')
* grid_ok=igok.eq.ga_nnodes()
* if(.not.grid_ok) call grid_cleanup(.true.)
* if (movecs_converged(rtdb, ao_bas_han,theory, movecs_in).
* & and.grid_ok) then
c
c check if we can get the grid from a file
c
* if(theory.eq.'sodft') then
* if (.not. rtdb_get(rtdb, 'sodft:energy', mt_dbl, 1,
* & energy)) call errquit(
* . 'sodft_aoints: failed getting converged energy',0,
* . RTDB_ERR)
* else
* if (.not. rtdb_get(rtdb, 'dft:energy', mt_dbl, 1,
* & energy)) call errquit(
* . 'sodft_aoints: failed getting converged energy',0,
* & RTDB_ERR)
* endif
* if (iproc.eq.0 .and. oprint)then
* write(LuOut,1011)energy
* call util_flush(Luout)
* endif
* converged = .true.
* sodft_aoints = .true.
* goto 1101
* endif
* endif
* endif
*1011 format(/' The DFT is already converged '//,
* & ' Total DFT energy =', f20.12/)
c
c Initialize the integral object
c
if(XCFIT.and.CDFIT) then
i3 = 3
bases(1) = AO_bas_han
bases(2) = CD_bas_han
bases(3) = XC_bas_han
elseif((.not.XCFIT).and.CDFIT) then
i3 = 2
bases(1) = AO_bas_han
bases(2) = CD_bas_han
elseif((.not.CDFIT).and.XCFIT) then
i3 = 2
bases(1) = AO_bas_han
bases(2) = XC_bas_han
else
i3 = 1
bases(1) = AO_bas_han
endif
c
cgk debug
* write(*,*)'gk: !!!!!!!!!!!!!!calling int_init from sodft_aoints'
* write(*,*)'gk: i3 = ', i3
cgk end
call int_init(rtdb, i3, bases)
cgk debug
* write(*,*)'gk: !!!!!!!!!!!!!!back from int_init from sodft_aoints'
cgk end
c
c get 1-el in cache
c
cgk debug
* write(*,*)'gk: !!!!!!!!!!!calling schwartz_init from sodft_aoints'
cgk end
c do we want Schwarz screening?
call schwarz_init(geom, AO_bas_han)
cgk debug
* write(*,*)'gk: !!!!!!!!!back from schwartz_init from sodft_aoints'
cgk end
c
l1ecache=.false.
* if ((nbf_ao*nbf_ao*3*8)/ga_nnodes() .lt. 8000000) then! Cache if < 8MB/node
* if (iproc.eq.0 .and. oprint)then
* write(LuOut,*) ' Caching 1-el integrals '
* call util_flush(Luout)
* endif
* l1ecache=.true.
* call int_1e_cache_ga(AO_bas_han, oskel)
* endif
c
c Build the grid.
c
grid_written=.false.
cgk debug
* write(*,*)'gk: calling get_setquad in sodft_aoints'
* write(*,*)'gk: building the grid'
cgk end
call grid_setquad(rtdb)
cgk debug
* write(*,*)'gk: back from get_setquad in sodft_aoints'
cgk end
c
c Analyze set parameters for consistency
c
cgk debug
* write(*,*)'gk: calling dft_inpana in sodft_aoints'
* write(*,*)'gk: analysing dft parameters '
cgk end
call dft_inpana(rtdb)
cgk debug
* write(*,*)'gk: back from dft_inpana in sodft_aoints'
cgk end
c
c
c--> determine system type (molecule, polymer,surface,crystal)
c if no system type exists in db --> default = molecular case
c
if(.not. geom_systype_get(geom,itype))
& itype = 0
c
if(itype.eq.0) then
#ifdef GA_TRACE
call trace_init(100000) ! initialize trace
#endif
if (inp_compare(.false., 'aoints', operation)) then
cgk debug
* write(*,*)'gk: calling wrt_dft_aoints from sodft_aoints'
cgk end
sodft_aoints = wrt_dft_aoints(rtdb)
cgk debug
* write(*,*)'gk: back from from sodft_aoints'
cgk end
else if(inp_compare(.false., 'grdtrc', operation)) then
call int_terminate
call intd_init(rtdb, i3, bases)
sodft_aoints = nadct_trace(rtdb,.false.,.false.)
else if(inp_compare(.false., 'nadtrc', operation)) then
call int_terminate
call intd_init(rtdb, i3, bases)
sodft_aoints = nadct_trace(rtdb,.true.,.true.)
#ifdef GA_TRACE
call trace_end(ga_nodeid()) ! end trace
#endif
else
call errquit('sodft_aoints: task not defined',0,
& INPUT_ERR)
endif
else
call errquit('sodft_aoints: problem with system type',0,
& INPUT_ERR)
endif
c
c terminate integrals
c
cgk debug
* write(*,*)'gk: starting dft cleanup in sodft_aoints'
cgk end
if(l1ecache) call int_1e_uncache_ga()
cgk doing this in nadctrc, etc.
c call schwarz_tidy
c call int_terminate
cgk end
c
c Go here if already converged
c
1101 continue
c
c Clean up cosmo stuff if necessary
c
if (cosmo_on) then
call cosmo_tidy()
cosmo_on = .false.
cosmo_phase = 1
endif
c
c destroy BS handles and geom object
c
XCFIT = .TRUE.
if( XC_bas_han.eq.-99999 ) XCFIT = .FALSE.
CDFIT = .TRUE.
if( CD_bas_han.eq.-99999 ) CDFIT = .FALSE.
if(XCFIT)then
if(.not.bas_destroy(XC_bas_han))
& call errquit ('dft:destroying XC basis handle',0, BASIS_ERR)
endif
if(CDFIT)then
if(.not.bas_destroy(CD_bas_han))
& call errquit ('dft:destroying CD basis handle',0, BASIS_ERR)
endif
if (.not.(
& bas_destroy(AO_bas_han)
& .and. geom_destroy(geom)))
& call errquit ('dft:destroying geom and basis handles',0,
& BASIS_ERR)
c
c--> Flush ouput buffer.
c
if(iproc.eq.0) call util_flush(LuOut)
call dft_pstat_print
call util_print_pop
cgk debug
* write(*,*)'gk: leaving sodft_aoints'
cgk end
c
return
c
end

180
src/NWints/aoints/dint_block_trc.F Normal file
View file

@ -0,0 +1,180 @@
subroutine dint_block_trc(blklen, d2blk, ilo, ihi, jlo, jhi,
& klo, khi, llo, lhi, nshell, geri, idatom, nat, frc_2el, fact,
& shoff, s, t, u, v, nbft, b2s, dbg)
integer blklen
double precision d2blk(*)
integer ilo, ihi, jlo, jhi, klo, khi, llo, lhi
integer nshell
double precision geri(llo:lhi, klo:khi, jlo:jhi, ilo:ihi, 3, 4)
integer idatom(4)
integer nat
double precision frc_2el(3,nat)
double precision fact
integer shoff
integer s, t, u, v
integer nbft
integer b2s(nbft)
integer slen, tlen, ulen, vlen, is, js, ks, ls
integer st, uv, a, b, c, d, ab, cd, abcd, stlen, ab1, cd1
integer i, j, k, l, ijkl
integer icart, iat
double precision dE
integer dbg(4,*)
integer i0, j0, k0, l0
logical index_err
character*2 deriv(3)
data deriv/'dx','dy','dz'/
double precision cntr
st=(s*(s-1))/2+t
uv=(u*(u-1))/2+v
slen=ihi-ilo+1
tlen=jhi-jlo+1
ulen=khi-klo+1
vlen=lhi-llo+1
* do iat=1,nat
* do icart=1,3
* frc_2el(icart,iat)=0d0
* enddo
* enddo
* write(6,2321) 'slen=',slen,' tlen=',tlen,' ulen=',ulen,' vlen=',
* & vlen
2321 format(4(a,i2))
* write(6,2312) '(st)=', st, ' (uv)=',uv
* write(6,'(a,4(i2,x))')'idatom: ', (idatom(iat), iat=1,4)
do iat = 1, 4
if(idatom(iat).ge.1) then
do icart = 1, 3
* write(*,'(a,i3)') '!! atom # = ', idatom(iat)
* write(*,'(a,i3)') '$$ icart = ', icart
dE = 0d0
do i=ilo, ihi
is=i-ilo+1
do j=jlo, jhi
js=j-jlo+1
* s=b2s(a)
* t=b2s(b)
* st=nshell*(t-1)-((t-1)*t)/2+s
if(s.eq.t) then
a=max(is,js)
b=min(is,js)
ab=((a-1)*a)/2+b
stlen=((slen+1)*slen)/2
else
a=is
b=js
* ab=(b-1)*slen+a
ab=(a-1)*tlen+b
stlen=slen*tlen
endif
do k=klo, khi
ks=k-klo+1
do l=llo, lhi
ls=l-llo+1
if(u.eq.v) then
c=max(ks,ls)
d=min(ks,ls)
cd=((c-1)*c)/2+d
else
c=ks
d=ls
* cd=(d-1)*ulen+c
cd=(c-1)*vlen+d
endif
if(st.eq.uv) then
ab1=max(ab,cd)
cd1=min(ab,cd)
abcd=stlen*(cd1-1)-((cd1-1)*cd1)/2+ab1
else
abcd=stlen*(cd-1)+ab
endif
ijkl=shoff+abcd
cgk debug
if(i.lt.j) then
i0=j
j0=i
else
i0=i
j0=j
endif
if(k.lt.l) then
k0=l
l0=k
else
k0=k
l0=l
endif
index_err=.false.
if(dbg(1,ijkl).ne.i0) index_err=.true.
if(dbg(2,ijkl).ne.j0) index_err=.true.
if(dbg(3,ijkl).ne.k0) index_err=.true.
if(dbg(4,ijkl).ne.l0) index_err=.true.
cgk I do not know how to switch bra and kets to make the error trap
cgk absolute
* if(index_err) then
* write(6,'(a)')'!!!WARNING Index error'
* write(6,2322) 'normal ordered [',i0,',',j0,';',k0,',',l0,']'
* write(6,2322) 'debug index[',dbg(1,ijkl),',',dbg(2,ijkl),';',
* & dbg(3,ijkl),',',dbg(4,ijkl),']'
* endif
2322 format(4(a,i2),a1)
* write(6,2322) '[ab;cd] = [',a,',',b,';',c,',',d,']'
* write(6,2312) 'ab=', ab, ' cd=',cd
2312 format(2(a,i2))
* write(6,'(a,i6)') ' abcd=', abcd
* write(6,'((a,i6))') 'ijkl=',ijkl
* write(6,2324) 'd2(',i,',',j,';',k,',',l,')=',d2blk(ijkl)
2324 format(4(a,i2),a,f20.12)
* write(6,1200) idatom(iat),deriv(icart),'geri=',
* & geri(l,k,j,i,icart,iat)
1200 format(i2,a2,a,f20.12)
cgk end
* df=1d0
* if(k.nq.l) df=df*2d0
* if(i.nq.j) df=df*2d0
* if(i.ne.j.and.k.ne.l) df=df*2d0
cntr=d2blk(ijkl)*geri(l,k,j,i,icart,iat)
dE = dE + cntr
if(icart.eq.1.and.abs(cntr).gt.1d-10) then
* write(27,2324) 'd2(',i,',',j,';',k,',',l,')=',d2blk(ijkl)
* write(27,1200) idatom(iat),deriv(icart),'geri=',
* & geri(l,k,j,i,icart,iat)
* write(27,*)'cntr = ', cntr
endif
enddo
enddo
enddo
enddo
frc_2el(icart,idatom(iat)) = frc_2el(icart,idatom(iat)) +
& dE*fact
enddo
endif
enddo
* if (st.lt.uv) then
* write(*,*)'st < uv'
* endif
* write(6,2200) '2 electron energy gradient',' ',' ',
* $ ((frc_2el(i,j),i=1,3),j=1,nat)
*2200 format(A,A,A/,1000(3(1x,F20.10),/))
cgk debug
* if(idatom(1).gt.0.or.idatom(2).gt.0.or.idatom(3).gt.0.or.
* & idatom(4).gt.0) then
* write(27,2222)'gk: shells: ', s, t, u, v
*2222 format(a,4i2)
* WRITE(27, 4020) ((frc_2el(i,j),i=1,3),j=1,nat)
*4020 FORMAT (6F12.6)
* endif
cgk end
return
end

710
src/NWints/aoints/egrad_trace.F Normal file
View file

@ -0,0 +1,710 @@
subroutine egrad_trace(ibas, aodens, ninfo, info, nbft, g_force,
& frc_nuc, frc_kin, frc_eff, frc_2el, frc_sow, frc_sox, frc_soy,
& frc_soz, frc_tot, nat, rtdb, nadct, aodensa, infoa, trans )
implicit none
#include "errquit.fh"
* #include "cint1cache.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "bas.fh"
#include "geom.fh"
#include "rtdb.fh"
#include "msgids.fh"
* #include "inp.fh"
* #include "apiP.fh"
* #include "cscfps.fh"
* #include "sym.fh"
* #include "cdft.fh"
#include "stdio.fh"
#include "cdft.fh"
integer ibas
integer*4 aodens
integer*4 ninfo
integer*4 info(ninfo) ! not to be confused with NWChem info
integer nbft
integer g_force, g_d1, g_eff
integer g_sod(3)
integer lforce
double precision frc_nuc(3,nat)
double precision frc_kin(3,nat)
double precision frc_eff(3,nat)
double precision frc_2el(3,nat)
double precision frc_sow(3,nat)
double precision frc_sox(3,nat)
double precision frc_soy(3,nat)
double precision frc_soz(3,nat)
double precision frc_tot(3,nat)
integer nat
integer rtdb
logical nadct
integer ga_create_atom_blocked
external ga_create_atom_blocked
! integer geom
! logical oskel
integer nshell
integer ishell, mem1, max1e, lrecal, n1mx
integer max2e, mscratch_1e, mscratch_2e, lbuf1, lscratch,
& max_at_bf, lsqatom
integer l_buf, l_scr, l_ilab, l_jlab, l_sifbuf, l_d1, l_eff,
& l_sifval, l_slab, l_sm, l_sym, l_map, l_b2s,
& l_smin, l_smax, l_sodx, l_sody, l_sodz, l_buf2, l_scr2
integer k_buf, k_scr, k_ilab, k_jlab, k_sifbuf, k_d1, k_eff,
& k_sifval, k_slab, k_sm, k_sym, k_map, k_b2s,
& k_smin, k_smax, k_sodx, k_sody, k_sodz, k_buf2, k_scr2
integer sizeof_int, szlabs, numd1, sz_bf_list
integer i, j, ij, k
integer me, nproc
integer max2e4c, mscr2e4c, lbuf2, lscr2
c
c sifs parameters
c
integer*4 l1rec, n1max, l2rec, n2max, ierr
integer*4 nsym, nmap
parameter(nsym=1, nmap=0)
* integer*4 ietype(mxenrgy)
integer*4 kntin(nsym), nmpsy(nsym)
* real*8 energy(nenrgy)
integer*4 itypea, itypeb
integer*4 numd1s, nbfs
integer*4 aodens2
logical status
INTEGER :: g_ad1, g_hd, l_asym, k_asym, l_hd, k_hd, nadx
DOUBLE PRECISION :: nadxfl(nat,3)
INTEGER :: hdol ! Handle or half-drv ovlp GA
INTEGER*4 :: aodensa, infoa(ninfo)
LOGICAL :: trans
INTEGER :: dims(3),chunk(3),iend(2),jstrt(3),jend(3)
INTEGER, PARAMETER :: istrt(2) = (/1,1/)
me = ga_nodeid()
nproc = ga_nnodes()
c
cgk debug
* write(*,*)'gk: entered egrad_trace'
* write(*,*)'gk: nbft = ', nbft
* write(*,*)'gk: nadct = ',nadct
cgk end
c
c Get info about the basis sets
c
cgk debug
* write(*,*)'gk: ibas = ', ibas
cgk end
if (.not. bas_numcont(ibas, nshell)) call errquit
$ ('egrad_trace: bas_numcont failed for ibas',
& ibas, BASIS_ERR)
cgk debug
* write(*,*)'gk: nshell = ', nshell
cgk end
l1rec = info(2)
n1max = info(3)
c allocate necessary local temporary arrays on the stack
c get memory for sifs 1e labels
sizeof_int=MA_sizeof(MT_INT, 1, MT_BYTE)
if(sizeof_int.eq.4) then
szlabs=(n1max+1)/2
sz_bf_list=(nbft+1)/2
elseif(sizeof_int.eq.8) then
szlabs=n1max
sz_bf_list=nbft
else
call errquit('wrt_dft_aoints: unsupporte integer type size=',
& sizeof_int, MA_ERR)
endif
c sifs integral label buffer (1 electron)
if (.not. MA_push_get(MT_Int,2*szlabs,'wrt_dft_aoints:ilab',
& l_slab, k_slab))
$ call errquit('wrt_dft_aoints: ma failed for sif labs',
& 2*szlabs, MA_ERR)
lrecal=l1rec
c get memory for sifs work buffer space
if (.not. MA_push_get(MT_DBL, lrecal,'egrad_trace:sifbuf',
& l_sifbuf, k_sifbuf))
$ call errquit('egrad_trace: ma failed getting sifbuf', lrecal,
& MA_ERR)
n1mx=n1max
c sifs value buffer
if (.not. MA_push_get(MT_DBL, n1mx,'egrad_trace:sifval',
& l_sifval, k_sifval))
$ call errquit('egrad_trace: ma failed getting sifval', n1mx,
& MA_ERR)
numd1=nbft*(nbft+1)/2
numd1s=numd1
cgk debug
* write(*,*)'gk: numd1=',numd1,' nbft=',nbft
cgk end
c sifs total effective density
if (.not. MA_push_get(MT_DBL, numd1,'egrad_trace:d1',
& l_sm, k_sm))
$ call errquit('egrad_trace: ma failed getting d1', numd1,
& MA_ERR)
c memory for symmetry labels (all one)
if (.not. MA_push_get(MT_Int, sz_bf_list,'egrad_trace:sym',
& l_sym, k_sym))
$ call errquit('egrad_trace: ma failed getting sym',
& sz_bf_list, MA_ERR)
c memory for maps
if (.not. MA_push_get(MT_Int, sz_bf_list,'egrad_trace:map',
& l_map, k_map))
$ call errquit('egrad_trace: ma failed getting map',
& sz_bf_list, MA_ERR)
nmpsy(1)=nbft
nbfs=nbft
if (.not. bas_geom(ibas, geom)) call errquit
$ ('egrad_trace: basis corrupt?', 0, BASIS_ERR)
c global 1e columbus matrices
c effective total non spin orbit density
g_d1 = ga_create_atom_blocked(geom, ibas,
& 'total eff density matrix')
c effective non spin orbit fock matrix
g_eff = ga_create_atom_blocked(geom, ibas,
& 'effective fock matrix')
C LB
c antisymmetric transition density
g_ad1 = ga_create_atom_blocked(geom, ibas,
& 'antisymmteric density matrix')
c half-derivative overlap, single derivative scratch
g_hd = ga_create_atom_blocked(geom, ibas,
& 'half-derivative overlap matrix')
c half-derivative overlap, all derivatives (3*# atoms)
dims(1) = 3*nat
dims(2) = nbft
dims(3) = nbft
chunk(1)= dims(1)
chunk(2)= -1
chunk(3)= -1
IF(.NOT. nga_create(MT_DBL,3,dims,
& 'half-derivative overlap GA',chunk,hdol))
$ call errquit('egrad_trace: failed to create hdol GA',
& 3*nat*nbft*nbft,GA_ERR)
CALL ga_zero(hdol)
CALL ga_print(hdol)
C LB
c Do this for node 0 only to avoid conflics reading aodens
if(me.eq.0) then
c read one electron density
itypea=0
itypeb=7
cgk debug
* write(*,*)'gk: getting effective density'
* write(*,*)'gk: calling rd1mat from egrad_trace'
* WRITE(*,*)"LB egrad_trace, l1rec=",l1rec,"n1max=",n1max,
* & "numd1s=",numd1s,"nsym=",nsym,"nmpsy=",nmpsy,"kntin=",kntin,
* & "nbfs=",nbfs
cgkd end
call rd1mat(aodens, l1rec, n1max, info, dbl_mb(k_sifbuf),
& int_mb(k_slab), dbl_mb(k_sifval), dbl_mb(k_sm), numd1s,
& int_mb(k_sym), int_mb(k_map), nsym, nmpsy, kntin, itypea,
& itypeb, nbfs)
cgk debug
* write(*,*)'gk: back from rd1mat'
* write(*,*)'gk: total effective density matrix:'
* ij=0
* do i=1,nbft
* do j=1,i
* write(*,'(2i4,f20.12)') i,j,dbl_mb(k_sm+ij)
* ij=ij+1
* enddo
* enddo
cgk end
c put effective density in global arrays
call sif2ga(g_d1,dbl_mb(k_sm),ibas,nbft)
c get sifs effective fock matrix
itypea=0
itypeb=8
call rd1mat(aodens, l1rec, n1max, info, dbl_mb(k_sifbuf),
& int_mb(k_slab), dbl_mb(k_sifval), dbl_mb(k_sm), numd1s,
& int_mb(k_sym), int_mb(k_map), nsym, nmpsy, kntin, itypea,
& itypeb, nbfs)
cgk debug
* write(*,*)'gk: back from rd1mat'
* write(*,*)'gk: total effective fock matrix:'
* ij=0
* do i=1,nbft
* do j=1,i
* write(*,'(2i4,f20.12)') i,j,dbl_mb(k_sm+ij)
* ij=ij+1
* enddo
* enddo
cgk end
c put effective density in global arrays
call sif2ga(g_eff,dbl_mb(k_sm),ibas,nbft)
C LB
c get antisymmetric density matrix
! USING MEMORY ALLOCATIONS FROM AODENS -- will this always work?
IF(trans) THEN
itypea=2
itypeb=9
call rd1mat(aodensa,infoa(2),infoa(3), info, dbl_mb(k_sifbuf),
& int_mb(k_slab), dbl_mb(k_sifval), dbl_mb(k_sm), numd1s,
& int_mb(k_sym), int_mb(k_map), nsym, nmpsy, kntin, itypea,
& itypeb, nbfs)
! WRITE(*,*)"LB egrad_trace, back from rd1mat for asym dens"
c put antisymmetric density in global arrays
call asif2ga(g_ad1,dbl_mb(k_sm),ibas,nbft)
* write(*,*)'gk: antisymmteric density matrix global array'
call ga_print(g_ad1)
ENDIF
* write(*,*)'gk: effective densitiy matrix global array'
call ga_print(g_d1)
* write(*,*)'gk: effective fock matrix global array'
call ga_print(g_eff)
C LB
end if
call ga_sync()
c
c using ~/src/gradients/grad_force.F as guide.
c
call schwarz_tidy()
call intd_terminate()
call int_init(rtdb, 1, ibas)
call schwarz_init(geom, ibas)
call int_terminate()
call intd_init(rtdb, 1, ibas)
max1e=0
max2e=0
mscratch_1e=0
mscratch_2e=0
cgk debug
call int_mem_1e(max1e,mscratch_1e)
* write(*,*)'gk: max1e=',max1e
* write(*,*)'gk: mscratch_1e=',mscratch_1e
cgk end
call int_mem(max1e, max2e, mscratch_1e, mscratch_2e)
cgk debug
* write(*,*)'gk: max1e=',max1e,' max2e=',max2e
* write(*,*)'gk: mscratch_1e=',mscratch_1e,' mscratch_2e=',
* & mscratch_2e
cgk end
* max2e = max(max2e,1296*100) ! 100 D quartets
lbuf1 = max(max1e, max2e)
lbuf1 = 3*2*lbuf1
lscratch = max(mscratch_1e,mscratch_2e)
cgk debug
* write(*,*)'gk: in egrad_trace before meme alloc'
* write(*,*)'gk: max2e = ', max2e
* write(*,*)'gk: lbuf1 = ', lbuf1
* write(*,*)'gk: lscratch = ', lscratch
cgk end
c buffers for one electron integral derivatives
if (.not. ma_push_get(mt_dbl,lbuf1,'deriv buffer',l_buf,k_buf))
$ call errquit('egrad_trace:could not allocate buffer',lbuf1,
& MA_ERR)
* call ma_summarize_allocated_blocks()
cgk debug
* write(*,*)'gk: l_scr=',l_scr,' k_scr=',k_scr,' lscratch=',
* & lscratch
cgk end
if (.not. ma_push_get(mt_dbl,lscratch,'deriv scratch',
$ l_scr, k_scr))
& call errquit('egrad_trace: scratch alloc failed',
$ lscratch, MA_ERR)
c local density matrix block
if (.not. bas_nbf_ce_max(ibas,max_at_bf))
$ call errquit('egrad_trace: could not get max_at_bf',0,
$ BASIS_ERR)
lsqatom = max_at_bf * max_at_bf
if (.not. ma_push_get(mt_dbl,lsqatom,'local_density',
$ l_d1,k_d1))
$ call errquit('egrad_trace:could not allocate l_d1',lsqatom,
& MA_ERR)
if (.not. ma_push_get(mt_dbl,lsqatom,'local_eff_fock',l_eff,
$ k_eff)) call errquit
& ('egrad_trace;could not allocate l_eff',lsqatom, MA_ERR)
oskel=.false.
c g_force does not appear to be used
call grad1(dbl_mb(k_buf), lbuf1, dbl_mb(k_scr), lscratch,
$ dbl_mb(k_d1), dbl_mb(k_eff), frc_nuc,
$ frc_kin, frc_eff, g_force,
$ g_d1, g_eff, ibas, geom, nproc, nat, max_at_bf,
$ rtdb, oskel)
if(me.eq.0) then
write(6,2200) 'kinetic energy gradient',' ',' ',
$ ((frc_kin(i,j),i=1,3),j=1,nat)
if(.not.nadct)
$ write(6,2200) 'nuclear repulsion gradient',' ',' ',
$ ((frc_nuc(i,j),i=1,3),j=1,nat)
write(6,2200) 'weighted density gradient',' ',' ',
$ ((frc_eff(i,j),i=1,3),j=1,nat)
2200 format(A,A,A/,1000(3(1x,F16.10),/))
endif
* goto 1000
C LB calculate half-derivative overlap matrices
IF (trans) THEN
! ENDIF
if (.not. bas_numbf(AO_bas_han, nbf_ao)) then
call errquit('Exiting from nadct_trace',1, BASIS_ERR)
endif
CALL hdoverlap(dbl_mb(k_buf), lbuf1, lscratch, AO_bas_han, geom,
& nproc, nat, max1e, mem1, hdol, nbft)
nadxfl=0D0
iend(1) =nbft
iend(2) =nbft
jstrt(2)=1
jstrt(3)=1
jend(2) =nbft
jend(3) =nbft
!CALL ga_print(hdol)
DO i = 0, nat-1 !atom number
DO j = 0, 2 !derivative direction
jstrt(1)=i*3+j+1
jend(1) =i*3+j+1
!WRITE(*,*)"In egrad_trace, jstrt=",jstrt(1)
CALL nga_copy_patch('N',hdol,jstrt,jend,g_hd,istrt,iend)
CALL ga_print(g_hd)
nadxfl(i+1,j+1)=-2*ga_ddot(g_hd,g_ad1)
ENDDO
ENDDO
!CALL ga_print(hdol)
WRITE(LUOUT,*)"nadxfl :"
CALL output(nadxfl,1,nat,1,3,nat,3,1)
OPEN(nadx,file='nadxfl',form='formatted')
WRITE(nadx,1010) (nadxfl(i,:), i=1,nat)
1010 FORMAT (3d15.6)
ENDIF
C LB
c spin orbit contributions
do i=1, 3
status = ga_create(mt_dbl,nbft,nbft,'so density matrix', 0, 0,
& g_sod(i))
if (.not. status) then
if (ga_nodeid() .eq. 0) then
write(6,*) ' egrad_trace: ', 'so density matrix'
call util_flush(6)
endif
call ga_sync()
call errquit('egrad_trace: ga_create ', 0, GA_ERR)
endif
call ga_zero(g_sod(i))
end do
call ga_zero(g_eff) ! reuse effective density matrix
if (.not. ma_push_get(mt_dbl,lsqatom,'local sox den',l_sodx,
$ k_sodx)) call errquit
& ('egrad_trace;could not allocate l_sodx',lsqatom, MA_ERR)
if (.not. ma_push_get(mt_dbl,lsqatom,'local sox den',l_sody,
$ k_sody)) call errquit
& ('egrad_trace;could not allocate l_sody',lsqatom, MA_ERR)
if (.not. ma_push_get(mt_dbl,lsqatom,'local soz den',l_sodz,
$ k_sodz)) call errquit
& ('egrad_trace;could not allocate l_sodz',lsqatom, MA_ERR)
if(me.eq.0) then
itypea=2
itypeb=11
cgk debug
* write(*,*)'gk: x spin orbit density'
cgkd end
call rd1mat(aodens, l1rec, n1max, info, dbl_mb(k_sifbuf),
& int_mb(k_slab), dbl_mb(k_sifval), dbl_mb(k_sm), numd1s,
& int_mb(k_sym), int_mb(k_map), nsym, nmpsy, kntin, itypea,
& itypeb, nbfs)
cgk debug
* write(*,*)'gk: x spin orbit density matrix matrix:'
* ij=0
* do i=1,nbft
* do j=1,i
* write(*,'(2i4,f20.12)') i,j,dbl_mb(k_sm+ij)
* ij=ij+1
* enddo
* enddo
cgk end
c put x so effective density in global arrays
call asif2ga(g_sod(1),dbl_mb(k_sm),ibas,nbft)
write(*,*)'x so effective density global array'
call ga_print(g_sod(1))
itypea=2
itypeb=12
cgk debug
* write(*,*)'gk: y spin orbit density'
cgkd end
call rd1mat(aodens, l1rec, n1max, info, dbl_mb(k_sifbuf),
& int_mb(k_slab), dbl_mb(k_sifval), dbl_mb(k_sm), numd1s,
& int_mb(k_sym), int_mb(k_map), nsym, nmpsy, kntin, itypea,
& itypeb, nbfs)
cgk debug
* write(*,*)'y spin orbit density matrix matrix:'
* ij=0
* do i=1,nbft
* do j=1,i
* write(*,'(2i4,f20.12)') i,j,dbl_mb(k_sm+ij)
* ij=ij+1
* enddo
* enddo
cgk end
c put y so effective density in global arrays
call asif2ga(g_sod(2),dbl_mb(k_sm),ibas,nbft)
write(*,*)'y so effective density global array'
call ga_print(g_sod(2))
itypea=2
itypeb=13
cgk debug
* write(*,*)'gk: z spin orbit density'
cgkd end
call rd1mat(aodens, l1rec, n1max, info, dbl_mb(k_sifbuf),
& int_mb(k_slab), dbl_mb(k_sifval), dbl_mb(k_sm), numd1s,
& int_mb(k_sym), int_mb(k_map), nsym, nmpsy, kntin, itypea,
& itypeb, nbfs)
cgk debug
* write(*,*)'gk: z spin orbit dinsity matrix matrix:'
* ij=0
* do i=1,nbft
* do j=1,i
* write(*,'(2i4,f20.12)') i,j,dbl_mb(k_sm+ij)
* ij=ij+1
* enddo
* enddo
cgk end
c put z so effective density in global arrays
call asif2ga(g_sod(3),dbl_mb(k_sm),ibas,nbft)
write(*,*)'z so effective density global array'
call ga_print(g_sod(3))
c get spin orbit effective fock matrix
itypea=1
itypeb=35
call rd1mat(aodens, l1rec, n1max, info, dbl_mb(k_sifbuf),
& int_mb(k_slab), dbl_mb(k_sifval), dbl_mb(k_sm), numd1s,
& int_mb(k_sym), int_mb(k_map), nsym, nmpsy, kntin, itypea,
& itypeb, nbfs)
cgk debug
* write(*,*)'gk: total spin orbit effective fock matrix:'
* ij=0
* do i=1,nbft
* do j=1,i
* write(*,'(2i4,f20.12)') i,j,dbl_mb(k_sm+ij)
* ij=ij+1
* enddo
* enddo
cgk end
c put effective density in global arrays
call sif2ga(g_eff,dbl_mb(k_sm),ibas,nbft)
cgk debug
* write(*,*)'spin orbit effective fock matrix global array'
* call ga_print(g_eff)
cgk end
end if
cgk debug
* write(*,*)'gk:calling sopgrdtrc'
cgk end
call sopgrdtrc(dbl_mb(k_buf),lbuf1,dbl_mb(k_scr),lscratch,
& dbl_mb(k_eff),frc_sow,g_sod,g_eff,ibas,geom,nproc,nat,max_at_bf,
& oskel,frc_sox,frc_soy,frc_soz, dbl_mb(k_sodx),dbl_mb(k_sody),
& dbl_mb(k_sodz))
cgk debug
* write(*,*)'gk: back from sopgrdtrc'
cgk end
if(me.eq.0) then
write(6,2201) 'SOx gradient',' ',' ',
$ ((frc_sox(i,j),i=1,3),j=1,nat)
write(6,2201) 'SOy gradient',' ',' ',
$ ((frc_soy(i,j),i=1,3),j=1,nat)
write(6,2201) 'SOz gradient',' ',' ',
$ ((frc_soz(i,j),i=1,3),j=1,nat)
write(6,2201) 'SO F.dS gradient',' ',' ',
$ ((frc_sow(i,j),i=1,3),j=1,nat)
endif
2201 format(A,A,A/,1000(3(1x,E18.6),/))
do j=1, nat
do i=1,3
frc_sow(i,j)=-frc_sow(i,j)+frc_sox(i,j)+frc_soy(i,j)+
& frc_soz(i,j)
enddo
enddo
if(me.eq.0) then
write(6,2201) 'total SO gradient',' ',' ',
$ ((frc_sow(i,j),i=1,3),j=1,nat)
endif
c chop stack at first item allocated
if (.not. MA_pop_stack(l_sodz)) call errquit
$ ('egrad_trace: pop failed at l_sodz', 0, GA_ERR)
if (.not. MA_pop_stack(l_sody)) call errquit
$ ('egrad_trace: pop failed at l_sody', 0, GA_ERR)
if (.not. MA_pop_stack(l_sodx)) call errquit
$ ('egrad_trace: pop failed at l_sodx', 0, GA_ERR)
c jump to here of no spin orbit
1000 continue
call schwarz_tidy()
call intd_terminate()
call int_init(rtdb, 1, ibas)
call schwarz_init(geom, ibas)
call int_terminate()
call intd_init(rtdb, 1, ibas)
if(me.eq.0) then
call sifo2f(aodens, aodens, 'aodens', info, aodens2, ierr)
endif
l2rec = info(4)
n2max = info(5)
c memory for basis function to shell mapping array
if (.not. ma_push_get(mt_dbl,nbft,'bf2shell_map',l_b2s,
$ k_b2s)) call errquit
& ('egrad_trace;could not allocate l_b2s',nbft, MA_ERR)
if (.not. ma_push_get(mt_dbl,nshell,'shell_min',l_smin,
$ k_smin)) call errquit
& ('egrad_trace;could not allocate l_smin',nshell, MA_ERR)
if (.not. ma_push_get(mt_dbl,nshell,'shell_min',l_smax,
$ k_smax)) call errquit
& ('egrad_trace;could not allocate l_smax',nshell, MA_ERR)
call int_mem_2e4c(max2e4c, mscr2e4c) ! blocking algorithm
* write(*,*)'gk: max2e4c = ', max2e4c
* write(*,*)'gk: mscr2e4c = ', mscr2e4c
lbuf2=3*4*max2e4c
lscr2=mscr2e4c
if (.not. ma_push_get(mt_dbl,lbuf2,'buf2',l_buf2,
$ k_buf2)) call errquit
& ('egrad_trace;could not allocate l_b2s',lbuf2, MA_ERR)
if (.not. ma_push_get(mt_dbl,lscr2,'buf2',l_scr2,
$ k_scr2)) call errquit
& ('egrad_trace;could not allocate l_b2s',lscr2, MA_ERR)
c two electron trace
call sif2da_2e_sort(aodens, l2rec, n2max, info, ibas, nbft,
& nshell, dbl_mb(k_b2s), dbl_mb(k_smin), dbl_mb(k_smax), frc_2el,
& nat, lbuf2, dbl_mb(k_buf2), lscr2, dbl_mb(k_scr2) )
call schwarz_tidy()
call intd_terminate()
if(nadct) then
do j=1, nat
do i=1,3
frc_tot(i,j)=frc_kin(i,j)+frc_eff(i,j)+
& +frc_sow(i,j)+frc_2el(i,j)
enddo
enddo
else
do j=1, nat
do i=1,3
frc_tot(i,j)=frc_kin(i,j)+frc_nuc(i,j)+frc_eff(i,j)+
& +frc_sow(i,j)+frc_2el(i,j)
enddo
enddo
endif
if(me.eq.0) then
write(6,2200) '2 electron energy gradient',' ',' ',
$ ((frc_2el(i,j),i=1,3),j=1,nat)
write(6,2200) 'total gradient',' ',' ',
$ ((frc_tot(i,j),i=1,3),j=1,nat)
write(6,2201) 'sci total gradient',' ',' ',
$ ((frc_tot(i,j),i=1,3),j=1,nat)
endif
if (.not. MA_pop_stack(l_scr2)) call errquit
$ ('egrad_trace: pop failed at l_scr2', 0, GA_ERR)
if (.not. MA_pop_stack(l_buf2)) call errquit
$ ('egrad_trace: pop failed at l_buf2', 0, GA_ERR)
if (.not. MA_pop_stack(l_smax)) call errquit
$ ('egrad_trace: pop failed at l_smax', 0, GA_ERR)
if (.not. MA_pop_stack(l_smin)) call errquit
$ ('egrad_trace: pop failed at l_smin', 0, GA_ERR)
if (.not. MA_pop_stack(l_b2s)) call errquit
$ ('egrad_trace: pop failed at l_b2s', 0, GA_ERR)
if (.not. MA_pop_stack(l_eff)) call errquit
$ ('egrad_trace: pop failed at l_eff', 0, GA_ERR)
if (.not. MA_pop_stack(l_d1)) call errquit
$ ('egrad_trace: pop failed at l_d1', 0, GA_ERR)
if (.not. MA_pop_stack(l_scr)) call errquit
$ ('egrad_trace: pop failed at l_scr', 0, GA_ERR)
if (.not. MA_pop_stack(l_buf)) call errquit
$ ('egrad_trace: pop failed at l_buf', 0, GA_ERR)
if (.not. MA_pop_stack(l_map)) call errquit
$ ('egrad_trace: pop failed at l_map', 0, GA_ERR)
if (.not. MA_pop_stack(l_sym)) call errquit
$ ('egrad_trace: pop failed at l_sym', 0, GA_ERR)
if (.not. MA_pop_stack(l_sm)) call errquit
$ ('egrad_trace: pop failed at l_sm', 0, GA_ERR)
if (.not. MA_pop_stack(l_sifval)) call errquit
$ ('egrad_trace: pop failed at l_sifval', 0, GA_ERR)
if (.not. MA_pop_stack(l_sifbuf)) call errquit
$ ('egrad_trace: pop failed at l_sifbuf', 0, GA_ERR)
if (.not. MA_pop_stack(l_slab)) call errquit
$ ('egrad_trace: pop failed at l_slab', 0, GA_ERR)
* if (.not. MA_pop_stack(l_bfnlab)) call errquit
* if (.not. MA_pop_stack(l_bfnlab)) call errquit
* $ ('egrad_trace: pop failed at l_bfnlab', 0, GA_ERR)
* if (.not. MA_pop_stack(l_info)) call errquit
* $ ('egrad_trace: pop failed at l_info', 0, GA_ERR)
* if (.not. MA_pop_stack(l_jlab)) call errquit
* $ ('egrad_trace: pop failed at l_jlab', 0, GA_ERR)
* if (.not. MA_pop_stack(l_ilab)) call errquit
* $ ('egrad_trace: pop failed at l_ilab', 0, GA_ERR)
* if (.not. MA_pop_stack(l_scr)) call errquit
* $ ('egrad_trace: pop failed at l_scr', 0, GA_ERR)
* if (.not. MA_pop_stack(l_buf)) call errquit
* $ ('egrad_trace: pop failed at l_buf', 0, GA_ERR)
c
c ----- restore ocache -----
c
call ga_sync() ! So that no nasty races can result
c
cgk debug
* write(*,*)'gk: leaving egrad_trace at end of routine'
cgk end
end

424
src/NWints/aoints/grad1.F Normal file
View file

@ -0,0 +1,424 @@
subroutine grad1 ( H, lbuf, scr, lscr, dens, wdens, frc_nuc,
$ frc_kin, frc_wgh, g_force,
$ g_dens, g_wdens, basis, geom, nproc, nat,
$ max_at_bf, rtdb, oskel )
c$Id: grad1.F,v 1.17 2008/04/22 22:00:19 niri Exp $
C one electron contribution to RHF, ROHF and UHF gradients
C now also UMP2
implicit none
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "bas.fh"
#include "rtdb.fh"
#include "sym.fh"
#include "bq_params.fh"
C-------------------------parameters--------------------------------
integer lbuf, lscr,
$ g_dens, ! density matrix (summed if ROHF, UHF)
$ g_wdens, ! weighted density (Lagrangian)
$ g_force, ! global force array
$ basis, geom, nproc, nat, max_at_bf, rtdb
double precision H, ! integral derivatives
$ scr,
$ dens, ! local density block
$ wdens, ! local weighted density block
$ frc_nuc, frc_kin, frc_wgh ! forces arrays
dimension H ( lbuf ), frc_nuc(3, nat), frc_kin(3, nat),
$ frc_wgh(3, nat), scr(lscr),
$ dens(max_at_bf,max_at_bf), wdens(max_at_bf,max_at_bf)
logical oskel ! symmetry?
C-------------------------local variables--------------------------
integer ijatom, next, iat1, iat2, iat3, ish1, ish2,
$ iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l,
$ if1, il1, if2, il2,
$ icart, ic, nint, ip1, ip2
double precision crd1, crd2 ! atomic coordinates
dimension crd1(3), crd2(3)
integer idatom
dimension idatom(2)
double precision dE, dx, dy, dz, qfac, fact, q1, q2
logical status, pointforce,dobq
character*16 name
integer bq_ncent
integer i_qbq,i_cbq
double precision r12
c
integer nxtask, task_size
external nxtask
task_size = 1
status = rtdb_parallel(.true.) ! Broadcast reads to all processes
pointforce = geom_include_bqbq(geom)
dobq = geom_extbq_on()
call hf_print_set(1)
cgk debug
* write(*,*)'gk: in grad1, dobq = ', dobq
* write(*,*)'gk: in grad1, pointforce = ', pointforce
cgk end
ijatom = -1
next = nxtask(nproc,task_size)
do 90, iat1 = 1, nat
do 80, iat2 = 1, iat1
ijatom = ijatom + 1
if ( ijatom .eq. next ) then
status = bas_ce2bfr(basis,iat1,iab1f,iab1l)
status = bas_ce2bfr(basis,iat2,iab2f,iab2l)
if (iab1f.le.0 .or. iab2f.le.0) then
c
c At least one center has no functions on it ... next atom
c
goto 1010
endif
if (oskel) then
if (.not. sym_atom_pair(geom, iat1, iat2, qfac))
$ goto 1010
else
qfac = 1.0d0
endif
status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
call ga_get (g_dens, iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
call ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf)
do 70, ish1 = iac1f, iac1l
if ( iat1.eq.iat2 ) iac2l = ish1
do 60, ish2 = iac2f, iac2l
C shell block in atomic (D/Dw)-matrix block
status = bas_cn2bfr(basis,ish1,if1,il1)
if1 = if1 - iab1f + 1
il1 = il1 - iab1f + 1
status = bas_cn2bfr(basis,ish2,if2,il2)
if2 = if2 - iab2f + 1
il2 = il2 - iab2f + 1
nint = ( il1 - if1 + 1 ) * ( il2 - if2 + 1 )
C overlap derivatives
call intd_1eov(basis,ish1,basis,ish2,lscr,scr,
& lbuf,H,idatom)
C Dw x S
if ( idatom(1) .ge. 1 ) then
C idatom(1).ge.0 <=> idatom(2).ge.0 (no check necessary)
ic = 1
do 28, icart = 1, 3
de = 0.D0
do 22, ip1 = if1, il1
do 20, ip2 = if2, il2
dE = dE + wdens(ip1,ip2) * H(ic)
ic = ic + 1
20 continue
22 continue
dE = dE * qfac
frc_wgh(icart,idatom(1)) = frc_wgh(icart,idatom(1))
$ - dE - dE
frc_wgh(icart,idatom(2)) = frc_wgh(icart,idatom(2))
$ + dE + dE
28 continue
endif
C 1el. derivatives
if(.not.dobq) then
call intd_1eh1(basis,ish1,basis,ish2,lscr,scr,
& lbuf,H)
else
call intd_1epot(basis,ish1,basis,ish2,lscr,scr,
& lbuf,H)
end if
C D x H
ic=1
do 50, iat3 = 1, nat
do 40, icart = 1, 3
dE = 0.D0
do 31, ip1 = if1, il1
do 30, ip2 = if2, il2
dE = dE + dens(ip1,ip2) * H(ic)
ic = ic + 1
30 continue
31 continue
if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
dE = dE * qfac
frc_kin(icart,iat3) = frc_kin(icart,iat3) + dE
40 continue
50 continue
60 continue
70 continue
1010 continue
C Vnn
if ( iat1 .NE. iat2 ) then
if (iab1f.ne.0 .or. iab2f.ne.0 .or. pointforce ) then
C no forces between point charges (for John Nicholas)
status = geom_cent_get (geom, iat1, name, crd1, q1)
status = geom_cent_get (geom, iat2, name, crd2, q2)
dx = crd2(1) - crd1(1)
dy = crd2(2) - crd1(2)
dz = crd2(3) - crd1(3)
fact = q1 * q2 / SQRT ( dx*dx + dy*dy + dz*dz ) **3
dE = dx * fact
frc_nuc(1,iat1) = frc_nuc(1,iat1) + dE
frc_nuc(1,iat2) = frc_nuc(1,iat2) - dE
dE = dy * fact
frc_nuc(2,iat1) = frc_nuc(2,iat1) + dE
frc_nuc(2,iat2) = frc_nuc(2,iat2) - dE
dE = dz * fact
frc_nuc(3,iat1) = frc_nuc(3,iat1) + dE
frc_nuc(3,iat2) = frc_nuc(3,iat2) - dE
endif
endif
next = nxtask(nproc,task_size)
endif
80 continue
90 continue
next = nxtask(-nproc,task_size)
c
c continue with forces if there are external charges (MV)
c --------------------------------------------------
if(.not.dobq) return
cgk debug
* write(*,*)'gk: continuing with forces for external charges'
cgk end
bq_ncent = geom_extbq_ncenter()
i_cbq = geom_extbq_coord()
i_qbq = geom_extbq_charge()
ijatom = -1
next = nxtask(nproc,task_size)
do 91, iat1 = 1, nat
do 81, iat2 = 1, bq_ncent
ijatom = ijatom + 1
if ( ijatom .eq. next ) then
status = geom_cent_get (geom, iat1, name, crd1, q1)
q2 = dbl_mb(i_qbq+iat2-1)
crd2(1) = dbl_mb(i_cbq+(iat2-1)*3)
crd2(2) = dbl_mb(i_cbq+(iat2-1)*3+1)
crd2(3) = dbl_mb(i_cbq+(iat2-1)*3+2)
dx = crd2(1) - crd1(1)
dy = crd2(2) - crd1(2)
dz = crd2(3) - crd1(3)
c == check for small distances (the regions could overlap) ==
r12 = sqrt(dx*dx + dy*dy + dz*dz)
if (r12.gt.bq_smalldist) then
fact = q1 * q2 / SQRT ( dx*dx + dy*dy + dz*dz ) **3
dE = dx * fact
frc_nuc(1,iat1) = frc_nuc(1,iat1) + dE
dE = dy * fact
frc_nuc(2,iat1) = frc_nuc(2,iat1) + dE
dE = dz * fact
frc_nuc(3,iat1) = frc_nuc(3,iat1) + dE
end if
c
next = nxtask(nproc,task_size)
endif
81 continue
91 continue
next = nxtask(-nproc,task_size)
cgk debug
* write(*,*)'gk: leaving grad1'
cgk end
return
end
subroutine grad_hnd_cos ( H, lbuf, scr, lscr,
$ dens, wdens, frc_nuc,
$ frc_kin, frc_wgh, g_force,
$ g_dens, g_wdens, basis, geom, nproc, nat,
$ max_at_bf, rtdb, oskel )
c$Id: grad1.F,v 1.17 2008/04/22 22:00:19 niri Exp $
C one electron contribution to RHF, ROHF and UHF gradients
C now also UMP2
implicit none
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "bas.fh"
#include "rtdb.fh"
#include "sym.fh"
#include "stdio.fh"
C-------------------------parameters--------------------------------
integer lbuf, lscr,
$ g_dens, ! density matrix (summed if ROHF, UHF)
$ g_wdens, ! weighted density (Lagrangian)
$ g_force, ! global force array
$ basis, geom, nproc, nat, max_at_bf, rtdb
double precision H, ! integral derivatives
$ scr,
$ dens, ! local density block
$ wdens, ! local weighted density block
$ frc_nuc, frc_kin, frc_wgh ! forces arrays
dimension H ( lbuf ), frc_nuc(3, nat), frc_kin(3, nat),
$ frc_wgh(3, nat), scr(lscr),
$ dens(max_at_bf,max_at_bf), wdens(max_at_bf,max_at_bf)
logical oskel ! symmetry?
C-------------------------local variables--------------------------
integer ijatom, next, iat1, iat2, iat3, ish1, ish2,
$ iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l,
$ if1, il1, if2, il2,
$ icart, ic, nint, ip1, ip2
double precision dE, qfac
logical status, pointforce
integer nxtask, task_size
external nxtask
c ---- -cosmo- gradient term -----
logical odbug
odbug=.true.
if(odbug) then
write(Luout,*) 'in -grad1_hnd_cos- ...'
endif
c
task_size = 1
status = rtdb_parallel(.true.) ! Broadcast reads to all processes
pointforce = geom_include_bqbq(geom)
call hf_print_set(1)
ijatom = -1
next = nxtask(nproc,task_size)
do 90, iat1 = 1, nat
do 80, iat2 = 1, iat1
ijatom = ijatom + 1
if ( ijatom .eq. next ) then
status = bas_ce2bfr(basis,iat1,iab1f,iab1l)
status = bas_ce2bfr(basis,iat2,iab2f,iab2l)
if (iab1f.le.0 .or. iab2f.le.0) then
c
c At least one center has no functions on it ... next atom
c
goto 1010
endif
if (oskel) then
if (.not. sym_atom_pair(geom, iat1, iat2, qfac))
$ goto 1010
else
qfac = 1.0d0
endif
status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
call ga_get (g_dens, iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
call ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf)
do 70, ish1 = iac1f, iac1l
if ( iat1.eq.iat2 ) iac2l = ish1
do 60, ish2 = iac2f, iac2l
C shell block in atomic (D/Dw)-matrix block
status = bas_cn2bfr(basis,ish1,if1,il1)
if1 = if1 - iab1f + 1
il1 = il1 - iab1f + 1
status = bas_cn2bfr(basis,ish2,if2,il2)
if2 = if2 - iab2f + 1
il2 = il2 - iab2f + 1
nint = ( il1 - if1 + 1 ) * ( il2 - if2 + 1 )
ic=1
do iat3 = 1, nat
do icart = 1, 3
do ip1 = if1, il1
do ip2 = if2, il2
H(ic)=0.D0
ic = ic + 1
enddo
enddo
enddo
enddo
C 1el. -cosmo- derivatives
c- call intd_1eh1(basis,ish1,basis,ish2,lscr,scr,
c- & lbuf,H)
C D x H
ic=1
do 50, iat3 = 1, nat
do 40, icart = 1, 3
dE = 0.D0
do 31, ip1 = if1, il1
do 30, ip2 = if2, il2
dE = dE + dens(ip1,ip2) * H(ic)
ic = ic + 1
30 continue
31 continue
if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
dE = dE * qfac
frc_kin(icart,iat3) = frc_kin(icart,iat3) + dE
40 continue
50 continue
60 continue
70 continue
1010 continue
next = nxtask(nproc,task_size)
endif
80 continue
90 continue
next = nxtask(-nproc,task_size)
return
end

281
src/NWints/aoints/hdoverlap.F Normal file
View file

@ -0,0 +1,281 @@
subroutine hdoverlap ( H, l_buf, l_scr, basis, geom, nproc,
& nat, max1e, mem1, hdol, nmot)
C LB based loosely on grad1. Program to calculate half-derivative overlap
C integrals, 1/2[<i|j^x>-<i^x|j>]. This capability already exists
C within intd_1eov, it just needed some teasing out.
implicit none
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "bas.fh"
#include "rtdb.fh"
#include "sym.fh"
#include "bq_params.fh"
!#include "cdft.fh"
#include "stdio.fh"
#include "errquit.fh"
#include "nwc_const.fh"
#include "basP.fh"
C-------------------------parameters--------------------------------
INTEGER :: basis, lbus, max1e, mem1, nproc, nat, nbf
INTEGER :: nmot
INTEGER :: l_buf, l_scr, geom
INTEGER :: hdol !handle for hdol GA
double precision :: H ! integral derivatives
DOUBLE PRECISION :: scr ! scratch
DOUBLE PRECISION :: SQHDOL (3*nat,nmot,nmot)
dimension H ( l_buf ), scr(l_scr)
logical oskel ! symmetry?
C-------------------------local variables--------------------------
integer ijatom, next, iat1, iat2, iat3, ish1, ish2,
$ iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l,
$ if1, il1, if2, il2,
$ icart, ic, nint, ip1, ip2
double precision crd1, crd2 ! atomic coordinates
dimension crd1(3), crd2(3)
integer idatom
dimension idatom(2)
double precision qfac
logical status, pointforce,dobq
character*16 name
integer bq_ncent
integer i_qbq,i_cbq
double precision r12
c
integer nxtask, task_size
external nxtask
LOGICAL :: hlfdrv !true = calculate half-derivative overlap integrals
INTEGER :: atmidx1, atmidx2, nmss2, idx3 ! indices the SQHDOL matrices
INTEGER :: contf, contl, shlf, shll ! more indices
INTEGER :: orbstrt(nat) ! gives the orbital label where each atom
! center starts. There may be a better
! built-in way to achieve this
INTEGER :: i,j,k ! counters
INTEGER :: itype, jtype ! angular momentum of contraction
INTEGER :: nintsph ! nint after conversion to spherical
INTEGER :: dum1, dum2, dum3 ! dummy variables
INTEGER :: dims(3)
INTEGER, PARAMETER :: istrt(3) = (/1,1,1/)
c LB
! WRITE(LUOUT,*)"LB Debug, now in hdoverlap"
SQHDOL=0.0D0
! WRITE(LUOUT,*)"LB Debug, iab1f=",iab1f,"iab2f=",iab2f
! WRITE(LUOUT,*)"nmot=",nmot
nbf=nmot ! convert to INTEGER*8
orbstrt=0! construct a vector containing the orbital
! index at which each center starts. There
! could be an easier way to do this. MAKE SURE THIS IS
! COUNTING BASED ON SPHERICAL, NOT CARTESIAN!
orbstrt(1)=1
DO i = 1, nat-1
orbstrt(i+1)=orbstrt(i)
status = bas_ce2cnr(basis,i,contf,contl)
DO j = contf, contl
status = bas_cn2bfr(basis,j,shlf,shll)
orbstrt(i+1)=orbstrt(i+1)+shll-shlf+1
ENDDO
ENDDO
C LB
task_size = 1
status = rtdb_parallel(.true.) ! Broadcast reads to all processes
pointforce = geom_include_bqbq(geom)
dobq = geom_extbq_on()
call hf_print_set(1)
ijatom = -1
next = nxtask(nproc,task_size)
! Loop Through Bra atom
do 90, iat1 = 1, nat
! Loop through ket atom
do 80, iat2 = 1, iat1
ijatom = ijatom + 1
C LB
! WRITE(LUOUT,*)"LB Debug, iat1=",iat1,"iat2=",iat2
! & ,"ijatom=",ijatom
C LB
if ( ijatom .eq. next ) then
status = bas_ce2bfr(basis,iat1,iab1f,iab1l)
status = bas_ce2bfr(basis,iat2,iab2f,iab2l)
if (iab1f.le.0 .or. iab2f.le.0) then
c
c At least one center has no functions on it ... next atom
c
goto 1010
endif
if (oskel) then
if (.not. sym_atom_pair(geom, iat1, iat2, qfac))
$ goto 1010
else
qfac = 1.0d0
endif
status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
! call ga_get (g_dens, iab1f,iab1l,iab2f,iab2l,dens,max_at_bf)
! call ga_get(g_wdens,iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf)
! Loop through bra atom shells
do 70, ish1 = iac1f, iac1l
if ( iat1.eq.iat2 ) iac2l = ish1
! Loop through ket atom shells
do 60, ish2 = iac2f, iac2l
C LB
! WRITE(LUOUT,*)"LB Debug, ishl=",ish1,"ish2=",ish2
C LB
C shell block in atomic (D/Dw)-matrix block
status = bas_cn2bfr(basis,ish1,if1,il1)
C LB
! WRITE(LUOUT,*)"LB hdoverlap, if1=",if1,"il1=",il1
C LB
if1 = if1 - iab1f + 1
il1 = il1 - iab1f + 1
status = bas_cn2bfr(basis,ish2,if2,il2)
C LB
! WRITE(LUOUT,*)"LB hdoverlap, if2=",if2,"il2=",il2
C LB
if2 = if2 - iab2f + 1
il2 = il2 - iab2f + 1
nint = ( il1 - if1 + 1 ) * ( il2 - if2 + 1 )
C half-derivative overlap derivatives
C LB
hlfdrv=.true.
! WRITE(LUOUT,*)"LB hd, iat1=",iat1,"iat2=",iat2
! WRITE(LUOUT,*)"LB hd, ish1=",ish1,"ish2=",ish2
C LB
call intd_1eovQ(basis,ish1,basis,ish2,l_scr,scr,
& l_buf,H,idatom,hlfdrv)
C LB
! WRITE(LUOUT,*)"LB Debug, back from intd_1eovQ"
status=bas_continfo(basis,ish1,itype,dum1,dum2,dum3)
status=bas_continfo(basis,ish2,jtype,dum1,dum2,dum3)
nintsph=(2*itype+1)*(2*jtype+1)
! WRITE(LUOUT,*)"nintsph=",nintsph
! WRITE(LUOUT,*)"Ova="
! WRITE(LUOUT,*)H(1:nintsph*6)
C LB We now have H (BO2i) which includes the integrals organized by
C subshell and dx dy dz wrt both atom centers involved. We now
C need to get these into a matrix in the AO space, organized by
C atom center derivative. I am not going to do this by symmetry
C block right now since the initial SONADCT capability will only
C be available in C1 symmetry.
! WRITE(LUOUT,*)"LB hdoverlap, orbstrt=",orbstrt
! WRITE(LUOUT,*)"Bra Atom=",iat1,"shell=",ish1
! WRITE(LUOUT,*)"Ket Atom=",iat2,"shell=",ish2
atmidx1=iat1*3-2 !index of atom center derivative
atmidx2=iat2*3-2 !index of atom center derivative
! WRITE(LUOUT,*)"if1=",if1,"il1=",il1
! WRITE(LUOUT,*)"if2=",if2,"il2=",il2
! WRITE(LUOUT,*)"iac1f=",iac1f,"iac1l=",iac1l
! WRITE(LUOUT,*)"iac2f=",iac2f,"iac2l=",iac2l
j=0
IF (bas_spherical(basis)) THEN
il1=if1+2*itype
il2=if2+2*jtype
nint=nintsph
! WRITE(LUOUT,*)"in spherical region,
! & il1=",il1,"il2=",il2
ENDIF
DO i = if1, il1
j=j+1
nmss2=il2-if2+1
idx3=(j-1)*nmss2
IF (iat1.EQ.iat2) THEN
!basis functions from same atom
SQHDOL(atmidx1, orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& -H(idx3+1:idx3+nmss2+1 )!x1
SQHDOL(atmidx1+1,orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& -H(idx3+nint+1:idx3+nint+nmss2 )!y1
SQHDOL(atmidx1+2,orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& -H(idx3+2*nint+1:idx3+2*nint+nmss2)!z1
ELSE
!basis functions from different atoms
SQHDOL(atmidx1, orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& H(idx3+1:idx3+nmss2+1 )!x1
SQHDOL(atmidx1+1,orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& H(idx3+nint+1:idx3+nint+nmss2 )!y1
SQHDOL(atmidx1+2,orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& H(idx3+2*nint+1:idx3+2*nint+nmss2)!z1
SQHDOL(atmidx2 ,orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& H(idx3+1:idx3+nmss2+1)!x2
SQHDOL(atmidx2+1,orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& H(idx3+nint+1:idx3+nint+nmss2)!y2
SQHDOL(atmidx2+2,orbstrt(iat1)+i-1,
& orbstrt(iat2)+if2-1:orbstrt(iat2)+il2-1)=
& H(idx3+2*nint+1:idx3+2*nint+nmss2)!z2
ENDIF
ENDDO
! DO i=1,3*nat
! WRITE(LUOUT,*)"SQHDOL",i
! CALL output(SQHDOL(i,:,:),1,nmot,1,nmot,nmot,nmot,1)
! ENDDO
60 continue
70 continue
!
1010 continue
next = nxtask(nproc,task_size)
endif
!
80 continue
90 continue
next = nxtask(-nproc,task_size)
!copy temp SQHDOL matrix too hdol global array
dims(1) = 3*nat
dims(2) = nmot
dims(3) = nmot
CALL nga_put(hdol,istrt,dims,SQHDOL,dims(1:2))
cgk debug
! write(*,*)'gk: leaving hdoverlap'
cgk end
return
end

471
src/NWints/aoints/int_1e_sifs.F Normal file
View file

@ -0,0 +1,471 @@
subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft,
& nmap,map,imtype,ifmt1,ibvtyp,ibitv,l1rec,n1max,
& clab, ninfo, info)
implicit none
#include "errquit.fh"
#include "cint1cache.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "bas.fh"
#include "geom.fh"
#include "rtdb.fh"
#include "inp.fh"
#include "apiP.fh"
#include "cscfps.fh"
#include "sym.fh"
#include "cdft.fh"
c
c Compute the desired type of integrals (kinetic, potential, overlap)
c
c Oskel indicates that the skeleton (petite-list symmetry) matrix should be
c built ...
c
c arguments
c
integer ibas
character*(10) integ_type ! [input] Name of integrals to compute
* logical oskel ! [input] If true generate symmetry unique list
c
c local variables
c
integer nshell
integer ishell, jshell, mem1, max1e, lrecal, n1mx
integer ijshell, ilo, ihi, idim
integer l_buf, l_scr, l_ilab, l_jlab, l_info, l_bfnlab, l_sifbuf,
& l_sifval
integer k_buf, k_scr, k_ilab, k_jlab, k_info, k_bfnlab, k_sifbuf,
& k_sifval
integer type
logical odoit
double precision q2
external block_int1e ! For T3D
integer i, noffset,g_loc, j, ijmap
c
logical ocache_save
c
logical odbug, Lresult
logical osome
character*8 blabtmp
integer nbft
c
c sifs parameters
c
integer*4 aoints
integer*4 ifmt1, ibvtyp, n1max, ierr
INTEGER l1rec, ntitle, ibuf, nsym
integer*4 ibitv
integer*4 mxbf
parameter (mxbf=1000)
c header 1
integer*4 nbas, ninfo, mxenrgy, nenrgy, nmap
parameter(ntitle=1,nsym=1,mxenrgy=1)
c header 2
integer*4 info(ninfo) ! not to be confused with NWChem info
integer*4 ietype(mxenrgy)
integer*4 imtype(nmap) ! map will be taken from dynamic memory
integer*4 map(nbft,nmap)
integer*4 clab(2,*)
real*8 energy(nenrgy)
character*80 stitle(ntitle)
character*4 slabel(nsym)
integer*4 nbpsy(nsym)
character*1 shtypes(-1:7)
data shtypes/'l','s', 'p', 'd', 'f', 'g', 'h', 'i', 'k'/
integer*4 shmap(-1:7)
data shmap / 0, 1, 2, 4, 6, 9, 12, 16, 0 /
integer otype, nprim, ngen, sphcart, iat, shdim, igen, ibf
integer numtot
integer*4 msame,nmsame,nomore
parameter(msame=0, nmsame=1, nomore= 2)
integer*4 itypea, itypeb, nrec, num, last
real*8 fcore
logical basok
cgk provisional
* change thresh to the appropriate user supplied zero tolerance
real*8 thresh
parameter (thresh=1d-12)
cgk end
integer numints
c
cgk debug
! write(*,*)'gk: **** entered int_1e_sifs'
! WRITE(*,*)'LB: info=',info
cgk end
basok=.false.
odbug=.true.
osome=.true.
osome=osome.or.odbug
odbug=odbug.and.(ga_nodeid().eq.0)
osome=osome.and.(ga_nodeid().eq.0)
call ga_sync()
if (oscfps) call pstat_on(ps_int_1e)
c
c
c
c ----- save ocache logical variable -----
c
if(type.eq.6 .or. type.eq.7 ) then
ocache_save=ocache
ocache =.false.
endif
c
c Get info about the basis sets
c
if (.not. bas_numcont(ibas, nshell)) call errquit
$ ('int_1e_sifs: bas_numcont failed for ibas', ibas,
& BASIS_ERR)
if (nbft.gt.mxbf) call errquit
$ ('int_1e_sifs: nbft gt maximum aoints basis functions', nbft,
& BASIS_ERR)
c
c allocate necessary local temporary arrays on the stack
c
c k_* are the offsets corrsponding to the l_* handles
c
call int_mem_1e(max1e, mem1)
cgk debug
! write(*,*)'gk: in int_mem_1e, max1e=', max1e, ' mem1=',mem1
cgk end
c get memory for integral buffer
if (.not. MA_push_get(MT_DBL,max1e,'int_1e_sifs:buf',l_buf,k_buf))
$ call errquit('int_1e_sifs: ma failed', max1e, MA_ERR)
c get memory for scratch space
if (.not. MA_push_get(MT_DBL, mem1,'int_1e_sifs:scr',l_scr,k_scr))
$ call errquit('int_1e_sifs: ma failed', mem1, MA_ERR)
c get memory for bra labels
if (.not. MA_push_get(MT_Int,max1e,'int_1e_sifs:ilab',l_ilab,
& k_ilab))
$ call errquit('int_1e_sifs: ma failed for ilab', max1e,
& MA_ERR)
c get memory for ket labels
if (.not. MA_push_get(MT_Int,max1e,'int_1e_sifs:jlab',l_jlab,
& k_jlab))
$ call errquit('int_1e_sifs: ma failed for jlab', max1e,
& MA_ERR)
c get memory for basis info (character*16)
if (.not. ma_push_get(mt_byte,nbft*16,'int_1e_sifs:info',
$ l_info,k_info))
$ call errquit('int_1e_sifs: error getting info mem',
& nbft*16, MA_ERR)
c get memory for sifs bfnlab (character*8)
if (.not. ma_push_get(mt_byte,nbft*8,'int_1e_sifs:bfnlab',
$ l_bfnlab,k_bfnlab))
$ call errquit('int_1e_sifs: error getting bfnlab mem',
& nbft*8, MA_ERR)
c get memory for sifs buffer space
lrecal=l1rec
if (.not. MA_push_get(MT_DBL, lrecal,'int_1e_sifs:sifbuf',
& l_sifbuf, k_sifbuf))
$ call errquit('int_1e_sifs: ma failed getting sifbuf', lrecal,
& MA_ERR)
n1mx=n1max
if (.not. MA_push_get(MT_DBL, n1mx,'int_1e_sifs:sifval',
& l_sifval, k_sifval))
$ call errquit('int_1e_sifs: ma failed getting sifval', n1mx,
& MA_ERR)
c
c Set up SIFS header iformation
c
call bas_vec_info(ibas, byte_mb(k_info))
cgk debug
* write(*,*)'gk: basis set info'
* do i=0,nbft-1
* write(6,'(i4,2x,16a1)') i+1, (byte_mb(k_info+i*16+j), j=0,15)
* enddo
cgk end
c gk: instead of doing this, use 16 byte labels in SIFS/COLUMBUS.
c gk: consider this a temporary hack. 8 characters is not enough.
do i=0,nbft-1
do j=0,7
write(byte_mb(k_bfnlab+i*8+j),'(a1)')
& byte_mb(k_info+i*16+j)
enddo
enddo
cgk debug
! write(*,*)'gk: bfnlab:'
! do i=0,nbft-1
! write(6,'(i4,2x,8a1)') i+1, (byte_mb(k_bfnlab+i*8+j), j=0,7)
! enddo
cgk end
c put effective nuclear repulsion energy in energy(1)
ietype(1) = -1 ! nuclear repulsion=1 ... for now, but core later
LResult = geom_nuc_rep_energy(geom, energy(1))
imtype(1)=3 ! bfn-to-center map vector
imtype(2)=4 ! bfn-to-orbital_type map vector
fcore=0d0
cgk debug
! write(*,*)'gk: i, iat, shmap(otype)'
! write(*,*)'gk: nmap=',nmap
cgk end
i=0
do ishell=1,nshell
basok=bas_continfo(ibas, ishell, otype, nprim, ngen, sphcart)
basok=bas_cn2bfr(ibas, ishell, ilo, ihi)
basok=bas_cn2ce(ibas, ishell, iat)
if (otype.gt.6 .or. otype.lt.0) call errquit
$ ('int_1e_sifs: unsupported sifs basis otype?',otype,
& BASIS_ERR)
c
shdim = (ihi-ilo+1)
if (otype.ge.0) then
shdim = shdim / ngen
else ! Fudge for SP shells
ngen = 1
endif
if (otype .le. 1) sphcart = 0
do igen = 1, ngen
do ibf = 1, shdim
i = i + 1
if (i.gt.nbft) call errquit
$ ('int_1e_sifs: error counting basis functions',nbft,
& BASIS_ERR)
map(i,1)=iat
map(i,2)=shmap(otype)
cgk debug
* write(*,'(a,i,i,i)')'gk: ', i, iat, shmap(otype)
* write(*,'(a,i,i,i)')'gk: ', i, map(i,1), map(i,2)
cgk end
* int_mb(k_map+i-1) = iat
* int_mb(k_map+nbft+i-1) = shmap(otype)
enddo
enddo
enddo
slabel(1)='a1'
stitle(1)='AO integrals from NWChem'
nbpsy(1)=nbft
* call sifwh( aoints, ntitle, nsym, nbft,
* & ninfo, nenrgy, nmap, stitle, nbpsy, slabel, info,
* & bfnlab, ietype, energy, imtype, map, ierr )
ijshell = 0
q2 = 1.0d0
if (0 .eq. ga_nodeid()) then
cgk debug
! write(*,*)'gk: calling sifwh with aoints=',aoints
! write(*,*)'gk: nenrgy=',nenrgy
! write(*,*)'gk: ntitle=',ntitle
! write(*,*)'gk: ninfo=',ninfo
! WRITE(*,*)'LB: l2rec=',info(4)
cgk end
call sifwh( aoints, ntitle, nsym, nbft,
& ninfo, nenrgy, nmap, stitle, nbpsy, slabel, info,
& byte_mb(k_bfnlab), ietype, energy, imtype, map,
& ierr )
cgk debug
! write(*,*)'gk: back from sifwh with ierr=',ierr
cgk end
c start loop over integral types here.
do type=1,3
nrec=0
num=0
numtot=0
ibuf=0
last=msame
if (type.eq.1) then
integ_type='kinetic'
itypea=0
itypeb=1
else if (type.eq.2) then
integ_type='potential'
itypea=0
itypeb=2
else if (type.eq.3) then
integ_type='overlap'
itypea=0
itypeb=0
end if
cgk debug
! write(*,*)'gk: processing integral type: ', integ_type
cgk end
do jshell = 1, nshell
do ishell = 1, nshell
cgk debug
* write(*,*)'gk: ishell=',ishell,' jshell=',jshell
cgk end
odoit = .true.
if (oskel)
$ odoit = sym_shell_pair(ibas, ishell, jshell, q2)
if (odoit) then
if (.not. bas_cn2bfr(ibas, ishell, ilo, ihi))
$ call errquit('int_1e_ga: bas_cn2bfr ?', ibas,
& BASIS_ERR)
idim = ihi - ilo + 1
call ifill(max1e,0,int_mb(k_ilab),1)
call ifill(max1e,0,int_mb(k_jlab),1)
if (type .eq. 1) then
call int_l1eke (ibas, jshell, ibas, ishell,
& thresh, int_mb(k_ilab),int_mb(k_jlab), max1e,
$ dbl_mb(k_buf), mem1, dbl_mb(k_scr), numints)
else if (type .eq. 2) then
cgk debug
* write(*,*)'gk: calling int_l1epe'
* write(*,*)'gk: ibas = ', ibas
* write(*,*)'gk: ishell = ', ishell
* write(*,*)'gk: jshell = ', jshell
* write(*,*)'gk: thresh = ', thresh
* write(*,*)'gk: k_ilab = ', k_ilab
* write(*,*)'gk: k_jlab = ', k_jlab
* write(*,*)'gk: max1e = ', max1e
* write(*,*)'gk: k_buf = ', k_buf
* write(*,*)'gk: mem1 = ', mem1
* write(*,*)'gk: k_scr = ', k_scr
* write(*,*)'gk: numints = ', numints
cgk end
call int_l1epe (ibas, jshell, ibas, ishell,
& thresh, int_mb(k_ilab),int_mb(k_jlab), max1e,
$ dbl_mb(k_buf), mem1, dbl_mb(k_scr), numints)
cgk debug
* write(*,*)'gk: back from int_l1epe'
cgk end
else if (type .eq. 3) then
c ECP is summed here, but COLUMBUS should not care.
call int_l1eov (ibas, ishell, ibas, jshell, thresh,
& int_mb(k_ilab),int_mb(k_jlab), max1e,
& dbl_mb(k_buf), mem1, dbl_mb(k_scr), numints)
end if
c
cgk debug
* write(*,*)'gk: this batch, numints = ', numints
cgk end
do i=1,numints
if(ibuf.eq.n1max) then
numtot = numtot + ibuf
C LB
! WRITE(*,*)"int_1e_sifs calling sifew1 1,
! &last=",last,"ibvtyp=",ibvtyp
! WRITE(*,*)"sizeof(dbl_mb)=",sizeof(dbl_mb(k_sifval)),
! & "sizeof(buf)=",sizeof(dbl_mb(k_sifbuf))
C LB
call sifew1(aoints, info, 2, ibuf, last,
& itypea, itypeb, ifmt1, ibvtyp,
& dbl_mb(k_sifval), clab, fcore, ibitv,
& dbl_mb(k_sifbuf), nrec, ierr)
c ibuf on return has the number of unwritten
c integrals. dbl_mb(k_sifval+0:(ibuf-1))
c contains the unwritten values, ditto for labels
numtot = numtot - ibuf
endif
ibuf=ibuf+1
clab(1,ibuf)=int_mb(k_ilab+i-1)
clab(2,ibuf)=int_mb(k_jlab+i-1)
dbl_mb(k_sifval+ibuf-1)=dbl_mb(k_buf+i-1)
cgk debug
* write(*,*)'gk: ibuf=', ibuf
* write(*,'(a,2i4,f20.12)')
* & 'gk: ilab,jlab,dbl_mb:',
* & int_mb(k_ilab+i-1),int_mb(k_jlab+i-1),
* & dbl_mb(k_buf+i-1)
cgk end
enddo
end if
ijshell = ijshell + 1
end do
end do
last=nmsame
numtot=numtot+ibuf
cgk debug
* write(*,*)'gk: nrec=', nrec
* write(*,*)'gk: ibuf=', ibuf
* write(*,*)'gk: last=', last,"ibvtyp=",ibvtyp
* write(*,*)'gk: numtot=', numtot
* WRITE(*,*)"int_mom_sifs calling sifew1 2"
cgk end
call sifew1(aoints, info, 2, ibuf, last,
& itypea, itypeb, ifmt1, ibvtyp,
& dbl_mb(k_sifval), clab, fcore, ibitv,
& dbl_mb(k_sifbuf), nrec, ierr)
write(6,'(a,i10,1x,a,a)') 'Wrote ',numtot, integ_type,
& ' integrals to aoints'
enddo
endif
call ga_sync() ! So that no nasty races can result
cgk debug
! write(*,*)'gk: calling int_so_sifs from int_1e_sifs'
* write(*,*)'gk: mem1=', mem1
cgk end
call int_so_sifs(ibas, oskel, aoints, nbft, max1e, mem1, l1rec,
& n1max, dbl_mb(k_sifbuf), dbl_mb(k_sifval), ninfo, info, clab,
& fcore, ibvtyp, ibitv, ifmt1)
cgk debug
! write(*,*)'gk: calling int_mom_sifs from int_1e_sifs'
* write(*,*)'gk: mem1=', mem1
cgk end
call int_mom_sifs(ibas, oskel, aoints, nbft, l1rec,
& n1max, dbl_mb(k_sifbuf), dbl_mb(k_sifval), ninfo, info,
& fcore, ibvtyp, ibitv, ifmt1)
cgk debug
! write(*,*)'gk: back from int_so_sifs in int_1e_sifs'
cgk end
c
c chop stack at first item allocated
c
if (.not. MA_pop_stack(l_sifval)) call errquit
$ ('int_1e_sifs: pop failed at l_bfnlab', 0, GA_ERR)
if (.not. MA_pop_stack(l_sifbuf)) call errquit
$ ('int_1e_sifs: pop failed at l_bfnlab', 0, GA_ERR)
if (.not. MA_pop_stack(l_bfnlab)) call errquit
$ ('int_1e_sifs: pop failed at l_bfnlab', 0, GA_ERR)
if (.not. MA_pop_stack(l_info)) call errquit
$ ('int_1e_sifs: pop failed at l_info', 0, GA_ERR)
if (.not. MA_pop_stack(l_jlab)) call errquit
$ ('int_1e_sifs: pop failed at l_jlab', 0, GA_ERR)
if (.not. MA_pop_stack(l_ilab)) call errquit
$ ('int_1e_sifs: pop failed at l_ilab', 0, GA_ERR)
if (.not. MA_pop_stack(l_scr)) call errquit
$ ('int_1e_sifs: pop failed at l_scr', 0, GA_ERR)
if (.not. MA_pop_stack(l_buf)) call errquit
$ ('int_1e_sifs: pop failed at l_buf', 0, GA_ERR)
c
c ----- restore ocache -----
c
if(type.eq.6 .or. type.eq.7) then
ocache=ocache_save
endif
c
* if (.not.inp_compare(.false., integ_type, 'so')) then
* call ga_dadd(1.0d0, g_loc, 1.0d0, g, g)
* if (.not. ga_destroy(g_loc)) call errquit('i1ega: gad?',0,
* & GA_ERR)
* endif
call ga_sync() ! So that no nasty races can result
c
if (oscfps) call pstat_off(ps_int_1e)
c
cgk debug
! write(*,*)'gk: leaving int_1e_sifs at end of routine'
cgk end
end

221
src/NWints/aoints/int_2e_sifs.F Normal file
View file

@ -0,0 +1,221 @@
subroutine int_2e_sifs( geom, basis, tol2e, oskel, aoint2, ninfo,
& info, ifmt, ibvtyp, ibitv, clabs)
c
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "bas.fh"
#include "geom.fh"
#include "schwarz.fh"
#include "cscfps.fh"
#include "cfock.fh"
#include "util.fh"
c
c driver routine to write 2e integrals to sifs file based on
c ao_fock_2e
c
c arguments
c
integer geom, basis ! [input] parameter handles
integer*4 aoint2
integer*4 ninfo, ifmt, ibvtyp
integer*4 ibitv
integer*4 info(ninfo)
double precision tol2e ! [input] integral selection threshold
logical oskel ! [input] toggle skeleton Fock matrix
c
c local variables
c
integer natoms, ablklen, nshb, nbf, maxg, lenscr, nsh
integer i, mem_rep
integer l_dij, l_dik, l_dli, l_djk, l_dlj, l_dlk
integer l_fij, l_fik, l_fli, l_fjk, l_flj, l_flk
integer k_dij, k_dik, k_dli, k_djk, k_dlj, k_dlk
integer k_fij, k_fik, k_fli, k_fjk, k_flj, k_flk
integer l_atmp, k_atmp
integer l_block, k_block, nblock, maxblock
logical status
logical int2e_test_mem, int2e_buf_write, int2e_file_rewind
external int2e_test_mem, int2e_buf_write, int2e_file_rewind
integer ncomputed, nwritten, num
common/sif2e/ncomputed, nwritten, num
integer*4 max2e, l2rec, n2max, ierr
integer*4 clabs(4,*)
c
cgk debug
* write(*,*)'gk: entered int_2e_sifs'
cgk end
if (.not. bas_numbf(basis, nbf)) call errquit
$ ('ao_replicated: bas_numbf failed', 0, BASIS_ERR)
cgk debug
* write(*,*)'gk: nbf= ', nbf
cgk end
if (.not. bas_numcont(basis, nsh)) call errquit
$ ('ao_replicated: bas_numcont failed', 0, BASIS_ERR)
cgk debug
* write(*,*)'gk: nsh= ', nsh
cgk end
if (.not. geom_ncent(geom, natoms)) call errquit
$ ('ao_fock_2e: geom_ncent failed', 0, GEOM_ERR)
cgk debug
* write(*,*)'gk: natoms= ', natoms
cgk end
c
c Default is replicated data. This can be overridden by setting
c 'fock:replicated' as false. It is also overriden here when the
c replicated data algorithm would exceed memory requirements
c
* call intb_mem_2e4c(maxg, lenscr)
* mem_rep = nbf*nbf*2*nfock + nsh*nsh
* if (owritefile) mem_rep = mem_rep + max(0,memsize)
* if (owritefile .or. (.not. (owritefile .or. oreadfile)))
* $ mem_rep = mem_rep + maxeri*5 + lenscr
* if (oreplicated) then
* oreplicated = (mem_rep .lt. ma_inquire_avail(mt_dbl))
* endif
c
* if (oreplicated) then
* call ao_replicated( geom, basis, nfock, jfac, kfac,
* $ tol2e, oskel, vg_dens, vg_fock, asym )
* return
* endif
c
c allocate necessary local temporary arrays on the stack
c
c l_scr ... workspace for integral routines
c l_d** ... ** block of density matrix
c l_f** ... ** block of fock matrix
c
c k_* are the offsets corrsponding to the l_* handles
c
c If reading/writing to/from cache/file then rewind
c
* if (oreadfile .or. owritefile) then
* if (.not. int2e_file_rewind()) call errquit
* $ ('ao_fock_2e: failed rewinding integral file',0,
* & UNKNOWN_ERR)
* endif
c
c Halve exchange factor to conform to internal definition
c
* call dscal(nfock, 0.5d0, kfac, 1)
c
c Determine appropriate task chunking and max no. of bf in a
c block of the density/fock matrix
c
if (.not. geom_ncent(geom, natoms)) call errquit
$ ('ao_fock_2e: geom_ncent failed', 0, GEOM_ERR)
if (.not. ma_push_get(MT_INT, 2*natoms, 'fock2e:block',
$ l_block, k_block))call errquit('fock_2e: ma failed',2*natoms,
& MA_ERR)
cgk debug
* write(*,*)'gk: calling fock_2e_block_atoms'
cgk end
call fock_2e_block_atoms(basis, oskel, tol2e,
$ int_mb(k_block), nblock, maxblock)
cgk debug
* write(*,*)'gk: back from fock_2e_block_atoms'
cgk end
c
ablklen = maxblock**2
* maxd = nfock*ablklen
c
status = .true.
* status = status .and. ma_push_get(MT_DBL, maxd, 'dij',
* $ l_dij, k_dij)
* status = status .and. ma_push_get(MT_DBL, maxd, 'dik',
* $ l_dik, k_dik)
* status = status .and. ma_push_get(MT_DBL, maxd, 'dli',
* $ l_dli, k_dli)
* status = status .and. ma_push_get(MT_DBL, maxd, 'djk',
* $ l_djk, k_djk)
* status = status .and. ma_push_get(MT_DBL, maxd, 'dlj',
* $ l_dlj, k_dlj)
* status = status .and. ma_push_get(MT_DBL, maxd, 'dlk',
* $ l_dlk, k_dlk)
* status = status .and. ma_push_get(MT_DBL, maxd, 'fij',
* $ l_fij, k_fij)
* status = status .and. ma_push_get(MT_DBL, maxd, 'fik',
* $ l_fik, k_fik)
* status = status .and. ma_push_get(MT_DBL, maxd, 'fli',
* $ l_fli, k_fli)
* status = status .and. ma_push_get(MT_DBL, maxd, 'fjk',
* $ l_fjk, k_fjk)
* status = status .and. ma_push_get(MT_DBL, maxd, 'flj',
* $ l_flj, k_flj)
* status = status .and. ma_push_get(MT_DBL, maxd, 'flk',
* $ l_flk, k_flk)
status = ma_push_get(MT_DBL, ablklen, 'atmp',
$ l_atmp, k_atmp)
if (.not. status) call errquit('fock_2e: ablklen failed', ablklen,
& MA_ERR)
c
c If we can hold all of the integrals in core then we should
c do this ... i.e., the default should be semi-direct, though
c right now we must store all of the integrals.
c
cgk debug
* write(*,*)'gk: before int2e_test_mem oreadfile = ', oreadfile
* write(*,*)'gk: before int2e_test_mem owritefile = ', owritefile
cgk end
* if (.not. (oreadfile .or. owritefile))
* $ owritefile = int2e_test_mem(geom, basis, tol2e*0.01d0)
cgk debug
* write(*,*)'gk: after int2e_test_mem owritefile = ', owritefile
cgk end
c
* if (oreadfile) then
* call fock_2e_from_file( geom, basis, nfock, ablklen,
* $ jfac, kfac, tol2e, oskel,
* $ dbl_mb(k_dij), dbl_mb(k_dik), dbl_mb(k_dli),
* $ dbl_mb(k_djk), dbl_mb(k_dlj), dbl_mb(k_dlk),
* $ dbl_mb(k_fij), dbl_mb(k_fik), dbl_mb(k_fli),
* $ dbl_mb(k_fjk), dbl_mb(k_flj), dbl_mb(k_flk),
* $ dbl_mb(k_atmp), vg_dens, vg_fock)
* else
c
c tell the screening statistics we're starting another fock build
c
* iscreen(6) = iscreen(6) + 1
c
cgk debug
* write(*,*)'gk: calling int_2e_sifs_a'
cgk end
call int_2e_sifs_a( geom, basis, ablklen, tol2e, oskel,
$ dbl_mb(k_atmp), int_mb(k_block), nblock, aoint2, ninfo,
& info, ifmt, ibvtyp, ibitv, clabs)
cgk debug
* write(*,*)'gk: back from int_2e_sifs_a'
cgk end
c
c print the integral screening statistics for this fock build
c schwarz_print must be called by all nodes.
c
if (util_print('screening statistics', print_debug)) then
if ( .not. bas_numcont(basis, nshb) ) call errquit(
$ 'ao_fock_2e: problem with call to bas_numcont', basis,
& BASIS_ERR)
cgk debug
* write(*,*)'gk: calling schwarz_print'
cgk end
call schwarz_print(natoms, nshb)
cgk debug
* write(*,*)'gk: back from schwarz_print'
cgk end
endif
C
* endif
c
status = .true.
status = status .and. ma_pop_stack(l_atmp)
status = status .and. ma_pop_stack(l_block)
if (.not. status) call errquit('fock_2e: ma_pop?', 0, MA_ERR)
c
call ga_sync()
write(*,*) 'number of 2e integrals computed = ', ncomputed
write(*,*) 'number of 2e integrals written to aoints = ', nwritten
end

280
src/NWints/aoints/int_2e_sifs_a.F Normal file
View file

@ -0,0 +1,280 @@
subroutine int_2e_sifs_a( geom, basis, ablklen, tol2e, oskel,
$ tmp, blocks, nblock, aoint2, ninfo, info, ifmt, ibvtyp,
& ibitv, clabs)
implicit none
#include "errquit.fh"
c
#include "tcgmsg.fh"
#include "global.fh"
#include "schwarz.fh"
#include "geom.fh"
#include "cfock.fh"
#include "sym.fh"
#include "util.fh"
#include "cscfps.fh"
#include "nwc_const.fh"
#include "mafdecls.fh"
c
integer geom, basis
integer ablklen
double precision tol2e
logical oskel
* double precision dij(nfock*ablklen),dik(nfock*ablklen)
* double precision dli(nfock*ablklen),djk(nfock*ablklen)
* double precision dlj(nfock*ablklen),dlk(nfock*ablklen)
* double precision fij(nfock*ablklen),fik(nfock*ablklen)
* double precision fli(nfock*ablklen),fjk(nfock*ablklen)
* double precision flj(nfock*ablklen),flk(nfock*ablklen)
double precision tmp(ablklen)
* integer vg_dens(nfock)
* integer vg_fock(nfock)
integer blocks(2,*)
integer nblock
c
integer ncomputed, nwritten, num
common/sif2e/ncomputed, nwritten, num
c
integer natoms, nproc, ijkl, next
integer ib, jb, kb, lb, lbhi, i
integer iatlo, jatlo, katlo, latlo
integer iathi, jathi, kathi, lathi
integer ijk_prev(3,2) ! (i/j/k, lo/hi)
integer blklen
double precision sij, smax, skl ! schwarz info
double precision denmax, dentol, temp, integ_acc
logical otest
double precision schwarz_atom_block
external schwarz_atom_block
logical sym_atom_block_quartet
external sym_atom_block_quartet
c
integer nxtask
external nxtask
integer l_sifval, l_sifbuf, l2rec, n2max
integer k_sifval, k_sifbuf
integer*4 ninfo, ifmt, ibvtyp, ibitv, aoint2, nipv, itypea,
& itypeb, wait, nrec, reqnum, ierr, iwait
integer*4 info(ninfo)
integer*4 clabs(4,*)
integer*4 msame,nmsame,nomore
parameter(msame=0, nmsame=1, nomore= 2)
integer l_i, l_j, l_k, l_l, l_g
integer k_i, k_j, k_k, k_l, k_g
logical status
integer max2e, mem2
c
c Figure out the threshold on which to screen just integrals
c by dividing by maximum density matrix element. 10x for safety
c since really wanted max over atomic blocks
c
* denmax = 1d-10 ! To avoid a hard zero
* do i = 1, nfock
* call ga_maxelt(vg_dens(i), temp)
* denmax = max(temp, denmax)
* enddo
* dentol = min(dentolmax,tol2e/(10.0d0*denmax)) ! For screen integs only
c
c SEE AO_REPLICATED.F to make the next parameters consistent
c
c NOTE THAT ONLY USER SPECIFIED INTEGRAL ACCURACY IS NOW USED
c
* if (intacc .eq. 0.0d0) then
* integ_acc = min(1d-10,max(0.01d0*dentol,1d-20)) ! Variable
* call int_acc_set(integ_acc)
* else
cgk debug
* write(*,*)'gk: entered int_2e_sifs_a'
cgk end
integ_acc = intacc ! User controlled
cgk debug
* write(*,*)'gk: intacc = ', intacc
cgk end
* call int_acc_set(integ_acc)
integ_acc = 1d-10
cgk debug
* write(*,*)'gk: integ_acc = ', integ_acc
cgk end
* endif
c
c How many atoms are there?
c
ncomputed = 0
nwritten=0
nrec=0
if (.not. geom_ncent(geom, natoms)) call errquit
$ ('fock_2e: geom_ncent failed', 0, GEOM_ERR)
cgk debug
* write(*,*)'gk: natoms = ', natoms
cgk end
if (natoms .gt. nw_max_atom) call errquit
$ ('fock_2e: hard dimension failed', natoms, UNKNOWN_ERR)
cgk debug
* write(*,*)'gk: nw_max_atom = ', nw_max_atom
cgk end
c Get sifs memory here.
c values
n2max=info(5)
if (.not. MA_push_get(MT_DBL, n2max, 'int_2e_sifs_a:sifbuf',
& l_sifval, k_sifval))
$ call errquit('int_2e_sifs_a: ma failed getting sifbuf',
& n2max, MA_ERR)
c buffer
l2rec=info(4)
if (.not. MA_push_get(MT_DBL, l2rec, 'int_2e_sifs_a:sifbuf',
& l_sifbuf, k_sifbuf))
$ call errquit('int_2e_sifs_a: ma failed getting sifbuf',
& l2rec, MA_ERR)
call int_mem_2e4c(max2e, mem2)
mem2 = max(mem2,max2e)
status = .true.
status = status .and.
$ ma_push_get(MT_DBL, max2e, 'int_2e_sifs_b: buf', l_g, k_g)
status = status .and.
$ ma_push_get(MT_INT, max2e, 'int_2e_sifs_b: buf', l_i, k_i)
status = status .and.
$ ma_push_get(MT_INT, max2e, 'int_2e_sifs_b: buf', l_j, k_j)
status = status .and.
$ ma_push_get(MT_INT, max2e, 'int_2e_sifs_b: buf', l_k, k_k)
status = status .and.
$ ma_push_get(MT_INT, max2e, 'int_2e_sifs_b: buf', l_l, k_l)
if (.not. status) call errquit('fock_2e_b: not enuf memory',0,
& MA_ERR)
c
c How many processors are there?
c
nproc = ga_nnodes()
cgk debug
* write(*,*)'gk: nproc = ', nproc
cgk end
c
ijk_prev(1,1) = -1
ijk_prev(2,1) = -1
ijk_prev(3,1) = -1
ijk_prev(1,2) = -1
ijk_prev(2,2) = -1
ijk_prev(3,2) = -1
c
blklen = ablklen
c
ijkl = 0 ! this should force work to node zero
next = nxtask(nproc, 1)
c
c Loop thru blocked atomic quartets
c
smax = schwarz_max()
cgk debug
* write(*,*)'gk: starting loop through atomic quartets'
* write(*,*)'gk: smax = ', smax
* write(*,*)'gk: next = ', next
* write(*,*)'gk: nblock = ', nblock
* write(*,*)'gk: dentol = ', dentol
cgk end
do ib = nblock, 1, -1
iatlo = blocks(1,ib)
iathi = blocks(2,ib)
do jb = 1, ib
jatlo = blocks(1,jb)
jathi = blocks(2,jb)
* sij = schwarz_atom_block(iatlo,iathi,jatlo,jathi)
*
* Cannot use symmetry of pairs since the blocking of atoms means
* that we may have (kl) > (ij).
*
* if (sij*smax .ge. dentol) then
do kb = ib, 1, -1
katlo = blocks(1,kb)
kathi = blocks(2,kb)
lbhi = kb
if (ib .eq. kb) lbhi = jb
do lb = 1, lbhi
latlo = blocks(1,lb)
lathi = blocks(2,lb)
* skl = schwarz_atom_block(katlo,kathi,latlo,lathi)
* if (oskel) otest=sym_atom_block_quartet
* $ (geom, iatlo, iathi, jatlo, jathi,
* $ katlo, kathi, latlo, lathi)
otest=.true.
c
c Load balance over non-zero interactions
c
cgk but we need to need to send blocks of integrals to node zero
cgk to write out, so is it necessary to do this in parallel?
cgk debug
* write(*,*)'gk: otest = ', otest
* write(*,*)'gk: ijkl = ', ijkl
cgk end
if (otest .and. (ijkl .eq. next)) then
cgk debug
* write(*,*)'gk: calling sifs_2e_task'
cgk end
call sifs_2e_task(
$ geom, basis, oskel, iatlo, jatlo, katlo,
& latlo, iathi, jathi, kathi, lathi,
& ijk_prev, tmp, tol2e, aoint2, ninfo, info,
& ifmt, ibvtyp, dbl_mb(k_sifval), ibitv,
& dbl_mb(k_sifbuf), clabs, int_mb(k_i),
& int_mb(k_j), int_mb(k_k), int_mb(k_l),
& dbl_mb(k_g))
cgk debug
* write(*,*)'gk: back from sifs_2e_task'
cgk end
c
next = nxtask(nproc, 1)
end if
if (otest) ijkl = ijkl + 1
end do
enddo
* endif
end do
end do
nipv=4
itypea=3
itypeb=0
iwait=1
if (0 .eq. ga_nodeid()) then
cgk debug
* write(*,*)'gk: before sifew2'
* write(*,*)'gk: nwritten = ', nwritten
* write(*,*)'gk: num = ', num
cgk end
nwritten=nwritten+num
call sifew2( aoint2, info, nipv, num,
& nomore, itypea, itypeb, ifmt,
& ibvtyp, dbl_mb(k_sifval), clabs, ibitv,
& dbl_mb(k_sifbuf), iwait, nrec, reqnum,
& ierr )
cgk debug
* write(*,*)'gk: after sifew2'
* write(*,*)'gk: nwritten = ', nwritten
* write(*,*)'gk: num = ', num
cgk end
nwritten=nwritten+num
endif
next = nxtask(-nproc, 1)
call ga_sync()
status = status .and. ma_pop_stack(l_l)
status = status .and. ma_pop_stack(l_k)
status = status .and. ma_pop_stack(l_j)
status = status .and. ma_pop_stack(l_i)
status = status .and. ma_pop_stack(l_g)
if (.not. status) call errquit
$ ('int_2e_sifs_b: pop stack failed?', 0, MA_ERR)
if (.not. MA_pop_stack(l_sifbuf)) call errquit
$ ('int_2e_sifs_a: pop failed at l_bfnlab', 0, GA_ERR)
if (.not. MA_pop_stack(l_sifval)) call errquit
$ ('int_2e_sifs_a: pop failed at l_bfnlab', 0, GA_ERR)
cgk debug
* write(*,*)'gk: leaving int_2e_sifs_a'
cgk end
c
end

190
src/NWints/aoints/int_2e_sifs_b.F Normal file
View file

@ -0,0 +1,190 @@
subroutine int_2e_sifs_b( basis, tol2e, q4, iat, jat, kat, lat,
$ ilo, jlo, klo, llo, ihi, jhi, khi, lhi, aoint2, ninfo,
$ info, ifmt, ibvtyp, values, ibitv, buffer, clabs,
& ilabs, jlabs, klabs, llabs, eris )
c
implicit none
#include "errquit.fh"
#include "mafdecls.fh"
#include "schwarz.fh"
#include "bas.fh"
#include "cfock.fh"
#include "util.fh"
#include "cscfps.fh"
#include "stdio.fh"
integer ilabs(*), jlabs(*), klabs(*), llabs(*)
double precision eris(*)
integer basis
integer nfock
double precision tol2e, dentol
double precision sijkl
double precision q4
integer iat, jat, kat, lat
integer ilo, jlo, klo, llo
integer ihi, jhi, khi, lhi
integer ncomputed, nwritten, num
common/sif2e/ncomputed, nwritten, num
logical oij, okl, oikjl, status
integer ish, jsh, ksh, lsh, mem2, max2e
integer ishlo, ishhi, jshlo, jshhi, kshlo, kshhi, lshlo, lshhi
integer ibflo, jbflo, kbflo, lbflo, ibfhi, jbfhi, kbfhi, lbfhi
integer ishbflo, ishbfhi, jshbflo, jshbfhi
integer kshbflo, kshbfhi, lshbflo, lshbfhi
integer jshtop, kshtop, lshtop
integer l_scr, k_scr
double precision sij, smax, denmax
integer neri
integer i
integer*4 aoint2, ninfo, ifmt, ibvtyp, last, itypea, itypeb,
& iwait, nrec, reqnum, ierr, ibitv, nipv
parameter(nipv=4)
integer*4 clabs(nipv,*)
integer*4 info(ninfo)
double precision values(*), buffer(*)
integer*4 msame,nmsame,nomore
parameter(msame=0, nmsame=1, nomore= 2)
integer n2max
c iat, jat, kat, lat externally loop in canonical order over
c atoms.
c
c scratch space needed for 2-electron routines
c
c l_scr ... integral scratch space
c l_g ... buffer for eris
c l_i, j, k, l ... integral labels
c
c Shell and function ranges for each atom
c
cgk debug
* write(*,*)'gk: entered int_2e_sifs_b'
cgk end
* oerilabel=.true.
status = .true.
status = status .and. bas_ce2cnr(basis, iat, ishlo, ishhi)
status = status .and. bas_ce2cnr(basis, jat, jshlo, jshhi)
status = status .and. bas_ce2cnr(basis, kat, kshlo, kshhi)
status = status .and. bas_ce2cnr(basis, lat, lshlo, lshhi)
if (.not. owritefile) then
status = status .and. bas_ce2bfr(basis, iat, ibflo, ibfhi)
status = status .and. bas_ce2bfr(basis, jat, jbflo, jbfhi)
status = status .and. bas_ce2bfr(basis, kat, kbflo, kbfhi)
status = status .and. bas_ce2bfr(basis, lat, lbflo, lbfhi)
endif
if (.not. status) call errquit
$ ('int_2e_sifs_b: bas_ce to cnr/bfr failed', 0, BASIS_ERR)
call int_mem_2e4c(max2e, mem2)
mem2 = max(mem2,max2e)
status = .true.
status = status .and.
& ma_push_get(MT_DBL, mem2, 'int_2e_sifs_b: scr', l_scr, k_scr)
if (.not. status) call errquit('int_2e_sifs_b: not enuf memory',0,
& MA_ERR)
c
c Permutation symmetry boolean
c
oij = iat .eq. jat
okl = kat .eq. lat
oikjl = iat.eq.kat .and. jat.eq.lat
c
c loop over shells on each atom
c
smax = schwarz_max()
n2max=info(5)
iwait=1
itypea=3
itypeb=0
cgk debug
* write(*,*)'gk: at top of loop over shell quartets'
cgk end
do ish = ishlo, ishhi
if (.not. bas_cn2bfr(basis, ish, ishbflo, ishbfhi))
$ call errquit('int_2e_sifs_b: bas_cn2bfr', ish, BASIS_ERR)
jshtop = jshhi
if (oij) jshtop = ish
do jsh = jshlo, jshtop
sij = schwarz_shell(ish,jsh)
* if (sij*smax.ge.tol2e) then
if (.not. bas_cn2bfr(basis, jsh, jshbflo, jshbfhi))
$ call errquit('int_2e_sifs_b: bas_cn2bfr', jsh,
& BASIS_ERR)
kshtop = kshhi
if (oikjl) kshtop = ish
do ksh = kshlo, kshtop
if (.not. bas_cn2bfr(basis, ksh, kshbflo, kshbfhi))
$ call errquit('int_2e_sifs_b: bas_cn2bfr', ksh,
& BASIS_ERR)
lshtop = lshhi
if (okl) lshtop = ksh
if (oikjl .and. ksh.eq.ish) lshtop = jsh
do lsh = lshlo, lshtop
* if (sij*schwarz_shell(ksh,lsh)
* $ .ge. tol2e) then
if (.not. bas_cn2bfr(basis, lsh,
$ lshbflo, lshbfhi)) call errquit
$ ('int_2e_sifs_b: bas_cn2bfr', lsh,
& BASIS_ERR)
c
c Compute the shell block of integrals with labels
c
cgk debug
* write(*,*)'gk: calling int_2e4c'
cgk end
call int_2e4c(basis, ish, jsh, basis, ksh, lsh,
$ mem2, dbl_mb(k_scr), max2e, eris)
cgk debug
* write(*,*)'gk: oerilable = ', oerilabel
cgk end
call fock_2e_shell_to_labels( tol2e,
$ ishbflo, ishbfhi, jshbflo, jshbfhi,
$ kshbflo, kshbfhi, lshbflo, lshbfhi,
& ilabs, jlabs, klabs, llabs, eris,
$ dbl_mb(k_scr), neri)
ncomputed = ncomputed + neri
cgk debug
* write(*,*)'neri = ', neri
do i=1,neri
cgk debug
* write(*,2223) ilabs(i), jlabs(i), klabs(i), llabs(i), eris(i)
2223 format(4i4,f20.12)
cgk end
if(num.eq.n2max) then
nwritten=nwritten+num
call sifew2( aoint2, info, nipv, num,
& msame, itypea, itypeb, ifmt,
& ibvtyp, values, clabs, ibitv,
& buffer, iwait, nrec, reqnum,
& ierr )
nwritten=nwritten-num
endif
num=num+1
clabs(1,num)=ilabs(i)
clabs(2,num)=jlabs(i)
clabs(3,num)=klabs(i)
clabs(4,num)=llabs(i)
values(num)=eris(i)
enddo
* end if
end do
end do
* end if
end do
end do
c
status = ma_pop_stack(l_scr)
if (.not. status) call errquit
$ ('int_2e_sifs_b: pop stack failed?', 0, MA_ERR)
cgk debug
* write(*,*)'gk: nwritten = ', nwritten
* write(*,*)'gk: num = ', num
* write(*,*)'gk: leaving int_2e_sifs_b'
cgk end
c
end

367
src/NWints/aoints/int_mom_sifs.F Normal file
View file

@ -0,0 +1,367 @@
subroutine int_mom_sifs(ibas, oskel, aoints, nbft,
& l1rec, n1max,
& sifbuf, sifval, ninfo, info, fcore, ibvtyp, ibitv,
& ifmt1)
implicit none
#include "errquit.fh"
#include "cint1cache.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "rtdb.fh"
#include "inp.fh"
#include "apiP.fh"
#include "bas.fh"
#include "cscfps.fh"
#include "sym.fh"
c
c Oskel indicates that the skeleton (petite-list symmetry) matrix should be
c built ...
c
c arguments
c
integer ibas
logical oskel ! [input] If true generate symmetry unique list
double precision sifval(n1max)
double precision sifbuf(l1rec)
integer*4 ninfo
integer*4 info(ninfo) ! not to be confused with NWChem info
integer*4 clab(2,n1max)
real*8 fcore
c
c local variables
c
integer nshell
integer ishell, jshell, mem1, max1e, lrecal, n1mx, maxint
integer ilo, ihi, jlo, jhi, max1mom, mem1mom
integer l_buf, l_scr,
& l_momx, l_momy, l_momz, l_imom, l_jmom
integer k_buf, k_scr,
& k_momx, k_momy, k_momz, k_imom, k_jmom
integer noffsetx, noffsety, noffsetz
logical odoit
double precision q2
external block_int1e ! For T3D
integer i, noffset,g_loc, j, ijmap, symmap, symints
c
logical ocache_save
c
logical odbug
logical osome
character*8 blabtmp
integer nbft
c
c sifs parameters
c
integer*4 aoints
integer*4 ifmt1, l1rec, n1max, ierr
integer*4 ibvtyp
integer*4 ibitv
integer*4 mxbf
c make sure this is compatible with columbus
parameter (mxbf=1000)
c header 1
integer*4 ntitle, nsym, nbas, mxenrgy
parameter(ntitle=1,nsym=1,mxenrgy=1)
c header 2
integer otype, nprim, ngen, sphcart, iat, shdim, igen, ibf
integer numtot
integer*4 msame,nmsame,nomore
parameter(msame=0, nmsame=1, nomore= 2)
integer*4 nrec, num
integer itypea, itypeb, last
integer*4 ibuf
logical basok
real*8 momval
cgk provisional
* change thresh to the appropriate user supplied zero tolerance
real*8 thresh
parameter (thresh=1d-12)
cgk end
integer numints, icart, kcart, int_type
c
cgk debug
* write(*,*)'gk: **** entered int_mom_sifs'
cgk end
basok=.false.
odbug=.true.
osome=.true.
osome=osome.or.odbug
odbug=odbug.and.(ga_nodeid().eq.0)
osome=osome.and.(ga_nodeid().eq.0)
c
call ga_sync()
if (oscfps) call pstat_on(ps_int_1e)
c
c Get info about the basis sets
c
if (.not. bas_numcont(ibas, nshell)) call errquit
$ ('int_mom_sifs: bas_numcont failed for ibas', ibas,
& BASIS_ERR)
if (nbft.gt.mxbf) call errquit
$ ('int_mom_sifs: nbft gt maximum aoints basis functions',
& nbft, BASIS_ERR)
c
c allocate necessary local temporary arrays on the stack
c
c l_buf ... buffer to hold shell block of matrix
c l_s ... buffer to hold shell block of matrix
c l_scr ... workspace for integral routines
c
c k_* are the offsets corrsponding to the l_* handles
c
cgk debug
* write(*,*)'gk: max1e=', max1e
* write(*,*)'gk: mem1=', mem1
* write(*,*)'gk: nbft=', nbft
cgk end
call int_init_dip(max1e,mem1,ibas)
max1mom=3*max1e
mem1mom=3*mem1
mem1mom=max(mem1mom,max1mom)
maxint=nbft*nbft
cgk debug
* write(*,*)'gk: max1mom=', max1mom
* write(*,*)'gk: mem1mom=', mem1mom
* write(*,*)'gk: maxint=', maxint
cgk end
c get memory for integral buffer
if (.not. MA_push_get(MT_DBL,max1mom,'int_mom_sifs:buf',l_buf,
& k_buf))
$ call errquit('int_mom_sifs: ma failed', max1mom, MA_ERR)
c get memory for scratch space
if (.not. MA_push_get(MT_DBL, mem1mom,'int_mom_sifs:scr',l_scr,
& k_scr))
$ call errquit('int_mom_sifs: ma failed', mem1mom, MA_ERR)
c get memory for x
if (.not. MA_push_get(MT_DBL,maxint,'int_mom_sifs:momx',l_momx,
& k_momx))
$ call errquit('int_mom_sifs: ma failed for mom x', maxint,
& MA_ERR)
c get memory for y
if (.not. MA_push_get(MT_DBL,maxint,'int_mom_sifs:momy',l_momy,
& k_momy))
$ call errquit('int_mom_sifs: ma failed for mom y', maxint,
& MA_ERR)
c get memory for z
if (.not. MA_push_get(MT_DBL,maxint,'int_mom_sifs:momz',l_momz,
& k_momz))
$ call errquit('int_mom_sifs: ma failed for mom z', maxint,
& MA_ERR)
c get memory for full bra label array
if (.not. MA_push_get(MT_Int,maxint,'int_mom_sifs:imom',l_imom,
& k_imom))
$ call errquit('int_mom_sifs: ma failed for imom', maxint,
& MA_ERR)
c get memory for full ket label array
if (.not. MA_push_get(MT_Int,maxint,'int_mom_sifs:jmom',l_jmom,
& k_jmom))
$ call errquit('int_mom_sifs: ma failed for imom', maxint,
& MA_ERR)
c start loop over integral types here
c int_type = 0, spin-orbit potential
c int_type = 1, linear momentum
c int_type = 2, angular momentum
do int_type=1,2
if(int_type.eq.1.or.int_type.eq.2) itypea=2 ! gonna do dipole
if (0 .eq. ga_nodeid()) then
cgk debug
* write(*,*)'gk: int_type = ', int_type
cgk end
c start loop over integrals here.
numints=0
do jshell = 1, nshell
do ishell = 1, nshell
cgk debug
* write(*,*)'gk: ishell=',ishell,' jshell=',jshell
cgk end
odoit = .true.
if (oskel)
$ odoit = sym_shell_pair(ibas, ishell, jshell, q2)
if (odoit) then
if (.not. bas_cn2bfr(ibas, ishell, ilo, ihi))
$ call errquit('int_1e_ga: bas_cn2bfr ishell',
& ishell, BASIS_ERR)
if (.not. bas_cn2bfr(ibas, jshell, jlo, jhi))
$ call errquit('int_1e_ga: bas_cn2bfr jshell',
& jshell, BASIS_ERR)
if(int_type.eq.1) then
cgk debug
* write(*,*)'gk: calling int_veloc'
cgk end
C LB
! call int_veloc (ibas, jshell, ibas, ishell,
! $ mem1mom, dbl_mb(k_scr), max1mom, dbl_mb(k_buf))
call int_veloc (ibas, ishell, ibas, jshell,
$ mem1mom, dbl_mb(k_scr), max1mom, dbl_mb(k_buf))
elseif(int_type.eq.2) then
C LB
cgk debug
* write(*,*)'gk: calling int_angmom'
cgk end
C LB
! call int_angmom (ibas, jshell, ibas, ishell,
! $ mem1mom, dbl_mb(k_scr), max1mom, dbl_mb(k_buf))
call int_angmom (ibas, ishell, ibas, jshell,
$ mem1mom, dbl_mb(k_scr), max1mom, dbl_mb(k_buf))
C LB
endif
noffsetz = 0
noffsety = (ihi-ilo+1)*(jhi-jlo+1)
noffsetx = (ihi-ilo+1)*(jhi-jlo+1)*2
do j=jlo,jhi
do i=ilo,ihi
cgk debug
* write(*,'(a,i4,a,i4,a,i4,a,i4)')
* & 'gk: ilo=',ilo,' ihi=',ihi,' jlo=',jlo,' jhi=',jhi
* write(*,*)'gk: numints=',numints
cgk end
ijmap=(j-jlo)*(ihi-ilo+1)+(i-ilo)
if(i.ge.j) then
symmap=nbft*(j-1)-((j-1)*j)/2+i
* write(*,*)'gk: symmap=',symmap
dbl_mb(k_momx+symmap-1)=
& dbl_mb(k_buf+noffsetx+ijmap)
dbl_mb(k_momy+symmap-1)=
& dbl_mb(k_buf+noffsety+ijmap)
dbl_mb(k_momz+symmap-1)=
& dbl_mb(k_buf+noffsetz+ijmap)
int_mb(k_imom+symmap-1)=i
int_mb(k_jmom+symmap-1)=j
cgk debug
* write(*,'(a,3i3,f18.12)')
* & 'gk: i,j,ij,dbl_mb(mom_x):',
* & i,j,numints+ijmap,dbl_mb(k_momx+symmap-1)
* write(*,'(a,3i3,f18.12)')
* & 'gk: i,j,ij,dbl_mb(mom_y):',
* & i,j,numints+ijmap,dbl_mb(k_momy+symmap-1)
* write(*,'(a,3i3,f18.12)')
* & 'gk: i,j,ij,dbl_mb(mom_z):',
* & i,j,numints+ijmap,dbl_mb(k_momz+symmap-1)
cgk end
endif
enddo
enddo
numints=numints+ijmap+1
end if
end do
end do
c call sifs
endif
call ga_sync() ! So that no nasty races can result
cgk debug
* write(*,*)'gk: numints = ', numints
* write(*,*)'gk: maxint = ', maxint
* symmap=0
* do i=1,nbft
* do j = 1, i
* write(*,'(a,2i8,3f18.12)')
* & 'gk: i,j,so_x,so_y,so_z:', int_mb(k_imom+symmap),
* & int_mb(k_jmom+symmap),dbl_mb(k_momx+symmap),
* & dbl_mb(k_momy+symmap), dbl_mb(k_momz+symmap)
* symmap=symmap+1
* enddo
* enddo
do icart = 1, 3
if(icart.eq.1) then
kcart=k_momz
itypeb=0+int_type*3
elseif(icart.eq.2) then
kcart=k_momy
itypeb=1+int_type*3
elseif(icart.eq.3) then
kcart=k_momx
itypeb=2+int_type*3
endif
symmap=0
numtot=0
ibuf=0
do i=1,nbft
do j = 1, i
if(ibuf.eq.n1max) then
numtot = numtot + ibuf
C LB
* WRITE(*,*)"LB: in int_mom_sifs"
* WRITE(*,*)"LB: calling sifew1 1 last=",last,"ibvtyp=",ibvtyp
C LB
call sifew1(aoints, info, 2, ibuf, msame,
& itypea, itypeb, ifmt1, ibvtyp,
& sifval, clab, fcore, ibitv,
& sifbuf, nrec, ierr)
numtot = numtot - ibuf
endif
momval=dbl_mb(kcart+symmap)
if(abs(momval).gt.thresh) then
ibuf=ibuf+1
clab(1,ibuf)=int_mb(k_imom+symmap)
clab(2,ibuf)=int_mb(k_jmom+symmap)
sifval(ibuf)=momval
endif
symmap=symmap+1
enddo
enddo
if (icart.eq.3.and.int_type.eq.2) then
last=nomore
else
last=nmsame
endif
numtot=numtot+ibuf
cgk debug
* write(*,*)'gk: nrec=', nrec
* write(*,*)'gk: ibuf=', ibuf
* write(*,*)'gk: numtot=', numtot
cgk end
C LB
! WRITE(*,*)"LB: in int_mom_sifs"
! WRITE(*,*)"LB: calling sifew1 2 last=",last,"ibvtyp=",ibvtyp
C LB
call sifew1(aoints, info, 2, ibuf, last,
& itypea, itypeb, ifmt1, ibvtyp,
& sifval, clab, fcore, ibitv,
& sifbuf, nrec, ierr)
enddo
enddo
cgk end
cgk debug
* write(*,*)'gk: l_jmom = ', l_jmom
* write(*,*)'gk: k_jmom = ', k_jmom
cgk end
c
c chop stack at first item allocated
c
if (.not. MA_pop_stack(l_jmom)) call errquit
$ ('int_mom_sifs: pop failed at l_jmom', 0, GA_ERR)
if (.not. MA_pop_stack(l_imom)) call errquit
$ ('int_mom_sifs: pop failed at l_imom', 0, GA_ERR)
if (.not. MA_pop_stack(l_momz)) call errquit
$ ('int_mom_sifs: pop failed at l_momz', 0, GA_ERR)
if (.not. MA_pop_stack(l_momy)) call errquit
$ ('int_mom_sifs: pop failed at l_momy', 0, GA_ERR)
if (.not. MA_pop_stack(l_momx)) call errquit
$ ('int_mom_sifs: pop failed at l_momx', 0, GA_ERR)
if (.not. MA_pop_stack(l_scr)) call errquit
$ ('int_mom_sifs: pop failed at l_scr', 0, GA_ERR)
if (.not. MA_pop_stack(l_buf)) call errquit
$ ('int_mom_sifs: pop failed at l_buf', 0, GA_ERR)
call ga_sync() ! So that no nasty races can result
if (oscfps) call pstat_off(ps_int_1e)
c
cgk debug
* write(*,*)'gk: leaving int_so_sifs at end of routine'
cgk end
end

356
src/NWints/aoints/int_so_sifs.F Normal file
View file

@ -0,0 +1,356 @@
subroutine int_so_sifs(ibas, oskel, aoints, nbft, max1e, mem1,
& l1rec, n1max,
& sifbuf, sifval, ninfo, info, clab, fcore, ibvtyp, ibitv,
& ifmt1)
implicit none
#include "errquit.fh"
#include "cint1cache.fh"
#include "mafdecls.fh"
#include "global.fh"
#include "rtdb.fh"
#include "inp.fh"
#include "apiP.fh"
#include "bas.fh"
#include "cscfps.fh"
#include "sym.fh"
c
c Oskel indicates that the skeleton (petite-list symmetry) matrix should be
c built ...
c
c arguments
c
integer ibas
logical oskel ! [input] If true generate symmetry unique list
double precision sifval(n1max)
double precision sifbuf(l1rec)
integer*4 ninfo
integer*4 info(ninfo) ! not to be confused with NWChem info
integer*4 clab(2,max1e)
real*8 fcore
c
c local variables
c
integer nshell
integer ishell, jshell, mem1, max1e, lrecal, n1mx, maxint
integer ilo, ihi, jlo, jhi, max1eso, mem1so
integer l_buf, l_scr,
& l_sox, l_soy, l_soz, l_iso, l_jso
integer k_buf, k_scr, k_ilab, k_jlab,
& k_sox, k_soy, k_soz, k_iso, k_jso
integer type, noffsetx, noffsety, noffsetz
logical odoit
double precision q2
external block_int1e ! For T3D
integer i, noffset,g_loc, j, ijmap, symmap, symints
c
logical ocache_save
c
logical odbug
logical osome
character*8 blabtmp
integer nbft
c
c sifs parameters
c
integer*4 aoints
integer*4 ifmt1, l1rec, n1max, ierr
integer*4 ibvtyp
integer*4 ibitv
integer*4 mxbf
c make sure this is compatible with columbus
parameter (mxbf=1000)
c header 1
integer*4 ntitle, nsym, nbas, mxenrgy
parameter(ntitle=1,nsym=1,mxenrgy=1)
c header 2
integer otype, nprim, ngen, sphcart, iat, shdim, igen, ibf
integer numtot
integer*4 msame,nmsame,nomore
parameter(msame=0, nmsame=1, nomore= 2)
integer*4 nrec, num
integer itypea, itypeb, last
integer*4 ibuf
logical basok
* character*3 sotype(3)
* data sotype /'sox','soy','soz'/
real*8 soval
cgk provisional
* change thresh to the appropriate user supplied zero tolerance
real*8 thresh
parameter (thresh=1d-12)
cgk end
integer numints, icart, kcart
c
cgk debug
* write(*,*)'gk: **** entered int_so_sifs'
cgk end
basok=.false.
odbug=.true.
osome=.true.
osome=osome.or.odbug
odbug=odbug.and.(ga_nodeid().eq.0)
osome=osome.and.(ga_nodeid().eq.0)
c
call ga_sync()
if (oscfps) call pstat_on(ps_int_1e)
c
c
c
c ----- save ocache logical variable -----
c
if(type.eq.6 .or. type.eq.7 ) then
ocache_save=ocache
ocache =.false.
endif
c
c Get info about the basis sets
c
if (.not. bas_numcont(ibas, nshell)) call errquit
$ ('int_so_sifs: bas_numcont failed for ibas', ibas,
& BASIS_ERR)
if (nbft.gt.mxbf) call errquit
$ ('int_so_sifs: nbft gt maximum aoints basis functions', nbft,
& BASIS_ERR)
c
c allocate necessary local temporary arrays on the stack
c
c l_buf ... buffer to hold shell block of matrix
c l_s ... buffer to hold shell block of matrix
c l_scr ... workspace for integral routines
c
c k_* are the offsets corrsponding to the l_* handles
c
cgk debug
* write(*,*)'gk: max1e=', max1e
* write(*,*)'gk: mem1=', mem1
* write(*,*)'gk: nbft=', nbft
cgk end
max1eso=3*max1e
mem1so=3*mem1
maxint=nbft*(nbft+1)/2
cgk debug
* write(*,*)'gk: max1eso=', max1eso
* write(*,*)'gk: mem1so=', mem1so
* write(*,*)'gk: maxint=', maxint
* WRITE(*,*)"LB: start MA_push_get"
cgk end
C LB
* WRITE(*,*)"LB: integral buffer"
C LB
c get memory for integral buffer
if (.not. MA_push_get(MT_DBL,max1eso,'int_so_sifs:buf',l_buf,
& k_buf))
$ call errquit('int_so_sifs: ma failed', max1eso, MA_ERR)
c get memory for scratch space
if (.not. MA_push_get(MT_DBL, mem1so,'int_so_sifs:scr',l_scr,
& k_scr))
$ call errquit('int_so_sifs: ma failed', mem1so, MA_ERR)
c get memory for x spin orbit
if (.not. MA_push_get(MT_DBL,maxint,'int_so_sifs:sox',l_sox,
& k_sox))
$ call errquit('int_so_sifs: ma failed for so x', maxint,
& MA_ERR)
c get memory for y spin orbit
if (.not. MA_push_get(MT_DBL,maxint,'int_so_sifs:soy',l_soy,
& k_soy))
$ call errquit('int_so_sifs: ma failed for so y', maxint,
& MA_ERR)
c get memory for z spin orbit
if (.not. MA_push_get(MT_DBL,maxint,'int_so_sifs:soz',l_soz,
& k_soz))
$ call errquit('int_so_sifs: ma failed for so z', maxint,
& MA_ERR)
c get memory for full bra label array
if (.not. MA_push_get(MT_Int,maxint,'int_so_sifs:iso',l_iso,
& k_iso))
$ call errquit('int_so_sifs: ma failed for iso', maxint,
& MA_ERR)
c get memory for full ket label array
if (.not. MA_push_get(MT_Int,maxint,'int_so_sifs:jso',l_jso,
& k_jso))
$ call errquit('int_so_sifs: ma failed for iso', maxint,
& MA_ERR)
if (0 .eq. ga_nodeid()) then
c start loop over so integrals here.
numints=0
do jshell = 1, nshell
do ishell = 1, nshell
cgk debug
* write(*,*)'gk: ishell=',ishell,' jshell=',jshell
cgk end
odoit = .true.
if (oskel)
$ odoit = sym_shell_pair(ibas, ishell, jshell, q2)
if (odoit) then
if (.not. bas_cn2bfr(ibas, ishell, ilo, ihi))
$ call errquit('int_1e_ga: bas_cn2bfr ishell',
& ishell, BASIS_ERR)
if (.not. bas_cn2bfr(ibas, jshell, jlo, jhi))
$ call errquit('int_1e_ga: bas_cn2bfr jshell',
& jshell, BASIS_ERR)
call intso_1e (ibas, jshell, ibas, ishell,
$ mem1, dbl_mb(k_scr), max1eso, dbl_mb(k_buf))
noffsetz = 0
noffsety = (ihi-ilo+1)*(jhi-jlo+1)
noffsetx = (ihi-ilo+1)*(jhi-jlo+1)*2
do j=jlo,jhi
do i=ilo,ihi
cgk debug
* write(*,'(a,i4,a,i4,a,i4,a,i4)')
* & 'gk: ilo=',ilo,' ihi=',ihi,' jlo=',jlo,' jhi=',jhi
* write(*,*)'gk: numints=',numints
cgk end
ijmap=(j-jlo)*(ihi-ilo+1)+(i-ilo)
if(i.ge.j) then
symmap=nbft*(j-1)-((j-1)*j)/2+i
* write(*,*)'gk: symmap=',symmap
dbl_mb(k_sox+symmap-1)=
& dbl_mb(k_buf+noffsetx+ijmap)
dbl_mb(k_soy+symmap-1)=
& dbl_mb(k_buf+noffsety+ijmap)
dbl_mb(k_soz+symmap-1)=
& dbl_mb(k_buf+noffsetz+ijmap)
int_mb(k_iso+symmap-1)=i
int_mb(k_jso+symmap-1)=j
cgk debug
* write(*,'(a,3i3,f18.12)')
* & 'gk: i,j,ij,dbl_mb(so_x):',
* & i,j,numints+ijmap,dbl_mb(k_sox+symmap-1)
* write(*,'(a,3i3,f18.12)')
* & 'gk: i,j,ij,dbl_mb(so_y):',
* & i,j,numints+ijmap,dbl_mb(k_soy+symmap-1)
* write(*,'(a,3i3,f18.12)')
* & 'gk: i,j,ij,dbl_mb(so_z):',
* & i,j,numints+ijmap,dbl_mb(k_soz+symmap-1)
cgk end
endif
enddo
enddo
numints=numints+ijmap+1
end if
end do
end do
c call sifs
endif
call ga_sync() ! So that no nasty races can result
cgk debug
* write(*,*)'gk: numints = ', numints
* write(*,*)'gk: maxint = ', maxint
* symmap=0
* do i=1,nbft
* do j = 1, i
* write(*,'(a,2i8,3f18.12)')
* & 'gk: i,j,so_x,so_y,so_z:', int_mb(k_iso+symmap),
* & int_mb(k_jso+symmap),dbl_mb(k_sox+symmap),
* & dbl_mb(k_soy+symmap), dbl_mb(k_soz+symmap)
* symmap=symmap+1
* enddo
* enddo
do icart = 1, 3
itypea=2
if(icart.eq.1) then
kcart=k_soz
itypeb=2
elseif(icart.eq.2) then
kcart=k_soy
itypeb=1
elseif(icart.eq.3) then
kcart=k_sox
itypeb=0
endif
symmap=0
numtot=0
ibuf=0
do i=1,nbft
do j = 1, i
if(ibuf.eq.n1max) then
numtot = numtot + ibuf
cgk debug
* WRITE(*,*)"int_so_sifs calling sifew1 1"
* write(*,*)'gk: nrec=', nrec
* write(*,*)'gk: ibuf=', ibuf
* write(*,*)'gk: numtot=', numtot
* WRITE(*,*)"last=",last,"ibvtyp=",ibvtyp
cgk end
call sifew1(aoints, info, 2, ibuf, msame,
& itypea, itypeb, ifmt1, ibvtyp,
& sifval, clab, fcore, ibitv,
& sifbuf, nrec, ierr)
numtot = numtot - ibuf
endif
soval=dbl_mb(kcart+symmap)
if(abs(soval).gt.thresh) then
ibuf=ibuf+1
clab(1,ibuf)=int_mb(k_iso+symmap)
clab(2,ibuf)=int_mb(k_jso+symmap)
sifval(ibuf)=soval
endif
symmap=symmap+1
enddo
enddo
* if (icart.eq.3) then
* last=nomore
* else
last=nmsame
* endif
numtot=numtot+ibuf
cgk debug
* WRITE(*,*)"int_so_sifs calling sifew1 2"
* write(*,*)'gk: nrec=', nrec
* write(*,*)'gk: ibuf=', ibuf
* write(*,*)'gk: numtot=', numtot
* WRITE(*,*)"last=",last,"ibvtyp=",ibvtyp
cgk end
call sifew1(aoints, info, 2, ibuf, last,
& itypea, itypeb, ifmt1, ibvtyp,
& sifval, clab, fcore, ibitv,
& sifbuf, nrec, ierr)
* write(6,'(a,i10,1x,a,a)') 'Wrote ',numtot, sotype(icart),
* & ' integrals to aoints'
enddo
cgk end
cgk debug
* write(*,*)'gk: l_jso = ', l_jso
* write(*,*)'gk: k_jso = ', k_jso
cgk end
c
c chop stack at first item allocated
c
* if (.not. MA_pop_stack(l_sifval)) call errquit
* $ ('int_so_sifs: pop failed at l_sifval', 0, GA_ERR)
* if (.not. MA_pop_stack(l_sifbuf)) call errquit
* $ ('int_so_sifs: pop failed at l_sifbuf', 0, GA_ERR)
if (.not. MA_pop_stack(l_jso)) call errquit
$ ('int_so_sifs: pop failed at l_jso', 0, GA_ERR)
if (.not. MA_pop_stack(l_iso)) call errquit
$ ('int_so_sifs: pop failed at l_iso', 0, GA_ERR)
if (.not. MA_pop_stack(l_soz)) call errquit
$ ('int_so_sifs: pop failed at l_soz', 0, GA_ERR)
if (.not. MA_pop_stack(l_soy)) call errquit
$ ('int_so_sifs: pop failed at l_soy', 0, GA_ERR)
if (.not. MA_pop_stack(l_sox)) call errquit
$ ('int_so_sifs: pop failed at l_sox', 0, GA_ERR)
if (.not. MA_pop_stack(l_scr)) call errquit
$ ('int_so_sifs: pop failed at l_scr', 0, GA_ERR)
if (.not. MA_pop_stack(l_buf)) call errquit
$ ('int_so_sifs: pop failed at l_buf', 0, GA_ERR)
call ga_sync() ! So that no nasty races can result
if (oscfps) call pstat_off(ps_int_1e)
c
cgk debug
* write(*,*)'gk: leaving int_so_sifs at end of routine'
cgk end
end

496
src/NWints/aoints/nadct_trace.F Normal file
View file

@ -0,0 +1,496 @@
logical function nadct_trace(rtdb, nadct, trans)
c
c Compute the terms that go into the nonadiabatic derivative coupling
c terms
c f_I,J(R)^x = <I|d/dR_x J>
c
c Each term is computed by a trace of data from COLUMBUS with
c gradient integrals from NWChem.
c
c This is also computing gradients apparently
implicit none
#include "errquit.fh"
#include "bas.fh"
#include "geom.fh"
#include "rtdb.fh"
#include "mafdecls.fh"
#include "stdio.fh"
#include "cdft.fh"
#include "global.fh"
#include "util.fh"
#include "schwarz.fh"
#include "dftps.fh"
#include "zora.fh"
#include "rel_consts.fh"
#include "case.fh"
#include "dft_scaleig.fh"
#ifdef SCALAPACK
integer ga_cholesky, ga_llt_i
external ga_cholesky, ga_llt_i
#endif
c
integer rtdb
logical nadct
integer me, nproc, cinfo
logical IOLGC, oprint, oprint_time, oprint_ioinfo
logical ldirect
integer iVxc_opt, lmiss, ivcoul_opt
integer g_xcinv, g_oi, g_tmpb
integer lsvals, isvals, n_dep
double precision toll_s
integer nat
integer g_force
integer lforce
integer l_force, l_frc_nuc, l_frc_kin, l_frc_fds, l_frc_2el,
& l_frc_sow, l_frc_sox, l_frc_soy, l_frc_soz, l_frc_atr
integer k_force, k_frc_nuc, k_frc_kin, k_frc_fds, k_frc_2el,
& k_frc_sow, k_frc_sox, k_frc_soy, k_frc_soz, k_frc_atr
integer*4 aodens, aodensa
parameter(aodens=24, aodensa=25)
integer*4 ninfo, ninfmx, ntitmx, nsymx
parameter(ninfo=5,ninfmx=10,ntitmx=30,nsymx=8)
integer*4 infocd(ninfmx)
character*80 siftits(ntitmx)
integer*4 nsym, nmot, nmbpsy(nsymx), nsift
integer*4 infocda(ninfmx)
character*80 siftitsa(ntitmx)
integer*4 nsifta
LOGICAL :: trans ! is this a transition property run?
cgk debug
* write(*,*)'gk: in nadct_trace'
cgk end
nadct_trace=.true.
c
me = ga_nodeid()
nproc = ga_nnodes()
oprint = util_print('dft cinfo', print_default)
oprint_time = util_print('dft timings', print_high)
oprint_ioinfo = util_print('io info', print_debug)
if(ga_nodeid().eq.0.and.
. util_print('geometry', print_high)) then
if(.not.geom_print(geom)) call errquit(
. ' nadct_trace: geom print failed',0, GEOM_ERR)
endif
IOLGC = .TRUE.
if (noio.eq.1) IOLGC = .FALSE.
!#ifdef SCALAPACK
!#if defined(PARALLEL_DIAG)
! call slexit
! call slexit2
! if (.not. rtdb_cget(rtdb, 'dft:scaladiag',1,yev))
! & yev='d'
! if (.not. rtdb_get(rtdb, 'dft:scalanb',mt_int,1,nb))
! & nb=0
!#endif
!#endif
c
c Determine the characteristics of the AO Gaussian basis set.
c This uses the old NWdft data structures, rather than tthe rtdb
c
if (.not. bas_numbf(AO_bas_han, nbf_ao)) then
call errquit('Exiting from nadct_trace',1, BASIS_ERR)
endif
cgk debug
* write(*,*)'in nadct_trac, AO_bas_han = ', AO_bas_han
cgk end
* iter = 0
if(nproc.gt.1) then
write(6,'(a,a)')'WARNING: Writing aoints, calculating 1e ints ',
& 'with 1 CPU'
write(6,'(a,a)')'WARNING: we are currently supporting 1 ',
& 'process for aoints'
endif
c
c Check zora options
do_zora_get_correction=.false.
if (rtdb_get(rtdb,'zora:calc',mt_log,1,
& do_zora_get_correction))
& do_zora_get_correction=.true.
c
do_zora = .false.
if (rtdb_get(rtdb,'zora:on',mt_log,1,do_zora))
& do_zora = .true.
c
c == Get the grid cutoff or set default ==
if(.not.rtdb_get(rtdb,'zora:cutoff',mt_dbl,1,zoracutoff))
& zoracutoff = 1.0d-08
c
c Get the light speed and other constants, if specified
call rel_init_con(rtdb)
clight_au = cau
c
*c Get and check the Coulomb Attenuation Method (CAM) parameters
* if (.not.rtdb_get(rtdb, 'dft:cam_exch', mt_log, 1,
* & cam_exch)) cam_exch=.false.
* if (.not.rtdb_get(rtdb, 'dft:cam_omega', mt_dbl, 1,
* & cam_omega)) cam_omega=0.d0
* if (.not.rtdb_get(rtdb, 'dft:cam_alpha', mt_dbl, 1,
* & cam_alpha)) cam_alpha=0.d0
* if (.not.rtdb_get(rtdb, 'dft:cam_beta', mt_dbl, 1,
* & cam_beta)) cam_beta=0.d0
*c
* ldirect = .false.
* if (cam_exch .and. (.not.CDFIT) ) then
* if (rtdb_get(rtdb, 'dft:direct', mt_log, 1, ldirect)) then
* if (.not.ldirect) then
* if (ga_nodeid().eq.0) then
* write(luout,*) "CAM type calculations without CDFIT are only"
* write(luout,*) "available with the direct option"
* end if
* call errquit('Please use the direct option',4, INPUT_ERR)
* end if
* end if
* end if
c
c Select method to compute the XC potential.
if (XC_bas_han.eq.-99999)then
iVxc_opt = 0
XCFIT = .FALSE.
else
call errquit('nadct_trace: XC basis fitting not implemented',
& 0,CAPMIS_ERR)
endif
* if (XCFIT)then
*c
*c Determine the characteristics of the XC Gaussian basis set.
*c
* if (.not. bas_numbf(XC_bas_han, nbf_xc)) then
* call errquit('Exiting from scfdriver',4, BASIS_ERR)
* endif
*c
*c Compute the matrix inverse of the 2-ctr OIs.
*c
**ga:1:0
* if (.not. ga_create(mt_dbl, nbf_xc, nbf_xc, 'XCinv',
* & nbf_xc, 0, g_xcinv))
* & call errquit('nadct_trace: error creating g_xcinv',
* & 0,GA_ERR)
* call ga_zero(g_xcinv)
**ga:1:0
* if (.not. ga_create(mt_dbl, nbf_xc, nbf_xc, 'g_oi',
* & nbf_xc, 0, g_oi))
* & call errquit('nadct_trace: error creating g_oi',0,
* & GA_ERR)
* call ga_zero(g_oi)
* call int_1e_ga(xc_bas_han, xc_bas_han,
* . g_oi, 'overlap', oskel)
* if (oskel)call sym_symmetrize(
* ( geom, xc_bas_han, .false., g_oi)
*#if defined(PARALLEL_DIAG)
* call ga_chol(g_oi, g_xcinv, cinfo)
*#else
* call ga_chol_seq(g_oi, g_xcinv, cinfo)
*#endif
* if (cinfo.ne.0)then
* if (me.eq.0)then
* write(LuOut,*)' Problem in performing a Choleski '
* write(LuOut,*)' decomposition of the 2-ctr overlap '
* write(LuOut,*)' matrix using XC fitting basis. '
* write(LuOut,*)' Attempting a diag/inverse. '
* endif
* endif
* if (cinfo.eq.0) then
*#if defined(PARALLEL_DIAG)
* call ga_inverse(g_xcinv, g_oi)
* if (odftps) call pstat_on(ps_bdgemm)
* call ga_dgemm('T', 'N', nbf_xc, nbf_xc, nbf_xc, 1.d0,
* & g_oi, g_oi, 0.d0, g_xcinv)
* if (odftps) call pstat_off(ps_bdgemm)
*#else
* call ga_copy(g_xcinv, g_oi)
* call ga_inv_seq(g_oi, g_xcinv)
*#endif
* else
* call dft_invdiag(g_oi, g_xcinv, nbf_xc)
* endif
* if (.not. ga_destroy(g_oi)) call errquit
* & ('nadct_trace: could not destroy g_oi', 0, GA_ERR)
* if (me.eq.0.and.IOLGC) then
* lmiss = 0
* call dft_invio('XC', g_xcinv, nbf_xc, 'WRITE', lmiss)
* if (lmiss.ne.0)call errquit
* & (' nadct_trace; dft_invio - abnormal write of XC ', 0,
* & DISK_ERR)
* endif
* if (IOLGC)then
* if (.not. ga_destroy(g_xcinv)) call errquit
* & ('nadct_trace: could not destroy g_xcinv', 0, GA_ERR)
* endif
* endif
*c
c Done with XC setup stuff....
c Select method to compute Coulomb potential.
c
if (CD_bas_han.eq.-99999)then
iVcoul_opt = 0
CDFIT = .FALSE.
else
call errquit('nadct_trace: charge density fitting not implement
&ed', 0,CAPMIS_ERR)
endif
c
c
* if (CDFIT)then
*c
*c Determine the characteristics of the CD Gaussian basis set.
*c
* if (odftps) call pstat_on(ps_2cer)
* if (.not. bas_numbf(CD_bas_han, nbf_cd)) then
* call errquit('Exiting from scfdriver',5, BASIS_ERR)
* endif
*c
*c Compute the matrix inverse of the CD 2-ctr ERIs.
*c
* if (.not. ga_create(mt_dbl, nbf_cd, nbf_cd, 'CD 2cERI',
* & 0, nbf_cd, g_2ceri))
* & call errquit('nadct_trace: error creating g_2ceri',0,
* & GA_ERR)
* call ga_zero(g_2ceri)
* call dft_get2eri(CD_bas_han, g_2ceri,oskel)
* if (odftps) call pstat_off(ps_2cer)
* if (oskel)call
* . sym_symmetrize(geom,cd_bas_han,.false.,g_2ceri)
* call ga_sync()
* if (.not. ga_duplicate(g_2ceri, g_cdinv, 'CD 2cERInv'))
* & call errquit('nadct_trace: error creating g_cdinv',0,
* & GA_ERR)
* if (odftps) call pstat_on(ps_inv)
* call ga_zero(g_cdinv)
* lmiss = 1
*c if (lmiss.eq.1) then
* call ga_zero(g_cdinv)
* cinfo = 0
*#if defined(PARALLEL_DIAG)
*#ifdef SCALAPACK
* call ga_copy(g_2ceri, g_cdinv)
* call ga_sync()
* cinfo= ga_cholesky('U',g_cdinv)
*#else
* call ga_chol(g_2ceri, g_cdinv, cinfo)
*#endif
*#else
* call ga_chol_seq(g_2ceri, g_cdinv, cinfo)
*#endif
* if (cinfo.ne.0)then
* if (me.eq.0)then
* write(LuOut,*)' Problem in performing a Choleski '
* write(LuOut,*)' decomposition of the 2-ctr ERI '
* write(LuOut,*)' matrix using CD fitting basis. '
* write(LuOut,*)' Attempting a diag/inverse. '
* endif
* endif
* if (cinfo.eq.0) then
* g_tmpb = g_2ceri
*#if defined(PARALLEL_DIAG)
*#ifdef SCALAPACK
* cinfo = ga_llt_i('U',g_cdinv,-1)
* if (cinfo.ne.0)then
* if (me.eq.0)then
* write(LuOut,*)' Problem in performing a Invers. '
* write(LuOut,*)' of the 2-ctr ERI '
* endif
* call ga_sync
* call errquit('nadct_trace: inverse failed ',0,0)
* endif
*
*#else
* call ga_inverse(g_cdinv, g_tmpb)
*#endif
* if (odftps) call pstat_on(ps_bdgemm)
*#ifndef SCALAPACK
* call ga_dgemm('T', 'N', nbf_cd, nbf_cd, nbf_cd, 1.d0,
* & g_tmpb, g_tmpb, 0.d0, g_cdinv)
*#endif
* if (odftps) call pstat_off(ps_bdgemm)
*#else
* call ga_copy(g_cdinv, g_tmpb)
* call ga_inv_seq(g_tmpb, g_cdinv)
*#endif
* else
* call dft_invdiag(g_2ceri, g_cdinv,
* & nbf_cd)
* endif
* if (odftps) call pstat_off(ps_inv)
*#ifndef SCALAPACK
*c
*c second build of g_2ceri needed becuase previous calls destroyed it
*c
* if (odftps) call pstat_on(ps_2cer)
* call ga_zero(g_2ceri)
* call dft_get2eri(CD_bas_han, g_2ceri,oskel)
* if (oskel)call
* . sym_symmetrize(geom,cd_bas_han,.false.,g_2ceri)
* if (odftps) call pstat_off(ps_2cer)
*#endif
* if (IOLGC.and.(me.eq.0)) then
* lmiss = 0
* call dft_invio('CDI', g_cdinv, nbf_cd, 'WRITE', lmiss)
* if (lmiss.ne.0)call errquit
* & (' nadct_trace; dft_invio - abnormal write of CDI ', 0
* & ,DISK_ERR)
* lmiss = 0
* call dft_invio('CD', g_2ceri, nbf_cd, 'WRITE', lmiss)
* if (lmiss.ne.0)call errquit
* & (' nadct_trace; dft_invio - abnormal write of CD ', 0,
* & DISK_ERR)
* endif
* if (oprint_ioinfo)then
* call ga_print(g_cdinv)
* call ga_print(g_2ceri)
* endif
*c endif
* if (IOLGC) then
* if (.not. ga_destroy(g_cdinv)) call errquit
* & ('nadct_trace: could not destroy g_xcinv', 0, GA_ERR)
* if (.not. ga_destroy(g_2ceri)) call errquit
* & ('nadct_trace: could not destroy g_xcinv', 0, GA_ERR)
* endif
* endif
*c
c Done with CD setup stuff
c
if (.not.MA_Push_Get(MT_Dbl, nbf_ao, 'ovl eig vals', lsvals,
& isvals))
& call errquit('nadct_trace: cannot allocate ovl eig vals', 0,
& MA_ERR)
call dfill(nbf_ao, 0.0d0, dbl_mb(isvals), 1)
c
c Set orbital overlap tolerance
c
call util_tolls(rtdb,.false.,toll_s,n_dep,dbl_mb(isvals),nbf_ao)
c
c open aodens
c
if(me.eq.0) then
if(ga_nodeid().eq.0) then
open(unit=aodens,status='old',form='unformatted',
& file='aodens')
rewind aodens
endif
cgk debug
* write(*,*)'gk: calling rdhcid'
cgk end
call rdhcid(aodens, infocd, nsift, siftits, nsym, nmot, nmbpsy)
cgk debug
* write(*,*)'gk: back from rdhcid'
* write(*,*)'gk: nsym = ', nsym
cgk end
endif
if (nmot.ne.nbf_ao) then
call errquit('nadct_trace: aodens number of orbitals does not m
&atch NWChem orbital number', nmot, BASIS_ERR)
endif
if (nsym.ne.1) then
call errquit('nadct_trace: only c1 symmetry implented...invalid
& aodens, nsym=', nsym, BASIS_ERR)
endif
c open aodensa
c
IF (trans) THEN
if(me.eq.0) then
if(ga_nodeid().eq.0) then
open(unit=aodensa,status='old',form='unformatted',
& file='aodens.asym')
rewind aodens
endif
call rdhcid(aodensa, infocda,nsifta,siftitsa,nsym,nmot,nmbpsy)
endif
if (nmot.ne.nbf_ao) then
call errquit('nadct_trace: aodens number of orbitals does not m
&atch NWChem orbital number', nmot, BASIS_ERR)
endif
if (nsym.ne.1) then
call errquit('nadct_trace: only c1 symmetry implented...invalid
& aodens, nsym=', nsym, BASIS_ERR)
endif
ENDIF
c set up force arrays
if (.not. geom_ncent(geom,nat))
$ call errquit('nadct_trace: could not get natoms',0, GEOM_ERR)
c global force array
if (.not. ga_create(mt_dbl, 3, nat, 'forces', 3, 0, g_force))
$ call errquit('nadct_trace: failed to create force GA',0,
& GA_ERR)
call ga_zero(g_force)
c local replication for separate pieces
lforce = nat * 3
if (.not. ma_push_get(mt_dbl,lforce,'forces',l_force,k_force))
$ call errquit('nadct_trace:could not allocate l_force',lforce,
& MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_nuc,
$ k_frc_nuc)) call errquit
& ('nadct_trace:could not allocate l_frc_nuc',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_kin,
$ k_frc_kin)) call errquit
& ('nadct_trace:could not allocate l_frc_kin',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_fds,
$ k_frc_fds)) call errquit
& ('nadct_trace:could not allocate l_frc_fds',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_2el,
$ k_frc_2el)) call errquit
& ('nadct_trace:could not allocate l_frc_2el',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_sow,
$ k_frc_sow)) call errquit
& ('nadct_trace:could not allocate l_frc_sow',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_sox,
$ k_frc_sox)) call errquit
& ('nadct_trace:could not allocate l_frc_sox',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_soy,
$ k_frc_soy)) call errquit
& ('nadct_trace:could not allocate l_frc_soy',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_soz,
$ k_frc_soz)) call errquit
& ('nadct_trace:could not allocate l_frc_soz',lforce, MA_ERR)
if (.not. ma_push_get(mt_dbl, lforce, 'forces', l_frc_atr,
$ k_frc_atr)) call errquit
& ('nadct_trace:could not allocate l_frc_atr',lforce, MA_ERR)
call dfill(lforce, 0.0D0, dbl_mb(k_force), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_nuc), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_kin), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_fds), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_2el), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_sow), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_sox), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_soy), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_soz), 1)
call dfill(lforce, 0.0D0, dbl_mb(k_frc_atr), 1)
c
c Calculate non spin-orbit energy gradient terms
c trace( del(h)*D + del(g)*d - del(S)*F )
c
cgk debug
* write(*,*)'gk: calling egrad_trace'
* write(*,*)'gk: nmot = ', nmot
* write(*,*)'gk: nbf_ao = ', nbf_ao
cgk end
call egrad_trace(AO_bas_han, aodens, ninfo, infocd, nbf_ao,
& g_force, dbl_mb(k_frc_nuc), dbl_mb(k_frc_kin),
& dbl_mb(k_frc_fds), dbl_mb(k_frc_2el), dbl_mb(k_frc_sow),
& dbl_mb(k_frc_sox), dbl_mb(k_frc_soy), dbl_mb(k_frc_soz),
& dbl_mb(k_force), nat, rtdb, nadct, aodensa, infocda, trans)
cgk debug
* write(*,*)'gk: back from egrad_trace'
cgk end
close(unit=aodens)
IF(trans) CLOSE(UNIT=aodensa)
1000 continue
cgk debug
* write(*,*)'gk: leaving nadct_traces'
cgk end
return
end

29
src/NWints/aoints/print_dint_block.F Normal file
View file

@ -0,0 +1,29 @@
subroutine print_dint_block(ilo, ihi, jlo, jhi,
& klo, khi, llo, lhi, geri, idatom)
integer ilo, ihi, jlo, jhi, klo, khi, llo, lhi
double precision geri(llo:lhi, klo:khi, jlo:jhi, ilo:ihi, 3, 4)
integer idatom(4)
character*4 dcart(3)
data dcart/'dx','dy','dz'/
integer i, j, k, l, icart, iat
do iat = 1, 4
if(idatom(iat).ge.1) then
do i=ilo, ihi
do j=jlo, jhi
do k=klo, khi
do l=llo, lhi
do icart = 1, 3
write(6,2000) idatom(iat),dcart(icart),'[',i,',',j,';',
& k,',',l,'] = ', geri(l,k,j,i,icart,iat)
2000 format(i2,a2,4(a1,i2),a,f20.10)
enddo
enddo
enddo
enddo
enddo
endif
enddo
return
end

18
src/NWints/aoints/print_soblock.F Normal file
View file

@ -0,0 +1,18 @@
subroutine print_soblock(H, ilo, ihi, jlo, jhi, iat, icart,
& socomp, nat)
integer ilo, ihi, jlo, jhi, iat, icart, socomp, nat
double precision H(jlo:jhi, ilo:ihi, 3, 3, nat)
character*2 dcart(3)
data dcart/'dx','dy','dz'/
integer i, j
do i=ilo, ihi
do j=jlo, jhi
write(6,2000) iat,dcart(icart),'<',i,'|Hso',socomp,'|',j,'>=',
& H(j,i,socomp,icart,iat)
2000 format(i2,a2,3(a,i2),a,f20.10)
enddo
enddo
return
end

100
src/NWints/aoints/rd1mat.F Normal file
View file

@ -0,0 +1,100 @@
subroutine rd1mat(ntape, lbuf, nipbuf, infomo, buf, ilab, val,
& d1, nnmt, sym, map, nsym, nmpsy, kntin, itypea, itypeb, nbf)
c
c read a 1-e density matrix
c
c input:
c ntape = input unit number.
c lbuf = buffer length.
c nipbuf = maximum number of integrals in a buffer.
c infomo = info array
c buf(*) = buffer for integrals and packed labels.
c ilab(*)= unpacked labels.
c val(*) = unpacked values.
c sym(*) = symmetry of each orbital
c nnmt = total size of symmetric, symmetry-packed, lower-triangular-
c packed matrix.
c map(*) = orbital index mapping vector.
c
c output:
c d1(*) = 1-e density matrix
c
implicit none
c
integer*4 lbuf, nipbuf, ntape, nnmt, nsym, nbf
integer*4 itypea, itypeb
integer*4 nipv, msame, nmsame, nomore
parameter(nipv=2, msame=0,nmsame = 1, nomore = 2 )
integer*4 last, nrec, ierr
real*8 buf(lbuf),val(nipbuf),fcore
double precision d1(nnmt)
integer*4 ilab(nipv,nipbuf),sym(*),map(*)
integer*4 nmpsy(*)
integer*4 infomo(*),kntin(*)
integer*4 iretbv, symoff(1),lasta,lastb
parameter(iretbv = 0)
real*8 zero
parameter(zero = 0.0d0)
integer*4 wrnerr, nfterr, faterr
parameter(wrnerr=0,nfterr=1,faterr=2)
integer*4 i
integer numd
c
cgk debug
* write(*,*)'gk: calling sifr1x from rd1mat'
* write(*,*)'gk: ntape = ', ntape
* write(*,*)'gk: itypea = ', itypea
* write(*,*)'gk: itypeb = ', itypeb
* write(*,*)'gk: nipbuf = ', nipbuf
* write(*,*)'gk: lbuf = ', lbuf
* write(*,*)'gk: nsym = ', nsym
* write(*,*)'gk: nnmt = ', nnmt
* write(*,*)'gk: nbf = ', nbf
* write(*,*)'gk: nmpsy(*) = ', (nmpsy(i),i=1,nsym)
* write(*,*)'gk: infomo(1) = ', infomo(1)
* write(*,*)'gk: infomo(2) = ', infomo(2)
* write(*,*)'gk: infomo(3) = ', infomo(3)
* write(*,*)'gk: infomo(4) = ', infomo(4)
* write(*,*)'gk: infomo(5) = ', infomo(5)
* write(*,*)'gk: symoff(1) = ', symoff(1)
cgk end
fcore= zero
symoff(1)=0
c
cgk debug
* write(*,*)'gk:calling dfill'
* write(*,*)' nnmt = ', nnmt
* call flush(6)
cgk end
numd=nnmt
call dfill(numd,0d0,d1,1)
cgk debug
* write(*,*)'gk:back from dfill'
* call flush(6)
cgk end
do i=1,nbf
map(i)=i
enddo
cgk debug
* write(*,*)'gk: calling sifr1x'
* WRITE(*,*)"LB: itypeb=",itypeb
cgk end
call sifr1x(
1 ntape, infomo, itypea, itypeb,
2 nsym, nmpsy, symoff, buf,
3 val, ilab, map, sym,
4 nnmt, d1, fcore, kntin,
5 lasta,lastb, last, nrec,
6 ierr)
cgk debug
* write(*,*)'gk: back from sifr1x in rd1mat'
* write(*,*)'gk: nmpsy(1)=', nmpsy(1)
cgk end
if ( ierr.ne. 0) then
call bummer('rd1mat: from sifr1x, ierr=',ierr,faterr)
elseif (fcore .ne. zero) then
write(6,*)' frozen core orbitals are not yet allowed'
call bummer('rd1mat: from sifr1x(), fcore=',fcore,wrnerr)
endif
return
end

67
src/NWints/aoints/rd_d2bl.F Normal file
View file

@ -0,0 +1,67 @@
c
c***********************************************************************
c
*deck indast
subroutine rd_d2bl(
& vals, labs, nbuf, buffer,
& last, sfile, info)
c
c
implicit none
c
c ## paramter & common block section
c
integer*4 nipv,iretbv
parameter (nipv=4,iretbv=0)
c
integer*4 ninfmx
parameter(ninfmx=10)
c
c # bummer error types.
integer*4 wrnerr, nfterr, faterr
parameter(wrnerr=0,nfterr=1,faterr=2)
c
c ## integer*4 section
c
integer*4 info(ninfmx),itypea,itypeb,ifmt,ibvtyp,ierr,ibitv(1)
integer*4 last,labs(4,*)
integer*4 nbuf
integer*4 sfile
c
c ## real*8
c
real*8 buffer(*)
real*8 vals(*)
c
c
c-----------------------------------------------------
c
cgk debug
! write(*,*)'gk: before sifrd2'
! write(*,*)'gk: sfile = ', sfile
! write(*,*)'gk: info(1):', info(1)
! write(*,*)'gk: info(2):', info(2)
! write(*,*)'gk: info(3):', info(3)
! write(*,*)'gk: info(4):', info(4)
! write(*,*)'gk: info(5):', info(5)
! write(*,*)'gk: last = ', last
cgk end
call sifrd2(sfile, info, nipv, iretbv,
& buffer, nbuf, last, itypea,
& itypeb, ifmt, ibvtyp, vals,
& labs, ibitv, ierr)
cgk debug
! write(*,*)'gk: after sifrd2'
! write(*,*)'in rd_d2bl, itypea = ', itypea
! write(*,*)'in rd_d2bl, itypeb = ', itypeb
! write(*,*)'in rd_d2bl, last = ', last
! write(*,*)'in rd_d2bl, labs(1,1) = ', labs(1,1)
! write(*,*)'in rd_d2bl, labs(2,1) = ', labs(2,1)
! write(*,*)'in rd_d2bl, labs(3,1) = ', labs(3,1)
! write(*,*)'in rd_d2bl, labs(4,1) = ', labs(4,1)
! write(*,*)' vals(1) = ', vals(1)
cgk end
c
if (ierr.gt.0) call bummer("error in sifrd2",0,faterr)
return
end

75
src/NWints/aoints/rdhcid.F Normal file
View file

@ -0,0 +1,75 @@
c
c**********************************************************************
c
*deck rdhcid
subroutine rdhcid(denfl,infocd,ntitle,title,nsym,nmot,
& nbpsy)
c
c read the header info from the ci density matrix file.
c
c written by ron shepard.
c
implicit none
c
character*80 title(*)
integer*4 infocd(*)
integer*4 ntitle,nsym,nmot,ninfo,nenrgy,ierr,nmap
integer*4 nmotx
parameter (nmotx=255)
integer*4 i, ntitmx, denfl,idummy, nsymx, nengmx
parameter( ntitmx = 30, nsymx=1, nengmx=20)
real*8 cienrg(nengmx)
character*4 slabel(nsymx)
character*8 cidlab(nmotx)
integer*4 nbpsy(*),cietyp(nengmx)
c
real*8 one, zero
parameter(one=1d0, zero=0d0)
c
integer*4 nmotn, nsymn
c
c bummer error types.
integer wrnerr, nfterr, faterr
parameter(wrnerr=0,nfterr=1,faterr=2)
c
rewind denfl
call sifrh1(
1 denfl, ntitle, nsym, nmot,
2 ninfo, nenrgy, nmap, ierr)
c
cgk debug
* write(*,*)'gk: in rdhcid, nsym=',nsym
* write(*,*)'gk: nmot = ', nmot
* write(*,*)'gk: nmap = ', nmap
* write(*,*)'gk: nenrgy = ', nenrgy
cgk end
if (ierr .ne. 0) then
call bummer('rdhcid: from sifrh1, ierr=',ierr,faterr)
elseif (nenrgy .gt. nengmx) then
call bummer('rdhcid: from sifrh1, nenrgy=',ntitle,faterr)
elseif (nsym .gt. nsymx) then
call bummer('rdhcid: from sifrh1, nsymx=',nsymx,faterr)
elseif (ntitle .gt. ntitmx) then
call bummer('rdhcid: from sifrh1, ntitle=',ntitle,faterr)
elseif (ninfo .gt. 10) then
call bummer('rdhcid: from sifrh1, ninfo=',ninfo,wrnerr)
endif
c #ignore map
nmap = 0
call sifrh2(
1 denfl, ntitle, nsym, nmot,
2 ninfo, nenrgy, nmap, title,
3 nbpsy, slabel, infocd, cidlab,
4 cietyp, cienrg, idummy, idummy,
5 ierr)
write(*,*)'ntitle=',ntitle
c
20 continue
c
write(*,6020)(title(i),i=1,ntitle)
6020 format(/' ci density file header information:'/(1x,a))
write(*,6040)(i,nbpsy(i),i=1,nsym)
6040 format(' (isym:nbpsy)',8(i3,':',i3))
c
return
end

154
src/NWints/aoints/sequ.F Normal file
View file

@ -0,0 +1,154 @@
integer function indx(a,b,c,d,lena,lenb,lenc,lend,s,t,u,v,dbg,
& iimn, jjmn, kkmn, llmn)
c safe but inefficient brute force shell index computer
c replace with direct computation of index
integer a, b, c, d
integer lena, lenb, lenc, lend
integer s, t, u, v
integer iimn, jjmn, kkmn, llmn
integer ii, jj, kk, ll
integer i, j, k, l, imn, imx, jmn, jmx, lmx, kmx
logical dbg
integer ij, kl
integer st, uv
* write(26,*) 'in indx'
st=(s*(s-1))/2+t
uv=(u*(u-1))/2+v
if(dbg) then
write(26,100)'calculating local shell quartet index for ',
& a,b,c,d
write(26,100)' shell sizes: ', lena, lenb, lenc, lend
endif
100 format(a,4i2)
if(st.lt.uv) then
write(26,*) 'u == v and d > c'
indx=-1
goto 1000
endif
if(u.eq.v.and.d.gt.c) then
write(26,*) 'u == v and d > c'
indx=-1
goto 1000
endif
* if(t.eq.v.and.d.gt.b) then
* write(26,*) 't == v and d > b'
* indx=-1
* goto 1000
* endif
* if(s.eq.u.and.c.gt.a) then
* write(26,*) 's == u and c > a'
* indx=-1
* goto 1000
* endif
if(s.eq.t.and.b.gt.a) then
write(26,*) 's == t and b > a'
indx=-1
goto 1000
endif
if(t.gt.s) then
write(26,*) 't > s'
indx=-1
goto 1000
endif
if(v.gt.u) then
write(26,*) 'v > u'
indx=-1
goto 1000
endif
indx=0
* i=1
* j=1
* k=1
* l=1
do k=1,c
if (v.eq.u) then
* if(k.lt.c.or.j.lt.b.or.i.lt.a) then
if(k.lt.c) then
lmx=k
else
lmx=d
endif
else
* if(k.lt.c.or.j.lt.b.or.i.lt.a) then
if(k.lt.c) then
lmx=lend
else
lmx=d
endif
endif
kk=k+kkmn-1
do l=1,lmx
* if (u.eq.s) then
* imn=k
* else
imn=1
* endif
if(l.lt.d.or.k.lt.c) then
imx=lena
else
imx=a
endif
ll=l+llmn-1
kl=((kk-1)*kk)/2+ll
do i=imn,imx
* write(*,*)'i = ', i, ' b = ', b, ' lenb = ', lenb
if(s.eq.t) then
if(l.lt.d.or.k.lt.c.or.i.lt.a) then
jmx=i
else
jmx=b
endif
else
if(l.lt.d.or.k.lt.c.or.i.lt.a) then
jmx=lenb
else
jmx=b
endif
endif
jmn=1
* write(*,*)'jmn=',jmn,' jmx=',jmx
ii=iimn+i-1
do j=jmn,jmx
jj=jjmn+j-1
ij=((ii-1)*ii)/2+jj
* write(26,2323) '[',i,',',j,';',k,',',l,']'
2323 format(4(a,i2),a1)
* write(26,'(2(a,i3))') '(ij)=',ij,' (kl)=',kl
if(st.eq.uv) then
if(ij.ge.kl) then
indx=indx+1
if(dbg) then
write(26,'(4i2,a,i3)') i,j,k,l, ' - ', indx
endif
endif
else
indx=indx+1
if(dbg) then
write(26,'(4i2,a,i3)') i,j,k,l, ' - ', indx
endif
endif
enddo
enddo
enddo
enddo
if(dbg) then
write(26,'(a,i4)') 'local index=', indx
endif
* write(26,*) 'leaving indx'
1000 continue
return
end

182
src/NWints/aoints/sif2arr.F Normal file
View file

@ -0,0 +1,182 @@
subroutine sif2arr(aodens, d2, mx_d2, valbuf, nipv, labbuf, buf,
& infocd, b2s, sminb, smaxb, shoff, nsp, nshell, dbg)
#include "errquit.fh"
integer*4 aodens
double precision d2(mx_d2)
integer mx_d2
double precision valbuf(*)
integer nipv
integer*4 labbuf(nipv,*)
double precision buf(*)
integer*4 infocd(*)
integer b2s(*)
integer sminb(*)
integer smaxb(*)
integer shoff(nsp,nsp)
integer nsp
integer nshell
integer dbg(4,*)
integer*4 nbuf, last
integer i, j, k, l, ij, kl
integer i0, j0, k0, l0
integer s, t, u, v, st, uv
integer s0, t0, u0, v0, st0, uv0
integer is, js, ks, ls, a, b, c, d, ab, cd, abcd
integer imin, jmin, kmin, lmin, imax, jmax
integer slen, tlen, stlen, ulen, vlen
integer ints, ijkl, ab1, cd1
logical sort_err
integer indx
external indx
integer idbg
integer imn, jmn, kmn, lmn
cgk debug
* write(*,*) 'gk: entered sif2arr'
cgk end
do ints=1,mx_d2
d2(ints)=0d0
enddo
100 continue
cgk debug
* write(*,*)'gk: calling rd_d2bl'
cgk end
call rd_d2bl(valbuf, labbuf, nbuf, buf, last, aodens, infocd)
cgk debug
* write(*,*)'gk: back from rd_d2bl'
* write(*,*)'gk: nbuf = ', nbuf
cgk end
c loop over integrals in buffer and place in shell block order
do ints=1,nbuf
c create mapping from [ij;kl] in AO basis order to [ab;cd] shell
c blocked order
sort_err=.false.
i0=labbuf(1,ints)
j0=labbuf(2,ints)
k0=labbuf(3,ints)
l0=labbuf(4,ints)
* ij=((i-1)*i)/2+j
* kl=((k-1)*k)/2+l
u0=b2s(k0) ! shell for AO k
v0=b2s(l0) ! shell for AO l
uv0=(u0*(u0-1))/2+v0
s0=b2s(i0) ! shell for AO i
t0=b2s(j0) ! shell for AO j
st0=(s0*(s0-1))/2+t0
* write(26,*)'gk:# start index computation'
if(uv0.gt.st0) then ! switch shell pairs
i=k0
j=l0
k=i0
l=j0
s=u0
t=v0
u=s0
v=t0
st=uv0
uv=st0
else
i=i0
j=j0
k=k0
l=l0
s=s0
t=t0
u=u0
v=v0
st=st0
uv=uv0
endif
kmn=sminb(u)
ks=k-kmn+1 ! min AO index in shell u
lmn=sminb(v)
ls=l-lmn+1 ! min AO index in shell v
ulen=smaxb(u)-kmn+1
vlen=smaxb(v)-lmn+1
c construct the (cd) for indexing within the uv shell pair
* if(u.eq.v) then
* c=max(ks,ls)
* d=min(ks,ls)
* cd=((c-1)*c)/2+d ! linearized index within uv shell pair
* else
* c=ks
* d=js
* cd=(ls-1)*ulen+ks ! linearized index within uv shell pair
* endif
imn=sminb(s)
is=i-imn+1 ! min index for shell s
jmn=sminb(t)
js=j-jmn+1 ! min index of shell t
slen=smaxb(s)-imn+1
tlen=smaxb(t)-jmn+1
abcd=indx(is,js,ks,ls,slen,tlen,ulen,vlen,s,t,u,v,.false.,
& imn,jmn,kmn,lmn)
ijkl=shoff(st,uv)+abcd ! sorted index
if(dbg(1,ijkl).ne.i) sort_err=.true.
if(dbg(2,ijkl).ne.j) sort_err=.true.
if(dbg(3,ijkl).ne.k) sort_err=.true.
if(dbg(4,ijkl).ne.l) sort_err=.true.
d2(ijkl)=valbuf(ints)
cgk debug
* write(26,2321) 's=',s,' t=',t,' u=',u,' v=', v
* write(26,2321) 'sl=',slen,' tl=',tlen,' ul=',ulen,' vl=', vlen
2321 format(4(a,i2))
* write(26,2312) '(st)=', st, ' (uv)=',uv
* if (st .lt. uv) write(26,*) 'st < uv'
2312 format(2(a,i2))
* write(26,2322) '[ab;cd] = [',is,',',js,';',ks,',',ls,']'
* write(26,2312) 'ab=', ab, ' cd=',cd
2322 format(4(a,i2),a1)
* write(26,'(2(a,i6))') 'shoff=',shoff(st,uv), ' abcd=', abcd
* write(26,'((a,i6))') 'ijkl=',ijkl
* write(26,2322) '[',i,',',j,';',k,',',l,']'
* write(26,2324) 'd2(',i,',',j,';',k,',',l,')=',d2(ijkl)
* write(26,2322) 'debug(',dbg(1,ijkl),',',dbg(2,ijkl),';',
* & dbg(3,ijkl),',',dbg(4,ijkl),')'
2324 format(4(a,i2),a,f20.12)
cgk end
if(sort_err) then
* write(26,*)'is = ', is
* write(26,*)'s = ', s
* write(26,*)'sminb(s) = ', sminb(s)
* write(26,*)'is = ', is
do idbg=1,nshell
write(26,'(a,i2,a,i3)')'sminb(',idbg,') = ', sminb(idbg)
enddo
abcd=indx(slen,tlen,ulen,vlen,slen,tlen,ulen,vlen,s,t,u,v,
& .true., imn,jmn,kmn,lmn)
abcd=indx(is,js,ks,ls,slen,tlen,ulen,vlen,s,t,u,v,.true.,
& imn,jmn, kmn,lmn)
call errquit ('sif2arr: sort error', 0, GA_ERR)
endif
enddo
if (last.eq.2) goto 1000
go to 100
1000 continue
cgk debug
* write(*,*)'gk: leaving sif2arr'
cgk end
return
end

276
src/NWints/aoints/sif2da_2e_sort.F Normal file
View file

@ -0,0 +1,276 @@
subroutine sif2da_2e_sort(aodens, l2rec, n2max, info, ibas, nbft,
& nshell, b2s, sminb, smaxb, frc_2el, nat, lbuf, buf, lscr,
& scr)
#include "bas.fh"
#include "errquit.fh"
#include "mafdecls.fh"
integer*4 aodens
integer*4 l2rec
integer*4 n2max
integer*4 info(*)
integer ibas
integer nbft
integer nshell
integer b2s(nbft) ! basis number to shell mapping array (output)
integer sminb(nshell) ! min basis in shell (output)
integer smaxb(nshell) ! max basis in shell (output)
double precision frc_2el(3,nat)
integer nat
integer lbuf
double precision buf(lbuf)
integer lscr
double precision scr(lscr)
c local
integer ish, jsh, ksh, lsh
integer ilo, ihi, jlo, jhi, klo, khi, llo, lhi
integer i, j, k, l, ij, kl
integer ijkl, mx_d2, sh_ijkl, ijsh, klsh, sh_ijkl_len
integer ij_max, kl_min, mxp
integer sizeof_int, szlabs
integer l_sq, l_d2, l_vb, l_lb, l_buf, l_dbg
integer k_sq, k_d2, k_vb, k_lb, k_buf, k_dbg
integer cijkl
logical status
integer nipv
parameter(nipv=4)
integer s, t, u, v, st, uv, stuv, shoff, sqindx
integer nsp, is, js, ks, ls, a, b, c, d, ab, cd, abcd
integer imin, jmin, kmin, lmin, imax, jmax, kmax, lmax
integer slen, tlen, ulen, vlen, stlen
integer i8n2max, i8l2rec
integer avail
integer iln, kln
cgk debug
* write(*,*)'gk: entered sif2da_2e_sort'
cgk end
c maximimum number of unique density matrix elements
mx_d2 = (nbft*(nbft+1)+2)*nbft*(nbft+1)/8
c maximum pairs
mxp = (nbft*(nbft+1))/2
nsp=((nshell+1)*nshell)/2
cgk debug
* write(*,*)'gk: nsp = ', nsp
cgk end
c generate local data to avoid overhead
do lsh=1,nshell
status=bas_cn2bfr(ibas,lsh,llo,lhi)
do l=llo,lhi
b2s(l)=lsh
enddo
sminb(lsh)=llo
smaxb(lsh)=lhi
enddo
cgk debug
* write(*,*)'gk: nbft = ', nbft
* write(*,*)'gk: b2s:'
* do l=1,nbft
* write(*,*)l, b2s(l)
* enddo
* write(*,*)'gk: sminb, smaxb:'
* do lsh=1,nshell
* write(*,*)lsh, sminb(lsh), smaxb(lsh)
* enddo
cgk end
c Get memory for shell block off set. Diminsioned for
c number of shell pairs by number of shell pairs
if (.not. ma_push_get(mt_int,nsp*nsp,'num_shell_quart',
$ l_sq, k_sq)) call errquit
& ('sif2da_2e_sort;could not allocate l_sq',nsp*nsp,
& MA_ERR)
cgk debug
c Get memory for debuggin array
if (.not. ma_push_get(mt_int,4*mx_d2,'debug_array',
$ l_dbg, k_dbg)) call errquit
& ('sif2da_2e_sort;could not allocate dbg array',nsp*nsp,
& MA_ERR)
cgk end
call ifill(4*mx_d2,0,int_mb(k_dbg),1)
call ifill(nsp*nsp, 0, int_mb(k_sq), 1)
c go order similar to aodens and figure shell quartet offests
ijkl = 0
sh_ijkl=0
sh_ijkl_len=0
shoff=0
do ksh=1,nshell
klo=sminb(ksh)
khi=smaxb(ksh)
kln=khi-klo+1
do lsh=1,ksh
* klsh=nshell*(lsh-1)-((lsh-1)*lsh)/2+ksh
klsh=((ksh-1)*ksh)/2+lsh
llo=sminb(lsh)
lhi=smaxb(lsh)
do ish=ksh,nshell
ilo=sminb(ish)
ihi=smaxb(ish)
iln=ihi-ilo+1
do jsh=1,ish
* ijsh=nshell*(jsh-1)-((jsh-1)*jsh)/2+ish
ijsh=((ish-1)*ish)/2+jsh
jlo=sminb(jsh)
jhi=smaxb(jsh)
if(ijsh.ge.klsh) then
sh_ijkl=sh_ijkl+1
sqindx=(klsh-1)*nsp+ijsh
int_mb(k_sq-1+sqindx)=shoff
cgk debug
* write(6,2323) '(',ish,',',jsh,';',ksh,',',lsh,')'
* write(6,'(2(a,i4))')'st=',ijsh, ' uv=', klsh
* write(6,'(2(a,i5))') 'sqindx = ', sqindx,
* & ' shoff=', int_mb(k_sq-1+sqindx)
cgk end
sh_ijkl_len=0
do k=klo,khi
if(lsh.eq.ksh) lhi=k
do l=llo,lhi
kl=((k-1)*k)/2+l
do i=ilo,ihi
if(jsh.eq.ish) jhi=i
do j=jlo,jhi
ij=((i-1)*i)/2+j
if(ijsh.eq.klsh) then
if(ij.ge.kl) then
sh_ijkl_len=sh_ijkl_len+1
ijkl=ijkl+1
cgk debug
int_mb(k_dbg-1+4*(ijkl-1)+1)=i
int_mb(k_dbg-1+4*(ijkl-1)+2)=j
int_mb(k_dbg-1+4*(ijkl-1)+3)=k
int_mb(k_dbg-1+4*(ijkl-1)+4)=l
cgk end
cgk debug
* write(6,'(a,i4)') '(ijkl)=', ijkl
* write(6,2323) '(',ish,',',jsh,';',ksh,',',lsh,')'
* write(6,2323) '[',i,',',j,';',k,',',l,']'
cgk end
endif
else
sh_ijkl_len=sh_ijkl_len+1
ijkl=ijkl+1
cgk debug
int_mb(k_dbg-1+4*(ijkl-1)+1)=i
int_mb(k_dbg-1+4*(ijkl-1)+2)=j
int_mb(k_dbg-1+4*(ijkl-1)+3)=k
int_mb(k_dbg-1+4*(ijkl-1)+4)=l
cgk end
cgk debug
2322 format(2(a,i2))
* write(6,'(a,i4)') '(ijkl)=', ijkl
* write(6,2323) '(',ish,',',jsh,';',ksh,',',lsh,')'
* write(6,2323) '[',i,',',j,';',k,',',l,']'
2323 format(4(a1,i2),a1)
endif
enddo
enddo
enddo
enddo
shoff=shoff+sh_ijkl_len
* write(*,*)'gk: shell length=',sh_ijkl_len
endif
enddo
enddo
enddo
enddo
cgk debug
* write(*,*)'gk: ijkl = ', ijkl
* write(*,*)'gk: mx_d2 = ', mx_d2
* call MA_summarize_allocated_blocks
* avail = MA_inquire_avail(mt_dbl)
* write(*,*)'gk: available mem in sif2da before d2 is allocated = ',
* & avail
cgk end
c memory for the *entire* 2 electron density array
if (.not. ma_push_get(mt_dbl, mx_d2,'d2 array',l_d2,
$ k_d2)) call errquit
& ('sif2da_2e_sort;could not allocate l_d2',mx_d2, MA_ERR)
c memory for the 2 electron sifs buffer
i8n2max=n2max
if (.not. ma_push_get(mt_dbl,i8n2max,'vals buf array',l_vb,
$ k_vb)) call errquit
& ('sif2da_2e_sort;could not allocate l_vb',i8n2max, MA_ERR)
c get memory for sifs 2e labels
sizeof_int=MA_sizeof(MT_INT, 1, MT_BYTE)
if(sizeof_int.eq.4) then
szlabs=(i8n2max+1)/2
elseif(sizeof_int.eq.8) then
szlabs=i8n2max
else
call errquit('wrt_dft_aoints: unsupporte integer type size=',
& sizeof_int, MA_ERR)
endif
if (.not. MA_push_get(MT_Int,nipv*szlabs,'wrt_dft_aoints:dlab',
& l_lb, k_lb))
$ call errquit('wrt_dft_aoints: ma failed for sif labs',
& nipv*szlabs, MA_ERR)
c memory for the 2 electron sifs buffer
i8l2rec=l2rec
if (.not. ma_push_get(mt_dbl,i8l2rec,'vals buf array',l_buf,
$ k_buf)) call errquit
& ('sif2da_2e_sort;could not allocate l_buf',l2rec, MA_ERR)
c
c generate array of integrals from sifs file
c
cgk debug
* write(*,*)'gk: calling sif2arr'
* open(unit=26,form='formatted')
cgk end
call sif2arr(aodens, dbl_mb(k_d2), mx_d2, dbl_mb(k_vb), nipv,
& int_mb(k_lb), dbl_mb(k_buf), info, b2s, sminb, smaxb,
& int_mb(k_sq), nsp, nshell, int_mb(k_dbg))
cgk debug
* write(*,*)'gk: back from sif2arr'
* close(26)
cgk end
c dump sif 2e buffers
if (.not. MA_pop_stack(l_buf)) call errquit
$ ('sif2da_2e: pop failed at l_buf', 0, GA_ERR)
if (.not. MA_pop_stack(l_lb)) call errquit
$ ('sif2da_2e: pop failed at l_lb', 0, GA_ERR)
if (.not. MA_pop_stack(l_vb)) call errquit
$ ('sif2da_2e: pop failed at l_vb', 0, GA_ERR)
c get memory for the gradient integrals
cgk debug
* write(*,*)'gk: in sif2da, lbuf = ', lbuf
* write(*,*)'gk: in sif2da, lscr = ', lscr
cgk end
call d2geri_trace(mx_d2, dbl_mb(k_d2), lscr, scr, lbuf, buf,
& nbft, b2s, nshell, sminb, smaxb, nsp, int_mb(k_sq), nat,
& frc_2el, ibas, int_mb(k_dbg))
cgk debug
* write(*,*)'gk: in sif2da_2e_sort back for d2geri_trace'
cgk end
if (.not. MA_pop_stack(l_d2)) call errquit
$ ('sif2da_2e: pop failed at l_d2', 0, GA_ERR)
cgk debug
if (.not. MA_pop_stack(l_dbg)) call errquit
$ ('sif2da_2e: pop failed at l_dbg', 0, GA_ERR)
cgk end
if (.not. MA_pop_stack(l_sq)) call errquit
$ ('sif2da_2e: pop failed at l_sq', 0, GA_ERR)
return
end

84
src/NWints/aoints/sif2ga.F Normal file
View file

@ -0,0 +1,84 @@
subroutine sif2ga(g,locar,basis,nbft)
#include "global.fh"
#include "bas.fh"
#include "nwc_const.fh"
#include "errquit.fh"
#include "mafdecls.fh"
integer g ! atom blocked global array
double precision locar(*) ! local array from sifs with all values
integer basis
integer nbft
integer nshell
integer ishell, jshell
integer ijshell, ilo, ihi, jlo, jhi, blkdim
integer max1e, mem1
integer l_buf
integer k_buf
integer iloc, i, j, ijmap, il, jl
logical shells_ok
logical int_chk_sh
external int_chk_sh
cgk debug
* write(*,*)'gk: entered sif2ga'
cgk end
if (.not. bas_numcont(basis, nshell)) call errquit
& ('sif2ga: bas_numcont failed for basis', basis,
& BASIS_ERR)
call int_mem_1e(max1e, mem1)
cgk debug
* write(*,*)'gk: max1e = ', max1e
cgk end
if (.not. MA_push_get(MT_DBL,max1e,'sif2ga:buf',l_buf,k_buf))
$ call errquit('sif2ga: ma failed', max1e, MA_ERR)
call ga_zero(g)
cgk debug
* write(*,*)'gk: nshell = ', nshell
cgk end
ijmap = 0
do jshell = 1, nshell
* write(*,*)'gk: jshell = ', jshell
if (.not. bas_cn2bfr(basis, jshell, jlo, jhi))
$ call errquit('sif2ga: bas_cn2bfr jshell', jshell,
& BASIS_ERR)
* write(*,*)'gk: jlo = ', jlo, ' jhi = ', jhi
do ishell = 1, nshell
* write(*,*)'gk: ishell = ', ishell
if (.not. bas_cn2bfr(basis, ishell, ilo, ihi))
$ call errquit('sif2ga: bas_cn2bfr ishell', ishell,
& BASIS_ERR)
* write(*,*)'gk: ilo = ', ilo, ' ihi = ', ihi
call dfill(max1e, 0.0d0, dbl_mb(k_buf), 1)
blkdim=ihi-ilo+1
iloc=0
do j = jlo, jhi
do i = ilo, ihi
jl=max(i,j)
il=min(i,j)
ijmap=(jl-1)*jl/2+il
* write(*,2000) 'gk: j=',j,' i=',i,' ijmap=',ijmap,
* & ' iloc=',iloc
2000 format(4(a,i4))
* write(*,2010) 'gk: ',locar(ijmap)
2010 format(a,f20.12)
dbl_mb(k_buf+iloc)=locar(ijmap)
iloc=iloc+1
enddo
enddo
call ga_put(g,ilo,ihi,jlo,jhi,dbl_mb(k_buf),blkdim)
enddo
enddo
if (.not. MA_pop_stack(l_buf)) call errquit
$ ('sif2ga: pop failed', 0, GA_ERR)
return
end

113
src/NWints/aoints/sifs_2e_task.F Normal file
View file

@ -0,0 +1,113 @@
cgk be sure int_acc_set is called properly.
subroutine sifs_2e_task( geom, basis, oskel, iatlo, jatlo,
& katlo, latlo, iathi, jathi, kathi, lathi, ijk_prev,
$ tmp, tol2e, aoint2, ninfo, info, ifmt, ibvtyp, values,
& ibitv, buffer, clabs, ilabs, jlabs, klabs, llabs, eris )
implicit none
#include "errquit.fh"
#include "schwarz.fh"
#include "bas.fh"
#include "cfock.fh"
#include "cscfps.fh"
#include "sym.fh"
c
c
integer geom
integer basis
logical oskel
integer iatlo, jatlo, katlo, latlo
integer iathi, jathi, kathi, lathi
integer ijk_prev(3,2)
double precision tmp(*)
double precision tol2e
integer ilabs(*), jlabs(*), klabs(*), llabs(*)
double precision eris(*)
c
c Given a task which comprises computing the contributions
c arising from a quartet of blocks of atoms, do everything necessary
c
double precision sij, sijkl, q4, smax
logical otest, otestij
integer ilo, ihi, jlo, jhi, klo, khi, llo, lhi
integer iat, jat, kat, lat, lattop, kattop, jattop
logical oij, okl, oikjl
integer junk
logical int2e_set_bf_range
external int2e_set_bf_range
integer*4 clabs(4,*)
integer*4 aoint2, ninfo, ifmt, ibvtyp, ibitv
integer*4 info(ninfo)
double precision values(*), buffer(*)
cgk debug
* write(*,*)'gk: entered sifs_2e_task'
* write(*,*)'gk: oskel = ', oskel
cgk end
c
c With caching get density and flush necessary fock blocks
c and return magnitude of density blocks for screening
c
if (.not. bas_ce2bfr(basis, iatlo, ilo, junk))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
if (.not. bas_ce2bfr(basis, iathi, junk, ihi))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
c
if (.not. bas_ce2bfr(basis, jatlo, jlo, junk))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
if (.not. bas_ce2bfr(basis, jathi, junk, jhi))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
c
if (.not. bas_ce2bfr(basis, katlo, klo, junk))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
if (.not. bas_ce2bfr(basis, kathi, junk, khi))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
c
if (.not. bas_ce2bfr(basis, latlo, llo, junk))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
if (.not. bas_ce2bfr(basis, lathi, junk, lhi))
$ call errquit('fock_2e_task: ce2bfr failed', 0, BASIS_ERR)
c
if (owritefile) then
if (.not. int2e_set_bf_range(
$ ilo, ihi, jlo, jhi, klo, khi, llo, lhi))
$ call errquit('should not happen', 0, UNKNOWN_ERR)
end if
c
q4 = 1.0d0
oij = iatlo.eq.jatlo
okl = katlo.eq.latlo
oikjl = (iatlo.eq.katlo) .and. (jatlo.eq.latlo)
c
do iat = iatlo, iathi
jattop = jathi
if (oij) jattop = iat
do jat = jatlo, jattop
kattop = kathi
if (oikjl) kattop = iat
do kat = katlo, kattop
lattop = lathi
if (okl) lattop = kat
if (oikjl .and. kat.eq.iat) lattop = jat
do lat = latlo,lattop
if (oskel) otest = sym_atom_quartet
$ (geom, iat, jat, kat, lat, q4)
* if (otest) then
cgk debug
* write(*,'(a,4i3)')'gk: iat,jat,kat,lat: ',iat,jat,kat,lat
cgk end
call int_2e_sifs_b(basis, tol2e, q4, iat, jat,
& kat, lat, ilo, jlo, klo, llo, ihi, jhi, khi,
& lhi, aoint2, ninfo, info, ifmt, ibvtyp,
& values, ibitv, buffer, clabs,
& ilabs, jlabs, klabs, llabs, eris)
* end if
end do
end do
end do
end do
cgk debug
* write(*,*)'gk: leaving sifs_2e_task'
cgk end
c
end

315
src/NWints/aoints/sopgrdtrc.F Normal file
View file

@ -0,0 +1,315 @@
subroutine sopgrdtrc
$ ( H, lbuf, scr, lscr, wdens,
$ frc_sow,
$ g_dens, g_wdens, basis, geom, nproc, nat,
$ max_at_bf, oskel,
& frc_sox, frc_soy, frc_soz, densx, densy, densz)
c based on grad1_so in nwdft/sodft -gsk
c keep the parallel stuff. Sure why not?
C spin orbit potential gradient using effective densities from
c columbus
implicit none
#include "mafdecls.fh"
#include "global.fh"
#include "geom.fh"
#include "bas.fh"
#include "rtdb.fh"
#include "sym.fh"
C-------------------------parameters--------------------------------
integer lbuf, lscr,
$ g_dens(3), ! density matrix (summed if ROHF, UHF)
$ g_wdens, ! weighted density (Lagrangian)
$ basis, geom, nproc, nat, max_at_bf
double precision H, ! integral derivatives
$ scr,
$ densx, densy, densz, ! local density block
$ wdens, ! local weighted density block
$ frc_sow, frc_sox, frc_soy, frc_soz ! forces arrays
dimension H ( lbuf ),
$ frc_sow(3, nat), frc_sox(3, nat), frc_soy(3,nat),
$ frc_soz(3,nat), scr(lscr),
$ wdens(max_at_bf,max_at_bf),
$ densx(max_at_bf,max_at_bf), densy(max_at_bf,max_at_bf),
$ densz(max_at_bf,max_at_bf)
logical oskel ! symmetry?
C-------------------------local variables--------------------------
integer ijatom, next, iat1, iat2, iat3, ish1, ish2,
$ iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l,
$ if1, il1, if2, il2,
$ icart, ic, ip1, ip2
integer i, j
double precision crd1, crd2 ! atomic coordinates
dimension crd1(3), crd2(3)
integer idatom
dimension idatom(2)
double precision dE, dx, dy, dz, qfac, fact, q1, q2
logical status, pointforce
character*16 name
integer nxtask, task_size
external nxtask
cgk debug
character*2 deriv(3)
data deriv/'dx','dy','dz'/
cgk end
task_size = 1
status = rtdb_parallel(.true.) ! Broadcast reads to all processes
pointforce = geom_include_bqbq(geom)
call hf_print_set(1)
do j=1, nat
do i=1,3
frc_sox(i,j)=0d0
frc_soy(i,j)=0d0
frc_soz(i,j)=0d0
frc_sow(i,j)=0d0
enddo
enddo
cgk debug
* write(*,*)'gk: printing so blocks'
cgk end
ijatom = -1
next = nxtask(nproc,task_size)
do iat1=1, nat
do iat2 = 1, nat
ijatom = ijatom + 1
if ( ijatom .eq. next ) then
status = bas_ce2bfr(basis,iat1,iab1f,iab1l)
status = bas_ce2bfr(basis,iat2,iab2f,iab2l)
if (iab1f.le.0 .or. iab2f.le.0) goto 999
status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
call ga_get (g_dens(1),
& iab1f,iab1l,iab2f,iab2l,densx,max_at_bf)
call ga_get (g_dens(2),
& iab1f,iab1l,iab2f,iab2l,densy,max_at_bf)
call ga_get (g_dens(3),
& iab1f,iab1l,iab2f,iab2l,densz,max_at_bf)
do ish1 = iac1f, iac1l
do ish2 = iac2f, iac2l
status = bas_cn2bfr(basis,ish1,if1,il1)
status = bas_cn2bfr(basis,ish2,if2,il2)
cgk debug
* write(*,'(a,2i4)')'gk: int block for shells, ', ish1, ish2
cgk end
call intd_1eso(basis,ish1,basis,ish2,lscr,scr,
& lbuf,H)
do iat3 = 1, nat
do icart = 1, 3
* call print_soblock(H,if1,il1,if2,il2,iat3,icart,1,
* & nat)
call trc_soblock(H,if1,il1,if2,il2,iat3,icart,1,
& nat,iab1f,iab2f,max_at_bf,frc_soz,densz)
* call print_soblock(H,if1,il1,if2,il2,iat3,icart,2,
* & nat)
call trc_soblock(H,if1,il1,if2,il2,iat3,icart,2,
& nat,iab1f,iab2f,max_at_bf,frc_soy,densy)
* call print_soblock(H,if1,il1,if2,il2,iat3,icart,3,
* & nat)
call trc_soblock(H,if1,il1,if2,il2,iat3,icart,3,
& nat,iab1f,iab2f,max_at_bf,frc_sox,densx)
enddo
enddo
enddo
enddo
999 continue
next = nxtask(nproc,task_size)
endif
enddo
enddo
* write(*,'("forces sow",9e13.5)')((frc_sow(i,j),i=1,3),j=1,nat)
* write(*,'("forces sox",9e13.5)')((frc_sox(i,j),i=1,3),j=1,nat)
* write(*,'("forces soy",9e13.5)')((frc_soy(i,j),i=1,3),j=1,nat)
* write(*,'("forces soz",9e13.5)')((frc_soz(i,j),i=1,3),j=1,nat)
next = nxtask(-nproc, task_size)
do j=1, nat
do i=1,3
* frc_sox(i,j)=0d0
* frc_soy(i,j)=0d0
* frc_soz(i,j)=0d0
frc_sow(i,j)=0d0
enddo
enddo
ijatom = -1
next = nxtask(nproc,task_size)
do 90, iat1 = 1, nat
do 80, iat2 = 1, iat1
ijatom = ijatom + 1
if ( ijatom .eq. next ) then
status = bas_ce2bfr(basis,iat1,iab1f,iab1l)
status = bas_ce2bfr(basis,iat2,iab2f,iab2l)
if (iab1f.le.0 .or. iab2f.le.0) then
c
c At least one center has no functions on it ... next atom
c
goto 1010
endif
if (oskel) then
if (.not. sym_atom_pair(geom, iat1, iat2, qfac))
$ goto 1010
else
qfac = 1.0d0
endif
status = bas_ce2cnr(basis,iat1,iac1f,iac1l)
status = bas_ce2cnr(basis,iat2,iac2f,iac2l)
call ga_get(g_wdens,
& iab1f,iab1l,iab2f,iab2l,wdens,max_at_bf)
* call ga_get (g_dens(1),
* & iab1f,iab1l,iab2f,iab2l,densx,max_at_bf)
* call ga_get (g_dens(2),
* & iab1f,iab1l,iab2f,iab2l,densy,max_at_bf)
* call ga_get (g_dens(3),
* & iab1f,iab1l,iab2f,iab2l,densz,max_at_bf)
do 70, ish1 = iac1f, iac1l
if ( iat1.eq.iat2 ) iac2l = ish1
do 60, ish2 = iac2f, iac2l
C shell block in atomic (D/Dw)-matrix block
status = bas_cn2bfr(basis,ish1,if1,il1)
if1 = if1 - iab1f + 1
il1 = il1 - iab1f + 1
status = bas_cn2bfr(basis,ish2,if2,il2)
if2 = if2 - iab2f + 1
il2 = il2 - iab2f + 1
C overlap derivatives
call intd_1eov(basis,ish1,basis,ish2,lscr,scr,
& lbuf,H,idatom)
C Dsow x S
if ( idatom(1) .ge. 1 ) then
C idatom(1).ge.0 <=> idatom(2).ge.0 (no check necessary)
ic = 1
do 28, icart = 1, 3
de = 0.D0
do 22, ip1 = if1, il1
do 20, ip2 = if2, il2
dE = dE + wdens(ip1,ip2) * H(ic)
ic = ic + 1
20 continue
22 continue
dE = dE * qfac
frc_sow(icart,idatom(1)) = frc_sow(icart,idatom(1))
$ - dE - dE
frc_sow(icart,idatom(2)) = frc_sow(icart,idatom(2))
$ + dE + dE
28 continue
endif
C 1el. so. derivatives
* call intd_1eso(basis,ish1,basis,ish2,lscr,scr,
* & lbuf,H)
C Dso x Hso
* ic=1
* do 150, iat3 = 1, nat
* do 140, icart = 1, 3
c z componet
cgk debug
* write(*,*)'gk: z so trace'
* write(*,2010)'deriv: ',iat3,deriv(icart)
2010 format(a,i1,x,a2)
cgk end
* dE = 0.D0
* do 131, ip1 = if1, il1
* do 128, ip2 = if2, il2
cgk debug
* write(*,2000)'i=',ip1+iab1f-1,' j=',ip2+iab2f-1,' den=',
* & densz(ip1,ip2),' gint=',H(ic)
2000 format(a,i3,a,i3,a,f20.10,a,f20.10)
cgk end
* dE = dE + densz(ip1,ip2)*H(ic)
* ic = ic + 1
*128 continue
*131 continue
* if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
* dE = dE * qfac
* frc_soz(icart,iat3) = frc_soz(icart,iat3) + dE
c y componet
cgk debug
* write(*,*)'gk: y so trace'
* write(*,2010)'deriv: ',iat3,deriv(icart)
cgk end
* dE = 0.D0
* do 230, ip1 = if1, il1
* do 231, ip2 = if2, il2
cgk debug
* write(*,2000)'i=',ip1+iab1f-1,' j=',ip2+iab2f-1,' den=',
* & densz(ip1,ip2),' gint=',H(ic)
cgk end
* dE = dE + densy(ip1,ip2)*H(ic)
* ic = ic + 1
*231 continue
*230 continue
* if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
* dE = dE * qfac
* frc_soy(icart,iat3) = frc_soy(icart,iat3) + dE
c x component
cgk debug
* write(*,*)'gk: x so trace'
* write(*,2010)'deriv: ',iat3,deriv(icart)
cgk end
* dE = 0.D0
* do 250, ip1 = if1, il1
* do 251, ip2 = if2, il2
cgk debug
* write(*,2000)'i=',ip1+iab1f-1,' j=',ip2+iab2f-1,' den=',
* & densz(ip1,ip2),' gint=',H(ic)
cgk end
* dE = dE + densx(ip1,ip2)*H(ic)
* ic = ic + 1
*251 continue
*250 continue
* if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE
* dE = dE * qfac
* frc_sox(icart,iat3) = frc_sox(icart,iat3) + dE
*140 continue
*150 continue
60 continue
70 continue
1010 continue
next = nxtask(nproc,task_size)
endif
80 continue
90 continue
next = nxtask(-nproc,task_size)
write(*,'("forces sow",9e13.5)')((frc_sow(i,j),i=1,3),j=1,nat)
write(*,'("forces sox",9e13.5)')((frc_sox(i,j),i=1,3),j=1,nat)
write(*,'("forces soy",9e13.5)')((frc_soy(i,j),i=1,3),j=1,nat)
write(*,'("forces soz",9e13.5)')((frc_soz(i,j),i=1,3),j=1,nat)
return
end

28
src/NWints/aoints/trc_soblock.F Normal file
View file

@ -0,0 +1,28 @@
subroutine trc_soblock(H, ilo, ihi, jlo, jhi, iat, icart,
& socomp, nat, iabf, jabf, max_at_bf, frc, dblk)
integer ilo, ihi, jlo, jhi, iat, icart, socomp, nat, iabf, jabf,
& max_at_bf
double precision H(jlo:jhi, ilo:ihi, 3, 3, nat)
double precision frc(3,nat)
double precision dblk(max_at_bf,max_at_bf)
character*2 dcart(3)
data dcart/'dx','dy','dz'/
integer i, j, ip, jp
double precision de
de = 0d0
do i=ilo, ihi
ip=i-iabf+1
do j=jlo, jhi
jp=j-jabf+1
de=de+H(j,i,socomp,icart,iat)*dblk(ip,jp)
* write(6,2000) iat,dcart(icart),'<',i,'|Hso',socomp,'|',j,'>=',
* & H(j,i,socomp,icart,iat), ' den=',dblk(ip,jp)
2000 format(i2,a2,3(a,i2),2(a,f20.10))
enddo
enddo
frc(icart,iat)=frc(icart,iat)+de
* write(*,'("forces",i2,9e13.5)')socomp,((frc(i,j),i=1,3),j=1,nat)
return
end

562
src/NWints/aoints/wrt_dft_aoints.F Normal file
View file

@ -0,0 +1,562 @@
c
c /* Main DFT driver */
c
logical function wrt_dft_aoints(rtdb)
c
c $Id:
c
implicit none
#include "errquit.fh"
c
integer rtdb ! rtdb handle [input]
Logical IOLGC, oprint, oprint_time, oprint_ioinfo
c
#include "bas.fh"
#include "geom.fh"
#include "rtdb.fh"
#include "mafdecls.fh"
#include "stdio.fh"
#include "cdft.fh"
#include "global.fh"
#include "util.fh"
#include "tcgmsg.fh"
#include "schwarz.fh"
#include "msgids.fh"
#include "dftps.fh"
#include "bgj.fh"
#include "dft_scaleig.fh"
!#ifdef SCALAPACK
! integer ga_cholesky,ga_llt_i
! external ga_cholesky,ga_llt_i
!#endif
c
#include "zora.fh"
#include "rel_consts.fh"
c
#include "case.fh"
#include "geomP.fh"
c
integer ga_create_atom_blocked
external ga_create_atom_blocked
c
integer g_xcinv, g_oi
* integer g_2ceri, gcdinv
integer g_tmpb, g_tmp
integer nvxc
integer me, nproc, cinfo
integer iter, nconstr
integer iVxc_opt, lmiss, ivcoul_opt
integer lsvals, isvals, idum,n_dep
double precision toll_s
logical dft_scf,atmblk, lcdft
external dft_scf
logical dft_scf_so
external dft_scf_so
logical status
logical ldirect
c sifs
integer sizeof_int, nmap, l_map, k_map, szintarr, szlabs
integer k_lab
integer l_lab
integer nmax
parameter(nmap=2)
integer*4 aoints
parameter(aoints=23)
integer*4 ifmt1, ibvtyp, n1max, ierr, ifmt, n2max
INTEGER l1rec, l2rec
integer lrecal
integer*4 ibitv
integer*4 mxengy
parameter (mxengy=1)
real*8 energy(mxengy)
integer*4 ninfo, iunit2, aoint2
parameter(ninfo=5)
integer*4 info(ninfo) ! not to be confused with NWChem info
real*8 tol2e
parameter(tol2e=1d-10) ! integrate with user parameters...
cgk debug
integer mem1, max1e, mem2, max2e
cgk end
ierr=0 !initialize
wrt_dft_aoints = .false.
c
me = ga_nodeid()
nproc = ga_nnodes()
oprint = util_print('dft cinfo', print_default)
oprint_time = util_print('dft timings', print_high)
oprint_ioinfo = util_print('io info', print_debug)
if(ga_nodeid().eq.0.and.
. util_print('geometry', print_high)) then
if(.not.geom_print(geom)) call errquit(
. ' wrt_dft_aoints: geom print failed',0, GEOM_ERR)
endif
IOLGC = .TRUE.
if (noio.eq.1) IOLGC = .FALSE.
!#ifdef SCALAPACK
!#if defined(PARALLEL_DIAG)
! call slexit
! call slexit2
! if (.not. rtdb_cget(rtdb, 'dft:scaladiag',1,yev))
! & yev='d'
! if (.not. rtdb_get(rtdb, 'dft:scalanb',mt_int,1,nb))
! & nb=0
!#endif
!#endif
c
c Determine the characteristics of the AO Gaussian basis set.
c This uses the old NWdft data structures, rather than tthe rtdb
c
if (.not. bas_numbf(AO_bas_han, nbf_ao)) then
call errquit('Exiting from wrt_dft_aoints',1, BASIS_ERR)
endif
* iter = 0
if(nproc.gt.1) then
write(6,'(a,a)')'WARNING: Writing aoints, calculating 1e ints ',
& 'with 1 CPU'
write(6,'(a,a)')'WARNING: we are currently supporting 1 ',
& 'process for aoints'
endif
c
c open the sifs file
c
if(ga_nodeid().eq.0) then
open(unit=aoints,status='unknown',form='unformatted',
& file='aoints')
rewind aoints
* open(unit=23,status='unknown',form='formatted',
* & file='testfile')
* write(23,*) 'testing testfile'
endif
c
c Check zora options
do_zora_get_correction=.false.
if (rtdb_get(rtdb,'zora:calc',mt_log,1,
& do_zora_get_correction))
& do_zora_get_correction=.true.
c
do_zora = .false.
if (rtdb_get(rtdb,'zora:on',mt_log,1,do_zora))
& do_zora = .true.
c
c == Get the grid cutoff or set default ==
if(.not.rtdb_get(rtdb,'zora:cutoff',mt_dbl,1,zoracutoff))
& zoracutoff = 1.0d-08
c
c Get the light speed and other constants, if specified
call rel_init_con(rtdb)
clight_au = cau
c
*c Get and check the Coulomb Attenuation Method (CAM) parameters
* if (.not.rtdb_get(rtdb, 'dft:cam_exch', mt_log, 1,
* & cam_exch)) cam_exch=.false.
* if (.not.rtdb_get(rtdb, 'dft:cam_omega', mt_dbl, 1,
* & cam_omega)) cam_omega=0.d0
* if (.not.rtdb_get(rtdb, 'dft:cam_alpha', mt_dbl, 1,
* & cam_alpha)) cam_alpha=0.d0
* if (.not.rtdb_get(rtdb, 'dft:cam_beta', mt_dbl, 1,
* & cam_beta)) cam_beta=0.d0
*c
* ldirect = .false.
* if (cam_exch .and. (.not.CDFIT) ) then
* if (rtdb_get(rtdb, 'dft:direct', mt_log, 1, ldirect)) then
* if (.not.ldirect) then
* if (ga_nodeid().eq.0) then
* write(luout,*) "CAM type calculations without CDFIT are only"
* write(luout,*) "available with the direct option"
* end if
* call errquit('Please use the direct option',4, INPUT_ERR)
* end if
* end if
* end if
c
c Select method to compute the XC potential.
if (XC_bas_han.eq.-99999)then
iVxc_opt = 0
XCFIT = .FALSE.
else
call errquit('wrt_dft_aoints: XC basis fitting not implemented'
&, 0,CAPMIS_ERR)
endif
* if (XCFIT)then
*c
*c Determine the characteristics of the XC Gaussian basis set.
*c
* if (.not. bas_numbf(XC_bas_han, nbf_xc)) then
* call errquit('Exiting from scfdriver',4, BASIS_ERR)
* endif
*c
*c Compute the matrix inverse of the 2-ctr OIs.
*c
**ga:1:0
* if (.not. ga_create(mt_dbl, nbf_xc, nbf_xc, 'XCinv',
* & nbf_xc, 0, g_xcinv))
* & call errquit('wrt_dft_aoints: error creating g_xcinv',
* & 0,GA_ERR)
* call ga_zero(g_xcinv)
**ga:1:0
* if (.not. ga_create(mt_dbl, nbf_xc, nbf_xc, 'g_oi',
* & nbf_xc, 0, g_oi))
* & call errquit('wrt_dft_aoints: error creating g_oi',0,
* & GA_ERR)
* call ga_zero(g_oi)
* call int_1e_ga(xc_bas_han, xc_bas_han,
* . g_oi, 'overlap', oskel)
* if (oskel)call sym_symmetrize(
* ( geom, xc_bas_han, .false., g_oi)
*#if defined(PARALLEL_DIAG)
* call ga_chol(g_oi, g_xcinv, cinfo)
*#else
* call ga_chol_seq(g_oi, g_xcinv, cinfo)
*#endif
* if (cinfo.ne.0)then
* if (me.eq.0)then
* write(LuOut,*)' Problem in performing a Choleski '
* write(LuOut,*)' decomposition of the 2-ctr overlap '
* write(LuOut,*)' matrix using XC fitting basis. '
* write(LuOut,*)' Attempting a diag/inverse. '
* endif
* endif
* if (cinfo.eq.0) then
*#if defined(PARALLEL_DIAG)
* call ga_inverse(g_xcinv, g_oi)
* if (odftps) call pstat_on(ps_bdgemm)
* call ga_dgemm('T', 'N', nbf_xc, nbf_xc, nbf_xc, 1.d0,
* & g_oi, g_oi, 0.d0, g_xcinv)
* if (odftps) call pstat_off(ps_bdgemm)
*#else
* call ga_copy(g_xcinv, g_oi)
* call ga_inv_seq(g_oi, g_xcinv)
*#endif
* else
* call dft_invdiag(g_oi, g_xcinv, nbf_xc)
* endif
* if (.not. ga_destroy(g_oi)) call errquit
* & ('wrt_dft_aoints: could not destroy g_oi', 0, GA_ERR)
* if (me.eq.0.and.IOLGC) then
* lmiss = 0
* call dft_invio('XC', g_xcinv, nbf_xc, 'WRITE', lmiss)
* if (lmiss.ne.0)call errquit
* & (' wrt_dft_aoints; dft_invio - abnormal write of XC ', 0,
* & DISK_ERR)
* endif
* if (IOLGC)then
* if (.not. ga_destroy(g_xcinv)) call errquit
* & ('wrt_dft_aoints: could not destroy g_xcinv', 0, GA_ERR)
* endif
* endif
*c
c Done with XC setup stuff....
c Select method to compute Coulomb potential.
c
if (CD_bas_han.eq.-99999)then
iVcoul_opt = 0
CDFIT = .FALSE.
else
call errquit('wrt_dft_aoints: charge density fitting not implem
&ented', 0,CAPMIS_ERR)
endif
c
c
* if (CDFIT)then
*c
*c Determine the characteristics of the CD Gaussian basis set.
*c
* if (odftps) call pstat_on(ps_2cer)
* if (.not. bas_numbf(CD_bas_han, nbf_cd)) then
* call errquit('Exiting from scfdriver',5, BASIS_ERR)
* endif
*c
*c Compute the matrix inverse of the CD 2-ctr ERIs.
*c
* if (.not. ga_create(mt_dbl, nbf_cd, nbf_cd, 'CD 2cERI',
* & 0, nbf_cd, g_2ceri))
* & call errquit('wrt_dft_aoints: error creating g_2ceri',0,
* & GA_ERR)
* call ga_zero(g_2ceri)
* call dft_get2eri(CD_bas_han, g_2ceri,oskel)
* if (odftps) call pstat_off(ps_2cer)
* if (oskel)call
* . sym_symmetrize(geom,cd_bas_han,.false.,g_2ceri)
* call ga_sync()
* if (.not. ga_duplicate(g_2ceri, g_cdinv, 'CD 2cERInv'))
* & call errquit('wrt_dft_aoints: error creating g_cdinv',0,
* & GA_ERR)
* if (odftps) call pstat_on(ps_inv)
* call ga_zero(g_cdinv)
* lmiss = 1
*c if (lmiss.eq.1) then
* call ga_zero(g_cdinv)
* cinfo = 0
*#if defined(PARALLEL_DIAG)
*#ifdef SCALAPACK
* call ga_copy(g_2ceri, g_cdinv)
* call ga_sync()
* cinfo= ga_cholesky('U',g_cdinv)
*#else
* call ga_chol(g_2ceri, g_cdinv, cinfo)
*#endif
*#else
* call ga_chol_seq(g_2ceri, g_cdinv, cinfo)
*#endif
* if (cinfo.ne.0)then
* if (me.eq.0)then
* write(LuOut,*)' Problem in performing a Choleski '
* write(LuOut,*)' decomposition of the 2-ctr ERI '
* write(LuOut,*)' matrix using CD fitting basis. '
* write(LuOut,*)' Attempting a diag/inverse. '
* endif
* endif
* if (cinfo.eq.0) then
* g_tmpb = g_2ceri
*#if defined(PARALLEL_DIAG)
*#ifdef SCALAPACK
* cinfo = ga_llt_i('U',g_cdinv,-1)
* if (cinfo.ne.0)then
* if (me.eq.0)then
* write(LuOut,*)' Problem in performing a Invers. '
* write(LuOut,*)' of the 2-ctr ERI '
* endif
* call ga_sync
* call errquit('wrt_dft_aoints: inverse failed ',0,0)
* endif
*
*#else
* call ga_inverse(g_cdinv, g_tmpb)
*#endif
* if (odftps) call pstat_on(ps_bdgemm)
*#ifndef SCALAPACK
* call ga_dgemm('T', 'N', nbf_cd, nbf_cd, nbf_cd, 1.d0,
* & g_tmpb, g_tmpb, 0.d0, g_cdinv)
*#endif
* if (odftps) call pstat_off(ps_bdgemm)
*#else
* call ga_copy(g_cdinv, g_tmpb)
* call ga_inv_seq(g_tmpb, g_cdinv)
*#endif
* else
* call dft_invdiag(g_2ceri, g_cdinv,
* & nbf_cd)
* endif
* if (odftps) call pstat_off(ps_inv)
*#ifndef SCALAPACK
*c
*c second build of g_2ceri needed becuase previous calls destroyed it
*c
* if (odftps) call pstat_on(ps_2cer)
* call ga_zero(g_2ceri)
* call dft_get2eri(CD_bas_han, g_2ceri,oskel)
* if (oskel)call
* . sym_symmetrize(geom,cd_bas_han,.false.,g_2ceri)
* if (odftps) call pstat_off(ps_2cer)
*#endif
* if (IOLGC.and.(me.eq.0)) then
* lmiss = 0
* call dft_invio('CDI', g_cdinv, nbf_cd, 'WRITE', lmiss)
* if (lmiss.ne.0)call errquit
* & (' wrt_dft_aoints; dft_invio - abnormal write of CDI ', 0
* & ,DISK_ERR)
* lmiss = 0
* call dft_invio('CD', g_2ceri, nbf_cd, 'WRITE', lmiss)
* if (lmiss.ne.0)call errquit
* & (' wrt_dft_aoints; dft_invio - abnormal write of CD ', 0,
* & DISK_ERR)
* endif
* if (oprint_ioinfo)then
* call ga_print(g_cdinv)
* call ga_print(g_2ceri)
* endif
*c endif
* if (IOLGC) then
* if (.not. ga_destroy(g_cdinv)) call errquit
* & ('wrt_dft_aoints: could not destroy g_xcinv', 0, GA_ERR)
* if (.not. ga_destroy(g_2ceri)) call errquit
* & ('wrt_dft_aoints: could not destroy g_xcinv', 0, GA_ERR)
* endif
* endif
*c
c Done with CD setup stuff
c
if (.not.MA_Push_Get(MT_Dbl, nbf_ao, 'ovl eig vals', lsvals,
& isvals))
& call errquit('wrt_dft_aoints: cannot allocate ovl eig vals', 0,
& MA_ERR)
call dfill(nbf_ao, 0.0d0, dbl_mb(isvals), 1)
c
c Set orbital overlap tolerance
c
call util_tolls(rtdb,.false.,toll_s,n_dep,dbl_mb(isvals),nbf_ao)
c
c set up SIFS buffer and other data for 1e ints.
c
ibvtyp=0
lrecal = 4096
C LB
c write(*,*)'LB: before sifcfg(1), ierr=',ierr
c WRITE(*,*)'n1max=',n1max,'n2max=',n2max,'ifmt1=',ifmt1,'ifmt',ifmt
c WRITE(*,*)"LB, sizeof l1rec=",sizeof(l1rec)
c WRITE(*,*)"LB, sizeof lrecal=",sizeof(lrecal)
C LB
call sifcfg(1, lrecal, nbf_ao, ibvtyp, ifmt1, l1rec, n1max, ierr)
lrecal = 4096
C LB
c write(*,*)'LB: after sifcfg(1), ierr=',ierr
c WRITE(*,*)'l1rec=',l1rec
C LB
call sifcfg(2, lrecal, nbf, ibvtyp, ifmt, l2rec, n2max, ierr)
C LB
c write(*,*)'LB: after sifcfg(2), ierr=',ierr
c WRITE(*,*)'l2rec=',l2rec
C LB
info(1) = 1 ! =FSPLIT; do not split aoints file
info(2) = l1rec
info(3) = n1max
info(4) = l2rec ! supposed be something else but not ready yet
info(5) = n2max ! supposed be something else but not ready yet
nmax=max(n1max,n2max)
cgk debug
c write(*,*)'gk: n1max = ', n1max
c write(*,*)'gk: l1rec = ', l1rec
c write(*,*)'gk: n2max = ', n1max
c write(*,*)'gk: l2rec = ', l2rec
c write(*,*)'gk: ifmt1 = ', ifmt1
c write(*,*)'gk: ifmt = ', ifmt
cgk end
if(ierr.ne.0)
& call errquit('wrt_dft_aoints: error in sifcfg', ierr,
& MA_ERR)
* if(n1max.le.max1e)
* & call errquit('wrt_dft_aoints: sifcfg does not match max1e',
* & max1e, MA_ERR)
c
c Allocate memory for i4 arrays used in sifs routines in wrt_dft_aoints
c
c get memory for integer map(npbf,2)
sizeof_int=MA_sizeof(MT_INT, 1, MT_BYTE)
cgk debug
* write (*,*) 'gk: MA_sizeof(MT_INT, 1, MT_BYTE) = ',sizeof_int
cgk end
if(sizeof_int.eq.4) then
szintarr=(nbf_ao+1)/2
szlabs=(nmax+1)/2
elseif(sizeof_int.eq.8) then
szintarr=nbf_ao
szlabs=nmax
else
call errquit('wrt_dft_aoints: unsupporte integer type size=',
& sizeof_int, MA_ERR)
endif
if (.not. MA_push_get(MT_Int,(szintarr+1)*nmap,
& 'wrt_dft_aoints:map', l_map, k_map))
$ call errquit('wrt_dft_aoints: ma failed for map',
& szintarr*nmap, MA_ERR)
c get memory for labels
if (.not. MA_push_get(MT_Int,4*szlabs,'wrt_dft_aoints:ilab',l_lab,
& k_lab))
$ call errquit('wrt_dft_aoints: ma failed for lab', 2*szlabs,
& MA_ERR)
call schwarz_tidy()
call int_terminate()
* call int_init(rtdb, 1, AO_bas_han)
* call schwarz_init(geom, AO_bas_han)
* call int_terminate()
* call intd_init(rtdb, 1, AO_bas_han)
cgk debug
* write(*,*)'gk: calling int_init'
cgk end
call int_init(rtdb, 1, AO_bas_han)
cgk debug
* write(*,*)'gk: back from int_init'
cgk end
call int_mem_1e(max1e, mem1)
* call int_mem(max1e, max2e, mem1, mem2)
cgk debug
* write(*,*)'gk: before int_1e_sifs, mem1 = ', mem1
* write(*,*)'gk: mem2 = ', mem2
* write(*,*)'gk: max1e = ', max1e
* write(*,*)'gk: max2e = ', max2e
* call int_mem_print()
cgk end
cgk debug
* write(*,*)'gk: calling int_1e_sifs'
* write(*,*)'gk: mxengy=', mxengy
* write(*,*)'gk: info=', info
cgk write
call int_1e_sifs(AO_bas_han, aoints, energy, mxengy,
& nbf_ao, nmap, int_mb(k_map),int_mb(k_map+szintarr*nmap+1),
& ifmt1, ibvtyp, ibitv, l1rec, n1max, int_mb(k_lab),
& ninfo, info)
cgk debug
* write(*,*)'gk: back from int_1e_sifs'
* write(*,*)'gk: nbf_ao=',nbf_ao
cgk write
cgk debug
c write(*,*)'gk: done with one electron integrals'
c write(*,*)'gk: now let us do 2e ints.'
c write(*,*)'LB: aoints=',aoints,'iunit2=',iunit2,'aoint2=',aoint2
c &,'ierr=',ierr
c WRITE(*,*)'LB: info =',info(1),info(2),info(3),info(4)
cgk end
c open 2 electron itegrals file
iunit2=aoints
call sifo2f(aoints,iunit2,'aoints',info,aoint2,ierr)
if(ierr.ne.0)
& call errquit('wrt_dft_aoints: error in sifo2f', ierr,
& MA_ERR)
cgk debug
* write(*,*)'gk: writing 2e ints to unit, aoint2 = ', aoint2
* write(*,*)'gk: iunit2 = ', iunit2
cgk end
call int_2e_sifs(geom, AO_bas_han, tol2e, oskel, aoint2, ninfo,
& info, ifmt, ibvtyp, ibitv, int_mb(k_lab))
c
c done with memory for sifs files
c
if (.not. MA_pop_stack(l_lab)) call errquit
$ ('wrt_dft_aoints: pop failed at l_lab', 0, GA_ERR)
if (.not. MA_pop_stack(l_map)) call errquit
$ ('wrt_dft_aoints: pop failed at l_map', 0, GA_ERR)
if(ga_nodeid().eq.0) then
close(unit=aoints)
endif
wrt_dft_aoints=.true.
! if (me.eq.0) then
! if (CDFIT.and.IOLGC) then
! call dft_invio('CD', idum, 1, 'KILL', lmiss)
! call dft_invio('CDI', idum, 1, 'KILL', lmiss)
! endif
! if (XCFIT.and.IOLGC) then
! call dft_invio('XC', idum, 1, 'KILL', lmiss)
! endif
! endif
! if (CDFIT.and.(.not.IOLGC)) then
! if (.not. ga_destroy(g_cdinv)) call errquit
! & ('wrt_dft_aoints: could not destroy g_cdinv', 0, GA_ERR)
! if (.not. ga_destroy(g_2ceri)) call errquit
! & ('wrt_dft_aoints: could not destroy g_2ceri', 0, GA_ERR)
! endif
! if (XCFIT.and.(.not.IOLGC))then
! if (.not. ga_destroy(g_xcinv)) call errquit
! & ('wrt_dft_aoints: could not destroy g_xcinv', 0, GA_ERR)
! if (.not. ga_destroy(g_oi)) call errquit
! & ('wrt_dft_aoints: could not destroy g_oi', 0, GA_ERR)
! endif
return
end

61
src/NWints/api/intd_1eov.F
View file

@ -1,3 +1,40 @@
C> Old intd_1eov is now intd_1eovQ (below). This intd_1eov is a wrapper
C> for backward compatibility. intd_1eovQ has a new input argument,
C> hlfdrv, required for computing the half-derivative overlap matrix used
C> by COLUMBUS to compute non-adiabatic derivative couplings. Routines
C> which need to calculate the HDOM should call intd_1eovQ directly,
C> specifying hlfdrv=.true.
C>
C> Lachlan Belcher 14 Feb 2023
subroutine intd_1eov(i_basis,ish,j_basis,jsh,lscr,scr,
& lOva,Ova,idatom)
c $Id$
implicit none
integer i_basis !< [Input] ish basis set handle
integer ish !< [Input] "i" contraction index
integer j_basis !< [Input] jsh basis set handle
integer jsh !< [Input] "j" contraction index
integer lscr !< [Input] length of scratch space "scr"
integer lOva !< [Input] number of overlap integral derivatives in shells ish and jsh
double precision scr(lscr) !< [Scratch] Work space
double precision Ova(lOva) !< [Output] The overlap integral derivatives
c ! NOTE: nint*3 integral derivatives returned per unique center
integer idatom(*) !< [Output] array identifying centers for derivatives
c ! e.g., the first nint*3 derivatives go to center idatom(1)
c ! the second nint*3 derivatives go to center idatom(2)
logical hlfdrv ! [input] if true, calculate same atom as well
! as different atom integrals. This is necessary
! when calculating the half-derivative overlap
! matrix for NADCTs
c
hlfdrv=.false.
CALL intd_1eovQ(i_basis,ish,j_basis,jsh,lscr,scr,
& lOva,Ova,idatom,hlfdrv)
end
C--------------------------------------------------------------------------------------------
C> \ingroup nwint
C> @{
C>
@ -14,8 +51,8 @@ C> of integrals, `ncoord` are the Cartesian coordinates x, y and z, and
C> is 2 for the number of expansion centers on which the basis functions are
C> positioned.
C>
subroutine intd_1eov(i_basis,ish,j_basis,jsh,lscr,scr,
& lOva,Ova,idatom)
subroutine intd_1eovQ(i_basis,ish,j_basis,jsh,lscr,scr,
& lOva,Ova,idatom,hlfdrv)
c $Id$
implicit none
#include "stdio.fh"
@ -34,6 +71,10 @@ c ! NOTE: nint*3 integral derivatives returned per unique
integer idatom(*) !< [Output] array identifying centers for derivatives
c ! e.g., the first nint*3 derivatives go to center idatom(1)
c ! the second nint*3 derivatives go to center idatom(2)
logical hlfdrv ! [input] if true, calculate same atom as well
! as different atom integrals. This is necessary
! when calculating the half-derivative overlap
! matrix for NADCTs
c
c Order is... nint*3*2 (3=> xyz, 2=atoms)
c
@ -74,11 +115,11 @@ c d[idatom(2),z]/
& INT_ERR)
endif
call intd_1eovP(i_basis,ish,j_basis,jsh,lscr,scr,nint,
& OVa,idatom)
& OVa,idatom,hlfdrv)
end
subroutine intd_1eovP(i_basis,ish,j_basis,jsh,lscr,scr,
& nint,Ova,idatom)
& nint,Ova,idatom,hlfdrv)
implicit none
c
c returns 1e overlap derivative integrals for shells/contractions ish,jsh
@ -117,6 +158,10 @@ c ! NOTE: nint*3 integral derivatives returned per unique
integer idatom(*) ! [output] array identifying centers for derivatives
c ! e.g., the first nint*3 derivatives go to center idatom(1)
c ! the second nint*3 derivatives go to center idatom(2)
logical hlfdrv ! [input] if true, calculate same atom as well
! as different atom integrals. This is necessary
! when calculating the half-derivative overlap
! matrix for NADCTs
c
c Order is... nint*3*2 (3=> xyz, 2=atoms)
c
@ -180,7 +225,7 @@ c
c
iatom = (sf_ibs_cn2ce(ish,ibas))
jatom = (sf_ibs_cn2ce(jsh,jbas))
if (iatom.eq.jatom) then
if ((iatom.eq.jatom).and.(.not.hlfdrv))then
idatom(1) = 0
idatom(2) = 0
return
@ -287,7 +332,7 @@ c........................ doS doT doV canAB
c........... nonrel dryrun
& ohnd_ok,oNR,.false.,scr,lscr,rel_dbg,rel_typ)
else if (onw_ok.or.osim_ok) then
call hf1d(
call hf1dQ(
& coords(1,iatom,igeom),
& dbl_mb(mb_exndcf(iexp,ibas)),
& dbl_mb(mb_exndcf(icf,ibas)),
@ -303,7 +348,7 @@ c
& Ova,scr,scr,nint,
c...........overlap, k-e, pot-e, canab, dryrun
& .true., .false., .false., .false., .false.,
& scr,lscr)
& scr,lscr,hlfdrv)
elseif (ohnd_ok) then
call hnd_stvintd(
& coords(1,iatom,igeom),
@ -363,7 +408,7 @@ c
call spcart_tran1e(Ova(1,Zyx,Zatom),scr,
& j_nbf_x,i_nbf_x,0,jgen,
& j_nbf_s,i_nbf_s,Itype,igen,
& .false.)
& .true.)
enddo
enddo
else if (bas_spherical(jbas)) then

83
src/NWints/deriv/hf1d.F
View file

@ -1,4 +1,14 @@
Subroutine hf1d(
C> Old hf1d is now hf1dQ (below). This hf1d is a wrapper
C> for backward compatibility. hf1dQ has a new input argument,
C> hlfdrv, required for computing the half-derivative overlap matrix used
C> by COLUMBUS to compute non-adiabatic derivative couplings. Routines
C> which need to calculate the HDOM should call hf1dQ directly,
C> specifying hlfdrv=.true.
C>
C> Lachlan Belcher 14 Feb 2023
SUBROUTINE hf1d(
& Axyz,Aprims,Acoefs,NPA,NCA,La,ictrA,
& Bxyz,Bprims,Bcoefs,NPB,NCB,Lb,ictrB,
& Cxyz,zan,exinv,ncenters,
@ -8,8 +18,68 @@ c $Id$
Implicit real*8 (a-h,o-z)
Implicit integer (i-n)
logical hlfdrv ! [input] if true, calculate same atom as well
! as different atom integrals. This is necessary
! when calculating the half-derivative overlap
! matrix for NADCTs
integer ictrA, ictrB !NB already used implicitly
Logical dryrun
Logical O2I,KEI,NAI,canAB,GenCon
c--> Cartesian Coordinates, Primitives & Contraction Coefficients
Dimension Axyz(3),Aprims(NPA),Acoefs(NPA,NCA)
Dimension Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB)
c--> Nuclear Cartesian Coordinates, Charges & Inverse Exponents
Dimension Cxyz(3,ncenters),zan(ncenters),exinv(ncenters)
c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals
Dimension bO2I(Nint,*),bKEI(Nint,*),bNAI(Nint,*)
hlfdrv=.false.
CALL hf1dQ(
& Axyz,Aprims,Acoefs,NPA,NCA,La,ictrA,
& Bxyz,Bprims,Bcoefs,NPB,NCB,Lb,ictrB,
& Cxyz,zan,exinv,ncenters,
& bO2I,bKEI,bNAI,Nint,
& O2I,KEI,NAI,canAB,dryrun,W0,maxW0,
& hlfdrv)
end
C-------------------------------------------------------------------------
Subroutine hf1dQ(
& Axyz,Aprims,Acoefs,NPA,NCA,La,ictrA,
& Bxyz,Bprims,Bcoefs,NPB,NCB,Lb,ictrB,
& Cxyz,zan,exinv,ncenters,
& bO2I,bKEI,bNAI,Nint,
& O2I,KEI,NAI,canAB,dryrun,W0,maxW0,
& hlfdrv)
c $Id$
Implicit real*8 (a-h,o-z)
Implicit integer (i-n)
#include "errquit.fh"
#include "stdio.fh"
LOGICAL smatm ! if true, both integrals use the same atom
! center. Necessary for knowing when to
! multiply by the 1/2 factor for half-derivative
! integrals
logical hlfdrv ! [input] if true, calculate same atom as well
! as different atom integrals. This is necessary
! when calculating the half-derivative overlap
! matrix for NADCTs
integer ictrA, ictrB !NB already used implicitly
Logical dryrun
Logical O2I,KEI,NAI,canAB,GenCon
@ -66,6 +136,8 @@ c******************************************************************************
call hf_print_set(0)
#endif
MXD = 1
smatm=.FALSE.
IF(ictrA.EQ.ictrB ) smatm = .TRUE.
c Determine whether general or segmented contraction is used.
@ -188,6 +260,13 @@ c Compute the 2-center overlap integrals, <a|S|b>.
& canAB)
do 200 n = 1,Nint
! WRITE(*,*)"LB in hf1d, smatm=",smatm,"hlfdrv=",hlfdrv
! WRITE(*,*)"ictrA=",ictrA,"ictrB=",ictrB
IF ((.NOT. smatm) .AND. (hlfdrv)) THEN
bO2I(n,1)=-5D-01*bO2I(n,1)
bO2I(n,2)=-5D-01*bO2I(n,2)
bO2I(n,3)=-5D-01*bO2I(n,3)
ENDIF
bO2I(n,4) = -bO2I(n,1)
bO2I(n,5) = -bO2I(n,2)
bO2I(n,6) = -bO2I(n,3)
@ -356,7 +435,7 @@ c (1) the derivative wrt to each of the nuclear attraction centers,
c (2) the derivative wrt to the centers of the Gaussian functions.
if( ictrA.eq.ictrB )then
call daxpy(Nint*3,1.D0,W0(i_VA),1,bNAI(1,((ictrA-1)*3+1)),1)
call daxpy(Nint*3,1.D0,W0(i_VA),1,bNAI(1,((ictrA-1)*3+1)),1)
else
call daxpy(Nint*3,1.D0,W0(i_VA),1,bNAI(1,((ictrA-1)*3+1)),1)
call daxpy(Nint*3,1.D0,W0(i_VB),1,bNAI(1,((ictrB-1)*3+1)),1)

45
src/NWints/sifs/GNUmakefile
View file

@ -1,38 +1,47 @@
#
# $Id$
#
USE_I4FLAGS=1
include ../../config/makefile.h
ifdef NO_SIFS
OBJ = sifs_stubs.o
else
OBJ_OPTIMIZE = sifrd1.o sifr2.o sifr1x.o sifr1n.o sifpre.o \
sifo2f.o sifmge.o siffw2.o siffw1.o siffr2.o siffr1.o sifew2.o \
sifew1.o sife2.o sife1.o sifd2.o sifd1.o sifcfg.o sifc2f.o \
sif2w8.o sifzwh.o sifzh2.o sifwh.o sifw2.o sifw1x.o siftyp.o \
siftdy.o sifstv.o sifskh.o sifsk1.o sifsce.o sifrsh.o sifrh2.o \
sifrh1.o sifrd2.o
OBJ_OPTIMIZE = sifrd1.o sifr2.o sifr1x.o sifr1n.o \
sifo2f.o siffw2.o siffw1.o sifew2.o \
sifew1.o sife2.o sife1.o sifd2.o sifd1.o sifcfg.o \
sif2w8.o sifzwh.o sifzh2.o sifwh.o sifw2.o \
sifskh.o sifsk1.o sifrh2.o \
sifrh1.o sifrd2.o aiopen.o aiwait.o airead.o \
aiwrit.o seqwbf.o seqrbf.o plab1.o plab8.o plab16.o plab32.o \
ulab1.o ulab8.o ulab16.o ulab32.o cdcopy.o izero.o
OBJ = bummer.o
endif
LIBRARY = libnwints.a
ifndef NO_SIFS
USES_BLAS = sifd1.F sifd2.F sife1.F sife2.F
endif
LIB_DEFINES =
LIB_INCLUDES = -I../../nwdft/grid -I../../nwdft/include \
-I../../ddscf
LIB_DEFINES =
include ../../config/makefile.h
include ../../config/makelib.h
bummer.F sif2w8.F sifc2f.F sifcfg.F \
bummer.F \
aiopen.F aiwait.F airead.F aiwrit.F \
seqwbf.F seqrbf.F \
plab1.F plab8.F plab16.F plab32.F \
ulab1.F ulab8.F ulab16.F ulab32.F \
sif2w8.F sifcfg.F \
sifd1.F sifd2.F sife1.F sife2.F \
sifew1.F sifew2.F siffr1.F siffr2.F \
siffw1.F siffw2.F sifmge.F sifo2f.F \
sifpre.F sifr1n.F sifr1x.F sifr2.F \
sifrd1.F sifrd2.F sifrh1.F sifrh2.F sifrsh.F \
sifsce.F sifsk1.F sifskh.F sifstv.F \
siftdy.F siftyp.F sifw1x.F sifw2.F \
sifew1.F sifew2.F \
siffw1.F siffw2.F sifo2f.F \
sifr1n.F sifr2.F \
sifrd1.F sifrd2.F sifrh1.F sifrh2.F \
sifsk1.F sifskh.F \
sifw2.F izero.F\
sifwh.F sifzh2.F sifzwh.F:
./build_sifs.sh

24
src/NWints/sifs/build_sifs.sh
View file

@ -1,6 +1,7 @@
#!/usr/bin/env bash
TGZ=columbus-master.tar.gz
[ -f "$TGZ" ] && gunzip -t "$TGZ" > /dev/null && USE_TGZ=1
USE_TGZ=0
if [[ "$USE_TGZ" == 1 ]]; then
echo "using existing" "$TGZ"
else
@ -8,6 +9,22 @@ rm -rf columbu* *.F sifs.patched
tries=0 ; until [ "$tries" -ge 10 ] ; do \
curl > "$TGZ" https://gitlab.com/api/v4/projects/36816383/repository/archive?path=Columbus/source/gcfci/colib/sifs \
&& curl > bummer.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/humanio/bummer.F?inline=false\
&& curl > bummer.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/humanio/bummer.F?inline=false\
&& curl > aiopen.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/basicio/aiopen.F?inline=false\
&& curl > airead.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/basicio/airead.F?inline=false\
&& curl > aiwait.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/basicio/aiwait.F?inline=false\
&& curl > aiwrit.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/basicio/aiwrit.F?inline=false\
&& curl > seqwbf.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/basicio/seqwbf.F?inline=false\
&& curl > seqrbf.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/basicio/seqrbf.F?inline=false\
&& curl > plab1.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/plab1.F?inline=false\
&& curl > plab8.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/plab8.F?inline=false\
&& curl > plab16.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/plab16.F?inline=false\
&& curl > plab32.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/plab32.F?inline=false\
&& curl > ulab1.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/ulab1.F?inline=false\
&& curl > ulab8.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/ulab8.F?inline=false\
&& curl > ulab16.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/ulab16.F?inline=false\
&& curl > ulab32.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/pack/ulab32.F?inline=false\
&& curl > izero.F https://gitlab.com/columbus-program-system/columbus/-/raw/master/Columbus/source/gcfci/colib/linear_algebra/izero.F?inline=false\
&& break ;\
tries=$((tries+1)) ; echo attempt no. $tries ; sleep 2 ; done
fi
@ -60,3 +77,10 @@ patch -p0 -s -N < siftdy.patch
rm -f siftdy.patch
echo yes > sifs.patched
fi
# dcopy is cursed, won't work with COLUMBUS
sed -i -e 's/dcopy/cdcopy/g' ./sif[d,e][1,2].F
# SO effective density matrix has btype 35. Needed until this change is made in COLUMBUS
sed -i -e 's/ibtypmx=34/ibtypmx=40/g' ./sifr1n.F
sed -i -e 's/btypmx=20/btypmx=40/g' ./sifr1x.F

38
src/NWints/sifs/cdcopy.f Normal file
View file

@ -0,0 +1,38 @@
subroutine cdcopy(n,dx,incx,dy,incy)
c
c copies a vector, x, to a vector, y.
c uses rolled loops for increments equal to one.
c written by ron shepard, based on cdcopy written by
c jack dongarra, linpack, 3/11/78.
c need to avoid blas version for integer rep incompatibility
c with columbus
real*8 dx(*),dy(*)
integer i,incx,incy,ix,iy,m,mp1,n
c
if ( n.le.0 ) return
if ( incx.eq.1 .and. incy.eq.1 ) go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if ( incx.lt.0 ) ix = (-n+1)*incx + 1
if ( incy.lt.0 ) iy = (-n+1)*incy + 1
do 10 i = 1, n
dy(iy) = dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
20 continue
do 50 i = 1, n
dy(i) = dx(i)
50 continue
return
end

3
src/task/GNUmakefile
View file

@ -9,7 +9,8 @@
task_save_state.o task.o task_input.o task_shell_input.o task_num_grad.o \
task_vscf.o task_property.o task_dntmc.o task_bsse.o \
task_jefftce.o task_ncc.o task_rfcf.o \
task_gwmol.o
task_gwmol.o \
task_aoints.o task_nadtrc.o
USES_BLAS = task_num_grad.F task_bsse.F

15
src/task/task.F
View file

@ -154,6 +154,15 @@ c
logical task_gwmol
external task_gwmol
logical task_aoints
external task_aoints
logical task_grdtrc
external task_grdtrc
logical task_nadtrc
external task_nadtrc
c logical embed_optimize
c external embed_optimize
c
@ -377,6 +386,12 @@ c
status = task_energy(rtdb)
else if (inp_compare(.false., 'gradient', operation)) then
status = task_gradient(rtdb)
else if (inp_compare(.false., 'aoints', operation)) then
status = task_aoints(rtdb)
else if (inp_compare(.false., 'grdtrc', operation)) then
status = task_nadtrc(rtdb)
else if (inp_compare(.false., 'nadtrc', operation)) then
status = task_nadtrc(rtdb)
else if (inp_compare(.false., 'property', operation)) then
status = task_property(rtdb)
else if (inp_compare(.false., 'hessian', operation)) then

426
src/task/task_aoints.F Normal file
View file

@ -0,0 +1,426 @@
logical function task_aoints(rtdb)
implicit none
#include "errquit.fh"
#include "rtdb.fh"
#include "mafdecls.fh"
#include "inp.fh"
#include "util.fh"
#include "global.fh"
#include "stdio.fh"
integer rtdb
c
logical status, ignore
logical ocons
character*32 theory
double precision cpu, wall
logical task_aoints_doit
external task_aoints_doit
* logical task_qmmm
* logical task_qmmm_energy
* external task_qmmm_energy
* logical cons_active
* external cons_active
c
c logical task_embed
c logical embed_energy
c external embed_energy
cc
* logical task_bsse
* logical bsse_energy
* external bsse_energy
logical oecce
cc
cgk debug
* write(*,*)'gk: entered task_aoints'
cgk end
oecce = .false.
cpu = util_cpusec()
wall = util_wallsec()
c
ignore = rtdb_delete(rtdb, 'task:aoints')
if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, .false.))
$ call errquit('task_aoints: failed to invalidate status',0,
$ RTDB_ERR)
c
if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory))
$ call errquit('task:aoints: theory not specified',0,
$ RTDB_ERR)
if (theory .ne. 'sodft') then
call errquit('task_aoints: theory not supported for aoints'
& ,0,RTDB_ERR)
else
status=task_aoints_doit(rtdb,theory)
end if
c
* if(task_qmmm) then
* status = task_qmmm_energy(rtdb)
* if (.not. rtdb_get(rtdb, 'qmmm:energy', mt_dbl, 1, energy))
* $ call errquit('task_energy: no energy in rtdb?',0,
* & RTDB_ERR)
* else if (task_embed) then
c status = embed_energy(rtdb)
c if (.not. rtdb_get(rtdb, 'embed:energy', mt_dbl, 1, energy))
c $ call errquit('task_energy: no energy in rtdb?',0,
c & RTDB_ERR)
* else if (theory .ne. 'oniom') then
cgk debug
* write(*,*)'gk: calling task_energy_doit'
cgk end
* status = task_energy_doit(rtdb,theory,energy)
cgk debug
* write(*,*)'gk: back from task_energy_doit'
cgk end
cc
* if(task_bsse) then
* status = bsse_energy(rtdb,theory,energy)
* if(.not.status)
* $ call errquit('task_errr: error call bsse energy',911)
* endif
cc
* else
* status = oniom_energy(rtdb)
* if (status) then
* if (.not. rtdb_get(rtdb, 'oniom:energy', mt_dbl, 1, energy))
* $ call errquit('task_energy: no energy in rtdb?',0,
* & RTDB_ERR)
* if (rtdb_get(rtdb, 'oniom:dipole', mt_dbl, 3, dipole)) then
* if (.not. rtdb_put(rtdb,'task:dipole',mt_dbl,3,dipole))
* $ call errquit('task_energy: failed storing dipole',0,
* & RTDB_ERR)
* else
* ignore = rtdb_delete(rtdb, 'task:dipole')
* endif
* endif
* endif
* ocons = cons_active(rtdb)
* if(.not.rtdb_get(rtdb, 'cons:active', mt_log, 1, ocons)) then
c ocons=.false.
c end if
* if(ocons) call cons_add_energy(rtdb,energy)
* if (.not. rtdb_put(rtdb, 'task:energy', mt_dbl, 1, energy))
* $ call errquit('task_energy: failed storing energy',0,
* & RTDB_ERR)
c
cpu = util_cpusec() - cpu
wall = util_wallsec() - wall
if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, status))
$ call errquit('task_energy: failed to store status',0,
& RTDB_ERR)
if (.not. rtdb_put(rtdb, 'task:cputime', mt_dbl, 1, cpu))
$ call errquit('task_energy: failed storing cputime',0,
& RTDB_ERR)
if (.not. rtdb_put(rtdb, 'task:walltime', mt_dbl, 1, wall))
$ call errquit('task_energy: failed storing walltime',0,
& RTDB_ERR)
c
call ecce_print1('cpu time', mt_dbl, cpu, 1)
call ecce_print1('wall time', mt_dbl, wall, 1)
* if (theory .eq. 'md') oecce=.false.
* if (theory .eq. 'smd') oecce=.false.
* if (status) then
* if (oecce) call geom_print_rtdb_ecce(rtdb)
* call ecce_print1('total energy', mt_dbl, energy, 1)
* call ecce_print_module_exit('task energy', 'ok')
* else
* call ecce_print_module_exit('task energy', 'failed')
* endif
c
task_aoints = status
cgk debug
* write(*,*)'gk: leaving task_aoints'
cgk end
c
end
logical function task_aoints_doit(rtdb,theory)
*
* $Id: task_energy.F,v 1.53 2009/02/04 01:17:17 niri Exp $
*
implicit none
#include "errquit.fh"
#include "rtdb.fh"
#include "mafdecls.fh"
#include "inp.fh"
#include "util.fh"
#include "global.fh"
integer rtdb
character*(*) theory
double precision energy
c
c Generic NWChem interface to compute the ao integrals in SIFS format.
c Currently only the sodft theory is supported.
c
c RTDB input parameters
c ---------------------
c task:theory (string) - name of (QM) level of theory to use
c
c RTDB output parameters
c ----------------------
c task:status (logical)- T/F for success/failure
c if (status) then
* . task:energy (real) - total energy
* . task:dipole(real(3)) - total dipole moment if available
c . task:cputime (real) - cpu time to execute the task
c . task:walltime (real) - wall time to execute the task
c
c Also returns status through the function value
c
* logical scf, rimp2, oimp2, direct_mp2, ccsd, mcscf, selci
* external scf, rimp2, oimp2, direct_mp2, ccsd, mcscf, selci
* logical semi_dir_mp2, rimp2g,nwmd,pspw_energy,band_energy
* external semi_dir_mp2, rimp2g,nwmd,pspw_energy,band_energy
* logical dft_energy, sodft_energy, tddft_energy
* external dft_energy, sodft_energy, tddft_energy
* logical tce_energy_fragment
* external tce_energy_fragment
* logical paw_energy,ofpw_energy
* external paw_energy,ofpw_energy
* logical task_smd_energy
* external task_smd_energy
* logical dc_energy
* external dc_energy
* logical leps_energy
* external leps_energy
logical sodft_aoints
external sodft_energy
c
logical status
* double precision dipole(3)
character*80 key, prefix
c
* logical dftmp2
* integer which_mp2
c
cgk debug
* write(*,*)'gk: in task_aoints_doit'
cgk end
call ecce_print_module_entry('task energy')
c
task_aoints_doit = .false.
c
c Set the theory in the database so that all is consistent for
c the lower level modules. Should not really have to do this,
c but this value is used (appropriately and inappropriately) and
c things are not yet disentangled.
c
if (.not. rtdb_cput(rtdb, 'task:theory', 1, theory))
$ call errquit('task_energy: setting theory?',0,
& RTDB_ERR)
c
c Actually do the deed
c
prefix = theory ! Most common scenario
* if (theory .eq. 'scf') then
* status = scf(rtdb)
* else if (theory .eq. 'smd') then
* status = task_smd_energy(rtdb)
* else if (theory .eq. 'dft') then
* if (.not. rtdb_get(rtdb,'dft:dftmp2',mt_log,1,dftmp2))
* 1 dftmp2 = .false.
* if (dftmp2) then
* if (.not.rtdb_get(rtdb,'dft:which_mp2',mt_int,1,which_mp2))
* 1 which_mp2 = 1
* if (which_mp2.eq.1) then ! semidirect MP2
* prefix = 'mp2'
* if(.not. rtdb_cput(rtdb,'sdmp2:task',1,'energy')) then
* call errquit('task_energy: rtdb_put failed',37,
* 1 RTDB_ERR)
* endif
* status = semi_dir_mp2(rtdb)
* elseif (which_mp2.eq.2) then ! direct MP2
* status = direct_mp2(rtdb)
* elseif (which_mp2.eq.3) then ! RI-MP2
* status = rimp2(rtdb)
* else
* write(6,*) 'task_energy: problem with DFT+MP2'
* endif
* else
* status = dft_energy(rtdb)
* endif
* else if (theory .eq. 'sodft')then
*gk debug
* write(*,*)'gk: calling sodft_energy from task_energy_doit'
*gk end
* status = sodft_aoints(rtdb)
*gk debug
* write(*,*)'gk: back from sodft_energy from task_energy_doit'
*gk end
*gk newcode
* else if (theory .eq. 'integrals') then
* status = integrals_energy(rtdb)
cgk end
* else if (theory .eq. 'tddft') then
* status = tddft_energy(rtdb)
* else if (theory .eq. 'tce') then
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'uccd') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'ccd'))
* 1 call errquit('task_energy: failed writing to rtdb',0,
* & RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ulccd') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'lccd'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'uccsd') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'ccsd'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ulccsd') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'lccsd'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'uqcisd') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'qcisd'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'uccsdt') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'ccsdt'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'uccsdtq') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'ccsdtq'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ucisd') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'cisd'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ucisdt') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'cisdt'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ucisdtq') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'cisdtq'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ump2') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'mbpt2'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ump3') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'mbpt3'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'ump4') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'mbpt4'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'umbpt2') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'mbpt2'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'umbpt3') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'mbpt3'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'umbpt4') then
* prefix = 'tce'
* if (.not.rtdb_cput(rtdb,'tce:model',1,'mbpt4'))
* 1 call errquit('task_energy: failed writing to rtdb',0, RTDB_ERR)
* status = tce_energy_fragment(rtdb)
* else if (theory .eq. 'rimp2') then
* status = rimp2(rtdb)
* else if (theory .eq. 'rimp2_grad') then
* status = rimp2g(rtdb)
* else if (theory .eq. 'oimp2') then
* status = oimp2(rtdb)
* else if (theory .eq. 'direct_mp2') then
*edo prefix = 'mp2'
* status = direct_mp2(rtdb)
* else if (theory .eq. 'ccsd') then
* status = ccsd(rtdb)
* else if (theory .eq. 'ccsd(t)') then
* prefix = 'ccsd'
* status = ccsd(rtdb)
* else if (theory .eq. 'ccsd+t(ccsd)') then
* prefix = 'ccsd'
* status = ccsd(rtdb)
* else if (theory .eq. 'mcscf') then
* status = mcscf(rtdb)
* else if (theory .eq. 'selci') then
* status = selci(rtdb)
* else if (theory .eq. 'mp2' .or. theory.eq.'semi_dir_mp2') then
* prefix = 'mp2'
* if(.not. rtdb_cput(rtdb,'sdmp2:task',1,'energy'))
* $ call errquit('task_energy: rtdb_put failed',0, RTDB_ERR)
* status = semi_dir_mp2(rtdb)
* else if (theory .eq. 'pspw') then
* status = pspw_energy(rtdb)
* else if (theory .eq. 'band') then
* status = band_energy(rtdb)
* else if (theory .eq. 'paw') then
* status = paw_energy(rtdb)
* else if (theory .eq. 'ofpw') then
* status = ofpw_energy(rtdb)
* else if (theory .eq. 'md') then
* status = nwmd(rtdb)
* else if (theory .eq. 'qmmm') then
* if(ga_nodeid().eq.0) then
* write(6,*)'task_energy: qmmm should never be the theory '
* write(6,*)'task_energy: you probably reversed the order ',
* & 'of the task input '
* write(6,*)'task [qmmm] <string theory> [<string operatio',
* & 'n = energy>] [numerical || analytic] [ignore]'
* endif
* status = .false.
* else if (theory .eq. 'embed') then
* if(ga_nodeid().eq.0) then
* write(6,*)'task [embed] <string theory> [<string operatio',
* & 'n = energy>] [numerical || analytic] [ignore]'
* endif
* status = .false.
* else if (theory .eq. 'leps') then
* status = leps_energy(rtdb)
* else if (theory .eq. 'dangchang') then
* status = dc_energy(rtdb)
* hack below
if (theory .eq. 'sodft') then
status = sodft_aoints(rtdb)
else
call errquit('task_aoints: unknown theory',0, INPUT_ERR)
endif
c
* if (status) then
*
* ENERGY is returned via the subroutine interface.
* DIPOLE is currently quietly stored.
*
* key = ' '
* write(key,'(a,'':energy'')') prefix(1:inp_strlen(prefix))
* if (.not. rtdb_get(rtdb, key, mt_dbl, 1, energy))
* $ call errquit('task_energy: no energy in rtdb?',0,
* & RTDB_ERR)
*
* key = ' '
* write(key,'(a,'':dipole'')') prefix(1:inp_strlen(prefix))
* if (rtdb_get(rtdb, key, mt_dbl, 3, dipole)) then
* if (.not. rtdb_put(rtdb, 'task:dipole', mt_dbl, 3, dipole))
* $ call errquit('task_energy: failed storing dipole',0,
* & RTDB_ERR)
* endif
* endif
*
task_aoints_doit = status
c
cgk debug
* write(*,*)'gk: leaving task_aoints_doit'
cgk end
end

9
src/task/task_input.F
View file

@ -231,6 +231,15 @@ c oembed = .true.
> 'rdf')) then
operation='rdf'
else if (inp_compare(.false.,test,
> 'aoints')) then
operation='aoints'
else if (inp_compare(.false.,test,
> 'grdtrc')) then
operation='grdtrc'
else if (inp_compare(.false.,test,
> 'nadtrc')) then
operation='nadtrc'
else if (inp_compare(.false.,test,
> 'lcao')) then
operation='lcao'
else if (inp_compare(.false.,test,

63
src/task/task_nadtrc.F Normal file
View file

@ -0,0 +1,63 @@
logical function task_nadtrc(rtdb)
implicit none
#include "errquit.fh"
#include "rtdb.fh"
#include "mafdecls.fh"
#include "inp.fh"
#include "util.fh"
#include "global.fh"
#include "stdio.fh"
integer rtdb
c
logical status, ignore
logical ocons
character*32 theory
double precision cpu, wall
logical task_aoints_doit
external task_aoints_doit
logical oecce
cc
cgk debug
* write(*,*)'gk: entered task_nadtrc'
cgk end
oecce = .false.
cpu = util_cpusec()
wall = util_wallsec()
c
ignore = rtdb_delete(rtdb, 'task:nadtrc')
if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, .false.))
$ call errquit('task_nadtrc: failed to invalidate status',0,
$ RTDB_ERR)
c
if (.not. rtdb_cget(rtdb, 'task:theory', 1, theory))
$ call errquit('task:nadtrc: theory not specified',0,
$ RTDB_ERR)
if (theory .ne. 'sodft') then
call errquit('task_nadtrc: theory not supported for nadtrc'
& ,0,RTDB_ERR)
else
status=task_aoints_doit(rtdb,theory)
end if
c
cpu = util_cpusec() - cpu
wall = util_wallsec() - wall
if (.not. rtdb_put(rtdb, 'task:status', mt_log, 1, status))
$ call errquit('task_energy: failed to store status',0,
& RTDB_ERR)
if (.not. rtdb_put(rtdb, 'task:cputime', mt_dbl, 1, cpu))
$ call errquit('task_energy: failed storing cputime',0,
& RTDB_ERR)
if (.not. rtdb_put(rtdb, 'task:walltime', mt_dbl, 1, wall))
$ call errquit('task_energy: failed storing walltime',0,
& RTDB_ERR)
c
call ecce_print1('cpu time', mt_dbl, cpu, 1)
call ecce_print1('wall time', mt_dbl, wall, 1)
c
task_nadtrc = status
cgk debug
* write(*,*)'gk: leaving task_nadtrc'
cgk end
c
end