      Program ROTMOL
C**********************************************************************
C* January 1998                                                       *
C*                      Program ROTMOL                                *
C*                                                                    *
C*                      by M. Nardelli                                *
C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica  *
C* Chimica Fisica della Universita' di Parma, Centro di Studio per la *
C* Strutturistica Diffrattometrica del CNR.                           *
C*           Viale delle Scienze 78, I-43100 Parma, Italy.            *
C*                                                                    *
C* This routine calculates the coordinates of the atoms of a molecule *
C* or a molecular fragment after rotation about a direction defined   *
C* by the first two atoms. It can be used in connection with problems *
C* of rotational disorder.                                            *
C*                                                                    *
C*                         INPUT                                      *
C* 1. TITLE (A80)                                                     *
C* 2. a, b, c, alpha, beta, gamma (free FORMAT)                       *
C* 3. N (total number of atoms, max.200), FI (rotation angle, deg.)   *
C* 4. Format of coordinates [e.g.(6A1,3X,3F8.5)]                      *
C* 5. Atom, x, y, z  (max number of atoms: 200)                       *
C*                                                                    *
C*                         OUTPUT                                     *
C* Title                                                              *
C* Crystal data                                                       *
C* Orthogonalization matrix                                           *
C* A Table giving the atomic coordinates of the atoms of the un-      *
C* rotated and rotated molecule                                       *
C*                                                                    *
C*                      SUBROUTINES                                   *
C* ORTHOG(C,A,RD)                                                     *
C* NORM(X)                                                            *
C* ROTFRG(N,NG,X,FI,XO)                                               *
C* INVMAT3(B,BI)                                                      *
C**********************************************************************
      CHARACTER*1 TIT(80),FRT(80),AT(200,6)
      CHARACTER*14 FILIN,FILOUT
      DIMENSION C(3),A(3),X(200,3),XO(200,3),XRO(200,3),XR(200,3),
     1 O(3,3),OI(3,3)
      COMMON/GRP1/AT
     1      /GRP2/O,OI
      RD=0.0174532925
C-----Assign input and output files
      WRITE(*,500)
500   FORMAT(' Key in the input file'/' ?>')
      READ(*,'(A14)')FILIN
      WRITE(*,501)
501   FORMAT(' Key in the output file'/' ?>')
      READ(*,'(A14)')FILOUT
      OPEN(UNIT=13,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=13)
      OPEN(UNIT=15,FILE=FILOUT,FORM='FORMATTED')
      REWIND(UNIT=15)
C-----Read title, codes and parameters
      READ(13,100)TIT
100   FORMAT(80A1)
      WRITE(15,101)TIT
101   FORMAT(1X,80A1/)
      READ(13,*)C,A
      WRITE(15,102)(C(I),A(I),I=1,3)
102   FORMAT(' Crystal data'/' a = ',F7.4,25X,'alpha= ',F6.2/
     1 ' b = ',F7.4,25X,'beta = ',F6.2/' c = ',F7.4,25X,'gamma= ',F6.2)
      CALL ORTHOG(C,A,RD)
      WRITE(15,103)((O(I,J),J=1,3),I=1,3)
103   FORMAT(/' Orthogonalization matrix'/3(1X,3F10.5/))
      READ(13,*)N,FI
      FI=FI*RD
      READ(13,100) FRT
      WRITE(15,104)
104   FORMAT(/' Atomic coordinates'//'Atom ',5X,'X/a',6X,'Y/b',6X,
     1 'Z/c',12X,'X',8X,'Y',8X,'Z'//' before rotation:')
      DO 2 I=1,N
      READ(13,FRT)(AT(I,J),J=1,6),(X(I,J),J=1,3)
C-----Orthogonal coordinates
      DO 1 J=1,3
      XO(I,J)=0.0
      DO 1 K=1,3
1     XO(I,J)=XO(I,J)+O(J,K)*X(I,K)
2     WRITE(15,105)(AT(I,L),L=1,6),(X(I,J),J=1,3),(XO(I,K),K=1,3)
105   FORMAT(6A1,1X,3F9.5,5X,3F9.4)
      CALL ROTFRG(N,XO,FI,XRO)
      DO 3 I=1,N
      DO 3 J=1,3
      XR(I,J)=0.0
      DO 3 K=1,3
3     XR(I,J)=XR(I,J)+OI(J,K)*XRO(I,K)
      FI=FI/RD
      WRITE(15,106)(AT(1,L),L=1,6),(AT(2,L),L=1,6),FI
106   FORMAT(/' after rotation about ',6A1,'-',6A1,' by',F6.1,' degr.')
      DO 7 I=1,N
 
      DO 5 L=3,6
      IF(AT(I,L).EQ.' ')THEN
      AT(I,L)=''''
      GO TO 6
      ENDIF
5     CONTINUE
      WRITE(15,107)(AT(I,L),L=1,6),(XR(I,K),K=1,3),(XRO(I,K),K=1,3)
107   FORMAT(6A1,'''',1X,3F9.5,5X,3F9.5)
      GO TO 7
6     WRITE(15,108)(AT(I,L),L=1,6),(XR(I,K),K=1,3),(XRO(I,K),K=1,3)
108   FORMAT(6A1,2X,3F9.5,5X,3F9.5)
7     CONTINUE
      STOP
      END
 
      SUBROUTINE ORTHOG(C,A,RD)
C-----This routine calculates the orthogonalization matrix O(3,3)
      DIMENSION C(3),A(3),SN(3),CS(3),O(3,3),OI(3,3)
      COMMON/GRP2/O,OI
      DO 1 I=1,3
      AR=A(I)*RD
      SN(I)=SIN(AR)
1     CS(I)=COS(AR)
      CSAR=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3))
      SNAR=SQRT(1.0-CSAR**2)
