C C CALCULATE BASE PLANE NORMALS, AND ROLL AND TILT ANGLES C C THIS PROGRAM USES THE OUTPUT COORDINATE LISTING FROM JOHN C ROSENBERG'S HELIX PROGRAM (AHLXRD.OUT) TO CALCULATE DIRECTION C COSINES AND CORRESPONDING ANGLES FOR THE NORMALS TO ALL BASE C PLANES AND TO THE BEST PLANE THROUGH BOTH BASES OF A PAIR. C IT THEN CALCULATES HELIX PARAMETERS THAT DEPEND ON BASE PLANE C NORMALS, OR ON THE ORIENTATION OF THE C6(PYRIMIDINE)-TO- C (C8)PURINE LINE THAT DEFINES THE LONG AXIS OF EACH BASE PAIR. C C ***NOTE: BUCKLE AND PROPELLER TWIST CORRECTIONS MADE 29.9.87. RED. C ***NOTE: NEW CAMBRIDGE NAMES AND SIGN CONVENTIONS ADDED 30.3.89 RED. C C PARAMETERS CALCULATED INCLUDE: C C TIP AND INCLINATION ANGLES FOR INDIVIDUAL BASES OR BASE PAIRS C ROLL AND TILT ANGLES BETWEEN ADJACENT BASES OR BASE PAIRS C (FOR ANALYTICAL EXPRESSIONS, SEE FRATINI, KOPKA, DREW AND C DICKERSON (1982) J.BIOL.CHEM. 257:14686-14707.) C PROPELLER TWIST BETWEEN BASES OF A PAIR, NEGATIVE FOR CLOCKWISE C ROTATION OF THE NEARER BASE IN A VIEW DOWN THE LONG AXIS C OF THE BASE PAIR. C BUCKLE = DIHEDRAL ANGLE BETWEEN BASES ALONG THEIR SHORT AXIS, C AFTER PROPELLER TWIST HAS BEEN ROTATED BACK TO ZERO. BUCKLE C IS POSITIVE IF THE BASE PAIR DOMES IN THE 5'-TO-3' DIRECTION C OF STRAND 1, CONVEX TOWARD THE 3' END. C C6/C8 = VECTOR LENGTH BETWEEN C6 AND C8 ATOMS. C INCL* = ANGLE OF INCLINATION OF C6-C8 VECTOR AWAY FROM C A PLANE PERPENDICULAR TO THE HELIX AXIS. C X DSP = PERPENDICULAR DISPLACEMENT OF THE C6-C8 VECTOR AWAY C FROM THE HELIX AXIS, VIEWED IN PROJECTION DOWN THE AXIS. C Y DSP = IN A PROJECTION ALONG THE HELIX AXIS, DROP A PERPENDI- C CULAR FROM THE HELIX CENTER TO THE C6/C8 LINE. Y DSP THEN C MEASURES THE DISTANCE OF THE MIDPOINT OF THE C6/C8 LINE C FROM THE INTERSECTION OF THE HELIX PERPENDICULAR. Y DSP C INDICATES THE EXTENT TO WHICH A SINGLE BASE PAIR IS PULLED C OUT OF THE HELIX STACK ALONG ITS LONG AXIS, VIEWED DOWN THE C HELIX AXIS. POSITIVE Y DSP IS FOR THE MOTION THAT PULLS THE C STRAND 1 BASE OUT OF THE STACK. C SLIDE = THE RELATIVE DISPLACEMENT OF MIDPOINTS OF THE C6-C8 LINES C OF TWO ADJACENT BASE PAIRS, VIEWED IN PROJECTION ON A PLANE C COMMON TO THE TWO BASE PAIRS. IT HENCE MEASURES LATERAL C DISPLACEMENT OF ADJACENT BASE PAIRS RELATIVE TO ONE ANOTHER, AND C IS TOTALLY INDEPENDENT OF THE CHOICE OF HELIX AXIS. C C (FOR FURTHER CLARIFICATION, AND AN ANALYTICAL EXPRESSION C FOR SLIDE, SEE THE APPENDIX TO "STRUCTURES OF BIOLOGICAL C MOLECULES AND ASSEMBLIES", EDITED BY A. MCPHERSON AND F. JURNAK, C WILEY, 1984.) C C INPUT FILE 7 (AHLXRD.OUT) IS THE ROSENBERG COORDINATE LIST OR A C SIMULATION THEREOF. THE HELIX AXIS IS ASSUMED TO BE ALONG Z. C TO OBTAIN CORRECT ROLL AND TILT SIGNS, STRAND 1 MUST ASCEND C THE Z AXIS TO MORE POSITIVE Z VALUES. IF NECESSARY, REVERSE THE C SIGNS OF Z, AND ALSO OF ONE OTHER AXIS TO MAINTAIN A RIGHT-HANDED C COORDINATE SYSTEM. C THE FILE BEGINS WITH THREE DISCRETIONARY TITLE CARDS. (ROSENBERG'S C PROGRAM EMITS A BLANK LINE, ONE TITLE LINE AND ANOTHER BLANK LINE.) C THEN FOLLOW ALL OF THE COORDINATE CARDS CONTAINING: X, Y, Z, ATOM C TYPE (C=1, N=2, O=3, P=4, OTHER=5), IN FORMAT(3F10.5,11X,A4). C C INPUT FILE 8 (BROLL.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. C6 OR C8 ON STRAND 1, C8 OR C6 ON STRAND 2 (2I5) C 3. NUMBERS OF ALL ATOMS IN BASE PLANE ON STRAND 1 (20I5) C 4. NUMBERS OF ALL ATOMS IN BASE PLANE ON STRAND 2 (20I5) C 5.....REPEAT CARDS 2-4, ONCE FOR EACH BASE PLANE IN THE HELIX. C C OUTPUT FILE 6 (BROLL.REC) CHRONICLES THE COURSE OF ITERATION OF C EIGENVECTORS, BASE PAIR BY BASE PAIR. ALTHOUGH INTERESTING DURING C PROGRAM DEVELOPMENT, IT HAS PROVEN OF LITTLE SUBSEQUENT INTEREST. C C OUTPUT FILE 9 (BROLL.OUT) CONTAINS THE HELIX PARAMETERS OF INTEREST: C FIRST A TABLE OF DIRECTION COSINES AND ANGLES, AND THEN A TABLE OF C HELIX DATA AS DESCRIBED ABOVE. "STRAND 3" MEANS BEST PLANE THROUGH C BOTH BASES OF EACH PAIR. C DIMENSION A(3,3),EIG(20,3),ATM(500,3),NAT(500),NTP(500),NA1(15), *NA2(15),XMIN(30,3),XMAJ(30,3),PLA(20,3,3),NC(20,2),TITL(19), *TITM(19),TITN(19),ANG(3),PHI(20,3,2),THE(20,3,2),TIL(20), *SV(3),RV(2,3),PV(2,3),PM(2),PTW(20),ALB(2),BKL(20), *SLI(20),DIS(20),HYP(20),ROT12(20),AM(3),BM(3),SLJ(20) C READ(8,1)NATM,NPLA 1 FORMAT(2I5) WRITE(6,7) WRITE(9,7) 7 FORMAT(' BASE ROLL AND TILT OUTPUT, CAMBRIDGE NOMENCLATURE') WRITE(6,8)NATM,NPLA WRITE(9,8)NATM,NPLA 8 FORMAT(I5,' ATOMS AND ',I5,' BASE PAIRS',I5) WRITE(6,200) WRITE(9,200) 200 FORMAT(1H ) READ(7,2)TITL WRITE(6,3)TITL WRITE(9,3)TITL 2 FORMAT(19A4) 3 FORMAT(19A4) READ(7,2)TITM WRITE(6,3)TITM WRITE(9,3)TITM READ(7,2)TITN WRITE(6,3)TITN WRITE(9,3)TITN WRITE(9,200) DO 115 I=1,20 DO 115 J=1,3 DO 115 K=1,2 115 THE(I,J,K)=0.0 C DO 4 I=1,NATM READ(7,5)(ATM(I,J),J=1,3),NTP(I),NAT(I) 5 FORMAT(3F10.5,11X,A4,5X,I5) 4 CONTINUE C DO 666 M=1,NPLA WRITE(6,200) WRITE(6,14)M 14 FORMAT(' BASE PAIR',I5) READ(8,30)NC(M,1),NC(M,2) 30 FORMAT(2I5) WRITE(6,31)NC(M,1),NC(M,2) 31 FORMAT(2I5) READ(8,9)(NA1(J),J=1,15) READ(8,9)(NA2(J),J=1,15) 9 FORMAT(15I5) WRITE(6,9)(NA1(J),J=1,15) WRITE(6,9)(NA2(J),J=1,15) C WRITE(6,15) 15 FORMAT(' STRAND 1') DO 10 K=1,15 NTEST=NA1(K) IF(NTEST)13,13,11 11 DO 12 J=1,3 XMIN(K,J)=ATM(NTEST,J) 12 CONTINUE 10 CONTINUE K=K+1 13 NPL=K-1 CALL MATRIX(XMIN,A,NPL) CALL EIGEN(A,EIG) DO 32 J=1,3 PLA(M,1,J)=EIG(20,J) 32 CONTINUE C WRITE(6,16) 16 FORMAT(' STRAND 2') DO 17 K=1,15 NTEST=NA2(K) IF(NTEST)20,20,18 18 DO 19 J=1,3 XMIN(K,J)=ATM(NTEST,J) 19 CONTINUE 17 CONTINUE K=K+1 20 NPL=K-1 CALL MATRIX(XMIN,A,NPL) CALL EIGEN(A,EIG) DO 33 J=1,3 PLA(M,2,J)=EIG(20,J) 33 CONTINUE C WRITE(6,21) 21 FORMAT(' BOTH STRANDS') DO 22 K=1,15 NTEST=NA1(K) IF(NTEST)23,23,24 24 DO 25 J=1,3 XMIN(K,J)=ATM(NTEST,J) 25 CONTINUE 22 CONTINUE K=K+1 23 K=K-1 DO 26 L=1,15 NTEST=NA2(L) IF(NTEST)27,27,28 28 DO 29 J=1,3 LI=K+L XMIN(LI,J)=ATM(NTEST,J) 29 CONTINUE 26 CONTINUE L=L+1 27 NPL=K+L-1 CALL MATRIX(XMIN,A,NPL) CALL EIGEN(A,EIG) DO 34 J=1,3 PLA(M,3,J)=EIG(20,J) 34 CONTINUE 666 CONTINUE 36 FORMAT(1H1) DO 37 J=1,3 WRITE(9,38)J 38 FORMAT(' STRAND',I4,' BASE NORMAL COSINES AND ANGLES') DO 39 M=1,NPLA DO 50 K=1,3 CSN=PLA(M,J,K) ANG(K)=57.29578*ACOS(CSN) 50 CONTINUE WRITE(9,40)(PLA(M,J,K),K=1,3),(ANG(K),K=1,3) 40 FORMAT(3F10.5,5X,3F10.2) 39 CONTINUE 37 CONTINUE C C CALCULATE AND STORE PHI ANGLES C DO 100 I=1,NPLA NC1=NC(I,1) NC2=NC(I,2) DELX=ATM(NC1,1)-ATM(NC2,1) DELY=ATM(NC1,2)-ATM(NC2,2) DENO=SQRT(DELX**2+DELY**2) SUMX=(ATM(NC1,1)+ATM(NC2,1))*0.5 SUMY=(ATM(NC1,2)+ATM(NC2,2))*0.5 HVEC=SQRT(SUMX**2+SUMY**2) SLI(I)=-(DELX*SUMX+DELY*SUMY)/DENO DIS(I)=(SUMX*DELY-DELX*SUMY)/DENO HYP(I)=HVEC ROT12(I)=DENO DELZ=ATM(NC1,3)-ATM(NC2,3) SV(1)=DELX SV(2)=DELY SV(3)=DELZ DO 130 K=1,3 RV(1,K)=PLA(I,1,K) RV(2,K)=PLA(I,2,K) 130 CONTINUE SV2=SV(1)**2+SV(2)**2+SV(3)**2 DO 131 J=1,2 RDS=RV(J,1)*SV(1)+RV(J,2)*SV(2)+RV(J,3)*SV(3) AKO=RDS/SV2 CALB=AKO*SQRT(SV2) ALB(J)=57.29578*ACOS(CALB) DO 132 K=1,3 PV(J,K)=RV(J,K)-AKO*SV(K) 132 CONTINUE PM(J)=SQRT(PV(J,1)**2+PV(J,2)**2+PV(J,3)**2) 131 CONTINUE BKL(I)=ALB(2)-ALB(1) COPT=PV(1,1)*PV(2,1)+PV(1,2)*PV(2,2)+PV(1,3)*PV(2,3) COPU=COPT/(PM(1)*PM(2)) ANSUP=57.29578*ACOS(COPU) PTW(I)=-ANSUP TIL(I)=57.29578*ATAN2(DELZ,DENO) DO 101 J=1,3 SINR=(PLA(I,J,1)*DELY-PLA(I,J,2)*DELX)/DENO PHI(I,J,1)=57.29578*ASIN(SINR) SINT=(PLA(I,J,1)*DELX+PLA(I,J,2)*DELY)/DENO PHI(I,J,2)=-57.29578*ASIN(SINT) 101 CONTINUE 100 CONTINUE C C CALCULATE AND STORE THETA ANGLES C NST=NPLA-1 DO 102 I=1,NST I1=I+1 NC1=NC(I,1) NC2=NC(I,2) ND1=NC(I1,1) ND2=NC(I1,2) DX=ATM(ND1,1)-ATM(ND2,1)-ATM(NC1,1)+ATM(NC2,1) DY=ATM(ND1,2)-ATM(ND2,2)-ATM(NC1,2)+ATM(NC2,2) DENO=SQRT(DX**2+DY**2) DO 103 J=1,3 DA=PLA(I1,J,1)-PLA(I,J,1) DB=PLA(I1,J,2)-PLA(I,J,2) SINR=(DA*DX+DB*DY)/DENO THE(I,J,1)=-57.29578*ASIN(SINR) SINT=(DB*DX-DA*DY)/DENO THE(I,J,2)=57.29578*ASIN(SINT) C C INSERTED COMMANDS BEGIN DO 250 K=1,3 AM(K)=(ATM(NC1,K)-ATM(NC2,K)+ATM(ND1,K)-ATM(ND2,K))/2. BM(K)=(ATM(NC1,K)+ATM(NC2,K)-ATM(ND1,K)-ATM(ND2,K))/2. 250 CONTINUE AMAG=SQRT(AM(1)**2+AM(2)**2+AM(3)**2) ADB=AM(1)*BM(1)+AM(2)*BM(2)+AM(3)*BM(3) SLJ(I)=-ADB/AMAG SLJ(I1)=0.0 C INSERTED COMMANDS END C 103 CONTINUE 102 CONTINUE C C PRINT OUT ROLL AND TILT ANGLES C WRITE(9,36) WRITE(9,3)TITL WRITE(9,3)TITM WRITE(9,3)TITN WRITE(9,200) DO 110 J=2,3 IF(J-2)150,150,151 150 WRITE(9,111) 111 FORMAT(' STRAND 1 ROLL AND TILT ANGLES ST *RAND 2 ROLL AND TILT ANGLES') WRITE(9,153) 153 FORMAT(' TIP INCL ROLL TILT TIP INCL * ROLL TILT') GO TO 154 151 WRITE(9,175) 175 FORMAT(1H ) WRITE(9,152) 152 FORMAT(' BEST PLANE THROUGH BOTH BASES') WRITE(9,112) 112 FORMAT(' TIP INCL ROLL TILT INCL* PR * TW BUCKLE SLIDE X DSP Y DSP C6/C8') 154 DO 113 M=1,NPLA IF(J-2)160,160,161 160 K=J-1 WRITE(9,174)PHI(M,K,1),PHI(M,K,2),THE(M,K,1),THE(M,K,2), *PHI(M,J,1),PHI(M,J,2),THE(M,J,1),THE(M,J,2) 174 FORMAT(4F7.2,7X,4F7.2) GO TO 162 161 WRITE(9,114)PHI(M,J,1),PHI(M,J,2),THE(M,J,1),THE(M,J,2),TIL(M),PTW *(M),BKL(M),SLJ(M),DIS(M),SLI(M),ROT12(M) 114 FORMAT(11F7.2) 162 CONTINUE 113 CONTINUE 110 CONTINUE WRITE(9,121) 121 FORMAT(1H ) WRITE(9,120) 120 FORMAT(' NOTE: ANGLES ARE CALCULATED FROM 5" END TO 3" END OF STR *AND 1, AND SIGNS OF') WRITE(9,122) 122 FORMAT(' ANGLES ALSO ARE CALCULATED WITH RESPECT TO STR *AND 1. TO EXAMINE INDIVIDUAL') WRITE(9,124) 124 FORMAT(' STRAND 2 BASES WITH RESPECT TO THEIR OWN STR *AND, REVERSE THE SIGNS OF ') WRITE(9,176) 176 FORMAT(' TIP AND TILT. FOR Z-DNA, REVERSE SIGNS OF INCLINATION, * X DSP AND Y DSP.') STOP END SUBROUTINE MATRIX(XMIN,A,NPL) C C BUILT THE A MATRIX FROM A LIST OF ATOM COORDINATES C DIMENSION XMIN(30,3),A(3,3) DIMENSION XMAJ(30,3),XAV(3) C DO 1 J=1,3 XAV(J)=0.0 1 CONTINUE DO 2 K=1,NPL DO 2 J=1,3 XAV(J)=XAV(J)+XMIN(K,J) 2 CONTINUE DO 3 J=1,3 XAV(J)=XAV(J)/NPL 3 CONTINUE DO 4 K=1,NPL DO 4 J=1,3 XMAJ(K,J)=XMIN(K,J)-XAV(J) 4 CONTINUE DO 5 I=1,3 DO 5 J=1,3 A(I,J)=0.0 5 CONTINUE DO 6 K=1,NPL DO 6 I=1,3 DO 6 J=1,3 A(I,J)=A(I,J)+XMIN(K,I)*XMAJ(K,J) 6 CONTINUE RETURN END SUBROUTINE EIGEN(A,EIG) DIMENSION A(3,3),EIG(20,3) DIMENSION AX(6,6),B(3,3) C C FIND THE LEAST EIGENVECTOR OF EQUATION AE=LE, WHERE A IS A 3X3 C MATRIX, L IS THE EIGENVALUE, AND E IS THE EIGENVECTOR. START C WITH E OR EIG=(0,0,1) AND ITERATE BY THE SWMB METHOD C (ACTA CRYST 14:185,1959) BY MULTIPLYING EIG BY B, WHERE C IS THE ADJOINT MATRIX OF A. HALT WHEN EIG CONVERGES (LESS THAN C 0.1% CHANGE BETWEEN CYCLES), OR AFTER 20 ITERATIONS. WRITE OUT C EACH EIGENVECTOR FOR INSPECTION. C DO 1 I=1,3 DO 1 J=1,3 I3=I+3 J3=J+3 AX(I,J)=A(I,J) AX(I3,J)=A(I,J) AX(I,J3)=A(I,J) AX(I3,J3)=A(I,J) 1 CONTINUE C DO 2 I=1,3 DO 2 J=1,3 I1=I+1 I2=I+2 J1=J+1 J2=J+2 B(J,I)=AX(I1,J1)*AX(I2,J2)-AX(I1,J2)*AX(I2,J1) 2 CONTINUE C DO 3 I=2,20 DO 3 J=1,3 EIG(I,J)=0.0 3 CONTINUE EIG(1,1)=0.0 EIG(1,2)=0.0 EIG(1,3)=1.0 WRITE(6,5)(EIG(1,J),J=1,3) 5 FORMAT(' EIGENVECTOR=',3F10.5) DO 4 K=1,19 K1=K+1 EX=B(1,1)*EIG(K,1)+B(1,2)*EIG(K,2)+B(1,3)*EIG(K,3) EY=B(2,1)*EIG(K,1)+B(2,2)*EIG(K,2)+B(2,3)*EIG(K,3) EZ=B(3,1)*EIG(K,1)+B(3,2)*EIG(K,2)+B(3,3)*EIG(K,3) EN=SQRT(EX**2+EY**2+EZ**2) EIG(K1,1)=EX/EN EIG(K1,2)=EY/EN EIG(K1,3)=EZ/EN WRITE(6,5)(EIG(K1,J),J=1,3) DO 6 J=1,3 DIF=ABS(EIG(K1,J)-EIG(K,J)) ARR=DIF/EIG(K1,J) ERR=ABS(ARR) IF(ERR-0.001)6,6,4 6 CONTINUE EIGZ=EIG(K1,3) SIGN=EIGZ/ABS(EIGZ) DO 8 J=1,3 EIG(20,J)=EIG(K1,J)*SIGN 8 CONTINUE GO TO 9 4 CONTINUE C C TEST THE EIGENVECTOR BY CARRYING OUT AE=LE AND EXAMINING L C 9 EX=A(1,1)*EIG(20,1)+A(1,2)*EIG(20,2)+A(1,3)*EIG(20,3) EY=A(2,1)*EIG(20,1)+A(2,2)*EIG(20,2)+A(2,3)*EIG(20,3) EZ=A(3,1)*EIG(20,1)+A(3,2)*EIG(20,2)+A(3,3)*EIG(20,3) ELAMX=EX/EIG(20,1) ELAMY=EY/EIG(20,2) ELAMZ=EZ/EIG(20,3) WRITE(6,10)ELAMX,ELAMY,ELAMZ 10 FORMAT(' LEAST EIGENVALUES= ',3F12.8) RETURN END