fpocket: tpocket.c Source File

00001 #include "../headers/tpocket.h" 
00002 
00003 /*
00004 
00005 ## GENERAL INFORMATION
00006 ##
00007 ## FILE                                         tpocket.c
00008 ## AUTHORS                                      P. Schmidtke and V. Le Guilloux
00009 ## LAST MODIFIED                        01-04-08
00010 ##
00011 ## SPECIFICATIONS
00012 ##
00013 ##  This file contains all routines used to execute tpocket and test the fpocket
00014 ##  performance, given a set of complexe-apo-ligand structure stored in a single
00015 ##  text file.
00016 ##
00017 ## MODIFICATIONS HISTORY
00018 ##
00019 ##      10-03-09        (v)  Mean number of atom per pocket + ligand volume added
00020 ##      05-03-09        (v)  Mean pocket volume added to tpocket output
00021 ##      19-01-09        (v)  Minor modif (status output printed on the same line)
00022 ##      28-11-08        (v)  Relooking + comments UTD
00023 ##      27-11-08        (v)  Write fpocket output if asked + relooking
00024 ##      15-10-08        (v)  Fixed bad count when PDB file not found.
00025 ##      01-04-08        (v)  Added template for comments and creation of history
00026 ##      01-01-08        (vp) Created (random date...)
00027 ##      
00028 ## TODO or SUGGESTIONS
00029 ##
00030 ##      (v) CLEAN THE CODE! It's such a mess here...
00031 ##
00032 
00033 */
00034 
00035 /*
00036     COPYRIGHT DISCLAIMER
00037 
00038     Vincent Le Guilloux, Peter Schmidtke and Pierre Tuffery, hereby
00039         disclaim all copyright interest in the program “fpocket” (which
00040         performs protein cavity detection) written by Vincent Le Guilloux and Peter
00041         Schmidtke.
00042 
00043     Vincent Le Guilloux  28 November 2008
00044     Peter Schmidtke      28 November 2008
00045     Pierre Tuffery       28 November 2008
00046 
00047     GNU GPL
00048 
00049     This file is part of the fpocket package.
00050 
00051     fpocket is free software: you can redistribute it and/or modify
00052     it under the terms of the GNU General Public License as published by
00053     the Free Software Foundation, either version 3 of the License, or
00054     (at your option) any later version.
00055 
00056     fpocket is distributed in the hope that it will be useful,
00057     but WITHOUT ANY WARRANTY; without even the implied warranty of
00058     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00059     GNU General Public License for more details.
00060 
00061     You should have received a copy of the GNU General Public License
00062     along with fpocket.  If not, see <http://www.gnu.org/licenses/>.
00063 
00064 **/
00065 
00066 /**
00067    ## FUNCTION: 
00068         test_fpocket
00069   
00070    ## SPECIFICATION: 
00071         Test fpocket for a set of pdb files. The output is writen in this function,
00072         and consists of two files (see documentation for details)
00073   
00074    ## PARAMETRES:
00075         @ s_tparams *par: Parameters, contain the fpocket parameters, and the list
00076                                           of data set information.
00077   
00078    ## RETURN:
00079         void
00080   
00081 */
00082 void test_fpocket(s_tparams *par)
00083 {
00084         if(! par || par->nfiles <= 0)  return ;
00085         
00086         int nranks = 13,  novlp = 5, novlp2 = 9,novlp3 = 8,
00087                 i, j, k, nok = 0  ;
00088         int ranks [] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 50} ;
00089         float ovlp[] = {50.0, 60.0, 70.0, 80.0,90.0} ;
00090         float ovlp2[] = {0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9, 1.0} ;
00091         float ovlp3[] = {0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7} ;
00092         
00093         float mean_ov1 = 0.0, mean_ov2 = 0.0, mean_dst = 0.0, mean_ov4 = 0.0, mean_ov5 = 0.0,
00094                   mean_ovr1 = 0.0, mean_ovr2 = 0.0, mean_ovr3 = 0.0, mean_ovr4 = 0.0, mean_ovr5 = 0.0, mean_ovr6 = 0.0,
00095                   mean_pvol3 = 0.0, mean_pvol6 = 0.0 ;
00096 
00097         float mean_lvol = 0.0 ;
00098         int mean_nbatm3 = 0, mean_nbatm6 = 0 ;
00099         
00100         int n1, n2, n3, n4, n5, n6, N = 0 ;
00101 
00102         /* Store statistics for each files: */
00103         int status[par->nfiles] ;
00104         float ddata [par->nfiles][M_NDDATA];
00105         int idata [par->nfiles][M_NIDATA];
00106 
00107         /* Test all files */
00108         for(i = 0 ; i < par->nfiles ; i++) {
00109                 strcpy(par->fpar->pdb_path, par->fapo[i]) ;
00110                 status [i] = test_set(par, i, ddata, idata) ;
00111                 fprintf(stdout, "> %3d : %s output code %d", i+1, par->fapo[i],
00112                                                                                                            status [i]) ;
00113                 if(i == par->nfiles - 1) fprintf(stdout, "\n") ;
00114                 else fprintf(stdout, "\r") ;
00115                 fflush(stdout) ;
00116 
00117                 if(status[i] == M_OK) N++ ;
00118         }
00119 
00120         /* Printing each complexe statistics */
00121         FILE *fp = fopen(par->p_output, "w") ;
00122         if(fp) {
00123                 n1 = 0, n2 = 0, n3 = 0, n4 = 0, n5 = 0, n6 = 0 ;
00124                 fprintf(fp, "LIG | COMPLEXE | APO | NB_PCK | CRIT1 | CRIT2 | CRIT3 | CRIT4 | CRIT5 | CRIT6 | POS1 | POS2 | POS3 | POS4 | POS5 | POS6 | REL_OVLP1 | REL_OVLP2 | REL_OVLP3 | REL_OVLP4 | REL_OVLP5 | REL_OVLP6 | LIGMASS | LIGVOL | PVOL3 | NATM3 | PVOL6 | NATM6 \n") ;
00125                 for(i = 0 ; i < par->nfiles ; i++) {
00126                         remove_path(par->fcomplex[i]) ;
00127                         remove_path(par->fapo[i]) ;
00128                         if(status[i] == M_OK) {
00129                                 fprintf(fp, "%s %s %s %5d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %4d %4d %4d %4d %4d %4d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %9.2f %9.2f %12.2f %4d %12.2f %4d\n",
00130                                                 par->fligan[i], par->fcomplex[i], par->fapo[i], idata[i][M_NPOCKET],
00131                                                 ddata[i][M_MAXPCT1], ddata[i][M_MAXPCT2], ddata[i][M_MINDST], ddata[i][M_CRIT4], ddata[i][M_CRIT5], ddata[i][M_CRIT6],
00132                                                 idata[i][M_POS1], idata[i][M_POS2], idata[i][M_POS3], idata[i][M_POS4], idata[i][M_POS5], idata[i][M_POS6],
00133                                                 ddata[i][M_OREL1], ddata[i][M_OREL2], ddata[i][M_OREL3], ddata[i][M_OREL4], ddata[i][M_OREL5], ddata[i][M_OREL6],
00134                                                 ddata[i][M_LIGMASS], ddata[i][M_LIGVOL], ddata[i][M_POCKETVOL_C3], idata[i][M_NATM3], ddata[i][M_POCKETVOL_C6], idata[i][M_NATM6]) ;
00135 
00136                                 if(idata[i][M_POS1] > 0) {
00137                                         mean_ov1 += ddata[i][M_MAXPCT1] ;
00138                                         mean_ovr1 += ddata[i][M_OREL1];
00139                                         n1++ ;
00140                                 }
00141 
00142                                 if(idata[i][M_POS2] > 0) {
00143                                         mean_ov2 += ddata[i][M_MAXPCT2] ;
00144                                         mean_ovr2 += ddata[i][M_OREL2];
00145                                         n2++ ;
00146                                 }
00147 
00148                                 if(idata[i][M_POS3] > 0) {
00149                                         mean_dst += ddata[i][M_MINDST];
00150                                         mean_ovr3 += ddata[i][M_OREL3];
00151                                         mean_pvol3 += ddata[i][M_POCKETVOL_C3] ;
00152                                         mean_nbatm3 += idata[i][M_NATM3] ;
00153                                         n3 ++ ;
00154                                 }
00155 
00156                                 if(idata[i][M_POS4] > 0) {
00157                                         mean_ov4 += ddata[i][M_CRIT4];
00158                                         mean_ovr4 += ddata[i][M_OREL4];
00159                                         n4 ++ ;
00160                                 }
00161 
00162                                 if(idata[i][M_POS5] > 0) {
00163                                         mean_ov5 += ddata[i][M_CRIT5];
00164                                         mean_ovr5 += ddata[i][M_OREL5];
00165                                         n5 ++ ;
00166                                 }
00167 
00168                                 if(idata[i][M_POS6] > 0) {
00169                                         mean_ovr6 += ddata[i][M_OREL5];
00170                                         mean_pvol6 += ddata[i][M_POCKETVOL_C6] ;
00171                                         mean_nbatm6 += idata[i][M_NATM6] ;
00172                                         n6 ++ ;
00173                                 }
00174                         }
00175                         else {
00176                                 fprintf(fp, "%s %s %s %5d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %4d %4d %4d %4d %4d %4d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %9.2f %9.2f %12.2f %4d %12.2f %4d \n",
00177                                                 par->fligan[i], par->fcomplex[i], par->fapo[i], -1, 
00178                                                 -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
00179                                                 -1, -1, -1, -1, -1, -1,
00180                                                 -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
00181                                                 -1.0, -1.0, -1.0, -1, -1.0, -1) ;
00182                         }
00183                         mean_lvol += ddata[i][M_LIGVOL] ;
00184                 }
00185 
00186                 mean_ov1 /= (float) n1 ; mean_ovr1 /= (float) n1 ;
00187                 mean_ov2 /= (float) n2 ; mean_ovr2 /= (float) n2 ;
00188 
00189                 mean_dst /= (float) n3 ; mean_ovr3 /= (float) n3 ;
00190                 mean_pvol3 /= (float) n3 ;
00191                 mean_nbatm3 /= (float) n3 ;
00192 
00193                 mean_ov4 /= (float) n4 ; mean_ovr4 /= (float) n4 ;
00194                 mean_ov5 /= (float) n5 ; mean_ovr5 /= (float) n5 ;
00195                 
00196                 mean_ovr6 /= (float) n6 ;
00197                 mean_pvol6 /= (float) n6 ;
00198                 mean_nbatm6 /= (float) n6 ;
00199 
00200                 mean_lvol /= (float) par->nfiles ;
00201 
00202                 fclose(fp) ;
00203         }
00204         else fprintf(stdout, "The file %s could not be opened\n", par->p_output) ;
00205 
00206         /* Printing global statistics */
00207         FILE *fg = fopen(par->g_output, "w") ;
00208         if(fg) {
00209                 /* Write the first criteria statistics */
00210                 fprintf(fg, "===================== General statistics on all complexes =======================\n") ;
00211 
00212 
00213                 /* Write the first criteria statistics */
00214                 fprintf(fg, "\n\t--\n") ;
00215                 fprintf(fg,   "\t-                       _ Distance criteria _                      -\n") ;
00216                 fprintf(fg,   "\t--\n\n") ;
00217                 fprintf(fg, "   Ratio of good predictions (dist = 4A) \n") ;
00218                 fprintf(fg, "\n") ;
00219 
00220                 for(i = 0 ; i < nranks ; i++) {
00221                         nok = 0 ;
00222                         for(k = 0 ; k < par->nfiles ; k++) {
00223                                 if(status[k] == M_OK && idata[k][M_POS3] <= ranks[i]
00224                                 && idata[k][M_POS3] > 0) {
00225                                         nok ++ ;
00226                                 }
00227                         }
00228                         //printf("%d NOK: %d, %d, %f\n", ranks[i], nok, N, ((float)nok) / ((float) N)) ;
00229                         fprintf(fg, "Rank <= %2d  :\t\t%6.2f\n", ranks[i],
00230                                         ((float)nok) / ((float) N)) ;
00231                 }
00232 
00233                 fprintf(fg, "-\n") ;
00234                 fprintf(fg, "Mean distance                   : %8.2f\n", mean_dst) ;
00235                 fprintf(fg, "Mean relative overlap           : %8.2f\n", mean_ovr3) ;
00236                 fprintf(fg, "Mean pocket volume (estimation) : %8.2f\n", mean_pvol3) ;
00237                 fprintf(fg, "Mean ligand volume (estimation) : %8.2f\n", mean_lvol) ;
00238                 fprintf(fg, "Mean number of pocket atom      : %4d\n", mean_nbatm3) ;
00239 
00240 
00241                 /* Write the 2nd criteria statistics */
00242                 fprintf(fg, "\n\t--\n") ;
00243                 fprintf(fg,   "\t- _ 1st overlap criteria (use of ligand's alpha sphere neighbors)_ -\n") ;
00244                 fprintf(fg,   "\t--\n\n") ;
00245                 fprintf(fg,"           :");
00246                 for( j = 0 ; j < novlp ; j++) {
00247                                 fprintf(fg, "  >%5.2f  :", ovlp[j]) ;
00248                 }
00249                 fprintf(fg,"\n");
00250                 fprintf(fg,"--");
00251                 for( j = 0 ; j < novlp ; j++) {
00252                         fprintf(fg, "-") ;
00253                 }
00254                 fprintf(fg,"\n");
00255                 
00256                 for(i = 0 ; i < nranks ; i++) {
00257                         fprintf(fg, "Rank <= %2d :", ranks[i]) ;
00258                         for( j = 0 ; j < novlp ; j++) {
00259                                 nok = 0 ;
00260                                 for(k = 0 ; k < par->nfiles ; k++) {
00261                                         if( status[k] == M_OK &&
00262                                                 ddata[k][M_MAXPCT2] >= ovlp[j] &&
00263                                                 idata[k][M_POS2] <= ranks[i] &&
00264                                                 idata[k][M_POS2] > 0) {
00265                                                 nok ++ ;
00266                                         }
00267                                 }
00268                                 fprintf(fg, "  %6.2f  :", ((float)nok) / ((float) N)) ;
00269                         }
00270                         fprintf(fg, "\n") ;
00271                 }
00272                 fprintf(fg, "-\n") ;
00273                 fprintf(fg, "Mean overlap          : %8.2f\n", mean_ov2) ;
00274                 fprintf(fg, "Mean relative overlap : %8.2f\n", mean_ovr2) ;
00275 
00276                 /* Write the 3rd criteria statistics */
00277                 
00278                 fprintf(fg, "\n\t--\n") ;
00279                 fprintf(fg,   "\t-        _ 2nd overlap criteria (simple distance criteria) _       -\n") ;
00280                 fprintf(fg,   "\t--\n\n") ;
00281                 fprintf(fg,"           :");
00282 
00283                 for( j = 0 ; j < novlp ; j++) {
00284                                 fprintf(fg, "  >%5.2f  :", ovlp[j]) ;
00285                 }
00286 
00287                 fprintf(fg,"\n");
00288                 fprintf(fg,"---");
00289 
00290                 for( j = 0 ; j < novlp ; j++) {
00291                         fprintf(fg, "-") ;
00292                 }
00293                 fprintf(fg,"\n");
00294 
00295                 for(i = 0 ; i < nranks ; i++) {
00296                         fprintf(fg, "Rank <= %2d :", ranks[i]) ;
00297                         for( j = 0 ; j < novlp ; j++) {
00298                                 nok = 0 ;
00299                                 for(k = 0 ; k < par->nfiles ; k++) {
00300                                         if( status[k] == M_OK
00301                                                 && ddata[k][M_MAXPCT1] >= ovlp[j]
00302                                                 && idata[k][M_POS1] <= ranks[i]
00303                                                 && idata[k][M_POS1] > 0) {
00304                                                 nok ++ ;
00305                                         }
00306                                 }
00307                                 fprintf(fg, "  %6.2f  :", ((float)nok) / ((float) N)) ;
00308                         }
00309                         fprintf(fg, "\n") ;
00310                 }
00311                 fprintf(fg, "-\n") ;
00312                 fprintf(fg, "Mean overlap          : %f\n", mean_ov1) ;
00313                 fprintf(fg, "Mean relative overlap : %f\n", mean_ovr1) ;
00314 
00315                 /* Write the 4th criteria statistics */
00316                 
00317                 fprintf(fg, "\n\t--\n") ;
00318                 fprintf(fg,   "\t-          _ 4th overlap criteria (alpha sphere overlap) _         -\n") ;
00319                 fprintf(fg,   "\t--\n\n") ;
00320                 fprintf(fg,"           :");
00321                 for( j = 0 ; j < novlp2 ; j++) {
00322                                 fprintf(fg, "  >%5.2f  :", ovlp2[j]) ;
00323                 }
00324                 fprintf(fg,"\n");
00325                 fprintf(fg,"----");
00326                 for( j = 0 ; j < novlp2 ; j++) {
00327                         fprintf(fg, "-") ;
00328                 }
00329                 fprintf(fg,"\n");
00330 
00331                 for(i = 0 ; i < nranks ; i++) {
00332                         fprintf(fg, "Rank <= %2d :", ranks[i]) ;
00333                         for( j = 0 ; j < novlp2 ; j++) {
00334                                 nok = 0 ;
00335                                 for(k = 0 ; k < par->nfiles ; k++) {
00336                                         if( status[k] == M_OK &&
00337                                                 ddata[k][M_CRIT4] >= ovlp2[j] &&
00338                                                 idata[k][M_POS4] <= ranks[i] &&
00339                                                 idata[k][M_POS4] > 0) {
00340                                                 nok ++ ;
00341                                         }
00342                                 }
00343                                 fprintf(fg, "  %6.2f  :", ((float)nok) / ((float) N)) ;
00344                         }
00345                         fprintf(fg, "\n") ;
00346                 }
00347                 fprintf(fg, "-\n") ;
00348                 fprintf(fg, "Mean overlap          : %f\n", mean_ov4) ;
00349                 fprintf(fg, "Mean relative overlap : %f\n", mean_ovr4) ;
00350 
00351                 /* Write the 4th criteria statistics */
00352                 
00353                 fprintf(fg, "\n\t--\n") ;
00354                 fprintf(fg,   "\t-          _ 5th overlap criteria (alpha sphere overlap) _         -\n") ;
00355                 fprintf(fg,   "\t--\n\n") ;
00356                 fprintf(fg,"           :");
00357                 for( j = 0 ; j < novlp3 ; j++) {
00358                                 fprintf(fg, "  >%5.2f  :", ovlp3[j]) ;
00359                 }
00360                 fprintf(fg,"\n");
00361                 fprintf(fg,"----");
00362                 for( j = 0 ; j < novlp3 ; j++) {
00363                         fprintf(fg, "-") ;
00364                 }
00365                 fprintf(fg,"\n");
00366 
00367                 for(i = 0 ; i < nranks ; i++) {
00368                         fprintf(fg, "Rank <= %2d :", ranks[i]) ;
00369                         for( j = 0 ; j < novlp3 ; j++) {
00370                                 nok = 0 ;
00371                                 for(k = 0 ; k < par->nfiles ; k++) {
00372                                         if( status[k] == M_OK &&
00373                                                 ddata[k][M_CRIT5] >= ovlp3[j] &&
00374                                                 idata[k][M_POS5] <= ranks[i] &&
00375                                                 idata[k][M_POS5] > 0) {
00376                                                 nok ++ ;
00377                                         }
00378                                 }
00379                                 fprintf(fg, "  %6.2f  :", ((float)nok) / ((float) N)) ;
00380                         }
00381                         fprintf(fg, "\n") ;
00382                 }
00383                 fprintf(fg, "-\n") ;
00384                 fprintf(fg, "Mean overlap          : %f\n", mean_ov5) ;
00385                 fprintf(fg, "Mean relative overlap : %f\n", mean_ovr5) ;
00386 
00387 
00388                 fprintf(fg, "\n\t--\n") ;
00389                 fprintf(fg,   "\t-      _ Concensus overlap criteria (alpha sphere overlap) _       -\n") ;
00390                 fprintf(fg,   "\t--\n\n") ;
00391                 fprintf(fg, "   Ratio of good predictions (dist = 3A) \n") ;
00392                 fprintf(fg, "\n") ;
00393 
00394                 for(i = 0 ; i < nranks ; i++) {
00395                         nok = 0 ;
00396                         for(k = 0 ; k < par->nfiles ; k++) {
00397                                 if(status[k] == M_OK && idata[k][M_POS6] <= ranks[i]
00398                                 && idata[k][M_POS6] > 0) {
00399                                         nok ++ ;
00400                                 }
00401                         }
00402                         //printf("%d NOK: %d, %d, %f\n", ranks[i], nok, N, ((float)nok) / ((float) N)) ;
00403                         fprintf(fg, "Rank <= %2d  :\t\t%6.2f\n", ranks[i],
00404                                         ((float)nok) / ((float) N)) ;
00405                 }
00406 
00407                 fprintf(fg, "-\n") ;
00408                 fprintf(fg, "Mean relative overlap           : %7.2f\n", mean_ovr6) ;
00409                 fprintf(fg, "Mean pocket volume (estimation) : %10.2f\n", mean_pvol6) ;
00410                 fprintf(fg, "Mean number of pocket atom      : %4d\n", mean_nbatm6) ;
00411 
00412                 fclose(fg) ;
00413         }
00414         else fprintf(stdout, "The file %s could not be opened\n", par->g_output) ;
00415 }
00416 
00417 /**
00418    ## FUNCTION: 
00419         void test_set(s_tparams *par)
00420   
00421    ## SPECIFICATION: 
00422         Test fpocket for a single set of apo/complexe/ligand.
00423   
00424    ## PARAMETRES:
00425         @ s_tparams *par  : Parameters, contain the fpocket parameters, and the list
00426                                             of data set to test. 
00427     @ int i           : ID if the current test set. Used to fill the two following
00428                                             statistic arrays
00429     @ float ddata[][] : Used to store statistics for the evaluation (float)
00430     @ float idata[][] : Used to store statistics for the evaluation (integers)
00431   
00432    ## RETURN:
00433         int - A flag saying whether or not the evaluation is successfull. Values:
00434                 M_LIGNOTFOUND -2
00435                 M_PDBOPENFAILED -1
00436                 M_OK 0 
00437                 M_NOPOCKETFOUND 1
00438   
00439 */
00440 int test_set(s_tparams *par, int i, float ddata [][M_NDDATA], int idata [][M_NIDATA]) 
00441 {
00442         s_atm **accpck = NULL,
00443                   **accpck2 = NULL ;
00444         
00445         char *fa = par->fapo[i], 
00446                  *fc = par->fcomplex[i],
00447                  *lig = par->fligan[i] ;
00448         
00449         int naccpck = 0,
00450                 naccpck2 = 0 ;
00451 
00452         /** Load the apo form, the complexe form with ligand, and the complexe form 
00453           * without ligand. The complexe without ligand will allow us to get exact 
00454           * atoms contacted by voronoi vertices detected by fpocket on the complexe 
00455           * form and then detect which atoms are involved in the interaction.
00456          **/
00457         s_pdb *apdb = rpdb_open(fa, lig, M_DONT_KEEP_LIG) ;
00458         s_pdb *cpdb = rpdb_open(fc, lig, M_KEEP_LIG);
00459         s_pdb *cpdb_nolig = rpdb_open(fc, lig, M_DONT_KEEP_LIG) ;
00460         
00461         /* Check that both pdb files are opened */
00462         if(! apdb || !cpdb) {
00463                 fprintf(stderr, "!! PDB loading failed for %s-%s ligand %s... %p %p\n", 
00464                                                 fc, fa, lig, apdb, cpdb) ;
00465                 if(cpdb) free_pdb_atoms(cpdb) ;
00466                 if(apdb) free_pdb_atoms(apdb) ;
00467                 
00468                 return M_LIGNOTFOUND ;
00469         }
00470         
00471         /* Check now if the ligand have been found */
00472         if(cpdb->natm_lig <= 0) {
00473                 
00474                 fprintf(stderr, "! Ligand '%s' not found in the complex %s...\n",
00475                                                 lig, fc) ;
00476                 
00477                 return M_LIGNOTFOUND ;
00478         }
00479         
00480         
00481         /* PDBs and ligand OK, now read PDB and launch fpocket */
00482         
00483         /* Get rid of HETATM in the apo structure*/
00484         rpdb_read(apdb, lig, M_DONT_KEEP_LIG) ;
00485         rpdb_read(cpdb, lig, M_KEEP_LIG) ;
00486         rpdb_read(cpdb_nolig, lig, M_DONT_KEEP_LIG) ;
00487 
00488         c_lst_pockets *pockets = search_pocket(apdb, par->fpar,cpdb);
00489         set_pockets_bary(pockets) ;
00490         
00491         /* Check that at least one pocket have been found */
00492         if(!pockets || pockets->n_pockets <= 0) {
00493                 if(pockets) c_lst_pocket_free(pockets) ;
00494                 return M_NOPOCKETFOUND ;
00495         }
00496         
00497         /* Everything is OK, so now perform the evaluation*/
00498         
00499         /* We found pocket on the apo form get volume and mass of the ligand */
00500         idata[i][M_NPOCKET] = pockets->n_pockets ;
00501         ddata[i][M_LIGMASS] = get_mol_mass_ptr(cpdb->latm_lig, cpdb->natm_lig) ;
00502         ddata[i][M_LIGVOL]  = get_mol_volume_ptr(cpdb->latm_lig, cpdb->natm_lig, 
00503                                                                                          par->fpar->nb_mcv_iter);
00504 
00505         /* Get atoms involved in the actual pocket */
00506         accpck = get_actual_pocket_DEPRECATED(cpdb, M_CRIT2_D, &naccpck) ;
00507         accpck2 = get_actual_pocket(cpdb, cpdb_nolig, i, par, &naccpck2) ;
00508 
00509         fflush(stdout) ;
00510         if (naccpck2 <= 0) {
00511                 fprintf(stdout, "! Warning: actual pocket has 0 atoms!! %s %d\n", fa,
00512                                                 naccpck2) ;
00513         }
00514 
00515         /* Calculate evaluation criterias */
00516         check_pockets(pockets, accpck2, naccpck2, cpdb->latm_lig, cpdb->natm_lig, 
00517                                   accpck, naccpck, ddata, idata, i) ;
00518 
00519         if(par->keep_fpout != 0) {
00520                 write_out_fpocket(pockets, apdb, par->fapo[i]) ;
00521         }
00522         /* write_pdb_com(cpdb, par->fcomplex[i]) ; */
00523         
00524         /* Free memory */
00525         c_lst_pocket_free(pockets) ;
00526         my_free(accpck) ;
00527         free_pdb_atoms(apdb) ;
00528         free_pdb_atoms(cpdb) ;
00529         free_pdb_atoms(cpdb_nolig) ;
00530                 
00531         return M_OK ;
00532 }
00533 
00534 /**
00535    ## FUNCTION: 
00536         get_actual_pocket
00537   
00538    ## SPECIFICATION: 
00539         Get atoms contacted by voronoi vertices present in the actual pocket and 
00540         situated at a distance of M_DST_CRIT from each ligand atoms.
00541         This way of finding exact atoms contacted by the ligand is more accurate,
00542         but need an extra search by fpocket to find voronoi vertices in the complexe
00543         form without the ligand...
00544   
00545    ## PARAMETRES:
00546         @ s_pdb *cpdb       : The pdb structure of the complexe.
00547         @ s_pdb *cpdb_nolig : The pdb structure of the complexe WITHOUT the ligand
00548         @ int i             : ID of the current test set
00549         @ s_tparams *par    : Parameters
00550         @ int *nb_atm       : Number of atom (this an output, it should be set here)
00551   
00552    ## RETURN:
00553         s_atm **: A 2D table containing pointers to atoms which defines the pocket.
00554   
00555 */
00556 s_atm** get_actual_pocket(s_pdb *cpdb, s_pdb *cpdb_nolig, int i, s_tparams *par, int *nb_atm) 
00557 {
00558         s_atm **neigh = NULL ;
00559 
00560         c_lst_pockets *pockets = search_pocket(cpdb_nolig, par->fpar,cpdb);
00561         if(pockets && pockets->n_pockets > 0) {
00562                 /* */
00563                 /* Get the list of atoms contacted by the vertices near the ligand. */
00564                 
00565                 neigh = get_mol_ctd_atm_neigh(cpdb->latm_lig, cpdb->natm_lig,
00566                                                                           pockets->vertices->pvertices,
00567                                                                           pockets->vertices->nvert, M_CRIT1_D,
00568                                                                           M_INTERFACE_SEARCH, nb_atm) ;
00569 
00570                 c_lst_pocket_free(pockets) ;
00571         }
00572         else {
00573                 fprintf(stderr, "! No pocket found for apo %s and complex %s ligand %s...\n", 
00574                                                 par->fapo[i], par->fcomplex[i], par->fligan[i]) ;
00575 
00576                 if(pockets) c_lst_pocket_free(pockets) ;
00577         }
00578 
00579         return neigh ;
00580 }
00581 
00582 /**
00583    ## FUNCTION: 
00584         get_actual_pocket
00585   
00586    ## SPECIFICATION: 
00587         Get atoms situated at a distance of lig_dist_crit from each ligand atoms.
00588     This is a quite inacurate way of defining the pocket -> DEPRECATED
00589     If the distance criteria used does not allow any identification, we will
00590     test a larger distance until 6A. 
00591   
00592    ## PARAMETRES:
00593         @ s_pdb *cpdb         : The pdb structure of the complexe.
00594         @ float lig_dist_crit : The distance criteria
00595         @ int *nb_atm         : OUTPUT: Number of atom found in the pocket
00596   
00597    ## RETURN:
00598         s_atm **: A 2D table containing pointers to atoms which defines the pocket.
00599   
00600 */
00601 s_atm** get_actual_pocket_DEPRECATED(s_pdb *cpdb, float lig_dist_crit, int *nb_atm) 
00602 {
00603         /* Getting the ligan's neighbors. */
00604         s_atm **alneigh = get_mol_atm_neigh(cpdb->latm_lig, cpdb->natm_lig,
00605                                                                                 cpdb->latoms_p, cpdb->natoms,
00606                                                                                 lig_dist_crit, nb_atm) ;
00607         if(*nb_atm <= 0) {
00608                 if(alneigh) my_free(alneigh) ;
00609                 alneigh = get_mol_atm_neigh(cpdb->latm_lig, cpdb->natm_lig,
00610                                                                         cpdb->latoms_p, cpdb->natoms,
00611                                                                         lig_dist_crit+1.0, nb_atm) ;
00612 
00613                 if(*nb_atm <= 0) {
00614                         if(alneigh) my_free(alneigh) ;
00615                         alneigh = get_mol_atm_neigh(cpdb->latm_lig, cpdb->natm_lig,
00616                                                                                 cpdb->latoms_p, cpdb->natoms,
00617                                                                                 lig_dist_crit+2.0, nb_atm) ;
00618                         if(*nb_atm <= 0) {
00619                                 if(alneigh) my_free(alneigh) ;
00620                                 alneigh = NULL ;
00621                         }
00622                 }
00623         }
00624         
00625         return alneigh ;
00626 }
00627 
00628 /**
00629    ## FUNCTION: 
00630         check_pockets
00631   
00632    ## SPECIFICATION: 
00633         Checking the pocket prediction for one data set. (apo+complex+ligand)
00634         This function will calculate three criterias:
00635                 1 - Overlap using the first definition of the pocket (atoms contacted
00636                         by vertices next to the ligand)
00637                 2 - Overlap using the second definition of the pocket (atoms situated
00638                         at N A of the ligand)
00639                 3 - Smallest distance of each ligand atom from the barycenter of the
00640                         pocket. Cutoff: 4A
00641   
00642    The two arrays ddata and idata given in paramet are used as output.
00643   
00644    ## PARAMETRES:
00645         @ c_lst_pockets *pockets: The pockets found for this set.
00646     @ s_atm accpck          : The actual pocket using first definition
00647     @ int naccpckk          : Number of atoms in the actual pocket definition 1
00648     @ s_atm lig             : Ligand atoms
00649     @ int nalig             : Number of ligand atoms
00650     @ s_atm alneigh         : The actual pocket using second definition
00651     @ int nlneigh           : Number of atoms in the actual pocket definition 2
00652     @ float ddata[][]       : Used to store statistics for the evaluation (float)
00653     @ float idata[][]       : Used to store statistics for the evaluation (integers)
00654     @ int i                 : ID if the current test set. Used to fill the two 
00655                                                           following statistic arrays
00656   
00657    ## RETURN:
00658         void
00659   
00660 */
00661 void check_pockets(c_lst_pockets *pockets, s_atm **accpck, int naccpck, s_atm **lig, 
00662                                    int nalig, s_atm **alneigh, int nlneigh, 
00663                                    float ddata [][M_NDDATA], int idata [][M_NIDATA], int i)
00664 {
00665         int found [] = {0, 0, 0, 0, 0, 0} ;
00666         int npneigh = 0, pos, j ;
00667         float ov1, ov2, ov3, ov4, dst = 0.0 ;
00668         /* float ovol ; */
00669         
00670         s_atm **pneigh = NULL ;
00671         s_vvertice **pvert = NULL ;
00672         
00673         node_pocket *ncur = NULL ;
00674         s_pocket *pcur = NULL ;
00675         
00676         /* Check the correspondance for each pocket found */
00677         ncur = pockets->first ;
00678 
00679         int n_lig_molecules=1;
00680         char chain_tmp[2];
00681         int resnumber_tmp;
00682         strcpy(chain_tmp,lig[0]->chain);
00683         resnumber_tmp = lig[0]->res_id;
00684         
00685         for (j = 1 ; j < nalig ; j++) {
00686             if(strcmp(chain_tmp,lig[j]->chain) !=0 || resnumber_tmp!=lig[j]->res_id){
00687                 n_lig_molecules++;
00688                 strcpy(chain_tmp,lig[j]->chain);
00689                 resnumber_tmp =lig[j]->res_id;
00690             }
00691         }
00692         
00693         pos = 0 ;
00694         while(ncur) {
00695                 pos ++ ;
00696                 pcur = ncur->pocket ;
00697 
00698                 /*
00699                         TEST THE FIRST CRITERIA (DISTANCE TO BARYCENTER)
00700                  */
00701                 if(! found[0]) {
00702                         for (j = 0 ; j < nalig ; j++) {
00703                                 dst = dist(lig[j]->x, lig[j]->y, lig[j]->z, 
00704                                                         pcur->bary[0], pcur->bary[1], pcur->bary[2]) ;
00705                                 if(dst < M_CRIT3_VAL) {
00706                                 /* Criteria 3 OK  */
00707                                         ddata[i][M_MINDST] = dst ;
00708                                         ddata[i][M_OREL3] = (naccpck <= 0.0)?-1.0 :(float)npneigh/(float)naccpck*100.0 ;
00709                                         ddata[i][M_POCKETVOL_C3] = pcur->pdesc->volume ;
00710                                         idata[i][M_POS3] = pos ;
00711                                         idata[i][M_NATM3] = count_pocket_contacted_atms(pcur) ;
00712                                         found[0] = 1 ; break ;
00713                                 }
00714                         }
00715                 }
00716 
00717                 /*
00718                         TEST THE 2nd and 3rd CRITERIA (ATOMIC OVERLAP)
00719                  */
00720                 if(!found[1] || !found[2]) {
00721                         pneigh = get_pocket_contacted_atms(pcur, &npneigh) ;
00722                         if(! found[1]) {
00723                                 ov1 = atm_corsp(alneigh, nlneigh, pneigh, npneigh) ;
00724                                 if(ov1 > M_CRIT1_VAL) {
00725                                         idata[i][M_POS1] = pos ;
00726                                         ddata[i][M_MAXPCT1] = ov1 ;
00727                                         ddata[i][M_OREL1] = (nlneigh <= 0.0)? -1.0 : (float)npneigh/(float)nlneigh*100.0 ;
00728                                         found[1] = 1 ;
00729                                 }
00730                         }
00731 
00732                         if(! found[2]) {
00733                                 ov2 = atm_corsp(accpck, naccpck, pneigh, npneigh) ;
00734                                 if(ov2 > M_CRIT2_VAL) {
00735                                         idata[i][M_POS2] = pos ;
00736                                         ddata[i][M_MAXPCT2] = ov2 ;
00737                                         ddata[i][M_OREL2] = (naccpck <= 0.0)?-1.0 : (float)npneigh/(float)naccpck*100.0 ;
00738                                         found[2] = 1 ;
00739                                 }
00740                         }
00741                         my_free(pneigh) ;
00742                 }
00743 
00744                 /*
00745                         TEST THE 4th and 5th CRITERIA (VERTICE OVERLAPS)
00746                  */
00747                 if(!found[3] || !found[4] || !found[5]) {
00748                         pvert = get_pocket_pvertices(pcur) ;
00749                         if(!found[3]){
00750                                 /* Calculate proportion of ligand atom that lies within 3A
00751                                    of at least one vertice */
00752                             
00753                                 ov3 = count_atm_prop_vert_neigh( lig, nalig,pvert, pcur->size, M_CRIT4_D,n_lig_molecules) ;
00754                                 if(ov3 > M_CRIT4_VAL) {
00755                                         idata[i][M_POS4] = pos ;
00756                                         ddata[i][M_CRIT4] = ov3 ;
00757                                         ddata[i][M_OREL4] = ov3*100 ;
00758                                         found[3] = 1 ;
00759                                 }
00760                         }
00761 
00762                         if(!found[4]){
00763                                 /* Calculate proportion of pocket vertice that lies within 3A
00764                                    of at least one ligand atom */
00765                                 ov4 = count_pocket_lig_vert_ovlp(lig, nalig,
00766                                                                                                  pvert, pcur->size, M_CRIT5_D) ;
00767                                 if(ov4 > M_CRIT5_VAL) {
00768                                         idata[i][M_POS5] = pos ;
00769                                         ddata[i][M_CRIT5] = ov4 ;
00770                                         ddata[i][M_OREL5] = ov4*100 ;
00771                                         found[4] = 1 ;
00772                                 }
00773                         }
00774 
00775                         if(!found[5]){
00776                                 ov3 = count_atm_prop_vert_neigh( lig, nalig,
00777                                                                                          pvert, pcur->size, M_CRIT4_D,n_lig_molecules );
00778                                 ov4 = count_pocket_lig_vert_ovlp(lig, nalig,
00779                                                                                                  pvert, pcur->size, M_CRIT5_D) ;
00780                                 if(ov4 > M_CRIT5_VAL && ov3 > M_CRIT4_VAL) {
00781                                         idata[i][M_POS6] = pos ;
00782                                         ddata[i][M_CRIT6] = 1.0 ;
00783                                         ddata[i][M_POCKETVOL_C6] = pcur->pdesc->volume ;
00784                                         ddata[i][M_OREL6] = ov4*100 ;
00785                                         idata[i][M_NATM6] = count_pocket_contacted_atms(pcur) ;
00786                                         found[5] = 1 ;
00787                                 }
00788                         }
00789                         my_free(pvert) ;
00790                 }
00791                 
00792                 ncur = ncur->next ;
00793 
00794                 /* Break the loop if all criteria are OK */
00795                 if(found[0] && found[1] && found[2] && found[3] && found[4] && found[5]) break ;
00796                 
00797         }
00798 
00799         if (! found[0]) {
00800                 ddata[i][M_MINDST] = 0.0 ;
00801                 ddata[i][M_OREL3] = 0.0 ;
00802                 idata[i][M_POS3] = 0 ;
00803                 idata[i][M_NATM3] = -1 ;
00804                 ddata[i][M_POCKETVOL_C3] = 0.0 ;
00805         }
00806 
00807         if (! found[1]) {
00808                 ddata[i][M_MAXPCT1] = 0.0 ;
00809                 ddata[i][M_OREL1] = 0.0 ;
00810                 idata[i][M_POS1] = 0 ;
00811         }
00812 
00813         if (! found[2]) {
00814                 ddata[i][M_MAXPCT2] = 0.0 ;
00815                 ddata[i][M_OREL2] = 0.0 ;
00816                 idata[i][M_POS2] = 0 ;
00817         }
00818 
00819         if (! found[3]) {
00820                 ddata[i][M_CRIT4] = 0.0 ;
00821                 ddata[i][M_OREL4] = 0.0 ;
00822                 idata[i][M_POS4] = 0 ;
00823         }
00824 
00825         if (! found[4]) {
00826                 ddata[i][M_CRIT5] = 0.0 ;
00827                 ddata[i][M_OREL5] = 0.0 ;
00828                 idata[i][M_POS5] = 0 ;
00829         }
00830         
00831         if (! found[5]) {
00832                 ddata[i][M_CRIT6] = 0.0 ;
00833                 ddata[i][M_OREL6] = 0.0 ;
00834                 idata[i][M_POS6] = 0 ;
00835                 idata[i][M_NATM6] = -1 ;
00836                 ddata[i][M_POCKETVOL_C6] = 0.0 ;
00837         }
00838 }
00839 
00840 /**
00841         --
00842         --
00843 
00844         ALL THE FOPLLOWING FUNCTIONS SHOULD BE USED TO CALCULATE AN OVERLAP VOLUME
00845         BETWEEN A GIVEN POPCKET AND THE REAL POCKET.
00846 
00847         THE FUNCTION set_overlap_volumes WILL SET AN OVERLAP VOLUME FOR EACH POCKET
00848         IN THE LIST GIVEN IN ARGUMENT, USING MONTE-CARLO OR BASIC DISCRETE ALGORITHM 
00849         DEPENDING ON THE PARAMETERS OF THE PROGRAM.
00850 
00851         BOTH ALGORITHM CAN BE USED INDEPENDENLTLY WITHOUT THE NEED OF A PARAMETER
00852         STRUCTURE (s_fparams)
00853 
00854         --
00855         --
00856    
00857 */
00858 
00859 /**
00860    ## FUNCTION: 
00861         void set_overlap_volumes
00862   
00863    ## SPECIFICATION: 
00864         This function will set an overlap volume between a pocket, defined here by
00865     its voronoi vertices, and the ligand.
00866         Depending on the parameters, the method used will be either Monte Carlo 
00867         approximation, or discret basic approximation.
00868   
00869    ## PARAMETRES:
00870         @ s_pocket *pockets : The pocket
00871         @ s_atm **lig       : The ligand atoms.
00872         @ int natoms natoms : Number of atoms of the ligand.
00873         @ float lig_vol     : Volume of the ligand
00874         @ s_fparams *params : Parameters, giving the algorithm to use and the 
00875                                               corresponding parameter.
00876   
00877    ## RETURN:
00878         float: The overlap volume found
00879   
00880 */
00881 float set_overlap_volumes(s_pocket *pocket, s_atm **lig, int natoms, float lig_vol, 
00882                                                   s_fparams *params) 
00883 {
00884         /* The function to use for calculation */
00885         float (*pf_set_vol)(s_atm**, int, float, s_pocket*, int ) = NULL ;      
00886         int crit = 0,
00887                 method = 0 ;
00888 
00889         if(params->basic_volume_div <= 0) {
00890                 pf_set_vol = set_mc_overlap_volume ;
00891                 crit = params->nb_mcv_iter ;
00892         }
00893         else {
00894                 pf_set_vol = set_basic_overlap_volume ;
00895                 crit = params->basic_volume_div ;
00896                 method = 1 ;
00897         }
00898 
00899         if(pocket->v_lst->n_vertices > 150) {
00900                 if(method == 0 && crit < 20000) crit = 20000 ;
00901         }
00902         
00903         return pf_set_vol(lig, natoms, lig_vol, pocket, crit) ;
00904 }
00905 
00906 /**
00907    ## FUNCTION:
00908         set_mc_overlap_volume
00909   
00910    ## SPECIFICATION: 
00911         This function will set an overlap volume between each pocket and the ligand.
00912         Method: monte Carlo aproximation.
00913   
00914    ## PARAMETRES:
00915         @ s_atm **lig       : The ligand atoms.
00916         @ int natoms natoms : Number of atoms of the ligand.
00917         @ float lig_vol     : Volume of the ligand
00918         @ s_pocket *pockets : The pocket
00919         @ int niter         : Number of monte-carlo iterations
00920   
00921    ## RETURN:
00922   
00923 */
00924 float set_mc_overlap_volume(s_atm **lig, int natoms, float lig_vol, 
00925                                                         s_pocket *pocket, int niter)
00926 {
00927         c_lst_vertices *vertices = pocket->v_lst ;
00928         s_vvertice *vcur = NULL ;
00929 
00930         int i = 0, j = 0,
00931                 nb_in = 0,
00932                 found = 0 ;
00933 
00934         float xmin = 0.0, xmax = 0.0,
00935                   ymin = 0.0, ymax = 0.0,
00936                   zmin = 0.0, zmax = 0.0,
00937                   xtmp = 0.0, ytmp = 0.0, ztmp = 0.0,
00938                   xr   = 0.0, yr   = 0.0, zr   = 0.0,
00939                   vbox = 0.0, full_vol = 0.0;
00940 
00941         s_atm *acur = NULL ;
00942 
00943         /* First, search extrems coordinates in the ligan */
00944         for(i = 0 ; i < natoms ; i++) {
00945                 acur = lig[i] ;
00946                 if(i == 0) {
00947                         xmin = acur->x - acur->radius ; xmax = acur->x + acur->radius ;
00948                         ymin = acur->y - acur->radius ; ymax = acur->y + acur->radius ;
00949                         zmin = acur->z - acur->radius ; zmax = acur->z + acur->radius ;
00950                 }
00951                 else {
00952                 /* Update the minimum and maximum extreme point */
00953                         if(xmin > (xtmp = acur->x - acur->radius)) xmin = xtmp ;
00954                         else if(xmax < (xtmp = acur->x + acur->radius)) xmax = xtmp ;
00955         
00956                         if(ymin > (ytmp = acur->y - acur->radius)) ymin = ytmp ;
00957                         else if(ymax < (ytmp = acur->y + acur->radius)) ymax = ytmp ;
00958         
00959                         if(zmin > (ztmp = acur->z - acur->radius)) zmin = ztmp ;
00960                         else if(zmax < (ztmp = acur->z + acur->radius)) zmax = ztmp ;
00961                 }
00962         }
00963 
00964         node_vertice *cur = vertices->first ;
00965         /* Continue with the pocket */
00966         while(cur) {
00967                 vcur = cur->vertice ;
00968                 if(xmin > (xtmp = vcur->x - vcur->ray)) xmin = xtmp ;
00969                 else if(xmax < (xtmp = vcur->x + vcur->ray)) xmax = xtmp ;
00970         
00971                 if(ymin > (ytmp = vcur->y - vcur->ray)) ymin = ytmp ;
00972                 else if(ymax < (ytmp = vcur->y + vcur->ray)) ymax = ytmp ;
00973         
00974                 if(zmin > (ztmp = vcur->z - vcur->ray)) zmin = ztmp ;
00975                 else if(zmax < (ztmp = vcur->z + vcur->ray)) zmax = ztmp ;
00976 
00977                 cur = cur->next ;
00978         }
00979 
00980 
00981         /* Next calculate the box volume */
00982         vbox = (xmax - xmin)*(ymax - ymin)*(zmax - zmin) ;
00983 
00984         /* Then apply monte carlo approximation of the volume */
00985         for(i = 0 ; i < niter ; i++) {
00986                 found = 0 ;
00987                 xr = rand_uniform(xmin, xmax) ;
00988                 yr = rand_uniform(ymin, ymax) ;
00989                 zr = rand_uniform(zmin, zmax) ;
00990 
00991                 cur = vertices->first ;
00992                 while(cur) {
00993                         vcur = cur->vertice ;
00994                         xtmp = vcur->x - xr ;
00995                         ytmp = vcur->y - yr ;
00996                         ztmp = vcur->z - zr ;
00997 
00998                 /* Compare r^2 and dist(center, random_point)^2 */
00999                         if((vcur->ray*vcur->ray) > (xtmp*xtmp + ytmp*ytmp + ztmp*ztmp)) {
01000                         /* The point is inside one of the vertice!!*/ 
01001                                 nb_in ++ ; found = 1 ; break ;
01002                         }
01003                         cur = cur->next ;
01004                 }
01005 
01006                 if(!found) {
01007                         for(j = 0 ; j < natoms ; j++) {
01008                                 acur = lig[j] ;
01009                                 xtmp = acur->x - xr ;
01010                                 ytmp = acur->y - yr ;
01011                                 ztmp = acur->z - zr ;
01012         
01013                         /* Compare r^2 and dist(center, random_point)^2 */
01014                                 if((acur->radius*acur->radius) > (xtmp*xtmp + ytmp*ytmp + ztmp*ztmp)) {
01015                                 /* The point is inside one of the vertice!! */
01016                                         nb_in ++ ; break ;
01017                                 }
01018                         }
01019                 }
01020         }
01021 
01022         full_vol = ((float)nb_in)/((float)niter)*vbox ; 
01023         pocket->vol_corresp = ((pocket->pdesc->volume + lig_vol) - full_vol)/lig_vol ;
01024 
01025         return pocket->vol_corresp ;
01026 }
01027 
01028 /**
01029    ## FUNCTION: 
01030         set_basic_overlap_volume
01031   
01032    ## SPECIFICATION: 
01033         This function will set an overlap volume between each pocket and the ligand.
01034         Discrete aproximation. (algo in idiscret*idiscret*idiscret operation minimum)
01035   
01036    ## PARAMETRES:
01037         @ s_atm **lig       : The ligand atoms.
01038         @ int natoms natoms : Number of atoms of the ligand.
01039         @ float lig_vol     : Volume of the ligand
01040         @ s_pocket *pockets : The pocket
01041         @ int idiscret      : Discretisation of the space: splitt the space in N*N*N
01042                                                   cubes.
01043   
01044    ## RETURN:
01045   
01046 */
01047 float set_basic_overlap_volume(s_atm **lig, int natoms, float lig_vol,
01048                                                            s_pocket *pocket, int idiscret)
01049 {
01050         c_lst_vertices *vertices = pocket->v_lst ;
01051         s_vvertice *vcur = NULL ;
01052 
01053         int i = 0, j = 0,
01054                 nb_in = 0,
01055                 found = 0,
01056                 niter = 0 ;
01057 
01058         float discret = 1.0/(float)idiscret ;
01059 
01060         float x = 0.0, y = 0.0, z = 0.0,
01061                   xstep = 0.0, ystep = 0.0, zstep = 0.0 ;
01062 
01063         float xmin = 0.0, xmax = 0.0,
01064                   ymin = 0.0, ymax = 0.0,
01065                   zmin = 0.0, zmax = 0.0,
01066                   xtmp = 0.0, ytmp = 0.0, ztmp = 0.0,
01067                   vbox = 0.0, full_vol = 0.0;
01068 
01069         s_atm *acur = NULL ;
01070 
01071         /* First, search extrems coordinates in the ligan */
01072         for(i = 0 ; i < natoms ; i++) {
01073                 acur = lig[i] ;
01074                 if(i == 0) {
01075                         xmin = acur->x - acur->radius ; xmax = acur->x + acur->radius ;
01076                         ymin = acur->y - acur->radius ; ymax = acur->y + acur->radius ;
01077                         zmin = acur->z - acur->radius ; zmax = acur->z + acur->radius ;
01078                 }
01079                 else {
01080                 /* Update the minimum and maximum extreme point */
01081                         if(xmin > (xtmp = acur->x - acur->radius)) xmin = xtmp ;
01082                         else if(xmax < (xtmp = acur->x + acur->radius)) xmax = xtmp ;
01083         
01084                         if(ymin > (ytmp = acur->y - acur->radius)) ymin = ytmp ;
01085                         else if(ymax < (ytmp = acur->y + acur->radius)) ymax = ytmp ;
01086         
01087                         if(zmin > (ztmp = acur->z - acur->radius)) zmin = ztmp ;
01088                         else if(zmax < (ztmp = acur->z + acur->radius)) zmax = ztmp ;
01089                 }
01090         }
01091 
01092         node_vertice *cur = vertices->first ;
01093         /* Continue with the pocket */
01094         while(cur) {
01095                 vcur = cur->vertice ;
01096                 if(xmin > (xtmp = vcur->x - vcur->ray)) xmin = xtmp ;
01097                 else if(xmax < (xtmp = vcur->x + vcur->ray)) xmax = xtmp ;
01098         
01099                 if(ymin > (ytmp = vcur->y - vcur->ray)) ymin = ytmp ;
01100                 else if(ymax < (ytmp = vcur->y + vcur->ray)) ymax = ytmp ;
01101         
01102                 if(zmin > (ztmp = vcur->z - vcur->ray)) zmin = ztmp ;
01103                 else if(zmax < (ztmp = vcur->z + vcur->ray)) zmax = ztmp ;
01104 
01105                 cur = cur->next ;
01106         }
01107 
01108 
01109         /* Next calculate the box volume */
01110         vbox = (xmax - xmin)*(ymax - ymin)*(zmax - zmin) ;
01111 
01112         xstep = discret * (xmax - xmin) ;
01113         ystep = discret * (ymax - ymin) ;
01114         zstep = discret * (zmax - zmin) ;
01115 
01116         /* Then apply monte carlo approximation of the volume */
01117         for(x = xmin ; x < xmax ; x += xstep) {
01118                 for(y = ymin ; y < ymax ; y += ystep) { 
01119                         for(z = zmin ; z < zmax ; z += zstep) {
01120                                 found = 0 ;
01121                                 cur = vertices->first ;
01122                                 while(cur) {
01123                                         vcur = cur->vertice ;
01124                                         xtmp = vcur->x - x ;
01125                                         ytmp = vcur->y - y ;
01126                                         ztmp = vcur->z - z ;
01127                 
01128                                 /* Compare r^2 and dist(center, random_point)^2 */
01129                                         if((vcur->ray*vcur->ray) > (xtmp*xtmp + ytmp*ytmp + ztmp*ztmp)) {
01130                                         /* The point is inside one of the vertice!! */
01131                                                 nb_in ++ ; found = 1 ; break ;
01132                                         }
01133                                         cur = cur->next ;
01134                                 }
01135                 
01136                                 if(!found) {
01137                                         for(j = 0 ; j < natoms ; j++) {
01138                                                 acur = lig[j] ;
01139                                                 xtmp = acur->x - x ;
01140                                                 ytmp = acur->y - y ;
01141                                                 ztmp = acur->z - z ;
01142                         
01143                                         /* Compare r^2 and dist(center, random_point)^2 */
01144                                                 if((acur->radius*acur->radius) > (xtmp*xtmp + ytmp*ytmp + ztmp*ztmp)) {
01145                                                 /* The point is inside one of the vertice!! */
01146                                                         nb_in ++ ; break ;
01147                                                 }
01148                                         }
01149                                 }
01150                                 niter++ ;
01151                         }
01152                 }
01153         }
01154 
01155         full_vol = ((float)nb_in)/((float)niter)*vbox ; 
01156         pocket->vol_corresp = ((pocket->pdesc->volume + lig_vol) - full_vol)/lig_vol ;
01157 
01158         /* Ok lets just return the volume Vpok = Nb_in/Niter*Vbox */
01159         return pocket->vol_corresp ;
01160 }