Laboratory for Computation and Visualization in Mathematics and Mechanics (LCVM^{2})

Institut des mathématiques Bernoulli (IMRI) - Section de mathématiques (SMA) Faculté des Sciences de Base (SB)

Computational physics project Final report

C_{α} circles

Christian Schlatter^{(1)}

11, avenue de Mont-Goulin CH-1008 Prilly Switzerland

July, 2002

Made with L^{A}TEX 2ε

(1) student in physics, christian.schlatter@epfl.ch

Contents

1 Introduction 1

2 Theoretical background 2 2.1 Proteins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 Global curvature, thickness and f-ratio . . . . . . . . . . . . . . . . . . . . . . . 4 2.3 The Protein Data Bank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 The algorithm 7 3.1 EXTRACTOR (Release 1.31) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.2 CALCULATOR (Release 1.62) . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.2.1 Fitting circles through three points . . . . . . . . . . . . . . . . . . . . . 9 3.2.2 Fitting spheres through four points . . . . . . . . . . . . . . . . . . . . . 9 3.2.3 Data saving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.3 Sourcecode details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4 Discussion 11 4.1 Program testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4.2.1 Global and local radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.2.2 Dislocation of circles and spheres . . . . . . . . . . . . . . . . . . . . . . 14 4.2.3 Optimality of protein shaping . . . . . . . . . . . . . . . . . . . . . . . . 15 4.2.4 Thickness of proteins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4.3 Runtime analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4.4.1 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

A Sample data files 22 A.1 C_{α} input data file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 A.2 C_{α} output data file (circles) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 A.3 C_{α} output data file (spheres) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

B Sourcecode of the algorithm 25 B.1 EXTRACTOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 B.2 CALCULATOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Projet de simulation numérique 2002  i  Christian Schlatter

List of Figures

2.1 The backbone of a protein. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 The bovine pancreatic trypsin inhibitor (BPTI) protein. . . . . . . . . . . . . . 3

4.1 Local and global circular radius of curvature. . . . . . . . . . . . . . . . . . . . . 12 4.2 Local and global spherical radius of curvature. . . . . . . . . . . . . . . . . . . . 12 4.3 Distribution of the smallest global circular radii of curvature. . . . . . . . . . . . 13 4.4 Distribution of the smallest global spherical radii of curvature. . . . . . . . . . . 13 4.5 Behavior of the nine smallest circle radii for each C_{α} atom. . . . . . . . . . . . . 14 4.6 Behavior of the nine smallest sphere radii for each C_{α} atom. . . . . . . . . . . . 15 4.7 Dislocation number ✠ for circles fitted through BPTI. . . . . . . . . . . . . . . . 16 4.8 Dislocation number ✠ for spheres fitted through BPTI. . . . . . . . . . . . . . . 16 4.9 Packing of α-helix in the chicken cartilage matrix protein. . . . . . . . . . . . . 17 4.10 Optimality of packing (circular). . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.11 Optimality of packing (spherical). . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.12 Runtime analysis for sphere calculation. . . . . . . . . . . . . . . . . . . . . . . 19

Projet de simulation numérique 2002  ii  Christian Schlatter

List of Tables

2.1 Description of the bovine pancreatic trypsin inhibitor in the Protein Data Bank (PDB). Parts of the header and the description of the backbone are shown. … 6

4.1 Circular and spherical thickness △ for some processed proteins. Note that the number of atoms corresponds to the number of extracted C_{α} atoms. The thickness doesn’t depend on the length of the proteins backbone. . . . . . . . . . . . 17

Projet de simulation numérique 2002  iii  Christian Schlatter

Chapter 1

Introduction

This computational physics project makes part of the two-semestrial lecture ”Simulation numérique de systèmes physiques” taught by professor Alfredo Pasquarello to students in physics. I had the pleasure to get a project at the Bernoulli institute from professor John Maddocks (Chaire d’analyse appliquée).

Proteins show a high packing density of residues. This packing has not been completely characterized^{(1)} yet. The molecule packing effects on functional dynamics. Very recently, Bagci, Jernigan and Bahar studied the packing of amino acids of different sizes and shapes [10] and [11]. Bagci and al. found that the residue of coarse-grained^{(2)} protein structures show a quite uniform density. Too, their analysis of the packing geometry leads to the conclusion, that the idealization of the molecule as homopolymeric chain (formed of equivalent residues) packs as spheres^{(3)} do.

Other recent work [7], [8] suggested that the packing of 3-dimensional structures is related to the concept of global curvature [9]. Using this approach, Maritan and al. [7] studied the packing of the most densely packed core regions, which are found to be those where the backbone of the protein forms special tertiary structures like helices. Such compact strings are found to be optimally shaped.

The algorithm developed within this project serves to analyze the presence and organization of several 3-dimensional protein configurations (determined by experiment) using the concept of global curvature. The whole backbone (or only sequences forming helices) of proteins can be treated using the coarse grained idealization of geometrical structure. The algorithm permits to study the behavior of the global curvature based on the fitting of circles and spheres through three or four residues.

The next chapter lights up some facts about proteins, the concept of the global radius and the Protein Data Bank. In chapter 3 the organization and content of the code is presented. The last chapter gives the results of the calculations and opens the interpretation phase. The appendix lists the program source code and some sample data files.

(1) Analyzing the genome could be simplified if packing metrics would be known.

(2) One space point per residue.

(3)

Closest packing of spheres corresponds to a fcc-lattice, as suggested already by Kepler and definitely

proven in 1998 only! [12].

Projet de simulation numérique 2002  1  Christian Schlatter

Theoretical background

2.1 Proteins

A protein molecule is a linear polymer that is made from a long chain of amino acids, each linked to its neighbor through a covalent peptide bond. The repeating sequence of atoms along the chain is referred to as the polypeptide backbone (see figure 2.1).

Different^{(1)} amino acid side chains are attached to one of the central carbon atoms (called the C_{α} atom) of each backbone amino acid. The side chains are these portions of the amino acids that are not involved in making a peptide bond and which gave each amino acid its unique properties.

Figure 2.1: The backbone of a protein. Four residues (blue blocks) are shown. From the left to the right a methionine, an aspartic acid, a leucine and a tyrosine are bonded together to form the backbone of a protein. The two ends of a polypeptide chain are always different. The beginning group consists of an NH_{2} - the end group of COOH. Proteins are always read from the N to the C direction.

(1) Out of 20 kinds.

Projet de simulation numérique 2002  2  Christian Schlatter

Proteins are extremely varied and of complicated shape. The most simple protein has about 30 amino acids, the most complicated about 10’000. Its backbones and side chains are folded into a conformation of lowest energy. The over all conformation of each protein is unique, but locally, two regular folding patterns are often found: the α-helix (the backbone turns around itself to form a cylinder) and the β-sheet (several parallel aligned backbones sequences).

The bovine pancreatic trypsin inhibitor (BPTI) is a protein which has been largely studied in chemistry, biology, physics and numerical experiments. BPTI is present in many tissues of the human body where it inhibits the proteolytic action of the enzyme trypsin. BPTI is a monomeric polypeptide containing 58 amino acid residues. Figure 2.2 shows BPTI and is colored by secondary structure.

Figure 2.2: Bovine pancreatic trypsin inhibitor (BPTI) protein. The coil (in grey) is colored by secondary structure. In red the two α-helices, in blue the β-sheets. The three disulfide bonds are shown in yellow.

