PROGRAM CCTF
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MX=600,MXB=10,IOX=20)
      DIMENSION FF(3*MX,3*MX),R(3,MX),IND(MX),IB(MX,MXB),ND(MX),
     1NDS(MX),RS(3,MX),FFS(3*MX,3*MX),IS(MX,MXB),AM(MX,3,3),
     2INDBIG(IOX,MX),IBONDS(4,MX),IBONDSB(4,MX),
     3NBO(MX),NBOB(MX),
     6zq(MX),zqs(MX),e(3*MX),es(3*MX)
      LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT,
     1LAPT,LALF,LRDO
c
      OPEN(3,FILE='CCTF.OUT')
      WRITE(3,3000)
      WRITE(6,3000)
3000  FORMAT(/,
     1       ' Compares 2 S-matrices for similar molecules'
     8,/,/,  '           By Petr Bour, Prague 2000',/,/,
     9       ' INPUT FILES:    BIG.X    coordinates of the target',
     1/,     '                BIGF.INP  S-matrix of the target',/,
     3       '               SMALL.X    coordinates of the source',/,
     2       '              SMALLF.INP  S-matrix of the source',/,
     4       '                CTTF.INP  overlap atom assignment',/,/,
     5       ' OUTPUT     :   CTTF.OUT  checking output file',/)
c
      OPEN(2,FILE='BIG.X',STATUS='OLD')
      CALL READRX(MX,NAT,R,ND,IBONDSB,NBOB)
      CLOSE(2)
      WRITE(3,*)NAT, ' atoms in the big molecule'
      WRITE(6,*)NAT, ' atoms in the big molecule'
      OPEN(2,FILE='SMALL.X',STATUS='OLD')
      CALL READRX(MX,NATS,RS,NDS,IBONDS,NBO)
      CLOSE(2)
      WRITE(3,*)NATS, ' atoms in the small molecule'
      WRITE(6,*)NATS, ' atoms in the small molecule'
c
      OPEN(2,FILE='CCTF.INP',STATUS='OLD')
      CALL INPUTIN(NS,IB,IS,IND,MX,MXB,INDBIG,IOX,NO,NAT,ND,NDS,
     1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LSTRICT,LINV,LSELECT,LAPT,
     2LALF,LRDO)
      CLOSE(2)
c
      OPEN(2,FILE='BIGF.INP',STATUS='OLD')
      call READS(3*MX,FF,nat,r,zq,ni,e)
      close(2)
      WRITE(3,*)NAT, ' atoms in big S'
      WRITE(6,*)NAT, ' atoms in big S'
      OPEN(2,FILE='SMALLF.INP',STATUS='OLD')
      call READS(3*MX,FFS,nats,rs,zqs,nis,es)
      WRITE(6,*)NATs, ' atoms in small S'
      WRITE(3,*)NATs, ' atoms in small S'
      close(2)
C
c     loop over modes in small molecule:
      write(6,6002)
      write(3,6002)
6002  format('mode  small   best fit in big     agreement',/,80(1h-))
      do 600 i=1,nis
c
c     look for all big modes:
      ibest=1
      i2best=1
      abest=1.0d20
      a2best=1.0d20
      abestn=1.0d20
      a2bestn=1.0d20
      do 601 iqb=1,ni
c
c     try plus and minus phase
      CALL IMPROVE(FF,FFS,R,RS,NAT,NATS,NS,IND,IB,IS,
     1AM,INDBIG,NO,IBONDSB,NBOB,LWR,LABS,LVCD,LOFF,LDIA,
     2IWG,LROA,LSTRICT,LINV,LSELECT,iqb,an,i,1.0d0,ans,anb)
      CALL IMPROVE(FF,FFS,R,RS,NAT,NATS,NS,IND,IB,IS,
     1AM,INDBIG,NO,IBONDSB,NBOB,LWR,LABS,LVCD,LOFF,LDIA,
     2IWG,LROA,LSTRICT,LINV,LSELECT,iqb,anminus,i,-1.0d0,ans,anb)
      if(anminus.lt.an)an=anminus
      if(an.lt.abest)then
       i2best=ibest
       a2best=abest
       a2bestn=abestn
       abest=an
       ibest=iqb
       abestn=anb
      endif
      if(an.lt.a2best.and.iqb.ne.ibest)then
       a2best=an
       i2best=iqb
       a2bestn=anb
      endif
601   continue
c
      proc=100.0d0*(1.0d0-abest/(ans+abestn)/2.0d0)
      proc2=100.0d0*(1.0d0-a2best/(ans+a2bestn)/2.0d0)
      write(3,6001)i,es(i),ibest, e(ibest), proc,
     1                     i2best,e(i2best),proc2,ans,abestn,a2bestn
