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