      subroutine surfac(density,scale,c,rad,cn,rn,cpts,potc,natomtype,
     &                  nptsvec,iatomtype,nconlimit,ncentr,maxpts,npts,
     &                  nespts)

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)

c  origins of this routine:
c
c  this subroutine was lifted from michael connolly's surface program
c  for ucsf graphics systems by u. chandra singh and peter a. kollman
c  and modified for use in quest
c
c  k. m. merz adapted and modified this program for use in ampac and
c  mopac
c
c  ian r. gould further modified this code for use in gaussian 90

      dimension vandwl(100),ciatom(3),centre(3)

      dimension c(3,*),rad(*),cn(3,*),rn(*),cpts(3,*),potc(maxpts,*)

      dimension natomtype(*),iatomtype(*),nptsvec(*)

      logical collid

      data vandwl/1.20d0,1.20d0,1.37d0,1.45d0,1.45d0,1.50d0,
     &            1.50d0,1.40d0,1.35d0,1.30d0,1.57d0,1.36d0,
     &            1.24d0,1.17d0,1.80d0,1.75d0,1.70d0,
     &            17*2.5d0,2.5d0,65*2.5d0/ 

      data zero,bohr,pi,four/0.0d0,0.529177d0,3.1415926535898d0,4.0d0/

      npts=0

      do 1001 iatom=1,ncentr

      iatomtype(iatom)=0

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

      rad(iatom)=scale*vandwl(natomtype(iatom))/bohr

      if(rad(iatom).eq.zero) then

          write(6,1002) iatom
 1002     format(' van der waals radius for atom',i5,' is set to 2.5',/,
     &           ' you may want to supply another value in the',
     &           ' subroutine surfac',/)

                             endif

      iatomtype(iatom)=1

 1001 continue

      do 1003 iatom=1,ncentr

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

      riatom=rad(iatom)

      ciatom(1)=c(1,iatom)
      ciatom(2)=c(2,iatom)
      ciatom(3)=c(3,iatom)

      nnbr=0

      do 1004 jatom=1,ncentr

      if((iatom.eq.jatom).or.(iatomtype(jatom).eq.0)) goto 1004

      if(dist(ciatom,c(1,jatom)).ge.(riatom+rad(jatom))) goto 1004

      nnbr=nnbr+1

      rn(nnbr)=rad(jatom)

      cn(1,nnbr)=c(1,jatom)
      cn(2,nnbr)=c(2,jatom)
      cn(3,nnbr)=c(3,jatom)

 1004 continue

      ncon=(four*pi*(riatom**2))*density
      if(ncon.gt.nconlimit) ncon=nconlimit

      call genun(cpts,ncon,number)

      do 1005 i=1,number

      centre(1)=ciatom(1)+riatom*cpts(1,i)
      centre(2)=ciatom(2)+riatom*cpts(2,i)
      centre(3)=ciatom(3)+riatom*cpts(3,i)

      if(collid(centre,zero,cn,rn,nnbr)) goto 1005

      npts=npts+1

      if(npts.gt.maxpts) then
                               stop 'too many points in surfac'
                         else
                               potc(npts,1)=ciatom(1)+riatom*cpts(1,i)
                               potc(npts,2)=ciatom(2)+riatom*cpts(2,i)
                               potc(npts,3)=ciatom(3)+riatom*cpts(3,i)
                         endif

      nptsvec(nespts+npts)=iatom

 1005 continue

 1003 continue

      write(66) npts

      call fastwr(66,potc(1,1),npts)
      call fastwr(66,potc(1,2),npts)
      call fastwr(66,potc(1,3),npts)

      return
      end
