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