600   write(6,6001)i,es(i),ibest ,e(ibest) ,proc,
     1                     i2best,e(i2best),proc2
6001  format(i4,f8.1,i5,f8.1,f6.2,'%',i5,f8.1,f6.2,'%',3f9.2)
      CLOSE(3)
      WRITE(6,*)' Program finished OK.'
      STOP
      END

      SUBROUTINE IMPROVE(FFB,FFS,RB,RS,NAT,NATS,NS,IND,IB,IS,AM,INDBIG,
     1NO,IBONDS,NBOB,LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,
     2LSTRICT,LINV,LSELECT,iqb,an,iq,phase,ans,anb)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MX=600,MXB=10,IOX=20)
      DIMENSION FFB(3*MX,3*MX),FFS(3*MX,3*MX),RB(3,NAT),RS(3,NATS),
     1IND(MX),TS(3),TB(3),IS(MX,MXB),IB(MX,MXB),AM(MX,3,3),
     3B(3,3),INDBIG(IOX,MX),TIJ(3),OTIJ(3),
     4IBONDS(4,MX),NBOB(MX),
     5IOL(IOX),WF(IOX)
      LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT
C     Takes big force field FFB and implants some constants from FFS
      IF(NO.EQ.0)call report('NO=0')
      IF(LSTRICT.and.lwr)WRITE(3,*)' Only transformed atoms in the beds'
      IWG=0
      IF(LINV.and.lwr)WRITE(3,*)' Smaller fragment as optical antipode' 
C
      an=0.0d0
      anb=0.0d0
      ans=0.0d0
c
      LOFF=.false.
C     Look at all atom pairs i-j:
      DO 20 IA= 1,NAT
C
C     Calculate center of i/j radiusvectors:
      DO 22 IX=1,3
22    TIJ(IX)=RB(IX,IA)
C
C     Look at all overlaps io, find that which fits best the i-j pair:
      DISTM=100000.0d0
      NF=0
      DO 21 IO=1,NO
C
C     Check, if io contains i:
      IF(INDBIG(IO,IA).GT.0)THEN
       NF=NF+1
       IOL(NF)=IO
C
C      Calculate center of the overlap
       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
       DO 26 IX=1,3
26     OTIJ(IX)=OTIJ(IX)/REAL(NATIO)
C      Calculate distance from i-j center
       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)=SQRT(DIST)
      ENDIF
21    CONTINUE
      IF(NF.EQ.0)THEN
3002  FORMAT(' Pair ',2I4,': no overlap found')
       GOTO 20
      ELSE
       if(LWR)WRITE(3,3003)NF,IA,IA,IOIJ
3003   FORMAT(I5,' overlaps found for pair ',2I4,'; best is # ',I4)
C 
C      Calculate the weighting factors:
       SUM=0.0d0
       DO 43 I=1,NF
43     IF(WF(I).GT.0.000000001d0)SUM=SUM+1/WF(I)
       DO 42 I=1,NF
       IF(WF(I).LE.0.000000001)THEN
        WF(1)=1.0d0
        NF=1
        GOTO 44
       ENDIF     
       WF(I)=1.0d0/WF(I)/SUM
42     CONTINUE
44     CONTINUE
C
C      Loop over possible overlaps:
       ISUM=0
       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

C
       IPASS=1
C      Try two passes, firstly only connected atoms taken:
666    ITU=0
       IF(IPASS.EQ.2.and.lwr) WRITE(3,*)' Second pass ... '
C
       DO 27 II=1,NAT
       IF(INDBIG(IPO,II).EQ.0)GOTO 27
C
C      Take IA and 
       IF(II.EQ.IA)GOTO 277
       DO 29 IBA=1,NBOB(II)
C
C      Take atoms connected to IA or 
       IF(IBONDS(IBA,II).EQ.IA)GOTO 277
29     CONTINUE
C
       IF(IPASS.EQ.2)THEN
        DO 31 IBA=1,NBOB(II)
        IIA=IBONDS(IBA,II)
        DO 31 JB=1,NBOB(IIA)
        IF(IBONDS(JB,IIA).EQ.IA)GOTO 277
31      CONTINUE
       ENDIF
       GOTO 27
C
C      If Lstrict, don't take non-transformed atoms:
277    IF(INDBIG(IPO,II).LT.0.AND.LSTRICT)GOTO 27
       ITU=ITU+1
       DO 28 IX=1,3
       XS(ITU,IX)=RS(IX,ABS(INDBIG(IPO,II)))
