C**************************************************************************
C     Program to calculate differential rms-displacements. Uses input
C     from an edited ortep output file(UNIT 1) containing cartesian 
C     coordinates of atoms, principal rmsd's(Ui) and direction cosines. 
C     The U-matrix is trasformed from diagonal form to cartesian frame. 
C     The rmsd along the bond(from the first atom in the list) is calculated.
C     The difference between the rmsd-values for the two atoms forming 
C     the bond is printed.        
C*****************************mvr/december '96*****************************
C     SC = SCALE FACTOR (overall scale = ....inch/angstrom)
C     IF SC < 0 MORE OUTPUT IN UNIT 2
C     MN = CENTRAL (FIRST) ATOM NAME; AN = ATOM NAME 
C     U(cart.) = T'U(prin.)T ;    1          X
C                                 2  =  (T)  Y
C                                 3          Z
C                            (Principal   (Cartesian
C                              frame)        frame)
C     D = V'U(cart.)V , V' = unit row vector along the M-A bond in cartesian.
C         PM or PA = U(prin.) 
C         UM or UA = U(cart.) 
C         Q = U(prin.)T
C     OUTPUT TABLE IN FILE "DELU"
C*****
	CHARACTER*6 MN,AN
	DIMENSION TM(3,3),Q(3,3),PM(3,3),UM(3,3),PA(3,3),UA(3,3),V(3),
     1  TA(3,3),AM(3),AA(3)
        CON=57.29578
        OPEN(3,FILE='DELU')
        WRITE(3,60)
        READ(1,*) SC
        WRITE(2,50) SC
	READ(1,10) MN,XM,YM,ZM 
	WRITE(2,40) MN,XM,YM,ZM 
	DO 11 I=1,3
	DO 11 J=1,3
   11	PM(I,J)=0.
	READ(1,*)   ((TM(I,J),J=1,3),PM(I,I),I=1,3)
	WRITE(2,20) ((TM(I,J),J=1,3),PM(I,I),I=1,3)
        DO 16 I=1,3
   16   PM(I,I)=PM(I,I)*PM(I,I)
	CALL MULT(PM,TM,Q)
        IF(SC.LT.0)CALL OUTM(Q)
	CALL TRAN(TM)
        IF(SC.LT.0)CALL OUTM(TM)
	CALL MULT(TM,Q,UM)
	CALL OUTM(UM)
   12	READ(1,10,END=100) AN,XA,YA,ZA 
	WRITE(2,40) AN,XA,YA,ZA 
	DO 13 I=1,3
	DO 13 J=1,3
   13   PA(I,J)=0.
	READ(1,*)   ((TA(I,J),J=1,3),PA(I,I),I=1,3)
	WRITE(2,20) ((TA(I,J),J=1,3),PA(I,I),I=1,3)
        DO 17 I=1,3
   17   PA(I,I)=PA(I,I)*PA(I,I)
	CALL MULT(PA,TA,Q)
        IF(SC.LT.0)CALL OUTM(Q)
	CALL TRAN(TA)
        IF(SC.LT.0)CALL OUTM(TA)
	CALL MULT(TA,Q,UA)
	CALL OUTM(UA)
	V1=XM-XA
	V2=YM-YA
	V3=ZM-ZA
	D=SQRT(V1*V1+V2*V2+V3*V3)
	V(1)=V1/D
	V(2)=V2/D
	V(3)=V3/D
        DI=D/ABS(SC)
        WRITE(2,50) DI
        CALL OUTV(V)
        DO 14 I=1,3
        AM(I)=CON*ACOS(TM(I,1)*V(1)+TM(I,2)*V(2)+TM(I,3)*V(3))        
        AA(I)=CON*ACOS(TA(I,1)*V(1)+TA(I,2)*V(2)+TA(I,3)*V(3))
   14   CONTINUE
        WRITE(2,80)MN,(PM(I,I),AM(I),I=1,3)
        WRITE(2,80)AN,(PA(I,I),AA(I),I=1,3)
	CALL UBON(V,UM,DM)
	CALL UBON(V,UA,DA)
	DIFF=DA-DM
        DIF=SQRT(ABS(DIFF))
        IF (DIFF.LT.0.)DIF=-DIF
	WRITE(2,30) DM,DA,DIF
        WRITE(3,70) MN,DM,AN,DA,DIF,DI,(INT(AM(I)),I=1,3),
     1  (INT(AA(I)),I=1,3)
	GOTO 12
   10 FORMAT(/11X,A6,33X,3F10.6)
   20 FORMAT(5X,3F10.4,2X,F10.6)
   30 FORMAT(5X,'U-m = ',F10.4,1X,
     1'U-a = ',F10.4,2X,'DIFD = ',F10.4) 
   40 FORMAT(5X,A6,3F10.6)
   50 FORMAT(5X,F10.6)
   60 FORMAT(11X,' U/A2  ',8X,' U/A2   ', 'DEL/A', '    DIST/A',
     1'         ANGLES'/3X,'-------',
     2'---------------------------------------------------------------')  
   70 FORMAT(3X,A6,F7.4,2X,A6,F7.4,1X,F7.4,2X,F8.4,6I4)
   80 FORMAT(5X,A6,6F10.4)
  100 STOP
      END   
C*****
      SUBROUTINE MULT(A,B,C)
	DIMENSION A(3,3),B(3,3),C(3,3)
	N=3
	DO 30 I=1,N
	DO 20 J=1,N
	S=0.
	DO 10 K=1,N
	S=S+A(I,K)*B(K,J)
   10 CONTINUE
      C(I,J)=S
   20 CONTINUE
   30 CONTINUE		
      RETURN
	END
C*****
      SUBROUTINE TRAN(A)
	DIMENSION A(3,3)
	N=3
	DO 20 I=1,N
	DO 10 J=I+1,N
	T=A(I,J)
	A(I,J)=A(J,I)
	A(J,I)=T
   10 CONTINUE
   20 CONTINUE
      RETURN
	END
C*****
      SUBROUTINE UBON(V,A,S)
	DIMENSION V(3),A(3,3),T(3)
	N=3
	DO 20 I=1,N
	S=0.
	DO 10 K=1,N
	S=S+V(K)*A(K,I)
   10 CONTINUE
      T(I)=S
   20 CONTINUE		 		
      S=0.
	DO 30 I=1,N
	S=S+T(I)*V(I)
   30 CONTINUE
      RETURN
	END
C*****		
      SUBROUTINE OUTM(A)
	DIMENSION A(3,3)
	N=3
	DO 11 I=1,N
	WRITE(2,10) (A(I,J),J=1,N)
   11 CONTINUE	
   10 FORMAT(5X,3(F10.4,1X))	   
      RETURN
	END
C*****
      SUBROUTINE OUTV(A)
      DIMENSION A(3)
      N=3
      WRITE(2,10) (A(I),I=1,N)
   10 FORMAT(5X,3(F10.4,1X))
      RETURN
      END