program correcttdc
      implicit none
      integer*4 nat,N,i,ix,jx,ii,j,jj,n1,n2,no,io,it,itn
      real*8 bohr,d1,d2,d5,rij(3),uij,uir,ujr,epsr
      real*8,allocatable::r(:),p(:,:),f(:,:),fr(:,:),fn(:,:)
      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' FILE.TEN for atomic pairs not found in CCT.INP',/,/,
     1'      Input:     FILE.FC',/,
     1'                 FILE.TEN',/,
     1'                 FILE.X  ',/,
     1'                 CCT.INP, optional',/,
     1'                 (eps) as argument, optional',/,/,
     1'     Output:      NEW.FC',/,
     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),p(N,3),F(N,N),fr(N,N),fn(N,N),ind(N),ij(N,N))
      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,N
 3    r(i)=r(i)/bohr
      write(6,*)'Geometry read'

      close(9)
      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')

      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

      do 11 i=1,N
      do 11 j=1,N
11    fn(i,j)=F(i,j)

      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 TDC/APT:
       do 13 ix=1,3
13     rij(ix)=r(ix+3*(i-1))-r(ix+3*(j-1))
       d2=rij(1)**2+rij(2)**2+rij(3)**2
       d1=dsqrt(d2)
       d5=d2*d2*d1
       do 10 ix=1,3
       do 10 jx=1,3
       ii=ix+3*(i-1)
       jj=jx+3*(j-1)
       uij=p(ii,1)*p(jj,1)+p(ii,2)*p(jj,2)+p(ii,3)*p(jj,3)
       uir=p(ii,1)* rij(1)+p(ii,2)* rij(2)+p(ii,3)* rij(3)
       ujr=p(jj,1)* rij(1)+p(jj,2)* rij(2)+p(jj,3)* rij(3)
       fn(ii,jj)=fn(ii,jj)+(uij*d2-3.0d0*uir*ujr)/d5/epsr
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'

      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