Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-04-11 08:04:39 

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