1C> \ingroup task 2C> @{ 3 logical function task_ecp_print_integrals(rtdb) 4 implicit none 5#include "errquit.fh" 6* $Id$ 7c 8#include "stdio.fh" 9#include "mafdecls.fh" 10#include "geom.fh" 11#include "bas.fh" 12c 13 logical int_normalize, raktask26_a 14 external int_normalize, raktask26_a 15c 16c 17 integer rtdb 18c 19 logical status 20 integer basis, geom 21 integer nbf, nat, nshell 22 integer maxg1, maxs1 23 integer hbuf, hscr, hint 24 integer kbuf, kscr, kint 25c 26 task_ecp_print_integrals = .false. 27c 28 if (.not.geom_create(geom,'geometry')) call errquit 29 & ('geom create failed',911, GEOM_ERR) 30 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit 31 & ('geom_rtdb_load failed',911, RTDB_ERR) 32c 33 if (.not.bas_create(basis,'ao basis')) call errquit 34 & ('bas_create failed',911, BASIS_ERR) 35 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) call errquit 36 & ('bas_rtdb_load failed',911, RTDB_ERR) 37c 38 write(6,*)' geom/basis loaded' 39c 40 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 41c 42 if (.not. bas_print(basis)) 43 $ call errquit(' basis print failed', 0, BASIS_ERR) 44c 45 if (.not.bas_numbf(basis,nbf)) call errquit 46 & ('numbf failed',911, BASIS_ERR) 47 if (.not.bas_numcont(basis,nshell)) call errquit 48 & ('numcont failed',911, BASIS_ERR) 49c 50 if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe' 51 write(6,*) 'number of atoms ', nat 52c 53 call int_init(rtdb, 1, basis) 54 call int_mem_1e(maxg1, maxs1) 55 write(luout,*)' maxg1 = ',maxg1 56 write(luout,*)' maxs1 = ',maxs1 57 maxs1 = max(maxs1,(nbf*nbf)) 58 write(luout,*)' maxs1 = ',maxs1, ' after max for copy ' 59 status = .true. 60 status = status .and. 61 & ma_alloc_get(mt_dbl,maxg1,'int buffer' ,hbuf,kbuf) 62 status = status .and. 63 & ma_alloc_get(mt_dbl,maxs1,'scr buffer' ,hscr,kscr) 64 status = status .and. 65 & ma_alloc_get(mt_dbl,(nbf*nbf),'ints' ,hint,kint) 66c 67 task_ecp_print_integrals = raktask26_a(rtdb, 68 & geom, basis, nbf, nat, nshell, maxg1, maxs1, 69 & dbl_mb(kint), 70 & dbl_mb(kbuf), 71 & dbl_mb(kscr)) 72 status = .true. 73 status = status.and.ma_free_heap(hint) 74 status = status.and.ma_free_heap(hbuf) 75 status = status.and.ma_free_heap(hscr) 76 status = status.and.bas_destroy(basis) 77 status = status.and.geom_destroy(geom) 78 task_ecp_print_integrals = task_ecp_print_integrals.and.status 79 call int_terminate() 80 end 81C> @} 82 logical function raktask26_a(rtdb, 83 & geom, basis, nbf, nat, nshell, szb, szs, 84 & zint,buf,scr) 85 implicit none 86#include "errquit.fh" 87#include "nwc_const.fh" 88#include "mafdecls.fh" 89#include "geom.fh" 90#include "geomP.fh" 91#include "basdeclsP.fh" 92#include "basP.fh" 93#include "bas.fh" 94#include "stdio.fh" 95#include "geobasmapP.fh" 96#include "inp.fh" 97#include "bas_exndcf_dec.fh" 98#include "bas_ibs_dec.fh" 99c 100 integer rtdb, geom, basis, nbf, nat, nshell, szb, szs 101 double precision zint(nbf,nbf) 102 double precision buf(szb) 103 double precision scr(szs) 104c 105 double precision val 106 integer nshell_use 107 integer ishell, ilo, ihi, nbfshi 108 integer jshell, jlo, jhi, nbfshj 109 integer ii_np, ii_gen, ii_exp, ii_cf, ii_type, ii_atom 110 integer jj_np, jj_gen, jj_exp, jj_cf, jj_type, jj_atom 111 integer nbfsh, ucont, xbas, cnt, i, j, lu 112 logical does_it_exist 113c 114 character*255 filename 115c 116#include "bas_exndcf_sfn.fh" 117#include "bas_ibs_sfn.fh" 118c 119 call dfill (nbf*nbf,0.0d00,zint,1) 120 call dfill (szb,0.0d00,buf,1) 121 call dfill (szs,0.0d00,scr,1) 122c 123 xbas = basis + BASIS_HANDLE_OFFSET 124 nshell_use = nshell 125c 126 do ishell = 1, nshell_use 127 write(6,*)' fd: ishell = ',ishell,' of ',nshell_use 128 call util_flush(6) 129 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 130 & stop 'cn2bfr error i' 131 nbfshi = ihi - ilo + 1 132 ucont = (sf_ibs_cn2ucn(ishell,xbas)) 133 ii_np = infbs_cont(CONT_NPRIM,ucont,xbas) 134 ii_gen = infbs_cont(CONT_NGEN,ucont,xbas) 135 ii_exp = infbs_cont(CONT_IEXP,ucont,xbas) 136 ii_cf = infbs_cont(CONT_ICFP,ucont,xbas) 137 ii_type = infbs_cont(CONT_TYPE,ucont,xbas) 138 ii_atom = (sf_ibs_cn2ce(ishell,xbas)) 139 do jshell = 1, ishell 140 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 141 & stop 'cn2bfr error j' 142 nbfshj = jhi - jlo + 1 143 nbfsh = nbfshi*nbfshj 144* write(6,*)' fd: jshell = ',jshell,' size =',nbfsh 145 ucont = (sf_ibs_cn2ucn(jshell,xbas)) 146 jj_np = infbs_cont(CONT_NPRIM,ucont,xbas) 147 jj_gen = infbs_cont(CONT_NGEN,ucont,xbas) 148 jj_exp = infbs_cont(CONT_IEXP,ucont,xbas) 149 jj_cf = infbs_cont(CONT_ICFP,ucont,xbas) 150 jj_type = infbs_cont(CONT_TYPE,ucont,xbas) 151 jj_atom = (sf_ibs_cn2ce(jshell,xbas)) 152* 153 call dfill (szb,0.0d00,buf,1) 154 call dfill (szs,0.0d00,scr,1) 155 call int_ecp_hf1( 156 & coords(1,ii_atom,geom), 157 & dbl_mb(mb_exndcf(ii_exp,xbas)), 158 & dbl_mb(mb_exndcf(ii_cf,xbas)), 159 & ii_np, ii_gen, ii_type, ii_atom, 160 161 & coords(1,jj_atom,geom), 162 & dbl_mb(mb_exndcf(jj_exp,xbas)), 163 & dbl_mb(mb_exndcf(jj_cf,xbas)), 164 & jj_np, jj_gen, jj_type, jj_atom, 165 166 & buf,nbfsh,scr,szs,.false.) 167*-------- 168 cnt = 1 169 do i = ilo,ihi 170 do j = jlo, jhi 171 zint(i,j) = buf(cnt) 172 zint(j,i) = buf(cnt) 173 enddo 174 enddo 175*-------- 176 enddo 177 enddo 178* now have full ecp integral matrix print it out 179c 180 call util_file_name('ecp_integrals',.false.,.false.,filename) 181c 182 lu = 69 183 does_it_exist = .false. 184 inquire(file=filename,exist=does_it_exist) 185 if (does_it_exist) then 186 write(luout,*) 187 & 'rak26: overwrite of existing file', 188 & filename(1:inp_strlen(filename)) 189 call util_file_unlink(filename) 190 endif 191* 192 does_it_exist = .false. 193 inquire(file=filename,exist=does_it_exist) 194 if (does_it_exist) then 195 write(luout,*) 196 & 'rak26: file not unlinked: ', 197 & filename(1:inp_strlen(filename)) 198 call errquit('rak26: fatal error ',911, DISK_ERR) 199 endif 200* 201 open(unit=lu,file=filename, 202 & form='formatted', 203 & access='sequential', 204 & status='new') 205 cnt = 0 206 do i = 1,nbf 207 do j = 1,i 208 cnt = cnt + 1 209 val = zint(i,j) 210 if (abs(val).gt.1.0d-8) then 211 write(lu,10000)cnt,i,j,val 212 endif 213 enddo 214 enddo 215 close(unit=lu,status='keep') 216 raktask26_a = .true. 21710000 format(1x,i10,i5,i5,1pd20.10) 218 end 219 220 221