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