program correctchar
      implicit none
      integer*4 nat,N,i,ix,jx,ii,j,jj,k,n1,n2,no,io,it,itn,iargc
      real*8 bohr,epsr,rki(3),
     1dki,dki2,dki3,dki5,qiqk
      real*8,allocatable::r(:),f(:,:),fr(:,:),fn(:,:),q(:),
     1fq(:,:)
      integer*4,allocatable::ind(:),ij(:,:)
      logical lref,lnames,lcct
      character*10 s10
      bohr=0.529177d0
      write(6,800)
800   format(/,' Corrects FILE.FC by TDC interactions based on',/,
     1' partial atomic charges',/,/,
     1'      Input:     FILE.FC',/,
     1'                 FILE.X with the charges  ',/,
     1'                 CCT.INP, optional',/,
     1'                 (eps) as argument, optional',/,/,
     1'     Output:      NEW.FC new FF',/,
     1'                  CHAR.FC charge FF',/,
     1'                  REF.TXT (if REF.FC present)',/)

      if(iargc().eq.0)then
       epsr=1.0d0
      else
       call getarg(1,s10)
       read(s10,*)epsr
      endif
      write(6,609)epsr
609   format(/,' Epsr = ',f6.2,/)

      open(9,file='FILE.X',status='old')
      read(9,*)
      read(9,*)nat
      N=3*nat
      allocate(r(N),F(N,N),fr(N,N),fn(N,N),ind(N),ij(N,N),
     1fq(N,N),q(nat))
      do 1 i=1,nat
1     read(9,*)r(1+3*(i-1)),(r(ix+3*(i-1)),ix=1,3),(q(i),ix=1,8)
      do 3 i=1,N
 3    r(i)=r(i)/bohr
      write(6,*)'Geometry read'
      close(9)

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

      do 8 i=1,nat
      do 8 j=1,nat
8     ij(i,j)=0

      no=0
      lnames=.false.
      inquire(file='CCT.INP',exist=lcct)
      if(lcct)then
       open(9,file='CCT.INP',status='old')
4      read(9,90,end=99,err=99)s10
90     format(a10)
       if(s10(1:6).eq.'LNAMES')read(9,*)lnames
       if(s10(1:7).eq.'POLYMER')then
        read(9,*)no
        do 5 io=1,no
        read(9,*)
        if(lnames)read(9,*)
        n1=1
6       n2=min(n1+19,nat)
        read(9,220)(ind(i),i=n1,n2)
 220    format(20i3)
        if(n2.lt.nat)then
         n1=n2+1
         goto 6
        endif
        n1=1
7       n2=min(n1+19,nat)
        read(9,220)(ind(i),i=n1,n2)
        if(n2.lt.nat)then
         n1=n2+1
         goto 7
        endif
        do 5 i=1,nat
        do 5 j=1,nat
5       if(ind(i).ne.0.and.ind(j).ne.0)ij(i,j)=ij(i,j)+1
        goto 99
       endif
       goto 4
99     close(9)
       write(6,*)'CCT.INP ',no,' overlaps'
      endif

      inquire(file='REF.FC',exist=lref)
      if(lref)then
       CALL READFF(N,fr,'REF.FC')
       open(20,file='REF.TXT')
      endif

      fn=F

      fq=0.0d0
      do 91 k=1,nat
      do 91 i=1,nat
      if(i.ne.k)then
       qiqk=q(i)*q(k)
       do 113 ix=1,3
113    rki(ix)=r(ix+3*(k-1))-r(ix+3*(i-1))
       dki2=rki(1)**2+rki(2)**2+rki(3)**2
       dki=dsqrt(dki2)
       dki3=dki2*dki
       dki5=dki3*dki2
       do 92 ix=1,3
       ii=ix+3*(k-1)
c      d^2E/drk^2:
       fq(ii,ii)=fq(ii,ii)-qiqk/dki3
c      d^2E/drk/dri(i<>k):
       jj=ix+3*(i-1)
       fq(ii,jj)=qiqk/dki3
       do 92 jx=1,3
       jj=jx+3*(k-1)
c      d^2E/drk^2:
       fq(ii,jj)=fq(ii,jj)+3.0d0*qiqk*rki(ix)*rki(jx)/dki5
c      d^2E/drk/dri(i<>k):
       jj=jx+3*(i-1)
92     fq(ii,jj)=-3.0d0*qiqk*rki(ix)*rki(jx)/dki5
      endif
91    continue

      fq=fq/epsr

      it=0
      itn=0
      do 9 i=1,nat
      do 9 j=i+1,nat
      it=it+1
      if((lcct.and.ij(i,j).eq.0).or.(.not.lcct.and.
     1 dabs(f(i,j)).lt.1.0d-10))then
       itn=itn+1
c      Fij was not defined or is zero, apply correction:
       do 10 ix=1,3
       do 10 jx=1,3
       ii=ix+3*(i-1)
       jj=jx+3*(j-1)
       fn(ii,jj)=f(ii,j)+fq(ii,jj)
10     fn(jj,ii)=fn(ii,jj)
       if(lref)then
        do 12 ix=1,3
        do 12 jx=1,3
        ii=ix+3*(i-1)
        jj=jx+3*(j-1)
12      write(20,200)i,j,ix,jx,f(ii,jj),fn(ii,jj),fr(ii,jj)
200     format(4i4,3g12.4)
       endif
      endif
9     continue
      write(6,*)itn,' pairs of ',it,' corrected'

      if(lref)then
       close(20)
       write(6,*)'REF.TXT written'
      endif

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

      OPEN(20,FILE='CHAR.FC')
      CALL WRITEFF(N,N,fq)
      CLOSE(20)
      WRITE(*,*)'Charge FF written into CHAR.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