1C> \ingroup task
2C> @{
3      logical function task_pderiv(rtdb)
4*
5* $Id$
6*
7***********************************************************************
8* This is about the nuclear coordinate derivatives of multicenter
9* integrals that we need.  Notation:
10*  i,j,p  compound basis function (in practice, shell) indices,
11*       including, specification of which center the function is on
12*
13*  Xi, Xj   orbital basis functions (chi), position of center not shown
14*           explicitly here
15*
16*  Gp      density basis function, ditto
17*
18*  R, R2   crystal lattice translation vectors
19*
20*  rnuc    position of some nucleus
21*
22*   1,2     electronic coordinates to be integrated (vol. elements d1,d2)
23*
24* I will just display the integrals whose derivatives we need, not the
25*  differentiation itself.
26*
27* 1. orbital overlap   Int d1 Xi(1) Xj(1-R)   (i think we have this)  eq. 78
28*
29* 2. Kinetic energy    Int d1 Xi(1) [-del^2/2] Xj(1-R)     (ditto)      eq. 30
30*
31* 3. 1-center attraction  -Sum(nuc)Z(nuc) Int d1 Gp(1-R) /|1-rnuc|      eq. 60
32*
33* 4. 2-center repulsion  Int d1 d2 Gp(1) (1/r12) Gp(2-R)  [r12 obvious,
34*     p = some other p index] eq. 52
35*
36* 5. 2-center attraction -Sum(nuc)Z(nuc) Int d1 Xi(1+R2)Xj(1+R2-R)/|1-rnuc|  eq.42
37*
38* 6. 3-center repulsion Int d1 d2 Xi(1)Xj(1-R)(1/r12)Gp(2-R2)         eq. 41
39*
40* All equations refer to our gross paper JCP 105, 10983 (1995),
41* doi:10.1063/1.472866 of which I believe you have a copy of, if not I will
42* bring one.
43*
44* Probably you will need some more details on how we want this, or we can
45* talk to you to understand how we get it, then we can reshuffle things as
46* necessary.  For example do we get the derivative of any integral with
47* repect to any arbitrary nucleus that we specify, or with respect to all
48* nuclei in a single block of data, or only for those nuclei which are
49* relevant (nonzero deriv.) for a given integral?  Call or message me as
50* needed and Ill come over.)
51*
52* Later
53* John
54***********************************************************************
55*
56
57      implicit none
58#include "errquit.fh"
59#include "stdio.fh"
60#include "mafdecls.fh"
61#include "geom.fh"
62#include "bas.fh"
63#include "util.fh"
64c
65      logical int_normalize
66      external int_normalize
67      logical pderiv_compute_1e1cpe
68      external pderiv_compute_1e1cpe
69      logical pderiv_compute_1epe
70      external pderiv_compute_1epe
71      logical pderiv_compute_1eov
72      external pderiv_compute_1eov
73      logical pderiv_compute_1eke
74      external pderiv_compute_1eke
75      logical pderiv_compute_2e2c
76      external pderiv_compute_2e2c
77      logical pderiv_compute_2e3c
78      external pderiv_compute_2e3c
79      logical pderiv_compute_mpole
80      external pderiv_compute_mpole
81      logical pderiv_compute_3ov
82      external pderiv_compute_3ov
83      logical pderiv_compute_2e3ct, pderiv_compute_d2e3ct
84      external pderiv_compute_2e3ct, pderiv_compute_d2e3ct
85c
86      integer rtdb
87c
88      logical status
89      integer basis, geom
90      integer nat, nat3, size, maxg1, maxs1, maxg2, maxs2, maxg, maxs
91      integer nbf, nbfsq
92      integer hbuf, hbufp, hbufm, hscr, h1ec, h1efd, hxyz
93      integer kbuf, kbufp, kbufm, kscr, k1ec, k1efd, kxyz
94c
95      task_pderiv = .false.
96c
97      if (.not.geom_create(geom,'geometry')) call errquit
98     &    ('geom create failed',911, GEOM_ERR)
99      if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit
100     &    ('geom_rtdb_load failed',911, RTDB_ERR)
101c
102      if (.not.bas_create(basis,'ao basis')) call errquit
103     &    ('bas_create failed',911, BASIS_ERR)
104      if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) call errquit
105     &    ('bas_rtdb_load failed',911, RTDB_ERR)
106c
107      write(6,*)' geom/basis loaded'
108c
109      if (.not.int_normalize(rtdb,basis)) stop ' norm error 1'
110c
111      if (.not. bas_print(basis))
112     $    call errquit(' basis print failed', 0, BASIS_ERR)
113c
114      if (.not.bas_numbf(basis,nbf)) call errquit
115     &    ('numbf failed',911, BASIS_ERR)
116c
117      nbfsq = nbf*nbf
118      if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe'
119      write(6,*) 'number of atoms ', nat
120      nat3 = 3*nat
121      size = max(nbf,56)  ! slmax = 35 for mpole
122      size = nat3*nbfsq*size
123      size = 2*size
124c
125      call intd_init(rtdb,1,basis)
126      call int_mem_print()
127      call int_mem_1e(maxg1,maxs1)
128      call int_mem_2e3c(maxg2,maxs2)
129      maxg2 = maxg2*3*4
130      maxg1 = maxg1*35  ! mpole
131      maxg = max(maxg1,maxg2)
132      maxs = max(maxs1,maxs2)
133      maxg = maxg + maxg/10
134      maxs = maxs + maxs/10
135      maxs = 2*maxs
136      maxg = 2*maxg
137      write(6,*)' normal maxg1 ',maxg1
138      write(6,*)' normal maxs1 ',maxs1
139      write(6,*)' normal maxg2 ',maxg2
140      write(6,*)' normal maxs2 ',maxs2
141      write(6,*)' normal maxg  ',maxg
142      write(6,*)' normal maxs  ',maxs
143      status = ma_alloc_get(mt_dbl,maxg,'int buffer' ,hbuf,kbuf)
144      status = status .and.
145     &    ma_alloc_get(mt_dbl,maxg,'int buffer plus',hbufp,kbufp)
146      status = status .and.
147     &    ma_alloc_get(mt_dbl,maxg,'int buffer minus',hbufm,kbufm)
148      status = status .and.
149     &    ma_alloc_get(mt_dbl,maxs,'int scratch',hscr,kscr)
150      status = status .and.
151     &    ma_alloc_get(mt_dbl,size,'block c',h1ec,k1ec)
152      status = status .and.
153     &    ma_alloc_get(mt_dbl,size,'block fd',h1efd,k1efd)
154      status = status .and.
155     &    ma_alloc_get(mt_dbl,3*nat,'my coords',hxyz,kxyz)
156      if (.not.status) stop ' memory alloc failed rak20 (1)'
157c
158      task_pderiv =  pderiv_compute_1e1cpe(
159     &    geom,basis,nbf,nat,
160     &    maxg,maxs,
161     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
162     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
163     &    dbl_mb(kxyz))
164c
165      task_pderiv = task_pderiv .and.
166     &    pderiv_compute_1eov(
167     &    geom,basis,nbf,nat,
168     &    maxg,maxs,
169     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
170     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
171     &    dbl_mb(kxyz))
172c
173      task_pderiv = task_pderiv .and.
174     &    pderiv_compute_1eke(
175     &    geom,basis,nbf,nat,
176     &    maxg,maxs,
177     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
178     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
179     &    dbl_mb(kxyz))
180c
181      task_pderiv = task_pderiv .and.
182     &    pderiv_compute_1epe(
183     &    geom,basis,nbf,nat,
184     &    maxg,maxs,
185     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
186     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
187     &    dbl_mb(kxyz))
188c
189      task_pderiv = task_pderiv .and.
190     &    pderiv_compute_2e2c(
191     &    geom,basis,nbf,nat,
192     &    maxg,maxs,
193     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
194     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
195     &    dbl_mb(kxyz))
196c
197      task_pderiv = task_pderiv .and.
198     &    pderiv_compute_2e3c(
199     &    geom,basis,nbf,nat,
200     &    maxg,maxs,
201     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
202     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
203     &    dbl_mb(kxyz))
204
205      task_pderiv = task_pderiv .and.
206     &    pderiv_compute_3ov(
207     &    geom,basis,nbf,nat,
208     &    maxg,maxs,
209     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
210     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
211     &    dbl_mb(kxyz))
212
213      task_pderiv = task_pderiv .and.
214     &    pderiv_compute_mpole(
215     &    geom,basis,nbf,nat,
216     &    maxg,maxs,
217     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
218     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
219     &    dbl_mb(kxyz))
220
221      task_pderiv = task_pderiv .and.
222     &    pderiv_compute_2e3ct(
223     &    geom,basis,nbf,nat,
224     &    maxg,maxs,
225     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
226     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
227     &    dbl_mb(kxyz))
228
229      task_pderiv = task_pderiv .and.
230     &    pderiv_compute_d2e3ct(
231     &    geom,basis,nbf,nat,
232     &    maxg,maxs,
233     &    dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm),
234     &    dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd),
235     &    dbl_mb(kxyz))
236
237
238
239      call intd_terminate()
240      if (.not.bas_destroy(basis)) stop ' bas destroy fail'
241      if (.not.geom_destroy(geom)) stop ' geom destroy fail'
242      status = .true.
243      status = status.and.ma_free_heap(hbuf)
244      status = status.and.ma_free_heap(hbufp)
245      status = status.and.ma_free_heap(hbufm)
246      status = status.and.ma_free_heap(hscr)
247      status = status.and.ma_free_heap(h1ec)
248      status = status.and.ma_free_heap(h1efd)
249      status = status.and.ma_free_heap(hxyz)
250      task_pderiv = status.and.task_pderiv
251c
252      end
253C> @}
254      logical function pderiv_compute_1e1cpe(
255     &    geom,basis,nbf,nat,lbuf,lscr,
256     &    buf,bufp,bufm,scr,p1c,p1fd,xyz)
257      implicit none
258#include "mafdecls.fh"
259#include "nwc_const.fh"
260#include "geomP.fh"
261#include "bas.fh"
262#include "stdio.fh"
263c
264      double precision ddot
265      external ddot
266c
267      integer geom
268      integer basis
269      integer nbf
270      integer nat
271      integer lbuf
272      integer lscr
273      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
274      double precision scr(lscr)
275      double precision p1c(nbf,3,nat)
276      double precision p1fd(nbf,3,nat)
277      double precision xyz(3,nat)
278c
279      integer nat3
280      integer atom, ixyz
281      double precision delta, factor, thresh, norm
282      double precision R(3)
283      integer nzero1, nzero2
284      integer ishell, nshell, ilo, ihi, nbfsh, cnt, i
285      integer IR
286c
287      call dfill(3,0.0d00,R,1)
288      call dfill((3*nat),0.0d00,xyz,1)
289      call dfill((nbf*3*nat),0.0d00,p1c,1)
290      call dfill((nbf*3*nat),0.0d00,p1fd,1)
291      call dfill(lbuf,0.0d00,buf,1)
292      call dfill(lbuf,0.0d00,bufp,1)
293      call dfill(lbuf,0.0d00,bufm,1)
294      call dfill(lscr,0.0d00,scr,1)
295      do IR = 1,2
296        if (IR.eq.1) call dfill(3,0.0d00,R,1)
297        if (IR.eq.2) then
298          R(1) = 1.0d00
299          R(2) = 2.0d00
300          R(3) = 3.0d00
301        endif
302        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 1'
303        thresh = 1.0d-10
304        delta = 0.0001d00
305        factor = 1.0d00/(2.0d00*delta)
306        nat3 = 3*nat
307* store original coordintates
308        call dcopy(nat3,coords(1,1,geom),1,xyz,1)
309*
310        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
311
312        do ishell = 1, nshell
313          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
314     &        stop 'cn2bfr error i'
315          nbfsh = ihi - ilo + 1
316          do atom = 1,nat
317            do ixyz = 1,3
318              call dcopy(nat3,xyz,1,coords(1,1,geom),1)
319              coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta
320              call dfill(lscr,0.0d00,scr,1)
321              call dfill(lbuf,0.0d00,buf,1)
322              call dfill(lbuf,0.0d00,bufp,1)
323              call intp_1e1cpe(basis,ishell,R,lscr,scr,lbuf,bufp)
324*
325              call dcopy(nat3,xyz,1,coords(1,1,geom),1)
326              coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta
327              call dfill(lscr,0.0d00,scr,1)
328              call dfill(lbuf,0.0d00,bufm,1)
329              call intp_1e1cpe(basis,ishell,R,lscr,scr,lbuf,bufm)
330*
331              call dcopy(nbfsh,bufp,1,buf,1)
332              call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
333              call dscal(nbfsh,factor,buf,1)
334              cnt = 1
335              do i = ilo,ihi
336                p1fd(i,ixyz,atom) = buf(cnt)
337                cnt = cnt + 1
338              enddo
339            enddo
340          enddo
341        enddo
342c
343        write(6,*)' fd: full list '
344        nzero1 = 0
345        do atom = 1,nat
346          do ixyz = 1,3
347            do i = 1,nbf
348              if (abs(p1fd(i,ixyz,atom)).le.thresh) nzero1 = nzero1 + 1
349*              write(6,10000)i,ixyz,atom,p1fd(i,ixyz,atom)
350            enddo
351          enddo
352        enddo
353        write(6,*)' fd: non-zero list '
354        nzero2 = 0
355        do atom = 1,nat
356          do ixyz = 1,3
357            do i = 1,nbf
358              if (abs(p1fd(i,ixyz,atom)).gt.thresh) then
359*               write(6,10000)i,ixyz,atom,p1fd(i,ixyz,atom)
360                continue
361              else
362                nzero2 = nzero2 + 1
363              endif
364            enddo
365          enddo
366        enddo
367        write(6,*)' fd: num zeros :1:', nzero1
368        write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
369     &      (nzero2-nzero1)
370        write(6,*)' fd: non-zero  :1:', ((nbf*3*nat)-nzero1)
371        write(6,*)' fd: non-zero  :2:', ((nbf*3*nat)-nzero2)
372        write(6,*)' fd: possible  : :', (nbf*3*nat)
373c
374        call dcopy(nat3,xyz,1,coords(1,1,geom),1)
375        do ishell = 1,nshell
376          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
377     &        stop 'cn2bfr error i'
378          call dfill(lscr,0.0d00,scr,1)
379          call dfill(lbuf,0.0d00,buf,1)
380          call intpd_1e1cpe(basis,ishell,R,lscr,scr,lbuf,buf)
381          cnt = 1
382          do atom=1,nat
383            do ixyz = 1,3
384              do i = ilo, ihi
385                p1c(i,ixyz,atom) = buf(cnt)
386                cnt = cnt + 1
387              enddo
388            enddo
389          enddo
390        enddo
391        write(6,*)'  c: full list '
392        nzero1 = 0
393        do atom = 1,nat
394          do ixyz = 1,3
395            do i = 1,nbf
396              if (abs(p1c(i,ixyz,atom)).le.thresh) nzero1 = nzero1 + 1
397*              write(6,10001)i,ixyz,atom,p1c(i,ixyz,atom)
398            enddo
399          enddo
400        enddo
401        write(6,*)'  c: non-zero list '
402        nzero2 = 0
403        do atom = 1,nat
404          do ixyz = 1,3
405            do i = 1,nbf
406              if (abs(p1c(i,ixyz,atom)).gt.thresh) then
407*                write(6,10001)i,ixyz,atom,p1c(i,ixyz,atom)
408                continue
409              else
410                nzero2 = nzero2 + 1
411              endif
412            enddo
413          enddo
414        enddo
415        write(6,*)'  c: num zeros :1:', nzero1
416        write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
417     &      (nzero2-nzero1)
418        write(6,*)'  c: non-zero  :1:', ((nbf*3*nat)-nzero1)
419        write(6,*)'  c: non-zero  :2:', ((nbf*3*nat)-nzero2)
420        write(6,*)'  c: possible  : :', (nbf*3*nat)
421c
42210000   format(1x,'p1fd(',i3,',',i2,',',i3,') = ',1pd20.10)
42310001   format(1x,' p1c(',i3,',',i2,',',i3,') = ',1pd20.10)
424c
425        call daxpy((nbf*3*nat),-1.0d00,p1fd,1,p1c,1)
426        norm = ddot((nbf*3*nat),p1c,1,p1c,1)
427        write(luout,*)' 1e1cpe difference norm:',ir,' ',norm
428c
429        pderiv_compute_1e1cpe = norm.lt.thresh
430      enddo
431      end
432      logical function pderiv_compute_1eov(
433     &    geom,basis,nbf,nat,lbuf,lscr,
434     &    buf,bufp,bufm,scr,ovc,ovfd,xyz)
435      implicit none
436#include "mafdecls.fh"
437#include "nwc_const.fh"
438#include "geomP.fh"
439#include "bas.fh"
440#include "stdio.fh"
441c
442      double precision ddot
443      external ddot
444c
445      integer geom
446      integer basis
447      integer nbf
448      integer nat
449      integer lbuf
450      integer lscr
451      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
452      double precision scr(lscr)
453      double precision ovc(nbf,nbf,3,nat)
454      double precision ovfd(nbf,nbf,3,nat)
455      double precision xyz(3,nat)
456c
457      integer nat3
458      integer atom, ixyz, IR
459      double precision delta, factor, thresh, norm
460      double precision R(3)
461      integer nzero1, nzero2
462      integer nshell, cnt, nbfsh, watom(2)
463      integer ishell, ilo, ihi, nbfshi, i
464      integer jshell, jlo, jhi, nbfshj, j
465c
466      call dfill(3,0.0d00,R,1)
467      call dfill((3*nat),0.0d00,xyz,1)
468      call dfill((nbf*nbf*3*nat),0.0d00,ovc,1)
469      call dfill((nbf*nbf*3*nat),0.0d00,ovfd,1)
470      call dfill(lbuf,0.0d00,buf,1)
471      call dfill(lbuf,0.0d00,bufp,1)
472      call dfill(lbuf,0.0d00,bufm,1)
473      call dfill(lscr,0.0d00,scr,1)
474      do IR = 1,2
475        if (IR.eq.1) call dfill(3,0.0d00,R,1)
476        if (IR.eq.2) then
477          R(1) = 1.0d00
478          R(2) = 2.0d00
479          R(3) = 3.0d00
480        endif
481        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 2'
482        thresh = 1.0d-10
483        delta = 0.0001d00
484        factor = 1.0d00/(2.0d00*delta)
485        nat3 = 3*nat
486* store original coordintates
487        call dcopy(nat3,coords(1,1,geom),1,xyz,1)
488*
489        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
490
491        do ishell = 1, nshell
492          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
493     &        stop 'cn2bfr error i'
494          nbfshi = ihi - ilo + 1
495          do jshell = 1, ishell
496            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
497     &          stop 'cn2bfr error j'
498            nbfshj = jhi - jlo + 1
499            nbfsh = nbfshi*nbfshj
500            do atom = 1,nat
501              do ixyz = 1,3
502                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
503                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta
504                call dfill(lscr,0.0d00,scr,1)
505                call dfill(lbuf,0.0d00,buf,1)
506                call dfill(lbuf,0.0d00,bufp,1)
507                call intp_1eov(basis,ishell,basis,jshell,
508     &              R,lscr,scr,lbuf,bufp)
509*
510                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
511                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta
512                call dfill(lscr,0.0d00,scr,1)
513                call dfill(lbuf,0.0d00,bufm,1)
514                call intp_1eov(basis,ishell,basis,jshell,
515     &              R,lscr,scr,lbuf,bufm)
516*
517                call dcopy(nbfsh,bufp,1,buf,1)
518                call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
519                call dscal(nbfsh,factor,buf,1)
520                cnt = 1
521                do i = ilo,ihi
522                  do j = jlo, jhi
523                    ovfd(i,j,ixyz,atom) = buf(cnt)
524                    ovfd(j,i,ixyz,atom) = buf(cnt)
525                    cnt = cnt + 1
526                  enddo
527                enddo
528                if (.not.ma_verify_allocator_stuff())
529     &              stop ' ma broke 3'
530              enddo
531            enddo
532          enddo
533        enddo
534c
535        write(6,*)' fd: full list i,j'
536        nzero1 = 0
537        do atom = 1,nat
538          do ixyz = 1,3
539            do i = 1,nbf
540              do j = 1,nbf
541                if (abs(ovfd(i,j,ixyz,atom)).le.thresh)
542     &              nzero1 = nzero1 + 1
543*                write(6,10000)i,j,ixyz,atom,ovfd(i,j,ixyz,atom)
544              enddo
545            enddo
546          enddo
547        enddo
548        write(6,*)' fd: non-zero list '
549        nzero2 = 0
550        do atom = 1,nat
551          do ixyz = 1,3
552            do i = 1,nbf
553              do j = 1,nbf
554                if (abs(ovfd(i,j,ixyz,atom)).gt.thresh) then
555*                  write(6,10000)i,j,ixyz,atom,ovfd(i,j,ixyz,atom)
556                  continue
557                else
558                  nzero2 = nzero2 + 1
559                endif
560              enddo
561            enddo
562          enddo
563        enddo
564        write(6,*)' fd: num zeros :1:', nzero1
565        write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
566     &      (nzero2-nzero1)
567        write(6,*)' fd: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
568        write(6,*)' fd: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
569        write(6,*)' fd: possible  : :', (nbf*nbf*3*nat)
570        nzero1 = 0
571        nzero2 = 0
572c
573        call dcopy(nat3,xyz,1,coords(1,1,geom),1)
574        do ishell = 1,nshell
575          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
576     &        stop 'cn2bfr error i'
577          do jshell = 1,ishell
578            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
579     &          stop 'cn2bfr error j'
580            call dfill(lscr,0.0d00,scr,1)
581            call dfill(lbuf,0.0d00,buf,1)
582            call intpd_1eov(basis,ishell,basis,jshell,R,lscr,scr,
583     &          lbuf,buf,watom)
584*            write(6,*)'watom 1eov',watom
585            cnt = 1
586            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 4'
587            do atom=1,2
588              if (watom(atom).gt.0) then
589                do ixyz = 1,3
590                  do i = ilo, ihi
591                    do j = jlo, jhi
592*buffer is <jlo:jhi, ilo:ihi, 1:3, 2>
593                      ovc(i,j,ixyz,watom(atom)) = buf(cnt)
594                      ovc(j,i,ixyz,watom(atom)) = buf(cnt)
595                      cnt = cnt + 1
596                    enddo
597                  enddo
598                enddo
599              endif
600            enddo
601            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 5'
602          enddo
603        enddo
604        write(6,*)'  c: full list i,j'
605        nzero1 = 0
606        do atom = 1,nat
607          do ixyz = 1,3
608            do i = 1,nbf
609              do j = 1,nbf
610                if (abs(ovc(i,j,ixyz,atom)).le.thresh)
611     &              nzero1 = nzero1 + 1
612*                write(6,10001)i,j,ixyz,atom,ovc(i,j,ixyz,atom)
613              enddo
614            enddo
615          enddo
616        enddo
617        write(6,*)'  c: non-zero list '
618        nzero2 = 0
619        do atom = 1,nat
620          do ixyz = 1,3
621            do i = 1,nbf
622              do j = 1,nbf
623                if (abs(ovc(i,j,ixyz,atom)).gt.thresh) then
624*                  write(6,10001)i,j,ixyz,atom,ovc(i,j,ixyz,atom)
625                  continue
626                else
627                  nzero2 = nzero2 + 1
628                endif
629              enddo
630            enddo
631          enddo
632        enddo
633        write(6,*)'  c: num zeros :1:', nzero1
634        write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
635     &      (nzero2-nzero1)
636        write(6,*)'  c: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
637        write(6,*)'  c: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
638        write(6,*)'  c: possible  : :', (nbf*nbf*3*nat)
639c
64010000   format(1x,'ovfd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
64110001   format(1x,' ovc(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
642c
643        call daxpy((nbf*nbf*3*nat),-1.0d00,ovfd,1,ovc,1)
644        norm = ddot((nbf*nbf*3*nat),ovc,1,ovc,1)
645        write(luout,*)' 1eov difference norm:',ir,' ',norm
646c
647        pderiv_compute_1eov = norm.lt.thresh
648      enddo
649      end
650      logical function pderiv_compute_1eke(
651     &    geom,basis,nbf,nat,lbuf,lscr,
652     &    buf,bufp,bufm,scr,kec,kefd,xyz)
653      implicit none
654#include "mafdecls.fh"
655#include "nwc_const.fh"
656#include "geomP.fh"
657#include "bas.fh"
658#include "stdio.fh"
659c
660      double precision ddot
661      external ddot
662c
663      integer geom
664      integer basis
665      integer nbf
666      integer nat
667      integer lbuf
668      integer lscr
669      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
670      double precision scr(lscr)
671      double precision kec(nbf,nbf,3,nat)
672      double precision kefd(nbf,nbf,3,nat)
673      double precision xyz(3,nat)
674c
675      integer nat3
676      integer atom, ixyz, IR
677      double precision delta, factor, thresh, norm
678      double precision R(3)
679      integer nzero1, nzero2
680      integer nshell, cnt, nbfsh, watom(2)
681      integer ishell, ilo, ihi, nbfshi, i
682      integer jshell, jlo, jhi, nbfshj, j
683c
684      call dfill(3,0.0d00,R,1)
685      call dfill((3*nat),0.0d00,xyz,1)
686      call dfill((nbf*nbf*3*nat),0.0d00,kec,1)
687      call dfill((nbf*nbf*3*nat),0.0d00,kefd,1)
688      call dfill(lbuf,0.0d00,buf,1)
689      call dfill(lbuf,0.0d00,bufp,1)
690      call dfill(lbuf,0.0d00,bufm,1)
691      call dfill(lscr,0.0d00,scr,1)
692      do IR = 1,2
693        if (IR.eq.1) call dfill(3,0.0d00,R,1)
694        if (IR.eq.2) then
695          R(1) = 1.0d00
696          R(2) = 2.0d00
697          R(3) = 3.0d00
698        endif
699        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 6'
700        thresh = 1.0d-10
701        delta = 0.0001d00
702        factor = 1.0d00/(2.0d00*delta)
703        nat3 = 3*nat
704* store original coordintates
705        call dcopy(nat3,coords(1,1,geom),1,xyz,1)
706*
707        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
708
709        do ishell = 1, nshell
710          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
711     &        stop 'cn2bfr error i'
712          nbfshi = ihi - ilo + 1
713          do jshell = 1, ishell
714            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
715     &          stop 'cn2bfr error j'
716            nbfshj = jhi - jlo + 1
717            nbfsh = nbfshi*nbfshj
718            do atom = 1,nat
719              do ixyz = 1,3
720                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
721                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta
722                call dfill(lscr,0.0d00,scr,1)
723                call dfill(lbuf,0.0d00,buf,1)
724                call dfill(lbuf,0.0d00,bufp,1)
725                call intp_1eke(basis,ishell,basis,jshell,
726     &              R,lscr,scr,lbuf,bufp)
727*
728                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
729                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta
730                call dfill(lscr,0.0d00,scr,1)
731                call dfill(lbuf,0.0d00,bufm,1)
732                call intp_1eke(basis,ishell,basis,jshell,
733     &              R,lscr,scr,lbuf,bufm)
734*
735                call dcopy(nbfsh,bufp,1,buf,1)
736                call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
737                call dscal(nbfsh,factor,buf,1)
738                cnt = 1
739                do i = ilo,ihi
740                  do j = jlo, jhi
741                    kefd(i,j,ixyz,atom) = buf(cnt)
742                    kefd(j,i,ixyz,atom) = buf(cnt)
743                    cnt = cnt + 1
744                  enddo
745                enddo
746                if (.not.ma_verify_allocator_stuff())
747     &              stop ' ma broke 7'
748              enddo
749            enddo
750          enddo
751        enddo
752c
753        write(6,*)' fd: full list i,j'
754        nzero1 = 0
755        do atom = 1,nat
756          do ixyz = 1,3
757            do i = 1,nbf
758              do j = 1,nbf
759                if (abs(kefd(i,j,ixyz,atom)).le.thresh)
760     &              nzero1 = nzero1 + 1
761*                write(6,10000)i,j,ixyz,atom,kefd(i,j,ixyz,atom)
762              enddo
763            enddo
764          enddo
765        enddo
766        write(6,*)' fd: non-zero list '
767        nzero2 = 0
768        do atom = 1,nat
769          do ixyz = 1,3
770            do i = 1,nbf
771              do j = 1,nbf
772                if (abs(kefd(i,j,ixyz,atom)).gt.thresh) then
773*                  write(6,10000)i,j,ixyz,atom,kefd(i,j,ixyz,atom)
774                  continue
775                else
776                  nzero2 = nzero2 + 1
777                endif
778              enddo
779            enddo
780          enddo
781        enddo
782        write(6,*)' fd: num zeros :1:', nzero1
783        write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
784     &      (nzero2-nzero1)
785        write(6,*)' fd: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
786        write(6,*)' fd: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
787        write(6,*)' fd: possible  : :', (nbf*nbf*3*nat)
788        nzero1 = 0
789        nzero2 = 0
790c
791        call dcopy(nat3,xyz,1,coords(1,1,geom),1)
792        do ishell = 1,nshell
793          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
794     &        stop 'cn2bfr error i'
795          do jshell = 1,ishell
796            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
797     &          stop 'cn2bfr error j'
798            call dfill(lscr,0.0d00,scr,1)
799            call dfill(lbuf,0.0d00,buf,1)
800            call intpd_1eke(basis,ishell,basis,jshell,R,lscr,scr,
801     &          lbuf,buf,watom)
802*            write(6,*)'watom 1eke',watom
803            cnt = 1
804            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 8'
805            do atom=1,2
806              if (watom(atom).gt.0) then
807                do ixyz = 1,3
808                  do i = ilo, ihi
809                    do j = jlo, jhi
810                      kec(i,j,ixyz,watom(atom)) = buf(cnt)
811                      kec(j,i,ixyz,watom(atom)) = buf(cnt)
812                      cnt = cnt + 1
813                    enddo
814                  enddo
815                enddo
816              endif
817            enddo
818            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 9'
819          enddo
820        enddo
821        write(6,*)'  c: full list i,j'
822        nzero1 = 0
823        do atom = 1,nat
824          do ixyz = 1,3
825            do i = 1,nbf
826              do j = 1,nbf
827                if (abs(kec(i,j,ixyz,atom)).le.thresh)
828     &              nzero1 = nzero1 + 1
829*                write(6,10001)i,j,ixyz,atom,kec(i,j,ixyz,atom)
830              enddo
831            enddo
832          enddo
833        enddo
834        write(6,*)'  c: non-zero list '
835        nzero2 = 0
836        do atom = 1,nat
837          do ixyz = 1,3
838            do i = 1,nbf
839              do j = 1,nbf
840                if (abs(kec(i,j,ixyz,atom)).gt.thresh) then
841*                  write(6,10001)i,j,ixyz,atom,kec(i,j,ixyz,atom)
842                  continue
843                else
844                  nzero2 = nzero2 + 1
845                endif
846              enddo
847            enddo
848          enddo
849        enddo
850        write(6,*)'  c: num zeros :1:', nzero1
851        write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
852     &      (nzero2-nzero1)
853        write(6,*)'  c: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
854        write(6,*)'  c: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
855        write(6,*)'  c: possible  : :', (nbf*nbf*3*nat)
856c
85710000   format(1x,'kefd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
85810001   format(1x,' kec(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
859c
860        call daxpy((nbf*nbf*3*nat),-1.0d00,kefd,1,kec,1)
861        norm = ddot((nbf*nbf*3*nat),kec,1,kec,1)
862        write(luout,*)' 1eke difference norm:',ir,' ',norm
863c
864        pderiv_compute_1eke = norm.lt.thresh
865      enddo
866      end
867      logical function pderiv_compute_1epe(
868     &    geom,basis,nbf,nat,lbuf,lscr,
869     &    buf,bufp,bufm,scr,pec,pefd,xyz)
870      implicit none
871#include "mafdecls.fh"
872#include "nwc_const.fh"
873#include "geomP.fh"
874#include "bas.fh"
875#include "stdio.fh"
876c
877      double precision ddot
878      external ddot
879c
880      integer geom
881      integer basis
882      integer nbf
883      integer nat
884      integer lbuf
885      integer lscr
886      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
887      double precision scr(lscr)
888      double precision pec(nbf,nbf,3,nat)
889      double precision pefd(nbf,nbf,3,nat)
890      double precision xyz(3,nat)
891c
892      integer nat3
893      integer atom, ixyz, IR
894      double precision delta, factor, thresh, norm
895      double precision R1(3), R2(3)
896      integer nzero1, nzero2
897      integer nshell, cnt, nbfsh
898      integer ishell, ilo, ihi, nbfshi, i
899      integer jshell, jlo, jhi, nbfshj, j
900c
901      call dfill(3,0.0d00,R1,1)
902      call dfill(3,0.0d00,R2,1)
903      call dfill((3*nat),0.0d00,xyz,1)
904      call dfill((nbf*nbf*3*nat),0.0d00,pec,1)
905      call dfill((nbf*nbf*3*nat),0.0d00,pefd,1)
906      call dfill(lbuf,0.0d00,buf,1)
907      call dfill(lbuf,0.0d00,bufp,1)
908      call dfill(lbuf,0.0d00,bufm,1)
909      call dfill(lscr,0.0d00,scr,1)
910      do IR = 1,3
911        if (IR.eq.1) then
912          call dfill(3,0.0d00,R1,1)
913          call dfill(3,0.0d00,R2,1)
914        elseif (IR.eq.2) then
915          call dfill(3,0.0d00,R2,1)
916          R1(1) = 1.0d00
917          R1(2) = 2.0d00
918          R1(3) = 3.0d00
919        elseif (IR.eq.3) then
920          R1(1) = 1.0d00
921          R1(2) = 2.0d00
922          R1(3) = 3.0d00
923          R2(1) = 3.0d00
924          R2(2) = 4.0d00
925          R2(3) = 5.0d00
926        else
927          stop ' how did IR get here'
928        endif
929        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 10'
930        thresh = 1.0d-10
931        delta = 0.0001d00
932        factor = 1.0d00/(2.0d00*delta)
933        nat3 = 3*nat
934* store original coordintates
935        call dcopy(nat3,coords(1,1,geom),1,xyz,1)
936*
937        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
938
939        do ishell = 1, nshell
940          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
941     &        stop 'cn2bfr error i'
942          nbfshi = ihi - ilo + 1
943          do jshell = 1, ishell
944            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
945     &          stop 'cn2bfr error j'
946            nbfshj = jhi - jlo + 1
947            nbfsh = nbfshi*nbfshj
948            do atom = 1,nat
949              do ixyz = 1,3
950                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
951                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta
952                call dfill(lscr,0.0d00,scr,1)
953                call dfill(lbuf,0.0d00,buf,1)
954                call dfill(lbuf,0.0d00,bufp,1)
955                call intp_1epe(basis,ishell,R1,basis,jshell,
956     &              R2,lscr,scr,lbuf,bufp)
957*
958                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
959                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta
960                call dfill(lscr,0.0d00,scr,1)
961                call dfill(lbuf,0.0d00,bufm,1)
962                call intp_1epe(basis,ishell,R1,basis,jshell,
963     &              R2,lscr,scr,lbuf,bufm)
964*
965                call dcopy(nbfsh,bufp,1,buf,1)
966                call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
967                call dscal(nbfsh,factor,buf,1)
968                cnt = 1
969                do i = ilo,ihi
970                  do j = jlo, jhi
971                    pefd(i,j,ixyz,atom) = buf(cnt)
972                    pefd(j,i,ixyz,atom) = buf(cnt)
973                    cnt = cnt + 1
974                  enddo
975                enddo
976                if (.not.ma_verify_allocator_stuff())
977     &              stop ' ma broke 11'
978              enddo
979            enddo
980          enddo
981        enddo
982c
983        write(6,*)' fd: full list i,j'
984        nzero1 = 0
985        do atom = 1,nat
986          do ixyz = 1,3
987            do i = 1,nbf
988              do j = 1,nbf
989                if (abs(pefd(i,j,ixyz,atom)).le.thresh)
990     &              nzero1 = nzero1 + 1
991*                write(6,10000)i,j,ixyz,atom,pefd(i,j,ixyz,atom)
992              enddo
993            enddo
994          enddo
995        enddo
996        write(6,*)' fd: non-zero list '
997        nzero2 = 0
998        do atom = 1,nat
999          do ixyz = 1,3
1000            do i = 1,nbf
1001              do j = 1,nbf
1002                if (abs(pefd(i,j,ixyz,atom)).gt.thresh) then
1003*                  write(6,10000)i,j,ixyz,atom,pefd(i,j,ixyz,atom)
1004                  continue
1005                else
1006                  nzero2 = nzero2 + 1
1007                endif
1008              enddo
1009            enddo
1010          enddo
1011        enddo
1012        write(6,*)' fd: num zeros :1:', nzero1
1013        write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
1014     &      (nzero2-nzero1)
1015        write(6,*)' fd: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
1016        write(6,*)' fd: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
1017        write(6,*)' fd: possible  : :', (nbf*nbf*3*nat)
1018        nzero1 = 0
1019        nzero2 = 0
1020c
1021        call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1022        do ishell = 1,nshell
1023          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1024     &        stop 'cn2bfr error i'
1025          do jshell = 1,ishell
1026            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1027     &          stop 'cn2bfr error j'
1028            call dfill(lscr,0.0d00,scr,1)
1029            call dfill(lbuf,0.0d00,buf,1)
1030            call intpd_1epe(basis,ishell,R1,basis,jshell,R2,
1031     &          lscr,scr,lbuf,buf)
1032            cnt = 1
1033            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 12'
1034            do atom=1,nat
1035              do ixyz = 1,3
1036                do i = ilo, ihi
1037                  do j = jlo, jhi
1038                    pec(i,j,ixyz,atom) = buf(cnt)
1039                    pec(j,i,ixyz,atom) = buf(cnt)
1040                    cnt = cnt + 1
1041                  enddo
1042                enddo
1043              enddo
1044            enddo
1045            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 13'
1046          enddo
1047        enddo
1048        write(6,*)'  c: full list i,j'
1049        nzero1 = 0
1050        do atom = 1,nat
1051          do ixyz = 1,3
1052            do i = 1,nbf
1053              do j = 1,nbf
1054                if (abs(pec(i,j,ixyz,atom)).le.thresh)
1055     &              nzero1 = nzero1 + 1
1056*                write(6,10001)i,j,ixyz,atom,pec(i,j,ixyz,atom)
1057              enddo
1058            enddo
1059          enddo
1060        enddo
1061        write(6,*)'  c: non-zero list '
1062        nzero2 = 0
1063        do atom = 1,nat
1064          do ixyz = 1,3
1065            do i = 1,nbf
1066              do j = 1,nbf
1067                if (abs(pec(i,j,ixyz,atom)).gt.thresh) then
1068*                  write(6,10001)i,j,ixyz,atom,pec(i,j,ixyz,atom)
1069                  continue
1070                else
1071                  nzero2 = nzero2 + 1
1072                endif
1073              enddo
1074            enddo
1075          enddo
1076        enddo
1077        write(6,*)'  c: num zeros :1:', nzero1
1078        write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
1079     &      (nzero2-nzero1)
1080        write(6,*)'  c: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
1081        write(6,*)'  c: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
1082        write(6,*)'  c: possible  : :', (nbf*nbf*3*nat)
1083c
108410000   format(1x,'pefd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
108510001   format(1x,' pec(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
1086c
1087        call daxpy((nbf*nbf*3*nat),-1.0d00,pefd,1,pec,1)
1088        norm = ddot((nbf*nbf*3*nat),pec,1,pec,1)
1089        write(luout,*)' 1epe difference norm:',IR,' ',norm
1090c
1091        pderiv_compute_1epe = norm.lt.thresh
1092      enddo
1093      end
1094      logical function pderiv_compute_2e2c(
1095     &    geom,basis,nbf,nat,lbuf,lscr,
1096     &    buf,bufp,bufm,scr,eri2c,eri2fd,xyz)
1097      implicit none
1098#include "mafdecls.fh"
1099#include "nwc_const.fh"
1100#include "geomP.fh"
1101#include "bas.fh"
1102#include "stdio.fh"
1103c
1104      double precision ddot
1105      external ddot
1106c
1107      integer geom
1108      integer basis
1109      integer nbf
1110      integer nat
1111      integer lbuf
1112      integer lscr
1113      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
1114      double precision scr(lscr)
1115      double precision eri2c(nbf,nbf,3,nat)
1116      double precision eri2fd(nbf,nbf,3,nat)
1117      double precision xyz(3,nat)
1118c
1119      integer nat3
1120      integer atom, ixyz, IR
1121      double precision delta, factor, thresh, norm
1122      double precision R(3)
1123      integer nzero1, nzero2
1124      integer nshell, cnt, nbfsh, watom(2)
1125      integer ishell, ilo, ihi, nbfshi, i
1126      integer jshell, jlo, jhi, nbfshj, j
1127c
1128      call dfill(3,0.0d00,R,1)
1129      call dfill((3*nat),0.0d00,xyz,1)
1130      call dfill((nbf*nbf*3*nat),0.0d00,eri2c,1)
1131      call dfill((nbf*nbf*3*nat),0.0d00,eri2fd,1)
1132      call dfill(lbuf,0.0d00,buf,1)
1133      call dfill(lbuf,0.0d00,bufp,1)
1134      call dfill(lbuf,0.0d00,bufm,1)
1135      call dfill(lscr,0.0d00,scr,1)
1136      do IR = 1,2
1137        if (IR.eq.1) call dfill(3,0.0d00,R,1)
1138        if (IR.eq.2) then
1139          R(1) = 1.0d00
1140          R(2) = 2.0d00
1141          R(3) = 3.0d00
1142        endif
1143        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 14'
1144        thresh = 1.0d-10
1145        delta = 0.0001d00
1146        factor = 1.0d00/(2.0d00*delta)
1147        nat3 = 3*nat
1148* store original coordintates
1149        call dcopy(nat3,coords(1,1,geom),1,xyz,1)
1150*
1151        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
1152
1153        do ishell = 1, nshell
1154          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1155     &        stop 'cn2bfr error i'
1156          nbfshi = ihi - ilo + 1
1157          do jshell = 1, ishell
1158            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1159     &          stop 'cn2bfr error j'
1160            nbfshj = jhi - jlo + 1
1161            nbfsh = nbfshi*nbfshj
1162            do atom = 1,nat
1163              do ixyz = 1,3
1164                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1165                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta
1166                call dfill(lscr,0.0d00,scr,1)
1167                call dfill(lbuf,0.0d00,buf,1)
1168                call dfill(lbuf,0.0d00,bufp,1)
1169                call intp_2e2c(basis,ishell,basis,jshell,
1170     &              R,lscr,scr,lbuf,bufp)
1171*
1172                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1173                coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta
1174                call dfill(lscr,0.0d00,scr,1)
1175                call dfill(lbuf,0.0d00,bufm,1)
1176                call intp_2e2c(basis,ishell,basis,jshell,
1177     &              R,lscr,scr,lbuf,bufm)
1178*
1179                call dcopy(nbfsh,bufp,1,buf,1)
1180                call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
1181                call dscal(nbfsh,factor,buf,1)
1182                cnt = 1
1183                do i = ilo,ihi
1184                  do j = jlo, jhi
1185                    eri2fd(i,j,ixyz,atom) = buf(cnt)
1186                    eri2fd(j,i,ixyz,atom) = buf(cnt)
1187                    cnt = cnt + 1
1188                  enddo
1189                enddo
1190                if (.not.ma_verify_allocator_stuff())
1191     &              stop ' ma broke 15'
1192              enddo
1193            enddo
1194          enddo
1195        enddo
1196c
1197        write(6,*)' fd: full list i,j'
1198        nzero1 = 0
1199        do atom = 1,nat
1200          do ixyz = 1,3
1201            do i = 1,nbf
1202              do j = 1,nbf
1203                if (abs(eri2fd(i,j,ixyz,atom)).le.thresh)
1204     &              nzero1 = nzero1 + 1
1205*                write(6,10000)i,j,ixyz,atom,eri2fd(i,j,ixyz,atom)
1206              enddo
1207            enddo
1208          enddo
1209        enddo
1210        write(6,*)' fd: non-zero list '
1211        nzero2 = 0
1212        do atom = 1,nat
1213          do ixyz = 1,3
1214            do i = 1,nbf
1215              do j = 1,nbf
1216                if (abs(eri2fd(i,j,ixyz,atom)).gt.thresh) then
1217*                  write(6,10000)i,j,ixyz,atom,eri2fd(i,j,ixyz,atom)
1218                  continue
1219                else
1220                  nzero2 = nzero2 + 1
1221                endif
1222              enddo
1223            enddo
1224          enddo
1225        enddo
1226        write(6,*)' fd: num zeros :1:', nzero1
1227        write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
1228     &      (nzero2-nzero1)
1229        write(6,*)' fd: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
1230        write(6,*)' fd: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
1231        write(6,*)' fd: possible  : :', (nbf*nbf*3*nat)
1232        nzero1 = 0
1233        nzero2 = 0
1234c
1235        call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1236        do ishell = 1,nshell
1237          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1238     &        stop 'cn2bfr error i'
1239          do jshell = 1,ishell
1240            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1241     &          stop 'cn2bfr error j'
1242            call dfill(lscr,0.0d00,scr,1)
1243            call dfill(lbuf,0.0d00,buf,1)
1244*            write(6,*)'rak20: ishell, jshell ',ishell,jshell
1245            call intpd_2e2c(basis,ishell,basis,jshell,R,lscr,scr,
1246     &          lbuf,buf,watom)
1247*            write(6,*)'watom 2e2c',watom
1248            cnt = 1
1249            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 16'
1250            do atom=1,2
1251              if (watom(atom).gt.0) then
1252                do ixyz = 1,3
1253                  do i = ilo, ihi
1254                    do j = jlo, jhi
1255                      eri2c(i,j,ixyz,watom(atom)) = buf(cnt)
1256                      eri2c(j,i,ixyz,watom(atom)) = buf(cnt)
1257                      cnt = cnt + 1
1258                    enddo
1259                  enddo
1260                enddo
1261              endif
1262            enddo
1263            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 17'
1264          enddo
1265        enddo
1266        write(6,*)'  c: full list i,j'
1267        nzero1 = 0
1268        do atom = 1,nat
1269          do ixyz = 1,3
1270            do i = 1,nbf
1271              do j = 1,nbf
1272                if (abs(eri2c(i,j,ixyz,atom)).le.thresh)
1273     &              nzero1 = nzero1 + 1
1274*                write(6,10001)i,j,ixyz,atom,eri2c(i,j,ixyz,atom)
1275              enddo
1276            enddo
1277          enddo
1278        enddo
1279        write(6,*)'  c: non-zero list '
1280        nzero2 = 0
1281        do atom = 1,nat
1282          do ixyz = 1,3
1283            do i = 1,nbf
1284              do j = 1,nbf
1285                if (abs(eri2c(i,j,ixyz,atom)).gt.thresh) then
1286*                  write(6,10001)i,j,ixyz,atom,eri2c(i,j,ixyz,atom)
1287                  continue
1288                else
1289                  nzero2 = nzero2 + 1
1290                endif
1291              enddo
1292            enddo
1293          enddo
1294        enddo
1295        write(6,*)'  c: num zeros :1:', nzero1
1296        write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
1297     &      (nzero2-nzero1)
1298        write(6,*)'  c: non-zero  :1:', ((nbf*nbf*3*nat)-nzero1)
1299        write(6,*)'  c: non-zero  :2:', ((nbf*nbf*3*nat)-nzero2)
1300        write(6,*)'  c: possible  : :', (nbf*nbf*3*nat)
1301c
130210000   format(1x,'eri2fd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
130310001   format(1x,' eri2c(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10)
1304c
1305        call daxpy((nbf*nbf*3*nat),-1.0d00,eri2fd,1,eri2c,1)
1306        norm = ddot((nbf*nbf*3*nat),eri2c,1,eri2c,1)
1307        write(luout,*)' 2e2c difference norm:',ir,' ',norm
1308c
1309        pderiv_compute_2e2c = norm.lt.thresh
1310      enddo
1311      end
1312      logical function pderiv_compute_2e3c(
1313     &    geom,basis,nbf,nat,lbuf,lscr,
1314     &    buf,bufp,bufm,scr,eri3c,eri3fd,xyz)
1315      implicit none
1316#include "mafdecls.fh"
1317#include "nwc_const.fh"
1318#include "geomP.fh"
1319#include "bas.fh"
1320#include "stdio.fh"
1321c
1322      double precision ddot
1323      external ddot
1324c
1325      integer geom
1326      integer basis
1327      integer nbf
1328      integer nat
1329      integer lbuf
1330      integer lscr
1331      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
1332      double precision scr(lscr)
1333      double precision eri3c(nbf,nbf,nbf,3,nat)
1334      double precision eri3fd(nbf,nbf,nbf,3,nat)
1335      double precision xyz(3,nat)
1336c
1337      integer nat3, nintz
1338      integer atom, ixyz, IR
1339      double precision delta, factor, thresh, norm
1340      double precision R1(3), R2(3)
1341      integer nzero1, nzero2
1342      integer nshell, cnt, nbfsh, watom(4)
1343      integer ishell, ilo, ihi, nbfshi, i
1344      integer jshell, jlo, jhi, nbfshj, j
1345      integer kshell, klo, khi, nbfshk, k
1346c
1347      call dfill(3,0.0d00,R1,1)
1348      call dfill(3,0.0d00,R2,1)
1349      call dfill((3*nat),0.0d00,xyz,1)
1350      call dfill((nbf*nbf*nbf*3*nat),0.0d00,eri3c,1)
1351      call dfill((nbf*nbf*nbf*3*nat),0.0d00,eri3fd,1)
1352      call dfill(lbuf,0.0d00,buf,1)
1353      call dfill(lbuf,0.0d00,bufp,1)
1354      call dfill(lbuf,0.0d00,bufm,1)
1355      call dfill(lscr,0.0d00,scr,1)
1356      do IR = 1,3
1357        if (IR.eq.1) then
1358          call dfill(3,0.0d00,R1,1)
1359          call dfill(3,0.0d00,R2,1)
1360        elseif (IR.eq.2) then
1361          call dfill(3,0.0d00,R2,1)
1362          R1(1) = 1.0d00
1363          R1(2) = 2.0d00
1364          R1(3) = 3.0d00
1365        elseif (IR.eq.3) then
1366          R1(1) = 1.0d00
1367          R1(2) = 2.0d00
1368          R1(3) = 3.0d00
1369          R2(1) = 3.0d00
1370          R2(2) = 4.0d00
1371          R2(3) = 5.0d00
1372        else
1373          stop ' how did IR get here'
1374        endif
1375        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18'
1376        thresh = 1.0d-10
1377        delta = 0.0001d00
1378        factor = 1.0d00/(2.0d00*delta)
1379        nat3 = 3*nat
1380* store original coordintates
1381        call dcopy(nat3,coords(1,1,geom),1,xyz,1)
1382*
1383        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
1384
1385        do ishell = 1, nshell
1386          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1387     &        stop 'cn2bfr error i'
1388          nbfshi = ihi - ilo + 1
1389          do jshell = 1, nshell
1390            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1391     &          stop 'cn2bfr error j'
1392            nbfshj = jhi - jlo + 1
1393            do kshell = 1, jshell
1394              if (.not.bas_cn2bfr(basis,kshell,klo,khi))
1395     &            stop 'cn2bfr error k'
1396              nbfshk = khi - klo + 1
1397              nbfsh = nbfshi*nbfshj*nbfshk
1398              do atom = 1,nat
1399                do ixyz = 1,3
1400                  call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1401                  coords(ixyz,atom,geom) =
1402     &                coords(ixyz,atom,geom) + delta
1403                  call dfill(lscr,0.0d00,scr,1)
1404                  call dfill(lbuf,0.0d00,buf,1)
1405                  call dfill(lbuf,0.0d00,bufp,1)
1406                  call intp_2e3c(basis,ishell,basis,jshell,kshell,
1407     &              R1,R2,lscr,scr,lbuf,bufp)
1408*
1409                  call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1410                  coords(ixyz,atom,geom) =
1411     &                coords(ixyz,atom,geom) - delta
1412                  call dfill(lscr,0.0d00,scr,1)
1413                  call dfill(lbuf,0.0d00,bufm,1)
1414                  call intp_2e3c(basis,ishell,basis,jshell,kshell,
1415     &              R1,R2,lscr,scr,lbuf,bufm)
1416*
1417                  call dcopy(nbfsh,bufp,1,buf,1)
1418                  call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
1419                  call dscal(nbfsh,factor,buf,1)
1420                  cnt = 1
1421                  do i = ilo,ihi
1422                    do j = jlo, jhi
1423                      do k = klo, khi
1424                        eri3fd(i,j,k,ixyz,atom) = buf(cnt)
1425                        eri3fd(i,k,j,ixyz,atom) = buf(cnt)
1426                        cnt = cnt + 1
1427                      enddo
1428                    enddo
1429                  enddo
1430                  if (.not.ma_verify_allocator_stuff())
1431     &                stop ' ma broke 19'
1432                enddo
1433              enddo
1434            enddo
1435          enddo
1436        enddo
1437c
1438        write(6,*)' fd: full list i,j'
1439        nzero1 = 0
1440        do atom = 1,nat
1441          do ixyz = 1,3
1442            do i = 1,nbf
1443              do j = 1,nbf
1444                do k = 1,nbf
1445                  if (abs(eri3fd(i,j,k,ixyz,atom)).le.thresh)
1446     &                nzero1 = nzero1 + 1
1447*                  write(6,10000)i,j,k,ixyz,atom,eri3fd(i,j,k,ixyz,atom)
1448                enddo
1449              enddo
1450            enddo
1451          enddo
1452        enddo
1453        write(6,*)' fd: non-zero list '
1454        nzero2 = 0
1455        do atom = 1,nat
1456          do ixyz = 1,3
1457            do i = 1,nbf
1458              do j = 1,nbf
1459                do k = 1,nbf
1460                  if (abs(eri3fd(i,j,k,ixyz,atom)).gt.thresh) then
1461*                    write(6,10000)i,j,k,ixyz,atom,
1462*     &                  eri3fd(i,j,k,ixyz,atom)
1463                    continue
1464                  else
1465                    nzero2 = nzero2 + 1
1466                  endif
1467                enddo
1468              enddo
1469            enddo
1470          enddo
1471        enddo
1472        write(6,*)' fd: num zeros :1:', nzero1
1473        write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
1474     &      (nzero2-nzero1)
1475        write(6,*)' fd: non-zero  :1:', ((nbf*nbf*nbf*3*nat)-nzero1)
1476        write(6,*)' fd: non-zero  :2:', ((nbf*nbf*nbf*3*nat)-nzero2)
1477        write(6,*)' fd: possible  : :', (nbf*nbf*nbf*3*nat)
1478        nzero1 = 0
1479        nzero2 = 0
1480c
1481        call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1482        do ishell = 1,nshell
1483          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1484     &        stop 'cn2bfr error i'
1485          do jshell = 1,nshell
1486            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1487     &          stop 'cn2bfr error j'
1488            do kshell = 1, jshell
1489              if (.not.bas_cn2bfr(basis,kshell,klo,khi))
1490     &            stop 'cn2bfr error k'
1491              call dfill(lscr,0.0d00,scr,1)
1492              call dfill(lbuf,0.0d00,buf,1)
1493*              write(6,*)'rak20: ishell, jshell, kshell ',
1494*     &            ishell,jshell,kshell
1495              call intpd_2e3c(basis,ishell,basis,jshell,kshell,R1,R2,
1496     &            lscr,scr,lbuf,buf,watom)
1497*              write(6,*)'watom 2e3c',watom
1498              if (.not.ma_verify_allocator_stuff()) stop ' ma broke 20'
1499              nintz = (ihi-ilo+1)*(jhi-jlo+1)*(khi-klo+1)
1500              do atom=1,4
1501                if (watom(atom).gt.0) then
1502                  cnt = ((atom-1)*nintz*3) + 1
1503                  do ixyz = 1,3
1504                    do i = ilo, ihi
1505                      do j = jlo, jhi
1506                        do k = klo, khi
1507                          eri3c(i,j,k,ixyz,watom(atom)) = buf(cnt)
1508                          eri3c(i,k,j,ixyz,watom(atom)) = buf(cnt)
1509                          cnt = cnt + 1
1510                        enddo
1511                      enddo
1512                    enddo
1513                  enddo
1514                endif
1515              enddo
1516              if (.not.ma_verify_allocator_stuff()) stop ' ma broke 21'
1517            enddo
1518          enddo
1519        enddo
1520        write(6,*)'  c: full list i,j'
1521        nzero1 = 0
1522        do atom = 1,nat
1523          do ixyz = 1,3
1524            do i = 1,nbf
1525              do j = 1,nbf
1526                do k = 1,nbf
1527                  if (abs(eri3c(i,j,k,ixyz,atom)).le.thresh)
1528     &                nzero1 = nzero1 + 1
1529*                  write(6,10001)i,j,k,ixyz,atom,eri3c(i,j,k,ixyz,atom)
1530                enddo
1531              enddo
1532            enddo
1533          enddo
1534        enddo
1535        write(6,*)'  c: non-zero list '
1536        nzero2 = 0
1537        do atom = 1,nat
1538          do ixyz = 1,3
1539            do i = 1,nbf
1540              do j = 1,nbf
1541                do k = 1,nbf
1542                  if (abs(eri3c(i,j,k,ixyz,atom)).gt.thresh) then
1543*                    write(6,10001)i,j,k,ixyz,atom,
1544*     &                  eri3c(i,j,k,ixyz,atom)
1545                    continue
1546                  else
1547                    nzero2 = nzero2 + 1
1548                  endif
1549                enddo
1550              enddo
1551            enddo
1552          enddo
1553        enddo
1554        write(6,*)'  c: num zeros :1:', nzero1
1555        write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
1556     &      (nzero2-nzero1)
1557        write(6,*)'  c: non-zero  :1:', ((nbf*nbf*nbf*3*nat)-nzero1)
1558        write(6,*)'  c: non-zero  :2:', ((nbf*nbf*nbf*3*nat)-nzero2)
1559        write(6,*)'  c: possible  : :', (nbf*nbf*nbf*3*nat)
1560c
156110000   format(1x,'eri3fd(',i3,',',i3,',',i3,',',i2,',',i3,
1562     &      ') = ',1pd20.10)
156310001   format(1x,' eri3c(',i3,',',i3,',',i3,',',i2,',',i3,
1564     &      ') = ',1pd20.10)
1565c
1566        call daxpy((nbf*nbf*nbf*3*nat),-1.0d00,eri3fd,1,eri3c,1)
1567        norm = ddot((nbf*nbf*nbf*3*nat),eri3c,1,eri3c,1)
1568        write(luout,*)' 2e3c difference norm:',ir,' ',norm
1569c
1570        pderiv_compute_2e3c = norm.lt.thresh
1571      enddo
1572      end
1573      logical function pderiv_compute_3ov(
1574     &    geom,basis,nbf,nat,lbuf,lscr,
1575     &    buf,bufp,bufm,scr,ov3c,ov3fd,xyz)
1576      implicit none
1577#include "mafdecls.fh"
1578#include "nwc_const.fh"
1579#include "geomP.fh"
1580#include "bas.fh"
1581#include "stdio.fh"
1582c
1583      double precision ddot
1584      external ddot
1585c
1586      integer geom
1587      integer basis
1588      integer nbf
1589      integer nat
1590      integer lbuf
1591      integer lscr
1592      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
1593      double precision scr(lscr)
1594      double precision ov3c(nbf,nbf,nbf,3,nat)
1595      double precision ov3fd(nbf,nbf,nbf,3,nat)
1596      double precision xyz(3,nat)
1597c
1598      integer nat3, nintz
1599      integer atom, ixyz
1600      double precision delta, factor, thresh, norm
1601*rak:      double precision zeta(50)
1602*rak:      integer ztype, znp, zng, zsph
1603      integer nzero1, nzero2
1604      integer nshell, cnt, nbfsh, watom(3)
1605      integer ishell, ilo, ihi, nbfshi, i
1606      integer jshell, jlo, jhi, nbfshj, j
1607      integer kshell, klo, khi, nbfshk, k
1608c
1609      call dfill((3*nat),0.0d00,xyz,1)
1610      call dfill((nbf*nbf*nbf*3*nat),0.0d00,ov3c,1)
1611      call dfill((nbf*nbf*nbf*3*nat),0.0d00,ov3fd,1)
1612      call dfill(lbuf,0.0d00,buf,1)
1613      call dfill(lbuf,0.0d00,bufp,1)
1614      call dfill(lbuf,0.0d00,bufm,1)
1615      call dfill(lscr,0.0d00,scr,1)
1616      if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18'
1617      thresh = 1.0d-10
1618      delta = 0.0001d00
1619      factor = 1.0d00/(2.0d00*delta)
1620      nat3 = 3*nat
1621* store original coordintates
1622      call dcopy(nat3,coords(1,1,geom),1,xyz,1)
1623*
1624      if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
1625*
1626*rak:      if (.not.bas_continfo(basis,1,ztype,znp,zng,zsph)) stop 'e1'
1627*rak:      if (znp.gt.50) stop ' zeta dimension too small '
1628*rak:      call dfill(znp,0.0d00,zeta,1)
1629*rak:      if (.not.bas_get_exponent(basis,1,zeta)) stop 'bas_get_e failed'
1630*rak:      zeta(1) = 1.0d-50
1631*rak:      if (.not.bas_set_exponent(basis,1,zeta,znp))
1632*rak:     &    stop 'bas_set_e failed'
1633*
1634      do ishell = 1, nshell
1635        if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1636     &      stop 'cn2bfr error i'
1637        nbfshi = ihi - ilo + 1
1638        do jshell = 1, nshell
1639          if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1640     &        stop 'cn2bfr error j'
1641          nbfshj = jhi - jlo + 1
1642          do kshell = 1, nshell
1643            if (.not.bas_cn2bfr(basis,kshell,klo,khi))
1644     &          stop 'cn2bfr error k'
1645            nbfshk = khi - klo + 1
1646            nbfsh = nbfshi*nbfshj*nbfshk
1647            do atom = 1,nat
1648              do ixyz = 1,3
1649                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1650                coords(ixyz,atom,geom) =
1651     &              coords(ixyz,atom,geom) + delta
1652                call dfill(lscr,0.0d00,scr,1)
1653                call dfill(lbuf,0.0d00,buf,1)
1654                call dfill(lbuf,0.0d00,bufp,1)
1655
1656                call int_1e3ov(
1657     &              basis,ishell,basis,jshell,basis,kshell,
1658     &              lscr,scr,lbuf,bufp)
1659*
1660                call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1661                coords(ixyz,atom,geom) =
1662     &              coords(ixyz,atom,geom) - delta
1663                call dfill(lscr,0.0d00,scr,1)
1664                call dfill(lbuf,0.0d00,bufm,1)
1665                call int_1e3ov(
1666     &              basis,ishell,basis,jshell,basis,kshell,
1667     &              lscr,scr,lbuf,bufm)
1668*
1669                call dcopy(nbfsh,bufp,1,buf,1)
1670                call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
1671                call dscal(nbfsh,factor,buf,1)
1672                cnt = 1
1673                do i = ilo,ihi
1674                  do j = jlo, jhi
1675                    do k = klo, khi
1676                      ov3fd(i,j,k,ixyz,atom) = buf(cnt)
1677                      cnt = cnt + 1
1678                    enddo
1679                  enddo
1680                enddo
1681                if (.not.ma_verify_allocator_stuff())
1682     &              stop ' ma broke 19'
1683              enddo
1684            enddo
1685          enddo
1686        enddo
1687      enddo
1688c
1689      write(6,*)' fd: full list i,j'
1690      nzero1 = 0
1691      do atom = 1,nat
1692        do ixyz = 1,3
1693          do i = 1,nbf
1694            do j = 1,nbf
1695              do k = 1,nbf
1696                if (abs(ov3fd(i,j,k,ixyz,atom)).le.thresh)
1697     &              nzero1 = nzero1 + 1
1698*                  write(6,10000)i,j,k,ixyz,atom,ov3fd(i,j,k,ixyz,atom)
1699              enddo
1700            enddo
1701          enddo
1702        enddo
1703      enddo
1704      write(6,*)' fd: non-zero list '
1705      nzero2 = 0
1706      do atom = 1,nat
1707        do ixyz = 1,3
1708          do i = 1,nbf
1709            do j = 1,nbf
1710              do k = 1,nbf
1711                if (abs(ov3fd(i,j,k,ixyz,atom)).gt.thresh) then
1712*                    write(6,10000)i,j,k,ixyz,atom,
1713*     &                  ov3fd(i,j,k,ixyz,atom)
1714                  continue
1715                else
1716                  nzero2 = nzero2 + 1
1717                endif
1718              enddo
1719            enddo
1720          enddo
1721        enddo
1722      enddo
1723      write(6,*)' fd: num zeros :1:', nzero1
1724      write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
1725     &    (nzero2-nzero1)
1726      write(6,*)' fd: non-zero  :1:', ((nbf*nbf*nbf*3*nat)-nzero1)
1727      write(6,*)' fd: non-zero  :2:', ((nbf*nbf*nbf*3*nat)-nzero2)
1728      write(6,*)' fd: possible  : :', (nbf*nbf*nbf*3*nat)
1729      nzero1 = 0
1730      nzero2 = 0
1731c
1732      call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1733      do ishell = 1,nshell
1734        if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1735     &      stop 'cn2bfr error i'
1736        do jshell = 1,nshell
1737          if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1738     &        stop 'cn2bfr error j'
1739          do kshell = 1, jshell
1740            if (.not.bas_cn2bfr(basis,kshell,klo,khi))
1741     &          stop 'cn2bfr error k'
1742            call dfill(lscr,0.0d00,scr,1)
1743            call dfill(lbuf,0.0d00,buf,1)
1744*              write(6,*)'rak20: ishell, jshell, kshell ',
1745*     &            ishell,jshell,kshell
1746            call intd_1e3ov(
1747     &          basis,ishell,basis,jshell,basis,kshell,
1748     &          lscr,scr,lbuf,buf,watom)
1749*              write(6,*)'watom 3ov',watom
1750            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 20'
1751            nintz = (ihi-ilo+1)*(jhi-jlo+1)*(khi-klo+1)
1752            do atom=1,3
1753              if (watom(atom).gt.0) then
1754                cnt = ((atom-1)*nintz*3) + 1
1755                do ixyz = 1,3
1756                  do i = ilo, ihi
1757                    do j = jlo, jhi
1758                      do k = klo, khi
1759                        ov3c(i,j,k,ixyz,watom(atom)) = buf(cnt)
1760                        ov3c(i,k,j,ixyz,watom(atom)) = buf(cnt)
1761                        cnt = cnt + 1
1762                      enddo
1763                    enddo
1764                  enddo
1765                enddo
1766              endif
1767            enddo
1768            if (.not.ma_verify_allocator_stuff()) stop ' ma broke 21'
1769          enddo
1770        enddo
1771      enddo
1772      write(6,*)'  c: full list i,j'
1773      nzero1 = 0
1774      do atom = 1,nat
1775        do ixyz = 1,3
1776          do i = 1,nbf
1777            do j = 1,nbf
1778              do k = 1,nbf
1779                if (abs(ov3c(i,j,k,ixyz,atom)).le.thresh)
1780     &              nzero1 = nzero1 + 1
1781*                  write(6,10001)i,j,k,ixyz,atom,ov3c(i,j,k,ixyz,atom)
1782              enddo
1783            enddo
1784          enddo
1785        enddo
1786      enddo
1787      write(6,*)'  c: non-zero list '
1788      nzero2 = 0
1789      do atom = 1,nat
1790        do ixyz = 1,3
1791          do i = 1,nbf
1792            do j = 1,nbf
1793              do k = 1,nbf
1794                if (abs(ov3c(i,j,k,ixyz,atom)).gt.thresh) then
1795*                    write(6,10001)i,j,k,ixyz,atom,
1796*     &                  ov3c(i,j,k,ixyz,atom)
1797                  continue
1798                else
1799                  nzero2 = nzero2 + 1
1800                endif
1801              enddo
1802            enddo
1803          enddo
1804        enddo
1805      enddo
1806      write(6,*)'  c: num zeros :1:', nzero1
1807      write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
1808     &    (nzero2-nzero1)
1809      write(6,*)'  c: non-zero  :1:', ((nbf*nbf*nbf*3*nat)-nzero1)
1810      write(6,*)'  c: non-zero  :2:', ((nbf*nbf*nbf*3*nat)-nzero2)
1811      write(6,*)'  c: possible  : :', (nbf*nbf*nbf*3*nat)
1812c
181310000 format(1x,'ov3fd(',i3,',',i3,',',i3,',',i2,',',i3,
1814     &    ') = ',1pd20.10)
181510001 format(1x,' ov3c(',i3,',',i3,',',i3,',',i2,',',i3,
1816     &    ') = ',1pd20.10)
1817c
1818      call daxpy((nbf*nbf*nbf*3*nat),-1.0d00,ov3fd,1,ov3c,1)
1819      norm = ddot((nbf*nbf*nbf*3*nat),ov3c,1,ov3c,1)
1820      write(luout,*)' 3ov difference norm: ',norm
1821c
1822      pderiv_compute_3ov = norm.lt.thresh
1823      end
1824      logical function pderiv_compute_mpole(
1825     &    geom,basis,nbf,nat,lbuf,lscr,
1826     &    buf,bufp,bufm,scr,mp3c,mp3fd,xyz)
1827      implicit none
1828#include "mafdecls.fh"
1829#include "nwc_const.fh"
1830#include "geomP.fh"
1831#include "bas.fh"
1832#include "stdio.fh"
1833c
1834      double precision ddot
1835      external ddot
1836c
1837      integer geom
1838      integer basis
1839      integer nbf
1840      integer nat
1841      integer lbuf
1842      integer lscr
1843      integer lval_max
1844      parameter (lval_max = 5)
1845      integer slmax
1846*      parameter (slmax = 1)  ! for lval_max = 0
1847*      parameter (slmax = 4)  ! for lval_max = 1
1848*      parameter (slmax = 35) ! ! for lval_max = 4
1849      parameter (slmax = 56) ! 1 +3+6+10+15+21 = 35+21 = 56
1850      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
1851      double precision scr(lscr)
1852      double precision mp3c(slmax,nbf,nbf,3,(nat+1))   ! add one for the multipole center
1853      double precision mp3fd(slmax,nbf,nbf,3,(nat+1))  ! ditto
1854      double precision xyz(3,nat)
1855c
1856      integer nat3, nintz
1857      integer atom, ixyz, IR, lval, nbflval
1858      double precision delta, factor, thresh, norm
1859      double precision R1(3), l_center(3), l_orgcenter(3)
1860      integer nzero1, nzero2
1861      integer nshell, cnt, nbfsh, watom(3), nshell_bas
1862      integer ishell, ilo, ihi, nbfshi, i
1863      integer jshell, jlo, jhi, nbfshj, j
1864      integer k, klo, khi, num_ret
1865      integer h_diff, k_diff
1866      integer imp_cent
1867c
1868      integer pole_index
1869      pole_index(lval) = lval*((lval+3)*(lval+3)+2)/6 + 1
1870c
1871      call dfill(3,0.0d00,R1,1)
1872      call dfill((3*nat),0.0d00,xyz,1)
1873      call dfill((slmax*nbf*nbf*3*(nat+1)),0.0d00,mp3c,1)
1874      call dfill((slmax*nbf*nbf*3*(nat+1)),0.0d00,mp3fd,1)
1875      call dfill(lbuf,0.0d00,buf,1)
1876      call dfill(lbuf,0.0d00,bufp,1)
1877      call dfill(lbuf,0.0d00,bufm,1)
1878      call dfill(lscr,0.0d00,scr,1)
1879      do imp_cent = 1, 2
1880        if (imp_cent.eq.1) then
1881          call dfill(3,0.0d00,l_center,1) ! start with multipole center at origin
1882        else if (imp_cent.eq.2) then
1883          call dcopy(3,coords(1,1,geom),1,l_center,1)
1884        endif
1885        call dcopy(3,l_center,1,l_orgcenter,1)
1886
1887      do IR = 1, 3
1888        if (IR.eq.1) then
1889          call dfill(3,0.0d00,R1,1)
1890        elseif (IR.eq.2) then
1891          R1(1) = 1.0d00
1892          R1(2) = 2.0d00
1893          R1(3) = 3.0d00
1894        elseif (IR.eq.3) then
1895          R1(1) = 1.0d00
1896          R1(2) = 1.0d00
1897          R1(3) = 1.0d00
1898        else
1899          stop ' how did IR get here'
1900        endif
1901        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18'
1902        thresh = 1.0d-10
1903        delta = 0.0001d00
1904        factor = 1.0d00/(2.0d00*delta)
1905*rak:        write(6,*)' factor = ',factor
1906        nat3 = 3*nat
1907* store original coordintates
1908        call dcopy(nat3,coords(1,1,geom),1,xyz,1)
1909*rak:        write(6,*)' original coords '
1910*rak:        call output(coords(1,1,geom),1,3,1,nat,3,nat,1)
1911*rak:        write(6,*)' l_center',l_center
1912*rak:        write(6,*)' l_orgcenter',l_orgcenter
1913*
1914        if (.not.bas_numcont(basis,nshell_bas)) stop 'bas_numcont'
1915        nshell = nshell_bas
1916        nshell = 1 ! tmp change
1917*rak:        write(6,*)' nshell = ',nshell
1918
1919        do ishell = 1, nshell
1920          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
1921     &        stop 'cn2bfr error i'
1922          nbfshi = ihi - ilo + 1
1923          do jshell = 1, nshell
1924            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
1925     &          stop 'cn2bfr error j'
1926            nbfshj = jhi - jlo + 1
1927            do lval = 0,lval_max
1928*rak:              write(6,*)'................................1',
1929*rak:     &            ishell,jshell,lval
1930              nbflval = (lval+1)*(lval+2)/2
1931              nbfsh = nbfshj*nbfshi*nbflval
1932              klo = pole_index((lval-1)) + 1
1933              klo = max(klo,1)  ! minumum value of klo is 1
1934              khi = pole_index(lval)
1935*rak:              write(6,*)' klo, khi',klo,khi
1936              do atom = 1,nat
1937                do ixyz = 1,3
1938                  call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1939                  coords(ixyz,atom,geom) =
1940     &                coords(ixyz,atom,geom) + delta
1941                  call dfill(lscr,0.0d00,scr,1)
1942                  call dfill(lbuf,0.0d00,buf,1)
1943                  call dfill(lbuf,0.0d00,bufp,1)
1944*rak:                  write(6,*)' coords intp_mpolel(1)',atom,ixyz
1945*rak:                  call output(coords(1,1,geom),1,3,1,nat,3,nat,1)
1946*rak:                  write(6,*)' l_center',l_center
1947*rak:                  write(6,*)' l_orgcenter',l_orgcenter,' -*'
1948                  call intp_mpolel(basis,ishell,basis,jshell,
1949     &              R1,lval,l_orgcenter,lscr,scr,lbuf,bufp,num_ret)
1950*                  call intp_1eov(basis,ishell,basis,jshell,
1951*     &                R1,lscr,scr,lbuf,buf)
1952*                  cnt = 0
1953*                  write(6,*)
1954*     &                ' i   j  cnt  mpole  ovl   diff + atom ixyz'
1955*                  do i = ilo,ihi
1956*                    do j = jlo,jhi
1957*                      cnt = cnt + 1
1958*                      norm = abs(bufp(cnt)-buf(cnt))
1959*                      write(6,*)i,j,cnt,bufp(cnt),buf(cnt),norm,
1960*     &                    atom, ixyz
1961*                    enddo
1962*                  enddo
1963*                  norm = 0.0d00
1964*                  call dfill(lbuf,0.0d00,buf,1)
1965*
1966                  call dcopy(nat3,xyz,1,coords(1,1,geom),1)
1967                  coords(ixyz,atom,geom) =
1968     &                coords(ixyz,atom,geom) - delta
1969                  call dfill(lscr,0.0d00,scr,1)
1970                  call dfill(lbuf,0.0d00,bufm,1)
1971*rak:                  write(6,*)' coords intp_mpolel(2)',atom,ixyz
1972*rak:                  call output(coords(1,1,geom),1,3,1,nat,3,nat,1)
1973*rak:                  write(6,*)' l_center',l_center
1974*rak:                  write(6,*)' l_orgcenter',l_orgcenter,' -*'
1975                  call intp_mpolel(basis,ishell,basis,jshell,
1976     &              R1,lval,l_orgcenter,lscr,scr,lbuf,bufm,num_ret)
1977*rak:                  call intp_1eov(basis,ishell,basis,jshell,
1978*rak:     &                R1,lscr,scr,lbuf,buf)
1979*rak:                  cnt = 0
1980*rak:                  write(6,*)
1981*rak:     &                ' i   j  cnt  mpole  ovl   diff - atom ixyz'
1982*rak:                  do i = ilo,ihi
1983*rak:                    do j = jlo,jhi
1984*rak:                      cnt = cnt + 1
1985*rak:                      norm = abs(bufp(cnt)-buf(cnt))
1986*rak:                      write(6,*)i,j,cnt,bufp(cnt),buf(cnt),norm,
1987*rak:     &                    atom,ixyz
1988*rak:                    enddo
1989*rak:                  enddo
1990*rak:                  norm = 0.0d00
1991*rak:                  call dfill(lbuf,0.0d00,buf,1)
1992*
1993*rak:                  call prak3mat(bufp,nbfsh,nbfshi,nbfshj,nbflval,
1994*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
1995*rak:     &                'bufp in atom move',atom,ixyz)
1996*rak:                  call prak3mat(bufm,nbfsh,nbfshi,nbfshj,nbflval,
1997*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
1998*rak:     &                'bufm in atom move',atom,ixyz)
1999                  call dcopy(nbfsh,bufp,1,buf,1)
2000                  call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
2001*rak:                  call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval,
2002*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
2003*rak:     &                'buf before scale in atom move',atom,ixyz)
2004                  call dscal(nbfsh,factor,buf,1)
2005*rak:                  call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval,
2006*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
2007*rak:     &                'buf in atom move',atom,ixyz)
2008                  cnt = 1
2009                  do i = ilo,ihi
2010                    do j = jlo, jhi
2011                      do k = klo, khi
2012                        mp3fd(k,j,i,ixyz,atom) = buf(cnt)
2013*rak:                        write(70,*)
2014*rak:     &                      'fd<',k,j,i,ixyz,atom,'>=',buf(cnt)
2015                        cnt = cnt + 1
2016                      enddo
2017                    enddo
2018                  enddo
2019*rak:                  write(6,*)' nbfsh ',nbfsh,' cnt ',cnt
2020                  if (.not.ma_verify_allocator_stuff())
2021     &                stop ' ma broke 19'
2022                enddo
2023              enddo
2024*             now do multipole center move
2025              call dcopy(nat3,xyz,1,coords(1,1,geom),1)  ! restore original atom coords
2026              do ixyz = 1,3
2027                  call dcopy(3,l_orgcenter,1,l_center,1)
2028                  l_center(ixyz) = l_center(ixyz) + delta
2029                  call dfill(lscr,0.0d00,scr,1)
2030                  call dfill(lbuf,0.0d00,buf,1)
2031                  call dfill(lbuf,0.0d00,bufp,1)
2032*rak:                  write(6,*)' coords intp_mpolel(3)',atom,ixyz
2033*rak:                  call output(coords(1,1,geom),1,3,1,nat,3,nat,1)
2034*rak:                  write(6,*)' l_center',l_center,' -*'
2035*rak:                  write(6,*)' l_orgcenter',l_orgcenter
2036                  call intp_mpolel(basis,ishell,basis,jshell,
2037     &              R1,lval,l_center,lscr,scr,lbuf,bufp,num_ret)
2038*
2039                  call dcopy(3,l_orgcenter,1,l_center,1)
2040                  l_center(ixyz) = l_center(ixyz) - delta
2041                  call dfill(lscr,0.0d00,scr,1)
2042                  call dfill(lbuf,0.0d00,bufm,1)
2043*rak:                  write(6,*)' coords intp_mpolel(4)',atom,ixyz
2044*rak:                  call output(coords(1,1,geom),1,3,1,nat,3,nat,1)
2045*rak:                  write(6,*)' l_center',l_center,' -*'
2046*rak:                  write(6,*)' l_orgcenter',l_orgcenter
2047                  call intp_mpolel(basis,ishell,basis,jshell,
2048     &              R1,lval,l_center,lscr,scr,lbuf,bufm,num_ret)
2049*
2050*rak:                  call prak3mat(bufp,nbfsh,nbfshi,nbfshj,nbflval,
2051*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
2052*rak:     &                'bufp',(nat+1),ixyz)
2053*rak:                  call prak3mat(bufm,nbfsh,nbfshi,nbfshj,nbflval,
2054*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
2055*rak:     &                'bufm',(nat+1),ixyz)
2056                  call dcopy(nbfsh,bufp,1,buf,1)
2057                  call daxpy(nbfsh,-1.0d00,bufm,1,buf,1)
2058*rak:                  call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval,
2059*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
2060*rak:     &                'buf b4 scale in mp center move',(nat+1),ixyz)
2061                  call dscal(nbfsh,factor,buf,1)
2062*rak:                  call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval,
2063*rak:     &                ilo,ihi,jlo,jhi,klo,khi,
2064*rak:     &                'buf in mp center move',(nat+1),ixyz)
2065                  cnt = 1
2066                  do i = ilo,ihi
2067                    do j = jlo, jhi
2068                      do k = klo, khi
2069                        mp3fd(k,j,i,ixyz,nat+1) = buf(cnt)
2070*rak:                        write(70,*)
2071*rak:     &                      'fd<',k,j,i,ixyz,(nat+1),'>=',buf(cnt)
2072                        cnt = cnt + 1
2073                      enddo
2074                    enddo
2075                  enddo
2076                  if (.not.ma_verify_allocator_stuff())
2077     &                stop ' ma broke 19'
2078
2079              enddo
2080            enddo
2081          enddo
2082        enddo
2083c
2084        write(6,*)' fd: full list i,j'
2085        nzero1 = 0
2086        do atom = 1,nat+1
2087          do ixyz = 1,3
2088            do i = 1,nbf
2089              do j = 1,nbf
2090                do k = 1,slmax
2091                  if (abs(mp3fd(k,j,i,ixyz,atom)).le.thresh)
2092     &                nzero1 = nzero1 + 1
2093*                  write(6,10000)k,j,i,ixyz,atom,mp3fd(k,j,i,ixyz,atom)
2094                enddo
2095              enddo
2096            enddo
2097          enddo
2098        enddo
2099        write(6,*)' fd: non-zero list '
2100        nzero2 = 0
2101        do atom = 1,nat+1
2102          do ixyz = 1,3
2103            do i = 1,nbf
2104              do j = 1,nbf
2105                do k = 1,slmax
2106                  if (abs(mp3fd(k,j,i,ixyz,atom)).gt.thresh) then
2107*                    write(6,10000)k,j,i,ixyz,atom,
2108*     &                  mp3fd(k,j,i,ixyz,atom)
2109                    continue
2110                  else
2111                    nzero2 = nzero2 + 1
2112                  endif
2113                enddo
2114              enddo
2115            enddo
2116          enddo
2117        enddo
2118        write(6,*)' fd: num zeros :1:', nzero1
2119        write(6,*)' fd: num zeros :2:', nzero2, ':delta:',
2120     &      (nzero2-nzero1)
2121        write(6,*)' fd: non-zero  :1:',
2122     &      ((nbf*nbf*slmax*3*(nat+1))-nzero1)
2123        write(6,*)' fd: non-zero  :2:',
2124     &      ((nbf*nbf*slmax*3*(nat+1))-nzero2)
2125        write(6,*)' fd: possible  : :',
2126     &      (nbf*nbf*slmax*3*(nat+1))
2127        nzero1 = 0
2128        nzero2 = 0
2129c
2130        call dcopy(nat3,xyz,1,coords(1,1,geom),1)
2131        do ishell = 1,nshell
2132          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
2133     &        stop 'cn2bfr error i'
2134          do jshell = 1,nshell
2135            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
2136     &          stop 'cn2bfr error j'
2137            do lval = 0, lval_max
2138*rak:              write(6,*)'................................2',
2139*rak:     &            ishell,jshell,lval
2140              klo = pole_index((lval-1)) + 1
2141              klo = max(klo,1)  ! minumum value of klo is 1
2142              khi = pole_index(lval)
2143*rak:              write(6,*)' klo, khi',klo,khi
2144c
2145              call dfill(lscr,0.0d00,scr,1)
2146              call dfill(lbuf,0.0d00,buf,1)
2147              call dfill(lbuf,0.0d00,bufp,1)
2148*rak:              write(6,*)'rak20: ishell, jshell,lval ',
2149*rak:     &            ishell,jshell,lval
2150*rak:              write(6,*)' lscr ', lscr
2151*rak:              call intpd_1eov(basis,ishell,basis,jshell,R1,lscr,scr,
2152*rak:     &            lbuf,bufp,watom)
2153*rak:              write(6,*)' coords intpd_mpolel'
2154*rak:              call output(coords(1,1,geom),1,3,1,nat,3,nat,1)
2155*rak:              write(6,*)' l_center',l_center
2156*rak:              write(6,*)' l_orgcenter',l_orgcenter,' -*'
2157              call intpd_mpolel(basis,ishell,basis,jshell,R1,
2158     &            lval,l_orgcenter,
2159     &            lscr,scr,lbuf,buf,num_ret,watom)
2160*rak:              cnt = 0
2161*rak:              write(6,*) ' atom ixyz  i   j  cnt  mpole  ovl   diff '
2162*rak:              do atom = 1,nat+1
2163*rak:                do ixyz = 1,3
2164*rak:                  do i = ilo,ihi
2165*rak:                    do j = jlo,jhi
2166*rak:                      cnt = cnt + 1
2167*rak:                      norm = abs(bufp(cnt)-buf(cnt))
2168*rak:                      write(6,*)
2169*rak:     &                    atom,ixyz,i,j,cnt,buf(cnt),bufp(cnt),norm
2170*rak:                    enddo
2171*rak:                  enddo
2172*rak:                enddo
2173*rak:              enddo
2174*rak:              call dfill(lbuf,0.0d00,bufp,1)
2175*rak:              write(6,*)'watom mpole',watom
2176              if (.not.ma_verify_allocator_stuff()) stop ' ma broke 20'
2177              nintz = (ihi-ilo+1)*(jhi-jlo+1)*((lval+1)*(lval+2)/2)
2178*rak:              write(6,*)'nintz ',nintz
2179*rak:              write(6,*)'ilo/hi',ilo,ihi
2180*rak:              write(6,*)'jlo/hi',jlo,jhi
2181*rak:              write(6,*)'klo/hi',klo,khi
2182              do atom=1,3
2183*rak:                write(6,*)' atom,watom(atom)',atom,watom(atom)
2184                if (watom(atom).gt.0) then
2185                  cnt = ((atom-1)*nintz*3) + 1
2186                  do ixyz = 1,3
2187*rak:                    k_diff = cnt + (ixyz-1)*nintz
2188*rak:                    call prak3mat(buf(k_diff),nbfsh,
2189*rak:     &                  nbfshi,nbfshj,nbflval,
2190*rak:     &                  ilo,ihi,jlo,jhi,klo,khi,
2191*rak:     &                  'buf in derivs call(1)',watom(atom),ixyz)
2192                    do i = ilo, ihi
2193                      do j = jlo, jhi
2194                        do k = klo, khi
2195                          mp3c(k,j,i,ixyz,watom(atom)) = buf(cnt)
2196*rak:                          write(71,*)
2197*rak:     &                        ' c<',k,j,i,ixyz,watom(atom),'>=',buf(cnt)
2198                          cnt = cnt + 1
2199                        enddo
2200                      enddo
2201                    enddo
2202                  enddo
2203                endif
2204                if (watom(atom).eq.-3) then
2205                  cnt = ((atom-1)*nintz*3) + 1
2206                  do ixyz = 1,3
2207*rak:                    k_diff = cnt + (ixyz-1)*nintz
2208*rak:                    call prak3mat(buf(k_diff),nbfsh,
2209*rak:     &                  nbfshi,nbfshj,nbflval,
2210*rak:     &                  ilo,ihi,jlo,jhi,klo,khi,
2211*rak:     &                  'buf in derivs call(2)',watom(atom),ixyz)
2212                    do i = ilo, ihi
2213                      do j = jlo, jhi
2214                        do k = klo, khi
2215                          mp3c(k,j,i,ixyz,nat+1) = buf(cnt)
2216*rak:                          write(71,*)
2217*rak:     &                        ' c<',k,j,i,ixyz,(nat+1),'>=',buf(cnt)
2218                          cnt = cnt + 1
2219                        enddo
2220                      enddo
2221                    enddo
2222                  enddo
2223                endif
2224              enddo
2225              if (.not.ma_verify_allocator_stuff()) stop ' ma broke 21'
2226            enddo
2227          enddo
2228        enddo
2229        write(6,*)'  c: full list i,j'
2230        nzero1 = 0
2231        do atom = 1,nat+1
2232          do ixyz = 1,3
2233            do i = 1,nbf
2234              do j = 1,nbf
2235                do k = 1,slmax
2236                  if (abs(mp3c(k,j,i,ixyz,atom)).le.thresh)
2237     &                nzero1 = nzero1 + 1
2238*                  write(6,10001)k,j,i,ixyz,atom,mp3c(k,j,i,ixyz,atom)
2239                enddo
2240              enddo
2241            enddo
2242          enddo
2243        enddo
2244        write(6,*)'  c: non-zero list '
2245        nzero2 = 0
2246        do atom = 1,nat+1
2247          do ixyz = 1,3
2248            do i = 1,nbf
2249              do j = 1,nbf
2250                do k = 1,slmax
2251                  if (abs(mp3c(k,j,i,ixyz,atom)).gt.thresh) then
2252*                    write(6,10001)k,j,i,ixyz,atom,
2253*     &                  mp3c(k,j,i,ixyz,atom)
2254                    continue
2255                  else
2256                    nzero2 = nzero2 + 1
2257                  endif
2258                enddo
2259              enddo
2260            enddo
2261          enddo
2262        enddo
2263        write(6,*)'  c: num zeros :1:', nzero1
2264        write(6,*)'  c: num zeros :2:', nzero2, ':delta:',
2265     &      (nzero2-nzero1)
2266        write(6,*)'  c: non-zero  :1:',
2267     &      ((nbf*nbf*slmax*3*(nat+1))-nzero1)
2268        write(6,*)'  c: non-zero  :2:',
2269     &      ((nbf*nbf*slmax*3*(nat+1))-nzero2)
2270        write(6,*)'  c: possible  : :',
2271     &      (nbf*nbf*slmax*3*(nat+1))
2272c
227310000   format(1x,'mp3fd(',i3,',',i3,',',i3,',',i2,',',i3,
2274     &      ') = ',1pd20.10)
227510001   format(1x,' mp3c(',i3,',',i3,',',i3,',',i2,',',i3,
2276     &      ') = ',1pd20.10)
2277c
2278        if (.not.ma_alloc_get(mt_dbl,
2279     &      (nbf*nbf*slmax*3*(nat+1)),
2280     &      'diff buffer',
2281     &      h_diff,
2282     &      k_diff)) stop 'mpole: ma alloc failed'
2283        call dcopy ((nbf*nbf*slmax*3*(nat+1)),mp3c,1,dbl_mb(k_diff),1)
2284        call  daxpy((nbf*nbf*slmax*3*(nat+1)),
2285     &      -1.0d00,mp3fd,1,dbl_mb(k_diff),1)
2286        norm = ddot((nbf*nbf*slmax*3*(nat+1)),
2287     &      dbl_mb(k_diff),1,dbl_mb(k_diff),1)
2288        write(luout,*)' mpole difference norm:ir=',ir,
2289     &      ':imp_cent= ',imp_cent,': ',norm
2290c
2291        pderiv_compute_mpole = norm.lt.thresh
2292        if (.not.pderiv_compute_mpole) then
2293          do atom = 1,(nat+1)
2294            do ixyz = 1,3
2295              nintz = (atom-1)*3*slmax*nbf*nbf
2296              nintz = nintz + (ixyz-1)*slmax*nbf*nbf
2297              call prak3mat(dbl_mb(k_diff + nintz),
2298     &            (nbf*nbf*slmax),
2299     &            nbf,nbf,slmax,
2300     &            1,nbf,1,nbf,1,slmax,
2301     &            'diff matrix',atom,ixyz)
2302            enddo
2303          enddo
2304        endif
2305        if (.not.ma_free_heap(h_diff)) stop 'mpole: ma free failed '
2306      enddo
2307      enddo
2308      end
2309      subroutine prak3mat(buf,nbf,nbfi,nbfj,nbfk,
2310     &    ilo,ihi,jlo,jhi,klo,khi,msg,atom,ixyz)
2311      implicit none
2312      integer atom
2313      integer ixyz
2314      integer nbf, nbfi, nbfj, nbfk
2315      integer ilo,ihi,jlo,jhi,klo,khi
2316      character*(*) msg
2317      double precision buf(nbfk,nbfj,nbfi)
2318c
2319      integer i, j, k
2320c
2321      write(6,*)' '
2322      write(6,*)' prak3mat<',msg,'> atom=',atom,' ixyz=',ixyz
2323c
2324      if (nbf.ne.(nbfk*nbfj*nbfi)) then
2325        write(6,*)' nbf error '
2326        write(6,*)' nbf  = ',nbf
2327        write(6,*)' nbfi = ',nbfi
2328        write(6,*)' nbfj = ',nbfj
2329        write(6,*)' nbfk = ',nbfk
2330      endif
2331      do i = ilo,ihi
2332        do j = jlo,jhi
2333          do k = klo,khi
2334            if (abs(buf(k,j,i)).gt.1.0d-07) then
2335              write(6,10000)k,j,i,buf(k,j,i)
2336            endif
2337          enddo
2338        enddo
2339      enddo
2340      write(6,*)' '
234110000 format(1x,'buf(',i3,',',i3,',',i3,')=',1pd20.10)
2342      end
2343      logical function pderiv_compute_2e3ct(
2344     &    geom,basis,nbf,nat,lbuf,lscr,
2345     &    buf,bufp,bufm,scr,eri3,eri3t,xyz)
2346      implicit none
2347#include "mafdecls.fh"
2348#include "nwc_const.fh"
2349#include "geomP.fh"
2350#include "bas.fh"
2351#include "stdio.fh"
2352c
2353      double precision ddot
2354      external ddot
2355c
2356      integer geom
2357      integer basis
2358      integer nbf
2359      integer nat
2360      integer lbuf
2361      integer lscr
2362      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
2363      double precision scr(lscr)
2364      double precision eri3t(nbf,nbf,nbf)
2365      double precision eri3(nbf,nbf,nbf)
2366      double precision xyz(3,nat)
2367c
2368      integer IR
2369      double precision thresh, norm
2370      double precision RJ(3), RK(3)
2371      integer nshell, cnt, nbfsh
2372      integer ishell, ilo, ihi, nbfshi, i
2373      integer jshell, jlo, jhi, nbfshj, j
2374      integer kshell, klo, khi, nbfshk, k
2375      integer nzero, nzerot, npos, npost
2376      integer h_diff, k_diff
2377c
2378      call dfill(3,0.0d00,RJ,1)
2379      call dfill(3,0.0d00,RK,1)
2380      call dfill((3*nat),0.0d00,xyz,1)
2381      call dfill((nbf*nbf*nbf),0.0d00,eri3t,1)
2382      call dfill((nbf*nbf*nbf),0.0d00,eri3,1)
2383      call dfill(lbuf,0.0d00,buf,1)
2384      call dfill(lbuf,0.0d00,bufp,1)
2385      call dfill(lbuf,0.0d00,bufm,1)
2386      call dfill(lscr,0.0d00,scr,1)
2387      do IR = 1,3
2388        if (IR.eq.1) then
2389          call dfill(3,0.0d00,RJ,1)
2390          call dfill(3,0.0d00,RK,1)
2391        elseif (IR.eq.2) then
2392          call dfill(3,0.0d00,RK,1)
2393          RJ(1) = 1.0d00
2394          RJ(2) = 2.0d00
2395          RJ(3) = 3.0d00
2396        elseif (IR.eq.3) then
2397          RJ(1) = 1.0d00
2398          RJ(2) = 2.0d00
2399          RJ(3) = 3.0d00
2400          RK(1) = 3.0d00
2401          RK(2) = 4.0d00
2402          RK(3) = 5.0d00
2403        else
2404          stop ' how did IR get here'
2405        endif
2406        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18'
2407        thresh = 1.0d-10
2408*
2409        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
2410
2411        npos   = 0
2412        nzero  = 0
2413        nzerot = 0
2414        npost  = 0
2415        do ishell = 1, nshell ! basis functon NOT translated
2416          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
2417     &        stop 'cn2bfr error i'
2418          nbfshi = ihi - ilo + 1
2419          do jshell = 1, nshell ! basis functon translated
2420            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
2421     &          stop 'cn2bfr error j'
2422            nbfshj = jhi - jlo + 1
2423            do kshell = 1, nshell   ! fitting function
2424              if (.not.bas_cn2bfr(basis,kshell,klo,khi))
2425     &            stop 'cn2bfr error k'
2426              nbfshk = khi - klo + 1
2427              nbfsh = nbfshi*nbfshj*nbfshk
2428              call dfill(lscr,0.0d00,scr,1)
2429              call dfill(lbuf,0.0d00,buf,1)
2430              call intp_2e3c(basis,kshell,basis,jshell,ishell,
2431     &              Rk,Rj,lscr,scr,lbuf,buf)
2432
2433              cnt = 1
2434              do k = klo, khi
2435                do j = jlo, jhi
2436                  do i = ilo, ihi
2437*                    write(70,*)i,j,k,buf(cnt)
2438                    npos = npos + 1
2439                    if (abs(buf(cnt)).lt.thresh) nzero = nzero+1
2440                    eri3(i,j,k) = buf(cnt)
2441                    cnt = cnt + 1
2442                  enddo
2443                enddo
2444              enddo
2445
2446              call dfill(lscr,0.0d00,scr,1)
2447              call dfill(lbuf,0.0d00,buf,1)
2448              call intp_2e3ct(basis,ishell,jshell,basis,kshell,
2449     &            RJ,RK,lscr,scr,lbuf,buf)
2450*
2451              cnt = 1
2452              do i = ilo, ihi
2453                do j = jlo, jhi
2454                  do k = klo, khi
2455*                    write(71,*)i,j,k,buf(cnt)
2456                    npost = npost + 1
2457                    if (abs(buf(cnt)).lt.thresh) nzerot=nzerot+1
2458                    eri3t(i,j,k) = buf(cnt)
2459                    cnt = cnt + 1
2460                  enddo
2461                enddo
2462              enddo
2463c
2464            enddo
2465          enddo
2466        enddo
2467c
2468        if (.not.ma_alloc_get(mt_dbl,
2469     &      (nbf*nbf*nbf),
2470     &      'diff buffer',
2471     &      h_diff,
2472     &      k_diff)) stop '2e3ct: ma alloc failed'
2473        call dcopy((nbf*nbf*nbf),eri3,1,dbl_mb(k_diff),1)
2474        call daxpy((nbf*nbf*nbf),-1.0d00,eri3t,1,dbl_mb(k_diff),1)
2475        norm = ddot((nbf*nbf*nbf),dbl_mb(k_diff),1,dbl_mb(k_diff),1)
2476        write(luout,*)' 2e3ct info:'
2477        write(luout,*)' number possible in regular:',npos
2478        write(luout,*)' number zero     in regular:',nzero
2479        write(luout,*)' number non-zero in regular:',(npos-nzero)
2480        write(luout,*)' number possible in transpo:',npost
2481        write(luout,*)' number zero     in transpo:',nzerot
2482        write(luout,*)' number non-zero in transpo:',(npost-nzerot)
2483        write(luout,*)' 2e3ct difference norm:',ir,' ',norm
2484c
2485        pderiv_compute_2e3ct = norm.lt.thresh
2486        if (.not.pderiv_compute_2e3ct) then
2487          call prak2e3c(nbf,eri3,eri3t,dbl_mb(k_diff))
2488        endif
2489        if (.not.ma_free_heap(h_diff)) stop '2e3ct: ma free failed '
2490      enddo
2491      end
2492      subroutine prak2e3c(nbf,eri3,eri3t,diff)
2493      implicit none
2494      integer nbf
2495      double precision eri3 (nbf,nbf,nbf)
2496      double precision eri3t(nbf,nbf,nbf)
2497      double precision diff (nbf,nbf,nbf)
2498c
2499      integer i,j,k
2500c
2501      do i = 1,nbf
2502        do j = 1,nbf
2503          do k = 1,nbf
2504            if (abs(diff(i,j,k)).gt.1.0d-07) then
2505              write(6,10000)i,j,k,
2506     &            eri3(i,j,k),
2507     &            eri3t(i,j,k),
2508     &            diff(i,j,k)
2509            endif
2510          enddo
2511        enddo
2512      enddo
251310000 format(1x,'<',i3,',',i3,',',i3,'> =',3(1pd20.10))
2514      end
2515      logical function pderiv_compute_d2e3ct(
2516     &    geom,basis,nbf,nat,lbuf,lscr,
2517     &    buf,bufp,bufm,scr,eri3,eri3t,xyz)
2518      implicit none
2519#include "mafdecls.fh"
2520#include "nwc_const.fh"
2521#include "geomP.fh"
2522#include "bas.fh"
2523#include "stdio.fh"
2524c
2525      double precision ddot
2526      external ddot
2527c
2528      integer geom
2529      integer basis
2530      integer nbf
2531      integer nat
2532      integer lbuf
2533      integer lscr
2534      double precision buf(lbuf), bufp(lbuf), bufm(lbuf)
2535      double precision scr(lscr)
2536      double precision eri3t(nbf,nbf,nbf,3,nat)
2537      double precision eri3(nbf,nbf,nbf,3,nat)
2538      double precision xyz(3,nat)
2539c
2540      integer nat3, nintz
2541      integer atom, ixyz, IR
2542      double precision thresh, norm
2543      double precision RJ(3), RK(3)
2544      integer nshell, cnt, nbfsh, watom(4)
2545      integer ishell, ilo, ihi, nbfshi, i
2546      integer jshell, jlo, jhi, nbfshj, j
2547      integer kshell, klo, khi, nbfshk, k
2548      integer nzero, nzerot, npos, npost
2549      integer h_diff, k_diff
2550c
2551      nat3 = 3*nat
2552      call dfill(3,0.0d00,RJ,1)
2553      call dfill(3,0.0d00,RK,1)
2554      call dfill((nat3),0.0d00,xyz,1)
2555      call dfill((nbf*nbf*nbf*nat3),0.0d00,eri3t,1)
2556      call dfill((nbf*nbf*nbf*nat3),0.0d00,eri3,1)
2557      call dfill(lbuf,0.0d00,buf,1)
2558      call dfill(lbuf,0.0d00,bufp,1)
2559      call dfill(lbuf,0.0d00,bufm,1)
2560      call dfill(lscr,0.0d00,scr,1)
2561      do IR = 1,3
2562        if (IR.eq.1) then
2563          call dfill(3,0.0d00,RJ,1)
2564          call dfill(3,0.0d00,RK,1)
2565        elseif (IR.eq.2) then
2566          call dfill(3,0.0d00,RK,1)
2567          RJ(1) = 1.0d00
2568          RJ(2) = 2.0d00
2569          RJ(3) = 3.0d00
2570        elseif (IR.eq.3) then
2571          RJ(1) = 1.0d00
2572          RJ(2) = 2.0d00
2573          RJ(3) = 3.0d00
2574          RK(1) = 3.0d00
2575          RK(2) = 4.0d00
2576          RK(3) = 5.0d00
2577        else
2578          stop ' how did IR get here'
2579        endif
2580        if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18'
2581        thresh = 1.0d-10
2582*
2583        if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont'
2584
2585        npos   = 0
2586        nzero  = 0
2587        nzerot = 0
2588        npost  = 0
2589        do ishell = 1, nshell ! basis functon NOT translated
2590          if (.not.bas_cn2bfr(basis,ishell,ilo,ihi))
2591     &        stop 'cn2bfr error i'
2592          nbfshi = ihi - ilo + 1
2593          do jshell = 1, nshell ! basis functon translated
2594            if (.not.bas_cn2bfr(basis,jshell,jlo,jhi))
2595     &          stop 'cn2bfr error j'
2596            nbfshj = jhi - jlo + 1
2597            do kshell = 1, nshell   ! fitting function
2598              if (.not.bas_cn2bfr(basis,kshell,klo,khi))
2599     &            stop 'cn2bfr error k'
2600              nbfshk = khi - klo + 1
2601              nbfsh = nbfshi*nbfshj*nbfshk
2602              call dfill(lscr,0.0d00,scr,1)
2603              call dfill(lbuf,0.0d00,buf,1)
2604              call intpd_2e3c(basis,kshell,basis,jshell,ishell,
2605     &              Rk,Rj,lscr,scr,lbuf,buf,watom)
2606
2607              nintz = nbfsh
2608              do atom = 1,4
2609                if (watom(atom).gt.0) then
2610                  cnt = ((atom-1)*nintz*3) + 1
2611                  do ixyz = 1,3
2612                    do k = klo, khi
2613                      do j = jlo, jhi
2614                        do i = ilo, ihi
2615                          npos = npos + 1
2616                          if (abs(buf(cnt)).lt.thresh) nzero = nzero+1
2617*                          write(70,*)i,j,k,ixyz,watom(atom),buf(cnt)
2618                          eri3(i,j,k,ixyz,watom(atom)) = buf(cnt)
2619                          cnt = cnt + 1
2620                        enddo
2621                      enddo
2622                    enddo
2623                  enddo
2624                endif
2625              enddo
2626
2627              call dfill(lscr,0.0d00,scr,1)
2628              call dfill(lbuf,0.0d00,buf,1)
2629              call intpd_2e3ct(basis,ishell,jshell,basis,kshell,
2630     &            RJ,RK,lscr,scr,lbuf,buf,watom)
2631*
2632              do atom = 1,4
2633                if (watom(atom).gt.0) then
2634                  cnt = ((atom-1)*nintz*3) + 1
2635                  do ixyz = 1,3
2636                    do i = ilo, ihi
2637                      do j = jlo, jhi
2638                        do k = klo, khi
2639                          npost = npost + 1
2640                          if (abs(buf(cnt)).lt.thresh) nzerot=nzerot+1
2641*                          write(71,*)i,j,k,ixyz,watom(atom),buf(cnt)
2642                          eri3t(i,j,k,ixyz,watom(atom)) = buf(cnt)
2643                          cnt = cnt + 1
2644                        enddo
2645                      enddo
2646                    enddo
2647                  enddo
2648                endif
2649              enddo
2650c
2651            enddo
2652          enddo
2653        enddo
2654c
2655        if (.not.ma_alloc_get(mt_dbl,
2656     &      (nbf*nbf*nbf*nat3),
2657     &      'diff buffer',
2658     &      h_diff,
2659     &      k_diff)) stop 'd2e3ct: ma alloc failed'
2660        call dcopy((nbf*nbf*nbf*nat3),eri3,1,dbl_mb(k_diff),1)
2661        call daxpy((nbf*nbf*nbf*nat3),-1.0d00,eri3t,1,dbl_mb(k_diff),1)
2662        norm =
2663     &      ddot((nbf*nbf*nbf*nat3),dbl_mb(k_diff),1,dbl_mb(k_diff),1)
2664        write(luout,*)' d2e3ct info:'
2665        write(luout,*)' number possible in regular:',npos
2666        write(luout,*)' number zero     in regular:',nzero
2667        write(luout,*)' number non-zero in regular:',(npos-nzero)
2668        write(luout,*)' number possible in transpo:',npost
2669        write(luout,*)' number zero     in transpo:',nzerot
2670        write(luout,*)' number non-zero in transpo:',(npost-nzerot)
2671        write(luout,*)' d2e3ct difference norm:',ir,' ',norm
2672c
2673        pderiv_compute_d2e3ct = norm.lt.thresh
2674        if (.not.pderiv_compute_d2e3ct) then
2675          call prakd2e3c(nbf,nat,eri3,eri3t,dbl_mb(k_diff))
2676        endif
2677        if (.not.ma_free_heap(h_diff)) stop 'd2e3ct: ma free failed '
2678      enddo
2679      end
2680      subroutine prakd2e3c(nbf,nat,eri3,eri3t,diff)
2681      implicit none
2682      integer nbf
2683      integer nat
2684      double precision eri3(nbf,nbf,nbf,3,nat)
2685      double precision eri3t(nbf,nbf,nbf,3,nat)
2686      double precision diff(nbf,nbf,nbf,3,nat)
2687c
2688      integer atom,ixyz,i,j,k
2689c
2690      do atom = 1,nat
2691        do ixyz = 1,3
2692          do i = 1,nbf
2693            do j = 1,nbf
2694              do k = 1,nbf
2695                if (abs(diff(i,j,k,ixyz,atom)).gt.1.0d-07) then
2696                  write(6,10000)i,j,k,ixyz,atom,
2697     &            eri3(i,j,k,ixyz,atom),
2698     &            eri3t(i,j,k,ixyz,atom),
2699     &            diff(i,j,k,ixyz,atom)
2700                endif
2701              enddo
2702            enddo
2703          enddo
2704        enddo
2705      enddo
270610000 format(1x,'<',i3,',',i3,',',i3,',',i3,',',i3,'> =',3(1pd20.10))
2707      end
2708