program qdiff
c------------------------------------------------------------------
	include 'qdiffpar4.h'
c------------------------------------------------------------------
      character*8 hour
      character*9 day
      character*60 toplbl,comment
      character*10 nxtlbl
	character*9 bclab(4)
	character*80 line,filnam
	character*6 head
	character*3 tres,sres
	character*1 schn
	character*6 tatm,satm
	character*4 snum
	character*5 atnum 
	character*24 crdstr
	character*13 radstr
	character*8 testchar
	logical iper,iconc,isite,iatout,ibios,isph
	logical itobig,itest2,phiwrt,diff
c------------------------------------------------------------------
c------------------------------------------------------------------
	dimension iper(3)
	dimension cmin(3)
	dimension cmax(3)
	dimension cran(3)
	dimension xo(3)
	dimension xn(3),xn2(3,natmax),ibgrd(4,nsp),dbval(0:1,0:6,0:1)
	integer*4 gchrgp(3,ngcrg)
	logical   it1(0:6)
	dimension rad3(natmax),chrgv2(4,ncrgmx),idpos(nsp),db(6,nsp)
	dimension ismin(3),sf1(nhgp),sf2(nhgp),fc(0:1,0:1,0:1)
	dimension igi(0:1,0:1,0:1),qval(ngcrg),iqpos(ngcrg),iepsv(nsp)
	dimension ismax(3),ioff(3,12000),iceps(5),sfd(5,0:1)
	dimension sq1(-15:15),sq2(-15:15),sq3(-15:15)
	dimension rad2a1(-15:15),rad2a2(-15:15),rad2a3(-15:15)
	dimension xyzmid(3),chrgv4(natmax),fc2(0:1,0:1,0:1)
	dimension offset(3),iepsv2(nsp),gchrg(ngcrg)
	dimension cqplus(3),cqmin(3),summid(3),tempmid(3)
	dimension tomid(3,3)   
	character enc(4)
	logical	  logg,logs,logc,loga,ipdbwrt,ipdbrd,ifrcrd,ifrcwrt
	logical   igraph,ipoten,imem,iautocon
c	offset of midlines from grid points
c-----------------------------------------------------------------
c------------------------------------------------------------------
c
c conversion to kT/charge at 25 celsius
c
      data epkt,fpi /561.0,12.566/
c
c converts between ionic st (M) and debye lenght
c
	data dfact / 3.047/	
c
c boundary condition types
c
	data bclab / 'zero     ','dipolar  ',
     &             'focussing','coulombic' /
c
c defines offset to midlines along three axes
c
      data tomid /0.5,0.0,0.0, 0.0,0.5,0.0, 0.0,0.0,0.5/
c------------------------------------------------------------------
c------------------------------------------------------------------
	write(6,*)'   '
	write(6,*)'--------------------QDIFFXS------------------------'
	write(6,*)'*                                                 *'
	write(6,*)'* a program to solve the P-B equation             *'
	write(6,*)'* in 3D, using non-linear form, incorporating     *'
	write(6,*)'* 2 dielectric regions, ionic strength, periodic  *'
	write(6,*)'* and focussing boundary conditions, utilizing    *'
	write(6,*)'* stripped optimum successive over-relaxation     *'
	write(6,*)'* and surface charge position correction          *'
	write(6,*)'*                                                 *'
	write(6,*)'--------------------QDIFFXS------------------------'
	write(6,*)'  '
c------------------------------------------------------------------
	call rdline
c------------------------------------------------------------------
c
c initialization
c
	start = cputime(0.0)
	do 9000 i = 1,nrlist
	  irnumb(i) = 0
	  irlink(i) = 0
9000	continue
	do 9001 i = 1,nclist
	  icnumb(i) = 0
	  iclink(i) = 0
9001	continue
c
c
204	format(a80)
205	format(3f8.3)
c
	call rdprm(iautocon,perfil,offset,epsin,epsout,rionst,exrad
     &        ,radprb,ibctyp,iper,nlit,nnit,iconc,ibios,isite,iatout
     &        ,toplbl,diff,isph,ipdbwrt,ipdbrd,ifrcwrt,ifrcrd,ipoten
     &        ,igraph,imem,icon1,icon2,phiwrt,logg,logs,logc,loga)
c
c convert ionic strength to debye length
c
	if(rionst.ne.0.) then
	  deblen = dfact/sqrt(rionst)
	else
	  deblen = 1.e6
	end if
c
c print title, description
c
      call date(day)
c     call time(hour)
      hour='not working'
	write(6,*)'  '
      write(6,*)' program run on ',day
      write(6,*)'             at ',hour
	write(6,*)'  '
	write(6,*)'parameters read from file'
	write(6,*)filnam
	write(6,*)'grid size                  :',igrid
	if(perfil.ne.0) then
	write(6,*)'percent of box to be filled:',perfil
	else
	write(6,*)'scale,in grids/A, set to be:',scale 
	end if
	write(6,*)'object centred at (gu)     :',offset
	write(6,*)'inner,outer dielectrics    :',epsin,epsout
	write(6,*)'ionic strength (M)         :',rionst
	write(6,*)'debye length (A)           :',deblen
	write(6,*)'ion exclusion radius (A)   :',exrad
	write(6,*)'probe radius (A)           :',radprb
	write(6,*)'boundary conditions        : ',bclab(ibctyp)
	write(6,*)'x,y,z periodic bc. flags   :',iper
	if(iautocon) then
	if(gten.eq.0.) write(6,*)'# of linear iterations     : automatic'
	if(gten.ne.0.) write(6,*)'convergence by grid energy :',gten,' kt'
	else 
	write(6,*)'# of linear iterations     :',nlit
	end if 
	write(6,*)'# of non-linear iterations :',nnit
	write(6,*)'concentration map output   :',iconc
	write(6,*)'spherical charge distbn.   :',isph
	write(6,*)'INSIGHT format output      :',ibios
	write(6,*)'site potential output      :',isite
	write(6,*)'modified atom file output  :',iatout
	write(6,*)'map file label             :'
	write(6,201)toplbl
