| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // This example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708
30 // Delage et al. PDB4DNA: implementation of DN 30 // Delage et al. PDB4DNA: implementation of DNA geometry from the Protein Data
31 // Bank (PDB) description for 31 // Bank (PDB) description for Geant4-DNA Monte-Carlo
32 // simulations (submitted to 32 // simulations (submitted to Comput. Phys. Commun.)
33 // The Geant4-DNA web site is available at htt 33 // The Geant4-DNA web site is available at http://geant4-dna.org
34 // 34 //
35 // << 35 //
36 /// \file PDBlib.cc 36 /// \file PDBlib.cc
37 /// \brief Implementation file for PDBlib clas 37 /// \brief Implementation file for PDBlib class
38 38
39 #include "PDBlib.hh" 39 #include "PDBlib.hh"
40 40
41 // define if the program is running with Geant << 41 //define if the program is running with Geant4
42 #define GEANT4 42 #define GEANT4
43 43
44 #ifdef GEANT4 44 #ifdef GEANT4
45 // Specific to Geant4, globals.hh is used for << 45 //Specific to Geant4, globals.hh is used for G4cout
46 # include "globals.hh" << 46 #include "globals.hh"
47 #else 47 #else
48 # define G4cout std::cout << 48 #define G4cout std::cout
49 # define G4cerr std::cerr << 49 #define G4cerr std::cerr
50 # define G4endl std::endl << 50 #define G4endl std::endl
51 # define G4String std::string // string incl << 51 #define G4String std::string // string included in PDBatom.hh
52 # include <cfloat> << 52 #include <cfloat>
53 #endif 53 #endif
54 #include <cmath> <<
55 #include <fstream> 54 #include <fstream>
56 #include <iostream> 55 #include <iostream>
57 #include <limits> 56 #include <limits>
>> 57 #include <cmath>
58 #include <sstream> 58 #include <sstream>
59 #include <stdlib.h> 59 #include <stdlib.h>
60 60
61 using namespace std; 61 using namespace std;
62 62
63 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 64
65 PDBlib::PDBlib() : fNbNucleotidsPerStrand(0) { << 65 PDBlib::PDBlib() :
>> 66 fNbNucleotidsPerStrand(0)
>> 67 {
>> 68 }
66 69
67 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68 71
69 /** 72 /**
70 * \brief Load a PDB file into memory 73 * \brief Load a PDB file into memory
71 * \details molecule (polymer,?), << 74 * \details molecule (polymer,?),
72 * read line, key words 75 * read line, key words
73 * \param filename G4String for filename 76 * \param filename G4String for filename
74 * \param isDNA 77 * \param isDNA
75 * \param verbose 78 * \param verbose
76 * \return List of molecules 79 * \return List of molecules
77 */ 80 */
78 81
79 Molecule* PDBlib::Load(const std::string& file << 82 Molecule * PDBlib::Load(const std::string& filename,
80 unsigned short int verb << 83 unsigned short int& isDNA,
>> 84 unsigned short int verbose = 0)
81 { 85 {
82 G4String sLine = ""; 86 G4String sLine = "";
83 ifstream infile; 87 ifstream infile;
84 infile.open(filename.c_str()); 88 infile.open(filename.c_str());
85 if (!infile) { << 89 if(!infile)
86 G4cout << "PDBlib::load >> file " << filen << 90 {
>> 91 G4cout << "PDBlib::load >> file " << filename << " not found !!!!"
>> 92 << G4endl;
87 } 93 }
88 else { << 94 else
89 int nbAtom = 0; // total number of atoms << 95 {
90 int numAtomInRes = 0; // number of an ato << 96 int nbAtom = 0; // total number of atoms in a residue
91 int nbResidue = 0; // total number of res << 97 int numAtomInRes = 0; // number of an atom in a residue
92 int nbMolecule = 0; // total number of mo << 98 int nbResidue = 0; // total number of residues
93 << 99 int nbMolecule = 0; // total number of molecule
94 Atom* AtomicOld = nullptr; << 100
95 Atom* AtomicNext = nullptr; << 101 Atom * AtomicOld = nullptr;
96 << 102 Atom * AtomicNext = nullptr;
97 int serial; // Atom serial number << 103
98 G4String atomName; // Atom name << 104 int serial; //Atom serial number
99 G4String element; // Element Symbol << 105 G4String atomName; //Atom name
100 G4String resName; // Residue name for thi << 106 G4String element; //Element Symbol
101 double x, y, z; // Orthogonal coordinates << 107 G4String resName; //Residue name for this atom
102 double occupancy; // Occupancy << 108 double x, y, z; //Orthogonal coordinates in Angstroms
103 << 109 double occupancy; //Occupancy
104 Residue* residueOld = nullptr; << 110
105 Residue* residueFirst = nullptr; << 111 Residue * residueOld = nullptr;
106 Residue* residueNext = nullptr; << 112 Residue * residueFirst = nullptr;
107 << 113 Residue * residueNext = nullptr;
108 Molecule* moleculeFirst = nullptr; << 114
109 Molecule* moleculeOld = nullptr; << 115 Molecule * moleculeFirst = nullptr;
110 Molecule* moleculeNext = nullptr; << 116 Molecule * moleculeOld = nullptr;
>> 117 Molecule * moleculeNext = nullptr;
111 118
112 ///////////////////////////////// 119 /////////////////////////////////
113 // NEW variable to draw a fitting cylinder << 120 //NEW variable to draw a fitting cylinder if z oriented
114 //=> fitting box 121 //=> fitting box
115 double minGlobZ, maxGlobZ; 122 double minGlobZ, maxGlobZ;
116 double minGlobX, maxGlobX; 123 double minGlobX, maxGlobX;
117 double minGlobY, maxGlobY; 124 double minGlobY, maxGlobY;
118 double minX, maxX, minY, maxY, minZ, maxZ; << 125 double minX, maxX, minY, maxY, minZ, maxZ; //Sort of 'mother volume' box
119 126
120 minGlobZ = -DBL_MAX; 127 minGlobZ = -DBL_MAX;
121 minGlobX = -DBL_MAX; 128 minGlobX = -DBL_MAX;
122 minGlobY = -DBL_MAX; 129 minGlobY = -DBL_MAX;
123 maxGlobZ = DBL_MAX; 130 maxGlobZ = DBL_MAX;
124 maxGlobX = DBL_MAX; 131 maxGlobX = DBL_MAX;
125 maxGlobY = DBL_MAX; 132 maxGlobY = DBL_MAX;
126 minX = -DBL_MAX; 133 minX = -DBL_MAX;
127 minY = -DBL_MAX; 134 minY = -DBL_MAX;
128 minZ = -DBL_MAX; 135 minZ = -DBL_MAX;
129 maxX = DBL_MAX; 136 maxX = DBL_MAX;
130 maxY = DBL_MAX; 137 maxY = DBL_MAX;
131 maxZ = DBL_MAX; 138 maxZ = DBL_MAX;
132 139
133 int lastResSeq = -1; 140 int lastResSeq = -1;
134 int resSeq = 0; 141 int resSeq = 0;
135 142
136 G4String nameMol; 143 G4String nameMol;
137 144
138 unsigned short int numModel = 0; 145 unsigned short int numModel = 0;
139 unsigned short int model = 0; 146 unsigned short int model = 0;
140 147
141 // Number of TER << 148 //Number of TER
142 int ter = 0; 149 int ter = 0;
143 // Number of TER (chain) max for this file << 150 //Number of TER (chain) max for this file
144 151
145 int terMax = INT_MAX; 152 int terMax = INT_MAX;
146 153
147 if (!infile.eof()) { << 154 if(!infile.eof())
>> 155 {
148 getline(infile, sLine); 156 getline(infile, sLine);
149 std::size_t found = sLine.find("DNA"); 157 std::size_t found = sLine.find("DNA");
150 if (found != G4String::npos) { << 158 if(found != G4String::npos)
>> 159 {
151 terMax = 2; 160 terMax = 2;
152 isDNA = 1; 161 isDNA = 1;
153 } 162 }
154 else << 163 else isDNA = 0;
155 isDNA = 0; << 164 //If PDB file have not a header line
156 // If PDB file have not a header line <<
157 found = sLine.find("HEADER"); 165 found = sLine.find("HEADER");
158 if (found == G4String::npos) { << 166 if(found == G4String::npos)
>> 167 {
159 infile.close(); 168 infile.close();
160 infile.open(filename.c_str()); 169 infile.open(filename.c_str());
161 170
162 G4cout << "PDBlib::load >> No header f 171 G4cout << "PDBlib::load >> No header found !!!!" << G4endl;
163 } 172 }
164 } 173 }
165 174
166 while (!infile.eof()) { << 175 while(!infile.eof())
>> 176 {
167 getline(infile, sLine); 177 getline(infile, sLine);
168 // In the case of several models << 178 //In the case of several models
169 if ((sLine.substr(0, 6)).compare("NUMMDL << 179 if((sLine.substr(0, 6)).compare("NUMMDL") == 0)
>> 180 {
170 istringstream((sLine.substr(10, 4))) > 181 istringstream((sLine.substr(10, 4))) >> numModel;
171 } 182 }
172 if ((numModel > 0) && ((sLine.substr(0, << 183 if((numModel > 0) && ((sLine.substr(0, 6)).compare("MODEL ") == 0))
>> 184 {
173 istringstream((sLine.substr(10, 4))) > 185 istringstream((sLine.substr(10, 4))) >> model;
174 ////////////////////////////////////// 186 //////////////////////////////////////////
175 if (model > 1) break; << 187 if(model > 1) break;
176 ////////////////////////////////////// 188 //////////////////////////////////////////
177 } 189 }
178 // For verification of residue sequence << 190 //For verification of residue sequence
179 if ((sLine.substr(0, 6)).compare("SEQRES << 191 if((sLine.substr(0, 6)).compare("SEQRES") == 0)
180 // Create list of molecule here << 192 {
181 if (verbose > 1) G4cout << sLine << G4 << 193 //Create list of molecule here
>> 194 if(verbose > 1) G4cout << sLine << G4endl;
182 } 195 }
183 196
184 // Coordinate section << 197 //Coordinate section
185 if ((numModel > 0) && ((sLine.substr(0, << 198 if((numModel > 0) && ((sLine.substr(0, 6)).compare("ENDMDL") == 0))
>> 199 {
186 ; 200 ;
187 } 201 }
188 else if ((sLine.substr(0, 6)).compare("T << 202 else if((sLine.substr(0, 6)).compare("TER ") == 0) //3 spaces
189 { 203 {
190 ////////////////////////////////////// 204 //////////////////////////////////////////
191 // Currently retrieve only the two fir << 205 //Currently retrieve only the two first chains(TER) => two DNA strands
192 ter++; 206 ter++;
193 if (ter > terMax) break; << 207 if(ter > terMax) break;
194 ////////////////////////////////////// 208 //////////////////////////////////////////
195 209
196 // if (verbose>1) G4cout << sLine << G << 210 //if (verbose>1) G4cout << sLine << G4endl;
197 /************ Begin TER ************** 211 /************ Begin TER ******************/
198 lastResSeq = -1; 212 lastResSeq = -1;
199 resSeq = 0; 213 resSeq = 0;
200 214
201 AtomicOld->SetNext(NULL); 215 AtomicOld->SetNext(NULL);
202 residueOld->SetNext(NULL); 216 residueOld->SetNext(NULL);
203 residueOld->fNbAtom = nbAtom; 217 residueOld->fNbAtom = nbAtom;
204 218
205 // Molecule creation: << 219 //Molecule creation:
206 if (moleculeOld == NULL) { << 220 if(moleculeOld == NULL)
207 nameMol = filename; //+numModel << 221 {
>> 222 nameMol = filename; //+numModel
208 moleculeOld = new Molecule(nameMol, 223 moleculeOld = new Molecule(nameMol, nbMolecule);
209 moleculeOld->SetFirst(residueFirst); 224 moleculeOld->SetFirst(residueFirst);
210 moleculeOld->fNbResidue = nbResidue; 225 moleculeOld->fNbResidue = nbResidue;
211 moleculeFirst = moleculeOld; 226 moleculeFirst = moleculeOld;
212 } 227 }
213 else { << 228 else
>> 229 {
214 moleculeNext = new Molecule(nameMol, 230 moleculeNext = new Molecule(nameMol, nbMolecule);
215 moleculeOld->SetNext(moleculeNext); 231 moleculeOld->SetNext(moleculeNext);
216 moleculeNext->SetFirst(residueFirst) 232 moleculeNext->SetFirst(residueFirst);
217 moleculeNext->fNbResidue = nbResidue 233 moleculeNext->fNbResidue = nbResidue;
218 moleculeOld = moleculeNext; 234 moleculeOld = moleculeNext;
219 } 235 }
220 236
221 nbMolecule++; 237 nbMolecule++;
222 moleculeOld->SetNext(NULL); 238 moleculeOld->SetNext(NULL);
223 moleculeOld->fCenterX = (int)((minX + << 239 moleculeOld->fCenterX = (int) ((minX + maxX) / 2.);
224 moleculeOld->fCenterY = (int)((minY + << 240 moleculeOld->fCenterY = (int) ((minY + maxY) / 2.);
225 moleculeOld->fCenterZ = (int)((minZ + << 241 moleculeOld->fCenterZ = (int) ((minZ + maxZ) / 2.);
226 moleculeOld->fMaxGlobZ = maxGlobZ; 242 moleculeOld->fMaxGlobZ = maxGlobZ;
227 moleculeOld->fMinGlobZ = minGlobZ; 243 moleculeOld->fMinGlobZ = minGlobZ;
228 moleculeOld->fMaxGlobX = maxGlobX; 244 moleculeOld->fMaxGlobX = maxGlobX;
229 moleculeOld->fMinGlobX = minGlobX; 245 moleculeOld->fMinGlobX = minGlobX;
230 moleculeOld->fMaxGlobY = maxGlobY; 246 moleculeOld->fMaxGlobY = maxGlobY;
231 moleculeOld->fMinGlobY = minGlobY; 247 moleculeOld->fMinGlobY = minGlobY;
232 248
233 minGlobZ = -DBL_MAX; 249 minGlobZ = -DBL_MAX;
234 minGlobX = -DBL_MAX; 250 minGlobX = -DBL_MAX;
235 minGlobY = -DBL_MAX; 251 minGlobY = -DBL_MAX;
236 maxGlobZ = DBL_MAX; 252 maxGlobZ = DBL_MAX;
237 maxGlobX = DBL_MAX; 253 maxGlobX = DBL_MAX;
238 maxGlobY = DBL_MAX; 254 maxGlobY = DBL_MAX;
239 minX = -DBL_MAX; 255 minX = -DBL_MAX;
240 minY = -DBL_MAX; 256 minY = -DBL_MAX;
241 minZ = -DBL_MAX; 257 minZ = -DBL_MAX;
242 maxX = DBL_MAX; 258 maxX = DBL_MAX;
243 maxY = DBL_MAX; 259 maxY = DBL_MAX;
244 maxZ = DBL_MAX; 260 maxZ = DBL_MAX;
245 261
246 nbAtom = 0; 262 nbAtom = 0;
247 numAtomInRes = 0; 263 numAtomInRes = 0;
248 nbResidue = 0; 264 nbResidue = 0;
249 AtomicOld = NULL; 265 AtomicOld = NULL;
250 AtomicNext = NULL; 266 AtomicNext = NULL;
251 residueOld = NULL; 267 residueOld = NULL;
252 residueFirst = NULL; 268 residueFirst = NULL;
253 residueNext = NULL; 269 residueNext = NULL;
254 270
255 ///////////// End TER //////////////// 271 ///////////// End TER ///////////////////
256 } 272 }
257 else if ((sLine.substr(0, 6)).compare("A << 273 else if((sLine.substr(0, 6)).compare("ATOM ") == 0)
>> 274 {
>> 275
258 /************ Begin ATOM ************* 276 /************ Begin ATOM ******************/
259 // serial << 277 //serial
260 istringstream((sLine.substr(6, 5))) >> 278 istringstream((sLine.substr(6, 5))) >> serial;
261 279
262 // To be improved << 280 //To be improved
263 // atomName : << 281 //atomName :
264 atomName = sLine.substr(12, 4); 282 atomName = sLine.substr(12, 4);
265 if (atomName.substr(0, 1).compare(" ") << 283 if(atomName.substr(0, 1).compare(" ") == 0) element = sLine.substr(13,
266 element = sLine.substr(13, 1); << 284 1);
267 else << 285 else element = sLine.substr(12, 1);
268 element = sLine.substr(12, 1); <<
269 286
270 // set Van der Waals radius expressed 287 // set Van der Waals radius expressed in Angstrom
271 double vdwRadius = -1.; 288 double vdwRadius = -1.;
272 if (element == "H") { << 289 if(element == "H")
>> 290 {
273 vdwRadius = 1.2; 291 vdwRadius = 1.2;
274 } 292 }
275 else if (element == "C") { << 293 else if(element == "C")
>> 294 {
276 vdwRadius = 1.7; 295 vdwRadius = 1.7;
277 } 296 }
278 else if (element == "N") { << 297 else if(element == "N")
>> 298 {
279 vdwRadius = 1.55; 299 vdwRadius = 1.55;
280 } 300 }
281 else if (element == "O") { << 301 else if(element == "O")
>> 302 {
282 vdwRadius = 1.52; 303 vdwRadius = 1.52;
283 } 304 }
284 else if (element == "P") { << 305 else if(element == "P")
>> 306 {
285 vdwRadius = 1.8; 307 vdwRadius = 1.8;
286 } 308 }
287 else if (element == "S") { << 309 else if(element == "S")
>> 310 {
288 vdwRadius = 1.8; 311 vdwRadius = 1.8;
289 } 312 }
290 else { << 313 else
291 #ifndef GEANT4 << 314 {
>> 315 #ifndef GEANT4
292 G4cerr << "Element not recognized : 316 G4cerr << "Element not recognized : " << element << G4endl;
293 G4cerr << "Stop now" << G4endl; 317 G4cerr << "Stop now" << G4endl;
294 exit(1); 318 exit(1);
295 #else 319 #else
296 G4ExceptionDescription errMsg; 320 G4ExceptionDescription errMsg;
297 errMsg << "Element not recognized : 321 errMsg << "Element not recognized : " << element << G4endl;
298 322
299 G4Exception("PDBlib::Load", "ELEM_NO << 323 G4Exception("PDBlib::Load",
300 #endif << 324 "ELEM_NOT_RECOGNIZED",
>> 325 FatalException,
>> 326 errMsg);
>> 327 #endif
301 } 328 }
302 329
303 { 330 {
304 // resName : << 331 //resName :
305 resName = sLine.substr(17, 3); 332 resName = sLine.substr(17, 3);
306 // resSeq : << 333 //resSeq :
307 istringstream((sLine.substr(22, 4))) 334 istringstream((sLine.substr(22, 4))) >> resSeq;
308 // x,y,z : << 335 //x,y,z :
309 stringstream((sLine.substr(30, 8))) 336 stringstream((sLine.substr(30, 8))) >> x;
310 stringstream((sLine.substr(38, 8))) 337 stringstream((sLine.substr(38, 8))) >> y;
311 stringstream((sLine.substr(46, 8))) 338 stringstream((sLine.substr(46, 8))) >> z;
312 // occupancy : << 339 //occupancy :
313 occupancy = 1.; 340 occupancy = 1.;
314 341
315 if (minGlobZ < z) minGlobZ = z; << 342 if(minGlobZ < z) minGlobZ = z;
316 if (maxGlobZ > z) maxGlobZ = z; << 343 if(maxGlobZ > z) maxGlobZ = z;
317 if (minGlobX < x) minGlobX = x; << 344 if(minGlobX < x) minGlobX = x;
318 if (maxGlobX > x) maxGlobX = x; << 345 if(maxGlobX > x) maxGlobX = x;
319 if (minGlobY < y) minGlobY = y; << 346 if(minGlobY < y) minGlobY = y;
320 if (maxGlobY > y) maxGlobY = y; << 347 if(maxGlobY > y) maxGlobY = y;
321 if (minX > x) minX = x; << 348 if(minX > x) minX = x;
322 if (maxX < x) maxX = x; << 349 if(maxX < x) maxX = x;
323 if (minY > y) minY = y; << 350 if(minY > y) minY = y;
324 if (maxY < y) maxY = y; << 351 if(maxY < y) maxY = y;
325 if (minZ > z) minZ = z; << 352 if(minZ > z) minZ = z;
326 if (maxZ < z) maxZ = z; << 353 if(maxZ < z) maxZ = z;
327 << 354
328 // treatment for Atoms: << 355 //treatment for Atoms:
329 if (AtomicOld == NULL) { << 356 if(AtomicOld == NULL)
330 AtomicOld = << 357 {
331 new Atom(serial, atomName, "", 0 << 358 AtomicOld = new Atom(serial,
332 AtomicOld->SetNext(NULL); // If o << 359 atomName,
333 } << 360 "",
334 else { << 361 0,
335 AtomicNext = << 362 0,
336 new Atom(serial, atomName, "", 0 << 363 x,
>> 364 y,
>> 365 z,
>> 366 vdwRadius,
>> 367 occupancy,
>> 368 0,
>> 369 element);
>> 370 AtomicOld->SetNext(NULL); //If only one Atom inside residue
>> 371 }
>> 372 else
>> 373 {
>> 374 AtomicNext = new Atom(serial,
>> 375 atomName,
>> 376 "",
>> 377 0,
>> 378 0,
>> 379 x,
>> 380 y,
>> 381 z,
>> 382 vdwRadius,
>> 383 occupancy,
>> 384 0,
>> 385 element);
337 AtomicOld->SetNext(AtomicNext); 386 AtomicOld->SetNext(AtomicNext);
338 AtomicOld = AtomicNext; 387 AtomicOld = AtomicNext;
339 } 388 }
340 nbAtom++; 389 nbAtom++;
341 } // END if (element!="H") << 390 } //END if (element!="H")
342 /****************************Begin Res 391 /****************************Begin Residue************************/
343 // treatment for residues: << 392 //treatment for residues:
344 if (residueOld == NULL) { << 393 if(residueOld == NULL)
345 if (verbose > 2) G4cout << "residueO << 394 {
>> 395 if(verbose > 2) G4cout << "residueOld == NULL" << G4endl;
346 396
347 AtomicOld->fNumInRes = 0; 397 AtomicOld->fNumInRes = 0;
348 residueOld = new Residue(resName, re 398 residueOld = new Residue(resName, resSeq);
349 residueOld->SetFirst(AtomicOld); 399 residueOld->SetFirst(AtomicOld);
350 residueOld->SetNext(NULL); 400 residueOld->SetNext(NULL);
351 residueFirst = residueOld; 401 residueFirst = residueOld;
352 lastResSeq = resSeq; 402 lastResSeq = resSeq;
353 nbResidue++; 403 nbResidue++;
354 } 404 }
355 else { << 405 else
356 if (lastResSeq == resSeq) { << 406 {
>> 407 if(lastResSeq == resSeq)
>> 408 {
357 numAtomInRes++; 409 numAtomInRes++;
358 AtomicOld->fNumInRes = numAtomInRe 410 AtomicOld->fNumInRes = numAtomInRes;
359 } 411 }
360 else { << 412 else
>> 413 {
361 numAtomInRes = 0; 414 numAtomInRes = 0;
362 AtomicOld->fNumInRes = numAtomInRe 415 AtomicOld->fNumInRes = numAtomInRes;
363 416
364 residueNext = new Residue(resName, 417 residueNext = new Residue(resName, resSeq);
365 residueNext->SetFirst(AtomicOld); 418 residueNext->SetFirst(AtomicOld);
366 residueOld->SetNext(residueNext); 419 residueOld->SetNext(residueNext);
367 residueOld->fNbAtom = nbAtom - 1; 420 residueOld->fNbAtom = nbAtom - 1;
368 421
369 nbAtom = 1; 422 nbAtom = 1;
370 residueOld = residueNext; 423 residueOld = residueNext;
371 lastResSeq = resSeq; 424 lastResSeq = resSeq;
372 nbResidue++; 425 nbResidue++;
373 } 426 }
374 } 427 }
375 /////////////////////////End Residue// 428 /////////////////////////End Residue////////////
376 ///////////// End ATOM /////////////// 429 ///////////// End ATOM ///////////////////
377 } // END if Atom << 430 } //END if Atom
378 } // END while (!infile.eof()) << 431 } //END while (!infile.eof())
379 432
380 infile.close(); 433 infile.close();
381 return moleculeFirst; 434 return moleculeFirst;
382 } // END else if (!infile) << 435 } //END else if (!infile)
383 return NULL; 436 return NULL;
384 } 437 }
385 438
386 //....oooOO0OOooo........oooOO0OOooo........oo 439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
387 440
388 /** 441 /**
389 * \brief Compute barycenters 442 * \brief Compute barycenters
390 * \details Compute barycenters and its coord << 443 * \details Compute barycenters and its coordinate
391 * for nucleotides 444 * for nucleotides
392 * \param moleculeListTemp * moleculeLis 445 * \param moleculeListTemp * moleculeList
393 * \return Barycenter * 446 * \return Barycenter *
394 */ 447 */
395 Barycenter* PDBlib::ComputeNucleotideBarycente << 448 Barycenter * PDBlib::ComputeNucleotideBarycenters(Molecule * moleculeListTemp)
396 { 449 {
397 //////////////////////////////////////////// 450 ///////////////////////////////////////////////////////
398 // Placement and physical volume constructio << 451 //Placement and physical volume construction from memory
399 Barycenter* BarycenterFirst = NULL; << 452 Barycenter * BarycenterFirst = NULL;
400 Barycenter* BarycenterOld = NULL; << 453 Barycenter * BarycenterOld = NULL;
401 Barycenter* BarycenterNext = NULL; << 454 Barycenter * BarycenterNext = NULL;
402 << 455
403 // Residue (Base, Phosphate,sugar) list << 456 //Residue (Base, Phosphate,sugar) list
404 Residue* residueListTemp; << 457 Residue *residueListTemp;
405 // Atom list inside a residu << 458 //Atom list inside a residu
406 Atom* AtomTemp; << 459 Atom *AtomTemp;
407 460
408 int k = 0; 461 int k = 0;
409 int j_old = 0; 462 int j_old = 0;
410 463
411 while (moleculeListTemp) { << 464 while(moleculeListTemp)
>> 465 {
412 residueListTemp = moleculeListTemp->GetFir 466 residueListTemp = moleculeListTemp->GetFirst();
413 467
414 k++; 468 k++;
415 int j = 0; 469 int j = 0;
416 470
417 // Check numerotation style (1->n per stra << 471 //Check numerotation style (1->n per strand or 1->2n for two strand)
418 int correctNumerotationNumber = 0; 472 int correctNumerotationNumber = 0;
419 if (k == 2 && residueListTemp->fResSeq > 1 << 473 if(k == 2 && residueListTemp->fResSeq > 1)
>> 474 {
420 correctNumerotationNumber = residueListT 475 correctNumerotationNumber = residueListTemp->fResSeq;
421 } 476 }
422 477
423 while (residueListTemp) { << 478 while(residueListTemp)
>> 479 {
424 AtomTemp = residueListTemp->GetFirst(); 480 AtomTemp = residueListTemp->GetFirst();
425 j++; 481 j++;
426 482
427 // Correction consequently to numerotati << 483 //Correction consequently to numerotation check
428 if (correctNumerotationNumber != 0) { << 484 if(correctNumerotationNumber != 0)
429 residueListTemp->fResSeq = residueList << 485 {
>> 486 residueListTemp->fResSeq = residueListTemp->fResSeq
>> 487 - correctNumerotationNumber
>> 488 + 1;
430 } 489 }
431 490
432 // Barycenter computation << 491 //Barycenter computation
433 double baryX = 0., baryY = 0., baryZ = 0 492 double baryX = 0., baryY = 0., baryZ = 0.;
434 double baryBaseX = 0., baryBaseY = 0., b 493 double baryBaseX = 0., baryBaseY = 0., baryBaseZ = 0.;
435 double barySugX = 0., barySugY = 0., bar 494 double barySugX = 0., barySugY = 0., barySugZ = 0.;
436 double baryPhosX = 0., baryPhosY = 0., b 495 double baryPhosX = 0., baryPhosY = 0., baryPhosZ = 0.;
437 unsigned short int nbAtomInBase = 0; 496 unsigned short int nbAtomInBase = 0;
438 497
439 for (int i = 0; i < residueListTemp->fNb << 498 for(int i = 0; i < residueListTemp->fNbAtom; i++)
440 // Compute barycenter of the nucletoti << 499 {
>> 500 //Compute barycenter of the nucletotide
441 baryX += AtomTemp->fX; 501 baryX += AtomTemp->fX;
442 baryY += AtomTemp->fY; 502 baryY += AtomTemp->fY;
443 baryZ += AtomTemp->fZ; 503 baryZ += AtomTemp->fZ;
444 // Compute barycenters for Base Sugar << 504 //Compute barycenters for Base Sugar Phosphat
445 if (residueListTemp->fResSeq == 1) { << 505 if(residueListTemp->fResSeq == 1)
446 if (i == 0) { << 506 {
>> 507 if(i == 0)
>> 508 {
447 baryPhosX += AtomTemp->fX; 509 baryPhosX += AtomTemp->fX;
448 baryPhosY += AtomTemp->fY; 510 baryPhosY += AtomTemp->fY;
449 baryPhosZ += AtomTemp->fZ; 511 baryPhosZ += AtomTemp->fZ;
450 } 512 }
451 else if (i < 8) { << 513 else if(i < 8)
>> 514 {
452 barySugX += AtomTemp->fX; 515 barySugX += AtomTemp->fX;
453 barySugY += AtomTemp->fY; 516 barySugY += AtomTemp->fY;
454 barySugZ += AtomTemp->fZ; 517 barySugZ += AtomTemp->fZ;
455 } 518 }
456 else { << 519 else
457 // hydrogen are placed at the end << 520 {
458 // We don't want them for this cal << 521 //hydrogen are placed at the end of the residue in a PDB file
459 if (AtomTemp->fElement != "H") { << 522 //We don't want them for this calculation
>> 523 if(AtomTemp->fElement != "H")
>> 524 {
460 baryBaseX += AtomTemp->fX; 525 baryBaseX += AtomTemp->fX;
461 baryBaseY += AtomTemp->fY; 526 baryBaseY += AtomTemp->fY;
462 baryBaseZ += AtomTemp->fZ; 527 baryBaseZ += AtomTemp->fZ;
463 nbAtomInBase++; 528 nbAtomInBase++;
464 } 529 }
465 } 530 }
466 } 531 }
467 else { << 532 else
468 if (i < 4) { << 533 {
>> 534 if(i < 4)
>> 535 {
469 baryPhosX += AtomTemp->fX; 536 baryPhosX += AtomTemp->fX;
470 baryPhosY += AtomTemp->fY; 537 baryPhosY += AtomTemp->fY;
471 baryPhosZ += AtomTemp->fZ; 538 baryPhosZ += AtomTemp->fZ;
472 } 539 }
473 else if (i < 11) { << 540 else if(i < 11)
>> 541 {
474 barySugX += AtomTemp->fX; 542 barySugX += AtomTemp->fX;
475 barySugY += AtomTemp->fY; 543 barySugY += AtomTemp->fY;
476 barySugZ += AtomTemp->fZ; 544 barySugZ += AtomTemp->fZ;
477 } 545 }
478 else { << 546 else
479 // hydrogen are placed at the end << 547 {
480 // We don't want them for this cal << 548 //hydrogen are placed at the end of the residue in a PDB file
481 if (AtomTemp->fElement != "H") { << 549 //We don't want them for this calculation
>> 550 if(AtomTemp->fElement != "H")
>> 551 { // break;
482 baryBaseX += AtomTemp->fX; 552 baryBaseX += AtomTemp->fX;
483 baryBaseY += AtomTemp->fY; 553 baryBaseY += AtomTemp->fY;
484 baryBaseZ += AtomTemp->fZ; 554 baryBaseZ += AtomTemp->fZ;
485 nbAtomInBase++; 555 nbAtomInBase++;
486 } 556 }
487 } 557 }
488 } 558 }
489 AtomTemp = AtomTemp->GetNext(); 559 AtomTemp = AtomTemp->GetNext();
490 } // end of for ( i=0 ; i < residueLis << 560 } //end of for ( i=0 ; i < residueListTemp->nbAtom ; i++)
491 561
492 baryX = baryX / (double)residueListTemp- << 562 baryX = baryX / (double) residueListTemp->fNbAtom;
493 baryY = baryY / (double)residueListTemp- << 563 baryY = baryY / (double) residueListTemp->fNbAtom;
494 baryZ = baryZ / (double)residueListTemp- << 564 baryZ = baryZ / (double) residueListTemp->fNbAtom;
495 565
496 if (residueListTemp->fResSeq != 1) // S << 566 if(residueListTemp->fResSeq != 1) //Special case first Phosphate
497 { 567 {
498 baryPhosX = baryPhosX / 4.; 568 baryPhosX = baryPhosX / 4.;
499 baryPhosY = baryPhosY / 4.; 569 baryPhosY = baryPhosY / 4.;
500 baryPhosZ = baryPhosZ / 4.; 570 baryPhosZ = baryPhosZ / 4.;
501 } 571 }
502 barySugX = barySugX / 7.; 572 barySugX = barySugX / 7.;
503 barySugY = barySugY / 7.; 573 barySugY = barySugY / 7.;
504 barySugZ = barySugZ / 7.; 574 barySugZ = barySugZ / 7.;
505 baryBaseX = baryBaseX / (double)nbAtomIn << 575 baryBaseX = baryBaseX / (double) nbAtomInBase;
506 baryBaseY = baryBaseY / (double)nbAtomIn << 576 baryBaseY = baryBaseY / (double) nbAtomInBase;
507 baryBaseZ = baryBaseZ / (double)nbAtomIn << 577 baryBaseZ = baryBaseZ / (double) nbAtomInBase;
508 << 578
509 // Barycenter creation: << 579 //Barycenter creation:
510 if (BarycenterOld == NULL) { << 580 if(BarycenterOld == NULL)
511 BarycenterOld = new Barycenter(j + j_o << 581 {
512 baryBas << 582 BarycenterOld = new Barycenter(j + j_old, baryX, baryY, baryZ, //j [1..n]
513 barySug << 583 baryBaseX,
>> 584 baryBaseY,
>> 585 baryBaseZ,
>> 586 barySugX,
>> 587 barySugY,
>> 588 barySugZ,
>> 589 baryPhosX,
>> 590 baryPhosY,
>> 591 baryPhosZ);
514 BarycenterFirst = BarycenterOld; 592 BarycenterFirst = BarycenterOld;
515 } 593 }
516 else { << 594 else
517 BarycenterNext = << 595 {
518 new Barycenter(j + j_old, baryX, bar << 596 BarycenterNext = new Barycenter(j + j_old,
519 barySugY, barySugZ, b << 597 baryX,
>> 598 baryY,
>> 599 baryZ,
>> 600 baryBaseX,
>> 601 baryBaseY,
>> 602 baryBaseZ,
>> 603 barySugX,
>> 604 barySugY,
>> 605 barySugZ,
>> 606 baryPhosX,
>> 607 baryPhosY,
>> 608 baryPhosZ);
520 BarycenterOld->SetNext(BarycenterNext) 609 BarycenterOld->SetNext(BarycenterNext);
521 BarycenterOld = BarycenterNext; 610 BarycenterOld = BarycenterNext;
522 } 611 }
523 612
524 //////////////////////////////////////// 613 /////////////////////////////////////////////////
525 // distance computation between all atom << 614 //distance computation between all atoms inside
526 // a residue and the barycenter << 615 //a residue and the barycenter
527 AtomTemp = residueListTemp->GetFirst(); 616 AtomTemp = residueListTemp->GetFirst();
528 double dT3Dp; 617 double dT3Dp;
529 double max = 0.; 618 double max = 0.;
530 for (int ii = 0; ii < residueListTemp->f << 619 for(int ii = 0; ii < residueListTemp->fNbAtom; ii++)
531 dT3Dp = DistanceTwo3Dpoints(AtomTemp-> << 620 {
532 Barycenter << 621 dT3Dp = DistanceTwo3Dpoints(AtomTemp->fX,
>> 622 BarycenterOld->fCenterX,
>> 623 AtomTemp->fY,
>> 624 BarycenterOld->fCenterY,
>> 625 AtomTemp->fZ,
>> 626 BarycenterOld->fCenterZ);
533 BarycenterOld->SetDistance(ii, dT3Dp); 627 BarycenterOld->SetDistance(ii, dT3Dp);
534 if (dT3Dp > max) max = dT3Dp; << 628 if(dT3Dp > max) max = dT3Dp;
535 AtomTemp = AtomTemp->GetNext(); 629 AtomTemp = AtomTemp->GetNext();
536 } // end of for ( i=0 ; i < residueLis << 630 } //end of for ( i=0 ; i < residueListTemp->nbAtom ; i++)
537 631
538 BarycenterOld->SetRadius(max + 1.8); 632 BarycenterOld->SetRadius(max + 1.8);
539 residueListTemp = residueListTemp->GetNe 633 residueListTemp = residueListTemp->GetNext();
540 634
541 } // end of while sur residueListTemp << 635 } //end of while sur residueListTemp
542 636
543 j_old += j; 637 j_old += j;
544 638
545 /// molecs->push_back(*moleculeListTemp); << 639 ///molecs->push_back(*moleculeListTemp);
546 moleculeListTemp = moleculeListTemp->GetNe 640 moleculeListTemp = moleculeListTemp->GetNext();
547 } // end of while sur moleculeListTemp << 641 } //end of while sur moleculeListTemp
548 642
549 if (BarycenterNext != NULL) { << 643 if(BarycenterNext != NULL)
>> 644 {
550 BarycenterNext->SetNext(NULL); 645 BarycenterNext->SetNext(NULL);
551 } 646 }
552 647
553 return BarycenterFirst; 648 return BarycenterFirst;
554 } 649 }
555 650
556 //....oooOO0OOooo........oooOO0OOooo........oo 651 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
557 652
558 /** 653 /**
559 * \brief the corresponding bounding volume 654 * \brief the corresponding bounding volume parameters
560 * \details the corresponding bounding volume 655 * \details the corresponding bounding volume parameters
561 * to build a box from atoms coordi 656 * to build a box from atoms coordinates
562 */ 657 */
563 void PDBlib::ComputeBoundingVolumeParams(Molec << 658 void PDBlib::ComputeBoundingVolumeParams(Molecule *moleculeListTemp,
564 doubl << 659 double &dX,
565 doubl << 660 double &dY,
566 doubl << 661 double &dZ, //Dimensions for bounding volume
>> 662 double &tX,
>> 663 double &tY,
>> 664 double &tZ) //Translation for bounding volume
567 { 665 {
568 double minminX, minminY, minminZ; // minimu << 666 double minminX, minminY, minminZ; //minimum minimorum
569 double maxmaxX, maxmaxY, maxmaxZ; // maximu << 667 double maxmaxX, maxmaxY, maxmaxZ; //maximum maximorum
570 668
571 minminX = DBL_MAX; 669 minminX = DBL_MAX;
572 minminY = DBL_MAX; 670 minminY = DBL_MAX;
573 minminZ = DBL_MAX; 671 minminZ = DBL_MAX;
574 maxmaxX = -DBL_MAX; 672 maxmaxX = -DBL_MAX;
575 maxmaxY = -DBL_MAX; 673 maxmaxY = -DBL_MAX;
576 maxmaxZ = -DBL_MAX; 674 maxmaxZ = -DBL_MAX;
577 675
578 while (moleculeListTemp) { << 676 while(moleculeListTemp)
579 if (minminX > moleculeListTemp->fMaxGlobX) << 677 {
>> 678 if(minminX > moleculeListTemp->fMaxGlobX)
>> 679 {
580 minminX = moleculeListTemp->fMaxGlobX; 680 minminX = moleculeListTemp->fMaxGlobX;
581 } 681 }
582 if (minminY > moleculeListTemp->fMaxGlobY) << 682 if(minminY > moleculeListTemp->fMaxGlobY)
>> 683 {
583 minminY = moleculeListTemp->fMaxGlobY; 684 minminY = moleculeListTemp->fMaxGlobY;
584 } 685 }
585 if (minminZ > moleculeListTemp->fMaxGlobZ) << 686 if(minminZ > moleculeListTemp->fMaxGlobZ)
>> 687 {
586 minminZ = moleculeListTemp->fMaxGlobZ; 688 minminZ = moleculeListTemp->fMaxGlobZ;
587 } 689 }
588 if (maxmaxX < moleculeListTemp->fMinGlobX) << 690 if(maxmaxX < moleculeListTemp->fMinGlobX)
>> 691 {
589 maxmaxX = moleculeListTemp->fMinGlobX; 692 maxmaxX = moleculeListTemp->fMinGlobX;
590 } 693 }
591 if (maxmaxY < moleculeListTemp->fMinGlobY) << 694 if(maxmaxY < moleculeListTemp->fMinGlobY)
>> 695 {
592 maxmaxY = moleculeListTemp->fMinGlobY; 696 maxmaxY = moleculeListTemp->fMinGlobY;
593 } 697 }
594 if (maxmaxZ < moleculeListTemp->fMinGlobZ) << 698 if(maxmaxZ < moleculeListTemp->fMinGlobZ)
>> 699 {
595 maxmaxZ = moleculeListTemp->fMinGlobZ; 700 maxmaxZ = moleculeListTemp->fMinGlobZ;
596 } 701 }
597 702
598 moleculeListTemp = moleculeListTemp->GetNe 703 moleculeListTemp = moleculeListTemp->GetNext();
599 } // end of while sur moleculeListTemp << 704 } //end of while sur moleculeListTemp
600 705
601 dX = (maxmaxX - minminX) / 2. + 1.8; // 1.8 << 706 dX = (maxmaxX - minminX) / 2. + 1.8; //1.8 => size of biggest radius for atoms
602 dY = (maxmaxY - minminY) / 2. + 1.8; 707 dY = (maxmaxY - minminY) / 2. + 1.8;
603 dZ = (maxmaxZ - minminZ) / 2. + 1.8; 708 dZ = (maxmaxZ - minminZ) / 2. + 1.8;
604 709
605 tX = minminX + (maxmaxX - minminX) / 2.; 710 tX = minminX + (maxmaxX - minminX) / 2.;
606 tY = minminY + (maxmaxY - minminY) / 2.; 711 tY = minminY + (maxmaxY - minminY) / 2.;
607 tZ = minminZ + (maxmaxZ - minminZ) / 2.; 712 tZ = minminZ + (maxmaxZ - minminZ) / 2.;
608 } 713 }
609 714
610 //....oooOO0OOooo........oooOO0OOooo........oo 715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
611 716
612 /** 717 /**
613 * \brief Compute number of nucleotide per 718 * \brief Compute number of nucleotide per strand
614 * \details Compute number of nucleotide per << 719 * \details Compute number of nucleotide per strand
615 * for DNA 720 * for DNA
616 */ 721 */
617 void PDBlib::ComputeNbNucleotidsPerStrand(Mole << 722 void PDBlib::ComputeNbNucleotidsPerStrand(Molecule * moleculeListTemp)
618 { 723 {
619 Residue* residueListTemp; << 724 Residue *residueListTemp;
620 725
621 int j_old = 0; 726 int j_old = 0;
622 727
623 while (moleculeListTemp) { << 728 while(moleculeListTemp)
>> 729 {
624 residueListTemp = moleculeListTemp->GetFir 730 residueListTemp = moleculeListTemp->GetFirst();
625 731
626 int j = 0; 732 int j = 0;
627 733
628 while (residueListTemp) { << 734 while(residueListTemp)
>> 735 {
629 j++; 736 j++;
630 residueListTemp = residueListTemp->GetNe 737 residueListTemp = residueListTemp->GetNext();
631 } // end of while sur residueListTemp << 738 } //end of while sur residueListTemp
632 739
633 j_old += j; 740 j_old += j;
634 moleculeListTemp = moleculeListTemp->GetNe 741 moleculeListTemp = moleculeListTemp->GetNext();
635 } // end of while sur moleculeListTemp << 742 } //end of while sur moleculeListTemp
636 743
637 fNbNucleotidsPerStrand = j_old / 2; 744 fNbNucleotidsPerStrand = j_old / 2;
638 } 745 }
639 746
640 //....oooOO0OOooo........oooOO0OOooo........oo 747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
641 748
642 /** 749 /**
643 * \brief Compute barycenters 750 * \brief Compute barycenters
644 * \details Compute barycenters and its coord << 751 * \details Compute barycenters and its coordinate
645 * for nucleotides 752 * for nucleotides
646 */ 753 */
647 unsigned short int PDBlib::ComputeMatchEdepDNA << 754 unsigned short int PDBlib::ComputeMatchEdepDNA(Barycenter *BarycenterList,
648 << 755 Molecule *moleculeListTemp,
649 << 756 double x,
650 << 757 double y,
>> 758 double z,
>> 759 int &numStrand,
>> 760 int &numNucleotid,
>> 761 int &codeResidue)
651 { 762 {
652 unsigned short int matchFound = 0; 763 unsigned short int matchFound = 0;
653 764
654 Molecule* mLTsavedPointer = moleculeListTemp << 765 Molecule *mLTsavedPointer = moleculeListTemp;
655 Barycenter* BLsavedPointer = BarycenterList; << 766 Barycenter *BLsavedPointer = BarycenterList;
656 767
657 short int strandNum = 0; // Strand number << 768 short int strandNum = 0; //Strand number
658 int residueNum = 1; // Residue (nucleotide) << 769 int residueNum = 1; //Residue (nucleotide) number
659 G4String baseName; // Base name [A,C,T,G] << 770 G4String baseName; //Base name [A,C,T,G]
660 unsigned short int BSP = 2; // Base (defaul << 771 unsigned short int BSP = 2; //Base (default value), Sugar, Phosphat
661 772
662 double smallestDist; // smallest dist Atom << 773 double smallestDist; //smallest dist Atom <-> edep coordinates
663 double distEdepDNA; 774 double distEdepDNA;
664 double distEdepAtom; 775 double distEdepAtom;
665 776
666 // Residue (Base, Phosphate,suggar) list << 777 //Residue (Base, Phosphate,suggar) list
667 Residue* residueListTemp; << 778 Residue *residueListTemp;
668 // Atom list inside a residue << 779 //Atom list inside a residue
669 Atom* AtomTemp; << 780 Atom *AtomTemp;
670 781
671 int k = 0; // Molecule number << 782 int k = 0; //Molecule number
672 moleculeListTemp = mLTsavedPointer; 783 moleculeListTemp = mLTsavedPointer;
673 BarycenterList = BLsavedPointer; 784 BarycenterList = BLsavedPointer;
674 785
675 smallestDist = 33.0; // Sufficiently large << 786 smallestDist = 33.0; //Sufficiently large value
676 while (moleculeListTemp) { << 787 while(moleculeListTemp)
>> 788 {
677 k++; 789 k++;
678 residueListTemp = moleculeListTemp->GetFir 790 residueListTemp = moleculeListTemp->GetFirst();
679 791
680 int j = 0; // Residue number << 792 int j = 0; //Residue number
681 793
682 int j_save = INT_MAX; // Saved res. numbe << 794 int j_save = INT_MAX; //Saved res. number if match found
683 795
684 while (residueListTemp) { << 796 while(residueListTemp)
>> 797 {
685 j++; 798 j++;
686 799
687 if (j - j_save > 2) break; << 800 if(j - j_save > 2) break;
688 801
689 distEdepDNA = DistanceTwo3Dpoints(x, Bar << 802 distEdepDNA = DistanceTwo3Dpoints(x,
>> 803 BarycenterList->fCenterX,
>> 804 y,
>> 805 BarycenterList->fCenterY,
>> 806 z,
690 Baryce 807 BarycenterList->fCenterZ);
691 if (distEdepDNA < BarycenterList->GetRad << 808 if(distEdepDNA < BarycenterList->GetRadius())
692 // Find the closest atom << 809 {
693 // Compute distance between energy dep << 810 //Find the closest atom
694 // if distance <1.8 then match OK but << 811 //Compute distance between energy deposited and atoms for a residue
>> 812 //if distance <1.8 then match OK but search inside 2 next residues
695 AtomTemp = residueListTemp->GetFirst() 813 AtomTemp = residueListTemp->GetFirst();
696 for (int iii = 0; iii < residueListTem << 814 for(int iii = 0; iii < residueListTemp->fNbAtom; iii++)
697 distEdepAtom = << 815 {
698 DistanceTwo3Dpoints(x, AtomTemp->G <<
699 816
700 if ((distEdepAtom < AtomTemp->GetVan << 817 distEdepAtom = DistanceTwo3Dpoints(x,
>> 818 AtomTemp->GetX(),
>> 819 y,
>> 820 AtomTemp->GetY(),
>> 821 z,
>> 822 AtomTemp->GetZ());
>> 823
>> 824 if((distEdepAtom < AtomTemp->GetVanDerWaalsRadius()) && (smallestDist
>> 825 > distEdepAtom))
>> 826 {
701 strandNum = k; 827 strandNum = k;
702 828
703 if (k == 2) { << 829 if(k == 2)
704 residueNum = fNbNucleotidsPerStr << 830 {
>> 831 residueNum = fNbNucleotidsPerStrand + 1
>> 832 - residueListTemp->fResSeq;
705 } 833 }
706 else { << 834 else
>> 835 {
707 residueNum = residueListTemp->fR 836 residueNum = residueListTemp->fResSeq;
708 } 837 }
709 838
710 baseName = (residueListTemp->fResN 839 baseName = (residueListTemp->fResName)[2];
711 if (residueListTemp->fResSeq == 1) << 840 if(residueListTemp->fResSeq == 1)
712 if (iii == 0) << 841 {
713 BSP = 0; //"Phosphate" << 842 if(iii == 0) BSP = 0; //"Phosphate"
714 else if (iii < 8) << 843 else if(iii < 8) BSP = 1; //"Sugar"
715 BSP = 1; //"Sugar" << 844 else BSP = 2; //"Base"
716 else <<
717 BSP = 2; //"Base" <<
718 } 845 }
719 else { << 846 else
720 if (iii < 4) << 847 {
721 BSP = 0; //"Phosphate" << 848 if(iii < 4) BSP = 0; //"Phosphate"
722 else if (iii < 11) << 849 else if(iii < 11) BSP = 1; //"Sugar"
723 BSP = 1; //"Sugar" << 850 else BSP = 2; //"Base"
724 else <<
725 BSP = 2; //"Base" <<
726 } 851 }
727 852
728 smallestDist = distEdepAtom; 853 smallestDist = distEdepAtom;
729 854
730 int j_max_value = INT_MAX; 855 int j_max_value = INT_MAX;
731 856
732 if (j_save == j_max_value) j_save << 857 if(j_save == j_max_value) j_save = j;
733 matchFound = 1; 858 matchFound = 1;
734 } 859 }
735 AtomTemp = AtomTemp->GetNext(); 860 AtomTemp = AtomTemp->GetNext();
736 } // end of for ( iii=0 ; iii < resi << 861 } //end of for ( iii=0 ; iii < residueListTemp->nbAtom ; iii++)
737 } // end for if (distEdepDNA < Barycent << 862 } //end for if (distEdepDNA < BarycenterList->GetRadius())
738 BarycenterList = BarycenterList->GetNext 863 BarycenterList = BarycenterList->GetNext();
739 residueListTemp = residueListTemp->GetNe 864 residueListTemp = residueListTemp->GetNext();
740 } // end of while sur residueListTemp << 865 } //end of while sur residueListTemp
741 moleculeListTemp = moleculeListTemp->GetNe 866 moleculeListTemp = moleculeListTemp->GetNext();
742 } // end of while sur moleculeListTemp << 867 } //end of while sur moleculeListTemp
743 868
744 numStrand = strandNum; 869 numStrand = strandNum;
745 numNucleotid = residueNum; 870 numNucleotid = residueNum;
746 codeResidue = BSP; 871 codeResidue = BSP;
747 872
748 return matchFound; 873 return matchFound;
749 } 874 }
750 875
751 //....oooOO0OOooo........oooOO0OOooo........oo 876 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
752 877
753 double PDBlib::DistanceTwo3Dpoints(double xA, << 878 double PDBlib::DistanceTwo3Dpoints(double xA,
>> 879 double xB,
>> 880 double yA,
>> 881 double yB,
>> 882 double zA,
>> 883 double zB)
754 { 884 {
755 return sqrt((xA - xB) * (xA - xB) + (yA - yB << 885 return sqrt((xA - xB) * (xA - xB) + (yA - yB) * (yA - yB)
>> 886 + (zA - zB) * (zA - zB));
756 } 887 }
757 888