From 60c46833411e9131d7c5f71c18bc51fdac8b0af6 Mon Sep 17 00:00:00 2001 From: "lachlan.belcher" Date: Tue, 8 Oct 2024 10:54:43 -0300 Subject: [PATCH 1/2] addition of files sym_adapt.F and sym_adapt.fh to begin the process of symmetry adaptation for Columbus. So far, this is the calculation of the AO to SO matrix. --- src/include/sym_adapt.fh | 59 +++ src/nwc_columbus/aoints/GNUmakefile | 2 +- src/nwc_columbus/aoints/dft_aoints.F | 11 +- src/nwc_columbus/aoints/int_1e_sifs.F | 129 ++++-- src/nwc_columbus/aoints/sym_adapt.F | 518 +++++++++++++++++++++++ src/nwc_columbus/aoints/wrt_dft_aoints.F | 4 + src/task/task_aoints.F | 10 +- 7 files changed, 695 insertions(+), 38 deletions(-) create mode 100644 src/include/sym_adapt.fh create mode 100644 src/nwc_columbus/aoints/sym_adapt.F diff --git a/src/include/sym_adapt.fh b/src/include/sym_adapt.fh new file mode 100644 index 0000000000..ec55dde653 --- /dev/null +++ b/src/include/sym_adapt.fh @@ -0,0 +1,59 @@ +#include "nwc_const.fh" +#include "geomP.fh" + 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 :: IS(0:7) + INTEGER :: parbit(0:7) + INTEGER :: itran(mxbf,mxbf) + INTEGER :: irrepmap(0:mxsym-1) !map from Dalton to NWChem irrep order + INTEGER :: opmap(0:mxsym-1) !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 operators other than E + 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 + + +! DATA IS/0,1,1,2,1,2,2,3/ + DATA parbit/1,-1,-1,1,-1,1,1,-1/ + DATA irrepmap/0,1,2,3,4,5,6,7/ ! will be changed in sym_adapt depending on group + DATA opmap/0,1,2,3,4,5,6,7/ ! will be changed in sym_adapt depending on group + + COMMON/sym_ad/ isymax,sh_op_map,bf_op_map, + & 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,bf_per_ir, + & bf_per_ir_cum,bf_so_ir,so_uq_to_lab, + & so_lab1,so_lab2,bf_phase,char_tab, +! & stabilizer, ngen,is,parbit,itran, irrepmap, + & stabilizer, ngen,parbit,itran, irrepmap, + & wt,nbpsy, nop, nir diff --git a/src/nwc_columbus/aoints/GNUmakefile b/src/nwc_columbus/aoints/GNUmakefile index 65e58d5b24..4a4ae1d627 100644 --- a/src/nwc_columbus/aoints/GNUmakefile +++ b/src/nwc_columbus/aoints/GNUmakefile @@ -11,7 +11,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 diff --git a/src/nwc_columbus/aoints/dft_aoints.F b/src/nwc_columbus/aoints/dft_aoints.F index 3bf7a71cd6..026a427428 100644 --- a/src/nwc_columbus/aoints/dft_aoints.F +++ b/src/nwc_columbus/aoints/dft_aoints.F @@ -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 @@ logical function sodft_aoints(rtdb) #include "util.fh" #include "cgridfile.fh" #include "cosmo.fh" + +#include "geomP.fh" c c local declarations c @@ -67,7 +69,7 @@ logical function sodft_aoints(rtdb) 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 @@ logical function sodft_aoints(rtdb) 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 @@ logical function sodft_aoints(rtdb) 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 @@ logical function sodft_aoints(rtdb) 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' diff --git a/src/nwc_columbus/aoints/int_1e_sifs.F b/src/nwc_columbus/aoints/int_1e_sifs.F index 92383d6827..31fb843a27 100644 --- a/src/nwc_columbus/aoints/int_1e_sifs.F +++ b/src/nwc_columbus/aoints/int_1e_sifs.F @@ -14,6 +14,8 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, #include "cscfps.fh" #include "sym.fh" #include "cdft.fh" +!#include "geomP.fh" +#include "sym_adapt.fh" c c Compute the desired type of integrals (kinetic, potential, overlap) c @@ -52,14 +54,15 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, c integer aoints integer n1max - INTEGER l1rec, ntitle, ibuf, nsym, ibvtyp, ierr + INTEGER l1rec, ntitle, ibuf, nsym, ibvtyp, ierr!,mxsym 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) @@ -68,8 +71,8 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, integer clab(2,*) real*8 energy(nenrgy) character*80 stitle(ntitle) - character*4 slabel(nsym) - integer nbpsy(nsym) + character*4 slabel(mxsym) + !integer nbpsy(mxsym) character*1 shtypes(-1:7) data shtypes/'l','s', 'p', 'd', 'f', 'g', 'h', 'i', 'k'/ integer shmap(-1:7) @@ -83,6 +86,32 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, real*8 fcore logical basok +C LB + !logical status + !character*20 zname + !character*8 zir(mxsym), zclass(mxsym) + !integer iclass(mxsym), nop, nir + !double precision chars(mxsym*mxsym) + !character*8 grp_name + + !INTEGER, PARAMETER :: mxshel=120 + !integer sh_op_map(8,mxshel), bf_op_map(8,mxbf) + !double precision bf_phase(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) + !double precision char_tab(1:8,0:7) + !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) + !logical oprint + +C LB cgk provisional * change thresh to the appropriate user supplied zero tolerance @@ -125,6 +154,7 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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 @@ -209,12 +239,22 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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) +! WRITE(*,*)"C LB, ibas=",ibas +! WRITE(*,*)"C LB, ishell=",ishell +! WRITE(*,*)"C LB, otype=",otype +! WRITE(*,*)"C LB, nprim=",nprim +! WRITE(*,*)"C LB, ngen=",ngen +! 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) @@ -249,9 +289,24 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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) + nsym=nir + 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,7 +327,8 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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 ) cgk debug @@ -302,18 +358,20 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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 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,32 +382,46 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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 + else if (type .eq. 2) then !Potential Energy 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 + 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 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: numints = ', numints ! write(*,*)'gk: back from int_l1epe' cgk end - else if (type .eq. 3) then + else if (type .eq. 3) then !Overlap +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 +cgk end + 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, @@ -357,6 +429,7 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, end if c cgk debug + write(*,*)'gk: numints = ', numints * write(*,*)'gk: this batch, numints = ', numints cgk end diff --git a/src/nwc_columbus/aoints/sym_adapt.F b/src/nwc_columbus/aoints/sym_adapt.F new file mode 100644 index 0000000000..feac1149ae --- /dev/null +++ b/src/nwc_columbus/aoints/sym_adapt.F @@ -0,0 +1,518 @@ + SUBROUTINE sym_adapt(ibas) + ! Routine to adapt the AOs to symmetrized orbitals + IMPLICIT NONE +#include "errquit.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 :: nop, nir, 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) + 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) + + ! 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+1))/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)) + itran(nsorb,jkb)= + & bf_op_map(mm,bf_uq(naorbd))*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 ireplm(L,M) + ! Taken from DALTON (hersol.F) + ! Symmetry of RLM integrals for spherical basis + + IMPLICIT NONE +#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 ! ireplm + + + + + +!------------------------------------------------------------------ + SUBROUTINE lmnval(nhkta,khkta,l,m,n) + ! Taken from DALTON (hergam.F and carpow.F) + ! determine l, m, n for cartesian basis + + IMPLICIT NONE +#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 + + IMPLICIT NONE +#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 diff --git a/src/nwc_columbus/aoints/wrt_dft_aoints.F b/src/nwc_columbus/aoints/wrt_dft_aoints.F index 226407a2f2..70fc121f01 100644 --- a/src/nwc_columbus/aoints/wrt_dft_aoints.F +++ b/src/nwc_columbus/aoints/wrt_dft_aoints.F @@ -469,6 +469,10 @@ logical function wrt_dft_aoints(rtdb) * 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 diff --git a/src/task/task_aoints.F b/src/task/task_aoints.F index 71e1ce2f81..b7b615d7d0 100644 --- a/src/task/task_aoints.F +++ b/src/task/task_aoints.F @@ -31,7 +31,7 @@ logical function task_aoints(rtdb) 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 @@ logical function task_aoints(rtdb) $ 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 @@ logical function task_aoints_doit(rtdb,theory) * 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 @@ logical function task_aoints_doit(rtdb,theory) * 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 From 987b029c42f83f3863ccaf972838e7de9ee8850b Mon Sep 17 00:00:00 2001 From: "lachlan.belcher" Date: Mon, 11 Nov 2024 22:37:16 -0300 Subject: [PATCH 2/2] corrected indexing error for dipole integrals; starting to work towards abelian symmetry; currenty only identifies the linear combinations that make up the symmetry orbitals, but does not implement. CONTINUE TO USE ONLY C1 SYMMETRY! --- src/nwc_columbus/aoints/GNUmakefile | 3 +- src/nwc_columbus/aoints/egrad_trace.F | 3 +- src/nwc_columbus/aoints/int_1e_sifs.F | 188 ++++------- src/nwc_columbus/aoints/int_mom_sifs.F | 100 ++++-- src/nwc_columbus/aoints/int_so_sifs.F | 21 +- .../aoints/nwc_sym_mod.F} | 35 +- src/nwc_columbus/aoints/sym_adapt.F | 319 +++++++++++++----- src/nwc_columbus/aoints/wrt_dft_aoints.F | 17 +- 8 files changed, 411 insertions(+), 275 deletions(-) rename src/{include/sym_adapt.fh => nwc_columbus/aoints/nwc_sym_mod.F} (62%) diff --git a/src/nwc_columbus/aoints/GNUmakefile b/src/nwc_columbus/aoints/GNUmakefile index 4a4ae1d627..2398d1011f 100644 --- a/src/nwc_columbus/aoints/GNUmakefile +++ b/src/nwc_columbus/aoints/GNUmakefile @@ -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 \ diff --git a/src/nwc_columbus/aoints/egrad_trace.F b/src/nwc_columbus/aoints/egrad_trace.F index cfb48cd235..039d8477a3 100644 --- a/src/nwc_columbus/aoints/egrad_trace.F +++ b/src/nwc_columbus/aoints/egrad_trace.F @@ -415,7 +415,8 @@ subroutine egrad_trace(ibas, aodens, ninfo, info, nbft, g_force, !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) diff --git a/src/nwc_columbus/aoints/int_1e_sifs.F b/src/nwc_columbus/aoints/int_1e_sifs.F index 31fb843a27..ef3db35ede 100644 --- a/src/nwc_columbus/aoints/int_1e_sifs.F +++ b/src/nwc_columbus/aoints/int_1e_sifs.F @@ -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" @@ -14,8 +15,6 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, #include "cscfps.fh" #include "sym.fh" #include "cdft.fh" -!#include "geomP.fh" -#include "sym_adapt.fh" c c Compute the desired type of integrals (kinetic, potential, overlap) c @@ -34,14 +33,14 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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,10 +53,10 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, c integer aoints integer n1max - INTEGER l1rec, ntitle, ibuf, nsym, ibvtyp, ierr!,mxsym + 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 @@ -69,15 +68,15 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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(mxsym) - !integer nbpsy(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) @@ -87,30 +86,7 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, logical basok C LB - !logical status - !character*20 zname - !character*8 zir(mxsym), zclass(mxsym) - !integer iclass(mxsym), nop, nir - !double precision chars(mxsym*mxsym) - !character*8 grp_name - - !INTEGER, PARAMETER :: mxshel=120 - !integer sh_op_map(8,mxshel), bf_op_map(8,mxbf) - !double precision bf_phase(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) - !double precision char_tab(1:8,0:7) - !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) - !logical oprint - + INTEGER :: nsoints ! number of symmetry orbital integrals C LB cgk provisional @@ -123,7 +99,6 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, c cgk debug write(*,*)'gk: **** entered int_1e_sifs' -! WRITE(*,*)'LB: info=',info(1:5) cgk end basok=.false. odbug=.true. @@ -201,7 +176,11 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, & 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 @@ -243,12 +222,12 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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, ngen=",ngen +! WRITE(*,*)"C LB, ngeno=",ngeno ! WRITE(*,*)"C LB, sphcart=",sphcart basok=bas_cn2bfr(ibas, ishell, ilo, ihi) ! WRITE(*,*)"C LB, ilo=",ilo @@ -261,13 +240,13 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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 @@ -291,19 +270,23 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, c slabel(1)='a1' C LB: Add symmetry labels to aoints files - !status = sym_char_table(grp_name,nop,nir,iclass,zir,zclass,chars) +! 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 +! 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 +! 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 @@ -331,6 +314,7 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, & 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) @@ -362,14 +346,14 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, cgk end do jshell = 1, nshell - do ishell = 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) +! & odoit = sym_shell_pair(ibas, ishell, jshell, q2) ! WRITE(*,*)"C LB, odoit= ", odoit ! WRITE(*,*)"C LB, q2= ", q2 @@ -382,71 +366,40 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, call ifill(max1e,0,int_mb(k_ilab),1) call ifill(max1e,0,int_mb(k_jlab),1) - if (type .eq. 1) then !Kinetic Energy + 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 !Potential Energy -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 -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: numints = ', numints -! write(*,*)'gk: back from int_l1epe' -cgk end - else if (type .eq. 3) then !Overlap -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 -cgk end + 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: numints = ', numints -* 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)) @@ -458,32 +411,23 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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, @@ -500,7 +444,7 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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' @@ -508,13 +452,15 @@ subroutine int_1e_sifs(ibas, aoints, energy, nenrgy, nbft, 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 diff --git a/src/nwc_columbus/aoints/int_mom_sifs.F b/src/nwc_columbus/aoints/int_mom_sifs.F index ad8f5d1575..3bc9536d52 100644 --- a/src/nwc_columbus/aoints/int_mom_sifs.F +++ b/src/nwc_columbus/aoints/int_mom_sifs.F @@ -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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, #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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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 @@ subroutine int_mom_sifs(ibas, oskel, aoints, nbft, 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' diff --git a/src/nwc_columbus/aoints/int_so_sifs.F b/src/nwc_columbus/aoints/int_so_sifs.F index bb62a93cb8..750129ef51 100644 --- a/src/nwc_columbus/aoints/int_so_sifs.F +++ b/src/nwc_columbus/aoints/int_so_sifs.F @@ -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 @@ subroutine int_so_sifs(ibas, oskel, aoints, nbft, max1e, mem1, 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 @@ subroutine int_so_sifs(ibas, oskel, aoints, nbft, max1e, mem1, 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 @@ subroutine int_so_sifs(ibas, oskel, aoints, nbft, max1e, mem1, 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 @@ subroutine int_so_sifs(ibas, oskel, aoints, nbft, max1e, mem1, * 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' diff --git a/src/include/sym_adapt.fh b/src/nwc_columbus/aoints/nwc_sym_mod.F similarity index 62% rename from src/include/sym_adapt.fh rename to src/nwc_columbus/aoints/nwc_sym_mod.F index ec55dde653..d0a42c3ff5 100644 --- a/src/include/sym_adapt.fh +++ b/src/nwc_columbus/aoints/nwc_sym_mod.F @@ -1,5 +1,11 @@ -#include "nwc_const.fh" -#include "geomP.fh" + 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 @@ -24,14 +30,13 @@ INTEGER :: isymao(mxqn,mxaqn) INTEGER :: stabilizer(max_cent) !stabilizer for each center INTEGER :: ngen ! number of generating operators -! INTEGER :: IS(0:7) - INTEGER :: parbit(0:7) + 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) !map from Dalton to NWChem irrep order - INTEGER :: opmap(0:mxsym-1) !map from Dalton to NWChem operator order + 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 operators other than E + INTEGER :: nop !number of symmetry operators INTEGER :: nir !number of irreps LOGICAL :: oprint @@ -43,17 +48,5 @@ CHARACTER*8 :: zir(mxsym) ! irrep names -! DATA IS/0,1,1,2,1,2,2,3/ - DATA parbit/1,-1,-1,1,-1,1,1,-1/ - DATA irrepmap/0,1,2,3,4,5,6,7/ ! will be changed in sym_adapt depending on group - DATA opmap/0,1,2,3,4,5,6,7/ ! will be changed in sym_adapt depending on group - - COMMON/sym_ad/ isymax,sh_op_map,bf_op_map, - & 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,bf_per_ir, - & bf_per_ir_cum,bf_so_ir,so_uq_to_lab, - & so_lab1,so_lab2,bf_phase,char_tab, -! & stabilizer, ngen,is,parbit,itran, irrepmap, - & stabilizer, ngen,parbit,itran, irrepmap, - & wt,nbpsy, nop, nir + + END MODULE nwc_sym diff --git a/src/nwc_columbus/aoints/sym_adapt.F b/src/nwc_columbus/aoints/sym_adapt.F index feac1149ae..e1b5a09081 100644 --- a/src/nwc_columbus/aoints/sym_adapt.F +++ b/src/nwc_columbus/aoints/sym_adapt.F @@ -1,8 +1,10 @@ SUBROUTINE sym_adapt(ibas) ! Routine to adapt the AOs to symmetrized orbitals + USE nwc_sym IMPLICIT NONE #include "errquit.fh" -#include "sym_adapt.fh" +!#include "nwc_const.fh" +!#include "sym_adapt.fh" #include "sym.fh" #include "geom.fh" #include "bas.fh" @@ -16,12 +18,13 @@ SUBROUTINE sym_adapt(ibas) INTEGER :: l(mxaqn),m(mxaqn),n(mxaqn) INTEGER :: Lvar, Mvar, Nvar INTEGER :: ang_max !maximum L for basis ibas - INTEGER :: nop, nir, iclass(mxsym),otype + 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) @@ -50,6 +53,7 @@ SUBROUTINE sym_adapt(ibas) ! 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)) @@ -70,92 +74,92 @@ SUBROUTINE sym_adapt(ibas) $ 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 +! 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+1))/LOG(2.0)) + 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 @@ -357,7 +361,7 @@ SUBROUTINE sym_adapt(ibas) ! WRITE(*,*)" C LB, orbital contributes!" nsorb = nsorb + 1 ! WRITE(*,'(a,1i4)') -! &" C LB, this will be symmetrized orbital",nsorb +! &" C LB, this will be symmetrized orbital",nsorb jkb = 0 ! WRITE(*,'(a,1i4)')"C LB, naorbd=",naorbd @@ -372,10 +376,16 @@ SUBROUTINE sym_adapt(ibas) ! 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(mm,bf_uq(naorbd))*parbit(IAND(kk,ivarb)) + & 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) + ENDDO !mm (operators) wt(nsorb,1)=nsorb wt(nsorb,2)=jkb @@ -407,14 +417,45 @@ SUBROUTINE sym_adapt(ibas) +!----------------------------------------------------------------- + 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 "sym_adapt.fh" +!#include "nwc_const.fh" +!#include "sym_adapt.fh" ! IN INTEGER, INTENT(IN) :: L, M @@ -425,7 +466,7 @@ INTEGER FUNCTION ireplm(L,M) IF (M.LT.0) ireplm=IEOR(ireplm,IEOR(isymax(1,1),isymax(2,1))) RETURN - END ! ireplm + END !FUNCTION ireplm @@ -436,8 +477,10 @@ 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 "sym_adapt.fh" +!#include "nwc_const.fh" +!#include "sym_adapt.fh" ! IN INTEGER, INTENT(IN) :: nhkta, khkta @@ -476,8 +519,10 @@ 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 "sym_adapt.fh" +!#include "nwc_const.fh" +!#include "sym_adapt.fh" ! IN INTEGER, INTENT(IN) :: ncent ! number of unique atom centers @@ -516,3 +561,91 @@ SUBROUTINE stab_cent(geom,ncent) 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 diff --git a/src/nwc_columbus/aoints/wrt_dft_aoints.F b/src/nwc_columbus/aoints/wrt_dft_aoints.F index 70fc121f01..468b76e49d 100644 --- a/src/nwc_columbus/aoints/wrt_dft_aoints.F +++ b/src/nwc_columbus/aoints/wrt_dft_aoints.F @@ -79,6 +79,9 @@ logical function wrt_dft_aoints(rtdb) 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 @@ logical function wrt_dft_aoints(rtdb) $ 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() @@ -483,7 +492,7 @@ logical function wrt_dft_aoints(rtdb) 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 @@ -493,7 +502,7 @@ logical function wrt_dft_aoints(rtdb) 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 @@ -501,7 +510,7 @@ logical function wrt_dft_aoints(rtdb) 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' @@ -535,6 +544,8 @@ logical function wrt_dft_aoints(rtdb) 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