c 
	if(ipdbrd.or.ipdbwrt) then
	  if(ipdbrd) write(6,*) 'set to  read unformatted pdb file'
	if(ipdbrd.and.ipdbwrt) then
	write(6,*) ' '
	write(6,*) ' WARNING: can not write an unformatted pdb'
	write(6,*) ' file, while reading in an unformatted pdb file'
	write(6,*) ' Therefore the write option is disabled'
	else
	if(ipdbwrt) write(6,*) ' set to write unformatted pdb file' 
	end if 
	end if 
c 
	if(ifrcrd.or.ifrcwrt) then 
	if(ifrcrd.and.isite) write(6,*) 'set to read unformatted frc.pdb file'
	if(ifrcrd.and.ifrcwrt) then
	write(6,*) ' '
	write(6,*) ' WARNING: can not write an unformatted frc'
	write(6,*) ' file, while reading in an unformatted frc file'
	write(6,*) ' Therefore the write option is disabled'
	else
	if(ifrcwrt) write(6,*) ' set to write unformatted frc.pdb file'
	end if 
	end if 
c
	if(.not.igraph) write(6,*) ' convergence graph turned off'
	if(.not.ipoten) write(6,*) ' potential listings turned off'
	if((icon1.ne.10).or.(icon2.ne.1)) then
	write(6,*) 'convergence test interval is every',icon1,'loops'
	write(6,*) 'testing',100/icon2,'%'
	end if 
	write(6,*)' '
	
c
c scale 4.pi factor to grid
c
c
c scale dielectric so that potentials turn out as kT when 
c charges are in e's and distances are in angstroms
c
	epsin = epsin/epkT
	epsout = epsout/epkT
c
c set eps and deb map values to their exterior values
c prior to resetting interior values using atom coords and radii
c
	write(6,*)'   '
	write(6,*)'initializing dielectric and debye maps...'
	write(6,*)'   '
	  do 9008 k = 1,igrid
	    do 9009 j = 1,igrid
		do 9010 i = 1,igrid
		  iepsmp(i,j,k,1) = 0
		  iepsmp(i,j,k,2) = 0
		  iepsmp(i,j,k,3) = 0
		  iepsmp2(i,j,k,1) = 0
		  iepsmp2(i,j,k,2) = 0
		  iepsmp2(i,j,k,3) = 0
		  idebmap(i,j,k)  = 1
9010		continue
9009	    continue
9008	  continue
c
c jump point for shortened version
c
c
c initialize a few things
c
	qnet  = 0.0
	qplus = 0.0
	qmin  = 0.0
	do 9014 k  = 1,3
	  cqplus(k) = 0.0
	  cqmin(k) = 0.0
9014	continue
	nqass = 0
	natom = 0
c
c
c read radius file and store in hash table
c
	open(unit=11,status='unknown',err=901)
	filnam = ' '
	inquire(11,name = filnam)
	write(6,*)'atom radii read from file'
	write(6,*)filnam
	write(6,*)' '
c
c skip/print comments (start with !) and one column header line 
c
105	read(11,201,end=901)comment
	if(comment(1:1).eq.'!') then
	  write(6,*)comment
	  goto 105
	end if
201   format(a)
	nrdrec = 0
100   continue
	  nrdrec = nrdrec + 1
	  if(nrdrec.gt.nrmax) then
	    write(6,*)' maximum # of radius records exceeded'
	    write(6,*)' increase nrmax'
	    stop 
	  end if
	  read(11,200,err=904,end=300)atm,res,rad
	  call up(atm,6)
	  call elb(atm,6)
	  call up(res,3)
	  call elb(res,3)
	  atnam(nrdrec) = atm
	  rnam(nrdrec)  = res
	  radt(nrdrec)  = rad
	  call rent(atm,res,nrdrec)
	goto 100
300   continue
	close(11)
200   format(A6,A3,F8.3)
	nrdrec = nrdrec - 1
	write(6,*)'# of radius parameter records:',nrdrec
c
c read charge parameter file and store in hash table
c
	open(unit=12,status='unknown',err=902)
	filnam = ' '
	inquire(12,name = filnam)
	write(6,*)' '
	write(6,*)'atomic charges read from file'
	write(6,*)filnam
	write(6,*)' '
c
c skip/print comments (start with !) and one column header line 
c
106	read(12,201,end=901)comment
	if(comment(1:1).eq.'!') then
	  write(6,*)comment
	  goto 106
	end if
	nchrec = 0
