Atoms.cpp

Go to the documentation of this file.
00001 //
00002 // Utilitaires  pour manipuler les fichiers d'atomes
00003 //
00004 #include "Atoms.h"
00005 //
00006 #include <algorithm>
00007 #include <iostream>
00008 #include <iomanip>
00009 #include <sstream>
00010 #include <stdlib.h>
00011 //
00012 setOfAtoms::setOfAtoms() : _numberOfAtoms(0),_numberOfAA(0),_atoms(NULL)
00013 {}
00014 
00015 setOfAtoms::setOfAtoms(const int size) : _numberOfAtoms(size),
00016                                          _atoms(new atomElement[size])
00017 {}
00018 setOfAtoms::~setOfAtoms() 
00019 {
00020   if(_atoms != NULL)
00021     {delete [] _atoms ;_atoms= NULL ;}
00022   _numberOfAtoms = 0 ; 
00023   
00024 }
00025 void setOfAtoms::reset() 
00026 {
00027   if(_atoms != NULL)
00028     {delete [] _atoms ;_atoms= NULL ;}
00029   _numberOfAtoms = 0 ; 
00030 }
00031 void setOfAtoms::setNumberOfAtoms(const int size) {
00032   _numberOfAtoms = size ; 
00033   _atoms         = new atomElement[size] ;
00034 }
00035 void setOfAtoms::findMinMax(double *atomMin , double * atomMax  ){
00036   //
00037   atomMax[0] = atomMax[1] = atomMax[2] = -10000.0 ;
00038   atomMin[0] = atomMin[1] = atomMin[2] =  10000.0 ;
00039   //
00040   for(int i = 0 ; i < this->_numberOfAtoms ; ++i ){
00041     atomMax[0] = std::max(atomMax[0], this->_atoms[i]._xyz[0]) ;
00042     atomMax[1] = std::max(atomMax[1], this->_atoms[i]._xyz[1]) ;
00043     atomMax[2] = std::max(atomMax[2], this->_atoms[i]._xyz[2]) ;
00044     atomMin[0] = std::min(atomMin[0], this->_atoms[i]._xyz[0]) ;
00045     atomMin[1] = std::min(atomMin[1], this->_atoms[i]._xyz[1]) ;
00046     atomMin[2] = std::min(atomMin[2], this->_atoms[i]._xyz[2]) ;
00047   }
00048 }
00049 void setOfAtoms::buildGravityCenter(double *G  ) {
00050 
00051   G[0] = G[1] = G[2] = 0.0 ;
00052   for(int i = 0 ; i < this->_numberOfAtoms ; ++i ){
00053     G[0] += this->_atoms[i]._xyz[0] ;
00054     G[1] += this->_atoms[i]._xyz[1] ;
00055     G[2] += this->_atoms[i]._xyz[2] ;
00056   }
00057   G[0] /= this->_numberOfAtoms ;
00058   G[1] /= this->_numberOfAtoms ;
00059   G[2] /= this->_numberOfAtoms ;
00060 }
00061 // Lecture des atomes dans un fichier .qco
00062 //
00063 //
00064 //
00065 int  setOfAtoms::findNumberOfAtoms(std::ifstream & data){
00066   //
00067   const std::string stringToFind("QCPP NB ATOMS") ;
00068   int numberOfAtoms = 0 ;
00069   bool isFind = false ;
00070   std::string line(" ") ;
00071   std::stringstream sline(" ") ;
00072   while ( ! isFind ){
00073     getline(data,line) ; 
00074     isFind = ! line.find(stringToFind);
00075   }
00076   if(isFind) { data >> numberOfAtoms; getline(data,line) ; }
00077   return numberOfAtoms ;
00078 }
00079 //
00080 bool setOfAtoms::readQCIFile(std::ifstream & data ){
00081   //
00082   //
00083   // il faut sauter les commentaires et lire le nombre d'atomes
00084   //
00085   this->_numberOfAtoms = findNumberOfAtoms(data) ;
00086   //
00087   if (this->_numberOfAtoms <1 ) {
00088     std::cerr <<" Error in setOfAtoms::readQCIFile, the number Of atoms is 0" <<std::endl ;
00089     exit(1) ;
00090     }
00091   _atoms = new atomElement[this->_numberOfAtoms] ;
00092   //
00093   bool isFind = true , useFRAG = false;
00094   std::string line, comment("//"), frag  ;
00095   std:: string AA("   "), oldAA ; 
00096   int numberMin = -1  ;
00097   //
00098   while ( isFind ){
00099     getline(data,line) ; 
00100     isFind = ! line.find(comment);
00101   }
00102   std::stringstream sline ; sline <<line ;
00103   this->_atoms[0].read(sline); 
00104   numberMin     = this->_atoms[0].number() ;
00105   //
00106   for (int i = 1 ; i < this->_numberOfAtoms ; ++i){
00107     std::stringstream sline ;
00108      getline(data,line) ; 
00109      while( line.size() ==0 ) {getline(data,line);    }  // autorise des lignes vides
00110     // getline(data,line);     
00111     sline << line ;
00112     this->_atoms[i].read(sline); 
00113     numberMin = std::min(numberMin,this->_atoms[i].number()) ;
00114   }
00115   //
00116   //  C numerotation for atoms (0,..this->_numberOfAtoms-1
00117   //
00118   if(numberMin == 1){
00119     for (int i = 0 ; i < this->_numberOfAtoms ; ++i){
00120       this->_atoms[i]._num -= 1 ;
00121     }
00122   }
00123   else if( numberMin == -1){
00124     for (int i = 0 ; i < this->_numberOfAtoms ; ++i){
00125       this->_atoms[i]._num = i  ;
00126     }
00127   }
00128   //
00129   for (int i = 0 ; i < this->_numberOfAtoms ; ++i){
00130     //
00131     //  Donne le nombre amines acides
00132     //
00133     frag = this->_atoms[i].fragmentName() ;
00134     useFRAG = ( frag != "DefaultGroup" ) ;
00135     if( useFRAG ){
00136       oldAA = AA ;
00137       AA.assign(frag,0,frag.find(":"));
00138       if( oldAA != AA) 
00139         { _numberOfAA++ ;}
00140     }
00141   }
00142   this->getNumberOfFragments() ;
00143   //
00144   return true ;
00145 }
00146 void setOfAtoms::readXYZFile(std::ifstream & data){
00147   //
00148   // write header
00149   //
00150   std::string header;
00151   data >>  this->_numberOfAtoms ;
00152   data >> header;
00153   // 
00154   _atoms = new atomElement[this->_numberOfAtoms] ;
00155   // read atoms
00156   for( int i = 0 ; i < this->_numberOfAtoms ; ++i){
00157     data >> this->_atoms[i]._name >> this->_atoms[i]._xyz[0] 
00158         >> this->_atoms[i]._xyz[1] >> this->_atoms[i]._xyz[2]  ;
00159     this->_atoms[i]._num = i  ;
00160    }
00161 }
00162 //
00163 //
00164 //
00165 //
00166 // void sortWithNumerotation(const int numberOfAtoms,  atomsStruct *pos , int * newNumerotation){
00167 //   //
00168 //   int newi, j, current;
00169 //   for(int i = 0 ; i < numberOfAtoms ; ++i) {
00170 //     newi = newNumerotation[i] ;
00171 //     while ( i != newi){
00172 //       swap(pos[i],pos[newi]); 
00173 //       newNumerotation[i]    = newNumerotation[newi] ; 
00174 //       newNumerotation[newi] = newi ; 
00175 //       // on continue
00176 //       newi = newNumerotation[i] ;
00177 //     }
00178 //   }
00179 // }
00180 
00181 //
00182 //
00183 void setOfAtoms::writeVMDFile(std::ofstream & out, const std::vector< std::string > &color ) {
00184   //
00185   // write header
00186   //
00187   out << "# partition file for VMD " <<std::endl 
00188       << "# Generated by atoms::writeVMDFile " <<std::endl
00189       << "draw materials off " <<std::endl;
00190   std::string previousColor(color[0]);
00191   out <<"draw color "<< color[0] <<std::endl ;
00192   for( int i = 0 ; i < this->_numberOfAtoms ; ++i){
00193     if( previousColor !=  color[i]  ) {
00194       out <<"draw color "<< color[i] <<std::endl ;
00195       previousColor =  color[i] ;
00196     }
00197     out <<"draw sphere {"
00198         << this->_atoms[i]._xyz[0]<< " "
00199         << this->_atoms[i]._xyz[1]<< " "
00200         << this->_atoms[i]._xyz[2]<< " } " << "radius 0.2 " << " resolution 6 " <<std::endl;
00201    }
00202   int atomLinked , previousNum = 0 ;
00203   for( int i = this->_numberOfAtoms -1; i > -1   ; --i){
00204     for( int j = 0 ; j < static_cast<int>(this->_atoms[i]._connect.size()) ; ++j ){
00205       atomLinked = (this->_atoms[i])._connect[j] - 1;
00206       if( i < atomLinked ){
00207       // il y a une liaison
00208         if(color[i] == color[atomLinked]){
00209           // meme partition
00210           if( previousColor !=  color[i]  ) {
00211             out <<"draw color "<< color[i] <<std::endl ;
00212             previousColor =  color[i] ;
00213           }
00214           out <<"draw cylinder  {"
00215               << this->_atoms[i]._xyz[0]<< " "
00216               << this->_atoms[i]._xyz[1]<< " "
00217               << this->_atoms[i]._xyz[2]<< " }  {" 
00218               << this->_atoms[atomLinked]._xyz[0]<< " "
00219               << this->_atoms[atomLinked]._xyz[1]<< " "
00220               << this->_atoms[atomLinked]._xyz[2]<< " } " 
00221               << "radius 0.05 " <<std::endl;
00222         }
00223         else {
00224           if( previousColor !=  color[i]  ) {
00225             out <<"draw color "<< color[i] <<std::endl ;
00226           }
00227           //      out << "draw color "<< color[i] <<std::endl;
00228           out <<"draw cylinder  {"
00229               << this->_atoms[i]._xyz[0]<< " "
00230               << this->_atoms[i]._xyz[1]<< " "
00231               << this->_atoms[i]._xyz[2]<< " }  {" 
00232               << 0.5 * (this->_atoms[atomLinked]._xyz[0]  +  this->_atoms[i]._xyz[0] )<< " "
00233               << 0.5 * (this->_atoms[atomLinked]._xyz[1]  +  this->_atoms[i]._xyz[1] )<< " "
00234               << 0.5 * (this->_atoms[atomLinked]._xyz[2]  +  this->_atoms[i]._xyz[2] )<< " } " 
00235               << "radius 0.05 " <<std::endl;
00236           out << "draw color "<< color[atomLinked] <<std::endl;
00237           out <<"draw cylinder  {"
00238               << this->_atoms[atomLinked]._xyz[0]<< " "
00239               << this->_atoms[atomLinked]._xyz[1]<< " "
00240               << this->_atoms[atomLinked]._xyz[2]<< " }  {" 
00241               << 0.5 * (this->_atoms[atomLinked]._xyz[0]  + this->_atoms[i]._xyz[0] )<< " "
00242               << 0.5 * ( this->_atoms[atomLinked]._xyz[1] + this->_atoms[i]._xyz[1] )<< " "
00243               << 0.5 * (this->_atoms[atomLinked]._xyz[2]  + this->_atoms[i]._xyz[2] )<< " } " 
00244               << "radius 0.05 " <<std::endl;
00245           previousColor =  color[atomLinked] ; 
00246           previousNum   = atomLinked ;
00247         }
00248       }
00249     }          
00250   }
00251   //
00252 }
00253 
00254 void setOfAtoms::writeQCIFile(std::ofstream & out){
00255   //
00256   // write header
00257   //
00258   out << "//"<<std::endl <<"// Generated by QC++::setOfAtoms::writeQCIFile "<<std::endl
00259       << "//La charge totale du systeme"<<std::endl <<"QCPP SYSTEM CHARGE"<<std::endl <<"0"<<std::endl
00260       << "//La multiplicite de spin    "<<std::endl <<"QCPP SPIN MULTIPLICITY"<<std::endl <<"1"<<std::endl
00261       << "//"<<std::endl <<"//Le nombre d'atomes"<<std::endl <<"QCPP NB ATOMS"<<std::endl ;
00262   out <<  this->_numberOfAtoms <<std::endl  ;
00263   // write atoms
00264   for( int i = 0 ; i < this->_numberOfAtoms ; ++i){
00265     out <<"   "<<this->_atoms[i]._name
00266         << std::setw(15) << this->_atoms[i]._xyz[0] 
00267         << std::setw(15) << this->_atoms[i]._xyz[1] 
00268         << std::setw(15) << this->_atoms[i]._xyz[2] 
00269         << std::setw(7) << this->_atoms[i]._num  + 1  // + 1 a cause de la numerotation C
00270         << std::setw(20) <<_atoms[i]._fragName << "  " ;
00271     //
00272     out << this->_atoms[i]._connect.size() ;
00273     for( unsigned int j = 0 ; j < this->_atoms[i]._connect.size() ; ++j ){
00274       out << "  " <<  (this->_atoms[i])._connect[j] ;
00275     }          
00276     out << std::endl;
00277   }
00278 }void setOfAtoms::writeQCIFile(std::ofstream & out, std::vector<int>& index){
00279   //
00280   // write header
00281   //
00282   std::cout << "WARNING PB with the connectivity"<<std::endl;
00283   out << "//"<<std::endl <<"// Generated by QC++::setOfAtoms::writeQCIFile "<<std::endl
00284       << "//La charge totale du systeme"<<std::endl <<"QCPP SYSTEM CHARGE"<<std::endl <<"0"<<std::endl
00285       << "//La multiplicite de spin    "<<std::endl <<"QCPP SPIN MULTIPLICITY"<<std::endl <<"1"<<std::endl
00286       << "//"<<std::endl <<"//Le nombre d'atomes"<<std::endl <<"QCPP NB ATOMS"<<std::endl ;
00287   out <<  index.size() <<std::endl     << "//  " <<std::endl;
00288   // write atoms
00289   int i ;
00290   for( unsigned int k = 0 ; k < index.size() ; ++k){
00291     i = index[k];
00292     out <<"   "<<this->_atoms[i]._name
00293         << std::setw(15) << this->_atoms[i]._xyz[0] 
00294         << std::setw(15) << this->_atoms[i]._xyz[1] 
00295         << std::setw(15) << this->_atoms[i]._xyz[2] 
00296         << std::setw(7)  <<  this->_atoms[i]._num  + 1  // + 1 a cause de la numerotation C
00297         << std::setw(20) <<_atoms[i]._fragName << "  " ;
00298     //
00299     out << this->_atoms[i]._connect.size() ;
00300     for( unsigned int j = 0 ; j < this->_atoms[i]._connect.size() ; ++j ){
00301       out << "  " <<  (this->_atoms[i])._connect[j] ;
00302     }          
00303     out << std::endl;
00304   }
00305 }
00306 void setOfAtoms::writeXYZFile(std::ofstream & out, std::vector<int>& index){
00307   //
00308   double H[3] ;
00309   H[0] = H[1] = H[2] = 0.0 ;
00310   this->writeXYZFile(out, H ,index) ;
00311 }
00312 void setOfAtoms::writeXYZFile(std::ofstream & out){
00313   //
00314   double H[3] ;
00315   H[0] = H[1] = H[2] = 0.0 ;
00316   this->writeXYZFile(out, H ) ;
00317 }
00318 void setOfAtoms::writeXYZFile(std::ofstream & out,double *H,const  std::vector<int>& index){
00319   //
00320   // write header
00321   //
00322   out <<  index.size() <<std::endl 
00323       << "Generated by QC++::writeXYZFile(...,index) " <<std::endl;
00324   // write atoms
00325   int i ;
00326   for( unsigned int j = 0 ; j < index.size() ; ++j){
00327     i = index[j] ;
00328     out <<"   "<<this->_atoms[i]._name
00329         << std::setw(15) << this->_atoms[i]._xyz[0] + H[0]
00330         << std::setw(15) << this->_atoms[i]._xyz[1] + H[1]
00331         << std::setw(15) << this->_atoms[i]._xyz[2] + H[2]
00332         << std::endl;
00333    }
00334 }
00335 void setOfAtoms::writeXYZPosition(std::ofstream & out,double *H){
00336   //
00337   // write header
00338   //
00339   for( int i = 0 ; i < this->_numberOfAtoms  ; ++i){
00340     out <<"   "<<this->_atoms[i]._name
00341         << std::setw(15) << this->_atoms[i]._xyz[0] + H[0]
00342         << std::setw(15) << this->_atoms[i]._xyz[1] + H[1]
00343         << std::setw(15) << this->_atoms[i]._xyz[2] + H[2]
00344         << std::endl;
00345    }
00346 }
00347 void setOfAtoms::writeXYZFile(std::ofstream & out,double *H){
00348   //
00349   // write header
00350   //
00351   out <<  this->_numberOfAtoms <<std::endl 
00352       << "Generated by QC++::writeXYZFile " <<std::endl;
00353   // write atoms
00354   for( int i = 0 ; i < this->_numberOfAtoms ; ++i){
00355     out <<"   "<<this->_atoms[i]._name
00356         << std::setw(15) << this->_atoms[i]._xyz[0] + H[0]
00357         << std::setw(15) << this->_atoms[i]._xyz[1] + H[1]
00358         << std::setw(15) << this->_atoms[i]._xyz[2] + H[2]
00359         << std::endl;
00360    }
00361 }
00362 
00363 int setOfAtoms::getNumberOfFragments()
00364 {
00365   //  
00366   std::string blanc(" ");
00367   
00368   std::string frag = _atoms[_numberOfAtoms - 1].getFragmentName() ;
00369   int pos = frag.find_last_of("-") + 1;
00370   for(int i = 0 ; i < pos ; ++i){
00371     frag[i] =  blanc[0] ;
00372   }
00373   std::stringstream  s ; s << frag ; 
00374   int number ;   s >> number ;
00375   _numberOfFragments = number;
00376   return number ;
00377 }
00378 
00379 void setOfAtoms::copyFrom(const int beginAtom, const int  EndAtom, const setOfAtoms & atoms) {
00380   if(_atoms == NULL || _numberOfAtoms < EndAtom - beginAtom ){
00381     std::cerr << "Error in setOfAtoms::copy(beginAtom,EndAtom,atoms) in " << __FILE__ <<std::endl
00382               << "                  " <<  _numberOfAtoms << " must be < " << EndAtom - beginAtom
00383               << "          or      " << _atoms << " must be allocated "<<std::endl ;
00384     exit(1) ;
00385   }
00386   for(int i = 0 ; i < EndAtom - beginAtom ; ++i){
00387     _atoms[i].copy(atoms.getElement(beginAtom + i) ) ; 
00388   }
00389 }
00390   //
00391 void setOfAtoms::swap(const int i, const int  j){
00392   if( i <_numberOfAtoms && j <_numberOfAtoms){
00393     atomElement tmp(_atoms[i])  ;
00394     _atoms[i].copy(_atoms[j]) ;
00395     _atoms[j].copy(tmp) ;
00396   }
00397   else{
00398     std::cerr << "void setOfAtoms::swap" <<std::endl;
00399     exit(1) ;
00400   }
00401 }
00402 //----------------------------------------------------------------------------//
00403 // Private Method: InertiaMatrixCompute
00404 //
00405 // Description: Calcul la matrice d'inertie d'un système de
00406 //              points (en 3D), et retourne le vecteur
00407 //              propre orthonormé associé à la plus petite
00408 //              valeur propre.
00409 //----------------------------------------------------------------------------//
00410 bool  setOfAtoms::buildIinertialMatrix(double *coordG) {
00411   std::cout << " Nothing to do here BUG " <<std::endl<<std::endl ;
00412   exit(1) ;
00413   return false;
00414 }
00415 
00416 //
00417 void setOfAtoms::fixConnectivity(const std::vector<int> & reOrder){
00418   //
00419   for (int i = 0 ; i < this->_numberOfAtoms ; ++i){
00420     for ( unsigned int j = 0 ; j < this->_atoms[i]._connect.size() ; ++j){
00421       this->_atoms[i]._connect[j] =  reOrder[this->_atoms[i]._connect[j]-1] + 1 ;
00422     }  
00423   }
00424   //
00425 }

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