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