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     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     for (i = 0; (i < egNR+egSP); i++)
 184     {
 185         egrp_use[i] = FALSE;
 186     }
 187     egrp_use[egCOULSR] = bCoulSR;
 188     egrp_use[egLJSR]   = bLJSR;
 189     egrp_use[egBHAMSR] = bBhamSR;
 190     egrp_use[egCOULLR] = bCoulLR;
 191     egrp_use[egLJLR]   = bLJLR;
 192     egrp_use[egBHAMLR] = bBhamLR;
 193     egrp_use[egCOUL14] = bCoul14;
 194     egrp_use[egLJ14]   = bLJ14;
 195     egrp_use[egTotal]  = TRUE;
 196
 197     bRef = opt2bSet("-eref", NFILE, fnm);
 198     in   = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
 199     do_enxnms(in, &nre, &enm);
 200
 201     if (nre == 0)
 202     {
 203         gmx_fatal(FARGS, "No energies!\n");
 204     }
 205
 206     bCutmax = opt2parg_bSet("-max", asize(pa), pa);
 207     bCutmin = opt2parg_bSet("-min", asize(pa), pa);
 208
 209     nenergy = 0;
 210
 211     /* Read groupnames from input file and construct selection of
 212        energy groups from it*/
 213
 214     fprintf(stderr, "Will read groupnames from inputfile\n");
 215     ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
 216     fprintf(stderr, "Read %d groups\n", ngroups);
 217     snew(set, sqr(ngroups)*egNR/2);
 218     n     = 0;
 219     prevk = 0;
 220     for (i = 0; (i < ngroups); i++)
 221     {
 222         fprintf(stderr, "\rgroup %d", i);
 223         for (j = i; (j < ngroups); j++)
 224         {
 225             for (m = 0; (m < egNR); m++)
 226             {
 227                 if (egrp_use[m])
 228                 {
 229                     sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
 230 #ifdef DEBUG
 231                     fprintf(stderr, "\r%-15s %5d", groupname, n);
 232 #endif
 233                     for (k = prevk; (k < prevk+nre); k++)
 234                     {
 235                         if (strcmp(enm[k%nre].name, groupname) == 0)
 236                         {
 237                             set[n++] = k;
 238                             break;
 239                         }
 240                     }
 241                     if (k == prevk+nre)
 242                     {
 243                         fprintf(stderr, "WARNING! could not find group %s (%d,%d)"
 244                                 "in energy file\n", groupname, i, j);
 245                     }
 246                     else
 247                     {
 248                         prevk = k;
 249                     }
 250                 }
 251             }
 252         }
 253     }
 254     fprintf(stderr, "\n");
 255     nset = n;
 256     snew(eneset, nset+1);
 257     fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
 258
 259     /* Start reading energy frames */
 260     snew(fr, 1);
 261     do
 262     {
 263         do
 264         {
 265             bCont = do_enx(in, fr);
 266             if (bCont)
 267             {
 268                 timecheck = check_times(fr->t);
 269             }
 270         }
 271         while (bCont && (timecheck < 0));
 272
 273         if (timecheck == 0)
 274         {
 275 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
 276
 277             if (bCont)
 278             {
 279                 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
 280
 281                 if ((nenergy % 1000) == 0)
 282                 {
 283                     srenew(time, nenergy+1000);
 284                     for (i = 0; (i <= nset); i++)
 285                     {
 286                         srenew(eneset[i], nenergy+1000);
 287                     }
 288                 }
 289                 time[nenergy] = fr->t;
 290                 sum           = 0;
 291                 for (i = 0; (i < nset); i++)
 292                 {
 293                     eneset[i][nenergy] = fr->ener[set[i]].e;
 294                     sum               += fr->ener[set[i]].e;
 295                 }
 296                 if (bSum)
 297                 {
 298                     eneset[nset][nenergy] = sum;
 299                 }
 300                 nenergy++;
 301             }
 302             teller++;
 303         }
 304     }
 305     while (bCont && (timecheck == 0));
 306
 307     fprintf(stderr, "\n");
 308
 309     fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
 310             "over %d frames\n", ngroups, nset, nenergy);
 311
 312     snew(emat, egNR+egSP);
 313     for (j = 0; (j < egNR+egSP); j++)
 314     {
 315         if (egrp_use[m])
 316         {
 317             snew(emat[j], ngroups);
 318             for (i = 0; (i < ngroups); i++)
 319             {
 320                 snew(emat[j][i], ngroups);
 321             }
 322         }
 323     }
 324     snew(groupnr, ngroups);
 325     for (i = 0; (i < ngroups); i++)
 326     {
 327         groupnr[i] = i+1;
 328     }
 329     rlo.r  = 1.0, rlo.g  = 0.0, rlo.b  = 0.0;
 330     rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
 331     rhi.r  = 0.0, rhi.g  = 0.0, rhi.b  = 1.0;
 332     if (bMeanEmtx)
 333     {
 334         snew(e, ngroups);
 335         for (i = 0; (i < ngroups); i++)
 336         {
 337             snew(e[i], nenergy);
 338         }
 339         n = 0;
 340         for (i = 0; (i < ngroups); i++)
 341         {
 342             for (j = i; (j < ngroups); j++)
 343             {
 344                 for (m = 0; (m < egNR); m++)
 345                 {
 346                     if (egrp_use[m])
 347                     {
 348                         for (k = 0; (k < nenergy); k++)
 349                         {
 350                             emat[m][i][j] += eneset[n][k];
 351                             e[i][k]       += eneset[n][k]; /* *0.5; */
 352                             e[j][k]       += eneset[n][k]; /* *0.5; */
 353                         }
 354                         n++;
 355                         emat[egTotal][i][j] += emat[m][i][j];
 356                         emat[m][i][j]       /= nenergy;
 357                         emat[m][j][i]        = emat[m][i][j];
 358                     }
 359                 }
 360                 emat[egTotal][i][j] /= nenergy;
 361                 emat[egTotal][j][i]  = emat[egTotal][i][j];
 362             }
 363         }
 364         if (bFree)
 365         {
 366             if (bRef)
 367             {
 368                 fprintf(stderr, "Will read reference energies from inputfile\n");
 369                 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
 370                 fprintf(stderr, "Read %d reference energies\n", neref);
 371                 snew(eref, neref);
 372                 snew(erefres, neref);
 373                 for (i = 0; (i < neref); i++)
 374                 {
 375                     snew(erefres[i], 5);
 376                     sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
 377                     eref[i] = edum;
 378                 }
 379             }
 380             snew(eaver, ngroups);
 381             for (i = 0; (i < ngroups); i++)
 382             {
 383                 for (k = 0; (k < nenergy); k++)
 384                 {
 385                     eaver[i] += e[i][k];
 386                 }
 387                 eaver[i] /= nenergy;
 388             }
 389             beta = 1.0/(BOLTZ*reftemp);
 390             snew(efree, ngroups);
 391             snew(edif, ngroups);
 392             for (i = 0; (i < ngroups); i++)
 393             {
 394                 expE = 0;
 395                 for (k = 0; (k < nenergy); k++)
 396                 {
 397                     expE += exp(beta*(e[i][k]-eaver[i]));
 398                 }
 399                 efree[i] = log(expE/nenergy)/beta + eaver[i];
 400                 if (bRef)
 401                 {
 402                     n = search_str2(neref, erefres, groups[i]);
 403                     if (n != -1)
 404                     {
 405                         edif[i] = efree[i]-eref[n];
 406                     }
 407                     else
 408                     {
 409                         edif[i] = efree[i];
 410                         fprintf(stderr, "WARNING: group %s not found "
 411                                 "in reference energies.\n", groups[i]);
 412                     }
 413                 }
 414                 else
 415                 {
 416                     edif[i] = 0;
 417                 }
 418             }
 419         }
 420
 421         emid             = 0.0; /*(emin+emax)*0.5;*/
 422         egrp_nm[egTotal] = "total";
 423         for (m = 0; (m < egNR+egSP); m++)
 424         {
 425             if (egrp_use[m])
 426             {
 427                 emin = 1e10;
 428                 emax = -1e10;
 429                 for (i = 0; (i < ngroups); i++)
 430                 {
 431                     for (j = i; (j < ngroups); j++)
 432                     {
 433                         if (emat[m][i][j] > emax)
 434                         {
 435                             emax = emat[m][i][j];
 436                         }
 437                         else if (emat[m][i][j] < emin)
 438                         {
 439                             emin = emat[m][i][j];
 440                         }
 441                     }
 442                 }
 443                 if (emax == emin)
 444                 {
 445                     fprintf(stderr, "Matrix of %s energy is uniform at %f "
 446                             "(will not produce output).\n", egrp_nm[m], emax);
 447                 }
 448                 else
 449                 {
 450                     fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
 451                             egrp_nm[m], emin, emax);
 452                     if ((bCutmax) || (emax > cutmax))
 453                     {
 454                         emax = cutmax;
 455                     }
 456                     if ((bCutmin) || (emin < cutmin))
 457                     {
 458                         emin = cutmin;
 459                     }
 460                     if ((emax == cutmax) || (emin == cutmin))
 461                     {
 462                         fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
 463                     }
 464
 465                     sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
 466                     sprintf(label, "%s Interaction Energies", egrp_nm[m]);
 467                     out = ffopen(fn, "w");
 468                     if (emin >= emid)
 469                     {
 470                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 471                                   "Residue Index", "Residue Index",
 472                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 473                                   emid, emax, rmid, rhi, &nlevels);
 474                     }
 475                     else if (emax <= emid)
 476                     {
 477                         write_xpm(out, 0, label, "Energy (kJ/mol)",
 478                                   "Residue Index", "Residue Index",
 479                                   ngroups, ngroups, groupnr, groupnr, emat[m],
 480                                   emin, emid, rlo, rmid, &nlevels);
 481                     }
 482                     else
 483                     {
 484                         write_xpm3(out, 0, label, "Energy (kJ/mol)",
 485                                    "Residue Index", "Residue Index",
 486                                    ngroups, ngroups, groupnr, groupnr, emat[m],
 487                                    emin, emid, emax, rlo, rmid, rhi, &nlevels);
 488                     }
 489                     ffclose(out);
 490                 }
 491             }
 492         }
 493         snew(etot, egNR+egSP);
 494         for (m = 0; (m < egNR+egSP); m++)
 495         {
 496             snew(etot[m], ngroups);
 497             for (i = 0; (i < ngroups); i++)
 498             {
 499                 for (j = 0; (j < ngroups); j++)
 500                 {
 501                     etot[m][i] += emat[m][i][j];
 502                 }
 503             }
 504         }
 505
 506         out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
 507                        oenv);
 508         xvgr_legend(out, 0, NULL, oenv);
 509         j = 0;
 510         for (m = 0; (m < egNR+egSP); m++)
 511         {
 512             if (egrp_use[m])
 513             {
 514                 fprintf(out, "@ legend string %d \"%s\"\n", j++, egrp_nm[m]);
 515             }
 516         }
 517         if (bFree)
 518         {
 519             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Free");
 520         }
 521         if (bFree)
 522         {
 523             fprintf(out, "@ legend string %d \"%s\"\n", j++, "Diff");
 524         }
 525         fprintf(out, "@TYPE xy\n");
 526         fprintf(out, "#%3s", "grp");
 527         for (m = 0; (m < egNR+egSP); m++)
 528         {
 529             if (egrp_use[m])
 530             {
 531                 fprintf(out, " %9s", egrp_nm[m]);
 532             }
 533         }
 534         if (bFree)
 535         {
 536             fprintf(out, " %9s", "Free");
 537         }
 538         if (bFree)
 539         {
 540             fprintf(out, " %9s", "Diff");
 541         }
 542         fprintf(out, "\n");
 543         for (i = 0; (i < ngroups); i++)
 544         {
 545             fprintf(out, "%3.0f", groupnr[i]);
 546             for (m = 0; (m < egNR+egSP); m++)
 547             {
 548                 if (egrp_use[m])
 549                 {
 550                     fprintf(out, " %9.5g", etot[m][i]);
 551                 }
 552             }
 553             if (bFree)
 554             {
 555                 fprintf(out, " %9.5g", efree[i]);
 556             }
 557             if (bRef)
 558             {
 559                 fprintf(out, " %9.5g", edif[i]);
 560             }
 561             fprintf(out, "\n");
 562         }
 563         ffclose(out);
 564     }
 565     else
 566     {
 567         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
 568                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
 569 /*
 570     out=ftp2FILE(efMAT,NFILE,fnm,"w");
 571     n=0;
 572     emin=emax=0.0;
 573     for (k=0; (k<nenergy); k++) {
 574       for (i=0; (i<ngroups); i++)
 575     for (j=i+1; (j<ngroups); j++)
 576       emat[i][j]=eneset[n][k];
 577       sprintf(label,"t=%.0f ps",time[k]);
 578       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
 579       n++;
 580     }
 581     ffclose(out);
 582  */
 583     }
 584     close_enx(in);
 585
 586     return 0;
 587 }