      subroutine delphi(iswitch,iteration,ncenters,nptcharges,ngrid,
     &                  ngridby2,natomtype,nboundary,epsilon,coord,
     &                  charges,coordptcharge,ptcharges,phinew,phiold,
     &                  gridcharge,dielectric,avgdielectric,
     &                  coordneighbour,radneighbour)

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)

      logical connol

      dimension vanderwaal(18),cubepts(3,6)

      dimension natomtype(*),nboundary(ngrid,ngrid,*)

      dimension coord(3,*),charges(*),coordptcharge(3,*),ptcharges(*),
     &          phinew(ngrid,ngrid,*),phiold(ngrid,ngrid,*),
     &          gridcharge(ngrid,ngrid,*),dielectric(6,ngrid,ngrid,*),
     &          avgdielectric(ngrid,ngrid,*),coordneighbour(3,*),
     &          radneighbour(*)

      data zero,half,one/0.0d0,0.5d0,1.0d0/
      data pi,six/3.1415926535898d0,6.0d0/
      data bohr,efactor1,efactor2/0.529177d0,27.212d0,23.062d0/
      data hstep,tolerance,probe/0.50d0,0.0001d0,1.40d0/

      data vanderwaal/
     &     1.40d0,                                               1.62d0,
     &     2.02d0,0.00d0,     0.00d0,1.92d0,1.77d0,1.62d0,1.57d0,1.76d0,
     &     2.52d0,1.92d0,     0.00d0,2.32d0,2.07d0,2.07d0,2.02d0,2.14d0/

      rewind 49
      rewind 50

      fourpi=pi/half/half

      proberadius=probe/bohr

      do 1001 k=1,ngrid
      do 1001 j=1,ngrid
      do 1001 i=1,ngrid
      nboundary(i,j,k)=0
      gridcharge(i,j,k)=zero
 1001 phiold(i,j,k)=zero

      do 1002 i=1,ncenters

      if(natomtype(i).eq.0) goto 1002

      ix=int((coord(1,i)+real64(ngridby2)*hstep)/hstep)
      iy=int((coord(2,i)+real64(ngridby2)*hstep)/hstep)
      iz=int((coord(3,i)+real64(ngridby2)*hstep)/hstep)

      ix1=ix
      iy1=iy
      iz1=iz

      ix2=ix
      iy2=iy
      iz2=iz+1

      ix3=ix
      iy3=iy+1
      iz3=iz

      ix4=ix
      iy4=iy+1
      iz4=iz+1

      ix5=ix+1
      iy5=iy
      iz5=iz

      ix6=ix+1
      iy6=iy
      iz6=iz+1

      ix7=ix+1
      iy7=iy+1
      iz7=iz

      ix8=ix+1
      iy8=iy+1
      iz8=iz+1

      x1=real64(ix1-ngridby2)*hstep
      y1=real64(iy1-ngridby2)*hstep
      z1=real64(iz1-ngridby2)*hstep

      x2=real64(ix2-ngridby2)*hstep
      y2=real64(iy2-ngridby2)*hstep
      z2=real64(iz2-ngridby2)*hstep

      x3=real64(ix3-ngridby2)*hstep
      y3=real64(iy3-ngridby2)*hstep
      z3=real64(iz3-ngridby2)*hstep

      x4=real64(ix4-ngridby2)*hstep
      y4=real64(iy4-ngridby2)*hstep
      z4=real64(iz4-ngridby2)*hstep

      x5=real64(ix5-ngridby2)*hstep
      y5=real64(iy5-ngridby2)*hstep
      z5=real64(iz5-ngridby2)*hstep

      x6=real64(ix6-ngridby2)*hstep
      y6=real64(iy6-ngridby2)*hstep
      z6=real64(iz6-ngridby2)*hstep

      x7=real64(ix7-ngridby2)*hstep
      y7=real64(iy7-ngridby2)*hstep
      z7=real64(iz7-ngridby2)*hstep

      x8=real64(ix8-ngridby2)*hstep
      y8=real64(iy8-ngridby2)*hstep
      z8=real64(iz8-ngridby2)*hstep

      x1=abs(x1-coord(1,i))/hstep
      y1=abs(y1-coord(2,i))/hstep
      z1=abs(z1-coord(3,i))/hstep

      x2=abs(x2-coord(1,i))/hstep
      y2=abs(y2-coord(2,i))/hstep
      z2=abs(z2-coord(3,i))/hstep

      x3=abs(x3-coord(1,i))/hstep
      y3=abs(y3-coord(2,i))/hstep
      z3=abs(z3-coord(3,i))/hstep

      x4=abs(x4-coord(1,i))/hstep
      y4=abs(y4-coord(2,i))/hstep
      z4=abs(z4-coord(3,i))/hstep

      x5=abs(x5-coord(1,i))/hstep
      y5=abs(y5-coord(2,i))/hstep
      z5=abs(z5-coord(3,i))/hstep

      x6=abs(x6-coord(1,i))/hstep
      y6=abs(y6-coord(2,i))/hstep
      z6=abs(z6-coord(3,i))/hstep

      x7=abs(x7-coord(1,i))/hstep
      y7=abs(y7-coord(2,i))/hstep
      z7=abs(z7-coord(3,i))/hstep

      x8=abs(x8-coord(1,i))/hstep
      y8=abs(y8-coord(2,i))/hstep
      z8=abs(z8-coord(3,i))/hstep

      gridcharge(ix1,iy1,iz1)=gridcharge(ix1,iy1,iz1)
     &                       +(one-x1)*(one-y1)*(one-z1)*charges(i)

      gridcharge(ix2,iy2,iz2)=gridcharge(ix2,iy2,iz2)
     &                       +(one-x2)*(one-y2)*(one-z2)*charges(i)

      gridcharge(ix3,iy3,iz3)=gridcharge(ix3,iy3,iz3)
     &                       +(one-x3)*(one-y3)*(one-z3)*charges(i)

      gridcharge(ix4,iy4,iz4)=gridcharge(ix4,iy4,iz4)
     &                       +(one-x4)*(one-y4)*(one-z4)*charges(i)

      gridcharge(ix5,iy5,iz5)=gridcharge(ix5,iy5,iz5)
     &                       +(one-x5)*(one-y5)*(one-z5)*charges(i)

      gridcharge(ix6,iy6,iz6)=gridcharge(ix6,iy6,iz6)
     &                       +(one-x6)*(one-y6)*(one-z6)*charges(i)

      gridcharge(ix7,iy7,iz7)=gridcharge(ix7,iy7,iz7)
     &                       +(one-x7)*(one-y7)*(one-z7)*charges(i)

      gridcharge(ix8,iy8,iz8)=gridcharge(ix8,iy8,iz8)
     &                       +(one-x8)*(one-y8)*(one-z8)*charges(i)

 1002 continue

      if(iteration.gt.1) call fastrd(50,phiold,ngrid*ngrid*ngrid)

      da=real64(1-ngridby2)*hstep
      db=real64(ngrid-ngridby2)*hstep

      do 1003 j=1,ngrid
      do 1003 i=1,ngrid

      phiold(i,j,1)=zero
      phiold(i,1,j)=zero
      phiold(1,i,j)=zero

      phiold(i,j,ngrid)=zero
      phiold(i,ngrid,j)=zero
      phiold(ngrid,i,j)=zero

 1003 continue

      do 1004 k=1,ncenters

      if(natomtype(k).ne.0) then

               xx=coord(1,k)
               yy=coord(2,k)
               zz=coord(3,k)

               do 1005 j=2,ngrid-1

               dj=real64(j-ngridby2)*hstep

               do 1005 i=2,ngrid-1

               di=real64(i-ngridby2)*hstep

               phiold(i,j,1)=phiold(i,j,1)+charges(k)/epsilon/
     &                       sqrt((di-xx)**2+(dj-yy)**2+(da-zz)**2)

               phiold(i,1,j)=phiold(i,1,j)+charges(k)/epsilon/
     &                       sqrt((di-xx)**2+(da-yy)**2+(dj-zz)**2)

               phiold(1,i,j)=phiold(1,i,j)+charges(k)/epsilon/
     &                       sqrt((da-xx)**2+(di-yy)**2+(dj-zz)**2)

               phiold(i,j,ngrid)=phiold(i,j,ngrid)+charges(k)/epsilon/
     &                           sqrt((di-xx)**2+(dj-yy)**2+(db-zz)**2)

               phiold(i,ngrid,j)=phiold(i,ngrid,j)+charges(k)/epsilon/
     &                           sqrt((di-xx)**2+(db-yy)**2+(dj-zz)**2)

               phiold(ngrid,i,j)=phiold(ngrid,i,j)+charges(k)/epsilon/
     &                           sqrt((db-xx)**2+(di-yy)**2+(dj-zz)**2)

 1005          continue

                            endif
 1004 continue

      do 1006 j=1,ngrid
      do 1006 i=1,ngrid
      phinew(i,j,1)=phiold(i,j,1)
      phinew(i,1,j)=phiold(i,1,j)
      phinew(1,i,j)=phiold(1,i,j)
      phinew(i,j,ngrid)=phiold(i,j,ngrid)
      phinew(i,ngrid,j)=phiold(i,ngrid,j)
 1006 phinew(ngrid,i,j)=phiold(ngrid,i,j)

      if(iteration.gt.1) then

                call fastrd(49,dielectric,6*ngrid*ngrid*ngrid)

                         else

                do 1007 k=1,ngrid
                do 1007 j=1,ngrid
                do 1007 i=1,ngrid
                dielectric(1,i,j,k)=epsilon
                dielectric(2,i,j,k)=epsilon
                dielectric(3,i,j,k)=epsilon
                dielectric(4,i,j,k)=epsilon
                dielectric(5,i,j,k)=epsilon
 1007           dielectric(6,i,j,k)=epsilon

                do 1008 l=1,ncenters

                if(natomtype(l).eq.0) goto 1008 

                if(natomtype(l).eq.04) stop 'cannot do delphi'
                if(natomtype(l).eq.05) stop 'cannot do delphi'
                if(natomtype(l).eq.13) stop 'cannot do delphi'
                if(natomtype(l).gt.18) stop 'cannot do delphi'

                radius=vanderwaal(natomtype(l))/bohr

                neighbour=0

                do 1009 i=1,ncenters
                if((i.ne.l).and.(natomtype(i).ne.0)) then

                   neighbour=neighbour+1

                   coordneighbour(1,neighbour)=coord(1,i)
                   coordneighbour(2,neighbour)=coord(2,i)
                   coordneighbour(3,neighbour)=coord(3,i)

                   radneighbour(neighbour)=vanderwaal(natomtype(i))/bohr

                                                     endif
 1009           continue

                do 1010 k=2,ngrid-1
                z=real64(k-ngridby2)*hstep

                do 1010 j=2,ngrid-1
                y=real64(j-ngridby2)*hstep

                do 1010 i=2,ngrid-1
                x=real64(i-ngridby2)*hstep

                cubepts(1,1)=x+hstep*half
                cubepts(2,1)=y
                cubepts(3,1)=z

                cubepts(1,2)=x-hstep*half
                cubepts(2,2)=y
                cubepts(3,2)=z

                cubepts(1,3)=x
                cubepts(2,3)=y+hstep*half
                cubepts(3,3)=z

                cubepts(1,4)=x
                cubepts(2,4)=y-hstep*half
                cubepts(3,4)=z

                cubepts(1,5)=x
                cubepts(2,5)=y
                cubepts(3,5)=z+hstep*half

                cubepts(1,6)=x
                cubepts(2,6)=y
                cubepts(3,6)=z-hstep*half

       do 1011 m=1,6

       if(dist(cubepts(1,m),coord(1,l)).le.radius)
     &                                  dielectric(m,i,j,k)=one

       if((dist(cubepts(1,m),coord(1,l)).gt.radius).and.
     &    (dist(cubepts(1,m),coord(1,l)).le.(radius+proberadius))) then

       if(connol(cubepts(1,m),proberadius,coordneighbour,radneighbour,
     &    coord(1,l),radius,neighbour)) dielectric(m,i,j,k)=one

                                                                   endif

 1011  continue

 1010           continue

 1008           continue

                rewind 49

                call fastwr(49,dielectric,6*ngrid*ngrid*ngrid)

                         endif

      do 1012 k=2,ngrid-1
      do 1012 j=2,ngrid-1
      do 1012 i=2,ngrid-1
 1012 avgdielectric(i,j,k)=dielectric(1,i,j,k)+dielectric(2,i,j,k)
     &                    +dielectric(3,i,j,k)+dielectric(4,i,j,k)
     &                    +dielectric(5,i,j,k)+dielectric(6,i,j,k)

      do 1013 k=2,ngrid-1
      do 1013 j=2,ngrid-1
      do 1013 i=2,ngrid-1
      a1=abs(six-avgdielectric(i,j,k))
      a2=abs(six-avgdielectric(i,j,k)/epsilon)
 1013 if(min(a1,a2).gt.tolerance) nboundary(i,j,k)=1

      do 1014 n=1,150

      do 1015 k=2,ngrid-1
      do 1015 j=2,ngrid-1
      do 1015 i=2,ngrid-1

      phinew(i,j,k)=phiold(i+1,j,k)*dielectric(1,i,j,k)
     &             +phiold(i-1,j,k)*dielectric(2,i,j,k)
     &             +phiold(i,j+1,k)*dielectric(3,i,j,k)
     &             +phiold(i,j-1,k)*dielectric(4,i,j,k)
     &             +phiold(i,j,k+1)*dielectric(5,i,j,k)
     &             +phiold(i,j,k-1)*dielectric(6,i,j,k)
     &             +fourpi*gridcharge(i,j,k)/hstep

 1015 phinew(i,j,k)=phinew(i,j,k)/avgdielectric(i,j,k)

      do 1016 k=2,ngrid-1
      do 1016 j=2,ngrid-1
      do 1016 i=2,ngrid-1
 1016 phiold(i,j,k)=phinew(i,j,k)

 1014 continue

      rewind 50

      call fastwr(50,phiold,ngrid*ngrid*ngrid)

      do 1017 k=2,ngrid-1
      do 1017 j=2,ngrid-1
      do 1017 i=2,ngrid-1
 1017 phiold(i,j,k)=(hstep/fourpi)*(six*phinew(i,j,k)
     &                             -phinew(i+1,j,k)-phinew(i-1,j,k)
     &                             -phinew(i,j+1,k)-phinew(i,j-1,k)
     &                             -phinew(i,j,k+1)-phinew(i,j,k-1))
     &                             -gridcharge(i,j,k)

      nptcharges=0

      if(iswitch.eq.1) write(6,1018)
 1018 format(/,' surface induced charges',
     &         ' (can be displayed with MIDAS): ',/)

      do 1019 k=2,ngrid-1
      do 1019 j=2,ngrid-1
      do 1019 i=2,ngrid-1

      if(nboundary(i,j,k).eq.0) goto 1019

      nptcharges=nptcharges+1

      x=real64(i-ngridby2)*hstep
      y=real64(j-ngridby2)*hstep
      z=real64(k-ngridby2)*hstep

      coordptcharge(1,nptcharges)=x
      coordptcharge(2,nptcharges)=y
      coordptcharge(3,nptcharges)=z

      ptcharges(nptcharges)=phiold(i,j,k)

      if(iswitch.eq.1) then

         if(phiold(i,j,k).lt.zero) then

             write(6,1020) nptcharges,nptcharges+1,x,y,z,phiold(i,j,k)
 1020        format('ATOM',3x,i4,'  CL- CIO ',i5,4x,3f8.3,f12.9)

                                   else

             write(6,1021) nptcharges,nptcharges+1,x,y,z,phiold(i,j,k)
 1021        format('ATOM',3x,i4,'  NA+ CIO ',i5,4x,3f8.3,f12.9)

                                   endif

                       endif
 1019 continue

      deltag=zero

      do 1022 j=1,ncenters
      do 1022 i=1,nptcharges
 1022 deltag=deltag+half*ptcharges(i)*charges(j)/
     &                               dist(coordptcharge(1,i),coord(1,j))

      write(6,1023) nptcharges
 1023 format(/,' there are ',i5,' surface induced charges',/)

      write(6,1024) deltag*efactor1*efactor2
 1024 format(' electrostatic contribution to free energy of solvation',
     &       ' = ',f9.4,' kcal/mole',////)

      return
      end
