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