program correctpolar
      implicit none
      integer*4 nat,N,i,ix,jx,ii,j,jj,a,b,c,k,nats,Ns,nw,iargc
      real*8 epsr,alpha,s
      real*8 ri(3),rj(3),rw(3),rwi(3),rwj(3),ti(3,3),tj(3,3),di,dj
      real*8,allocatable::r(:),p(:,:),f(:,:),rs(:)
      character*10 s10
      write(6,800)
800   format(/,' Corrects FILE.FC by polarizabilities of a solvent',/,
     1' FILE.TEN for atomic pairs not found in CCT.INP',/,/,
     1'      Input:     FILE.FC',/,
     1'                 FILE.TEN',/,
     1'                 FILE.X    geometry without solvent',/,
     1'                 FILES.X   geometry with solvent',/,
     1'     Output:      NEW.FC',/)

      if(iargc().ne.2)call report('Usage: correctpolar <eps> <a>')
      call getarg(1,s10)
      read(s10,*)epsr
      call getarg(2,s10)
      read(s10,*)alpha
      write(6,609)epsr,alpha
609   format(/,' Epsr = ',f6.2,' Alpha = ',f6.2,/)

      call rx(0,'FILE.X',nat,r)
      N=3*nat
      allocate(r(N),p(N,3),F(N,N))
      call rx(1,'FILE.X',nat,r)
      write(6,*)'Geometry read'

      call rx(0,'FILES.X',nats,rs)
      Ns=3*nats
      allocate(rs(Ns))
      call rx(1,'FILES.X',nats,rs)
      write(6,*)'Geometry with solvent read'

      open(9,file='FILE.TEN',status='old')
      read(9,*)
      do 2 i=1,N
2     read(9,*)(p(i,j),j=1,3)
      close(9)
      write(6,*)'APT read'

      CALL READFF(N,F,'FILE.FC')

      nw=(nats-nat)/3
      write(6,600)1,nat,nat+1,nats,nw
600   format(' Supposing atoms ',i6,' .. ',i6,' of solute',/,
     1       '                 ',i6,' .. ',i6,' of solvent',/,
     2       ' (',i6,' waters)')

      do 9 i=1,nat
      ri(1)=r(1+3*(i-1))
      ri(2)=r(2+3*(i-1))
      ri(3)=r(3+3*(i-1))
      do 9 j=i+1,nat
      rj(1)=r(1+3*(j-1))
      rj(2)=r(2+3*(j-1))
      rj(3)=r(3+3*(j-1))
      do 9 k=nat+1,nats,3
c     water coordinates:
      do 91 ix=1,3
      rw(ix)=(rs(ix+3*(k-1))+rs(ix+3*k)+rs(ix+3*(k+1)))/3.0d0
      rwi(ix)=rw(ix)-ri(ix)
91    rwj(ix)=rw(ix)-rj(ix)
      di=dsqrt(rwi(1)**2+rwi(2)**2+rwi(3)**2)
      dj=dsqrt(rwj(1)**2+rwj(2)**2+rwj(3)**2)
      do 92 ix=1,3
      do 92 jx=1,3
      ti(ix,jx)=3.0d0*rwi(ix)*rwi(jx)
      tj(ix,jx)=3.0d0*rwj(ix)*rwj(jx)
      if(ix.eq.jx)ti(ix,jx)=ti(ix,jx)-di**2
      if(ix.eq.jx)tj(ix,jx)=tj(ix,jx)-dj**2
      ti(ix,jx)=ti(ix,jx)/di**5
 92   tj(ix,jx)=tj(ix,jx)/dj**5

      do 9 ix=1,3
      do 9 jx=1,3
      ii=ix+3*(i-1)
      jj=jx+3*(j-1)
      s=0.0d0
      do 11 a=1,3
      do 11 b=1,3
      do 11 c=1,3
11    s=s+p(ii,a)*ti(a,b)*tj(b,c)*p(jj,c)
      s=s*2.0d0*alpha/epsr**2
      f(ii,jj)=f(ii,jj)+s
9     f(jj,ii)=f(ii,jj)

      OPEN(20,FILE='NEW.FC')
      CALL WRITEFF(N,N,f)
      CLOSE(20)
      WRITE(*,*)'FF written into NEW.FC'

      end

C
      SUBROUTINE READFF(N,FCAR,s)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      character*(*) s
      OPEN(20,FILE=s,STATUS='OLD')
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      CLOSE(20)
      write(6,*)s//' read '
      RETURN
      END

C
C
 
      SUBROUTINE WRITEFF(MX3,N,FCAR)
      IMPLICIT INTEGER*4 (I-N) 
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION FCAR(MX3,N)
C     CONST=4.359828/0.5291772**2
      CONST=1.0d0
      DO 6 I=1,N
      DO 6 J=1,N
6     FCAR(I,J)=FCAR(I,J)/CONST
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(I4,5D14.6)
      RETURN
      END

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      subroutine rx(ic,s,nat,r)
      implicit none
      integer*4 ic,nat,i,ix
      real*8 r(*),bohr
      character*(*) s
      bohr=0.529177d0
      open(9,file=s,status='old')
      read(9,*)
      read(9,*)nat
      if(ic.eq.1)then
       do 1 i=1,nat
1      read(9,*)r(1+3*(i-1)),(r(ix+3*(i-1)),ix=1,3)
       do 3 i=1,3*nat
 3     r(i)=r(i)/bohr
      endif
      close(9)
      return
      end