comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 14:15:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Merge pull request #1041 from lachlanbelcher/nwc_sym

Browse source 
Nwc sym
This commit is contained in:
Edoardo Aprà 2024-11-12 10:58:26 -08:00 • committed by GitHub
commit 0779ca2e41
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
10 changed files with 921 additions and 128 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/nwc_columbus/aoints/GNUmakefile
View file

@ -3,7 +3,8 @@
include ../../config/makefile.h
OBJ_OPTIMIZE = dft_aoints.o wrt_dft_aoints.o int_1e_sifs.o \
OBJ_OPTIMIZE = nwc_sym_mod.o \
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 \
@ -11,7 +12,7 @@ include ../../config/makefile.h
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
int_mom_sifs.o hdoverlap.o sym_adapt.o
LIBRARY = libnwc_columbus.a

11
src/nwc_columbus/aoints/dft_aoints.F
View file

@ -18,7 +18,7 @@
* call grid_cleanup(.true.)
*c
* end
logical function sodft_aoints(rtdb)
logical function sodft_aoints(rtdb,theory)
c
c>>> driver
c
@ -50,6 +50,8 @@ c
#include "util.fh"
#include "cgridfile.fh"
#include "cosmo.fh"
#include "geomP.fh"
c
c local declarations
c
@ -67,7 +69,7 @@ c
external dft_main0d, movecs_converged,grid_reopen,xc_gotxc
logical grid_ok,l1ecache
integer igok
character*80 theory
character*32 theory
character*30 operation
c
cgk debug
@ -87,7 +89,7 @@ c
if (.not. rtdb_cget(rtdb, 'task:operation', 1, operation))
$ operation = ' '
if(.not.rtdb_cput(rtdb,'dft:theory', 1, 'sodft'))
if(.not.rtdb_cput(rtdb,'dft:theory', 1, theory))
& call errquit('dft_aoints: can not put theory on rtdb',
& 0,RTDB_ERR)
c
@ -110,7 +112,7 @@ cgk end
cgk debug
!write(6,*)'gk: back from dft_fdist_init from sodft_aoints'
cgk end
c
c
c Check for aoints...probably not neccessary now.
c
* if (.not. rtdb_cget(rtdb, 'task:operation',1,operation))
@ -140,6 +142,7 @@ c
cgk debug
!write(6,*)'gk: calling dft_rdinput from sodft_aoints'
cgk end
call dft_rdinput(rtdb)
cgk debug
!write(6,*)'gk: back from dft_rdinput from sodft_aoints'

3
src/nwc_columbus/aoints/egrad_trace.F
View file

@ -415,7 +415,8 @@ C LB calculate half-derivative overlap matrices
!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)
!nadxfl(i+1,j+1)=-2*ga_ddot(g_hd,g_ad1)
nadxfl(i+1,j+1)=2*ga_ddot(g_hd,g_ad1)
ENDDO
ENDDO
!CALL ga_print(hdol)

175
src/nwc_columbus/aoints/int_1e_sifs.F
View file