101   continue
	  nchrec = nchrec + 1
	  if(nchrec.gt.ncmax) then
	    write(6,*)' maximum # of charge records exceeded'
	    write(6,*)' - increase ncmax'
	    stop 
	  end if
	  read(12,202,err=905,end=301)atm,res,rnum,chn,chrgv
	  call up(atm,6)
	  call elb(atm,6)
	  call up(res,3)
	  call elb(res,3)
	  call up(rnum,4)
	  call elb(rnum,4)
	  call up(chn,1)
	  call elb(chn,1)
	  catnam(nchrec) = atm
	  crnam(nchrec)  = res
	  crnum(nchrec)  = rnum
	  cchn(nchrec)   = chn
	  chrgvt(nchrec) = chrgv
	  call cent(atm,res,rnum,chn,nchrec)
	goto 101
301   continue
	close(12)
202   format(A6,A3,A4,A1,F8.3)
	nchrec = nchrec - 1 
	write(6,*)'# of charge parameter records:',nchrec
c
c beginning of unformatted pdb file read
c 
	if(ipdbrd) then
c 
	open(unit=13,status='unknown',form='unformatted',err=903)
c 
	i=1
1711	continue
	read(13,end=2222) xo,rad,chrgv
	xn2(1,i)=xo(1)
	xn2(2,i)=xo(2)
	xn2(3,i)=xo(3)
	rad3(i)=rad
	chrgv4(i)=chrgv
	i=i+1
	goto 1711
2222	continue
	close(13)
	natom=i-1
c
	else
	call setrc(xn2,rad3,chrgv4,iatout,ipdbwrt,natom)
	end if
c
c find extrema and calculate scale according to them and
c to the percent box fill
c 
	cmin1=6000
	cmin2=6000
	cmin3=6000
	cmax1=-6000
	cmax2=-6000
	cmax3=-6000
	do 2103 ix=1,natom
	cmin1=amin1(cmin1,xn2(1,ix)-rad3(ix))
	cmin2=amin1(cmin2,xn2(2,ix)-rad3(ix))
	cmin3=amin1(cmin3,xn2(3,ix)-rad3(ix))
	cmax1=amax1(cmax1,xn2(1,ix)+rad3(ix))
	cmax2=amax1(cmax2,xn2(2,ix)+rad3(ix))
	cmax3=amax1(cmax3,xn2(3,ix)+rad3(ix))
2103	continue 
	oldmid(1)=(cmax1+cmin1)/2
	oldmid(2)=(cmax2+cmin2)/2
	oldmid(3)=(cmax3+cmin3)/2
	cran1=cmax1-cmin1
	cran2=cmax2-cmin2
	cran3=cmax3-cmin3 
	bran=amax1(cran1,cran2)
	bran=amax1(bran,cran3)
	midg=(igrid+1)/2
	rmidg=midg
	gran=igrid-1
c
c scale is:
c 
	if(perfil.ne.0.) scale=gran*perfil/(100.*bran) 
	if(perfil.eq.0.) perfil=100.*bran*scale/gran
c 
c scale factors and oldmid ( latter with offsets)
c 
	fpoh = fpi*scale
	do 2108 i = 1,3
	  oldmid(i) = oldmid(i) - offset(i)/scale
2108	continue
c
        if((offset(1).eq.999).or.(offset(1).eq.777)) then
	if(offset(1).eq.999) write(6,*)
     &       'modifying midpoints using frc input file'
	if(offset(1).eq.777) write(6,*)
     &       'modifying midpoints using fort.27'
	summid(1)=0.0
	summid(2)=0.0
	summid(3)=0.0
	frcnum=0.0
        if(offset(1).eq.999) open(unit=15,status='unknown',err=903)
        if(offset(1).eq.777) open(unit=27,status='unknown',err=903)
1124    read(15,204,end=828)line
        head = line(1:6)
	  call up(head,6)
c
c skip header lines
c
          if((head.ne.'ATOM  ').and.(head.ne.'HETATM')) then
	    go to 1124
	  end if
	  frcnum=frcnum + 1.0
	  crdstr = line(31:54)
        decode(24,205,crdstr)tempmid 
	summid(1)=summid(1)+tempmid(1)
	summid(2)=summid(2)+tempmid(2)
	summid(3)=summid(3)+tempmid(3)
	if((offset(2).eq.999).or.(offset(2).eq.777)) goto 828
	goto 1124
828	continue 
	if(frcnum.eq.0.0) then
	if(offset(1).eq.999) write(6,*) 'frc file empty of atoms for'
	if(offset(1).eq.777) write(6,*) 'unit 27 empty of atoms for'
	write(6,*) 'midpoint determination therefore assuming zero offsets'
	oldmid(1)=oldmid(1) + offset(1)/scale 
	oldmid(2)=oldmid(2) + offset(2)/scale 
	oldmid(3)=oldmid(3) + offset(3)/scale 
	else 
	oldmid(1)=summid(1)/frcnum
	oldmid(2)=summid(2)/frcnum
	oldmid(3)=summid(3)/frcnum
	end if 
c
	if(offset(1).eq.999) close(15)
	if(offset(1).eq.777) close(27)
        end if 
c
c 
c scale coordinates to grid
c 
c	addit1=rmidg-(oldmid(1)*scale)
c	addit2=rmidg-(oldmid(2)*scale)
c	addit3=rmidg-(oldmid(3)*scale)
c	do 2104 ix=1,natom
c	xn2(1,ix)=xn2(1,ix)*scale+addit1
c	xn2(2,ix)=xn2(2,ix)*scale+addit2
c	xn2(3,ix)=xn2(3,ix)*scale+addit3
c2104	continue 
c
	do 3104 ix=1,natom
	xo(1)=xn2(1,ix)
	xo(2)=xn2(2,ix)
	xo(3)=xn2(3,ix)
	call ctog(xo,xn)
	xn2(1,ix)=xn(1)
	xn2(2,ix)=xn(2)
	xn2(3,ix)=xn(3)
