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