tools.cpp

Go to the documentation of this file.
00001 
00002 #include <iostream>
00003 
00004 #include "tools.h"
00005 #include "commonDef.h"
00006 #include "pdbstruct.h"
00007 
00008 #if defined (_GNU_GSL_)
00009 #include <gsl/gsl_math.h>
00010 #include <gsl/gsl_eigen.h>
00011 #endif
00012 
00013 int  giveNumberOfOA(const std::string & atom) {
00014   int number = 0 , pos = 0 ;
00015   const std::string atomH("H");
00016 // const std::string heavyAtoms1("C,F,I,O,P,N,S,") , atomH("H");
00017 //  const std::string heavyAtoms2("Al,Be,Bi,Br,Cd,In,Mg,Pb,Hg,Sb,Si,Sn,Cl,Te,Tl,Zn,Ga,Ge,As,Se,") ;
00018   if( atom[0] >= '1' &&  atom[0] <= '9' ) 
00019     {pos++;}
00020   if( atom[pos] == atomH[0]){
00021     number = 1 ;}
00022   else{
00023     number = 4 ;
00024   }
00025   //  std::cout << "    atomName : " << atom << " the number of Orbital is   " << number << std::endl;
00026 
00027   return number ;
00028 }
00029 bool skipString(std::ifstream &data, const std::string &keyword) {
00030   //
00031   std::string line, key ;
00032   bool isFind = false ;
00033   //
00034   while ( !data.eof() ){
00035     std::stringstream sline ;
00036     std::getline(data,line) ; 
00037     sline << line ; 
00038     sline >> key ;
00039     if( keyword == key) 
00040       { isFind = true ; break; }
00041   } 
00042   //
00043   return isFind;
00044 } 
00045 
00046 void renumLinear(setOfAtoms   & atomsSet, std::vector<int>& oldNumerotation ){
00047   //
00048   //  First step we build the inner reper of the molecule
00049   //
00050   double gravityCenter[3],  axe1[3], axe2[3], axe3[3] ;
00051 
00052   buildInnerReperOfMolecule(atomsSet, gravityCenter, axe1, axe2, axe3 ) ;
00053   //
00054    //
00055   //  Trier les atomes dans ce nouveau repere
00056   //
00057   for(int i = 0 ; i < static_cast<int>(oldNumerotation.size()) ; ++i) {
00058     oldNumerotation[i] = i;
00059   }//for i
00060   //
00061   projectAtomsOnReper(atomsSet, gravityCenter, axe1, axe2, axe3);
00062   //
00063   sortInFirstDimension(atomsSet, oldNumerotation) ;
00064   //
00065 }
00066 template <typename SETOFATOMS>
00067 void buildInnerReperOfMolecule(SETOFATOMS   & atomsSet, double *gravityCenter,
00068                                double *axe1, double *axe2, double *axe3) {
00069   std::cout << "   buildInnerReperOfMolecule : " <<std::endl
00070             << "      nbAtoms "<< atomsSet.numberOfAtoms() <<std::endl;
00071   double inertialMatrix[3][3] ;
00072   //
00073   atomsSet.buildGravityCenter(gravityCenter) ;
00074   //
00075   std::cout<<"Centre de Gravite : ("<<gravityCenter[0]<<" ; "<<gravityCenter[1]<<" ; "
00076            <<gravityCenter[2]<<")"<<std::endl;
00077   //
00078   int i ;
00079   double Ixx, Iyy, Izz, Ixy, Ixz, Iyz ;
00080   Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0 ;
00081   for(i = 0 ; i < atomsSet.numberOfAtoms() ; ++i)
00082     {
00083       //      atomElement & atoms = atomsSet.getElement(i) ;
00084       //      std::cout << i  << " pos " <<  atoms._name<< " "<<   atoms._xyz[1] <<endl ;
00085       Ixx += ( atomsSet.getElement(i)._xyz[1] - gravityCenter[1])*( atomsSet.getElement(i)._xyz[1] - gravityCenter[1]) +
00086              ( atomsSet.getElement(i)._xyz[2] - gravityCenter[2])*( atomsSet.getElement(i)._xyz[2] - gravityCenter[2]);
00087       Iyy += ( atomsSet.getElement(i)._xyz[0] - gravityCenter[0])*( atomsSet.getElement(i)._xyz[0] - gravityCenter[0]) +
00088              ( atomsSet.getElement(i)._xyz[2] - gravityCenter[2])*( atomsSet.getElement(i)._xyz[2] - gravityCenter[2]);
00089       Izz += ( atomsSet.getElement(i)._xyz[0] - gravityCenter[0])*( atomsSet.getElement(i)._xyz[0] - gravityCenter[0]) +
00090              ( atomsSet.getElement(i)._xyz[1] - gravityCenter[1])*( atomsSet.getElement(i)._xyz[1] - gravityCenter[1]);
00091       Ixy -= ( atomsSet.getElement(i)._xyz[0] - gravityCenter[0])*( atomsSet.getElement(i)._xyz[1] - gravityCenter[1]);
00092       Ixz -= ( atomsSet.getElement(i)._xyz[0] - gravityCenter[0])*( atomsSet.getElement(i)._xyz[2] - gravityCenter[2]);
00093       Iyz -= ( atomsSet.getElement(i)._xyz[2] - gravityCenter[2])*( atomsSet.getElement(i)._xyz[1] - gravityCenter[1]);
00094     }//for i
00095   //
00096   inertialMatrix[0][0] = Ixx ;
00097   inertialMatrix[1][1] = Iyy ;
00098   inertialMatrix[2][2] = Izz ;
00099   inertialMatrix[0][1] = inertialMatrix[1][0] = Ixy ;
00100   inertialMatrix[0][2] = inertialMatrix[2][0] = Ixz ;
00101   inertialMatrix[1][2] = inertialMatrix[2][1] = Iyz ;
00102 //
00103   double eigenValues[3] ;
00104   int nbRoots = 3 ;
00105 
00106 #if defined (_GNU_GSL_)
00107   gsl_matrix_view m 
00108     = gsl_matrix_view_array(&inertialMatrix[0][0], 3, 3);
00109   gsl_vector *eval = gsl_vector_alloc (3);
00110   gsl_matrix *evec = gsl_matrix_alloc (3, 3);
00111 
00112   gsl_eigen_symmv_workspace * w = 
00113     gsl_eigen_symmv_alloc (3);
00114   
00115   gsl_eigen_symmv (&m.matrix, eval, evec, w);
00116 
00117   gsl_eigen_symmv_free(w);
00118 
00119   gsl_eigen_symmv_sort (eval, evec, 
00120                         GSL_EIGEN_SORT_VAL_ASC);
00121    for(i = 0 ; i < 3 ; ++i ){
00122      eigenValues[i] = gsl_vector_get(eval, i) ;
00123    }
00124    int j ;
00125   gsl_vector_view evec_i = gsl_matrix_column(evec, 0);
00126   for ( j = 0 ; j < 3 ; ++j){
00127     axe1[j] = gsl_vector_get(&evec_i.vector, j) ;
00128     //    _axe1[j] = gsl_matrix_get(m, i, j); gsl_vector_get(evec_i, j) ;
00129   }
00130    evec_i = gsl_matrix_column(evec, 1);
00131    for ( j = 0 ; j < 3 ; ++j){
00132      axe2[j] = gsl_vector_get(&evec_i.vector, j) ;
00133    }
00134    evec_i = gsl_matrix_column(evec, 2);
00135    for( j = 0 ; j < 3 ; ++j){
00136      axe3[j] = gsl_vector_get(&evec_i.vector, j) ;
00137    }
00138 #else
00139    //
00140    std::cerr <<"void buildInnerReperOfMolecule must be compiled with GSL"<<std::endl;
00141    exit(1) ;
00142 
00143 #endif 
00144   //
00145   std::cout << "Les valeurs Propres et Vecteurs propres : "<< nbRoots<<std::endl
00146             << "    "<< eigenValues[0]<< "  "<<axe1[0]<< "  "<< axe1[1]<< "  "<< axe1[2]<<std::endl
00147             << "    "<< eigenValues[1]<< "  "<<axe2[0]<< "  "<< axe2[1]<< "  "<< axe2[2]<<std::endl
00148             << "    "<< eigenValues[2]<< "  "<<axe3[0]<< "  "<< axe3[1]<< "  "<< axe3[2]<<std::endl
00149             << std::endl;
00150   //
00151   //  Orientation du repere
00152   //
00153 //   double dot ;
00154 //   double vect[3] ;
00155 //   dot      = axe1[0] * axe2[0] + axe1[1] * axe2[1] + axe1[2] * axe2[2] ;
00156 //   std::cout << "axe1 . axe2 : " << dot <<std::endl ;
00157 //   dot      = axe1[0] * axe3[0] + axe1[1] * axe3[1] + axe1[2] * axe3[2] ;
00158 //   std::cout << "axe1 . axe3 : " << dot <<std::endl ;
00159 //   dot      = axe3[0] * axe2[0] + axe3[1] * axe2[1] + axe3[2] * axe2[2] ;
00160 //   std::cout << "axe3 . axe2 : " << dot <<std::endl ;
00161 //   //
00162 //   vect[0]  = axe1[1] * axe2[2] - axe1[2] * axe2[1] ;
00163 //   vect[1]  = axe1[2] * axe2[0] - axe1[0] * axe2[2] ;
00164 //   vect[2]  = axe1[0] * axe2[1] - axe1[1] * axe2[0] ;
00165 //   //
00166 //   dot      = vect[0] * axe3[0] + vect[1] * axe3[1] + vect[2] * axe3[2] ;
00167 //   std::cout << "vect . axe3 : " << dot <<std::endl ;
00168 //   std::cout << "New Reper "<<std::endl
00169 //          << "  i    "<<axe1[0]<< "  "<< axe1[1]<< "  "<< axe1[2]<<std::endl
00170 //          << "  j    "<<axe2[0]<< "  "<< axe2[1]<< "  "<< axe2[2]<<std::endl
00171 //          << "  k    "<<axe3[0]<< "  "<< axe3[1]<< "  "<< axe3[2]<<std::endl
00172 //          << "  k    "<<vect[0]<< "  "<< vect[1]<< "  "<< vect[2]<<std::endl
00173 //          << std::endl;
00174   
00175 }
00176 template <typename SETOFATOMS>
00177 void projectAtomsOnReper(SETOFATOMS& atoms, double* gravityCenter, double* axe1, double* axe2, double* axe3){
00178 //
00179   double s1,s2,s3 ;
00180 //   s1 =  (atoms.getElement(0)._xyz[0]- gravityCenter[0])*axe1[0] +  
00181 //         (atoms.getElement(0)._xyz[1]- gravityCenter[1])*axe1[1] +  
00182 //         (atoms.getElement(0)._xyz[2]- gravityCenter[2])*axe1[2] ;
00183 //   if(s1 > 0) 
00184 //     { 
00185 //       std::cout << "Il y a un probleme le repere est mal oriente" << std::endl;
00186 //       axe1[0] = -axe1[0] ;axe1[1] = -axe1[1] ;axe1[2] = -axe1[2] ;
00187 //     }
00188   for(int i = 0 ; i < atoms.numberOfAtoms() ; ++i)
00189     {
00190       //      atomElement & atom = atoms.getElement(i) ;
00191       atoms.getElement(i)._xyz[0] -= gravityCenter[0];
00192       atoms.getElement(i)._xyz[1] -= gravityCenter[1];
00193       atoms.getElement(i)._xyz[2] -= gravityCenter[2];
00194       //
00195       s1            = atoms.getElement(i)._xyz[0]*axe1[0] +  atoms.getElement(i)._xyz[1]*axe1[1] +  atoms.getElement(i)._xyz[2]*axe1[2] ;
00196       s2            = atoms.getElement(i)._xyz[0]*axe2[0] +  atoms.getElement(i)._xyz[1]*axe2[1] +  atoms.getElement(i)._xyz[2]*axe2[2] ;
00197       s3            = atoms.getElement(i)._xyz[0]*axe3[0] +  atoms.getElement(i)._xyz[1]*axe3[1] +  atoms.getElement(i)._xyz[2]*axe3[2] ;
00198       atoms.getElement(i)._xyz[0]  = s1 ;
00199       atoms.getElement(i)._xyz[1]  = s2 ;
00200       atoms.getElement(i)._xyz[2]  = s3 ;
00201     }
00202 }
00203 void sortInFirstDimension(setOfAtoms& atoms, std::vector<int>&  oldNum) {
00204 
00205   // A noter que même si la variable nbAtoms n'a pas besoin d'être un
00206   // argument d'entrée, nous la mettons comme telle, car le partitionneur
00207   // ne restera pas forcément dans la classe CQCAtoms (où la variable
00208   // nbAtoms est connue).
00209 
00210   int      i, j;
00211   int      temp1;
00212   double   temp3;
00213 
00214   for(i = 1 ; i <  atoms.numberOfAtoms() ; ++i) {
00215     temp3 = atoms.getElement(i)._xyz[0] ;
00216     for(j = 0 ; j < i ; ++j) {
00217 
00218       if (temp3 < (atoms.getElement(j))._xyz[0] ){
00219         atoms.swap(i,j) ;
00220         temp1     = oldNum[j];
00221         oldNum[j] = oldNum[i];
00222         oldNum[i] = temp1;
00223       }//if temp3
00224 
00225     }//for j
00226   }//for i
00227   //
00228   //  Il faut remettre la connectivity avec les bons numeros ;
00229   //
00230   std::vector<int>  inverseNewNum(oldNum.size());
00231   for(i = 0 ; i <  static_cast<int>(oldNum.size()) ; ++i) {
00232     inverseNewNum[oldNum[i]] = i ;
00233   }
00234   //
00235   //
00236   //
00237   //  oldNum             : l'ancienne position dans le tableau
00238   //  inverseNewNum      :  
00239   //
00240    for(i = 0 ; i <  static_cast<int>(oldNum.size()) ; ++i) {
00241      std::cout << i << "   "<<  (atoms.getElement(i))._num << "   " 
00242               << oldNum[i]<< "   " << inverseNewNum[i] <<std::endl ;
00243    }
00244   //
00245   atoms.fixConnectivity(inverseNewNum) ;
00246 }
00247 void readPartitionAndSetColor(std::ifstream& part ,
00248                               std::vector<std::string>& color){
00249 
00250   std::vector<std::string> colorRef(10) ;
00251   std::string sharedZoneColor("gray") ;
00252   //
00253   colorRef[0] = "red" ;
00254   colorRef[1] = "blue" ;
00255   colorRef[2] = "green" ;
00256   colorRef[3] = "yellow" ;
00257   colorRef[4] = "gold" ;
00258   colorRef[5] = "purple" ;
00259   colorRef[6] = "orchid" ;
00260   colorRef[7] = "pink" ;
00261   colorRef[8] = "khaki" ;
00262   colorRef[9] = "white" ;
00263   // ouverture du fichier
00264   //
00265   part.clear() ;
00266   part.seekg(0,std::ios::beg) ;
00267   //
00268   std::string line ;
00269   int _nbOfbSubDomains=0 , _nbAtoms =0;
00270   //
00271   skipString(part,"NUMBER") ;
00272   part >> _nbOfbSubDomains ;  
00273   skipString(part,"NUMBER") ;
00274   part >> _nbAtoms ;  
00275   std::cout << "NUMBER OF SUBDOMAINS : " << _nbOfbSubDomains << std::endl
00276             << "NUMBER OF ATOMS      : " << _nbAtoms << "  " << std::endl;
00277   if(static_cast<int>(color.size()) != _nbAtoms )
00278     { exit(1) ;}
00279   skipString(part,"PARTITION") ;
00280   skipString(part,"PARTITION") ;
00281   int numAt, nbSD, numSD,zoneIJ, localOrder;
00282   std::ofstream out ;
00283   out.open("colorAtoms.tmp");
00284   bool debug ;
00285    for ( int i = 0 ; i < _nbAtoms ; ++i) {
00286      debug = false ;
00287      part >> numAt >>  nbSD  ;
00288      for(int j = 0 ; j < nbSD ; ++j) {
00289         part >> numSD >> zoneIJ >> localOrder ;
00290         if(  zoneIJ ==  Overlap_Zone ){
00291           color[i] = sharedZoneColor ; std::getline(part,line);
00292           debug =  true ;
00293           break ;
00294         } 
00295         else if (  zoneIJ ==  Core_Zone ){
00296           debug =  true ;
00297           color[i] = colorRef[numSD] ;  
00298 //        std::getline(part,line);
00299 //        break ;
00300         }
00301       }
00302      if ( ! debug){
00303        std::cout << "BUG COLOR is not affected " << i << "  "<<  numSD << "  "<< zoneIJ <<"  "<< localOrder <<std::endl;
00304      }
00305      out << i <<"     -"<< color[i]<<"-"<< std::endl ;
00306    }
00307    out.close();
00308 }
00309 
00310 // double  distanceBoxToBox(double * minBox, double * maxBox1, double * minBox2, double * maxBox2){
00311 //   //
00312 //   double dist = 0.0 ;
00313 //   double length1[3],length2[3]; 
00314 //   int ind1[3], ind2[3]
00315 //   int i, j ;
00316 //   //
00317 //   for(i = 0 ; i < 3 ; ++i) 
00318 //     {
00319 //       length1[i] = maxBox1[i] -  minBox1[i] ;
00320 //       length2[i] = maxBox2[i] -  minBox2[i] ;
00321 //       ind1[i]    = ind2[i] = -1 ;
00322 //     }
00323 //   for(i = 0 ; i < 3 ; ++i) {
00324 //     maxBox1[i] -  minBox2[i] ;
00325 //   }
00326 //   //
00327 //   return dist ;
00328 // }
00329 template  void buildInnerReperOfMolecule< setOfAtoms   >(setOfAtoms  & atomsSet, double *gravityCenter,
00330                                                        double *axe1, double *axe2, double *axe3);
00331 
00332 
00333 template  void buildInnerReperOfMolecule< PDBstructure >(PDBstructure  & atomsSet, double *gravityCenter,
00334                                                          double *axe1, double *axe2, double *axe3);
00335 //
00336 template void projectAtomsOnReper< setOfAtoms   >( setOfAtoms & atoms, double* gravityCenter, 
00337                                                  double* axe1, double* axe2, double* axe3);
00338 
00339 template void projectAtomsOnReper< PDBstructure >( PDBstructure & atoms, double* gravityCenter, 
00340                                                    double* axe1, double* axe2, double* axe3);

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