program HELANAL C**************************************************************************** C * C PROGRAM TO CHARACTERISE THE GEOMETRIES OF ALPHA HELICES IN PROTEINS. * C * C AUTHORS: Sandeep Kumar and Prof. Manju Bansal, * C MBU, Indian Institute of Science, Bangalore 560012, India * C * C e-mail address: mb@mbu.iisc.ernet.in * C * C Method outlined in the following paper: * C * C Kumar, S. and Bansal, M. (1996). Structural and Sequence Characteristics * C of Long Alpha Helices in Globular Proteins. Biophysical J.,71, 1574-1586. * C * C If any Bug is found, please report it to the authors. * C * C SUMMARY OF THE ALGORITHM * C * C Geometry of an alpha helix is characterised in terms of the angles between* C local helix axes and the path traced by the local helix origins, * C calculated using the procedure of Sugeta and Miyazawa (1967), for every * C set of four contiguous C-alpha atoms, and sliding this window over the * C length of the helix in steps of one C-alpha atom. * C Matrix M(I, J) contains the bending angles between local helix axes I & J * C which are used to characterize the overall geometry of the helix. * C The local helix origins trace out the path described by the helix in 3-D * C space. These origins are reoriented in X-Y plane and the reoriented * C points are used to fit a circle and a line by least squares method. * C Unit twist and unit height of the alpha helix are also calculated. * C * C A maximum of 5000 helices, each with 100 or less number of residues can * C be analysed at one time. * C * C In order to analyse a larger number of helices, increase the value of 'i' * C in the following statement: * C * C parameter (i=5000) in the main program. * C * C In order to analyse helices with lengths greater than 100 residues, * C increase the dimensions of the appropriate variables in the dimension * C statements of the main program and subroutines. * C * C This program has several options for input and is fully interactive. * C * C INPUT FILES : * C * C Input files to this program can be different depending upon the options * C chosen. This program can directly read HELIX records in one or more PDB * C files. In order to analyse helices found in helix records of a PDB file, * C give the PDB file name as input to the program. In order to analyse the * C helices found in the HELIX records of more than one PDB files, give the * C PDB file names sequentially when prompted by the program or input the * C name of the file containing the PDB files names in format (5x,a11) as the * C input to the program. * C * C HELIX records can be read from PDB files or can be read from a file * C containing information about the helix start/end residues written in * C the same format as PDB HELIX records or in a different format, in * C which case the format has to be keyed in when running the program. * C In the last case, files other than the PDB files or files containing * C C-alpha coordinates in the PDB format, can also be given as input to * C the program, upon specifying their format. * C * C * C All PDB files and other input files (if any) should be in the same * C directory. * C * C OUTPUT FILES : * C * C For each run of HELANAL on file(s) containing alpha helices with length * C greater than or equal to 9 residues, the following output files are * C created. * C * C RUN.ANS contains the questions and their answers during a run of HELANAL. * C * C HELINFO.OUT file created only when the HELIX records in the PDB files are * C used for information on helix start/end residues, contains these helix * C records. * C * C HELCA.OUT contains Coordinates of the C-alpha atoms constituting the * C helices. * C * C AXES.OUT contains the local helix axes fitted to 4 consecutive C-alpha * C atoms along with a matrix M(I, J) whose elements are the angles between * C local helix axes I and J. * C * C ANGLE.OUT contains the angle between successive local helix axes. It also * C lists mean bending angle and maximum bending angle for each helix. * C * C ORIGIN.OUT contains the local helix origins for the helix along with the * C statistics obtained by fitting least square plane, circle and line to the * C local helix origins. * C * C NH.OUT contains unit height and unit twist for every turn of the helix as * C well as average unit height and average unit twist for the whole helix. * C * C ***.PRM contains summary of various parameters obtained by HELANAL. * C * C ***.TAB contains the parameters from ***.PRM file in a tabular form along * C the overall geometry assignment. * C * C**************************************************************************** c parameter (i=5000) character*80 line, prmfile,forma,formacor,tabfile character*40 inpfile character*11 fnm,pdbfile character*3 rb,re,helnm,res character*4 atom character*1 cb,ce,ans1,ans2,ans3,ans4,ans5,chain character*1 geom c dimension rb(i),cb(i),nrb(i) dimension re(i),ce(i),nre(i),fnm(i) dimension helnm(i),pdbfile(i) dimension xca(100),yca(100),zca(100) dimension ox(100),oy(100),oz(100) dimension twist(100),rnou(100),angl(100) dimension vec12(3),vec23(3),vec34(3),dv13(3),dv24(3) dimension uloc1(100),uloc2(100),uloc3(100),kount(100) dimension rad(3),ibendangl(100,100),height(100) c c FORMATS USED c 1 format (/1x,'Do You wish to use HELIX records in PDB files?') 2 format (/1x,'Key in name of the INPUT file containing PDB file n -ames in (5x,a11) format',/1x,'e.g. ***.fil :') 3 format (/1x,'Key in name of the INPUT file containing helix st -art/end information',/1x,'e.g. ***.inp :') 4 format (1x,'C-ALPHA ATOMS',/) 5 format (30x,3f8.3) 6 format (1x,'LOCAL HELIX AXES FOR CA ATOMS (I - (I+3))',3x, - 'L',7x,'M',7x,'N',/) 7 format (5x,2(i5,5x),15x,3f8.3) 8 format (5x,3(f8.3)) 9 format (1x,'LOCAL HELIX ORIGINS FOR CA ATOMS (I - (I+3))' - ,2x,'X',7x,'Y',7x,'Z',/) 10 format(/1x,'Is the information in the same format as HELIX recor -ds in PDB files?') 11 format (5x,'O ',3(i5,2x),15x,3f8.3) 12 format (1x,'LOCAL BENDING ANGLES CALCULATED FROM LOCAL HELIX -AXES FITTED TO CA ATOMS'/1x,'((I-3) - I) AND (I - (I+3))'/) 13 format (4(i5,2x),5x,f8.3) 14 format (1x,'LOCAL HELIX PARAMETERS CALCULATED FOR CA ATOMS (I - - (I+3))') 16 format (1x,'MATRIX M(I,J) FOR ANGLES BETWEEN LOCAL HELIX AXES - I AND J') 17 format (i3,1x,40i3) 18 format (4x,40i3) 19 format(1x,'AVERAGE TWIST (DEG.)',6x,f8.2,2x,'S.D.',f8.2, - 5x,'ADV',f8.2) 20 format(1x,'AVERAGE N',17x,f8.2,2x,'S.D.',f8.2,5x,'ADV' - ,f8.2) 21 format(1x,'AVERAGE UNIT HEIGHT (ANG.)',f8.2,2x,'S.D.',f8.2, - 5x,'ADV',f8.2) 22 format(/1x,'If not, key in the input data format:'/1x,'(Specify fi -elds for coordinate file name, helix name, name, chain and number' -/1x,'for first and last residues. Use three letter code for residu -e names.)'/) 23 format(a80) 24 format (/1x,'Do you wish to type in the PDB file names?') 25 format(5x,'N------>C'/) 26 format(1x,'MAXIMUM BENDING ANGLE',5x,f8.3) 27 format(1x,'MEAN BENDING ANGLE',8x,f8.3,2x,'S.D.',f8.3, - 5x,'ADV',f8.3) 28 format(1x,'SUMMARY OF PARAMETERS OBTAINED BY HELANAL'/) 29 format(1x,'MEAN BENDING ANGLE (DEG.)',15x,f8.3,2x,'S.D.', - f8.3,2x,'ADV',f8.3) 30 format(1x,'MAXIMUM BENDING ANGLE (DEG.)',12x,f8.3) 31 format (42x,'TWIST',6x, 'N',4x,'HEIGHT'/) 32 format(a80) 33 format(13x,a3,1x,a3,1x,a1,1x,i3,4x,3f8.3) 34 format(11x,a3,1x,2(a3,1x,a1,1x,i4,2x),33x,a4) 35 format(/1x,a11,2x,'HELIX',2x,a3,1x,2(a3,1x,a1,1x,i3,4x)) 36 format(1x,'AVERAGE N',31x,f8.3,2x,'S.D.',f8.3,2x,'ADV' - ,f8.3) 37 format(1x,'AVERAGE UNIT HEIGHT (ANG.)',14x,f8.3,2x,'S.D.', - f8.3,2x,'ADV',f8.3) 38 format(1x,'LEAST SQUARES PLANE, CIRCLE AND LINE FITTED TO LOC -AL HELIX ORIGINS') 39 format(5x,a11) 40 format(/1x,'Key in a PDB file name :') 42 format(/1x,'QUESTIONS AND THEIR ANSWERS'/) 43 format(/1x,'All questions and their answers are written into the - file "RUN.ANS". '/) 44 format(/1x,'Is',1x,a11,1x,'a PDB file?') 45 format(/1x,'If not, key in the input data format:',/1x,'(Specify - fields for atom name (4 characters), res. name, chain, res. numbe -r',/1x,'and XYZ coordinates. Use three letter code for residue nam -es.)'/) 46 format(a40) 47 format(1x,'HELICES TAKEN FROM :',1x,a40) 48 format(/1x,'WRITING THE RESULTS TO :',1x,a40,/,1x,'AND IN TABULA -R FORM TO :',1x,a40) 49 format(1x,a40) 50 format(1x,'HELIX IGNORED : IT CONTAINS LESS THAN 9 RESIDUES.') 51 format(/1x,'Do you wish to use another PDB file?') 52 format(/1x,a11,' DOES NOT CONTAIN HELIX RECORDS') 53 format(3x,'File Name',6x,'Helix',5x,'n',4x,'h',3x,'Aver',2x, - 'Max',2x,'Radius',2x,'rmsC',2x,'rmsL',2x,'r2',1x,'Geometry', -/38x,'BA',3x,'BA',/) 54 format(1x, 'Helix : Chain identifier, number and name of fi -rst and last residues'/13x, 'of the helix.') 55 format (1x,'n : Average number of residues per turn of -the helix.') 56 format(1x,'h : Average unit height in the helix (Angstroms -).') 57 format(1x,'Aver BA : Average Bending angle between successive - local helix axes (Deg.).') 58 format(1x,'Max BA : Maximum Bending angle between succesive lo -cal helix axes (Deg.).') 59 format(1x,'Radius : Radius of least squares circle fitted to t -he local helix '/13x,'origins (Angstroms).') 60 format(1x,'rmsC : Root Mean Square Deviation for least squares - circle fitted to'/13x,'the local helix origins (Angstroms).') 61 format(1x,'rmsL : Root Mean Square Deviation for least squares - line fitted to the'/13x,'local helix origins (Angstroms). -') 62 format(1x,'r2 : Square of linear correlation coefficient -for least squares line'/13x,'fittted to the local helix 1origins. -') 63 format(1x,'Geometry : Overall geometry of the helix, Linear (L) -, Curved (C),'/13x,'Kinked (K) or unassigned (*).') 64 format(1x,'Std. Dev. : Standard Deviations of the average parame -ters for the helix'/13x,'in the previous line.') 66 format(/1x,'DESCRIPTION OF THE TABLE'/) 67 format(/'******************************************************* -************************') 68 format(/,5x,'NL=',i5,5x,'NC=',i5,5x,'NK=',i5,5x,'NA=',i5 -,5x,'NH=',i5) 69 format(1x,'NL, NC, NK: Number of helices with overall linear, curv -ed, kinked,',/1x,'NA, NH',6x,'unassigned geometry and total number - of helices.') nlinear=0 ncurved=0 nkinked=0 nogeom=0 c C OPENING INPUT AND OUTPUT FILES c open (unit=1,file='helca.out') open (unit=2,file='axes.out') open (unit=3,file='origin.out') open (unit=7,file='RUN.ANS') open (unit=8,file='angle.out') open (unit=9,file='nh.out') write(1,4) write(2,6) write(3,9) write(7,42) write(8,12) write(9,14) write(9,31) write(6,43) c ASKING THE FIRST QUESTION : USE HELIX RECORDS IN PDB FILE? 120 write (6,1) write (7,1) read(5,'(a1)')ans1 write(7,'(1x,a1)')ans1 if ((ans1.ne.'n').and.(ans1.ne.'N').and.(ans1.ne.'y') - .and.(ans1.ne.'Y')) then goto 120 endif c c READING HELIX RECORDS FROM PDB FILES c if ((ans1.eq.'y').or.(ans1.eq.'Y')) then c ASKING SECOND QUESTION : TYPE IN THE PDB FILE NAMES ? 119 write(6,24) write(7,24) read(5,'(a1)') ans2 write(7,'(1x,a1)') ans2 if ((ans2.ne.'n').and.(ans2.ne.'N').and.(ans2.ne.'y') - .and.(ans2.ne.'Y')) then goto 119 endif if((ans2.eq.'n').or.(ans2.eq.'N')) then 121 write(6,2) write(7,2) read(5,46) inpfile write(7,49) inpfile kkk = index(inpfile(1:40),':') if (kkk.eq.0) then kkk=1 else kkk=kkk+1 endif kkkk=index(inpfile(kkk:40),'.') kkkk=kkkk+kkk prmfile=inpfile(kkk:kkkk-1)//'prm' tabfile=inpfile(kkk:kkkk-1)//'tab' if((inpfile(kkkk:kkkk+2).eq.'ent').or.(inpfile(kkkk:kkkk+2) -.eq.'pdb').or.(inpfile(1:3).eq.'pdb')) goto 121 open (unit=12,file=inpfile,status='OLD',err=121) open (unit=11,file='helinfo.out') open (unit=10,file=prmfile) open (unit=13,file=tabfile) write(13,53) nfiles=0 do 501 j=1,i read(12,39,end=129,err=121) pdbfile(j) nfiles=nfiles+1 501 continue 129 continue do 502 l=1,nfiles open (unit=15,file=pdbfile(l),status='OLD',err=121) do 503 j=1,i read(15,32,end=131) line if (line(1:7).eq.'HELIX ') then write(11,32) line endif 503 continue 131 continue close (15) 502 continue rewind(11) write(1,47) inpfile write(2,47) inpfile write(3,47) inpfile write(8,47) inpfile write(9,47) inpfile write(10,28) write(10,47) inpfile write(6,48) prmfile,tabfile write(7,48) prmfile,tabfile goto 100 endif endif 155 if((ans1.eq.'y').or.(ans1.eq.'Y')) then if((ans2.eq.'y').or.(ans2.eq.'Y')) then 133 write(6,40) write(7,40) read(5,46) inpfile write(7,49) inpfile kkk = index(inpfile(1:40),':') if (kkk.eq.0) then kkk=1 else kkk=kkk+1 endif kkkk=index(inpfile(kkk:40),'.') kkkk=kkkk+kkk if((inpfile(kkkk:kkkk+2).ne.'ent').and.(inpfile(kkkk:kkkk+2) -.ne.'pdb')) goto 133 prmfile=inpfile(kkk:kkkk-1)//'prm' tabfile=inpfile(kkk:kkkk-1)//'tab' open (unit=12,file=inpfile,status='OLD',err=133) open (unit=11,file='helinfo.out') open (unit=10,file=prmfile) open (unit=13,file=tabfile) write(13,53) nlines = 0 do 497 j=1,i read(12,32,end=132) line if (line(1:7).eq.'HELIX ') then write(11,32) line nlines=nlines+1 endif 497 continue 132 continue if(nlines.eq.0) then write(6,52) inpfile write(7,52) inpfile goto 156 endif rewind(11) close (12) write(1,47) inpfile write(2,47) inpfile write(3,47) inpfile write(8,47) inpfile write(9,47) inpfile write(10,28) write(10,47) inpfile write(6,48) prmfile,tabfile write(7,48) prmfile,tabfile go to 100 endif endif c READING HELIX INFORMATION FROM A SEPARATE FILE if ((ans1.eq.'n').or.(ans1.eq.'N')) then 122 write (6,3) write (7,3) read(5,46) inpfile write(7,49) inpfile kkk = index(inpfile(1:40),':') if (kkk.eq.0) then kkk=1 else kkk=kkk+1 endif kkkk=index(inpfile(kkk:40),'.') kkkk=kkkk+kkk prmfile=inpfile(kkk:kkkk-1)//'prm' tabfile=inpfile(kkk:kkkk-1)//'tab' if((inpfile(kkkk:kkkk+2).eq.'ent').or.(inpfile(kkkk:kkkk+2) -.eq.'pdb').or.(inpfile(1:3).eq.'pdb')) goto 122 open (unit=11,file=inpfile,status='OLD',err=122) open (unit=10,file=prmfile) open (unit=13,file=tabfile) write(13,53) write(1,47) inpfile write(2,47) inpfile write(3,47) inpfile write(8,47) inpfile write(9,47) inpfile write(10,28) write(10,47) inpfile C ASKING THE THIRD QUESTION : IS THE FORMAT SAME AS HELIX RECORDS IN PDB FILES? c 123 write(6,10) write(7,10) read(5,'(a1)') ans3 write(7,'(1x,a1)') ans3 if ((ans3.ne.'n').and.(ans3.ne.'N').and.(ans3.ne.'y') - .and.(ans3.ne.'Y')) then goto 123 endif if ((ans3.eq.'N').or.(ans3.eq.'n')) then 127 rewind(11) write(6,22) write(7,22) read(5,'(a80)') forma write(7,'(1x,a79)') forma indexhel=0 do 498 j=1,i read(11,forma,end=126,err=127) fnm(j),helnm(j),rb(j),cb(j), 1 nrb(j),re(j),ce(j),nre(j) indexhel=indexhel+1 499 continue 498 continue 126 continue go to 125 endif c if((ans3.eq.'Y').or.(ans3.eq.'y')) then goto 100 endif endif c c READING HELIX INFORMATION c 100 indexhel=0 do 504 j=1,i read(11,34,end=99) helnm(j),rb(j),cb(j),nrb(j), 1 re(j),ce(j),nre(j),fnm(j)(4:7) do 505 k=5,7 if(fnm(j)(k:k).eq.'A') fnm(j)(k:k)='a' if(fnm(j)(k:k).eq.'B') fnm(j)(k:k)='b' if(fnm(j)(k:k).eq.'C') fnm(j)(k:k)='c' If(fnm(j)(k:k).eq.'D') fnm(j)(k:k)='d' if(fnm(j)(k:k).eq.'E') fnm(j)(k:k)='e' if(fnm(j)(k:k).eq.'F') fnm(j)(k:k)='f' if(fnm(j)(k:k).eq.'G') fnm(j)(k:k)='g' if(fnm(j)(k:k).eq.'H') fnm(j)(k:k)='h' if(fnm(j)(k:k).eq.'I') fnm(j)(k:k)='i' if(fnm(j)(k:k).eq.'J') fnm(j)(k:k)='j' if(fnm(j)(k:k).eq.'K') fnm(j)(k:k)='k' if(fnm(j)(k:k).eq.'L') fnm(j)(k:k)='l' if(fnm(j)(k:k).eq.'M') fnm(j)(k:k)='m' if(fnm(j)(k:k).eq.'N') fnm(j)(k:k)='n' if(fnm(j)(k:k).eq.'O') fnm(j)(k:k)='o' if(fnm(j)(k:k).eq.'P') fnm(j)(k:k)='p' if(fnm(j)(k:k).eq.'Q') fnm(j)(k:k)='q' if(fnm(j)(k:k).eq.'R') fnm(j)(k:k)='r' if(fnm(j)(k:k).eq.'S') fnm(j)(k:k)='s' if(fnm(j)(k:k).eq.'T') fnm(j)(k:k)='t' if(fnm(j)(k:k).eq.'U') fnm(j)(k:k)='u' if(fnm(j)(k:k).eq.'V') fnm(j)(k:k)='v' if(fnm(j)(k:k).eq.'W') fnm(j)(k:k)='w' if(fnm(j)(k:k).eq.'X') fnm(j)(k:k)='x' if(fnm(j)(k:k).eq.'Y') fnm(j)(k:k)='y' if(fnm(j)(k:k).eq.'Z') fnm(j)(k:k)='z' 505 continue fnm(j)(1:3)='pdb' fnm(j)(8:11)='.ent' indexhel=indexhel+1 504 continue 99 continue close (11) c WRITING HELIX INFORMATION ON OUTPUT FILES c 125 do 506 k1=1,indexhel len = nre(k1)-nrb(k1)+1 if (len.le.8) then write(6,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) write(7,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) write(6,50) write(7,50) else open(unit=4,file=fnm(k1)) write(1,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) write(2,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) write(3,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) write(8,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) write(9,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) if (fnm(k1).ne.fnm(k1-1)) write(10,67) write(10,35) fnm(k1),helnm(k1),rb(k1),cb(k1),nrb(k1), 1 re(k1),ce(k1),nre(k1) c c IS IT A PDB FILE? : FOURTH QUESTION c if ((ans1.eq.'n').or.(ans1.eq.'N')) then if(fnm(k1).ne.fnm(k1-1)) then 134 write(6,44) fnm(k1) write(7,44) fnm(k1) read(5,'(a1)') ans4 write(7,'(1x,a1)') ans4 if ((ans4.ne.'n').and.(ans4.ne.'N').and.(ans4.ne.'y') 1 .and.(ans4.ne.'Y')) then go to 134 endif if((ans4.eq.'y').or.(ans4.eq.'y')) then write(6,48) prmfile,tabfile write(7,48) prmfile,tabfile go to 135 endif if ((ans4.eq.'N').or.(ans4.eq.'n')) then 139 rewind(4) write(6,45) write(7,45) read(5,'(a80)') formacor write(7,'(1x,a79)') formacor j=1 137 continue read(4,32,end=136) line read(line,formacor,err=139) atom,res,chain,nores,xx,yy,zz if(chain.eq.cb(k1)) then if ((atom.eq.'CA ').or.(atom.eq.' CA').or. - (atom.eq.' CA ')) then if((nores.ge.nrb(k1)).and.(nores.le.nre(k1))) then kk=nores-nrb(k1)+1 xca(kk)=xx yca(kk)=yy zca(kk)=zz write(1,32) line endif endif endif j=j+1 go to 137 136 continue write(6,48) prmfile,tabfile write(7,48) prmfile,tabfile go to 138 endif endif endif c C EXTRACTING C-ALPHA ATOMS FROM PDB FILES c 135 read(4,32,end=999) line if ((line(1:4).eq.'ATOM').and.(line(14:15).eq.'CA')) then read(line,33) atom,res,chain,nores,xx,yy,zz if(chain.eq.cb(k1)) then if((nores.ge.nrb(k1)).and.(nores.le.nre(k1))) then kk=nores-nrb(k1)+1 xca(kk)=xx yca(kk)=yy zca(kk)=zz write(1,32) line endif endif endif go to 135 999 continue c c TAKING 4 CONSECUTIVE CA ATOMS c 138 kk3=kk-3 pi=180.0/acos(-1.0) do 508 l=1,kk3 kount(l)=l c c VECTORS JOINING CA ATOMS c vec12(1)=xca(l+1)-xca(l) vec12(2)=yca(l+1)-yca(l) vec12(3)=zca(l+1)-zca(l) c vec23(1)=xca(l+2)-xca(l+1) vec23(2)=yca(l+2)-yca(l+1) vec23(3)=zca(l+2)-zca(l+1) c vec34(1)=xca(l+3)-xca(l+2) vec34(2)=yca(l+3)-yca(l+2) vec34(3)=zca(l+3)-zca(l+2) c c DIFFERENCE OF VECTORS JOINING CA ATOMS c dv13(1)=vec12(1)-vec23(1) dv13(2)=vec12(2)-vec23(2) dv13(3)=vec12(3)-vec23(3) c dv24(1)=vec23(1)-vec34(1) dv24(2)=vec23(2)-vec34(2) dv24(3)=vec23(3)-vec34(3) c c DIRECTION OF THE LOCAL HELIX AXIS c call cross(dv13,dv24,el,em,en) write(2,7) l,l+3, el,em,en uloc1(l)=el uloc2(l)=em uloc3(l)=en c dmag=dot(dv13,dv13) dmag=sqrtf(dmag) emag=dot(dv24,dv24) emag=sqrtf(emag) costheta=dot(dv13,dv24)/(dmag*emag) c c TWIST OF THE LOCAL HELIX, and thus number of units per turn c twist(l)=acos(costheta)*pi rnou(l)=360.0/twist(l) c c RADIUS OF LOCAL HELIX CYLINDER c costheta1=1.0-costheta radmag=sqrtf(dmag*emag)/(2.0*costheta1) c c UNIT HEIGHT OF THE LOCAL HELIX c height(l)=vec23(1)*uloc1(l)+vec23(2)*uloc2(l)+vec23(3)*uloc3(l) write(9,7) l,l+3, twist(l),rnou(l),height(l) c c LOCAL HELIX ORIGIN c dv13(1)=dv13(1)/dmag dv13(2)=dv13(2)/dmag dv13(3)=dv13(3)/dmag c rad(1)=radmag*dv13(1) rad(2)=radmag*dv13(2) rad(3)=radmag*dv13(3) c ox(l)=xca(l+1)-rad(1) oy(l)=yca(l+1)-rad(2) oz(l)=zca(l+1)-rad(3) c dv24(1)=dv24(1)/emag dv24(2)=dv24(2)/emag dv24(3)=dv24(3)/emag c rad(1)=radmag*dv24(1) rad(2)=radmag*dv24(2) rad(3)=radmag*dv24(3) c ox(l+1)=xca(l+2)-rad(1) oy(l+1)=yca(l+2)-rad(2) oz(l+1)=zca(l+2)-rad(3) 508 continue c do 509 k=1,kk-2 if(abs(ox(k)).le.1.0E-04) ox(k) = 0.0000 if(abs(oy(k)).le.1.0E-04) oy(k) = 0.0000 if(abs(oz(k)).le.1.0E-04) oz(k) = 0.0000 write(3,11) k,k+1,k+2,ox(k),oy(k),oz(k) 509 continue write(3,38) c c LOCAL (k--k+3, k+3--k+6) BENDING ANGLES c kk6=kk3-3 do 510 k=1,kk6 angle=uloc1(k)*uloc1(k+3)+uloc2(k)*uloc2(k+3)+ - uloc3(k)*uloc3(k+3) if(abs(angle-1.0).le.1.0E-06) angle = 1.0 angl(k)=acos(angle)*pi write(8,13) k,k+3,K+3,k+6,angl(k) 510 continue c c MEAN BENDING ANGLE FOR THE WHOLE HELIX c call average(angl,kk6,avang,sdang,Advang) write(8,27) avang,sdang,Advang write(10,29) avang,sdang,Advang c c MAXIMUM BENDING ANGLE IN THE WHOLE HELIX c call max(angl,kk6,anglmax) write(8,26) anglmax write(10,30) anglmax c c LOCAL BENDING ANGLES MATRIX c do 511 j=1,kk3 do 512 k=1,kk3 proj=uloc1(j)*uloc1(k)+uloc2(j)*uloc2(k)+uloc3(j)*uloc3(k) if(abs(proj-1).le.1.0E-06) proj=1.0 bendangl=acos(proj)*pi ibend=int(bendangl) bend2=bendangl-float(ibend) c if(bend2.lt.0.5) then ibendangl(j,k)=int(bendangl) else ibendangl(j,k)=int(bendangl+1.0) endif c 512 continue 511 continue write(2,16) write(2,25) write(2,18) (kount(j),j=1,kk3) do 513 j=1,kk3 write(2,17) j,(ibendangl(j,k),k=1,kk3) 513 continue c c AVERAGE HELICAL PARAMETERS c c AVERAGE TWIST c call average(twist,kk3,av,sd,Adv) write(9,19) av,sd,Adv c c AVERAGE N c call average(rnou,kk3,avn,sdn,Advn) write(9,20) avn,sdn,Advn write(10,36) avn,sdn,Advn c c AVERAGE UNIT HEIGHT c call average(height,kk3,avh,sdh,Advh) write(9,21) avh,sdh,Advh write(10,37) avh,sdh,Advh c c FITTING LEAST SQUARES PLANE, LINE AND CIRCLE TO THE LOCAL HELIX ORIGINS c call fit(ox,oy,oz,kk-2,radc,rmsdc,rmsdl,r2) c c WRITING RESULTS IN TABULAR FORM c call table(fnm(k1),cb(k1),rb(k1),nrb(k1),re(k1),nre(k1), - avn,sdn,avh,sdh,avang,sdang,anglmax,radc,rmsdc,rmsdl,r2, - geom) if(geom.eq.'L') nlinear = nlinear+1 if(geom.eq.'C') ncurved = ncurved+1 if(geom.eq.'K') nkinked = nkinked+1 if(geom.eq.'*') nogeom = nogeom+1 close (4) endif 506 continue nhelices = nlinear+ncurved+nkinked+nogeom write(13,68) nlinear,ncurved,nkinked,nogeom,nhelices c C WRITING SOME TEXT TO TABLULAR FORM OF THE RESULTS. C write (13,66) write (13,54) write (13,55) write (13,56) write (13,57) write (13,58) write (13,59) write (13,60) write (13,61) write (13,62) write (13,63) write (13,64) write(13,69) C write(10,67) c C FIFTH QUESTION : DO YOU WISH TO USE ANOTHER PDB FILE? c 156 if((ans1.eq.'y').or.(ans1.eq.'Y')) then if((ans2.eq.'y').or.(ans2.eq.'Y')) then 153 write(6,51) write(7,51) read(5,'(a1)') ans5 write(7,'(1x,a1)') ans5 if ((ans5.ne.'n').and.(ans5.ne.'N').and.(ans5.ne.'y') - .and.(ans5.ne.'Y')) then goto 153 endif if ((ans5.eq.'y').or.(ans5.eq.'Y')) then go to 155 endif if ((ans5.eq.'n').or.(ans5.eq.'N')) then go to 154 endif endif endif 154 continue stop end c function dot (x,y) c c COMPUTES AND RETURNS THE DOT PRODUCT OF THE VECTORS X AND Y c dimension x(3),y(3) dot = 0.0 do 1000 i=1,3 dot = dot + x(i) * y(i) 1000 continue return end c c SUBROUTINE TO CALCULATE CROSS PRODUCTS OF VECTORS A and B. c subroutine cross(a,b,el,em,en) dimension a(3), b(3), c(3) c(1) = a(2) * b(3) - a(3) * b(2) c(2) = a(3) * b(1) - a(1) * b(3) c(3) = a(1) * b(2) - a(2) * b(1) fac1 = c(1) * c(1) + c(2) * c(2) + c(3) * c(3) fac1 = sqrtf(fac1) if(fac1.le.0.00001) go to 100 fac1 = 1.0/fac1 100 continue el = c(1) * fac1 em = c(2) * fac1 en = c(3) * fac1 return end c c SUBROUTINE TO CALCULATE AVERAGES c subroutine average(b,npts,av,sd,Adv) c CALCULATES AVERAGE, STANDARD DEVIATION AND MEAN ABSOLUTE DEVIATION dimension b(100) sum = 0.0 do 514 j = 1, npts sum = sum + b(j) 514 continue av= sum/float(npts) sum1=0.0 sum2=0.0 do 515 j=1,npts sum1=sum1+(b(j)-av)*(b(j)-av) sum2=sum2+abs(b(j)-av) 515 continue sum1=sum1/float(npts-1) sd=sqrtf(sum1) Adv=sum2/float(npts) return end c c SUBROUTINE TO FIT PLANE,CIRCLE AND LINE TO LOCAL HELIX ORIGINS subroutine fit(ox,oy,oz,np,radc,rmsdc,rmsdl,r2) c dimension ox(100),oy(100),oz(100) dimension rx(100),ry(100),rz(100) dimension rmp(100),rmc(100),rml(100) dimension matp(3,3),matc(3,3),matl(2,2),rotmt(3,3) dimension pmat(3,3),cmat(3,3),lmat(2,2) dimension v(3),ap(3),ac(3),al(2) c real matp,matc,matl,lmat double precision matp,pmat,matc,cmat,matl,lmat c 22 format(1x,'R.M.S. DEVIATION FROM BEST PLANE',8x,f8.3,2x,'ANG.') 23 format(1x,'L, M AND N OF THE BEST PLANE',11x,3(1x,f8.3)) 24 format(1x,'R.M.S. DEVIATION FROM BEST CIRCLE',7x,f8.3,2x,'ANG.') 25 format(1x,'RADIUS OF THE BEST CIRCLE',13x,f10.3,2x,'ANG.') 26 format(1x,'CENTRE OF THE BEST CIRCLE',13x,2(f10.3,1x)) 27 format(1x,'R.M.S. DEVIATION FROM BEST LINE',9x,f8.3,2x,'ANG.') 28 format(1x,'SLOPE OF THE BEST LINE',18x,f8.3) 29 format(1x,'INTERCEPT OF THE BEST LINE',14x,f8.3) 30 format(1x,'SQUARE OF LINEAR CORRELATION COEFFICIENT',f8.3) 31 format(1x,'FIT TO THE CIRCLE IS NOT GOOD') 32 format(1x,'FIT TO THE LINE IS NOT GOOD') 33 format(1x,'FIT TO THE PLANE IS NOT GOOD') 34 format(1x,'PERPENDICULAR DISTANCE OF BEST PLANE',4x,f8.3, 1 2x,'ANG.') c c FITTING LEAST SQUARES PLANE TO LOCAL HELIX ORIGINS c x2=0.0000000 y2=0.0000000 z2=0.0000000 xy=0.0000000 xz=0.0000000 yz=0.0000000 x=0.0000000 y=0.0000000 z=0.0000000 c do 516 j=1,np x2=ox(j)*ox(j)+x2 y2=oy(j)*oy(j)+y2 z2=oz(j)*oz(j)+z2 xy=ox(j)*oy(j)+xy xz=ox(j)*oz(j)+xz yz=oy(j)*oz(j)+yz x=ox(j)+x y=oy(j)+y z=oz(j)+z 516 continue c do 100 l=1,3 do 101 m=1,3 matp(l,m)=0.000000 101 continue 100 continue matp(1,1)=x2+matp(1,1) matp(1,2)=xy+matp(1,2) matp(1,3)=xz+matp(1,3) matp(2,1)=xy+matp(2,1) matp(2,2)=y2+matp(2,2) matp(2,3)=yz+matp(2,3) matp(3,1)=xz+matp(3,1) matp(3,2)=yz+matp(3,2) matp(3,3)=z2+matp(3,3) call matinv(matp,pmat) c bp1=x bp2=y bp3=z c do 517 j=1,3 ap(j)=pmat(j,1)*bp1+pmat(j,2)*bp2+pmat(j,3)*bp3 517 continue c sump=0.0 do 518 j=1,np rmp(j)=ap(1)*ox(j)+ap(2)*oy(j)+ap(3)*oz(j)-1 sump=sump+rmp(j)*rmp(j) 518 continue c if(sump.lt.0.0) then write(3,33) else rmsdp=(sqrtf(sump/np)) write(3,22) rmsdp endif c c CONVERTING TO THE NORMAL FORM OF THE PLANE c deno is denominator sqrt(ap(1)*ap(1)+ap(2)*ap(2)+ap(3)*ap(3)) c pp is the distance of plane from XY plane c deno=sqrtf(ap(1)*ap(1)+ap(2)*ap(2)+ap(3)*ap(3)) pp=1.0/deno pl=ap(1)*pp pm=ap(2)*pp pn=ap(3)*pp write(3,23) pl,pm,pn write(3,34) pp c c TO REORIENT THE POINTS SO THAT THE BEST FIT PLANE COINCIDES WITH c (X-Y) PLANE c The l,m,n of the normal to (X-Y) plane are (0,0,1) c Angle between the normal to best fit plane & the (X-Y) plane. c c c VECTOR NORMAL TO THE L,M,N OF THE BEST FIT & (X-Y) PLANE c do 102 j=1,3 v(j) = 0.0000000 102 continue v(1)=pm+v(1) v(2)=-pl+v(2) v(3)=0.0+v(3) vmag=sqrtf(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) v(1)=v(1)/vmag v(2)=v(2)/vmag v(3)=v(3)/vmag c THE ROTATION MATRIX REQUIRED TO MAKE l,m,n as 0,0,1 c costheta=pn costheta1=1.0-pn theta=acos(pn) sintheta=sin(theta) a1=v(1)*sintheta a2=v(2)*sintheta a3=v(3)*sintheta b1=v(2)*v(3)*costheta1 b2=v(3)*v(1)*costheta1 b3=v(1)*v(2)*costheta1 do 103 j=1,3 do 104 k=1,3 rotmt(j,k)=0.0000000 104 continue 103 continue rotmt(1,1)=costheta+v(1)*v(1)*costheta1+rotmt(1,1) rotmt(2,2)=costheta+v(2)*v(2)*costheta1+rotmt(2,2) rotmt(3,3)=costheta+v(3)*v(3)*costheta1+rotmt(3,3) rotmt(1,2)=b3-a3+rotmt(1,2) rotmt(2,1)=b3+a3+rotmt(2,1) rotmt(3,1)=b2-a2+rotmt(3,1) rotmt(1,3)=b2+a2+rotmt(1,3) rotmt(2,3)=b1-a1+rotmt(2,3) rotmt(3,2)=b1+a1+rotmt(3,2) c c ROTATING THE POINTS SO THAT THEY LIE IN (X-Y) PLANE c do 105 j=1,np rx(j) = 0.00000000 ry(j) = 0.00000000 rz(j) = 0.00000000 105 continue do 519 j=1,np rx(j)=rotmt(1,1)*ox(j)+rotmt(1,2)*oy(j)+rotmt(1,3)*oz(j)+rx(j) ry(j)=rotmt(2,1)*ox(j)+rotmt(2,2)*oy(j)+rotmt(2,3)*oz(j)+ry(j) rz(j)=rotmt(3,1)*ox(j)+rotmt(3,2)*oy(j)+rotmt(3,3)*oz(j)+rz(j) if(abs(rx(j)).le.1.0E-04) rx(j)=0.0000 if(abs(ry(j)).le.1.0E-04) ry(j)=0.0000 if(abs(rz(j)).le.1.0E-04) rz(j)=0.0000 519 continue c c TRANSLATING THE POINTS SO THAT THEY LIE IN X-Y PLANE OF THE FORM c ax+by=0 c do 520 j=1,np rx(j)=rx(j) ry(j)=ry(j) rz(j)=rz(j)-pp 520 continue c c FITTING CIRCLE TO THESE REORIENTED POINTS c c Equation used is (x+a)^2+(y+b)^2=r^2 c center of circle (-a,-b) ,c=r^2-a^2-b^2 c delta=x^2+y^2+2ax+2by-c c circle finds out radius of curvature and coordinates of center c from a set of given points using LEAST SQUARE FIT method. c x2=0.0000000 y2=0.0000000 x3=0.0000000 y3=0.0000000 xy2=0.0000000 x2y=0.0000000 xy=0.0000000 x=0.0000000 y=0.0000000 c do 521 j=1,np x3=rx(j)*rx(j)*rx(j)+x3 y3=ry(j)*ry(j)*ry(j)+y3 x2y=rx(j)*rx(j)*ry(j)+x2y xy2=rx(j)*ry(j)*ry(j)+xy2 x2=rx(j)*rx(j)+x2 y2=ry(j)*ry(j)+y2 xy=rx(j)*ry(j)+xy x=rx(j)+x y=ry(j)+y 521 continue c do 106 j=1,3 do 107 k=1,3 matc(j,k)=0.00000000 107 continue 106 continue matc(1,1)=x2+matc(1,1) matc(1,2)=xy+matc(1,2) matc(1,3)=-x+matc(1,3) matc(2,1)=xy+matc(2,1) matc(2,2)=y2+matc(2,2) matc(2,3)=-y+matc(2,3) matc(3,1)=x+matc(3,1) matc(3,2)=y+matc(3,2) matc(3,3)=-np+matc(3,3) c call matinv(matc,cmat) c bc1=-x3-xy2 bc2=-x2y-y3 bc3=-x2-y2 do 108 l=1,3 ac(l) = 0.00000000 108 continue do 522 j=1,3 ac(j)=cmat(j,1)*bc1+cmat(j,2)*bc2+cmat(j,3)*bc3+ac(j) 522 continue ac(1)=ac(1)/2.0 ac(2)=ac(2)/2.0 ac(3)=ac(3) radcsq=ac(1)*ac(1)+ac(2)*ac(2)+ac(3) if (radcsq.gt.0.0) then radc=sqrtf(radcsq) write(3,25) radc write(10,25) radc write(3,26) -ac(1), -ac(2) else radc = 0.0 write(10,31) endif c sumc=0.00000000 c c dev=sqrt((x+a)^2+(y+b)^2)-r c rms dev=sqrt(dev^2/n) c do 523 j=1,np rem=(rx(j)+ac(1))*(rx(j)+ac(1))+(ry(j)+ac(2))*(ry(j)+ac(2)) rmc(j)=sqrtf(rem)-radc sumc=sumc+rmc(j)*rmc(j) 523 continue if(sumc.lt.0.0) then write(3,31) write(10,31) else rmsdc=sqrtf(sumc/np) write(3,24) rmsdc write(10,24) rmsdc endif c c FITTING LINE TO THE REORIENTED POINTS c EQUATION OF LINE Y=MX+C c c do 109 l=1,2 do 110 m=1,2 matl(l,m) = 0.00000000 110 continue 109 continue matl(1,1)=x2+matl(1,1) matl(1,2)=x+matl(1,2) matl(2,1)=x+matl(2,1) matl(2,2)=np+matl(2,2) call matinv2(matl,lmat) c bl1=xy bl2=y c do 111 k=1,2 al(k) = 0.0000000 111 continue do 524 j=1,2 al(j)=lmat(j,1)*bl1+lmat(j,2)*bl2+al(j) 524 continue c write(3,28) al(1) write(3,29) al(2) c suml=0.00000000 c do 525 j=1,np rml(j)=ry(j)-al(1)*rx(j)-al(2) suml=suml+rml(j)*rml(j) 525 continue if(suml.lt.0.0) then write(3,32) write(10,32) else rmsdl=(sqrtf(suml/np)) write(3,27) rmsdl write(10,27) rmsdl endif c c LINEAR CORRELATION COEFFICIENT c r=(n*(xy)-x*y)/sqrt((n*x2-x*x)*(n*y2-y*y)) c rnum=np*xy-x*y rdeno=(np*x2-x*x)*(np*y2-y*y) r=rnum/sqrtf(rdeno) r2=r*r write(3,30) r2 write(10,30) r2 return end C FUNCTION SQRTF(X) 1 FORMAT(2X,'SQUARE ROOT OF A NEGATIVE QUANTITY?',E16.8) IF(X.LT.0.0) GO TO 101 SQRTF = SQRT(X) RETURN 101 CONTINUE WRITE(*,1) X SQRTF = 0.0 RETURN END c c SUBROUTINE TO COMPUTE INVERSE OF A 3 X 3 MATRIX c subroutine matinv(h,s) c CALCULATES INVERSE OF A 3x3 MATRIX dimension h(3,3),s(3,3) double precision c11,c12,c13,c21,c22,c23,c31,c32,c33 double precision deth1,deth2,deth3,deth,s,h do 112 j = 1,3 do 113 k = 1,3 s(j,k) = 0.0000000 113 continue 112 continue c11=h(2,2)*h(3,3)-h(3,2)*h(2,3) c12=h(2,1)*h(3,3)-h(3,1)*h(2,3) c13=h(2,1)*h(3,2)-h(3,1)*h(2,2) c21=h(1,2)*h(3,3)-h(1,3)*h(3,2) c22=h(1,1)*h(3,3)-h(1,3)*h(3,1) c23=h(1,1)*h(3,2)-h(3,1)*h(1,2) c31=h(1,2)*h(2,3)-h(1,3)*h(2,2) c32=h(1,1)*h(2,3)-h(1,3)*h(2,1) c33=h(1,1)*h(2,2)-h(1,2)*h(2,1) deth1=h(1,1)*c11 deth2=-h(1,2)*c12 deth3=h(1,3)*c13 deth=deth1+deth2+deth3 if(deth.ne.0.0) then s(1,1)=c11/deth+s(1,1) s(1,2)=-c21/deth+s(1,2) s(1,3)=c31/deth+s(1,3) s(2,1)=-c12/deth+s(2,1) s(2,2)=c22/deth+s(2,2) s(2,3)=-c32/deth+s(2,3) s(3,1)=c13/deth+s(3,1) s(3,2)=-c23/deth+s(3,2) s(3,3)=c33/deth+s(3,3) else write (*,*) 'MATRIX IS SINGULAR' endif return end c c SUBROUTINE TO COMPUTE INVERSE OF A 2 X 2 MATRIX c subroutine matinv2(p,q) C INVERSE OF A 2x2 MATRIX dimension p(2,2),q(2,2) double precision c11,c12,c21,c22 double precision detp1,detp2,detp,p,q do 114 j = 1,2 do 115 k = 1,2 q(j,k) = 0.00000000 115 continue 114 continue c11=p(2,2) c12=p(2,1) c21=p(1,2) c22=p(1,1) detp1=p(1,1)*c11 detp2=-p(1,2)*c12 detp=detp1+detp2 if(detp.ne.0.0) then q(1,1)=c11/detp+q(1,1) q(1,2)=-c21/detp+q(1,2) q(2,1)=-c12/detp+q(2,1) q(2,2)=c22/detp+q(2,2) else write(*,*) 'MATRIX IS SINGULAR' endif return end c c SUBROUTINE TO FIND LARGEST NUMBER subroutine max(x,no,xmax) c FINDS MAXIMUM NO. dimension x(100) xmax=x(1) do 525 j=1,no if (x(j).gt.xmax) xmax=x(j) 525 continue return end subroutine table(code,cbg,rbg,nrbg,red,nred,turnn,sdn, - height,sdh,bav,sdb,bam,rad,rmsc,rmsl,rsq,geom) c SUBROUTINE TO WRITE HELIX ANALYSIS RESULTS IN TABULAR FORM. c character*11 code character*3 rbg,red character*1 cbg,rb1,re1 character*1 geom 9 format(1x,a11,1x,a1,1x,i3,a1,' - ',i3,a1,1x,f4.2,1x, - f4.2,1x,f4.1,1x,f4.1,2x,i4,2x,2(f5.2,1x),f4.2,3x,a1) 10 format(16x,'Std. Dev.',2x,f4.2,1x,f4.1,1x,f4.1) if(rbg.eq.'ALA') rb1='A' if(rbg.eq.'CYS') rb1='C' if(rbg.eq.'ASP') rb1='D' if(rbg.eq.'GLU') rb1='E' if(rbg.eq.'PHE') rb1='F' if(rbg.eq.'GLY') rb1='G' if(rbg.eq.'HIS') rb1='H' if(rbg.eq.'ILE') rb1='I' if(rbg.eq.'LYS') rb1='K' if(rbg.eq.'LEU') rb1='L' if(rbg.eq.'MET') rb1='M' if(rbg.eq.'ASN') rb1='N' if(rbg.eq.'PRO') rb1='P' if(rbg.eq.'GLN') rb1='Q' if(rbg.eq.'ARG') rb1='R' if(rbg.eq.'SER') rb1='S' if(rbg.eq.'THR') rb1='T' if(rbg.eq.'VAL') rb1='V' if(rbg.eq.'TRP') rb1='W' if(rbg.eq.'TYR') rb1='Y' if(red.eq.'ALA') re1='A' if(red.eq.'CYS') re1='C' if(red.eq.'ASP') re1='D' if(red.eq.'GLU') re1='E' if(red.eq.'PHE') re1='F' if(red.eq.'GLY') re1='G' if(red.eq.'HIS') re1='H' if(red.eq.'ILE') re1='I' if(red.eq.'LYS') re1='K' if(red.eq.'LEU') re1='L' if(red.eq.'MET') re1='M' if(red.eq.'ASN') re1='N' if(red.eq.'PRO') re1='P' if(red.eq.'GLN') re1='Q' if(red.eq.'ARG') re1='R' if(red.eq.'SER') re1='S' if(red.eq.'THR') re1='T' if(red.eq.'VAL') re1='V' if(red.eq.'TRP') re1='W' if(red.eq.'TYR') re1='Y' if(bam.ge.20.0) geom='K' if (bam.lt.20.0) then if ((rmsc.gt.1.0).and.(rmsl.gt.1.0)) then geom = '*' else ratio=rmsl/rmsc if (ratio.gt.1.0) then geom='C' endif if (ratio.le.0.7) then if (rsq.ge.0.8) then geom='L' else geom='*' endif endif if ((ratio.gt.0.7).and.(ratio.le.1.0)) then if(rsq.le.0.5) then geom='C' endif if(rsq.ge.0.8) then geom='L' endif if ((rsq.gt.0.5).and.(rsq.lt.0.8)) then geom = '*' endif endif endif endif if ((rad-int(rad)).lt.0.5) then irad=int(rad) else irad=int(rad+0.5) endif write(13,9) code,cbg,nrbg,rb1,nred,re1,turnn,height,bav, - bam,irad,rmsc, rmsl,rsq,geom write(13,10) sdn, sdh, sdb return end