// // Make Ramachandran plot for the first Structure Cluster // // (C) Copyright 1996, Carlsberg Laboratory, Dept. of Chemistry // // PAGE DIGITS 6 DIM tmp$ OF 1 PRINT "Input structure cluster number: ", INPUT stc_id% DIM stc. OF stc_rec DIM atc. OF atc_rec DIM a1. OF atc_rec DIM a2. OF atc_rec DIM a3. OF atc_rec DIM a4. OF atc_rec TRAP CAT_STC_GET stc_id%,stc. HANDLER PRINT "Could not get Structure Cluster ",stc_id% END ENDTRAP natoms%:=stc.stc_natoms% // Count no. of residues: resmin%:=999 resmax%:=-999 st_id%:=CAT_STC_ST_FIRST(stc_id%) atc_id%:=CAT_ST_ATC_FIRST(st_id%) PRINT st_id%;atc_id% FOR i%:=1 TO natoms% DO CAT_ATC_GET atc_id%,atc. IF atc.atc_resno%resmax% THEN resmax%:=atc.atc_resno% atc_id%:=CAT_ST_ATC_NEXT#(atc_id%,st_id%) NEXT i% PRINT "resmin: ",resmin%," resmax: ",resmax% // Initialize graphics: PDEVICE "ps" PINIT "/indy2/mk/rama/nkt.ps",devxmin,devxmax,devymin,devymax IF devxmax>devymax THEN devxmax:=devymax IF devymax>devxmax THEN devymax:=devxmax xlow:=-180 xhigh:=180 ylow:=-180 yhigh:=180 // set up the coordinate system. charsize:=0.3 msize:=charsize*0.3 ntick%:=10 // plotting viewport, in cm: vxmin:=devxmin vymin:=devymin vxmax:=devxmax vymax:=devymax //title: vymax:=vymax-2*charsize //labels: vymin:=vymin+2.1*charsize vxmin:=vxmin+4*charsize //depends upon no. of digits PNOCLIP //border PCOLOR 1 PDENSITY 1 PTHICK 1 PMOVE vxmin,vymin PDRAW vxmax,vymin PDRAW vxmax,vymax PDRAW vxmin,vymax PDRAW vxmin,vymin //axes scale2(xlow,xhigh,ntick%,xlowa,xhigha,nxaxis,xtick,xndec%) scale2(ylow,yhigh,ntick%,ylowa,yhigha,nyaxis,ytick,yndec%) PRINT "xlow: ",xlow," xhigh: ",xhigh PRINT "xlowa: ",xlowa," xhigha: ",xhigha," xtick: ",xtick PRINT "ylow: ",ylow," yhigh: ",yhigh PRINT "ylowa: ",ylowa," yhigha: ",yhigha," ytick: ",ytick PCOLOR 1 PDENSITY 1 PTHICK 1 PTEXTSTYLE charsize,6 x:=xlowa WHILE x<=xhigha DO IF x>=xlow AND x<=xhigh THEN PTEXT xtrans(x),vymin-0.5*charsize,STR$(x) PMOVE xtrans(x),vymin PDRAW xtrans(x),vymin+charsize/2 PMOVE xtrans(x),vymax PDRAW xtrans(x),vymax-charsize/2 ENDIF x:=x+xtick ENDWHILE PTEXTSTYLE charsize,8 y:=ylowa WHILE y<=yhigha DO IF y>=ylow AND y<=yhigh THEN PTEXT vxmin-1*charsize,ytrans(y),STR$(y) PMOVE vxmin,ytrans(y) PDRAW vxmin+charsize/2,ytrans(y) PMOVE vxmax,ytrans(y) PDRAW vxmax-charsize/2,ytrans(y) ENDIF y:=y+ytick ENDWHILE PCLIP vxmin,vymin,vxmax,vymax // Plot allowable regions: PCOLOR 2 PMOVE xtrans(-156.6),ytrans(91.3) PDRAW xtrans(-70.4),ytrans(91.3) PDRAW xtrans(-54.7),ytrans(112.8) PDRAW xtrans(-54.7),ytrans(173.2) PDRAW xtrans(-136.95),ytrans(173.2) PDRAW xtrans(-136.9),ytrans(155.8) PDRAW xtrans(-156.5),ytrans(135.6) PDRAW xtrans(-156.5),ytrans(91.3) PMOVE xtrans(-180),ytrans(42.9) PDRAW xtrans(-140.8),ytrans(16.1) PDRAW xtrans(-86),ytrans(16.1) PDRAW xtrans(-74.3),ytrans(45.6) PDRAW xtrans(-74.3),ytrans(72.5) PDRAW xtrans(-44.3),ytrans(102) PDRAW xtrans(-44.3),ytrans(161.1) PDRAW xtrans(-46.9),ytrans(179.9) PMOVE xtrans(-180),ytrans(-34.9) PDRAW xtrans(-164.3),ytrans(-42.9) PDRAW xtrans(-133),ytrans(-42.9) PDRAW xtrans(-109.5),ytrans(-32.2) PDRAW xtrans(-106.9),ytrans(-21.4) PDRAW xtrans(-44.3),ytrans(-21.4) PDRAW xtrans(-44.3),ytrans(-71.1) PDRAW xtrans(-180),ytrans(-71.1) PMOVE xtrans(-156.5),ytrans(-60.4) PDRAW xtrans(-54.7),ytrans(-60.4) PDRAW xtrans(-54.7),ytrans(-40.2) PDRAW xtrans(-100.4),ytrans(-40.2) PDRAW xtrans(-123.9),ytrans(-51) PDRAW xtrans(-156.5),ytrans(-51) PDRAW xtrans(-156.5),ytrans(-60.4) PMOVE xtrans(-180),ytrans(-163.8) PDRAW xtrans(-75.6),ytrans(-163.8) PDRAW xtrans(-46.9),ytrans(-180) PMOVE xtrans(62.6),ytrans(14.7) PDRAW xtrans(62.6),ytrans(96.7) PDRAW xtrans(45.6),ytrans(79.2) PDRAW xtrans(45.6),ytrans(26.8) PDRAW xtrans(62.6),ytrans(14.7) PCOLOR 1 PDENSITY 1 PTHICK 0 PRECTANGLESTYLE 2,1 DIM id_c%(resmin%:resmax%) DIM id_n%(resmin%:resmax%) DIM id_ca%(resmin%:resmax%) DIM gly%(resmin%:resmax%) OPEN FILE 1,"/indy2/mk/rama/pils.lst",WRITE st_id%:=CAT_STC_ST_FIRST(stc_id%) // use icolor to change colors between substructures //icolor%=1 WHILE st_id%<>-1 DO //pcolor icolor% //print "Using icolor: ",icolor% //icolor%=icolor%+1 FOR i%:=resmin% TO resmax% DO id_c%(i%):=-1 id_n%(i%):=-1 id_ca%(i%):=-1 NEXT i% atc_id%:=CAT_ST_ATC_FIRST(st_id%) FOR i%:=1 TO natoms% DO CAT_ATC_GET atc_id%,atc. j%:=atc.atc_resno% gly%(j%):=atc.atc_resname!="GLY" IF atc.atc_aname!="C" THEN IF id_c%(j%)<>-1 THEN PRINT "Twice C: ",id_c%(j%),atc_id% ELSE id_c%(j%):=atc_id% ENDIF ELIF atc.atc_aname!="N" THEN IF id_n%(j%)<>-1 THEN PRINT "Twice N: ",id_n%(j%),atc_id% ELSE id_n%(j%):=atc_id% ENDIF ELIF atc.atc_aname!="CA" THEN IF id_ca%(j%)<>-1 THEN PRINT "Twice CA: ",id_ca%(j%),atc_id% ELSE id_ca%(j%):=atc_id% ENDIF ENDIF atc_id%:=CAT_ST_ATC_NEXT#(atc_id%,st_id%) NEXT i% FOR i%:=resmin%+1 TO resmax%-1 DO // PRINT i%;id_ca%(i%);id_c%(i%);id_n%(i%) // phi: c(i-1),n(i),ca(i),c(i) IF id_c%(i%-1)<>-1 AND id_n%(i%)<>-1 AND id_ca%(i%)<>-1 AND id_c%(i%)<>-1 THEN CAT_ATC_GET id_c%(i%-1),a1. CAT_ATC_GET id_n%(i%),a2. CAT_ATC_GET id_ca%(i%),a3. CAT_ATC_GET id_c%(i%),a4. phi:=angle(a1.,a2.,a3.,a4.) ELSE phi:=999 ENDIF // psi: n(i),ca(i),c(i),n(i+1) IF id_n%(i%)<>-1 AND id_ca%(i%)<>-1 AND id_c%(i%)<>-1 AND id_n%(i%+1)<>-1 THEN CAT_ATC_GET id_n%(i%),a1. CAT_ATC_GET id_ca%(i%),a2. CAT_ATC_GET id_c%(i%),a3. CAT_ATC_GET id_n%(i%+1),a4. psi:=angle(a1.,a2.,a3.,a4.) ELSE psi:=999 ENDIF PRINT FILE 1: st_id%;i%;gly%(i%);phi;psi IF phi<>999 AND psi<>999 THEN // plot the point: x:=xtrans(phi) y:=ytrans(psi) IF gly%(i%) THEN PRINT "x", PMOVE x-msize,y-msize PDRAW x+msize,y+msize PMOVE x-msize,y+msize PDRAW x+msize,y-msize PTEXTSTYLE charsize*0.4,2 PTEXT x+msize,y,STR$(i%) ELSE PRINT "o", PMOVE x-msize,y-msize PDRAW x+msize,y-msize PDRAW x+msize,y+msize PDRAW x-msize,y+msize PDRAW x-msize,y-msize PTEXTSTYLE charsize*0.4,5 PTEXT x,y,STR$(i%) ENDIF ENDIF NEXT i% PRINT st_id%:=CAT_STC_ST_NEXT#(st_id%,stc_id%) ENDWHILE CLOSE FILE 1 PNOCLIP PRINT "Press return to finish: ", INPUT tmp$ PPAGE PFINISH FUNC angle(REF a1.,REF a2.,REF a3.,REF a4.) CLOSED xi:=a1.atc_x# xj:=a2.atc_x# xk:=a3.atc_x# xl:=a4.atc_x# yi:=a1.atc_y# yj:=a2.atc_y# yk:=a3.atc_y# yl:=a4.atc_y# zi:=a1.atc_z# zj:=a2.atc_z# zk:=a3.atc_z# zl:=a4.atc_z# rijx:=xi-xj rijy:=yi-yj rijz:=zi-zj rjkx:=xj-xk rjky:=yj-yk rjkz:=zj-zk rklx:=xk-xl rkly:=yk-yl rklz:=zk-zl ax:=rijy*rjkz-rijz*rjky ay:=rijz*rjkx-rijx*rjkz az:=rijx*rjky-rijy*rjkx bx:=rjky*rklz-rkly*rjkz by:=rjkz*rklx-rklz*rjkx bz:=rjkx*rkly-rklx*rjky cx:=rjky*az-rjkz*ay cy:=rjkz*ax-rjkx*az cz:=rjkx*ay-rjky*ax rar:=ax*ax+ay*ay+az*az IF rar<0.0001 THEN rar:=0.0001 rar:=1/SQR(rar) rbr:=bx*bx+by*by+bz*bz IF rbr<0.0001 THEN rbr:=0.0001 rbr:=1/SQR(rbr) rcr:=cx*cx+cy*cy+cz*cz IF rcr<0.0001 THEN rcr:=0.0001 rcr:=1/SQR(rcr) ax:=ax*rar ay:=ay*rar az:=az*rar bx:=bx*rbr by:=by*rbr bz:=bz*rbr cx:=cx*rcr cy:=cy*rcr cz:=cz*rcr cp:=ax*bx+ay*by+az*bz sp:=cx*bx+cy*by+cz*bz IF cp<-1 THEN cp:=-1 IF cp>1 THEN cp:=1 phi:=ACS(cp) IF sp>=0 THEN phi:=-phi RETURN phi/PI*180 ENDFUNC angle PROC scale2(min1,max1,ntick1%,REF min2,REF max2,REF ntick2,REF tick,REF ndec%) d:=(max1-min1)/ntick1% ld:=LOG(d)/LOG(10) ild:=INT(ld) ndec%:=0 IF ild<0 THEN ndec%:=-ild fld:=10^(ld-ild) tick:=10^ild IF fld>5 THEN tick:=tick*10 ndec%:=ndec%-1 IF ndec%<0 THEN ndec%:=0 ELIF fld>2 THEN tick:=tick*5 ELIF fld>1 THEN tick:=tick*2 ENDIF min2:=FLR(min1/tick)*tick max2:=FLR(max1/tick+0.99999)*tick ENDPROC scale2 FUNC xtrans(x) RETURN vxmin+(x-xlow)/(xhigh-xlow)*(vxmax-vxmin) ENDFUNC xtrans FUNC ytrans(y) RETURN vymin+(y-ylow)/(yhigh-ylow)*(vymax-vymin) ENDFUNC ytrans