src/gromacs/gmxana/gmx_enemat.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <string.h>
  39 #include <math.h>
  40
  41 #include "string2.h"
  42 #include "typedefs.h"
  43 #include "gmx_fatal.h"
  44 #include "vec.h"
  45 #include "smalloc.h"
  46 #include "gromacs/fileio/enxio.h"
  47 #include "statutil.h"
  48 #include "names.h"
  49 #include "macros.h"
  50 #include "xvgr.h"
  51 #include "gstat.h"
  52 #include "physics.h"
  53 #include "gromacs/fileio/matio.h"
  54 #include "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 = 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                     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         for (m = 0; (m < egNR+egSP); m++)
 515         {
 516             if (egrp_use[m])
 517             {
 518                 fprintf(out, "@ legend string %d \"%s\"\n", j++, egrp_nm[m]);
 519             }
 520         }
 521         if (bFree)
 522         {
 523             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Free");
 524         }
 525         if (bFree)
 526         {
 527             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Diff");
 528         }
 529         fprintf(out, "@TYPE xy\n");
 530         fprintf(out, "#%3s", "grp");
 531         for (m = 0; (m < egNR+egSP); m++)
 532         {
 533             if (egrp_use[m])
 534             {
 535                 fprintf(out, " %9s", egrp_nm[m]);
 536             }
 537         }
 538         if (bFree)
 539         {
 540             fprintf(out, " %9s", "Free");
 541         }
 542         if (bFree)
 543         {
 544             fprintf(out, " %9s", "Diff");
 545         }
 546         fprintf(out, "\n");
 547         for (i = 0; (i < ngroups); i++)
 548         {
 549             fprintf(out, "%3.0f", groupnr[i]);
 550             for (m = 0; (m < egNR+egSP); m++)
 551             {
 552                 if (egrp_use[m])
 553                 {
 554                     fprintf(out, " %9.5g", etot[m][i]);
 555                 }
 556             }
 557             if (bFree)
 558             {
 559                 fprintf(out, " %9.5g", efree[i]);
 560             }
 561             if (bRef)
 562             {
 563                 fprintf(out, " %9.5g", edif[i]);
 564             }
 565             fprintf(out, "\n");
 566         }
 567         ffclose(out);
 568     }
 569     else
 570     {
 571         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
 572                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
 573 /*
 574     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 575     n=0;
 576     emin=emax=0.0;
 577     for (k=0; (k<nenergy); k++) {
 578       for (i=0; (i<ngroups); i++)
 579     for (j=i+1; (j<ngroups); j++)
 580       emat[i][j]=eneset[n][k];
 581       sprintf(label,"t=%.0f ps",time[k]);
 582       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 583       n++;
 584     }
 585     ffclose(out);
 586  */
 587     }
 588     close_enx(in);
 589
 590     return 0;
 591 }