The study of the spatial configuration of the backbone at a coarse grained level can be done by considering the coordinates of only one repetitive backbone amino acid atom, as for example the C_{α} carbon atom (hence the name of this project.

2.2 Global curvature, thickness and f-ratio

The n atoms of a chain-like molecule (as the backbone of a protein) can be represented as discrete curve C_{n} by a list of distinct points {q_{1}, .., q_{n}} in three-dimensional space. The points are given by a repetitive atom type figuring in the backbone (coarse-grain level consideration).

Let q_{i}, q_{j}, q_{k} ∈ ℝ^{3}, with i, j, k ∈ {1, .., n} be three such non-collinear points. An unique circle passing through each of these points can be construed whose radius r(q_{i}, q_{j}, q_{k}) is given by

r(q_{i}, q_{j}, q_{k})^{2} = \frac{−D_{ij}^{2}D_{ik}^{2}D_{jk}^{2}}{D_{ij}^{4} + D_{ik}^{4} + D_{jk}^{4} − 2(D_{ij}^{2}D_{ik}^{2} + D_{ik}^{2}D_{jk}^{2} + D_{ij}^{2}D_{jk}^{2})}  (2.1)

where D_{ij} = |q_{i} − q_{j}| is the euclidian distance between points i and j.

The local radius of curvature ρ at point q_{i} of curve C_{n} corresponds to the circle which locally best approximates the curve. In the discrete case it is simply defined by

ρ(q_{i}) = r(q_{i−1}, q_{i}, q_{i+1}), 1 < i < n  (2.2)

The global radius of curvature ρ_{G} at point q_{i} of C_{n} generalizes the concept of local curvature. It is the radius of the smallest circle containing q_{i} and any two other distinct points q_{j} and q_{k} i.e.

ρ_{G}(q_{i}) = min_{\begin{matrix}1⩽j,k⩽n \\ i≠j≠k≠i \end{matrix}} r(q_{i}, q_{j}, q_{k})  (2.3)

The global radius of curvature ρ_{G} takes into account both, the local bending of the curve and the non-local proximity of different parts of the curve and is therefore predestined for a characterization of its packing.

Finally the thickness Δ[C_{n}] associated to C_{n}, i.e.

Δ[C_{n}] = min_{1⩽i⩽n} ρ_{G}(q_{i})  (2.4)

denotes the maximum radius if a uniform tube with the curve passing through its axis, beyond which the tube either ceases to be smooth, owing to tight local bends, or it self-intersects.

To characterize optimal strings (such as helices) in terms of the sphere packing problem, Maritan and al. [7] introduced the f-ratio. Optimal helices show equality of the local radius of curvature ρ and the radius that is associated with a suitable triplet of non-consecutive points lying in two successive turns of the helix. Let us call the non-local radius of this special non-local triplet ρ_{N}, defined by

ρ_{N}(q_{i}) = min_{\begin{matrix}1⩽i<j<k⩽n \\ (k−j)+(j−i)>2 \end{matrix}} r(q_{i}, q_{j}, q_{k})  (2.5)

The local f-ratio f(q_{i}) is then defined by

f(q_{i}) = \frac{ρ_{N}}{ρ}  (2.6)

For f(q_{i}) > 1 the global radius of curvature ρ_{G}(q_{i}) corresponds to the local radius of curvature ρ(q_{i}) (local regime), for f(q_{i}) < 1 the global radius of curvature is strictly lower than the local radius, and the thickness of the helix is determined basically by the distance between two consecutive helix turns. Optimal packing selects the very special helix with f(q_{i}) = 1.

The global f-ratio F [C_{n}], averaged over the whole curve (or an interval of the curve) measures then the optimality of its packing.

All the above notions can be extended to the case where spheres are fitted through four points of the backbone. Let q_{i}, q_{j}, q_{k}, q_{l} ∈ ℝ^{3}, with i, j, k, l ∈ {1, .., n} be four non-coplanar points. An unique sphere passing through each of these points can be construed whose radius r(q_{i}, q_{j}, q_{k}, q_{l}) is given by

r(q_{i}, q_{j}, q_{k}, q_{l})^{2} = \frac{− | \begin{matrix}0 & D_{ij}^{2} & D_{ik}^{2} & D_{il}^{2} \\ D_{ij}^{2} & 0 & D_{jk}^{2} & D_{jl}^{2} \\ D_{ik}^{2} & D_{jk}^{2} & 0 & D_{kl}^{2} \\ D_{il}^{2} & D_{jl}^{2} & D_{kl}^{2} & 0 \end{matrix} |}{2 | \begin{matrix}0 & D_{ij}^{2} & D_{ik}^{2} & D_{il}^{2} & 1 \\ D_{ij}^{2} & 0 & D_{jk}^{2} & D_{jl}^{2} & 1 \\ D_{ik}^{2} & D_{jk}^{2} & 0 & D_{jl}^{2} & 1 \\ D_{il}^{2} & D_{jkl}^{2} & D_{kl}^{2} & 0 & 1 \\ 1 & 1 & 1 & 1 & 0 \end{matrix} |}  (2.7)

where D_{ij} = |q_{i} − q_{j}| is the euclidian distance between points i and j.

The local radius of curvature ρ at point q_{i} and q_{i+1} is logically

ρ(q_{i}) = ρ(q_{i+1}) = r(q_{i−1}, q_{i}, q_{i+1}, q_{i+2}), 1 < i < (n − 1)  (2.8)

To be consistent, the global radius of curvature ρ_{G} at point q_{i} is

ρ_{G}(q_{i}) = min_{\begin{matrix}1⩽j,k⩽n \\ |u−v|≠0, \\ u,v∈\{i,j,k,l\} \end{matrix}} r(q_{i}, q_{j}, q_{k})  (2.9)

and the non-local radius ρ_{N} will then be

ρ_{N}(q_{i}) = min_{\begin{matrix}1⩽i<j<k<l⩽n \\ (l−k)+(k−j)+(j−i)>3 \end{matrix}} r(q_{i}, q_{j}, q_{k}, q_{l})  (2.10)

The thickness and the f-ratio are calculated accordingly.

2.3 The Protein Data Bank  6

2.3 The Protein Data Bank

The algorithm uses real protein structure data which has been determined by experiment. Protein data can been downloaded from the huge Protein Data Bank (PDB) database^{(2)}, which is hosting almost 20’000 biological macromolecular structures.

The Protein Data Bank uses a proprietary data format to save many information about protein. Each single protein is documented in so called PDB-files. Table 2.1 shows some parts of the file describing the BPTI protein.

A comprehensive description of the PDB file format can be found in [14]. The algorithm at issue conforms to release 2.2 of the PDB file format description.

(2) http://www.rcsb.org/pdb/

HEADER PROTEINASE INHIBITOR (TRYPSIN)  13−MAY−87 6PTI  6PTI 3
COMPND BOVINE PANCREATIC TRYPSIN INHIBITOR ( / BPTI$\$$ ,CRYSTAL FORM  6PTI 4

…

REMARK 2 RESOLUTION. 1 . 7 ANGSTROMS.  6PTI 36

…

HELIX 1 H1 PRO  2 GLU  7 5 ALL DONORS,ACCEPTORS INCLUDED 6PTI 58
HELIX 2 H2 SER  47 GLY  56 1 ALL DONORS,ACCEPTORS INCLUDED 6PTI 59

…

ATOM  1 N ARG  1  1 2 . 5 7 9 −6.640 1 1 . 3 7 6 1 . 0 0 2 4 . 0 4  6PTI 73
ATOM  2 CA ARG  1  1 2 . 1 3 9 −5.568 1 0 . 4 6 4 1 . 0 0 2 1 . 5 2  6PTI 74
ATOM  3 C ARG  1  1 1 . 3 0 4 −6.259 9 . 3 6 6 1 . 0 0 1 8 . 3 4  6PTI 75
ATOM  4 O ARG  1  1 1 . 7 2 4 −7.230 8 . 7 5 7 1 . 0 0 1 7 . 2 2  6PTI 76
ATOM  5 CB ARG  1  1 3 . 1 9 0 −4.672 9 . 8 7 8 1 . 0 0 2 2 . 8 5  6PTI 77
ATOM  6 CG ARG  1  1 3 . 8 2 0 −4.993 8 . 5 6 3 1 . 0 0 2 4 . 9 0  6PTI 78
ATOM  7 CD ARG  1  1 4 . 4 7 8 −3.884 7 . 8 5 2 1 . 0 0 2 6 . 6 2  6PTI 79
ATOM  8 NE ARG  1  1 5 . 3 4 7 −3.042 8 . 5 9 5 1 . 0 0 2 9 . 0 3  6PTI 80
ATOM  9 CZ ARG  1  1 6 . 5 5 5 −3.097 9 . 0 9 1 1 . 0 0 2 9 . 0 5  6PTI 81
ATOM 10 NH1 ARG  1  1 7 . 4 0 0 −4.129 9 . 0 0 5 1 . 0 0 2 9 . 3 6  6PTI 82
ATOM 11 NH2 ARG  1  1 7 . 0 4 2 −2.028 9 . 7 5 0 1 . 0 0 2 9 . 9 4  6PTI 83
ATOM 12 N PRO  2  1 0 . 0 9 9 −5.724 9 . 1 9 3 1 . 0 0 1 5 . 9 0  6PTI 84
ATOM 13 CA PRO  2  9 . 1 6 6 −6.247 8 . 1 8 1 1 . 0 0 1 4 . 6 3  6PTI 85
ATOM 14 C PRO  2  9 . 7 0 0 −6.139 6 . 7 7 4 1 . 0 0 1 2 . 1 4  6PTI 86
ATOM 15 O PRO  2  1 0 . 4 2 8 −5.188 6 . 3 8 9 1 . 0 0 1 1 . 6 4  6PTI 87
ATOM 16 CB PRO  2  7 . 9 1 5 −5.381 8 . 3 9 1 1 . 0 0 1 4 . 7 7  6PTI 88
ATOM 17 CG PRO  2  8 . 0 7 6 −4.722 9 . 7 0 8 1 . 0 0 1 5 . 3 2  6PTI 89
ATOM 18 CD PRO  2  9 . 5 7 8 −4.563 9 . 9 1 4 1 . 0 0 1 5 . 6 8  6PTI 90
ATOM 19 N ASP  3  9 . 3 2 4 −7.061 5 . 8 9 9 1 . 0 0 1 2 . 4 0  6PTI 91
ATOM 20 CA ASP  3  9 . 8 0 5 −7.020 4 . 4 9 9 1 . 0 0 1 2 . 6 7  6PTI 92
ATOM 21 C ASP  3  9 . 4 4 3 −5.726 3 . 7 8 1 1 . 0 0 9 . 9 5  6PTI 93
ATOM 22 O ASP  3  1 0 . 0 7 7 −5.313 2 . 8 0 2 1 . 0 0 9 . 2 5  6PTI 94
ATOM 23 CB ASP  3  9 . 1 9 3 −8.193 3 . 6 9 9 1 . 0 0 1 8 . 9 8  6PTI 95
ATOM 24 CG ASP  3  9 . 8 1 4 −9.527 4 . 0 3 1 1 . 0 0 2 6 . 0 4  6PTI 96
ATOM 25 OD1 ASP  3  9.236 −10.570 3 . 6 3 8 1 . 0 0 2 9 . 2 8  6PTI 97
ATOM 26 OD2 ASP  3  1 0 . 8 9 2 −9.510 4 . 6 9 1 1 . 0 0 2 6 . 0 7  6PTI 98

Table 2.1: Description of the bovine pancreatic trypsin inhibitor in the Protein Data Bank (PDB). Parts of the header and the description of the backbone are shown.

The algorithm

The present algorithm consists mainly of two programs, called EXTRACTOR and CAL-CULATOR. EXTRACTOR has to be used to generate the data files which will be read by CALCULATOR. Its input are data files downloaded from PDB. EXTRACTOR loads and analyzes these files and extracts the information we are interested in. The output are proprietary data files in the C_{α} source format (see appendix A, C_{α} input data file format). Input and output files are read and written from and to a subdirectory tree on the local disk. This tree structure can be created using the batch file DIRECTORIES.BAT.

CALCULATOR does only accept the C_{α} data files as input, which will be used for the ’circling’ and ’sphereing’ of the protein. The results of the calculation are stored into a subdirectory tree on disk in the proprietary C_{α} result format (see appendix A, C_{α} output data file format).

3.1 EXTRACTOR (Release 1.31)

When calling the EXTRACTOR the user has to provide the filename of the PDB file to be processed. Alternatively, several PDB files can be listed^{(1)} in a file called Files.ini. EXTRACTOR will then process file by file given by this list. All PDB files have to be stored in the subfolder PDB.

Next, the user has to provide the type of atoms to extract. For alpha-carbons (C_{α}) the corresponding identifier is CA (the same identifier as in PDB files is used).

The third step is to select if the whole backbone of the protein should get extracted or only atoms which are located in secondary structure α-helices. The latter case has been implemented for comparison purposes with existing results done by [7] (see figure 4.9 in section 4.2).

The PDB files are processed accordingly to the PDB format description guide [14] from December 20, 1996. PDB ASCII files are read and processed line by line by the EXTRACTOR. Each line starts with a record name followed by data fields of different lengths (see table 2.1). These fields are addressed by dividing the line string up in substrings of the corresponding length.

(1) One complete filename (including the file extension) by line, no blank lines.

Projet de simulation numérique 2002  7  Christian Schlatter

From PDB files the HEADER, REMARK 2, SETQRES, HELIX and ATOM records are interpreted. The HEADER entry contains classification details about the protein, the unique PDB ID, the date of acquisition of the PDB data and sometimes experimental techniques. The REMARK 2 entry may^{(2)} contain information about the highest resolution (in Å) that was used in building the model. The resolution is extracted for error calculation purposes, but the release at issue doesn’t make use of it (if present its written to the header of the output file). The SETQRES record lists the number of residues in the backbone. In case the whole backbone is extracted, this number is used for initialization of variables. By contrast, when extracting is limited to α-helices only, the HELIX record is used to determine the location and number of the concerned atoms. There may be several α-helices in the protein, both types - right- and left-handed - are considered. Finally, the atom coordinates and the corresponding chain identifier are provided from the ATOM record. If multiple chains are present, all are getting extracted. If the PDB file contains more than one model^{(3)}, only the first is taken into consideration.

Please be aware that PDB files may contain errors. If the number of residues given in the SETQRES or HELIX record are not correct, then the algorithm shrinks or expands the variables by re-initialization to extract the data correctly.

Rarely neighboring Cα atoms have been observed, what is obviously impossible. In these cases only the first atom is considered.

If for a given file there were no helices, the file is ignored.

The gathered data is then stored to disk in the folder source. For an example of the proprietary file format see appendix A. The output filename are the same as the input filenames. The file extension is changed to bck for whole backbone extractions and hel for helices only.

3.2 CALCULATOR (Release 1.62)

CALCULATOR follows the philosophy introduced with EXTRACTOR. It can be run in batch or single file mode. In batch mode the files have to be listed^{(4)} in the Source.ini file; in single mode the user provides the filename. CALCULATOR only accepts input files in the C_{α} input file format. Therefore, any PDB file has previously to be processed by EXTRACTOR.

The user has the choice between fitting circles only or circles and spheres. Fitting spheres does take much more time than fitting circles, see section 4.3 for more information.

(2) For diffraction experiments only.

(3) Using different experimental techniques or mutants.

(4) One complete filename (including the file extension) by line, no blank lines.

Projet de simulation numérique 2002  Christian Schlatter

3.2.1 Fitting circles through three points

The radius of a circle going through 3 points is calculated according to formula 2.1.

When three points are collinear one can’t draw a circle through them of course. The denominator in equation 2.1 vanishes for three collinear points, so this criterion is used to prevent collinearity.

The calculated radii are classed in the following three categories :

- local radii : three neighboring atoms i, j and k with i < j = i + 1 < k = j + 1 define the local radius for the atom j in the middle.

- partly adjacent radii : two of the three radii are neighbors.

- non-adjacent radii : no atoms are neighbors.

For each atom the smallest radii are ordered with increasing size, the number of memorized radii is limited by the flag memory.

When all combinations of circles have been calculated the global radius, the non-local radius, the f-number and the thickness are determined.

To know the distribution of the smallest circles the thickness of the protein has to be known (it is impossible to store all combinations of circles and doing statistics on the whole set of radii). To do that, a second call of the circling routine is initiated and the number of smallest radii with a resolution of 1 % is counted. The upper limit radius is given with the flag limit.

3.2.2 Fitting spheres through four points

The radius of a sphere going through 4 points is calculated according to formula 2.7.

Again, when the four points are coplanar, no sphere can be built with them. For four coplanar points the determinant in the denominator of equation 2.7 vanishes, what can be used as criterion.

The determinant is calculated using the partial pivoting Gaussian elimination scheme [1].

The calculated radii are classed in the following two categories :

- local radii : four neighboring atoms i, j, k and l with i < j = i + 1 < k = j + 1 < l = k + 1 define the local radius for the two atoms j and k in the middle.

- non-local radii : all other combinations of atoms.

Afterwards, the global radii, the thickness and the f-number are determined. A second run determines the distribution of the smallest spheres.

3.2.3 Data saving

The results are stored to the folder results. If the input file contained α-helices the data is saved to results/helix, backbone data to results/whole. The radii are always saved twice : one file is written to the subfolder human whose data is formatted to be easy readable by the user (see appendix A), another to the subfolder machine, which contains raw data assigned for further processing, e.g. for plotting with Gnuplot.

Please be aware that uninitialized values, incalculable values and infinity are noted by 0.0 (zero).

The content of these output files is determined by their file extensions, which are :

- cir : local and global circles for each atom, circular thickness, circular f-ratio.

- sph : local and global sphere for each atom, spherical thickness, spherical f-ratio.

- cdt : distribution of the circles.

- sdt : distribution of the spheres.

3.3 Sourcecode details

The algorithm has been developed under Microsoft Windows XP Professional using Compaq’s Visual Fortran 6.6 Development Environment (a fully Fortran 95 compiler, the former Digital compiler). The code at issue has been made compatible with SUN’s Fortran 90 compiler to run the program on the SUN workstation in the LCVM^{2} calculation laboratory.

The algorithm has been written as flexible as possible. Most variables are allocatables and the heaviest variables are pointers. The final code is a compromise between memory demands and execution speed.

The source code is structured using modules and subroutines. Adjustable parameters are grouped in a special module parameters just at the beginning of the files extractor.f90 and calculator.f90. The data type of the main variable (alphacarbons) which is used in both programs is also defined (data type carbons) in a single module var.

EXTRACTOR calls the subroutines helices (or backbones) which are located in the module reading to extract the coordinates of the atoms in the helix (or in the whole backbone).

CALCULATOR reads its input files using the subroutine import in the module loading. The circles (or spheres) are then calculated using the subroutines circles (or spheres) in module distance. The subroutine spheres makes use of the function determinant in module mathe to calculate the determinant of the expressions in formula 2.7.

For further details please refer to the source code of the algorithm in appendix B.

Discussion

4.1 Program testing

During development, the algorithm has been tested on several PDB files. All mutants or variants^{(1)} of the BPTI protein have been successfully processed. From April to August 2002 a bunch of as different as possible proteins have been processed. At present different output files are ready for interpretation. The plausibility of the resulting global radii has been confirmed by comparing to previous work, see section 4.2.3

4.2 Results

The following figures illustrate the data which is produced using CALCULATOR. All the calculations have been done using C_{α} backbone atoms only.

4.2.1 Global and local radius

Figures 4.1 and 4.2 show local and global radii of curvature along the backbone of the bovine pancreatic trypsin inhibitor protein (BPTI) consisting of 58 residues. PDB file 1PTI has been used. The local radius shows quite strong oscillations along the backbone. The circular thickness is achieved with residues 11, 35 and 36 and equals 2.46 Å. The spherical thickness is achieved with residues 49 to 52 and equals 2.84.

During the second run of CALCULATOR (i.e. when the thickness is known), the statistics of the smallest circles and spheres can be revealed. Figure 4.3 and 4.4 show the distribution regarding to the thickness for BPTI. The abscissa is given by the normalized global radius ρ_{G}/△. The behavior of the distribution of the radii depends on the number of atoms of the protein. A gap as it can be seen on the distribution of the smallest radii of the circular global radii (figure 4.3) is not always observed.

A more illustrative plot of the behavior of the global radii is given with figures 4.5 and 4.6. They show the smallest found circles and spheres along the backbone. The lowermost curve corresponds to the global radius of curvature. The polygons rely points which belongs to the same position in the hierarchy of increasing global radii.

(1) The Protein Data Bank returns at present 59 entries on a query for ’BPTI’.

Projet de simulation numérique 2002  11  Christian Schlatter

\underline{4.2 Results 12}

Figure 4.1: Local and global circular radius of curvature for BPTI. Blue squares show the local radii, the red triangles the global radii. The vertical axis show normalized values, i.e. all radii has been divided by the circular thickness of BPTI.

Figure 4.2: Local and global spherical radius of curvature for BPTI. Blue squares show the local radii, the red triangles the global radii. The vertical axis show normalized values, i.e. all radii has been divided by the spherical thickness of BPTI.

Figure 4.3: Distribution of the smallest global circular radii of curvature for BPTI. The resolution of the normalized global radius is 1 %. There are no other circles up to 5 % of the thickness.

Figure 4.4: Distribution of the smallest global spherical radii of curvature for BPTI. The resolution of the scaling factor is 1 %. Two spheres have radius in the interval {1.00, 1.01} times the thickness.

Figure 4.5: Behavior of the nine smallest circle radii for each C_{α} atom. The lowermost curve shows the global circular radii along the backbone. The curve just above shows the circles which have radius most next to those realizing the global radius, one could call them second global radii, and so on.

4.2.2 Dislocation of circles and spheres

An interesting feature available in the output files is the dislocation number ✠ which is given by

✠ = (k − j) + (j − i)  (4.1)

with i, j, k in ascending order. It is some kind of measure of how un-local the global radius is.

For ✠ = 2 the global circle is a local circle. Figure 4.7 shows circular ✠ for BPTI.

For spheres the dislocation number becomes

✠ = (l − k) + (k − j) + (j − i)  (4.2)

with i, j, k, l in ascending order. For ✠ = 3 the global sphere is a local sphere. Figure 4.8 shows spherical ✠ for BPTI.

Please note that circles are not that much scattered than spheres.

Figure 4.6: Behavior of the nine smallest sphere radii for each C_{α} atom for BPTI. The lowermost curve shows the global spherical radii along the backbone. Did you remark the strange gap between forth and fifth smallest circle for residue 26?.

4.2.3 Optimality of protein shaping

Maritan and al. calculated the local circular f-ratio for the helix in the the chicken cartilage protein. They used the PDB entry 1AQ5 and published a plot in their publication^{(2)}. A good way to test the developed algorithm was to reproduce this plot. Figure 4.9 shows the f-ratio calculated with CALCULATOR. The plots coincide and Maritan’s and al. global f-ratio for residues 14 to 42 stating the optimality of helices is also found : F = 1.01.

The f-ratio is now also calculated for the whole backbone of proteins. Figure 4.10 shows the local f-ratio calculated along BPTI. The concept of optimality measured by the f-ratio has been extended to spherical radii. Figure 4.11 shows the spherical f-ratio for BPTI. Please note the stability f ≈ 1.0 for the residues forming the α-helix.

4.2.4 Thickness of proteins

Table 4.1 shows the circular and spherical thickness for some processed proteins. The absence of the spherical thickness for larger proteins is simply due to the fact that it took too much time to calculate them (see next section).

(2) Figure 3 on page 289 in [7] Projet de simulation numérique 2002  Christian Schlatter

Figure 4.7: Dislocation number ✠ for circles fitted through BPTI. ✠ = 2 means that the global circle corresponds to the local circle.

Figure 4.8: Dislocation number ✠ for circles fitted through BPTI. ✠ = 3 means that the global sphere corresponds to the local sphere.

Figure 4.9: Packing of α-helix in the chicken cartilage matrix protein. The local f-ratio is is plotted for residues 14 to 42 which shape the helix.

Protein # atoms circular △ [Å] spherical △ [Å] Description
4INS  21  2.60  2.83 Insulin
1PPB  36  2.50  2.80 Thrombin
1PTI  58  2.46  2.84 Proteinase inhibitor
7HVP  97  2.39  2.78 HIV-1 Protease
1IGT  214  2.30  2.79 Antibodies
1FBB  224  2.46  2.70 Bacteriorhodopsin
2PLV  288  2.38  2.89 Poliovirus
1YVN  372  2.51  2.75 Actin
1HTB  374  2.31  2.77 Alcohol Dehydrogenase
1ASZ  423  2.44  2.72 Aminoacyl-tRNA Synthetases
1N2C  478  2.35  2.76 Nitrogenase
1PRH  522  2.33  - Cyclooxygenase
1ACC  665  2.37  - Anthrax Toxin
1JB0  740  2.30  - Photosystem I
1B7T  766  2.41  - Myosin
6GPB  828  2.34  - Gycogen Phosphorylase

Table 4.1: Circular and spherical thickness △ for some processed proteins. Note that the number of atoms corresponds to the number of extracted C_{α} atoms. The thickness doesn’t depend on the length of the proteins backbone.

Figure 4.10: Optimality of packing. Circular f-ratio for BPTI. Two helices are present: a right-handed 3_{10} helix from residue 3 to 7, a right-handed α-helix from residue 47 to 56.

Figure 4.11: Optimality of packing. Spherical f-ratio for BPTI. Two helices are present: a right-handed 3_{10} helix from residue 3 to 7, a right-handed α-helix from residue 47 to 56.

4.3 Runtime analysis

Large macromolecules require an astronomical number of radii to be calculated. Using the Fortran package StopWatch [2], the execution time for different proteins has been analyzed. The measurements have been done on a SUN workstation ?? running Solaris 8 ??

Figure 4.12 shows the number of hours necessary to process several bigger proteins. The time dependence t is linear in the number of radii to calculate, i.e.

t = A · \binom{n}{k}  (4.3)

where n is the number of extracted atoms (i.e. the number of residues in the protein), k = 3 for circles or k = 4 for spheres.

For the 1N2C file (a complex of nitrogenase proteins containing 478 C_{α} atoms) the calculation of the spheres takes up to twelve (!) days, whilst for its circles, only 30 seconds were necessary. On LCVM^{2}’s computer the coefficient A equals 4.95·10^{−4} for spheres if t is measured in seconds.

Please bear in mind that the Sun system is not very fast. For comparison purposes some files have also been processed on an AMD Athlon running at 1.2 GHz under Microsoft Windows XP Professional, where all the jobs were fulfilled more than 6 times faster.

Figure 4.12: Runtime analysis for sphere calculation. Execution time is linear in the number of combinations of spheres to calculate. During the twelve days necessary to process PDB file 1N2C, 2’148’006’525 spheres were calculated.

4.4 Conclusion

The algorithm has been tested and data for several proteins has been produced. Further data can easily be produced and the gathered data is open for interpretation.

The present implementation of CALCULATOR works fine and delivers results for smaller proteins within short time. For longer proteins as for example DNA, a precalculation to select neighboring atoms should be used.

4.4.1 Acknowledgements

I would like to thank professor John H. Maddocks for accepting me as project student and for taking counsel with me during the project period. A big thank goes to Jana Smutny, doctoral candidate at LCVM^{2}. She was providing me with the formulae 2.1 and 2.7 [6] and assisted my test runs on LCVM^{2} workstations helping me to make my algorithms compatible with the Sun Fortran 90 compiler. Thanks also to Mister Caussignac, for the help solving computer problems, hard disk crashes and launching remote jobs on Sun OS in the computer laboratory at the institute.

Prilly, August 2002

Christian Schlatter

References

[1] Jacques Rappaz and Marco Picasso; Introduction à l’analyse numérique, Presses polytechniques et universitaires romandes, Lausanne 1998.

[2] William F. Mitchell, StopWatch: A Module for Portable Measurement of Execution Time, Fortran Journal Vol. 9, No. 1, 1997.

[3] Bruce Alberts, Dennis Bray et al.; Essential Cell Biology, An introduction to the molecular biology of the cell; Garland Publishing, New York 1998.

[4] J. C. Adams, W. S. Brainerd et al.; The Fortran 95 Handbook, MIT Press, Cambridge 1997.

[5] M. Posternak, Alfonso Baldereschi; Introduction aux moyens informatiques, EPFL 2000.

[6] Jana Smutny, Radii and centres of spheres and circles, EPFL 2001, internal document.

[7] Amos Maritan, Cristian Micheletti, Antonio Trovato and Jayanth R. Banavar; Optimal shapes of compact strings, Nature, vol. 406, pp. 287-290, London, july 2000.

[8] Andrzej Stasiak and John H. Maddocks; Best packing in proteins and DNA, Nature, vol. 406, pp 251-253, London, july 2000.

[9] Oscar Gonzalez and John H. Maddocks; Global curvature, thickness, and the ideal shape of knots, Proc. Natl. Acad. Sci. USA, vol. 96, pp. 4769-4773, april 1999.

[10] Zerrin Bagci, Robert L. Jernigan and Ivet Bahar, Residue packing in proteins: Uniform distribution on a coarse-grained scale, Journal of chemical physics, vol. 116, nr. 5, pp 2269-2276, february 2002.

[11] Zerrin Bagci, Robert L. Jernigan and Ivet Bahar, Residue coordination in proteins conforms to the closest packing of spheres, polymer, vol. 43, pp. 451-459, 2002.

[12] N. J. A. Sloane, Kepler’s concecture confirmed, Nature, vol. 395, pp. 435-436, London, 1998.

[13] H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I.N. Shindyalov, P.E. Bourne; The Protein Data Bank. Nucleic Acids Research, vol. 28, pp. 235-242, 2000

[14] The Protein Data Bank Contents Guide: Atomic Coordinate Entry Format Description, http://www.rcsb.org/pdb/docs/format/pdbguide2.2/guide2.2 frame.html

Appendix A

Sample data files

The following sections give illustrations of the proprietary C_{α} data formats. They show generated files for the PDB file 4INS, describing the insulin hormone from pigs. The C_{α} atoms coordinates of chain B from molecule I has been processed.

The three first lines of C_{α} files always identify what is inside the file. The release of the corresponding algorithms is listed there for compatibility reasons (the data format has been modified throughout several releases).

The fourth line states the atom type (CA for C_{α}) to which the coordinates in the body of the file applies to. Line 5 gives the resolution used to determine the coordinates by experiment. If the value is 0.0 then the resolution information was not present in the PDB file. The entry column only serializes the atoms (the residue numbering sometimes doesn’t start with 1 or a C_{α} atom is missing in a residue).

A.1 C_{α} input data file

The C_{α} input data file is generated by EXTRACTOR by translating the PDB file (here 4INS.ent). If more chains are present they are identified by the chain entry. All coordinates are given in Å.

Atom c a r t e s i a n c o o r d i n a t e s from PDB
C−alpha−c i r c l e p r o j e c t . File created using Extractor .
Release\begin{matrix}c \\ ◯ \end{matrix} 1 . 3 1 ( 2 0 0 2 / 0 4 / 2 8 ) . 2002 EPFL Lausanne
Atom type : CA

ResolutionÅ [ ] : 1 . 5 0

Entry  xÅ [ ]  yÅ [ ]  zÅ [ ] r e s i d u e chain
1 −20.374  1 . 3 6 8  4 . 0 5 3  1  B
2 −18.754  0 . 5 9 8  7 . 4 0 6  2  B
3 −16.872  2 . 6 9 1  9 . 9 5 0  3  B
4 −13.137  2 . 5 4 2  9 . 5 7 1  4  B
5 −10.083  3 . 0 0 0  1 1 . 7 7 9  5  B
6  −7.468  4 . 7 0 9  9 . 6 1 1  6  B
7  −4.058  5 . 8 3 5  1 0 . 5 6 9  7  B
8  −1.070  7 . 2 6 1  8 . 6 3 2  8  B
9  −0.863  6 . 5 4 4  4 . 9 8 0  9  B
10  −3.940  4 . 3 7 9  5 . 1 8 8  10  B

A.2 C_{α} output data file (circles)  23

11  −6.071  7 . 5 1 8  5 . 6 1 7  11  B
12  −4.579  9 . 0 5 7  2 . 4 9 5  12  B
13  −5.432  5 . 9 4 6  0 . 5 4 2  13  B
14  −9.061  6 . 1 6 4  1 . 7 1 9  14  B
15  −9.224  9 . 8 5 4  0 . 7 1 7  15  B
16  −7.838  8 . 9 6 1  −2.686  16  B
17 −10.566  6 . 2 9 0  −3.167  17  B
18 −13.351  8 . 5 4 5  −1.933  18  B
19 −12.433  1 1 . 7 5 8  −3.756  19  B
20 −10.685  1 0 . 3 5 9  −6.851  20  B
21 −10.398  1 3 . 0 4 3  −9.501  21  B

A.2 C_{α} output data file (circles)

On line 6 and 7 of the file the calculated thickness in Åand the entry numbers of the atoms realizing this radius are given. The memory flag on line 8 corresponds to the number of smallest global radii which are listed - together with local radius, f-ratio - in the body of the file. The dislocation number ✠ is also given there (see section 4.2.2 for its definition).

The illustration below shows the human readable version of the file. The machine readable version contains the same information, but without formatting (as parenthesis and additional explanatorily text).

C i r c l i n g of p r o t e i n s
C−alpha−c i r c l e p r o j e c t . File created using C a l c u l a t o r .
Release\begin{matrix}c \\ ◯ \end{matrix} 1 . 6 1 ( 2 0 0 2 / 0 7 / 0 7 ) 2002 EPFL Lausanne
Atom type : CA

ResolutionÅ [ ] : 1 . 5 0 C i r c u l a r t h i c k n e s s Å [ ] : 2 . 6 0

Thickness combination : 1 6 , 17 and 18
Memory : 2
Entry  Local  f−r a t i o Dis  1 . g l o b a l r a d i u s  2 . g l o b a l r a d i u s
1  0 . 0 0 0  0 . 0 0 0  2 ( 1 . 1 1 4 ; 1 2 3 ) ( 1 . 1 3 5 ; 1 2 4)
2  1 . 1 1 4  0 . 9 7 3  2 ( 1 . 0 5 3 ; 2 3 4 ) ( 1 . 0 8 4 ; 2 4 5)
3  1 . 0 5 3  1 . 0 4 2  2 ( 1 . 0 5 3 ; 2 3 4 ) ( 1 . 0 7 6 ; 3 4 5)
4  1 . 0 7 6  1 . 0 0 8  2 ( 1 . 0 5 3 ; 2 3 4 ) ( 1 . 0 6 9 ; 4 5 6)
5  1 . 0 6 9  1 . 0 1 4  2 ( 1 . 0 4 5 ; 5 6 7 ) ( 1 . 0 6 9 ; 4 5 6)
6  1 . 0 4 5  0 . 9 6 2  5 ( 1 . 0 0 5 ; 6 10 1 1 ) ( 1 . 0 4 5 ; 5 6 7)
7  1 . 1 6 6  0 . 9 4 5  2 ( 1 . 0 4 5 ; 5 6 7 ) ( 1 . 0 9 0 ; 7 8 9)
8  1 . 0 9 0  0 . 9 9 5  4 ( 1 . 0 8 4 ; 5 8 9 ) ( 1 . 0 9 0 ; 7 8 9)
9  1 . 6 8 2  0 . 6 4 4  4 ( 1 . 0 8 2 ; 6 9 1 0 ) ( 1 . 0 8 4 ; 5 8 9)
10  1 . 8 4 8  0 . 5 4 4  5 ( 1 . 0 0 5 ; 6 10 1 1 ) ( 1 . 0 8 2 ; 6 9 10)
11  1 . 9 6 4  0 . 5 1 2  5 ( 1 . 0 0 5 ; 6 10 1 1 ) ( 1 . 2 6 4 ; 11 12 15)
12  1 . 7 4 8  0 . 6 2 5  2 ( 1 . 0 1 6 ; 12 13 1 4 ) ( 1 . 0 9 3 ; 12 14 15)
13  1 . 0 1 6  1 . 0 5 8  2 ( 1 . 0 1 6 ; 12 13 1 4 ) ( 1 . 0 3 3 ; 13 14 15)
14  1 . 0 3 3  1 . 0 4 4  2 ( 1 . 0 1 6 ; 12 13 1 4 ) ( 1 . 0 3 3 ; 13 14 15)
15  1 . 0 4 1  1 . 0 3 2  2 ( 1 . 0 3 3 ; 13 14 1 5 ) ( 1 . 0 3 5 ; 15 16 17)
16  1 . 0 3 5  1 . 0 3 8  2 ( 1 . 0 0 0 ; 16 17 1 8 ) ( 1 . 0 3 5 ; 15 16 17)
17  1 . 0 0 0  1 . 0 9 2  2 ( 1 . 0 0 0 ; 16 17 1 8 ) ( 1 . 0 2 2 ; 17 18 19)
18  1 . 0 2 2  1 . 0 9 4  2 ( 1 . 0 0 0 ; 16 17 1 8 ) ( 1 . 0 2 2 ; 17 18 19)
19  1 . 1 0 6  1 . 0 0 5  2 ( 1 . 0 2 2 ; 17 18 1 9 ) ( 1 . 1 0 6 ; 18 19 20)
20  1 . 4 3 1  0 . 7 7 7  2 ( 1 . 1 0 6 ; 18 19 2 0 ) ( 1 . 1 1 2 ; 17 19 20)
21  0 . 0 0 0  0 . 0 0 0  2 ( 1 . 4 3 1 ; 19 20 2 1 ) ( 1 . 6 5 5 ; 18 19 21)

A.3 C_{α} output data file (spheres)  24

A.3 C_{α} output data file (spheres)

The same remarks as in the previous section are valid.

Sphering of p r o t e i n s
C−alpha−c i r c l e p r o j e c t . File created using C a l c u l a t o r .
Release\begin{matrix}c \\ ◯ \end{matrix} 1 . 6 1 ( 2 0 0 2 / 0 7 / 0 7 ) . 2002 EPFL Lausanne
Atom type : CA

ResolutionÅ [ ] : 1 . 5 0 S p h e r i c a l t h i c k n e s s Å [ ] : 2 . 8 3

Thickness combination : 1 3 , 1 4 , 15 and 16
Memory : 2
Entry  Local  f−r a t i o Dis  1 . g l o b a l r a d i u s  2 . g l o b a l r a d i u s
1  0 . 0 0 0  0 . 0 0 0  4 ( 1 . 0 6 0 ; 1 2 3 4 ) ( 1 . 0 8 8 ; 1 2 4 5)
2  1 . 0 6 0  1 . 0 2 7  3 ( 1 . 0 0 9 ; 2 3 4 5 ) ( 1 . 0 6 0 ; 1 2 3 4)
3  1 . 0 0 9  1 . 0 8 2  3 ( 1 . 0 0 9 ; 2 3 4 5 ) ( 1 . 0 5 0 ; 3 4 5 6)
4  1 . 0 0 9  2 . 2 7 4  3 ( 1 . 0 0 9 ; 2 3 4 5 ) ( 1 . 0 2 7 ; 4 5 6 7)
5  1 . 0 2 7  1 . 0 4 8  3 ( 1 . 0 0 9 ; 2 3 4 5 ) ( 1 . 0 2 7 ; 4 5 6 7)
6  1 . 0 2 7  1 . 0 6 1  3 ( 1 . 0 2 7 ; 4 5 6 7 ) ( 1 . 0 5 0 ; 3 4 5 6)
7  1 . 0 7 0  1 . 0 0 6  3 ( 1 . 0 2 7 ; 4 5 6 7 ) ( 1 . 0 7 0 ; 5 6 7 8)
8  1 . 0 7 8  0 . 9 9 9  4 ( 1 . 0 7 0 ; 5 6 7 8 ) ( 1 . 0 7 7 ; 4 5 7 8)
9  1 . 9 9 7  0 . 5 4 6  4 ( 1 . 0 7 8 ; 6 7 8 9 ) ( 1 . 0 9 0 ; 5 6 8 9)
10  1 . 9 9 7  0 . 5 6 4  5 ( 1 . 1 2 6 ; 5 6 9 1 0 ) ( 1 . 1 7 7 ; 6 7 9 10)
11  3 . 6 9 0  0 . 3 3 9  6 ( 1 . 2 5 2 ; 5 6 10 1 1 ) ( 1 . 3 5 0 ; 5 6 9 11)
12  2 . 1 2 9  0 . 5 5 7  4 ( 1 . 0 1 7 ; 12 13 14 1 5 ) ( 1 . 1 8 6 ; 12 13 15 16)
13  1 . 0 1 7  1 . 1 0 0  3 ( 1 . 0 0 0 ; 13 14 15 1 6 ) ( 1 . 0 1 7 ; 12 13 14 15)
14  1 . 0 0 0  1 . 1 1 9  3 ( 1 . 0 0 0 ; 13 14 15 1 6 ) ( 1 . 0 1 4 ; 14 15 16 17)
15  1 . 0 0 0  1 . 1 8 6  3 ( 1 . 0 0 0 ; 13 14 15 1 6 ) ( 1 . 0 1 4 ; 14 15 16 17)
16  1 . 0 1 4  1 . 1 0 4  3 ( 1 . 0 0 0 ; 13 14 15 1 6 ) ( 1 . 0 1 4 ; 14 15 16 17)
17  1 . 0 3 7  1 . 0 8 0  3 ( 1 . 0 1 4 ; 14 15 16 1 7 ) ( 1 . 0 3 7 ; 16 17 18 19)
18  1 . 0 3 7  1 . 2 7 6  3 ( 1 . 0 3 7 ; 16 17 18 1 9 ) ( 1 . 0 4 6 ; 15 16 17 18)
19  1 . 0 5 2  1 . 3 1 9  3 ( 1 . 0 3 7 ; 16 17 18 1 9 ) ( 1 . 0 5 2 ; 17 18 19 20)
20  1 . 7 1 6  0 . 0 0 0  4 ( 1 . 0 5 2 ; 17 18 19 2 0 ) ( 1 . 3 8 7 ; 16 18 19 20)
21  0 . 0 0 0  0 . 0 0 0  4 ( 1 . 6 2 5 ; 17 18 19 2 1 ) ( 1 . 7 1 6 ; 18 19 20 21)

Appendix B Sourcecode of the algorithm

B.1 EXTRACTOR

! Projet de s i m u l a t i o n énumrique ! student : C h r i s t i a n S c h l a t t e r , Physique 4 e éanne EPFL

! p r o f e s s e u r : John Maddocks , DMA EPFL

!

! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗  5

!

! This r e l e a s e meets the v e r s i o n 2 . 2 o f the PDB format d e s c r i p t i o n guide
! which i s a v a i l a b l e on
! http : / /www. rcsb . org /pdb/ docs / format / pdbguide2 . 2 / guide2 . 2 frame . html
! Future m o d i f i c a t i o n o f the f i e l d l e n g t h s ( as proposed f o r the atom  10
! entry w i l l r e q u i r e a m o d i f i c a t i o n o f t h i s code .
! I ’m using Visual Fortran 6 . 6A from Compaq ( formerly D i g i t a l ) ,
! f u l l y Fortran 9 5 compatible compiler .
! I t should a l s o be runnnable on a f u l l y Fortran 9 0 compiler .  15
! For r e l e a s e information see module r e l e a s e
module r e l e a s e

20

! Release information
character (LEN=25) , parameter : : r e l e a s e i n f o = ’ Release 1 . 3 1 ( 2 0 0 2 / 0 4 / 2 8 ) ’
end module r e l e a s e

25

module var
! Data type d e f i n i t i o n s
type carbons

30

integer : : r e s i d u e
real , dimension ( 3 ) : : c o o r d i n a t e s
character (LEN=1) : : c h a i n i d e n t
end type carbons  35
end module var
module reading

40

Projet de simulation numérique 2002  25  Christian Schlatter

B.1 EXTRACTOR  26

! Subroutines
implicit none

contains

45

subroutine h e l i c e s ( filename , r e s o l u t i o n , alphacarbons , talkmode , atomtype )
! Extracting o f alpha h e l i c e s
use var  50
implicit none
! V a r i a b l e s
integer : : i , j , k , sequencenumber , 55
integer  : : io , oldatom , newatom
integer  : : numbercarbon ,
integer  : : beginsequence ,
integer  : : endsequence , numberhelix
integer , allocatable , dimension ( : , : ) : : chainseq 60
! l i n e s : found h e l i c e s
! begin at sequence end at sequence
character (LEN=1)  : : test3 , ident
65 character (LEN=2)  : : h e l i x t y p e
character (LEN=4)  : : atomname , seq , beginseq ,
character (LEN=4)  : : endseq , test6 , atomtype
character (LEN=5)  : : res , t e s t 4
character (LEN=6)  : : atom , h e l i x
70 character (LEN=8)  : : x , y , z , t e s t 2
character (LEN=9)  : : t e s t 1
character (LEN=10)  : : t e s t 5
character (LEN=22)  : : t e s t
character (LEN=32)  : : t e s t 3 2
75 character (LEN=60)  : : molecule
character ∗ ( ∗) , intent ( in )  : : filename
character (LEN=1) , allocatable , dimension ( : , : ) : : chain
! l i n e s : found h e l i c e s
! rows : l o c a t i o n , o r i e n t a t i o n

80

real  : : xpos , ypos , zpos
real , intent ( out )  : : r e s o l u t i o n
type ( carbons ) , pointer , dimension ( : )  : : alphacarbons , dummy

85

logical  : : f l a g , talkmode

OPEN (UNIT=1 , FILE=filename , FORM=’FORMATTED’ , STATUS=’OLD’ )

90

rewind 1
! Reading o f the PDB f i l e header l i n e
95 read ( 1 , ’ ( a10 , a60 ) ’ , ADVANCE=’YES’ ) test5 , molecule
i f ( talkmode ) then

     write ( ∗ , ’ ( "    The following protein is going to be analysed:" ) ’ )

write ( ∗ , ∗)  100
write ( ∗ , ’ ( "    " , a60 ) ’ ) molecule
end i f
105 read ( 1 , ’ ( a10 , a60 ) ’ , ADVANCE=’YES’ ) test5 , molecule
i f ( talkmode ) then
write ( ∗ , ’ ( "    " , a60 ) ’ ) molecule
write ( ∗ , ∗)  110
end i f

do

115

! Reading o f the experimental data r e s o l u t i o n value
read ( 1 , ’ ( a22 , a5 ) ’ , ADVANCE=’YES’ ) t e s t , r e s
120 i f ( t e s t == ’REMARK 2 RESOLUTION. ’ ) then

exit

endif

125

end do
i f ( r e s == ’ NOT ’ ) then
i f ( talkmode ) then  130

     write ( ∗ , ’ ( "    No resolution specified !!! " ) ’ )

write ( ∗ , ∗)
end i f  135

else

read ( res , ’ ( f 5 . 0 ) ’ ) r e s o l u t i o n

140

i f ( talkmode ) then

     write ( ∗ , ’ ( "    The structure data has a resolution of " , f 4 . 2 , "

angstroems" ) ’ ) r e s o l u t i o n
write ( ∗ , ∗)

145

end i f
end i f
! Determination o f the number o f present alpha h e l i c e s  150
rewind 1
numberhelix = 0

155

do

read ( 1 , ’ ( a6 , a32 , a2 ) ’ , ADVANCE=’YES’ , IOSTAT = i o ) h e l i x , test32 ,

h e l i x t y p e

i f ( i o < 0) then  160
i f ( ( talkmode ) . and . ( numberhelix == 0) ) then

    write ( ∗ , ’ ( "    No alpha helices in this file !!" ) ’ )

write ( ∗ , ∗)  165

return

end i f

170 exit

end i f
i f ( ( h e l i x == ’HELIX ’ ) . and . ( ( h e l i x t y p e == ’ 6 ’) . or . ( h e l i x t y p e
== ’ 1 ’) ) ) then

175

! Alpha h e l i x types : 1 means right −handed ; 6 means l e f t −handed h e l i x
numberhelix = numberhelix + 1
end i f  180
end do
allocate ( chainseq ( numberhelix , 2 ) )
allocate ( chain ( numberhelix , 2 ) )  185
! Determination o f l o c a t i o n o f the alpha h e l i x
rewind 1

190

i = 1
chain = X
f l a g = . f a l s e .
HEL : do  195
i f ( i > numberhelix ) then
! All alpha h e l i c e s have been read

200

exit

end i f
205 read ( 1 , ’ ( a6 , a13 , a1 , a1 , a4 , a8 , a4 , a1 , a2 ) ’ , ADVANCE=’YES’ , &
IOSTAT = i o ) h e l i x , test1 , ident , test3 , beginseq , &
test2 , endseq , test3 , h e l i x t y p e
i f ( ( h e l i x == ’HELIX ’ ) . and . ( ( h e l i x t y p e == ’ 6 ’) . or . ( h e l i x t y p e
== ’ 1 ’) ) ) then

210

chain ( i , 1 ) = ident

read ( beginseq , ’ ( i 4 ) ’ ) chainseq ( i , 1 ) ! Helix s t a r t i n g sequence

number

read ( endseq , ’ ( i 4 ) ’ ) chainseq ( i , 2 )  ! Helix stopping sequence

number

215

f l a g = . true .
i f ( h e l i x t y p e == ’ 1 ’) then
i f ( talkmode ) then  220

  write ( ∗ , ’ ( "    Right-handed alpha helix found in chain " , a1 ) ’ )

chain ( i , 1 )
end i f

225

chain ( i , 2 ) = ’R’
i = i + 1

else

230

i f ( talkmode ) then

  write ( ∗ , ’ ( "    Left-handed alpha helix foundin chain " , a1 ) ’ )

chain ( i , 1 )
end i f  235
chain ( i , 2 ) = ’L ’
i = i + 1

240 endif

endif

i f ( ( h e l i x /= ’HELIX ’ ) . and . ( f l a g ) ) then

245

exit HEL ! End o f h e l i x s p e c i f i c a t o r l i s t reached

endif

end do HEL  250
i f ( talkmode ) then
write ( ∗ , ∗)
write ( ∗ , ’ ( "    " , i3 , " helices found." ) ’ ) i − 1  255
write ( ∗ , ∗)
end i f
j = 1  260
do k = 1 , i − 1
numbercarbon = chainseq (k , 2 ) − chainseq (k , 1 ) + 1

265

i f ( a s s o c i a t e d ( alphacarbons ) ) then

! Extending the v a r i a b l e ’ alphacarbons ’ f o r the new h e l i x
allocate (dummy( size ( alphacarbons ) ) )  270
dummy = alphacarbons
deallocate ( alphacarbons )
allocate ( alphacarbons ( size (dummy) + numbercarbon ) )
alphacarbons ( 1 : j −1) = dummy ( : )
deallocate (dummy)  275

else

! A l l o c a t i n g the v a r i a b l e ’ alphacarbons ’ f o r the f i r s t h e l i x

280

allocate ( alphacarbons ( numbercarbon ) )
end i f
i f ( talkmode ) then  285

write ( ∗ , ’ ( "    There are " , i10 , " carbon atoms located in the &
                  alpha helix of chain " , a1 ) ’ ) &

numbercarbon , chain (k , 1 )
write ( ∗ , ∗)  290

write ( ∗ , ’ ( "    The carbon atoms have the following coordinates:" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( "        #    SEQ    x [A]    y [A]    z [A]" ) ’ )

write ( ∗ , ∗)

295

end i f
rewind 1
sequencenumber = 0
oldatom = −10  300
ATM : do
! Reading o f the c o o r d i n a t e s o f the atoms l o c a t e d in the alpha h e l i c e s

305

read ( 1 , ’ ( a6 , a5 , a5 , a5 , a1 , a4 , a4 , a8 , a8 , a8 ) ’ , ADVANCE=’YES’ ) &
atom , test4 , atomname , test5 , ident , seq , test6 , x , y , z
i f ( ( atom == ’ATOM ’ ) . and . ( ident == chain (k , 1 ) ) ) then

310

read ( test4 , ’ ( i 5 ) ’ ) newatom
i f ( ( atomname == atomtype ) . and . ( newatom == oldatom + 1) ) then
! Error in PDB f i l e or strange mutation : The same atom found as  315

neighbours

oldatom = newatom

cycle

end i f  320
i f ( atomname == atomtype ) then
oldatom = newatom

325

end i f
read ( seq , ’ ( i 4 ) ’ ) sequencenumber
read (x , ’ ( f 8 . 0 ) ’ ) xpos
read (y , ’ ( f 8 . 0 ) ’ ) ypos  330
read ( z , ’ ( f 8 . 0 ) ’ ) zpos

endif

335 i f ( ( atom == ’ATOM ’ ) . and . ( atomname == atomtype ) . and . &
& ( sequencenumber >= chainseq (k , 1 ) ) . and . ( sequencenumber <=
chainseq (k , 2 ) ) ) then
! Atom o f type ’ atomtype ’ found i n s i d e a h e l i x
i f ( talkmode ) then  340
write ( ∗ , ’ ( "  " , i5 , i6 , f 8 . 3 , f 8 . 3 , f 8 . 3 ) ’ ) j , sequencenumber ,
xpos , ypos , zpos
end i f
alphacarbons ( j )%r e s i d u e = sequencenumber
345 alphacarbons ( j )%c o o r d i n a t e s ( 1 ) = xpos
alphacarbons ( j )%c o o r d i n a t e s ( 2 ) = ypos
alphacarbons ( j )%c o o r d i n a t e s ( 3 ) = zpos
alphacarbons ( j )%c h a i n i d e n t = ident
j = j + 1

350

endif

i f ( sequencenumber >= chainseq (k , 2 ) . and . ( atomname == atomtype ) ) then
! End o f h e l i x reached  355
i f ( talkmode ) then
write ( ∗ , ∗)

360

end i f
exit ATM

365 endif

end do ATM
end do

370

CLOSE (UNIT=1 , STATUS=’KEEP’ )
deallocate ( chainseq , chain )
end subroutine h e l i c e s  375
subroutine backbones ( filename , r e s o l u t i o n , alphacarbons , talkmode , atomtype )
! Extracting the backbone o f the p r o t e i n

380

use var

implicit none
! V a r i a b l e s  385
integer , parameter  : : maxhelices = 10
! maximum p o s s i b l e number o f h e l i c e s
integer  : : i , j , k , sequencenumber , io ,  390
integer  : : oldatom , newatom
integer  : : oldsequencenumber
integer  : : numbercarbon , beginsequence ,
integer  : : endsequence , r e s i d u e s
integer , dimension ( maxhelices , 2 ) : : chainseq 395
! l i n e s : found h e l i c e s
! begin at sequence end at sequence
character (LEN=1)  : : ident
400 character (LEN=4)  : : atomname , seq , test6 , r e s i
character (LEN=4)  : : atomtype
character (LEN=5)  : : res , t e s t 4
character (LEN=6)  : : atom , h e l i x , s e q r e s
character (LEN=7)  : : t e s t 1
405 character (LEN=8)  : : x , y , z
character (LEN=10)  : : t e s t 5
character (LEN=22)  : : t e s t
character (LEN=60)  : : molecule
character ∗ ( ∗) , intent ( in )  : : filename
410 character (LEN=1) , dimension ( maxhelices , 2 ) : : chain
! l i n e s : found h e l i c e s
! rows : l o c a t i o n , o r i e n t a t i o n
real  : : xpos , ypos , zpos
real , intent ( out )  : : r e s o l u t i o n  415
type ( carbons ) , pointer , dimension ( : )  : : alphacarbons , dummy
logical  : : f l a g , talkmode

420

OPEN (UNIT=1 , FILE=filename , FORM=’FORMATTED’ , STATUS=’OLD’ )
! Reading o f the PDB f i l e header l i n e

425

rewind 1
read ( 1 , ’ ( a10 , a60 ) ’ , ADVANCE=’YES’ ) test5 , molecule
i f ( talkmode ) then  430

     write ( ∗ , ’ ( "    The following protein is going to be analysed:" ) ’ )

write ( ∗ , ∗)
write ( ∗ , ’ ( "    " , a60 ) ’ ) molecule

435

end i f
read ( 1 , ’ ( a10 , a60 ) ’ , ADVANCE=’YES’ ) test5 , molecule
i f ( talkmode ) then  440

write ( ∗ , ’ ( "    " , a60 ) ’ ) molecule
write ( ∗ , ∗)
end i f  445

do

! Reading ot the experimental data r e s o l u t i o n value

450

read ( 1 , ’ ( a22 , a5 ) ’ , ADVANCE=’YES’ ) t e s t , r e s
i f ( t e s t == ’REMARK 2 RESOLUTION. ’ ) then

455 exit endif

end do

460

i f ( r e s == ’ NOT ’ ) then
i f ( talkmode ) then

     write ( ∗ , ’ ( "    No resolution specified !!! " ) ’ )  465

write ( ∗ , ∗)
end i f

470 else

read ( res , ’ ( f 5 . 0 ) ’ ) r e s o l u t i o n
i f ( talkmode ) then

475

     write ( ∗ , ’ ( "    The structure data has a resolution of " , f 4 . 2 , "

angstroems" ) ’ ) r e s o l u t i o n
write ( ∗ , ∗)
end i f

480

end i f

do

! Determination o f the number o f r e s i d u e s in the backbone  485
read ( 1 , ’ ( a6 , a7 , a4 ) ’ , ADVANCE=’YES’ ) seqres , test1 , r e s i
i f ( s e q r e s == ’SEQRES’ ) then

490

exit endif

end do  495
read ( r e s i , ’ ( i 4 ) ’ ) r e s i d u e s
j = 1
Projet de simulation numérique 2002  Christian Schlatter

500

allocate ( alphacarbons ( r e s i d u e s ) )
i f ( talkmode ) then

write ( ∗ , ’ ( "    There are " , i10 , " residues in the backbone." ) ’ )  505

r e s i d u e s

write ( ∗ , ∗)

write ( ∗ , ’ ( "    The carbon atoms have the following coordinates:" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( "        #    SEQ    x [A]    y [A]    z [A]" ) ’ )

write ( ∗ , ∗)  510
end i f
rewind 1
sequencenumber = 0  515
oldatom = −10
ATM : do
oldsequencenumber = sequencenumber  520
! Reading o f the c o o r d i n a t e s o f the atoms
read ( 1 , ’ ( a6 , a5 , a5 , a5 , a1 , a4 , a4 , a8 , a8 , a8 ) ’ , ADVANCE=’YES’ , &
IOSTAT = i o ) atom , test4 , atomname , test5 , ident , seq , &  525
test6 , x , y , z
i f ( atom == ’ATOM ’ ) then
read ( test4 , ’ ( i 5 ) ’ ) newatom  530
i f ( ( atomname == atomtype ) . and . ( newatom == oldatom + 1) ) then
! Error in PDB f i l e or strange mutation : The same atom found as

neighbours

535

oldatom = newatom

cycle

end i f

540

i f ( atomname == atomtype ) then
oldatom = newatom
end i f  545
read ( seq , ’ ( i 4 ) ’ ) sequencenumber
read (x , ’ ( f 8 . 0 ) ’ ) xpos
read (y , ’ ( f 8 . 0 ) ’ ) ypos
read ( z , ’ ( f 8 . 0 ) ’ ) zpos  550

endif

i f ( ( sequencenumber < oldsequencenumber ) . or . ( i o < 0) ) then

555

! I f there were more than one atom o f the chosen type per

! residue , the ’ alphacarbons ’ v a r i a b l e has to be extended
allocate (dummy( oldsequencenumber ) )
dummy = alphacarbons ( 1 : oldsequencenumber )  560
deallocate ( alphacarbons )
allocate ( alphacarbons ( size (dummy) ) )
alphacarbons = dummy
deallocate (dummy)

565

write ( ∗ , ∗)
exit ATM
end i f

570

i f ( ( atom == ’ATOM ’ ) . and . ( atomname == atomtype ) ) then
! An atom has been found
i f ( j > r e s i d u e s ) then  575
! Extending the v a r i a b l e ’ alphacarbons ’ i f the number o f
! r e s i d u e s given in the PDB f i l e header was not c o r r e c t .
allocate (dummy( size ( alphacarbons ) ) )  580
dummy = alphacarbons
deallocate ( alphacarbons )
allocate ( alphacarbons ( size (dummy) + 1) )
alphacarbons ( 1 : j −1) = dummy ( : )
deallocate (dummy)  585
end i f
alphacarbons ( j )%r e s i d u e = sequencenumber
590 alphacarbons ( j )%c o o r d i n a t e s ( 1 ) = xpos
alphacarbons ( j )%c o o r d i n a t e s ( 2 ) = ypos
alphacarbons ( j )%c o o r d i n a t e s ( 3 ) = zpos
alphacarbons ( j )%c h a i n i d e n t = ident
i f ( talkmode ) then  595
write ( ∗ , ’ ( "  " , i5 , i6 , f 8 . 3 , f 8 . 3 , f 8 . 3 ) ’ ) j , sequencenumber ,
xpos , ypos , zpos
end i f

600

j = j + 1

endif

i f ( ( ( sequencenumber == r e s i d u e s ) . and . ( atomname == atomtype ) ) . or . (  605
atom == ’TER ’ ) ) then
! All atoms have been processed .
i f ( talkmode ) then

610

write ( ∗ , ∗)
end i f

i f ( j < r e s i d u e s ) then  615

! Correction i f number o f found alphacarbons i s s m a l l e r
! than the s i z e o f i t s a l l o c a t e d v a r i a b l e

620

allocate (dummy( j −1) )
dummy = alphacarbons ( 1 : j −1)
deallocate ( alphacarbons )
allocate ( alphacarbons ( size (dummy) ) )
alphacarbons ( : ) = dummy ( : )  625
deallocate (dummy)
end i f
exit ATM  630

endif

end do ATM

635

CLOSE (UNIT=1 , STATUS=’KEEP’ )
end subroutine backbones
end module reading  640
! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
program e x t r a c t

645

! Main program
use reading
use var
use r e l e a s e  650
implicit none
! V a r i a b l e s
655 character (LEN=12)  : : export , dummy
character (LEN=3)  : : subfolder , input
character (LEN=4)  : : extension , atomtype
character (LEN=6)  : : e x p o r t f o l d e r
character (LEN=12) , allocatable , dimension ( : ) : : filename
character , dimension ( : ) , pointer  : : name  660
real  : : r e s o l u t i o n , min
real ( 8 )  : : i n i t i a l , f i n a l
real , allocatable , dimension ( : , : )  : : a l l a l p h a c a r b o n s , oldcarbons
real , pointer , dimension ( : , : ) : : a l p h a c i r c l e s 665
real , allocatable , dimension ( : )  : : minradius
type ( carbons ) , pointer , dimension ( : )  : : alphacarbons
integer  : : i , j , k , numberfiles , io ,  670
integer  : : position , length
integer  : : what , r e s o u r c e f l a g

integer , allocatable , dimension ( : , : )  : : t r i p l e
logical  : : talkmode  675
logical , dimension ( : ) , allocatable  : : g o o d f i l e s
! r e s o u r c e f l a g = 0  ! e i t h e r 0 : a s i n g l e PDB− f i l e s p e c i f i e d by filename
! or  1 : a bunch o f PDB−f i l e s l i s t e d in F i l e s . i n i
!  i s loaded .  680
s u b f o l d e r = ’PDB’  ! name o f the s u b d i r e c t o r y where the PDB−r e s o u r c e
! f i l e s are s t o r e d in ( with e x t e n s i o n . ent or . pdb )
e x p o r t f o l d e r = ’ source ’ ! name o f the s u b f o l d e r f o r ASCII export f i l e s  685
talkmode = . true .  ! i f true , d e t a i l e d comments are printed on s c r e e n
i f ( talkmode ) then

690

write ( ∗ , ∗)

  write ( ∗ , ’ ( " Computational physics project" ) ’ )
  write ( ∗ , ’ ( " Christian Schlatter, student of physics, 4th year" ) ’ )

write ( ∗ , ∗) r e l e a s e i n f o
write ( ∗ , ’ ( "  695

       ***************************************************************" ) ’ )

write ( ∗ , ∗)
end i f
! PDB input f i l e c o l l e c t i o n  700

write ( ∗ , ’ ( " PDB File coordinates extractor" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " How many PDB files would you like to process ?" ) ’ )

write ( ∗ , ∗)  705

write ( ∗ , ’ ( " -> Type 1 for a single file or" ) ’ )
write ( ∗ , ’ ( " -> Type 2 for a list of files provided in Files.ini" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " Please choose your option now : " ) ’ , advance = ’ no ’ )

710

read ( ∗ , ∗) r e s o u r c e f l a g
write ( ∗ , ∗)
i f ( r e s o u r c e f l a g == 2) then ! Batch mode using F i l e s . i n i  715
write ( ∗ , ’ ( " Opening Files.ini... -> " ) ’ , advance = ’ no ’ )
open (UNIT=1 , FILE=’ F i l e s . i n i ’ , FORM=’FORMATTED’ , STATUS=’OLD’ )

720

rewind 1
n u m b e r f i l e s = 0

725 do

! Determination o f the number o f PDB f i l e s
read ( 1 , ’ ( a12 ) ’ , IOSTAT = i o ) dummy

730

i f ( i o == 0) then
n u m b e r f i l e s = n u m b e r f i l e s + 1

735 else

exit

end i f

740

end do
allocate ( filename ( n u m b e r f i l e s ) )
rewind 1  745
do i = 1 , n u m b e r f i l e s
! Determination o f the PDB f i l e names

750

read ( 1 , ’ ( a12 ) ’ ) filename ( i )
end do
755 close (UNIT=1 , STATUS=’KEEP’ )
write ( ∗ , ’ ( i10 , " files to process" ) ’ ) n u m b e r f i l e s
write ( ∗ , ∗)

760 else

allocate ( filename ( 1 ) )
! User provided PDB f i l e name

765

  write ( ∗ , ’ ( " Please provide the filename of the PDB file : " ) ’ , advance = ’

no ’ )
read ( ∗ , ∗) filename ( 1 )
write ( ∗ , ∗)
n u m b e r f i l e s = 1  770

endif

allocate ( g o o d f i l e s ( n u m b e r f i l e s ) )
g o o d f i l e s = . true .  775
! Atom type to be e x t r a c t e d
input = ’  ’

780

write ( ∗ , ∗)

write ( ∗ , ’ ( " Which type of atoms would you like to extract ?" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " Example : For alpha carbons type the two characters CA" ) ’ )

write ( ∗ , ∗)  785

write ( ∗ , ’ ( " Please choose your atom now : " ) ’ , advance = ’ no ’ )

read ( ∗ , ∗) input
Projet de simulation numérique 2002  Christian Schlatter

atomtype ( 1 : 1 ) = " "  790
do i = 1 , 3
atomtype ( i +1: i +1) = input ( i : i )

795

end do

! S p e c i a l t e r t i a r y c o n f i g u r a t i o n s e l e c t i o n

800

write ( ∗ , ∗)

write ( ∗ , ’ ( " What kind of information would you like to extract ?" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " -> Type 1 for coordinates of atoms located in alpha helices." ) ’ )
write ( ∗ , ’ ( " -> Type 2 for coordinates of atoms of the whole backbone of the  805
protein." ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " Please choose your option now : " ) ’ , advance = ’ no ’ )

read ( ∗ , ∗) what

810

do i = 1 , n u m b e r f i l e s
! Processing f i l e by f i l e
i f ( talkmode ) then  815
write ( ∗ , ∗)
write ( ∗ , ∗)

write ( ∗ , ’ ( " Extracting structure information from file " , a12 ) ’ )

filename ( i )

write ( ∗ , ’ ( " *******************************************************" ) ’ )  820

write ( ∗ , ∗)
write ( ∗ , ’ ( " Processing now file " , a12 ) ’ ) filename ( i )
write ( ∗ , ∗)
end i f  825
position = scan ( filename ( i ) , ’ . ’ )
select case ( what )

830

case ( 1 )
! Extracting h e l i x subsequences only
c a l l h e l i c e s ( s u b f o l d e r // ’/ ’// filename ( i ) , r e s o l u t i o n , alphacarbons ,  835
talkmode , atomtype )
filename ( i ) ( position +1: position+1) = ’h ’
filename ( i ) ( position +2: position+2) = ’ e ’
filename ( i ) ( position +3: position+3) = ’ l ’

840

case ( 2 )
! Extracting the whole backbone o f the p r o t e i n

c a l l backbones ( s u b f o l d e r // ’/ ’// filename ( i ) , r e s o l u t i o n , alphacarbons ,  845
talkmode , atomtype )
filename ( i ) ( position +1: position+1) = ’b ’
filename ( i ) ( position +2: position+2) = ’ c ’
filename ( i ) ( position +3: position+3) = ’ k ’

850

case default
end select
i f ( . not . ( a s s o c i a t e d ( alphacarbons ) ) ) then  855
! This may occur when no alpha h e l i c e s were present in the s e l e c t e d

p r o t e i n

i f ( talkmode ) then

860

     write ( ∗ , ’ ( " This file will be ignored. Please check if it is
        applicable." ) ’ )

write ( ∗ , ∗)

pause 865 else

     write ( ∗ , ’ ( "File " , a8 , " has not been created. No helices present." ) ’ )

filename ( i )
end i f

870

g o o d f i l e s ( i ) = . f a l s e .

cycle

end i f  875
! Saving the gathered i n f o r m a t i o n s to disk
write ( ∗ , ∗)

write ( ∗ , ’ ( " Exporting C-alpha coordinates to file " , a12 ) ’ ) filename ( i )  880

i f ( talkmode ) then

  write ( ∗ , ’ ( " **************************************************" ) ’ )

write ( ∗ , ∗)  885
end i f
open (UNIT=2 , FILE=e x p o r t f o l d e r // ’/ ’// filename ( i ) , FORM=’FORMATTED’ , STATUS
=’REPLACE’ )

890

! D e f i n i t i o n o f the C−alpha c i r c l e input f i l e header

write ( 2 , ’ ( "Atom cartesian coordinates from PDB" ) ’ )
write ( 2 , ’ ( "C-alpha-circle project. File created using Extractor." ) ’ )

write ( 2 , ’ ( a25 ) ’ , advance=’no ’ ) r e l e a s e i n f o  895

write ( 2 , ’ ( "\begin{matrix}c \\ ◯ \end{matrix}.  2002 EPFL Lausanne" ) ’ )

write ( 2 , ’ ( "Atom type : " , a4 ) ’ ) atomtype

write ( 2 , ’ ( "ResolutionÅ [] : " , f 4 . 2 ) ’ ) r e s o l u t i o n

Projet de simulation numérique 2002  Christian Schlatter

write ( 2 , ∗)

900 write ( 2 , ’ ( "Entry      xÅ []      yÅ []      zÅ []  residue   chain" ) ’ )

write ( 2 , ∗)
do j = 1 , size ( alphacarbons )
write ( 2 , ’ ( i5 , 3 f10 . 3 , i9 , a7 ) ’ ) j , alphacarbons ( j )%c o o r d i n a t e s ( 1 ) , &
905 alphacarbons ( j )%c o o r d i n a t e s ( 2 ) , &
alphacarbons ( j )%c o o r d i n a t e s ( 3 ) , &
alphacarbons ( j )%residue , &
alphacarbons ( j )%c h a i n i d e n t
end do

910

close (UNIT=2 , STATUS=’KEEP’ )
deallocate ( alphacarbons )
end do  915
! P r o t o c o l l i n g o f which f i l e s were processed s u c e s s f u l l y
open (UNIT=3 , FILE = ’ g o o d f i l e s . txt ’ , FORM=’FORMATTED’ , STATUS=’REPLACE’ )

920

do i = 1 , n u m b e r f i l e s
i f ( g o o d f i l e s ( i ) ) then
write ( 3 , ∗) filename ( i )  925
end i f
end do

930

close (UNIT=3 , STATUS=’KEEP’ )
deallocate ( filename )
write ( ∗ , ∗)  935
write ( ∗ , ’ ( " DONE ! " ) ’ )
write ( ∗ , ∗)
end program e x t r a c t

B.2 CALCULATOR  42

B.2 CALCULATOR

! Projet de s i m u l a t i o n énumrique

! student : C h r i s t i a n S c h l a t t e r , SB , s e c t i o n physique , EPFL

! p r o f e s s e u r : John Maddocks , SB , s e c t i o n émathmatique , EPFL

! This r e l e a s e meets the v e r s i o n 2 . 2 o f the PDB format d e s c r i p t i o n guide  5
! which i s a v a i l a b l e on
! http : / /www. rcsb . org /pdb/ docs / format / pdbguide2 . 2 / guide2 . 2 frame . html
! Future m o d i f i c a t i o n s o f the f i e l d l e n g t h s ( as proposed f o r the atom entry )
! w i l l r e q u i r e a m o d i f i c a t i o n o f t h i s code .

10

module parameters
! User a d j u s t a b l e parameter module
integer , parameter : : memory  =  12 ! Memory l i m i t o f s m a l l e s t r a d i i
integer , parameter : : l i m i t  = 200  ! Upper l i m i t f o r d i s t r i b u t i o n  15

vector

logical , parameter : : talkmode = . f a l s e . ! Comments on s c r e e n
logical , parameter : : secondrun = . true . ! S t a t i s t i c s run to c l a s s i f y
d i s t r i b u t i o n o f r a d i i
20 character (LEN=25) , parameter : : r e l e a s e i n f o = ’ Release 1 . 6 1 ( 2 0 0 2 / 0 7 / 0 7 ) ’
end module parameters
! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗

25

module var
! Data type d e f i n i t i o n s
type carbons

30

integer  : : r e s i d u e
real , dimension ( 3 ) : : c o o r d i n a t e s
character (LEN=1) : : c h a i n i d e n t
end type carbons  35
end module var
! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗

40

module mathe
! Mathematical r o u t i n e s

contains

45

function crossproduct ( a , b)
! C a l c u l a t e s the c r o s s product o f two v e c t o r s
implicit none

50

real , dimension ( 3 ) , intent ( in ) : : a , b
real , dimension ( 3 ) : : crossproduct
crossproduct ( 1 ) = a ( 2 ) ∗b ( 3 ) − a ( 3 ) ∗b ( 2 )
Projet de simulation numérique 2002  Christian Schlatter

B.2 CALCULATOR  43

crossproduct ( 2 ) = a ( 3 ) ∗b ( 1 ) − a ( 1 ) ∗b ( 3 )  55
crossproduct ( 3 ) = a ( 1 ) ∗b ( 2 ) − a ( 2 ) ∗b ( 1 )
end function crossproduct
function determinant (A)  60
! C a l c u l a t e s the determinant o f a matrix using
! the p a r t i a l p i v o t i n g Gaussian e l i m i n a t i o n scheme .
implicit none  65
real , pointer , dimension ( : , : ) : : A
real , allocatable , dimension ( : , : ) : : B
real , allocatable , dimension ( : ) : : C

70

integer , allocatable , dimension ( : ) : : indexer
integer : : i , j , k , n , tempindex , msgn
real : : s c a l e r , pivot , p i v o t i , p i v o t j
real : : determinant  75
n = size (A( : , 1 ) )
allocate (B(n , n) )
B = A  80
! index r e c o r d s the p i v o t i n g order
allocate ( indexer (n) )

85

! A i s the o r i g i n a l matrix in the input . B i s A transformed
! i n t o the matrix plus the p i v o t i n g element r a t i o s below
! the diagonal in the output .
do i = 1 , n  90
indexer ( i ) = i
end do

95

! Find the r e s c a l i n g f a c t o r s , one from each row
allocate (C(n) )
do i = 1 , n  100
s c a l e r = 0.0
do j = 1 , n

105

s c a l e r = amax1( s c a l e r , abs (B( i , j ) ) )
end do
C( i ) = s c a l e r  110
end do

! Search the p i v o t i n g ( l a r g e s t ) element from each column

115

do j = 1 , n−1
pivot = 0.0
do i = j , n  120
p i v o t i = abs (B( indexer ( i ) , j ) ) / C( indexer ( i ) )
i f ( p i v o t i .GT. pivot ) then

125

pivot = p i v o t i
k = i

endif

130

end do
! Interchange the rows via indexer (n) to record p i v o t i n g order
tempindex = indexer ( j )  135
indexer ( j ) = indexer ( k )
indexer ( k ) = tempindex
do i = j +1 , n

140

p i v o t j = B( indexer ( i ) , j ) /B( indexer ( j ) , j )
! Record p i v o t i n g r a t i o s below the diagonal
B( indexer ( i ) , j ) = p i v o t j  145
! Modify other elements a c c o r d i n g l y
do k = j +1 , n

150

B( indexer ( i ) , k ) = B( indexer ( i ) , k )−p i v o t j ∗B( indexer ( j ) , k )
end do
end do  155
end do
determinant = 1.0

160

do i = 1 , n
determinant = determinant ∗B( indexer ( i ) , i )
end do  165
msgn = 1
do i = 1 , n

170

do while ( i .NE. indexer ( i ) )

msgn = −msgn
j = indexer ( i )
indexer ( i ) = indexer ( j )  175
indexer ( j ) = j
end do
end do  180
determinant = msgn∗ determinant
deallocate ( indexer )
deallocate (B)  185
deallocate (C)
end function determinant
end module mathe  190
! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
module d i s t a n c e

195

implicit none

contains

subroutine c i r c l e s ( alphacarbons , f , g , t r u e l o c a l , g l o b a l r a d i u s , &  200
g l o b a l l o c a t i o n , s m a l l e s t r a d i u s , s m a l l e s t l o c a t i o n , &
d i s t r i b u t i o n , t h i c k n e s s )
! f i t s c i r c l e s through the macromolecule
use var  205
use mathe
use parameters
implicit none

210

type ( carbons ) , pointer , dimension ( : ) : : alphacarbons
real  : : r2 , A, B, C, t h i c k n e s s
real , dimension ( 3 )  : : D, f i r s t , second , t h i r d
real , dimension ( 3 ) : : crossed , r a d i i 215
real , pointer , dimension ( : )  : : f , g , l o c a l r a d i u s , nonadjacentradius
real , pointer , dimension ( : )  : : g l o b a l r a d i u s , n o n l o c a l r a d i u s
real , pointer , dimension ( : )  : : p a r t l y a d j a c e n t r a d i u s , t r u e l o c a l
real , pointer , dimension ( : , : )  : : s m a l l e s t r a d i u s

220

integer , pointer , dimension ( : )  : : f i r s t l o c a l
integer , pointer , dimension ( : , : )  : : n o n a d j a c e n t l o c a t i o n
integer , pointer , dimension ( : , : )  : : p a r t l y a d j a c e n t l o c a t i o n
integer , pointer , dimension ( : , : )  : : g l o b a l l o c a t i o n
integer , pointer , dimension ( : , : , : ) : : s m a l l e s t l o c a t i o n 225
integer , dimension ( 3 )  : : combination
integer  : : numbercarbon , i , j , k , l , m, q , r
integer  : : mem, index , pos , l o c a t i o n
integer , dimension ( 1 )  : : s e l e c t o r
! due to c o m p a t i b i l i t y i s s u e s with Fortran 90  230

integer ( 4 ) , dimension ( l i m i t )  : : d i s t r i b u t i o n
numbercarbon = size ( alphacarbons )

235

allocate ( l o c a l r a d i u s ( numbercarbon ) )
allocate ( t r u e l o c a l ( numbercarbon ) )
allocate ( nonadjacentradius ( numbercarbon ) )
allocate ( f i r s t l o c a l ( numbercarbon ) )
allocate ( n o n a d j a c e n t l o c a t i o n ( numbercarbon , 3 ) )  240
allocate ( p a r t l y a d j a c e n t r a d i u s ( numbercarbon ) )
allocate ( p a r t l y a d j a c e n t l o c a t i o n ( numbercarbon , 3 ) )
i f ( memory . gt . numbercarbon ) then

245

! I can ’ t remember more than present atoms !
mem = numbercarbon

else

250

mem = memory
end i f
l o c a l r a d i u s = 0.0  255
t r u e l o c a l = 0.0
nonadjacentradius = 0.0
p a r t l y a d j a c e n t r a d i u s = 0.0
s m a l l e s t r a d i u s = 0.0
index = 0  260
i f ( t h i c k n e s s . ne . 0 . 0 ) then
! S t a t i s t i c a l c l a s s i f i c a t i o n o f the r a d i i p o s s i b l e now

  write ( ∗ , ’ ( " Second run. Statistics running..." ) ’ )  265

write ( ∗ , ∗)
write ( ∗ , ∗)
d i s t r i b u t i o n = 0.0

270

end i f
! Checking a l l p o s s i b l e combination o f 3 p o i n t s
do i = 1 , numbercarbon − 2

275

do j = i + 1 , numbercarbon − 1
do k = j + 1 , numbercarbon
280 ! c o l i n e a r i t y check − > ((b−a ) \ times ( c−a ) ) ˆ 2 = 0
f i r s t = alphacarbons ( i )%c o o r d i n a t e s
second = alphacarbons ( j )%c o o r d i n a t e s
t h i r d = alphacarbons ( k )%c o o r d i n a t e s

285

c r o s s e d = crossproduct ( ( second−f i r s t ) , ( third −f i r s t ) )
! r a d i i c a l c u l a t i o n
i f ( dot product ( crossed , c r o s s e d ) == 0) then  290

r2 = 0.0

else

295

! Radius formula when three p o i n t s are given
D( 1 ) = (sum ( ( alphacarbons ( i )%c o o r d i n a t e s − &
alphacarbons ( j )%c o o r d i n a t e s ) ∗∗2) )
D( 2 ) = (sum ( ( alphacarbons ( i )%c o o r d i n a t e s − &
300 alphacarbons ( k )%c o o r d i n a t e s ) ∗∗2) )
D( 3 ) = (sum ( ( alphacarbons ( j )%c o o r d i n a t e s − &
alphacarbons ( k )%c o o r d i n a t e s ) ∗∗2) )
A = sum(D∗∗2)
305 B = 2∗((D( 1 ) ∗D( 2 ) )+(D( 1 ) ∗D( 3 ) )+(D( 2 ) ∗D( 3 ) ) )
i f ( (A−B) == 0) then
C = −1 ! The p o i n t s are c o l i n e a r

310

else

C = ( − product (D) / (A − B) )
end i f  315
i f (C < 0) then
! This may a l s o occur i f the p o i n t s are almost c o l i n e a r ,
! so that the r a d i u s o f the c i r c l e i s q u i t e l a r g e  320

   write ( ∗ , ’ ( "    Problem with alpha carbons : " , i10 , ", " , i10 , ", "

, i 1 0 ) ’ ) i , j , k

   write ( ∗ , ’ ( "    Colinearity check dot product equals: " , f10 . 3 ) ’ )

dot product ( crossed , c r o s s e d )
write ( ∗ , ∗)

325

r2 = 0.0

else

r2 = s q r t (C)  330
i f ( t h i c k n e s s . ne . 0 . 0 ) then
! This i s the second run :
! C l a s s i f i c a t i o n o f the found r a d i u s  335
pos = f l o o r ( ( ( r2 / t h i c k n e s s ) − 1) ∗ 1 0 0 ) + 1
i f ( pos . gt . l i m i t ) then

340

! Found r a d i u s i s b i g g e r than the given l i m i t
pos = l i m i t
end i f  345
d i s t r i b u t i o n ( pos ) = d i s t r i b u t i o n ( pos ) + 1

end i f

350

combination ( 1 ) = i  ! This i s due to Fortran 9 0 c o m p a t i b i l i t y
combination ( 2 ) = j
combination ( 3 ) = k
do m = 1 , 3  355
! Completing f o r a l l three p o i n t s
l = combination (m)

360

do q = 1 , mem
i f ( r2 . l t . s m a l l e s t r a d i u s ( l , q ) ) then
! A s m a l l e r r a d i u s has been found  365
do r = mem, q+1, −1
! S h i f t i n g the r a d i u s vector f o r i n s e r t i o n o f the
! new r a d i u s  370
s m a l l e s t r a d i u s ( l , r ) = s m a l l e s t r a d i u s ( l , r −1)
s m a l l e s t l o c a t i o n ( l , r , 1 ) = s m a l l e s t l o c a t i o n ( l , r −1 ,1)
s m a l l e s t l o c a t i o n ( l , r , 2 ) = s m a l l e s t l o c a t i o n ( l , r −1 ,2)
s m a l l e s t l o c a t i o n ( l , r , 3 ) = s m a l l e s t l o c a t i o n ( l , r −1 ,3)  375
end do
! Saving the new r a d i u s

380

s m a l l e s t r a d i u s ( l , q ) = r2
s m a l l e s t l o c a t i o n ( l , q , 1 ) = i
s m a l l e s t l o c a t i o n ( l , q , 2 ) = j
s m a l l e s t l o c a t i o n ( l , q , 3 ) = k

385

exit

end i f
i f ( s m a l l e s t r a d i u s ( l , q ) == 0.0) then  390
! This i s the f i r s t r a d i u s f o r the given point
s m a l l e s t r a d i u s ( l , q ) = r2
s m a l l e s t l o c a t i o n ( l , q , 1 ) = i  395
s m a l l e s t l o c a t i o n ( l , q , 2 ) = j
s m a l l e s t l o c a t i o n ( l , q , 3 ) = k

exit

400

end i f
end do
end do  405

end i f
end i f

410

i f ( ( ( k−j ) + ( j−i ) ) == 2) then
! The found r a d i u s i s l o c a l
combination ( 1 ) = i  415
combination ( 2 ) = j
combination ( 3 ) = k
i f ( ( t r u e l o c a l ( j ) == 0) . or . ( t r u e l o c a l ( j ) > r2 ) ) then

420

! Only f o r the middle point j ( i < j < k ) t h i s i s a true
! l o c a l r a d i u s
t r u e l o c a l ( j ) = r2

425

end i f
do m = 1 , 3
430 l = combination (m)
! Completion f o r a l l three p o i n t s
i f ( ( l o c a l r a d i u s ( l ) == 0) . or . ( l o c a l r a d i u s ( l ) > r2 ) ) then

435

l o c a l r a d i u s ( l ) = r2
f i r s t l o c a l ( l ) = i  ! I n d i c a t e s the f i r s t point
end i f

440

end do
e l s e i f ( ( ( k−j ) . gt . 1 ) . and . ( ( j−i ) . gt . 1 ) ) then
! None o f the atoms are neighbors  445
combination ( 1 ) = i
combination ( 2 ) = j
combination ( 3 ) = k

450

do m = 1 , 3
l = combination (m)
i f ( ( nonadjacentradius ( l ) == 0) . or . ( nonadjacentradius ( l ) > r2 ) )  455

then

nonadjacentradius ( l ) = r2
n o n a d j a c e n t l o c a t i o n ( l , 1 ) = i
n o n a d j a c e n t l o c a t i o n ( l , 2 ) = j
n o n a d j a c e n t l o c a t i o n ( l , 3 ) = k  460
end i f
end do

465

else

! Two o f the three atoms are d i r e c t neighbors
combination ( 1 ) = i  470
combination ( 2 ) = j
combination ( 3 ) = k
do m = 1 , 3

475

l = combination (m)
i f ( ( ( ( k−j ) . ne . 1 ) . and . ( ( j−i ) . ne . 1 ) ) . or . ( ( ( k−j ) == 1) . and
. ( ( j−i ) == 1) ) ) then
write ( ∗ , ’ ( "Error !!!" ) ’ )  480
end i f
i f ( ( p a r t l y a d j a c e n t r a d i u s ( l ) == 0) . or . ( p a r t l y a d j a c e n t r a d i u s ( l )
> r2 ) ) then

485

p a r t l y a d j a c e n t r a d i u s ( l ) = r2
p a r t l y a d j a c e n t l o c a t i o n ( l , 1 ) = i
p a r t l y a d j a c e n t l o c a t i o n ( l , 2 ) = j
p a r t l y a d j a c e n t l o c a t i o n ( l , 3 ) = k

490

end i f
end do
end i f  495
index = index + 1
end do

500

end do
end do

505

i f ( talkmode ) then

write ( ∗ , ’ ( " " , i10 , " circles have been calculated" ) ’ ) index

write ( ∗ , ∗)

510

end i f
! F i n a l i z i n g
allocate ( g l o b a l r a d i u s ( numbercarbon ) )  515
allocate ( g l o b a l l o c a t i o n ( numbercarbon , 3 ) )
allocate ( n o n l o c a l r a d i u s ( numbercarbon ) )
do i = 1 , numbercarbon

520

! Which c l a s s o f r a d i i provides the s m a l l e s t r a d i u s ?

r a d i i ( 1 ) = l o c a l r a d i u s ( i )
r a d i i ( 2 ) = nonadjacentradius ( i )
r a d i i ( 3 ) = p a r t l y a d j a c e n t r a d i u s ( i )  525
s e l e c t o r = minloc ( r a d i i , mask = r a d i i . ne . 0 . 0 )
l o c a t i o n = s e l e c t o r ( 1 )
! Global r a d i u s l i s t s the s m a l l e s t c i r c l e f o r each point  530
g l o b a l r a d i u s ( i ) = r a d i i ( l o c a t i o n )
select case ( l o c a t i o n )
! Indexes o f the p o i n t s r e a l i z i n g the s m a l l e s t c i r c l e  535
case ( 1 )
g l o b a l l o c a t i o n ( i , 1 ) = f i r s t l o c a l ( i )
g l o b a l l o c a t i o n ( i , 2 ) = f i r s t l o c a l ( i ) + 1
g l o b a l l o c a t i o n ( i , 3 ) = f i r s t l o c a l ( i ) + 2  540
case ( 2 )
g l o b a l l o c a t i o n ( i , : ) = n o n a d j a c e n t l o c a t i o n ( i , : )

545

case ( 3 )
g l o b a l l o c a t i o n ( i , : ) = p a r t l y a d j a c e n t l o c a t i o n ( i , : )
end select  550
! Doing the same but by i g n o r i n g the l o c a l r a d i i
r a d i i ( 1 ) = 0.0
r a d i i ( 2 ) = nonadjacentradius ( i )  555
r a d i i ( 3 ) = p a r t l y a d j a c e n t r a d i u s ( i )
s e l e c t o r = minloc ( r a d i i , mask = r a d i i . ne . 0 . 0 )
l o c a t i o n = s e l e c t o r ( 1 )

560

n o n l o c a l r a d i u s ( i ) = r a d i i ( l o c a t i o n )
end do
allocate ( f ( numbercarbon ) )  565
allocate ( g ( numbercarbon ) )
f = 0.0
g = 0.0

570

do i = 2 , numbercarbon − 1
! C a l c u l a t i o n o f the f and g f a c t o r
f ( i ) = n o n l o c a l r a d i u s ( i ) / t r u e l o c a l ( i )  575
g ( i ) = g l o b a l r a d i u s ( i ) / t r u e l o c a l ( i )
end do
deallocate ( nonadjacentradius , n o n a d j a c e n t l o c a t i o n )  580

deallocate ( p a r t l y a d j a c e n t r a d i u s , p a r t l y a d j a c e n t l o c a t i o n )
deallocate ( n o n l o c a l r a d i u s , f i r s t l o c a l , l o c a l r a d i u s )
end subroutine c i r c l e s

585

! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
subroutine spheres ( alphacarbons , f , g , l o c a l r a d i u s , g l o b a l r a d i u s , &
g l o b a l l o c a t i o n , s m a l l e s t r a d i u s , s m a l l e s t l o c a t i o n , &
d i s t r i b u t i o n , t h i c k n e s s )  590
! f i t s spheres through the macromolecule
use var
use mathe
use parameters  595
implicit none
type ( carbons ) , pointer , dimension ( : ) : : alphacarbons

600

real  : : A, B, maxradius , rquadrat
real  : : r3 , thickness , detn , detz
real  : : prevtime , time
real , dimension ( 4 , 4 )  : : D, E , H
real , dimension ( 3 ) : : C, f i r s t , second , third , fourth 605
real , dimension ( 2 )  : : r a d i i
real , dimension ( 4 )  : : minima
real , pointer , dimension ( : )  : : f , g , l o c a l r a d i u s , n o n l o c a l r a d i u s
real , pointer , dimension ( : )  : : g l o b a l r a d i u s
real , pointer , dimension ( : , : ) : : s m a l l e s t r a d i u s 610
real , pointer , dimension ( : , : )  : : coplanarity , N, Z
integer  : : numbercarbon , i , j , k , l , m, o
integer  : : q , r , s , index , e r r o r
integer : : l o c a t i o n , pos , mem 615
integer , pointer , dimension ( : )  : : f i r s t l o c a l
integer , pointer , dimension ( : , : )  : : n o n l o c a l l o c a t i o n , g l o b a l l o c a t i o n
integer , pointer , dimension ( : , : , : )  : : s m a l l e s t l o c a t i o n
integer , dimension ( 4 )  : : combination
integer , dimension ( 1 ) : : s e l e c t o r 620
integer ( 4 ) , dimension ( l i m i t )  : : d i s t r i b u t i o n
numbercarbon = size ( alphacarbons )
allocate ( l o c a l r a d i u s ( numbercarbon ) )  625
allocate ( f i r s t l o c a l ( numbercarbon ) )
allocate ( n o n l o c a l r a d i u s ( numbercarbon ) )
allocate ( n o n l o c a l l o c a t i o n ( numbercarbon , 4 ) )
i f ( memory . gt . numbercarbon ) then  630
! I can ’ t remember more than present atoms !
mem = numbercarbon

635 else

mem = memory
end i f

640

index = 0
e r r o r = 0
l o c a l r a d i u s = 0.0
n o n l o c a l r a d i u s = 0.0
s m a l l e s t r a d i u s = 0.0  645
allocate ( c o p l a n a r i t y ( 4 , 4 ) )
allocate (N( 4 , 4 ) )
allocate (Z( 5 , 5 ) )

650

i f ( t h i c k n e s s . ne . 0 . 0 ) then
! S t a t i s t i c a l c l a s s i f i c a t i o n o f the r a d i i p o s s i b l e now

write ( ∗ , ’ ( " Second run. Statistics running..." ) ’ )

write ( ∗ , ∗)  655
write ( ∗ , ∗)
d i s t r i b u t i o n = 0.0
end i f  660
! Checking a l l p o s s i b l e combination o f 4 p o i n t s
do i = 1 , numbercarbon − 3
do j = i + 1 , numbercarbon − 2  665
do k = j + 1 , numbercarbon − 1
do l = k + 1 , numbercarbon

670

! c o p l a n a r i t y check −> (d−a ) . ( ( b−a ) \ times ( c−a ) ) = 0
c o p l a n a r i t y = reshape ( ( / alphacarbons ( i )%coordinates , 1 . 0 , &
alphacarbons ( j )%coordinates , 1 . 0 , &
675 alphacarbons ( k )%coordinates , 1 . 0 , &
alphacarbons ( l )%coordinates , 1 . 0 / ) , ( / 4 , 4 / )

)

! r a d i i c a l c u l a t i o n
D( 1 : 4 , 1 : 4 ) = 0.0  680
D( 1 , 2 ) = (sum ( ( alphacarbons ( i )%c o o r d i n a t e s − &
alphacarbons ( j )%c o o r d i n a t e s ) ∗∗2) )
D( 1 , 3 ) = (sum ( ( alphacarbons ( i )%c o o r d i n a t e s − &
alphacarbons ( k )%c o o r d i n a t e s ) ∗∗2) )
685 D( 1 , 4 ) = (sum ( ( alphacarbons ( i )%c o o r d i n a t e s − &
alphacarbons ( l )%c o o r d i n a t e s ) ∗∗2) )
D( 2 , 3 ) = (sum ( ( alphacarbons ( j )%c o o r d i n a t e s − &
alphacarbons ( k )%c o o r d i n a t e s ) ∗∗2) )
D( 2 , 4 ) = (sum ( ( alphacarbons ( j )%c o o r d i n a t e s − &
690 alphacarbons ( l )%c o o r d i n a t e s ) ∗∗2) )
D( 3 , 4 ) = (sum ( ( alphacarbons ( k )%c o o r d i n a t e s − &
alphacarbons ( l )%c o o r d i n a t e s ) ∗∗2) )
N( 1 : 4 , 1 : 4 ) = 0.0
Z ( 1 : 5 , 1 : 5 ) = 0.0

695

N = D + transpose (D)
Z ( 1 : 4 , 1 : 4 ) = N

Z ( 1 : 4 , 5 ) = 1.0
Z ( 5 , 1 : 4 ) = 1.0

700

detn = determinant (N)
detz = determinant (Z)
i f ( detz == 0) then

705

! The p o i n t s are coplanar
rquadrat = −1

else

710

rquadrat = ( − detn / ( 2 ∗ detz ) )
end i f
i f ( ( abs ( determinant ( c o p l a n a r i t y ) ) == 0) .OR. ( rquadrat <0) ) then  715
! This may a l s o occur i f the p o i n t s are almost coplanar , so that
the r a d i u s o f
! the sphere i s q u i t e l a r g e
i f ( talkmode ) then  720

  write ( ∗ , ’ ( " Problem with alpha carbons : " , i10 , ", " , i10 , ", "

, i10 , ", " , i 1 0 ) ’ ) i , j , k , l

  write ( ∗ , ’ ( " Determinant of Z equals: " , f10 . 3 ) ’ ) detz
  write ( ∗ , ’ ( " Coplanarity check determinant equals: " , f10 . 3 ) ’ )

determinant ( c o p l a n a r i t y )
write ( ∗ , ∗)  725
end i f
r3 = 0.0

730

e r r o r = e r r o r + 1

else

r3 = s q r t ( rquadrat )  735
i f ( t h i c k n e s s . ne . 0 . 0 ) then
! This i s the second run :
! C l a s s i f i c a t i o n o f the found r a d i u s  740
pos = f l o o r ( ( ( r3 / t h i c k n e s s ) − 1) ∗ 1 0 0 ) + 1
i f ( pos . gt . l i m i t ) then

745

! Found r a d i u s i s b i g g e r than the given l i m i t
pos = l i m i t
end i f  750
d i s t r i b u t i o n ( pos ) = d i s t r i b u t i o n ( pos ) + 1

end i f

755

combination ( 1 ) = i ! This i s due to Fortran 9 0 c o m p a t i b i l i t y
combination ( 2 ) = j
combination ( 3 ) = k
combination ( 4 ) = l

760

do m = 1 , 4
! Completing f o r a l l four p o i n t s
765 s = combination (m)
do q = 1 , mem
i f ( r3 . l t . s m a l l e s t r a d i u s ( s , q ) ) then

770

! A s m a l l e r r a d i u s has been found
do r = mem, q+1, −1
! S h i f t i n g the r a d i u s vector f o r i n s e r t i o n o f the  775
! new r a d i u s
s m a l l e s t r a d i u s ( s , r ) = s m a l l e s t r a d i u s ( s , r −1)
s m a l l e s t l o c a t i o n ( s , r , 1 ) = s m a l l e s t l o c a t i o n ( s , r −1 ,1)
s m a l l e s t l o c a t i o n ( s , r , 2 ) = s m a l l e s t l o c a t i o n ( s , r −1 ,2)  780
s m a l l e s t l o c a t i o n ( s , r , 3 ) = s m a l l e s t l o c a t i o n ( s , r −1 ,3)
s m a l l e s t l o c a t i o n ( s , r , 4 ) = s m a l l e s t l o c a t i o n ( s , r −1 ,4)
end do

785

! Saving the new r a d i u s
s m a l l e s t r a d i u s ( s , q ) = r3
s m a l l e s t l o c a t i o n ( s , q , 1 ) = i
s m a l l e s t l o c a t i o n ( s , q , 2 ) = j  790
s m a l l e s t l o c a t i o n ( s , q , 3 ) = k
s m a l l e s t l o c a t i o n ( s , q , 4 ) = l

exit

795

end i f
i f ( s m a l l e s t r a d i u s ( s , q ) == 0.0) then
! This i s the f i r s t r a d i u s f o r the given point  800
s m a l l e s t r a d i u s ( s , q ) = r3
s m a l l e s t l o c a t i o n ( s , q , 1 ) = i
s m a l l e s t l o c a t i o n ( s , q , 2 ) = j
s m a l l e s t l o c a t i o n ( s , q , 3 ) = k  805
s m a l l e s t l o c a t i o n ( s , q , 4 ) = l

exit

end i f  810
end do

end do

815

end i f
i f ( ( ( l−k ) + ( k−j ) + ( j−i ) ) == 3) then
! The found r a d i u s i s l o c a l  820
i f ( ( l o c a l r a d i u s ( j ) > r3 ) . or . ( l o c a l r a d i u s ( j ) == 0.0) ) then
l o c a l r a d i u s ( j ) = r3

825

end i f
! The l o c a l r a d i u s i s a t t r i b u t e d to the both inner p o i n t s
l o c a l r a d i u s ( k ) = r3

830

else

! The r a d i u s i s not an l o c a l one
combination ( 1 ) = i  835
combination ( 2 ) = j
combination ( 3 ) = k
combination ( 4 ) = l
do m = 1 , 4  840
o = combination (m)
i f ( ( n o n l o c a l r a d i u s ( o ) == 0) . or . ( n o n l o c a l r a d i u s ( o ) > r3 ) ) then

845

n o n l o c a l r a d i u s ( o ) = r3
n o n l o c a l l o c a t i o n ( o , 1 ) = i
n o n l o c a l l o c a t i o n ( o , 2 ) = j
n o n l o c a l l o c a t i o n ( o , 3 ) = k
n o n l o c a l l o c a t i o n ( o , 4 ) = l  850
end i f
end do

855

end i f
index = index + 1
end do  860
end do
end do

865

end do

i f ( talkmode ) then

870

write ( ∗ , ’ ( " " , i10 , " spheres have been calculated." ) ’ ) index

write ( ∗ , ’ ( " For " , i10 , " spheres the determination of its radius was
    not possible." ) ’ ) e r r o r

write ( ∗ , ∗)
end i f  875
deallocate ( c o p l a n a r i t y )
deallocate (N)
deallocate (Z)

880

allocate ( g l o b a l r a d i u s ( numbercarbon ) )
allocate ( g l o b a l l o c a t i o n ( numbercarbon , 4 ) )
! F i n a l i z i n g

885

do i = 1 , numbercarbon
! Which c l a s s o f r a d i i provides the s m a l l e s t r a d i u s ?
r a d i i ( 1 ) = l o c a l r a d i u s ( i )  890
r a d i i ( 2 ) = n o n l o c a l r a d i u s ( i )
s e l e c t o r = minloc ( r a d i i , mask = r a d i i . ne . 0 . 0 )
l o c a t i o n = s e l e c t o r ( 1 )

895

! Global r a d i u s l i s t s the s m a l l e s t c i r c l e f o r each point
g l o b a l r a d i u s ( i ) = r a d i i ( l o c a t i o n )
select case ( l o c a t i o n )

900

! Indexes o f the p o i n t s r e a l i z i n g the s m a l l e s t c i r c l e
case ( 1 )
g l o b a l l o c a t i o n ( i , 1 ) = i − 1
g l o b a l l o c a t i o n ( i , 2 ) = i  905
g l o b a l l o c a t i o n ( i , 3 ) = i + 1
g l o b a l l o c a t i o n ( i , 4 ) = i + 2
case ( 2 )

910

g l o b a l l o c a t i o n ( i , : ) = n o n l o c a l l o c a t i o n ( i , : )
end select
end do  915
allocate ( f ( numbercarbon ) )
allocate ( g ( numbercarbon ) )
f = 0.0  920
g = 0.0
do i = 2 , numbercarbon − 2
! C a l c u l a t i o n o f the f and g f a c t o r  925
f ( i ) = n o n l o c a l r a d i u s ( i ) / l o c a l r a d i u s ( i )
g ( i ) = g l o b a l r a d i u s ( i ) / l o c a l r a d i u s ( i )

end do  930
deallocate ( n o n l o c a l r a d i u s , n o n l o c a l l o c a t i o n )
deallocate ( f i r s t l o c a l )
end subroutine spheres  935
end module d i s t a n c e
! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗

940

module loading
implicit none

945 contains

subroutine import ( i n p u t f i l e , alphacarbons , atomtype , r e s o l u t i o n )
! Loads the preprocessed input data f i l e s

950

use var
implicit none
real  : : r e s o l u t i o n  955
integer  : : i , io , numbercarbons , temp
character (LEN=18) , intent ( in )  : : i n p u t f i l e
! Compatibility with Fortran 90  960
character (LEN=4)  : : atomtype
character (LEN=12)  : : dummy
character (LEN=13)  : : dummy2
965 character (LEN=17)  : : dummy3
character (LEN=35)  : : t e s t
type ( carbons ) , pointer , dimension ( : ) : : alphacarbons
970 OPEN (UNIT=2 , FILE=i n p u t f i l e , FORM=’FORMATTED’ , STATUS=’OLD’ )
! Reading o f the header l i n e and the r e s o l u t i o n value
rewind 2

975

read ( 2 , ’ ( a35 ) ’ ) t e s t
i f ( t e s t . ne . ’ Atom c a r t e s i a n c o o r d i n a t e s from PDB’ ) then

  write ( ∗ , ’ ( "This is not an atom coordinates file" ) ’ )  980

return

end i f
read ( 2 , ∗)  985
read ( 2 , ∗)
read ( 2 , ’ ( a13 , a3 ) ’ ) dummy2 , atomtype

read ( 2 , ’ ( a17 , f 4 . 0 ) ’ ) dummy3 , r e s o l u t i o n  990
read ( 2 , ∗)
read ( 2 , ∗)
read ( 2 , ∗)

995

numbercarbons = 0

do

! Counting the numbers o f atoms  1000
read ( 2 , ’ ( a12 ) ’ , IOSTAT = i o ) dummy
i f ( i o == 0) then

1005

numbercarbons = numbercarbons + 1

else

1010 exit

end i f
end do

1015

allocate ( alphacarbons ( numbercarbons ) )
rewind 2
do i = 1 , 8  1020
read ( 2 , ∗)
end do

1025

do i = 1 , numbercarbons
! Reading the c o o r d i n a t e s
read ( 2 , ’ ( i5 , 3 f10 . 3 , i9 , i6 , a1 ) ’ ) temp , &  1030
alphacarbons ( i )%c o o r d i n a t e s ( 1 ) , &
alphacarbons ( i )%c o o r d i n a t e s ( 2 ) , &
alphacarbons ( i )%c o o r d i n a t e s ( 3 ) , &
alphacarbons ( i )%residue , temp , &
1035 alphacarbons ( i )%c h a i n i d e n t
end do
CLOSE (UNIT=2 , STATUS=’KEEP’ )
end subroutine import  1040
end module loading
! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗

1045

program c a l c u l a t o r

! Main program
use loading  1050
use var
use d i s t a n c e
use parameters
implicit none  1055
character (LEN=12)  : : export , dummy
character (LEN=7)  : : e x p o r t f o l d e r , machinefolder
character (LEN=5)  : : subfolder , humanfolder
1060 character (LEN=6)  : : i m p o r t f o l d e r
character (LEN=12) , allocatable , dimension ( : ) : : filename
character , dimension ( : ) , pointer  : : name
character (LEN=4)  : : extension , atomtype
real  : : r e s o l u t i o n , t h i c k n e s s  1065
real ( 8 )  : : i n i t i a l , f i n a l
real , pointer , dimension ( : )  : : f , g , l o c a l r a d i u s
real , pointer , dimension ( : )  : : g l o b a l r a d i u s
real , pointer , dimension ( : , : )  : : a l p h a c i r c l e s , s m a l l e s t r a d i u s

1070

type ( carbons ) , pointer , dimension ( : )  : : alphacarbons
type ( carbons ) , allocatable , dimension ( : )  : : a l l a l p h a c a r b o n s , oldcarbons
integer  : : i , j , k , numberfiles , i o
integer : : position , length , d i s l o c a t i o n 1075
integer  : : a , b , c , d , what , s o u r c e f l a g
integer  : : f i t f l a g , numbercarbon , mem
integer , allocatable , dimension ( : , : )  : : comb
integer , pointer , dimension ( : , : )  : : g l o b a l l o c a t i o n
integer , pointer , dimension ( : , : , : ) : : s m a l l e s t l o c a t i o n 1080
integer ( 4 ) , dimension ( l i m i t )  : : d i s t r i b u t i o n
logical , dimension ( : ) , allocatable  : : g o o d f i l e s
e x p o r t f o l d e r  = ’ r e s u l t s ’  1085
! name o f the s u b d i r e c t o r y where the r e s u l t s are s t o r e d to
! Adjust the above f i e l d length o f the s t r i n g ! ! !
i m p o r t f o l d e r  = ’ source ’
! name o f the s u b f o l d e r f o r ASCII export f i l e s  1090
! Adjust the above f i e l d length o f the s t r i n g ! ! !
humanfolder  = ’human ’
machinefolder  = ’ machine ’

1095

i f ( talkmode ) then
write ( ∗ , ∗)

  write ( ∗ , ’ ( " Computational physics project" ) ’ )
  write ( ∗ , ’ ( " Christian Schlatter, student, physics, 4th year" ) ’ )  1100

write ( ∗ , ∗) r e l e a s e i n f o
write ( ∗ , ’ ( "

       ****************************************************************" ) ’ )

write ( ∗ , ∗)
end i f  1105

! s o u r c e f l a g = 1  ! > TEMPORARY
! a l l o c a t e ( filename ( 1 ) )  ! > TEMPORARY
! filename ( 1 ) = ’1 acc . hel ’  ! > TEMPORARY
! n u m b e r f i l e s = 1  ! > TEMPORARY  1110

write ( ∗ , ’ ( " Global curvature analyser" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " How many proteins would you like to process ?" ) ’ )

write ( ∗ , ∗)  1115

write ( ∗ , ’ ( " -> Type 1 for a single input file or" ) ’ )
write ( ∗ , ’ ( " -> Type 2 for a list of files provided in Source.ini" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " Please choose your option now : " ) ’ , advance = ’ no ’ )

1120

s o u r c e f l a g = 2  ! < TEMPORARY
! read ( ∗ , ∗) s o u r c e f l a g  ! > TEMPORARY
write ( ∗ , ∗)  1125
i f ( s o u r c e f l a g == 2) then ! Batch mode using Source . i n i

  write ( ∗ , ’ ( " Opening Source.ini..." ) ’ )

write ( ∗ , ∗)  1130
open (UNIT=1 , FILE=’ Source . i n i ’ , FORM=’FORMATTED’ , STATUS=’OLD’ )
rewind 1

1135

n u m b e r f i l e s = 0

do

! Determination o f the number o f input f i l e s  1140
read ( 1 , ’ ( a12 ) ’ , IOSTAT = i o ) dummy
i f ( i o == 0) then

1145

n u m b e r f i l e s = n u m b e r f i l e s + 1

else

1150 exit

end i f
end do

1155

allocate ( filename ( n u m b e r f i l e s ) )
rewind 1
do i = 1 , n u m b e r f i l e s  1160
! Determination o f the input f i l e names
read ( 1 , ’ ( a12 ) ’ ) filename ( i )
Projet de simulation numérique 2002  Christian Schlatter

write ( ∗ , ∗) filename ( i )  1165
end do
close (UNIT=1 , STATUS=’KEEP’ )

1170

write ( ∗ , ∗)

write ( ∗ , ’ ( " You have chosen " , i3 , " files to process" ) ’ ) n u m b e r f i l e s

write ( ∗ , ∗)

1175 else

allocate ( filename ( 1 ) )
! User provided input f i l e name

1180

write ( ∗ , ’ ( " Please provide the filename of the import file: " ) ’ , advance

= ’ no ’ )
read ( ∗ , ∗) filename ( 1 )
write ( ∗ , ∗)  1185
n u m b e r f i l e s = 1

endif

1190

write ( ∗ , ’ ( " What would you like to fit into the protein ?" ) ’ )

write ( ∗ , ∗)

write ( ∗ , ’ ( " -> Type 1 for fitting circles only" ) ’ )
write ( ∗ , ’ ( " -> Type 2 for fitting both - circles and spheres" ) ’ )

write ( ∗ , ∗)  1195

write ( ∗ , ’ ( " Please choose your option now : " ) ’ , advance = ’ no ’ )

f i t f l a g = 2  ! < TEMPORARY
! read ( ∗ , ∗) f i t f l a g  ! > TEMPORARY  1200
write ( ∗ , ∗)
allocate ( g o o d f i l e s ( n u m b e r f i l e s ) )
g o o d f i l e s = . true .

1205

do i = 1 , n u m b e r f i l e s
i f ( talkmode ) then
write ( ∗ , ∗)  1210

write ( ∗ , ’ ( " Extracting structure information from file " , a12 ) ’ )

filename ( i )
write ( ∗ , ’ ( "

     *****************************************************************" ) ’ )

write ( ∗ , ∗)
end i f  1215
write ( ∗ , ’ ( " Processing now file " , a12 ) ’ ) filename ( i )
write ( ∗ , ∗)
nullify ( alphacarbons )  1220

i f ( a s s o c i a t e d ( alphacarbons ) ) then
deallocate ( alphacarbons )

1225

end i f
! Reading the input f i l e
c a l l import ( i m p o r t f o l d e r // ’/ ’// filename ( i ) , alphacarbons , atomtype ,  1230
r e s o l u t i o n )
numbercarbon = size ( alphacarbons )
i f ( memory . gt . numbercarbon ) then

1235

mem = numbercarbon

else

mem = memory  1240
end i f
allocate ( s m a l l e s t r a d i u s ( numbercarbon , mem) )
allocate ( s m a l l e s t l o c a t i o n ( numbercarbon ,mem, 3 ) )  1245
i f ( talkmode ) then

  write ( ∗ , ’ ( " Fitting circles through the alpha carbons" ) ’ )

write ( ∗ , ’ ( "  1250

       *****************************************************************" ) ’ )

write ( ∗ , ∗)
end i f
t h i c k n e s s = 0.0  1255
! F i t t i n g c i r c l e s through three atoms
c a l l c i r c l e s ( alphacarbons , f , g , l o c a l r a d i u s , g l o b a l r a d i u s , &
g l o b a l l o c a t i o n , s m a l l e s t r a d i u s , s m a l l e s t l o c a t i o n , &  1260
d i s t r i b u t i o n , t h i c k n e s s )
! C i r c u l a r t h i c k n e s s o f the macromolecule
t h i c k n e s s = minval ( g l o b a l r a d i u s )  1265
i f ( secondrun ) then
! Doing a second run f o r s t a t i s t i c s on the s m a l l e s t c i r c l e s

1270

deallocate ( l o c a l r a d i u s , g l o b a l r a d i u s , f , g , g l o b a l l o c a t i o n )
c a l l c i r c l e s ( alphacarbons , f , g , l o c a l r a d i u s , g l o b a l r a d i u s , &
g l o b a l l o c a t i o n , s m a l l e s t r a d i u s , s m a l l e s t l o c a t i o n , &
d i s t r i b u t i o n , t h i c k n e s s )  1275
end i f

position = scan ( filename ( i ) , ’ . ’ )

1280

i f ( filename ( i ) ( position +1: position+1) == ’h ’ ) then
! An alpha h e l i x input f i l e has been used
s u b f o l d e r = ’ h e l i x ’

1285

else

! An backbone input f i l e has been used
s u b f o l d e r = ’ whole ’

1290

end i f
! Saving the c i r c l e s to disk
filename ( i ) ( position +1: position+1) = ’ c ’  1295
filename ( i ) ( position +2: position+2) = ’ i ’
filename ( i ) ( position +3: position+3) = ’ r ’
open (UNIT=3 , FILE=e x p o r t f o l d e r // ’/ ’// s u b f o l d e r // ’/ ’// humanfolder // ’/ ’//
filename ( i ) , FORM=’FORMATTED’ , STATUS=’REPLACE’ )
1300 open (UNIT=4 , FILE=e x p o r t f o l d e r // ’/ ’// s u b f o l d e r // ’/ ’// machinefolder // ’/ ’//
filename ( i ) , FORM=’FORMATTED’ , STATUS=’REPLACE’ )
! Making the output f i l e header

write ( 3 , ’ ( "Circling of proteins" ) ’ )

write ( 4 , ’ ( "Circles" ) ’ )  1305

write ( 3 , ’ ( "C-alpha-circle project. File created using Calculator." ) ’ )

write ( 3 , ’ ( a25 ) ’ , advance=’no ’ ) r e l e a s e i n f o
write ( 4 , ’ ( a25 ) ’ ) r e l e a s e i n f o

write ( 3 , ’ ( "\begin{matrix}c \\ ◯ \end{matrix} 2002 EPFL Lausanne" ) ’ )

write ( 3 , ’ ( "Atom type : " , a4 ) ’ ) atomtype  1310
write ( 4 , ’ ( a4 ) ’ ) atomtype

write ( 3 , ’ ( "ResolutionÅ [] : " , f 4 . 2 ) ’ ) r e s o l u t i o n

write ( 4 , ’ ( f 4 . 2 ) ’ ) r e s o l u t i o n

1315

write ( 3 , ’ ( "Circular thicknessÅ [] : " , f 4 . 2 ) ’ ) t h i c k n e s s

write ( 4 , ’ ( f 4 . 2 ) ’ ) t h i c k n e s s

write ( 3 , ’ ( "Thickness combination : " , i3 , ", " , i3 , " and " , i 3 ) ’ ) &

g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,1) , &  1320
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,2) , &
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,3)
write ( 4 , ’ ( 3 i 4 ) ’ ) g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,1) , &
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,2) , &  1325
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,3)
write ( 3 , ’ ( "Memory : " , i 3 ) ’ ) mem
write ( 4 , ’ ( i 3 ) ’ ) mem

1330

i f ( secondrun ) then
! Saving s t a t i s t i c s r e s u l t s to disk

filename ( i ) ( position +1: position+1) = ’ c ’  1335
filename ( i ) ( position +2: position+2) = ’d ’
filename ( i ) ( position +3: position+3) = ’ t ’
open (UNIT=5 , FILE=e x p o r t f o l d e r // ’/ ’// s u b f o l d e r // ’/ ’// machinefolder // ’/ ’//
filename ( i ) , FORM=’FORMATTED’ , STATUS=’REPLACE’ )

1340

do j = 1 , l i m i t
! D i s t r i b u t i o n o f the s m a l l e s t r a d i i
write ( 5 , ’ ( f10 . 2 , i 1 0 ) ’ ) ( ( j / 1 0 0 . 0 ) + 1 . 0 ) , d i s t r i b u t i o n ( j )  1345
end do
close (UNIT=5 , STATUS=’KEEP’ )

1350

end i f
write ( 3 , ∗)
write ( 4 , ∗)
write ( 3 , ’ ( "Entry      Local    f-ratio   Dis" ) ’ , advance = ’ no ’ )  1355
do j = 1 , mem
write ( 3 , ’ ( "     " , i 2 ) ’ , advance = ’ no ’ ) j
write ( 3 , ’ ( ". global radius " ) ’ , advance = ’ no ’ )  1360
end do
write ( 3 , ’ ( " Residue" ) ’ )
write ( 3 , ∗)  1365
! Writing data to f i l e
do j = 1 , ( size ( alphacarbons ) )

1370

a = minval ( g l o b a l l o c a t i o n ( j , : ) )
c = maxval ( g l o b a l l o c a t i o n ( j , : ) )
b = minval ( g l o b a l l o c a t i o n ( j , : ) , MASK = g l o b a l l o c a t i o n ( j , : ) .NE. a )
! Spacing o f the indexes  1375
d i s l o c a t i o n = ( c − b) + ( b − a )
write ( 3 , ’ ( i5 , 2 f10 . 3 , i 5 ) ’ , advance = ’ no ’ ) j , &
l o c a l r a d i u s ( j ) / thickness , f ( j ) , d i s l o c a t i o n  1380
write ( 4 , ’ ( i5 , 2 f10 . 3 , i 5 ) ’ , advance = ’ no ’ ) j , &
l o c a l r a d i u s ( j ) / thickness , f ( j ) , d i s l o c a t i o n
do k = 1 , mem  1385
write ( 3 , ’ ( " (" , f 5 . 3 , ";" , 3 i4 , ") " ) ’ , advance = ’ no ’ ) &
s m a l l e s t r a d i u s ( j , k ) / thickness , &
s m a l l e s t l o c a t i o n ( j , k , 1 ) , &
s m a l l e s t l o c a t i o n ( j , k , 2 ) , &  1390
s m a l l e s t l o c a t i o n ( j , k , 3 )

write ( 4 , ’ ( f 7 . 3 ) ’ , advance = ’ no ’ ) s m a l l e s t r a d i u s ( j , k ) / t h i c k n e s s
end do  1395
write ( 3 , ’ ( i 8 ) ’ ) alphacarbons ( j )%r e s i d u e
write ( 4 , ’ ( i 8 ) ’ ) alphacarbons ( j )%r e s i d u e
end do  1400

write ( ∗ , ’ ( " -> This protein has a thickness of circularily " , f 5 . 3 , "

Angstroem." ) ’ ) t h i c k n e s s
write ( ∗ , ∗)
1405 close (UNIT=3 , STATUS=’KEEP’ )
close (UNIT=4 , STATUS=’KEEP’ )
deallocate ( l o c a l r a d i u s , g l o b a l r a d i u s , f , g , g l o b a l l o c a t i o n , s m a l l e s t r a d i u s )
deallocate ( s m a l l e s t l o c a t i o n )

1410

i f ( f i t f l a g == 2) then
! F i t t i n g spheres through 4 atoms o f the macromolecule
allocate ( s m a l l e s t r a d i u s ( numbercarbon , mem) )  1415
allocate ( s m a l l e s t l o c a t i o n ( numbercarbon ,mem, 4 ) )
i f ( talkmode ) then

write ( ∗ , ’ ( " Fitting spheres through the alpha carbons" ) ’ )  1420

write ( ∗ , ’ ( "

     *****************************************************************" )

’ )

write ( ∗ , ∗)
end i f

1425

t h i c k n e s s = 0.0
c a l l spheres ( alphacarbons , f , g , l o c a l r a d i u s , g l o b a l r a d i u s , &
g l o b a l l o c a t i o n , s m a l l e s t r a d i u s , s m a l l e s t l o c a t i o n , &
d i s t r i b u t i o n , t h i c k n e s s )  1430
! S p h e r i c a l t h i c k n e s s o f the macromolecule

write ( ∗ , ’ ( "First run successfully completed..." ) ’ )

1435

t h i c k n e s s = minval ( g l o b a l r a d i u s )
i f ( secondrun ) then
! Doing a second run f o r s t a t i s t i c s on the s m a l l e s t spheres  1440
deallocate ( l o c a l r a d i u s , g l o b a l r a d i u s , f , g , g l o b a l l o c a t i o n )
c a l l spheres ( alphacarbons , f , g , l o c a l r a d i u s , g l o b a l r a d i u s , &
g l o b a l l o c a t i o n , s m a l l e s t r a d i u s , s m a l l e s t l o c a t i o n , &  1445
d i s t r i b u t i o n , t h i c k n e s s )
end i f

position = scan ( filename ( i ) , ’ . ’ )  1450
! Saving the spheres to disk
filename ( i ) ( position +1: position+1) = ’ s ’
filename ( i ) ( position +2: position+2) = ’p ’  1455
filename ( i ) ( position +3: position+3) = ’h ’
open (UNIT=3 , FILE=e x p o r t f o l d e r // ’/ ’// s u b f o l d e r // ’/ ’// humanfolder // ’/ ’//
filename ( i ) , FORM=’FORMATTED’ , STATUS=’REPLACE’ )
open (UNIT=4 , FILE=e x p o r t f o l d e r // ’/ ’// s u b f o l d e r // ’/ ’// machinefolder // ’/ ’//
filename ( i ) , FORM=’FORMATTED’ , STATUS=’REPLACE’ )

1460

! Making the output f i l e header

write ( 3 , ’ ( "Sphering of proteins" ) ’ )

write ( 4 , ’ ( "Spheres" ) ’ )

write ( 3 , ’ ( "C-alpha-circle project. File created using Calculator." ) ’ )  1465

write ( 3 , ’ ( a25 ) ’ , advance=’no ’ ) r e l e a s e i n f o
write ( 4 , ’ ( a25 ) ’ ) r e l e a s e i n f o

write ( 3 , ’ ( "\begin{matrix}c \\ ◯ \end{matrix}.  2002 EPFL Lausanne" ) ’ )

write ( 3 , ’ ( "Atom type : " , a4 ) ’ ) atomtype
write ( 4 , ’ ( a4 ) ’ ) atomtype  1470

write ( 3 , ’ ( "ResolutionÅ [] : " , f 4 . 2 ) ’ ) r e s o l u t i o n

write ( 4 , ’ ( f 4 . 2 ) ’ ) r e s o l u t i o n

1475 write ( 3 , ’ ( "Spherical thicknessÅ [] : " , f 4 . 2 ) ’ ) t h i c k n e s s

write ( 4 , ’ ( f 4 . 2 ) ’ ) t h i c k n e s s

write ( 3 , ’ ( "Thickness combination : " , i3 , ", " , i3 , ", " , i3 , " and " , &

i 3 ) ’ ) g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,1) , &
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,2) , &  1480
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,3) , &
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,4)
write ( 4 , ’ ( 4 i 4 ) ’ ) g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,1) , &
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,2) , &  1485
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,3) , &
g l o b a l l o c a t i o n ( minloc ( g l o b a l r a d i u s ) ,4)
write ( 3 , ’ ( "Memory : " , i 3 ) ’ ) mem
write ( 4 , ’ ( i 3 ) ’ ) mem  1490
i f ( secondrun ) then
! Saving s t a t i s t i c s r e s u l t s to disk

1495

filename ( i ) ( position +1: position+1) = ’ s ’
filename ( i ) ( position +2: position+2) = ’d ’
filename ( i ) ( position +3: position+3) = ’ t ’
1500 open (UNIT=5 , FILE=e x p o r t f o l d e r // ’/ ’// s u b f o l d e r // ’/ ’// machinefolder
// ’/ ’// filename ( i ) , FORM=’FORMATTED’ , STATUS=’REPLACE’ )
do j = 1 , l i m i t
! D i s t r i b u t i o n o f the s m a l l e s t r a d i i

1505

write ( 5 , ’ ( f10 . 2 , i 1 0 ) ’ ) ( ( j / 1 0 0 . 0 ) + 1 . 0 ) , d i s t r i b u t i o n ( j )
end do
1510 close (UNIT=5 , STATUS=’KEEP’ )
end i f
write ( 3 , ∗)
write ( 4 , ∗)  1515
write ( 3 , ’ ( "Entry      Local    f-ratio   Dis" ) ’ , advance = ’ no ’ )
do j = 1 , mem
write ( 3 , ’ ( "         " , i 2 ) ’ , advance = ’ no ’ ) j  1520
write ( 3 , ’ ( ". global radius " ) ’ , advance = ’ no ’ )
end do
write ( 3 , ’ ( " Residue" ) ’ )  1525
write ( 3 , ∗)
! Writing data to f i l e
do j = 1 , ( size ( alphacarbons ) )  1530
a = minval ( g l o b a l l o c a t i o n ( j , : ) )
d = maxval ( g l o b a l l o c a t i o n ( j , : ) )
b = minval ( g l o b a l l o c a t i o n ( j , : ) , MASK = g l o b a l l o c a t i o n ( j , : ) .NE. a )
c = maxval ( g l o b a l l o c a t i o n ( j , : ) , MASK = g l o b a l l o c a t i o n ( j , : ) .NE. d)  1535
! Spacing o f the indexes
d i s l o c a t i o n = ( d − c ) + ( c − b) + ( b − a )

1540

write ( 3 , ’ ( i5 , 2 f10 . 3 , i 5 ) ’ , advance = ’ no ’ ) j , &
l o c a l r a d i u s ( j ) / thickness , f ( j ) , d i s l o c a t i o n
write ( 4 , ’ ( i5 , 2 f10 . 3 , i 5 ) ’ , advance = ’ no ’ ) j , &
l o c a l r a d i u s ( j ) / thickness , f ( j ) , d i s l o c a t i o n  1545
do k = 1 , mem
write ( 3 , ’ ( " (" , f 5 . 3 , ";" , 4 i4 , ") " ) ’ , advance = ’ no ’ ) &
s m a l l e s t r a d i u s ( j , k ) / thickness , &  1550
s m a l l e s t l o c a t i o n ( j , k , 1 ) , &
s m a l l e s t l o c a t i o n ( j , k , 2 ) , &
s m a l l e s t l o c a t i o n ( j , k , 3 ) , &
s m a l l e s t l o c a t i o n ( j , k , 4 )

1555

write ( 4 , ’ ( f 7 . 3 ) ’ , advance = ’ no ’ ) s m a l l e s t r a d i u s ( j , k ) / t h i c k n e s s
end do
1560 write ( 3 , ’ ( i 8 ) ’ ) alphacarbons ( j )%r e s i d u e
write ( 4 , ’ ( i 8 ) ’ ) alphacarbons ( j )%r e s i d u e
end do

write ( ∗ , ’ ( " -> This protein has a thickness of spherically " , f 5 . 3 , "  1565

Angstroem." ) ’ ) t h i c k n e s s
write ( ∗ , ∗)
close (UNIT=3 , STATUS=’KEEP’ )
close (UNIT=4 , STATUS=’KEEP’ )

1570

deallocate ( l o c a l r a d i u s , g l o b a l r a d i u s , f , g , g l o b a l l o c a t i o n , s m a l l e s t r a d i u s

)

deallocate ( s m a l l e s t l o c a t i o n )
end i f

1575

end do
deallocate ( filename )
deallocate ( g o o d f i l e s )

1580

write ( ∗ , ’ ( "Program execution finished" ) ’ )

write ( ∗ , ∗)
end program c a l c u l a t o r