Updating documentation #1
Conversation
|
Hello @fiona-naughton! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:
Comment last updated at 2024-11-01 02:00:39 UTC |
Welcome to Codecov 🎉Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests. Thanks for integrating Codecov - We've got you covered ☂️ |
|
@lilyminium as someone who has previously dealt with helanal, would you be able to have a quick look through and make sure there's no egregious errors etc? |
There was a problem hiding this comment.
Sorry, partial review -- I'm not great at geometry at the best of times and I need some time to actually focus on this instead of trying to dip in and out while doing other work. And my brain is a bit dead at the moment :( so sorry if my comments are a bit terse too
| * - ``local_origins``: the projected centre of each 4-atom window, in line | ||
| with atom :math:`c_{i+1}`. | ||
|
|
||
| Calculated as the approximate intersection of :math:`\mathbf{D_1}` and |
There was a problem hiding this comment.
The "approximate" and "assuming ~even spacing of atoms" is confusing me here since I'm not sure what that means. Looking at the code it has some convoluted logic about cylinders and radii -- does this shake out to be the intersection of the D1/D2 vectors? I might have to come back when my brain is less fuzzy.
There was a problem hiding this comment.
Yeah, no worries - this is one of the bits I had to spend a while wrapping my head around and trying to figure out how to phrase.
To just walk through my understanding here, if it helps (or in case you notice anything wrong with my logic 😅 ) . Apologies in advance for the terrible formatting:
-
By 'even spacing' I'm meaning the atoms are spaced evenly along a helical path, so that:
a) all B (vectors between neighboring atoms) have the same length
b) looking down the ‘local axis’ of a window, the atoms are evenly spaced around a circle (with radius r and origin O, which is what we’re trying to find)
c) there is the same ‘rise’ component (along the local axis) to each B, so D1 = B2-B1, D2=B3-B2 both have ‘rise’ 0, i.e. they’re both in the plane when looking down the helix axis. -
A vector D1 = B2-B1 is the diagonal of a parallelogram formed by the vectors B2 and B1. (Drawing it as this diagonal, it points to atom c1. In basically all the diagrams, we slide it to point out from c1 mostly because it makes it less messy that way).
By properties of parallelograms, the diagonals bisect the angles if it’s a rhombus – which assuming (1a), is the case here. I’m not actually sure how relevant this fact specifically is but the original code made a deal of it. -
More usefully, it seems apparent (there's probably a rigorous proof but I've not done it) that when we consider the view looking down the helix axis so the atoms are evenly spaced around a circle (1b), the Ds lie along ‘radial lines’. This wouldn’t be the case if say c1 and c2 were much closer than c1 and c0 (see Fig 1. below)
Thus the radius line R1 from the ‘local origin’ O to c1 is parallel to D1, and one R2 from O to c2 is parallel to D2. R1 and R2 both have some length R. D1 and D2 are both in plane (1b), so the angle between R1 and R2 is the angle theta between D1 and D2, which we have calculated from the cross product. -
Another property of parallelograms is the two diagonals bisect each other. D1 is one diagonal of the c0/c1/c2 parallelogram, and a line from c0 to c2 is the other. So if the diagonals intersect at X, the length of X -> c1 is half the length of D1. (Fig 2 below)
-
Again assuming (1a), c0/c1/c2 form a rhombus, so the diagonals intersect at right angles – ie. X-c2 and D1 (or R1) are perpendicular, and (O, c2, X) form a right triangle with hypotenuse O-c2 = R2 and angle c2-O-X = theta. Thus, OX has length Rcos(theta).
-
The radius vector R1 (O-c1) can be broken into OX and X-c1, so it’s length = R can also be calculated as Rcos(theta) + |D1|/2. Rearranging to find R gives |D1|/(2-2cos(theta)).
The code substitutes |D1| -> sqrt(|D1|*|D2|) here; best I can tell this is a concession to the fact things aren't perfectly even so some of the difference in spacing will be accounted for; but I couldn't figure out if there was a rigorous reason for this substitution (or it's even possible this is actually a more accurate derivation that avoids some of the assumptions I'm making?) -
The location of O is then calculated by starting at c2, and moving 'in' by R1 (ie moving along the direction of D1 by a distance of (minus) R, that is O = c2 – R D1/|D1|, with R calculated as above
So this all seems to make sense when we have the ‘even spacing assumption’ - D1, D2 will intersect at the origin and so on. But when we don’t have exact spacing I’m not sure exactly what is calculating – only that, given all the assumptions, I doesn't seem like it would be as straightforward as "the intersection of D1/D2", though it should be close enough most of the time – hence my cautious using of approximately.
I'm trusting that there's a reason for these what seem like very convoluted calculations. I couldn't find equations in the referenced papers, so I assume this comes from the original fortran code version of helanal, which doesn't appear to exist online anymore...
There was a problem hiding this comment.
Here's helanal.f, thank to the Internet Archive: https://web.archive.org/web/20090226192455/http://www.ccrnp.ncifcrf.gov/users/kumarsan/HELANAL/helanal.f
There was a problem hiding this comment.
(Maybe we can link the archived version in the docs?)
There was a problem hiding this comment.
And for posterity:
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
endThere was a problem hiding this comment.
Ah, thanks!! Ill have a look through that at some stage (though I'm no expert in fortran). Good idea to to add a link, I might do that now.
There was a problem hiding this comment.
Thanks so much for walking through that with me @fiona-naughton, your diagrams here were super helpful here too. I think I get it, although I do agree that if you have to do so much assuming that atoms are equally spaced, you could probably have fit a circle to the points. That being said the bottom of figure 1b also makes me think that the approximation still works relatively well -- if your helix is contorted like that, there will be an impact on the origin too.
I'll leave this unresolved in case you wanted a reminder to look at the fortran, but the explanation really cleared things up for me.
| (or :math:`\mathbf{D_2}`) for each atom :math:`c_{i+1}` (or | ||
| :math:`c_{i+2}`). | ||
|
|
||
| Assuming ~even spacing of the atoms, this vector will bisect the angle |
There was a problem hiding this comment.
I'm not sure the "Assuming ~even spacing" would be helpful wording to me as a user trying to figure out what this is, although I have a serious case of brain-fuzz at the moment. Or maybe if there was an additional footnote explaining why it's only approximately bisecting it because the B_x vectors aren't normalised would be helpful.
There was a problem hiding this comment.
Very fair - I've sort of addressed this above.
As said above, I don't think the fact these are bisectors is actually specifically relevant. I'd kept the references because the code is still strongly linked to it (e.g. with the variable for the D vectors being called bisectors) and I didn't want to mess with those, and thought it might help people visualise what's going on - but with all my cautious it's only approximate it may be adding confusion more than helping.
Perhaps for now, I could simplify to only saying they'll point from the local origin to the atom, and maybe in the future add some notes later expanding on the exact calculation, something like the above?
There was a problem hiding this comment.
Maybe saying both would be clearest? So something like "this vector points from the local origin to the atom. Assuming ~even spacing of the atoms, ..."
some notes later expanding on the exact calculation
This would be amazing so we don't keep revisiting this 😅 And for any user who really really wants to know what's happening without reading the paper or the code...
|
By the way, I really like the table and visual depiction of each property!! Huge improvement over the previous hacky helanal docs. They are just making me actually think about each property though, which is why I've been slow to review -- sorry about that. |
Co-authored-by: Lily Wang <31115101+lilyminium@users.noreply.github.com>
Co-authored-by: Lily Wang <31115101+lilyminium@users.noreply.github.com>
Co-authored-by: Lily Wang <31115101+lilyminium@users.noreply.github.com>
Co-authored-by: Lily Wang <31115101+lilyminium@users.noreply.github.com>
|
Thanks @lilyminium for the review so far! Apologies for throwing this on you, I do definitely sympathise it's a pain to wrap one's head around |
|
Everything else LGTM -- sorry, maybe I was just running sick that day when I was confused about the properties! Something went wrong with conda for the unix tests -- I'll restart CI and if that doesn't work we might need to fix the mamba action |
|
Thanks @lilyminium! I restarted the tests again and everything seems to be passing now |
|
Thanks @fiona-naughton! |
In the process of preparing helanal for an MDAKit I got a little lost exactly what the code was doing, and ended up re-writing parts of the documentation to (hopefully) be clearer, including some new diagrams.