PROGRAM DODCG
      IMPLICIT none
      integer*4 MX,iper,nt0,it,i,N,J,ii,nc,ns,jx,MXB,n6,npar0,
     1ig,is,ix,NG,nat,ncs,ierr,ITER,NATB,nwl,nwl0,n4,n40,
     2ivel,imag,ntr
      PARAMETER (MX=4600,nt0=550,MXB=15,npar0=14,n40=6,nwl0=nt0)
c     Limits:
c     MX -# of atoms in peptide
c     nt0 -# of transitions on one chromophore
c     MXB -# of overlaped atoms
c     npar0 -# of parameters to follow
c     n40 -# of potential sites
c     nwl0 -# of fitted transitions
      integer*4 ic(MX),ics(MX),ttype(MX),IP(npar0)
      REAL*8 R(3*MX),rc(3),wc(MX),rn(3*MX),mag(nt0,3),elel(nt0,3),
     1elev(nt0,3),en(nt0),sum,w,A,XS,XB,RTOL,um(3),ue(3),dco,
     1PG(nt0,npar0),zim(n40),dco0,
     1buvcd(nwl0,n40),fi(n40),F(npar0),enc(nt0),
     1dcn,dcn0,dhn,dhn0,DPCX,DPCY,DPCZ,elelc(nt0,3),rcs(3)
      character*80 s80
      logical*4 lper,lfdp(4)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NATB,ITER
c
c     takes electronic dipoles from G.OUT, spreads them on large peptide
c
c     Input: G.OUT       --   TDDFT G Output
c            tinker.xyz  --- large molecule
c
c     Read in transition energies and dipoles from gaussian:
      open(30,file='G.OUT',status='old')
10    read(30,300,end=949,err=949)s80
300   format(a80)
      if(s80(2:44).eq.'Ground to excited state Transition electric')
     1call rd(nt0,it,elel)
      if(s80(2:44).eq.'Ground to excited state transition velocity')
     1call rd(nt0,it,elev)
      if(s80(2:44).eq.'Ground to excited state transition magnetic')
     1call rd(nt0,it,mag)
      if(s80(2:14).eq.'Excited State')call re(en)
      if(s80(26:46).eq.'Standard orientation:')call ro(rn,nat)
      goto 10
949   close(30)
      write(6,*)it,' transitions in G.OUT'
      
c     read-in geometry
      OPEN(2,FILE='tinker.xyz',STATUS='OLD')
      read(2,*)N
      if(N.gt.MX)call report('Too many atoms.')
      do 1 i=1,N
1     READ(2,2000)(R(J+3*(i-1)),J=1,3),ttype(i)
2000  format(13x,3f12.6,i6)
      close(2)
      write(6,*)N, ' atoms in tinker.xyz'

      open(4,file='DODC.CON')
      read(4,*)ivel,imag
      if(ivel.eq.1)then
       write(6,*)' velocity dipoles used'
       do 48 i=1,it
       do 48 ix=1,3
48     elel(i,ix)=elev(i,ix)
      else
       write(6,*)' radial dipoles used'
c
c      gaussians gives them with wrong sign!:
       do 481 i=1,it
       do 481 ix=1,3
481    elel(i,ix)=-elel(i,ix)
      endif
c
c     transform aus to debyes, debyes*A
      do 49 i=1,it
      do 49 ix=1,3
      mag(i,ix)=mag(i,ix)*2.541765d0*0.529177d0
49    elel(i,ix)=elel(i,ix)*2.541765d0


      if(imag.eq.1)then
       write(6,*)' magnetic moments added'
      else
       write(6,*)' magnetic moments skipped'
       do 47 i=1,it
       do 47 ix=1,3
47     mag(i,ix)=0.0d0
      endif
      open(3,file='DC.SC')
c     nma chromophore:   H N C O
      read(4,*)ncs,(ics(i),i=1,ncs)
      dco0=dsqrt((rn(1+3*(ics(4)-1))-rn(1+3*(ics(3)-1)))**2+
     1           (rn(2+3*(ics(4)-1))-rn(2+3*(ics(3)-1)))**2+
     1           (rn(3+3*(ics(4)-1))-rn(3+3*(ics(3)-1)))**2)
      dcn0=dsqrt((rn(1+3*(ics(2)-1))-rn(1+3*(ics(3)-1)))**2+
     1           (rn(2+3*(ics(2)-1))-rn(2+3*(ics(3)-1)))**2+
     1           (rn(3+3*(ics(2)-1))-rn(3+3*(ics(3)-1)))**2)
      dhn0=dsqrt((rn(1+3*(ics(2)-1))-rn(1+3*(ics(1)-1)))**2+
     1           (rn(2+3*(ics(2)-1))-rn(2+3*(ics(1)-1)))**2+
     1           (rn(3+3*(ics(2)-1))-rn(3+3*(ics(1)-1)))**2)


      read(4,*)NG
      write(3,*)NG,' groups'
      write(6,*)NG,' groups'