3104	continue 
c assign charges for boundary conditions
c 
	ic1=0
	do 2101 ix=1,natom
	if(chrgv4(ix).ne.0.) then
	ic1=ic1+1
	atmcrg(1,ic1)=xn2(1,ix)
	atmcrg(2,ic1)=xn2(2,ix)
	atmcrg(3,ic1)=xn2(3,ix)
	atmcrg(4,ic1)=chrgv4(ix)
	end if 
2101	continue 
c 
c ic = number of charges
c 
c find charge moments for dipole approximation
c 
	qnet=0.0
	qplus=0.0
	qmin=0.0
	do 2106 ix=1,3
	cqplus(ix)=0.0
	cqmin(ix)=0.0
2106	continue 
	do 2105 ix=1,ic1
	chrg=atmcrg(4,ix) 
	qnet=qnet + chrg  
	if(chrg.gt.0.) then
	qplus=qplus + chrg
	cqplus(1)=cqplus(1) + chrg*atmcrg(1,ix)
	cqplus(2)=cqplus(2) + chrg*atmcrg(2,ix)
	cqplus(3)=cqplus(3) + chrg*atmcrg(3,ix)
	else
	qmin=qmin + chrg
	cqmin(1)=cqmin(1) + chrg*atmcrg(1,ix)
	cqmin(2)=cqmin(2) + chrg*atmcrg(2,ix)
	cqmin(3)=cqmin(3) + chrg*atmcrg(3,ix)
	end if
2105    continue 
c
c divide by charge totals
c 
	if(qplus.ne.0.0) then
	  do 2110 k = 1,3
	    cqplus(k) = cqplus(k)/qplus
2110	  continue
	end if
	if(qmin.ne.0.0) then
	  do 2111 k = 1,3
	    cqmin(k) = cqmin(k)/qmin
2111	  continue
	end if
c
c select those charges which will be charging the grid
c 
	    nqass=ic1
	    rgrid=igrid
	    ic2=0
	do 2102 ix=1,ic1
	    if((atmcrg(1,ix).gt.1.).and.(atmcrg(1,ix).lt.rgrid)) then
	    if((atmcrg(2,ix).gt.1.).and.(atmcrg(2,ix).lt.rgrid)) then
	    if((atmcrg(3,ix).gt.1.).and.(atmcrg(3,ix).lt.rgrid)) then
	      ic2=ic2+1
	      chrgv2(1,ic2)=atmcrg(1,ix)
	      chrgv2(2,ic2)=atmcrg(2,ix)
	      chrgv2(3,ic2)=atmcrg(3,ix)
	      chrgv2(4,ic2)=atmcrg(4,ix)
	    end if
	    end if
	    end if
2102	continue 
	nqgrd=ic2 
c
c write details
c
	write(6,*)'  '
	write(6,*)'box fill      (%):',perfil
	write(6,*)'xmin,xmax     (A):',cmin1,cmax1
	write(6,*)'ymin,ymax     (A):',cmin2,cmax2
	write(6,*)'zmin,zma      (A):',cmin3,cmax3
	write(6,*)'x,y,z range   (A): ',cran1,cran2,cran3
	write(6,*)'scale   (grids/A): ',scale
	write(6,*)'object centre (A): ',oldmid
	write(6,*)'  '
c
c
c
c
	finish=cputime(start)
	write(6,*) 'time to read in and/or assign rad/chrg=',finish 
c
c
	write(6,*)'   '
	write(6,*)'number of atom coordinates read  : ',natom
	write(6,*)'total number of charged atoms    : ',nqass
	write(6,*)'net assigned charge              : ',qnet
	write(6,*)'assigned positive charge         : ',qplus
	write(6,*)'centred at (gu) :',cqplus
	write(6,*)'assigned negative charge         : ',qmin
	write(6,*)'centred at (gu) :',cqmin
	write(6,*)'   '
c
	call setout(xn2,rad3,natom,exrad,radprb)
c
c finish off dielectric map, producing solvent accessible
c volume defined by prob radius
c
	finish=cputime(start)
	write(6,*) 'time to turn everything out is',finish 
c check for boundary elements
c
	do i=1,3
	do iz=1,igrid
	do iy=1,igrid
	do ix=1,igrid
	if(iepsmp2(ix,iy,iz,i).ne.0) iepsmp(ix,iy,iz,i)=1
	end do
	end do
	end do
	end do
c
	it1(0)= .false.
	it1(6)= .false.
	do 487 ix=1,5
	it1(ix)= .true.
487	continue
	ibnum=0
	do 671 k=2,igrid-1
	  do 672 j=2,igrid-1
	    do 673 i=2,igrid-1
              ieps=   iepsmp(i,j,k,1) +
     &   iepsmp(i,j,k,2) +
     &   iepsmp(i,j,k,3) +
     &   iepsmp(i-1,j,k,1)+ 
     &   iepsmp(i,j-1,k,2) +
     &   iepsmp(i,j,k-1,3) 
	if(it1(ieps)) then
		ibnum=ibnum+1
		ibgrd(1,ibnum)=i
		ibgrd(2,ibnum)=j
		ibgrd(3,ibnum)=k
		ibgrd(4,ibnum)=iatmap(i,j,k)
	end if
