src/gromacs/gmxana/gmx_enemat.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <string.h>
  39 #include <math.h>
  40
  41 #include "string2.h"
  42 #include "typedefs.h"
  43 #include "gmx_fatal.h"
  44 #include "vec.h"
  45 #include "smalloc.h"
  46 #include "enxio.h"
  47 #include "statutil.h"
  48 #include "names.h"
  49 #include "macros.h"
  50 #include "xvgr.h"
  51 #include "gstat.h"
  52 #include "physics.h"
  53 #include "matio.h"
  54 #include "strdb.h"
  55 #include "gmx_ana.h"
  56
  57
  58 static int search_str2(int nstr, char **str, char *key)
  59 {
  60     int  i, n;
  61     int  keylen = strlen(key);
  62     /* Linear search */
  63     n = 0;
  64     while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
  65     {
  66         n++;
  67     }
  68     for (i = 0; (i < nstr); i++)
  69     {
  70         if (gmx_strncasecmp(str[i], key, n) == 0)
  71         {
  72             return i;
  73         }
  74     }
  75
  76     return -1;
  77 }
  78
  79 int gmx_enemat(int argc, char *argv[])
  80 {
  81     const char     *desc[] = {
  82         "[TT]g_enemat[tt] extracts an energy matrix from the energy file ([TT]-f[tt]).",
  83         "With [TT]-groups[tt] a file must be supplied with on each",
  84         "line a group of atoms to be used. For these groups matrix of",
  85         "interaction energies will be extracted from the energy file",
  86         "by looking for energy groups with names corresponding to pairs",
  87         "of groups of atoms, e.g. if your [TT]-groups[tt] file contains:[BR]",
  88         "[TT]2[tt][BR]",
  89         "[TT]Protein[tt][BR]",
  90         "[TT]SOL[tt][BR]",
  91         "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
  92         "'LJ:Protein-SOL' are expected in the energy file (although",
  93         "[TT]g_enemat[tt] is most useful if many groups are analyzed",
  94         "simultaneously). Matrices for different energy types are written",
  95         "out separately, as controlled by the",
  96         "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
  97         "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
  98         "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
  99         "Finally, the total interaction energy energy per group can be ",
 100         "calculated ([TT]-etot[tt]).[PAR]",
 101
 102         "An approximation of the free energy can be calculated using:",
 103         "[MATH]E[SUB]free[sub] = E[SUB]0[sub] + kT [LOG][CHEVRON][EXP](E-E[SUB]0[sub])/kT[exp][chevron][log][math], where '[MATH][CHEVRON][chevron][math]'",
 104         "stands for time-average. A file with reference free energies",
 105         "can be supplied to calculate the free energy difference",
 106         "with some reference state. Group names (e.g. residue names)",
 107         "in the reference file should correspond to the group names",
 108         "as used in the [TT]-groups[tt] file, but a appended number",
 109         "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
 110         "in the comparison."
 111     };
 112     static gmx_bool bSum      = FALSE;
 113     static gmx_bool bMeanEmtx = TRUE;
 114     static int      skip      = 0, nlevels = 20;
 115     static real     cutmax    = 1e20, cutmin = -1e20, reftemp = 300.0;
 116     static gmx_bool bCoulSR   = TRUE, bCoulLR = FALSE, bCoul14 = FALSE;
 117     static gmx_bool bLJSR     = TRUE, bLJLR = FALSE, bLJ14 = FALSE, bBhamSR = FALSE, bBhamLR = FALSE,
 118                     bFree     = TRUE;
 119     t_pargs         pa[]      = {
 120         { "-sum",  FALSE, etBOOL, {&bSum},
 121           "Sum the energy terms selected rather than display them all" },
 122         { "-skip", FALSE, etINT,  {&skip},
 123           "Skip number of frames between data points" },
 124         { "-mean", FALSE, etBOOL, {&bMeanEmtx},
 125           "with [TT]-groups[tt] extracts matrix of mean energies instead of "
 126           "matrix for each timestep" },
 127         { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
 128         { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
 129         { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
 130         { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
 131         { "-coullr", FALSE, etBOOL, {&bCoulLR}, "extract Coulomb LR energies"},
 132         { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
 133         { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
 134         { "-ljlr", FALSE, etBOOL, {&bLJLR}, "extract Lennard-Jones LR energies"},
 135         { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
 136         { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
 137         { "-bhamlr", FALSE, etBOOL, {&bBhamLR}, "extract Buckingham LR energies"},
 138         { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
 139         { "-temp", FALSE, etREAL, {&reftemp},
 140           "reference temperature for free energy calculation"}
 141     };
 142     /* We will define egSP more energy-groups:
 143        egTotal (total energy) */
 144 #define egTotal egNR
 145 #define egSP 1
 146     gmx_bool       egrp_use[egNR+egSP];
 147     ener_file_t    in;
 148     FILE          *out;
 149     int            timecheck = 0;
 150     gmx_enxnm_t   *enm       = NULL;
 151     t_enxframe    *fr;
 152     int            teller = 0;
 153     real           sum;
 154     gmx_bool       bCont, bRef;
 155     gmx_bool       bCutmax, bCutmin;
 156     real         **eneset, *time = NULL;
 157     int           *set, i, j, k, prevk, m = 0, n, nre, nset, nenergy;
 158     char         **groups = NULL;
 159     char           groupname[255], fn[255];
 160     int            ngroups;
 161     t_rgb          rlo, rhi, rmid;
 162     real           emax, emid, emin;
 163     real        ***emat, **etot, *groupnr;
 164     double         beta, expE, **e, *eaver, *efree = NULL, edum;
 165     char           label[234];
 166     char         **ereflines, **erefres = NULL;
 167     real          *eref  = NULL, *edif = NULL;
 168     int            neref = 0;
 169     output_env_t   oenv;
 170
 171     t_filenm       fnm[] = {
 172         { efEDR, "-f", NULL, ffOPTRD },
 173         { efDAT, "-groups", "groups.dat", ffREAD },
 174         { efDAT, "-eref",   "eref.dat", ffOPTRD },
 175         { efXPM, "-emat",   "emat", ffWRITE },
 176         { efXVG, "-etot",   "energy", ffWRITE }
 177     };
 178 #define NFILE asize(fnm)
 179
 180     parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 181                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
 182
 183     egrp_use[egCOULSR] = bCoulSR;
 184     egrp_use[egLJSR]   = bLJSR;
 185     egrp_use[egBHAMSR] = bBhamSR;
 186     egrp_use[egCOULLR] = bCoulLR;
 187     egrp_use[egLJLR]   = bLJLR;
 188     egrp_use[egBHAMLR] = bBhamLR;
 189     egrp_use[egCOUL14] = bCoul14;
 190     egrp_use[egLJ14]   = bLJ14;
 191     egrp_use[egTotal]  = TRUE;
 192
 193     bRef = opt2bSet("-eref", NFILE, fnm);
 194     in   = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
 195     do_enxnms(in, &nre, &enm);
 196
 197     if (nre == 0)
 198     {
 199         gmx_fatal(FARGS, "No energies!\n");
 200     }
 201
 202     bCutmax = opt2parg_bSet("-max", asize(pa), pa);
 203     bCutmin = opt2parg_bSet("-min", asize(pa), pa);
 204
 205     nenergy = 0;
 206
 207     /* Read groupnames from input file and construct selection of
 208        energy groups from it*/
 209
 210     fprintf(stderr, "Will read groupnames from inputfile\n");
 211     ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
 212     fprintf(stderr, "Read %d groups\n", ngroups);
 213     snew(set, sqr(ngroups)*egNR/2);
 214     n     = 0;
 215     prevk = 0;
 216     for (i = 0; (i < ngroups); i++)
 217     {
 218         fprintf(stderr, "\rgroup %d", i);
 219         for (j = i; (j < ngroups); j++)
 220         {
 221             for (m = 0; (m < egNR); m++)
 222             {
 223                 if (egrp_use[m])
 224                 {
 225                     sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
 226 #ifdef DEBUG
 227                     fprintf(stderr, "\r%-15s %5d", groupname, n);
 228 #endif
 229                     for (k = prevk; (k < prevk+nre); k++)
 230                     {
 231                         if (strcmp(enm[k%nre].name, groupname) == 0)
 232                         {
 233                             set[n++] = k;
 234                             break;
 235                         }
 236                     }
 237                     if (k == prevk+nre)
 238                     {
 239                         fprintf(stderr, "WARNING! could not find group %s (%d,%d)"
 240                                 "in energy file\n", groupname, i, j);
 241                     }
 242                     else
 243                     {
 244                         prevk = k;
 245                     }
 246                 }
 247             }
 248         }
 249     }
 250     fprintf(stderr, "\n");
 251     nset = n;
 252     snew(eneset, nset+1);
 253     fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
 254
 255     /* Start reading energy frames */
 256     snew(fr, 1);
 257     do
 258     {
 259         do
 260         {
 261             bCont = do_enx(in, fr);
 262             if (bCont)
 263             {
 264                 timecheck = check_times(fr->t);
 265             }
 266         }
 267         while (bCont && (timecheck < 0));
 268
 269         if (timecheck == 0)
 270         {
 271 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
 272
 273             if (bCont)
 274             {
 275                 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
 276
 277                 if ((nenergy % 1000) == 0)
 278                 {
 279                     srenew(time, nenergy+1000);
 280                     for (i = 0; (i <= nset); i++)
 281                     {
 282                         srenew(eneset[i], nenergy+1000);
 283                     }
 284                 }
 285                 time[nenergy] = fr->t;
 286                 sum           = 0;
 287                 for (i = 0; (i < nset); i++)
 288                 {
 289                     eneset[i][nenergy] = fr->ener[set[i]].e;
 290                     sum               += fr->ener[set[i]].e;
 291                 }
 292                 if (bSum)
 293                 {
 294                     eneset[nset][nenergy] = sum;
 295                 }
 296                 nenergy++;
 297             }
 298             teller++;
 299         }
 300     }
 301     while (bCont && (timecheck == 0));
 302
 303     fprintf(stderr, "\n");
 304
 305     fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
 306             "over %d frames\n", ngroups, nset, nenergy);
 307
 308     snew(emat, egNR+egSP);
 309     for (j = 0; (j < egNR+egSP); j++)
 310     {
 311         if (egrp_use[m])
 312         {
 313             snew(emat[j], ngroups);
 314             for (i = 0; (i < ngroups); i++)
 315             {
 316                 snew(emat[j][i], ngroups);
 317             }
 318         }
 319     }
 320     snew(groupnr, ngroups);
 321     for (i = 0; (i < ngroups); i++)
 322     {
 323         groupnr[i] = i+1;
 324     }
 325     rlo.r  = 1.0, rlo.g  = 0.0, rlo.b  = 0.0;
 326     rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
 327     rhi.r  = 0.0, rhi.g  = 0.0, rhi.b  = 1.0;
 328     if (bMeanEmtx)
 329     {
 330         snew(e, ngroups);
 331         for (i = 0; (i < ngroups); i++)
 332         {
 333             snew(e[i], nenergy);
 334         }
 335         n = 0;
 336         for (i = 0; (i < ngroups); i++)
 337         {
 338             for (j = i; (j < ngroups); j++)
 339             {
 340                 for (m = 0; (m < egNR); m++)
 341                 {
 342                     if (egrp_use[m])
 343                     {
 344                         for (k = 0; (k < nenergy); k++)
 345                         {
 346                             emat[m][i][j] += eneset[n][k];
 347                             e[i][k]       += eneset[n][k]; /* *0.5; */
 348                             e[j][k]       += eneset[n][k]; /* *0.5; */
 349                         }
 350                         n++;
 351                         emat[egTotal][i][j] += emat[m][i][j];
 352                         emat[m][i][j]       /= nenergy;
 353                         emat[m][j][i]        = emat[m][i][j];
 354                     }
 355                 }
 356                 emat[egTotal][i][j] /= nenergy;
 357                 emat[egTotal][j][i]  = emat[egTotal][i][j];
 358             }
 359         }
 360         if (bFree)
 361         {
 362             if (bRef)
 363             {
 364                 fprintf(stderr, "Will read reference energies from inputfile\n");
 365                 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
 366                 fprintf(stderr, "Read %d reference energies\n", neref);
 367                 snew(eref, neref);
 368                 snew(erefres, neref);
 369                 for (i = 0; (i < neref); i++)
 370                 {
 371                     snew(erefres[i], 5);
 372                     sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
 373                     eref[i] = edum;
 374                 }
 375             }
 376             snew(eaver, ngroups);
 377             for (i = 0; (i < ngroups); i++)
 378             {
 379                 for (k = 0; (k < nenergy); k++)
 380                 {
 381                     eaver[i] += e[i][k];
 382                 }
 383                 eaver[i] /= nenergy;
 384             }
 385             beta = 1.0/(BOLTZ*reftemp);
 386             snew(efree, ngroups);
 387             snew(edif, ngroups);
 388             for (i = 0; (i < ngroups); i++)
 389             {
 390                 expE = 0;
 391                 for (k = 0; (k < nenergy); k++)
 392                 {
 393                     expE += exp(beta*(e[i][k]-eaver[i]));
 394                 }
 395                 efree[i] = log(expE/nenergy)/beta + eaver[i];
 396                 if (bRef)
 397                 {
 398                     n = search_str2(neref, erefres, groups[i]);
 399                     if (n != -1)
 400                     {
 401                         edif[i] = efree[i]-eref[n];
 402                     }
 403                     else
 404                     {
 405                         edif[i] = efree[i];
 406                         fprintf(stderr, "WARNING: group %s not found "
 407                                 "in reference energies.\n", groups[i]);
 408                     }
 409                 }
 410                 else
 411                 {
 412                     edif[i] = 0;
 413                 }
 414             }
 415         }
 416
 417         emid             = 0.0; /*(emin+emax)*0.5;*/
 418         egrp_nm[egTotal] = "total";
 419         for (m = 0; (m < egNR+egSP); m++)
 420         {
 421             if (egrp_use[m])
 422             {
 423                 emin = 1e10;
 424                 emax = -1e10;
 425                 for (i = 0; (i < ngroups); i++)
 426                 {
 427                     for (j = i; (j < ngroups); j++)
 428                     {
 429                         if (emat[m][i][j] > emax)
 430                         {
 431                             emax = emat[m][i][j];
 432                         }
 433                         else if (emat[m][i][j] < emin)
 434                         {
 435                             emin = emat[m][i][j];
 436                         }
 437                     }
 438                 }
 439                 if (emax == emin)
 440                 {
 441                     fprintf(stderr, "Matrix of %s energy is uniform at %f "
 442                             "(will not produce output).\n", egrp_nm[m], emax);
 443                 }
 444                 else
 445                 {
 446                     fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
 447                             egrp_nm[m], emin, emax);
 448                     if ((bCutmax) || (emax > cutmax))
 449                     {
 450                         emax = cutmax;
 451                     }
 452                     if ((bCutmin) || (emin < cutmin))
 453                     {
 454                         emin = cutmin;
 455                     }
 456                     if ((emax == cutmax) || (emin == cutmin))
 457                     {
 458                         fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
 459                     }
 460
 461                     sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
 462                     sprintf(label, "%s Interaction Energies", egrp_nm[m]);
 463                     out = ffopen(fn, "w");
 464                     if (emin >= emid)
 465                     {
 466                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 467                                   "Residue Index", "Residue Index",
 468                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 469                                   emid, emax, rmid, rhi, &nlevels);
 470                     }
 471                     else if (emax <= emid)
 472                     {
 473                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 474                                   "Residue Index", "Residue Index",
 475                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 476                                   emin, emid, rlo, rmid, &nlevels);
 477                     }
 478                     else
 479                     {
 480                         write_xpm3(out, 0, label, "Energy (kJ/mol)",
 481                                    "Residue Index", "Residue Index",
 482                                    ngroups, ngroups, groupnr, groupnr, emat[m],
 483                                    emin, emid, emax, rlo, rmid, rhi, &nlevels);
 484                     }
 485                     ffclose(out);
 486                 }
 487             }
 488         }
 489         snew(etot, egNR+egSP);
 490         for (m = 0; (m < egNR+egSP); m++)
 491         {
 492             snew(etot[m], ngroups);
 493             for (i = 0; (i < ngroups); i++)
 494             {
 495                 for (j = 0; (j < ngroups); j++)
 496                 {
 497                     etot[m][i] += emat[m][i][j];
 498                 }
 499             }
 500         }
 501
 502         out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
 503                        oenv);
 504         xvgr_legend(out, 0, NULL, oenv);
 505         j = 0;
 506         for (m = 0; (m < egNR+egSP); m++)
 507         {
 508             if (egrp_use[m])
 509             {
 510                 fprintf(out, "@ legend string %d \"%s\"\n", j++, egrp_nm[m]);
 511             }
 512         }
 513         if (bFree)
 514         {
 515             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Free");
 516         }
 517         if (bFree)
 518         {
 519             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Diff");
 520         }
 521         fprintf(out, "@TYPE xy\n");
 522         fprintf(out, "#%3s", "grp");
 523         for (m = 0; (m < egNR+egSP); m++)
 524         {
 525             if (egrp_use[m])
 526             {
 527                 fprintf(out, " %9s", egrp_nm[m]);
 528             }
 529         }
 530         if (bFree)
 531         {
 532             fprintf(out, " %9s", "Free");
 533         }
 534         if (bFree)
 535         {
 536             fprintf(out, " %9s", "Diff");
 537         }
 538         fprintf(out, "\n");
 539         for (i = 0; (i < ngroups); i++)
 540         {
 541             fprintf(out, "%3.0f", groupnr[i]);
 542             for (m = 0; (m < egNR+egSP); m++)
 543             {
 544                 if (egrp_use[m])
 545                 {
 546                     fprintf(out, " %9.5g", etot[m][i]);
 547                 }
 548             }
 549             if (bFree)
 550             {
 551                 fprintf(out, " %9.5g", efree[i]);
 552             }
 553             if (bRef)
 554             {
 555                 fprintf(out, " %9.5g", edif[i]);
 556             }
 557             fprintf(out, "\n");
 558         }
 559         ffclose(out);
 560     }
 561     else
 562     {
 563         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
 564                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
 565 /*
 566     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 567     n=0;
 568     emin=emax=0.0;
 569     for (k=0; (k<nenergy); k++) {
 570       for (i=0; (i<ngroups); i++)
 571     for (j=i+1; (j<ngroups); j++)
 572       emat[i][j]=eneset[n][k];
 573       sprintf(label,"t=%.0f ps",time[k]);
 574       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 575       n++;
 576     }
 577     ffclose(out);
 578  */
 579     }
 580     close_enx(in);
 581
 582     return 0;
 583 }