PROGRAM CCTS
      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),INDBIGT(IOX,MX),
     3NBO(MX),NBOB(MX),PS(MX,3,3),AS(MX,3,3),PB(MX,3,3),AB(MX,3,3),
     3PS2(MX,3,3),AS2(MX,3,3),PT(MX,3,3),AT(MX,3,3),NBO2(MX),
     4ALPHA(3*MX,3,3),ALPHAS(3*MX,3,3),AA(3*MX,3,3,3),
     5AAS(3*MX,3,3,3),G(3*MX,3,3),GS(3*MX,3,3),FFS2(3*MX,3*MX),
     5GS2(3*MX,3,3),AAS2(3*MX,3,3,3),ALPHAS2(3*MX,3,3),
     5GT(3*MX,3,3),AAT(3*MX,3,3,3),ALPHAT(3*MX,3,3),FFT(3*MX,3*MX),
     6RS2(3,MX),IBONDS2(4,MX),NDS2(MX)
      LOGICAL OPT,LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT,
     1LAPT,LALF,LRDO
      logical lxfile
c
      OPEN(3,FILE='CCTS.OUT')
      WRITE(3,3000)
      WRITE(6,3000)
3000  FORMAT(/,
     1       ' Perturbational Scaling of Molecular Tensors ',/,
     8       '           By Petr Bour, Prague 1999',/,/,
     9       ' TARGET MOLECULE: BIG.X    coordinates',/,
     1       '                  BIG. FC  force field',/,
     1       '                  BIG.TEN  APTs and AATs',/,
     3       '                  BIG.TTT  ROA tensors',/,/,
     3       ' TWO MOLECULES DEFINING THE PERTURBATION:',/,
     3       ' WORSE MOLECULE: ',/,
     3       ' SMALL1.X, SMALL1.FC, SMALL1.TEN, SMALL1.TTT',/, 
     3       ' BETTER MOLECULE: ',/,
     3       ' SMALL2.X, SMALL2.FC, SMALL2.TEN, SMALL2.TTT',/,/, 
     4       '                CCTS.INP  atom assignment',/,/,
     5       ' OUTPUT     :',/,
     7       ' FILE.X, FILE. FC, FILE.TEN, FILE.TTT',/,
     8       '                CCTS.OUT  control output file',/)
c
      inquire(FILE='BIG.X',EXIST=lxfile)
      if(lxfile)then
       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'
      else
       call report('BIG.X not found')
      endif
      inquire(FILE='SMALL1.X',EXIST=lxfile)
      if(lxfile)then
       OPEN(2,FILE='SMALL1.X',STATUS='OLD')
       CALL READRX(MX,NATS,RS,NDS,IBONDS,NBO)
       CLOSE(2)
       WRITE(3,*)NATS, ' atoms in the small molecule'
      else 
       call report('SMALL1.X not found')
      endif
      inquire(FILE='SMALL2.X',EXIST=lxfile)
      if(lxfile)then
       OPEN(2,FILE='SMALL2.X',STATUS='OLD')
       CALL READRX(MX,NATS2,RS2,NDS2,IBONDS2,NBO2)
       CLOSE(2)
       WRITE(3,*)NATS2, ' atoms in the second small molecule'
       if(NATS2.NE.NATS)call report('small mols must have same # of at')
      else 
       call report('SMALL2.X not found')
      endif
