74 unsigned short int verbose=0)
78 infile.open(filename.c_str());
82 G4cout<<
"PDBlib::load >> file "<<filename<<
" not found !!!!"<<
G4endl;
84 cout <<
"PDBlib::load >> file "<<filename<<
" not found !!!!"<<endl;
95 Atom * AtomicOld = NULL;
96 Atom * AtomicNext = NULL;
116 double minGlobZ,maxGlobZ;
117 double minGlobX,maxGlobX;
118 double minGlobY,maxGlobY;
119 double minX,maxX,minY,maxY,minZ,maxZ;
122 minGlobZ=-2000000000;
123 minGlobX=-2000000000;
124 minGlobY=-2000000000;
154 unsigned short int numModel=0;
155 unsigned short int model=0;
161 int terMax=2000000000;
167 getline(infile, sLine);
168 std::size_t found = sLine.find(
"DNA");
169 if (found!=std::string::npos)
177 found = sLine.find(
"HEADER");
178 if (found==std::string::npos)
181 infile.open(filename.c_str());
186 cout<<
"PDBlib::load >> No header found !!!!"<<endl;
192 while (!infile.eof())
194 getline(infile, sLine);
196 if ((sLine.substr(0,6)).compare(
"NUMMDL") == 0)
198 istringstream ((sLine.substr(10,4))) >> numModel;
200 if ((numModel > 0) && ( (sLine.substr(0,6)).compare(
"MODEL ") == 0))
202 istringstream ((sLine.substr(10,4))) >> model;
204 if (model > 1)
break;
208 if ((sLine.substr(0,6)).compare(
"SEQRES") == 0)
214 if (verbose > 1) cout << sLine << endl;
219 if ((numModel > 0) && ( (sLine.substr(0,6)).compare(
"ENDMDL") == 0))
222 else if ((sLine.substr(0,6)).compare(
"TER ") == 0)
227 if (ter > terMax)
break;
240 if (moleculeOld == NULL)
243 moleculeOld =
new Molecule(nameMol,nbMolecule);
244 moleculeOld->
SetFirst(residueFirst);
246 moleculeFirst = moleculeOld;
250 moleculeNext =
new Molecule(nameMol,nbMolecule);
251 moleculeOld->
SetNext(moleculeNext);
252 moleculeNext ->
SetFirst(residueFirst);
254 moleculeOld = moleculeNext;
259 moleculeOld->
fCenterX = (int)((minX + maxX) /2.);
260 moleculeOld->
fCenterY = (int)((minY + maxY) /2.);
261 moleculeOld->
fCenterZ = (int)((minZ + maxZ) /2.);
270 minGlobZ=-2000000000;
271 minGlobX=-2000000000;
272 minGlobY=-2000000000;
308 else if ((sLine.substr(0,6)).compare(
"ATOM ") == 0)
313 istringstream ((sLine.substr(6,5))) >> serial;
317 atomName=sLine.substr(12,4);
318 if (atomName.substr(0,1).compare(
" ") == 0 )
319 element=sLine.substr(13,1);
321 element=sLine.substr(12,1);
329 else if (element==
"C")
333 else if (element==
"N")
337 else if (element==
"O")
341 else if (element==
"P")
345 else if (element==
"S")
352 G4cout <<
"Element not recognized : " << element <<
G4endl;
355 cout <<
"Element not recognized : " << element << endl;
356 cout <<
"Stop now" << endl;
365 resName=sLine.substr(17,3);
367 istringstream ((sLine.substr(22,4))) >> resSeq;
369 stringstream ((sLine.substr(30,8))) >>
x;
370 stringstream ((sLine.substr(38,8))) >> y;
371 stringstream ((sLine.substr(46,8))) >>
z;
375 if (minGlobZ < z) minGlobZ=
z;
376 if (maxGlobZ > z) maxGlobZ=
z;
377 if (minGlobX < x) minGlobX=
x;
378 if (maxGlobX > x) maxGlobX=
x;
379 if (minGlobY < y) minGlobY=y;
380 if (maxGlobY > y) maxGlobY=y;
381 if (minX > x) minX=
x;
382 if (maxX < x) maxX=
x;
383 if (minY > y) minY=y;
384 if (maxY < y) maxY=y;
385 if (minZ > z) minZ=
z;
386 if (maxZ < z) maxZ=
z;
389 if (AtomicOld == NULL)
391 AtomicOld =
new Atom(serial,
405 AtomicNext =
new Atom(serial,
415 AtomicOld->
SetNext(AtomicNext);
416 AtomicOld = AtomicNext;
422 if (residueOld == NULL)
427 if (verbose>2) cout <<
"residueOld == NULL"<<endl;
430 residueOld =
new Residue(resName,resSeq);
433 residueFirst = residueOld;
439 if (lastResSeq==resSeq)
449 residueNext =
new Residue(resName,resSeq);
451 residueOld->
SetNext(residueNext);
455 residueOld = residueNext;
466 return moleculeFirst;
496 while (moleculeListTemp)
498 residueListTemp = moleculeListTemp->
GetFirst();
504 int correctNumerotationNumber=0;
505 if (k==2 && residueListTemp->
fResSeq > 1)
507 correctNumerotationNumber=residueListTemp->
fResSeq;
510 while (residueListTemp)
512 AtomTemp=residueListTemp->
GetFirst();
516 if (correctNumerotationNumber!=0)
519 correctNumerotationNumber+1;
523 double baryX=0.,baryY=0.,baryZ=0.;
524 double baryBaseX=0.,baryBaseY=0.,baryBaseZ=0.;
525 double barySugX=0.,barySugY=0.,barySugZ=0.;
526 double baryPhosX=0.,baryPhosY=0.,baryPhosZ=0.;
527 unsigned short int nbAtomInBase=0;
529 for (
int i=0 ; i < residueListTemp->
fNbAtom ; i++)
536 if (residueListTemp->
fResSeq == 1)
540 baryPhosX+=AtomTemp->
fX;
541 baryPhosY+=AtomTemp->
fY;
542 baryPhosZ+=AtomTemp->
fZ;
546 barySugX+=AtomTemp->
fX;
547 barySugY+=AtomTemp->
fY;
548 barySugZ+=AtomTemp->
fZ;
555 baryBaseX+=AtomTemp->
fX;
556 baryBaseY+=AtomTemp->
fY;
557 baryBaseZ+=AtomTemp->
fZ;
565 baryPhosX+=AtomTemp->
fX;
566 baryPhosY+=AtomTemp->
fY;
567 baryPhosZ+=AtomTemp->
fZ;
571 barySugX+=AtomTemp->
fX;
572 barySugY+=AtomTemp->
fY;
573 barySugZ+=AtomTemp->
fZ;
580 baryBaseX+=AtomTemp->
fX;
581 baryBaseY+=AtomTemp->
fY;
582 baryBaseZ+=AtomTemp->
fZ;
589 baryX = baryX / (double)residueListTemp->
fNbAtom;
590 baryY = baryY / (
double)residueListTemp->
fNbAtom;
591 baryZ = baryZ / (double)residueListTemp->
fNbAtom;
593 if (residueListTemp->
fResSeq != 1)
595 baryPhosX = baryPhosX / 4.;
596 baryPhosY = baryPhosY / 4.;
597 baryPhosZ = baryPhosZ / 4.;
599 barySugX = barySugX / 7.;
600 barySugY = barySugY / 7.;
601 barySugZ = barySugZ / 7.;
602 baryBaseX = baryBaseX / (double)nbAtomInBase;
603 baryBaseY = baryBaseY / (double)nbAtomInBase;
604 baryBaseZ = baryBaseZ / (double)nbAtomInBase;
607 if (BarycenterOld == NULL)
609 BarycenterOld =
new Barycenter(j+j_old,baryX,baryY,baryZ,
610 baryBaseX,baryBaseY,baryBaseZ,
611 barySugX,barySugY,barySugZ,
612 baryPhosX,baryPhosY,baryPhosZ);
613 BarycenterFirst = BarycenterOld;
617 BarycenterNext =
new Barycenter(j+j_old,baryX,baryY,baryZ,
618 baryBaseX,baryBaseY,baryBaseZ,
619 barySugX,barySugY,barySugZ,
620 baryPhosX,baryPhosY,baryPhosZ);
621 BarycenterOld -> SetNext(BarycenterNext);
622 BarycenterOld = BarycenterNext;
628 AtomTemp=residueListTemp->
GetFirst();
631 for (
int ii=0 ; ii < residueListTemp->
fNbAtom ; ii++)
637 if (dT3Dp>max) max=dT3Dp;
642 residueListTemp=residueListTemp->
GetNext();
649 moleculeListTemp=moleculeListTemp->
GetNext();
652 if(BarycenterNext!=NULL)
654 BarycenterNext -> SetNext(NULL);
657 return BarycenterFirst;
668 double &dX,
double &dY,
double &dZ,
669 double &tX,
double &tY,
double &tZ)
671 double minminX,minminY,minminZ;
672 double maxmaxX,maxmaxY,maxmaxZ;
690 while (moleculeListTemp)
692 if (minminX > moleculeListTemp->
fMaxGlobX)
696 if (minminY > moleculeListTemp->
fMaxGlobY)
700 if (minminZ > moleculeListTemp->
fMaxGlobZ)
704 if (maxmaxX < moleculeListTemp->fMinGlobX)
708 if (maxmaxY < moleculeListTemp->fMinGlobY)
712 if (maxmaxZ < moleculeListTemp->fMinGlobZ)
717 moleculeListTemp=moleculeListTemp->
GetNext();
720 dX=(maxmaxX-minminX)/2.+1.8;
721 dY=(maxmaxY-minminY)/2.+1.8;
722 dZ=(maxmaxZ-minminZ)/2.+1.8;
724 tX=minminX+(maxmaxX-minminX)/2.;
725 tY=minminY+(maxmaxY-minminY)/2.;
726 tZ=minminZ+(maxmaxZ-minminZ)/2.;
743 while (moleculeListTemp)
745 residueListTemp = moleculeListTemp->
GetFirst();
750 while (residueListTemp)
753 residueListTemp=residueListTemp->
GetNext();
757 moleculeListTemp=moleculeListTemp->
GetNext();
774 double x,
double y,
double z,
775 int &numStrand,
int &numNucleotid,
int &codeResidue)
777 unsigned short int matchFound=0;
779 Molecule *mLTsavedPointer = moleculeListTemp;
782 short int strandNum=0;
785 unsigned short int BSP=2;
797 moleculeListTemp = mLTsavedPointer;
798 BarycenterList = BLsavedPointer;
801 while (moleculeListTemp)
804 residueListTemp = moleculeListTemp->
GetFirst();
809 int j_save=2000000000;
814 while (residueListTemp)
818 if (j - j_save > 2 )
break;
823 if (distEdepDNA < BarycenterList->GetRadius())
828 AtomTemp=residueListTemp->
GetFirst();
829 for (
int iii=0 ; iii < residueListTemp->
fNbAtom ; iii++)
836 if ((distEdepAtom < AtomTemp->GetVanDerWaalsRadius())
837 && (smallestDist > distEdepAtom))
847 residueNum = residueListTemp->
fResSeq;
850 baseName = (residueListTemp->
fResName)[2];
851 if (residueListTemp->
fResSeq == 1)
853 if (iii == 0) BSP = 0;
854 else if (iii < 8) BSP = 1;
859 if (iii < 4) BSP = 0;
860 else if (iii < 11) BSP = 1;
864 smallestDist=distEdepAtom;
867 int j_max_value=2000000000;
871 if (j_save == j_max_value) j_save = j;
877 BarycenterList=BarycenterList->
GetNext();
878 residueListTemp=residueListTemp->
GetNext();
880 moleculeListTemp=moleculeListTemp->
GetNext();
883 numStrand = strandNum;
884 numNucleotid = residueNum;
895 return std::sqrt ( (xA-xB)*(xA-xB) + (yA-yB)*(yA-yB) + (zA-zB)*(zA-zB) );
int fCenterZ
"Z center" of this Molecule (for rotation...)
void SetRadius(double)
Set the distance between the farther atom and nucleotide barycenter.
double fCenterY
"Y coordinate" of this nucelotide Barycenter
double fCenterX
"X coordinate" of this nucelotide Barycenter
void SetNext(Molecule *)
Set the next Molecule.
int fNbAtom
Number of atoms into a residue (usage before vector)
Barycenter * ComputeNucleotideBarycenters(Molecule *moleculeListTemp)
Compute nucleotide barycenter from memory.
Residue * GetFirst()
Get the first Residue.
Residue * GetNext()
Get the next residue.
double fZ
Z orthogonal coordinates in Angstroms.
void SetFirst(Residue *)
Set the first Residue.
double GetX()
Return the X position for the Atom.
void SetNext(Atom *)
Set the next atom.
int fCenterX
"X center" of this Molecule (for rotation...)
double DistanceTwo3Dpoints(double xA, double xB, double yA, double yB, double zA, double zB)
return distance between two 3D points
void SetDistance(int i, double)
Set the distance between atom i and nucleotide barycenter.
int fNumInRes
its number in residue sequence
G4GLOB_DLL std::ostream G4cout
double GetZ()
Return the Z position for the Atom.
int fCenterY
"Y center" of this Molecule (for rotation...)
string fElement
Element symbol extracted from 'atom name'.
void ComputeNbNucleotidsPerStrand(Molecule *moleculeListTemp)
Compute number of nucleotide per strand.
unsigned short int ComputeMatchEdepDNA(Barycenter *, Molecule *, double x, double y, double z, int &numStrand, int &numNucleotid, int &codeResidue)
Compute if energy is deposited in per atom.
Molecule * GetNext()
information about molecule (not implemented)
Barycenter * GetNext()
Get the next Barycenter.
Molecule * Load(const string &filename, unsigned short int &isDNA, unsigned short int verbose)
Load PDB file into memory.
int fNbResidue
Number of residue inside the molecule.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
string fResName
Residue name.
Atom * GetFirst()
Get the first atom.
const G4double x[NPOINTSGL]
void SetNext(Residue *)
Set the next residue.
int fNbNucleotidsPerStrand
Number of nucleotid per strand.
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void ComputeBoundingVolumeParams(Molecule *moleculeListTemp, double &dX, double &dY, double &dZ, double &tX, double &tY, double &tZ)
Compute the corresponding bounding volume parameters.
double fX
X orthogonal coordinates in Angstroms.
Atom * GetNext()
Returns the next Atom.
double fY
Y orthogonal coordinates in Angstroms.
int fResSeq
Residue sequence number.
double fCenterZ
"Z coordinate" of this nucelotide Barycenter
PDBlib()
First constructor.
double GetY()
Return the Y position for the Atom.
Definition of the PDBlib class.
void SetFirst(Atom *)
Set the first Atom of the residue.