src/gromacs/gmxana/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  * Copyright (c) 2013,2014, 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 <string.h>
  40 #include <math.h>
  41
  42 #include "gromacs/utility/cstringutil.h"
  43 #include "gromacs/legacyheaders/typedefs.h"
  44 #include "gromacs/utility/fatalerror.h"
  45 #include "gromacs/math/vec.h"
  46 #include "gromacs/utility/smalloc.h"
  47 #include "gromacs/fileio/enxio.h"
  48 #include "gromacs/commandline/pargs.h"
  49 #include "gromacs/legacyheaders/names.h"
  50 #include "gromacs/legacyheaders/macros.h"
  51 #include "gromacs/fileio/xvgr.h"
  52 #include "gstat.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/fileio/matio.h"
  55 #include "gromacs/fileio/strdb.h"
  56 #include "gmx_ana.h"
  57 #include "gromacs/fileio/trxio.h"
  58
  59
  60 static int search_str2(int nstr, char **str, char *key)
  61 {
  62     int  i, n;
  63     int  keylen = strlen(key);
  64     /* Linear search */
  65     n = 0;
  66     while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
  67     {
  68         n++;
  69     }
  70     for (i = 0; (i < nstr); i++)
  71     {
  72         if (gmx_strncasecmp(str[i], key, n) == 0)
  73         {
  74             return i;
  75         }
  76     }
  77
  78     return -1;
  79 }
  80
  81 int gmx_enemat(int argc, char *argv[])
  82 {
  83     const char     *desc[] = {
  84         "[THISMODULE] extracts an energy matrix from the energy file ([TT]-f[tt]).",
  85         "With [TT]-groups[tt] a file must be supplied with on each",
  86         "line a group of atoms to be used. For these groups matrix of",
  87         "interaction energies will be extracted from the energy file",
  88         "by looking for energy groups with names corresponding to pairs",
  89         "of groups of atoms, e.g. if your [TT]-groups[tt] file contains:[BR]",
  90         "[TT]2[tt][BR]",
  91         "[TT]Protein[tt][BR]",
  92         "[TT]SOL[tt][BR]",
  93         "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
  94         "'LJ:Protein-SOL' are expected in the energy file (although",
  95         "[THISMODULE] is most useful if many groups are analyzed",
  96         "simultaneously). Matrices for different energy types are written",
  97         "out separately, as controlled by the",
  98         "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
  99         "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
 100         "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
 101         "Finally, the total interaction energy energy per group can be ",
 102         "calculated ([TT]-etot[tt]).[PAR]",
 103
 104         "An approximation of the free energy can be calculated using:",
 105         "[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]'",
 106         "stands for time-average. A file with reference free energies",
 107         "can be supplied to calculate the free energy difference",
 108         "with some reference state. Group names (e.g. residue names)",
 109         "in the reference file should correspond to the group names",
 110         "as used in the [TT]-groups[tt] file, but a appended number",
 111         "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
 112         "in the comparison."
 113     };
 114     static gmx_bool bSum      = FALSE;
 115     static gmx_bool bMeanEmtx = TRUE;
 116     static int      skip      = 0, nlevels = 20;
 117     static real     cutmax    = 1e20, cutmin = -1e20, reftemp = 300.0;
 118     static gmx_bool bCoulSR   = TRUE, bCoulLR = FALSE, bCoul14 = FALSE;
 119     static gmx_bool bLJSR     = TRUE, bLJLR = FALSE, bLJ14 = FALSE, bBhamSR = FALSE, bBhamLR = FALSE,
 120                     bFree     = TRUE;
 121     t_pargs         pa[]      = {
 122         { "-sum",  FALSE, etBOOL, {&bSum},
 123           "Sum the energy terms selected rather than display them all" },
 124         { "-skip", FALSE, etINT,  {&skip},
 125           "Skip number of frames between data points" },
 126         { "-mean", FALSE, etBOOL, {&bMeanEmtx},
 127           "with [TT]-groups[tt] extracts matrix of mean energies instead of "
 128           "matrix for each timestep" },
 129         { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
 130         { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
 131         { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
 132         { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
 133         { "-coullr", FALSE, etBOOL, {&bCoulLR}, "extract Coulomb LR energies"},
 134         { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
 135         { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
 136         { "-ljlr", FALSE, etBOOL, {&bLJLR}, "extract Lennard-Jones LR energies"},
 137         { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
 138         { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
 139         { "-bhamlr", FALSE, etBOOL, {&bBhamLR}, "extract Buckingham LR energies"},
 140         { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
 141         { "-temp", FALSE, etREAL, {&reftemp},
 142           "reference temperature for free energy calculation"}
 143     };
 144     /* We will define egSP more energy-groups:
 145        egTotal (total energy) */
 146 #define egTotal egNR
 147 #define egSP 1
 148     gmx_bool       egrp_use[egNR+egSP];
 149     ener_file_t    in;
 150     FILE          *out;
 151     int            timecheck = 0;
 152     gmx_enxnm_t   *enm       = NULL;
 153     t_enxframe    *fr;
 154     int            teller = 0;
 155     real           sum;
 156     gmx_bool       bCont, bRef;
 157     gmx_bool       bCutmax, bCutmin;
 158     real         **eneset, *time = NULL;
 159     int           *set, i, j, k, prevk, m = 0, n, nre, nset, nenergy;
 160     char         **groups = NULL;
 161     char           groupname[255], fn[255];
 162     int            ngroups;
 163     t_rgb          rlo, rhi, rmid;
 164     real           emax, emid, emin;
 165     real        ***emat, **etot, *groupnr;
 166     double         beta, expE, **e, *eaver, *efree = NULL, edum;
 167     char           label[234];
 168     char         **ereflines, **erefres = NULL;
 169     real          *eref  = NULL, *edif = NULL;
 170     int            neref = 0;
 171     output_env_t   oenv;
 172
 173     t_filenm       fnm[] = {
 174         { efEDR, "-f", NULL, ffOPTRD },
 175         { efDAT, "-groups", "groups", ffREAD },
 176         { efDAT, "-eref",   "eref",   ffOPTRD },
 177         { efXPM, "-emat",   "emat",   ffWRITE },
 178         { efXVG, "-etot",   "energy", ffWRITE }
 179     };
 180 #define NFILE asize(fnm)
 181
 182     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 183                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 184     {
 185         return 0;
 186     }
 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 = gmx_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                     gmx_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         if (output_env_get_print_xvgr_codes(oenv))
 516         {
 517             char str1[STRLEN], str2[STRLEN];
 518             if (output_env_get_xvg_format(oenv) == exvgXMGR)
 519             {
 520                 sprintf(str1, "@ legend string ");
 521                 sprintf(str2, " ");
 522             }
 523             else
 524             {
 525                 sprintf(str1, "@ s");
 526                 sprintf(str2, " legend ");
 527             }
 528
 529             for (m = 0; (m < egNR+egSP); m++)
 530             {
 531                 if (egrp_use[m])
 532                 {
 533                     fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
 534                 }
 535             }
 536             if (bFree)
 537             {
 538                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
 539             }
 540             if (bFree)
 541             {
 542                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
 543             }
 544             fprintf(out, "@TYPE xy\n");
 545             fprintf(out, "#%3s", "grp");
 546
 547             for (m = 0; (m < egNR+egSP); m++)
 548             {
 549                 if (egrp_use[m])
 550                 {
 551                     fprintf(out, " %9s", egrp_nm[m]);
 552                 }
 553             }
 554             if (bFree)
 555             {
 556                 fprintf(out, " %9s", "Free");
 557             }
 558             if (bFree)
 559             {
 560                 fprintf(out, " %9s", "Diff");
 561             }
 562             fprintf(out, "\n");
 563         }
 564         for (i = 0; (i < ngroups); i++)
 565         {
 566             fprintf(out, "%3.0f", groupnr[i]);
 567             for (m = 0; (m < egNR+egSP); m++)
 568             {
 569                 if (egrp_use[m])
 570                 {
 571                     fprintf(out, " %9.5g", etot[m][i]);
 572                 }
 573             }
 574             if (bFree)
 575             {
 576                 fprintf(out, " %9.5g", efree[i]);
 577             }
 578             if (bRef)
 579             {
 580                 fprintf(out, " %9.5g", edif[i]);
 581             }
 582             fprintf(out, "\n");
 583         }
 584         gmx_ffclose(out);
 585     }
 586     else
 587     {
 588         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
 589                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
 590 /*
 591     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 592     n=0;
 593     emin=emax=0.0;
 594     for (k=0; (k<nenergy); k++) {
 595       for (i=0; (i<ngroups); i++)
 596     for (j=i+1; (j<ngroups); j++)
 597       emat[i][j]=eneset[n][k];
 598       sprintf(label,"t=%.0f ps",time[k]);
 599       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 600       n++;
 601     }
 602     gmx_ffclose(out);
 603  */
 604     }
 605     close_enx(in);
 606
 607     return 0;
 608 }