1      subroutine esp_esp(ene,occ,dns,x,q,grid,val)
2c
3c $Id$
4c
5      implicit none
6c
7#include "esp_params.fh"
8#include "esp_common.fh"
9#include "global.fh"
10#include "bas.fh"
11#include "rtdb.fh"
12#include "geom.fh"
13#include "mafdecls.fh"
14#include "msgids.fh"
15#include "util.fh"
16#include "schwarz.fh"
17#include "stdio.fh"
18c
19      integer ga_create_atom_blocked
20      external ga_create_atom_blocked
21c
22      real*8 ene(2*nbf),occ(2*nbf),dns(mbf,mbf)
23      real*8 x(3,natoms),q(natoms)
24      real*8 grid(3,mxgrid),val(mxgrid)
25c
26      integer i,j,ish,jsh,ix
27      integer ilo,ihi,jlo,jhi
28      real*8 dist,vt
29      integer imin,nval,jshi
30      real*8 dmin,fact
31      character*10 today,now
32c
33      integer maxbuf,maxscr,avail
34      integer batch,no_batches,ngrid_batches
35      integer iptr
36      logical ldummy(1)
37c
38c     get the density matrix, occupation and energies
39c
40      call esp_denmat(occ,ene)
41c
42      if (lscrn.ne.0) call schwarz_init(igeom,ibasis)
43c
44c     get electrostatic potential on the grid points
45c
46      do i=1,ngrid
47         val(i)=0.0d0
48      end do
49c
50c     estimate mem. req.
51c
52      call int_init_1eelec(maxbuf,maxscr,ibasis,0,ngrid)
53      avail=ma_inquire_avail(mt_dbl)*8/10
54      no_batches=nint(dble(maxbuf+maxscr)/dble(avail))
55      no_batches=max(no_batches,1)
56      ngrid_batches=ngrid/no_batches
57#ifdef DEBUG
58      if(me.eq.0) then
59      write(6,*) ' avail ',avail
60      write(6,*) ' no_batches ',no_batches
61      write(6,*) ' ngrid_batches ',ngrid_batches
62      endif
63#endif
64      iptr=1
65      do batch=1,no_batches
66         if(batch.eq.no_batches) ngrid_batches=ngrid-iptr+1
67         call hnd_elfcon_schw(ibasis,igeom,lg_d,
68     C        grid(1,iptr),ngrid_batches,val(iptr),0)
69         iptr=iptr+ngrid_batches
70      enddo
71      if(me.eq.0) then
72      do 8 i=1,ngrid
73         imin=0
74         val(i)=-val(i)
75         vt=val(i)
76         do 9 ix=1,natoms
77            dist=sqrt((grid(1,i)-x(1,ix))*(grid(1,i)-x(1,ix))+
78     +          (grid(2,i)-x(2,ix))*(grid(2,i)-x(2,ix))+
79     +          (grid(3,i)-x(3,ix))*(grid(3,i)-x(3,ix)))
80            if(imin.eq.0) then
81               dmin=dist
82               imin=ix
83            endif
84            if(dmin.gt.dist) then
85               dmin=dist
86               imin=ix
87            endif
88            val(i)=val(i)+q(ix)/dist
89    9    continue
90    8 continue
91      endif
92c
93      if(lscrn.ne.0) call schwarz_tidy()
94c
95      if(np.gt.1) then
96        call ga_brdcst(mre_002,val,ngrid*ma_sizeof(mt_dbl,1,mt_byte),0)
97      endif
98c
99c     deallocate memory for the integrals
100c
101      if(.not.ga_destroy(lg_d))
102     + call md_abort('esp_denmat: ga_destroy lg_d failed',me)
103c
104c     integral termination
105c
106      call int_terminate()
107c
108      if(me.eq.0) then
109      call util_file_name('grid',.false.,.false.,grdfil)
110      open(unit=lfngrd,file=grdfil,form='formatted',status='unknown')
111      rewind(lfngrd)
112      write(lfngrd,'(i10,f20.10)') ngrid,charge
113      do i=1,ngrid
114        write(lfngrd,'(4f20.10)') (grid(j,i),j=1,3),val(i)
115      enddo
116      close(unit=lfngrd)
117      endif
118c
119      return
120      end
121