C************************************************************************
C* Release September 1996                                               *
C*                                                                      *
C*                     Program CHANGDAT                                 *
C*                                                                      *
C*                      by M. Nardelli                                  *
C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica    *
C*     e Chimica Fisica della Universita' degli Studi di Parma,         *
C* Centro di Studio per la Strutturistica Diffrattometrica del CNR,     *
C*          Viale delle Scienze 78, I-43100 Parma, Italy                *
C*                                                                      *
C* This program modifies the lattice parameters, reflection indices,    *
C* atomic coordinates and anisotropic displacement parameters according *
C* to given rotation matrices and translation vector.                   *
C*    Three options are available:                                      *
C*     (1) the reflection indices in the F^2^ (or F) SHELX input file;  *
C* the cell constants, the atomic coordinates and the anisotropic       *
C* displacement parameters:                                             *
C*     (2) in the SHELXL-93 name.ins input file,                        *
C*     (3) in the PARST input file.                                     *
C* The matrices and vector are keyed in on request.                     *
C*                                                                      *
C* N.B. - If the lattice parameters are not to be changed, the hkl      *
C*        transformation matrix must be (1 0 0 /0 1 0/0 0 1)            *
C*      - If the transformation matrix of coordinates, RX, is unknown,  *
C*        but the transformation matrix of cell parameters, RH, has     *
C*        been given, then give for RX a matrix of zeros.               *
C************************************************************************
      PROGRAM CHANGDAT
      COMMON/IOU/IN,IO
     1      /GRP1/RH(3,3)
     2      /GRP2/RX(3,3),TX(3)
      DOUBLE PRECISION A(3,3),B(3,3)
      IN=15
      IO=16
      CALL INOU
      WRITE(*,100)
100   FORMAT(' According to your choice, key in:'/
     1 ' 0 if you want to change the intensity (or F) file'/
     2 ' 1 if you want to change the SHELXL-93 name.ins file'/
     3 ' 2 if you want to change the PARST input file'/' ?>')
      READ(*,*)KDC
      WRITE(*,101)
101   FORMAT(' Key in the hkl transformation matrix by rows ',
     1 '(FREE FORMAT)'/' ?>')
      READ(*,*)((RH(I,J),J=1,3),I=1,3)
      IF(KDC.EQ.0)THEN
      CALL NEWHKL
      ELSE
      WRITE(*,102)
102   FORMAT(' Key in the xyz transformation matrix by rows ',
     1 '(FREE FORMAT)'/' If you do not know this matrix, then give',
     2 ' a matrix of zeros'/' ?>')
      READ(*,*)((RX(I,J),J=1,3),I=1,3)
      DO 1 I=1,3
      DO 1 J=1,3
      IF(ABS(RX(I,J)).GT.0.01)GO TO 4
1     CONTINUE
      DO 2 I=1,3
      DO 2 J=1,3
2     A(I,J)=RH(I,J)
      CALL INVMAT3(A,B,IER)
      DO 3 I=1,3
      DO 3 J=1,3
3     RX(J,I)=B(I,J)
4     WRITE(*,103)
103   FORMAT(' Key in the Tx translation vector (FREE FORMAT)'/' ?>')
      READ(*,*)(TX(I),I=1,3)
      ENDIF
      IF(KDC.EQ.1)THEN
      CALL NEWSHELX
      WRITE(IO,104)
104   FORMAT(/' WARNING: If the number of molecules and the SYMM ',
     1 'and UNIT instructions'/' =======  have to be changed, these ',
     2 'changes must be done by editor.')
      ELSEIF(KDC.EQ.2)THEN
      CALL NEWPARST
      WRITE(IO,105)
105   FORMAT(/' WARNING: If the Bravais symbol, the number of ',
     1 'molecules, the centric (-1)'/' =======  or acentric (1) code',
     2 ' and equivalent positions have to be changed,'/ 10X,'these ',
     3 'changes must be done by editor.')
      ENDIF
      CLOSE (IN)
      CLOSE (IO)
      STOP
      END
 
      SUBROUTINE NEWHKL
C-----Trasforma gli indici del file di input nome.hkl per SHELXL-93
C-----secondo una matrice richiesta in interattivo
      DIMENSION AM(3,3),IND(3),AN(3),IN1(3)
      COMMON/IOU/IN,IO
     1      /GRP1/RH(3,3)
