subroutine chrg1up(nqgrd,chrgv2,qval,iqpos,fpoh,
     &            sixeps,difeps,debfct,gchrgp,gchrg,icount1a,icount1b)
c	Mike Gilson's fancy charge distribution embedded into qdiffx.
c	A combination of Kim Sharp's chrgit1 and Anthony Nicholls' chrgup.
c       Grafting by Richard Friedman.
c	Latest version:6/2/89. First Version:4/25/89.
	include 'qdiffpar4.h'
c 
	logical isum(3*ngrid),flag
	dimension gchrg(ngcrg),gchrgd(ngcrg),chrgv2(4,ncrgmx)
	dimension qval(ngcrg),iqpos(ngcrg)
	dimension ddist2(500), rsave(500)
	dimension clist(4,50,500), nlist(500)
        dimension kb(3),dg(3)      
c	integer*2 gchrgp(3,ngcrg)
	integer gchrgp(3,ngcrg)
	integer*4 gchrg2(ngcrg)
	parameter (crit = .0001)
	data radmin,radmax,radstep / 1.2, 2.0, .005 /
c
c assign fractional charges to grid points, using phimap
c as a dummy array (set to zero again later)
        do 9010 ig=1,nqgrd
        do 9001 i = 1,3
c
c trucate to nearest grid point
c
          kb(i) = chrgv2(i,ig)
c
c find position of charge in box in terms
c of fractional distance along edges
c
          dg(i) = chrgv2(i,ig) - kb(i)
c        end do
9001	  continue
c
c loop over increasing radius to find suitable one
c
	  irmax = 0
	  do 9002 rmax = radmin, radmax, radstep
	    rmax2 = rmax*rmax
	    irmax = irmax + 1
	    rsave (irmax) = rmax
	    ilist = 0
	    ir = rmax + dg(1) + 1
c
c loop over all grid points within range of rmax:
c place base point of grid box containing the charge, at the origin, for now.
c
	    do 9003 i = -ir, ir
	      xdist = i - dg(1)
	      do 9004 j = -ir, ir
	        ydist = j - dg(2)
	        do 9005 k = -ir, ir
	          zdist = k - dg(3)
c
c calculate distance from current grid point to charge:
c
	          dist2 =  xdist*xdist + ydist*ydist + zdist*zdist
	          dist = sqrt(dist2)
c
c if grid point is closer than R, store the point and its distance:
c
		    if (dist2.le. Rmax2) then
	     		ilist = ilist + 1
	     		clist(1,ilist, irmax)  = i
	     		clist(2,ilist, irmax) = j
	     		clist(3,ilist, irmax) = k
	     		clist(4,ilist, irmax) = dist
	          endif
c	        end do
c	      end do
c	    end do
9005		  continue
9004        continue
9003      continue
	    nlist(irmax) = ilist
c
c generate weighting factors for this rmax:
c sum normalizes the weighting to one:
c
	    sum = 0
	    do 9006 il = 1, nlist(irmax)
	      clist(4,il, irmax) = Rmax - clist(4,il, irmax)
	      sum = sum + clist(4,il, irmax)
c	    end do
9006	    continue
c
c normalize the weighting factors to sum to one:
c
	    do 9007 il = 1, nlist(irmax)
	      clist(4,il, irmax) = clist(4,il, irmax)/sum
c	    end do
9007	    continue
c
c calculate center of charge for this rmax:
c
	    xsum = 0
	    ysum = 0
	    zsum = 0
	    do 9008 il = 1, nlist(irmax)
	      xsum = xsum + clist(4,il, irmax) * clist(1,il, irmax)
	      ysum = ysum + clist(4,il, irmax) * clist(2,il, irmax)
	      zsum = zsum + clist(4,il, irmax) * clist(3,il, irmax)
c	    end do
9008	    continue
c	  print *, 'center of charge in map:' , xsum,ysum, zsum
c	  print *, 'actual location:        ', dg(1), dg(2), dg(3)
c
c check whether criterion is satisfied, and if so, exit rmax loop.
c
	    ddx = dg(1) - xsum
	    ddy = dg(2) - ysum
	    ddz = dg(3) - zsum
	    ddist2(irmax) = ddx * ddx + ddy*ddy + ddz * ddz
	    if (ddist2(irmax) .le.crit) goto 1000
c
c otherwise, try another cutoff radius:		     
c
c	  end do
9002    continue
c
c if loop gets finished without a radius yielding good enough
c results, print warning,and use the best cutoff radius found:
c	print *, 'Criterion not satisfied for charge #, ', iq
c	print *, ' Using best cutoff radius found...'
	  call minimum (ddist2, irmax, ichoice, dmin)
	  dmin = sqrt (dmin)
	  rmax = rsave(ichoice)
	  goto 1500

1000	  continue
	  ichoice = irmax
	  dmin = sqrt(ddist2(ichoice))

1500	  continue

