1      subroutine dim_elfcon(basis,geom,g_dens,points,npt,elfval,nder)
2c
3c $Id: hnd_elfcon.F 19707 2010-10-29 17:59:36Z d3y133 $
4c
5c     This routine calculates the electronic contribution of the
6c     electronic integral defined by nder for a given density at
7c     the grid points defined in points.
8c
9c     It returns an array (max(nder*3,1),npts) which holds all
10c     max(nder*3,1) components for each grid point
11c
12      implicit none
13#include "nwc_const.fh"
14#include "errquit.fh"
15#include "global.fh"
16#include "bas.fh"
17#include "mafdecls.fh"
18#include "geom.fh"
19#include "stdio.fh"
20#include "msgids.fh"
21c
22      integer basis    ! [input] basis set
23      integer geom     ! [input] geometry
24      integer g_dens   ! [input] GA with density
25      integer npt      ! [input] number of coord points
26      integer nder     ! [input] electronic integral type
27      double precision points(3,npt) ! [input] coordinates for points
28      double precision elfval(*)     ! [output] efg values for each coord
29c
30      integer ishell, jshell, ijshell, nshell, nbf_max, me, nproc
31      integer ilo, ihi, jlo, jhi, idim, jdim, nint
32      integer l_dens, k_dens, l_scr, k_scr, l_buf, k_buf
33      integer maxbuf, maxscr, i
34c
35      me = ga_nodeid()
36      nproc = ga_nnodes()
37c
38c     ----- calculate buffer and scratch space -----
39c           buffer = (lmax*(lmax+1)/2)^2 * (max(nder*3,1) * ngridpoints
40c           scratch = see hnd_elfder wrapper routine
41c
42      call int_init_1eelec(maxbuf,maxscr,basis,nder,npt)
43c
44      if (.not. bas_geom(basis, geom)) call errquit
45     $   ('hnd_elfcon: bad basis', 555, BASIS_ERR)
46      if (.not. bas_numcont(basis, nshell)) call errquit
47     $   ('hnd_elfcon: bas_numcont failed for basis', basis, BASIS_ERR)
48      if (.not. bas_nbf_cn_max(basis,nbf_max)) call errquit
49     &   ('hnd_elfcon: bas_nbf_cn_max failed',555, BASIS_ERR)
50c
51      if (.not. ma_push_get(mt_dbl,nbf_max*nbf_max,'dens patch',l_dens,
52     &    k_dens)) call errquit('hnd_elfcon: ma 1 failed',911,MA_ERR)
53      if (.not. ma_push_get(mt_dbl,maxscr,'scratch',l_scr,k_scr))
54     &    call errquit('hnd_elfcon: ma 2 failed',911,MA_ERR)
55      if (.not. ma_push_get(mt_dbl,maxbuf,'int buf',l_buf,k_buf))
56     &    call errquit('hnd_elfcon: ma 3 failed',911,MA_ERR)
57c
58c     Zero elfval result array
59c
60      call dcopy(max(nder*3,1)*npt,0.0d0,0,elfval,1)
61c
62c     ----- calculate electronic integral component(s) at all points -----
63c
64
65      ijshell = 0
66      do ishell = 1, nshell
67c
68c     get basis info
69c
70         if (.not. bas_cn2bfr(basis, ishell, ilo, ihi)) call errquit
71     &      ('hnd_elfcon: bas_cn2bfr failed for basis',basis,BASIS_ERR)
72         idim = ihi - ilo + 1
73
74         do jshell = 1, nshell
75            ijshell = ijshell + 1
76            if (mod(ijshell,nproc) .eq. me) then
77c
78c     get basis info
79c
80               if (.not. bas_cn2bfr(basis, jshell, jlo, jhi)) call
81     &            errquit('hnd_elfcon: bas_cn2bfr',basis,BASIS_ERR)
82               jdim = jhi - jlo + 1
83               nint = idim * jdim
84c
85c     Get the density patch, make the integrals and contract
86c
87               call ga_get(g_dens, ilo, ihi, jlo, jhi,
88     $                     dbl_mb(k_dens), idim)
89c
90               call dim_1eelec(basis,ishell,basis,jshell,maxscr,
91     &                         dbl_mb(k_scr),nint,dbl_mb(k_buf),
92     &                         nder,points,npt)
93c
94               call dim_multi_reduce(dbl_mb(k_buf),dbl_mb(k_dens),
95     &                           elfval,idim,jdim,npt*(max(nder*3,1)))
96            end if  ! mod parallel loop
97         end do   ! jshell
98      end do    ! ishell
99c
100c     Collect components from all the nodes for all points
101c
102      call ga_sync()
103      call ga_dgop(msg_efgs_col,elfval,npt*(max(nder*3,1)),'+')
104c
105c     Clean up MA data blocks
106c
107      if (.not.ma_pop_stack(l_buf)) call errquit
108     &   ('hnd_elfcon, ma_pop_stack of l_buf failed',911,MA_ERR)
109      if (.not.ma_pop_stack(l_scr)) call errquit
110     &   ('hnd_elfcon, ma_pop_stack of l_scr failed',911,MA_ERR)
111      if (.not.ma_pop_stack(l_dens)) call errquit
112     &   ('hnd_elfcon, ma_pop_stack of l_dens failed',911,MA_ERR)
113      return
114      end
115c
116      subroutine dim_multi_reduce(mblock,block,rblock,idim,jdim,nblock)
117c
118      implicit none
119#include "stdio.fh"
120      integer idim,jdim,nblock
121      double precision mblock(idim,jdim,nblock), block(idim,jdim)
122      double precision rblock(nblock)
123c
124      integer iblock,i,j
125c
126c      write(luout,*) "integral patch:"
127c      write(luout,*) mblock
128      do iblock = 1, nblock
129         do i = 1, idim
130            do j = 1, jdim
131               rblock(iblock)=rblock(iblock)+mblock(i,j,iblock)*
132     &                                       block(i,j)
133            enddo
134         enddo
135      enddo
136c
137      return
138      end
139