Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // This example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration 27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software 28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication: 29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708 30 // Delage et al. PDB4DNA: implementation of DN 30 // Delage et al. PDB4DNA: implementation of DNA geometry from the Protein Data 31 // Bank (PDB) description for 31 // Bank (PDB) description for Geant4-DNA Monte-Carlo 32 // simulations (submitted to 32 // simulations (submitted to Comput. Phys. Commun.) 33 // The Geant4-DNA web site is available at htt 33 // The Geant4-DNA web site is available at http://geant4-dna.org 34 // 34 // 35 // << 35 // 36 /// \file PDBlib.cc 36 /// \file PDBlib.cc 37 /// \brief Implementation file for PDBlib clas 37 /// \brief Implementation file for PDBlib class 38 38 39 #include "PDBlib.hh" 39 #include "PDBlib.hh" 40 40 41 // define if the program is running with Geant << 41 //define if the program is running with Geant4 42 #define GEANT4 42 #define GEANT4 43 43 44 #ifdef GEANT4 44 #ifdef GEANT4 45 // Specific to Geant4, globals.hh is used for << 45 //Specific to Geant4, globals.hh is used for G4cout 46 # include "globals.hh" << 46 #include "globals.hh" 47 #else 47 #else 48 # define G4cout std::cout << 48 #define G4cout std::cout 49 # define G4cerr std::cerr << 49 #define G4cerr std::cerr 50 # define G4endl std::endl << 50 #define G4endl std::endl 51 # define G4String std::string // string incl << 51 #define G4String std::string // string included in PDBatom.hh 52 # include <cfloat> << 52 #include <cfloat> 53 #endif 53 #endif 54 #include <cmath> << 55 #include <fstream> 54 #include <fstream> 56 #include <iostream> 55 #include <iostream> 57 #include <limits> 56 #include <limits> >> 57 #include <cmath> 58 #include <sstream> 58 #include <sstream> 59 #include <stdlib.h> 59 #include <stdlib.h> 60 60 61 using namespace std; 61 using namespace std; 62 62 63 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 64 64 65 PDBlib::PDBlib() : fNbNucleotidsPerStrand(0) { << 65 PDBlib::PDBlib() : >> 66 fNbNucleotidsPerStrand(0) >> 67 { >> 68 } 66 69 67 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 68 71 69 /** 72 /** 70 * \brief Load a PDB file into memory 73 * \brief Load a PDB file into memory 71 * \details molecule (polymer,?), << 74 * \details molecule (polymer,?), 72 * read line, key words 75 * read line, key words 73 * \param filename G4String for filename 76 * \param filename G4String for filename 74 * \param isDNA 77 * \param isDNA 75 * \param verbose 78 * \param verbose 76 * \return List of molecules 79 * \return List of molecules 77 */ 80 */ 78 81 79 Molecule* PDBlib::Load(const std::string& file << 82 Molecule * PDBlib::Load(const std::string& filename, 80 unsigned short int verb << 83 unsigned short int& isDNA, >> 84 unsigned short int verbose = 0) 81 { 85 { 82 G4String sLine = ""; 86 G4String sLine = ""; 83 ifstream infile; 87 ifstream infile; 84 infile.open(filename.c_str()); 88 infile.open(filename.c_str()); 85 if (!infile) { << 89 if(!infile) 86 G4cout << "PDBlib::load >> file " << filen << 90 { >> 91 G4cout << "PDBlib::load >> file " << filename << " not found !!!!" >> 92 << G4endl; 87 } 93 } 88 else { << 94 else 89 int nbAtom = 0; // total number of atoms << 95 { 90 int numAtomInRes = 0; // number of an ato << 96 int nbAtomTot = 0; // total number of atoms 91 int nbResidue = 0; // total number of res << 97 int nbAtom = 0; // total number of atoms in a residue 92 int nbMolecule = 0; // total number of mo << 98 int numAtomInRes = 0; // number of an atom in a residue 93 << 99 int nbResidue = 0; // total number of residues 94 Atom* AtomicOld = nullptr; << 100 int nbMolecule = 0; // total number of molecule 95 Atom* AtomicNext = nullptr; << 101 96 << 102 Atom * AtomicOld = nullptr; 97 int serial; // Atom serial number << 103 Atom * AtomicNext = nullptr; 98 G4String atomName; // Atom name << 104 99 G4String element; // Element Symbol << 105 int serial; //Atom serial number 100 G4String resName; // Residue name for thi << 106 G4String atomName; //Atom name 101 double x, y, z; // Orthogonal coordinates << 107 G4String element; //Element Symbol 102 double occupancy; // Occupancy << 108 G4String resName; //Residue name for this atom 103 << 109 double x, y, z; //Orthogonal coordinates in Angstroms 104 Residue* residueOld = nullptr; << 110 double occupancy; //Occupancy 105 Residue* residueFirst = nullptr; << 111 106 Residue* residueNext = nullptr; << 112 Residue * residueOld = nullptr; 107 << 113 Residue * residueFirst = nullptr; 108 Molecule* moleculeFirst = nullptr; << 114 Residue * residueNext = nullptr; 109 Molecule* moleculeOld = nullptr; << 115 110 Molecule* moleculeNext = nullptr; << 116 Molecule * moleculeFirst = nullptr; >> 117 Molecule * moleculeOld = nullptr; >> 118 Molecule * moleculeNext = nullptr; 111 119 112 ///////////////////////////////// 120 ///////////////////////////////// 113 // NEW variable to draw a fitting cylinder << 121 //NEW variable to draw a fitting cylinder if z oriented 114 //=> fitting box 122 //=> fitting box 115 double minGlobZ, maxGlobZ; 123 double minGlobZ, maxGlobZ; 116 double minGlobX, maxGlobX; 124 double minGlobX, maxGlobX; 117 double minGlobY, maxGlobY; 125 double minGlobY, maxGlobY; 118 double minX, maxX, minY, maxY, minZ, maxZ; << 126 double minX, maxX, minY, maxY, minZ, maxZ; //Sort of 'mother volume' box 119 127 120 minGlobZ = -DBL_MAX; 128 minGlobZ = -DBL_MAX; 121 minGlobX = -DBL_MAX; 129 minGlobX = -DBL_MAX; 122 minGlobY = -DBL_MAX; 130 minGlobY = -DBL_MAX; 123 maxGlobZ = DBL_MAX; 131 maxGlobZ = DBL_MAX; 124 maxGlobX = DBL_MAX; 132 maxGlobX = DBL_MAX; 125 maxGlobY = DBL_MAX; 133 maxGlobY = DBL_MAX; 126 minX = -DBL_MAX; 134 minX = -DBL_MAX; 127 minY = -DBL_MAX; 135 minY = -DBL_MAX; 128 minZ = -DBL_MAX; 136 minZ = -DBL_MAX; 129 maxX = DBL_MAX; 137 maxX = DBL_MAX; 130 maxY = DBL_MAX; 138 maxY = DBL_MAX; 131 maxZ = DBL_MAX; 139 maxZ = DBL_MAX; 132 140 133 int lastResSeq = -1; 141 int lastResSeq = -1; 134 int resSeq = 0; 142 int resSeq = 0; 135 143 136 G4String nameMol; 144 G4String nameMol; 137 145 138 unsigned short int numModel = 0; 146 unsigned short int numModel = 0; 139 unsigned short int model = 0; 147 unsigned short int model = 0; 140 148 141 // Number of TER << 149 //Number of TER 142 int ter = 0; 150 int ter = 0; 143 // Number of TER (chain) max for this file << 151 //Number of TER (chain) max for this file 144 152 145 int terMax = INT_MAX; 153 int terMax = INT_MAX; 146 154 147 if (!infile.eof()) { << 155 if(!infile.eof()) >> 156 { 148 getline(infile, sLine); 157 getline(infile, sLine); 149 std::size_t found = sLine.find("DNA"); 158 std::size_t found = sLine.find("DNA"); 150 if (found != G4String::npos) { << 159 if(found != G4String::npos) >> 160 { 151 terMax = 2; 161 terMax = 2; 152 isDNA = 1; 162 isDNA = 1; 153 } 163 } 154 else << 164 else isDNA = 0; 155 isDNA = 0; << 165 //If PDB file have not a header line 156 // If PDB file have not a header line << 157 found = sLine.find("HEADER"); 166 found = sLine.find("HEADER"); 158 if (found == G4String::npos) { << 167 if(found == G4String::npos) >> 168 { 159 infile.close(); 169 infile.close(); 160 infile.open(filename.c_str()); 170 infile.open(filename.c_str()); 161 171 162 G4cout << "PDBlib::load >> No header f 172 G4cout << "PDBlib::load >> No header found !!!!" << G4endl; 163 } 173 } 164 } 174 } 165 175 166 while (!infile.eof()) { << 176 while(!infile.eof()) >> 177 { 167 getline(infile, sLine); 178 getline(infile, sLine); 168 // In the case of several models << 179 //In the case of several models 169 if ((sLine.substr(0, 6)).compare("NUMMDL << 180 if((sLine.substr(0, 6)).compare("NUMMDL") == 0) >> 181 { 170 istringstream((sLine.substr(10, 4))) > 182 istringstream((sLine.substr(10, 4))) >> numModel; 171 } 183 } 172 if ((numModel > 0) && ((sLine.substr(0, << 184 if((numModel > 0) && ((sLine.substr(0, 6)).compare("MODEL ") == 0)) >> 185 { 173 istringstream((sLine.substr(10, 4))) > 186 istringstream((sLine.substr(10, 4))) >> model; 174 ////////////////////////////////////// 187 ////////////////////////////////////////// 175 if (model > 1) break; << 188 if(model > 1) break; 176 ////////////////////////////////////// 189 ////////////////////////////////////////// 177 } 190 } 178 // For verification of residue sequence << 191 //For verification of residue sequence 179 if ((sLine.substr(0, 6)).compare("SEQRES << 192 if((sLine.substr(0, 6)).compare("SEQRES") == 0) 180 // Create list of molecule here << 193 { 181 if (verbose > 1) G4cout << sLine << G4 << 194 //Create list of molecule here >> 195 if(verbose > 1) G4cout << sLine << G4endl; 182 } 196 } 183 197 184 // Coordinate section << 198 //Coordinate section 185 if ((numModel > 0) && ((sLine.substr(0, << 199 if((numModel > 0) && ((sLine.substr(0, 6)).compare("ENDMDL") == 0)) >> 200 { 186 ; 201 ; 187 } 202 } 188 else if ((sLine.substr(0, 6)).compare("T << 203 else if((sLine.substr(0, 6)).compare("TER ") == 0) //3 spaces 189 { 204 { 190 ////////////////////////////////////// 205 ////////////////////////////////////////// 191 // Currently retrieve only the two fir << 206 //Currently retrieve only the two first chains(TER) => two DNA strands 192 ter++; 207 ter++; 193 if (ter > terMax) break; << 208 if(ter > terMax) break; 194 ////////////////////////////////////// 209 ////////////////////////////////////////// 195 210 196 // if (verbose>1) G4cout << sLine << G << 211 //if (verbose>1) G4cout << sLine << G4endl; 197 /************ Begin TER ************** 212 /************ Begin TER ******************/ 198 lastResSeq = -1; 213 lastResSeq = -1; 199 resSeq = 0; 214 resSeq = 0; 200 215 201 AtomicOld->SetNext(NULL); 216 AtomicOld->SetNext(NULL); 202 residueOld->SetNext(NULL); 217 residueOld->SetNext(NULL); 203 residueOld->fNbAtom = nbAtom; 218 residueOld->fNbAtom = nbAtom; 204 219 205 // Molecule creation: << 220 //Molecule creation: 206 if (moleculeOld == NULL) { << 221 if(moleculeOld == NULL) 207 nameMol = filename; //+numModel << 222 { >> 223 nameMol = filename; //+numModel 208 moleculeOld = new Molecule(nameMol, 224 moleculeOld = new Molecule(nameMol, nbMolecule); 209 moleculeOld->SetFirst(residueFirst); 225 moleculeOld->SetFirst(residueFirst); 210 moleculeOld->fNbResidue = nbResidue; 226 moleculeOld->fNbResidue = nbResidue; 211 moleculeFirst = moleculeOld; 227 moleculeFirst = moleculeOld; 212 } 228 } 213 else { << 229 else >> 230 { 214 moleculeNext = new Molecule(nameMol, 231 moleculeNext = new Molecule(nameMol, nbMolecule); 215 moleculeOld->SetNext(moleculeNext); 232 moleculeOld->SetNext(moleculeNext); 216 moleculeNext->SetFirst(residueFirst) 233 moleculeNext->SetFirst(residueFirst); 217 moleculeNext->fNbResidue = nbResidue 234 moleculeNext->fNbResidue = nbResidue; 218 moleculeOld = moleculeNext; 235 moleculeOld = moleculeNext; 219 } 236 } 220 237 221 nbMolecule++; 238 nbMolecule++; 222 moleculeOld->SetNext(NULL); 239 moleculeOld->SetNext(NULL); 223 moleculeOld->fCenterX = (int)((minX + << 240 moleculeOld->fCenterX = (int) ((minX + maxX) / 2.); 224 moleculeOld->fCenterY = (int)((minY + << 241 moleculeOld->fCenterY = (int) ((minY + maxY) / 2.); 225 moleculeOld->fCenterZ = (int)((minZ + << 242 moleculeOld->fCenterZ = (int) ((minZ + maxZ) / 2.); 226 moleculeOld->fMaxGlobZ = maxGlobZ; 243 moleculeOld->fMaxGlobZ = maxGlobZ; 227 moleculeOld->fMinGlobZ = minGlobZ; 244 moleculeOld->fMinGlobZ = minGlobZ; 228 moleculeOld->fMaxGlobX = maxGlobX; 245 moleculeOld->fMaxGlobX = maxGlobX; 229 moleculeOld->fMinGlobX = minGlobX; 246 moleculeOld->fMinGlobX = minGlobX; 230 moleculeOld->fMaxGlobY = maxGlobY; 247 moleculeOld->fMaxGlobY = maxGlobY; 231 moleculeOld->fMinGlobY = minGlobY; 248 moleculeOld->fMinGlobY = minGlobY; 232 249 233 minGlobZ = -DBL_MAX; 250 minGlobZ = -DBL_MAX; 234 minGlobX = -DBL_MAX; 251 minGlobX = -DBL_MAX; 235 minGlobY = -DBL_MAX; 252 minGlobY = -DBL_MAX; 236 maxGlobZ = DBL_MAX; 253 maxGlobZ = DBL_MAX; 237 maxGlobX = DBL_MAX; 254 maxGlobX = DBL_MAX; 238 maxGlobY = DBL_MAX; 255 maxGlobY = DBL_MAX; 239 minX = -DBL_MAX; 256 minX = -DBL_MAX; 240 minY = -DBL_MAX; 257 minY = -DBL_MAX; 241 minZ = -DBL_MAX; 258 minZ = -DBL_MAX; 242 maxX = DBL_MAX; 259 maxX = DBL_MAX; 243 maxY = DBL_MAX; 260 maxY = DBL_MAX; 244 maxZ = DBL_MAX; 261 maxZ = DBL_MAX; 245 262 246 nbAtom = 0; 263 nbAtom = 0; 247 numAtomInRes = 0; 264 numAtomInRes = 0; 248 nbResidue = 0; 265 nbResidue = 0; 249 AtomicOld = NULL; 266 AtomicOld = NULL; 250 AtomicNext = NULL; 267 AtomicNext = NULL; 251 residueOld = NULL; 268 residueOld = NULL; 252 residueFirst = NULL; 269 residueFirst = NULL; 253 residueNext = NULL; 270 residueNext = NULL; 254 271 255 ///////////// End TER //////////////// 272 ///////////// End TER /////////////////// 256 } 273 } 257 else if ((sLine.substr(0, 6)).compare("A << 274 else if((sLine.substr(0, 6)).compare("ATOM ") == 0) >> 275 { >> 276 258 /************ Begin ATOM ************* 277 /************ Begin ATOM ******************/ 259 // serial << 278 //serial 260 istringstream((sLine.substr(6, 5))) >> 279 istringstream((sLine.substr(6, 5))) >> serial; 261 280 262 // To be improved << 281 //To be improved 263 // atomName : << 282 //atomName : 264 atomName = sLine.substr(12, 4); 283 atomName = sLine.substr(12, 4); 265 if (atomName.substr(0, 1).compare(" ") << 284 if(atomName.substr(0, 1).compare(" ") == 0) element = sLine.substr(13, 266 element = sLine.substr(13, 1); << 285 1); 267 else << 286 else element = sLine.substr(12, 1); 268 element = sLine.substr(12, 1); << 269 287 270 // set Van der Waals radius expressed 288 // set Van der Waals radius expressed in Angstrom 271 double vdwRadius = -1.; 289 double vdwRadius = -1.; 272 if (element == "H") { << 290 if(element == "H") >> 291 { 273 vdwRadius = 1.2; 292 vdwRadius = 1.2; 274 } 293 } 275 else if (element == "C") { << 294 else if(element == "C") >> 295 { 276 vdwRadius = 1.7; 296 vdwRadius = 1.7; 277 } 297 } 278 else if (element == "N") { << 298 else if(element == "N") >> 299 { 279 vdwRadius = 1.55; 300 vdwRadius = 1.55; 280 } 301 } 281 else if (element == "O") { << 302 else if(element == "O") >> 303 { 282 vdwRadius = 1.52; 304 vdwRadius = 1.52; 283 } 305 } 284 else if (element == "P") { << 306 else if(element == "P") >> 307 { 285 vdwRadius = 1.8; 308 vdwRadius = 1.8; 286 } 309 } 287 else if (element == "S") { << 310 else if(element == "S") >> 311 { 288 vdwRadius = 1.8; 312 vdwRadius = 1.8; 289 } 313 } 290 else { << 314 else 291 #ifndef GEANT4 << 315 { >> 316 #ifndef GEANT4 292 G4cerr << "Element not recognized : 317 G4cerr << "Element not recognized : " << element << G4endl; 293 G4cerr << "Stop now" << G4endl; 318 G4cerr << "Stop now" << G4endl; 294 exit(1); 319 exit(1); 295 #else 320 #else 296 G4ExceptionDescription errMsg; 321 G4ExceptionDescription errMsg; 297 errMsg << "Element not recognized : 322 errMsg << "Element not recognized : " << element << G4endl; 298 323 299 G4Exception("PDBlib::Load", "ELEM_NO << 324 G4Exception("PDBlib::Load", 300 #endif << 325 "ELEM_NOT_RECOGNIZED", >> 326 FatalException, >> 327 errMsg); >> 328 #endif 301 } 329 } 302 330 303 { 331 { 304 // resName : << 332 nbAtomTot++; >> 333 >> 334 //resName : 305 resName = sLine.substr(17, 3); 335 resName = sLine.substr(17, 3); 306 // resSeq : << 336 //resSeq : 307 istringstream((sLine.substr(22, 4))) 337 istringstream((sLine.substr(22, 4))) >> resSeq; 308 // x,y,z : << 338 //x,y,z : 309 stringstream((sLine.substr(30, 8))) 339 stringstream((sLine.substr(30, 8))) >> x; 310 stringstream((sLine.substr(38, 8))) 340 stringstream((sLine.substr(38, 8))) >> y; 311 stringstream((sLine.substr(46, 8))) 341 stringstream((sLine.substr(46, 8))) >> z; 312 // occupancy : << 342 //occupancy : 313 occupancy = 1.; 343 occupancy = 1.; 314 344 315 if (minGlobZ < z) minGlobZ = z; << 345 if(minGlobZ < z) minGlobZ = z; 316 if (maxGlobZ > z) maxGlobZ = z; << 346 if(maxGlobZ > z) maxGlobZ = z; 317 if (minGlobX < x) minGlobX = x; << 347 if(minGlobX < x) minGlobX = x; 318 if (maxGlobX > x) maxGlobX = x; << 348 if(maxGlobX > x) maxGlobX = x; 319 if (minGlobY < y) minGlobY = y; << 349 if(minGlobY < y) minGlobY = y; 320 if (maxGlobY > y) maxGlobY = y; << 350 if(maxGlobY > y) maxGlobY = y; 321 if (minX > x) minX = x; << 351 if(minX > x) minX = x; 322 if (maxX < x) maxX = x; << 352 if(maxX < x) maxX = x; 323 if (minY > y) minY = y; << 353 if(minY > y) minY = y; 324 if (maxY < y) maxY = y; << 354 if(maxY < y) maxY = y; 325 if (minZ > z) minZ = z; << 355 if(minZ > z) minZ = z; 326 if (maxZ < z) maxZ = z; << 356 if(maxZ < z) maxZ = z; 327 << 357 328 // treatment for Atoms: << 358 //treatment for Atoms: 329 if (AtomicOld == NULL) { << 359 if(AtomicOld == NULL) 330 AtomicOld = << 360 { 331 new Atom(serial, atomName, "", 0 << 361 AtomicOld = new Atom(serial, 332 AtomicOld->SetNext(NULL); // If o << 362 atomName, 333 } << 363 "", 334 else { << 364 0, 335 AtomicNext = << 365 0, 336 new Atom(serial, atomName, "", 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); 337 AtomicOld->SetNext(AtomicNext); 389 AtomicOld->SetNext(AtomicNext); 338 AtomicOld = AtomicNext; 390 AtomicOld = AtomicNext; 339 } 391 } 340 nbAtom++; 392 nbAtom++; 341 } // END if (element!="H") << 393 } //END if (element!="H") 342 /****************************Begin Res 394 /****************************Begin Residue************************/ 343 // treatment for residues: << 395 //treatment for residues: 344 if (residueOld == NULL) { << 396 if(residueOld == NULL) 345 if (verbose > 2) G4cout << "residueO << 397 { >> 398 if(verbose > 2) G4cout << "residueOld == NULL" << G4endl; 346 399 347 AtomicOld->fNumInRes = 0; 400 AtomicOld->fNumInRes = 0; 348 residueOld = new Residue(resName, re 401 residueOld = new Residue(resName, resSeq); 349 residueOld->SetFirst(AtomicOld); 402 residueOld->SetFirst(AtomicOld); 350 residueOld->SetNext(NULL); 403 residueOld->SetNext(NULL); 351 residueFirst = residueOld; 404 residueFirst = residueOld; 352 lastResSeq = resSeq; 405 lastResSeq = resSeq; 353 nbResidue++; 406 nbResidue++; 354 } 407 } 355 else { << 408 else 356 if (lastResSeq == resSeq) { << 409 { >> 410 if(lastResSeq == resSeq) >> 411 { 357 numAtomInRes++; 412 numAtomInRes++; 358 AtomicOld->fNumInRes = numAtomInRe 413 AtomicOld->fNumInRes = numAtomInRes; 359 } 414 } 360 else { << 415 else >> 416 { 361 numAtomInRes = 0; 417 numAtomInRes = 0; 362 AtomicOld->fNumInRes = numAtomInRe 418 AtomicOld->fNumInRes = numAtomInRes; 363 419 364 residueNext = new Residue(resName, 420 residueNext = new Residue(resName, resSeq); 365 residueNext->SetFirst(AtomicOld); 421 residueNext->SetFirst(AtomicOld); 366 residueOld->SetNext(residueNext); 422 residueOld->SetNext(residueNext); 367 residueOld->fNbAtom = nbAtom - 1; 423 residueOld->fNbAtom = nbAtom - 1; 368 424 369 nbAtom = 1; 425 nbAtom = 1; 370 residueOld = residueNext; 426 residueOld = residueNext; 371 lastResSeq = resSeq; 427 lastResSeq = resSeq; 372 nbResidue++; 428 nbResidue++; 373 } 429 } 374 } 430 } 375 /////////////////////////End Residue// 431 /////////////////////////End Residue//////////// 376 ///////////// End ATOM /////////////// 432 ///////////// End ATOM /////////////////// 377 } // END if Atom << 433 } //END if Atom 378 } // END while (!infile.eof()) << 434 } //END while (!infile.eof()) 379 435 380 infile.close(); 436 infile.close(); 381 return moleculeFirst; 437 return moleculeFirst; 382 } // END else if (!infile) << 438 } //END else if (!infile) 383 return NULL; 439 return NULL; 384 } 440 } 385 441 386 //....oooOO0OOooo........oooOO0OOooo........oo 442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 387 443 388 /** 444 /** 389 * \brief Compute barycenters 445 * \brief Compute barycenters 390 * \details Compute barycenters and its coord << 446 * \details Compute barycenters and its coordinate 391 * for nucleotides 447 * for nucleotides 392 * \param moleculeListTemp * moleculeLis 448 * \param moleculeListTemp * moleculeList 393 * \return Barycenter * 449 * \return Barycenter * 394 */ 450 */ 395 Barycenter* PDBlib::ComputeNucleotideBarycente << 451 Barycenter * PDBlib::ComputeNucleotideBarycenters(Molecule * moleculeListTemp) 396 { 452 { 397 //////////////////////////////////////////// 453 /////////////////////////////////////////////////////// 398 // Placement and physical volume constructio << 454 //Placement and physical volume construction from memory 399 Barycenter* BarycenterFirst = NULL; << 455 Barycenter * BarycenterFirst = NULL; 400 Barycenter* BarycenterOld = NULL; << 456 Barycenter * BarycenterOld = NULL; 401 Barycenter* BarycenterNext = NULL; << 457 Barycenter * BarycenterNext = NULL; 402 << 458 403 // Residue (Base, Phosphate,sugar) list << 459 //Residue (Base, Phosphate,sugar) list 404 Residue* residueListTemp; << 460 Residue *residueListTemp; 405 // Atom list inside a residu << 461 //Atom list inside a residu 406 Atom* AtomTemp; << 462 Atom *AtomTemp; 407 463 408 int k = 0; 464 int k = 0; 409 int j_old = 0; 465 int j_old = 0; 410 466 411 while (moleculeListTemp) { << 467 while(moleculeListTemp) >> 468 { 412 residueListTemp = moleculeListTemp->GetFir 469 residueListTemp = moleculeListTemp->GetFirst(); 413 470 414 k++; 471 k++; 415 int j = 0; 472 int j = 0; 416 473 417 // Check numerotation style (1->n per stra << 474 //Check numerotation style (1->n per strand or 1->2n for two strand) 418 int correctNumerotationNumber = 0; 475 int correctNumerotationNumber = 0; 419 if (k == 2 && residueListTemp->fResSeq > 1 << 476 if(k == 2 && residueListTemp->fResSeq > 1) >> 477 { 420 correctNumerotationNumber = residueListT 478 correctNumerotationNumber = residueListTemp->fResSeq; 421 } 479 } 422 480 423 while (residueListTemp) { << 481 while(residueListTemp) >> 482 { 424 AtomTemp = residueListTemp->GetFirst(); 483 AtomTemp = residueListTemp->GetFirst(); 425 j++; 484 j++; 426 485 427 // Correction consequently to numerotati << 486 //Correction consequently to numerotation check 428 if (correctNumerotationNumber != 0) { << 487 if(correctNumerotationNumber != 0) 429 residueListTemp->fResSeq = residueList << 488 { >> 489 residueListTemp->fResSeq = residueListTemp->fResSeq >> 490 - correctNumerotationNumber >> 491 + 1; 430 } 492 } 431 493 432 // Barycenter computation << 494 //Barycenter computation 433 double baryX = 0., baryY = 0., baryZ = 0 495 double baryX = 0., baryY = 0., baryZ = 0.; 434 double baryBaseX = 0., baryBaseY = 0., b 496 double baryBaseX = 0., baryBaseY = 0., baryBaseZ = 0.; 435 double barySugX = 0., barySugY = 0., bar 497 double barySugX = 0., barySugY = 0., barySugZ = 0.; 436 double baryPhosX = 0., baryPhosY = 0., b 498 double baryPhosX = 0., baryPhosY = 0., baryPhosZ = 0.; 437 unsigned short int nbAtomInBase = 0; 499 unsigned short int nbAtomInBase = 0; 438 500 439 for (int i = 0; i < residueListTemp->fNb << 501 for(int i = 0; i < residueListTemp->fNbAtom; i++) 440 // Compute barycenter of the nucletoti << 502 { >> 503 //Compute barycenter of the nucletotide 441 baryX += AtomTemp->fX; 504 baryX += AtomTemp->fX; 442 baryY += AtomTemp->fY; 505 baryY += AtomTemp->fY; 443 baryZ += AtomTemp->fZ; 506 baryZ += AtomTemp->fZ; 444 // Compute barycenters for Base Sugar << 507 //Compute barycenters for Base Sugar Phosphat 445 if (residueListTemp->fResSeq == 1) { << 508 if(residueListTemp->fResSeq == 1) 446 if (i == 0) { << 509 { >> 510 if(i == 0) >> 511 { 447 baryPhosX += AtomTemp->fX; 512 baryPhosX += AtomTemp->fX; 448 baryPhosY += AtomTemp->fY; 513 baryPhosY += AtomTemp->fY; 449 baryPhosZ += AtomTemp->fZ; 514 baryPhosZ += AtomTemp->fZ; 450 } 515 } 451 else if (i < 8) { << 516 else if(i < 8) >> 517 { 452 barySugX += AtomTemp->fX; 518 barySugX += AtomTemp->fX; 453 barySugY += AtomTemp->fY; 519 barySugY += AtomTemp->fY; 454 barySugZ += AtomTemp->fZ; 520 barySugZ += AtomTemp->fZ; 455 } 521 } 456 else { << 522 else 457 // hydrogen are placed at the end << 523 { 458 // We don't want them for this cal << 524 //hydrogen are placed at the end of the residue in a PDB file 459 if (AtomTemp->fElement != "H") { << 525 //We don't want them for this calculation >> 526 if(AtomTemp->fElement != "H") >> 527 { 460 baryBaseX += AtomTemp->fX; 528 baryBaseX += AtomTemp->fX; 461 baryBaseY += AtomTemp->fY; 529 baryBaseY += AtomTemp->fY; 462 baryBaseZ += AtomTemp->fZ; 530 baryBaseZ += AtomTemp->fZ; 463 nbAtomInBase++; 531 nbAtomInBase++; 464 } 532 } 465 } 533 } 466 } 534 } 467 else { << 535 else 468 if (i < 4) { << 536 { >> 537 if(i < 4) >> 538 { 469 baryPhosX += AtomTemp->fX; 539 baryPhosX += AtomTemp->fX; 470 baryPhosY += AtomTemp->fY; 540 baryPhosY += AtomTemp->fY; 471 baryPhosZ += AtomTemp->fZ; 541 baryPhosZ += AtomTemp->fZ; 472 } 542 } 473 else if (i < 11) { << 543 else if(i < 11) >> 544 { 474 barySugX += AtomTemp->fX; 545 barySugX += AtomTemp->fX; 475 barySugY += AtomTemp->fY; 546 barySugY += AtomTemp->fY; 476 barySugZ += AtomTemp->fZ; 547 barySugZ += AtomTemp->fZ; 477 } 548 } 478 else { << 549 else 479 // hydrogen are placed at the end << 550 { 480 // We don't want them for this cal << 551 //hydrogen are placed at the end of the residue in a PDB file 481 if (AtomTemp->fElement != "H") { << 552 //We don't want them for this calculation >> 553 if(AtomTemp->fElement != "H") >> 554 { // break; 482 baryBaseX += AtomTemp->fX; 555 baryBaseX += AtomTemp->fX; 483 baryBaseY += AtomTemp->fY; 556 baryBaseY += AtomTemp->fY; 484 baryBaseZ += AtomTemp->fZ; 557 baryBaseZ += AtomTemp->fZ; 485 nbAtomInBase++; 558 nbAtomInBase++; 486 } 559 } 487 } 560 } 488 } 561 } 489 AtomTemp = AtomTemp->GetNext(); 562 AtomTemp = AtomTemp->GetNext(); 490 } // end of for ( i=0 ; i < residueLis << 563 } //end of for ( i=0 ; i < residueListTemp->nbAtom ; i++) 491 564 492 baryX = baryX / (double)residueListTemp- << 565 baryX = baryX / (double) residueListTemp->fNbAtom; 493 baryY = baryY / (double)residueListTemp- << 566 baryY = baryY / (double) residueListTemp->fNbAtom; 494 baryZ = baryZ / (double)residueListTemp- << 567 baryZ = baryZ / (double) residueListTemp->fNbAtom; 495 568 496 if (residueListTemp->fResSeq != 1) // S << 569 if(residueListTemp->fResSeq != 1) //Special case first Phosphate 497 { 570 { 498 baryPhosX = baryPhosX / 4.; 571 baryPhosX = baryPhosX / 4.; 499 baryPhosY = baryPhosY / 4.; 572 baryPhosY = baryPhosY / 4.; 500 baryPhosZ = baryPhosZ / 4.; 573 baryPhosZ = baryPhosZ / 4.; 501 } 574 } 502 barySugX = barySugX / 7.; 575 barySugX = barySugX / 7.; 503 barySugY = barySugY / 7.; 576 barySugY = barySugY / 7.; 504 barySugZ = barySugZ / 7.; 577 barySugZ = barySugZ / 7.; 505 baryBaseX = baryBaseX / (double)nbAtomIn << 578 baryBaseX = baryBaseX / (double) nbAtomInBase; 506 baryBaseY = baryBaseY / (double)nbAtomIn << 579 baryBaseY = baryBaseY / (double) nbAtomInBase; 507 baryBaseZ = baryBaseZ / (double)nbAtomIn << 580 baryBaseZ = baryBaseZ / (double) nbAtomInBase; 508 << 581 509 // Barycenter creation: << 582 //Barycenter creation: 510 if (BarycenterOld == NULL) { << 583 if(BarycenterOld == NULL) 511 BarycenterOld = new Barycenter(j + j_o << 584 { 512 baryBas << 585 BarycenterOld = new Barycenter(j + j_old, baryX, baryY, baryZ, //j [1..n] 513 barySug << 586 baryBaseX, >> 587 baryBaseY, >> 588 baryBaseZ, >> 589 barySugX, >> 590 barySugY, >> 591 barySugZ, >> 592 baryPhosX, >> 593 baryPhosY, >> 594 baryPhosZ); 514 BarycenterFirst = BarycenterOld; 595 BarycenterFirst = BarycenterOld; 515 } 596 } 516 else { << 597 else 517 BarycenterNext = << 598 { 518 new Barycenter(j + j_old, baryX, bar << 599 BarycenterNext = new Barycenter(j + j_old, 519 barySugY, barySugZ, b << 600 baryX, >> 601 baryY, >> 602 baryZ, >> 603 baryBaseX, >> 604 baryBaseY, >> 605 baryBaseZ, >> 606 barySugX, >> 607 barySugY, >> 608 barySugZ, >> 609 baryPhosX, >> 610 baryPhosY, >> 611 baryPhosZ); 520 BarycenterOld->SetNext(BarycenterNext) 612 BarycenterOld->SetNext(BarycenterNext); 521 BarycenterOld = BarycenterNext; 613 BarycenterOld = BarycenterNext; 522 } 614 } 523 615 524 //////////////////////////////////////// 616 ///////////////////////////////////////////////// 525 // distance computation between all atom << 617 //distance computation between all atoms inside 526 // a residue and the barycenter << 618 //a residue and the barycenter 527 AtomTemp = residueListTemp->GetFirst(); 619 AtomTemp = residueListTemp->GetFirst(); 528 double dT3Dp; 620 double dT3Dp; 529 double max = 0.; 621 double max = 0.; 530 for (int ii = 0; ii < residueListTemp->f << 622 for(int ii = 0; ii < residueListTemp->fNbAtom; ii++) 531 dT3Dp = DistanceTwo3Dpoints(AtomTemp-> << 623 { 532 Barycenter << 624 dT3Dp = DistanceTwo3Dpoints(AtomTemp->fX, >> 625 BarycenterOld->fCenterX, >> 626 AtomTemp->fY, >> 627 BarycenterOld->fCenterY, >> 628 AtomTemp->fZ, >> 629 BarycenterOld->fCenterZ); 533 BarycenterOld->SetDistance(ii, dT3Dp); 630 BarycenterOld->SetDistance(ii, dT3Dp); 534 if (dT3Dp > max) max = dT3Dp; << 631 if(dT3Dp > max) max = dT3Dp; 535 AtomTemp = AtomTemp->GetNext(); 632 AtomTemp = AtomTemp->GetNext(); 536 } // end of for ( i=0 ; i < residueLis << 633 } //end of for ( i=0 ; i < residueListTemp->nbAtom ; i++) 537 634 538 BarycenterOld->SetRadius(max + 1.8); 635 BarycenterOld->SetRadius(max + 1.8); 539 residueListTemp = residueListTemp->GetNe 636 residueListTemp = residueListTemp->GetNext(); 540 637 541 } // end of while sur residueListTemp << 638 } //end of while sur residueListTemp 542 639 543 j_old += j; 640 j_old += j; 544 641 545 /// molecs->push_back(*moleculeListTemp); << 642 ///molecs->push_back(*moleculeListTemp); 546 moleculeListTemp = moleculeListTemp->GetNe 643 moleculeListTemp = moleculeListTemp->GetNext(); 547 } // end of while sur moleculeListTemp << 644 } //end of while sur moleculeListTemp 548 645 549 if (BarycenterNext != NULL) { << 646 if(BarycenterNext != NULL) >> 647 { 550 BarycenterNext->SetNext(NULL); 648 BarycenterNext->SetNext(NULL); 551 } 649 } 552 650 553 return BarycenterFirst; 651 return BarycenterFirst; 554 } 652 } 555 653 556 //....oooOO0OOooo........oooOO0OOooo........oo 654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 557 655 558 /** 656 /** 559 * \brief the corresponding bounding volume 657 * \brief the corresponding bounding volume parameters 560 * \details the corresponding bounding volume 658 * \details the corresponding bounding volume parameters 561 * to build a box from atoms coordi 659 * to build a box from atoms coordinates 562 */ 660 */ 563 void PDBlib::ComputeBoundingVolumeParams(Molec << 661 void PDBlib::ComputeBoundingVolumeParams(Molecule *moleculeListTemp, 564 doubl << 662 double &dX, 565 doubl << 663 double &dY, 566 doubl << 664 double &dZ, //Dimensions for bounding volume >> 665 double &tX, >> 666 double &tY, >> 667 double &tZ) //Translation for bounding volume 567 { 668 { 568 double minminX, minminY, minminZ; // minimu << 669 double minminX, minminY, minminZ; //minimum minimorum 569 double maxmaxX, maxmaxY, maxmaxZ; // maximu << 670 double maxmaxX, maxmaxY, maxmaxZ; //maximum maximorum 570 671 571 minminX = DBL_MAX; 672 minminX = DBL_MAX; 572 minminY = DBL_MAX; 673 minminY = DBL_MAX; 573 minminZ = DBL_MAX; 674 minminZ = DBL_MAX; 574 maxmaxX = -DBL_MAX; 675 maxmaxX = -DBL_MAX; 575 maxmaxY = -DBL_MAX; 676 maxmaxY = -DBL_MAX; 576 maxmaxZ = -DBL_MAX; 677 maxmaxZ = -DBL_MAX; 577 678 578 while (moleculeListTemp) { << 679 while(moleculeListTemp) 579 if (minminX > moleculeListTemp->fMaxGlobX) << 680 { >> 681 if(minminX > moleculeListTemp->fMaxGlobX) >> 682 { 580 minminX = moleculeListTemp->fMaxGlobX; 683 minminX = moleculeListTemp->fMaxGlobX; 581 } 684 } 582 if (minminY > moleculeListTemp->fMaxGlobY) << 685 if(minminY > moleculeListTemp->fMaxGlobY) >> 686 { 583 minminY = moleculeListTemp->fMaxGlobY; 687 minminY = moleculeListTemp->fMaxGlobY; 584 } 688 } 585 if (minminZ > moleculeListTemp->fMaxGlobZ) << 689 if(minminZ > moleculeListTemp->fMaxGlobZ) >> 690 { 586 minminZ = moleculeListTemp->fMaxGlobZ; 691 minminZ = moleculeListTemp->fMaxGlobZ; 587 } 692 } 588 if (maxmaxX < moleculeListTemp->fMinGlobX) << 693 if(maxmaxX < moleculeListTemp->fMinGlobX) >> 694 { 589 maxmaxX = moleculeListTemp->fMinGlobX; 695 maxmaxX = moleculeListTemp->fMinGlobX; 590 } 696 } 591 if (maxmaxY < moleculeListTemp->fMinGlobY) << 697 if(maxmaxY < moleculeListTemp->fMinGlobY) >> 698 { 592 maxmaxY = moleculeListTemp->fMinGlobY; 699 maxmaxY = moleculeListTemp->fMinGlobY; 593 } 700 } 594 if (maxmaxZ < moleculeListTemp->fMinGlobZ) << 701 if(maxmaxZ < moleculeListTemp->fMinGlobZ) >> 702 { 595 maxmaxZ = moleculeListTemp->fMinGlobZ; 703 maxmaxZ = moleculeListTemp->fMinGlobZ; 596 } 704 } 597 705 598 moleculeListTemp = moleculeListTemp->GetNe 706 moleculeListTemp = moleculeListTemp->GetNext(); 599 } // end of while sur moleculeListTemp << 707 } //end of while sur moleculeListTemp 600 708 601 dX = (maxmaxX - minminX) / 2. + 1.8; // 1.8 << 709 dX = (maxmaxX - minminX) / 2. + 1.8; //1.8 => size of biggest radius for atoms 602 dY = (maxmaxY - minminY) / 2. + 1.8; 710 dY = (maxmaxY - minminY) / 2. + 1.8; 603 dZ = (maxmaxZ - minminZ) / 2. + 1.8; 711 dZ = (maxmaxZ - minminZ) / 2. + 1.8; 604 712 605 tX = minminX + (maxmaxX - minminX) / 2.; 713 tX = minminX + (maxmaxX - minminX) / 2.; 606 tY = minminY + (maxmaxY - minminY) / 2.; 714 tY = minminY + (maxmaxY - minminY) / 2.; 607 tZ = minminZ + (maxmaxZ - minminZ) / 2.; 715 tZ = minminZ + (maxmaxZ - minminZ) / 2.; 608 } 716 } 609 717 610 //....oooOO0OOooo........oooOO0OOooo........oo 718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 611 719 612 /** 720 /** 613 * \brief Compute number of nucleotide per 721 * \brief Compute number of nucleotide per strand 614 * \details Compute number of nucleotide per << 722 * \details Compute number of nucleotide per strand 615 * for DNA 723 * for DNA 616 */ 724 */ 617 void PDBlib::ComputeNbNucleotidsPerStrand(Mole << 725 void PDBlib::ComputeNbNucleotidsPerStrand(Molecule * moleculeListTemp) 618 { 726 { 619 Residue* residueListTemp; << 727 Residue *residueListTemp; 620 728 >> 729 int k = 0; 621 int j_old = 0; 730 int j_old = 0; 622 731 623 while (moleculeListTemp) { << 732 while(moleculeListTemp) >> 733 { 624 residueListTemp = moleculeListTemp->GetFir 734 residueListTemp = moleculeListTemp->GetFirst(); 625 735 >> 736 k++; 626 int j = 0; 737 int j = 0; 627 738 628 while (residueListTemp) { << 739 while(residueListTemp) >> 740 { 629 j++; 741 j++; 630 residueListTemp = residueListTemp->GetNe 742 residueListTemp = residueListTemp->GetNext(); 631 } // end of while sur residueListTemp << 743 } //end of while sur residueListTemp 632 744 633 j_old += j; 745 j_old += j; 634 moleculeListTemp = moleculeListTemp->GetNe 746 moleculeListTemp = moleculeListTemp->GetNext(); 635 } // end of while sur moleculeListTemp << 747 } //end of while sur moleculeListTemp 636 748 637 fNbNucleotidsPerStrand = j_old / 2; 749 fNbNucleotidsPerStrand = j_old / 2; 638 } 750 } 639 751 640 //....oooOO0OOooo........oooOO0OOooo........oo 752 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 641 753 642 /** 754 /** 643 * \brief Compute barycenters 755 * \brief Compute barycenters 644 * \details Compute barycenters and its coord << 756 * \details Compute barycenters and its coordinate 645 * for nucleotides 757 * for nucleotides 646 */ 758 */ 647 unsigned short int PDBlib::ComputeMatchEdepDNA << 759 unsigned short int PDBlib::ComputeMatchEdepDNA(Barycenter *BarycenterList, 648 << 760 Molecule *moleculeListTemp, 649 << 761 double x, 650 << 762 double y, >> 763 double z, >> 764 int &numStrand, >> 765 int &numNucleotid, >> 766 int &codeResidue) 651 { 767 { 652 unsigned short int matchFound = 0; 768 unsigned short int matchFound = 0; 653 769 654 Molecule* mLTsavedPointer = moleculeListTemp << 770 Molecule *mLTsavedPointer = moleculeListTemp; 655 Barycenter* BLsavedPointer = BarycenterList; << 771 Barycenter *BLsavedPointer = BarycenterList; 656 772 657 short int strandNum = 0; // Strand number << 773 short int strandNum = 0; //Strand number 658 int residueNum = 1; // Residue (nucleotide) << 774 int residueNum = 1; //Residue (nucleotide) number 659 G4String baseName; // Base name [A,C,T,G] << 775 G4String baseName; //Base name [A,C,T,G] 660 unsigned short int BSP = 2; // Base (defaul << 776 unsigned short int BSP = 2; //Base (default value), Sugar, Phosphat 661 777 662 double smallestDist; // smallest dist Atom << 778 double smallestDist; //smallest dist Atom <-> edep coordinates 663 double distEdepDNA; 779 double distEdepDNA; 664 double distEdepAtom; 780 double distEdepAtom; 665 781 666 // Residue (Base, Phosphate,suggar) list << 782 //Residue (Base, Phosphate,suggar) list 667 Residue* residueListTemp; << 783 Residue *residueListTemp; 668 // Atom list inside a residue << 784 //Atom list inside a residue 669 Atom* AtomTemp; << 785 Atom *AtomTemp; 670 786 671 int k = 0; // Molecule number << 787 int k = 0; //Molecule number 672 moleculeListTemp = mLTsavedPointer; 788 moleculeListTemp = mLTsavedPointer; 673 BarycenterList = BLsavedPointer; 789 BarycenterList = BLsavedPointer; 674 790 675 smallestDist = 33.0; // Sufficiently large << 791 smallestDist = 33.0; //Sufficiently large value 676 while (moleculeListTemp) { << 792 while(moleculeListTemp) >> 793 { 677 k++; 794 k++; 678 residueListTemp = moleculeListTemp->GetFir 795 residueListTemp = moleculeListTemp->GetFirst(); 679 796 680 int j = 0; // Residue number << 797 int j = 0; //Residue number 681 798 682 int j_save = INT_MAX; // Saved res. numbe << 799 int j_save = INT_MAX; //Saved res. number if match found 683 800 684 while (residueListTemp) { << 801 while(residueListTemp) >> 802 { 685 j++; 803 j++; 686 804 687 if (j - j_save > 2) break; << 805 if(j - j_save > 2) break; 688 806 689 distEdepDNA = DistanceTwo3Dpoints(x, Bar << 807 distEdepDNA = DistanceTwo3Dpoints(x, >> 808 BarycenterList->fCenterX, >> 809 y, >> 810 BarycenterList->fCenterY, >> 811 z, 690 Baryce 812 BarycenterList->fCenterZ); 691 if (distEdepDNA < BarycenterList->GetRad << 813 if(distEdepDNA < BarycenterList->GetRadius()) 692 // Find the closest atom << 814 { 693 // Compute distance between energy dep << 815 //Find the closest atom 694 // if distance <1.8 then match OK but << 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 695 AtomTemp = residueListTemp->GetFirst() 818 AtomTemp = residueListTemp->GetFirst(); 696 for (int iii = 0; iii < residueListTem << 819 for(int iii = 0; iii < residueListTemp->fNbAtom; iii++) 697 distEdepAtom = << 820 { 698 DistanceTwo3Dpoints(x, AtomTemp->G << 699 821 700 if ((distEdepAtom < AtomTemp->GetVan << 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 { 701 strandNum = k; 832 strandNum = k; 702 833 703 if (k == 2) { << 834 if(k == 2) 704 residueNum = fNbNucleotidsPerStr << 835 { >> 836 residueNum = fNbNucleotidsPerStrand + 1 >> 837 - residueListTemp->fResSeq; 705 } 838 } 706 else { << 839 else >> 840 { 707 residueNum = residueListTemp->fR 841 residueNum = residueListTemp->fResSeq; 708 } 842 } 709 843 710 baseName = (residueListTemp->fResN 844 baseName = (residueListTemp->fResName)[2]; 711 if (residueListTemp->fResSeq == 1) << 845 if(residueListTemp->fResSeq == 1) 712 if (iii == 0) << 846 { 713 BSP = 0; //"Phosphate" << 847 if(iii == 0) BSP = 0; //"Phosphate" 714 else if (iii < 8) << 848 else if(iii < 8) BSP = 1; //"Sugar" 715 BSP = 1; //"Sugar" << 849 else BSP = 2; //"Base" 716 else << 717 BSP = 2; //"Base" << 718 } 850 } 719 else { << 851 else 720 if (iii < 4) << 852 { 721 BSP = 0; //"Phosphate" << 853 if(iii < 4) BSP = 0; //"Phosphate" 722 else if (iii < 11) << 854 else if(iii < 11) BSP = 1; //"Sugar" 723 BSP = 1; //"Sugar" << 855 else BSP = 2; //"Base" 724 else << 725 BSP = 2; //"Base" << 726 } 856 } 727 857 728 smallestDist = distEdepAtom; 858 smallestDist = distEdepAtom; 729 859 730 int j_max_value = INT_MAX; 860 int j_max_value = INT_MAX; 731 861 732 if (j_save == j_max_value) j_save << 862 if(j_save == j_max_value) j_save = j; 733 matchFound = 1; 863 matchFound = 1; 734 } 864 } 735 AtomTemp = AtomTemp->GetNext(); 865 AtomTemp = AtomTemp->GetNext(); 736 } // end of for ( iii=0 ; iii < resi << 866 } //end of for ( iii=0 ; iii < residueListTemp->nbAtom ; iii++) 737 } // end for if (distEdepDNA < Barycent << 867 } //end for if (distEdepDNA < BarycenterList->GetRadius()) 738 BarycenterList = BarycenterList->GetNext 868 BarycenterList = BarycenterList->GetNext(); 739 residueListTemp = residueListTemp->GetNe 869 residueListTemp = residueListTemp->GetNext(); 740 } // end of while sur residueListTemp << 870 } //end of while sur residueListTemp 741 moleculeListTemp = moleculeListTemp->GetNe 871 moleculeListTemp = moleculeListTemp->GetNext(); 742 } // end of while sur moleculeListTemp << 872 } //end of while sur moleculeListTemp 743 873 744 numStrand = strandNum; 874 numStrand = strandNum; 745 numNucleotid = residueNum; 875 numNucleotid = residueNum; 746 codeResidue = BSP; 876 codeResidue = BSP; 747 877 748 return matchFound; 878 return matchFound; 749 } 879 } 750 880 751 //....oooOO0OOooo........oooOO0OOooo........oo 881 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 752 882 753 double PDBlib::DistanceTwo3Dpoints(double xA, << 883 double PDBlib::DistanceTwo3Dpoints(double xA, >> 884 double xB, >> 885 double yA, >> 886 double yB, >> 887 double zA, >> 888 double zB) 754 { 889 { 755 return sqrt((xA - xB) * (xA - xB) + (yA - yB << 890 return sqrt((xA - xB) * (xA - xB) + (yA - yB) * (yA - yB) >> 891 + (zA - zB) * (zA - zB)); 756 } 892 } 757 893