673	continue
672	continue
671	continue
c
c reset iepsmap using these boundary elements
	write(6,*) "number of grid points on expanded surface= ",ibnum
c
	if(ihs .ne. 0) then
	write(6,*) "writing expanded surface data file: hsurf.dat"
	open(40,file="hsurf.dat")
	sahen=1.2079/(scale**2)
	do i=1,ibnum
	xo(1)=ibgrd(1,i)
	xo(2)=ibgrd(2,i)
	xo(3)=ibgrd(3,i)
	j=ibgrd(4,i)
	call gtoc(xo,xn)
	write(40,414) j,xn,sahen
	end do
414	format(i5,3f8.3,f8.3)
	close(40)
	end if
c
	call mkeps(radprb,epsout,ibnum,ibgrd,rad3,xn2)
c
c
c have used iepsmp2, now reassign iepsmp2
c
	finish=cputime(start)
	write(6,*) 'redo insides at ',finish 
c
	call setin(xn2,rad3,natom)
c
	finish=cputime(start)
	write(6,*) 'time to turn everything in is',finish 
c
c COMMENT OUT MEMBRANES FOR THE MOMENT
	 if(imem) call mem
c 
c recalculate boundary elements
c
	debfct = epsout/(deblen*scale)**2
	difeps=epsin-epsout
	sixeps=epsout*6.0
	sixth=1.0/6.0
	sixsalt=sixth*((1/(1+debfct/sixeps))-1.0)
	icgrid=igrid*igrid*igrid 
c
	if(rionst.ne.0.) then
	do 491 iz=0,1
	    do 493 ix=1,3
	    denom= sixeps + ix*difeps + iz*debfct
	    dbval(0,ix,iz)= 0.0
	    dbval(1,ix,iz)= difeps/denom
	    sfd(ix,iz)=epsout/denom
493	continue
491	continue
	do 492 iz=0,1
	    do 494 ix=4,5
	    denom= sixeps + ix*difeps + iz*debfct
	    dbval(1,ix,iz)= 0.0
	    dbval(0,ix,iz)= -difeps/denom
	    sfd(ix,iz)=epsin/denom
494	continue
492	continue
	else
	do 591 iz=0,1
	    do 593 ix=1,5
	    denom= sixeps + ix*difeps 
	    dbval(0,ix,iz)= (epsout/denom) -sixth 
	    dbval(1,ix,iz)= (epsin/denom) - sixth
593	continue
591	continue
	end if
c
c dbval for db
c 
	ibnum=0
	isgrid=igrid**2
	do 771 k=2,igrid-1
	  do 772 j=2,igrid-1
	    do 773 i=2,igrid-1
            ieps =  (iepsmp(i,j,k,1) + 
     &	                iepsmp(i,j,k,2) + 
     &	                iepsmp(i,j,k,3) + 
     &	                iepsmp(i-1,j,k,1) + 
     &	                iepsmp(i,j-1,k,2) + 
     &	                iepsmp(i,j,k-1,3)) 
	if(it1(ieps)) then
		ibnum=ibnum+1
		ibgrd(1,ibnum)=i
		ibgrd(2,ibnum)=j
		ibgrd(3,ibnum)=k
		ibgrd(4,ibnum)=iatmap2(i,j,k)
		iepsv2(ibnum)=ieps
	end if
773	continue
772	continue
771	continue
c
	ibnum1=0
	ibnum2=nsp/2
	do 774 ix=1,ibnum
	  i=ibgrd(1,ix)
	  j=ibgrd(2,ix)
	  k=ibgrd(3,ix)
	  ieps=iepsv2(ix)
	  iw=isgrid*(k-1) + igrid*(j-1) + i
	  iv=(iw+1)/2
	  deb=idebmap(i,j,k)
	  if(iw.ne.(2*iv)) then
		ibnum1=ibnum1+1
		ibnum3=ibnum1
	  else
		ibnum2=ibnum2+1
		ibnum3=ibnum2
	  end if
	  idpos(ibnum3) = iv
	  iepsv(ibnum3) = ieps
	  db(1,ibnum3)=dbval(iepsmp(i-1,j,k,1),ieps,int(deb))
	  db(2,ibnum3)=dbval(iepsmp(i,j,k,1),ieps,int(deb))
	  db(3,ibnum3)=dbval(iepsmp(i,j-1,k,2),ieps,int(deb))
	  db(4,ibnum3)=dbval(iepsmp(i,j,k,2),ieps,int(deb))
	  db(5,ibnum3)=dbval(iepsmp(i,j,k-1,3),ieps,int(deb))
	  db(6,ibnum3)=dbval(iepsmp(i,j,k,3),ieps,int(deb))
774	continue
c realign idpos and db
c
	icount2a=ibnum1
	icount2b=icount2a+ibnum2-(nsp/2)
	itemp=(nsp/2)
	do 781 ix=icount2a+1,icount2b
	itemp=itemp+1
	idpos(ix)=idpos(itemp)
	iepsv(ix)=iepsv(itemp)
	db(1,ix)=db(1,itemp)
	db(2,ix)=db(2,itemp)
	db(3,ix)=db(3,itemp)
	db(4,ix)=db(4,itemp)
	db(5,ix)=db(5,itemp)
	db(6,ix)=db(6,itemp)
781	continue
c
c
	write(6,*) 'number of dielectric boundary points', icount2b 
