PROGRAM CCTG
      IMPLICIT none
      integer*4 NAT,NO,NNO,IWG,NATS_max,NE,I,ia,ix,NATS
      real*8 CUTOFF,d,TE
      logical loldformat,lnew,LHALTONERROR,LSTRICT,LNAMES,LWR
      CHARACTER*80 basicname,xname
      real*8,allocatable::R(:,:),RS(:,:),rname(:,:,:),Re(:,:)
      integer*4,allocatable::ND(:),IBONDSB(:,:),NBOB(:),INDBIG(:,:),
     1NDS(:),NDNAME(:,:),NATNAME(:)
      CHARACTER*80,allocatable:: NAME(:)
c
      WRITE(6,3000)
3000  FORMAT(/,
     1       ' Transfer of Molecular Tensors in Cartesian Coordinates:'
     8,/,    ' adjusting interatomic distances according to fragments'
     8,/,/,  '           By Petr Bour, Prague 2021',/,/,
     9       '    INPUT FILES:',/,
     9       '      BIG.X coordinates',/,
     3       '      SMALL.X coordinates (or named .x)',/,
     4       '      CTT.INP  atom assignment',/,/,
     5       '    OUTPUT:',/,
     5       '      BIGN.X repaired',/)

c     Input geometry of the target molecule:
      OPEN(2,FILE='BIG.X',STATUS='OLD')
      read(2,*)
      read(2,*)NAT
      allocate(R(3,NAT),ND(NAT),IBONDSB(4,NAT),NBOB(NAT),Re(3,NAT))
      rewind 2
      CALL READRX(NAT,R,ND,IBONDSB,NBOB)
      CLOSE(2)
      WRITE(6,*)NAT, ' atoms in the big molecule'

      NO=NNO()
      if(NO.eq.0)call report('no overlaps')
      allocate(INDBIG(NO,NAT),NAME(NO),NATNAME(NO))
c
c     Open file with option and overlaps and read it in:
      CALL readopt(NAT,NO,INDBIG,NAME,loldformat,lnew,
     1LHALTONERROR,LSTRICT,LNAMES,IWG,CUTOFF,LWR,d,TE)
c
c     if tensors constructed from named fragments, load them all:
      if(LNAMES)then
       NATS_max=0
       DO 111 I=1,NO 
       call getname(NAME(I),basicname,NE)
       xname=basicname(1:NE)//'.x'
       open(2,file=xname,status='old')
       read(2,*)
       read(2,*)NATS
       CLOSE(2)
       write(6,*)' Fragment ',basicname(1:NE),': ',NATS,' atoms'
       NATNAME(I)=NATS
111    if(NATS.gt.NATS_max)NATS_max=NATS
       allocate(RS(3,NATS_max),NDS(NATS_max),NDNAME(NO,NATS_max),
     1 rname(NO,3,NATS_max))

       DO 11 I=1,NO 
       call getname(NAME(I),basicname,NE)
       xname=basicname(1:NE)//'.x'
       open(2,file=xname,status='old')
       CALL READRXS(NATS_max,RS,NDS)
       CLOSE(2)
       do 11 ia=1,NATS
       NDNAME(I,ia)=NDS(ia)
       do 11 ix=1,3
11     rname(I,ix,ia)=RS(ix,ia)
      else
       OPEN(2,FILE='SMALL.X',STATUS='OLD')
       read(2,*)
       read(2,*)NATS
       CLOSE(2)
       WRITE(6,*)NATS, ' atoms in SMALL.X'
       NATS_max=NATS
       allocate(RS(3,NATS),NDS(NATS),NDNAME(1,NATS),rname(1,3,NATS))
       OPEN(2,FILE='SMALL.X',STATUS='OLD')
       CALL READRXS(NATS,RS,NDS)
       CLOSE(2)
       NATNAME(1)=NATS
       do 131 ia=1,NATS
       NDNAME(1,ia)=NDS(ia)
       do 131 ix=1,3
131    rname(1,ix,ia)=RS(ix,ia)
      endif
c     names
c
c     check consistency in atomic types:
      do 22 i=1,NO
      if(LNAMES)then
       NATS=natname(i)
       do 23 ia=1,NATS
