00001
00002
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
00029 std::string atom ;
00030 sline >> _numInMolecule ;
00031 sline >> _atomName ;
00032
00033
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
00049 _numResiduInMolecule--;
00050
00051 _numInMolecule-- ;
00052 }
00053
00054 void PDBelement::readConnection(std::stringstream & sline)
00055 {
00056
00057
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
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
00145
00146
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 ;
00165 this->_atoms[at-1].readConnection(sline) ;
00166 }
00167 }
00168
00169
00170
00171 _numberOfResidus = this->_atoms[_numberOfAtoms-1].numberOfResidu() +1 ;
00172
00173 _numberOfFragments = _numberOfResidus ;
00174
00175 return true ;
00176 }
00177
00178 void PDBstructure::writeQCIFile(std::ofstream & out)
00179 {
00180
00181
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
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
00201 out << " "<< ((this->_atoms[i]._connection))[j] +1 ;
00202 }
00203 out << std::endl ;
00204 }
00205 }
00206 void PDBstructure::writeXYZFile(std::ofstream & out){
00207
00208
00209
00210 int at ;
00211 out << this->_numberOfAtoms <<std::endl
00212 << "Generated by PDBstructure::writeXYZFile " <<std::endl;
00213
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
00226
00227 int at, num ;
00228 out << this->_numberOfAtoms <<std::endl
00229 << "Generated by PDBstructure::writeXYZFile " <<std::endl;
00230
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"
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
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
00302 if( ! base.getFragmentNameFrom(_atoms[i].getResiduName(), _atoms[i].getAtomName(),tmp ) ){
00303
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
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
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
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
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
00360 --j;
00361 }
00362 }
00363 }
00364 else {
00365 dejaFait = false ;
00366 }
00367 }
00368
00369
00370
00371
00372
00373 nbFrags = 0 ; old =_atoms[0]._fragName ;
00374 int molWat = 0 ;
00375
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 }