c
      OPEN(2,FILE='CCTS.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
c     read force fields:
      if(LDIA.or.LOFF)then
       INQUIRE(FILE='BIG.FC',EXIST=OPT)
       IF(.NOT.OPT)call report('BIG.FC not found')
       OPEN(2,FILE='BIG.FC',STATUS='OLD')
       CALL READFF(MX,NAT,FF)
       CLOSE(2)
       WRITE(3,3001)
       WRITE(6,3001)
3001   FORMAT(' Force field of the big molecule found')
       INQUIRE(FILE='SMALL1.FC',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('SMALL1.FC not found')
       OPEN(2,FILE='SMALL1.FC',STATUS='OLD')
       CALL READFF(MX,NATS,FFS)
       CLOSE(2)
       INQUIRE(FILE='SMALL2.FC',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('SMALL2.FC not found')
       OPEN(2,FILE='SMALL2.FC',STATUS='OLD')
       CALL READFF(MX,NATS,FFS2)
       CLOSE(2)
       WRITE(3,3002)
       WRITE(6,3002)
3002   FORMAT(' Force fields of the small molecules found')
       do 1 i=1,3*NATS
       do 1 j=1,3*NATS
1      FFT(i,j)=FFS2(i,j)
      endif
      IF(LABS.OR.LVCD)THEN
       INQUIRE(FILE='SMALL1.TEN',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('SMALL1.TEN not found')
       INQUIRE(FILE='SMALL2.TEN',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('SMALL2.TEN not found')
       INQUIRE(FILE='BIG.TEN',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('BIG.TEN not found')
       OPEN(15,FILE='SMALL1.TEN',STATUS='OLD')
       CALL READTEN(MX,PS,AS,LAPT,RS,NATS)
       CLOSE(15)
       OPEN(15,FILE='SMALL2.TEN',STATUS='OLD')
       CALL READTEN(MX,PS2,AS2,LAPT,RS2,NATS)
       CLOSE(15)
       OPEN(15,FILE='BIG.TEN')
       CALL READTEN(MX,PB,AB,LAPT,R,NAT)
       CLOSE(15)
       WRITE(3,3005)
       WRITE(6,3005)
3005   FORMAT(' Tensor files read in')
       do 2 i=1,NATS
       do 2 j=1,3
       do 2 k=1,3
       AT(i,j,k)=AS2(i,j,k)
2      PT(i,j,k)=PS2(i,j,k)
      ENDIF
C
      IF(LROA)THEN
       INQUIRE(FILE='BIG.TTT',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('BIG.TTT not found')
       INQUIRE(FILE='SMALL2.TTT',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('SMALL2.TTT not found')
       INQUIRE(FILE='SMALL1.TTT',EXIST=OPT)
       IF(.NOT.OPT)CALL REPORT('SMALL1.TTT not found')
       OPEN(15,FILE='SMALL1.TTT',STATUS='OLD')
       CALL READTTT(MX,ALPHAS,AAS,GS,NATS)
       CLOSE(15)
       ICO=1 
       CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,ICO,RS)
       OPEN(15,FILE='SMALL2.TTT',STATUS='OLD')
       CALL READTTT(MX,ALPHAS2,AAS2,GS2,NATS)
       CLOSE(15)
       CALL TTTT(3*MX,ALPHAS2,AAS2,GS2,NATS,ICO,RS2)
       OPEN(15,FILE='BIG.TTT')
       CALL READTTT(MX,ALPHA,AA,G,NAT)
       CLOSE(15)
       CALL TTTT(3*MX,ALPHA,AA,G,NAT,ICO,R)
       do 3 i=1,NATS*3
       do 3 j=1,3
       do 3 k=1,3
       ALPHAT(i,j,k)=ALPHAS2(i,j,k)
       GT(i,j,k)=GS2(i,j,k)
       do 3 l=1,3
3      AAT(i,j,k,l)=AAS2(i,j,k,l)
      ENDIF
C
      WRITE(6,*)' Transferring the tensors SMALL1->SMALL2 ...'
      NOT=1
      do 4 i=1,NATS
4     INDBIGT(1,i)=i
      CALL IMPROVE(FFT,FFS,RS2,RS,NATS,NATS,NS,IND,IB,IS,
     1AM,INDBIGT,NOT,IBONDS2,NBO2,LWR,PT,AT,LABS,LVCD,AS,PS,LOFF,LDIA,
     2IWG,LROA,ALPHAT,AAT,GT,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT)
c     CALL IMPROVE(FF,FFS,R,RS,NAT,NATS,NS,IND,IB,IS,
c    1AM,INDBIG,NO,IBONDSB,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA,
c    2IWG,LROA,ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT)
      if(LDIA.or.LOFF)then
       OPEN(2,FILE='TEM.FC')
       CALL WRITEFF(MX,NATS,FFT)
       CLOSE(2)
      endif
      IF(LABS.OR.LVCD)THEN
       OPEN(15,FILE='TEM.TEN')
       CALL WRITETEN(MX,PT,AT,RS2,NATS)
       CLOSE(15)
      ENDIF
      IF(LROA)THEN
       ICO=2
       CALL TTTT(3*MX,ALPHAT,AAT,GT,NATS,ICO,RS2)
       OPEN(2,FILE='TEM.TTT')
       CALL WRITETTT(MX,NATS,ALPHAT,AAT,GT)
       CLOSE(2)
      ENDIF
C
      WRITE(6,*)' Correcting the BIG ..'
      CALL IMPROVEC(INDBIG,NO,IWG,
     1LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,
     3FF  ,PB ,AB ,AA  ,G  ,ALPHA  ,R  ,NAT  ,IBONDSB,NBOB,
     3FFS2,PS2,AS2,AAS2,GS2,ALPHAS2,RS2,NATS2,
     3FFT ,PT ,AT ,AAT ,GT ,ALPHAT )
c
      OPEN(2,FILE='FILE.X')
      CALL WRITERX(MX,NAT,R,ND,IBONDSB)
      WRITE(3,*)' File FILE.X written '
      WRITE(6,*)' File FILE.X written '
      CLOSE(2)
      if(LDIA.or.LOFF)then
       OPEN(2,FILE='FILE.FC')
       CALL WRITEFF(MX,NAT,FF)
       CLOSE(2)
       WRITE(3,*)' File FILE.FC written '
       WRITE(6,*)' File FILE.FC written '
      endif
      IF(LABS.OR.LVCD)THEN
       OPEN(15,FILE='FILE.TEN')
       CALL WRITETEN(MX,PB,AB,R,NAT)
       CLOSE(15)
       WRITE(3,*)' File FILE.TEN written '
       WRITE(6,*)' File FILE.TEN written '
      ENDIF
      IF(LROA)THEN
       ICO=2
       CALL TTTT(3*MX,ALPHA,AA,G,NAT,ICO,R)
       OPEN(2,FILE='FILE.TTT')
       CALL WRITETTT(MX,NAT,ALPHA,AA,G)
       CLOSE(2)
       WRITE(3,*)' File FILE.TTT written '
       WRITE(6,*)' File FILE.TTT written '
      ENDIF
      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,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA,IWG,LROA,
     2ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT)
      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),PB(MX,3,3),AB(MX,3,3),PS(MX,3,3),
     5AS(MX,3,3),IOL(IOX),WF(IOX),ALPHA(3*MX,3,3),ALPHAS(3*MX,3,3),
     6AA(3*MX,3,3,3),AAS(3*MX,3,3,3),G(3*MX,3,3),GS(3*MX,3,3)
      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)GOTO 19
      IF(LDIA)WRITE(3,*)' Diagonal force field will be changed'
      IF(LOFF)WRITE(3,*)' Off-diagonals will be changed'
      IF(LABS)WRITE(3,*)' Eletric derivatives will be transformed'
      IF(LVCD)WRITE(3,*)' Magnetic derivatives will be transformed'
      IF(LROA)WRITE(3,*)' ROA parameters will be transformed'
      IF(LSTRICT)WRITE(3,*)' Only transformed atoms in the beds'
      IF(IWG.EQ.1)WRITE(3,*)' Multiple overlaps will be weighted'
      IF(LINV)WRITE(3,*)' Smaller fragment as optical antipode' 
      IF(LSELECT)WRITE(3,*)' FC for non-transformed atoms set to 0 !'
C
C     Look at all atom pairs i-j:
      DO 20 IA= 1,NAT
      DO 20 JA=IA,NAT
      IF(.NOT.LOFF.AND.IA.NE.JA)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))
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 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
       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
C      WRITE(3,3002)IA,JA
3002  FORMAT(' Pair ',2I4,': no overlap found')
       GOTO 20
      ELSE
       WRITE(3,3003)NF,IA,JA,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.000000001)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      Find the atomic groups to be placed on together for IA-JA:
C
       IPASS=1
C      Try two passes, firstly only connected atoms taken:
666    ITU=0
       IF(IPASS.EQ.2) WRITE(3,*)' Second pass ... '
C
       DO 27 II=1,NAT
       IF(INDBIG(IPO,II).EQ.0)GOTO 27
C
C      Take IA and JA:
       IF(II.EQ.IA.OR.II.EQ.JA)GOTO 277
       DO 29 IBA=1,NBOB(II)
C
C      Take atoms connected to IA or JA:
       IF(IBONDS(IBA,II).EQ.IA.OR.IBONDS(IBA,II).EQ.JA)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.OR.IBONDS(JB,IIA).EQ.JA)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)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 !')
       WRITE(3,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU)
       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      Transform the force field:
       IF((LDIA.AND.IA.EQ.JA).OR.(LOFF.AND.IA.NE.JA))THEN
        DO 36 IX=1,3
        INDIB=3*IA-3+IX
        JSTART=1
        IF(IA.EQ.JA)JSTART=IX
        DO 36 JX=JSTART,3
        INDJB=3*JA-3+JX
        FOLD=FFB(INDIB,INDJB)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)FFB(INDIB,INDJB)=0.0d0
        FNEW=0.0d0
        DO 37 IIX=1,3
        INDIS=3*INDBIG(IPO,IA)-3+IIX
        DO 37 JJX=1,3
        INDJS=3*INDBIG(IPO,JA)-3+JJX
