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     for (i = 0; (i < egNR+egSP); i++)
 189     {
 190         egrp_use[i] = FALSE;
 191     }
 192     egrp_use[egCOULSR] = bCoulSR;
 193     egrp_use[egLJSR]   = bLJSR;
 194     egrp_use[egBHAMSR] = bBhamSR;
 195     egrp_use[egCOULLR] = bCoulLR;
 196     egrp_use[egLJLR]   = bLJLR;
 197     egrp_use[egBHAMLR] = bBhamLR;
 198     egrp_use[egCOUL14] = bCoul14;
 199     egrp_use[egLJ14]   = bLJ14;
 200     egrp_use[egTotal]  = TRUE;
 201
 202     bRef = opt2bSet("-eref", NFILE, fnm);
 203     in   = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
 204     do_enxnms(in, &nre, &enm);
 205
 206     if (nre == 0)
 207     {
 208         gmx_fatal(FARGS, "No energies!\n");
 209     }
 210
 211     bCutmax = opt2parg_bSet("-max", asize(pa), pa);
 212     bCutmin = opt2parg_bSet("-min", asize(pa), pa);
 213
 214     nenergy = 0;
 215
 216     /* Read groupnames from input file and construct selection of
 217        energy groups from it*/
 218
 219     fprintf(stderr, "Will read groupnames from inputfile\n");
 220     ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
 221     fprintf(stderr, "Read %d groups\n", ngroups);
 222     snew(set, sqr(ngroups)*egNR/2);
 223     n     = 0;
 224     prevk = 0;
 225     for (i = 0; (i < ngroups); i++)
 226     {
 227         fprintf(stderr, "\rgroup %d", i);
 228         for (j = i; (j < ngroups); j++)
 229         {
 230             for (m = 0; (m < egNR); m++)
 231             {
 232                 if (egrp_use[m])
 233                 {
 234                     sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
 235 #ifdef DEBUG
 236                     fprintf(stderr, "\r%-15s %5d", groupname, n);
 237 #endif
 238                     for (k = prevk; (k < prevk+nre); k++)
 239                     {
 240                         if (strcmp(enm[k%nre].name, groupname) == 0)
 241                         {
 242                             set[n++] = k;
 243                             break;
 244                         }
 245                     }
 246                     if (k == prevk+nre)
 247                     {
 248                         fprintf(stderr, "WARNING! could not find group %s (%d,%d)"
 249                                 "in energy file\n", groupname, i, j);
 250                     }
 251                     else
 252                     {
 253                         prevk = k;
 254                     }
 255                 }
 256             }
 257         }
 258     }
 259     fprintf(stderr, "\n");
 260     nset = n;
 261     snew(eneset, nset+1);
 262     fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
 263
 264     /* Start reading energy frames */
 265     snew(fr, 1);
 266     do
 267     {
 268         do
 269         {
 270             bCont = do_enx(in, fr);
 271             if (bCont)
 272             {
 273                 timecheck = check_times(fr->t);
 274             }
 275         }
 276         while (bCont && (timecheck < 0));
 277
 278         if (timecheck == 0)
 279         {
 280 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
 281
 282             if (bCont)
 283             {
 284                 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
 285
 286                 if ((nenergy % 1000) == 0)
 287                 {
 288                     srenew(time, nenergy+1000);
 289                     for (i = 0; (i <= nset); i++)
 290                     {
 291                         srenew(eneset[i], nenergy+1000);
 292                     }
 293                 }
 294                 time[nenergy] = fr->t;
 295                 sum           = 0;
 296                 for (i = 0; (i < nset); i++)
 297                 {
 298                     eneset[i][nenergy] = fr->ener[set[i]].e;
 299                     sum               += fr->ener[set[i]].e;
 300                 }
 301                 if (bSum)
 302                 {
 303                     eneset[nset][nenergy] = sum;
 304                 }
 305                 nenergy++;
 306             }
 307             teller++;
 308         }
 309     }
 310     while (bCont && (timecheck == 0));
 311
 312     fprintf(stderr, "\n");
 313
 314     fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
 315             "over %d frames\n", ngroups, nset, nenergy);
 316
 317     snew(emat, egNR+egSP);
 318     for (j = 0; (j < egNR+egSP); j++)
 319     {
 320         if (egrp_use[m])
 321         {
 322             snew(emat[j], ngroups);
 323             for (i = 0; (i < ngroups); i++)
 324             {
 325                 snew(emat[j][i], ngroups);
 326             }
 327         }
 328     }
 329     snew(groupnr, ngroups);
 330     for (i = 0; (i < ngroups); i++)
 331     {
 332         groupnr[i] = i+1;
 333     }
 334     rlo.r  = 1.0, rlo.g  = 0.0, rlo.b  = 0.0;
 335     rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
 336     rhi.r  = 0.0, rhi.g  = 0.0, rhi.b  = 1.0;
 337     if (bMeanEmtx)
 338     {
 339         snew(e, ngroups);
 340         for (i = 0; (i < ngroups); i++)
 341         {
 342             snew(e[i], nenergy);
 343         }
 344         n = 0;
 345         for (i = 0; (i < ngroups); i++)
 346         {
 347             for (j = i; (j < ngroups); j++)
 348             {
 349                 for (m = 0; (m < egNR); m++)
 350                 {
 351                     if (egrp_use[m])
 352                     {
 353                         for (k = 0; (k < nenergy); k++)
 354                         {
 355                             emat[m][i][j] += eneset[n][k];
 356                             e[i][k]       += eneset[n][k]; /* *0.5; */
 357                             e[j][k]       += eneset[n][k]; /* *0.5; */
 358                         }
 359                         n++;
 360                         emat[egTotal][i][j] += emat[m][i][j];
 361                         emat[m][i][j]       /= nenergy;
 362                         emat[m][j][i]        = emat[m][i][j];
 363                     }
 364                 }
 365                 emat[egTotal][i][j] /= nenergy;
 366                 emat[egTotal][j][i]  = emat[egTotal][i][j];
 367             }
 368         }
 369         if (bFree)
 370         {
 371             if (bRef)
 372             {
 373                 fprintf(stderr, "Will read reference energies from inputfile\n");
 374                 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
 375                 fprintf(stderr, "Read %d reference energies\n", neref);
 376                 snew(eref, neref);
 377                 snew(erefres, neref);
 378                 for (i = 0; (i < neref); i++)
 379                 {
 380                     snew(erefres[i], 5);
 381                     sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
 382                     eref[i] = edum;
 383                 }
 384             }
 385             snew(eaver, ngroups);
 386             for (i = 0; (i < ngroups); i++)
 387             {
 388                 for (k = 0; (k < nenergy); k++)
 389                 {
 390                     eaver[i] += e[i][k];
 391                 }
 392                 eaver[i] /= nenergy;
 393             }
 394             beta = 1.0/(BOLTZ*reftemp);
 395             snew(efree, ngroups);
 396             snew(edif, ngroups);
 397             for (i = 0; (i < ngroups); i++)
 398             {
 399                 expE = 0;
 400                 for (k = 0; (k < nenergy); k++)
 401                 {
 402                     expE += exp(beta*(e[i][k]-eaver[i]));
 403                 }
 404                 efree[i] = log(expE/nenergy)/beta + eaver[i];
 405                 if (bRef)
 406                 {
 407                     n = search_str2(neref, erefres, groups[i]);
 408                     if (n != -1)
 409                     {
 410                         edif[i] = efree[i]-eref[n];
 411                     }
 412                     else
 413                     {
 414                         edif[i] = efree[i];
 415                         fprintf(stderr, "WARNING: group %s not found "
 416                                 "in reference energies.\n", groups[i]);
 417                     }
 418                 }
 419                 else
 420                 {
 421                     edif[i] = 0;
 422                 }
 423             }
 424         }
 425
 426         emid             = 0.0; /*(emin+emax)*0.5;*/
 427         egrp_nm[egTotal] = "total";
 428         for (m = 0; (m < egNR+egSP); m++)
 429         {
 430             if (egrp_use[m])
 431             {
 432                 emin = 1e10;
 433                 emax = -1e10;
 434                 for (i = 0; (i < ngroups); i++)
 435                 {
 436                     for (j = i; (j < ngroups); j++)
 437                     {
 438                         if (emat[m][i][j] > emax)
 439                         {
 440                             emax = emat[m][i][j];
 441                         }
 442                         else if (emat[m][i][j] < emin)
 443                         {
 444                             emin = emat[m][i][j];
 445                         }
 446                     }
 447                 }
 448                 if (emax == emin)
 449                 {
 450                     fprintf(stderr, "Matrix of %s energy is uniform at %f "
 451                             "(will not produce output).\n", egrp_nm[m], emax);
 452                 }
 453                 else
 454                 {
 455                     fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
 456                             egrp_nm[m], emin, emax);
 457                     if ((bCutmax) || (emax > cutmax))
 458                     {
 459                         emax = cutmax;
 460                     }
 461                     if ((bCutmin) || (emin < cutmin))
 462                     {
 463                         emin = cutmin;
 464                     }
 465                     if ((emax == cutmax) || (emin == cutmin))
 466                     {
 467                         fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
 468                     }
 469
 470                     sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
 471                     sprintf(label, "%s Interaction Energies", egrp_nm[m]);
 472                     out = ffopen(fn, "w");
 473                     if (emin >= emid)
 474                     {
 475                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 476                                   "Residue Index", "Residue Index",
 477                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 478                                   emid, emax, rmid, rhi, &nlevels);
 479                     }
 480                     else if (emax <= emid)
 481                     {
 482                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 483                                   "Residue Index", "Residue Index",
 484                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 485                                   emin, emid, rlo, rmid, &nlevels);
 486                     }
 487                     else
 488                     {
 489                         write_xpm3(out, 0, label, "Energy (kJ/mol)",
 490                                    "Residue Index", "Residue Index",
 491                                    ngroups, ngroups, groupnr, groupnr, emat[m],
 492                                    emin, emid, emax, rlo, rmid, rhi, &nlevels);
 493                     }
 494                     ffclose(out);
 495                 }
 496             }
 497         }
 498         snew(etot, egNR+egSP);
 499         for (m = 0; (m < egNR+egSP); m++)
 500         {
 501             snew(etot[m], ngroups);
 502             for (i = 0; (i < ngroups); i++)
 503             {
 504                 for (j = 0; (j < ngroups); j++)
 505                 {
 506                     etot[m][i] += emat[m][i][j];
 507                 }
 508             }
 509         }
 510
 511         out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
 512                        oenv);
 513         xvgr_legend(out, 0, NULL, oenv);
 514         j = 0;
 515         for (m = 0; (m < egNR+egSP); m++)
 516         {
 517             if (egrp_use[m])
 518             {
 519                 fprintf(out, "@ legend string %d \"%s\"\n", j++, egrp_nm[m]);
 520             }
 521         }
 522         if (bFree)
 523         {
 524             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Free");
 525         }
 526         if (bFree)
 527         {
 528             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Diff");
 529         }
 530         fprintf(out, "@TYPE xy\n");
 531         fprintf(out, "#%3s", "grp");
 532         for (m = 0; (m < egNR+egSP); m++)
 533         {
 534             if (egrp_use[m])
 535             {
 536                 fprintf(out, " %9s", egrp_nm[m]);
 537             }
 538         }
 539         if (bFree)
 540         {
 541             fprintf(out, " %9s", "Free");
 542         }
 543         if (bFree)
 544         {
 545             fprintf(out, " %9s", "Diff");
 546         }
 547         fprintf(out, "\n");
 548         for (i = 0; (i < ngroups); i++)
 549         {
 550             fprintf(out, "%3.0f", groupnr[i]);
 551             for (m = 0; (m < egNR+egSP); m++)
 552             {
 553                 if (egrp_use[m])
 554                 {
 555                     fprintf(out, " %9.5g", etot[m][i]);
 556                 }
 557             }
 558             if (bFree)
 559             {
 560                 fprintf(out, " %9.5g", efree[i]);
 561             }
 562             if (bRef)
 563             {
 564                 fprintf(out, " %9.5g", edif[i]);
 565             }
 566             fprintf(out, "\n");
 567         }
 568         ffclose(out);
 569     }
 570     else
 571     {
 572         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
 573                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
 574 /*
 575     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 576     n=0;
 577     emin=emax=0.0;
 578     for (k=0; (k<nenergy); k++) {
 579       for (i=0; (i<ngroups); i++)
 580     for (j=i+1; (j<ngroups); j++)
 581       emat[i][j]=eneset[n][k];
 582       sprintf(label,"t=%.0f ps",time[k]);
 583       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 584       n++;
 585     }
 586     ffclose(out);
 587  */
 588     }
 589     close_enx(in);
 590
 591     thanx(stderr);
 592
 593     return 0;
 594 }