src/gromacs/gmxana/gmx_enemat.c

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