37      FNEW=FNEW+A(IX,IIX)*FFS(INDIS,INDJS)*A(JX,JJX)
        IF(LWR)WRITE(3,3001)IA,JA, IX,JX,FOLD,FNEW,IIT,TOL,WF(IP)
        FFB(INDIB,INDJB)=FNEW*WF(IP)+FFB(INDIB,INDJB)
        FFB(INDJB,INDIB)=FFB(INDIB,INDJB)
36      CONTINUE
       ENDIF
C
C      Transform the polarization tensors
       IF(IA.EQ.JA.AND.(LABS.OR.LVCD))THEN
        IAS=INDBIG(IPO,IA)
        DO 40 IX=1,3
        DO 40 JX=1,3
        POLD=PB(IA,IX,JX)
        AOLD=AB(IA,IX,JX)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)THEN
         PB(IA,IX,JX)=0.0d0
         AB(IA,IX,JX)=0.0d0
        ENDIF
        PNEW=0.0d0
        ANEW=0.0d0
        DO 41 II=1,3
        AD1=A(IX,II)
        DO 41 JJ=1,3
        AD2=AD1*A(JX,JJ)
        PNEW=PNEW+PS(IAS,II,JJ)*AD2
41      ANEW=ANEW+AS(IAS,II,JJ)*AD2
        PB(IA,IX,JX)=PNEW*WF(IP)+PB(IA,IX,JX)
        AB(IA,IX,JX)=ANEW*WF(IP)+AB(IA,IX,JX)
        IF (.NOT.LVCD)AB(IA,IX,JX)=AOLD        
        IF (.NOT.LVCD)ANEW=AOLD        
        IF (.NOT.LABS)PB(IA,IX,JX)=POLD        
        IF (.NOT.LABS)PNEW=POLD        
40      IF(LWR)WRITE(3,3006)IA,IX,JX,POLD,PNEW,AOLD,ANEW,WF(IP)
3006    FORMAT(3I3,F11.6,'/',F11.6,3X,F11.6,'/',F11.6,' (',F7.4,')')
       ENDIF
