      PROGRAM OHCOOR
C************************************************************************
C*  January 1998                                                        *
C*                           Program OHCOOR                             *
C*                                                                      *
C*                         by Mario NARDELLI                            *
C*  Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica,  * della Universita` degli Studi di
C*  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 calculates the coordinates of an atom H in a hydroxyl  *
C*  group A-B-O-H  (1-2-3-4) with the O-H bond forming a torsion angle  *
C*  Tau with respect to the A-B bond at a requested angle B-O-H (Theta) *
C*  and distance O-H (D34).                                             *
C*                               INPUT                                  *
C* - Title (80A1)                                                       *
C* - a,b,c,alpha,beta,gamma (free format)                               *
C* - 'ATOM(1)',x/a,y/b,z/c     (free format; ATOM(1):A4)                *
C* - 'ATOM(2)',x/a,y/b,z/c     (free format; ATOM(2):A4)                *
C* - 'ATOM(3)',x/a,y/b,z/c     (free format; ATOM(3):A4)                *
C* - 'ATOM(4)'                 (free format; ATOM(4):A4)                *
C* - 0,0,0  (atom1)                                                     *
C* - distance (atom1-atom2),0,0                                         *
C* - distance (atom2-atom3), bond angle (atom1-atom2-atom3),0           *
C* - distance (atom3-atom4), bond angle (atom2-atom3-atom4),torsion     *
C*   angle (atom1-atom2-atom3-atom4)                                    *
C*                                                                      *
C*                               OUTPUT                                 *
C*   Table of fractional coordinates                                    *
C*   Table of bond distances, bond angles and the torsion angle         *
C*                                                                      *
C* N.B. The sign of the torsion angle is not defined. To define the     *
C*      correct calculated coordinates of the H atom, the analisis of   *
C*      the orientation of the group OH with ROTENER1 must be carried   *
C*      out.
C*
C*   Subroutines requested:
C*                         ROTMAT(X,B)
C*                         INVMAT3(A,B)
C*                         DSTANG4(ATM,XO,RD)
C*                         OH(RD,XO)
C************************************************************************
      CHARACTER*1 TIT(80)
      CHARACTER*4 ATM(4)
      CHARACTER*14 FILIN,FILOUT
      DOUBLE PRECISION O(3,3),OI(3,3),R(3,3),RI(3,3)
      DIMENSION X(4,3),C(3),A(3),XO(4,3),SN(3),CS(3),DD1(3),X1(3),X2(3),
     1 DD2(3),XX(3),A2(3),XO1(4,3),XO2(4,3),D(2,3),X3(4,3),XO3(4,3),
     2 XO4(3,3),DX(3)
      RD=57.295779513
      WRITE(*,100)
100   FORMAT(' Key in the name of the input file'/' ?>')
      READ(*,'(A14)')FILIN
      WRITE(*,101)
101   FORMAT(' Key in the name of 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 data
      READ(13,102)TIT
102   FORMAT(80A1)
      WRITE(15,103)TIT
103   FORMAT(1X,80A1/)
      READ(13,*)C,A
      READ(13,*)(ATM(I),(X(I,K),K=1,3),I=1,3)
      READ(13,*)ATM(4)
      WRITE(15,105)(ATM(I),(X(I,K),K=1,3),I=1,3)
105   FORMAT(/1X,'Atom',9X,'x',9X,'y',9X,'z'/3(1X,A4,3X,3F10.5/))
C-----Calculation of cartesian local coordinates of the fragment
C-----from the internal coordinates
      CALL OH(RD,XO1)
C-----Calculation of the orthogonalization matrix
      DO 1 I=1,3
      AA=A(I)/RD
      SN(I)=SIN(AA)
1     CS(I)=COS(AA)
      CSAR=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3))
      SNAR=SQRT(1.0-CSAR**2)
      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-----Orthogonalization of atomic coordinates
      DO 2 I=1,3
      DO 2 J=1,3
      XO(I,J)=0.0
      DO 2 K=1,3
2     XO(I,J)=XO(I,J)+O(J,K)*X(I,K)
C-----Read the orthogonal coordinates in the local system of axes
C     READ(13,*)((XO1(I,K),K=1,3),I=1,4)
C-----Transfer of the origin
      DO 5 I=1,4
      DO 5 J=1,3
5     XO2(I,J)=XO1(I,J)+XO(1,J)
C-----Rotate the system of axes
      CALL ROTMAT(XO,R)
      DO 8 I=1,3
      DO 8 K=1,3
      XO4(I,K)=0.0
      DO 8 J=1,3
8     XO4(I,K)=XO4(I,K)+R(K,J)*XO(I,J)
      DO 3 I=1,3
3     DX(I)=XO4(1,I)-XO2(1,I)
      DO 4 I=1,4
      DO 4 J=1,3
4     XO2(I,J)=XO2(I,J)+DX(J)
      CALL INVMAT3(R,RI)
      DO 7 I=1,4
      DO 7 K=1,3
      XO3(I,K)=0.0
      DO 7 J=1,3
7     XO3(I,K)=XO3(I,K)+RI(K,J)*XO2(I,J)
C-----Calculation of fractional coordinates of H
      CALL INVMAT3(O,OI)
      DO 6 I=1,4
      DO 6 K=1,3
      X3(I,K)=0.0
      DO 6 J=1,3
6     X3(I,K)=X3(I,K)+OI(K,J)*XO3(I,J)
C-----Print fractional coordinates
      WRITE(15,104)(ATM(I),(X3(I,K),K=1,3),I=1,4)
104   FORMAT(/1X,'Atom',8X,'X/a',7X,'Y/b',7X,'Z/c'/
     1 4(1X,A4,3X,3F10.5/))
      CALL DSTANG4(ATM,XO3,RD)
      STOP
      END
 
      SUBROUTINE ROTMAT(X,B)
C-----This routine built up a matrix which transform the orthogonal
C-----coordinates of three atoms 1,2,3 in a new orthogonal system
C-----having atom 1 at the origin, x along the 1-2 direction,
C-----y orthogonal to x and in the 1,2,3 plane, z normal to the x
C-----and y axes
      DOUBLE PRECISION B(3,3),RIS
      DIMENSION X(4,3),X1(3,3),R2(3),R3(3)
      DO 1 I=1,3
      DO 1 J=1,3
1     X1(I,J)=X(I,J)-X(1,J)
      X1(2,1)=-X1(2,1)
      X1(2,2)=-X1(2,2)
      X1(2,3)=-X1(2,3)
      AA=0.0
      DO 2 K=1,3
2     AA=AA+X1(2,K)**2
      AA=SQRT(AA)
      DO 3 K=1,3
3     B(3,K)=X1(2,K)/AA
      V=0.0
      DO 4 K=1,3
4     V=V+X1(3,K)*B(3,K)
      DO 5 K=1,3
5     R2(K)=V*B(3,K)
      DO 6 K=1,3
6     R3(K)=X1(3,K)-R2(K)
      AA=0.0
      DO 7 K=1,3
7     AA=AA+R3(K)**2
      AA=SQRT(AA)
      DO 8 K=1,3
8     B(2,K)=R3(K)/AA
      B(1,1)=B(2,2)*B(3,3)-B(2,3)*B(3,2)
      B(1,2)=B(2,3)*B(3,1)-B(2,1)*B(3,3)
      B(1,3)=B(2,1)*B(3,2)-B(2,2)*B(3,1)
      DO 9 K=1,3
      RIS=-B(1,K)
      B(1,K)=B(3,K)
9     B(3,K)=RIS
      RETURN
      END
 
      SUBROUTINE INVMAT3(A,B)
C-----Inverte una matrice 3x3 con il metodo di partizione.
      DOUBLE PRECISION A(3,3),B(3,3),D,AM1,AM2,AM3,AM4,C,X
      D=A(1,1)*A(2,2)-A(1,2)*A(2,1)
      IF(DABS(D).LT.1D-10) GO TO 1
      AM1=A(1,2)*A(2,3)-A(1,3)*A(2,2)
      AM2=A(1,3)*A(2,1)-A(1,1)*A(2,3)
      AM3=A(2,1)*A(3,2)-A(2,2)*A(3,1)
      AM4=A(1,2)*A(3,1)-A(1,1)*A(3,2)
      X=A(3,3)+(AM1*A(3,1)+AM2*A(3,2))/D
      IF(DABS(X).LT.1D-10) GO TO 1
      B(3,3)=1D0/X
      X=X*D
      B(1,3)=AM1/X
      B(1,1)=(B(1,3)*AM3+A(2,2))/D
      B(1,2)=(B(1,3)*AM4-A(1,2))/D
      B(2,3)=AM2/X
      B(2,1)=(B(2,3)*AM3-A(2,1))/D
      B(2,2)=(B(2,3)*AM4+A(1,1))/D
      B(3,1)=AM3/X
      B(3,2)=AM4/X
      GO TO 2
1     WRITE(15,100)
100   FORMAT(' The matrix cannot be inverted')
      STOP
2     RETURN
      END
 
      SUBROUTINE DSTANG4(ATM,XO,RD)
C-----This routine calculates the bond distances, the bond angles
C-----and the torsion angle formed by four atoms in a A1-A2-A3-A4
C-----system
      CHARACTER*4 ATM(4)
      REAL M
      DIMENSION XO(4,3),DLT(3,3),DST(3),M(3,3)
      DO 1 J=1,3
      DLT(1,J)=XO(2,J)-XO(1,J)
      DLT(2,J)=XO(2,J)-XO(3,J)
1     DLT(3,J)=XO(4,J)-XO(3,J)
      DO 10 I=1,3
      DST(I)=0.0
      DO 2 J=1,3
2     DST(I)=DST(I)+DLT(I,J)**2
      DST(I)=SQRT(DST(I))
10    WRITE(15,100)ATM(I),ATM(I+1),DST(I)
100   FORMAT(1X,A4,'-',A4,11X,F7.3)
      DO 3 I=1,3
      DO 3 J=1,3
3     M(I,J)=DLT(I,J)/DST(I)
      CF1=M(1,1)*M(2,1)+M(1,2)*M(2,2)+M(1,3)*M(2,3)
      CF2=M(2,1)*M(3,1)+M(2,2)*M(3,2)+M(2,3)*M(3,3)
      SF1=SQRT(1.0-CF1**2)
      SF2=SQRT(1.0-CF2**2)
      P=SF1*SF2
      F1=ATAN2(SF1,CF1)*RD
      F2=ATAN2(SF2,CF2)*RD
      WRITE(15,101)ATM(1),ATM(2),ATM(3),F1,ATM(2),ATM(3),ATM(4),F2
101   FORMAT(1X,A4,'-',A4,'-',A4,3X,F10.2/1X,A4,'-',A4,'-',A4,3X,F10.2)
      IF(ABS(CF1).GT.0.98.OR.ABS(CF2).GT.0.98)GO TO 4
      IF(ABS(P).LE.1.E-5)GO TO 4
      STAU=(M(3,1)*(M(2,2)*M(1,3)-M(1,2)*M(2,3))+
     1 M(3,2)*(M(1,1)*M(2,3)-M(1,3)*M(2,1))+
     2 M(3,3)*(M(2,1)*M(1,2)-M(1,1)*M(2,2)))/P
      CTAU=((M(2,2)*M(1,3)-M(1,2)*M(2,3))*(M(3,2)*M(2,3)-M(2,2)*M(3,3))
     1 +(M(2,3)*M(1,1)-M(2,1)*M(1,3))*(M(3,3)*M(2,1)-M(3,1)*M(2,3))
     2 +(M(2,1)*M(1,2)-M(1,1)*M(2,2))*(M(3,1)*M(2,2)-M(2,1)*M(3,2)))/P
      TAU=ATAN2(STAU,CTAU)*RD
      WRITE(15,102)(ATM(I),I=1,4),TAU
102   FORMAT(1X,A4,'-',A4,'-',A4,'-',A4,F8.2)
      GO TO 5
4     WRITE(15,103)
103   FORMAT(' The torsion angle cannot be calculated')
5     RETURN
      END
 
      SUBROUTINE OH(RD,XO)
C-----This routine calculates the local orthogonal coordinates
C-----in a system of four atoms whose internal coordinates are
C-----given
      DIMENSION Q(0:2,0:100),A(0:8,0:100),XO(4,3)
      DO 1 K=1,7
1     A(K,0)=0.
      A(0,0)=1.
      A(4,0)=1.
      A(8,0)=1.
      Q(0,0)=0.
      Q(1,0)=0.
      Q(2,0)=0.
      DO 2 I=0,3
      READ(13,*)D,W,T
      N=I
      IF(I.EQ.0) GO TO 6
      IF(I.GT.2) GO TO 4
      IF(I.LT.2) GO TO 3
      T=0.0
      M=1
      GO TO 5
3     W=0.
      T=180.0-T
      M=0
      GO TO 5
4     M=N-1
5     N=M
      B0=-COS(W/RD)
      B1=-SIN(W/RD)
      B5=-SIN(T/RD)
      B8=COS(T/RD)
      B3=-B1*B8
      B4=B0*B8
      B6=B1*B5
      B7=-B0*B5
      B2=0.
      A(0,I)=A(0,M)*B0+A(1,M)*B3+A(2,M)*B6
      A(1,I)=A(0,M)*B1+A(1,M)*B4+A(2,M)*B7
      A(2,I)=A(0,M)*B2+A(1,M)*B5+A(2,M)*B8
      A(3,I)=A(3,M)*B0+A(4,M)*B3+A(5,M)*B6
      A(4,I)=A(3,M)*B1+A(4,M)*B4+A(5,M)*B7
      A(5,I)=A(3,M)*B2+A(4,M)*B5+A(5,M)*B8
      A(6,I)=A(6,M)*B0+A(7,M)*B3+A(8,M)*B6
      A(7,I)=A(6,M)*B1+A(7,M)*B4+A(8,M)*B7
      A(8,I)=A(6,M)*B2+A(7,M)*B5+A(8,M)*B8
      Q(0,I)=Q(0,M)+A(0,I)*D
      Q(1,I)=Q(1,M)+A(3,I)*D
      Q(2,I)=Q(2,M)+A(6,I)*D
6     DO 7 K=1,3
7     XO(I+1,K)=Q(K-1,I)
2     CONTINUE
      RETURN
      END
