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 static const char* egrp_nm[egNR + 1] = { "Coul-SR", "LJ-SR", "Buck-SR", "Coul-14", "LJ-14", nullptr };
  98
  99
 100 int gmx_enemat(int argc, char* argv[])
 101 {
 102     const char* desc[] = {
 103         "[THISMODULE] extracts an energy matrix from the energy file ([TT]-f[tt]).",
 104         "With [TT]-groups[tt] a file must be supplied with on each",
 105         "line a group of atoms to be used. For these groups matrix of",
 106         "interaction energies will be extracted from the energy file",
 107         "by looking for energy groups with names corresponding to pairs",
 108         "of groups of atoms, e.g. if your [TT]-groups[tt] file contains::",
 109         "",
 110         "    2",
 111         "    Protein",
 112         "    SOL",
 113         "",
 114         "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
 115         "'LJ:Protein-SOL' are expected in the energy file (although",
 116         "[THISMODULE] is most useful if many groups are analyzed",
 117         "simultaneously). Matrices for different energy types are written",
 118         "out separately, as controlled by the",
 119         "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
 120         "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
 121         "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
 122         "Finally, the total interaction energy energy per group can be ",
 123         "calculated ([TT]-etot[tt]).[PAR]",
 124
 125         "An approximation of the free energy can be calculated using:",
 126         "[MATH]E[SUB]free[sub] = E[SUB]0[sub] + kT ",
 127         "[LOG][CHEVRON][EXP](E-E[SUB]0[sub])/kT[exp][chevron][log][math], where ",
 128         "'[MATH][CHEVRON][chevron][math]'",
 129         "stands for time-average. A file with reference free energies",
 130         "can be supplied to calculate the free energy difference",
 131         "with some reference state. Group names (e.g. residue names)",
 132         "in the reference file should correspond to the group names",
 133         "as used in the [TT]-groups[tt] file, but a appended number",
 134         "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
 135         "in the comparison."
 136     };
 137     static gmx_bool bSum      = FALSE;
 138     static gmx_bool bMeanEmtx = TRUE;
 139     static int      skip = 0, nlevels = 20;
 140     static real     cutmax = 1e20, cutmin = -1e20, reftemp = 300.0;
 141     static gmx_bool bCoulSR = TRUE, bCoul14 = FALSE;
 142     static gmx_bool bLJSR = TRUE, bLJ14 = FALSE, bBhamSR = FALSE, bFree = TRUE;
 143     t_pargs         pa[] = {
 144         { "-sum",
 145           FALSE,
 146           etBOOL,
 147           { &bSum },
 148           "Sum the energy terms selected rather than display them all" },
 149         { "-skip", FALSE, etINT, { &skip }, "Skip number of frames between data points" },
 150         { "-mean",
 151           FALSE,
 152           etBOOL,
 153           { &bMeanEmtx },
 154           "with [TT]-groups[tt] extracts matrix of mean energies instead of "
 155           "matrix for each timestep" },
 156         { "-nlevels", FALSE, etINT, { &nlevels }, "number of levels for matrix colors" },
 157         { "-max", FALSE, etREAL, { &cutmax }, "max value for energies" },
 158         { "-min", FALSE, etREAL, { &cutmin }, "min value for energies" },
 159         { "-coulsr", FALSE, etBOOL, { &bCoulSR }, "extract Coulomb SR energies" },
 160         { "-coul14", FALSE, etBOOL, { &bCoul14 }, "extract Coulomb 1-4 energies" },
 161         { "-ljsr", FALSE, etBOOL, { &bLJSR }, "extract Lennard-Jones SR energies" },
 162         { "-lj14", FALSE, etBOOL, { &bLJ14 }, "extract Lennard-Jones 1-4 energies" },
 163         { "-bhamsr", FALSE, etBOOL, { &bBhamSR }, "extract Buckingham SR energies" },
 164         { "-free", FALSE, etBOOL, { &bFree }, "calculate free energy" },
 165         { "-temp",
 166           FALSE,
 167           etREAL,
 168           { &reftemp },
 169           "reference temperature for free energy calculation" }
 170     };
 171     /* We will define egSP more energy-groups:
 172        egTotal (total energy) */
 173 #define egTotal egNR
 174 #define egSP 1
 175     gmx_bool          egrp_use[egNR + egSP];
 176     ener_file_t       in;
 177     FILE*             out;
 178     int               timecheck = 0;
 179     gmx_enxnm_t*      enm       = nullptr;
 180     t_enxframe*       fr;
 181     int               teller = 0;
 182     real              sum;
 183     gmx_bool          bCont, bRef;
 184     gmx_bool          bCutmax, bCutmin;
 185     real **           eneset, *time = nullptr;
 186     int *             set, i, j, prevk, k, m = 0, n, nre, nset, nenergy;
 187     char**            groups = nullptr;
 188     char              groupname[255], fn[255];
 189     int               ngroups;
 190     t_rgb             rlo, rhi, rmid;
 191     real              emax, emid, emin;
 192     real ***          emat, **etot, *groupnr;
 193     double            beta, expE, **e, *eaver, *efree = nullptr, edum;
 194     char              label[234];
 195     char **           ereflines, **erefres = nullptr;
 196     real *            eref = nullptr, *edif = nullptr;
 197     int               neref = 0;
 198     gmx_output_env_t* oenv;
 199
 200     t_filenm fnm[] = { { efEDR, "-f", nullptr, ffOPTRD },
 201                        { efDAT, "-groups", "groups", ffREAD },
 202                        { efDAT, "-eref", "eref", ffOPTRD },
 203                        { efXPM, "-emat", "emat", ffWRITE },
 204                        { efXVG, "-etot", "energy", ffWRITE } };
 205 #define NFILE asize(fnm)
 206
 207     if (!parse_common_args(
 208                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 209     {
 210         return 0;
 211     }
 212
 213     for (i = 0; (i < egNR + egSP); i++)
 214     {
 215         egrp_use[i] = FALSE;
 216     }
 217     egrp_use[egCOULSR] = bCoulSR;
 218     egrp_use[egLJSR]   = bLJSR;
 219     egrp_use[egBHAMSR] = bBhamSR;
 220     egrp_use[egCOUL14] = bCoul14;
 221     egrp_use[egLJ14]   = bLJ14;
 222     egrp_use[egTotal]  = TRUE;
 223
 224     bRef = opt2bSet("-eref", NFILE, fnm);
 225     in   = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
 226     do_enxnms(in, &nre, &enm);
 227
 228     if (nre == 0)
 229     {
 230         gmx_fatal(FARGS, "No energies!\n");
 231     }
 232
 233     bCutmax = opt2parg_bSet("-max", asize(pa), pa);
 234     bCutmin = opt2parg_bSet("-min", asize(pa), pa);
 235
 236     nenergy = 0;
 237
 238     /* Read groupnames from input file and construct selection of
 239        energy groups from it*/
 240
 241     fprintf(stderr, "Will read groupnames from inputfile\n");
 242     ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
 243     fprintf(stderr, "Read %d groups\n", ngroups);
 244     snew(set, static_cast<size_t>(gmx::square(ngroups) * egNR / 2));
 245     n     = 0;
 246     prevk = 0;
 247     for (i = 0; (i < ngroups); i++)
 248     {
 249         for (j = i; (j < ngroups); j++)
 250         {
 251             for (m = 0; (m < egNR); m++)
 252             {
 253                 if (egrp_use[m])
 254                 {
 255                     sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
 256                     bool foundMatch = false;
 257                     for (k = prevk; (k < prevk + nre); k++)
 258                     {
 259                         if (std::strcmp(enm[k % nre].name, groupname) == 0)
 260                         {
 261                             set[n++]   = k;
 262                             foundMatch = true;
 263                             break;
 264                         }
 265                     }
 266                     if (!foundMatch)
 267                     {
 268                         fprintf(stderr,
 269                                 "WARNING! could not find group %s (%d,%d) "
 270                                 "in energy file\n",
 271                                 groupname,
 272                                 i,
 273                                 j);
 274                     }
 275                     else
 276                     {
 277                         prevk = k;
 278                     }
 279                 }
 280             }
 281         }
 282     }
 283     fprintf(stderr, "\n");
 284     if (n == 0)
 285     {
 286         // Return an error, can't do what the user asked for
 287         fprintf(stderr,
 288                 "None of the specified energy groups were found in this .edr file.\n"
 289                 "Perhaps you used the wrong groups, the wrong files, or didn't use a .tpr\n"
 290                 "that was made from an .mdp file that specified these energy groups.\n");
 291         return 1;
 292     }
 293     nset = n;
 294     snew(eneset, nset + 1);
 295     fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
 296
 297     /* Start reading energy frames */
 298     snew(fr, 1);
 299     do
 300     {
 301         do
 302         {
 303             bCont = do_enx(in, fr);
 304             if (bCont)
 305             {
 306                 timecheck = check_times(fr->t);
 307             }
 308         } while (bCont && (timecheck < 0));
 309
 310         if (timecheck == 0)
 311         {
 312 #define DONTSKIP(cnt) (skip) ? (((cnt) % skip) == 0) : TRUE
 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 }