      subroutine mespgen(nreal,ninteger,iptsflag,natoms,nesppts,
     &                   ncontacts,nauxiliary,level,ncenters,nbsystems,
     &                   nequivnb,nbneutrals,nalpha,nbeta,nconts,ncontp,
     &                   ncontd,ncds,ncdspd,nshels,nshelp,nsheld,
     &                   iprimvec,nprimvec,icfunc,mtloca,icdlfunc,
     &                   icdcfunc,natomtype,nequivvec,neutral,isys,
     &                   iswitch,nptsvec,iwkvec,coord,charge,dmata,
     &                   dmatb,alpha,coeff,alphacd,coeffscd,coeffpcd,
     &                   coeffdcd,cdfitc,potpt,qqa,qqb,radius,
     &                   cneighbour,rneighbour,contacts,wkvec)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

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

      dimension scale(4)

      dimension iprimvec(*),nprimvec(*),icfunc(*),mtloca(*),nshels(*),
     &          nshelp(*),nsheld(*),icdlfunc(*),icdcfunc(*),
     &          natomtype(*),iswitch(*),nptsvec(*),nequivvec(natoms,*),
     &          neutral(natoms,*),isys(natoms,*),iwkvec(*)

      dimension coord(3,*),charge(*),dmata(*),dmatb(*),alpha(*),
     &          coeff(*),alphacd(*),coeffscd(*),coeffpcd(*),coeffdcd(*),
     &          cdfitc(*),potpt(nesppts,*),qqa(2*natoms,*),qqb(*),
     &          radius(*),cneighbour(3,*),rneighbour(*),contacts(3,*),
     &          wkvec(*)

      data zero,bohr,one,debye,density/0.0d0,0.529177d0,
     &                                 1.0d0,2.54176568d0,1.0d0/
      data scale/1.40d0,1.60d0,1.80d0,2.00d0/

      rewind 66

      ncntrt=nconts+3*ncontp+6*ncontd
      ndim=ncntrt*(ncntrt+1)/2

      ncdfunc=ncds+10*ncdspd

      npoints=ncenters/level

      nespts=0

      dnetcharge=zero

      do 1001 i=1,ncenters
 1001 dnetcharge=dnetcharge+charge(i)

      dnetcharge=dnetcharge-real64(nalpha+nbeta)

      if(iptsflag.eq.0) then

                  nespts=0

                  rewind 90

 1002             nespts=nespts+1

                  read(90,1003,end=1004) wkvec(0*250000+nespts),
     &                                   wkvec(1*250000+nespts),
     &                                   wkvec(2*250000+nespts)
 1003             format(3f20.12)

                  wkvec(0*250000+nespts)=wkvec(0*250000+nespts)/bohr
                  wkvec(1*250000+nespts)=wkvec(1*250000+nespts)/bohr
                  wkvec(2*250000+nespts)=wkvec(2*250000+nespts)/bohr

                  goto 1002

 1004             continue

                  nespts=nespts-1

                  write(66) nespts

                  call fastwr(66,wkvec(000001),nespts)
                  call fastwr(66,wkvec(250001),nespts)
                  call fastwr(66,wkvec(500001),nespts)

                        else

                  call surfac(density,scale(1),coord,radius,cneighbour,
     &                        rneighbour,contacts,potpt,natomtype,
     &                        nptsvec,iswitch,ncontacts,ncenters,
     &                        nesppts,nesp1,nespts)

                  nespts=nespts+nesp1

                  call surfac(density,scale(2),coord,radius,cneighbour,
     &                        rneighbour,contacts,potpt,natomtype,
     &                        nptsvec,iswitch,ncontacts,ncenters,
     &                        nesppts,nesp2,nespts)

                  nespts=nespts+nesp2

                  call surfac(density,scale(3),coord,radius,cneighbour,
     &                        rneighbour,contacts,potpt,natomtype,
     &                        nptsvec,iswitch,ncontacts,ncenters,
     &                        nesppts,nesp3,nespts)

                  nespts=nespts+nesp3

                  call surfac(density,scale(4),coord,radius,cneighbour,
     &                        rneighbour,contacts,potpt,natomtype,
     &                        nptsvec,iswitch,ncontacts,ncenters,
     &                        nesppts,nesp4,nespts)

                  nespts=nespts+nesp4

                        endif

      rewind 66

      do 1005 i=1,ndim
 1005 wkvec(i)=dmata(i)+dmatb(i)

      do 1006 i=1,npoints
      qqa(i,npoints+1)=one
 1006 qqa(npoints+1,i)=one

      qqa(npoints+1,npoints+1)=zero

      if(nequivnb.gt.0) then
                              do 1007 i=1,nequivnb
                              qqa(nequivvec(i,1),npoints+1+i)=+one
                              qqa(nequivvec(i,2),npoints+1+i)=-one
                              qqa(npoints+1+i,nequivvec(i,1))=+one
 1007                         qqa(npoints+1+i,nequivvec(i,2))=-one
                        endif

      if(nbneutrals.gt.0) then

                     do 1008 i=1,nbneutrals
                     do 1008 j=1,neutral(i,1)
                     qqa(neutral(i,j+1),npoints+1+nequivnb+i)=+one
 1008                qqa(npoints+1+nequivnb+i,neutral(i,j+1))=+one

                          endif

      qqb(npoints+1)=dnetcharge

      iesp1=1
      iesp2=iesp1+ndim
      iesp3=iesp2+ndim
      iesp4=iesp3+nespts*ncenters
      iesp5=iesp4+nespts*7
      iesp6=iesp5+ncdfunc

      if(iesp6.gt.nreal) stop 'mesp evaluation'

      call mespmm(nreal-iesp4,ninteger,iptsflag,nauxiliary,npoints,
     &            nespts,ncenters,natoms,nbsystems,nconts,ncontp,ncontd,
     &            ncds,ncdspd,nshels,nshelp,nsheld,mtloca,iprimvec,
     &            nprimvec,icfunc,icdcfunc,icdlfunc,nptsvec,isys,iwkvec,
     &            coord,charge,alpha,coeff,alphacd,coeffscd,coeffpcd,
     &            coeffdcd,wkvec(iesp5),wkvec(iesp1),wkvec(iesp2),qqa,
     &            qqb,wkvec(iesp3),wkvec(iesp4),wkvec(iesp6))

      call matinv(qqa,qqb,npoints+1+nequivnb+nbneutrals,1,2*natoms)

      write(6,1009) nespts
 1009 format(//,' number of points used to fit esp = ',i6,/)

      write(6,1010)
 1010 format(' esp-fitted charges using the true density :',/)

      do 1011 i=1,npoints
 1011 write(6,1012) i,coord(1,i)*bohr,coord(2,i)*bohr,
     &                                coord(3,i)*bohr,qqb(i)
 1012 format(8x,'atom # ',i3,' at ',3f10.5,8x,' charge = ',f8.5)

      write(6,1013)
 1013 format()

      xcomp=zero
      ycomp=zero
      zcomp=zero

      do 1014 i=1,npoints
      xcomp=xcomp+coord(1,i)*qqb(i)
      ycomp=ycomp+coord(2,i)*qqb(i)
 1014 zcomp=zcomp+coord(3,i)*qqb(i)

      dipole=sqrt(xcomp*xcomp+ycomp*ycomp+zcomp*zcomp)*debye

      write(6,1015) dipole
 1015 format(8x,'dipole moment from point charges = ',f9.4,' debyes',//)

      return
      end