1     READ(IN,101,END=5)(IND(I),I=1,3),F2,SF2,IA
101   FORMAT(3I4,2F8.2,I4)
      IF(IND(1).EQ.0.AND.IND(2).EQ.0.AND.IND(3).EQ.0)GO TO 4
      DO 3 I=1,3
      AN(I)=0.0
      DO 2 J=1,3
2     AN(I)=AN(I)+RH(I,J)*IND(J)
3     IN1(I)=AN(I)
      WRITE(IO,101)(IN1(I),I=1,3),F2,SF2,IA
      GO TO 1
4     WRITE(IO,102)
102   FORMAT(3(3X,'0'),3(4X,'0.00'))
5     RETURN
      END
 
      SUBROUTINE ROTU(KKD,RI,U,U1)
C-----This routine calculates the U rotating matrix
      DOUBLE PRECISION A(3,3),B(3,3),U(3,3),UI(3,3),RI(3,3),
     1 RIT(3,3),RU(3,3),U1(3,3)
      COMMON/IOU/IN,IO
      CALL INVMAT3(U,UI,IER)
      IF(IER.EQ.1)RETURN
      DO 1 I=1,3
      DO 1 J=1,3
1     A(I,J)=0D0
      DO 2 I=1,3
      DO 2 J=1,3
      DO 2 K=1,3
2     A(I,J)=A(I,J)+RI(I,K)*UI(K,J)
      DO 3 I=1,3
      DO 3 J=1,3
3     B(I,J)=0D0
      DO 4 I=1,3
      DO 4 J=1,3
      DO 4 K=1,3
4     B(I,J)=B(I,J)+U(I,K)*A(K,J)
      DO 5 I=1,3
      DO 5 J=1,3
5     RIT(J,I)=RI(I,J)
      DO 6 I=1,3
      DO 6 J=1,3
6     RU(I,J)=0D0
      DO 7 I=1,3
      DO 7 J=1,3
      DO 7 K=1,3
7     RU(I,J)=RU(I,J)+RIT(I,K)*B(K,J)
      DO 8 I=1,3
      DO 8 J=1,3
8     U1(I,J)=0D0
      IF(KKD.EQ.0)THEN
      DO 9 I=1,3
      DO 9 J=1,3
      DO 9 K=1,3
9     U1(I,J)=U1(I,J)+RU(I,K)*U(K,J)
      ELSE
      DO 10 I=1,3
      DO 10 J=1,3
      DO 10 K=1,3
10    U1(I,J)=U1(I,J)+RU(I,K)*U(K,J)
      ENDIF
      RETURN
      END
 
      SUBROUTINE INVMAT3(A,B,IER)
C-----This routine inverts a 3x3 matrix
      DOUBLE PRECISION A(3,3),B(3,3),D
      COMMON/IOU/IN,IO
      IER=0
      D=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))-
     1  A(1,2)*(A(2,1)*A(3,3)-A(2,3)*A(3,1))+
     2  A(1,3)*(A(2,1)*A(3,2)-A(3,1)*A(2,2))
      IF(DABS(D).LT.1D-10) GO TO 1
      B(1,1)=(A(2,2)*A(3,3)-A(3,2)*A(2,3))/D
      B(2,1)=-(A(2,1)*A(3,3)-A(3,1)*A(2,3))/D
      B(3,1)=(A(2,1)*A(3,2)-A(3,1)*A(2,2))/D
      B(1,2)=-(A(1,2)*A(3,3)-A(3,2)*A(1,3))/D
      B(2,2)=(A(1,1)*A(3,3)-A(3,1)*A(1,3))/D
      B(3,2)=-(A(1,1)*A(3,2)-A(3,1)*A(1,2))/D
      B(1,3)=(A(1,2)*A(2,3)-A(2,2)*A(1,3))/D
      B(2,3)=-(A(1,1)*A(2,3)-A(2,1)*A(1,3))/D
      B(3,3)=(A(1,1)*A(2,2)-A(2,1)*A(1,2))/D
      GO TO 2
1     WRITE(IO,100)
100   FORMAT(' The matrix cannot be inverted')
      IER=1
2     RETURN
      END
 
      SUBROUTINE INOU
C-----Assign input and output files
      CHARACTER FILIN*14,STR1*23,STRI*2,STRO*3,STR2*8,STR3*2
      COMMON/IOU/IN,IO
      STR1='Key in the name of the '
      STRI='in'
      STRO='out'
      STR2='put file'
      STR3='?>'
      WRITE(*,100)STR1,STRI,STR2,STR3