28     XB(ITU,IX)=RB(IX,II)
       IND(ITU)=II
27     CONTINUE
       IF(ITU.LT.3.and.LWR)WRITE(3,*)' ITU < 3 !'
       IF(ITU.LT.3.AND.IPASS.EQ.1)THEN
        IPASS=2
        GOTO 666
       ENDIF 
       IF(ITU.LT.3)CALL REPORT(' ITU < 3 !')
       IF(ITU.GT.MXB)CALL REPORT('ITU > MXB !')
       if(LWR)WRITE(3,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU)
       if(LWR)WRITE(3,3005)(INDBIG(IPO,IND(I)),I=1,ITU)
3004   FORMAT('Overlap ',I2,': weighting factor: ',F7.4,';',
     1         I4,' atom bed:',/,
     2            '      Big  : ',20I3)
3005   FORMAT(    '      Small: ',20I3,/)
C     Calculate center coordinates to be rotated:
       DO 33 I=1,3
       TB(I) =0.0d0
33     TS(I)=0.0d0
       DO 34 ISAT=1,ITU
       DO 34 J=1,3
       TB(J)=TB(J)+XB(ISAT,J)
34     TS(J)=TS(J)+XS(ISAT,J)
       DO 38 J=1,3
       TB(J)=TB(J)/REAL(ITU)
38     TS(J)=TS(J)/REAL(ITU)
C      Subtract mass-center:
       DO 35 ISAT=1,ITU
       DO 35 J=1,3
       XS(ISAT,J)=XS(ISAT,J)-TS(J)
35     XB(ISAT,J)=XB(ISAT,J)-TB(J)
C
       IF(LINV)THEN
        DO 39 ISAT=1,ITU
        DO 39 J=1,3
39      XS(ISAT,J)=-XS(ISAT,J)
       ENDIF
C
C      Calculate the transformation matrix:
       CALL DOMATRIX(IERR)
       IF(IERR.EQ.1.AND.IPASS.EQ.1)THEN
        IPASS=2
        GOTO 666
       ENDIF
       IF(IERR.EQ.1)CALL REPORT('Program cannot continue')
C
       IF(LINV)THEN
        DO 50 I=1,3
        DO 50 J=1,3
50      A(I,J)=-A(I,J)
       ENDIF
C
C      Rotate small matrix and compare to big
       IAS=INDBIG(IPO,IA)
       DO 40 IX=1,3
       POLD=FFB(3*(IA-1)+IX,iqb)*PHASE
C      Zero out at the start of the summation
       PNEW=0.0
       DO 41 II=1,3
41     PNEW=PNEW+FFS(3*(IAS-1)+II,iq)*A(IX,II)
       an=an+(PNEW-POLD)**2
       ans=ans+PNEW**2
       anb=anb+POLD**2
40     IF(LWR)WRITE(3,3006)IA,IX,iq,iqb,POLD,PNEW
3006   FORMAT(4I4,3F11.6)
45     CONTINUE
      ENDIF
C
20    CONTINUE
C
      RETURN
C
      END

c
      SUBROUTINE READRX(N0,N,R,ND,IBONDS,NB)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(3,N0),ND(N0),IBONDS(4,N0),NB(N0)
      READ(2,*)
      read(2,*)N
      IF(N.GT.N0)CALL REPORT(' N > MX !')
      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 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=10)
      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=10)
      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

      SUBROUTINE INPUTIN(NS,IB,IS,IND,MXT,MXB,INDBIG,IOX,NO,NAT,ND,NDS,
     1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LSTRICT,LINV,LSELECT,LAPT,LALF,
     2LRDO)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION IB(MXT,MXB),IS(MXT,MXB),IND(MXT),INDBIG(IOX,MXT),
     1ND(MXT),NDS(MXT),NDD(20)
      CHARACTER*7 STR
      LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT,LAPT,
     1LALF,LRDO
      LRDO=.FALSE.
      LALF=.FALSE.
      LWR=.FALSE.
      LAPT=.FALSE.
      LABS=.FALSE.
      LVCD=.FALSE.
      LOFF=.TRUE.
      LDIA=.TRUE.
      LROA=.FALSE.
      LSTRICT=.TRUE.
      LINV=.FALSE.
      LSELECT=.FALSE.
      IWG=1
