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 #include "gmxpre.h"
  38
  39 #include <stdlib.h>
  40 #include <string.h>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "gromacs/fileio/xvgr.h"
  45 #include "gromacs/legacyheaders/copyrite.h"
  46 #include "gromacs/legacyheaders/macros.h"
  47 #include "gromacs/legacyheaders/typedefs.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/statistics/statistics.h"
  50 #include "gromacs/utility/cstringutil.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/futil.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 static const char *etitles[] = { "E-docked", "Free Energy" };
  56
  57 typedef struct {
  58     real     edocked, efree;
  59     int      index, cluster_id;
  60     t_atoms  atoms;
  61     rvec    *x;
  62     int      ePBC;
  63     matrix   box;
  64 } t_pdbfile;
  65
  66 static t_pdbfile *read_pdbf(const char *fn)
  67 {
  68     t_pdbfile *pdbf;
  69     double     e;
  70     char       buf[256], *ptr;
  71     int        natoms;
  72     FILE      *fp;
  73
  74     snew(pdbf, 1);
  75     get_stx_coordnum (fn, &natoms);
  76     init_t_atoms(&(pdbf->atoms), natoms, FALSE);
  77     snew(pdbf->x, natoms);
  78     read_stx_conf(fn, buf, &pdbf->atoms, pdbf->x, NULL, &pdbf->ePBC, pdbf->box);
  79     fp = gmx_ffopen(fn, "r");
  80     do
  81     {
  82         ptr = fgets2(buf, 255, fp);
  83         if (ptr)
  84         {
  85             if (strstr(buf, "Intermolecular") != NULL)
  86             {
  87                 ptr = strchr(buf, '=');
  88                 sscanf(ptr+1, "%lf", &e);
  89                 pdbf->edocked = e;
  90             }
  91             else if (strstr(buf, "Estimated Free") != NULL)
  92             {
  93                 ptr = strchr(buf, '=');
  94                 sscanf(ptr+1, "%lf", &e);
  95                 pdbf->efree = e;
  96             }
  97         }
  98     }
  99     while (ptr != NULL);
 100     gmx_ffclose(fp);
 101
 102     return pdbf;
 103 }
 104
 105 static t_pdbfile **read_em_all(const char *fn, int *npdbf)
 106 {
 107     t_pdbfile **pdbf = 0;
 108     int         i, maxpdbf;
 109     char        buf[256], name[256];
 110     gmx_bool    bExist;
 111
 112     strcpy(buf, fn);
 113     buf[strlen(buf)-4] = '\0';
 114     maxpdbf            = 100;
 115     snew(pdbf, maxpdbf);
 116     i = 0;
 117     do
 118     {
 119         sprintf(name, "%s_%d.pdb", buf, i+1);
 120         if ((bExist = gmx_fexist(name)) == TRUE)
 121         {
 122             pdbf[i]        = read_pdbf(name);
 123             pdbf[i]->index = i+1;
 124             i++;
 125             if (i >= maxpdbf)
 126             {
 127                 maxpdbf += 100;
 128                 srenew(pdbf, maxpdbf);
 129             }
 130         }
 131     }
 132     while (bExist);
 133
 134     *npdbf = i;
 135
 136     printf("Found %d pdb files\n", i);
 137
 138     return pdbf;
 139 }
 140
 141 static gmx_bool bFreeSort = FALSE;
 142
 143 static int pdbf_comp(const void *a, const void *b)
 144 {
 145     t_pdbfile *pa, *pb;
 146     real       x;
 147     int        dc;
 148
 149     pa = *(t_pdbfile **)a;
 150     pb = *(t_pdbfile **)b;
 151
 152     dc = pa->cluster_id - pb->cluster_id;
 153
 154     if (dc == 0)
 155     {
 156         if (bFreeSort)
 157         {
 158             x = pa->efree   - pb->efree;
 159         }
 160         else
 161         {
 162             x = pa->edocked - pb->edocked;
 163         }
 164
 165         if (x < 0)
 166         {
 167             return -1;
 168         }
 169         else if (x > 0)
 170         {
 171             return 1;
 172         }
 173         else
 174         {
 175             return 0;
 176         }
 177     }
 178     else
 179     {
 180         return dc;
 181     }
 182 }
 183
 184 static void analyse_em_all(int npdb, t_pdbfile *pdbf[], const char *edocked,
 185                            const char *efree, const output_env_t oenv)
 186 {
 187     FILE *fp;
 188     int   i;
 189
 190     for (bFreeSort = FALSE; (bFreeSort <= TRUE); bFreeSort++)
 191     {
 192         qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 193         fp = xvgropen(bFreeSort ? efree : edocked,
 194                       etitles[bFreeSort], "()", "E (kJ/mol)", oenv);
 195         for (i = 0; (i < npdb); i++)
 196         {
 197             fprintf(fp, "%12lf\n", bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
 198         }
 199         gmx_ffclose(fp);
 200     }
 201 }
 202
 203 static void clust_stat(FILE *fp, int start, int end, t_pdbfile *pdbf[])
 204 {
 205     int         i;
 206     gmx_stats_t ed, ef;
 207     real        aver, sigma;
 208
 209     ed = gmx_stats_init();
 210     ef = gmx_stats_init();
 211     for (i = start; (i < end); i++)
 212     {
 213         gmx_stats_add_point(ed, i-start, pdbf[i]->edocked, 0, 0);
 214         gmx_stats_add_point(ef, i-start, pdbf[i]->efree, 0, 0);
 215     }
 216     gmx_stats_get_ase(ed, &aver, &sigma, NULL);
 217     fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[FALSE], aver, sigma);
 218     gmx_stats_get_ase(ef, &aver, &sigma, NULL);
 219     fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[TRUE], aver, sigma);
 220     gmx_stats_done(ed);
 221     gmx_stats_done(ef);
 222     sfree(ed);
 223     sfree(ef);
 224 }
 225
 226 static real rmsd_dist(t_pdbfile *pa, t_pdbfile *pb, gmx_bool bRMSD)
 227 {
 228     int  i;
 229     real rmsd, dist;
 230     rvec acm, bcm, dx;
 231
 232     if (bRMSD)
 233     {
 234         rmsd = 0;
 235         for (i = 0; (i < pa->atoms.nr); i++)
 236         {
 237             rvec_sub(pa->x[i], pb->x[i], dx);
 238             rmsd += iprod(dx, dx);
 239         }
 240         rmsd = sqrt(rmsd/pa->atoms.nr);
 241     }
 242     else
 243     {
 244         dist = 0;
 245         clear_rvec(acm);
 246         clear_rvec(bcm);
 247         for (i = 0; (i < pa->atoms.nr); i++)
 248         {
 249             rvec_inc(acm, pa->x[i]);
 250             rvec_inc(bcm, pb->x[i]);
 251         }
 252         rvec_sub(acm, bcm, dx);
 253         for (i = 0; (i < DIM); i++)
 254         {
 255             dx[i] /= pa->atoms.nr;
 256         }
 257         rmsd = norm(dx);
 258     }
 259     return rmsd;
 260 }
 261
 262 static void line(FILE *fp)
 263 {
 264     fprintf(fp, "                   - - - - - - - - - -\n");
 265 }
 266
 267 static void cluster_em_all(FILE *fp, int npdb, t_pdbfile *pdbf[],
 268                            gmx_bool bFree, gmx_bool bRMSD, real cutoff)
 269 {
 270     int   i, j, k;
 271     int  *cndx, ncluster;
 272     real  rmsd;
 273
 274     bFreeSort = bFree;
 275     qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 276
 277     fprintf(fp, "Statistics over all structures:\n");
 278     clust_stat(fp, 0, npdb, pdbf);
 279     line(fp);
 280
 281     /* Index to first structure in a cluster */
 282     snew(cndx, npdb);
 283     ncluster = 1;
 284
 285     for (i = 1; (i < npdb); i++)
 286     {
 287         for (j = 0; (j < ncluster); j++)
 288         {
 289             rmsd = rmsd_dist(pdbf[cndx[j]], pdbf[i], bRMSD);
 290             if (rmsd <= cutoff)
 291             {
 292                 /* Structure i is in cluster j */
 293                 pdbf[i]->cluster_id = pdbf[cndx[j]]->cluster_id;
 294                 break;
 295             }
 296         }
 297         if (j == ncluster)
 298         {
 299             /* New cluster! Cool! */
 300             cndx[ncluster]      = i;
 301             pdbf[i]->cluster_id = ncluster;
 302             ncluster++;
 303         }
 304     }
 305     cndx[ncluster] = npdb;
 306     qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 307
 308     j         = 0;
 309     cndx[j++] = 0;
 310     for (i = 1; (i < npdb); i++)
 311     {
 312         if (pdbf[i]->cluster_id != pdbf[i-1]->cluster_id)
 313         {
 314             cndx[j++] = i;
 315         }
 316     }
 317     cndx[j] = npdb;
 318     if (j != ncluster)
 319     {
 320         gmx_fatal(FARGS, "Consistency error: j = %d, ncluster = %d", j, ncluster);
 321     }
 322
 323     fprintf(fp, "I found %d clusters based on %s and %s with a %.3f nm cut-off\n",
 324             ncluster, etitles[bFree], bRMSD ? "RMSD" : "distance", cutoff);
 325     line(fp);
 326     for (j = 0; (j < ncluster); j++)
 327     {
 328         fprintf(fp, "Cluster: %3d  %s: %10.5f kJ/mol %3d elements\n",
 329                 j, etitles[bFree],
 330                 bFree ? pdbf[cndx[j]]->efree : pdbf[cndx[j]]->edocked,
 331                 cndx[j+1]-cndx[j]);
 332         clust_stat(fp, cndx[j], cndx[j+1], pdbf);
 333         for (k = cndx[j]; (k < cndx[j+1]); k++)
 334         {
 335             fprintf(fp, "  %3d", pdbf[k]->index);
 336         }
 337         fprintf(fp, "\n");
 338         line(fp);
 339     }
 340     sfree(cndx);
 341 }
 342
 343 int gmx_anadock(int argc, char *argv[])
 344 {
 345     const char     *desc[] = {
 346         "[THISMODULE] analyses the results of an Autodock run and clusters the",
 347         "structures together, based on distance or RMSD. The docked energy",
 348         "and free energy estimates are analysed, and for each cluster the",
 349         "energy statistics are printed.[PAR]",
 350         "An alternative approach to this is to cluster the structures first",
 351         "using [gmx-cluster] and then sort the clusters on either lowest",
 352         "energy or average energy."
 353     };
 354     t_filenm        fnm[] = {
 355         { efPDB, "-f", NULL,       ffREAD  },
 356         { efXVG, "-od", "edocked", ffWRITE },
 357         { efXVG, "-of", "efree",   ffWRITE },
 358         { efLOG, "-g",  "anadock", ffWRITE }
 359     };
 360     output_env_t    oenv;
 361 #define NFILE asize(fnm)
 362     static gmx_bool bFree  = FALSE, bRMS = TRUE;
 363     static real     cutoff = 0.2;
 364     t_pargs         pa[]   = {
 365         { "-free",   FALSE, etBOOL, {&bFree},
 366           "Use Free energy estimate from autodock for sorting the classes" },
 367         { "-rms",    FALSE, etBOOL, {&bRMS},
 368           "Cluster on RMS or distance" },
 369         { "-cutoff", FALSE, etREAL, {&cutoff},
 370           "Maximum RMSD/distance for belonging to the same cluster" }
 371     };
 372 #define NPA asize(pa)
 373
 374     FILE       *fp;
 375     t_pdbfile **pdbf = NULL;
 376     int         npdbf;
 377
 378     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0,
 379                            NULL, &oenv))
 380     {
 381         return 0;
 382     }
 383
 384     fp = gmx_ffopen(opt2fn("-g", NFILE, fnm), "w");
 385     please_cite(stdout, "Hetenyi2002b");
 386     please_cite(fp, "Hetenyi2002b");
 387
 388     pdbf = read_em_all(opt2fn("-f", NFILE, fnm), &npdbf);
 389
 390     analyse_em_all(npdbf, pdbf, opt2fn("-od", NFILE, fnm), opt2fn("-of", NFILE, fnm),
 391                    oenv);
 392
 393     cluster_em_all(fp, npdbf, pdbf, bFree, bRMS, cutoff);
 394
 395     gmx_thanx(fp);
 396     gmx_ffclose(fp);
 397
 398     return 0;
 399 }