c
        if(isph)then
	   call chrg1up(nqgrd,chrgv2,qval,iqpos,fpoh,sixeps,
     &     difeps,debfct,gchrgp,gchrg,icount1a,icount1b)
	else
	call chrgup(nqgrd,chrgv2,qval,iqpos,fpoh,sixeps,difeps,
     &                 debfct,gchrgp,gchrg,icount1a,icount1b)
        endif
c 
c
c
	finish=cputime(start)
	write(6,*) 'iepsmp to db, and charging done at', finish 
c	
	write(6,*) 'number of grid points assigned charge', icount1b
c
c set saltmaps 1 and 2
	if(rionst.ne.0.) then 
	iw=1
	do 841 iz=1,igrid
	  do 842 iy=1,igrid
	    do 843 ix=1,igrid
		  phimap3(iw)=sixth + idebmap(ix,iy,iz)*sixsalt  
	      iw=iw+1
843	continue
842	continue
841	continue
	iy=0
	sf1((icgrid+1)/2)=phimap3(icgrid)
	do 850 ix=1,icgrid-2,2
	iy=iy+1
	sf1(iy)=phimap3(ix)
	sf2(iy)=phimap3(ix+1) 
850	continue 
	do 844 ix=1,icount2a
	  i=1
	  if(sf1(idpos(ix)).eq.sixth) i=0
	  sf1(idpos(ix))=sfd(iepsv(ix),i)
844	continue
	do 845 ix=icount2a+1,icount2b
	  i=1
	  if(sf2(idpos(ix)).eq.sixth) i=0
	  sf2(idpos(ix))=sfd(iepsv(ix),i)
845	continue
	end if 
c
c
c
c
c calculate boundary conditions
c
c
	write(6,*)'  '
	write(6,*)'setting boundary conditions...'
	write(6,*)'  '
c
c
	call setbc(ibctyp,iper,qplus,qmin,cqplus,cqmin,
     &          epsout,deblen,nqass)
c
c
c iterate
c
	call relfac(nlit,nnit,iper,idpos,db,sf1,sf2,
     &            icount2a,icount2b,rionst,spec)
c 
	noit=int(7.8/log(1.0 + sqrt(1-spec)))
	write(6,*) 'estimated iterations to convergence',noit
	if(iautocon) nlit = noit 
c
c 
	finish=cputime(start) 
	write(6,*)'  '
	write(6,*)'setup time was (sec) ',finish
	write(6,*)'  '
c
c
	if(nnit.eq.0.or.rionst.eq.0.0) then
	call itit(nlit,nnit,iper,idpos,db,sf1,sf2,iqpos,qval,
     &            icount2a,icount2b,icount1a,icount1b,rionst,spec,
     &            igraph,ipoten,icon1,icon2,diff,epsin,epsout)
c
	else 
c
	call nitit(nlit,nnit,iper,idpos,db,sf1,sf2,iqpos,qval,
     &            icount2a,icount2b,icount1a,icount1b,rionst,spec,
     &            igraph,ipoten,icon1,icon2,debfct,epsout)
	end if 
c
c
	if(loga)  then
	call anagrd(gchrgp,gchrg,icount1b,epsin*epkT,erga)
	write(6,*) 'analytic grid energy is        ',erga,' kt'
	end if
c
	if(logg) then
c
	ergg=0.0
	do 587 i=1,icount1b
	ix=gchrgp(1,i)
	iy=gchrgp(2,i)
	iz=gchrgp(3,i)
	ergg=ergg + phimap(ix,iy,iz)*gchrg(i)
587	continue
	ergg=ergg/2.0
	write(6,*) ' '
	write(6,*) 'total energy (including grid energy): ',ergg,' kt'
	end if 
c
	if(logg.and.loga) then
	write(6,*) 'difference energy, in kcals, is',(ergg-erga)*0.6
	end if
c 
	if(logs) then 
	if(diff) ibc=0
	epsval=epsin*epkT
	call react(ibgrd,nqass,icount2b,epsval,ergs,rad3,xn2)
	write(6,*) ' '
c	write(6,*) 'approximate solvation energy is',ergs,' kt'
	end if 
c
c
	if(logc) then 
	call clb(epsin,nqass,ergc)
	write(6,*) ' '
	write(6,*) 'coulombic energy:                     ',ergc,' kt'
	end if 
c 
	finish=cputime(start)
	write(6,*) 'energy calculations done at',finish 
	if(iconc) then
c
c convert potentials to concentrations
c
        nxtlbl = 'concentrat'
	  sixth = 1./6.
	  if(rionst.ne.0.0) then
	  write(6,*)'  '
	  write(6,*)'converting potentials to '
	  write(6,*)'net charge concentrations...'
	  write(6,*)'  '
	    do 9037 iz = 1,igrid
	      do 9038 iy = 1,igrid
		  do 9039 ix = 1,igrid
c
c use same number of terms in expansion
c of sinh as for iteration in itit.f
c
c		    temp = exp(-phimap(ix,iy,iz)) - exp(phimap(ix,iy,iz))
		    phi = phimap(ix,iy,iz)
		    phisq = phi**2
		    temp = phisq/120. + sixth
		    temp = phisq*temp + 1.
		    temp = 2.0*temp*phi
		    phimap(ix,iy,iz) = rionst*temp*idebmap(ix,iy,iz)
9039	    continue
9038	  continue
9037	continue
	  end if
	else
        nxtlbl = 'potential '
	end if