C      Transform the ROA tensors
       IF(IA.EQ.JA.AND.LROA)THEN
        IAS=INDBIG(IPO,IA)
        DO 46 IX=1,3
        IIND=3*(IA-1)+IX
        DO 46 I=1,3
        DO 46 J=1,3
        ALPHAOLD=ALPHA(IIND,I,J)
        GOLD=G(IIND,I,J)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)THEN
         ALPHA(IIND,I,J)=0.0d0
         G(IIND,I,J)=0.0d0
        ENDIF
        ALPHANEW=0.0d0
        GNEW=0.0d0
        DO 47 IXX=1,3
        AD0=A(IX,IXX)
        INDS=3*(IAS-1)+IXX
        DO 47 II=1,3
        AD1=AD0*A(I,II)
        DO 47 JJ=1,3
        AD2=AD1*A(J,JJ)
        ALPHANEW=ALPHANEW+ALPHAS(INDS,II,JJ)*AD2
47      GNEW=GNEW+GS(INDS,II,JJ)*AD2
        ALPHA(IIND,I,J)=ALPHANEW*WF(IP)+ALPHA(IIND,I,J)
        G(IIND,I,J)=GNEW*WF(IP)+G(IIND,I,J)
46      IF(LWR)WRITE(3,3006)IIND,I,J,ALPHAOLD,ALPHANEW,
     1   GOLD,GNEW,WF(IP)
        DO 48 IX=1,3
        IIND=3*(IA-1)+IX
        DO 48 I=1,3
        DO 48 J=1,3
        DO 48 K=1,3
        AAOLD=AA(IIND,I,J,K)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)AA(IIND,I,J,K)=0.0d0
        AANEW=0.0d0
        DO 49 IXX=1,3
        AD0=A(IX,IXX)
        INDS=3*(IAS-1)+IXX
        DO 49 II=1,3
        AD1=AD0*A(I,II)
        DO 49 JJ=1,3
        AD2=AD1*A(J,JJ)
        DO 49 KK=1,3
49      AANEW=AANEW+AAS(INDS,II,JJ,KK)*AD2*A(K,KK)
        AA(IIND,I,J,K)=AANEW*WF(IP)+AA(IIND,I,J,K)
48      IF(LWR)WRITE(3,3007)IIND,I,J,K,AAOLD,AANEW
3007    FORMAT('A',4I3,F12.6,'/',F12.6)
       ENDIF
45     CONTINUE
      ENDIF
C
20    CONTINUE
C
C     If only interaction among transfered atoms selcted, the other 
C     force constants will be set to zero in final force field:
      IF(LSELECT)THEN
       DO 51 IA=1,NAT
       DO 52 IO=1,NO
52     IF(INDBIG(IO,IA).GT.0)GOTO 51
       WRITE(3,*)' Atom ',IA,' taken out of the FF'
       DO 53 IX=1,3
       INDIB=3*IA-3+IX
       DO 53 JA=1,NAT
       DO 53 JX=1,3
       INDJB=3*JA-3+JX
       FFB(INDIB,INDJB)=0.0d0
53     FFB(INDJB,INDIB)=0.0d0
51     CONTINUE
      ENDIF
      RETURN
C
19    WRITE(3,3000)
3000  FORMAT(1X,'          Old matrix element   New element',/,
     1       1X,'   atoms    xi   xj      F    [a.u.]  F',5H' #it,/,
     2       1X,'-------------------------------------------------')
      DO 1 IA=1,NS
      DO 3 I=1,3
      TB(I) =0.0d0
3     TS(I)=0.0d0
      DO 4 ISAT=1,ABS(IND(IA))+1
      DO 4 J=1,3
      XS(ISAT,J)=RS(J,IS(IA,ISAT))
      XB(ISAT,J)=RB(J,IB(IA,ISAT))
      TB(J)=TB(J)+XB(ISAT,J)
4     TS(J)=TS(J)+XS(ISAT,J)
      DO 5 ISAT=1,ABS(IND(IA))+1
      DO 5 J=1,3
      XS(ISAT,J)=XS(ISAT,J)-TS(J)/REAL(ABS(IND(IA))+1.0d0)
5     XB(ISAT,J)=XB(ISAT,J)-TB(J)/REAL(ABS(IND(IA))+1.0d0)
      IF(IND(IA).LT.0)THEN
       WRITE(3,*)' Space inversion applied for vicinity of atom ',IA
       DO 9 ISAT=1,ABS(IND(IA))+1
       DO 9 J=1,3
9      XS(ISAT,J)=-XS(ISAT,J)
      ENDIF
      ITU=ABS(IND(IA))+1
      CALL DOMATRIX(IERR)
      IF(IERR.NE.0)CALL REPORT('Erroneous stop')
      DO 10 I=1,3
      DO 10 J=1,3
10    AM(IA,I,J)=A(I,J)
      DO 1 JA=1,IA
      DO 11 I=1,3
      DO 11 J=1,3
11    B(I,J)=AM(JA,I,J)
      DO 6 I=1,3
      INDI=3*IB(IA,1)-3+I
      JSTART=1
      IF(IA.EQ.JA)JSTART=I
      DO 6 J=JSTART,3
      INDJ=3*IB(JA,1)-3+J
      FOLD=FFB(INDI,INDJ)
      FFB(INDI,INDJ)=0.0d0
      DO 7 II=1,3
      IL=3*IS(IA,1)-3+II
      DO 7 JJ=1,3
      JS=3*IS(IA,1)-3+JJ
7     FFB(INDI,INDJ)=FFB(INDI,INDJ)+A(I,II)*FFS(IL,JS)*B(J,JJ)
      IF(I.EQ.1.AND.J.EQ.I.AND.JA.EQ.IA)THEN
       WRITE(3,3001)IA,JA,I,J,FOLD,FFB(INDI,INDJ),IIT
      ELSE
       WRITE(3,3001)IA,JA,I,J,FOLD,FFB(INDI,INDJ)
      ENDIF
      FFB(INDJ,INDI)=FFB(INDI,INDJ)