c
c     geometry and ES field corrections:
      read(4,*)(IP(i),i=1,npar0)
      write(6,6009)(IP(i),i=1,npar0)
6009  format(14i2)
      read(4,*)ntr
      do 6 it=1,ntr
6     read(4,*)(PG(it,i),i=1,npar0)
c
c     determine which reference value to use
      if(IP(1).lt.0)then
       write(6,*)' Ab initio value of d(C=O):'
      else
       dco0=1.231d0
       write(6,*)' MD value of d(C=O):'
      endif
      write(6,6070)dco0
 6070 format(f12.3,' A')
      if(IP(3).lt.0)then
       write(6,*)' Ab initio value of d(C=N):'
      else
       dcn0=1.335d0
       write(6,*)' MD value of d(C=N):'
      endif
      write(6,6070)dcn0
c
c     periodic boundary conditions
      read(4,*)iper,DPCX,DPCY,DPCZ
      lper=iper.eq.1
c
c     potential correction:
      if(IP(npar0).eq.1)then
       open(55,file='BIJ.SCR',status='old',form='unformatted')
       read(55)n4,n6,(zim(i),i=1,n4),(lfdp(i),i=1,4)
       write(6,*)'Zs:'
       write(6,6071)(zim(i),i=1,n4)
 6071  format(6f12.4)
       if(lfdp(1))read(55)
       if(lfdp(2))read(55)
       if(lfdp(3))read(55)
       if(lfdp(4))read(55)nwl,((buvcd(i,j),i=1,nwl),j=1,n6)
       close(55)
      endif

      do 2 IG=1,NG
      write(3,*)IG,' . group:'
      write(6,*)IG,' . group:'
c     Number of atoms determining center and their numbers:
      read(4,*)nc,(ic(i),i=1,nc)
      if(nc.ne.ncs)call report('nc<>ncs!')
c     center weights:'
      read(4,*)(wc(i),i=1,nc)
      sum=0.0d0
      do 321 ii=1,nc
321   sum=sum+wc(ii)
      do 3 i=1,3
      rc(i)=0.0d0
      rcs(i)=0.0d0
      do 32 ii=1,nc
      rcs(i)=rcs(i)+rn(i+3*(ics(ii)-1))*wc(ii)
32    rc(i)=rc(i)+R(i+3*(ic(ii)-1))*wc(ii)
      rcs(i)=rcs(i)/sum
3     rc(i)=rc(i)/sum
      write(3,1000)(rc(i)/10.0d0,i=1,3)
1000  format(3f15.6,' nm')
      
      ns=ntr+1
      write(3,*)ns,' states'
      write(3,*)
c
c     ground state:
      write(3,*)'0 0.0'
