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