1*
2* $Id$
3*
4      logical function raktask_intdd_3c(rtdb)
5      implicit none
6#include "errquit.fh"
7#include "stdio.fh"
8#include "bas.fh"
9#include "geom.fh"
10#include "mafdecls.fh"
11#include "global.fh"
12*
13*::functions
14      logical int_normalize
15      external int_normalize
16*::passed
17      integer rtdb
18*::local
19      integer basis, geom, nbf, cn_nbf_max, nshell
20      integer nat, nat3
21      integer size, sizesq, sizeg, scr_size
22      integer maxg1, maxg2, maxs1, maxs2
23      logical status
24      integer hbuf, kbuf, hscr, kscr
25      integer hfd, kfd, hfdsq, kfdsq, hxyz, kxyz
26      integer hgradp, kgradp
27      integer hgradm, kgradm
28      integer hbufsum, kbufsum
29*
30      raktask_intdd_3c = .false.
31*
32      if (.not.geom_create(geom,'geometry')) call errquit
33     &    ('geom create failed',911, GEOM_ERR)
34      if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit
35     &    ('geom_rtdb_load failed',911, RTDB_ERR)
36c
37      if (.not.bas_create(basis,'ao basis')) call errquit
38     &    ('bas_create failed',911, BASIS_ERR)
39      if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) call errquit
40     &    ('bas_rtdb_load failed',911, RTDB_ERR)
41c
42      write(luout,*)' geom/basis loaded'
43      call util_flush(luout)
44c
45      if (.not.int_normalize(rtdb,basis)) stop ' norm error 1'
46c
47      if (.not. bas_print(basis))
48     $    call errquit(' basis print failed', 0, BASIS_ERR)
49      if (.not. gbs_map_print(basis))
50     $    call errquit(' gbs_map_print failed', 0, UNKNOWN_ERR)
51c
52      if (.not.bas_numbf(basis,nbf)) call errquit
53     &    ('numbf failed',911, BASIS_ERR)
54c
55      if (.not.bas_numcont(basis,nshell)) call errquit
56     &    ('numbf failed',911, BASIS_ERR)
57c
58      if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe'
59      write(luout,*) 'number of atoms ', nat
60      nat3 = 3*nat
61c
62      if (.not.bas_nbf_cn_max(basis,cn_nbf_max))
63     &    stop 'bas_nbf_cn_max'
64c
65      size   = (cn_nbf_max**4)*78
66      sizesq = (cn_nbf_max**4)*12*12
67      sizeg  = (cn_nbf_max**4)*12
68      size   = size   + size  /10
69      sizesq = sizesq + sizesq/10
70      sizeg  = sizeg  + sizeg /10
71c
72      call intdd_init(rtdb,1,basis)
73      call int_mem_print()
74      call int_mem_1e(maxg1,maxs1)
75      call int_mem_2e4c(maxg2,maxs2)
76*
77      write(luout,*)' maxg1 :',maxg1
78      write(luout,*)' maxs1 :',maxs1
79      write(luout,*)' maxg2 :',maxg2
80      write(luout,*)' size  :',size
81      write(luout,*)' sizesq:',sizesq
82      write(luout,*)' sizeg :',sizeg
83      write(luout,*)' maxs2 :',maxs2
84      call util_flush(luout)
85*
86      scr_size = 2*maxs2
87      status =
88     &    ma_alloc_get(mt_dbl,size,'intdd buffer',hbuf,kbuf)
89      status = status.and.
90     &    ma_alloc_get(mt_dbl,sizesq,'intdd buffer summed',
91     &    hbufsum,kbufsum)
92      status = status.and.
93     &    ma_alloc_get(mt_dbl,scr_size,'scr buffer',hscr,kscr)
94      status = status.and.
95     &    ma_alloc_get(mt_dbl,size,'intdd fd buffer',hfd,kfd)
96      status = status.and.
97     &    ma_alloc_get(mt_dbl,sizesq,'intdd fd sq buffer',hfdsq,kfdsq)
98      status = status.and.
99     &    ma_alloc_get(mt_dbl,3*nat,'coords',hxyz,kxyz)
100      status = status.and.
101     &    ma_alloc_get(mt_dbl,sizeg,'grad +',hgradp,kgradp)
102      status = status.and.
103     &    ma_alloc_get(mt_dbl,sizeg,'grad -',hgradm,kgradm)
104      if (.not.status) stop ' memory alloc failed rak23 (1)'
105*
106      call raktask_intdd_3ca(geom,basis,nbf,nshell,cn_nbf_max,
107     &    size,dbl_mb(kbuf),
108     &    scr_size,dbl_mb(kscr),
109     &    dbl_mb(kgradp),
110     &    dbl_mb(kgradm),
111     &    dbl_mb(kfd),
112     &    dbl_mb(kfdsq),
113     &    nat,dbl_mb(kxyz),
114     &    size,sizesq,sizeg,
115     &    dbl_mb(kbufsum))
116*
117      call intdd_terminate()
118      raktask_intdd_3c = bas_destroy(basis)
119      raktask_intdd_3c = raktask_intdd_3c.and.
120     &    geom_destroy(geom)
121      raktask_intdd_3c = raktask_intdd_3c.and.
122     &    ma_free_heap(hscr)
123      raktask_intdd_3c = raktask_intdd_3c.and.
124     &    ma_free_heap(hbuf)
125      raktask_intdd_3c = raktask_intdd_3c.and.
126     &    ma_free_heap(hfd)
127      raktask_intdd_3c = raktask_intdd_3c.and.
128     &    ma_free_heap(hfdsq)
129      raktask_intdd_3c = raktask_intdd_3c.and.
130     &    ma_free_heap(hxyz)
131      raktask_intdd_3c = raktask_intdd_3c.and.
132     &    ma_free_heap(hgradp)
133      raktask_intdd_3c = raktask_intdd_3c.and.
134     &    ma_free_heap(hgradm)
135      raktask_intdd_3c = raktask_intdd_3c.and.
136     &    ma_free_heap(hbufsum)
137      end
138      subroutine raktask_intdd_3ca(geom,basis,nbf,nshell,cn_nbf_max,
139     &    lbuf,buf,lscr,scr,gradp, gradm,buffd,buffdsq,nat,xyz,
140     &    lbuffd, lbuffdsq, lgrad, bufsum)
141      implicit none
142#include "mafdecls.fh"
143#include "stdio.fh"
144#include "bas.fh"
145#include "nwc_const.fh"
146#include "geomP.fh"
147#include "inp.fh"
148c::functions
149      logical rakdd_checktrans
150      external rakdd_checktrans
151c::passed
152      integer geom
153      integer basis
154      integer nbf
155      integer nshell
156      integer cn_nbf_max
157      integer lbuf
158      integer lscr
159      integer nat
160      integer lbuffd, lbuffdsq, lgrad
161      double precision buffd(lbuffd), buffdsq(lbuffdsq)
162      double precision gradp(lgrad), gradm(lgrad)
163      double precision buf(lbuf)
164      double precision bufsum(lbuffdsq)
165      double precision scr(lscr)
166      double precision xyz(3,nat)
167*
168      integer hbufcp, kbufcp
169      integer nzero,ncount,count
170      integer iatom, jatom, katom, latom
171      integer ish, jsh, ksh
172      integer ilo, ihi, inbf
173      integer jlo, jhi, jnbf
174      integer klo, khi, knbf
175      integer nint, ninth, nintg
176      integer idatom(4)
177      integer idatoms(4)
178      integer idatomp(4)
179      integer idatomm(4)
180      integer atoms2move(4)
181      integer num_atoms2move, atom1, atom2
182      integer nat3
183      integer ixyz, zatom
184      double precision thresh, delta, scale, normmax, normmin, norm
185      character*99 string,strings, stringe
186*
187      normmax = 0.0d00
188      normmin = 1.2345678d05
189      call int_acc_high()
190* store original coordintates
191      nat3 = 3*nat
192      call dcopy(nat3,coords(1,1,geom),1,xyz,1)
193*
194      thresh = 1.0d-07
195      delta  = 0.001d00
196*
197      write(luout,*)'  ',nshell,' total shells '
198      call util_flush(luout)
199*      do ish = 1, nshell
200      do ish = 1, 1
201        if (.not.bas_cn2bfr(basis,ish,ilo,ihi))
202     &      stop 'cn2bfr error i'
203        inbf = ihi - ilo + 1
204        strings = ' '
205        call util_date(strings)
206        if (.not.bas_cn2ce(basis,ish,iatom))
207     &      stop 'bas_cn2ce error i'
208*        do jsh = 1, nshell
209        do jsh = 1, 1
210          if (.not.bas_cn2bfr(basis,jsh,jlo,jhi))
211     &          stop 'cn2bfr error j'
212          jnbf = jhi - jlo + 1
213          if (.not.bas_cn2ce(basis,jsh,jatom))
214     &          stop 'bas_cn2ce error j'
215*          do ksh =  1, nshell
216          do ksh =  1, 6
217            if (.not.bas_cn2bfr(basis,ksh,klo,khi))
218     &            stop 'cn2bfr error k'
219            knbf = khi - klo + 1
220            if (.not.bas_cn2ce(basis,ksh,katom))
221     &            stop 'bas_cn2ce error k'
222            if (.not.(iatom.eq.jatom.and.
223     &            jatom.eq.katom)) then
224              nint  = inbf*jnbf*knbf
225              nintg = nint*12
226              ninth = nint*78
227              call dfill (lbuf,0.0d00,buf,1)
228              call dcopy(nat3,xyz,1,coords(1,1,geom),1)
229              call intdd_2e3c(basis,ish,basis,jsh,ksh,
230     &              lscr,scr,lbuf,buf,idatom)
231              call rakdd_checkt_78(nint,buf)
232              write(luout,*)idatom
233              idatoms(1) = idatom(1)
234              idatoms(2) = idatom(2)
235              idatoms(3) = idatom(3)
236              idatoms(4) = idatom(4)
237              write(6,*)' idatom  1:', idatom
238              write(6,*)' idatoms 1:', idatoms
239              if (.not.ma_alloc_get(mt_dbl,lbuf,
240     &              'copy of buf',hbufcp,kbufcp)) stop ' ma alloc fail'
241              call dcopy(lbuf,buf,1,dbl_mb(kbufcp),1)
242              call rakdd_fill(12,nint,dbl_mb(kbufcp),bufsum,idatoms)
243              write(6,*)' idatom  2:', idatom
244              write(6,*)' idatoms 2:', idatoms
245              if (.not.ma_free_heap(hbufcp)) stop ' ma_free fail'
246              call rakdd_print_dd(nint,buf,idatom)
247              if (.not.rakdd_checktrans(nint,bufsum)) then
248                call rakdd_print_dd(nint,buf,idatom)
249                call rakdd_printsum(nint,bufsum,idatoms)
250              endif
251              nzero = 0
252              ncount = 0
253              do count = 1,ninth
254                if (abs(buf(count)).gt.thresh) then
255                  ncount = ncount + 1
256                else
257                  nzero = nzero + 1
258                endif
259              enddo
260              write(luout,*)nzero,'+',ncount,' != ',ninth
261              if ((nzero+ncount).ne.ninth)
262     &              write(luout,*)nzero,'+',ncount,' != ',ninth
263              atoms2move(1) = iatom
264              atoms2move(2) = iatom
265              atoms2move(3) = jatom
266              atoms2move(4) = katom
267              num_atoms2move = 0
268              do atom1 = 1,4
269                do atom2 = 2,4
270                  if (atom1.ne.atom2) then
271                    if (atoms2move(atom1).eq.atoms2move(atom2))
272     &                    atoms2move(atom2) = 0
273                  endif
274                enddo
275              enddo
276              num_atoms2move = 0
277              do atom1 = 1,4
278                if (atoms2move(atom1).gt.0)
279     &                num_atoms2move = num_atoms2move + 1
280              enddo
281*
282              call dfill (lbuffdsq,0.0d00,buffdsq,1)
283              write(6,*)' atoms2move ',atoms2move
284              do zatom = 1,4
285                if (atoms2move(zatom).gt.0) then
286                  do ixyz = 1,3
287                    nintg = nint*12
288                    call dcopy(nat3,xyz,1,coords(1,1,geom),1)
289                    coords(ixyz,atoms2move(zatom),geom) =
290     &                    coords(ixyz,atoms2move(zatom),geom) + delta
291                    call dfill(lgrad,0.0d00,gradp,1)
292                    call intd_2e3c(basis,ish,basis,jsh,ksh,
293     &                    lscr,scr,lgrad,gradp,idatomp)
294                    write(string,*)' grad +',ixyz,zatom
295*                      call rakdd_printgrad(nint,gradp,string)
296                    write(6,*)'idatomp:',idatomp
297                    call dcopy(nat3,xyz,1,coords(1,1,geom),1)
298                    coords(ixyz,atoms2move(zatom),geom) =
299     &                    coords(ixyz,atoms2move(zatom),geom) - delta
300                    call dfill(lgrad,0.0d00,gradm,1)
301                    call intd_2e3c(basis,ish,basis,jsh,ksh,
302     &                    lscr,scr,lgrad,gradm,idatomm)
303                    write(string,*)' grad -',ixyz,zatom
304*                      call rakdd_printgrad(nint,gradm,string)
305                    call dcopy(nat3,xyz,1,coords(1,1,geom),1)
306                    write(6,*)'idatomm:',idatomm
307                    call daxpy(nintg,-1.0d00,gradm,1,gradp,1)
308                    write(string,*)' grad diff',ixyz,zatom
309*                      call rakdd_printgrad(nint,gradp,string)
310                    scale = 1.0d00/(2.0d00*delta)
311                    call dscal(nintg,scale,gradp,1)
312                    write(string,*)' grad diff scaled',ixyz,zatom
313*                      call rakdd_printgrad(nint,gradp,string)
314                    call rakdd_grad_fill(gradp,idatomp,
315     &                    idatomm,ixyz,
316     &                    zatom,idatoms,
317     &                    atoms2move(zatom),buffdsq,nint)
318                  enddo
319                endif
320              enddo             ! zatom
321              call rakdd_fill_b(buffdsq,nint)  ! check if symmetric
322              call rakdd_printsq(buffdsq,nint)
323*                call rakdd_print_both(80,buffdsq,bufsum,nint)
324              call rakdd_compare3c(
325     &              ish,jsh,ksh,
326     &              (nint*12*12),buffdsq,bufsum,norm)
327              normmax = max(normmax,norm)
328              normmin = min(normmin,norm)
329            endif               ! 4 atoms the same
330          enddo                 ! ksh
331        enddo                   ! jsh
332        string  = ' '
333        stringe = ' '
334        call util_date(stringe)
335        write(string,'(1x,i5,a3,a27,a2,a27)')
336     &        ish,' ::',
337     &        strings(1:25),'::',
338     &        stringe(1:25)
339        do jsh = 1,inp_strlen(string)
340          ksh = ichar(string(jsh:jsh))
341          if (ksh.eq.10) string(jsh:jsh) = ' '
342        enddo
343        write(luout,'(1x,a)')string(1:inp_strlen(string))
344      enddo                     ! ish
345      write(luout,'(1x,a,1pd20.10)')
346     &      ' maximum difference norm over all quartets:',normmax
347      write(luout,'(1x,a,1pd20.10)')
348     &      ' minimum difference norm over all quartets:',normmin
349*
35010000 format(1x,a,4(i4),1x,a,4(i3),1x,a,4(1x,i5),1x,a,i6,/,
351     &      47x,a,4(1x,i5),1x,a,i2)
352      end
353      subroutine rakdd_compare3c(
354     &              ish,jsh,ksh,
355     &              len,buffd,buf,norm)
356      implicit none
357#include "stdio.fh"
358#include "mafdecls.fh"
359      integer ish,jsh,ksh
360      integer len
361      double precision buffd(len), buf(len)
362*
363      double precision ddot
364      external ddot
365*
366      double precision norm, thresh
367      integer h_diff, k_diff
368*
369      thresh = 1.0d-07
370      if (.not.ma_alloc_get(mt_dbl,len,
371     &    'diff buf',h_diff,k_diff)) stop ' ma alloc fail'
372
373      call dcopy(len,buf,1,dbl_mb(k_diff),1)
374      call daxpy(len,-1.0d00,buffd,1,dbl_mb(k_diff),1)
375      norm = ddot(len,dbl_mb(k_diff),1,dbl_mb(k_diff),1)
376*      if (norm.gt.thresh)
377*     &
378      write(luout,*)' ',ish,jsh,ksh,
379     &    ' norm = ',norm
380      if (.not.ma_free_heap(h_diff)) stop ' ma free fail'
381      end
382      subroutine rakdd_printsq(buffsq,nint)
383      implicit none
384      integer nint
385      double precision buffsq (3,4,3,4,nint)
386c
387      integer ia, ja, ix, jx, int
388      logical doit
389c
390      do int = 1,nint
391        do ia = 1,4
392          do ja = 1,4
393            if (ia.le.ja) then
394              do ix = 1,3
395                do jx = 1,3
396                  doit = .true.
397                  if (ia.eq.ja) doit = ix.le.jx
398                  if (doit) then
399                    write(6,10000)ix,ia,jx,ja,int,
400     &                    buffsq(ix,ia,jx,ja,int)
401                  endif
402                enddo
403              enddo
404            endif
405          enddo
406        enddo
407      enddo
40810000 format(1x,'dd(',4(i4,','),i4,') =',1pd20.10)
409      end
410