6     CONTINUE
3001  FORMAT(4I5,2F11.5,I5,F10.7,F7.4)
1     CONTINUE
      RETURN
      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

c
      SUBROUTINE WRITERX(N0,N,R,ND,IB)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION R(3,N0),ND(N0),IB(4,N0)
      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)
      RETURN
      END

      SUBROUTINE READFF(N0,NA,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(3*N0,3*N0)
      NAT3=3*NA
      N=NAT3
C     Read in the lower triangle of FF,
C     written in parts as      n1,n1
C                               .      .
C                              ln,n1   .     ln,ln
C                               .      .      .
C                               .      .      .        n3,n3
C                               .      .      .         .
C                               n,n1   .      .         n,n3
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(2,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
C
      DO 31 I=1,N
      DO 31 J=I,N
   31 FCAR(I,J)=FCAR(J,I)
      RETURN
      END

      SUBROUTINE WRITEFF(MX,NA,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(3*MX,3*MX)
      N=NA*3
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   WRITE(2,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 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.00000000001)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
      IF(STR(1:4).EQ.'POLY')THEN
       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
      ENDIF
      write(6,2000)STR
      write(3,2000)STR
      call report('unknown option')
      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
      SUBROUTINE READTEN(N0,P,A,LAPT,X,NAT)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MX=600)
      DIMENSION P(N0,3,3),A(N0,3,3),AI(MX,3,3),AJ(MX,3,3),PV(MX,3,3),
     1R(3,MX),X(3,NAT)
      LOGICAL LAPT
      BOHR=0.52917705993d0
      DO 3 L=1,NAT
      DO 3 I=1,3
3     R(I,L)=X(I,L)/BOHR
      READ (15,*) NOAT,NDUM,NDUM
      NAT=NOAT
      WRITE(3,*)NAT,' atoms'
      WRITE(*,*)NAT,' atoms'
C
C     Read dipole derivatives:
      DO 10 L=1,NAT
      DO 10 J=1,3                                                     
10    READ (15,*) (P(L,J,I),I=1,3)
      WRITE(6,*)' Dipole derivatives read - in '
      WRITE(3,*)' Dipole derivatives read - in '
C
C     If no axial tensor, then construct its polar part only:
      IF(LAPT)THEN
       DO 2203 L=1,NAT
       DO 2203 J=1,3
       DO 2203 I=1,3
       SUM=0.0d0
       DO 2204 ID=1,3
       DO 2204 IG=1,3
2204   SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID)
       AJ(L,J,I)=0.0d0
       PV(L,J,I)=P(L,J,I)
2203   AI(L,J,I)=SUM
       WRITE(6,*)' APT part of AAT constructed'
       WRITE(3,*)' APT part of AAT constructed'
       GOTO 1
      ENDIF
      DO 220 L=1,NOAT                                                 
      DO 220 J=1,3                                                    
220   READ (15,*) (AI(L,J,I),I=1,3)
      DO 230 L=1,NOAT                                                 
      DO 230 J=1,3                                                    
230   READ (15,*) (AJ(L,J,I),I=1,3)                                   
      DO 100 L=1,NOAT                                                 
      DO 100 J=1,3                                                    
100   READ (15,*) (PV(L,J,I),I=1,3)                                   
      WRITE(6,*)' VCD parameters read in ...'
      WRITE(3,*)' VCD parameters read in ...'
C     Distributed-origin gauge:
1     CALL VCDD0(AI,PV,P,R,NAT)   
      WRITE(6,*)'DOG used'
      WRITE(3,*)'DOG used'
C     
C     Total axial tensor = electronic + nuclear:
      DO 4 L=1,NAT
      DO 4 I=1,3
      DO 4 J=1,3
4     A(L,I,J)=AI(L,I,J)+AJ(L,I,J)
C
C     Put axial tensor in the origin on the atom L:
      DO 2201 L=1,NAT
      DO 2201 J=1,3
      DO 2201 I=1,3
      SUM=0.0d0
      DO 2202 ID=1,3
      DO 2202 IG=1,3
2202  SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID)
2201  A(L,J,I)=A(L,J,I)+SUM
      WRITE(3,*)' AAT shifted to atomic origin'
      WRITE(6,*)' AAT shifted to atomic origin'
      RETURN
      END
C
      SUBROUTINE WRITETEN(N0,P,A,X,NAT)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MX=600)
      DIMENSION P(N0,3,3),A(N0,3,3),R(3,MX),X(3,NAT),AT(N0,3,3)
      BOHR=0.52917705993d0
      DO 3 L=1,NAT
      DO 3 I=1,3
3     R(I,L)=X(I,L)/BOHR
      DO 2201 L=1,NAT
      DO 2201 J=1,3
      DO 2201 I=1,3
      SUM=0.0d0
      DO 2202 ID=1,3
      DO 2202 IG=1,3
2202  SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID)
2201  AT(L,J,I)=A(L,J,I)+SUM
      WRITE(3,*)' AAT shifted back to laboratory system'
      WRITE(6,*)' AAT shifted back to laboratory system'
      WRITE(15,1500) NAT,NAT-6,0