c
c write site potentials
c note if concentration option is selected, concentrations
c and conc. gradients are written
c
	if(ifrcrd.and.isite) goto 8888
c
c if site potentials required and unformatted read/write, skip
c 
      if(isite) then
c 
c during formatted frc file read/write can write unformatted frc.pdb
c 
      if(ifrcwrt) open(21,form='unformatted')
      pot=0.0
      ex=0.0
      ey=0.0
      ez=0.0
209	format(I5) 
	  write(6,*)'  '
	  write(6,*)'writing potentials at charge sites...'
	  write(6,*)'  '
        open(unit=15,status='unknown',err=903)
	  filnam = ' '
	  inquire(15,name = filnam)
	  write(6,*)'coordinates for potential output read from file'
	  write(6,*)filnam
	  write(6,*)'  '
        open(unit=16)
	  filnam = ' '
	  inquire(16,name = filnam)
	  write(6,*)'potentials written to file'
	  write(6,*)filnam
	  write(6,*)'  '
	  epsin = epsin*epkT
	  epsout = epsout*epkT
        write(16,*)'output from QDIFF   '
        write(16,*)'grid size,percent fill:',igrid,perfil
        write(16,*)'inner,outer dielectric:',epsin,epsout
        write(16,*)'ionic strength (M):',rionst
        write(16,*)'ion excl., probe radius:',exrad,radprb
        write(16,*)'linear, nolinear iterations:',nlit,nnit
        write(16,*)'boundary condition:',ibctyp
        write(16,*)'    '
        write(16,*)'title: '
        write(16,*)toplbl
	  if(iconc) then
          write(16,231)
          write(16,232)