@ -1,6 +1,7 @@
subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft,
& nmap,map,imtype,ibvtyp,ibitv,l1rec,n1max,
& clab, ninfo, info)
& clab, ninfo, info, csolab)
USE nwc_sym, ONLY: mxsym,mxbf,zir,nir,nbpsy
implicit none
#include "errquit.fh"
#include "cint1cache.fh"
@ -32,14 +33,14 @@ c
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
& l_sifval,l_SOval,l_SOlab
integer k_buf, k_scr, k_ilab, k_jlab, k_info, k_bfnlab, k_sifbuf,
& k_sifval
& k_sifval,k_SOval,k_SOlab
integer type
logical odoit
double precision q2
external block_int1e ! For T3D
integer i, noffset,g_loc, j, ijmap
integer i, noffset,g_loc, j, ijmap, k
c
logical ocache_save
c
@ -54,27 +55,28 @@ c
integer n1max
INTEGER l1rec, ntitle, ibuf, nsym, ibvtyp, ierr
integer ibitv
integer mxbf
parameter (mxbf=1000)
! integer mxbf
! parameter (mxbf=1000)
c header 1
integer nbas, mxenrgy, nenrgy, nmap
INTEGER ninfo
parameter(ntitle=1,nsym=1,mxenrgy=1)
!parameter(ntitle=1,nsym=1,mxenrgy=1)
parameter(ntitle=1,mxenrgy=1)
c header 2
integer info(ninfo) ! not to be confused with NWChem info
integer ietype(mxenrgy)
integer imtype(nmap) ! map will be taken from dynamic memory
integer map(nbft,nmap)
integer clab(2,*)
integer cSOlab(2,*)
real*8 energy(nenrgy)
character*80 stitle(ntitle)
character*4 slabel(nsym)
integer nbpsy(nsym)
character*4 slabel(mxsym)
character*1 shtypes(-1:7)
data shtypes/'l','s', 'p', 'd', 'f', 'g', 'h', 'i', 'k'/
integer shmap(-1:7)
data shmap / 0, 1, 2, 4, 6, 9, 12, 16, 0 /
integer otype, nprim, ngen, sphcart, iat, shdim, igen, ibf
integer otype, nprim, ngeno, sphcart, iat, shdim, igen, ibf
integer numtot
integer msame,nmsame,nomore
parameter(msame=0, nmsame=1, nomore= 2)
@ -83,6 +85,9 @@ c header 2
real*8 fcore
logical basok
C LB
INTEGER :: nsoints ! number of symmetry orbital integrals
C LB
cgk provisional
* change thresh to the appropriate user supplied zero tolerance
@ -94,7 +99,6 @@ cgk end
c
cgk debug
write(*,*)'gk: **** entered int_1e_sifs'
! WRITE(*,*)'LB: info=',info(1:5)
cgk end
basok=.false.
odbug=.true.
@ -125,6 +129,7 @@ c
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
@ -171,7 +176,11 @@ c get memory for sifs buffer space
& l_sifval, k_sifval))
$ call errquit('int_1e_sifs: ma failed getting sifval', n1mx,
& MA_ERR)
c get memory for SO buffer space
if (.not. MA_push_get(MT_DBL, n1mx,'int_1e_sifs:SOval',
& l_SOval, k_SOval))
$ call errquit('int_1e_sifs: ma failed getting SOval', n1mx,
& MA_ERR)
c
c Set up SIFS header iformation
@ -209,25 +218,35 @@ c put effective nuclear repulsion energy in energy(1)
cgk debug
! write(*,*)'gk: i, iat, shmap(otype)'
! write(*,*)'gk: nmap=',nmap
! WRITE(*,*)"C LB, nshell=",nshell
cgk end
i=0
do ishell=1,nshell
basok=bas_continfo(ibas, ishell, otype, nprim, ngen, sphcart)
basok=bas_continfo(ibas, ishell, otype, nprim, ngeno, sphcart)
! WRITE(*,*)"C LB, ibas=",ibas
! WRITE(*,*)"C LB, ishell=",ishell
! WRITE(*,*)"C LB, otype=",otype
! WRITE(*,*)"C LB, nprim=",nprim
! WRITE(*,*)"C LB, ngeno=",ngeno
! WRITE(*,*)"C LB, sphcart=",sphcart
basok=bas_cn2bfr(ibas, ishell, ilo, ihi)
! WRITE(*,*)"C LB, ilo=",ilo
! WRITE(*,*)"C LB, ihi=",ihi
basok=bas_cn2ce(ibas, ishell, iat)
! WRITE(*,*)"C LB, iat=",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
shdim = shdim / ngeno
else ! Fudge for SP shells
ngen = 1
ngeno = 1
endif
if (otype .le. 1) sphcart = 0
do igen = 1, ngen
do igen = 1, ngeno
do ibf = 1, shdim
i = i + 1
@ -249,9 +268,28 @@ cgk end
enddo
enddo
slabel(1)='a1'
c slabel(1)='a1'
C LB: Add symmetry labels to aoints files
! status = sym_char_table(grp_name,nop,nir,iclass,zir,zclass,chars)
! WRITE(*,*)"C LB, nir=",nir
! WRITE(*,*)"C LB, ibas=",ibas
! WRITE(*,*)"C LB, zir=",zir(1:nir)
nsym=nir
! WRITE(*,*)"C LB, nsym=",nsym
slabel(1:nsym)=zir(1:nsym)
! CALL sym_bas_irreps(ibas,.true.,nbpsy(1:nsym))
! WRITE(*,*)"----Geometry and Symmetry Information----"
! write(*,'(a,i3)')"geom= ",geom
! write(*,'(a,a)')"sym group= ",grp_name
! write(*,*)"slabel= ",slabel(1:nsym)
! WRITE(*,'(a,i3)')"ibas= ",ibas
! WRITE(*,'(a,i3)')"nbft= ",nbft
WRITE(*,'(a,8i3)')"nbpsy= ",nbpsy(1:nsym)
! WRITE(*,'(a,i3)')"sym_num_ops= ",sym_num_ops(geom)
! WRITE(*,'(a,i3)')"nop= ",nop
C LB
stitle(1)='AO integrals from NWChem'
nbpsy(1)=nbft
!nbpsy(1)=nbft
* call sifwh( aoints, ntitle, nsym, nbft,
* & ninfo, nenrgy, nmap, stitle, nbpsy, slabel, info,
@ -272,9 +310,11 @@ cgk debug
cgk end
call sifwh( aoints, ntitle, nsym, nbft,
& ninfo, nenrgy, nmap, stitle, nbpsy, slabel, info,
& ninfo, nenrgy, nmap, stitle,
& nbpsy, slabel(1:nsym), info,
& byte_mb(k_bfnlab), ietype, energy, imtype, map,
& ierr )
WRITE(*,*)"C LB, O"
cgk debug
! write(*,*)'gk: back from sifwh with ierr=',ierr
! WRITE(*,*)"LB after sifwh, info=",info(1:5)
@ -302,18 +342,20 @@ c start loop over integral types here.
end if
cgk debug
! write(*,*)'gk: processing integral type: ', integ_type
write(*,*)'gk: processing integral type: ', integ_type
cgk end
do jshell = 1, nshell
do ishell = 1, nshell
do ishell =1, nshell
cgk debug
! write(*,*)'gk: ishell=',ishell,' jshell=',jshell
write(*,*)'gk: ishell=',ishell,' jshell=',jshell
cgk end
odoit = .true.
if (oskel)
$ odoit = sym_shell_pair(ibas, ishell, jshell, q2)
! if (oskel)
! & odoit = sym_shell_pair(ibas, ishell, jshell, q2)
! WRITE(*,*)"C LB, odoit= ", odoit
! WRITE(*,*)"C LB, q2= ", q2
if (odoit) then
if (.not. bas_cn2bfr(ibas, ishell, ilo, ihi))
@ -324,56 +366,40 @@ cgk end
call ifill(max1e,0,int_mb(k_ilab),1)
call ifill(max1e,0,int_mb(k_jlab),1)
if (type .eq. 1) then
if (type .eq. 1) then !Kinetic Energy
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
else if (type .eq. 2) then !Potential Energy
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
else if (type .eq. 3) then !Overlap
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,
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))
WRITE(*,*)"C LB, calling sym_1int 1"
CALL sym_1int(ibuf,nsoints,
& dbl_mb(k_sifval), clab,
& dbl_mb(k_SOval), cSOlab)
WRITE(*,'(a,1i4)')"C LB, nsoints =",nsoints
WRITE(*,*)"C LB, cSOlab ="
WRITE(*,'(2i4)') cSOlab(1:2,1:nsoints)
C LB
call sifew1(aoints, info, 2, ibuf, last,
WRITE(*,*)"C LB, calling sifew1"
call sifew1(aoints, info, 2, nsoints, last,
& itypea, itypeb, ibvtyp,
& dbl_mb(k_sifval), clab, fcore, ibitv,
& dbl_mb(k_SOval), cSOlab, 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))
@ -385,32 +411,23 @@ c contains the unwritten values, ditto for labels
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
enddo ! i
end if
ijshell = ijshell + 1
end do
end do
end do !i_shell
end do !j_shell
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_1e_sifs calling sifew1 2"
cgk end
call sifew1(aoints, info, 2, ibuf, last,
WRITE(*,*)"C LB, calling sym_1int 2"
CALL sym_1int(ibuf,nsoints,
& dbl_mb(k_sifval), clab,
& dbl_mb(k_SOval), cSOlab)
WRITE(*,*)"C LB, calling sifew1"
call sifew1(aoints, info, 2, nsoints, last,
& itypea, itypeb, ibvtyp,
& dbl_mb(k_sifval), clab, fcore, ibitv,
& dbl_mb(k_SOval), cSOlab, fcore, ibitv,
& dbl_mb(k_sifbuf), nrec, ierr)
write(6,'(a,i10,1x,a,a)') 'Wrote ',numtot, integ_type,
@ -427,7 +444,7 @@ cgk debug
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)
& fcore, ibvtyp, ibitv,dbl_mb(k_SOval),cSOlab)
cgk debug
* write(*,*)'gk: calling int_mom_sifs from int_1e_sifs'
@ -435,13 +452,15 @@ cgk debug
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)
& fcore, ibvtyp, ibitv,dbl_mb(k_SOval),cSOlab)
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_SOval)) call errquit
$ ('int_1e_sifs: pop failed at l_SOval', 0, GA_ERR)
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

100
src/nwc_columbus/aoints/int_mom_sifs.F
View file