c
c     correct energies and dipoles according to
c     geometry and ES field:w=w0*(1+dd/d*p) dip=dip0*(1+dd/d*p')
c     first par-energy, second dipole
c      
      do 41 is=2,ns

      do 46 i=1,npar0
46    F(i)=0.0d0

      if(IP(1).ne.0)then
       dco=dsqrt((R(1+3*(ic(4)-1))-R(1+3*(ic(3)-1)))**2+
     1           (R(2+3*(ic(4)-1))-R(2+3*(ic(3)-1)))**2+
     1           (R(3+3*(ic(4)-1))-R(3+3*(ic(3)-1)))**2)
       F(1)=(dco-dco0)/dco0*PG(is-1,1)
       F(2)=(dco-dco0)/dco0*PG(is-1,2)
       if(is.eq.2)write(6,6001)1,dco,dco0
6001   format(i6,' d(C=O) = ',f10.5,'; D0 = ',f10.5)
      endif

c     C=N bond
      if(IP(3).ne.0)then
       dcn=dsqrt((R(1+3*(ic(2)-1))-R(1+3*(ic(3)-1)))**2+
     1           (R(2+3*(ic(2)-1))-R(2+3*(ic(3)-1)))**2+
     1           (R(3+3*(ic(2)-1))-R(3+3*(ic(3)-1)))**2)
       F(3)=(dcn-dcn0)/dcn0*PG(is-1,3)
       F(4)=(dcn-dcn0)/dcn0*PG(is-1,4)
       if(is.eq.2)write(6,6003)3,dcn,dcn0
6003   format(i6,' d(C-N) = ',f10.5,'; D0 = ',f10.5)
      endif

c     H-N bond
      if(IP(5).ne.0)then
       dhn=dsqrt((R(1+3*(ic(2)-1))-R(1+3*(ic(1)-1)))**2+
     1           (R(2+3*(ic(2)-1))-R(2+3*(ic(1)-1)))**2+
     1           (R(3+3*(ic(2)-1))-R(3+3*(ic(1)-1)))**2)
       F(5)=(dhn-dhn0)/dhn0*PG(is-1,5)
       F(6)=(dhn-dhn0)/dhn0*PG(is-1,6)
       if(is.eq.2)write(6,6004)5,dhn,dhn0
6004   format(i6,' d(N-H) = ',f10.5,'; D0 = ',f10.5)
      endif

      enc(is-1)=en(is-1)*(1.0d0+F(1)+F(3)+F(5))

      do 41 ix=1,3 
41    elelc(is-1,ix)=elel(is-1,ix)*(1.0d0+F(2)+F(4)+F(6))

      if(IP(npar0).ne.1)then
       call fw(n6,fi,N,ttype,R,ic(1),ic(2),ic(3),ic(4),lper,
     1 DPCX,DPCY,DPCZ)
       do 42 is=2,ns
       if(en(is-1).gt.200.0d0)then
	do 43 i=1,4
43      enc(is-1)=enc(is-1)+buvcd(1,i)*fi(i)
       else
	do 44 i=1,4
44      enc(is-1)=enc(is-1)+buvcd(2,i)*fi(i)
       endif
42     continue
      endif

      do 4 is=2,ns
      w=1.0d7/enc(is-1)
4     write(3,1001)is,w/8065.54093D0,w,enc(is-1)
1001  format(i3,f11.6,' eV =',f11.2,' cm-1 =',f11.2,' nm')
      write(3,*)
c
c     write group reference atoms into XS and corresponding ones in XB:
      if(ncs.gt.MXB)call report('ncs>MXB')
      do 5 i=1,ncs
      do 5 ix=1,3
      XB(i,ix)= R(ix+3*(ic(i)-1))-rc(ix)
5     XS(i,ix)=rn(ix+3*(ics(i)-1))-rcs(ix)
      NATB=ncs
c
c     overlap XS,XB, make the rotation matrix, store in A:
      call domatrix(ierr)
      if(ierr.ne.0)call report('rotation not successfull')
      write(6,*)' Transformation matrix'
      write(6,6767)A
6767  format(3(3f10.5,/))
c
c     write rotated dipole moments:
      do 2 is=2,ns

      do 8 ix=1,3
      ue(ix)=0.0d0
      um(ix)=0.0d0
      do 8 jx=1,3
      ue(ix)=ue(ix)+A(ix,jx)*elelc(is-1,jx)
8     um(ix)=um(ix)+A(ix,jx)* mag(is-1,jx)
c
2     write(3,1002)is,(ue(ix),ix=1,3),(um(ix),ix=1,3)
1002  format(i3,6f12.6)

      close(3)
      close(4)
      call report('DC.SC written')
      end
c
      SUBROUTINE REPORT(RS)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*(*) RS
      WRITE(6,3000)
3000  FORMAT(/,80(1H*))
      WRITE(6,*)RS
      WRITE(6,3001)
3001  FORMAT(80(1H*),/,/,'PROGRAM STOPPED')
      STOP
      END

      subroutine rd(nt0,it,e)
      implicit none
      integer*4 nt0,it,idumm,i
      character*10 s10
      real*8 e(nt0,3)
      logical nu
      it=0
      read(30,*)
31    read(30,300)s10
300   format(a10)
      if(nu(s10(10:10)))then
       it=it+1
       if(it.gt.nt0)then
        write(6,*)' too many transitions !'
	stop
       endif
       backspace 30
       read(30,*)idumm,(e(it,i),i=1,3)
       goto 31
      endif
      backspace 30
      return
      end

      subroutine ro(rn,nat)
      implicit none
      integer*4 dum,i,nat
      character*80 s80
      real*8 rn(*)
      logical nu

      nat=0
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(nu(s80(5:5)))then
       nat=nat+1
       read(s80,*)dum,dum,dum,(rn(i+3*(nat-1)),i=1,3)
       goto 31
      endif
      write(6,*)nat,' atoms in chromophore'
      return
      end

      subroutine re(e)
      implicit none
      integer*4 it
      character*80 s80
      real*8 e(*)
      backspace 30
      read(30,300)s80
300   format(a80)
      read(s80(15:18),*)it
      read(s80(47:55),*)e(it)
      return
      end

      function nu(c)
      character*1 c
      logical nu
      nu=c.eq.'1'.or.c.eq.'2'.or.c.eq.'3'.or.c.eq.'4'.or.c.eq.'5'.
     1or.c.eq.'6'.or.c.eq.'7'.or.c.eq.'8'.or.c.eq.'9'.or.c.eq.'0'
      return
      end

      SUBROUTINE DOMATRIX(IERR)
C     This is a Fortran version of Petr Malon's subroutine Amoeba
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=15)
      DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      IERR=0