23     NDS(ia)=ndname(i,ia)
      endif
      DO 22 ia=1,NAT
      IF(INDBIG(I,ia).NE.0)THEN
       IF(ND(ia).NE.NDS(ABS(INDBIG(I,ia))))THEN
        WRITE(6,3300)I,ia,ND(ia),NDS(ABS(INDBIG(I,ia)))
3300    FORMAT('Overlap ',I2,', atom ',I3,': ',2I4)
        CALL REPORT(' Different atomic types !')
       ENDIF
       IF(INDBIG(I,ia).gt.NATS)then
        write(6,*)'overlap ',i
        write(6,*)' NATS   ',NATS
        write(6,*)'IND     ',ia
        write(6,*)'INDB    ',INDBIG(I,ia)
        call report('Overflow NATS')
       ENDIF
      ENDIF
22    CONTINUE

      CALL IMPROVE(LSTRICT,IWG,LNAMES,LHALTONERROR,Re,
     1cutoff,NAT,R,LWR,NATS_max,RS,NO,rname,INDBIG,d)

      CALL WRITERX(NAT,Re,ND,IBONDSB,'BIGN.X')
      END

c
      SUBROUTINE IMPROVE(LSTRICT,IWG,LNAMES,LHALTONERROR,Re,
     1cutoff,NAT,RB,LWR,NATS_max,RS,NO,rname,INDBIG,d)
      implicit none
      logical LSTRICT,LNAMES,LHALTONERROR,LWR
      integer*4 IWG,NAT,IA,NF,IX,IO,NATS_max,NO,I,IAS,JAS,II,IIOII,
     1INDBIG(NO,NAT),IOIJ,IP,IPO,IPOOLD,ISUM,JA,NATIO
      real*8 cutoff,RB(3,NAT),xi,yi,zi,TIJ(3),DISTM,OTIJ(3),
     1RS(3,NATS_max),rname(NO,3,NATS_max),DIST,SUM,RMS,Re(3,NAT),RMSF,
     1d,di
      integer*4,allocatable::IOL(:)
      real*8,allocatable::WF(:),Rq(:,:,:)
      IF(LSTRICT)     WRITE(6,*)' Only transformed atoms in the beds'
      IF(IWG.EQ.1)    WRITE(6,*)' Multiple overlaps will be weighted'
      IF(LNAMES)      WRITE(6,*)' Named small fragments'
      IF(LHALTONERROR)WRITE(6,*)' Will stop on error'
      if(IWG.eq.0)    write(6,*)' IWG=0 best overlap only'
      if(IWG.eq.1)    write(6,*)' IWG=1 weighing of overlaps'
      WRITE(6,3999)CUTOFF
3999  format(' Distance cutoff ',f12.2,' A')
      IOIJ=0
      IPOOLD=0
      allocate(IOL(NO),WF(NO),Rq(NAT,NAT,3))

      DO 27 IA= 1,NAT
      DO 27 JA= 1,NAT
      Rq(IA,JA,1)=RB(1,IA)-RB(1,JA)
      Rq(IA,JA,2)=RB(2,IA)-RB(2,JA)
 27   Rq(IA,JA,3)=RB(3,IA)-RB(3,JA)
c
C     Look at all atom pairs i-j:
      DO 20 IA= 1,NAT
      DO 20 JA=IA+1,NAT
      di=dsqrt(Rq(IA,JA,1)**2+Rq(IA,JA,2)**2+Rq(IA,JA,3)**2)
      if(CUTOFF.gt.0.0d0.and.di.gt.CUTOFF)goto 20
C
C     Calculate center of i/j radiusvectors:
      DO 22 IX=1,3
22    TIJ(IX)=0.5d0*(RB(IX,IA)+RB(IX,JA))
      DISTM=1.0d5
      NF=0

C     Look at all overlaps io, find that which fits best the i-j pair:
      DO 21 IO=1,NO
C
C     Check, if io contains ij:
      IF(INDBIG(IO,IA).GT.0.AND.INDBIG(IO,JA).GT.0)THEN
       NF=NF+1
       IOL(NF)=IO
C
C      Calculate center of the overlap and distance from i-j center
c             i.........TIJ.............j
c      ----------OTIJ---------------------
       DO 23 IX=1,3