1500  FORMAT(3I5)
      DO 10 L=1,NAT
      DO 10 J=1,3
10    WRITE(15,1501) (P(L,J,I),I=1,3),L
1501  FORMAT(3F14.8,I5)
      DO 220 L=1,NAT                                                 
      DO 220 J=1,3                                                    
220   WRITE(15,1501) (AT(L,J,I),I=1,3),L
      Z0=0.0d0
      DO 230 L=1,NAT                                                 
      DO 230 J=1,3                                                    
230   WRITE(15,1501) (Z0,I=1,3),L
      DO 100 L=1,NAT
      DO 100 J=1,3
100   WRITE(15,1501) (P(L,J,I),I=1,3),L
      RETURN
      END
C                                                                      
      FUNCTION EPS(I,J,K)
      REAL*8 EPS
      INTEGER*4 I,J,K
      EPS=0.0d0
      IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0
      IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0
      IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0
      IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0
      IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0
      IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0
      RETURN                                                          
      END                                                             
C 
      SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT)
C     STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA                           
      PARAMETER (MX=600)
      DIMENSION VCD(MX,3,3),VEL(MX,3,3),DIP(MX,3,3),                  
     1C(3,NAT)
      DO 90 BETA =1,3                                                 
      DO 90 GAMMA=1,3                                                 
      DO 90 DELTA=1,3                                                 
      SKEW=0.25D0*EPS(BETA,GAMMA,DELTA)                               
      IF (DABS(SKEW).GT.1.0D-10)THEN                                          
      DO 80 ALPHA=1,3                                                 
      DO 80 LAMBDA=1,NAT                                              
80    VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA)                   
     1 -SKEW*C(GAMMA,LAMBDA)                                          
     2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA))             
      ENDIF                                                           
90    CONTINUE
      RETURN
      END
C
      SUBROUTINE READTTT(MX,ALPHA,A,G,NAT)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3)
      READ(15,*)
      READ(15,*)NAT
      READ(15,*)
      READ(15,*)
      DO 1 I=1,3
      READ(15,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IIND=3*(L-1)+IX
1     READ(15,*)LDUM,IXDUM,(ALPHA(IIND,I,J),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 2 I=1,3
      READ(15,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IIND=3*(L-1)+IX
2     READ(15,*)LDUM,IXDUM,(G(IIND,I,J),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(15,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IIND=3*(L-1)+IX
3     READ(15,*)LDUM,IXDUM,(A(IIND,I,J,K),K=1,3)
      CLOSE(15)
      RETURN
      END
C
      SUBROUTINE WRITETTT(MX,NAT,ALPHA,A,G)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3)
      WRITE(2,2000)NAT
2000  FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/,
     1' The electric-dipolar electric-dipolar polarizability:',/,
     2' Atom/x    jx           jy           jz')
      DO 1 I=1,3
      WRITE(2,2002)I
2002  FORMAT(' Alpha(',I1,',J):')
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IIND=3*(L-1)+IX
1     WRITE(2,2001)L,IX,(ALPHA(IIND,I,J),J=1,3)
2001  FORMAT(I5,1H ,I1,3F13.7)
      WRITE(2,2004)
2004  FORMAT(' The electric dipole magnetic dipole polarizability:',/,
     1       ' Atom/x    jx           jy           jz')
      DO 2 I=1,3
      WRITE(2,2003)I
2003  FORMAT(' G(',I1,',J):')
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IIND=3*(L-1)+IX
2     WRITE(2,2001)L,IX,(G(IIND,I,J),J=1,3)
      WRITE(2,2005)
2005  FORMAT(' The electric dipole electric quadrupole polarizability:',
     2/,     ' Atom/x    kx           ky           kz')
      DO 3 I=1,3
      DO 3 J=1,3
      WRITE(2,2006)I,J
2006  FORMAT(' A(',I1,',',I1,',K):')
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IIND=3*(L-1)+IX
3     WRITE(2,2001)L,IX,(A(IIND,I,J,K),K=1,3)
      RETURN
      END

      SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R)
C     Transforms A and G to local origins (ICO=1) or common origin (ICO=2);
C     sets A and G to zero for ICO=3.
C     R supposed to be passed in in A
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      PARAMETER (MX=600)
      DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),X(3,MX)
C
C     transfer into atomic coordinates:
      DO 5 I=1,3
      DO 5 J=1,NAT
5     X(I,J)=R(I,J)/0.529177D0
C
      IF(ICO.EQ.1)SIGN=1.0d0
      IF(ICO.EQ.2)SIGN=-1.0d0
      IF(ICO.EQ.3)SIGN=0.0d0
      W=1.0d0
C     supposedly the static limit G/w is calculated
C
      DO 2 IA=1,3
      DO 2 IB=1,3
      DO 2 L=1,NAT
      DO 2 IE=1,3
      IIND=3*(L-1)+IE
      SUM=0.0d0
      DO 1 IG=1,3
      DO 1 ID=1,3
1     SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID)
      G(IIND,IA,IB)=G(IIND,IA,IB)+SUM
      IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0
2     CONTINUE
      IF(ICO.EQ.1)WRITE(3,*)' G transformed into atomic origins'
      IF(ICO.EQ.2)WRITE(3,*)' G transformed into common origin'
      IF(ICO.EQ.3)WRITE(3,*)' G set to zero'
