ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PDBlib.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PDBlib.cc
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 //
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 <fstream>
55 #include <iostream>
56 #include <limits>
57 #include <cmath>
58 #include <sstream>
59 #include <stdlib.h>
60 
61 using namespace std;
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
65 PDBlib::PDBlib() :
66  fNbNucleotidsPerStrand(0)
67 {
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
82 Molecule * PDBlib::Load(const std::string& filename,
83  unsigned short int& isDNA,
84  unsigned short int verbose = 0)
85 {
86  G4String sLine = "";
87  ifstream infile;
88  infile.open(filename.c_str());
89  if(!infile)
90  {
91  G4cout << "PDBlib::load >> file " << filename << " not found !!!!"
92  << G4endl;
93  }
94  else
95  {
96  int nbAtomTot = 0; // total number of atoms
97  int nbAtom = 0; // total number of atoms in a residue
98  int numAtomInRes = 0; // number of an atom in a residue
99  int nbResidue = 0; // total number of residues
100  int nbMolecule = 0; // total number of molecule
101 
102  Atom * AtomicOld = nullptr;
103  Atom * AtomicNext = nullptr;
104 
105  int serial; //Atom serial number
106  G4String atomName; //Atom name
107  G4String element; //Element Symbol
108  G4String resName; //Residue name for this atom
109  double x, y, z; //Orthogonal coordinates in Angstroms
110  double occupancy; //Occupancy
111 
112  Residue * residueOld = nullptr;
113  Residue * residueFirst = nullptr;
114  Residue * residueNext = nullptr;
115 
116  Molecule * moleculeFirst = nullptr;
117  Molecule * moleculeOld = nullptr;
118  Molecule * moleculeNext = nullptr;
119 
121  //NEW variable to draw a fitting cylinder if z oriented
122  //=> fitting box
123  double minGlobZ, maxGlobZ;
124  double minGlobX, maxGlobX;
125  double minGlobY, maxGlobY;
126  double minX, maxX, minY, maxY, minZ, maxZ; //Sort of 'mother volume' box
127 
128  minGlobZ = -DBL_MAX;
129  minGlobX = -DBL_MAX;
130  minGlobY = -DBL_MAX;
131  maxGlobZ = DBL_MAX;
132  maxGlobX = DBL_MAX;
133  maxGlobY = DBL_MAX;
134  minX = -DBL_MAX;
135  minY = -DBL_MAX;
136  minZ = -DBL_MAX;
137  maxX = DBL_MAX;
138  maxY = DBL_MAX;
139  maxZ = DBL_MAX;
140 
141  int lastResSeq = -1;
142  int resSeq = 0;
143 
144  G4String nameMol;
145 
146  unsigned short int numModel = 0;
147  unsigned short int model = 0;
148 
149  //Number of TER
150  int ter = 0;
151  //Number of TER (chain) max for this file
152 
153  int terMax = INT_MAX;
154 
155  if(!infile.eof())
156  {
157  getline(infile, sLine);
158  std::size_t found = sLine.find("DNA");
159  if(found != G4String::npos)
160  {
161  terMax = 2;
162  isDNA = 1;
163  }
164  else isDNA = 0;
165  //If PDB file have not a header line
166  found = sLine.find("HEADER");
167  if(found == G4String::npos)
168  {
169  infile.close();
170  infile.open(filename.c_str());
171 
172  G4cout << "PDBlib::load >> No header found !!!!" << G4endl;
173  }
174  }
175 
176  while(!infile.eof())
177  {
178  getline(infile, sLine);
179  //In the case of several models
180  if((sLine.substr(0, 6)).compare("NUMMDL") == 0)
181  {
182  istringstream((sLine.substr(10, 4))) >> numModel;
183  }
184  if((numModel > 0) && ((sLine.substr(0, 6)).compare("MODEL ") == 0))
185  {
186  istringstream((sLine.substr(10, 4))) >> model;
188  if(model > 1) break;
190  }
191  //For verification of residue sequence
192  if((sLine.substr(0, 6)).compare("SEQRES") == 0)
193  {
194  //Create list of molecule here
195  if(verbose > 1) G4cout << sLine << G4endl;
196  }
197 
198  //Coordinate section
199  if((numModel > 0) && ((sLine.substr(0, 6)).compare("ENDMDL") == 0))
200  {
201  ;
202  }
203  else if((sLine.substr(0, 6)).compare("TER ") == 0) //3 spaces
204  {
206  //Currently retrieve only the two first chains(TER) => two DNA strands
207  ter++;
208  if(ter > terMax) break;
210 
211  //if (verbose>1) G4cout << sLine << G4endl;
212  /************ Begin TER ******************/
213  lastResSeq = -1;
214  resSeq = 0;
215 
216  AtomicOld->SetNext(NULL);
217  residueOld->SetNext(NULL);
218  residueOld->fNbAtom = nbAtom;
219 
220  //Molecule creation:
221  if(moleculeOld == NULL)
222  {
223  nameMol = filename; //+numModel
224  moleculeOld = new Molecule(nameMol, nbMolecule);
225  moleculeOld->SetFirst(residueFirst);
226  moleculeOld->fNbResidue = nbResidue;
227  moleculeFirst = moleculeOld;
228  }
229  else
230  {
231  moleculeNext = new Molecule(nameMol, nbMolecule);
232  moleculeOld->SetNext(moleculeNext);
233  moleculeNext->SetFirst(residueFirst);
234  moleculeNext->fNbResidue = nbResidue;
235  moleculeOld = moleculeNext;
236  }
237 
238  nbMolecule++;
239  moleculeOld->SetNext(NULL);
240  moleculeOld->fCenterX = (int) ((minX + maxX) / 2.);
241  moleculeOld->fCenterY = (int) ((minY + maxY) / 2.);
242  moleculeOld->fCenterZ = (int) ((minZ + maxZ) / 2.);
243  moleculeOld->fMaxGlobZ = maxGlobZ;
244  moleculeOld->fMinGlobZ = minGlobZ;
245  moleculeOld->fMaxGlobX = maxGlobX;
246  moleculeOld->fMinGlobX = minGlobX;
247  moleculeOld->fMaxGlobY = maxGlobY;
248  moleculeOld->fMinGlobY = minGlobY;
249 
250  minGlobZ = -DBL_MAX;
251  minGlobX = -DBL_MAX;
252  minGlobY = -DBL_MAX;
253  maxGlobZ = DBL_MAX;
254  maxGlobX = DBL_MAX;
255  maxGlobY = DBL_MAX;
256  minX = -DBL_MAX;
257  minY = -DBL_MAX;
258  minZ = -DBL_MAX;
259  maxX = DBL_MAX;
260  maxY = DBL_MAX;
261  maxZ = DBL_MAX;
262 
263  nbAtom = 0;
264  numAtomInRes = 0;
265  nbResidue = 0;
266  AtomicOld = NULL;
267  AtomicNext = NULL;
268  residueOld = NULL;
269  residueFirst = NULL;
270  residueNext = NULL;
271 
273  }
274  else if((sLine.substr(0, 6)).compare("ATOM ") == 0)
275  {
276 
277  /************ Begin ATOM ******************/
278  //serial
279  istringstream((sLine.substr(6, 5))) >> serial;
280 
281  //To be improved
282  //atomName :
283  atomName = sLine.substr(12, 4);
284  if(atomName.substr(0, 1).compare(" ") == 0) element = sLine.substr(13,
285  1);
286  else element = sLine.substr(12, 1);
287 
288  // set Van der Waals radius expressed in Angstrom
289  double vdwRadius = -1.;
290  if(element == "H")
291  {
292  vdwRadius = 1.2;
293  }
294  else if(element == "C")
295  {
296  vdwRadius = 1.7;
297  }
298  else if(element == "N")
299  {
300  vdwRadius = 1.55;
301  }
302  else if(element == "O")
303  {
304  vdwRadius = 1.52;
305  }
306  else if(element == "P")
307  {
308  vdwRadius = 1.8;
309  }
310  else if(element == "S")
311  {
312  vdwRadius = 1.8;
313  }
314  else
315  {
316 #ifndef GEANT4
317  G4cerr << "Element not recognized : " << element << G4endl;
318  G4cerr << "Stop now" << G4endl;
319  exit(1);
320 #else
321  G4ExceptionDescription errMsg;
322  errMsg << "Element not recognized : " << element << G4endl;
323 
324  G4Exception("PDBlib::Load",
325  "ELEM_NOT_RECOGNIZED",
326  FatalException,
327  errMsg);
328 #endif
329  }
330 
331  {
332  nbAtomTot++;
333 
334  //resName :
335  resName = sLine.substr(17, 3);
336  //resSeq :
337  istringstream((sLine.substr(22, 4))) >> resSeq;
338  //x,y,z :
339  stringstream((sLine.substr(30, 8))) >> x;
340  stringstream((sLine.substr(38, 8))) >> y;
341  stringstream((sLine.substr(46, 8))) >> z;
342  //occupancy :
343  occupancy = 1.;
344 
345  if(minGlobZ < z) minGlobZ = z;
346  if(maxGlobZ > z) maxGlobZ = z;
347  if(minGlobX < x) minGlobX = x;
348  if(maxGlobX > x) maxGlobX = x;
349  if(minGlobY < y) minGlobY = y;
350  if(maxGlobY > y) maxGlobY = y;
351  if(minX > x) minX = x;
352  if(maxX < x) maxX = x;
353  if(minY > y) minY = y;
354  if(maxY < y) maxY = y;
355  if(minZ > z) minZ = z;
356  if(maxZ < z) maxZ = z;
357 
358  //treatment for Atoms:
359  if(AtomicOld == NULL)
360  {
361  AtomicOld = new Atom(serial,
362  atomName,
363  "",
364  0,
365  0,
366  x,
367  y,
368  z,
369  vdwRadius,
370  occupancy,
371  0,
372  element);
373  AtomicOld->SetNext(NULL); //If only one Atom inside residue
374  }
375  else
376  {
377  AtomicNext = new Atom(serial,
378  atomName,
379  "",
380  0,
381  0,
382  x,
383  y,
384  z,
385  vdwRadius,
386  occupancy,
387  0,
388  element);
389  AtomicOld->SetNext(AtomicNext);
390  AtomicOld = AtomicNext;
391  }
392  nbAtom++;
393  } //END if (element!="H")
394  /****************************Begin Residue************************/
395  //treatment for residues:
396  if(residueOld == NULL)
397  {
398  if(verbose > 2) G4cout << "residueOld == NULL" << G4endl;
399 
400  AtomicOld->fNumInRes = 0;
401  residueOld = new Residue(resName, resSeq);
402  residueOld->SetFirst(AtomicOld);
403  residueOld->SetNext(NULL);
404  residueFirst = residueOld;
405  lastResSeq = resSeq;
406  nbResidue++;
407  }
408  else
409  {
410  if(lastResSeq == resSeq)
411  {
412  numAtomInRes++;
413  AtomicOld->fNumInRes = numAtomInRes;
414  }
415  else
416  {
417  numAtomInRes = 0;
418  AtomicOld->fNumInRes = numAtomInRes;
419 
420  residueNext = new Residue(resName, resSeq);
421  residueNext->SetFirst(AtomicOld);
422  residueOld->SetNext(residueNext);
423  residueOld->fNbAtom = nbAtom - 1;
424 
425  nbAtom = 1;
426  residueOld = residueNext;
427  lastResSeq = resSeq;
428  nbResidue++;
429  }
430  }
433  } //END if Atom
434  } //END while (!infile.eof())
435 
436  infile.close();
437  return moleculeFirst;
438  } //END else if (!infile)
439  return NULL;
440 }
441 
442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
443 
451 Barycenter * PDBlib::ComputeNucleotideBarycenters(Molecule * moleculeListTemp)
452 {
454  //Placement and physical volume construction from memory
455  Barycenter * BarycenterFirst = NULL;
456  Barycenter * BarycenterOld = NULL;
457  Barycenter * BarycenterNext = NULL;
458 
459  //Residue (Base, Phosphate,sugar) list
460  Residue *residueListTemp;
461  //Atom list inside a residu
462  Atom *AtomTemp;
463 
464  int k = 0;
465  int j_old = 0;
466 
467  while(moleculeListTemp)
468  {
469  residueListTemp = moleculeListTemp->GetFirst();
470 
471  k++;
472  int j = 0;
473 
474  //Check numerotation style (1->n per strand or 1->2n for two strand)
475  int correctNumerotationNumber = 0;
476  if(k == 2 && residueListTemp->fResSeq > 1)
477  {
478  correctNumerotationNumber = residueListTemp->fResSeq;
479  }
480 
481  while(residueListTemp)
482  {
483  AtomTemp = residueListTemp->GetFirst();
484  j++;
485 
486  //Correction consequently to numerotation check
487  if(correctNumerotationNumber != 0)
488  {
489  residueListTemp->fResSeq = residueListTemp->fResSeq
490  - correctNumerotationNumber
491  + 1;
492  }
493 
494  //Barycenter computation
495  double baryX = 0., baryY = 0., baryZ = 0.;
496  double baryBaseX = 0., baryBaseY = 0., baryBaseZ = 0.;
497  double barySugX = 0., barySugY = 0., barySugZ = 0.;
498  double baryPhosX = 0., baryPhosY = 0., baryPhosZ = 0.;
499  unsigned short int nbAtomInBase = 0;
500 
501  for(int i = 0; i < residueListTemp->fNbAtom; i++)
502  {
503  //Compute barycenter of the nucletotide
504  baryX += AtomTemp->fX;
505  baryY += AtomTemp->fY;
506  baryZ += AtomTemp->fZ;
507  //Compute barycenters for Base Sugar Phosphat
508  if(residueListTemp->fResSeq == 1)
509  {
510  if(i == 0)
511  {
512  baryPhosX += AtomTemp->fX;
513  baryPhosY += AtomTemp->fY;
514  baryPhosZ += AtomTemp->fZ;
515  }
516  else if(i < 8)
517  {
518  barySugX += AtomTemp->fX;
519  barySugY += AtomTemp->fY;
520  barySugZ += AtomTemp->fZ;
521  }
522  else
523  {
524  //hydrogen are placed at the end of the residue in a PDB file
525  //We don't want them for this calculation
526  if(AtomTemp->fElement != "H")
527  {
528  baryBaseX += AtomTemp->fX;
529  baryBaseY += AtomTemp->fY;
530  baryBaseZ += AtomTemp->fZ;
531  nbAtomInBase++;
532  }
533  }
534  }
535  else
536  {
537  if(i < 4)
538  {
539  baryPhosX += AtomTemp->fX;
540  baryPhosY += AtomTemp->fY;
541  baryPhosZ += AtomTemp->fZ;
542  }
543  else if(i < 11)
544  {
545  barySugX += AtomTemp->fX;
546  barySugY += AtomTemp->fY;
547  barySugZ += AtomTemp->fZ;
548  }
549  else
550  {
551  //hydrogen are placed at the end of the residue in a PDB file
552  //We don't want them for this calculation
553  if(AtomTemp->fElement != "H")
554  { // break;
555  baryBaseX += AtomTemp->fX;
556  baryBaseY += AtomTemp->fY;
557  baryBaseZ += AtomTemp->fZ;
558  nbAtomInBase++;
559  }
560  }
561  }
562  AtomTemp = AtomTemp->GetNext();
563  } //end of for ( i=0 ; i < residueListTemp->nbAtom ; i++)
564 
565  baryX = baryX / (double) residueListTemp->fNbAtom;
566  baryY = baryY / (double) residueListTemp->fNbAtom;
567  baryZ = baryZ / (double) residueListTemp->fNbAtom;
568 
569  if(residueListTemp->fResSeq != 1) //Special case first Phosphate
570  {
571  baryPhosX = baryPhosX / 4.;
572  baryPhosY = baryPhosY / 4.;
573  baryPhosZ = baryPhosZ / 4.;
574  }
575  barySugX = barySugX / 7.;
576  barySugY = barySugY / 7.;
577  barySugZ = barySugZ / 7.;
578  baryBaseX = baryBaseX / (double) nbAtomInBase;
579  baryBaseY = baryBaseY / (double) nbAtomInBase;
580  baryBaseZ = baryBaseZ / (double) nbAtomInBase;
581 
582  //Barycenter creation:
583  if(BarycenterOld == NULL)
584  {
585  BarycenterOld = new Barycenter(j + j_old, baryX, baryY, baryZ, //j [1..n]
586  baryBaseX,
587  baryBaseY,
588  baryBaseZ,
589  barySugX,
590  barySugY,
591  barySugZ,
592  baryPhosX,
593  baryPhosY,
594  baryPhosZ);
595  BarycenterFirst = BarycenterOld;
596  }
597  else
598  {
599  BarycenterNext = new Barycenter(j + j_old,
600  baryX,
601  baryY,
602  baryZ,
603  baryBaseX,
604  baryBaseY,
605  baryBaseZ,
606  barySugX,
607  barySugY,
608  barySugZ,
609  baryPhosX,
610  baryPhosY,
611  baryPhosZ);
612  BarycenterOld->SetNext(BarycenterNext);
613  BarycenterOld = BarycenterNext;
614  }
615 
617  //distance computation between all atoms inside
618  //a residue and the barycenter
619  AtomTemp = residueListTemp->GetFirst();
620  double dT3Dp;
621  double max = 0.;
622  for(int ii = 0; ii < residueListTemp->fNbAtom; ii++)
623  {
624  dT3Dp = DistanceTwo3Dpoints(AtomTemp->fX,
625  BarycenterOld->fCenterX,
626  AtomTemp->fY,
627  BarycenterOld->fCenterY,
628  AtomTemp->fZ,
629  BarycenterOld->fCenterZ);
630  BarycenterOld->SetDistance(ii, dT3Dp);
631  if(dT3Dp > max) max = dT3Dp;
632  AtomTemp = AtomTemp->GetNext();
633  } //end of for ( i=0 ; i < residueListTemp->nbAtom ; i++)
634 
635  BarycenterOld->SetRadius(max + 1.8);
636  residueListTemp = residueListTemp->GetNext();
637 
638  } //end of while sur residueListTemp
639 
640  j_old += j;
641 
643  moleculeListTemp = moleculeListTemp->GetNext();
644  } //end of while sur moleculeListTemp
645 
646  if(BarycenterNext != NULL)
647  {
648  BarycenterNext->SetNext(NULL);
649  }
650 
651  return BarycenterFirst;
652 }
653 
654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
655 
661 void PDBlib::ComputeBoundingVolumeParams(Molecule *moleculeListTemp,
662  double &dX,
663  double &dY,
664  double &dZ, //Dimensions for bounding volume
665  double &tX,
666  double &tY,
667  double &tZ) //Translation for bounding volume
668 {
669  double minminX, minminY, minminZ; //minimum minimorum
670  double maxmaxX, maxmaxY, maxmaxZ; //maximum maximorum
671 
672  minminX = DBL_MAX;
673  minminY = DBL_MAX;
674  minminZ = DBL_MAX;
675  maxmaxX = -DBL_MAX;
676  maxmaxY = -DBL_MAX;
677  maxmaxZ = -DBL_MAX;
678 
679  while(moleculeListTemp)
680  {
681  if(minminX > moleculeListTemp->fMaxGlobX)
682  {
683  minminX = moleculeListTemp->fMaxGlobX;
684  }
685  if(minminY > moleculeListTemp->fMaxGlobY)
686  {
687  minminY = moleculeListTemp->fMaxGlobY;
688  }
689  if(minminZ > moleculeListTemp->fMaxGlobZ)
690  {
691  minminZ = moleculeListTemp->fMaxGlobZ;
692  }
693  if(maxmaxX < moleculeListTemp->fMinGlobX)
694  {
695  maxmaxX = moleculeListTemp->fMinGlobX;
696  }
697  if(maxmaxY < moleculeListTemp->fMinGlobY)
698  {
699  maxmaxY = moleculeListTemp->fMinGlobY;
700  }
701  if(maxmaxZ < moleculeListTemp->fMinGlobZ)
702  {
703  maxmaxZ = moleculeListTemp->fMinGlobZ;
704  }
705 
706  moleculeListTemp = moleculeListTemp->GetNext();
707  } //end of while sur moleculeListTemp
708 
709  dX = (maxmaxX - minminX) / 2. + 1.8; //1.8 => size of biggest radius for atoms
710  dY = (maxmaxY - minminY) / 2. + 1.8;
711  dZ = (maxmaxZ - minminZ) / 2. + 1.8;
712 
713  tX = minminX + (maxmaxX - minminX) / 2.;
714  tY = minminY + (maxmaxY - minminY) / 2.;
715  tZ = minminZ + (maxmaxZ - minminZ) / 2.;
716 }
717 
718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
719 
725 void PDBlib::ComputeNbNucleotidsPerStrand(Molecule * moleculeListTemp)
726 {
727  Residue *residueListTemp;
728 
729  int k = 0;
730  int j_old = 0;
731 
732  while(moleculeListTemp)
733  {
734  residueListTemp = moleculeListTemp->GetFirst();
735 
736  k++;
737  int j = 0;
738 
739  while(residueListTemp)
740  {
741  j++;
742  residueListTemp = residueListTemp->GetNext();
743  } //end of while sur residueListTemp
744 
745  j_old += j;
746  moleculeListTemp = moleculeListTemp->GetNext();
747  } //end of while sur moleculeListTemp
748 
749  fNbNucleotidsPerStrand = j_old / 2;
750 }
751 
752 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
753 
759 unsigned short int PDBlib::ComputeMatchEdepDNA(Barycenter *BarycenterList,
760  Molecule *moleculeListTemp,
761  double x,
762  double y,
763  double z,
764  int &numStrand,
765  int &numNucleotid,
766  int &codeResidue)
767 {
768  unsigned short int matchFound = 0;
769 
770  Molecule *mLTsavedPointer = moleculeListTemp;
771  Barycenter *BLsavedPointer = BarycenterList;
772 
773  short int strandNum = 0; //Strand number
774  int residueNum = 1; //Residue (nucleotide) number
775  G4String baseName; //Base name [A,C,T,G]
776  unsigned short int BSP = 2; //Base (default value), Sugar, Phosphat
777 
778  double smallestDist; //smallest dist Atom <-> edep coordinates
779  double distEdepDNA;
780  double distEdepAtom;
781 
782  //Residue (Base, Phosphate,suggar) list
783  Residue *residueListTemp;
784  //Atom list inside a residue
785  Atom *AtomTemp;
786 
787  int k = 0; //Molecule number
788  moleculeListTemp = mLTsavedPointer;
789  BarycenterList = BLsavedPointer;
790 
791  smallestDist = 33.0; //Sufficiently large value
792  while(moleculeListTemp)
793  {
794  k++;
795  residueListTemp = moleculeListTemp->GetFirst();
796 
797  int j = 0; //Residue number
798 
799  int j_save = INT_MAX; //Saved res. number if match found
800 
801  while(residueListTemp)
802  {
803  j++;
804 
805  if(j - j_save > 2) break;
806 
807  distEdepDNA = DistanceTwo3Dpoints(x,
808  BarycenterList->fCenterX,
809  y,
810  BarycenterList->fCenterY,
811  z,
812  BarycenterList->fCenterZ);
813  if(distEdepDNA < BarycenterList->GetRadius())
814  {
815  //Find the closest atom
816  //Compute distance between energy deposited and atoms for a residue
817  //if distance <1.8 then match OK but search inside 2 next residues
818  AtomTemp = residueListTemp->GetFirst();
819  for(int iii = 0; iii < residueListTemp->fNbAtom; iii++)
820  {
821 
822  distEdepAtom = DistanceTwo3Dpoints(x,
823  AtomTemp->GetX(),
824  y,
825  AtomTemp->GetY(),
826  z,
827  AtomTemp->GetZ());
828 
829  if((distEdepAtom < AtomTemp->GetVanDerWaalsRadius()) && (smallestDist
830  > distEdepAtom))
831  {
832  strandNum = k;
833 
834  if(k == 2)
835  {
836  residueNum = fNbNucleotidsPerStrand + 1
837  - residueListTemp->fResSeq;
838  }
839  else
840  {
841  residueNum = residueListTemp->fResSeq;
842  }
843 
844  baseName = (residueListTemp->fResName)[2];
845  if(residueListTemp->fResSeq == 1)
846  {
847  if(iii == 0) BSP = 0; //"Phosphate"
848  else if(iii < 8) BSP = 1; //"Sugar"
849  else BSP = 2; //"Base"
850  }
851  else
852  {
853  if(iii < 4) BSP = 0; //"Phosphate"
854  else if(iii < 11) BSP = 1; //"Sugar"
855  else BSP = 2; //"Base"
856  }
857 
858  smallestDist = distEdepAtom;
859 
860  int j_max_value = INT_MAX;
861 
862  if(j_save == j_max_value) j_save = j;
863  matchFound = 1;
864  }
865  AtomTemp = AtomTemp->GetNext();
866  } //end of for ( iii=0 ; iii < residueListTemp->nbAtom ; iii++)
867  } //end for if (distEdepDNA < BarycenterList->GetRadius())
868  BarycenterList = BarycenterList->GetNext();
869  residueListTemp = residueListTemp->GetNext();
870  } //end of while sur residueListTemp
871  moleculeListTemp = moleculeListTemp->GetNext();
872  } //end of while sur moleculeListTemp
873 
874  numStrand = strandNum;
875  numNucleotid = residueNum;
876  codeResidue = BSP;
877 
878  return matchFound;
879 }
880 
881 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
882 
883 double PDBlib::DistanceTwo3Dpoints(double xA,
884  double xB,
885  double yA,
886  double yB,
887  double zA,
888  double zB)
889 {
890  return sqrt((xA - xB) * (xA - xB) + (yA - yB) * (yA - yB)
891  + (zA - zB) * (zA - zB));
892 }