100   FORMAT(1X,A22,1X,A2,A8/1X,A2)
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IN,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IN)
      WRITE(*,101)STR1,STRO,STR2,STR3
101   FORMAT(1X,A22,1X,A3,A8/1X,A2)
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IO,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IO)
      RETURN
      END
 
      SUBROUTINE NEWPARST
      CHARACTER*1 TIT(80),ATOM(6),UNIM(3,5),II(5)
      CHARACTER*80 FRT
      DIMENSION C(3),SC(3),A(3),SA(3),X(3),SX(3),B(6),SB(6),TAU(3),
     1 STAU(3),AN(3),SAN(3),X1(3),SX1(3)
      DOUBLE PRECISION U(3,3),SU(3,3),RHI(3,3),U1(3,3),SU1(3,3),
     1 RHD(3,3),RDXI(3,3)
      COMMON/IOU/IN,IO
     1      /GRP1/RH(3,3)
     2      /GRP2/RX(3,3),TX(3)
     3      /GRP3/C,SC,A,SA
      PG=3.14159265
      RD=57.295779513
C-----Read title, codes and parameters
      DO 1 I=1,2
      READ(IN,'(80A1)') TIT
1     WRITE(IO,'(80A1)') TIT
      READ(IN,*)N,D3,DM,NH1,NTN,LSP,LSL,NPR,NST,KL,KY,KC,KO,KB,KA,KT,
     1 KD,KE,KOR,NI
      BACKSPACE (IN)
      READ(IN,'(80A1)')TIT
      WRITE(IO,'(80A1)') TIT
      IF(KY.EQ.0) GO TO 6
      DO 2 I=1,2
      READ(IN,'(80A1)') TIT
2     WRITE(IO,'(80A1)') TIT
6     READ(IN,*)C,SC,A,SA
      DO 3 I=1,3
3     IF(SC(I).LE.1E-9)SC(I)=C(I)*1E-9
      DO 18 I=1,3
      DO 18 J=1,3
      IF(ABS(RH(I,I)-1.0).GT.0.01) GO TO 19
      IF(I.NE.J.AND.ABS(RH(I,J)).GT.0.01)GO TO 19
18    CONTINUE
      DO 21 I=1,3
      TAU(I)=C(I)
      STAU(I)=SC(I)
      AN(I)=A(I)
21    SAN(I)=SA(I)
      GO TO 20
19    CALL NEWCELL(TAU,STAU,AN,SAN)
20    WRITE(IO,'(3(F8.4,1X),3(F6.4,1X),3(F6.2,1X),3(F4.2,1X))')
     1 TAU,STAU,AN,SAN
      CALL CHANGRH(KKK)
C-----Read atom symbols,coordinates,thermal parameters & e.s.d.'s
      IF(KL.EQ.4) GO TO 7
      READ(IN,'(A80)')FRT
      WRITE(IO,'(A80)')FRT
7     DO 4 IK=1,N
      IF(KOR.EQ.0)THEN
      IF(KL.EQ.1.OR.KL.EQ.2)READ(IN,FRT)ATOM,X,SX
      IF(KL.EQ.3)READ(IN,FRT)ATOM,(X(J),SX(J),J=1,3)
      IF(KL.EQ.4)THEN
      READ(IN,100)ATOM,X,((UNIM(I,K),K=1,5),I=1,3)
100   FORMAT(7X,6A1,2X,3F8.5,12X,3(5A1,3X))
      DO 9 I=1,3
      DO 8 JJ=1,5
8     II(JJ)=UNIM(I,JJ)
      CALL CONVE(II,DH)
9     SX(I)=DH
      ENDIF
      ELSE
      IF(KL.EQ.1)READ(IN,FRT)ATOM,X,SX,B,SB
      IF(KL.EQ.2)READ(IN,FRT)ATOM,X,B,SX,SB
      IF(KL.EQ.3)READ(IN,FRT)ATOM,(X(J),SX(J),J=1,3),(B(K),
     1 SB(K),K=1,6)
      ENDIF
C-----Change the x,y,z coordinates
      DO 17 I=1,3
      X1(I)=0.0
      SX1(I)=0.0
      DO 10 J=1,3
      X1(I)=X1(I)+RX(I,J)*X(J)
10    SX1(I)=SX1(I)+RX(I,J)*SX(J)
      SX1(I)=ABS(SX1(I))
