00001
00002
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
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
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); }
00110
00111 sline << line ;
00112 this->_atoms[i].read(sline);
00113 numberMin = std::min(numberMin,this->_atoms[i].number()) ;
00114 }
00115
00116
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
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
00149
00150 std::string header;
00151 data >> this->_numberOfAtoms ;
00152 data >> header;
00153
00154 _atoms = new atomElement[this->_numberOfAtoms] ;
00155
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
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 void setOfAtoms::writeVMDFile(std::ofstream & out, const std::vector< std::string > &color ) {
00184
00185
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
00208 if(color[i] == color[atomLinked]){
00209
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
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
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
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
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
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
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
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
00321
00322 out << index.size() <<std::endl
00323 << "Generated by QC++::writeXYZFile(...,index) " <<std::endl;
00324
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
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
00350
00351 out << this->_numberOfAtoms <<std::endl
00352 << "Generated by QC++::writeXYZFile " <<std::endl;
00353
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
00404
00405
00406
00407
00408
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 }