1c $Id$
2      logical function kgdtest (rtdb)
3      implicit none
4#include "mafdecls.fh"
5#include "rtdb.fh"
6#include "context.fh"
7#include "global.fh"
8#include "stdio.fh"
9c::passed
10      integer rtdb          ! rtdb handle
11c::local
12      integer me
13      integer kgdtask, kgd_tmp
14c
15      call ga_sync()
16      call ga_sync()
17      kgdtask = 0
18      if (rtdb_get(rtdb,'kgdtask',MT_INT,1,kgd_tmp))
19     &    kgdtask = kgd_tmp
20c
21      call ga_sync()
22      call ga_sync()
23      me = ga_nodeid()
24      write (luout,*) 'in kgdtest'
25      call util_flush(luout)
26c
27      if (kgdtask.eq.0) then    !...................................   0
28        if (me.eq.0) then
29          write(luout,*)' default kgdtest task '
30          write(luout,*)' test use of kgdtest! '
31        endif
32        kgdtest = .true.
33      else if (kgdtask.eq.1) then !.................................   1
34        if (me.eq.0) write(luout,*)
35     &      ' kgdtest task 1: relativistic integral test'
36        call kgdtest_rel1e(rtdb)
37        kgdtest = .true.
38      else if (kgdtask.eq.2) then !.................................   2
39        if (me.eq.0) write(luout,*)
40     &      ' kgdtest task 2: general contraction test'
41        call kgdtest_gencon(rtdb)
42        kgdtest = .true.
43      else if (kgdtask.eq.3) then !.................................   3
44        if (me.eq.0) write(luout,*)
45     &      ' kgdtest task 3: so-ecp integral test'
46        call kgdtest_soecp(rtdb)
47        kgdtest = .true.
48      else if (kgdtask.eq.4) then !.................................   4
49        if (me.eq.0) write(luout,*)
50     &      ' kgdtest task 4: rel2e integral test'
51        call kgdtest_rel2e(rtdb)
52        kgdtest = .true.
53      else if (kgdtask.eq.5) then !.................................   5
54        if (me.eq.0) write(luout,*)
55     &      ' kgdtest task 5: ecp memory test'
56        call kgdtest_ecpmem(rtdb)
57        kgdtest = .true.
58      else if (kgdtask.eq.6) then !.................................   6
59        if (me.eq.0) write(luout,*)
60     &      ' kgdtest task 6: rel general contraction integral test'
61        call kgdtest_relgen(rtdb)
62        kgdtest = .true.
63      end if
64      end
65************************************************************************
66      subroutine kgdtest_rel1e(rtdb)
67      implicit none
68#include "errquit.fh"
69#include "mafdecls.fh"
70#include "rtdb.fh"
71#include "context.fh"
72#include "geom.fh"
73#include "bas.fh"
74#include "nwc_const.fh"
75#include "basP.fh"
76#include "basdeclsP.fh"
77#include "geomP.fh"
78#include "geobasmapP.fh"
79#include "bas_exndcf_dec.fh"
80#include "bas_ibs_dec.fh"
81c
82      logical int_normalize
83      external int_normalize
84c
85c test relativistic one-electron integrals
86c
87      integer rtdb
88      integer geom,basis, basis_id
89      integer nshell, memscr, membuf
90      integer h_scr, k_scr, h_buf, k_buf
91      integer ish, jsh, ucont
92      integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom
93      integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom
94C      integer nint_out
95      integer ihi,jhi
96      logical status
97      character*255 mo_basis, geom_name
98c
99#include "bas_exndcf_sfn.fh"
100#include "bas_ibs_sfn.fh"
101c
102      mo_basis = 'ao basis'
103      geom_name = 'geometry'
104c
105      if(.not.geom_create(geom,geom_name))call errquit
106     &    ('kgdtest_rel1e: geom create error',911, GEOM_ERR)
107      if(.not.bas_create(basis,mo_basis))call errquit
108     &    ('kgdtest_rel1e: basis create error',911, BASIS_ERR)
109c
110      if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit
111     &    ('kgdtest_rel1e: geom load ',911, RTDB_ERR)
112      if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit
113     &    ('kgdtest_rel1e: basis load ',911, RTDB_ERR)
114c
115      basis_id = basis + BASIS_HANDLE_OFFSET
116      nshell = ncont_tot_gb(basis_id)
117      if (.not.int_normalize(rtdb,basis)) call errquit
118     &    ('kgdtest_rel1e: error normalizing ',911, INT_ERR)
119c
120      call int_init(rtdb,1,basis)
121      memscr = 100 000
122      membuf = 1000
123      if (.not.ma_push_get(mt_dbl,memscr,' scratch ',
124     &    h_scr, k_scr)) call errquit
125     &    (' ma error 1',911, MA_ERR)
126      if (.not.ma_push_get(mt_dbl,membuf*3,' buf ',
127     &    h_buf, k_buf)) call errquit
128     &    (' ma error 2',911, MA_ERR)
129c
130      do ish = 1,nshell
131        do jsh = 1,ish
132          write(6,*)' ============= shells <',ish,'|',jsh,'>',
133     &        '==================== start =========='
134          write(6,*)' '
135
136          ucont = (sf_ibs_cn2ucn(ish,basis_id))
137          Li      = infbs_cont(CONT_TYPE ,ucont,basis_id)
138          i_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
139          i_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
140          i_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
141          i_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
142          i_cent  = (sf_ibs_cn2ce(ish,basis_id))
143          i_geom  = ibs_geom(basis_id)
144          ihi = i_gen*(Li+1)*(Li+2)/2
145c
146          ucont = (sf_ibs_cn2ucn(jsh,basis_id))
147          Lj      = infbs_cont(CONT_TYPE ,ucont,basis_id)
148          j_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
149          j_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
150          j_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
151          j_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
152          j_cent  = (sf_ibs_cn2ce(jsh,basis_id))
153          j_geom  = ibs_geom(basis_id)
154          jhi = j_gen*(Lj+1)*(Lj+2)/2
155*
156*   Calculate overlap and kinetic energy integrals
157*
158      call int_hf1sp(
159     &        coords(1,i_cent,i_geom),
160     &        dbl_mb(mb_exndcf(i_iexp,basis_id)),
161     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,i_cent,
162     &        coords(1,j_cent,j_geom),
163     &        dbl_mb(mb_exndcf(j_iexp,basis_id)),
164     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,j_cent,
165     &        coords(1,1,i_geom),charge(1,i_geom),ncenter(i_geom),
166     &        dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000),
167     &        .true.,.true.,.false.,.false.,
168     &        .false.,dbl_mb(k_scr),memscr,'kgdtest_rel1e')
169          call ecp_matpr (dbl_mb(k_buf),1,jhi,1,ihi,1,jhi,1,ihi,
170     &        'overlap integrals','E',120,8)
171          call ecp_matpr (dbl_mb(k_buf+1000),1,jhi,1,ihi,1,jhi,1,ihi,
172     &        'kinetic integrals','E',120,8)
173C     modified metric
174          call rel_onel (
175     &        coords(1,i_cent,i_geom),
176     &        dbl_mb(mb_exndcf(i_iexp,basis_id)),
177     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),
178     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,
179     &        coords(1,j_cent,j_geom),
180     &        dbl_mb(mb_exndcf(j_iexp,basis_id)),
181     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),
182     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,
183     &        coords(1,1,i_geom),charge(1,i_geom),
184     &        geom_invnucexp(1,i_geom),ncenter(i_geom),
185     &        dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000),
186     &        membuf,.true.,.false.,.false.,.false.,.true.,.false.,
187     &        .false.,.false.,dbl_mb(k_scr),memscr,0,1)
188          call ecp_matpr (dbl_mb(k_buf),1,jhi,1,ihi,1,jhi,1,ihi,
189     &        'modified overlap integrals','E',120,8)
190C     modified kinetic energy
191          call rel_onel (
192     &        coords(1,i_cent,i_geom),
193     &        dbl_mb(mb_exndcf(i_iexp,basis_id)),
194     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),
195     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,
196     &        coords(1,j_cent,j_geom),
197     &        dbl_mb(mb_exndcf(j_iexp,basis_id)),
198     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),
199     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,
200     &        coords(1,1,i_geom),charge(1,i_geom),
201     &        geom_invnucexp(1,i_geom),ncenter(i_geom),
202     &        dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000),
203     &        membuf,.false.,.true.,.false.,.false.,.true.,.false.,
204     &        .false.,.false.,dbl_mb(k_scr),memscr,0,1)
205          call ecp_matpr (dbl_mb(k_buf+1000),1,jhi,1,ihi,1,jhi,1,ihi,
206     &        'modified kinetic integrals','E',120,8)
207C     modified potential energy
208          call rel_onel (
209     &        coords(1,i_cent,i_geom),
210     &        dbl_mb(mb_exndcf(i_iexp,basis_id)),
211     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),
212     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,
213     &        coords(1,j_cent,j_geom),
214     &        dbl_mb(mb_exndcf(j_iexp,basis_id)),
215     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),
216     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,
217     &        coords(1,1,i_geom),charge(1,i_geom),
218     &        geom_invnucexp(1,i_geom),ncenter(i_geom),
219     &        dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000),
220     &        membuf,.false.,.false.,.true.,.false.,.true.,.false.,
221     &        .false.,.false.,dbl_mb(k_scr),memscr,0,1)
222          call ecp_matpr (dbl_mb(k_buf+2000),1,jhi,1,ihi,1,jhi,1,ihi,
223     &        'modified potential integrals','E',120,8)
224          write(6,*)' ============= shells <',ish,'|',jsh,'>',
225     &          '====================  end  =========='
226          write(6,*)' '
227        enddo
228      enddo
229c
230      call int_terminate()
231      status = ma_pop_stack(h_buf)
232      status = status.and.ma_pop_stack(h_scr)
233      if (.not.status) call errquit('pop failed',911, MA_ERR)
234      status = bas_destroy(basis)
235      status = status.and.geom_destroy(geom)
236      if (.not.status) call errquit('b/g destroy failed',911, GEOM_ERR)
237      return
238      end
239************************************************************************
240      subroutine kgdtest_gencon(rtdb)
241      implicit none
242#include "errquit.fh"
243#include "mafdecls.fh"
244#include "rtdb.fh"
245#include "context.fh"
246#include "geom.fh"
247#include "bas.fh"
248#include "nwc_const.fh"
249#include "basP.fh"
250#include "basdeclsP.fh"
251#include "geomP.fh"
252#include "geobasmapP.fh"
253#include "bas_exndcf_dec.fh"
254#include "bas_ibs_dec.fh"
255#include "stdio.fh"
256c
257      logical int_normalize
258      external int_normalize
259      integer idamax
260      external idamax
261c
262c test general contraction in McMD 2e integrals.
263c
264      integer rtdb
265      integer geom,basis, basis_id
266      integer nshell, memscr, membuf
267      integer h_scr, k_scr, h_buf, k_buf
268      integer ish, jsh, ksh,lsh, ucont
269      integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom
270      integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom
271      integer lk, k_prim, k_gen, k_iexp, k_icfp, k_cent, k_geom
272      integer ll, l_prim, l_gen, l_iexp, l_icfp, l_cent, l_geom
273      integer Nints
274      integer i_eri,i_c,j_c,k_c,l_c,ic,jc,kc,lc,i_seg
275      integer ihi,jhi,khi,lhi,i2,j2,k2,l2,n_cart,n_cont
276      character*255 mo_basis, geom_name
277      double precision difmax
278c
279#include "bas_exndcf_sfn.fh"
280#include "bas_ibs_sfn.fh"
281c
282      mo_basis = 'ao basis'
283      geom_name = 'geometry'
284c
285      if(.not.geom_create(geom,geom_name))call errquit
286     &    ('kgdtest_gencon: geom create error',911, GEOM_ERR)
287      if(.not.bas_create(basis,mo_basis))call errquit
288     &    ('kgdtest_gencon: basis create error',911, BASIS_ERR)
289c
290      if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit
291     &    ('kgdtest_gencon: geom load ',911, RTDB_ERR)
292      if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit
293     &    ('kgdtest_gencon: basis load ',911, RTDB_ERR)
294c
295      basis_id = basis + BASIS_HANDLE_OFFSET
296      nshell = ncont_tot_gb(basis_id)
297      if (.not.int_normalize(rtdb,basis)) call errquit
298     &    ('kgdtest_gencon: error normalizing ',911, INT_ERR)
299c
300      call int_init(rtdb,1,basis)
301      memscr = 1 000 000
302      membuf = 10 000
303      if (.not.ma_push_get(mt_dbl,memscr,' scratch ',
304     &    h_scr, k_scr)) call errquit
305     &    (' ma error 1',911, MA_ERR)
306      if (.not.ma_push_get(mt_dbl,membuf,' buf ',
307     &    h_buf, k_buf)) call errquit
308     &    (' ma error 2',911, MA_ERR)
309c
310      difmax = 0.0d0
311      do ish = 1,nshell
312        ucont = (sf_ibs_cn2ucn(ish,basis_id))
313        Li      = infbs_cont(CONT_TYPE ,ucont,basis_id)
314        i_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
315        i_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
316        i_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
317        i_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
318        i_cent  = (sf_ibs_cn2ce(ish,basis_id))
319        i_geom  = ibs_geom(basis_id)
320        i2 = (Li+1)*(Li+2)/2
321        ihi = i_gen*i2
322c
323        do jsh = 1,ish
324          ucont = (sf_ibs_cn2ucn(jsh,basis_id))
325          Lj      = infbs_cont(CONT_TYPE ,ucont,basis_id)
326          j_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
327          j_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
328          j_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
329          j_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
330          j_cent  = (sf_ibs_cn2ce(jsh,basis_id))
331          j_geom  = ibs_geom(basis_id)
332          j2 = (Lj+1)*(Lj+2)/2
333          jhi = j_gen*j2
334c
335          do ksh = 1,ish
336            ucont = (sf_ibs_cn2ucn(ksh,basis_id))
337            Lk      = infbs_cont(CONT_TYPE ,ucont,basis_id)
338            k_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
339            k_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
340            k_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
341            k_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
342            k_cent  = (sf_ibs_cn2ce(ish,basis_id))
343            k_geom  = ibs_geom(basis_id)
344            k2 = (Lk+1)*(Lk+2)/2
345            khi = k_gen*k2
346c
347            do lsh = 1,ksh
348              ucont = (sf_ibs_cn2ucn(lsh,basis_id))
349              Ll      = infbs_cont(CONT_TYPE ,ucont,basis_id)
350              l_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
351              l_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
352              l_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
353              l_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
354              l_cent  = (sf_ibs_cn2ce(lsh,basis_id))
355              l_geom  = ibs_geom(basis_id)
356              l2 = (Ll+1)*(Ll+2)/2
357              lhi = l_gen*l2
358
359              n_cart = i2*j2*k2*l2
360              n_cont = i_gen*j_gen*k_gen*l_gen
361              Nints = ihi*jhi*khi*lhi
362              write (LuOut,*) ish,Li,i_prim,i_gen,i_cent
363              write (LuOut,*) jsh,Lj,j_prim,j_gen,j_cent
364              write (LuOut,*) ksh,Lk,k_prim,k_gen,k_cent
365              write (LuOut,*) lsh,Ll,l_prim,l_gen,l_cent
366C              if (n_cont .gt. 1 .and. Li.ne.Lj .and. Lk.ne.Ll) then
367C                write(luout,*)' '
368C                write(luout,*)' ============= (',
369C     &              ish,':',Li,',',
370C     &              jsh,':',Lj,'|',
371C     &              ksh,':',Lk,',',
372C     &              lsh,':',Ll,')',
373C     &              '===================='
374C                write(luout,*)' '
375                write (luout,*) 'Nints',Nints
376
377                call hf2(
378     &            coords(1,i_cent,i_geom),
379     &              dbl_mb(mb_exndcf(i_iexp,basis_id)),
380     &              dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,
381     &              coords(1,j_cent,j_geom),
382     &              dbl_mb(mb_exndcf(j_iexp,basis_id)),
383     &              dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,
384     &              coords(1,k_cent,k_geom),
385     &              dbl_mb(mb_exndcf(k_iexp,basis_id)),
386     &              dbl_mb(mb_exndcf(k_icfp,basis_id)),k_prim,k_gen,Lk,
387     &              coords(1,l_cent,l_geom),
388     &              dbl_mb(mb_exndcf(l_iexp,basis_id)),
389     &              dbl_mb(mb_exndcf(l_icfp,basis_id)),l_prim,l_gen,Ll,
390     &              dbl_mb(k_buf),Nints,.false.,.false.,.false.,
391     &              .false.,dbl_mb(k_scr),memscr)
392
393                i_eri = k_buf+Nints
394                i_seg = i_eri
395                i_c = i_icfp
396                do ic = 1,i_gen
397                  j_c = j_icfp
398                  do jc = 1,j_gen
399                    k_c = k_icfp
400                    do kc = 1,k_gen
401                      l_c = l_icfp
402                      do lc = 1,l_gen
403C                        write (luout,*) 'ic,jc,kc,lc',ic,jc,kc,lc
404                        call hf2(
405     &                      coords(1,i_cent,i_geom),
406     &                      dbl_mb(mb_exndcf(i_iexp,basis_id)),
407     &                      dbl_mb(mb_exndcf(i_c,basis_id)),
408     &                      i_prim,1,Li,
409     &                      coords(1,j_cent,j_geom),
410     &                      dbl_mb(mb_exndcf(j_iexp,basis_id)),
411     &                      dbl_mb(mb_exndcf(j_c,basis_id)),
412     &                      j_prim,1,Lj,
413     &                      coords(1,k_cent,k_geom),
414     &                      dbl_mb(mb_exndcf(k_iexp,basis_id)),
415     &                      dbl_mb(mb_exndcf(k_c,basis_id)),
416     &                      k_prim,1,Lk,
417     &                      coords(1,l_cent,l_geom),
418     &                      dbl_mb(mb_exndcf(l_iexp,basis_id)),
419     &                      dbl_mb(mb_exndcf(l_c,basis_id)),
420     &                      l_prim,1,Ll,
421     &                      dbl_mb(i_eri),Nints,.false.,.false.,.false.,
422     &                      .false.,dbl_mb(k_scr),memscr)
423
424                        i_eri = i_eri+i2*j2*k2*l2
425                        l_c = l_c+l_prim
426                      end do
427                      k_c = k_c+k_prim
428                    end do
429                    j_c = j_c+j_prim
430                  end do
431                  i_c = i_c+i_prim
432                end do
433
434                call reorder_2eints(dbl_mb(i_eri),dbl_mb(i_seg),
435     &              l2,l_gen,k2,k_gen,j2,j_gen,i2,i_gen)
436                call dcopy (Nints,dbl_mb(i_eri),1,dbl_mb(i_seg),1)
437                call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1)
438                ic = idamax(Nints,dbl_mb(i_eri),1)-1
439                write (luout,*) 'Maximum difference',dbl_mb(i_eri+ic)
440                difmax = max(difmax,abs(dbl_mb(i_eri+ic)))
441                if (abs(dbl_mb(i_eri+ic)) .gt. 1.0d-12) then
442                  call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi,
443     &                1,khi*lhi,1,ihi*jhi,'General contraction',
444     &                'E',120,6)
445                  call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi,
446     &                1,khi*lhi,1,ihi*jhi,'Segmented contraction',
447     &                'E',120,6)
448                end if
449C              end if
450            end do
451          end do
452        end do
453      end do
454      write (luout,*) 'Maximum difference of all blocks',difmax
455
456      end
457************************************************************************
458      subroutine reorder_2eints(bnew,bold,lc,lg,kc,kg,jc,jg,ic,ig)
459      double precision bnew(lc,lg,kc,kg,jc,jg,ic*ig)
460      double precision bold(lc,kc,jc,ic,lg,kg,jg*ig)
461      integer lc,lg,kc,kg,jc,jg,ic,ig
462      integer i,j,k,l,ii,jj,kk,ll,inew,iold
463      iold = 0
464      inew = 0
465      do i = 1,ig
466        do j = 1,jg
467          iold = iold+1
468          do k = 1,kg
469            do l = 1,lg
470              do ii = 1,ic
471                do jj = 1,jc
472                  do kk = 1,kc
473                    do ll = 1,lc
474                      bnew(ll,l,kk,k,jj,j,inew+ii) =
475     &                    bold(ll,kk,jj,ii,l,k,iold)
476                    end do
477                  end do
478                end do
479              end do
480            end do
481          end do
482        end do
483        inew = inew+ic
484      end do
485      end
486
487************************************************************************
488      subroutine kgdtest_soecp(rtdb)
489      implicit none
490#include "errquit.fh"
491#include "mafdecls.fh"
492#include "rtdb.fh"
493#include "context.fh"
494#include "geom.fh"
495#include "bas.fh"
496#include "nwc_const.fh"
497#include "basP.fh"
498#include "basdeclsP.fh"
499#include "geomP.fh"
500#include "geobasmapP.fh"
501#include "bas_exndcf_dec.fh"
502#include "bas_ibs_dec.fh"
503#include "ecp_nwc.fh"
504#include "stdio.fh"
505c
506      logical int_normalize
507      external int_normalize
508      integer idamax
509      external idamax
510c
511c test spin-orbit integrals
512c
513      integer rtdb
514      integer geom,basis, basis_id
515      integer nshell, memscr, membuf
516      integer h_scr, k_scr, h_buf, k_buf
517      integer ish, jsh, ucont
518      integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom
519      integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom
520      integer ihi,jhi,i2,j2
521      integer ibug,n_blk
522      character*255 mo_basis, geom_name
523c
524#include "bas_exndcf_sfn.fh"
525#include "bas_ibs_sfn.fh"
526c
527      mo_basis = 'ao basis'
528      geom_name = 'geometry'
529c
530      if(.not.geom_create(geom,geom_name))call errquit
531     &    ('kgdtest_gencon: geom create error',911, GEOM_ERR)
532      if(.not.bas_create(basis,mo_basis))call errquit
533     &    ('kgdtest_gencon: basis create error',911, BASIS_ERR)
534c
535      if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit
536     &    ('kgdtest_gencon: geom load ',911, RTDB_ERR)
537      if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit
538     &    ('kgdtest_gencon: basis load ',911, RTDB_ERR)
539c
540      basis_id = basis + BASIS_HANDLE_OFFSET
541      nshell = ncont_tot_gb(basis_id)
542      if (.not.int_normalize(rtdb,basis)) call errquit
543     &    ('kgdtest_gencon: error normalizing ',911, INT_ERR)
544c
545      call int_init(rtdb,1,basis)
546      memscr = 1 000 000
547      membuf = 10 000
548      if (.not.ma_push_get(mt_dbl,memscr,' scratch ',
549     &    h_scr, k_scr)) call errquit
550     &    (' ma error 1',911, MA_ERR)
551      if (.not.ma_push_get(mt_dbl,membuf,' buf ',
552     &    h_buf, k_buf)) call errquit
553     &    (' ma error 2',911, MA_ERR)
554c
555      if (.not. rtdb_get(rtdb,'ecp:ibug',mt_int,1,ibug))
556     &    ibug = 3
557      if (.not. rtdb_get(rtdb,'ecp:n_blk',mt_int,1,n_blk))
558     &    n_blk = 3
559c
560      do ish = 1,nshell
561        ucont = (sf_ibs_cn2ucn(ish,basis_id))
562        Li      = infbs_cont(CONT_TYPE ,ucont,basis_id)
563        i_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
564        i_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
565        i_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
566        i_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
567        i_cent  = (sf_ibs_cn2ce(ish,basis_id))
568        i_geom  = ibs_geom(basis_id)
569        i2 = (Li+1)*(Li+2)/2
570        ihi = i_gen*i2
571c
572        do jsh = 1,ish
573          ucont = (sf_ibs_cn2ucn(jsh,basis_id))
574          Lj      = infbs_cont(CONT_TYPE ,ucont,basis_id)
575          j_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
576          j_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
577          j_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
578          j_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
579          j_cent  = (sf_ibs_cn2ce(jsh,basis_id))
580          j_geom  = ibs_geom(basis_id)
581          j2 = (Lj+1)*(Lj+2)/2
582          jhi = j_gen*j2
583
584          write(luout,*)' '
585          write(luout,*)' ============= ',ish,jsh,
586     &        '===================='
587          write(luout,*)' '
588
589          call ecp_integral (
590     &        coords(1,i_cent,i_geom),
591     &        dbl_mb(mb_exndcf(i_iexp,basis_id)),
592     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,i_cent,
593     &        coords(1,j_cent,j_geom),
594     &        dbl_mb(mb_exndcf(j_iexp,basis_id)),
595     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,j_cent,
596     &        dbl_mb(k_xyzecp),
597     &        dbl_mb(k_ecp_e),dbl_mb(k_ecp_c),
598     &        int_mb(k_ecp_nprim_c),
599     &        int_mb(k_ecp_ncoef_c), ! new name is n_colc_C
600     &        int_mb(k_ecp_ind_z),
601     &        int_mb(k_ecp_ind_c),
602     &        n_zeta_c,
603     &        n_zeta_c,
604     &        int_mb(k_ecp_l_c),
605     &        int_mb(k_ecp_lip),
606     &        n_ecp,l_ecp,
607     &        0,
608     &        dbl_mb(k_ecp_c2s),mem_c2s,
609     &        dbl_mb(k_buf),ihi*jhi,n_blk,
610     &        .false.,dbl_mb(k_scr),memscr,
611     &        ibug)
612
613        end do
614      end do
615      call int_terminate()
616      end
617************************************************************************
618      subroutine kgdtest_rel2e(rtdb)
619      implicit none
620#include "errquit.fh"
621#include "mafdecls.fh"
622#include "rtdb.fh"
623#include "context.fh"
624#include "geom.fh"
625#include "bas.fh"
626#include "nwc_const.fh"
627#include "basP.fh"
628#include "basdeclsP.fh"
629#include "geomP.fh"
630#include "geobasmapP.fh"
631#include "bas_exndcf_dec.fh"
632#include "bas_ibs_dec.fh"
633#include "apiP.fh"
634#include "rel_nwc.fh"
635#include "stdio.fh"
636c
637      logical int_normalize
638      external int_normalize
639      integer idamax
640      external idamax
641c
642c test relativistic 2e integrals for correct magnitude
643c
644      integer rtdb
645      integer geom, basis, basis_id
646      integer nshell, memscr, membuf
647      integer h_scr, k_scr, h_buf, k_buf
648      integer ish, jsh, ksh,lsh, ucont
649      integer li, ipr, igc, iex, icf, icfs, ice, igm
650      integer lj, jpr, jgc, jex, jcf, jcfs, jce, jgm
651      integer lk, kpr, kgc, kex, kcf, kcfs, kce, kgm
652      integer ll, lpr, lgc, lex, lcf, lcfs, lce, lgm
653      integer Nints
654      integer i_eri,i,abas,sbas
655      integer ihi,jhi,khi,lhi,i2,j2,k2,l2,n_cart,n_cont
656      character*255 mo_basis, geom_name
657      double precision difmax,errmax,dif,erimax
658c
659#include "bas_exndcf_sfn.fh"
660#include "bas_ibs_sfn.fh"
661c
662      mo_basis = 'ao basis'
663      geom_name = 'geometry'
664c
665      if(.not.geom_create(geom,geom_name))call errquit
666     &    ('kgdtest_rel2e: geom create error',911, GEOM_ERR)
667      if(.not.bas_create(basis,mo_basis))call errquit
668     &    ('kgdtest_rel2e: basis create error',911, BASIS_ERR)
669c
670      if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit
671     &    ('kgdtest_rel2e: geom load ',911, GEOM_ERR)
672      if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit
673     &    ('kgdtest_rel2e: basis load ',911, BASIS_ERR)
674c
675      basis_id = basis + BASIS_HANDLE_OFFSET
676      nshell = ncont_tot_gb(basis_id)
677      if (.not.int_normalize(rtdb,basis)) call errquit
678     &    ('kgdtest_rel2e: error normalizing ',911, INT_ERR)
679c
680      call int_init(rtdb,1,basis)
681      memscr = 5 000 000
682      membuf = 100 000
683      if (.not.ma_push_get(mt_dbl,memscr,' scratch ',
684     &    h_scr, k_scr)) call errquit
685     &    (' ma error 1',911, MA_ERR)
686      if (.not.ma_push_get(mt_dbl,membuf,' buf ',
687     &    h_buf, k_buf)) call errquit
688     &    (' ma error 2',911, MA_ERR)
689c
690      sbas = sc_bsh + BASIS_HANDLE_OFFSET
691      abas = ao_bsh + BASIS_HANDLE_OFFSET
692c
693      difmax = 0.0d0
694      errmax = 0.0d0
695      rel_dbg = 3
696C      do ish = 1,nshell
697C      do ish = 2,12,10
698      ish = 12
699        ucont = (sf_ibs_cn2ucn(ish,abas))
700        Li   = infbs_cont(CONT_TYPE ,ucont,abas)
701        ipr  = infbs_cont(CONT_NPRIM,ucont,abas)
702        igc  = infbs_cont(CONT_NGEN ,ucont,abas)
703        iex  = infbs_cont(CONT_IEXP ,ucont,abas)
704        icf  = infbs_cont(CONT_ICFP ,ucont,abas)
705        ice  = (sf_ibs_cn2ce(ish,abas))
706        igm  = ibs_geom(abas)
707        i2 = (Li+1)*(Li+2)/2
708        ihi = igc*i2
709        ucont = ao_to_ls(ucont)
710        icfs = infbs_cont(CONT_ICFP ,ucont,sbas)
711c
712C        do jsh = 1,ish
713        jsh = 11
714C        jsh = ish-1
715          ucont = (sf_ibs_cn2ucn(jsh,abas))
716          Lj   = infbs_cont(CONT_TYPE ,ucont,abas)
717          jpr  = infbs_cont(CONT_NPRIM,ucont,abas)
718          jgc  = infbs_cont(CONT_NGEN ,ucont,abas)
719          jex  = infbs_cont(CONT_IEXP ,ucont,abas)
720          jcf  = infbs_cont(CONT_ICFP ,ucont,abas)
721          jce  = (sf_ibs_cn2ce(jsh,abas))
722          jgm  = ibs_geom(abas)
723          j2 = (Lj+1)*(Lj+2)/2
724          jhi = jgc*j2
725          ucont = ao_to_ls(ucont)
726          jcfs = infbs_cont(CONT_ICFP ,ucont,sbas)
727c
728C          do ksh = 1,ish
729          ksh = 1
730            ucont = (sf_ibs_cn2ucn(ksh,abas))
731            Lk   = infbs_cont(CONT_TYPE ,ucont,abas)
732            kpr  = infbs_cont(CONT_NPRIM,ucont,abas)
733            kgc  = infbs_cont(CONT_NGEN ,ucont,abas)
734            kex  = infbs_cont(CONT_IEXP ,ucont,abas)
735            kcf  = infbs_cont(CONT_ICFP ,ucont,abas)
736            kce  = (sf_ibs_cn2ce(ksh,abas))
737            kgm  = ibs_geom(abas)
738            k2 = (Lk+1)*(Lk+2)/2
739            khi = kgc*k2
740            ucont = ao_to_ls(ucont)
741            kcfs = infbs_cont(CONT_ICFP ,ucont,sbas)
742c
743C            do lsh = 1,ksh
744            lsh = 1
745              ucont = (sf_ibs_cn2ucn(lsh,abas))
746              Ll   = infbs_cont(CONT_TYPE ,ucont,abas)
747              lpr  = infbs_cont(CONT_NPRIM,ucont,abas)
748              lgc  = infbs_cont(CONT_NGEN ,ucont,abas)
749              lex  = infbs_cont(CONT_IEXP ,ucont,abas)
750              lcf  = infbs_cont(CONT_ICFP ,ucont,abas)
751              lce  = (sf_ibs_cn2ce(lsh,abas))
752              lgm  = ibs_geom(abas)
753              l2 = (Ll+1)*(Ll+2)/2
754              lhi = lgc*l2
755              ucont = ao_to_ls(ucont)
756              lcfs = infbs_cont(CONT_ICFP ,ucont,sbas)
757
758              n_cart = i2*j2*k2*l2
759              n_cont = igc*jgc*kgc*lgc
760              Nints = ihi*jhi*khi*lhi
761              i_eri = k_buf+Nints
762C              if (n_cont .gt. 1 .and. Li.ne.Lj .and. Lk.ne.Ll) then
763                write(luout,*)' '
764                write(luout,*)' ============= (',
765     &              ish,':',Li,',',
766     &              jsh,':',Lj,'|',
767     &              ksh,':',Lk,',',
768     &              lsh,':',Ll,')',
769     &              '===================='
770                write(luout,*)' '
771C                write (luout,*) 'Nints',Nints
772
773              call hf2(
774     &            coords(1,ice,igm),
775     &            dbl_mb(mb_exndcf(iex,abas)),
776     &            dbl_mb(mb_exndcf(icf,abas)),ipr,igc,Li,
777     &            coords(1,jce,jgm),
778     &            dbl_mb(mb_exndcf(jex,abas)),
779     &            dbl_mb(mb_exndcf(jcf,abas)),jpr,jgc,Lj,
780     &            coords(1,kce,kgm),
781     &            dbl_mb(mb_exndcf(kex,abas)),
782     &            dbl_mb(mb_exndcf(kcf,abas)),kpr,kgc,Lk,
783     &            coords(1,lce,lgm),
784     &            dbl_mb(mb_exndcf(lex,abas)),
785     &            dbl_mb(mb_exndcf(lcf,abas)),lpr,lgc,Ll,
786     &            dbl_mb(k_buf),Nints,.false.,.false.,.false.,
787     &            .false.,dbl_mb(k_scr),memscr)
788                  call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi,
789     &                1,khi*lhi,1,ihi*jhi,'Large component only',
790     &                'E',120,6)
791C              write (luout,*) 'hf2 completed'
792C              call util_flush(6)
793              call rel_2e4c_sf (
794     &            coords(1,ice,igm),dbl_mb(mb_exndcf(iex,abas)),
795     &            dbl_mb(mb_exndcf(icf,abas)),
796     &            dbl_mb(mb_exndcf(icfs,sbas)),ipr,igc,Li,ice,
797     &            coords(1,jce,jgm),dbl_mb(mb_exndcf(jex,abas)),
798     &            dbl_mb(mb_exndcf(jcf,abas)),
799     &            dbl_mb(mb_exndcf(jcfs,sbas)),jpr,jgc,Lj,jce,
800     &            coords(1,kce,kgm),dbl_mb(mb_exndcf(kex,abas)),
801     &            dbl_mb(mb_exndcf(kcf,abas)),
802     &            dbl_mb(mb_exndcf(kcfs,sbas)),kpr,kgc,Lk,kce,
803     &            coords(1,lce,lgm),dbl_mb(mb_exndcf(lex,abas)),
804     &            dbl_mb(mb_exndcf(lcf,abas)),
805     &            dbl_mb(mb_exndcf(lcfs,sbas)),lpr,lgc,Ll,lce,
806c...................... canAB   canCD   canPQ   DryRun
807     &            dbl_mb(i_eri),Nints,.false.,.false.,.false.,.false.,
808     &            dbl_mb(k_scr),memscr,
809     &            .true.,.true.,ss_one_cent,do_ssss,rel_dbg)
810
811                  call ecp_matpr (dbl_mb(i_eri),1,khi*lhi,1,ihi*jhi,
812     &                1,khi*lhi,1,ihi*jhi,'Relativistic integrals',
813     &                'E',120,6)
814C              write (luout,*) 'rel_2e4c_sf completed'
815C              call util_flush(6)
816              call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1)
817              dif = 0.0d0
818              erimax = 0.0d0
819              do i = 0,Nints-1
820                erimax = max(erimax,abs(dbl_mb(k_buf+i)))
821                if (abs(dbl_mb(k_buf+i)) .gt. 1.0d-12) then
822                  dif = max(dif,abs(dbl_mb(i_eri+i)/dbl_mb(k_buf+i)))
823                else if (abs(dbl_mb(i_eri+i)) .gt. 1.0d-12) then
824                  if (abs(dbl_mb(k_buf+i)) .ne. 0.0d0) then
825                    dif = max(dif,abs(dbl_mb(i_eri+i)/dbl_mb(k_buf+i)))
826                  else
827                    errmax = max(errmax,abs(dbl_mb(i_eri+i)))
828                  end if
829                end if
830              end do
831              difmax = max(difmax,dif)
832              if (dif .gt. 0.1d0) then
833                write (LuOut,*) ish,Li,ipr,igc,ice
834                write (LuOut,*) jsh,Lj,jpr,jgc,jce
835                write (LuOut,*) ksh,Lk,kpr,kgc,kce
836                write (LuOut,*) lsh,Ll,lpr,lgc,lce
837                write (LuOut,*) dif,erimax
838                call util_flush(LuOut)
839              end if
840C            end do
841C          end do
842C        end do
843C      end do
844      write (luout,*) 'Maximum difference of all blocks',difmax
845      write (luout,*) 'Maximum difference from zero integrals',errmax
846      call int_terminate()
847
848      end
849************************************************************************
850      subroutine kgdtest_ecpmem(rtdb)
851      implicit none
852#include "errquit.fh"
853#include "mafdecls.fh"
854#include "rtdb.fh"
855#include "context.fh"
856#include "geom.fh"
857#include "bas.fh"
858#include "nwc_const.fh"
859#include "basP.fh"
860#include "basdeclsP.fh"
861#include "geomP.fh"
862#include "geobasmapP.fh"
863#include "bas_exndcf_dec.fh"
864#include "bas_ibs_dec.fh"
865#include "ecp_nwc.fh"
866#include "stdio.fh"
867c
868      logical int_normalize
869      external int_normalize
870      integer idamax
871      external idamax
872c
873c test spin-orbit integrals
874c
875      integer rtdb
876      integer geom,basis, basis_id
877      integer nshell, memscr, membuf, maxscr
878      integer h_scr, k_scr, h_buf, k_buf
879      integer ish, jsh, ucont
880      integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom
881      integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom
882      integer ihi,jhi,i2,j2
883      integer ibug,n_blk
884      character*255 mo_basis, geom_name
885c
886#include "bas_exndcf_sfn.fh"
887#include "bas_ibs_sfn.fh"
888c
889      mo_basis = 'ao basis'
890      geom_name = 'geometry'
891c
892      if(.not.geom_create(geom,geom_name))call errquit
893     &    ('kgdtest_gencon: geom create error',911, GEOM_ERR)
894      if(.not.bas_create(basis,mo_basis))call errquit
895     &    ('kgdtest_gencon: basis create error',911, BASIS_ERR)
896c
897      if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit
898     &    ('kgdtest_gencon: geom load ',911, RTDB_ERR)
899      if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit
900     &    ('kgdtest_gencon: basis load ',911, RTDB_ERR)
901c
902      basis_id = basis + BASIS_HANDLE_OFFSET
903      nshell = ncont_tot_gb(basis_id)
904      if (.not.int_normalize(rtdb,basis)) call errquit
905     &    ('kgdtest_gencon: error normalizing ',911, INT_ERR)
906c
907      call int_init(rtdb,1,basis)
908      memscr = 1 000 000
909      membuf = 10 000
910      if (.not.ma_push_get(mt_dbl,memscr,' scratch ',
911     &    h_scr, k_scr)) call errquit
912     &    (' ma error 1',911, MA_ERR)
913      if (.not.ma_push_get(mt_dbl,membuf,' buf ',
914     &    h_buf, k_buf)) call errquit
915     &    (' ma error 2',911, MA_ERR)
916c
917      if (.not. rtdb_get(rtdb,'ecp:ibug',mt_int,1,ibug))
918     &    ibug = 1
919      if (.not. rtdb_get(rtdb,'ecp:n_blk',mt_int,1,n_blk))
920     &    n_blk = 1
921c
922      maxscr = 0
923      do ish = 1,nshell
924        ucont = (sf_ibs_cn2ucn(ish,basis_id))
925        Li      = infbs_cont(CONT_TYPE ,ucont,basis_id)
926        i_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
927        i_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
928        i_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
929        i_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
930        i_cent  = (sf_ibs_cn2ce(ish,basis_id))
931        i_geom  = ibs_geom(basis_id)
932        i2 = (Li+1)*(Li+2)/2
933        ihi = i_gen*i2
934c
935        do jsh = 1,ish
936          ucont = (sf_ibs_cn2ucn(jsh,basis_id))
937          Lj      = infbs_cont(CONT_TYPE ,ucont,basis_id)
938          j_prim  = infbs_cont(CONT_NPRIM,ucont,basis_id)
939          j_gen   = infbs_cont(CONT_NGEN ,ucont,basis_id)
940          j_iexp  = infbs_cont(CONT_IEXP ,ucont,basis_id)
941          j_icfp  = infbs_cont(CONT_ICFP ,ucont,basis_id)
942          j_cent  = (sf_ibs_cn2ce(jsh,basis_id))
943          j_geom  = ibs_geom(basis_id)
944          j2 = (Lj+1)*(Lj+2)/2
945          jhi = j_gen*j2
946
947          write(luout,*)' '
948          write(luout,*)' ============= ',ish,jsh,
949     &        '===================='
950          write(luout,*)' '
951
952          call ecp_integral (
953     &        coords(1,i_cent,i_geom),
954     &        dbl_mb(mb_exndcf(i_iexp,basis_id)),
955     &        dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,i_cent,
956     &        coords(1,j_cent,j_geom),
957     &        dbl_mb(mb_exndcf(j_iexp,basis_id)),
958     &        dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,j_cent,
959     &        dbl_mb(k_xyzecp),
960     &        dbl_mb(k_ecp_e),dbl_mb(k_ecp_c),
961     &        int_mb(k_ecp_nprim_c),
962     &        int_mb(k_ecp_ncoef_c), ! new name is n_colc_C
963     &        int_mb(k_ecp_ind_z),
964     &        int_mb(k_ecp_ind_c),
965     &        n_zeta_c,
966     &        n_zeta_c,
967     &        int_mb(k_ecp_l_c),
968     &        int_mb(k_ecp_lip),
969     &        n_ecp,l_ecp,
970     &        0,
971     &        dbl_mb(k_ecp_c2s),mem_c2s,
972     &        dbl_mb(k_buf),ihi*jhi,n_blk,
973     &        .true.,dbl_mb(k_scr),memscr,
974     &        ibug)
975          maxscr = max(memscr,maxscr)
976
977        end do
978      end do
979      write (luout,*) 'Maximum scratch needed:',maxscr
980      call int_terminate()
981      end
982************************************************************************
983      subroutine kgdtest_relgen(rtdb)
984      implicit none
985#include "errquit.fh"
986#include "mafdecls.fh"
987#include "rtdb.fh"
988#include "context.fh"
989#include "geom.fh"
990#include "bas.fh"
991#include "nwc_const.fh"
992#include "basP.fh"
993#include "basdeclsP.fh"
994#include "geomP.fh"
995#include "geobasmapP.fh"
996#include "bas_exndcf_dec.fh"
997#include "bas_ibs_dec.fh"
998#include "apiP.fh"
999#include "rel_nwc.fh"
1000#include "stdio.fh"
1001c
1002      logical int_normalize
1003      external int_normalize
1004      integer idamax
1005      external idamax
1006c
1007c test general contraction in relativistic 2e integrals.
1008c
1009      integer rtdb
1010      integer geom, basis, basid
1011      integer nshell, memscr, membuf
1012      integer h_scr, k_scr, h_buf, k_buf
1013      integer ish, jsh, ksh,lsh, ucont
1014      integer li, ipr, igc, iex, icf, icfs, ice, ige
1015      integer lj, jpr, jgc, jex, jcf, jcfs, jce, jge
1016      integer lk, kpr, kgc, kex, kcf, kcfs, kce, kge
1017      integer ll, lpr, lgc, lex, lcf, lcfs, lce, lge
1018      integer Nints
1019      integer i_eri,i_c,j_c,k_c,l_c,ic,jc,kc,lc,i_seg
1020      integer ihi,jhi,khi,lhi,i2,j2,k2,l2,n_cart,n_cont
1021      integer abas,sbas,i_cs,j_cs,k_cs,l_cs
1022      character*255 basis_name, geom_name
1023      double precision difmax
1024c
1025#include "bas_exndcf_sfn.fh"
1026#include "bas_ibs_sfn.fh"
1027c
1028      write (luout,*) 'in kgdtest_relgen'
1029      basis_name = 'ao basis'
1030      geom_name = 'geometry'
1031c
1032      if(.not.geom_create(geom,geom_name))call errquit
1033     &    ('kgdtest_relgen: geom create error',911, GEOM_ERR)
1034      if(.not.bas_create(basis,basis_name))call errquit
1035     &    ('kgdtest_relgen: basis create error',911, BASIS_ERR)
1036c
1037      if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit
1038     &    ('kgdtest_relgen: geom load ',911, GEOM_ERR)
1039      if(.not.bas_rtdb_load(rtdb,geom,basis,basis_name)) call errquit
1040     &    ('kgdtest_relgen: basis load ',911, BASIS_ERR)
1041c
1042      basid = basis + BASIS_HANDLE_OFFSET
1043      nshell = ncont_tot_gb(basid)
1044      if (.not.int_normalize(rtdb,basis)) call errquit
1045     &    ('kgdtest_relgen: error normalizing ',911, INT_ERR)
1046c
1047      call int_init(rtdb,1,basis)
1048      memscr = 5 000 000
1049      membuf = 100 000
1050      if (.not.ma_push_get(mt_dbl,memscr,' scratch ',
1051     &    h_scr, k_scr)) call errquit
1052     &    (' ma error 1',911, MA_ERR)
1053      if (.not.ma_push_get(mt_dbl,membuf,' buf ',
1054     &    h_buf, k_buf)) call errquit
1055     &    (' ma error 2',911, MA_ERR)
1056c
1057      sbas = sc_bsh + BASIS_HANDLE_OFFSET
1058      abas = ao_bsh + BASIS_HANDLE_OFFSET
1059      write (luout,*) ' starting shell loops ...'
1060c
1061      difmax = 0.0d0
1062      do ish = 1,nshell
1063        ucont = (sf_ibs_cn2ucn(ish,abas))
1064        Li      = infbs_cont(CONT_TYPE ,ucont,abas)
1065        ipr  = infbs_cont(CONT_NPRIM,ucont,abas)
1066        igc   = infbs_cont(CONT_NGEN ,ucont,abas)
1067        iex  = infbs_cont(CONT_IEXP ,ucont,abas)
1068        icf  = infbs_cont(CONT_ICFP ,ucont,abas)
1069        ice  = (sf_ibs_cn2ce(ish,abas))
1070        ige  = ibs_geom(abas)
1071        i2 = (Li+1)*(Li+2)/2
1072        ihi = igc*i2
1073        ucont = ao_to_ls(ucont)
1074        icfs = infbs_cont(CONT_ICFP ,ucont,sbas)
1075c
1076        do jsh = 1,ish
1077          ucont = (sf_ibs_cn2ucn(jsh,abas))
1078          Lj      = infbs_cont(CONT_TYPE ,ucont,abas)
1079          jpr  = infbs_cont(CONT_NPRIM,ucont,abas)
1080          jgc   = infbs_cont(CONT_NGEN ,ucont,abas)
1081          jex  = infbs_cont(CONT_IEXP ,ucont,abas)
1082          jcf  = infbs_cont(CONT_ICFP ,ucont,abas)
1083          jce  = (sf_ibs_cn2ce(jsh,abas))
1084          jge  = ibs_geom(abas)
1085          j2 = (Lj+1)*(Lj+2)/2
1086          jhi = jgc*j2
1087          ucont = ao_to_ls(ucont)
1088          jcfs = infbs_cont(CONT_ICFP ,ucont,sbas)
1089c
1090          do ksh = 1,ish
1091            ucont = (sf_ibs_cn2ucn(ksh,abas))
1092            Lk      = infbs_cont(CONT_TYPE ,ucont,abas)
1093            kpr  = infbs_cont(CONT_NPRIM,ucont,abas)
1094            kgc   = infbs_cont(CONT_NGEN ,ucont,abas)
1095            kex  = infbs_cont(CONT_IEXP ,ucont,abas)
1096            kcf  = infbs_cont(CONT_ICFP ,ucont,abas)
1097            kce  = (sf_ibs_cn2ce(ish,abas))
1098            kge  = ibs_geom(abas)
1099            k2 = (Lk+1)*(Lk+2)/2
1100            khi = kgc*k2
1101            ucont = ao_to_ls(ucont)
1102            kcfs = infbs_cont(CONT_ICFP ,ucont,sbas)
1103c
1104            do lsh = 1,ksh
1105              ucont = (sf_ibs_cn2ucn(lsh,abas))
1106              Ll      = infbs_cont(CONT_TYPE ,ucont,abas)
1107              lpr  = infbs_cont(CONT_NPRIM,ucont,abas)
1108              lgc   = infbs_cont(CONT_NGEN ,ucont,abas)
1109              lex  = infbs_cont(CONT_IEXP ,ucont,abas)
1110              lcf  = infbs_cont(CONT_ICFP ,ucont,abas)
1111              lce  = (sf_ibs_cn2ce(lsh,abas))
1112              lge  = ibs_geom(abas)
1113              l2 = (Ll+1)*(Ll+2)/2
1114              lhi = lgc*l2
1115              ucont = ao_to_ls(ucont)
1116              lcfs = infbs_cont(CONT_ICFP ,ucont,sbas)
1117
1118              n_cart = i2*j2*k2*l2
1119              n_cont = igc*jgc*kgc*lgc
1120              Nints = ihi*jhi*khi*lhi
1121C              write (LuOut,*) ish,Li,ipr,igc,ice,ihi
1122C              write (LuOut,*) jsh,Lj,jpr,jgc,jce,jhi
1123C              write (LuOut,*) ksh,Lk,kpr,kgc,kce,khi
1124C              write (LuOut,*) lsh,Ll,lpr,lgc,lce,lhi
1125C              if (n_cont .gt. 1 .and. Li.ne.Lj .and. Lk.ne.Ll) then
1126C                write(luout,*)' '
1127C                write(luout,*)' ============= (',
1128C     &              ish,':',Li,',',
1129C     &              jsh,':',Lj,'|',
1130C     &              ksh,':',Lk,',',
1131C     &              lsh,':',Ll,')',
1132C     &              '===================='
1133C                write(luout,*)' '
1134C                write (luout,*) 'Nints',Nints
1135C                write (luout,*) 'calling rel_2e4c_sf gen'
1136C                call util_flush(6)
1137
1138                call rel_2e4c_sf (
1139     &              coords(1,ice,ige),dbl_mb(mb_exndcf(iex,abas)),
1140     &              dbl_mb(mb_exndcf(icf,abas)),
1141     &              dbl_mb(mb_exndcf(icfs,sbas)),ipr,igc,Li,ice,
1142     &              coords(1,jce,jge),dbl_mb(mb_exndcf(jex,abas)),
1143     &              dbl_mb(mb_exndcf(jcf,abas)),
1144     &              dbl_mb(mb_exndcf(jcfs,sbas)),jpr,jgc,Lj,jce,
1145     &              coords(1,kce,kge),dbl_mb(mb_exndcf(kex,abas)),
1146     &              dbl_mb(mb_exndcf(kcf,abas)),
1147     &              dbl_mb(mb_exndcf(kcfs,sbas)),kpr,kgc,Lk,kce,
1148     &              coords(1,lce,lge),dbl_mb(mb_exndcf(lex,abas)),
1149     &              dbl_mb(mb_exndcf(lcf,abas)),
1150     &              dbl_mb(mb_exndcf(lcfs,sbas)),lpr,lgc,Ll,lce,
1151     &              dbl_mb(k_buf),Nints,.false.,.false.,.false.,
1152     &              .false.,dbl_mb(k_scr),memscr,
1153     &              .true.,.true.,.false.,.true.,0)
1154C                  call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi,
1155C     &                1,khi*lhi,1,ihi*jhi,'General contraction',
1156C     &                'E',120,6)
1157C                write (luout,*) 'calling rel_2e4c_sf seg'
1158C                call util_flush(6)
1159
1160                i_eri = k_buf+Nints
1161                i_seg = i_eri
1162                i_c = icf
1163                i_cs = icfs
1164                do ic = 1,igc
1165                  j_c = jcf
1166                  j_cs = jcfs
1167                  do jc = 1,jgc
1168                    k_c = kcf
1169                    k_cs = kcfs
1170                    do kc = 1,kgc
1171                      l_c = lcf
1172                      l_cs = lcfs
1173                      do lc = 1,lgc
1174C                        write (luout,*) 'ic,jc,kc,lc',ic,jc,kc,lc
1175
1176                call rel_2e4c_sf (
1177     &              coords(1,ice,ige),dbl_mb(mb_exndcf(iex,abas)),
1178     &              dbl_mb(mb_exndcf(i_c,abas)),
1179     &              dbl_mb(mb_exndcf(i_cs,sbas)),ipr,1,Li,ice,
1180     &              coords(1,jce,jge),dbl_mb(mb_exndcf(jex,abas)),
1181     &              dbl_mb(mb_exndcf(j_c,abas)),
1182     &              dbl_mb(mb_exndcf(j_cs,sbas)),jpr,1,Lj,jce,
1183     &              coords(1,kce,kge),dbl_mb(mb_exndcf(kex,abas)),
1184     &              dbl_mb(mb_exndcf(k_c,abas)),
1185     &              dbl_mb(mb_exndcf(k_cs,sbas)),kpr,1,Lk,kce,
1186     &              coords(1,lce,lge),dbl_mb(mb_exndcf(lex,abas)),
1187     &              dbl_mb(mb_exndcf(l_c,abas)),
1188     &              dbl_mb(mb_exndcf(l_cs,sbas)),lpr,1,Ll,lce,
1189     &              dbl_mb(i_eri),Nints,.false.,.false.,.false.,
1190     &              .false.,dbl_mb(k_scr),memscr,
1191     &              .true.,.true.,.false.,.true.,0)
1192C                  call ecp_matpr (dbl_mb(k_buf),1,k2*l2,1,i2*j2,
1193C     &                1,k2*l2,1,i2*j2,'Segmented contraction',
1194C     &                'E',120,6)
1195
1196                        i_eri = i_eri+i2*j2*k2*l2
1197                        l_c = l_c+lpr
1198                        l_cs = l_cs+lpr
1199                      end do
1200                      k_c = k_c+kpr
1201                      k_cs = k_cs+kpr
1202                    end do
1203                    j_c = j_c+jpr
1204                    j_cs = j_cs+jpr
1205                  end do
1206                  i_c = i_c+ipr
1207                  i_cs = i_cs+ipr
1208                end do
1209
1210                call reorder_2eints(dbl_mb(i_eri),dbl_mb(i_seg),
1211     &              l2,lgc,k2,kgc,j2,jgc,i2,igc)
1212                call dcopy (Nints,dbl_mb(i_eri),1,dbl_mb(i_seg),1)
1213C                call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi,
1214C     &              1,khi*lhi,1,ihi*jhi,'Segmented contraction',
1215C     &              'E',120,6)
1216                call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1)
1217C                call ecp_matpr (dbl_mb(i_eri),1,khi*lhi,1,ihi*jhi,
1218C     &              1,khi*lhi,1,ihi*jhi,'Differences',
1219C     &              'E',120,6)
1220                ic = idamax(Nints,dbl_mb(i_eri),1)-1
1221C                write (luout,*) 'Maximum difference',dbl_mb(i_eri+ic)
1222                difmax = max(difmax,abs(dbl_mb(i_eri+ic)))
1223                if (abs(dbl_mb(i_eri+ic)) .gt. 1.0d-12) then
1224                  write (luout,*) ish,jsh,ksh,lsh
1225                  call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi,
1226     &                1,khi*lhi,1,ihi*jhi,'General contraction',
1227     &                'E',120,6)
1228                  call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi,
1229     &                1,khi*lhi,1,ihi*jhi,'Segmented contraction',
1230     &                'E',120,6)
1231                end if
1232C              end if
1233            end do
1234          end do
1235        end do
1236      end do
1237      write (luout,*) 'Maximum difference of all blocks',difmax
1238
1239      end
1240************************************************************************
1241      subroutine kgdtest_gen1d(rtdb)
1242      implicit none
1243#include "errquit.fh"
1244#include "mafdecls.fh"
1245#include "rtdb.fh"
1246#include "context.fh"
1247#include "geom.fh"
1248#include "bas.fh"
1249#include "nwc_const.fh"
1250#include "basP.fh"
1251#include "basdeclsP.fh"
1252#include "geomP.fh"
1253#include "geobasmapP.fh"
1254#include "bas_exndcf_dec.fh"
1255#include "bas_ibs_dec.fh"
1256#include "apiP.fh"
1257#include "rel_nwc.fh"
1258#include "stdio.fh"
1259c
1260      logical int_normalize
1261      external int_normalize
1262      integer idamax
1263      external idamax
1264c
1265c test general contraction in derivative 1e integrals
1266c
1267      integer rtdb
1268      integer geom, basis, basid
1269      integer nshell, memscr, membuf
1270      integer h_scr, k_scr, h_buf, k_buf
1271      integer ish, jsh, ucont
1272      integer li, ipr, igc, iex, icf, ice, ige
1273      integer lj, jpr, jgc, jex, jcf, jce, jge
1274      integer Nints
1275      integer i_o2i,i_kei,i_nai
1276      integer ihi,jhi,i2,j2,n_cart,n_cont
1277      character*255 basis_name, geom_name
1278      double precision difmax
1279c
1280#include "bas_exndcf_sfn.fh"
1281#include "bas_ibs_sfn.fh"
1282c
1283      write (luout,*) 'in kgdtest_relgen'
1284      basis_name = 'ao basis'
1285      geom_name = 'geometry'
1286c
1287      if(.not.geom_create(geom,geom_name))call errquit
1288     &    ('kgdtest_relgen: geom create error',911, GEOM_ERR)
1289      if(.not.bas_create(basis,basis_name))call errquit
1290     &    ('kgdtest_relgen: basis create error',911, BASIS_ERR)
1291c
1292      if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit
1293     &    ('kgdtest_relgen: geom load ',911, RTDB_ERR)
1294      if(.not.bas_rtdb_load(rtdb,geom,basis,basis_name)) call errquit
1295     &    ('kgdtest_relgen: basis load ',911, RTDB_ERR)
1296c
1297      basid = basis + BASIS_HANDLE_OFFSET
1298      nshell = ncont_tot_gb(basid)
1299      if (.not.int_normalize(rtdb,basis)) call errquit
1300     &    ('kgdtest_relgen: error normalizing ',911, INT_ERR)
1301c
1302      call int_init(rtdb,1,basis)
1303      memscr = 5 000 000
1304      membuf = 100 000
1305      if (.not.ma_push_get(mt_dbl,memscr,' scratch ',
1306     &    h_scr, k_scr)) call errquit
1307     &    (' ma error 1',911, MA_ERR)
1308      if (.not.ma_push_get(mt_dbl,membuf,' buf ',
1309     &    h_buf, k_buf)) call errquit
1310     &    (' ma error 2',911, MA_ERR)
1311c
1312      difmax = 0.0d0
1313      do ish = 1,nshell
1314        ucont = (sf_ibs_cn2ucn(ish,basid))
1315        Li      = infbs_cont(CONT_TYPE ,ucont,basid)
1316        ipr  = infbs_cont(CONT_NPRIM,ucont,basid)
1317        igc   = infbs_cont(CONT_NGEN ,ucont,basid)
1318        iex  = infbs_cont(CONT_IEXP ,ucont,basid)
1319        icf  = infbs_cont(CONT_ICFP ,ucont,basid)
1320        ice  = (sf_ibs_cn2ce(ish,basid))
1321        ige  = ibs_geom(basid)
1322        i2 = (Li+1)*(Li+2)/2
1323        ihi = igc*i2
1324c
1325        do jsh = 1,ish
1326          ucont = (sf_ibs_cn2ucn(jsh,basid))
1327          Lj      = infbs_cont(CONT_TYPE ,ucont,basid)
1328          jpr  = infbs_cont(CONT_NPRIM,ucont,basid)
1329          jgc   = infbs_cont(CONT_NGEN ,ucont,basid)
1330          jex  = infbs_cont(CONT_IEXP ,ucont,basid)
1331          jcf  = infbs_cont(CONT_ICFP ,ucont,basid)
1332          jce  = (sf_ibs_cn2ce(jsh,basid))
1333          jge  = ibs_geom(basid)
1334          j2 = (Lj+1)*(Lj+2)/2
1335          jhi = jgc*j2
1336
1337          n_cart = i2*j2
1338          n_cont = igc*jgc
1339          Nints = ihi*jhi
1340          write (LuOut,*) ish,Li,ipr,igc,ice,ihi
1341          write (LuOut,*) jsh,Lj,jpr,jgc,jce,jhi
1342          write (luout,*) 'Nints',Nints
1343          i_o2i = k_scr
1344          i_kei = i_o2i+Nints*6
1345          i_nai = i_kei+Nints*6
1346          call util_flush(6)
1347c
1348          call hf1d(
1349     &        coords(1,ice,ige),dbl_mb(mb_exndcf(iex,basid)),
1350     &        dbl_mb(mb_exndcf(icf,basid)),ipr,igc,Li,ice,
1351     &        coords(1,jce,jge),dbl_mb(mb_exndcf(jex,basid)),
1352     &        dbl_mb(mb_exndcf(jcf,basid)),jpr,jgc,Lj,jce,
1353     &        coords(1,1,ige),charge(1,ige),
1354     &        geom_invnucexp(1,ige),ncenter(ige),
1355     &        dbl_mb(i_o2i),dbl_mb(i_kei),dbl_mb(i_nai),Nints,
1356     &        .true.,.true.,.true.,.false.,.false.,
1357     &        dbl_mb(k_scr),memscr)
1358c
1359C          call hf1d_new(
1360C     &        coords(1,ice,ige),dbl_mb(mb_exndcf(iex,basid)),
1361C     &        dbl_mb(mb_exndcf(icf,basid)),ipr,igc,Li,ice,
1362C     &        coords(1,jce,jge),dbl_mb(mb_exndcf(jex,basid)),
1363C     &        dbl_mb(mb_exndcf(jcf,basid)),jpr,jgc,Lj,jce,
1364C     &        coords(1,1,ige),charge(1,ige),
1365C     &        geom_invnucexp(1,ige),ncenter(ige),
1366C     &        dbl_mb(i_o2i),dbl_mb(i_kei),dbl_mb(i_nai),Nints,
1367C     &        .true.,.true.,.true.,.false.,.false.,
1368C     &        dbl_mb(k_scr),memscr)
1369C
1370C                write (6,*) k_buf,i_seg,i_eri
1371C                call reorder_2eints(dbl_mb(i_eri),dbl_mb(i_seg),
1372C     &              l2,lgc,k2,kgc,j2,jgc,i2,igc)
1373C                call dcopy (Nints,dbl_mb(i_eri),1,dbl_mb(i_seg),1)
1374C                call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi,
1375C     &              1,khi*lhi,1,ihi*jhi,'Segmented contraction',
1376C     &              'E',120,6)
1377C                call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1)
1378C                call ecp_matpr (dbl_mb(i_eri),1,khi*lhi,1,ihi*jhi,
1379C     &              1,khi*lhi,1,ihi*jhi,'Differences',
1380C     &              'E',120,6)
1381C                ic = idamax(Nints,dbl_mb(i_eri),1)-1
1382C                write (luout,*) 'Maximum difference',dbl_mb(i_eri+ic)
1383C                difmax = max(difmax,abs(dbl_mb(i_eri+ic)))
1384C                if (abs(dbl_mb(i_eri+ic)) .gt. 1.0d-12) then
1385C                  call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi,
1386C     &                1,khi*lhi,1,ihi*jhi,'General contraction',
1387C     &                'E',120,6)
1388C                  call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi,
1389C     &                1,khi*lhi,1,ihi*jhi,'Segmented contraction',
1390C     &                'E',120,6)
1391C                end if
1392C            end do
1393C          end do
1394        end do
1395      end do
1396      write (luout,*) 'Maximum difference of all blocks',difmax
1397
1398      end
1399