23     OTIJ(IX)=0.0D0
       NATIO=0
       DO 24 II=1,NAT
       IIOII=INDBIG(IO,II)
C      If Lstrict, take only transformed atoms
       IF(LSTRICT.AND.IIOII.LT.0)IIOII=0
       IIOII=ABS(IIOII)
       IF(IIOII.NE.0)THEN
        NATIO=NATIO+1
        DO 25 IX=1,3
25      OTIJ(IX)=OTIJ(IX)+RB(IX,II)
       ENDIF
24     CONTINUE
c
       DO 26 IX=1,3
26     OTIJ(IX)=OTIJ(IX)/DBLE(NATIO)
       DIST=(OTIJ(1)-TIJ(1))**2+(OTIJ(2)-TIJ(2))**2+(OTIJ(3)-TIJ(3))**2
       IF(DIST.LT.DISTM)THEN
        DISTM=DIST
        IOIJ=IO
       ENDIF
       WF(NF)=DSQRT(DIST)
      ENDIF
21    CONTINUE
c
      if(LWR)then
       if(NF.EQ.0)then
        WRITE(6,3003)IA,JA,0,0
       else
        WRITE(6,3003)IA,JA,NF,IOIJ
3003    FORMAT(3I5,' overlaps found; best is # ',I4)
       endif
      endif
C 
C     Calculate the weighting factors:
      SUM=0.0d0
      DO 43 I=1,NF
43    IF(WF(I).GT.1.0d-9)SUM=SUM+1.0d0/WF(I)
      DO 42 I=1,NF
      IF(WF(I).LE.1.0d-9)THEN
       WF(1)=1.0d0
       NF=1
       GOTO 44
      ENDIF     
      WF(I)=1.0d0/WF(I)/SUM
42    CONTINUE
44    CONTINUE
C
c
C     Loop over possible overlaps:
      ISUM=0
c     IPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIP
      DO 45 IP=1,NF
      IPO=IOL(IP)
      IF(IWG.EQ.0)THEN
       IF(IPO.NE.IOIJ)GOTO 45
       WF(IP)=1.0d0
      ENDIF
      ISUM=ISUM+1
C
      if(LNAMES)then
       if(IPOOLD.EQ.0.or.IPOOLD.NE.IPO)then
        IPOOLD=IPO
       endif
      ENDIF

C     New distance
C     Zero out at the start of the summation
      IF(ISUM.EQ.1)then
       Rq(IA,JA,1)=0.0d0
       Rq(IA,JA,2)=0.0d0
       Rq(IA,JA,3)=0.0d0
      endif
      IAS=INDBIG(IPO,IA)
      JAS=INDBIG(IPO,JA)
      if(LNAMES)then
       xi=rname(IPO,1,IAS)-rname(IPO,1,JAS)
       yi=rname(IPO,2,IAS)-rname(IPO,2,JAS)
       zi=rname(IPO,3,IAS)-rname(IPO,3,JAS)
      else
       xi=RS(1,IAS)-RS(1,JAS)
       yi=RS(2,IAS)-RS(2,JAS)
       zi=RS(3,IAS)-RS(3,JAS)
      endif
      Rq(IA,JA,1)=Rq(IA,JA,1)+WF(IP)*xi
      Rq(IA,JA,2)=Rq(IA,JA,2)+WF(IP)*yi
      Rq(IA,JA,3)=Rq(IA,JA,3)+WF(IP)*zi
      Rq(JA,IA,1)=-Rq(IA,JA,1)
      Rq(JA,IA,2)=-Rq(IA,JA,2)
      Rq(JA,IA,3)=-Rq(IA,JA,3)
45    CONTINUE
c     IPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIP
C
20    CONTINUE
c     IA,JA

      RMS=RMSF(RB,Rq,NAT,d)
      write(6,601)RMS
