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