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/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, NULL, 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)) != NULL)
  81     {
  82         if (std::strstr(buf, "Intermolecular") != NULL)
  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") != NULL)
  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 = 0;
 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     t_pdbfile *pa, *pb;
 141     real       x;
 142     int        dc;
 143
 144     pa = *(t_pdbfile **)a;
 145     pb = *(t_pdbfile **)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 (bFreeSort = FALSE; (bFreeSort <= TRUE); bFreeSort++)
 186     {
 187         qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 188         fp = xvgropen(bFreeSort ? efree : edocked,
 189                       etitles[bFreeSort], "()", "E (kJ/mol)", oenv);
 190         for (i = 0; (i < npdb); i++)
 191         {
 192             fprintf(fp, "%12lf\n", bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
 193         }
 194         xvgrclose(fp);
 195     }
 196 }
 197
 198 static void clust_stat(FILE *fp, int start, int end, t_pdbfile *pdbf[])
 199 {
 200     int         i;
 201     gmx_stats_t ed, ef;
 202     real        aver, sigma;
 203
 204     ed = gmx_stats_init();
 205     ef = gmx_stats_init();
 206     for (i = start; (i < end); i++)
 207     {
 208         gmx_stats_add_point(ed, i-start, pdbf[i]->edocked, 0, 0);
 209         gmx_stats_add_point(ef, i-start, pdbf[i]->efree, 0, 0);
 210     }
 211     gmx_stats_get_ase(ed, &aver, &sigma, NULL);
 212     fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[FALSE], aver, sigma);
 213     gmx_stats_get_ase(ef, &aver, &sigma, NULL);
 214     fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[TRUE], aver, sigma);
 215     gmx_stats_free(ed);
 216     gmx_stats_free(ef);
 217 }
 218
 219 static real rmsd_dist(t_pdbfile *pa, t_pdbfile *pb, gmx_bool bRMSD)
 220 {
 221     int  i;
 222     real rmsd;
 223     rvec acm, bcm, dx;
 224
 225     if (bRMSD)
 226     {
 227         rmsd = 0;
 228         for (i = 0; (i < pa->atoms.nr); i++)
 229         {
 230             rvec_sub(pa->x[i], pb->x[i], dx);
 231             rmsd += iprod(dx, dx);
 232         }
 233         rmsd = std::sqrt(rmsd/pa->atoms.nr);
 234     }
 235     else
 236     {
 237         clear_rvec(acm);
 238         clear_rvec(bcm);
 239         for (i = 0; (i < pa->atoms.nr); i++)
 240         {
 241             rvec_inc(acm, pa->x[i]);
 242             rvec_inc(bcm, pb->x[i]);
 243         }
 244         rvec_sub(acm, bcm, dx);
 245         for (i = 0; (i < DIM); i++)
 246         {
 247             dx[i] /= pa->atoms.nr;
 248         }
 249         rmsd = norm(dx);
 250     }
 251     return rmsd;
 252 }
 253
 254 static void line(FILE *fp)
 255 {
 256     fprintf(fp, "                   - - - - - - - - - -\n");
 257 }
 258
 259 static void cluster_em_all(FILE *fp, int npdb, t_pdbfile *pdbf[],
 260                            gmx_bool bFree, gmx_bool bRMSD, real cutoff)
 261 {
 262     int   i, j, k;
 263     int  *cndx, ncluster;
 264     real  rmsd;
 265
 266     bFreeSort = bFree;
 267     qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 268
 269     fprintf(fp, "Statistics over all structures:\n");
 270     clust_stat(fp, 0, npdb, pdbf);
 271     line(fp);
 272
 273     /* Index to first structure in a cluster */
 274     snew(cndx, npdb);
 275     ncluster = 1;
 276
 277     for (i = 1; (i < npdb); i++)
 278     {
 279         for (j = 0; (j < ncluster); j++)
 280         {
 281             rmsd = rmsd_dist(pdbf[cndx[j]], pdbf[i], bRMSD);
 282             if (rmsd <= cutoff)
 283             {
 284                 /* Structure i is in cluster j */
 285                 pdbf[i]->cluster_id = pdbf[cndx[j]]->cluster_id;
 286                 break;
 287             }
 288         }
 289         if (j == ncluster)
 290         {
 291             /* New cluster! Cool! */
 292             cndx[ncluster]      = i;
 293             pdbf[i]->cluster_id = ncluster;
 294             ncluster++;
 295         }
 296     }
 297     cndx[ncluster] = npdb;
 298     qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
 299
 300     j         = 0;
 301     cndx[j++] = 0;
 302     for (i = 1; (i < npdb); i++)
 303     {
 304         if (pdbf[i]->cluster_id != pdbf[i-1]->cluster_id)
 305         {
 306             cndx[j++] = i;
 307         }
 308     }
 309     cndx[j] = npdb;
 310     if (j != ncluster)
 311     {
 312         gmx_fatal(FARGS, "Consistency error: j = %d, ncluster = %d", j, ncluster);
 313     }
 314
 315     fprintf(fp, "I found %d clusters based on %s and %s with a %.3f nm cut-off\n",
 316             ncluster, etitles[bFree], bRMSD ? "RMSD" : "distance", cutoff);
 317     line(fp);
 318     for (j = 0; (j < ncluster); j++)
 319     {
 320         fprintf(fp, "Cluster: %3d  %s: %10.5f kJ/mol %3d elements\n",
 321                 j, etitles[bFree],
 322                 bFree ? pdbf[cndx[j]]->efree : pdbf[cndx[j]]->edocked,
 323                 cndx[j+1]-cndx[j]);
 324         clust_stat(fp, cndx[j], cndx[j+1], pdbf);
 325         for (k = cndx[j]; (k < cndx[j+1]); k++)
 326         {
 327             fprintf(fp, "  %3d", pdbf[k]->index);
 328         }
 329         fprintf(fp, "\n");
 330         line(fp);
 331     }
 332     sfree(cndx);
 333 }
 334
 335 int gmx_anadock(int argc, char *argv[])
 336 {
 337     const char       *desc[] = {
 338         "[THISMODULE] analyses the results of an Autodock run and clusters the",
 339         "structures together, based on distance or RMSD. The docked energy",
 340         "and free energy estimates are analysed, and for each cluster the",
 341         "energy statistics are printed.[PAR]",
 342         "An alternative approach to this is to cluster the structures first",
 343         "using [gmx-cluster] and then sort the clusters on either lowest",
 344         "energy or average energy."
 345     };
 346     t_filenm          fnm[] = {
 347         { efPDB, "-f", NULL,       ffREAD  },
 348         { efXVG, "-od", "edocked", ffWRITE },
 349         { efXVG, "-of", "efree",   ffWRITE },
 350         { efLOG, "-g",  "anadock", ffWRITE }
 351     };
 352     gmx_output_env_t *oenv;
 353 #define NFILE asize(fnm)
 354     static gmx_bool   bFree  = FALSE, bRMS = TRUE;
 355     static real       cutoff = 0.2;
 356     t_pargs           pa[]   = {
 357         { "-free",   FALSE, etBOOL, {&bFree},
 358           "Use Free energy estimate from autodock for sorting the classes" },
 359         { "-rms",    FALSE, etBOOL, {&bRMS},
 360           "Cluster on RMS or distance" },
 361         { "-cutoff", FALSE, etREAL, {&cutoff},
 362           "Maximum RMSD/distance for belonging to the same cluster" }
 363     };
 364 #define NPA asize(pa)
 365
 366     FILE       *fp;
 367     t_pdbfile **pdbf = NULL;
 368     int         npdbf;
 369
 370     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0,
 371                            NULL, &oenv))
 372     {
 373         return 0;
 374     }
 375
 376     fp = gmx_ffopen(opt2fn("-g", NFILE, fnm), "w");
 377     please_cite(stdout, "Hetenyi2002b");
 378     please_cite(fp, "Hetenyi2002b");
 379
 380     pdbf = read_em_all(opt2fn("-f", NFILE, fnm), &npdbf);
 381
 382     analyse_em_all(npdbf, pdbf, opt2fn("-od", NFILE, fnm), opt2fn("-of", NFILE, fnm),
 383                    oenv);
 384
 385     cluster_em_all(fp, npdbf, pdbf, bFree, bRMS, cutoff);
 386
 387     gmx_ffclose(fp);
 388
 389     return 0;
 390 }