601   format(' RMS ',f12.3,' A')

      Re=RB
      call geg(NAT,d,RB,Re,Rq)
      RMS=RMSF(Re,Rq,NAT,d)
      write(6,601)RMS

      RETURN
      END

      subroutine geg(NAT,d,R,Re,Rq)
      implicit none
      integer*4 NAT,ka,kx,ja
      real*8 R(3,NAT),Re(3,NAT),fg,d,Rq(NAT,NAT,3),s1,s2
      do 1 ka=1,NAT
      do 1 kx=1,3
      s1=0.0d0
      s2=0.0d0
      do 2 ja=1,NAT
      if(ja.ne.ka)then
       s1=s1+(R(kx,ja)+Rq(ka,ja,kx))*fg(Rq(ka,ja,kx),d)
       s2=s2+fg(Rq(ka,ja,kx),d)
      endif
2     continue
1     Re(kx,ka)=s1/s2
      return
      end


      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')
      CLOSE(3)
      CLOSE(2)
      STOP
      END
c

c
      function NNO()
      implicit none
      integer*4 NNO,NO
      CHARACTER*7 STR
      NO=0
      OPEN(2,FILE='CCT.INP',STATUS='OLD')
9     READ(2,2000,end=99,err=99)STR
2000  FORMAT(A7)
      IF(STR.EQ.'POLYMER')READ(2,*)NO
      goto 9
   99 close(2)
      NNO=NO
      return
      end

      SUBROUTINE readopt(NAT,NO,INDBIG,NAME,loldformat,lnew,
     1LHALTONERROR,LSTRICT,LNAMES,IWG,CUTOFF,LWR,d,TE)
      implicit none
      INTEGER*4 NAT,NO,INDBIG(NO,NAT),IWG,IAB,I,IC,nt,IO,
     1ia,ja,ismall
      CHARACTER*7 STR
      CHARACTER*80 NAME(NO)
      LOGICAL loldformat,lnew,LHALTONERROR,LSTRICT,LNAMES,LWR
      real*8 CUTOFF,d,TE
      INTEGER*4 ,allocatable::bb(:,:)
      loldformat=.true.
      lnew=.false.
      LHALTONERROR=.TRUE.
      CUTOFF=-1.0d5
      LSTRICT=.TRUE.
      LWR=.FALSE.
      LNAMES=.FALSE.
      IWG=1
      d=3.0d0
      TE=1.0d-3
      OPEN(2,FILE='CCT.INP',STATUS='OLD')
9     READ(2,2000)STR
2000  FORMAT(A7)
      WRITE(6,2000)STR
c     20i3 format:
      IF(STR(1:7).EQ.'OLDFORM')READ(2,*)loldformat
c     IWG=0 best overlap, IWG=1 weighting:
      IF(STR(1:3).EQ.'IWG')    READ(2,*)IWG
c     named fragments:
      IF(STR(1:6).EQ.'LNAMES') READ(2,*)LNAMES
c     transfor atoms in beds only:
      IF(STR.EQ.'LSTRICT')     READ(2,*)LSTRICT
c     extensive output:
      IF(STR(1:5).EQ.'WRITE')  READ(2,*)LWR
c     New format of overlap map:
      IF(STR(1:4).EQ.'LNEW')   READ(2,*)lnew
c     Stop transfer if an error in atom maps occurs:
      IF(STR(1:7).EQ.'LHALTON')READ(2,*)LHALTONERROR
c     Cutoff for atomic distances considered:
      IF(STR(1:6).EQ.'CUTOFF') READ(2,*)CUTOFF
      IF(STR(1:5).EQ.'WIDTH') READ(2,*)d
      IF(STR(1:5).EQ.'TOLER') READ(2,*)TE
c     In polarization correction, use distances from 
      IF(STR.EQ.'POLYMER')THEN
       READ(2,*)NO
       DO 2 I=1,NO
       READ(2,*)IC
       IF(IC.NE.I)THEN
        WRITE(6,*)IC,I
        CALL REPORT('Error in overlap reading ')
       ENDIF
c      named fragments:
       IF(LNAMES)READ(2,2900)NAME(I)
2900   FORMAT(a80)
c      newest format
       if(lnew)then
        DO 52 IAB=1,NAT
52      INDBIG(I,IAB)=0
        read(2,*)nt
        allocate(bb(2,nt))
        backspace 2
        read(2,*)nt,((bb(IO,IAB),IO=1,2),IAB=1,nt)
        do 53 IAB=1,nt
