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