src/gromacs/gmxana/gmx_anadock.c

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