231       format('         coordinates            charge
     $   concentration   gradient (M/ang)   ')
232       format('     x         y         z       (e)
     &      (M)         Cx       Cy        Cz')
	  else
          write(16,233)
          write(16,234)
233       format('         coordinates            charge
     $   potential         field (kT/e/Ang.)')
234       format('     x         y         z       (e)
     &      (kT/e)      Ex       Ey        Ez')
	  end if
c      123456789 123456789 123456789 123456789 123456789 123456789 123456789
	  etot = 0.
c
c read atom coordinate file
c
	  natom = 0
104     continue
        read(15,204,end=304)line
        head = line(1:6)
	  call up(head,6)
c
c skip header lines
c
        if((head.ne.'ATOM  ').and.(head.ne.'HETATM')) then
	    go to 104
	  end if
	  natom = natom + 1
	  crdstr = line(31:54)
	  atnum  = line(7:11)
        decode(24,205,crdstr)xo
	decode(5,209,atnum)inum 
c        decode(24,205,line(31:54))xo
	  atm = line(12:16)
	  res = line(18:20)
	  rnum = line(23:26)
	  chn = line(22:22)
	  call up(atm,6)
	  call elb(atm,6)
	  call up(res,3)
	  call elb(res,3)
	  call up(rnum,4)
	  call elb(rnum,4)
	  call up(chn,1)
	  call elb(chn,1)
c
c scale atoms to grid space
c
	  call ctog(xo,xn)
c
c assign charge to atom, searching for decreasingly specific specification
c note if no charge record found, is assumed to be 0.0
c
c
c
	  call cfind(atm,res,rnum,chn,ifind,n)
	  if(ifind.eq.0) then
	    schn = chn
	    chn = ' '
	    call cfind(atm,res,rnum,chn,ifind,n)
	    if(ifind.eq.0) then
		chn = schn
		snum = rnum
		rnum = '    '
	      call cfind(atm,res,rnum,chn,ifind,n)
	      if(ifind.eq.0) then
	        schn = chn
	        chn = ' '
	        call cfind(atm,res,rnum,chn,ifind,n)
		  if(ifind.eq.0) then
		    chn = schn
		    rnum = snum
		    sres = res
		    res = '   '
	          call cfind(atm,res,rnum,chn,ifind,n)
	          if(ifind.eq.0) then
	            schn = chn
	            chn = ' '
	            call cfind(atm,res,rnum,chn,ifind,n)
	            if(ifind.eq.0) then
		        chn = schn
		        snum = rnum
		        rnum = '    '
	              call cfind(atm,res,rnum,chn,ifind,n)
	              if(ifind.eq.0) then
	                schn = chn
	                chn = ' '
	                call cfind(atm,res,rnum,chn,ifind,n)
			  endif
			endif
		    endif
		  end if
		end if
	    end if
	  end if
	  chrgv = 0.0
	  if(ifind.ne.0) then
	    chrgv = chrgvt(n)
	  end if
c
c
	if(ifrcwrt)	write(21) inum,xo,chrgv,pot,ex,ey,ez
c
c calculate potential field and energy, and output to file
c
        call phintp(xn,phiv)
	  xn(1) = xn(1) + 1.
        call phintp(xn,fxu)
	  xn(1) = xn(1) - 2.
        call phintp(xn,fxl)
	  xn(1) = xn(1) + 1.
	  xn(2) = xn(2) + 1.
        call phintp(xn,fyu)
	  xn(2) = xn(2) - 2.
        call phintp(xn,fyl)
	  xn(2) = xn(2) + 1.
	  xn(3) = xn(3) + 1.
        call phintp(xn,fzu)
	  xn(3) = xn(3) - 2.
        call phintp(xn,fzl)
        qphiv = chrgv*phiv
	  etot = etot + qphiv
	  fx = (fxu-fxl)*scale/2.
	  fy = (fyu-fyl)*scale/2.
	  fz = (fzu-fzl)*scale/2.
	  write(16,230)xo,chrgv,phiv,fx,fy,fz
230     format(8F10.3)
        go to 104
304     continue		
c	end of file
	  close(15)
	if(ifrcwrt)  close(21) 
c	  natom = natom - 1
	  write(6,*)'   '
	  write(6,*)'number of atom coordinates read  : ',natom
	  write(6,*)'   '
	  etot = etot/2.
	  write(16,*)'total energy = ',etot,' kt'
	  close(16)
c
c end of formatted frc read/write and unformatted frc write
c 
      end if 
c
c  alternate frc read 
c 
8888	if(ifrcrd.and.isite) then
	open(15,form='unformatted')
	open(22,form='unformatted')
320	read(15,end=310) inum,xo,chrgv,phiv,fx,fy,fz
	call ctog(xo,xn)
	etot=0.0
        call phintp(xn,phiv)
	  xn(1) = xn(1) + 1.
        call phintp(xn,fxu)
	  xn(1) = xn(1) - 2.
        call phintp(xn,fxl)
	  xn(1) = xn(1) + 1.
	  xn(2) = xn(2) + 1.
        call phintp(xn,fyu)
	  xn(2) = xn(2) - 2.
        call phintp(xn,fyl)
	  xn(2) = xn(2) + 1.
	  xn(3) = xn(3) + 1.
        call phintp(xn,fzu)
	  xn(3) = xn(3) - 2.
        call phintp(xn,fzl)
        qphiv = chrgv*phiv
	  etot = etot + qphiv
	  fx = (fxu-fxl)*scale/2.
	  fy = (fyu-fyl)*scale/2.
	  fz = (fzu-fzl)*scale/2.
	  write(22) inum,xo,chrgv,phiv,fx,fy,fz 
        goto 320
310	close(15)
	close(22) 
c 
c END of unformatted frc.pdb read and frc write 
c 
	end if 
c
	finish=cputime(start)
	write(6,*) 'frc stuff now done at',finish 
c
c
c	if phiwrt set true then write potential map
c 
	if(phiwrt) then
c
c 
	if(ibios) then
c
c write phimap in insight format
c
	  open(unit=14,form='unformatted')
	  filnam = ' '
	  inquire(14,name = filnam)
	  write(6,*)'potential map written in INSIGHT format to file'
	  write(6,*)filnam
	  write(6,*)'  '
	  ivary = 0
	  nbyte = 4
	  intdat = 0
	  xang = 90.
	  yang = 90.
	  zang = 90.
	  intx = igrid - 1
	  inty = igrid - 1
	  intz = igrid - 1
	  xmax = 0.
	  do 9040 k = 1,3
	    temp = abs(oldmid(k))
	    xmax = amax1(xmax,temp)
9040	  continue
	  range = (igrid-1.)/(2.*scale)
	  extent = range + xmax
	  xstart = (oldmid(1)-range)/extent
	  ystart = (oldmid(2)-range)/extent
	  zstart = (oldmid(3)-range)/extent
	  xend = (oldmid(1)+range)/extent
	  yend = (oldmid(2)+range)/extent
	  zend = (oldmid(3)+range)/extent
        write(14)toplbl
	  write(14)ivary,nbyte,intdat,extent,extent,extent,
     1  xang,yang,zang,xstart,xend,ystart,yend,zstart,
     1  zend,intx,inty,intz
	  do 9041 k = 1,igrid
	    do 9042 j = 1,igrid
		write(14)(phimap(i,j,k),i=1,igrid)
9042	  continue
9041	  continue
c
	  write(6,*)'  '
	  write(6,*)'expanding dielectric map'
	  write(6,*)'to full grid of ',ngrid,'...'
	  write(6,*)'  '
	  call expand(ibios)
	else
c
c expand potential map to full (ngrid) size
c and write in DELPHI format
c
	  write(6,*)'  '
	  write(6,*)'expanding potential and dielectric maps'
	  write(6,*)'to full grid of ',ngrid,'...'
	  write(6,*)'writing potential map in DELPHI format'
	  write(6,*)'  '
	  if(igrid.ne.ngrid) then 
	  call expand(ibios)
	  end if 
	  open(unit=14,form='unformatted')
	  filnam = ' '
	  inquire(14,name = filnam)
	  write(6,*)'potential map written to file'
	  write(6,*)filnam
	  write(6,*)'  '
        write(14)'now starting phimap '
        write(14)nxtlbl,toplbl
        write(14)phimap
        write(14)' end of phimap  '
        write(14)scale,oldmid
        close(14)
	end if
c
	end if 
c write dielectric map
c
	call wrteps

	goto 999
900	write(6,*) 'unexpected end or non-existence of parameter file'
	stop 
901	write(6,*) 'unexpected end or non-existence of radius file'
   	stop
902	write(6,*) 'unexpected end or non-existence of charge file'
   	stop
903	write(6,*) 'unexpected end or non-existence of atom file'
   	stop
904	write(6,*) 'error in reading radius file'
   	stop
905	write(6,*) 'error in reading charge file'
   	stop
999	continue
	finish = cputime(start)
c	finish = finish - start
	write(6,*)'  '
	write(6,*)'total cpu time was (sec) ',finish
	write(6,*)'  '
	write(6,*)'QDIFF QDIF. QDI.. QD... Q....'
	end