C
      DO 3 IA=1,3
      DO 3 IB=1,3
      DO 3 IC=1,3
      DO 3 L=1,NAT
      DO 3 IE=1,3
      IIND=3*(L-1)+IE
      SUM=0.0d0
      IF(IB.EQ.IC)THEN
       DO 4 ID=1,3
4      SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN
      ENDIF
      A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM
     1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB))
      IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0
3     CONTINUE
      IF(ICO.EQ.1)WRITE(3,*)' A transformed into atomic origins'
      IF(ICO.EQ.2)WRITE(3,*)' A transformed into common origin'
      IF(ICO.EQ.3)WRITE(3,*)' A set to zero'
      RETURN
      END
c     ===========================================================
      SUBROUTINE IMPROVEC(INDBIG,NO,IWG,
     1LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,
     3FFB ,PB ,AB ,AA  ,G  ,ALPHA  ,RB ,NAT  ,IBONDS ,NBOB,
     3FFS ,PS ,AS ,AAS ,GS ,ALPHAS ,RS ,NATS ,
     3FFT ,PT ,AT ,AAT ,GT ,ALPHAT )
c
c     FFB - big force field; tensors PB,AB,AA,G,ALPHA;coord R
c     FFS - the good small FF
c     FFT - the worse small FF
c
      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),   FFT(3*MX,3*MX),
     1          RB(3,NAT),       RS(3,NATS),       
     1          IBONDS(4,MX),    
     1          NBOB(MX),        
     1          PB(MX,3,3),      PS(MX,3,3),       PT(MX,3,3),
     1          AB(MX,3,3),      AS(MX,3,3),       AT(MX,3,3),
     1          ALPHA(3*MX,3,3), ALPHAS(3*MX,3,3), ALPHAT(3*MX,3,3),
     1          AA(3*MX,3,3,3),  AAS(3*MX,3,3,3),  AAT(3*MX,3,3,3),
     1          G(3*MX,3,3),     GS(3*MX,3,3),     GT(3*MX,3,3),
     1INDBIG(IOX,MX)
      DIMENSION TS(3),TB(3),TIJ(3),OTIJ(3),IOL(IOX),WF(IOX),IND(IOX)
      LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT
c
C     Takes big force field FFB and implants corrections with respect to
c     FFS and FFT:      "FFB=FFB+(FFS-FFT)"
c
      IF(LDIA)WRITE(3,*)' Diagonal force field will be changed'
      IF(LOFF)WRITE(3,*)' Off-diagonals will be changed'
      IF(LABS)WRITE(3,*)' Eletric derivatives will be transformed'
      IF(LVCD)WRITE(3,*)' Magnetic derivatives will be transformed'
      IF(LROA)WRITE(3,*)' ROA parameters will be transformed'
      IF(LSTRICT)WRITE(3,*)' Only transformed atoms in the beds'
      IF(IWG.EQ.1)WRITE(3,*)' Multiple overlaps will be weighted'
      IF(LINV)WRITE(3,*)' Smaller fragment as optical antipode' 
C
C     Look at all atom pairs i-j:
      DO 20 IA= 1,NAT
      DO 20 JA=IA,NAT
      IF(.NOT.LOFF.AND.IA.NE.JA)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))
C
C     Look at all overlaps io, find that which fits best the i-j pair:
      DISTM=100000.0d0
      NF=0
c
c     Lopp over all overlaps:
      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
       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
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
C      WRITE(3,3002)IA,JA
3002   FORMAT(' Pair ',2I4,': no overlap found')
       GOTO 20
      ELSE
       WRITE(3,3303)NF,IA,JA,IOIJ
3303   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.0d0/WF(I)
       DO 42 I=1,NF
       IF(WF(I).LE.0.000000001d0)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      Find the atomic groups to be placed on together for IA-JA:
C
       IPASS=1
C      Try two passes, firstly only connected atoms taken:
666    ITU=0
       IF(IPASS.EQ.2) WRITE(3,*)' Second pass ... '
C
       DO 27 II=1,NAT
       IF(INDBIG(IPO,II).EQ.0)GOTO 27
C
C      Take IA and JA:
       IF(II.EQ.IA.OR.II.EQ.JA)GOTO 277
       DO 29 IBA=1,NBOB(II)
C
C      Take atoms connected to IA or JA:
       IF(IBONDS(IBA,II).EQ.IA.OR.IBONDS(IBA,II).EQ.JA)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.OR.IBONDS(JB,IIA).EQ.JA)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)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 !')
       WRITE(3,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU)
       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')
       if(LWR)write(3,*)IIT,' iterations'
C
       IF(LINV)THEN
        DO 50 I=1,3
        DO 50 J=1,3
50      A(I,J)=-A(I,J)
       ENDIF
C
C      Correct the force field:
       IF((LDIA.AND.IA.EQ.JA).OR.(LOFF.AND.IA.NE.JA))THEN
        if(LWR)write(3,3003)
3003    format('ia    jb    Fiajb_old    Fiajb_new')
        DO 36 IX=1,3
        INDIB=3*IA-3+IX
        JSTART=1
        IF(IA.EQ.JA)JSTART=IX
        DO 36 JX=JSTART,3
        INDJB=3*JA-3+JX
        FOLD=FFB(INDIB,INDJB)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)FFB(INDIB,INDJB)=0.0D0
        DEL=0.0d0
        DO 37 IIX=1,3
        INDIS=3*INDBIG(IPO,IA)-3+IIX
        DO 37 JJX=1,3
        INDJS=3*INDBIG(IPO,JA)-3+JJX
