      subroutine gridgen(ngridtype,nscf,ncenters,nconts,ncontp,ncontd,
     &                   ilfunc,icfunc,ngauss,nptrs,nptrp,nptrd,
     &                   natomtype,nptsatom,iwkvec,cent,alpha,coeff,
     &                   xcoord,ycoord,zcoord,weight,dcp,dcc,wtbecke,
     &                   extent)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      include "mpif.h"

      dimension istat(mpi_status_size)

      dimension slater(105),gcroot(96),gcwght(96)

      dimension coordshell(3,194),wghtshell(194)

      dimension ilfunc(*),icfunc(*),ngauss(*),nptrs(*),nptrp(*),nptrd(*)
      dimension natomtype(*),nptsatom(*),iwkvec(*)

      dimension cent(3,*),alpha(*),coeff(*)

      dimension xcoord(*),ycoord(*),zcoord(*),weight(*)
      dimension dcp(ncenters,*),dcc(ncenters,*),wtbecke(*),extent(*)

      data slater/ 0.35d0, 2.00d0, 1.45d0, 1.05d0, 0.85d0, 0.70d0,
     &             0.65d0, 0.60d0, 0.50d0, 2.25d0, 1.80d0, 1.50d0,
     &             1.25d0, 1.10d0, 1.00d0, 1.00d0, 1.00d0, 2.50d0,
     &             2.20d0, 1.80d0, 1.60d0, 1.40d0, 1.35d0, 1.40d0,
     &             1.40d0, 1.40d0, 1.35d0, 1.35d0, 1.35d0, 1.35d0,
     &             1.30d0, 1.25d0, 1.15d0, 1.15d0, 1.15d0, 2.75d0,
     &             2.35d0, 2.00d0, 1.80d0, 1.55d0, 1.45d0, 1.45d0,
     &             1.35d0, 1.30d0, 1.35d0, 1.40d0, 1.60d0, 1.55d0,
     &             1.55d0, 1.45d0, 1.45d0, 1.40d0, 1.40d0, 3.00d0,
     &             2.60d0, 2.15d0, 1.95d0, 1.85d0, 1.85d0, 1.85d0,
     &             1.85d0, 1.85d0, 1.85d0, 1.80d0, 1.75d0, 1.75d0,
     &             1.75d0, 1.75d0, 1.75d0, 1.75d0, 1.75d0, 1.55d0,
     &             1.45d0, 1.35d0, 1.35d0, 1.30d0, 1.35d0, 1.35d0,
     &             1.35d0, 1.50d0, 1.90d0, 1.80d0, 1.60d0, 1.90d0,
     &             1.65d0, 3.25d0, 2.80d0, 2.15d0, 1.95d0, 1.80d0,
     &             1.80d0, 1.75d0, 1.75d0, 1.75d0, 1.75d0, 1.75d0,
     &             1.75d0, 1.75d0, 1.75d0, 1.75d0, 1.75d0, 1.75d0,
     &             1.75d0, 1.55d0, 1.55d0/

      data zero,half,one,two,four,five/0.0d0,0.5d0,1.0d0,
     &                                 2.0d0,4.0d0,5.0d0/
      data sixteen,twentyone,thirtyfive,one_o_5/16.0d0,21.0d0,
     &                                          35.0d0,105.0d0/
      data tolerance,bohr/0.0000000001d0,0.529177d0/
      data afactor,pi/0.64d0,3.1415926535898d0/
      data alpha1,alpha2,alpha3/0.40d0,1.20d0,2.00d0/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      nptstotal=0

      if(myid.eq.0) rewind 11

      if(nscf.eq.1) then
                          ncheck=mod(ngridtype,100)
                          if(ncheck.ge.10) nradpts=24
                          if(ncheck.ge.20) nradpts=32
                          if(ncheck.ge.30) nradpts=40
                          if(ncheck.ge.40) nradpts=48
                          if(ncheck.ge.50) nradpts=64
                          if(ncheck.ge.60) nradpts=80
                          if(ncheck.ge.70) nradpts=96
                          ncheck=mod(ngridtype,10)
                          if(ncheck.ge.1) nnangpts=026
                          if(ncheck.ge.2) nnangpts=050
                          if(ncheck.ge.3) nnangpts=110
                          if(ncheck.ge.4) nnangpts=194
                    else
                          ncheck=ngridtype
                          if(ncheck.ge.1000) nradpts=24
                          if(ncheck.ge.2000) nradpts=32
                          if(ncheck.ge.3000) nradpts=40
                          if(ncheck.ge.4000) nradpts=48
                          if(ncheck.ge.5000) nradpts=64
                          if(ncheck.ge.6000) nradpts=80
                          if(ncheck.ge.7000) nradpts=96
                          ncheck=mod(ngridtype,1000)
                          if(ncheck.ge.100) nnangpts=026
                          if(ncheck.ge.200) nnangpts=050
                          if(ncheck.ge.300) nnangpts=110
                          if(ncheck.ge.400) nnangpts=194
                    endif

      call diffuse1(ncenters,nconts,ncontp,ncontd,ilfunc,icfunc,ngauss,
     &              nptrs,nptrp,nptrd,tolerance,alpha,coeff,cent,extent)

      do 1001 jatom=1,ncenters
      do 1001 iatom=1,ncenters
 1001 dcc(iatom,jatom)=dist(cent(1,iatom),cent(1,jatom))

      if(myid.eq.0) then
                          isent=0

                          do 1002 islave=0,numprocs-1

                          ia=ncenters-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          call mpi_bsend(isent,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

 1002                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)

      call mpi_recv(isent,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)

      do 1003 iatom=isent+1,isent+nloop

      if(natomtype(iatom).eq.0) goto 1003

      smallest=one/tolerance

      do 1004 jatom=1,ncenters
      if((jatom.ne.iatom).and.(natomtype(jatom).ne.0)) then

              if(dcc(iatom,jatom).le.smallest) smallest=dcc(iatom,jatom)

                                                       endif
 1004 continue

      smallest=smallest*(one-afactor)*half

      npoints=0

      do 1005 nshell=1,nradpts

      gcroot(nshell)=-cos(real64(nshell)*pi/real64(nradpts+1))

      gcwght(nshell)=(pi/real64(nradpts+1))
      gcwght(nshell)=gcwght(nshell)
     &              *(sin(real64(nshell)*pi/real64(nradpts+1)))**2
      gcwght(nshell)=gcwght(nshell)/sqrt(one-gcroot(nshell)**2)

      halfslater=half*slater(natomtype(iatom))/bohr

      rfactor=(one+gcroot(nshell))/(one-gcroot(nshell))

      radius=halfslater*rfactor

      thickness=two*halfslater/((one-gcroot(nshell))**2)

      volume=four*pi*radius*radius*thickness*gcwght(nshell)

      nangpts=nnangpts

      if(rfactor.le.alpha3) then
                                  if(nnangpts.eq.194) nangpts=110
                                  if(nnangpts.eq.110) nangpts=050
                                  if(nnangpts.eq.050) nangpts=026
                                  if(nnangpts.eq.026) nangpts=012
                            endif

      if(rfactor.le.alpha2) then
                                  if(nnangpts.eq.194) nangpts=050
                                  if(nnangpts.eq.110) nangpts=026
                                  if(nnangpts.eq.050) nangpts=012
                                  if(nnangpts.eq.026) nangpts=012
                            endif

      if(rfactor.le.alpha1) then
                                  if(nnangpts.eq.194) nangpts=26
                                  if(nnangpts.eq.110) nangpts=12
                                  if(nnangpts.eq.050) nangpts=12
                                  if(nnangpts.eq.026) nangpts=12
                            endif

      call points(nangpts,radius,volume,cent(1,iatom),coordshell,
     &            wghtshell)

      do 1006 ipoint=1,nangpts
 1006 dcp(iatom,ipoint)=radius

      do 1005 ipoint=1,nangpts

      do 1007 icenter=1,ncenters
 1007 wtbecke(icenter)=zero

      if(radius.lt.smallest) then
                                   ihits=1
                                   iwkvec(1)=iatom
                                   wtbecke(iatom)=one
                                   goto 1008
                             endif

      ihits=0

      do 1009 icenter=1,ncenters
      rr=(cent(1,icenter)-coordshell(1,ipoint))**2
     &  +(cent(2,icenter)-coordshell(2,ipoint))**2
     &  +(cent(3,icenter)-coordshell(3,ipoint))**2
      if(rr.le.extent(icenter)) then
                                      ihits=ihits+1
                                      iwkvec(ihits)=icenter
                                endif
 1009 continue

      if(ihits.eq.0) goto 1005

      do 1010 icenter=1,ihits
      jcenter=iwkvec(icenter)
 1010 if(jcenter.ne.iatom) dcp(jcenter,ipoint)=dist(cent(1,jcenter),
     &                                         coordshell(1,ipoint))

      do 1011 icenter=1,ihits

      jcenter=iwkvec(icenter)

      if(natomtype(jcenter).eq.0) goto 1011 

      partition=one

      do 1012 kcenter=1,ihits

      lcenter=iwkvec(kcenter)

      if((lcenter.eq.jcenter).or.(natomtype(lcenter).eq.0)) goto 1012

      u=(dcp(jcenter,ipoint)-dcp(lcenter,ipoint))/dcc(jcenter,lcenter)

      if(u.le.(-afactor)) z=-one

      if(abs(u).lt.afactor) z=(thirtyfive*(u/afactor)
     &                        -thirtyfive*(u/afactor)**3
     &                        +twentyone*(u/afactor)**5
     &                        -five*(u/afactor)**7)/sixteen

      if(u.ge.(+afactor)) z=+one

      s=half*(one-z)

      partition=partition*s

 1012 continue

      wtbecke(jcenter)=partition

 1011 continue

 1008 continue

      totalwtbecke=zero

      do 1013 kcenter=1,ihits
      lcenter=iwkvec(kcenter)
 1013 totalwtbecke=totalwtbecke+wtbecke(lcenter)

      wghtshell(ipoint)=wghtshell(ipoint)*wtbecke(iatom)/totalwtbecke

      npoints=npoints+1

      xcoord(npoints)=coordshell(1,ipoint)
      ycoord(npoints)=coordshell(2,ipoint)
      zcoord(npoints)=coordshell(3,ipoint)

      weight(npoints)=wghtshell(ipoint)

 1005 continue

      if(numprocs.eq.1) then

                        nptsatom(iatom)=npoints

                        call fastwr(11,xcoord,nptsatom(iatom))
                        call fastwr(11,ycoord,nptsatom(iatom))
                        call fastwr(11,zcoord,nptsatom(iatom))
                        call fastwr(11,weight,nptsatom(iatom))

                        nptstotal=nptstotal+npoints

                        else

                        call mpi_bsend(npoints,1,mpi_integer,0,
     &                                 5*iatom+1,mpi_comm_world,ierr)

                        call mpi_bsend(xcoord,npoints,mpi_real8,0,
     &                                 5*iatom+2,mpi_comm_world,ierr)
                        call mpi_bsend(ycoord,npoints,mpi_real8,0,
     &                                 5*iatom+3,mpi_comm_world,ierr)
                        call mpi_bsend(zcoord,npoints,mpi_real8,0,
     &                                 5*iatom+4,mpi_comm_world,ierr)
                        call mpi_bsend(weight,npoints,mpi_real8,0,
     &                                 5*iatom+5,mpi_comm_world,ierr)

                        endif

 1003 continue

      if((myid.eq.0).and.(numprocs.gt.1)) then

               do 1014 iatom=1,ncenters

               call mpi_recv(npoints,1,mpi_integer,mpi_any_source,
     &                       5*iatom+1,mpi_comm_world,istat,ierr)

               call mpi_recv(xcoord,npoints,mpi_real8,mpi_any_source,
     &                       5*iatom+2,mpi_comm_world,istat,ierr)
               call mpi_recv(ycoord,npoints,mpi_real8,mpi_any_source,
     &                       5*iatom+3,mpi_comm_world,istat,ierr)
               call mpi_recv(zcoord,npoints,mpi_real8,mpi_any_source,
     &                       5*iatom+4,mpi_comm_world,istat,ierr)
               call mpi_recv(weight,npoints,mpi_real8,mpi_any_source,
     &                       5*iatom+5,mpi_comm_world,istat,ierr)

               nptsatom(iatom)=npoints

               call fastwr(11,xcoord,nptsatom(iatom))
               call fastwr(11,ycoord,nptsatom(iatom))
               call fastwr(11,zcoord,nptsatom(iatom))
               call fastwr(11,weight,nptsatom(iatom))

 1014          nptstotal=nptstotal+npoints

                                          endif

      if(myid.eq.0) write(6,1015) nptstotal
 1015 format(/,' total number of grid points = ',2i9,/)

      return
      end
