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