c
c now we know what set of grids we're distributing the charge over
c (clist(1-3, 1-ilist(ichoice), ichoice) for ichoice), and the weighting
c (clist(4,1-ilist, ichoice)...
c now, distribute the charge:
c
	  do 9009 ilist = 1, nlist(ichoice)
c
c get grid point by adding offset to it
c
	    kx = kb(1) + clist(1,ilist,ichoice)
	    ky = kb(2) + clist(2,ilist,ichoice)
	    kz = kb(3) + clist(3,ilist,ichoice)
c
c make sure grid point is within the big box (should be a problem only
c in cases where box edge cuts through
c or very near the protein):
c
	    flag = .false.
   	    if (kx.lt.1 ) then
	      kx = 1
	      flag = .true.
	    endif
   	    if (ky .lt.1) then
	      ky = 1
	      flag = .true.
	    endif
   	    if (kz .lt.1) then
	      kz = 1
	      flag = .true.
	    endif
   	    if (kx.gt.ngrid) then
	      kx = ngrid
	      flag = .true.
	    endif
   	    if (ky.gt.ngrid)then
	      ky = ngrid
	      flag = .true.
	    endif
   	    if (kz.gt.ngrid) then
	      kz = ngrid
	      flag = .true.
	    endif
  
	    if(flag)then
	      write(6,*)' problem for charge at', (chrgv2(j,ig),j=1,3)
	    end if

          phimap(kx,ky,kz)=phimap(kx,ky,kz)+clist(4,ilist,ichoice)*chrgv2(4,ig)

c	  end do
9009	  continue

 9010	continue
c set up odd/even logical array 
	do 10 i=1,3*ngrid,2
10	isum(i)=.true.
	do 20 i=2,(3*ngrid-1),2
20	isum(i)=.false. 
c
c find which grid points have charge assigned to them
c (will use this array later to calculate grid energy)
c 
	n=0
	do 100 k=2,igrid-1
	  do 110 j=2,igrid-1
	    do 120 i=2,igrid-1
		if(phimap(i,j,k).ne.0) then
		n=n+1
		gchrgp(1,n)=i
		gchrgp(2,n)=j
		gchrgp(3,n)=k
		gchrg(n)=phimap(i,j,k)
		phimap(i,j,k)=0.0
		end if
120	 continue
110	 continue
100	 continue 
	icount1b=n
c
c determine how many charged grid points are odd
c 
	icount1a=0
	do 200 i=1,n
	itemp=gchrgp(1,i)+gchrgp(2,i)+gchrgp(3,i)
	if(isum(itemp)) icount1a=icount1a+1
200	continue 
c
c set up odd/even pointer array, to be used in making qval
c and iqpos
c 
	i1=0
	i2=icount1a
	do 300 i=1,n
	itemp=gchrgp(1,i)+gchrgp(2,i)+gchrgp(3,i)
	if(isum(itemp)) then
	i1=i1+1
	gchrg2(i)=i1
	else
	i2=i2+1
	gchrg2(i)=i2
	end if
300	continue 
c
c determine denominator at all charged grid points
c 
	ib=0
	epsins6=sixeps+ 6.0*difeps 
	do 400 i=1,n
	iz=gchrgp(3,i)
	iy=gchrgp(2,i)
	ix=gchrgp(1,i)
	itemp1=iepsmp(ix,iy,iz,1)+iepsmp(ix-1,iy,iz,1)
	itemp2=iepsmp(ix,iy,iz,2)+iepsmp(ix,iy-1,iz,2)
	itemp3=iepsmp(ix,iy,iz,3)+iepsmp(ix,iy,iz-1,3)
	itemp=itemp1+itemp2+itemp3
	gchrgd(i)=itemp*difeps + debfct*idebmap(ix,iy,iz) + sixeps 
	if(itemp.ne.6) then
	ib=ib+1
	cgbp(1,ib)=ix 
	cgbp(2,ib)=iy
	cgbp(3,ib)=iz 
	cgbp(4,ib)=gchrg(i)*fpoh/(epsins6+debfct*idebmap(ix,iy,iz))
	cgbp(5,ib)=gchrg2(i)
	end if 
400	continue 
	ibc=ib 
	write(6,*) '# grid points charged and at boundary=',ib
c
c make qval, fpoh term so potentials will be in kt/e
c 
	do 500 i=1,n
	j=gchrg2(i)
	qval(j)=gchrg(i)*fpoh/gchrgd(i)
500	continue 
c
c make iqpos
c 
	isgrid=igrid**2
	do 600 i=1,n
	j=gchrg2(i)
	ix=gchrgp(1,i)
	iy=gchrgp(2,i)
	iz=gchrgp(3,i)
	iw=1+ix+igrid*(iy-1)+isgrid*(iz-1)
	iv=iw/2
	iqpos(j)=iv
600	continue 
c
c end of chrgup, return with qval,iqpos and gchrgp and gchrg
c also icount1a, icount1b 
	return
	end