17    X1(I)=X1(I)+TX(I)
      IF(KOR.EQ.0)GO TO 22
C-----Change the U thermal parameters
      EMIN=1E-5
      IF(B(2).LT.1E-5.AND.B(3).LT.1E-5.AND.B(4).LT.1E-5.AND.
     1 B(5).LT.1E-5.AND.B(6).LT.1E-5)GO TO 5
      DO 11 I=1,3
      U(I,I)=B(I)
11    SU(I,I)=SB(I)
      U(2,3)=B(4)
      U(3,2)=U(2,3)
      U(1,3)=B(5)
      U(3,1)=U(1,3)
      U(1,2)=B(6)
      U(2,1)=U(1,2)
      SU(2,3)=SB(4)
      SU(3,2)=SU(2,3)
      SU(1,3)=SB(5)
      SU(3,1)=SU(1,3)
      SU(1,2)=SB(6)
      SU(2,1)=SU(1,2)
      DO 16 I=1,3
      DO 16 J=1,3
16    RHD(I,J)=RH(I,J)
      CALL INVMAT3(RHD,RHI,IER)
      IF(IER.EQ.1)RETURN
      KKD=0
      CALL ROTU(KKD,RHI,U,U1)
      KKD=1
      CALL ROTU(KKD,RHI,SU,SU1)
      DO 12 I=1,3
      B(I)=U1(I,I)
12    SB(I)=ABS(SU1(I,I))
      B(4)=U1(3,2)
      B(5)=U1(3,1)
      B(6)=U1(2,1)
      SB(4)=ABS(SU1(3,2))
      SB(5)=ABS(SU1(3,1))
      SB(6)=ABS(SU1(2,1))
22    IF(KOR.EQ.0)THEN
      IF(KL.EQ.1.OR.KL.EQ.2)WRITE(IO,FRT)ATOM,X1,SX1
      IF(KL.EQ.3)WRITE(IO,FRT)ATOM,(X1(J),SX1(J),J=1,3)
      IF(KL.EQ.4)THEN
      DO 15 I=1,3
      SS=SX1(I)
      CALL CONVN(SS,II)
      DO 15 JJ=1,5
15    UNIM(I,JJ)=II(JJ)
      WRITE(IO,100)ATOM,X1,((UNIM(I,K),K=1,5),I=1,3)
      ENDIF
      ELSE
      IF(KL.EQ.1)WRITE(IO,FRT)ATOM,X1,SX1,B,SB
      IF(KL.EQ.2)WRITE(IO,FRT)ATOM,X1,B,SX1,SB
      IF(KL.EQ.3)WRITE(IO,FRT)ATOM,(X1(J),SX1(J),J=1,3),(B(K),
     1 SB(K),K=1,6)
      ENDIF
      GO TO 4
5     IF(KL.EQ.1)WRITE(IO,FRT)ATOM,X1,SX1,B,SB
      IF(KL.EQ.2)WRITE(IO,FRT)ATOM,X1,B,SX1,SB
      IF(KL.EQ.3)WRITE(IO,FRT)ATOM,(X1(J),SX1(J),J=1,3),(B(K),
     1 SB(K),K=1,6)
4     CONTINUE
13    CONTINUE
      READ(IN,'(80A1)',END=14)TIT
      WRITE(IO,'(80A1)')TIT
      GO TO 13
14    RETURN
      END
 
      SUBROUTINE CONVE(II,DX)
      DIMENSION IX(5)
      CHARACTER*1 II(5),IA(10)
      DATA IA/'0','1','2','3','4','5','6','7','8','9'/
      DO 2 I=1,5
      IX(I)=0
      DO 3 K=1,10
      IF(IA(K).EQ.II(I)) GO TO 1