53      INDBIG(I,bb(1,IAB))=bb(2,IAB)
        deallocate(bb)
       else
        if(loldformat)then
         READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT)
2001     FORMAT(20I3)
         DO 5 IAB=1,NAT
5        INDBIG(I,IAB)=0
         READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT)
        else
         READ(2,*)(INDBIG(I,IAB),IAB=1,NAT)
         DO 51 IAB=1,NAT
51       INDBIG(I,IAB)=0
         READ(2,*)(INDBIG(I,IAB),IAB=1,NAT)
        endif
       endif
2      CONTINUE
c
c      check, that atoms are not defined twice in the overlaps:
       do 1 io=1,no
       do 1 ia=1,nat
       ismall=indbig(io,ia)
       if(ismall.ne.0)then
        do 3 ja=ia+1,nat
        if(ismall.eq.indbig(io,ja))then
         write(6,*)'overlap ',io
         write(6,*)'big     ',ia,ja
         write(6,*)'small   ',ismall
         call report('umbiqous overlap')
        endif
 3      continue
       endif
 1     continue 
       
       goto 9999
      ENDIF

      goto 9
c
c
9999  close(2)
      RETURN
      END
c     ============================================================
      function fg(RN,d)
      real*8 fg,RN,d,e
      e=(RN/d)**2
      if(e.gt.20.0d0)then
       fg=0.0d0
      else
       fg=exp(-e)
      endif
      return
      end
c
      SUBROUTINE READRX(N,R,ND,IBONDS,NB)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(3,N),ND(N),IBONDS(4,N),NB(N)
      READ(2,*)
      read(2,*)N
      do 1 i=1,N
      READ(2,*)ND(i),R(1,i),R(2,i),R(3,i),(IBONDS(IB,i),IB=1,4)
      NB(i)=0
      DO 1 IB=1,4
1     IF(IBONDS(IB,i).NE.0)NB(i)=IB
      RETURN
      END

      SUBROUTINE READRXS(N0,R,ND)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(3,N0),ND(N0)
      READ(2,*)
      read(2,*)N
      do 1 i=1,N
1     READ(2,*)ND(i),R(1,i),R(2,i),R(3,i)
      RETURN
      END

      SUBROUTINE WRITERX(N,R,ND,IB,s)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(3,N),ND(N),IB(4,N)
      character*(*) s
      OPEN(2,FILE=s)
      WRITE(2,2000)
2000  FORMAT('FORCE FIELD IMPROVEMENT')
      write(2,2002)N
2002  format(I6)
      CHARGE=0.0d0
      DO 1 I=1,N
1     WRITE(2,2001)ND(I),R(1,I),R(2,I),R(3,I),(IB(J,I),J=1,4),
     10,0,0,CHARGE
2001  FORMAT(I5,3F12.6,7I5,F8.4)
      CLOSE(2)
      WRITE(6,*)
      WRITE(6,*)s//' written '
      RETURN
      END
      subroutine getname(NAME,basicname,NE)
      implicit none
      integer*4 NE
      character*80 NAME,basicname
      NE=80
      do 10 NE=80,1,-1
10    if(NAME(NE:NE).ne.' ')goto 12
12    basicname(1:NE)=NAME(1:NE)
      return
      end
c     ============================================================

      function RMSF(R,RQ,NAT,d)
      implicit none
      integer*4 NAT,IA,JA
      real*8 R(3,NAT),RQ(NAT,NAT,3),RMSF,RMS,d,fg,
     1xi,yi,zi,xj,yj,zj,RN
      RMS=0.0d0
      DO 2 IA= 1,NAT
      xi=R(1,IA)
      yi=R(2,IA)
      zi=R(3,IA)
      DO 2 JA=IA+1,NAT
      xj=xi-R(1,JA)-Rq(IA,JA,1)
      yj=yi-R(2,JA)-Rq(IA,JA,2)
      zj=zi-R(3,JA)-Rq(IA,JA,3)
      RN=dsqrt(Rq(IA,JA,1)**2+Rq(IA,JA,2)**2+Rq(IA,JA,3)**2)
2     RMS=RMS+(xj**2+yj**2+zj**2)*fg(RN,d)
      RMSF=RMS
      return
      end