9     READ(2,2000)STR
2000  FORMAT(A7)
      WRITE(3,2000)STR
      IF(STR(1:3).EQ.'IWG')THEN
       READ(2,*)IWG
       GOTO 9
      ENDIF
      IF(STR.EQ.'LSELECT')THEN
       READ(2,*)LSELECT
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LINV')THEN
       READ(2,*)LINV
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LAPT')THEN
       READ(2,*)LAPT
       GOTO 9
      ENDIF
      IF(STR.EQ.'LSTRICT')THEN
       READ(2,*)LSTRICT
       GOTO 9
      ENDIF
      IF(STR(1:3).EQ.'LWR')THEN
       READ(2,*)LWR
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LALF')THEN
       READ(2,*)LALF
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LRDO')THEN
       READ(2,*)LRDO
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LROA')THEN
       READ(2,*)LROA
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LOFF')THEN
       READ(2,*)LOFF
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LDIA')THEN
       READ(2,*)LDIA
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LABS')THEN
       READ(2,*)LABS
       GOTO 9
      ENDIF
      IF(STR(1:4).EQ.'LVCD')THEN
       READ(2,*)LVCD
       GOTO 9
      ENDIF
      NO=0
      IF(STR.NE.'POLYMER')THEN
      REWIND 2
      READ(2,*)NS
      DO 1 I=1,NS
      READ(2,*)IB(I,1),IS(I,1),IND(I)
      DO 1 J=1,ABS(IND(I))
1     READ(2,*)IB(I,J+1),IS(I,J+1)
C
C     NS atoms are assigned
C     IB(I) atom of the big molecule goes with the IS(I) atom of the smaller.
C     Each of them has IND(I) satelite atoms to define the reference  frame.
C     Reference atoms are given for big as IB and for small as IS strings.
C
      RETURN 
      ENDIF
      READ(2,*)NO
      WRITE(3,*)NO,' overlaps found, atoms assigned automatically ...'
      IF(NO.GT.IOX)CALL REPORT(' NO > IOX ! Change the dimensions.')
      DO 2 I=1,NO
      READ(2,*)IC
      IF(IC.NE.I)WRITE(3,*)IC,I
      IF(IC.NE.I)CALL REPORT('Error in reading ')
      READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT)
      DO 5 IAB=1,NAT
5     INDBIG(I,IAB)=0
      READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT)
      DO 3 IAB=1,NAT
      IF(INDBIG(I,IAB).EQ.0)GOTO 3
      IF(ND(IAB).NE.NDS(ABS(INDBIG(I,IAB))))THEN
       WRITE(3,3000)I,IAB,ND(IAB),NDS(ABS(INDBIG(I,IAB)))
3000  FORMAT('Overlap ',I2,', atom ',I3,': ',2I4)
       CALL REPORT(' Different atomic types !')
      ENDIF
3     CONTINUE
2     CONTINUE
2001  FORMAT(20I3)
      DO 4 IO=1,NO
      WRITE(3,2002)IO
2002  FORMAT(' Overlap number ',I5)
      IA=0
7     IA=IA+1
      IEND=IA+19
      IF(IEND.GT.NAT)IEND=NAT
      IF(IA.GT.NAT)GOTO 6
      DO 8 II=1,(IEND-IA+1)
8     NDD(II)=II+IA-1
      WRITE(3,2004)(NDD(II),II=1,IEND-IA+1)
      IA=IEND
      GOTO 7
6     WRITE(3,2004)(INDBIG(IO,IAB),IAB=1,NAT)
2004  FORMAT(20I3)
4     CONTINUE
      RETURN
      END
C
      SUBROUTINE REPORT(RS)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*(*) RS
      WRITE(3,3000)
      WRITE(6,3000)
3000  FORMAT(/,80(1H*))
      WRITE(3,*)RS
      WRITE(6,*)RS
      WRITE(3,3001)
      WRITE(6,3001)
3001  FORMAT(80(1H*),/,/,'PROGRAM STOPPED')
      CLOSE(3)
      CLOSE(2)
      STOP
      END
C
C                                                                      
c
      SUBROUTINE READS(N3,S,nat,x,z,nint,e)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION S(N3,N3),x(3,N3/3),z(N3/3),e(N3)
C     
      read(2,*)nint,ndumm,nat
      DO 1 I=1,NAT
1     read(2,*)z(I),(x(i1,I),i1=1,3)
      read(2,*)
      DO 2 I=1,NAT
      DO 2 J=1,NINT
2     read(2,*)kdumm,kdumm,s(3*(i-1)+1,j),
     1s(3*(i-1)+2,j),s(3*(i-1)+3,j)
      read(2,*)nint
      READ(2,4000)(E(I),I=1,nint)
4000  FORMAT(6F11.6)
      RETURN
      end