pdbstruct.cpp

Go to the documentation of this file.
00001 //
00002 // Utilitaires  pour manipuler les fichiers pdb
00003 //
00004 #include <iostream>
00005 #include <sstream>
00006 #include <iomanip>
00007 #include <algorithm>
00008 
00009 #include "pdbstruct.h"
00010 #include "tools.h"
00011 
00012 const std::string FRAG("FRAG") ;
00013 const std::string AA("AA") ;
00014 
00015 
00016 
00017 PDBelement::PDBelement()  : _numInMolecule(-1) ,_numInResidu(-1),_numInFragment(-1),_numResiduInMolecule(-1),
00018                             _connection(),_residuName("DefaultGroup"), _atomName(""),_fragName("DefaultGroup")
00019 {
00020   _xyz[0] = _xyz[1] = _xyz[2] = 0.0 ; 
00021 }
00022 PDBelement::~PDBelement() 
00023 {
00024 }
00025 void PDBelement::readAtom(std::stringstream & sline)
00026 {
00027   //
00028   // Attention il y a 1 blanc entre chaque element ;
00029   std::string atom ;
00030   sline >>  _numInMolecule ;
00031   sline >> _atomName ; 
00032   //
00033   // No blanc between atomName and AminoAcid Name
00034   //
00035   if(_atomName.size() <= 4){
00036     sline >> _residuName ;   
00037   } 
00038   else {
00039     std::string tmp ;
00040     tmp.assign(_atomName,0,4);
00041     _residuName.assign(_atomName,4,_atomName.size()-4);
00042     _atomName.assign(tmp) ;
00043   }
00044   //
00045   sline >>_numResiduInMolecule ;
00046   sline >>_xyz[0] >> _xyz[1] >>  _xyz[2] ;
00047   //
00048   // a cause du C (on commence a zero );
00049   _numResiduInMolecule--; 
00050 
00051   _numInMolecule-- ;
00052 }
00053 
00054 void PDBelement::readConnection(std::stringstream & sline)
00055 {
00056   //
00057   // It's Ok now for numbers of atoms bigger than 9999 
00058   int at , l ;
00059    std::string numAtoms ; 
00060    while(sline >> numAtoms ){
00061     l = numAtoms.size() ;
00062     while (l >5) {
00063       l -= 5 ;
00064       numAtoms.insert(l," ");
00065     }
00066     std::stringstream newsline ; newsline << numAtoms ;
00067     while (newsline >> at ){
00068       // a cause du C  at = at -1 (on commence a zero)
00069       _connection.push_back(--at); 
00070     }   
00071   }
00072 }
00073 void PDBelement::copy(const PDBelement & at) 
00074 {
00075   _xyz[0]   = at._xyz[0] ;
00076   _xyz[1]   = at._xyz[1] ;
00077   _xyz[2]   = at._xyz[2] ;
00078   //
00079   _numInMolecule       = at._numInMolecule    ;
00080   _numInResidu         = at._numInResidu    ;
00081   _numResiduInMolecule = at._numResiduInMolecule    ;
00082   _numInFragment       = at._numInFragment    ;
00083   //
00084   _atomName.assign(at._atomName) ;
00085   _residuName.assign(at._residuName) ;
00086   _fragName.assign(at._fragName) ;
00087   _connection.resize(at._connection.size()) ;
00088   for ( unsigned int i = 0 ; i < _connection.size(); ++i){
00089     _connection[i] = at._connection[i] ;
00090   }
00091 }
00092 
00093 
00094 PDBstructure::PDBstructure() :  _systemCharge(0),_spinMultiplicity(1),
00095                                  _numberOfAtoms(0), _numberOfResidus(0),
00096                                 _numberOfFragments(0),_atoms(NULL)
00097 {}
00098 PDBstructure::PDBstructure(const int size)  :  _systemCharge(0),_spinMultiplicity(1),
00099                                                _numberOfAtoms(size) , _numberOfResidus(0),
00100                                                _numberOfFragments(0),_atoms(new PDBelement[size])
00101 {}
00102 
00103 PDBstructure::~PDBstructure() 
00104 {
00105   if(_atoms != NULL)
00106     {delete [] _atoms ;_atoms= NULL ;}
00107   _numberOfAtoms = _numberOfResidus= _numberOfFragments =_systemCharge = 0 ; 
00108   _spinMultiplicity = 1;
00109   
00110 }
00111 bool PDBstructure::readPDBFile(std::ifstream & data) 
00112 {
00113   const std::string ATOM("ATOM") ;
00114   const std::string HETETOATOM("HETATM") ;
00115   const std::string CONECT("CONECT") ;
00116   const std::string BLANC(" ") ;
00117   std::string line, keyword ;
00118   this->_numberOfAtoms = 0 ;
00119   while ( !data.eof() )
00120     { 
00121       std::stringstream  sline ;
00122       std::getline(data,line) ; sline << line ; sline >>keyword ;
00123       //
00124       if( ( ! keyword.find(ATOM)) || ( ! keyword.find(HETETOATOM)) )
00125         {
00126           this->_numberOfAtoms++;keyword = "";
00127         }
00128     }
00129   std::cout << "  Number of atoms        "<< this->_numberOfAtoms <<std::endl;
00130   if (this->_numberOfAtoms == 0 ) 
00131     {std::cerr <<" Error in PDBstructure::readPDBFile, the number of atoms is 0" <<std::endl ;
00132     exit(1) ;
00133     }
00134   data.clear() ; 
00135   //
00136   int i = 0 , at = 0;
00137   data.seekg(0,std::ios::beg) ;
00138   this->_atoms = new PDBelement[this->_numberOfAtoms] ;
00139   while ( !data.eof() )
00140     { 
00141       std::stringstream  sline ;
00142       std::getline(data,line); sline << line ; sline >>line ;
00143       //
00144       // Tester si M est a la fin du mot line si oui pas
00145       //    de probleme si non, il faudait revenir en arriere pour
00146       //    lire le numero de l'atom.
00147       //
00148       unsigned int length = line.size() ;
00149       if( ( ! line.find(ATOM)) ){
00150         if ( length > ATOM.size() ) { 
00151           sline.seekg(ATOM.size(),std::ios::beg) ;
00152         }
00153         this->_atoms[i].readAtom(sline); i++;
00154       }
00155       else if  (! line.find(HETETOATOM)) 
00156         {
00157         if ( length > HETETOATOM.size() ) { 
00158           sline.seekg(HETETOATOM.size(),std::ios::beg) ;
00159         }
00160           this->_atoms[i].readAtom(sline); i++;
00161         }
00162       else if( ( line == CONECT) )
00163         {
00164           sline >> at ;  // a cause du C on commence a 0
00165           this->_atoms[at-1].readConnection(sline) ;
00166         }
00167     }
00168   //
00169   //  Fix the number of residus and fragments (here it is the same)
00170   //
00171   _numberOfResidus   = this->_atoms[_numberOfAtoms-1].numberOfResidu() +1 ; // a cause du C
00172   //
00173   _numberOfFragments = _numberOfResidus ;
00174   //
00175   return true ;
00176 }
00177 
00178 void PDBstructure::writeQCIFile(std::ofstream & out)
00179 {
00180   //
00181   // ecriture du header
00182   //
00183   int at ;
00184    out << "QCPP SYSTEM CHARGE"     << std::endl << this->_systemCharge     <<std::endl
00185        << "QCPP SPIN MULTIPLICITY" << std::endl << this->_spinMultiplicity << std::endl
00186        << "QCPP NB ATOMS"          << std::endl << this->_numberOfAtoms    << std::endl ;
00187    out << std::setprecision(4);
00188    //   out << setw(10) << std::ios::left;
00189    out<< setiosflags(std::ios::showpoint |std::ios::fixed ) ;
00190    for( int i = 0 ; i < this->_numberOfAtoms ; ++i){
00191      at = posAtomInString(this->_atoms[i]._atomName) ;
00192      out <<"   "<<this->_atoms[i]._atomName[at]
00193          << std::setw(13) << this->_atoms[i]._xyz[0]
00194          << std::setw(13) << this->_atoms[i]._xyz[1]
00195          << std::setw(13) << this->_atoms[i]._xyz[2] 
00196          << std::setw(8) << i + 1
00197          << "   " << this->_atoms[i]._fragName << "  ";
00198        out << this->_atoms[i]._connection.size() ;
00199        for(int j = 0 ; j < static_cast<int>(this->_atoms[i]._connection.size()) ; ++j){
00200          // +1 a cause du C 
00201          out << "  "<< ((this->_atoms[i]._connection))[j] +1 ;
00202        }
00203      out << std::endl ;
00204    }
00205 }
00206 void PDBstructure::writeXYZFile(std::ofstream & out){
00207   //
00208   // write header
00209   //
00210   int at ;
00211   out <<  this->_numberOfAtoms <<std::endl 
00212       << "Generated by PDBstructure::writeXYZFile " <<std::endl;
00213   // write atoms
00214   for( int i = 0 ; i < this->_numberOfAtoms ; ++i){
00215     at = posAtomInString(this->_atoms[i]._atomName) ;
00216     out <<"   "<<this->_atoms[i]._atomName[at]
00217         << std::setw(15) << this->_atoms[i]._xyz[0]
00218         << std::setw(15) << this->_atoms[i]._xyz[1]
00219         << std::setw(15) << this->_atoms[i]._xyz[2]
00220         << std::endl;
00221    }
00222 }
00223 void PDBstructure::writeXYZFileForAtoms(std::ofstream & out, std::vector<int> &index){
00224   //
00225   // write header
00226   //
00227   int at, num ;
00228   out <<  this->_numberOfAtoms <<std::endl 
00229       << "Generated by PDBstructure::writeXYZFile " <<std::endl;
00230   // write atoms
00231   for( unsigned int i = 0 ; i < index.size() ; ++i){
00232     num = index[i] ;
00233     at  = posAtomInString(this->_atoms[num]._atomName) ;
00234     out <<"   "<<this->_atoms[num]._atomName[at]
00235         << std::setw(15) << this->_atoms[num]._xyz[0]
00236         << std::setw(15) << this->_atoms[num]._xyz[1]
00237         << std::setw(15) << this->_atoms[num]._xyz[2]
00238         << std::endl;
00239    }
00240 }
00241 void PDBstructure::writePDBFileForAtoms(std::ofstream & out, std::vector<int> &index) 
00242 {
00243   int at, l, ideb ;
00244   std::string name("     "), un("1"), deux("2"), trois("3");
00245   for( unsigned int i = 0 ; i < index.size() ; ++i){
00246     ideb = 2 ;
00247     at   = index[i] ;
00248     l    = this->_atoms[at]._atomName.size() ;
00249     std::cout << name.size() << "   "<< l << std::endl;
00250     if ( l == 3 ){
00251       if ( this->_atoms[at]._atomName[0] == un[0] || 
00252            this->_atoms[at]._atomName[0] == deux[0]  || 
00253            this->_atoms[at]._atomName[0] == trois[0]  ) 
00254         { ideb = 1 ;}
00255       else
00256         {ideb = 2;}
00257     }
00258     else if( l == 4)
00259       {ideb = 1;}
00260     
00261     for(int j = 0 ; j < l ; ++j){
00262       name[ideb+j] = this->_atoms[at]._atomName[j] ;
00263     }
00264     //
00265     out << "ATOM " <<  std::setw(6) << i+1 << std::setw(5) <<  name<< " "<< this->_atoms[at]._residuName << std::setw(5);
00266     out <<this->_atoms[at]._numResiduInMolecule +1 << " \t" // +1 a cause du C
00267         << this->_atoms[at]._xyz[0] << " \t"<< this->_atoms[at]._xyz[1] <<" \t"<<  this->_atoms[at]._xyz[2] 
00268         << "   0.0 " <<std::endl;
00269     name = "     " ;
00270   }
00271   out <<"TER"<<std::endl;
00272 }
00273 std::string PDBstructure::getNumberInString(const int numero){
00274   int l , lori = 5;
00275   std::string num("00000") , s;
00276   std::stringstream ss ; ss << numero ; ss >> s; l = s.length() ;
00277   for(int j = 0 ;j < l ;++j)
00278     { num[lori-1-j] = s[l-1-j] ;}
00279   return num ;
00280 }
00281 bool PDBstructure::setFragName(FragmentsDataBase& base){
00282 
00283   return this->setFragName(base, FRAG); 
00284 
00285 }
00286 bool PDBstructure::setFragName(FragmentsDataBase& base, const std::string &type){
00287   //
00288   // std::cout <<"PDBstructure::setFragName(FragmentsDataBase& base, const std::string &type) " << std::endl;
00289   std::string atom, residu ;
00290   bool isInVector = false;
00291   std::vector<std::string> errorResidu ;
00292   bool isAA = (type == AA) ;
00293   for(int i = 0 ; i < this->_numberOfAtoms; ++i) 
00294     {
00295       std::string tmp ; 
00296       //
00297       if(isAA) {
00298         tmp = _atoms[i].getResiduName()  ;
00299       }
00300       else {
00301         //      tmp = base.findFragmentName(_atoms[i].getResiduName(), _atoms[i].getAtomName() ) ;
00302         if( ! base.getFragmentNameFrom(_atoms[i].getResiduName(), _atoms[i].getAtomName(),tmp ) ){
00303           //    The residu is not in the DataBase
00304           isInVector = false ;
00305           for( unsigned int j = 0 ; j <  errorResidu.size() ; ++j){
00306             std::cout << " ERROR : "<<_atoms[i].getResiduName()<< " :  " <<_atoms[i].getAtomName()<<std::endl;
00307             if( errorResidu[j] == _atoms[i].getResiduName()){
00308               isInVector = true ; break ;
00309             }
00310           }
00311           if(!isInVector) {
00312             errorResidu.push_back(_atoms[i].getResiduName()) ;}
00313         }
00314       }
00315       _atoms[i]._fragName.assign(tmp) ;
00316     }
00317   //
00318   // Gestion des erreurs :
00319   //
00320   if (errorResidu.size() > 0 ) {
00321     std::cerr << "The following residus are not in the Data Base : "<< base.nameOfDataBase() <<std::endl
00322               << "        " ;
00323     for( unsigned int j = 0 ; j <  errorResidu.size() ; ++j){
00324       std::cerr << errorResidu[j] << "  " ;
00325     }
00326     std::cerr << std::endl<<std::endl;
00327   }
00328   // On compte le nombre de residu ;
00329   std::string old ;
00330   old.assign( _atoms[0]._fragName) ;
00331   int  nbFrags = 1 , deb = 0 ;
00332   std::string tmp;
00333   if( !isAA){
00334     //
00335     //  Pour faire le Merge
00336     //
00337     bool dejaFait = false;
00338     int j ;
00339     std::string F1("F1"), TR("TR") ;
00340     //
00341     for(int i = 0 ; i < this->_numberOfAtoms; ++i) 
00342       {
00343         if( _atoms[i]._fragName == "F1" && _atoms[i]._residuName != "WAT" )
00344           { deb++;}
00345         else
00346           {break;}
00347       }
00348     for(int i = deb ; i < this->_numberOfAtoms; ++i) 
00349       {
00350         //      std::cout << i << " "<< _atoms[i]._fragName << std::endl ;
00351         if(_atoms[i]._fragName == "F1"   && _atoms[i]._residuName != "WAT") {
00352           _atoms[i]._fragName = TR ;
00353           if(! dejaFait){
00354             old  =_atoms[i-1]._fragName ; 
00355             dejaFait = true ;
00356             j = i-1 ;
00357             while ( old ==_atoms[j]._fragName) {
00358               _atoms[j]._fragName  = TR ;
00359               //            std::cout << "              "<< j << " "<< _atoms[j]._fragName << std::endl ;
00360               --j;
00361             }
00362           }
00363         }
00364         else {
00365           dejaFait = false ;
00366         }
00367       }
00368     //
00369     //  Affichage complet du nom du fragment
00370     //    nom-residu:nom_fragment-numero_fragment
00371     //    nom_residu-numero_residu:nom_fragment-numero_fragment
00372     //
00373     nbFrags = 0 ;  old  =_atoms[0]._fragName ; 
00374     int molWat = 0 ;
00375     //    tmp  = _atoms[0]._residuName + ":" +  _atoms[0]._fragName + "-" + this->getNumberInString(nbFrags)    ;
00376     
00377       tmp = _atoms[0]._residuName  + "-" +  this->getNumberInString(_atoms[0]._numResiduInMolecule) + 
00378         ":" +  _atoms[0]._fragName + "-" + this->getNumberInString(nbFrags) ;
00379 
00380     for(int i = 0 ; i < this->_numberOfAtoms; ++i) {
00381       if( _atoms[i]._residuName == "WAT" )
00382         { 
00383           if(molWat == 0 ){  nbFrags++ ;}
00384           molWat++;  
00385           if( molWat == 3){
00386             molWat = 0 ; 
00387           }
00388         }
00389       else if( old != _atoms[i]._fragName ){
00390         nbFrags++ ; old = _atoms[i]._fragName ;
00391       }
00392       //
00393       tmp = _atoms[i]._residuName  + "-" +  this->getNumberInString(_atoms[i]._numResiduInMolecule) + 
00394         ":" +  _atoms[i]._fragName + "-" + this->getNumberInString(nbFrags) ;
00395       _atoms[i]._fragName.assign(tmp) ;
00396     }
00397   }
00398   else {
00399     for(int i = 0 ; i < this->_numberOfAtoms; ++i) {
00400       tmp = _atoms[i]._residuName +  "-" +  
00401         this->getNumberInString(_atoms[i]._numResiduInMolecule) +":" +  _atoms[i]._fragName + "-" 
00402         + this->getNumberInString(_atoms[i]._numResiduInMolecule) ;
00403       nbFrags = std::max( _atoms[i]._numResiduInMolecule , nbFrags ) ;  
00404       _atoms[i]._fragName.assign(tmp) ;      
00405     }
00406     nbFrags++;
00407   }
00408     if(!isAA){_numberOfFragments = nbFrags + 1 ;}
00409   
00410   //
00411   std:: cout << "  Number of Residues     "<< _numberOfResidus   << std::endl ;
00412   std:: cout << "  Number of fragments    "<< _numberOfFragments << std::endl ;
00413   //
00414   return true ;
00415 }
00416 void PDBstructure::swap(const int i, const int  j){
00417   if( i <_numberOfAtoms && j <_numberOfAtoms){
00418     PDBelement tmp(_atoms[i])  ;
00419     _atoms[i].copy(_atoms[j]) ;
00420     _atoms[j].copy(tmp) ;
00421   }
00422   else{
00423     std::cerr << "void PDBstructure::swap" <<std::endl;
00424     exit(1) ;
00425   }
00426 }
00427 void PDBstructure::findMinMax(double *atomMin , double * atomMax  ){
00428   //
00429   atomMax[0] = atomMax[1] = atomMax[2] = -10000.0 ;
00430   atomMin[0] = atomMin[1] = atomMin[2] =  10000.0 ;
00431   //
00432   for(int i = 0 ; i < this->_numberOfAtoms ; ++i ){
00433     atomMax[0] = std::max(atomMax[0], this->_atoms[i]._xyz[0]) ;
00434     atomMax[1] = std::max(atomMax[1], this->_atoms[i]._xyz[1]) ;
00435     atomMax[2] = std::max(atomMax[2], this->_atoms[i]._xyz[2]) ;
00436     atomMin[0] = std::min(atomMin[0], this->_atoms[i]._xyz[0]) ;
00437     atomMin[1] = std::min(atomMin[1], this->_atoms[i]._xyz[1]) ;
00438     atomMin[2] = std::min(atomMin[2], this->_atoms[i]._xyz[2]) ;
00439   }
00440 }
00441 void PDBstructure::buildGravityCenter(double *G  ) {
00442 
00443   G[0] = G[1] = G[2] = 0.0 ;
00444   for(int i = 0 ; i < this->_numberOfAtoms ; ++i ){
00445     G[0] += this->_atoms[i]._xyz[0] ;
00446     G[1] += this->_atoms[i]._xyz[1] ;
00447     G[2] += this->_atoms[i]._xyz[2] ;
00448   }
00449   G[0] /= this->_numberOfAtoms ;
00450   G[1] /= this->_numberOfAtoms ;
00451   G[2] /= this->_numberOfAtoms ;
00452 }

Generated on Sat Jan 28 21:07:25 2006 for QC++ by  doxygen 1.4.4