3     CONTINUE
1     IX(I)=K-1
2     CONTINUE
      IF(II(5).NE.' ') DX=IX(1)*0.1+IX(2)*0.01+IX(3)*0.001+
     1 IX(4)*0.0001+IX(5)*0.00001
      IF(II(5).EQ.' ') DX=IX(1)*0.01+IX(2)*0.001+IX(3)*0.0001+
     1 IX(4)*0.00001
      IF(II(4).EQ.' ') DX=IX(1)*0.001+IX(2)*0.0001+IX(3)*0.00001
      IF(II(3).EQ.' ') DX=IX(1)*0.0001+IX(2)*0.00001
      IF(II(2).EQ.' ') DX=IX(1)*0.00001
      IF(II(1).EQ.' ') DX=0.0
      RETURN
      END
 
      SUBROUTINE CONVN(SX,NC)
      CHARACTER*1 NC(5),NUM(10)
      INTEGER*4 NA,IS,IA(5)
      DATA NUM/'0','1','2','3','4','5','6','7','8','9'/
      NA=SX*100000
      IA(1)=NA/10000
      IS=IA(1)*10
      IA(2)=NA/1000-IS
      IS=IA(1)*100+IA(2)*10
      IA(3)=NA/100-IS
      IS=IA(1)*1000+IA(2)*100+IA(3)*10
      IA(4)=NA/10-IS
      IS=IA(1)*10000+IA(2)*1000+IA(3)*100+IA(4)*10
      IA(5)=NA-IS
      DO 2 I=1,5
      DO 1 J=1,10
      IF(IA(I).NE.(J-1))GO TO 1
      NC(I)=NUM(J)
      GO TO 2
1     CONTINUE
2     CONTINUE
      RETURN
      END
 
      SUBROUTINE NEWCELL(TAU,STAU,AN,SAN)
      DIMENSION C(3),SC(3),A(3),SA(3),AN(3),SAN(3),
     1 CS(3),SN(3),TAU(3),STAU(3)
      REAL IND(2,3)
      COMMON/IOU/IN,IO
     1      /GRP1/RH(3,3)
     2      /GRP2/RX(3,3),TX(3)
     3      /GRP3/C,SC,A,SA
      DG=0.01745329
      DO 6 I=1,3
      CS(I)=COS(A(I)*DG)
      SN(I)=SIN(A(I)*DG)
6     SA(I)=SA(I)*DG
      DO 2 I=1,3
      TAU(I)=0.0
      STAU(I)=0.0
      DO 1 K=1,3
1     TAU(I)=TAU(I)+(RH(I,K)*C(K))**2
      TAU(I)=TAU(I)+2.0*(RH(I,1)*RH(I,2)*C(1)*C(2)*CS(3)+
     1 RH(I,2)*RH(I,3)*C(2)*C(3)*CS(1)+
     2 RH(I,1)*RH(I,3)*C(1)*C(3)*CS(2))
      TAU(I)=SQRT(ABS(TAU(I)))
      STAU(I)=((RH(I,1)**2*C(1)+RH(I,1)*RH(I,2)*C(2)*CS(3)+
     1 RH(I,1)*RH(I,3)*C(3)*CS(2))*SC(1))**2
      STAU(I)=STAU(I)+((RH(I,2)**2*C(2)+RH(I,1)*RH(I,2)*C(1)
     1 *CS(3)+RH(I,2)*RH(I,3)*C(3)*CS(1))*SC(2))**2
      STAU(I)=STAU(I)+((RH(I,3)**2*C(3)+RH(I,2)*RH(I,3)*C(2)*CS(1)
     1+RH(I,3)*RH(I,1)*C(1)*CS(2))*SC(3))**2
      STAU(I)=STAU(I)+(RH(I,2)*RH(I,3)*C(2)*C(3)*SN(1)*SA(1))**2
      STAU(I)=STAU(I)+(RH(I,3)*RH(I,1)*C(3)*C(1)*SN(2)*SA(2))**2
      STAU(I)=STAU(I)+(RH(I,1)*RH(I,2)*C(1)*C(2)*SN(3)*SA(3))**2
2     STAU(I)=SQRT(ABS(STAU(I)))/TAU(I)
      DO 3 I=1,3
      IND(1,I)=RH(2,I)
3     IND(2,I)=RH(3,I)
      T1=TAU(2)
      T2=TAU(3)
      CALL ANGLE(IND,C,SC,SN,CS,SA,T1,T2,TH,STH)
      AN(1)=TH
      SAN(1)=ABS(STH)
      DO 4 I=1,3
4     IND(1,I)=RH(1,I)
      T1=TAU(1)
      CALL ANGLE(IND,C,SC,SN,CS,SA,T1,T2,TH,STH)
      AN(2)=TH
      SAN(2)=ABS(STH)
      DO 5 I=1,3
5     IND(2,I)=RH(2,I)
      T2=TAU(2)
      CALL ANGLE(IND,C,SC,SN,CS,SA,T1,T2,TH,STH)
      AN(3)=TH
      SAN(3)=ABS(STH)
      RETURN
      END
 
      SUBROUTINE ANGLE(IND,C,SC,SN,CS,SA,T1,T2,TH,STH)
