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