@ -1,6 +1,8 @@
subroutine int_mom_sifs(ibas, oskel, aoints, nbft,
& l1rec, n1max,
& sifbuf, sifval, ninfo, info, fcore, ibvtyp, ibitv)
& sifbuf, sifval, ninfo, info, fcore, ibvtyp, ibitv,
& SOval, cSOlab)
USE nwc_sym, ONLY: mxsym,mxbf
implicit none
#include "errquit.fh"
#include "cint1cache.fh"
@ -12,6 +14,7 @@
#include "bas.fh"
#include "cscfps.fh"
#include "sym.fh"
#include "geom.fh"
c
c Oskel indicates that the skeleton (petite-list symmetry) matrix should be
c built ...
@ -25,6 +28,7 @@ c
integer ninfo
integer info(ninfo) ! not to be confused with NWChem info
integer clab(2,n1max)
integer cSOlab(2,n1max)
real*8 fcore
c
@ -55,6 +59,7 @@ C LB
integer lval !Order of dipole integrals
double precision centerl(3) !Center of dipole integrals
integer num_MPint !Number of dipole integrals
integer geom
C LB
c
c sifs parameters
@ -63,12 +68,13 @@ c
integer l1rec, n1max
integer ibvtyp
integer ibitv
integer mxbf
! integer mxbf
c make sure this is compatible with columbus
parameter (mxbf=1000)
! parameter (mxbf=1000)
c header 1
integer ntitle, nsym, nbas, mxenrgy
parameter(ntitle=1,nsym=1,mxenrgy=1)
! parameter(ntitle=1,nsym=1,mxenrgy=1)
parameter(ntitle=1,mxenrgy=1)
c header 2
integer otype, nprim, ngen, sphcart, iat, shdim, igen, ibf
integer numtot
@ -88,10 +94,14 @@ cgk provisional
parameter (thresh=1d-12)
cgk end
integer numints, icart, kcart, int_type
C LB
INTEGER :: nsoints ! number of symmetry orbital integrals
DOUBLE PRECISION :: SOval(n1max)
C LB
c
cgk debug
!write(*,*)'gk: **** entered int_mom_sifs'
write(*,*)'gk: **** entered int_mom_sifs'
cgk end
basok=.false.
odbug=.true.
@ -194,12 +204,12 @@ c int_type = 1, linear momentum
c int_type = 2, angular momentum
c int_type = 3, dipole moments
do int_type=1,3
!WRITE(*,*)"LB, starting int_type loop",int_type
WRITE(*,*)"LB, starting int_type loop",int_type
if(int_type.eq.1.or.int_type.eq.2) itypea=2 ! gonna do dipole
if(int_type.eq.3) itypea=1
!WRITE(*,*)"LB, starting int_type loop, itypea=",itypea
WRITE(*,*)" LB, starting int_type loop, itypea=",itypea
if (0 .eq. ga_nodeid()) then
cgk debug
@ -211,7 +221,7 @@ c start loop over integrals here.
do jshell = 1, nshell
do ishell = 1, nshell
cgk debug
* write(*,*)'gk: ishell=',ishell,' jshell=',jshell
write(*,'(a,i4,a,i4)')'gk: ishell=',ishell,' jshell=',jshell
cgk end
odoit = .true.
if (oskel)
@ -248,29 +258,53 @@ C LB
C LB Dipoles
elseif(int_type.eq.3) then
!write(*,*)'LB: calling int_mpolel'
lval=1;
centerl=0.0;!Should this be coordinate origin or
!atom center?
lval=1;! 1 for dipoles
IF (.NOT. bas_geom(ibas,geom)) CALL errquit
$ ('int_mom_sifs: bas_geom failed for ibas', ibas,
& BASIS_ERR)
IF (.NOT. geom_center_of_mass(geom, centerl)) CALL
& errquit ('int_mom_sifs: geom_center_of_mass failed'
& ,geom,GEOM_ERR)
WRITE(*,*)"LB, calling int_mpolel"
CALL int_mpolel (ibas, ishell, ibas, jshell,
$ lval, centerl,
$ mem1mom, dbl_mb(k_scr), max1mom, dbl_mb(k_buf),
$ Num_MPint)
WRITE(*,'(a,2i5)')"LB, ishell, jshell =", ishell, jshell
WRITE(*,'(a,1i5)')"LB, Num_MPint =", Num_MPint
WRITE(*,*)"poles="
WRITE(*,*)dbl_mb(k_buf:k_buf+Num_MPint-1)
C LB
endif
noffsetz = 0
noffsety = (ihi-ilo+1)*(jhi-jlo+1)
noffsetx = (ihi-ilo+1)*(jhi-jlo+1)*2
IF (int_type .lt. 3) THEN
noffsetz = 0
noffsety = (ihi-ilo+1)*(jhi-jlo+1)
noffsetx = (ihi-ilo+1)*(jhi-jlo+1)*2
ELSE ! dipoles
noffsetz = 0
noffsety = 1
noffsetx = 2
ENDIF
WRITE(*,'(a,3i4)')"LB, noffsetz,noffsety,noffsetx=",noffsetz,
& noffsety,noffsetx
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
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 (int_type .eq. 3) THEN! dipoles
!ijmap=ijmap*(ihi-ilo+1)
ijmap=(i-ilo)*(jhi-jlo+1)+(j-jlo)
ijmap=ijmap*3
ENDIF
if(i.ge.j) then
symmap=nbft*(j-1)-((j-1)*j)/2+i
WRITE(*,'(a,2i4)')"LB, ijmap, symmap=",ijmap,symmap
* write(*,*)'gk: symmap=',symmap
dbl_mb(k_momx+symmap-1)=
& dbl_mb(k_buf+noffsetx+ijmap)
@ -281,15 +315,15 @@ cgk end
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)
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 ! ilo
@ -360,7 +394,8 @@ cgk debug
numtot = numtot + ibuf
C LB
!WRITE(*,*)"LB: in int_mom_sifs"
!WRITE(*,*)"LB: itypea,itypeb=",itypea,itypeb
WRITE(*,*)"LB: in int_mom_sifs, ibuf=",ibuf
WRITE(*,*)"LB: itypea,itypeb=",itypea,itypeb
!WRITE(*,*)"LB: calling sifew1 1 last=",last,"ibvtyp=",ibvtyp
C LB
call sifew1(aoints, info, 2, ibuf, msame,
@ -395,12 +430,17 @@ cgk debug
cgk end
C LB
!WRITE(*,*)"LB: in int_mom_sifs"
!WRITE(*,*)"LB: itypea,itypeb=",itypea,itypeb
WRITE(*,*)"LB: in int_mom_sifs, ibuf=",ibuf
WRITE(*,*)"LB: itypea,itypeb=",itypea,itypeb
WRITE(*,*)"C LB, calling sym_1int 2"
CALL sym_1int(ibuf,nsoints,
& sifval, clab,
& SOval, cSOlab)
!WRITE(*,*)"LB: calling sifew1 2 last=",last,"ibvtyp=",ibvtyp
C LB
call sifew1(aoints, info, 2, ibuf, last,
call sifew1(aoints, info, 2, nsoints, last,
& itypea, itypeb, ibvtyp,
& sifval, clab, fcore, ibitv,
& SOval, cSOlab, fcore, ibitv,
& sifbuf, nrec, ierr)
write(6,'(a,i10,1x,a,a)') 'Wrote ',numtot, integ_type,
& ' integrals to aoints'

21
src/nwc_columbus/aoints/int_so_sifs.F
View file

@ -1,6 +1,8 @@
subroutine int_so_sifs(ibas, oskel, aoints, nbft, max1e, mem1,
& l1rec, n1max,
& sifbuf, sifval, ninfo, info, clab, fcore, ibvtyp, ibitv)
& sifbuf, sifval, ninfo, info, clab, fcore, ibvtyp, ibitv,
& SymOrbval, cSOlab)
USE nwc_sym, ONLY: mxsym,mxbf
implicit none
#include "errquit.fh"
#include "cint1cache.fh"
@ -25,6 +27,7 @@ c
integer ninfo
integer info(ninfo) ! not to be confused with NWChem info
integer clab(2,max1e)
integer cSOlab(2,max1e)
real*8 fcore
c
@ -56,9 +59,9 @@ c
integer l1rec, n1max
integer ibvtyp
integer ibitv
integer mxbf
! integer mxbf
c make sure this is compatible with columbus
parameter (mxbf=1000)
! parameter (mxbf=1000)
c header 1
integer ntitle, nsym, nbas, mxenrgy
parameter(ntitle=1,nsym=1,mxenrgy=1)
@ -83,6 +86,10 @@ cgk provisional
parameter (thresh=1d-12)
cgk end
integer numints, icart, kcart
C LB
INTEGER :: nsoints ! number of symmetry orbital integrals
DOUBLE PRECISION :: SymOrbval(n1max)
C LB
c
cgk debug
@ -310,15 +317,19 @@ cgk end
* endif
numtot=numtot+ibuf
cgk debug
WRITE(*,*)"C LB, calling sym_1int 2"
CALL sym_1int(ibuf,nsoints,
& sifval, clab,
& SymOrbval, cSOlab)
* 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,
call sifew1(aoints, info, 2, nsoints, last,
& itypea, itypeb, ibvtyp,
& sifval, clab, fcore, ibitv,
& SymOrbval, cSOlab, fcore, ibitv,
& sifbuf, nrec, ierr)
write(6,'(a,i10,1x,a,a)') 'Wrote ',numtot, integ_type,
& ' integrals to aoints'

52
src/nwc_columbus/aoints/nwc_sym_mod.F Normal file
View file

@ -0,0 +1,52 @@
MODULE nwc_sym
IMPLICIT NONE
SAVE
#include "nwc_const.fh" !nw_max_shells, nw_max_nbf
#include "geomP.fh" !max_cent, coords
INTEGER, PARAMETER :: mxsym = 8 !max number of symmetry groups
INTEGER, PARAMETER :: mxshel=nw_max_shells!max number of shells
INTEGER, PARAMETER :: mxbf=nw_max_nbf !max number of basis functions
INTEGER, PARAMETER :: mxqn=7
INTEGER, PARAMETER :: mxaqn=mxqn*(mxqn+1)/2
INTEGER, PARAMETER :: mxqnm=21
INTEGER, PARAMETER :: mxaqnm=mxqnm*(mxqnm+1)/2
INTEGER :: isymax(3,2) ! irreps of (c1) x,y,z, (c2) Rx,Ry,Rz
INTEGER :: sh_op_map(8,mxshel), bf_op_map(8,mxbf)
INTEGER :: sh_n_uq_op(mxshel), bf_n_uq_op(mxbf)
INTEGER :: sh_uq_op(8,mxshel), bf_uq_op(8,mxbf)
INTEGER :: sh_n_uq, bf_n_uq
INTEGER :: sh_uq(mxshel), bf_uq(mxbf)
INTEGER :: sh_nat(2,mxshel), bf_nat(2,mxbf), sh_uq_bf(4,mxshel)
INTEGER :: bf_per_ir(0:7)
INTEGER :: bf_per_ir_cum(0:7)
INTEGER :: bf_so_ir(0:7,mxbf)
INTEGER :: so_uq_to_lab(2,0:7,mxbf)
INTEGER :: so_lab1(3,mxbf)
INTEGER :: so_lab2(mxbf)
INTEGER :: isymao(mxqn,mxaqn)
INTEGER :: stabilizer(max_cent) !stabilizer for each center
INTEGER :: ngen ! number of generating operators
INTEGER :: parbit(0:7) = (/1,-1,-1,1,-1,1,1,-1/) !parity of 3-bit integer 000-111
INTEGER :: itran(mxbf,mxbf)
INTEGER :: irrepmap(0:mxsym-1)=(/0,1,2,3,4,5,6,7/) !map from Dalton to NWChem irrep order
INTEGER :: opmap(0:mxsym-1)=(/0,1,2,3,4,5,6,7/) !map from Dalton to NWChem operator order
INTEGER :: wt(mxbf,2) !How many AOs make up each SO
INTEGER :: nbpsy(1:mxsym) ! number of basis functions per symmetry
INTEGER :: nop !number of symmetry operators
INTEGER :: nir !number of irreps
LOGICAL :: oprint
DOUBLE PRECISION :: bf_phase(8,mxbf)
DOUBLE PRECISION :: char_tab(1:8,0:7)
CHARACTER*8 :: grp_name ! name of symmetry group
CHARACTER*8 :: zir(mxsym) ! irrep names
END MODULE nwc_sym

651
src/nwc_columbus/aoints/sym_adapt.F Normal file
View file

@ -0,0 +1,651 @@
SUBROUTINE sym_adapt(ibas)
! Routine to adapt the AOs to symmetrized orbitals
USE nwc_sym
IMPLICIT NONE
#include "errquit.fh"
!#include "nwc_const.fh"
!#include "sym_adapt.fh"
#include "sym.fh"
#include "geom.fh"
#include "bas.fh"
#include "basP.fh"
! IN
INTEGER, INTENT(IN) :: ibas ! basis handle
! LOCAL
INTEGER :: geom ! geometry handle
INTEGER :: i,j,k,ii,jj,kk,mm ! counters
INTEGER :: l(mxaqn),m(mxaqn),n(mxaqn)
INTEGER :: Lvar, Mvar, Nvar
INTEGER :: ang_max !maximum L for basis ibas
INTEGER :: iclass(mxsym),otype
INTEGER :: naorb, naorbd, nsorb, ivarb, norbpsh, ncent
INTEGER :: jkb, cent
INTEGER :: ireplm ! integer function
INTEGER :: nprim, ngeno, spcrt ! dummies for bas_continfo
INTEGER :: mult(0:7)
INTEGER :: lorb, lex_least
DOUBLE PRECISION :: chars(mxsym*mxsym)
CHARACTER*8 :: zclass(mxsym)
LOGICAL :: status
! Get geometry (geom)
IF (.NOT. bas_geom(ibas,geom)) CALL errquit
$ ('sym_adapt: bas_geom failed for ibas', ibas,
& BASIS_ERR)
! Get number of centers (ncent)
IF (.NOT. geom_ncent(geom,ncent)) CALL errquit
$ ('sym_adapt: geom_ncent failed for geom', geom,
& GEOM_ERR)
WRITE(*,*)"C LB, ncent=",ncent
! Get highest angular momentum for basis (ang_max)
status = bas_high_angular(ibas,ang_max)
! Get group name (grp_name)
CALL sym_group_name(geom,grp_name)
IF (.NOT. sym_abelian_group(geom)) CALL errquit
$ ('sym_adapt: Columbus can only process abelian point
& groups',geom, GEOM_ERR)
! Get group info A
status = sym_char_table(grp_name,nop,nir,iclass,zir,zclass,chars)
WRITE(*,*)"C LB, nop =",nop
! Get group info B
CALL sym_bas_irreps(ibas,.true.,nbpsy(1:nir))
! Get group info C
oprint=.TRUE.
CALL sym_abelian_bas_info(ibas,
$ sh_op_map, bf_op_map, bf_phase,
$ sh_n_uq_op, bf_n_uq_op,
$ sh_uq_op, bf_uq_op,
$ sh_n_uq, bf_n_uq,
$ sh_uq, bf_uq,
$ sh_nat, bf_nat, sh_uq_bf,
$ char_tab,
$ bf_per_ir,
$ bf_per_ir_cum,
$ bf_so_ir,
$ zir,
$ so_uq_to_lab, so_lab1, so_lab2, oprint)
! WRITE(*,'(a)')"C LB apply operator to shell"
! WRITE(*,'(a)')"C LB sh_op_map ="
! WRITE(*,'(8i4)')sh_op_map(:,1:6)
! WRITE(*,'(a)')"C LB apply operator to basis function"
! WRITE(*,'(a)')"C LB bf_op_map ="
! WRITE(*,'(8i4)')bf_op_map(:,1:10)
! WRITE(*,'(a)')"C LB apply operator to basis function (phase)"
! WRITE(*,'(a)')"C LB bf_phase ="
! WRITE(*,'(8f6.1)')bf_phase(:,1:10)
! WRITE(*,'(a)')"C LB number of operations that produce unique
! & shells"
! WRITE(*,'(a)')"C LB sh_n_uq_op ="
! WRITE(*,'(6i4)')sh_n_uq_op(1:6)
! WRITE(*,'(a)')"C LB number of operations that produce unique
! & basis functions"
! WRITE(*,'(a)')"C LB bf_n_uq_op ="
! WRITE(*,'(10i4)')bf_n_uq_op(1:10)
! WRITE(*,'(a)')"C LB list of operations that produce unique
! & shells"
! WRITE(*,'(a)')"C LB sh_uq_op ="
! WRITE(*,'(8i8)')sh_uq_op(:,1:6)
! WRITE(*,'(a)')"C LB list of operations that produce unique
! & bfns"
! WRITE(*,'(a)')"C LB bf_uq_op ="
! WRITE(*,'(8i4)')bf_uq_op(:,1:10)
! WRITE(*,'(a)')"C LB number of unique shells"
! WRITE(*,'(a)')"C LB sh_n_uq ="
! WRITE(*,'(i4)')sh_n_uq
! WRITE(*,'(a)')"C LB number of unique bfns"
! WRITE(*,'(a)')"C LB bf_n_uq ="
! WRITE(*,'(i4)')bf_n_uq
! WRITE(*,'(a)')"C LB lexically highest sym-related shell"
! WRITE(*,'(a)')"C LB sh_uq ="
! WRITE(*,'(6i4)')sh_uq(1:6)
! WRITE(*,'(a)')"C LB lexically highest sym-related bfn"
! WRITE(*,'(a)')"C LB bf_uq ="
! WRITE(*,'(10i4)')bf_uq(1:10)
! WRITE(*,'(a)')"C LB map shells from nwchem to natural AO order"
! WRITE(*,'(a)')"C LB sh_nat1 ="
! WRITE(*,'(6i4)')sh_nat(1,1:6)
! WRITE(*,'(a)')"C LB map shells from natural to nwchem AO order"
! WRITE(*,'(a)')"C LB sh_nat2 ="
! WRITE(*,'(6i4)')sh_nat(2,1:6)
! WRITE(*,'(a)')"C LB map bfn from nwchem to natural AO order"
! WRITE(*,'(a)')"C LB bf_nat ="
! WRITE(*,'(10i4)')bf_nat(1,1:10)
! WRITE(*,'(a)')"C LB map bfn from natural to nwchem AO order"
! WRITE(*,'(a)')"C LB bf_nat ="
! WRITE(*,'(10i4)')bf_nat(2,1:10)
! WRITE(*,'(a)')"1st natural bfn from unique shell"
! WRITE(*,'(a)')"C LB sh_uq_bf ="
! WRITE(*,'(6i4)')sh_uq_bf(1,1:6)
! WRITE(*,'(a)')"last natural bfn from unique shell"
! WRITE(*,'(a)')"C LB sh_uq_bf ="
! WRITE(*,'(6i4)')sh_uq_bf(2,1:6)
! WRITE(*,'(a)')"1st unique bfn from unique shell"
! WRITE(*,'(a)')"C LB sh_uq_bf ="
! WRITE(*,'(6i4)')sh_uq_bf(3,1:6)
! WRITE(*,'(a)')"last unique bfn from unique shell"
! WRITE(*,'(a)')"C LB sh_uq_bf ="
! WRITE(*,'(6i4)')sh_uq_bf(4,1:6)
! WRITE(*,'(a)')"C LB bf_per_ir ="
! WRITE(*,'(8i4)')bf_per_ir(:)
! WRITE(*,'(a)')"list of irreps generated from unique bfn"
! WRITE(*,'(a)')"C LB bf_so_ir ="
! WRITE(*,'(8i4)')bf_so_ir(:,1:10)
! WRITE(*,'(a)')"C LB so_uq_to_lab 1 ="
! WRITE(*,'(8i12)')so_uq_to_lab(1,:,1:10)
! WRITE(*,'(a)')"C LB so_uq_to_lab 2 ="
! WRITE(*,'(8i12)')so_uq_to_lab(2,:,1:10)
! WRITE(*,'(a)')"irrep of given SO"
! WRITE(*,'(a)')"C LB so_lab1 ="
! WRITE(*,'(10i4)')so_lab1(1,1:10)
! WRITE(*,'(a)')"map twixt labelling schemes"
! WRITE(*,'(a)')"C LB so_lab1 ="
! WRITE(*,'(10i4)')so_lab1(2,1:10)
! WRITE(*,'(a)')"unique bfn that is generator"
! WRITE(*,'(a)')"C LB so_lab1 ="
! WRITE(*,'(10i4)')so_lab1(3,1:10)
! WRITE(*,'(a)')"map from sym-blocked SO order to natural SO order"
! WRITE(*,'(a)')"C LB so_lab2 ="
! WRITE(*,'(10i4)')so_lab2(1:10)
! WRITE(*,*)"C LB, before isymax, grp_name=",grp_name
! number of generating operators
ngen = INT(LOG(DBLE(nop))/LOG(2.0))
! Initialize isymax, reinitialize irrep and operator maps
! column 1 of isymax is irreps of x, y, z
! column 2 of isymax is irreps of Rx,Ry,Rz
! Note that for the math here, the irreps are in DALTON order!
! This will be fixed with the irrep and operator maps.
!
! Group Irrep order Operator Order
! ----- ------------------------------ ---------------------------
! C1 A E
! Cs A' A'' E Oxy
! Ci Ag Au E i
! C2 A B E C2z
! D2 A B2 B1 B3 E C2y C2x C2z
! C2v A1 B1 B2 A2 E Oyz Oxz C2z
! C2h Ag Bu Au Bg E Oxy C2z i
! D2h Ag B3u B2u B1g B1u B2g B3g Au E Oyz Oxz C2z Oxy C2y C2x i
SELECT CASE (grp_name)
CASE ('C1')
ISYMAX(1,1) = 0
ISYMAX(2,1) = 0
ISYMAX(3,1) = 0
ISYMAX(1,2) = 0
ISYMAX(2,2) = 0
ISYMAX(3,2) = 0
CASE ('Cs')
ISYMAX(1,1) = 0
ISYMAX(2,1) = 0
ISYMAX(3,1) = 1
ISYMAX(1,2) = 1
ISYMAX(2,2) = 1
ISYMAX(3,2) = 0
CASE ('Ci')
ISYMAX(1,1) = 1
ISYMAX(2,1) = 1
ISYMAX(3,1) = 1
ISYMAX(1,2) = 0
ISYMAX(2,2) = 0
ISYMAX(3,2) = 0
CASE ('C2')
ISYMAX(1,1) = 1
ISYMAX(2,1) = 1
ISYMAX(3,1) = 0
ISYMAX(1,2) = 1
ISYMAX(2,2) = 1
ISYMAX(3,2) = 0
CASE ('D2')
ISYMAX(1,1) = 3
ISYMAX(2,1) = 1
ISYMAX(3,1) = 2
ISYMAX(1,2) = 3
ISYMAX(2,2) = 1
ISYMAX(3,2) = 2
irrepmap(1)=2
irrepmap(2)=1
irrepmap(3)=3
CASE ('C2v')
ISYMAX(1,1) = 1
ISYMAX(2,1) = 2
ISYMAX(3,1) = 0
ISYMAX(1,2) = 2
ISYMAX(2,2) = 1
ISYMAX(3,2) = 3
irrepmap(1)=3
irrepmap(2)=1
irrepmap(3)=2
opmap(1)=3
opmap(2)=2
opmap(3)=1
CASE ('C2h')
ISYMAX(1,1) = 1
ISYMAX(2,1) = 1
ISYMAX(3,1) = 2
ISYMAX(1,2) = 3
ISYMAX(2,2) = 3
ISYMAX(3,2) = 0
irrepmap(1)=2
irrepmap(2)=3
irrepmap(3)=1
opmap(1)=1
opmap(2)=3
opmap(3)=2
CASE ('D2h')
ISYMAX(1,1) = 1
ISYMAX(2,1) = 2
ISYMAX(3,1) = 4
ISYMAX(1,2) = 6
ISYMAX(2,2) = 5
ISYMAX(3,2) = 3
irrepmap(1)=7
irrepmap(2)=3
irrepmap(3)=4
irrepmap(4)=5
irrepmap(5)=2
irrepmap(6)=6
irrepmap(7)=1
opmap(1)=3
opmap(2)=5
opmap(3)=6
opmap(4)=7
opmap(5)=4
opmap(6)=2
opmap(7)=1
CASE DEFAULT
CALL errquit ('sym_adapt: symmetry not recognized',
& grp_name, GEOM_ERR)
END SELECT
! WRITE(*,*)"C LB, isymax="
! WRITE(*,'(3i4)')isymax(1:3,1:2)
! Initialize the stabilizer list
CALL stab_cent(geom,ncent)
! ! initialize mult
! DO i = 0,7
! mult(i)=2**MAX(0,ngen-is(i))
! ENDDO
! WRITE(*,*)"C LB, ibas =",ibas
! WRITE(*,*)"C LB, bas_is_spherical =",bas_is_spherical(ibas)
! WRITE(*,*)"C LB, ang_max=",ang_max
!
! Initialize isymao
IF (.NOT. bas_is_spherical(ibas)) THEN
! ----- using cartesian basis -----
DO i = 1, ang_max + 1
CALL lmnval(i,i*(i + 1)/2,l,m,n)
DO j = 1, i*(i + 1)/2
Lvar = MOD(l(j),2)*isymax(1,1)
Mvar = MOD(m(j),2)*isymax(2,1)
Nvar = MOD(n(j),2)*isymax(3,1)
isymao(i,j) = IEOR(Lvar,IEOR(Mvar,Nvar))
ENDDO
ENDDO
ELSE
! ----- using spherical basis -----
DO i = 0, ang_max
ii = i + 1
IF (i .EQ. 0) THEN
isymao(ii,1) = ireplm(0,0)
ELSEIF (i .EQ. 1) THEN
isymao(ii,1) = ireplm(1, 1)
isymao(ii,2) = ireplm(1,-1)
isymao(ii,3) = ireplm(1, 0)
ELSE
DO j = -i, i
jj = j + i + 1
isymao(ii,jj) = ireplm(i,j)
ENDDO
ENDIF
ENDDO
ENDIF
! ----- Calculate transformation matrix from AO to SO -----
nsorb = 0 ! number of symmetry orbitals
! --- Loop over irreps
DO ii = 0, nir-1
i=irrepmap(ii)
! WRITE(*,'(a,1i4)')" C LB, irrep= ",i
naorbd = 0
! --- Loop over unique shells
DO j = 1, sh_n_uq
jj=sh_uq(j)
IF (.NOT. bas_cn2ce(ibas,jj,cent)) CALL errquit
& ("sym_adapt: can't find center for bfn",jj,GEOM_ERR)
status = bas_continfo(ibas,jj,otype,nprim,ngeno,spcrt)
norbpsh = sh_uq_bf(4,j)-sh_uq_bf(3,j)+1
! --- Loop over orbitals
DO k = 1, norbpsh
naorbd = naorbd + 1
! WRITE(*,'(a,2i4)')" C LB, i=",i
! WRITE(*,'(a,2i4)')" C LB, shell,orbital=",jj,k
! WRITE(*,'(a,2i4)')" C LB, otype+1=",otype+1
! WRITE(*,'(a,2i4)')" C LB, isymao=",isymao(otype+1,k)
! naorb = naorb + 1
ivarb = IEOR(i,isymao(otype+1,k))
! WRITE(*,'(a,i4)')"ivarb =",ivarb
IF (IAND(stabilizer(cent),ivarb).EQ.0) THEN
! WRITE(*,*)" C LB, orbital contributes!"
nsorb = nsorb + 1
! WRITE(*,'(a,1i4)')
! &" C LB, this will be symmetrized orbital",nsorb
jkb = 0
! WRITE(*,'(a,1i4)')"C LB, naorbd=",naorbd
! --- Loop over symmetry operations
DO mm = 1, nop
kk = opmap(mm - 1)
! WRITE(*,'(a,2i4)')" C LB, mm=",mm
! WRITE(*,'(a,2i4)')"C LB, mm,bf_uq_op=",
! $ mm,bf_uq_op(mm,bf_uq(naorbd))
IF(bf_uq_op(mm,bf_uq(naorbd)) .EQ. 0) EXIT
jkb = jkb + 1
! WRITE(*,'(a,2i4)')"C LB, KK, IVARB =",kk,ivarb
! WRITE(*,'(a,2i4)')"C LB, IAND =",iand(kk,ivarb)
! WRITE(*,'(a,2i4)')"C LB, parbit =",parbit(iand(kk,ivarb))
lorb = lex_least(bf_uq(naorbd))
! WRITE(*,'(a,1i4)')"C LB, aorb =",bf_uq(naorbd)
! WRITE(*,'(a,1i4)')"C LB, lorb =",lorb
itran(nsorb,jkb)=
& bf_op_map(kk+1,lorb)*parbit(IAND(mm-1,ivarb))
! & bf_op_map(mm,bf_uq(naorbd))*parbit(IAND(kk,ivarb))
! & bf_op_map(mm,lorb)*bf_phase(mm,lorb)*parbit(IAND(kk,ivarb))
! WRITE(*,*)"C LB, itran =",itran(nsorb,jkb)
ENDDO !mm (operators)
wt(nsorb,1)=nsorb
wt(nsorb,2)=jkb
! ELSE
! naorbd = naorbd + mult(stabilizer(cent))
ENDIF
ENDDO !k (orbitals)
ENDDO !j (unique shells)
ENDDO !ii (irreps)
WRITE(*,'(a)')" "
WRITE(*,'(a)')"Symmetrized Orbitals"
WRITE(*,'(a)')"--------------------"
WRITE(*,'(a)')" Sym Orb Atomic Orbs"
k = 0
DO j = 1, nir
WRITE(*,'(a)')" "
WRITE(*,'(a,a)')"Group ",zir(j)
DO i = 1, nbpsy(j)
k = k + 1
WRITE(*,'(a,1i4,a,8i4)')' ',k,' ',
& itran(k,1:wt(k,2))
ENDDO ! i (sym orbitals)
ENDDO ! j (irreps)
RETURN
END !sym_adapt
!-----------------------------------------------------------------
INTEGER FUNCTION lex_least(orb)
! Find the lexically lowest equivalent orbital
USE nwc_sym, ONLY : mxbf, mxshel, bf_uq, bf_op_map, nop
IMPLICIT NONE
! IN
INTEGER, INTENT(IN) :: orb ! orbital in question
! LOCAL
INTEGER :: k ! counter
INTEGER :: lorb !marks where lowest eq orb is in the map
lorb = 1
DO k = 1, nop
IF (bf_op_map(k,orb) .EQ. 0) EXIT
IF (bf_op_map(k,orb) .LT. bf_op_map(lorb,orb)) THEN
lorb = k
ENDIF
ENDDO ! k
lex_least = bf_op_map(lorb,orb)
RETURN
END !FUNCTION lex_least
!------------------------------------------------------------------
INTEGER FUNCTION ireplm(L,M)
! Taken from DALTON (hersol.F)
! Symmetry of RLM integrals for spherical basis
USE nwc_sym, ONLY: isymax
IMPLICIT NONE
!#include "nwc_const.fh"
!#include "sym_adapt.fh"
! IN
INTEGER, INTENT(IN) :: L, M
ireplm = 0
IF (MOD(L + M,2) .EQ. 1) ireplm = isymax(3,1)
IF (MOD(ABS(M),2) .EQ. 1) ireplm = IEOR(ireplm,isymax(1,1))
IF (M.LT.0) ireplm=IEOR(ireplm,IEOR(isymax(1,1),isymax(2,1)))
RETURN
END !FUNCTION ireplm
!------------------------------------------------------------------
SUBROUTINE lmnval(nhkta,khkta,l,m,n)
! Taken from DALTON (hergam.F and carpow.F)
! determine l, m, n for cartesian basis
USE nwc_sym, ONLY : mxaqnm, mxqnm
IMPLICIT NONE
!#include "nwc_const.fh"
!#include "sym_adapt.fh"
! IN
INTEGER, INTENT(IN) :: nhkta, khkta
! OUT
INTEGER, INTENT(OUT) :: l(khkta),m(khkta),n(khkta)
! LOCAL
INTEGER :: icomp,istep(mxaqnm),mval(mxaqnm),nval(mxaqnm)
INTEGER :: i,j,ij ! counters
! Calculate cartesian powers
ij=0
DO i=1, mxqnm
DO j=1, i
ij=ij + 1
istep(ij) = i
mval(ij) = i - j
nval(ij) = j - 1
ENDDO
ENDDO
icomp = 0
DO i = 1, khkta
icomp = icomp + 1
l(icomp) = nhkta - istep(i)
m(icomp) = mval(i)
n(icomp) = nval(i)
ENDDO
RETURN
END ! lmnval
!------------------------------------------------------------------
SUBROUTINE stab_cent(geom,ncent)
! Taken from DALTON (herrdn.F)
! Map stabilizers of each atom center
USE nwc_sym, ONLY: ngen, isymax, stabilizer, coords
IMPLICIT NONE
!#include "nwc_const.fh"
!#include "sym_adapt.fh"
! IN
INTEGER, INTENT(IN) :: ncent ! number of unique atom centers
INTEGER, INTENT(IN) :: geom ! geom handle
! LOCAL
INTEGER :: i,l,ll,m ! counters
INTEGER :: mul
! --- Loop over atoms
DO i = 1, ncent
mul = 0
ll = 1
! --- Loop over generating operators
DO 10 l = 1, ngen
! WRITE(*,*)"C LB, L=",L
DO m = 1, 3
! WRITE(*,*)"C LB, M=",M
! WRITE(*,*)"C LB, LL=",LL
! WRITE(*,*)"C LB, isymax=",isymax(m,1)
IF(IAND(ll,isymax(m,1)).NE.0) THEN
! WRITE(*,*)"C LB, i,coords=",i,coords(m,i,geom)
IF(ABS(coords(m,i,geom)).GE.1D-6)
& GOTO 10
ENDIF
ENDDO
mul = mul + ll
! WRITE(*,*)"C LB, mul=",mul
10 ll = 2 * ll
stabilizer(i) = mul
! WRITE(*,*)"C LB, i, stabilizer=",i,stabilizer(i)
ENDDO
RETURN
END ! stab_cent
!-----------------------------------------------------------------
SUBROUTINE sym_1int(num,numints,intbuf,labbuf,val,lab)
! Calculate the SOs from the AOs
USE nwc_sym, ONLY: nir, bf_per_ir_cum, nbpsy, wt, itran
IMPLICIT NONE
!#include "nwc_const.fh"
!#include "sym_adapt.fh"
! IN
INTEGER :: num !# of values to place in packed buffer
DOUBLE PRECISION :: intbuf(1:num) ! AO integer buffer
INTEGER :: labbuf(1:2,1:num) ! AO label buffer
! LOCAL
INTEGER :: i,j,k,l,m !counters
INTEGER :: offset
INTEGER :: i_lab,j_lab ! SO integral labels
INTEGER :: ao_i_lab,ao_j_lab ! AO integral labels
INTEGER :: sgni, sgnj ! sign of AO
! OUT
INTEGER :: numints ! number of SO ints
DOUBLE PRECISION :: val(1:num) ! SO integrals
INTEGER :: lab(1:2,1:num) ! SO labels
WRITE(*,*)"C LB, now in sym_1int"
numints = 0
! irreps
DO i = 1, nir
WRITE(*,'(a,2i4)')"C LB, irrep =",i
offset = bf_per_ir_cum(i-1)
WRITE(*,'(a,2i4)')"C LB, offset =",offset
! SO integral labels
DO i_lab = offset+1, offset+nbpsy(i)
DO j_lab = i_lab,offset+nbpsy(i)
! WRITE(*,'(a,2i4)')" C LB, SO integral label = ",i_lab,j_lab
numints = numints + 1
val(numints) = 0.0D0
lab(1,numints) = i_lab
lab(2,numints) = j_lab
! AO linear combination for i_lab
DO j = 1, wt(i_lab,2)
ao_i_lab = itran(i_lab,j)
sgni = SIGN(1,ao_i_lab)
ao_i_lab = ABS(ao_i_lab)
! AO linear combination for j_lab
DO k = 1, wt(j_lab,2)
ao_j_lab = itran(j_lab,k)
sgnj = SIGN(1,ao_j_lab)
ao_j_lab = ABS(ao_j_lab)
! WRITE(*,'(a,2i4)')" C LB, AO integral label = ",
! & ao_i_lab,ao_j_lab
! Find the AO integral that corresponds to ao_i_lab
! ao_j_lab, multiply by sgni*sgnj, and add to val
DO l = 1, num
! WRITE(*,'(a,2i4)')" C LB, labbuf ="
! & ,labbuf(1:2,l)
IF (((labbuf(1,l) .EQ. ao_i_lab) .AND.
& (labbuf(2,l) .EQ. ao_j_lab)) .OR.
& ((labbuf(1,l) .EQ. ao_j_lab) .AND.
& (labbuf(2,l) .EQ. ao_i_lab))) THEN
! match found!
WRITE(*,'(a,1e15.3)')" C LB, integral found:",
& intbuf(l)
val(numints) = val(numints) + intbuf(l)*sgni*sgnj
ENDIF
ENDDO ! l
ENDDO !k
ENDDO !j
! WRITE(*,'(a,3i4,1e15.3)')" C LB, SO integral:",
! & numints,lab(1:2,numints),val(numints)
ENDDO !j_lab (SO lab j)
ENDDO !i_lab (SO lab i)
ENDDO !i (irreps)
RETURN
END ! sym_1int

21
src/nwc_columbus/aoints/wrt_dft_aoints.F
View file

@ -79,6 +79,9 @@ c sifs
cgk debug
integer mem1, max1e, mem2, max2e
cgk end
C LB
INTEGER :: k_SOlab, l_SOlab
C LB
ierr=0 !initialize
imtype=0
@ -458,10 +461,16 @@ cgk end
$ call errquit('wrt_dft_aoints: ma failed for map',
& szintarr*nmap, MA_ERR)
c get memory for labels
write(*,*)'LB before MA_push_get, szlabs=', szlabs
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)
c get memory for Sym Orb labels
if (.not. MA_push_get(MT_int,4*szlabs,'int_1e_sifs:SOlab',
$ l_SOlab,k_SOlab))
$ call errquit('int_1e_sifs: error getting SOlab',
& 2*szlabs, MA_ERR)
call schwarz_tidy()
call int_terminate()
@ -469,6 +478,10 @@ c get memory for labels
* call schwarz_init(geom, AO_bas_han)
* call int_terminate()
* call intd_init(rtdb, 1, AO_bas_han)
C LB
! Calculate the symmetrized orbital transition matrix
CALL sym_adapt(AO_bas_han)
C LB
cgk debug
* write(*,*)'gk: calling int_init'
cgk end
@ -479,7 +492,7 @@ 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: befor int_1e_sifs, mem1 = ', mem1
* write(*,*)'gk: mem2 = ', mem2
* write(*,*)'gk: max1e = ', max1e
* write(*,*)'gk: max2e = ', max2e
@ -489,7 +502,7 @@ cgk end
cgk debug
! write(*,*)'gk: calling int_1e_sifs'
! write(*,*)'gk: mxengy=', mxengy
! write(*,*)'LB before int_1e_sifs, info=', info(1:5)
write(*,*)'LB before int_1e_sifs, szlabs=', szlabs
! WRITE(*,*)"k_map=",k_map,"szintarr=",szintarr,"nmap=",nmap
cgk write
@ -497,7 +510,7 @@ cgk write
call int_1e_sifs(AO_bas_han, aoints, energy, mxengy,
& nbf_ao, nmap, int_mb(k_map),imtype,
& ibvtyp, ibitv, l1rec, n1max, int_mb(k_lab),
& ninfo, info)
& ninfo, info, int_mb(k_SOlab))
cgk debug
! write(*,*)'gk: back from int_1e_sifs'
@ -531,6 +544,8 @@ cgk end
c
c done with memory for sifs files
c
if (.not. MA_pop_stack(l_SOlab)) call errquit
$ ('int_1e_sifs: pop failed at l_SOlab', 0, GA_ERR)
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

