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
00017
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
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
00049
00050 double gravityCenter[3], axe1[3], axe2[3], axe3[3] ;
00051
00052 buildInnerReperOfMolecule(atomsSet, gravityCenter, axe1, axe2, axe3 ) ;
00053
00054
00055
00056
00057 for(int i = 0 ; i < static_cast<int>(oldNumerotation.size()) ; ++i) {
00058 oldNumerotation[i] = i;
00059 }
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
00084
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 }
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
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
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
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
00181
00182
00183
00184
00185
00186
00187
00188 for(int i = 0 ; i < atoms.numberOfAtoms() ; ++i)
00189 {
00190
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
00206
00207
00208
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 }
00224
00225 }
00226 }
00227
00228
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
00238
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
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
00299
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
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
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);