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