10
src/task/task_aoints.F
View file

@ -31,7 +31,7 @@ cc
logical oecce
cc
cgk debug
* write(*,*)'gk: entered task_aoints'
c write(*,*)'gk: entered task_aoints'
cgk end
oecce = .false.
cpu = util_cpusec()
@ -46,7 +46,7 @@ c
$ call errquit('task:aoints: theory not specified',0,
$ RTDB_ERR)
if (theory .ne. 'sodft') then
if ((theory .ne. 'sodft').and.(theory .ne. 'dft')) then
call errquit('task_aoints: theory not supported for aoints'
& ,0,RTDB_ERR)
else
@ -194,7 +194,7 @@ c
* integer which_mp2
c
cgk debug
* write(*,*)'gk: in task_aoints_doit'
* write(*,*)'gk: in task_aoints_doit, theory=',theory
cgk end
call ecce_print_module_entry('task energy')
c
@ -392,8 +392,8 @@ cgk end
* else if (theory .eq. 'dangchang') then
* status = dc_energy(rtdb)
* hack below
if (theory .eq. 'sodft') then
status = sodft_aoints(rtdb)
if ((theory .eq. 'sodft').OR.(theory .eq. 'dft')) then
status = sodft_aoints(rtdb,theory)
else
call errquit('task_aoints: unknown theory',0, INPUT_ERR)
endif