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