37      DEL=DEL+A(IX,IIX)*(FFS(INDIS,INDJS)-FFT(INDIS,INDJS))*A(JX,JJX)
        FNEW=FOLD+DEL
        IF(LWR)WRITE(3,3001)IA,IX,JA,JX,FOLD,FNEW
3001    format(2(i4,i2),2f10.5)
        FFB(INDIB,INDJB)=FNEW*WF(IP)+FFB(INDIB,INDJB)
        FFB(INDJB,INDIB)=FFB(INDIB,INDJB)
36      CONTINUE
       ENDIF
C
C      Correct the polarization tensors
       IF(IA.EQ.JA.AND.(LABS.OR.LVCD))THEN
        if(LWR)write(3,3307)
3307    format('iab    APTold   APTnew   AATold   AATnew')
        IAS=INDBIG(IPO,IA)
        DO 40 IX=1,3
        DO 40 JX=1,3
        POLD=PB(IA,IX,JX)
        AOLD=AB(IA,IX,JX)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)THEN
         PB(IA,IX,JX)=0.0d0
         AB(IA,IX,JX)=0.0d0
        ENDIF
        PDEL=0.0d0
        ADEL=0.0d0
        DO 41 II=1,3
        AD1=A(IX,II)
        DO 41 JJ=1,3
        AD2=AD1*A(JX,JJ)
        PDEL=PDEL+(PS(IAS,II,JJ)-PT(IAS,II,JJ))*AD2
41      ADEL=ADEL+(AS(IAS,II,JJ)-AT(IAS,II,JJ))*AD2
        PNEW=POLD+PDEL
        ANEW=AOLD+ADEL
        PB(IA,IX,JX)=PNEW*WF(IP)+PB(IA,IX,JX)
        AB(IA,IX,JX)=ANEW*WF(IP)+AB(IA,IX,JX)
        IF (.NOT.LVCD)AB(IA,IX,JX)=AOLD        
        IF (.NOT.LVCD)ANEW=AOLD        
        IF (.NOT.LABS)PB(IA,IX,JX)=POLD        
        IF (.NOT.LABS)PNEW=POLD        
40      IF(LWR)WRITE(3,3006)IA,IX,JX,POLD,PNEW,AOLD,ANEW
3006    FORMAT(i4,2I2,4F10.5)
       ENDIF
c
C      Correct the ROA tensors
       IF(IA.EQ.JA.AND.LROA)THEN
        if(LWR)write(3,3008)
3008    format('iab    ALPHAold   ALPHAnew   Gold   Gnew')
        IAS=INDBIG(IPO,IA)
        DO 46 IX=1,3
        IIND=3*(IA-1)+IX
        DO 46 I=1,3
        DO 46 J=1,3
        ALPHAOLD=ALPHA(IIND,I,J)
        GOLD=G(IIND,I,J)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)THEN
         ALPHA(IIND,I,J)=0.0d0
         G(IIND,I,J)=0.0d0
        ENDIF
        DG=0.0d0
        DA=0.0d0
        DO 47 IXX=1,3
        AD0=A(IX,IXX)
        INDS=3*(IAS-1)+IXX
        DO 47 II=1,3
        AD1=AD0*A(I,II)
        DO 47 JJ=1,3
        AD2=AD1*A(J,JJ)
        DA=DA+(ALPHAS(INDS,II,JJ)-ALPHAT(INDS,II,JJ))*AD2
47      DG=DG+(GS(INDS,II,JJ)    -GT(INDS,II,JJ)    )*AD2
        ALPHANEW=ALPHAOLD+DA
        GNEW=GOLD+DG
        ALPHA(IIND,I,J)=ALPHANEW*WF(IP)+ALPHA(IIND,I,J)
        G(IIND,I,J)=GNEW*WF(IP)+G(IIND,I,J)
46      IF(LWR)WRITE(3,3006)IIND,I,J,ALPHAOLD,ALPHANEW,
     1   GOLD,GNEW
        if(LWR)write(3,3009)
3009    format('iab    Aold   Anew')
        DO 48 IX=1,3
        IIND=3*(IA-1)+IX
        DO 48 I=1,3
        DO 48 J=1,3
        DO 48 K=1,3
        AAOLD=AA(IIND,I,J,K)
C       Zero out at the start of the summation
        IF(ISUM.EQ.1)AA(IIND,I,J,K)=0.0d0
        DA=0.0d0
        DO 49 IXX=1,3
        AD0=A(IX,IXX)
        INDS=3*(IAS-1)+IXX
        DO 49 II=1,3
        AD1=AD0*A(I,II)
        DO 49 JJ=1,3
        AD2=AD1*A(J,JJ)
        DO 49 KK=1,3
49      DA=DA+(AAS(INDS,II,JJ,KK)-AAT(INDS,II,JJ,KK))*AD2*A(K,KK)
        AANEW=AAOLD+DA
        AA(IIND,I,J,K)=AANEW*WF(IP)+AA(IIND,I,J,K)
48      IF(LWR)WRITE(3,3007)IA,IX,I,J,K,AAOLD,AANEW
3007    FORMAT(i4,4I2,2F10.5)
       ENDIF
45     CONTINUE
      ENDIF
C
20    CONTINUE
      RETURN
      END