program kuknauhly character*80 fn parameter (nat0=2000,iu0=200) dimension x(nat0),y(nat0),z(nat0),ilt(nat0),am(iu0,12) c list=angle dtabase kind+ijkl dimension list(nat0,5) character*7 angles(12) data angles/' alpha',' beta',' gamma',' delta','epsilon' 1,' zeta', ' v0',' v1',' v2',' v3',' v4' 2,' chi'/ character*1 symbol(15) data symbol(1)/'H'/,symbol(6)/'C'/,symbol(7)/'N'/,symbol(8)/'O'/, 1symbol(15)/'P'/ if(iargc().eq.0)then write(6,*)' Usage: kuknauhly ' stop endif c write(6,*)' Filename:' c read(5,*)fn call getarg(1,fn) open(2,file=fn) read(2,*) read(2,*) nat write(6,*) nat,' atoms' if(nat.gt.nat0)then write(6,*)'Too many atoms' close(2) stop endif itot=0 do 111 i=1,nat 111 read(2,*)ilt(i),x(i),y(i),z(i) close(2) do 1 i1=1,nat xp=x(i1) yp=y(i1) zp=z(i1) if(ilt(i1).eq.15)then c suppose that this is backbone phosphor, climbing to c4' do 2 i2=1,nat if(ilt(i2).eq.8)then xo=x(i2) yo=y(i2) zo=z(i2) d12=(xo-xp)**2+(yo-yp)**2+(zo-zp)**2 if(d12.gt.2.37.and.d12.lt.3.2)then do 3 i3=1,nat if(ilt(i3).eq.6)then xc=x(i3) yc=y(i3) zc=z(i3) d23=(xo-xc)**2+(yo-yc)**2+(zo-zc)**2 if(d23.gt.1.8.and.d23.lt.2.2)then nh=0 do i4=1,nat if(ilt(i4).eq.1)then dch=(x(i4)-xc)**2+(y(i4)-yc)**2+(z(i4)-zc)**2 if(dch.gt.0.8.and.dch.lt.1.2)nh=nh+1 endif enddo c Main Branch c4->c3: if(nh.eq.2)then c constructing angle alpha: do i4=1,nat if(ilt(i4).eq.8.and.i4.ne.i2)then dpo=(x(i4)-xp)**2+(y(i4)-yp)**2+(z(i4)-zp)**2 if(dpo.gt.2.37.and.dpo.lt.3.2)then call ada(i4,i1,i2,i3,1,itot,list) endif endif enddo c P-O-C(H2)-found, looking for C4' do 4 i4=1,nat x4=x(i4) y4=y(i4) z4=z(i4) if(ilt(i4).eq.6)then d34=(x4-xc)**2+(y4-yc)**2+(z4-zc)**2 if(d34.gt.2.04.and.d34.lt.2.54)then c P-O-C(H2)-C4' call ada(i1,i2,i3,i4,2,itot,list) do 5 i5=1,nat x5=x(i5) y5=y(i5) z5=z(i5) if(ilt(i5).eq.6.and.i5.ne.i3)then c Main Branch d45=(x4-x5)**2+(y4-y5)**2+(z4-z5)**2 if(d45.gt.2.14.and.d45.lt.2.54)then c P-O-C(H2)-C4'-C3' call ada(i2,i3,i4,i5,3,itot,list) do 6 i6=1,nat x6=x(i6) y6=y(i6) z6=z(i6) if(ilt(i6).eq.8)then d56=(x6-x5)**2+(y6-y5)**2+(z6-z5)**2 if(d56.gt.1.81.and.d56.lt.2.21)then c P-O-C(H2)-C4'-C3'-O call ada(i3,i4,i5,i6,4,itot,list) do 7 i7=1,nat x7=x(i7) y7=y(i7) z7=z(i7) if(ilt(i7).eq.15)then d67=(x6-x7)**2+(y6-y7)**2+(z6-z7)**2 if(d67.gt.2.37.and.d67.lt.2.76)then c P-O-C(H2)-C4'-C3'-O-P call ada(i4,i5,i6,i7,5,itot,list) do 8 i8=1,nat x8=x(i8) y8=y(i8) z8=z(i8) if(ilt(i8).eq.8.and.i8.ne.i6)then d78=(x8-x7)**2+(y8-y7)**2+(z8-z7)**2 if(d78.gt.2.37.and.d78.lt.2.76)then c P-O-C(H2)-C4'-C3'-O-P-O call ada(i5,i6,i7,i8,6,itot,list) do 9 i9=1,nat x9=x(i9) y9=y(i9) z9=z(i9) if(ilt(i9).eq.6)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.1.84.and.d89.lt.2.26)then c P-O-C(H2)-C4'-C3'-O-P-O-C call ada(i6,i7,i8,i9,1,itot,list) endif endif 9 continue endif endif 8 continue endif endif 7 continue endif endif 6 continue endif endif c end of main branch if(ilt(i5).eq.8)then c Sugar Branch d45=(x4-x5)**2+(y4-y5)**2+(z4-z5)**2 if(d45.gt.1.81.and.d45.lt.2.21)then c P-O-C(H2)-C4'-O do 61 i6=1,nat x6=x(i6) y6=y(i6) z6=z(i6) if(ilt(i6).eq.6.and.i4.ne.i6)then d56=(x6-x5)**2+(y6-y5)**2+(z6-z5)**2 if(d56.gt.1.81.and.d56.lt.2.21)then c P-O-C(H2)-C4'-O-C do 71 i7=1,nat x7=x(i7) y7=y(i7) z7=z(i7) if(ilt(i7).eq.7)then d67=(x6-x7)**2+(y6-y7)**2+(z6-z7)**2 if(d67.gt.1.99.and.d67.lt.2.39)then c P-O-C(H2)-C4'-O-C-N do 81 i8=1,nat x8=x(i8) y8=y(i8) z8=z(i8) if(ilt(i8).eq.6)then d78=(x8-x7)**2+(y8-y7)**2+(z8-z7)**2 if(d78.gt.1.70.and.d78.lt.2.10)then c P-O-C(H2)-C4'-O-C-N-C c is it c2 or c4 (with O or N) io=0 in=0 ih=0 do 91 i9=1,nat x9=x(i9) y9=y(i9) z9=z(i9) if(ilt(i9).eq.8)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.1.29.and.d89.lt.1.69)io=io+1 endif if(ilt(i9).eq.7)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.1.65.and.d89.lt.2.05)in=in+1 endif if(ilt(i9).eq.1)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.0.90.and.d89.lt.1.37)ih=ih+1 endif 91 continue c P-O-C(H2)-C4'-O-C-N-C2(O) c P-O-C(H2)-C4'-O-C-N-C4(N) if((io.eq.1.and.in.eq.2).or.(io.eq.0.and. 1 ih.eq.0.and. 1 in.eq.2))call ada(i5,i6,i7,i8,12,itot,list) endif endif 81 continue endif endif 71 continue endif endif 61 continue endif endif c end of sugar branch 5 continue endif endif 4 continue endif c end of main branch c4->c3 c Main Branch c3->c4: if(nh.eq.1)then c P-O-C3 c constructing angle zeta do i4=1,nat if(ilt(i4).eq.8.and.i4.ne.i2)then dpo=(x(i4)-xp)**2+(y(i4)-yp)**2+(z(i4)-zp)**2 if(dpo.gt.2.37.and.dpo.lt.3.2) 1 call ada(i4,i1,i2,i3,6,itot,list) endif enddo do 42 i4=1,nat x4=x(i4) y4=y(i4) z4=z(i4) if(ilt(i4).eq.6)then inh=0 do i5=1,nat if(ilt(i5).eq.1)then dch=(x(i5)-x4)**2+(y(i5)-y4)**2+(z(i5)-z4)**2 if(dch.gt.0.8.and.dch.lt.1.2)inh=inh+1 endif enddo d34=(x4-xc)**2+(y4-yc)**2+(z4-zc)**2 if(d34.gt.2.14.and.d34.lt.2.54.and.inh.eq.1)then c P-O-C3-C4' call ada(i1,i2,i3,i4,5,itot,list) do 52 i5=1,nat x5=x(i5) y5=y(i5) z5=z(i5) c Main branch if(ilt(i5).eq.6.and.i5.ne.i3)then d45=(x4-x5)**2+(y4-y5)**2+(z4-z5)**2 if(d45.gt.2.14.and.d45.lt.2.54)then c P-O-C3-C4-CH2 call ada(i2,i3,i4,i5,4,itot,list) do 62 i6=1,nat x6=x(i6) y6=y(i6) z6=z(i6) if(ilt(i6).eq.8)then d56=(x6-x5)**2+(y6-y5)**2+(z6-z5)**2 if(d56.gt.1.81.and.d56.lt.2.21)then c P-O-C3-C4-CH2-O call ada(i3,i4,i5,i6,3,itot,list) do 72 i7=1,nat x7=x(i7) y7=y(i7) z7=z(i7) if(ilt(i7).eq.15)then d67=(x6-x7)**2+(y6-y7)**2+(z6-z7)**2 if(d67.gt.2.37.and.d67.lt.2.76)then c P-O-_c3_C4-CH2-O-P call ada(i4,i5,i6,i7,2,itot,list) do 82 i8=1,nat x8=x(i8) y8=y(i8) z8=z(i8) if(ilt(i8).eq.8.and.i8.ne.i6)then d78=(x8-x7)**2+(y8-y7)**2+(z8-z7)**2 if(d78.gt.2.37.and.d78.lt.2.76)then c P-O-C3-C4-CH2-O-P-O call ada(i5,i6,i7,i8,1,itot,list) do 92 i9=1,nat x9=x(i9) y9=y(i9) z9=z(i9) if(ilt(i9).eq.6)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.1.84.and.d89.lt.2.26)then c P-O-C(H2)-C4'-C3'-O-P-O-C call ada(i6,i7,i8,i9,6,itot,list) endif endif 92 continue endif endif 82 continue endif endif 72 continue endif endif 62 continue endif endif if(ilt(i5).eq.8)then c Sugar Branch d45=(x4-x5)**2+(y4-y5)**2+(z4-z5)**2 if(d45.gt.1.81.and.d45.lt.2.21)then c P-O-C3-C4'-O do 63 i6=1,nat x6=x(i6) y6=y(i6) z6=z(i6) if(ilt(i6).eq.6.and.i4.ne.i6)then d56=(x6-x5)**2+(y6-y5)**2+(z6-z5)**2 if(d56.gt.1.81.and.d56.lt.2.21)then c P-O-C3-C4'-O-C do 73 i7=1,nat x7=x(i7) y7=y(i7) z7=z(i7) if(ilt(i7).eq.7)then d67=(x6-x7)**2+(y6-y7)**2+(z6-z7)**2 if(d67.gt.1.99.and.d67.lt.2.39)then c P-O-C3-C4'-O-C-N do 83 i8=1,nat x8=x(i8) y8=y(i8) z8=z(i8) if(ilt(i8).eq.6)then d78=(x8-x7)**2+(y8-y7)**2+(z8-z7)**2 if(d78.gt.1.70.and.d78.lt.2.10)then c P-O-C(H2)-C4'-O-C-N-C c is it c2 or c4 (with O or N) io=0 in=0 ih=0 do 93 i9=1,nat x9=x(i9) y9=y(i9) z9=z(i9) if(ilt(i9).eq.8)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.1.29.and.d89.lt.1.69)io=io+1 endif if(ilt(i9).eq.7)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.1.65.and.d89.lt.2.05)in=in+1 endif if(ilt(i9).eq.1)then d89=(x8-x9)**2+(y8-y9)**2+(z8-z9)**2 if(d89.gt.0.90.and.d89.lt.1.37)ih=ih+1 endif 93 continue c P-O-C(H2)-C4'-O-C-N-C2(O) c P-O-C(H2)-C4'-O-C-N-C4(N) if((io.eq.1.and.in.eq.2).or.(io.eq.0.and. 1 ih.eq.0.and. 1 in.eq.2))call ada(i5,i6,i7,i8,12,itot,list) endif endif 83 continue endif endif 73 continue endif endif 63 continue endif endif c end of sugar branch 52 continue endif endif 42 continue endif c end of main branch c4->c3 endif endif 3 continue endif endif 2 continue endif 1 continue iumax=0 do i=1,iu0 do j=1,12 am(i,j)=88888.8 enddo enddo do 991 ib=1,12 iu=0 do 99 i=1,itot i1=list(i,1) i2=list(i,2) i3=list(i,3) i4=list(i,4) m=list(i,5) if(m.eq.ib)then iu=iu+1 if(iu.gt.iumax)iumax=iu if(iu.gt.iu0)then write(6,*)'am matrix dimension overflow' stop endif call DBT(x(i1),y(i1),z(i1),x(i2),y(i2),z(i2),x(i3), 1 y(i3),z(i3),x(i4),y(i4),z(i4),AD,BANGLE,ANGLE) write(6,6001)iu,m,angles(m),i1,symbol(ilt(i1)),i2, 1 symbol(ilt(i2)),i3,symbol(ilt(i3)),i4,symbol(ilt(i4)),ANGLE 6001 format(i3,i3,1x,a7,4(i5,a1),f12.2) am(iu,ib)=ANGLE endif 99 continue if(iu.gt.0)write(6,*) 991 continue write(6,6003)(angles(j),j=1,6),angles(12) 6003 format(/,' res ',12(2x,a7)) do 992 i=1,iumax 992 write(6,6002)i,(am(i,j),j=1,6),am(i,12) 6002 format(i3,2x,12f9.1) end subroutine DBT(x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3, 1AD,BANGLE,ANGLE) dimension A(3),P(3),E0(3),XX(3),YY(3) PI=real(3.141592653589D0) A(1)=X2-X3 A(2)=Y2-Y3 A(3)=Z2-Z3 P(1)=X0-X1 P(2)=Y0-Y1 P(3)=Z0-Z1 E0(1)=X1-X2 E0(2)=Y1-Y2 E0(3)=Z1-Z2 SSE0E0=E0(1)*E0(1)+E0(2)*E0(2)+E0(3)*E0(3) DO 8 KK=1,3 8 E0(KK)=E0(KK)/SQRT(SSE0E0) SSPP=P(1)*P(1)+P(2)*P(2)+P(3)*P(3) AD=SQRT(SSPP) SSAE0=A(1)*E0(1)+A(2)*E0(2)+A(3)*E0(3) SSPE0=P(1)*E0(1)+P(2)*E0(2)+P(3)*E0(3) DO 9 KK=1,3 XX(KK)=A(KK)-E0(KK)*SSAE0 9 YY(KK)=P(KK)-E0(KK)*SSPE0 SSXXXX=XX(1)*XX(1)+XX(2)*XX(2)+XX(3)*XX(3) SSYYYY=YY(1)*YY(1)+YY(2)*YY(2)+YY(3)*YY(3) SSXXYY=XX(1)*YY(1)+XX(2)*YY(2)+XX(3)*YY(3) APX=XX(2)*YY(3)-XX(3)*YY(2) APY=XX(3)*YY(1)-XX(1)*YY(3) APZ=XX(1)*YY(2)-XX(2)*YY(1) SSE0AP=APX*E0(1)+APY*E0(2)+APZ*E0(3) SIGN=1.00 IF (SSE0AP.NE.0.00)SIGN=-SSE0AP/ABS(SSE0AP) CB=0.00 IF (AD.NE.0.00)CB=-SSPE0/AD CA=0.00 IF (SSXXXX.NE.0.00.AND.SSYYYY.NE.0.00) 1CA=-SSXXYY/SQRT(SSXXXX)/SQRT(SSYYYY) BANGLE=0.00 ANGLE=0.00 IF (CB.EQ.-1.00)BANGLE=PI IF (ABS(CB).LT.1.00)BANGLE=ACOS(CB) IF (CA.LE.-1.00) ANGLE=PI IF (ABS(CA).LT.1.00) ANGLE=ACOS(CA) ANGLE=180.00* ANGLE/PI BANGLE=180.00*BANGLE/PI ANGLE=SIGN*ANGLE if(ANGLE.lt.-180.0)ANGLE=360.0+ANGLE return end subroutine ada(i,j,k,l,m,itot,list) parameter (nat0=2000) dimension list(nat0,5) do 1 ii=1,itot 1 if( list(ii,1).eq.i.and.list(ii,2).eq.j.and.list(ii,3).eq.k. 1and.list(ii,4).eq.l.and.list(ii,5).eq.m)return do 2 ii=1,itot 2 if( list(ii,1).eq.l.and.list(ii,2).eq.k.and.list(ii,3).eq.j. 1and.list(ii,4).eq.i.and.list(ii,5).eq.m)return itot=itot+1 list(itot,1)=i list(itot,2)=j list(itot,3)=k list(itot,4)=l list(itot,5)=m return end