src/tools/gmx_anadock.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "confio.h"
  40 #include "copyrite.h"
  41 #include "futil.h"
  42 #include "gmx_fatal.h"
  43 #include "smalloc.h"
  44 #include "string2.h"
  45 #include "vec.h"
  46 #include "gmx_statistics.h"
  47 #include "statutil.h"
  48 #include "typedefs.h"
  49 #include "xvgr.h"
  50 #include "macros.h"
  51
  52 static const char *etitles[] = { "E-docked", "Free Energy" };
  53
  54 typedef struct {
  55     real     edocked, efree;
  56     int      index, cluster_id;
  57     t_atoms  atoms;
  58     rvec    *x;
  59     int      ePBC;
  60     matrix   box;
  61 } t_pdbfile;
  62
  63 static t_pdbfile *read_pdbf(const char *fn)
  64 {
  65     t_pdbfile *pdbf;
  66     double     e;
  67     char       buf[256], *ptr;
  68     int        natoms;
  69     FILE      *fp;
  70
  71     snew(pdbf, 1);
  72     get_stx_coordnum (fn, &natoms);
  73     init_t_atoms(&(pdbf->atoms), natoms, FALSE);
  74     snew(pdbf->x, natoms);
  75     read_stx_conf(fn, buf, &pdbf->atoms, pdbf->x, NULL, &pdbf->ePBC, pdbf->box);
  76     fp = ffopen(fn, "r");
  77     do
  78     {
  79         ptr = fgets2(buf, 255, fp);
  80         if (ptr)
  81         {
  82             if (strstr(buf, "Intermolecular") != NULL)
  83             {
  84                 ptr = strchr(buf, '=');
  85                 sscanf(ptr+1, "%lf", &e);
  86                 pdbf->edocked = e;
  87             }
  88             else if (strstr(buf, "Estimated Free") != NULL)
  89             {
  90                 ptr = strchr(buf, '=');
  91                 sscanf(ptr+1, "%lf", &e);
  92                 pdbf->efree = e;
  93             }
  94         }
  95     }
  96     while (ptr != NULL);
  97     ffclose(fp);
  98
  99     return pdbf;
 100 }
 101
 102 static t_pdbfile **read_em_all(const char *fn, int *npdbf)
 103 {
 104     t_pdbfile **pdbf = 0;
 105     int         i, maxpdbf;
 106     char        buf[256], name[256];
 107     gmx_bool    bExist;
 108
 109     strcpy(buf, fn);
 110     buf[strlen(buf)-4] = '\0';
 111     maxpdbf            = 100;
 112     snew(pdbf, maxpdbf);
 113     i = 0;
 114     do
 115     {
 116         sprintf(name, "%s_%d.pdb", buf, i+1);
 117         if ((bExist = gmx_fexist(name)) == TRUE)
 118         {
 119             pdbf[i]        = read_pdbf(name);
 120             pdbf[i]->index = i+1;
 121             i++;
 122             if (i >= maxpdbf)
 123             {
 124                 maxpdbf += 100;
 125                 srenew(pdbf, maxpdbf);
 126             }
 127         }
 128     }
 129     while (bExist);
 130
 131     *npdbf = i;
 132
 133     printf("Found %d pdb files\n", i);
 134
 135     return pdbf;
 136 }
 137
 138 static gmx_bool bFreeSort = FALSE;
 139
 140 static int pdbf_comp(const void *a, const void *b)
 141 {
 142     t_pdbfile *pa, *pb;
 143     real       x;
 144     int        dc;
 145
 146     pa = *(t_pdbfile **)a;
 147     pb = *(t_pdbfile **)b;
 148
 149     dc = pa->cluster_id - pb->cluster_id;
 150
 151     if (dc == 0)
 152     {
 153         if (bFreeSort)
 154         {
 155             x = pa->efree   - pb->efree;
 156         }
 157         else
 158         {
 159             x = pa->edocked - pb->edocked;
 160         }
 161
 162         if (x < 0)
 163         {
 164             return -1;
 165         }
 166         else if (x > 0)
 167         {
 168             return 1;
 169         }
 170         else
 171         {
 172             return 0;
 173         }
 174     }
 175     else
 176     {
 177         return dc;
 178     }
 179 }
 180
 181 static void analyse_em_all(int npdb, t_pdbfile *pdbf[], const char *edocked,
 182                            const char *efree, const output_env_t oenv)
 183 {
 184     FILE *fp;
 185     int   i;
 186
 187     for (bFreeSort = FALSE; (bFreeSort <= TRUE); bFreeSort++)
 188     {
 189         qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 190         fp = xvgropen(bFreeSort ? efree : edocked,
 191                       etitles[bFreeSort], "()", "E (kJ/mol)", oenv);
 192         for (i = 0; (i < npdb); i++)
 193         {
 194             fprintf(fp, "%12lf\n", bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
 195         }
 196         ffclose(fp);
 197     }
 198 }
 199
 200 static void clust_stat(FILE *fp, int start, int end, t_pdbfile *pdbf[])
 201 {
 202     int         i;
 203     gmx_stats_t ed, ef;
 204     real        aver, sigma;
 205
 206     ed = gmx_stats_init();
 207     ef = gmx_stats_init();
 208     for (i = start; (i < end); i++)
 209     {
 210         gmx_stats_add_point(ed, i-start, pdbf[i]->edocked, 0, 0);
 211         gmx_stats_add_point(ef, i-start, pdbf[i]->efree, 0, 0);
 212     }
 213     gmx_stats_get_ase(ed, &aver, &sigma, NULL);
 214     fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[FALSE], aver, sigma);
 215     gmx_stats_get_ase(ef, &aver, &sigma, NULL);
 216     fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[TRUE], aver, sigma);
 217     gmx_stats_done(ed);
 218     gmx_stats_done(ef);
 219     sfree(ed);
 220     sfree(ef);
 221 }
 222
 223 static real rmsd_dist(t_pdbfile *pa, t_pdbfile *pb, gmx_bool bRMSD)
 224 {
 225     int  i;
 226     real rmsd, dist;
 227     rvec acm, bcm, dx;
 228
 229     if (bRMSD)
 230     {
 231         rmsd = 0;
 232         for (i = 0; (i < pa->atoms.nr); i++)
 233         {
 234             rvec_sub(pa->x[i], pb->x[i], dx);
 235             rmsd += iprod(dx, dx);
 236         }
 237         rmsd = sqrt(rmsd/pa->atoms.nr);
 238     }
 239     else
 240     {
 241         dist = 0;
 242         clear_rvec(acm);
 243         clear_rvec(bcm);
 244         for (i = 0; (i < pa->atoms.nr); i++)
 245         {
 246             rvec_inc(acm, pa->x[i]);
 247             rvec_inc(bcm, pb->x[i]);
 248         }
 249         rvec_sub(acm, bcm, dx);
 250         for (i = 0; (i < DIM); i++)
 251         {
 252             dx[i] /= pa->atoms.nr;
 253         }
 254         rmsd = norm(dx);
 255     }
 256     return rmsd;
 257 }
 258
 259 static void line(FILE *fp)
 260 {
 261     fprintf(fp, "                   - - - - - - - - - -\n");
 262 }
 263
 264 static void cluster_em_all(FILE *fp, int npdb, t_pdbfile *pdbf[],
 265                            const char *pdbout, gmx_bool bFree, gmx_bool bRMSD, real cutoff)
 266 {
 267     int   i, j, k;
 268     int  *cndx, ncluster;
 269     real  rmsd;
 270
 271     bFreeSort = bFree;
 272     qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 273
 274     fprintf(fp, "Statistics over all structures:\n");
 275     clust_stat(fp, 0, npdb, pdbf);
 276     line(fp);
 277
 278     /* Index to first structure in a cluster */
 279     snew(cndx, npdb);
 280     ncluster = 1;
 281
 282     for (i = 1; (i < npdb); i++)
 283     {
 284         for (j = 0; (j < ncluster); j++)
 285         {
 286             rmsd = rmsd_dist(pdbf[cndx[j]], pdbf[i], bRMSD);
 287             if (rmsd <= cutoff)
 288             {
 289                 /* Structure i is in cluster j */
 290                 pdbf[i]->cluster_id = pdbf[cndx[j]]->cluster_id;
 291                 break;
 292             }
 293         }
 294         if (j == ncluster)
 295         {
 296             /* New cluster! Cool! */
 297             cndx[ncluster]      = i;
 298             pdbf[i]->cluster_id = ncluster;
 299             ncluster++;
 300         }
 301     }
 302     cndx[ncluster] = npdb;
 303     qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 304
 305     j         = 0;
 306     cndx[j++] = 0;
 307     for (i = 1; (i < npdb); i++)
 308     {
 309         if (pdbf[i]->cluster_id != pdbf[i-1]->cluster_id)
 310         {
 311             cndx[j++] = i;
 312         }
 313     }
 314     cndx[j] = npdb;
 315     if (j != ncluster)
 316     {
 317         gmx_fatal(FARGS, "Consistency error: j = %d, ncluster = %d", j, ncluster);
 318     }
 319
 320     fprintf(fp, "I found %d clusters based on %s and %s with a %.3f nm cut-off\n",
 321             ncluster, etitles[bFree], bRMSD ? "RMSD" : "distance", cutoff);
 322     line(fp);
 323     for (j = 0; (j < ncluster); j++)
 324     {
 325         fprintf(fp, "Cluster: %3d  %s: %10.5f kJ/mol %3d elements\n",
 326                 j, etitles[bFree],
 327                 bFree ? pdbf[cndx[j]]->efree : pdbf[cndx[j]]->edocked,
 328                 cndx[j+1]-cndx[j]);
 329         clust_stat(fp, cndx[j], cndx[j+1], pdbf);
 330         for (k = cndx[j]; (k < cndx[j+1]); k++)
 331         {
 332             fprintf(fp, "  %3d", pdbf[k]->index);
 333         }
 334         fprintf(fp, "\n");
 335         line(fp);
 336     }
 337     sfree(cndx);
 338 }
 339
 340 int gmx_anadock(int argc, char *argv[])
 341 {
 342     const char     *desc[] = {
 343         "[TT]g_anadock[tt] analyses the results of an Autodock run and clusters the",
 344         "structures together, based on distance or RMSD. The docked energy",
 345         "and free energy estimates are analysed, and for each cluster the",
 346         "energy statistics are printed.[PAR]",
 347         "An alternative approach to this is to cluster the structures first",
 348         "using [TT]g_cluster[tt] and then sort the clusters on either lowest",
 349         "energy or average energy."
 350     };
 351     t_filenm        fnm[] = {
 352         { efPDB, "-f", NULL,       ffREAD  },
 353         { efPDB, "-ox", "cluster", ffWRITE },
 354         { efXVG, "-od", "edocked", ffWRITE },
 355         { efXVG, "-of", "efree",   ffWRITE },
 356         { efLOG, "-g",  "anadock", ffWRITE }
 357     };
 358     output_env_t    oenv;
 359 #define NFILE asize(fnm)
 360     static gmx_bool bFree  = FALSE, bRMS = TRUE;
 361     static real     cutoff = 0.2;
 362     t_pargs         pa[]   = {
 363         { "-free",   FALSE, etBOOL, {&bFree},
 364           "Use Free energy estimate from autodock for sorting the classes" },
 365         { "-rms",    FALSE, etBOOL, {&bRMS},
 366           "Cluster on RMS or distance" },
 367         { "-cutoff", FALSE, etREAL, {&cutoff},
 368           "Maximum RMSD/distance for belonging to the same cluster" }
 369     };
 370 #define NPA asize(pa)
 371
 372     FILE       *fp;
 373     t_pdbfile **pdbf = NULL;
 374     int         npdbf;
 375
 376     parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0,
 377                       NULL, &oenv);
 378
 379     fp = ffopen(opt2fn("-g", NFILE, fnm), "w");
 380     please_cite(stdout, "Hetenyi2002b");
 381     please_cite(fp, "Hetenyi2002b");
 382
 383     pdbf = read_em_all(opt2fn("-f", NFILE, fnm), &npdbf);
 384
 385     analyse_em_all(npdbf, pdbf, opt2fn("-od", NFILE, fnm), opt2fn("-of", NFILE, fnm),
 386                    oenv);
 387
 388     cluster_em_all(fp, npdbf, pdbf, opt2fn("-ox", NFILE, fnm),
 389                    bFree, bRMS, cutoff);
 390
 391     thanx(fp);
 392     ffclose(fp);
 393
 394     thanx(stdout);
 395
 396     return 0;
 397 }