C-----Calculates the angle formed by two lattice vectors with e.s.d.
      DIMENSION C(3),SC(3),SN(3),CS(3),SA(3)
      REAL IND(2,3)
      CSTH=IND(1,1)*IND(2,1)*C(1)**2+IND(1,2)*IND(2,2)*C(2)**2
     1 +IND(1,3)*IND(2,3)*C(3)**2
      CSTH=CSTH+(IND(1,1)*IND(2,2)+IND(2,1)*IND(1,2))*C(1)*C(2)*CS(3)
      CSTH=CSTH+(IND(1,2)*IND(2,3)+IND(2,2)*IND(1,3))*C(2)*C(3)*CS(1)
      CSTH=CSTH+(IND(1,3)*IND(2,1)+IND(2,3)*IND(1,1))*C(3)*C(1)*CS(2)
      IF(ABS(CSTH).LT.1E-6)CSTH=0.0
      CSTH=CSTH/(T1*T2)
      SNTH=SQRT(ABS(1-CSTH**2))
      TH=ATAN2(SNTH,CSTH)
      STH1=2.0*IND(1,1)*IND(2,1)*C(1)+
     1 (IND(1,1)*IND(2,2)+IND(2,1)*IND(1,2))*C(2)*CS(3)
      STH1=((STH1+(IND(1,3)*IND(2,1)+IND(2,3)*IND(1,1))*C(3)*CS(2))*
     1 SC(1))**2
      STH2=2.0*IND(1,2)*IND(2,2)*C(2)+
     1 (IND(1,1)*IND(2,2)+IND(2,1)*IND(1,2))*C(1)*CS(3)
      STH2=((STH2+(IND(1,2)*IND(2,3)+IND(2,2)*IND(1,3))*C(3)*CS(1))*
     1 SC(2))**2
      STH3=2.0*IND(1,3)*IND(2,3)*C(3)+
     1 (IND(1,2)*IND(2,3)+IND(2,2)*IND(1,3))*C(2)*CS(1)
      STH3=((STH3+(IND(1,3)*IND(2,1)+IND(2,3)*IND(1,1))*C(1)*CS(2))*
     1 SC(3))**2
      STH=((IND(1,2)*IND(2,3)+IND(2,2)*IND(1,3))*C(2)*C(3)*SN(1)
     1 *SA(1))**2
      STH=STH+((IND(1,3)*IND(2,1)+IND(2,3)*IND(1,1))*C(3)*C(1)*SN(2)
     1 *SA(2))**2
      STH=STH+((IND(1,1)*IND(2,2)+IND(2,1)*IND(1,2))*C(1)*C(2)*SN(3)
     1 *SA(3))**2
      STH=SQRT(ABS(STH1+STH2+STH3+STH))/(T1*T2*SNTH)
      RD=57.29578
      TH=TH*RD
      STH=STH*RD
      RETURN
      END
 
      SUBROUTINE NEWSHELX
C-----This routine trasforms the input coordinates and thermal parameters
C-----for SHELXL-93 according to a rotation matrix and a translation
C-----vector
      CHARACTER*1 TIT(80),STR(10,10),MINUS(4)
      CHARACTER*4 CODE(61),CDD,CD1
      DIMENSION NC(10),A1(10)
      DIMENSION C(3),SC(3),A(3),SA(3),X(3),B(6),TAU(3),
     1 STAU(3),AN(3),SAN(3),X1(3)
      DOUBLE PRECISION U(3,3),RHI(3,3),U1(3,3),RHD(3,3),RDXI(3,3)
      COMMON/IOU/IN,IO
     1      /GRP1/RH(3,3)
     2      /GRP2/RX(3,3),TX(3)
     3      /GRP3/C,SC,A,SA
      EQUIVALENCE (TIT(1),CDD)
      DATA CODE/'ACTA','AFIX','ANIS','ATOM','BASF','BIND','BLOC',
     1 'BOND','CHIV','CONF','CONN','DISP','DAMP','DELU','DFIX',
     2 'EADP','EQIV','EXTI','END ','EXYZ','FEND','FLAT','FMAP','FRAG',
     3 'FREE','FVAR','GRID','HFIX','HKLF','ISOR','LATT','LAUE','LIST',
     9 'L.S.','MERG','MORE','MOVE','MPLA','MOLE','OMIT','PART','PLAN',
     5 'REM ','RESI','RTAB','SADI','SAME','SFAC','SHEL','SIMU','SIZE',
     6 'SUMP','SYMM','TEMP','TIME','TITL','TWIN','UNIT','VFIL','VOID',
     7 'WGHT'/
      KKK=0
      READ(IN,'(80A1)')TIT
      WRITE(IO,'(80A1)')TIT