C     STARTING VALUES FOR THE ITERATION:
      ANG(1)=0.2d0
      ANG(2)=0.0d0
      ANG(3)=0.0d0
      DO 1 I=1,3
1     P(1,I)=ANG(I)
      Y(1)=FU(ANG)
      ANG(2)=ANG(1)
      ANG(1)=0.0d0
      DO 2 I=1,3
2     P(2,I)=ANG(I)
      Y(2)=FU(ANG)
      ANG(3)=ANG(2)
      ANG(2)=0.0d0
      DO 3 I=1,3
3     P(3,I)=ANG(I)
      Y(3)=FU(ANG)
      ANG(1)=ANG(3)
      ANG(2)=ANG(3)
      ANG(3)=0.0d0
      DO 4 I=1,3
4     P(4,I)=ANG(I)
      Y(4)=FU(ANG)
      FTOL=0.0000001d0
      ITMAX=1500
      ITER=0
99999 ILO=1
      IF(Y(1).GT.Y(2))THEN
       IHI=1
       INHI=2
      ELSE
       IHI=2
       INHI=1
      ENDIF
      DO 5 I=1,4
      IF(Y(I).LT.Y(ILO))ILO=I
      IF(Y(I).GT.Y(IHI))THEN
       INHI=IHI
       IHI=I
      ELSE
       IF(Y(I).GT.Y(INHI))THEN
        IF(I.NE.IHI)INHI=I
       ENDIF
      ENDIF
5     CONTINUE
      IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.0.00000000001d0)THEN
       RTOL=0.0d0
      ELSE
       RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0
      ENDIF
      IF(RTOL.LT.FTOL)RETURN
      IF(ITER.EQ.ITMAX)THEN
       WRITE(3,*)' Rotation has not converged !'
       IERR=1
       RETURN
      ENDIF
      ITER=ITER+1
      DO 55 I=1,3
55    PBAR(I)=0.0d0
      DO 6 I=1,4
      IF(I.NE.IHI)THEN
       DO 7 J=1,3
7      PBAR(J)=PBAR(J)+P(I,J)
      ENDIF
6     CONTINUE
      DO 8 J=1,3
      PBAR(J)=PBAR(J)/3.0d0
8     PR(J)=2.0d0*PBAR(J)-P(IHI,J)
      YPR=FU(PR)
      IF(YPR.LE.Y(ILO))THEN
       DO 9 J=1,3
9      PRR(J)=2.0d0*PR(J)-PBAR(J)
       YPRR=FU(PRR)
       IF(YPRR.LT.Y(ILO))THEN
        DO 10 J=1,3
10      P(IHI,J)=PRR(J)
        Y(IHI)=YPRR
       ELSE
        DO 11 J=1,3
11      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ELSE
       IF(YPR.GE.Y(INHI))THEN
        IF(YPR.LT.Y(IHI))THEN
         DO 12 J=1,3
12       P(IHI,J)=PR(J)
         Y(IHI)=YPR
        ENDIF
        DO 13 J=1,3
13      PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J)
        YPRR=FU(PRR)
        IF(YPRR.LT.Y(IHI))THEN
         DO 14 J=1,3
14       P(IHI,J)=PRR(J)
         Y(IHI)=YPRR
        ELSE
         DO 15 I=1,4
         IF(I.NE.ILO)THEN
          DO 16 J=1,3
          PR(J)=0.5d0*(P(I,J)+P(ILO,J))
16        P(I,J)=PR(J)
          Y(I)=FU(PR)
         ENDIF
15       CONTINUE
        ENDIF
       ELSE
        DO 17 J=1,3
17      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ENDIF
      GOTO 99999
      END

      FUNCTION FU(ANG)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=15)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      DIMENSION ANG(3)
      S1=SIN(ANG(1))
      S2=SIN(ANG(2))
      S3=SIN(ANG(3))
      C1=COS(ANG(1))
      C2=COS(ANG(2))
      C3=COS(ANG(3))
