PROGRAM secpro IMPLICIT none character*80 fn,s80 real*8 p(100),av(3) integer*4 i,nat,j,iargc,ig,nnn(100),iagg,id,ii real*8,allocatable::r(:) integer*4,allocatable::iz(:),bt(:,:),nb(:),us(:), 1itpp(:,:) c WRITE(6,3000) 3000 FORMAT(' Analyzes fi psi omega in a protein (FILE.X)') fn='FILE.X' if(iargc().gt.0)call getarg(1,fn) open(33,file='CONF.PRN') do 3 i=1,100 p(i)=0.0d0 3 nnn(i)=0 av(1)=0.0d0 av(2)=0.0d0 av(3)=0.0d0 ig=0 open(2,file=fn) 1000 read(2,*,end=99,err=99) read(2,*,end=99,err=99)nat ig=ig+1 if(ig.eq.1)allocate(r(3*nat),iz(nat),nb(nat),bt(nat,7),us(nat), 1itpp(nat,8)) do 1 i=1,nat read(2,2000)s80 2000 format(a80) id=0 do 5 ii=1,len(s80) 5 if(s80(ii:ii).eq.'.')id=id+1 if(id.eq.3)then c only coordinates, make bond table later: read(s80,*)iz(i),(r(j+3*(i-1)),j=1,3) else if(id.eq.4)then c bond table present, read it read(s80,*)iz(i),(r(j+3*(i-1)),j=1,3),(bt(i,ii),ii=1,7) else call report('Unkown format') endif endif 1 continue if(ig.eq.1)then write(6,*)nat,' atoms' if(id.eq.3)call mkbonds(nat,r,bt,nb,iz) endif do 2 i=1,nat nb(i)=7 us(i)=1 do 21 ii=1,7 if(bt(i,ii).eq.0)then nb(i)=ii-1 goto 2 endif 21 continue 2 continue c determine amide groups: if(ig.eq.1)call getamides(nat,iz,us,nb,bt,r,iagg,itpp) c print the angles if(ig.eq.1)write(6,605) 605 format(' geometry fi psi w a-H R-aH beta', 1' 3_10 PPI PPII other') call anala(ig,nat,itpp,iagg,r,nnn,av) goto 1000 99 close(2) if(nnn(2).ne.0)then do 4 i=1,7 4 p(i)=dble(nnn(10+i))/dble(nnn(2))*100.0d0 endif if(ig.ne.0)then do 6 i=1,3 av(i)=av(i)/dble(ig) 6 call cper(av(i)) endif write(6,609)ig,(nnn(i),i=1,4), 1nnn(1+10),p(1), 1nnn(2+10),p(2), 1nnn(3+10),p(3), 1nnn(4+10),p(4), 1nnn(5+10),p(5), 1nnn(6+10),p(6), 1nnn(7+10),p(7), 1(av(i),i=1,3) 609 format(/,i9,' geometries',/, 1 i9,' analyzed amide groups',/, 1 i9,' analyzed fi psi omega angless',/, 1 i9,' cis amides',/, 1 i9,' trans amides',/,/, 1 i9,' alpha-helices (',f10.2,' %)',/, 1 i9,' left alpha-helices (',f10.2,' %)',/, 1 i9,' betas (',f10.2,' %)',/, 1 i9,' 3_10s (',f10.2,' %)',/, 1 i9,' PPI (',f10.2,' %)',/, 1 i9,' PPII (',f10.2,' %)',/, 1 i9,' others (',f10.2,' %)',/,/, 1 9x,' Average fi ',f10.2,/, 1 9x,' Average psi ',f10.2,/, 1 9x,' Average w ',f10.2,/) close(33) end subroutine report(s) character*(*) s write(6,*)s stop end subroutine mkbonds(nat,r,nbt,nb,iz) implicit none integer*4 nat,nbt(nat,7),nb(nat),iz(nat),i,j,ti,tj real*8 r(*),x1,y1,z1,x2,y2,z2,rt,rb c Du H He Li Be B C N O F c Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti c V Cr Mn Fe Co Ni Cu Zn Ga Ge As real bonding(88) data bonding/ 10.50,0.32,0.98,1.28,0.95,0.87,0.82,0.80,0.78,0.77, 10.76,1.59,1.41,1.23,1.16,1.11,1.20,1.04,1.03,2.08,1.79,1.49,1.37, 11.27,1.23,1.22,1.22,1.21,1.20,1.22,1.30,1.31,1.27,1.25,1.21,1.19, 11.17,2.21,1.96,1.67,1.50,1.39,1.35,1.32,1.30,1.30,1.33,1.39,1.53, 11.49,1.46,1.45,1.41,1.38,1.36,2.40,2.03,1.74,1.70,1.70,1.69,1.68, 11.67,1.90,1.66,1.64,1.64,1.63,1.62,1.61,1.75,1.61,1.49,1.39,1.35, 11.33,1.31,1.32,1.35,1.39,1.54,1.53,1.52,1.51,1.51,1.50,1.49,0.10/ do 1 i=1,nat nb(i)=0 do 1 j=1,7 1 nbt(i,j)=0 c do 2 i=1,nat ti=iz(i)+1 if(ti.lt.1.or.ti.gt.89)call report('type out of range') rb=dble(bonding(ti)) x1=r((i-1)*3+1) y1=r((i-1)*3+2) z1=r((i-1)*3+3) do 2 j=i+1,nat tj=iz(j)+1 if(ti.ne.2.or.tj.ne.2)then if(tj.lt.1.or.tj.gt.88)call report('type out of range') rt=(rb+dble(bonding(tj)))**2 x2=(x1-r((j-1)*3+1))**2 if(x2.lt.rt)then y2=(y1-r((j-1)*3+2))**2+x2 if(y2.lt.rt)then z2=(z1-r((j-1)*3+3))**2+y2 if(z2.lt.rt)then nb(i)=nb(i)+1 nb(j)=nb(j)+1 nbt(i,nb(i))=j nbt(j,nb(j))=i endif endif endif endif 2 continue write(6,*)'Bonds made' return end subroutine getamides(nat,ity,us,nt,bt,r,iagg,itpp) integer*4 ity(*),nat,ia,ia2,ia3,ia4,ia5,ia6,ia7,ia8,ia9, 1ix2,ix3,ix4,ix5,ix6,iag,iac,iah,iap,iax, 1us(*),nt(*),bt(nat,7),itpp(nat,8) logical lcc,lcon real*8 xi,yi,zi,r(*),xj,yj,zj,d2 integer*4,allocatable:: 1ipex(:,:),ippp(:,:),ipeh(:,:),ipep(:,:),ipec(:,:) c identify amide groups c number of ordinary trans amide groups: (ipep) iag=0 c number of ordinary cis amide groups: (ipec) iac=0 c number of terminal amide groups: (ipeh) iah=0 c number of proline amide groups: (ippp) iap=0 c number of strange amide groups: (ipex) iax=0 c c iagg - number of all recognized amide groups: (itpp) iagg=0 c ipt: c number of ordinary amide groups: iag 1 c number of terminal amide groups: iah 2 c number of proline amide groups: iap 3 c number of strange amide groups: iax 4 c number of cis amide groups: iac allocate(ipex(nat,8),ippp(nat,8),ipeh(nat,8),ipep(nat,8), 1ipec(nat,8)) do 9 ia=1,nat if(ity(ia).eq.8.and.us(ia).eq.1.and.nt(ia).eq.1)then do 124 ix2=1,nt(ia) ia2=bt(ia,ix2) if(ity(ia2).eq.6.and.us(ia2).eq.1)then do 125 ix3=1,nt(ia2) ia3=bt(ia2,ix3) if(ity(ia3).eq.7.and.us(ia3).eq.1)then c ia ia2 ia3 c O = C - N c c carbon on C: lcc=.false. do 135 ix6=1,nt(ia2) ia6=bt(ia2,ix6) if(ity(ia6).eq.6)then lcc=.true. goto 1003 endif 135 continue 1003 continue c do 126 ix4=1,nt(ia3) ia4=bt(ia3,ix4) if(ity(ia4).eq.1)then c ia ia2 ia3 ia4 c O = C - N - H xi=r(3*(ia -1)+1) yi=r(3*(ia -1)+2) zi=r(3*(ia -1)+3) xj=r(3*(ia4-1)+1) yj=r(3*(ia4-1)+2) zj=r(3*(ia4-1)+3) d2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 if(d2.gt.7.0d0)then c trans do 127 ix5=1,nt(ia3) ia5=bt(ia3,ix5) if(ity(ia5).eq.6.and.ia5.ne.ia2)then c O = C - N(-C) - H if(lcc)then c O = C(-C) - N(-C) - H call swallow(nat,us,iag,ipep, 1 ia,ia2,ia3,ia4,ia5,ia6,ia6,ia6) goto 9 endif endif if(ity(ia5).eq.1.and.ia5.ne.ia4)then c O = C - NH2 if(lcc)then c O = C(-C) - NH2 call swallow(nat,us,iah,ipeh, 1 ia,ia2,ia3,ia4,ia5,ia6,ia6,ia6) goto 9 endif endif 127 continue else c cis do 138 ix5=1,nt(ia3) ia5=bt(ia3,ix5) if(ity(ia5).eq.6.and.ia5.ne.ia2)then c O = C - N(-C) - H if(lcc)then c O = C(-C) - N(-C) - H call swallow(nat,us,iac,ipec, 1 ia,ia2,ia3,ia4,ia5,ia6,ia6,ia6) goto 9 endif endif 138 continue endif endif 126 continue do 129 ix4=1,nt(ia3) ia4=bt(ia3,ix4) if(ity(ia4).eq.6.and.ia4.ne.ia2)then c ia ia2 ia3 ia4 c O = C - N - C do 131 ix5=1,nt(ia3) ia5=bt(ia3,ix5) if(ity(ia5).eq.6.and.ia5.ne.ia2.and.ia5.ne.ia4)then c O = C - N(-C) - C c ensure that ia5 is the main chain continuation carbon: lcon=.false. do 136 ix7=1,nt(ia5) ia7=bt(ia5,ix7) if(ity(ia7).eq.6)then do 137 ix8=1,nt(ia7) ia8=bt(ia7,ix8) 137 if(ity(ia8).eq.8)lcon=.true. endif 136 continue if(.not.lcon)then ix7=ia4 ia4=ia5 ia5=ix7 endif if(lcc)then c O = C(-C) - N(-C) - C do 132 ix7=1,nt(ia5) ia7=bt(ia5,ix7) if(ity(ia7).eq.6)then do 133 ix8=1,nt(ia7) ia8=bt(ia7,ix8) if(ity(ia8).eq.6)then do 134 ix9=1,nt(ia8) ia9=bt(ia8,ix9) if(ia9.eq.ia4)then c proline amide call swallow(nat,us,iap,ippp, 1 ia,ia2,ia3,ia4,ia5,ia6,ia7,ia8) goto 9 endif 134 continue endif 133 continue endif 132 continue endif endif 131 continue endif 129 continue c unknown O = C - N group call swallow(nat,us,iax,ipex,ia,ia2,ia3,0,0,0,0,0) endif 125 continue endif 124 continue endif 9 continue if(iag.gt.0)then open(11,file='AMIDE.LST') write( 6,*)iag,' trans O = C(-C) - N(-C) - H amide groups:' write(11,*)iag,' trans O = C(-C) - N(-C) - H amide groups:' do 4 ia=1,iag write(11,6009)ia,(ipep(ia,ix),ix=1,8) 4 write(6,6009)ia,(ipep(ia,ix),ix=1,6) 6009 format(i4,4x,11i5) if(iac.eq.0)close(11) endif if(iac.gt.0)then open(11,file='AMIDE.LST') write( 6,*)iac,' cis O = C(-C) - N(-C) - H amide groups:' write(11,*)iac,' cis O = C(-C) - N(-C) - H amide groups:' do 71 ia=1,iac write(11,6009)ia,(ipec(ia,ix),ix=1,8) 71 write(6,6009)ia,(ipec(ia,ix),ix=1,6) close(11) endif if(iap.gt.0)then open(12,file='PROLINE.LST') write(6,*)iap,' proline O = C(-C) - N (amide) groups:' write(12,*)iap,' proline O = C(-C) - N (amide) groups:' do 11 ia=1,iap write(12,6009)ia,(ippp(ia,ix),ix=1,8) 11 write(6,6009)ia,(ippp(ia,ix),ix=1,3) close(12) endif if(iah.gt.0)then open(13,file='TERMINAL.LST') write(6,*)iah,' terminal O = C(-C) - NH2 amide groups:' write(13,*)iah,' terminal O = C(-C) - NH2 amide groups:' do 12 ia=1,iah write(13,6009)ia,(ipeh(ia,ix),ix=1,8) 12 write(6,6009)ia,(ipeh(ia,ix),ix=1,3) close(13) endif if(iax.gt.0)then open(14,file='UNKNOWN.LST') write(6,*)iax,' unknown O = C - N (amide) groups:' write(14,*)iax,' unknown O = C - N (amide) groups:' do 10 ia=1,iax write(14,6009)ia,(ipex(ia,ix),ix=1,8) 10 write(6,6009)ia,(ipex(ia,ix),ix=1,3) close(14) endif c put all amide groups together (skip strange): iagg=0 do 101 ia=1,iag iagg=iagg+1 do 101 ix=1,8 101 itpp(iagg,ix)=ipep(ia,ix) do 111 ia=1,iac iagg=iagg+1 do 111 ix=1,8 111 itpp(iagg,ix)=ipec(ia,ix) do 103 ia=1,iah iagg=iagg+1 do 103 ix=1,8 103 itpp(iagg,ix)=ipeh(ia,ix) do 105 ia=1,iap iagg=iagg+1 do 105 ix=1,8 105 itpp(iagg,ix)=ippp(ia,ix) write(6,*)iagg,' amides in total' return end subroutine anala(k,nat,p,iagg,r,nnn,av) implicit none integer*4 nat,p(nat,8),iagg,nnn(100),im,nc(10),i,j,ic,ia,k real*8 r(*),a(10,3),d,tol,df,dm,o(10),av(3),fia,psia,wa,w,ta, 1fi,psi c nnn( 1) ... number of analyzed amide groups c nnn( 2) ... number of analyzed fi, psi angles c nnn( 3) ... number of cis amide groups c nnn( 4) ... number of trans amide groups c nnn(11) ... number of alpha-helices c nnn(12) ... number of left alpha-helices c nnn(13) ... number of betas c nnn(14) ... number of 3_10 c nnn(15) ... number of PPI c nnn(16) ... number of PPII c c i aC(6) --- Co(2) --- Nh(3) --- aC(5) c j aC(6) --- Co(2) --- Nh(3) --- aC(5) c c alpha-helix: a(1,1)= -57.0d0 a(1,2)= -47.0d0 a(1,3)= 180.0d0 c a(1,1)= -61.0d0 c a(1,2)= -41.0d0 c a(1,3)= 180.0d0 c left alpha-helix a(2,1)= 57.0d0 a(2,2)= 47.0d0 a(2,3)= 180.0d0 c beta: a(3,1)=-139.0d0 a(3,2)= 120.0d0 a(3,3)= 180.0d0 c 310 helix: a(4,1)= -74.0d0 a(4,2)= -4.0d0 a(4,3)= 180.0d0 c polyproline I: a(5,1)= -83.0d0 a(5,2)= 152.0d0 a(5,3)= 0.0d0 c polyproline II: a(6,1)= -78.0d0 a(6,2)= 149.0d0 a(6,3)= 180.0d0 tol=100.0d0 if(iagg.eq.0)return c average angles in this snapshot: fia=0.0d0 psia=0.0d0 wa=0.0d0 c number of each conformations and percentage in this snapshot: do 2 i=1,10 o(i)=0.0d0 2 nc(i)=0 ia=0 do 1 i=1,iagg nnn(1)=nnn(1)+1 w=ta(p(i,6),p(i,2),p(i,3),p(i,5),r) call adda(wa,w,i) if(dabs(w).lt.90.0d0)then nnn(3)=nnn(3)+1 else nnn(4)=nnn(4)+1 endif do 1 j=1,iagg if(p(i,5).eq.p(j,6))then nnn(2)=nnn(2)+1 ia=ia+1 fi=ta(p(i,2),p(i,3),p(i,5),p(j,2),r) psi=ta(p(i,3),p(i,5),p(j,2),p(j,3),r) call adda(fia,fi,ia) call adda(psia,psi,ia) c determine closest conformation: dm=1.0d9 im=1 c six standard conformations: do 5 ic=1,6 c distance from standard conformation: d=df(ic,fi,psi,w,a) if(d.lt.dm)then dm=d im=ic endif 5 continue if(dm.lt.tol)then nc(im)=nc(im)+1 else c unknown conformation: nc(7)=nc(7)+1 endif endif 1 continue if(ia.ne.0)then fia=fia/dble(ia) psia=psia/dble(ia) do 3 i=1,7 3 o(i)=nc(i)/dble(ia)*100.0d0 endif wa=wa/dble(iagg) do 4 i=1,7 4 nnn(10+i)=nnn(10+i)+nc(i) call cper(fia) call cper(psia) call cper(wa) write(6,600)k,fia,psia,wa,(o(i),i=1,7) write(33,600)k,fia,psia,wa,(o(i),i=1,7) 600 format(i9,3f8.2,7f7.2) call adda(av(1),fia,k) call adda(av(2),psia,k) call adda(av(3),wa,k) return end function ta(l1,l2,l3,l4,r) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION XX(3),YY(3),A(3),P(3),E0(3),B0(3), 1AP(3),r(*) PI=3.141592653589D0 do 23 ix=1,3 A(ix) =r(ix+3*(l3-1))-r(ix+3*(l4-1)) P(ix) =r(ix+3*(l1-1))-r(ix+3*(l2-1)) 23 E0(ix)=r(ix+3*(l2-1))-r(ix+3*(l3-1)) call vp(B0,E0,A) SSB0B0=sp(B0,B0) SSE0E0=sp(E0,E0) DO 8 KK=1,3 B0(KK)=B0(KK)/DSQRT(SSB0B0) 8 E0(KK)=E0(KK)/DSQRT(SSE0E0) SSAE0=sp(A,E0) SSPE0=sp(P,E0) DO 9 KK=1,3 XX(KK)=A(KK)-E0(KK)*SSAE0 9 YY(KK)=P(KK)-E0(KK)*SSPE0 SSXXXX=sp(XX,XX) SSYYYY=sp(YY,YY) SSXXYY=sp(XX,YY) call vp(AP,XX,YY) SSE0AP=sp(AP,E0) SIGN=1.0D0 IF (dabs(SSE0AP).gt.1.0D-9)SIGN=-SSE0AP/ABS(SSE0AP) CA=0.0D0 IF (dabs(SSXXXX).gt.1.0D-9.AND.dabs(SSYYYY).gt.1.0D-9) 1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY) ANGLE=0.0D0 IF (CA.LE.-1.0D0) ANGLE=PI IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA) ANGLE=180.0D0* ANGLE/PI ANGLE=SIGN*ANGLE if(ANGLE.lt.-180.0d0)ANGLE= 360.0d0+ANGLE if(ANGLE.gt.+180.0d0)ANGLE=-360.0d0+ANGLE ta=ANGLE return end function df(ic,fi,psi,wi,a) integer*4 ic real*8 fi,psi,wi,f,p,w,df,a(10,3) f=fi if(dabs(f+360.0d0-a(ic,1)).lt.dabs(f-a(ic,1)))f=f+360.0d0 if(dabs(f-360.0d0-a(ic,1)).lt.dabs(f-a(ic,1)))f=f-360.0d0 p=psi if(dabs(p+360.0d0-a(ic,2)).lt.dabs(p-a(ic,2)))p=p+360.0d0 if(dabs(p-360.0d0-a(ic,2)).lt.dabs(p-a(ic,2)))p=p-360.0d0 w=wi if(dabs(w+360.0d0-a(ic,3)).lt.dabs(w-a(ic,3)))w=w+360.0d0 if(dabs(w-360.0d0-a(ic,3)).lt.dabs(w-a(ic,3)))w=w-360.0d0 df=(f-a(ic,1))**2+(p-a(ic,2))**2+(w-a(ic,3))**2 return end subroutine swallow(n,u,i,j,i1,i2,i3,i4,i5,i6,i7,i8) implicit none integer n,i,j(n,8),i1,i2,i3,i4,i5,i6,i7,i8,u(*),ii i=i+1 j(i,1)=i1 j(i,2)=i2 j(i,3)=i3 j(i,4)=i4 j(i,5)=i5 j(i,6)=i6 j(i,7)=i7 j(i,8)=i8 do 1 ii=1,8 1 u(j(i,ii))=0 return end function sp(a,b) IMPLICIT none real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine vp(A,X,Y) implicit none real*8 A(*),X(*),Y(*) A(1)=X(2)*Y(3)-X(3)*Y(2) A(2)=X(3)*Y(1)-X(1)*Y(3) A(3)=X(1)*Y(2)-X(2)*Y(1) return end subroutine adda(wa,w,i) real*8 wa,w,ta integer*4 i c temporary average: ta=0.0d0 if(i-1.gt.0)ta=wa/dble(i-1) if(dabs(w-360.0d0-ta).lt.dabs(w-ta))then wa=wa+w-360.0d0 else if(dabs(w+360.0d0-ta).lt.dabs(w-ta))then wa=wa+w+360.0d0 else wa=wa+w endif endif return end subroutine cper(w) real*8 w if(w.gt. 180.0d0)w=w-360.0 if(w.lt.-180.0d0)w=w+360.0 return end