src/gromacs/gmxana/gmx_anadock.cpp

   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,2015, 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 <cmath>
  40 #include <cstdlib>
  41 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/xvgr.h"
  46 #include "gromacs/gmxana/gmx_ana.h"
  47 #include "gromacs/legacyheaders/copyrite.h"
  48 #include "gromacs/legacyheaders/macros.h"
  49 #include "gromacs/legacyheaders/typedefs.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/statistics/statistics.h"
  52 #include "gromacs/utility/cstringutil.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/futil.h"
  55 #include "gromacs/utility/smalloc.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 (std::strstr(buf, "Intermolecular") != NULL)
  88             {
  89                 ptr = std::strchr(buf, '=');
  90                 sscanf(ptr+1, "%lf", &e);
  91                 pdbf->edocked = e;
  92             }
  93             else if (std::strstr(buf, "Estimated Free") != NULL)
  94             {
  95                 ptr = std::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     std::strcpy(buf, fn);
 115     buf[std::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         xvgrclose(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;
 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 = std::sqrt(rmsd/pa->atoms.nr);
 243     }
 244     else
 245     {
 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 }