C  CALCULATE TORSION ANGLES AND GLYCOSYL ANGLE
C    THIS PROGRAM USES THE OUTPUT COORDINATE LISTING FROM JOHN
C    ROSENBERG'S HELIX PROGRAM (AHLXRD.OUT) TO CALCULATE TORSION
C    ANGLES.  ANGLES ARE CURRENT IUB/IUPAC RECOMMENDATIONS:
C
C P-(ALPHA)-O5'-(BETA)-C5'-(GAMMA)-C4'-(DELTA)-C3'-(EPSILON)-O3'-(ZETA)-P
C
C    PYRIMIDINE:  C2-N1-(CHI)-C1'-O4'
C    PURINE:      C4-N9-(CHI)-C1'-04'
C
C    ANGLES ARE LISTED IN 5'-TO-3' DIRECTION ALONG EACH STRAND.  IN
C       STRAND 1 LISTING, DEL = DIFFERENCE IN DELTA ANGLES AT TWO 
C       ENDS OF BASE PAIR.  IN STRAND 2 LISTING, MEAN = MEAN OF DELTA
C       ANGLES AT TWO ENDS OF BASE PAIR.  (WARNING: NOTE REVERSAL OF
C       DIRECTION DOWN HELIX FOR STRANDS 1 AND 2.  THE VERY FIRST DEL
C       ENTRY CORRESPONDS TO THE VERY LAST MEAN ENTRY, ETC.)
C
C    INPUT FILE 7 (AHLXRD.OUT) IS THE ROSENBERG COORDINATE LIST AS IN
C       THE BROLL PROGRAM.
C
C    INPUT FILE 8 (DTORAN.INP) CONTAINS SEQUENTIAL NUMBERS OF CERTAIN ATOMS
C       IN THE FILE 7 LIST, AS FOLLOWS:
C
C       1. TOTAL NUMBER OF ATOMS, NUMBER OF BASE PAIRS IN HELIX (2I5)
C       2. NUMBERS OF ATOMS: O5', C5', C4', C3', O3', P IN FIRST NUCLEOTIDE
C           ALONG STRAND 1 (6I5).  
C       3..REPEAT LINE FOR EACH NUCLEOTIDE ALONG STRAND 1.
C       4. NUMBERS OF ATOMS: C2, N1, C1', O4' FOR PYRIMIDINES, OR
C           C4, N9, C1', 04' FOR PURINES, IN FIRST NUCLEOTIDE ALONG
C           STRAND 1 (4I5).
C       5..REPEAT LINE FOR EACH NUCLEOTIDE ALONG STRAND 1.
C    
C    THE PROGRAM AUTOMATICALLY GENERATES NUMBERS OF EQUIVALENT ATOMS IN
C       STRAND 2 FOR THE SAKE OF ECONOMY OF INPUT.  BUT THIS MEANS THAT
C       THE ORDERING OF ATOMS MUST BE THE SAME IN BOTH STRANDS (USUALLY
C       THE CASE).  
C      
C    ***NOTE: ANGLES HAVE BEEN CHANGED (29.3.89) TO LIE BETWEEN 0 AND 360
C       DEGREES (INSTEAD OF THE FORMER -180 TO +180).  THIS MAKES IT MUCH
C       EASIER TO AVERAGE OR OTHERWISE DEAL WITH TRANS ANGLES WHICH VARY 
C       TO EITHER SIDE OF 180.
C
C
      DIMENSION TITL(19),ATM(500,3),NATO(2,250),NCHI(2,25,4),NA(4),
     *ANG(300),ACHI(25),ANSW(50,8)
      WRITE(9,50)
