src/gromacs/gmxana/gmx_enemat.cpp

   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,2015,2016, 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 <cmath>
  40 #include <cstring>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/enxio.h"
  44 #include "gromacs/fileio/matio.h"
  45 #include "gromacs/fileio/trxio.h"
  46 #include "gromacs/fileio/xvgr.h"
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/gmxana/gstat.h"
  49 #include "gromacs/math/functions.h"
  50 #include "gromacs/math/units.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/mdtypes/forcerec.h"
  53 #include "gromacs/mdtypes/md_enums.h"
  54 #include "gromacs/utility/arraysize.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/smalloc.h"
  59 #include "gromacs/utility/strdb.h"
  60
  61
  62 static int search_str2(int nstr, char **str, char *key)
  63 {
  64     int  i, n;
  65     int  keylen = std::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         "[THISMODULE] 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::",
  92         "",
  93         "    2",
  94         "    Protein",
  95         "    SOL",
  96         "",
  97         "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
  98         "'LJ:Protein-SOL' are expected in the energy file (although",
  99         "[THISMODULE] is most useful if many groups are analyzed",
 100         "simultaneously). Matrices for different energy types are written",
 101         "out separately, as controlled by the",
 102         "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
 103         "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
 104         "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
 105         "Finally, the total interaction energy energy per group can be ",
 106         "calculated ([TT]-etot[tt]).[PAR]",
 107
 108         "An approximation of the free energy can be calculated using:",
 109         "[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]'",
 110         "stands for time-average. A file with reference free energies",
 111         "can be supplied to calculate the free energy difference",
 112         "with some reference state. Group names (e.g. residue names)",
 113         "in the reference file should correspond to the group names",
 114         "as used in the [TT]-groups[tt] file, but a appended number",
 115         "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
 116         "in the comparison."
 117     };
 118     static gmx_bool bSum      = FALSE;
 119     static gmx_bool bMeanEmtx = TRUE;
 120     static int      skip      = 0, nlevels = 20;
 121     static real     cutmax    = 1e20, cutmin = -1e20, reftemp = 300.0;
 122     static gmx_bool bCoulSR   = TRUE, bCoul14 = FALSE;
 123     static gmx_bool bLJSR     = TRUE, bLJ14 = FALSE, bBhamSR = FALSE,
 124                     bFree     = TRUE;
 125     t_pargs         pa[]      = {
 126         { "-sum",  FALSE, etBOOL, {&bSum},
 127           "Sum the energy terms selected rather than display them all" },
 128         { "-skip", FALSE, etINT,  {&skip},
 129           "Skip number of frames between data points" },
 130         { "-mean", FALSE, etBOOL, {&bMeanEmtx},
 131           "with [TT]-groups[tt] extracts matrix of mean energies instead of "
 132           "matrix for each timestep" },
 133         { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
 134         { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
 135         { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
 136         { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
 137         { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
 138         { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
 139         { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
 140         { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
 141         { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
 142         { "-temp", FALSE, etREAL, {&reftemp},
 143           "reference temperature for free energy calculation"}
 144     };
 145     /* We will define egSP more energy-groups:
 146        egTotal (total energy) */
 147 #define egTotal egNR
 148 #define egSP 1
 149     gmx_bool          egrp_use[egNR+egSP];
 150     ener_file_t       in;
 151     FILE             *out;
 152     int               timecheck = 0;
 153     gmx_enxnm_t      *enm       = NULL;
 154     t_enxframe       *fr;
 155     int               teller = 0;
 156     real              sum;
 157     gmx_bool          bCont, bRef;
 158     gmx_bool          bCutmax, bCutmin;
 159     real            **eneset, *time = NULL;
 160     int              *set, i, j, k, prevk, m = 0, n, nre, nset, nenergy;
 161     char            **groups = NULL;
 162     char              groupname[255], fn[255];
 163     int               ngroups;
 164     t_rgb             rlo, rhi, rmid;
 165     real              emax, emid, emin;
 166     real           ***emat, **etot, *groupnr;
 167     double            beta, expE, **e, *eaver, *efree = NULL, edum;
 168     char              label[234];
 169     char            **ereflines, **erefres = NULL;
 170     real             *eref  = NULL, *edif = NULL;
 171     int               neref = 0;
 172     gmx_output_env_t *oenv;
 173
 174     t_filenm          fnm[] = {
 175         { efEDR, "-f", NULL, ffOPTRD },
 176         { efDAT, "-groups", "groups", ffREAD },
 177         { efDAT, "-eref",   "eref",   ffOPTRD },
 178         { efXPM, "-emat",   "emat",   ffWRITE },
 179         { efXVG, "-etot",   "energy", ffWRITE }
 180     };
 181 #define NFILE asize(fnm)
 182
 183     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 184                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 185     {
 186         return 0;
 187     }
 188
 189     for (i = 0; (i < egNR+egSP); i++)
 190     {
 191         egrp_use[i] = FALSE;
 192     }
 193     egrp_use[egCOULSR] = bCoulSR;
 194     egrp_use[egLJSR]   = bLJSR;
 195     egrp_use[egBHAMSR] = bBhamSR;
 196     egrp_use[egCOUL14] = bCoul14;
 197     egrp_use[egLJ14]   = bLJ14;
 198     egrp_use[egTotal]  = TRUE;
 199
 200     bRef = opt2bSet("-eref", NFILE, fnm);
 201     in   = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
 202     do_enxnms(in, &nre, &enm);
 203
 204     if (nre == 0)
 205     {
 206         gmx_fatal(FARGS, "No energies!\n");
 207     }
 208
 209     bCutmax = opt2parg_bSet("-max", asize(pa), pa);
 210     bCutmin = opt2parg_bSet("-min", asize(pa), pa);
 211
 212     nenergy = 0;
 213
 214     /* Read groupnames from input file and construct selection of
 215        energy groups from it*/
 216
 217     fprintf(stderr, "Will read groupnames from inputfile\n");
 218     ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
 219     fprintf(stderr, "Read %d groups\n", ngroups);
 220     snew(set, static_cast<size_t>(gmx::square(ngroups)*egNR/2));
 221     n     = 0;
 222     prevk = 0;
 223     for (i = 0; (i < ngroups); i++)
 224     {
 225         fprintf(stderr, "\rgroup %d", i);
 226         fflush(stderr);
 227         for (j = i; (j < ngroups); j++)
 228         {
 229             for (m = 0; (m < egNR); m++)
 230             {
 231                 if (egrp_use[m])
 232                 {
 233                     sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
 234 #ifdef DEBUG
 235                     fprintf(stderr, "\r%-15s %5d", groupname, n);
 236                     fflush(stderr);
 237 #endif
 238                     for (k = prevk; (k < prevk+nre); k++)
 239                     {
 240                         if (std::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                 fflush(stderr);
 286
 287                 if ((nenergy % 1000) == 0)
 288                 {
 289                     srenew(time, nenergy+1000);
 290                     for (i = 0; (i <= nset); i++)
 291                     {
 292                         srenew(eneset[i], nenergy+1000);
 293                     }
 294                 }
 295                 time[nenergy] = fr->t;
 296                 sum           = 0;
 297                 for (i = 0; (i < nset); i++)
 298                 {
 299                     eneset[i][nenergy] = fr->ener[set[i]].e;
 300                     sum               += fr->ener[set[i]].e;
 301                 }
 302                 if (bSum)
 303                 {
 304                     eneset[nset][nenergy] = sum;
 305                 }
 306                 nenergy++;
 307             }
 308             teller++;
 309         }
 310     }
 311     while (bCont && (timecheck == 0));
 312
 313     fprintf(stderr, "\n");
 314
 315     fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
 316             "over %d frames\n", ngroups, nset, nenergy);
 317
 318     snew(emat, egNR+egSP);
 319     for (j = 0; (j < egNR+egSP); j++)
 320     {
 321         if (egrp_use[m])
 322         {
 323             snew(emat[j], ngroups);
 324             for (i = 0; (i < ngroups); i++)
 325             {
 326                 snew(emat[j][i], ngroups);
 327             }
 328         }
 329     }
 330     snew(groupnr, ngroups);
 331     for (i = 0; (i < ngroups); i++)
 332     {
 333         groupnr[i] = i+1;
 334     }
 335     rlo.r  = 1.0, rlo.g  = 0.0, rlo.b  = 0.0;
 336     rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
 337     rhi.r  = 0.0, rhi.g  = 0.0, rhi.b  = 1.0;
 338     if (bMeanEmtx)
 339     {
 340         snew(e, ngroups);
 341         for (i = 0; (i < ngroups); i++)
 342         {
 343             snew(e[i], nenergy);
 344         }
 345         n = 0;
 346         for (i = 0; (i < ngroups); i++)
 347         {
 348             for (j = i; (j < ngroups); j++)
 349             {
 350                 for (m = 0; (m < egNR); m++)
 351                 {
 352                     if (egrp_use[m])
 353                     {
 354                         for (k = 0; (k < nenergy); k++)
 355                         {
 356                             emat[m][i][j] += eneset[n][k];
 357                             e[i][k]       += eneset[n][k]; /* *0.5; */
 358                             e[j][k]       += eneset[n][k]; /* *0.5; */
 359                         }
 360                         n++;
 361                         emat[egTotal][i][j] += emat[m][i][j];
 362                         emat[m][i][j]       /= nenergy;
 363                         emat[m][j][i]        = emat[m][i][j];
 364                     }
 365                 }
 366                 emat[egTotal][i][j] /= nenergy;
 367                 emat[egTotal][j][i]  = emat[egTotal][i][j];
 368             }
 369         }
 370         if (bFree)
 371         {
 372             if (bRef)
 373             {
 374                 fprintf(stderr, "Will read reference energies from inputfile\n");
 375                 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
 376                 fprintf(stderr, "Read %d reference energies\n", neref);
 377                 snew(eref, neref);
 378                 snew(erefres, neref);
 379                 for (i = 0; (i < neref); i++)
 380                 {
 381                     snew(erefres[i], 5);
 382                     sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
 383                     eref[i] = edum;
 384                 }
 385             }
 386             snew(eaver, ngroups);
 387             for (i = 0; (i < ngroups); i++)
 388             {
 389                 for (k = 0; (k < nenergy); k++)
 390                 {
 391                     eaver[i] += e[i][k];
 392                 }
 393                 eaver[i] /= nenergy;
 394             }
 395             beta = 1.0/(BOLTZ*reftemp);
 396             snew(efree, ngroups);
 397             snew(edif, ngroups);
 398             for (i = 0; (i < ngroups); i++)
 399             {
 400                 expE = 0;
 401                 for (k = 0; (k < nenergy); k++)
 402                 {
 403                     expE += std::exp(beta*(e[i][k]-eaver[i]));
 404                 }
 405                 efree[i] = std::log(expE/nenergy)/beta + eaver[i];
 406                 if (bRef)
 407                 {
 408                     n = search_str2(neref, erefres, groups[i]);
 409                     if (n != -1)
 410                     {
 411                         edif[i] = efree[i]-eref[n];
 412                     }
 413                     else
 414                     {
 415                         edif[i] = efree[i];
 416                         fprintf(stderr, "WARNING: group %s not found "
 417                                 "in reference energies.\n", groups[i]);
 418                     }
 419                 }
 420                 else
 421                 {
 422                     edif[i] = 0;
 423                 }
 424             }
 425         }
 426
 427         emid             = 0.0; /*(emin+emax)*0.5;*/
 428         egrp_nm[egTotal] = "total";
 429         for (m = 0; (m < egNR+egSP); m++)
 430         {
 431             if (egrp_use[m])
 432             {
 433                 emin = 1e10;
 434                 emax = -1e10;
 435                 for (i = 0; (i < ngroups); i++)
 436                 {
 437                     for (j = i; (j < ngroups); j++)
 438                     {
 439                         if (emat[m][i][j] > emax)
 440                         {
 441                             emax = emat[m][i][j];
 442                         }
 443                         else if (emat[m][i][j] < emin)
 444                         {
 445                             emin = emat[m][i][j];
 446                         }
 447                     }
 448                 }
 449                 if (emax == emin)
 450                 {
 451                     fprintf(stderr, "Matrix of %s energy is uniform at %f "
 452                             "(will not produce output).\n", egrp_nm[m], emax);
 453                 }
 454                 else
 455                 {
 456                     fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
 457                             egrp_nm[m], emin, emax);
 458                     if ((bCutmax) || (emax > cutmax))
 459                     {
 460                         emax = cutmax;
 461                     }
 462                     if ((bCutmin) || (emin < cutmin))
 463                     {
 464                         emin = cutmin;
 465                     }
 466                     if ((emax == cutmax) || (emin == cutmin))
 467                     {
 468                         fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
 469                     }
 470
 471                     sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
 472                     sprintf(label, "%s Interaction Energies", egrp_nm[m]);
 473                     out = gmx_ffopen(fn, "w");
 474                     if (emin >= emid)
 475                     {
 476                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 477                                   "Residue Index", "Residue Index",
 478                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 479                                   emid, emax, rmid, rhi, &nlevels);
 480                     }
 481                     else if (emax <= emid)
 482                     {
 483                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 484                                   "Residue Index", "Residue Index",
 485                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 486                                   emin, emid, rlo, rmid, &nlevels);
 487                     }
 488                     else
 489                     {
 490                         write_xpm3(out, 0, label, "Energy (kJ/mol)",
 491                                    "Residue Index", "Residue Index",
 492                                    ngroups, ngroups, groupnr, groupnr, emat[m],
 493                                    emin, emid, emax, rlo, rmid, rhi, &nlevels);
 494                     }
 495                     gmx_ffclose(out);
 496                 }
 497             }
 498         }
 499         snew(etot, egNR+egSP);
 500         for (m = 0; (m < egNR+egSP); m++)
 501         {
 502             snew(etot[m], ngroups);
 503             for (i = 0; (i < ngroups); i++)
 504             {
 505                 for (j = 0; (j < ngroups); j++)
 506                 {
 507                     etot[m][i] += emat[m][i][j];
 508                 }
 509             }
 510         }
 511
 512         out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
 513                        oenv);
 514         xvgr_legend(out, 0, NULL, oenv);
 515         j = 0;
 516         if (output_env_get_print_xvgr_codes(oenv))
 517         {
 518             char str1[STRLEN], str2[STRLEN];
 519             if (output_env_get_xvg_format(oenv) == exvgXMGR)
 520             {
 521                 sprintf(str1, "@ legend string ");
 522                 sprintf(str2, " ");
 523             }
 524             else
 525             {
 526                 sprintf(str1, "@ s");
 527                 sprintf(str2, " legend ");
 528             }
 529
 530             for (m = 0; (m < egNR+egSP); m++)
 531             {
 532                 if (egrp_use[m])
 533                 {
 534                     fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
 535                 }
 536             }
 537             if (bFree)
 538             {
 539                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
 540             }
 541             if (bFree)
 542             {
 543                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
 544             }
 545             fprintf(out, "@TYPE xy\n");
 546             fprintf(out, "#%3s", "grp");
 547
 548             for (m = 0; (m < egNR+egSP); m++)
 549             {
 550                 if (egrp_use[m])
 551                 {
 552                     fprintf(out, " %9s", egrp_nm[m]);
 553                 }
 554             }
 555             if (bFree)
 556             {
 557                 fprintf(out, " %9s", "Free");
 558             }
 559             if (bFree)
 560             {
 561                 fprintf(out, " %9s", "Diff");
 562             }
 563             fprintf(out, "\n");
 564         }
 565         for (i = 0; (i < ngroups); i++)
 566         {
 567             fprintf(out, "%3.0f", groupnr[i]);
 568             for (m = 0; (m < egNR+egSP); m++)
 569             {
 570                 if (egrp_use[m])
 571                 {
 572                     fprintf(out, " %9.5g", etot[m][i]);
 573                 }
 574             }
 575             if (bFree)
 576             {
 577                 fprintf(out, " %9.5g", efree[i]);
 578             }
 579             if (bRef)
 580             {
 581                 fprintf(out, " %9.5g", edif[i]);
 582             }
 583             fprintf(out, "\n");
 584         }
 585         xvgrclose(out);
 586     }
 587     else
 588     {
 589         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
 590                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
 591 /*
 592     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 593     n=0;
 594     emin=emax=0.0;
 595     for (k=0; (k<nenergy); k++) {
 596       for (i=0; (i<ngroups); i++)
 597     for (j=i+1; (j<ngroups); j++)
 598       emat[i][j]=eneset[n][k];
 599       sprintf(label,"t=%.0f ps",time[k]);
 600       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 601       n++;
 602     }
 603     gmx_ffclose(out);
 604  */
 605     }
 606     close_enx(in);
 607
 608     return 0;
 609 }