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,2017,2018, 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/mdlib/mdebin.h"
  53 #include "gromacs/mdtypes/enerdata.h"
  54 #include "gromacs/mdtypes/md_enums.h"
  55 #include "gromacs/trajectory/energyframe.h"
  56 #include "gromacs/utility/arraysize.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/futil.h"
  60 #include "gromacs/utility/smalloc.h"
  61 #include "gromacs/utility/strdb.h"
  62
  63
  64 static int search_str2(int nstr, char **str, char *key)
  65 {
  66     int  i, n;
  67     int  keylen = std::strlen(key);
  68     /* Linear search */
  69     n = 0;
  70     while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
  71     {
  72         n++;
  73     }
  74     for (i = 0; (i < nstr); i++)
  75     {
  76         if (gmx_strncasecmp(str[i], key, n) == 0)
  77         {
  78             return i;
  79         }
  80     }
  81
  82     return -1;
  83 }
  84
  85 int gmx_enemat(int argc, char *argv[])
  86 {
  87     const char     *desc[] = {
  88         "[THISMODULE] extracts an energy matrix from the energy file ([TT]-f[tt]).",
  89         "With [TT]-groups[tt] a file must be supplied with on each",
  90         "line a group of atoms to be used. For these groups matrix of",
  91         "interaction energies will be extracted from the energy file",
  92         "by looking for energy groups with names corresponding to pairs",
  93         "of groups of atoms, e.g. if your [TT]-groups[tt] file contains::",
  94         "",
  95         "    2",
  96         "    Protein",
  97         "    SOL",
  98         "",
  99         "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
 100         "'LJ:Protein-SOL' are expected in the energy file (although",
 101         "[THISMODULE] is most useful if many groups are analyzed",
 102         "simultaneously). Matrices for different energy types are written",
 103         "out separately, as controlled by the",
 104         "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
 105         "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
 106         "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
 107         "Finally, the total interaction energy energy per group can be ",
 108         "calculated ([TT]-etot[tt]).[PAR]",
 109
 110         "An approximation of the free energy can be calculated using:",
 111         "[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]'",
 112         "stands for time-average. A file with reference free energies",
 113         "can be supplied to calculate the free energy difference",
 114         "with some reference state. Group names (e.g. residue names)",
 115         "in the reference file should correspond to the group names",
 116         "as used in the [TT]-groups[tt] file, but a appended number",
 117         "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
 118         "in the comparison."
 119     };
 120     static gmx_bool bSum      = FALSE;
 121     static gmx_bool bMeanEmtx = TRUE;
 122     static int      skip      = 0, nlevels = 20;
 123     static real     cutmax    = 1e20, cutmin = -1e20, reftemp = 300.0;
 124     static gmx_bool bCoulSR   = TRUE, bCoul14 = FALSE;
 125     static gmx_bool bLJSR     = TRUE, bLJ14 = FALSE, bBhamSR = FALSE,
 126                     bFree     = TRUE;
 127     t_pargs         pa[]      = {
 128         { "-sum",  FALSE, etBOOL, {&bSum},
 129           "Sum the energy terms selected rather than display them all" },
 130         { "-skip", FALSE, etINT,  {&skip},
 131           "Skip number of frames between data points" },
 132         { "-mean", FALSE, etBOOL, {&bMeanEmtx},
 133           "with [TT]-groups[tt] extracts matrix of mean energies instead of "
 134           "matrix for each timestep" },
 135         { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
 136         { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
 137         { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
 138         { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
 139         { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
 140         { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
 141         { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
 142         { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
 143         { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
 144         { "-temp", FALSE, etREAL, {&reftemp},
 145           "reference temperature for free energy calculation"}
 146     };
 147     /* We will define egSP more energy-groups:
 148        egTotal (total energy) */
 149 #define egTotal egNR
 150 #define egSP 1
 151     gmx_bool          egrp_use[egNR+egSP];
 152     ener_file_t       in;
 153     FILE             *out;
 154     int               timecheck = 0;
 155     gmx_enxnm_t      *enm       = nullptr;
 156     t_enxframe       *fr;
 157     int               teller = 0;
 158     real              sum;
 159     gmx_bool          bCont, bRef;
 160     gmx_bool          bCutmax, bCutmin;
 161     real            **eneset, *time = nullptr;
 162     int              *set, i, j, prevk, k, m = 0, n, nre, nset, nenergy;
 163     char            **groups = nullptr;
 164     char              groupname[255], fn[255];
 165     int               ngroups;
 166     t_rgb             rlo, rhi, rmid;
 167     real              emax, emid, emin;
 168     real           ***emat, **etot, *groupnr;
 169     double            beta, expE, **e, *eaver, *efree = nullptr, edum;
 170     char              label[234];
 171     char            **ereflines, **erefres = nullptr;
 172     real             *eref  = nullptr, *edif = nullptr;
 173     int               neref = 0;
 174     gmx_output_env_t *oenv;
 175
 176     t_filenm          fnm[] = {
 177         { efEDR, "-f", nullptr, ffOPTRD },
 178         { efDAT, "-groups", "groups", ffREAD },
 179         { efDAT, "-eref",   "eref",   ffOPTRD },
 180         { efXPM, "-emat",   "emat",   ffWRITE },
 181         { efXVG, "-etot",   "energy", ffWRITE }
 182     };
 183 #define NFILE asize(fnm)
 184
 185     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 186                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 187     {
 188         return 0;
 189     }
 190
 191     for (i = 0; (i < egNR+egSP); i++)
 192     {
 193         egrp_use[i] = FALSE;
 194     }
 195     egrp_use[egCOULSR] = bCoulSR;
 196     egrp_use[egLJSR]   = bLJSR;
 197     egrp_use[egBHAMSR] = bBhamSR;
 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, static_cast<size_t>(gmx::square(ngroups)*egNR/2));
 223     n     = 0;
 224     prevk = 0;
 225     for (i = 0; (i < ngroups); i++)
 226     {
 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                     bool foundMatch = false;
 235                     for (k = prevk; (k < prevk+nre); k++)
 236                     {
 237                         if (std::strcmp(enm[k % nre].name, groupname) == 0)
 238                         {
 239                             set[n++]   = k;
 240                             foundMatch = true;
 241                             break;
 242                         }
 243                     }
 244                     if (!foundMatch)
 245                     {
 246                         fprintf(stderr, "WARNING! could not find group %s (%d,%d) "
 247                                 "in energy file\n", groupname, i, j);
 248                     }
 249                     else
 250                     {
 251                         prevk = k;
 252                     }
 253                 }
 254             }
 255         }
 256     }
 257     fprintf(stderr, "\n");
 258     if (n == 0)
 259     {
 260         // Return an error, can't do what the user asked for
 261         fprintf(stderr,
 262                 "None of the specified energy groups were found in this .edr file.\n"
 263                 "Perhaps you used the wrong groups, the wrong files, or didn't use a .tpr\n"
 264                 "that was made from an .mdp file that specified these energy groups.\n");
 265         return 1;
 266     }
 267     nset = n;
 268     snew(eneset, nset+1);
 269     fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
 270
 271     /* Start reading energy frames */
 272     snew(fr, 1);
 273     do
 274     {
 275         do
 276         {
 277             bCont = do_enx(in, fr);
 278             if (bCont)
 279             {
 280                 timecheck = check_times(fr->t);
 281             }
 282         }
 283         while (bCont && (timecheck < 0));
 284
 285         if (timecheck == 0)
 286         {
 287 #define DONTSKIP(cnt) (skip) ? (((cnt) % skip) == 0) : TRUE
 288
 289             if (bCont)
 290             {
 291                 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
 292                 fflush(stderr);
 293
 294                 if ((nenergy % 1000) == 0)
 295                 {
 296                     srenew(time, nenergy+1000);
 297                     for (i = 0; (i <= nset); i++)
 298                     {
 299                         srenew(eneset[i], nenergy+1000);
 300                     }
 301                 }
 302                 time[nenergy] = fr->t;
 303                 sum           = 0;
 304                 for (i = 0; (i < nset); i++)
 305                 {
 306                     eneset[i][nenergy] = fr->ener[set[i]].e;
 307                     sum               += fr->ener[set[i]].e;
 308                 }
 309                 if (bSum)
 310                 {
 311                     eneset[nset][nenergy] = sum;
 312                 }
 313                 nenergy++;
 314             }
 315             teller++;
 316         }
 317     }
 318     while (bCont && (timecheck == 0));
 319
 320     fprintf(stderr, "\n");
 321
 322     fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
 323             "over %d frames\n", ngroups, nset, nenergy);
 324
 325     snew(emat, egNR+egSP);
 326     for (j = 0; (j < egNR+egSP); j++)
 327     {
 328         if (egrp_use[m])
 329         {
 330             snew(emat[j], ngroups);
 331             for (i = 0; (i < ngroups); i++)
 332             {
 333                 snew(emat[j][i], ngroups);
 334             }
 335         }
 336     }
 337     snew(groupnr, ngroups);
 338     for (i = 0; (i < ngroups); i++)
 339     {
 340         groupnr[i] = i+1;
 341     }
 342     rlo.r  = 1.0, rlo.g  = 0.0, rlo.b  = 0.0;
 343     rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
 344     rhi.r  = 0.0, rhi.g  = 0.0, rhi.b  = 1.0;
 345     if (bMeanEmtx)
 346     {
 347         snew(e, ngroups);
 348         for (i = 0; (i < ngroups); i++)
 349         {
 350             snew(e[i], nenergy);
 351         }
 352         n = 0;
 353         for (i = 0; (i < ngroups); i++)
 354         {
 355             for (j = i; (j < ngroups); j++)
 356             {
 357                 for (m = 0; (m < egNR); m++)
 358                 {
 359                     if (egrp_use[m])
 360                     {
 361                         for (k = 0; (k < nenergy); k++)
 362                         {
 363                             emat[m][i][j] += eneset[n][k];
 364                             e[i][k]       += eneset[n][k]; /* *0.5; */
 365                             e[j][k]       += eneset[n][k]; /* *0.5; */
 366                         }
 367                         n++;
 368                         emat[egTotal][i][j] += emat[m][i][j];
 369                         emat[m][i][j]       /= nenergy;
 370                         emat[m][j][i]        = emat[m][i][j];
 371                     }
 372                 }
 373                 emat[egTotal][i][j] /= nenergy;
 374                 emat[egTotal][j][i]  = emat[egTotal][i][j];
 375             }
 376         }
 377         if (bFree)
 378         {
 379             if (bRef)
 380             {
 381                 fprintf(stderr, "Will read reference energies from inputfile\n");
 382                 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
 383                 fprintf(stderr, "Read %d reference energies\n", neref);
 384                 snew(eref, neref);
 385                 snew(erefres, neref);
 386                 for (i = 0; (i < neref); i++)
 387                 {
 388                     snew(erefres[i], 5);
 389                     sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
 390                     eref[i] = edum;
 391                 }
 392             }
 393             snew(eaver, ngroups);
 394             for (i = 0; (i < ngroups); i++)
 395             {
 396                 for (k = 0; (k < nenergy); k++)
 397                 {
 398                     eaver[i] += e[i][k];
 399                 }
 400                 eaver[i] /= nenergy;
 401             }
 402             beta = 1.0/(BOLTZ*reftemp);
 403             snew(efree, ngroups);
 404             snew(edif, ngroups);
 405             for (i = 0; (i < ngroups); i++)
 406             {
 407                 expE = 0;
 408                 for (k = 0; (k < nenergy); k++)
 409                 {
 410                     expE += std::exp(beta*(e[i][k]-eaver[i]));
 411                 }
 412                 efree[i] = std::log(expE/nenergy)/beta + eaver[i];
 413                 if (bRef)
 414                 {
 415                     n = search_str2(neref, erefres, groups[i]);
 416                     if (n != -1)
 417                     {
 418                         edif[i] = efree[i]-eref[n];
 419                     }
 420                     else
 421                     {
 422                         edif[i] = efree[i];
 423                         fprintf(stderr, "WARNING: group %s not found "
 424                                 "in reference energies.\n", groups[i]);
 425                     }
 426                 }
 427                 else
 428                 {
 429                     edif[i] = 0;
 430                 }
 431             }
 432         }
 433
 434         emid             = 0.0; /*(emin+emax)*0.5;*/
 435         egrp_nm[egTotal] = "total";
 436         for (m = 0; (m < egNR+egSP); m++)
 437         {
 438             if (egrp_use[m])
 439             {
 440                 emin = 1e10;
 441                 emax = -1e10;
 442                 for (i = 0; (i < ngroups); i++)
 443                 {
 444                     for (j = i; (j < ngroups); j++)
 445                     {
 446                         if (emat[m][i][j] > emax)
 447                         {
 448                             emax = emat[m][i][j];
 449                         }
 450                         else if (emat[m][i][j] < emin)
 451                         {
 452                             emin = emat[m][i][j];
 453                         }
 454                     }
 455                 }
 456                 if (emax == emin)
 457                 {
 458                     fprintf(stderr, "Matrix of %s energy is uniform at %f "
 459                             "(will not produce output).\n", egrp_nm[m], emax);
 460                 }
 461                 else
 462                 {
 463                     fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
 464                             egrp_nm[m], emin, emax);
 465                     if ((bCutmax) || (emax > cutmax))
 466                     {
 467                         emax = cutmax;
 468                     }
 469                     if ((bCutmin) || (emin < cutmin))
 470                     {
 471                         emin = cutmin;
 472                     }
 473                     if ((emax == cutmax) || (emin == cutmin))
 474                     {
 475                         fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
 476                     }
 477
 478                     sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
 479                     sprintf(label, "%s Interaction Energies", egrp_nm[m]);
 480                     out = gmx_ffopen(fn, "w");
 481                     if (emin >= emid)
 482                     {
 483                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 484                                   "Residue Index", "Residue Index",
 485                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 486                                   emid, emax, rmid, rhi, &nlevels);
 487                     }
 488                     else if (emax <= emid)
 489                     {
 490                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 491                                   "Residue Index", "Residue Index",
 492                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 493                                   emin, emid, rlo, rmid, &nlevels);
 494                     }
 495                     else
 496                     {
 497                         write_xpm3(out, 0, label, "Energy (kJ/mol)",
 498                                    "Residue Index", "Residue Index",
 499                                    ngroups, ngroups, groupnr, groupnr, emat[m],
 500                                    emin, emid, emax, rlo, rmid, rhi, &nlevels);
 501                     }
 502                     gmx_ffclose(out);
 503                 }
 504             }
 505         }
 506         snew(etot, egNR+egSP);
 507         for (m = 0; (m < egNR+egSP); m++)
 508         {
 509             snew(etot[m], ngroups);
 510             for (i = 0; (i < ngroups); i++)
 511             {
 512                 for (j = 0; (j < ngroups); j++)
 513                 {
 514                     etot[m][i] += emat[m][i][j];
 515                 }
 516             }
 517         }
 518
 519         out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
 520                        oenv);
 521         xvgr_legend(out, 0, nullptr, oenv);
 522         j = 0;
 523         if (output_env_get_print_xvgr_codes(oenv))
 524         {
 525             char str1[STRLEN], str2[STRLEN];
 526             if (output_env_get_xvg_format(oenv) == exvgXMGR)
 527             {
 528                 sprintf(str1, "@ legend string ");
 529                 sprintf(str2, " ");
 530             }
 531             else
 532             {
 533                 sprintf(str1, "@ s");
 534                 sprintf(str2, " legend ");
 535             }
 536
 537             for (m = 0; (m < egNR+egSP); m++)
 538             {
 539                 if (egrp_use[m])
 540                 {
 541                     fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
 542                 }
 543             }
 544             if (bFree)
 545             {
 546                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
 547             }
 548             if (bFree)
 549             {
 550                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
 551             }
 552             fprintf(out, "@TYPE xy\n");
 553             fprintf(out, "#%3s", "grp");
 554
 555             for (m = 0; (m < egNR+egSP); m++)
 556             {
 557                 if (egrp_use[m])
 558                 {
 559                     fprintf(out, " %9s", egrp_nm[m]);
 560                 }
 561             }
 562             if (bFree)
 563             {
 564                 fprintf(out, " %9s", "Free");
 565             }
 566             if (bFree)
 567             {
 568                 fprintf(out, " %9s", "Diff");
 569             }
 570             fprintf(out, "\n");
 571         }
 572         for (i = 0; (i < ngroups); i++)
 573         {
 574             fprintf(out, "%3.0f", groupnr[i]);
 575             for (m = 0; (m < egNR+egSP); m++)
 576             {
 577                 if (egrp_use[m])
 578                 {
 579                     fprintf(out, " %9.5g", etot[m][i]);
 580                 }
 581             }
 582             if (bFree)
 583             {
 584                 fprintf(out, " %9.5g", efree[i]);
 585             }
 586             if (bRef)
 587             {
 588                 fprintf(out, " %9.5g", edif[i]);
 589             }
 590             fprintf(out, "\n");
 591         }
 592         xvgrclose(out);
 593     }
 594     else
 595     {
 596         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
 597                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
 598 /*
 599     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 600     n=0;
 601     emin=emax=0.0;
 602     for (k=0; (k<nenergy); k++) {
 603       for (i=0; (i<ngroups); i++)
 604     for (j=i+1; (j<ngroups); j++)
 605       emat[i][j]=eneset[n][k];
 606       sprintf(label,"t=%.0f ps",time[k]);
 607       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 608       n++;
 609     }
 610     gmx_ffclose(out);
 611  */
 612     }
 613     close_enx(in);
 614
 615     return 0;
 616 }