C     ang1 = theta
C     ang2 = phi
C     ang3 = cappa
      A(1,1)=C1*C2*C3-S2*S3
      A(1,2)=C1*S2*C3+C2*S3
      A(1,3)=-S1*C3
      A(2,1)=-C1*C2*S3-S2*C3
      A(2,2)=-C1*S2*S3+C2*C3
      A(2,3)=S1*S3
      A(3,1)=S1*C2
      A(3,2)=S1*S2
      A(3,3)=C1
      R=0.0d0
      DO 1 I=1,NAT
      TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3)
      TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3)
      TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3)
1     R=R+(TX-XB(I,1))*(TX-XB(I,1))+
     1    (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3))
      FU=R
      RETURN
      END
c
      subroutine fw(n6,fi,nat,ity,r,iih,iin,iic,iio,lper,
     1DPCX,DPCY,DPCZ)
c     find waters and their potentials
      implicit none
      integer*4 i,n6,jj,ity(*),ityo,ityh,j1,nh,kk,k1,iih,iin,iic,iio,
     1ii,nat
      logical*4 lper
      real*8 fi(*),xo,yo,zo,r(*),xh,yh,zh,d2,x1,y1,z1,
     1x2,y2,z2,qh,qo,DPCX,DPCY,DPCZ,dx2,dy2,dz2
      data ityo,ityh,qh,qo/649,650,0.41d0,-0.82d0/

      dx2=DPCX/2.0d0
      dy2=DPCY/2.0d0
      dz2=DPCZ/2.0d0
      do 4 i=1,n6
4     fi(i)=0.0d0
      do 6 jj=1,nat
      if(ity(jj).eq.ityo)then
       j1=3*(jj-1)
       xo=r(j1+1)
       yo=r(j1+2)
       zo=r(j1+3)
       nh=0
       do 7 kk=1,nat
       if(ity(kk).eq.ityh)then
        k1=3*(kk-1)
        xh=r(k1+1)
        yh=r(k1+2)
        zh=r(k1+3)
        d2=(xo-xh)**2+(yo-yh)**2+(zo-zh)**2
        if(d2.gt.0.64d0.and.d2.lt.1.44d0)then
         nh=nh+1
         if(nh.gt.2)call report('Too many Hs around one O')
         if(nh.eq.1)then
          x1=xh
          y1=yh
          z1=zh
         else 
          x2=xh
          y2=yh
          z2=zh
         endif
        endif
       endif
7      continue
       if(nh.eq.2)then
        call add(iih ,1,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
        call add(iin ,2,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
        call add(iic ,3,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
        call add(iio ,4,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
c       call add(iian,5,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
c    1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
c       call add(iiac,6,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
c    1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
       endif 
      endif
6     continue
      write(6,*)' Potentials:'
      write(6,1603)(fi(ii),ii=1,n6)
1603  format(8f9.3)
      return
      end
c
      subroutine add(ii,ic,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
      implicit none
      real*8 r(*),xo,yo,zo,x1,y1,z1,x2,y2,z2,fi(*),x,y,z,ro,r1,r2,qh,qo,
     1dx2,dy2,dz2,DPCX,DPCY,DPCZ
      integer*4 ii,ic
      logical lper
c      
      x=r(3*(ii-1)+1)
      y=r(3*(ii-1)+2)
      z=r(3*(ii-1)+3)
      if(lper)then
       if(x-xo.gt.dx2.and.x-x1.gt.dx2.and.x-x2.gt.dx2)x=x-DPCX
       if(y-yo.gt.dy2.and.y-y1.gt.dy2.and.y-y2.gt.dy2)y=y-DPCY
       if(z-zo.gt.dz2.and.z-z1.gt.dz2.and.z-z2.gt.dz2)z=z-DPCZ
       if(xo-x.gt.dx2.and.x1-x.gt.dx2.and.x2-x.gt.dx2)x=x+DPCX
       if(yo-y.gt.dy2.and.y1-y.gt.dy2.and.y2-y.gt.dy2)y=y+DPCY
       if(zo-z.gt.dz2.and.z1-z.gt.dz2.and.z2-z.gt.dz2)z=z+DPCZ
      endif 
      ro=dsqrt((xo-x)**2+(yo-y)**2+(zo-z)**2)
      r1=dsqrt((x1-x)**2+(y1-y)**2+(z1-z)**2)
      r2=dsqrt((x2-x)**2+(y2-y)**2+(z2-z)**2)
c      
      fi(ic)=fi(ic)+qh/r1+qh/r2+qo/ro
      return
      end