| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // This example is provided by the Geant4-DNA
27 // Any report or published results obtained us
28 // shall cite the following Geant4-DNA collabo
29 // Med. Phys. 37 (2010) 4692-4708
30 // Delage et al. PDB4DNA: implementation of DN
31 // Bank (PDB) description for
32 // simulations (submitted to
33 // The Geant4-DNA web site is available at htt
34 //
35 //
36 /// \file PDBlib.cc
37 /// \brief Implementation file for PDBlib clas
38
39 #include "PDBlib.hh"
40
41 // define if the program is running with Geant
42 #define GEANT4
43
44 #ifdef GEANT4
45 // Specific to Geant4, globals.hh is used for
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 incl
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........oo
64
65 PDBlib::PDBlib() : fNbNucleotidsPerStrand(0) {
66
67 //....oooOO0OOooo........oooOO0OOooo........oo
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& file
80 unsigned short int verb
81 {
82 G4String sLine = "";
83 ifstream infile;
84 infile.open(filename.c_str());
85 if (!infile) {
86 G4cout << "PDBlib::load >> file " << filen
87 }
88 else {
89 int nbAtom = 0; // total number of atoms
90 int numAtomInRes = 0; // number of an ato
91 int nbResidue = 0; // total number of res
92 int nbMolecule = 0; // total number of mo
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 thi
101 double x, y, z; // Orthogonal coordinates
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
114 //=> fitting box
115 double minGlobZ, maxGlobZ;
116 double minGlobX, maxGlobX;
117 double minGlobY, maxGlobY;
118 double minX, maxX, minY, maxY, minZ, maxZ;
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 f
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
170 istringstream((sLine.substr(10, 4))) >
171 }
172 if ((numModel > 0) && ((sLine.substr(0,
173 istringstream((sLine.substr(10, 4))) >
174 //////////////////////////////////////
175 if (model > 1) break;
176 //////////////////////////////////////
177 }
178 // For verification of residue sequence
179 if ((sLine.substr(0, 6)).compare("SEQRES
180 // Create list of molecule here
181 if (verbose > 1) G4cout << sLine << G4
182 }
183
184 // Coordinate section
185 if ((numModel > 0) && ((sLine.substr(0,
186 ;
187 }
188 else if ((sLine.substr(0, 6)).compare("T
189 {
190 //////////////////////////////////////
191 // Currently retrieve only the two fir
192 ter++;
193 if (ter > terMax) break;
194 //////////////////////////////////////
195
196 // if (verbose>1) G4cout << sLine << G
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,
209 moleculeOld->SetFirst(residueFirst);
210 moleculeOld->fNbResidue = nbResidue;
211 moleculeFirst = moleculeOld;
212 }
213 else {
214 moleculeNext = new Molecule(nameMol,
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 +
224 moleculeOld->fCenterY = (int)((minY +
225 moleculeOld->fCenterZ = (int)((minZ +
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("A
258 /************ Begin ATOM *************
259 // serial
260 istringstream((sLine.substr(6, 5))) >>
261
262 // To be improved
263 // atomName :
264 atomName = sLine.substr(12, 4);
265 if (atomName.substr(0, 1).compare(" ")
266 element = sLine.substr(13, 1);
267 else
268 element = sLine.substr(12, 1);
269
270 // set Van der Waals radius expressed
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 :
293 G4cerr << "Stop now" << G4endl;
294 exit(1);
295 #else
296 G4ExceptionDescription errMsg;
297 errMsg << "Element not recognized :
298
299 G4Exception("PDBlib::Load", "ELEM_NO
300 #endif
301 }
302
303 {
304 // resName :
305 resName = sLine.substr(17, 3);
306 // resSeq :
307 istringstream((sLine.substr(22, 4)))
308 // x,y,z :
309 stringstream((sLine.substr(30, 8)))
310 stringstream((sLine.substr(38, 8)))
311 stringstream((sLine.substr(46, 8)))
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
332 AtomicOld->SetNext(NULL); // If o
333 }
334 else {
335 AtomicNext =
336 new Atom(serial, atomName, "", 0
337 AtomicOld->SetNext(AtomicNext);
338 AtomicOld = AtomicNext;
339 }
340 nbAtom++;
341 } // END if (element!="H")
342 /****************************Begin Res
343 // treatment for residues:
344 if (residueOld == NULL) {
345 if (verbose > 2) G4cout << "residueO
346
347 AtomicOld->fNumInRes = 0;
348 residueOld = new Residue(resName, re
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 = numAtomInRe
359 }
360 else {
361 numAtomInRes = 0;
362 AtomicOld->fNumInRes = numAtomInRe
363
364 residueNext = new Residue(resName,
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........oo
387
388 /**
389 * \brief Compute barycenters
390 * \details Compute barycenters and its coord
391 * for nucleotides
392 * \param moleculeListTemp * moleculeLis
393 * \return Barycenter *
394 */
395 Barycenter* PDBlib::ComputeNucleotideBarycente
396 {
397 ////////////////////////////////////////////
398 // Placement and physical volume constructio
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->GetFir
413
414 k++;
415 int j = 0;
416
417 // Check numerotation style (1->n per stra
418 int correctNumerotationNumber = 0;
419 if (k == 2 && residueListTemp->fResSeq > 1
420 correctNumerotationNumber = residueListT
421 }
422
423 while (residueListTemp) {
424 AtomTemp = residueListTemp->GetFirst();
425 j++;
426
427 // Correction consequently to numerotati
428 if (correctNumerotationNumber != 0) {
429 residueListTemp->fResSeq = residueList
430 }
431
432 // Barycenter computation
433 double baryX = 0., baryY = 0., baryZ = 0
434 double baryBaseX = 0., baryBaseY = 0., b
435 double barySugX = 0., barySugY = 0., bar
436 double baryPhosX = 0., baryPhosY = 0., b
437 unsigned short int nbAtomInBase = 0;
438
439 for (int i = 0; i < residueListTemp->fNb
440 // Compute barycenter of the nucletoti
441 baryX += AtomTemp->fX;
442 baryY += AtomTemp->fY;
443 baryZ += AtomTemp->fZ;
444 // Compute barycenters for Base Sugar
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
458 // We don't want them for this cal
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
480 // We don't want them for this cal
481 if (AtomTemp->fElement != "H") {
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 < residueLis
491
492 baryX = baryX / (double)residueListTemp-
493 baryY = baryY / (double)residueListTemp-
494 baryZ = baryZ / (double)residueListTemp-
495
496 if (residueListTemp->fResSeq != 1) // S
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)nbAtomIn
506 baryBaseY = baryBaseY / (double)nbAtomIn
507 baryBaseZ = baryBaseZ / (double)nbAtomIn
508
509 // Barycenter creation:
510 if (BarycenterOld == NULL) {
511 BarycenterOld = new Barycenter(j + j_o
512 baryBas
513 barySug
514 BarycenterFirst = BarycenterOld;
515 }
516 else {
517 BarycenterNext =
518 new Barycenter(j + j_old, baryX, bar
519 barySugY, barySugZ, b
520 BarycenterOld->SetNext(BarycenterNext)
521 BarycenterOld = BarycenterNext;
522 }
523
524 ////////////////////////////////////////
525 // distance computation between all atom
526 // a residue and the barycenter
527 AtomTemp = residueListTemp->GetFirst();
528 double dT3Dp;
529 double max = 0.;
530 for (int ii = 0; ii < residueListTemp->f
531 dT3Dp = DistanceTwo3Dpoints(AtomTemp->
532 Barycenter
533 BarycenterOld->SetDistance(ii, dT3Dp);
534 if (dT3Dp > max) max = dT3Dp;
535 AtomTemp = AtomTemp->GetNext();
536 } // end of for ( i=0 ; i < residueLis
537
538 BarycenterOld->SetRadius(max + 1.8);
539 residueListTemp = residueListTemp->GetNe
540
541 } // end of while sur residueListTemp
542
543 j_old += j;
544
545 /// molecs->push_back(*moleculeListTemp);
546 moleculeListTemp = moleculeListTemp->GetNe
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........oo
557
558 /**
559 * \brief the corresponding bounding volume
560 * \details the corresponding bounding volume
561 * to build a box from atoms coordi
562 */
563 void PDBlib::ComputeBoundingVolumeParams(Molec
564 doubl
565 doubl
566 doubl
567 {
568 double minminX, minminY, minminZ; // minimu
569 double maxmaxX, maxmaxY, maxmaxZ; // maximu
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->GetNe
599 } // end of while sur moleculeListTemp
600
601 dX = (maxmaxX - minminX) / 2. + 1.8; // 1.8
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........oo
611
612 /**
613 * \brief Compute number of nucleotide per
614 * \details Compute number of nucleotide per
615 * for DNA
616 */
617 void PDBlib::ComputeNbNucleotidsPerStrand(Mole
618 {
619 Residue* residueListTemp;
620
621 int j_old = 0;
622
623 while (moleculeListTemp) {
624 residueListTemp = moleculeListTemp->GetFir
625
626 int j = 0;
627
628 while (residueListTemp) {
629 j++;
630 residueListTemp = residueListTemp->GetNe
631 } // end of while sur residueListTemp
632
633 j_old += j;
634 moleculeListTemp = moleculeListTemp->GetNe
635 } // end of while sur moleculeListTemp
636
637 fNbNucleotidsPerStrand = j_old / 2;
638 }
639
640 //....oooOO0OOooo........oooOO0OOooo........oo
641
642 /**
643 * \brief Compute barycenters
644 * \details Compute barycenters and its coord
645 * for nucleotides
646 */
647 unsigned short int PDBlib::ComputeMatchEdepDNA
648
649
650
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)
659 G4String baseName; // Base name [A,C,T,G]
660 unsigned short int BSP = 2; // Base (defaul
661
662 double smallestDist; // smallest dist Atom
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
676 while (moleculeListTemp) {
677 k++;
678 residueListTemp = moleculeListTemp->GetFir
679
680 int j = 0; // Residue number
681
682 int j_save = INT_MAX; // Saved res. numbe
683
684 while (residueListTemp) {
685 j++;
686
687 if (j - j_save > 2) break;
688
689 distEdepDNA = DistanceTwo3Dpoints(x, Bar
690 Baryce
691 if (distEdepDNA < BarycenterList->GetRad
692 // Find the closest atom
693 // Compute distance between energy dep
694 // if distance <1.8 then match OK but
695 AtomTemp = residueListTemp->GetFirst()
696 for (int iii = 0; iii < residueListTem
697 distEdepAtom =
698 DistanceTwo3Dpoints(x, AtomTemp->G
699
700 if ((distEdepAtom < AtomTemp->GetVan
701 strandNum = k;
702
703 if (k == 2) {
704 residueNum = fNbNucleotidsPerStr
705 }
706 else {
707 residueNum = residueListTemp->fR
708 }
709
710 baseName = (residueListTemp->fResN
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
733 matchFound = 1;
734 }
735 AtomTemp = AtomTemp->GetNext();
736 } // end of for ( iii=0 ; iii < resi
737 } // end for if (distEdepDNA < Barycent
738 BarycenterList = BarycenterList->GetNext
739 residueListTemp = residueListTemp->GetNe
740 } // end of while sur residueListTemp
741 moleculeListTemp = moleculeListTemp->GetNe
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........oo
752
753 double PDBlib::DistanceTwo3Dpoints(double xA,
754 {
755 return sqrt((xA - xB) * (xA - xB) + (yA - yB
756 }
757