7     READ(IN,'(80A1)',END=6)TIT
      DO 35 I=1,80
      IF(TIT(I).NE.' ')GO TO 40
35    CONTINUE
      WRITE(IO,'(80A1)')TIT
      GO TO 7
40    DO 24 I=1,4
24    MINUS(I)=TIT(I)
      NCS=4
      CALL MINUMA(NCS,MINUS)
      DO 25 I=1,4
25    TIT(I)=MINUS(I)
      IF(CDD.NE.'CELL')GO TO 13
      CALL DECODSTR(TIT,STR,NC,NS,A1,KEQ)
      AL=A1(1)
      DO 10 I=2,4
      C(I-1)=A1(I)
10    A(I-1)=A1(I+3)
      GO TO 7
13    IF(CDD.NE.'ZERR')GO TO 14
      CALL DECODSTR(TIT,STR,NC,NS,A1,KEQ)
      IZ=A1(1)
      DO 11 I=2,4
      SC(I-1)=A1(I)
11    SA(I-1)=A1(I+3)
      DO 12 I=1,3
12    IF(SC(I).LE.1E-9)SC(I)=C(I)*1E-9
      DO 36 I=1,3
      DO 36 J=1,3
      IF(ABS(RH(I,I)-1.0).GT.0.01) GO TO 38
      IF(I.NE.J.AND.ABS(RH(I,J)).GT.0.01)GO TO 38
36    CONTINUE
      DO 37 I=1,3
      TAU(I)=C(I)
      STAU(I)=SC(I)
      AN(I)=A(I)
37    SAN(I)=SA(I)
      GO TO 39
38    CALL NEWCELL(TAU,STAU,AN,SAN)
39    WRITE(IO,100)AL,TAU,AN,IZ,STAU,SAN
100   FORMAT('CELL',1X,F7.5,3(1X,F8.4),3(1X,F7.3)/'ZERR',1X,I3,
     1 3(1X,F6.4),3(1X,F5.3))
      GO TO 7
14    IF(KKK.EQ.1) CALL CHANGRH(KKK)
      DO 8 I=1,61
      IF(CDD.NE.CODE(I))GO TO 8
      WRITE(IO,'(80A1)')TIT
      GO TO 7
8     CONTINUE
      IF(CDD.EQ.'    ')GO TO 15
      CD1=CDD
      CALL DECODSTR(TIT,STR,NC,NS,A1,KEQ)
      IA=A1(1)
C-----Change the x,y,z coordinates
      DO 21 I=2,4
21    X(I-1)=A1(I)
      DO 17 I=1,3
      X1(I)=0.0
      DO 16 J=1,3
16    X1(I)=X1(I)+RX(I,J)*X(J)
17    X1(I)=X1(I)+TX(I)
      SOF=A1(5)
      B(1)= A1(6)
      IF(NS.GT.6)B(2)=A1(7)
      IF(KEQ.EQ.1)GO TO 22
      WRITE(IO,'(A4,1X,I4,5(1X,F9.5))')CDD,IA,X1,SOF,B(1)
      GO TO 7
22    B(2)=A1(7)
      GO TO 7
15    CALL DECODSTR(TIT,STR,NC,NS,A1,KEQ)
      DO 23 I=1,4
23    B(I+2)=A1(I)
C-----Change the U thermal parameters
      DO 18 I=1,3
18    U(I,I)=B(I)
      U(2,3)=B(4)
      U(3,2)=U(2,3)
      U(1,3)=B(5)
      U(3,1)=U(1,3)
      U(1,2)=B(6)
      U(2,1)=U(1,2)
      DO 19 I=1,3
      DO 19 J=1,3
19    RHD(I,J)=RH(I,J)
      CALL INVMAT3(RHD,RHI,IER)
      IF(IER.EQ.1)RETURN
      KKD=0
      CALL ROTU(KKD,RHI,U,U1)
      DO 20 I=1,3
20    B(I)=U1(I,I)
      B(4)=U1(3,2)
      B(5)=U1(3,1)
      B(6)=U1(2,1)
      WRITE(IO,'(A4,1X,I4,6(F10.5),1X,''=''/5X,4(F10.5))')CD1,IA,X1,
     1 SOF,B
