Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-06-24 07:53:28 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file PDBlib.cc
0027 /// \brief Implementation of the PDBlib class
0028 
0029 // This example is provided by the Geant4-DNA collaboration
0030 // Any report or published results obtained using the Geant4-DNA software
0031 // shall cite the following Geant4-DNA collaboration publication:
0032 // Med. Phys. 37 (2010) 4692-4708
0033 // Delage et al. PDB4DNA: implementation of DNA geometry from the Protein Data
0034 //                  Bank (PDB) description for Geant4-DNA Monte-Carlo
0035 //                  simulations (submitted to Comput. Phys. Commun.)
0036 // The Geant4-DNA web site is available at http://geant4-dna.org
0037 //
0038 //
0039 
0040 #include "PDBlib.hh"
0041 
0042 // define if the program is running with Geant4
0043 #define GEANT4
0044 
0045 #ifdef GEANT4
0046 // Specific to Geant4, globals.hh is used for G4cout
0047 #  include "globals.hh"
0048 #else
0049 #  define G4cout std::cout
0050 #  define G4cerr std::cerr
0051 #  define G4endl std::endl
0052 #  define G4String std::string  // string included in PDBatom.hh
0053 #  include <cfloat>
0054 #endif
0055 #include <cmath>
0056 #include <fstream>
0057 #include <iostream>
0058 #include <limits>
0059 #include <sstream>
0060 #include <stdlib.h>
0061 
0062 using namespace std;
0063 
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065 
0066 PDBlib::PDBlib() : fNbNucleotidsPerStrand(0) {}
0067 
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069 
0070 /**
0071  * \brief    Load a PDB file into memory
0072  * \details  molecule (polymer,?),
0073  *                  read line, key words
0074  * \param    filename    G4String for filename
0075  * \param    isDNA
0076  * \param    verbose
0077  * \return   List of molecules
0078  */
0079 
0080 Molecule* PDBlib::Load(const std::string& filename, unsigned short int& isDNA,
0081                        unsigned short int verbose = 0)
0082 {
0083   G4String sLine = "";
0084   ifstream infile;
0085   infile.open(filename.c_str());
0086   if (!infile) {
0087     G4cout << "PDBlib::load >> file " << filename << " not found !!!!" << G4endl;
0088   }
0089   else {
0090     int nbAtom = 0;  // total number of atoms in a residue
0091     int numAtomInRes = 0;  // number of an atom in a residue
0092     int nbResidue = 0;  // total number of residues
0093     int nbMolecule = 0;  // total number of molecule
0094 
0095     Atom* AtomicOld = nullptr;
0096     Atom* AtomicNext = nullptr;
0097 
0098     int serial;  // Atom serial number
0099     G4String atomName;  // Atom name
0100     G4String element;  // Element Symbol
0101     G4String resName;  // Residue name for this atom
0102     double x, y, z;  // Orthogonal coordinates in Angstroms
0103     double occupancy;  // Occupancy
0104 
0105     Residue* residueOld = nullptr;
0106     Residue* residueFirst = nullptr;
0107     Residue* residueNext = nullptr;
0108 
0109     Molecule* moleculeFirst = nullptr;
0110     Molecule* moleculeOld = nullptr;
0111     Molecule* moleculeNext = nullptr;
0112 
0113     /////////////////////////////////
0114     // NEW variable to draw a fitting cylinder if z oriented
0115     //=> fitting box
0116     double minGlobZ, maxGlobZ;
0117     double minGlobX, maxGlobX;
0118     double minGlobY, maxGlobY;
0119     double minX, maxX, minY, maxY, minZ, maxZ;  // Sort of 'mother volume' box
0120 
0121     minGlobZ = -DBL_MAX;
0122     minGlobX = -DBL_MAX;
0123     minGlobY = -DBL_MAX;
0124     maxGlobZ = DBL_MAX;
0125     maxGlobX = DBL_MAX;
0126     maxGlobY = DBL_MAX;
0127     minX = -DBL_MAX;
0128     minY = -DBL_MAX;
0129     minZ = -DBL_MAX;
0130     maxX = DBL_MAX;
0131     maxY = DBL_MAX;
0132     maxZ = DBL_MAX;
0133 
0134     int lastResSeq = -1;
0135     int resSeq = 0;
0136 
0137     G4String nameMol;
0138 
0139     unsigned short int numModel = 0;
0140     unsigned short int model = 0;
0141 
0142     // Number of TER
0143     int ter = 0;
0144     // Number of TER (chain) max for this file
0145 
0146     int terMax = INT_MAX;
0147 
0148     if (!infile.eof()) {
0149       getline(infile, sLine);
0150       std::size_t found = sLine.find("DNA");
0151       if (found != G4String::npos) {
0152         terMax = 2;
0153         isDNA = 1;
0154       }
0155       else
0156         isDNA = 0;
0157       // If PDB file have not a header line
0158       found = sLine.find("HEADER");
0159       if (found == G4String::npos) {
0160         infile.close();
0161         infile.open(filename.c_str());
0162 
0163         G4cout << "PDBlib::load >> No header found !!!!" << G4endl;
0164       }
0165     }
0166 
0167     while (!infile.eof()) {
0168       getline(infile, sLine);
0169       // In the case of several models
0170       if ((sLine.substr(0, 6)).compare("NUMMDL") == 0) {
0171         istringstream((sLine.substr(10, 4))) >> numModel;
0172       }
0173       if ((numModel > 0) && ((sLine.substr(0, 6)).compare("MODEL ") == 0)) {
0174         istringstream((sLine.substr(10, 4))) >> model;
0175         //////////////////////////////////////////
0176         if (model > 1) break;
0177         //////////////////////////////////////////
0178       }
0179       // For verification of residue sequence
0180       if ((sLine.substr(0, 6)).compare("SEQRES") == 0) {
0181         // Create list of molecule here
0182         if (verbose > 1) G4cout << sLine << G4endl;
0183       }
0184 
0185       // Coordinate section
0186       if ((numModel > 0) && ((sLine.substr(0, 6)).compare("ENDMDL") == 0)) {
0187         ;
0188       }
0189       else if ((sLine.substr(0, 6)).compare("TER   ") == 0)  // 3 spaces
0190       {
0191         //////////////////////////////////////////
0192         // Currently retrieve only the two first chains(TER) => two DNA strands
0193         ter++;
0194         if (ter > terMax) break;
0195         //////////////////////////////////////////
0196 
0197         // if (verbose>1) G4cout << sLine << G4endl;
0198         /************ Begin TER ******************/
0199         lastResSeq = -1;
0200         resSeq = 0;
0201 
0202         AtomicOld->SetNext(NULL);
0203         residueOld->SetNext(NULL);
0204         residueOld->fNbAtom = nbAtom;
0205 
0206         // Molecule creation:
0207         if (moleculeOld == NULL) {
0208           nameMol = filename;  //+numModel
0209           moleculeOld = new Molecule(nameMol, nbMolecule);
0210           moleculeOld->SetFirst(residueFirst);
0211           moleculeOld->fNbResidue = nbResidue;
0212           moleculeFirst = moleculeOld;
0213         }
0214         else {
0215           moleculeNext = new Molecule(nameMol, nbMolecule);
0216           moleculeOld->SetNext(moleculeNext);
0217           moleculeNext->SetFirst(residueFirst);
0218           moleculeNext->fNbResidue = nbResidue;
0219           moleculeOld = moleculeNext;
0220         }
0221 
0222         nbMolecule++;
0223         moleculeOld->SetNext(NULL);
0224         moleculeOld->fCenterX = (int)((minX + maxX) / 2.);
0225         moleculeOld->fCenterY = (int)((minY + maxY) / 2.);
0226         moleculeOld->fCenterZ = (int)((minZ + maxZ) / 2.);
0227         moleculeOld->fMaxGlobZ = maxGlobZ;
0228         moleculeOld->fMinGlobZ = minGlobZ;
0229         moleculeOld->fMaxGlobX = maxGlobX;
0230         moleculeOld->fMinGlobX = minGlobX;
0231         moleculeOld->fMaxGlobY = maxGlobY;
0232         moleculeOld->fMinGlobY = minGlobY;
0233 
0234         minGlobZ = -DBL_MAX;
0235         minGlobX = -DBL_MAX;
0236         minGlobY = -DBL_MAX;
0237         maxGlobZ = DBL_MAX;
0238         maxGlobX = DBL_MAX;
0239         maxGlobY = DBL_MAX;
0240         minX = -DBL_MAX;
0241         minY = -DBL_MAX;
0242         minZ = -DBL_MAX;
0243         maxX = DBL_MAX;
0244         maxY = DBL_MAX;
0245         maxZ = DBL_MAX;
0246 
0247         nbAtom = 0;
0248         numAtomInRes = 0;
0249         nbResidue = 0;
0250         AtomicOld = NULL;
0251         AtomicNext = NULL;
0252         residueOld = NULL;
0253         residueFirst = NULL;
0254         residueNext = NULL;
0255 
0256         ///////////// End TER ///////////////////
0257       }
0258       else if ((sLine.substr(0, 6)).compare("ATOM  ") == 0) {
0259         /************ Begin ATOM ******************/
0260         // serial
0261         istringstream((sLine.substr(6, 5))) >> serial;
0262 
0263         // To be improved
0264         // atomName :
0265         atomName = sLine.substr(12, 4);
0266         if (atomName.substr(0, 1).compare(" ") == 0)
0267           element = sLine.substr(13, 1);
0268         else
0269           element = sLine.substr(12, 1);
0270 
0271         // set Van der Waals radius expressed in Angstrom
0272         double vdwRadius = -1.;
0273         if (element == "H") {
0274           vdwRadius = 1.2;
0275         }
0276         else if (element == "C") {
0277           vdwRadius = 1.7;
0278         }
0279         else if (element == "N") {
0280           vdwRadius = 1.55;
0281         }
0282         else if (element == "O") {
0283           vdwRadius = 1.52;
0284         }
0285         else if (element == "P") {
0286           vdwRadius = 1.8;
0287         }
0288         else if (element == "S") {
0289           vdwRadius = 1.8;
0290         }
0291         else {
0292 #ifndef GEANT4
0293           G4cerr << "Element not recognized : " << element << G4endl;
0294           G4cerr << "Stop now" << G4endl;
0295           exit(1);
0296 #else
0297           G4ExceptionDescription errMsg;
0298           errMsg << "Element not recognized : " << element << G4endl;
0299 
0300           G4Exception("PDBlib::Load", "ELEM_NOT_RECOGNIZED", FatalException, errMsg);
0301 #endif
0302         }
0303 
0304         {
0305           // resName :
0306           resName = sLine.substr(17, 3);
0307           // resSeq :
0308           istringstream((sLine.substr(22, 4))) >> resSeq;
0309           // x,y,z :
0310           stringstream((sLine.substr(30, 8))) >> x;
0311           stringstream((sLine.substr(38, 8))) >> y;
0312           stringstream((sLine.substr(46, 8))) >> z;
0313           // occupancy :
0314           occupancy = 1.;
0315 
0316           if (minGlobZ < z) minGlobZ = z;
0317           if (maxGlobZ > z) maxGlobZ = z;
0318           if (minGlobX < x) minGlobX = x;
0319           if (maxGlobX > x) maxGlobX = x;
0320           if (minGlobY < y) minGlobY = y;
0321           if (maxGlobY > y) maxGlobY = y;
0322           if (minX > x) minX = x;
0323           if (maxX < x) maxX = x;
0324           if (minY > y) minY = y;
0325           if (maxY < y) maxY = y;
0326           if (minZ > z) minZ = z;
0327           if (maxZ < z) maxZ = z;
0328 
0329           // treatment for Atoms:
0330           if (AtomicOld == NULL) {
0331             AtomicOld =
0332               new Atom(serial, atomName, "", 0, 0, x, y, z, vdwRadius, occupancy, 0, element);
0333             AtomicOld->SetNext(NULL);  // If only one Atom inside residue
0334           }
0335           else {
0336             AtomicNext =
0337               new Atom(serial, atomName, "", 0, 0, x, y, z, vdwRadius, occupancy, 0, element);
0338             AtomicOld->SetNext(AtomicNext);
0339             AtomicOld = AtomicNext;
0340           }
0341           nbAtom++;
0342         }  // END if (element!="H")
0343         /****************************Begin Residue************************/
0344         // treatment for residues:
0345         if (residueOld == NULL) {
0346           if (verbose > 2) G4cout << "residueOld == NULL" << G4endl;
0347 
0348           AtomicOld->fNumInRes = 0;
0349           residueOld = new Residue(resName, resSeq);
0350           residueOld->SetFirst(AtomicOld);
0351           residueOld->SetNext(NULL);
0352           residueFirst = residueOld;
0353           lastResSeq = resSeq;
0354           nbResidue++;
0355         }
0356         else {
0357           if (lastResSeq == resSeq) {
0358             numAtomInRes++;
0359             AtomicOld->fNumInRes = numAtomInRes;
0360           }
0361           else {
0362             numAtomInRes = 0;
0363             AtomicOld->fNumInRes = numAtomInRes;
0364 
0365             residueNext = new Residue(resName, resSeq);
0366             residueNext->SetFirst(AtomicOld);
0367             residueOld->SetNext(residueNext);
0368             residueOld->fNbAtom = nbAtom - 1;
0369 
0370             nbAtom = 1;
0371             residueOld = residueNext;
0372             lastResSeq = resSeq;
0373             nbResidue++;
0374           }
0375         }
0376         /////////////////////////End Residue////////////
0377         ///////////// End ATOM ///////////////////
0378       }  // END if Atom
0379     }  // END while (!infile.eof())
0380 
0381     infile.close();
0382     return moleculeFirst;
0383   }  // END else if (!infile)
0384   return NULL;
0385 }
0386 
0387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0388 
0389 /**
0390  * \brief    Compute barycenters
0391  * \details  Compute barycenters and its coordinate
0392  *                  for nucleotides
0393  * \param    moleculeListTemp *    moleculeList
0394  * \return   Barycenter *
0395  */
0396 Barycenter* PDBlib::ComputeNucleotideBarycenters(Molecule* moleculeListTemp)
0397 {
0398   ///////////////////////////////////////////////////////
0399   // Placement and physical volume construction from memory
0400   Barycenter* BarycenterFirst = NULL;
0401   Barycenter* BarycenterOld = NULL;
0402   Barycenter* BarycenterNext = NULL;
0403 
0404   // Residue (Base, Phosphate,sugar) list
0405   Residue* residueListTemp;
0406   // Atom list inside a residu
0407   Atom* AtomTemp;
0408 
0409   int k = 0;
0410   int j_old = 0;
0411 
0412   while (moleculeListTemp) {
0413     residueListTemp = moleculeListTemp->GetFirst();
0414 
0415     k++;
0416     int j = 0;
0417 
0418     // Check numerotation style (1->n per strand or 1->2n for two strand)
0419     int correctNumerotationNumber = 0;
0420     if (k == 2 && residueListTemp->fResSeq > 1) {
0421       correctNumerotationNumber = residueListTemp->fResSeq;
0422     }
0423 
0424     while (residueListTemp) {
0425       AtomTemp = residueListTemp->GetFirst();
0426       j++;
0427 
0428       // Correction consequently to numerotation check
0429       if (correctNumerotationNumber != 0) {
0430         residueListTemp->fResSeq = residueListTemp->fResSeq - correctNumerotationNumber + 1;
0431       }
0432 
0433       // Barycenter computation
0434       double baryX = 0., baryY = 0., baryZ = 0.;
0435       double baryBaseX = 0., baryBaseY = 0., baryBaseZ = 0.;
0436       double barySugX = 0., barySugY = 0., barySugZ = 0.;
0437       double baryPhosX = 0., baryPhosY = 0., baryPhosZ = 0.;
0438       unsigned short int nbAtomInBase = 0;
0439 
0440       for (int i = 0; i < residueListTemp->fNbAtom; i++) {
0441         // Compute barycenter of the nucletotide
0442         baryX += AtomTemp->fX;
0443         baryY += AtomTemp->fY;
0444         baryZ += AtomTemp->fZ;
0445         // Compute barycenters for Base Sugar Phosphat
0446         if (residueListTemp->fResSeq == 1) {
0447           if (i == 0) {
0448             baryPhosX += AtomTemp->fX;
0449             baryPhosY += AtomTemp->fY;
0450             baryPhosZ += AtomTemp->fZ;
0451           }
0452           else if (i < 8) {
0453             barySugX += AtomTemp->fX;
0454             barySugY += AtomTemp->fY;
0455             barySugZ += AtomTemp->fZ;
0456           }
0457           else {
0458             // hydrogen are placed at the end of the residue in a PDB file
0459             // We don't want them for this calculation
0460             if (AtomTemp->fElement != "H") {
0461               baryBaseX += AtomTemp->fX;
0462               baryBaseY += AtomTemp->fY;
0463               baryBaseZ += AtomTemp->fZ;
0464               nbAtomInBase++;
0465             }
0466           }
0467         }
0468         else {
0469           if (i < 4) {
0470             baryPhosX += AtomTemp->fX;
0471             baryPhosY += AtomTemp->fY;
0472             baryPhosZ += AtomTemp->fZ;
0473           }
0474           else if (i < 11) {
0475             barySugX += AtomTemp->fX;
0476             barySugY += AtomTemp->fY;
0477             barySugZ += AtomTemp->fZ;
0478           }
0479           else {
0480             // hydrogen are placed at the end of the residue in a PDB file
0481             // We don't want them for this calculation
0482             if (AtomTemp->fElement != "H") {  // break;
0483               baryBaseX += AtomTemp->fX;
0484               baryBaseY += AtomTemp->fY;
0485               baryBaseZ += AtomTemp->fZ;
0486               nbAtomInBase++;
0487             }
0488           }
0489         }
0490         AtomTemp = AtomTemp->GetNext();
0491       }  // end of for (  i=0 ; i < residueListTemp->nbAtom ; i++)
0492 
0493       baryX = baryX / (double)residueListTemp->fNbAtom;
0494       baryY = baryY / (double)residueListTemp->fNbAtom;
0495       baryZ = baryZ / (double)residueListTemp->fNbAtom;
0496 
0497       if (residueListTemp->fResSeq != 1)  // Special case first Phosphate
0498       {
0499         baryPhosX = baryPhosX / 4.;
0500         baryPhosY = baryPhosY / 4.;
0501         baryPhosZ = baryPhosZ / 4.;
0502       }
0503       barySugX = barySugX / 7.;
0504       barySugY = barySugY / 7.;
0505       barySugZ = barySugZ / 7.;
0506       baryBaseX = baryBaseX / (double)nbAtomInBase;
0507       baryBaseY = baryBaseY / (double)nbAtomInBase;
0508       baryBaseZ = baryBaseZ / (double)nbAtomInBase;
0509 
0510       // Barycenter creation:
0511       if (BarycenterOld == NULL) {
0512         BarycenterOld = new Barycenter(j + j_old, baryX, baryY, baryZ,  // j [1..n]
0513                                        baryBaseX, baryBaseY, baryBaseZ, barySugX, barySugY,
0514                                        barySugZ, baryPhosX, baryPhosY, baryPhosZ);
0515         BarycenterFirst = BarycenterOld;
0516       }
0517       else {
0518         BarycenterNext =
0519           new Barycenter(j + j_old, baryX, baryY, baryZ, baryBaseX, baryBaseY, baryBaseZ, barySugX,
0520                          barySugY, barySugZ, baryPhosX, baryPhosY, baryPhosZ);
0521         BarycenterOld->SetNext(BarycenterNext);
0522         BarycenterOld = BarycenterNext;
0523       }
0524 
0525       /////////////////////////////////////////////////
0526       // distance computation between all atoms inside
0527       // a residue and the barycenter
0528       AtomTemp = residueListTemp->GetFirst();
0529       double dT3Dp;
0530       double max = 0.;
0531       for (int ii = 0; ii < residueListTemp->fNbAtom; ii++) {
0532         dT3Dp = DistanceTwo3Dpoints(AtomTemp->fX, BarycenterOld->fCenterX, AtomTemp->fY,
0533                                     BarycenterOld->fCenterY, AtomTemp->fZ, BarycenterOld->fCenterZ);
0534         BarycenterOld->SetDistance(ii, dT3Dp);
0535         if (dT3Dp > max) max = dT3Dp;
0536         AtomTemp = AtomTemp->GetNext();
0537       }  // end of for (  i=0 ; i < residueListTemp->nbAtom ; i++)
0538 
0539       BarycenterOld->SetRadius(max + 1.8);
0540       residueListTemp = residueListTemp->GetNext();
0541 
0542     }  // end of while sur residueListTemp
0543 
0544     j_old += j;
0545 
0546     /// molecs->push_back(*moleculeListTemp);
0547     moleculeListTemp = moleculeListTemp->GetNext();
0548   }  // end of while sur moleculeListTemp
0549 
0550   if (BarycenterNext != NULL) {
0551     BarycenterNext->SetNext(NULL);
0552   }
0553 
0554   return BarycenterFirst;
0555 }
0556 
0557 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0558 
0559 /**
0560  * \brief    the corresponding bounding volume parameters
0561  * \details  the corresponding bounding volume parameters
0562  *            to build a box from atoms coordinates
0563  */
0564 void PDBlib::ComputeBoundingVolumeParams(Molecule* moleculeListTemp, double& dX, double& dY,
0565                                          double& dZ,  // Dimensions for bounding volume
0566                                          double& tX, double& tY,
0567                                          double& tZ)  // Translation for bounding volume
0568 {
0569   double minminX, minminY, minminZ;  // minimum minimorum
0570   double maxmaxX, maxmaxY, maxmaxZ;  // maximum maximorum
0571 
0572   minminX = DBL_MAX;
0573   minminY = DBL_MAX;
0574   minminZ = DBL_MAX;
0575   maxmaxX = -DBL_MAX;
0576   maxmaxY = -DBL_MAX;
0577   maxmaxZ = -DBL_MAX;
0578 
0579   while (moleculeListTemp) {
0580     if (minminX > moleculeListTemp->fMaxGlobX) {
0581       minminX = moleculeListTemp->fMaxGlobX;
0582     }
0583     if (minminY > moleculeListTemp->fMaxGlobY) {
0584       minminY = moleculeListTemp->fMaxGlobY;
0585     }
0586     if (minminZ > moleculeListTemp->fMaxGlobZ) {
0587       minminZ = moleculeListTemp->fMaxGlobZ;
0588     }
0589     if (maxmaxX < moleculeListTemp->fMinGlobX) {
0590       maxmaxX = moleculeListTemp->fMinGlobX;
0591     }
0592     if (maxmaxY < moleculeListTemp->fMinGlobY) {
0593       maxmaxY = moleculeListTemp->fMinGlobY;
0594     }
0595     if (maxmaxZ < moleculeListTemp->fMinGlobZ) {
0596       maxmaxZ = moleculeListTemp->fMinGlobZ;
0597     }
0598 
0599     moleculeListTemp = moleculeListTemp->GetNext();
0600   }  // end of while sur moleculeListTemp
0601 
0602   dX = (maxmaxX - minminX) / 2. + 1.8;  // 1.8 => size of biggest radius for atoms
0603   dY = (maxmaxY - minminY) / 2. + 1.8;
0604   dZ = (maxmaxZ - minminZ) / 2. + 1.8;
0605 
0606   tX = minminX + (maxmaxX - minminX) / 2.;
0607   tY = minminY + (maxmaxY - minminY) / 2.;
0608   tZ = minminZ + (maxmaxZ - minminZ) / 2.;
0609 }
0610 
0611 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0612 
0613 /**
0614  * \brief    Compute number of nucleotide per strand
0615  * \details  Compute number of nucleotide per strand
0616  *                  for DNA
0617  */
0618 void PDBlib::ComputeNbNucleotidsPerStrand(Molecule* moleculeListTemp)
0619 {
0620   Residue* residueListTemp;
0621 
0622   int j_old = 0;
0623 
0624   while (moleculeListTemp) {
0625     residueListTemp = moleculeListTemp->GetFirst();
0626 
0627     int j = 0;
0628 
0629     while (residueListTemp) {
0630       j++;
0631       residueListTemp = residueListTemp->GetNext();
0632     }  // end of while sur residueListTemp
0633 
0634     j_old += j;
0635     moleculeListTemp = moleculeListTemp->GetNext();
0636   }  // end of while sur moleculeListTemp
0637 
0638   fNbNucleotidsPerStrand = j_old / 2;
0639 }
0640 
0641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0642 
0643 /**
0644  * \brief    Compute barycenters
0645  * \details  Compute barycenters and its coordinate
0646  *                  for nucleotides
0647  */
0648 unsigned short int PDBlib::ComputeMatchEdepDNA(Barycenter* BarycenterList,
0649                                                Molecule* moleculeListTemp, double x, double y,
0650                                                double z, int& numStrand, int& numNucleotid,
0651                                                int& codeResidue)
0652 {
0653   unsigned short int matchFound = 0;
0654 
0655   Molecule* mLTsavedPointer = moleculeListTemp;
0656   Barycenter* BLsavedPointer = BarycenterList;
0657 
0658   short int strandNum = 0;  // Strand number
0659   int residueNum = 1;  // Residue (nucleotide) number
0660   G4String baseName;  // Base name [A,C,T,G]
0661   unsigned short int BSP = 2;  // Base (default value), Sugar, Phosphat
0662 
0663   double smallestDist;  // smallest dist Atom <-> edep coordinates
0664   double distEdepDNA;
0665   double distEdepAtom;
0666 
0667   // Residue (Base, Phosphate,suggar) list
0668   Residue* residueListTemp;
0669   // Atom list inside a residue
0670   Atom* AtomTemp;
0671 
0672   int k = 0;  // Molecule number
0673   moleculeListTemp = mLTsavedPointer;
0674   BarycenterList = BLsavedPointer;
0675 
0676   smallestDist = 33.0;  // Sufficiently large value
0677   while (moleculeListTemp) {
0678     k++;
0679     residueListTemp = moleculeListTemp->GetFirst();
0680 
0681     int j = 0;  // Residue number
0682 
0683     int j_save = INT_MAX;  // Saved res. number if match found
0684 
0685     while (residueListTemp) {
0686       j++;
0687 
0688       if (j - j_save > 2) break;
0689 
0690       distEdepDNA = DistanceTwo3Dpoints(x, BarycenterList->fCenterX, y, BarycenterList->fCenterY, z,
0691                                         BarycenterList->fCenterZ);
0692       if (distEdepDNA < BarycenterList->GetRadius()) {
0693         // Find the closest atom
0694         // Compute distance between energy deposited and atoms for a residue
0695         // if distance <1.8 then match OK but search inside 2 next residues
0696         AtomTemp = residueListTemp->GetFirst();
0697         for (int iii = 0; iii < residueListTemp->fNbAtom; iii++) {
0698           distEdepAtom =
0699             DistanceTwo3Dpoints(x, AtomTemp->GetX(), y, AtomTemp->GetY(), z, AtomTemp->GetZ());
0700 
0701           if ((distEdepAtom < AtomTemp->GetVanDerWaalsRadius()) && (smallestDist > distEdepAtom)) {
0702             strandNum = k;
0703 
0704             if (k == 2) {
0705               residueNum = fNbNucleotidsPerStrand + 1 - residueListTemp->fResSeq;
0706             }
0707             else {
0708               residueNum = residueListTemp->fResSeq;
0709             }
0710 
0711             baseName = (residueListTemp->fResName)[2];
0712             if (residueListTemp->fResSeq == 1) {
0713               if (iii == 0)
0714                 BSP = 0;  //"Phosphate"
0715               else if (iii < 8)
0716                 BSP = 1;  //"Sugar"
0717               else
0718                 BSP = 2;  //"Base"
0719             }
0720             else {
0721               if (iii < 4)
0722                 BSP = 0;  //"Phosphate"
0723               else if (iii < 11)
0724                 BSP = 1;  //"Sugar"
0725               else
0726                 BSP = 2;  //"Base"
0727             }
0728 
0729             smallestDist = distEdepAtom;
0730 
0731             int j_max_value = INT_MAX;
0732 
0733             if (j_save == j_max_value) j_save = j;
0734             matchFound = 1;
0735           }
0736           AtomTemp = AtomTemp->GetNext();
0737         }  // end of for (  iii=0 ; iii < residueListTemp->nbAtom ; iii++)
0738       }  // end for if (distEdepDNA < BarycenterList->GetRadius())
0739       BarycenterList = BarycenterList->GetNext();
0740       residueListTemp = residueListTemp->GetNext();
0741     }  // end of while sur residueListTemp
0742     moleculeListTemp = moleculeListTemp->GetNext();
0743   }  // end of while sur moleculeListTemp
0744 
0745   numStrand = strandNum;
0746   numNucleotid = residueNum;
0747   codeResidue = BSP;
0748 
0749   return matchFound;
0750 }
0751 
0752 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0753 
0754 double PDBlib::DistanceTwo3Dpoints(double xA, double xB, double yA, double yB, double zA, double zB)
0755 {
0756   return sqrt((xA - xB) * (xA - xB) + (yA - yB) * (yA - yB) + (zA - zB) * (zA - zB));
0757 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!