2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/fileio/enxio.h"
44 #include "gromacs/fileio/matio.h"
45 #include "gromacs/fileio/strdb.h"
46 #include "gromacs/fileio/trxio.h"
47 #include "gromacs/fileio/xvgr.h"
48 #include "gromacs/gmxana/gmx_ana.h"
49 #include "gromacs/gmxana/gstat.h"
50 #include "gromacs/legacyheaders/macros.h"
51 #include "gromacs/legacyheaders/names.h"
52 #include "gromacs/legacyheaders/typedefs.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
61 static int search_str2(int nstr, char **str, char *key)
64 int keylen = std::strlen(key);
67 while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
71 for (i = 0; (i < nstr); i++)
73 if (gmx_strncasecmp(str[i], key, n) == 0)
82 int gmx_enemat(int argc, char *argv[])
84 const char *desc[] = {
85 "[THISMODULE] extracts an energy matrix from the energy file ([TT]-f[tt]).",
86 "With [TT]-groups[tt] a file must be supplied with on each",
87 "line a group of atoms to be used. For these groups matrix of",
88 "interaction energies will be extracted from the energy file",
89 "by looking for energy groups with names corresponding to pairs",
90 "of groups of atoms, e.g. if your [TT]-groups[tt] file contains::",
96 "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
97 "'LJ:Protein-SOL' are expected in the energy file (although",
98 "[THISMODULE] is most useful if many groups are analyzed",
99 "simultaneously). Matrices for different energy types are written",
100 "out separately, as controlled by the",
101 "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
102 "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
103 "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
104 "Finally, the total interaction energy energy per group can be ",
105 "calculated ([TT]-etot[tt]).[PAR]",
107 "An approximation of the free energy can be calculated using:",
108 "[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]'",
109 "stands for time-average. A file with reference free energies",
110 "can be supplied to calculate the free energy difference",
111 "with some reference state. Group names (e.g. residue names)",
112 "in the reference file should correspond to the group names",
113 "as used in the [TT]-groups[tt] file, but a appended number",
114 "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
117 static gmx_bool bSum = FALSE;
118 static gmx_bool bMeanEmtx = TRUE;
119 static int skip = 0, nlevels = 20;
120 static real cutmax = 1e20, cutmin = -1e20, reftemp = 300.0;
121 static gmx_bool bCoulSR = TRUE, bCoulLR = FALSE, bCoul14 = FALSE;
122 static gmx_bool bLJSR = TRUE, bLJLR = FALSE, bLJ14 = FALSE, bBhamSR = FALSE, bBhamLR = FALSE,
125 { "-sum", FALSE, etBOOL, {&bSum},
126 "Sum the energy terms selected rather than display them all" },
127 { "-skip", FALSE, etINT, {&skip},
128 "Skip number of frames between data points" },
129 { "-mean", FALSE, etBOOL, {&bMeanEmtx},
130 "with [TT]-groups[tt] extracts matrix of mean energies instead of "
131 "matrix for each timestep" },
132 { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
133 { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
134 { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
135 { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
136 { "-coullr", FALSE, etBOOL, {&bCoulLR}, "extract Coulomb LR energies"},
137 { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
138 { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
139 { "-ljlr", FALSE, etBOOL, {&bLJLR}, "extract Lennard-Jones LR energies"},
140 { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
141 { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
142 { "-bhamlr", FALSE, etBOOL, {&bBhamLR}, "extract Buckingham LR energies"},
143 { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
144 { "-temp", FALSE, etREAL, {&reftemp},
145 "reference temperature for free energy calculation"}
147 /* We will define egSP more energy-groups:
148 egTotal (total energy) */
151 gmx_bool egrp_use[egNR+egSP];
155 gmx_enxnm_t *enm = NULL;
159 gmx_bool bCont, bRef;
160 gmx_bool bCutmax, bCutmin;
161 real **eneset, *time = NULL;
162 int *set, i, j, k, prevk, m = 0, n, nre, nset, nenergy;
163 char **groups = NULL;
164 char groupname[255], fn[255];
166 t_rgb rlo, rhi, rmid;
167 real emax, emid, emin;
168 real ***emat, **etot, *groupnr;
169 double beta, expE, **e, *eaver, *efree = NULL, edum;
171 char **ereflines, **erefres = NULL;
172 real *eref = NULL, *edif = NULL;
177 { efEDR, "-f", NULL, ffOPTRD },
178 { efDAT, "-groups", "groups", ffREAD },
179 { efDAT, "-eref", "eref", ffOPTRD },
180 { efXPM, "-emat", "emat", ffWRITE },
181 { efXVG, "-etot", "energy", ffWRITE }
183 #define NFILE asize(fnm)
185 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
186 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
191 for (i = 0; (i < egNR+egSP); i++)
195 egrp_use[egCOULSR] = bCoulSR;
196 egrp_use[egLJSR] = bLJSR;
197 egrp_use[egBHAMSR] = bBhamSR;
198 egrp_use[egCOULLR] = bCoulLR;
199 egrp_use[egLJLR] = bLJLR;
200 egrp_use[egBHAMLR] = bBhamLR;
201 egrp_use[egCOUL14] = bCoul14;
202 egrp_use[egLJ14] = bLJ14;
203 egrp_use[egTotal] = TRUE;
205 bRef = opt2bSet("-eref", NFILE, fnm);
206 in = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
207 do_enxnms(in, &nre, &enm);
211 gmx_fatal(FARGS, "No energies!\n");
214 bCutmax = opt2parg_bSet("-max", asize(pa), pa);
215 bCutmin = opt2parg_bSet("-min", asize(pa), pa);
219 /* Read groupnames from input file and construct selection of
220 energy groups from it*/
222 fprintf(stderr, "Will read groupnames from inputfile\n");
223 ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
224 fprintf(stderr, "Read %d groups\n", ngroups);
225 snew(set, static_cast<size_t>(sqr(ngroups)*egNR/2));
228 for (i = 0; (i < ngroups); i++)
230 fprintf(stderr, "\rgroup %d", i);
231 for (j = i; (j < ngroups); j++)
233 for (m = 0; (m < egNR); m++)
237 sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
239 fprintf(stderr, "\r%-15s %5d", groupname, n);
241 for (k = prevk; (k < prevk+nre); k++)
243 if (std::strcmp(enm[k%nre].name, groupname) == 0)
251 fprintf(stderr, "WARNING! could not find group %s (%d,%d)"
252 "in energy file\n", groupname, i, j);
262 fprintf(stderr, "\n");
264 snew(eneset, nset+1);
265 fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
267 /* Start reading energy frames */
273 bCont = do_enx(in, fr);
276 timecheck = check_times(fr->t);
279 while (bCont && (timecheck < 0));
283 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
287 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
289 if ((nenergy % 1000) == 0)
291 srenew(time, nenergy+1000);
292 for (i = 0; (i <= nset); i++)
294 srenew(eneset[i], nenergy+1000);
297 time[nenergy] = fr->t;
299 for (i = 0; (i < nset); i++)
301 eneset[i][nenergy] = fr->ener[set[i]].e;
302 sum += fr->ener[set[i]].e;
306 eneset[nset][nenergy] = sum;
313 while (bCont && (timecheck == 0));
315 fprintf(stderr, "\n");
317 fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
318 "over %d frames\n", ngroups, nset, nenergy);
320 snew(emat, egNR+egSP);
321 for (j = 0; (j < egNR+egSP); j++)
325 snew(emat[j], ngroups);
326 for (i = 0; (i < ngroups); i++)
328 snew(emat[j][i], ngroups);
332 snew(groupnr, ngroups);
333 for (i = 0; (i < ngroups); i++)
337 rlo.r = 1.0, rlo.g = 0.0, rlo.b = 0.0;
338 rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
339 rhi.r = 0.0, rhi.g = 0.0, rhi.b = 1.0;
343 for (i = 0; (i < ngroups); i++)
348 for (i = 0; (i < ngroups); i++)
350 for (j = i; (j < ngroups); j++)
352 for (m = 0; (m < egNR); m++)
356 for (k = 0; (k < nenergy); k++)
358 emat[m][i][j] += eneset[n][k];
359 e[i][k] += eneset[n][k]; /* *0.5; */
360 e[j][k] += eneset[n][k]; /* *0.5; */
363 emat[egTotal][i][j] += emat[m][i][j];
364 emat[m][i][j] /= nenergy;
365 emat[m][j][i] = emat[m][i][j];
368 emat[egTotal][i][j] /= nenergy;
369 emat[egTotal][j][i] = emat[egTotal][i][j];
376 fprintf(stderr, "Will read reference energies from inputfile\n");
377 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
378 fprintf(stderr, "Read %d reference energies\n", neref);
380 snew(erefres, neref);
381 for (i = 0; (i < neref); i++)
384 sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
388 snew(eaver, ngroups);
389 for (i = 0; (i < ngroups); i++)
391 for (k = 0; (k < nenergy); k++)
397 beta = 1.0/(BOLTZ*reftemp);
398 snew(efree, ngroups);
400 for (i = 0; (i < ngroups); i++)
403 for (k = 0; (k < nenergy); k++)
405 expE += std::exp(beta*(e[i][k]-eaver[i]));
407 efree[i] = std::log(expE/nenergy)/beta + eaver[i];
410 n = search_str2(neref, erefres, groups[i]);
413 edif[i] = efree[i]-eref[n];
418 fprintf(stderr, "WARNING: group %s not found "
419 "in reference energies.\n", groups[i]);
429 emid = 0.0; /*(emin+emax)*0.5;*/
430 egrp_nm[egTotal] = "total";
431 for (m = 0; (m < egNR+egSP); m++)
437 for (i = 0; (i < ngroups); i++)
439 for (j = i; (j < ngroups); j++)
441 if (emat[m][i][j] > emax)
443 emax = emat[m][i][j];
445 else if (emat[m][i][j] < emin)
447 emin = emat[m][i][j];
453 fprintf(stderr, "Matrix of %s energy is uniform at %f "
454 "(will not produce output).\n", egrp_nm[m], emax);
458 fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
459 egrp_nm[m], emin, emax);
460 if ((bCutmax) || (emax > cutmax))
464 if ((bCutmin) || (emin < cutmin))
468 if ((emax == cutmax) || (emin == cutmin))
470 fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
473 sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
474 sprintf(label, "%s Interaction Energies", egrp_nm[m]);
475 out = gmx_ffopen(fn, "w");
478 write_xpm(out, 0, label, "Energy (kJ/mol)",
479 "Residue Index", "Residue Index",
480 ngroups, ngroups, groupnr, groupnr, emat[m],
481 emid, emax, rmid, rhi, &nlevels);
483 else if (emax <= emid)
485 write_xpm(out, 0, label, "Energy (kJ/mol)",
486 "Residue Index", "Residue Index",
487 ngroups, ngroups, groupnr, groupnr, emat[m],
488 emin, emid, rlo, rmid, &nlevels);
492 write_xpm3(out, 0, label, "Energy (kJ/mol)",
493 "Residue Index", "Residue Index",
494 ngroups, ngroups, groupnr, groupnr, emat[m],
495 emin, emid, emax, rlo, rmid, rhi, &nlevels);
501 snew(etot, egNR+egSP);
502 for (m = 0; (m < egNR+egSP); m++)
504 snew(etot[m], ngroups);
505 for (i = 0; (i < ngroups); i++)
507 for (j = 0; (j < ngroups); j++)
509 etot[m][i] += emat[m][i][j];
514 out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
516 xvgr_legend(out, 0, NULL, oenv);
518 if (output_env_get_print_xvgr_codes(oenv))
520 char str1[STRLEN], str2[STRLEN];
521 if (output_env_get_xvg_format(oenv) == exvgXMGR)
523 sprintf(str1, "@ legend string ");
528 sprintf(str1, "@ s");
529 sprintf(str2, " legend ");
532 for (m = 0; (m < egNR+egSP); m++)
536 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
541 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
545 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
547 fprintf(out, "@TYPE xy\n");
548 fprintf(out, "#%3s", "grp");
550 for (m = 0; (m < egNR+egSP); m++)
554 fprintf(out, " %9s", egrp_nm[m]);
559 fprintf(out, " %9s", "Free");
563 fprintf(out, " %9s", "Diff");
567 for (i = 0; (i < ngroups); i++)
569 fprintf(out, "%3.0f", groupnr[i]);
570 for (m = 0; (m < egNR+egSP); m++)
574 fprintf(out, " %9.5g", etot[m][i]);
579 fprintf(out, " %9.5g", efree[i]);
583 fprintf(out, " %9.5g", edif[i]);
591 fprintf(stderr, "While typing at your keyboard, suddenly...\n"
592 "...nothing happens.\nWARNING: Not Implemented Yet\n");
594 out=ftp2FILE(efMAT,NFILE,fnm,"w");
597 for (k=0; (k<nenergy); k++) {
598 for (i=0; (i<ngroups); i++)
599 for (j=i+1; (j<ngroups); j++)
600 emat[i][j]=eneset[n][k];
601 sprintf(label,"t=%.0f ps",time[k]);
602 write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);