9     GO TO 7
6     RETURN
      END
 
      SUBROUTINE DECODSTR(TIT,STR,NC,NS,A1,KEQ)
C-----Decode an alphanumeric string
      CHARACTER*1 TIT(80),STR(10,10),STR1(10)
      DIMENSION NC(10),A1(10)
      KD=0
      J=0
      KEQ=0
      DO 1 I=5,80
      IF(KD.EQ.0)NN=0
      IF(TIT(I).EQ.' '.AND.TIT(I+1).NE.' ')THEN
      J=J+1
      DO 3 K=1,10
3     STR(J,K)=' '
      ENDIF
      IF(TIT(I).EQ.' ')GO TO 2
      NN=NN+1
      STR(J,NN)=TIT(I)
      KD=1
      NC(J)=NN
      GO TO 1
2     KD=0
1     CONTINUE
      NS=J-1
      DO 5 I=1,NS
      DO 4 J=1,NC(I)
4     STR1(J)=STR(I,J)
      IF(STR1(1).EQ.'=')GO TO 6
      NC1=NC(I)
      CALL CHARNUM(NC1,STR1,AN)
5     A1(I)=AN
      GO TO 7
6     KEQ=1
7     RETURN
      END
 
      SUBROUTINE CHARNUM(NC,STR,AA)
C-----Trasforma una stringa alfanumerica in un numero decimale
      CHARACTER*1 STR(NC),NUM(10)
      DOUBLE PRECISION AN,AB
      DATA NUM/'0','1','2','3','4','5','6','7','8','9'/
      IF(NC.NE.1)GO TO 9
      DO 7 I=1,10
      IF(STR(1).NE.NUM(I))GO TO 7
      AA=I-1
      GO TO 8
7     CONTINUE
8     RETURN
9     IS=1
      N1=1
      IF(STR(1).EQ.'-')THEN
      IS=-1
      N1=2
      ENDIF
      DO 1 I=1,NC
      IF(STR(I).NE.'.')GO TO 1
      GO TO 2
1     CONTINUE
2     N=I
      AN=0.0
      DO 4 I=N1,N-1
      DO 3 J=1,10
      IF(STR(I).NE.NUM(J))GO TO 3
      IA=10**(N-I-1)
      A=(J-1)*IA
      AN=AN+A
      GO TO 4
3     CONTINUE
4     CONTINUE
      DO 6 I=N+1,NC
      DO 5 J=1,10
      IF(STR(I).NE.NUM(J))GO TO 5
      A=J-1
      B=10**(I-N)
      AB=A/B
      AN=AN+AB
      GO TO 6
5     CONTINUE
6     CONTINUE
      AN=AN*IS
      AA=AN
      RETURN
      END
 
      SUBROUTINE MINUMA(N,STR)
C-----This routine transforms the lower case characters into upper
C-----case characters, in a string of N characters
      CHARACTER*1 MINU(26),MAIU(26),STR(N)
      DATA MINU/'a','b','c','d','e','f','g','h','i','j','k','l','m',
     1 'n','o','p','q','r','s','t','u','v','w','x','y','z'/
      DATA MAIU/'A','B','C','D','E','F','G','H','I','J','K','L','M',
     1 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
      DO 2 I=1,N
      DO 1 J=1,26
      IF(STR(I).NE.MINU(J)) GO TO 1
      STR(I)=MAIU(J)
      GO TO 2
1     CONTINUE
2     CONTINUE
      RETURN
      END
 
      SUBROUTINE CHANGRH(KKK)
C-----This routine changes the RH matrix
      DOUBLE PRECISION RDX(3,3),RDXI(3,3)
      COMMON/IOU/IN,IO
     1      /GRP1/RH(3,3)
     2      /GRP2/RX(3,3),TX(3)
      KKK=0
      DO 1 I=1,3
      DO 1 J=1,3
      IF(ABS(RH(I,I)-1.0).GT.0.01) GO TO 4
      IF(I.NE.J.AND.ABS(RH(I,J)).GT.0.01)GO TO 4
1     CONTINUE
      DO 2 I=1,3
      DO 2 J=1,3
2     RDX(I,J)=RX(I,J)
      CALL INVMAT3(RDX,RDXI,IER)
      DO 3 I=1,3
      DO 3 J=1,3
3     RH(J,I)=RDXI(I,J)
4     KKK=1
      RETURN
      END