C-----Calculation of the orthogonalization matrix
      O(1,1)=C(1)
      O(1,2)=C(2)*CS(3)
      O(1,3)=C(3)*CS(2)
      O(2,1)=0.0
      O(2,2)=C(2)*SN(3)
      O(2,3)=-C(3)*SN(2)*CSAR
      O(3,1)=0.0
      O(3,2)=0.0
      O(3,3)=C(3)*SN(2)*SNAR
C-----Calculation of the inverted orthogonalization matrix
      OI(1,1)=1./O(1,1)
      OI(1,2)=-O(1,2)/(O(1,1)*O(2,2))
      OI(1,3)=(O(1,2)*O(2,3)-O(2,2)*O(1,3))/(O(1,1)*O(2,2)*O(3,3))
      OI(2,1)=0.0
      OI(2,2)=1./O(2,2)
      OI(2,3)=-O(2,3)/(O(2,2)*O(3,3))
      OI(3,1)=0.0
      OI(3,2)=0.0
      OI(3,3)=1./O(3,3)
      RETURN
      END
 
      SUBROUTINE NORM(X)
C-----This routine normalizes a 3-dimensional vector
      DIMENSION X(3)
      S=0.0
      DO 1 I=1,3
1     S=S+X(I)**2
      S=SQRT(S)
      DO 2 I=1,3
2     X(I)=X(I)/S
      RETURN
      END
 
      SUBROUTINE ROTFRG(N,X,FI,XO)
C-----This routine modifies the coordinates so to have them referred
C-----to a system of orthogonal axes with the origin at the first atom,
C-----the z axis directed from the first to the second atom, the y
C-----axis normal to z and lying in the plane defined by the first
C-----three atoms and the x axis normal to y and z.
      DIMENSION X(200,3),XO(200,3),XR(3),X1(200,3),X2(200,3),B(3,3),
     1 RT(3,3),BI(3,3)
      DO 1 I=1,N
      DO 1 J=1,3
1     X1(I,J)=X(I,J)-X(1,J)
C-----Building up matrix B
      XR(1)=X1(2,1)
      XR(2)=X1(2,2)
      XR(3)=X1(2,3)
      CALL NORM(XR)
      DO 2 I=1,3
2     B(3,I)=XR(I)
      V=0
      DO 3 K=1,3
3     V=V+X1(3,K)*B(3,K)
      DO 4 K=1,3
4     XR(K)=X1(3,K)-V*B(3,K)
      CALL NORM(XR)
      DO 5 I=1,3
5     B(2,I)=XR(I)
      XR(1)=B(2,2)*B(3,3)-B(3,2)*B(2,3)
      XR(2)=B(3,1)*B(2,3)-B(2,1)*B(3,3)
      XR(3)=B(2,1)*B(3,2)-B(3,1)*B(2,2)
      CALL NORM(XR)
      DO 6 I=1,3
6     B(1,I)=XR(I)
      DO 7 I=1,N
      DO 7 J=1,3
      X2(I,J)=0.0
      DO 7 K=1,3
7     X2(I,J)=X2(I,J)+B(J,K)*X1(I,K)
C-----Building up rotation matrix
      RT(1,1)=COS(FI)
      RT(1,2)=-SIN(FI)
      RT(1,3)=0.0
      RT(2,1)=SIN(FI)
      RT(2,2)=COS(FI)
      RT(2,3)=0.0
      RT(3,1)=0.0
      RT(3,2)=0.0
      RT(3,3)=1.0
C-----Applying rotation matrix
      DO 8 I=1,N
      DO 8 J=1,3
      XO(I,J)=0.0
      DO 8 K=1,3
8     XO(I,J)=XO(I,J)+RT(J,K)*X2(I,K)
      CALL INVMAT3(B,BI)
      DO 9 I=1,N
      DO 9 J=1,3
      X1(I,J)=0.0
      DO 9 K=1,3
9     X1(I,J)=X1(I,J)+BI(J,K)*XO(I,K)
      DO 10 I=1,N
      DO 10 J=1,3
10    XO(I,J)=X1(I,J)+X(1,J)
      RETURN
      END
 
      SUBROUTINE INVMAT3(B,BI)
C-----This routine inverts a 3x3 matrix
      DIMENSION B(3,3),BI(3,3)
      BI(1,1)=B(2,2)*B(3,3)-B(2,3)*B(3,2)
      BI(2,1)=B(3,1)*B(2,3)-B(2,1)*B(3,3)
      BI(3,1)=B(2,1)*B(3,2)-B(3,1)*B(2,2)
      BI(1,2)=B(3,2)*B(1,3)-B(1,2)*B(3,3)
      BI(2,2)=B(1,1)*B(3,3)-B(3,1)*B(1,3)
      BI(3,2)=B(3,1)*B(1,2)-B(1,1)*B(3,2)
      BI(1,3)=B(1,2)*B(2,3)-B(2,2)*B(1,3)
      BI(2,3)=B(2,1)*B(1,3)-B(1,1)*B(2,3)
      BI(3,3)=B(1,1)*B(2,2)-B(1,2)*B(2,1)
      RETURN
      END