50    FORMAT(' MAIN CHAIN AND GLYCOSYL TORSION ANGLE CAL
     *CULATION')
      READ(7,1)TITL
      WRITE(9,1)TITL
1     FORMAT(19A4)
      READ(7,1)TITL
      WRITE(9,1)TITL
      READ(7,1)TITL
      WRITE(9,1)TITL
4     FORMAT(2I5)
5     FORMAT(I5,' ATOMS AND',I5,' BASE PAIRS')
7     FORMAT(6I5)
9     FORMAT(4I5)
20    FORMAT(8F8.1)
21    FORMAT(' STRAND',I5)
22    FORMAT('   ALPHA    BETA   GAMMA   DELTA  EPSILON   
     1ZETA     CHI     DEL')
32    FORMAT('   ALPHA    BETA   GAMMA   DELTA  EPSILON   
     1ZETA     CHI    MEAN')
25    FORMAT(3F10.5)
666   FORMAT(1H )
C
C   READ STARTING PARAMETERS
C
      READ(8,4)NAT,NBP
      WRITE(9,5)NAT,NBP
C
C   READ MAIN CHAIN ATOM NUMBERS
C
      DO 6 I=1,NBP
      NLO=I*6-5
      NUP=I*6
6     READ(8,7)(NATO(1,J),J=NLO,NUP)
      NHA=NAT/2
      MAX=6*NBP-1
      DO 10 I=1,MAX
10    NATO(2,I)=NATO(1,I)+NHA
C
C   READ CHI ATOM NUMBERS
C
      DO 8 I=1,NBP
      READ(8,9)(NCHI(1,I,J),J=1,4)
      DO 12 J=1,4
12    NCHI(2,I,J)=NCHI(1,I,J)+NHA
8     CONTINUE
C
C   READ ATOM COORDINATES
C
      DO 23 I=1,NAT
23    READ(7,25)(ATM(I,J),J=1,3)
C
      DO 30 IC=1,50
      DO 30 JC=1,8
30    ANSW(IC,JC)=0.
C
      DO 13 M=1,2
C
C   CALCULATE MAIN CHAIN TORSION ANGLES
C
      ANG(1)=0.0
      ANG(2)=0.0
      NMA=6*NBP-4
      DO 14 I=1,NMA
      DO 24 J=1,4
      IU=I+J-1
24    NA(J)=NATO(M,IU)
      CALL TORSION(NA,TE,ATM)
      I2=I+2
14    ANG(I2)=TE
      NMB=NMA+3
      NMC=NMA+4
      ANG(NMB)=0.0
      ANG(NMC)=0.0
C
C   CALCULATE GLYCOSYL CHI ANGLES
C
      DO 15 I=1,NBP
      DO 16 J=1,4
16    NA(J)=NCHI(M,I,J)
      CALL TORSION(NA,TE,ATM)
15    ACHI(I)=TE
C
C   REARRANGE ANGLES AND PRINT OUT
C
      DO 17 I=1,NBP
      DO 18 J=1,6
      NO=6*I-6+J
      IM=NBP*(M-1)+I
18    ANSW(IM,J)=ANG(NO)
      ANSW(IM,7)=ACHI(I)
17    CONTINUE
13    CONTINUE
      M=1
      WRITE(9,666)
      WRITE(9,21)M
      WRITE(9,22)
      DO 31 I=1,NBP
      IM=NBP*2-I+1
      DEL=ANSW(I,4)-ANSW(IM,4)
      ANSW(I,8)=DEL
      WRITE(9,20)(ANSW(I,J),J=1,8)
31    CONTINUE
C
      M=2
      WRITE(9,666)
      WRITE(9,21)M
      WRITE(9,32)
      DO 33 I=1,NBP
      IM=I+NBP
      IN=NBP-I+1
      FMEAN=0.5*(ANSW(IM,4)+ANSW(IN,4))
      ANSW(IM,8)=FMEAN
      WRITE(9,20)(ANSW(IM,J),J=1,8)
33    CONTINUE
      STOP
      END
C
C   TORSION ANGLE CALCULATING SUBROUTINE
C
      SUBROUTINE TORSION(NA,TE,ATM)
      DIMENSION ATM(500,3),NA(4)
      DIMENSION AV(3),BV(3),CV(3),S(3),T(3)
      DO 1 J=1,3
      J1=J+1
      NS=NA(J)
      NE=NA(J1)
      AV(J)=ATM(NE,1)-ATM(NS,1)
      BV(J)=ATM(NE,2)-ATM(NS,2)
1     CV(J)=ATM(NE,3)-ATM(NS,3)
      S(1)=BV(1)*CV(2)-BV(2)*CV(1)
      S(2)=CV(1)*AV(2)-CV(2)*AV(1)
      S(3)=AV(1)*BV(2)-AV(2)*BV(1)
      T(1)=BV(2)*CV(3)-BV(3)*CV(2)
      T(2)=CV(2)*AV(3)-CV(3)*AV(2)
      T(3)=AV(2)*BV(3)-AV(3)*BV(2)
C
      RM=SQRT(AV(2)**2+BV(2)**2+CV(2)**2)
      SM=SQRT(S(1)**2+S(2)**2+S(3)**2)
      TM=SQRT(T(1)**2+T(2)**2+T(3)**2)
      TOP=AV(1)*(BV(2)*CV(3)-BV(3)*CV(2))+
     *AV(2)*(BV(3)*CV(1)-BV(1)*CV(3))+
     *AV(3)*(BV(1)*CV(2)-BV(2)*CV(1))
      SAL=(RM*TOP)/(SM*TM)
      CAL=(S(1)*T(1)+S(2)*T(2)+S(3)*T(3))/(SM*TM)
      TE=57.29578*ATAN2(SAL,CAL)
      IF(TE)667,668,668
667   TE=TE+360.
668   CONTINUE
      RETURN
      END