subroutine react(ibgrd,nqass,icount2b,epsval,ergs,rad3,xn2)
c
	include 'qdiffpar4.h'
c
	dimension schrg(nsp),sqs(nsp),ibgrd(4,nsp),rbgrd(3,nsp)
	dimension rad3(natmax),xn2(3,natmax),mv(6),sqa(natmax)
	dimension cgrid(ngrid,ngrid,ngrid),xo(3),xn(3),sen(nsp)
	fact=0.9549296586/2.0
	sixth=1.0/6.0
	en1=0.0
c
c calculate surface charge
c
	do 10 i=1,icount2b
 	  ix=ibgrd(1,i)
	  iy=ibgrd(2,i)
	  iz=ibgrd(3,i)
	  temp1=phimap(ix-1,iy,iz)+phimap(ix+1,iy,iz)
	  temp2=phimap(ix,iy-1,iz)+phimap(ix,iy+1,iz)
  	  temp3=phimap(ix,iy,iz-1)+phimap(ix,iy,iz+1)
	  schrg(i)=phimap(ix,iy,iz)-(temp1+temp2+temp3)*sixth 
	  schrg(i)=schrg(i)*fact
	  en1=en1+(schrg(i))
10	continue 
c
c	write(6,*) "total s.charge= ", en1/scale
c
c remove fixed surface charge
c
	do i=1,icount2b
	ix=ibgrd(1,i)
	iy=ibgrd(2,i)
	iz=ibgrd(3,i)
	cgrid(ix,iy,iz)=schrg(i)
	end do
c
	do i=1,ibc
	ix=cgbp(1,i)
	iy=cgbp(2,i)
	iz=cgbp(3,i)
	cgrid(ix,iy,iz)=cgrid(ix,iy,iz)-cgbp(4,i)*fact
	end do
c
	do i=1,icount2b
	ix=ibgrd(1,i)
	iy=ibgrd(2,i)
	iz=ibgrd(3,i)
	schrg(i)=cgrid(ix,iy,iz)
	end do
c
c NB real charge is schrg/scale
c 
c generate pseudo distances for surface points
c
	ivz=0
	imt=0
	do 20 i=1,icount2b
c
	ib1=ibgrd(1,i)
	ib2=ibgrd(2,i)
	ib3=ibgrd(3,i)
	rx=float(ib1)
	ry=float(ib2)
	rz=float(ib3)
c
	mv(1)=iepsmp2(ib1,ib2,ib3,1)
	mv(2)=iepsmp2(ib1,ib2,ib3,2)
	mv(3)=iepsmp2(ib1,ib2,ib3,3)
	mv(4)=iepsmp2(ib1-1,ib2,ib3,1)
	mv(5)=iepsmp2(ib1,ib2-1,ib3,2)
	mv(6)=iepsmp2(ib1,ib2,ib3-1,3)
c
	do j=1,6
	if(mv(j).lt.0) then
	if(mv(j).eq.-1) rz=rmmax
	if(mv(j).eq.-2) rz=rmmin
	imt=imt+1
	goto 44
	end if
	end do
c
	dism=1000.
	iat=100000
	do j=1,6
	m=mv(j)
	if((m.ne.0).and.(m.ne.iat)) then
	dist=(xn2(1,m)-rx)**2 + (xn2(2,m)-ry)**2 + (xn2(3,m)-rz)**2
	if(dist.lt.dism) then
	iat=m
	dism=dist
	end if
	end if
	end do
	if(iat.eq.100000) then
	ivz=ivz+1
	iv=0
	else
	iv=iat
	end if
c
	if(iv.ne.0) then
	rad=rad3(iv)*scale
	arad=(xn2(1,iv)-rx)**2 + (xn2(2,iv)-ry)**2 + (xn2(3,iv)-rz)**2
	if(arad.eq.0.) then
	write(6,*) "oops!!"
	else
	sfact=rad/sqrt(arad)
	end if
	rx=xn2(1,iv) + sfact*(rx-xn2(1,iv))
	ry=xn2(2,iv) + sfact*(ry-xn2(2,iv))
	rz=xn2(3,iv) + sfact*(rz-xn2(3,iv))
	end if
c
c
44	rbgrd(1,i)=rx
	rbgrd(2,i)=ry
	rbgrd(3,i)=rz
	sqs(i)=(rx**2 + ry**2 + rz**2)/2.0
20	continue 
c
	write(6,*) "number of unassigned boundary points= ",ivz
	write(6,*) "number of membrane points scaled= ",imt
c
	en=0.0
	en2=0.0
	en3=0.0
	sfact=2.0/(scale*scale)
c
c calculate back energies
c
	do 40 i=1,nqass 
	dist1=atmcrg(1,i)
	dist2=atmcrg(2,i) 
	dist3=atmcrg(3,i) 
	dist4=(dist1**2 + dist2**2 + dist3**2)/2.0
	chrg=atmcrg(4,i)/sqrt(2.0)
	en1=0.0
c
	do 50 j=1,icount2b
	prod=dist1*rbgrd(1,j)+dist2*rbgrd(2,j)+dist3*rbgrd(3,j)
	dist=dist4+sqs(j)-prod 
	temp=schrg(j)/sqrt(dist)
	en1=en1+temp
	sen(j)=sen(j) + chrg*temp
50	continue
c
	en=en+en1*chrg
40	continue 
	ergs=en/2.0
c
	if(isch .ge. 1) then
	write(6,*) "writing surface charge file: surfch.dat"
	open(41,file="surfch.dat")
	do i=1,icount2b
	xo(1)=rbgrd(1,i)
	xo(2)=rbgrd(2,i) 
	xo(3)=rbgrd(3,i) 
	call gtoc(xo,xn)
	en1=schrg(i)
	write(41,414) i,xn,en1
	end do
	close(41)
	end if
c
c
	if(isen .ge. 1) then
	write(6,*) "writing surface energy file: surfen.dat"
	open(41,file="surfen.dat")
	en2=0.0
	do i=1,icount2b
	xo(1)=rbgrd(1,i)
	xo(2)=rbgrd(2,i) 
	xo(3)=rbgrd(3,i) 
	call gtoc(xo,xn)
	en1=sen(i)/2.0
	write(41,414) i,xn,en1
	end do
	close(41)
	end if
c
c if isen, calculate surface energy positions and values
c
414	format(i5,3f8.3,f8.3)
c
c
	write(6,*) 'corrected reaction field energy:      ',ergs,' kt'
	return
	end